← Back to branch summary

~airpollution/fluidity/fluidity_airpollution

« back to all changes in this revision

Viewing changes to libadaptivity/load_balance/src/functional_2d.cpp

  • Committer: ziyouzhj
  • Date: 2013-12-09 16:51:29 UTC
  • Revision ID: ziyouzhj@gmail.com-20131209165129-ucoetc3u0atyy05c
airpolution

Show diffs side-by-side

added added

removed removed

Lines of Context:
libadaptivity/load_balance/src/functional_2d.cpp
 
1
/*
 
2
 *      Copyright (c) 2004-2006 by Gerard Gorman
 
3
 *      Copyright (c) 2006- Imperial College London
 
4
 *      See COPYING file for copying and redistribution conditions.
 
5
 *
 
6
 *      This program is free software; you can redistribute it and/or modify
 
7
 *      it under the terms of the GNU General Public License as published by
 
8
 *      the Free Software Foundation; version 2 of the License.
 
9
 *
 
10
 *      This program is distributed in the hope that it will be useful,
 
11
 *      but WITHOUT ANY WARRANTY; without even the implied warranty of
 
12
 *      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
13
 *      GNU General Public License for more details.
 
14
 *
 
15
 *      Contact info: gerard.j.gorman@gmail.com/g.gorman@imperial.ac.uk
 
16
 */
 
17
 
 
18
#include "functional_2d.h"
 
19
#include <iostream>
 
20
 
 
21
using namespace std;
 
22
 
 
23
functional_2d::functional_2d(): alpha(0.5/sqrt(3.0)){}
 
24
functional_2d::~functional_2d(){}
 
25
 
 
26
samfloat_t functional_2d::dot(samfloat_t x0, samfloat_t y0,
 
27
                       samfloat_t x1, samfloat_t y1){
 
28
  return x0*(x1*m[0]+y1*m[1]) + y0*(x1*m[1]+y1*m[2]);
 
29
}
 
30
 
 
31
samfloat_t functional_2d::len2(samfloat_t x0, samfloat_t y0,
 
32
                        samfloat_t x1, samfloat_t y1){
 
33
  samfloat_t vx = x0-x1;
 
34
  samfloat_t vy = y0-y1;
 
35
  return vx*vx*m[0] + 2*vx*vy*m[1] + vy*vy*m[2];
 
36
}
 
37
 
 
38
samfloat_t functional_2d::GetInRadius(samfloat_t x0, samfloat_t y0,
 
39
                               samfloat_t x1, samfloat_t y1,
 
40
                               samfloat_t x2, samfloat_t y2){
 
41
  samfloat_t r0 = sqrt(len2(x1, y1, x2, y2));
 
42
  samfloat_t r1 = sqrt(len2(x2, y2, x0, y0));
 
43
  samfloat_t r2 = sqrt(len2(x0, y0, x1, y1));
 
44
  
 
45
  samfloat_t s=0.5*(r0+r1+r2);
 
46
  return sqrt((s-r0)*(s-r1)*(s-r2)/s);  
 
47
}
 
48
 
 
49
samfloat_t functional_2d::GetArea(samfloat_t r0_x, samfloat_t r0_y,
 
50
                                          samfloat_t r1_x, samfloat_t r1_y,
 
51
                                          samfloat_t r2_x, samfloat_t r2_y){
 
52
  samfloat_t x1 = r1_x - r0_x;
 
53
  samfloat_t y1 = r1_y - r0_y;
 
54
  
 
55
  samfloat_t x2 = r2_x - r0_x;
 
56
  samfloat_t y2 = r2_y - r0_y;
 
57
  
 
58
  return -0.5*(x1*y2 - x2*y1);
 
59
}
 
60
 
 
61
 
 
62
samfloat_t functional_2d::GetShape(samfloat_t x0, samfloat_t y0,
 
63
                            samfloat_t x1, samfloat_t y1,
 
64
                            samfloat_t x2, samfloat_t y2){
 
65
  
 
66
  samfloat_t inradius = GetInRadius(x0, y0,
 
67
                                x1, y1,
 
68
                                x2, y2);
 
69
  
 
70
  samfloat_t shape=0.0;
 
71
  if(inradius>10e-6)
 
72
    shape = (alpha/inradius - 1)*(alpha/inradius - 1);
 
73
  else
 
74
    shape = 1.0e6;
 
75
  
 
76
  return shape;
 
77
}
 
78
 
 
79
samfloat_t functional_2d::GetSize(samfloat_t x0, samfloat_t y0,
 
80
                           samfloat_t x1, samfloat_t y1,
 
81
                           samfloat_t x2, samfloat_t y2){
 
82
  samfloat_t size=0.0;
 
83
  
 
84
  // edge a
 
85
  samfloat_t r0 = sqrt(len2(x0, y0, x1, y1));
 
86
  size += ((r0-1)*(r0-1));
 
87
  
 
88
  // edge b
 
89
  samfloat_t r1 = sqrt(len2(x0, y0, x2, y2));
 
90
  size += ((r1-1)*(r1-1));
 
91
  
 
92
  // edge c
 
93
  samfloat_t r2 = sqrt(len2(x1, y1, x2, y2));
 
94
  size += ((r2-1)*(r2-1));
 
95
  
 
96
  return size;
 
97
}
 
98
 
 
99
samfloat_t functional_2d::standard(samfloat_t x0, samfloat_t y0, const samfloat_t *m0,
 
100
                                   samfloat_t x1, samfloat_t y1, const samfloat_t *m1,
 
101
                                   samfloat_t x2, samfloat_t y2, const samfloat_t *m2){
 
102
  
 
103
  m[0] = (m0[0]+m1[0]+m2[0])/3.0;
 
104
  m[1] = (m0[1]+m1[1]+m2[1])/3.0;
 
105
  m[2] = (m0[2]+m1[2]+m2[2])/3.0;
 
106
  
 
107
  samfloat_t shape = GetShape(x0, y0,
 
108
                              x1, y1,
 
109
                              x2, y2);
 
110
  
 
111
  samfloat_t size = GetSize(x0, y0,
 
112
                            x1, y1,
 
113
                            x2, y2);
 
114
  
 
115
  samfloat_t functional = 10.0*size + shape;
 
116
  if(!finite(functional)){
 
117
    return 1.0e6;
 
118
  }
 
119
 
 
120
  return functional;
 
121
}
 
122
 
 
123
samfloat_t functional_2d::oddy(samfloat_t x0, samfloat_t y0, const samfloat_t *m0,
 
124
                        samfloat_t x1, samfloat_t y1, const samfloat_t *m1,
 
125
                        samfloat_t x2, samfloat_t y2, const samfloat_t *m2){
 
126
 
 
127
  m[0] = (m0[0]+m1[0]+m2[0])/3.0;
 
128
  m[1] = (m0[1]+m1[1]+m2[1])/3.0;
 
129
  m[2] = (m0[2]+m1[2]+m2[2])/3.0;
 
130
  
 
131
  if(GetArea(x0, y0, x1, y1, x2, y2)<0.0)
 
132
    return 1.0e6;
 
133
  
 
134
  samfloat_t l2[3];
 
135
  l2[0] = len2(x1, y1, x2, y2);
 
136
  l2[1] = len2(x2, y2, x0, y0);
 
137
  l2[2] = len2(x0, y0, x1, y1);
 
138
  
 
139
  samfloat_t b2[3];
 
140
  b2[0] = dot(x2-x1, y2-y1, x0-x2, y0-y2);
 
141
  b2[1] = dot(x0-x2, y0-y2, x1-x0, y1-y0);
 
142
  b2[2] = dot(x1-x0, y1-y0, x2-x1, y2-y1);
 
143
  
 
144
  samfloat_t g[3];
 
145
  for(size_t i=0;i<3;i++){
 
146
    b2[i]*=b2[i];
 
147
    g[i] = l2[i]*l2[(i+1)%3] - b2[i];
 
148
  }
 
149
  
 
150
  samfloat_t functional = 0;
 
151
  for(size_t i=0;i<3;i++){
 
152
    samfloat_t a = (l2[i]-l2[(i+1)%3]);
 
153
    functional += (a*a + 4*b2[i])/g[i];
 
154
    // cout<<"functional "<<functional<<" "<<a<<" "<<b2[i]<<" "<<g[i]<<endl;
 
155
  }
 
156
  
 
157
  return 0.5*functional;
 
158
}