1 // GetDP - Copyright (C) 1997-2021 P. Dular and C. Geuzaine, University of Liege
2 //
3 // See the LICENSE.txt file for license information. Please report all
4 // issues on https://gitlab.onelab.info/getdp/getdp/issues.
5 //
6 // Contributor(s):
7 // Ruth Sabariego
8 //
9
10 #include <math.h>
11 #include "GetDPConfig.h"
12 #include "ProData.h"
13 #include "ProDefine.h"
14 #include "BF.h"
15 #include "GF.h"
16 #include "Message.h"
17
18 #if defined(HAVE_KERNEL)
19 #include "GeoData.h"
20 #endif
21
22 #define SQU(a) ((a)*(a))
23 #define THESIGN(a) ((a)>=0 ? 1 : -1)
24 #define ONE_OVER_TWO_PI 1.5915494309189534E-01
25 #define ONE_OVER_FOUR_PI 7.9577471545947668E-02
26
27 #define MAX_NODES 6
28 #define EPSILON 1.e-8
29 #define EPSILON2 1.e-20
30 #define RADIUS 0.154797 /* this is a hack... */
31
32 /* ------------------------------------------------------------------------ */
33 /* G F _ L a p l a c e x F o r m */
34 /* ------------------------------------------------------------------------ */
35
GF_LaplacexForm(GF_ARGX)36 void GF_LaplacexForm(GF_ARGX)
37 {
38 #if !defined(HAVE_KERNEL)
39 Message::Error("GF_LaplacexForm requires Kernel");
40 #else
41 double xs[MAX_NODES], ys[MAX_NODES], zs[MAX_NODES], u[3], v[3], n[3];
42 double u2=0., v2=0., xl=0., yl=0., zl=0., zl_2=0. ;
43 double Area, m0[3], m1[3], m2[3] ;
44 int Type_Int=0, i, j = 1 ;
45 double a=0., b=0., c=0., d, e, f, i1, I1 = 0., Iua, Iva, r2;
46 double s0m=0., s0p=0., s1m=0., s1p=0., s2m=0., s2p=0., t00, t10, t20, t0m_2, t0p_2, t1p_2;
47 double r00_2=0., r10_2=0., r20_2=0., r00, r10, r20, r0p=0., r0m=0., r1p=0.;
48 double f20=0., f21=0., f22=0., B0, B1, B2 ;
49 double f30, f31, f32, N10, N20, N30 ;
50
51 Val->Val[MAX_DIM] = 0.0 ;
52
53 switch ((int)Fct->Para[0]) {
54
55 case DIM_2D :
56
57 switch (Element->ElementSource->Type) {
58
59 case POINT_ELEMENT :
60 xs[0] = Element->ElementSource->x[0] ; ys[0] = Element->ElementSource->y[0] ;
61
62 r2 = SQU(x-xs[0])+SQU(y-ys[0]) ;
63 if (r2 > SQU(RADIUS)){
64 Val->Type = SCALAR ;
65 Val->Val[0] = - ONE_OVER_FOUR_PI * log(r2) ;
66 }
67 else{
68 Val->Type = SCALAR ;
69 Val->Val[0] = - ONE_OVER_FOUR_PI * log(SQU(RADIUS)) ;
70 }
71 break ;
72
73 case LINE :
74 xs[0] = Element->ElementSource->x[0] ; ys[0] = Element->ElementSource->y[0] ;
75 xs[1] = Element->ElementSource->x[1] ; ys[1] = Element->ElementSource->y[1] ;
76
77 if(xFunctionBF == (void(*)())BF_Volume) {
78 a = SQU(xs[0]-xs[1]) + SQU(ys[0]-ys[1]) ;
79 b = 2. * ((x-xs[0])*(xs[0]-xs[1]) + (y-ys[0])*(ys[0]-ys[1])) ;
80 c = SQU(x-xs[0]) + SQU(y-ys[0]) ;
81 d = 0.5 * b / a ;
82 e = c / a ;
83 f = e - d*d ;
84
85 if (f > EPSILON) { Type_Int = 1; }
86 else if (fabs(f) < EPSILON){ Type_Int = 0; }
87 else { Type_Int = -1; f = -f; }
88 if (Element->Num == Element->ElementSource->Num) Type_Int = 2 ;
89 if ((c == 0) || ((b == -2*a) && (c == a))) Type_Int = 3 ;
90
91 switch (Type_Int) {
92 case -1 :
93 I1 = log(a) +
94 ( (d+1.) * log(SQU(d+1.) - f) - 2.*(d+1.) +
95 sqrt(f) * log((d+1.+sqrt(f))/(d+1.-sqrt(f))) ) -
96 ( d*log(d*d-f) - 2.*d + sqrt(f)*log((d+sqrt(f))/(d-sqrt(f))) ) ;
97 break ;
98 case 0 :
99 I1 = log(a) + (d+1.)*log(SQU(d+1.)) - d*log(SQU(d)) - 2. ;
100 break ;
101 case 1 :
102 I1 = log(a) +
103 ( (d+1.) * log(SQU(d+1.) + f) - 2.*(d+1.) +
104 2.*sqrt(f) * atan((d+1.)/sqrt(f)) ) -
105 ( d*log(d*d+f) - 2.*d + 2.*sqrt(f)*atan(d/sqrt(f)) ) ;
106 break ;
107 case 2 :
108 i1 = -b / (2.*a) ;
109 I1 = 2. * i1 * (log(i1) - 1.) +
110 2. * (1.-i1) * (log(1.-i1) - 1.) + log(a) ;
111 break ;
112 case 3 :
113 I1 = .5 * log(a) - 1. ;
114 break ;
115 }
116
117 Val->Type = SCALAR ;
118 Val->Val[0] = - ONE_OVER_FOUR_PI * I1 ;
119 }
120 else {
121 Message::Error("Unknown Basis Function Type for 'GF_LaplacexForm'");
122 }
123 break ;
124
125 default :
126 Message::Error("Unknown Element Type (%s) for 'GF_LaplacexForm'",
127 Get_StringForDefine(Element_Type, Element->ElementSource->Type));
128 }
129
130 break;
131
132 case DIM_3D :
133
134 switch (Element->ElementSource->Type) {
135
136 case LINE :
137 xs[0] = Element->ElementSource->x[0] ; ys[0] = Element->ElementSource->y[0] ;
138 zs[0] = Element->ElementSource->z[0] ;
139 xs[1] = Element->ElementSource->x[1] ; ys[1] = Element->ElementSource->y[1] ;
140 zs[1] = Element->ElementSource->z[1] ;
141
142 a = SQU(xs[0]-xs[1]) + SQU(ys[0]-ys[1]) + SQU(zs[0]-zs[1]) ;
143 b = 2. * ((x-xs[0])*(xs[0]-xs[1]) +
144 (y-ys[0])*(ys[0]-ys[1]) +
145 (z-zs[0])*(zs[0]-zs[1])) ;
146 c = SQU(x-xs[0]) + SQU(y-ys[0]) + SQU(z-zs[0]) + SQU(RADIUS) ;
147
148 Val->Val[0] = ONE_OVER_FOUR_PI *
149 log( ( 2.*sqrt(a*(a+b+c))+2.*a+b ) / ( 2.*sqrt(a*c)+b ) ) ;
150 Val->Type = SCALAR ;
151 break ;
152
153 case TRIANGLE :
154 case QUADRANGLE :
155 if(xFunctionBF == (void(*)())BF_Volume) Type_Int = 1 ;
156 if(xFunctionBF == (void(*)())BF_Node) Type_Int = 2 ;
157
158 xs[0] = Element->ElementSource->x[0] ; ys[0] = Element->ElementSource->y[0] ;
159 zs[0] = Element->ElementSource->z[0] ;
160 xs[1] = Element->ElementSource->x[1] ; ys[1] = Element->ElementSource->y[1] ;
161 zs[1] = Element->ElementSource->z[1] ;
162 xs[2] = Element->ElementSource->x[2] ; ys[2] = Element->ElementSource->y[2] ;
163 zs[2] = Element->ElementSource->z[2] ;
164
165 if (Element->ElementSource->Type == QUADRANGLE) {
166 xs[3] = Element->ElementSource->x[3] ; ys[3] = Element->ElementSource->y[3] ;
167 zs[3] = Element->ElementSource->z[3] ;
168 j = 0 ;
169 };
170
171 for(i = j; i < 2; i++){
172 /* triangle side lengths */
173 a = sqrt(SQU(xs[1]-xs[0]) + SQU(ys[1]-ys[0]) + SQU(zs[1]-zs[0]));
174 b = sqrt(SQU(xs[2]-xs[1]) + SQU(ys[2]-ys[1]) + SQU(zs[2]-zs[1]));
175 c = sqrt(SQU(xs[2]-xs[0]) + SQU(ys[2]-ys[0]) + SQU(zs[2]-zs[0]));
176
177 /* local system (u,v,w) centered at (xs[0],ys[0],zs[0]) */
178 u[0] = (xs[1]-xs[0])/a;
179 u[1] = (ys[1]-ys[0])/a;
180 u[2] = (zs[1]-zs[0])/a;
181
182 /* triangle normal */
183 Geo_CreateNormal(Element->ElementSource->Type,xs,ys,zs,n);
184
185 /* v = n /\ u */
186 v[0] = n[1]*u[2]-n[2]*u[1];
187 v[1] = n[2]*u[0]-n[0]*u[2];
188 v[2] = n[0]*u[1]-n[1]*u[0];
189
190 u2 = (xs[2]-xs[0])*u[0] + (ys[2]-ys[0])*u[1] + (zs[2]-zs[0])*u[2]; /* u2 coordinate */
191 v2 = (xs[2]-xs[0])*v[0] + (ys[2]-ys[0])*v[1] + (zs[2]-zs[0])*v[2]; /* triangle height, v2 coordinate */
192
193 /* local coordinates of the observation point (xl, yl, zl) */
194 xl = u[0] * (x-xs[0]) + u[1] * (y-ys[0]) + u[2] * (z-zs[0]);
195 yl = v[0] * (x-xs[0]) + v[1] * (y-ys[0]) + v[2] * (z-zs[0]);
196 zl = n[0] * (x-xs[0]) + n[1] * (y-ys[0]) + n[2] * (z-zs[0]);
197
198 s0m = -( (a-xl) * (a-u2) + yl*v2 ) / b;
199 s0p = s0m + b;
200 s1p = ( xl * u2 + yl * v2 ) / c;
201 s1m = s1p - c;
202 s2m = - xl;
203 s2p = a - xl;
204
205 /* distance observation point projection on triangle plane to triangle local vertices*/
206 /* t1m = t0p ; t2p = t0m ; t2m = t1p ; */
207 t00 = (yl * (u2-a) + v2 * (a-xl)) / b;
208 t10 = (xl * v2 - yl * u2) / c;
209 t20 = yl;
210
211 t0m_2 = (a-xl)*(a-xl) + yl*yl;
212 t0p_2 = (u2-xl)*(u2-xl) + (v2-yl)*(v2-yl);
213 t1p_2 = xl*xl + yl*yl;
214
215 /* minimum distances^2 from the observation point to each triangle side*/
216 zl_2 = SQU(zl) ;
217 r00_2 = SQU(t00) + zl_2 ;
218 r10_2 = SQU(t10) + zl_2 ;
219 r20_2 = SQU(t20) + zl_2 ;
220
221 /* distances from observation point to the vertices*/
222 r0p = sqrt(t0p_2 + zl_2);
223 r0m = sqrt(t0m_2 + zl_2);
224 r1p = sqrt(t1p_2 + zl_2);
225
226 r00 = sqrt(r00_2);
227 r10 = sqrt(r10_2);
228 r20 = sqrt(r20_2);
229
230 /* intermediate functions */
231 if(r00 <= EPSILON*(fabs(s0m)+fabs(s0p)) ){ f20 = log(s0m/s0p) ; B0 = 0; }
232 else{
233 if (!(r0m + s0m))
234 Message::Error("1/0 in GF_LaplacexForm (case _3D TRIANGLE) Num %d Obs %.15e %.15e %.15e",
235 Element->ElementSource->Num, x, y, z) ;
236 f20 = log((r0p + s0p) / (r0m + s0m));
237 B0 = atan(t00*s0p/(r00_2+fabs(zl)*r0p))-atan(t00*s0m/(r00_2+fabs(zl)*r0m));
238 }
239
240 if(r10 <= EPSILON*(fabs(s1m)+fabs(s1p)) ){ f21 = log(s1m/s1p); B1 = 0; }
241 else{
242 if(!(r0p + s1m))
243 Message::Error("1/0 in GF_LaplacexForm (case _3D TRIANGLE) Num %d Obs %.15e %.15e %.15e",
244 Element->ElementSource->Num, x, y, z) ;
245 f21 = log((r1p + s1p) / (r0p + s1m));
246 B1 = atan(t10*s1p/(r10_2+fabs(zl)*r1p))-atan(t10*s1m/(r10_2+fabs(zl)*r0p));
247 }
248
249 if(r20 <= EPSILON*(fabs(s2m)+fabs(s2p)) ){ f22 = log(s2m/s2p); B2 = 0; }
250 else{
251 if(!(r1p+s2m))
252 Message::Error("1/0 in GF_LaplacexForm (case _3D TRIANGLE) Num %d Obs %.15e %.15e %.15e",
253 Element->ElementSource->Num, x, y, z) ;
254 f22 = log((r0m + s2p) / (r1p + s2m));
255 B2 = atan(t20*s2p/(r20_2+fabs(zl)*r0m))-atan(t20*s2m/(r20_2+fabs(zl)*r1p));
256 }
257
258 I1 += -fabs(zl)*(B0+B1+B2) + t00*f20+t10*f21+t20*f22 ; /* 1/r integral solution*/
259
260 if (j == 0){ xs[1] = xs[2]; ys[1] = ys[2]; zs[1] = zs[2]; xs[2] = xs[3]; ys[2] = ys[3]; zs[2] = zs[3];}
261 }
262
263
264 switch ( Type_Int ){
265 case 1 : /* BF_Volume */
266 Area = a * v2/2 ;/* Triangle area */
267 Val->Val[0] = I1 /Area ;
268 break;
269
270 case 2 : /* BF_Node */
271 if (!v2)
272 Message::Error("1/0 in GF_LaplacexForm (case _3D TRIANGLE) v2 %e", v2);
273
274 f30 = (s0p*r0p-s0m*r0m) + r00_2 * f20 ; /* f3i */
275 f31 = (s1p*r1p-s1m*r0p) + r10_2 * f21 ;
276 f32 = (s2p*r0m-s2m*r1p) + r20_2 * f22 ;
277
278 m0[0] = ((ys[2] - ys[1]) * n[2] - (zs[2] - zs[1]) * n[1])*f30/b ;
279 m0[1] = ((zs[2] - zs[1]) * n[0] - (xs[2] - xs[1]) * n[2])*f30/b ;
280 m0[2] = ((xs[2] - xs[1]) * n[1] - (ys[2] - ys[1]) * n[0])*f30/b ;
281
282 m1[0] = ((ys[0] - ys[2]) * n[2] - (zs[0] - zs[2]) * n[1])*f31/c ;
283 m1[1] = ((zs[0] - zs[2]) * n[0] - (xs[0] - xs[2]) * n[2])*f31/c ;
284 m1[2] = ((xs[0] - xs[2]) * n[1] - (ys[0] - ys[2]) * n[0])*f31/c ;
285
286 m2[0] = (u[1] * n[2] - u[2]* n[1])*f32 ;
287 m2[1] = (u[2] * n[0] - u[0]* n[2])*f32 ;
288 m2[2] = (u[0] * n[1] - u[1]* n[0])*f32 ;
289
290
291 Iua = (u[0] * (m0[0] + m1[0] + m2[0]) +
292 u[1] * (m0[1] + m1[1] + m2[1]) +
293 u[2] * (m0[2] + m1[2] + m2[2]))/2 ;
294
295 Iva = (v[0] * (m0[0] + m1[0] + m2[0]) +
296 v[1] * (m0[1] + m1[1] + m2[1]) +
297 v[2] * (m0[2] + m1[2] + m2[2]))/2 ;
298
299 switch(EntityNum){
300 case 1 :
301 N10 = 1 - xl/a + (u2/a -1) * yl/v2 ;
302 Val->Val[0] = N10 * I1 - Iua/a + (u2/a-1) * Iva/v2 ;
303 break;
304 case 2 :
305 N20 = xl/a - u2/a * yl/v2 ;
306 Val->Val[0] = N20 * I1 + Iua/a - u2/a * Iva/v2 ;
307 break;
308 case 3 :
309 N30 = yl/v2 ;
310 Val->Val[0] = N30 * I1 + Iva/v2 ;
311 break;
312 }
313 break;
314 default :
315 Message::Error("Unknown Basis Function Type for 'GF_LaplacexForm'");
316 }
317
318 Val->Val[0] *= ONE_OVER_FOUR_PI ;
319 if (j == 0){ Val->Val[0] /= 2; }
320 Val->Type = SCALAR ;
321 break ;
322
323 default :
324 Message::Error("Unknown Element Type (%s) for 'GF_LaplacexForm'",
325 Get_StringForDefine(Element_Type, Element->ElementSource->Type));
326 }
327 break ;
328
329 default :
330 Message::Error("Unknown Dimension (%d) for 'GF_LaplacexForm'",
331 (int)Fct->Para[0]);
332
333 }
334 #endif
335 }
336
337 /* ------------------------------------------------------------------------ */
338 /* G F _ G r a d L a p l a c e x F o r m */
339 /* ------------------------------------------------------------------------ */
340
GF_GradLaplacexForm(GF_ARGX)341 void GF_GradLaplacexForm(GF_ARGX)
342 {
343 #if !defined(HAVE_KERNEL)
344 Message::Error("GF_GradLaplacexForm requires Kernel");
345 #else
346 double xs[MAX_NODES], ys[MAX_NODES], zs[MAX_NODES] ;
347 double xxs, yys, r2, EPS ;
348 double a, b, c, a2, I1, I2 ;
349 double f0[3], f1[3], f2[3], N10, N20, N30 ;
350 double m0[3], m1[3], m2[3], s0[3], s1[3] ;
351 double umf2i, us0, us1, us2, vmf2i, vs0, vs1, vs2 ;
352 double u[3], v[3], n[3], u2, v2, xl, yl, zl, zl_2 ;
353 double area, I[3], Iua[3], Iva[3] ;
354 double s0m, s0p, s1m, s1p, s2m, s2p, t00, t10, t20, t0m_2, t0p_2, t1p_2;
355 double r00_2, r10_2, r20_2, r00, r10, r20, r0p, r0m, r1p, f20, f21, f22, B0, B1, B2, B ;
356 int Type_Int=0 ;
357
358 Val->Val[MAX_DIM] = Val->Val[MAX_DIM + 1] = Val->Val[MAX_DIM + 2] = 0. ;
359
360 switch ((int)Fct->Para[0]) {
361
362 case DIM_2D :
363
364 switch (Element->ElementSource->Type) {
365
366 case POINT_ELEMENT :
367 Val->Type = VECTOR ;
368
369 if (Element->Num == Element->ElementSource->Num) {
370 Val->Val[0] = Val->Val[1] = Val->Val[2] = 0. ;
371 return ;
372 }
373
374 xxs = x - Element->ElementSource->x[0] ;
375 yys = y - Element->ElementSource->y[0] ;
376 r2 = SQU(xxs)+SQU(yys) ;
377
378 if (r2 > EPSILON2) {
379 Val->Val[0] = - ONE_OVER_TWO_PI * xxs / r2 ;
380 Val->Val[1] = - ONE_OVER_TWO_PI * yys / r2 ;
381 Val->Val[2] = 0. ;
382 }
383 else {
384 Val->Val[0] = Val->Val[1] = Val->Val[2] = 0. ;
385 }
386 break ;
387
388 default :
389 Message::Error("Unknown Element Type (%s) for 'GF_GradLaplacexForm'",
390 Get_StringForDefine(Element_Type, Element->ElementSource->Type));
391 }
392 break ;
393
394
395 case DIM_3D :
396
397 switch (Element->ElementSource->Type) {
398
399 case LINE :
400 Val->Type = VECTOR ;
401
402 xs[0] = Element->ElementSource->x[0] ; ys[0] = Element->ElementSource->y[0] ;
403 zs[0] = Element->ElementSource->z[0] ;
404 xs[1] = Element->ElementSource->x[1] ; ys[1] = Element->ElementSource->y[1] ;
405 zs[1] = Element->ElementSource->z[1] ;
406
407 a = SQU(xs[0]-xs[1]) + SQU(ys[0]-ys[1]) + SQU(zs[0]-zs[1]) ;
408 b = 2. * ((x-xs[0])*(xs[0]-xs[1]) +
409 (y-ys[0])*(ys[0]-ys[1]) +
410 (z-zs[0])*(zs[0]-zs[1])) ;
411 c = SQU(x-xs[0]) + SQU(y-ys[0]) + SQU(z-zs[0]) + SQU(RADIUS) ;
412
413 I1 = 2./(4.*a*c-b*b) * ( (2.*a+b)/sqrt(a+b+c) - b/sqrt(c) ) ;
414 I2 = 2./(-4.*a*c+b*b) * ( (2.*c+b)/sqrt(a+b+c) - 2.*sqrt(c) ) ;
415 a2 = sqrt(a) ;
416
417 Val->Val[0] = ONE_OVER_FOUR_PI * ( (xs[0]-x) * I1 + (xs[1]-xs[0]) * I2 ) * a2 ;
418 Val->Val[1] = ONE_OVER_FOUR_PI * ( (ys[0]-y) * I1 + (ys[1]-ys[0]) * I2 ) * a2 ;
419 Val->Val[2] = ONE_OVER_FOUR_PI * ( (zs[0]-z) * I1 + (zs[1]-zs[0]) * I2 ) * a2 ;
420 break ;
421
422 case TRIANGLE :
423 Val->Type = VECTOR ;
424
425 xs[0] = Element->ElementSource->x[0] ; ys[0] = Element->ElementSource->y[0] ;
426 zs[0] = Element->ElementSource->z[0] ;
427 xs[1] = Element->ElementSource->x[1] ; ys[1] = Element->ElementSource->y[1] ;
428 zs[1] = Element->ElementSource->z[1] ;
429 xs[2] = Element->ElementSource->x[2] ; ys[2] = Element->ElementSource->y[2] ;
430 zs[2] = Element->ElementSource->z[2] ;
431
432
433 if(xFunctionBF == (void(*)())BF_Volume) Type_Int = 1 ;
434 if(xFunctionBF == (void(*)())BF_Node) Type_Int = 2 ;
435
436 /* triangle side lengths */
437 a = sqrt(SQU(xs[1]-xs[0]) + SQU(ys[1]-ys[0]) + SQU(zs[1]-zs[0]));
438 b = sqrt(SQU(xs[2]-xs[1]) + SQU(ys[2]-ys[1]) + SQU(zs[2]-zs[1]));
439 c = sqrt(SQU(xs[2]-xs[0]) + SQU(ys[2]-ys[0]) + SQU(zs[2]-zs[0]));
440
441 /* local system (u,v,w) centered at (xs[0],ys[0],zs[0]) */
442 u[0] = (xs[1]-xs[0])/a;
443 u[1] = (ys[1]-ys[0])/a;
444 u[2] = (zs[1]-zs[0])/a;
445
446 /* triangle normal */
447 Geo_CreateNormal(Element->ElementSource->Type,xs,ys,zs,n);
448
449 v[0] = n[1]*u[2]-n[2]*u[1];
450 v[1] = n[2]*u[0]-n[0]*u[2];
451 v[2] = n[0]*u[1]-n[1]*u[0];
452
453 u2 = (xs[2]-xs[0])*u[0] + (ys[2]-ys[0])*u[1] + (zs[2]-zs[0])*u[2]; /* u2 coordinate */
454 v2 = (xs[2]-xs[0])*v[0] + (ys[2]-ys[0])*v[1] + (zs[2]-zs[0])*v[2]; /* triangle height, v2 coordinate*/
455
456 /* local coordinates of the observation point (xl, yl, zl)*/
457 xl = u[0] * (x-xs[0]) + u[1] * (y-ys[0]) + u[2] * (z-zs[0]);
458 yl = v[0] * (x-xs[0]) + v[1] * (y-ys[0]) + v[2] * (z-zs[0]);
459 zl = n[0] * (x-xs[0]) + n[1] * (y-ys[0]) + n[2] * (z-zs[0]);
460
461 area = a * v2/2 ;/* Triangle area */
462 if (!zl) zl = sqrt(area) * 1e-15 ;
463
464 s0m = -( (a-xl) * (a-u2) + yl*v2 ) / b;
465 s0p = s0m + b;
466 s1p = ( xl * u2 + yl * v2 ) / c;
467 s1m = s1p - c;
468 s2m = - xl;
469 s2p = a - xl;
470
471 /* distance observation point projection on triangle plane to triangle local vertices*/
472 t00 = (yl * (u2-a) + v2 * (a-xl)) / b;
473 t10 = (xl * v2 - yl * u2) / c;
474 t20 = yl;
475
476 t0m_2 = ((a-xl)*(a-xl) + yl*yl);
477 t0p_2 = ((u2-xl)*(u2-xl) + (v2-yl)*(v2-yl));
478 t1p_2 = (xl*xl + yl*yl);
479
480 /* minimum distances^2 from the observation point to each triangle side*/
481 zl_2 = SQU(zl) ;
482 r00_2 = SQU(t00) + zl_2 ;
483 r10_2 = SQU(t10) + zl_2 ;
484 r20_2 = SQU(t20) + zl_2 ;
485
486 r00 = sqrt(r00_2);
487 r10 = sqrt(r10_2);
488 r20 = sqrt(r20_2);
489
490 /* distances from observation point to the vertices*/
491 r0p = sqrt(t0p_2 + zl_2);
492 r0m = sqrt(t0m_2 + zl_2);
493 r1p = sqrt(t1p_2 + zl_2);
494
495 EPS = EPSILON*(fabs(s0m)+fabs(s0p));
496
497 B0 = (r00 <= EPS) ? 0. : atan(t00*s0p/(r00_2+fabs(zl)*r0p))-atan(t00*s0m/(r00_2+fabs(zl)*r0m));
498 f20 = ((r0m + s0m) <= EPS) ? log(s0m/s0p) : log((r0p + s0p) / (r0m + s0m)) ;
499
500 EPS = EPSILON*(fabs(s1m)+fabs(s1p)) ;
501 B1 = (r10 <=EPS) ? 0. : atan(t10*s1p/(r10_2+fabs(zl)*r1p))-atan(t10*s1m/(r10_2+fabs(zl)*r0p));
502 f21 = ((r0p + s1m)<=EPS) ? log(s1m/s1p) : log((r1p + s1p) / (r0p + s1m));
503
504 EPS = EPSILON*(fabs(s2m)+fabs(s2p)) ;
505 B2 = (r20 <= EPS) ? 0. : atan(t20*s2p/(r20_2+fabs(zl)*r0m))-atan(t20*s2m/(r20_2+fabs(zl)*r1p));
506 f22 = ((r1p + s2m)< EPS) ? log(s2m/s2p): log((r0m + s2p) / (r1p + s2m));
507
508 B = B0 + B1 + B2 ;
509
510 s0[0] = (xs[2] - xs[1])/b ;
511 s0[1] = (ys[2] - ys[1])/b ;
512 s0[2] = (zs[2] - zs[1])/b ;
513
514 s1[0] = (xs[0] - xs[2])/c ;
515 s1[1] = (ys[0] - ys[2])/c ;
516 s1[2] = (zs[0] - zs[2])/c ;
517
518 m0[0] = s0[1] * n[2] - s0[2]* n[1] ;
519 m0[1] = s0[2] * n[0] - s0[0]* n[2] ;
520 m0[2] = s0[0] * n[1] - s0[1]* n[0] ;
521
522 m1[0] = s1[1]* n[2] - s1[2] * n[1] ;
523 m1[1] = s1[2]* n[0] - s1[0] * n[2] ;
524 m1[2] = s1[0]* n[1] - s1[1] * n[0] ;
525
526 m2[0] = u[1] * n[2] - u[2]* n[1] ;
527 m2[1] = u[2] * n[0] - u[0]* n[2] ;
528 m2[2] = u[0] * n[1] - u[1]* n[0] ;
529
530 /* Grad(1/r) integral solution*/
531 I[0] = -n[0] * THESIGN(zl) * B - (m0[0]*f20 + m1[0]*f21 + m2[0]*f22) ;
532 I[1] = -n[1] * THESIGN(zl) * B - (m0[1]*f20 + m1[1]*f21 + m2[1]*f22) ;
533 I[2] = -n[2] * THESIGN(zl) * B - (m0[2]*f20 + m1[2]*f21 + m2[2]*f22) ;
534
535 switch ( Type_Int ){
536 case 1 : /* BF_Volume */
537 Val->Val[0] = I[0]/area ;
538 Val->Val[1] = I[1]/area ;
539 Val->Val[2] = I[2]/area ;
540 break;
541
542 case 2 : /* BF_Node */
543 if (!v2 )
544 Message::Error("1/0 in GF_LaplacexForm (case _3D TRIANGLE) v2 %e", v2);
545
546 f0[0] = s0[0] * t00 * f20 - m0[0]*(r0p-r0m) ; /* fi */
547 f0[1] = s0[1] * t00 * f20 - m0[1]*(r0p-r0m) ;
548 f0[2] = s0[2] * t00 * f20 - m0[2]*(r0p-r0m) ;
549
550 f1[0] = s1[0] * t10 * f21 - m1[0]*(r1p-r0p) ;
551 f1[1] = s1[1] * t10 * f21 - m1[1]*(r1p-r0p) ;
552 f1[2] = s1[2] * t10 * f21 - m1[2]*(r1p-r0p) ;
553
554 f2[0] = u[0] * t20 * f22 - m2[0]*(r0m-r1p) ;
555 f2[1] = u[1] * t20 * f22 - m2[1]*(r0m-r1p) ;
556 f2[2] = u[2] * t20 * f22 - m2[2]*(r0m-r1p) ;
557
558 umf2i = u[0]*(m0[0]*f20 + m1[0]*f21 + m2[0]*f22) +
559 u[1]*(m0[1]*f20 + m1[1]*f21 + m2[1]*f22) +
560 u[2]*(m0[2]*f20 + m1[2]*f21 + m2[2]*f22) ;
561
562 us0 = u[0] * s0[0] + u[1] * s0[1] + u[2] * s0[2] ;
563 us1 = u[0] * s1[0] + u[1] * s1[1] + u[2] * s1[2] ;
564 us2 = u[0] * u[0] + u[1] * u[1] + u[2] * u[2] ;
565
566 vmf2i = v[0]*(m0[0]*f20 + m1[0]*f21 + m2[0]*f22) +
567 v[1]*(m0[1]*f20 + m1[1]*f21 + m2[1]*f22) +
568 v[2]*(m0[2]*f20 + m1[2]*f21 + m2[2]*f22) ;
569
570 vs0 = v[0] * s0[0] + v[1] * s0[1] + v[2] * s0[2] ;
571 vs1 = v[0] * s1[0] + v[1] * s1[1] + v[2] * s1[2] ;
572 vs2 = v[0] * u[0] + v[1] * u[1] + v[2] * u[2] ;
573
574 B *= fabs(zl);
575 umf2i *= zl ;
576 vmf2i *= zl ;
577
578 Iua[0] = n[0] * umf2i - B * u[0] + f0[0] * us0 + f1[0] * us1 + f2[0] * us2 ;
579 Iua[1] = n[1] * umf2i - B * u[1] + f0[1] * us0 + f1[1] * us1 + f2[1] * us2 ;
580 Iua[2] = n[2] * umf2i - B * u[2] + f0[2] * us0 + f1[2] * us1 + f2[2] * us2 ;
581
582 Iva[0] = n[0] * vmf2i - B * v[0] + f0[0] * vs0 + f1[0] * vs1 + f2[0] * vs2 ;
583 Iva[1] = n[1] * vmf2i - B * v[1] + f0[1] * vs0 + f1[1] * vs1 + f2[1] * vs2 ;
584 Iva[2] = n[2] * vmf2i - B * v[2] + f0[2] * vs0 + f1[2] * vs1 + f2[2] * vs2 ;
585
586 switch(EntityNum){
587 case 1 :
588 N10 = 1 - xl/a + (u2/a -1) * yl/v2 ;
589 Val->Val[0] = N10 * I[0] - Iua[0]/a + (u2/a-1) * Iva[0]/v2 ;
590 Val->Val[1] = N10 * I[1] - Iua[1]/a + (u2/a-1) * Iva[1]/v2 ;
591 Val->Val[2] = N10 * I[2] - Iua[2]/a + (u2/a-1) * Iva[2]/v2 ;
592 break;
593 case 2 :
594 N20 = xl/a - u2/a * yl/v2 ;
595 Val->Val[0] = N20 * I[0] + Iua[0]/a - u2/a * Iva[0]/v2 ;
596 Val->Val[1] = N20 * I[1] + Iua[1]/a - u2/a * Iva[1]/v2 ;
597 Val->Val[2] = N20 * I[2] + Iua[2]/a - u2/a * Iva[2]/v2 ;
598 break;
599 case 3 :
600 N30 = yl/v2 ;
601 Val->Val[0] = N30 * I[0] + Iva[0]/v2 ;
602 Val->Val[1] = N30 * I[1] + Iva[1]/v2 ;
603 Val->Val[2] = N30 * I[2] + Iva[2]/v2 ;
604 break;
605 }
606 break;
607 }
608
609 Val->Val[0] *= ONE_OVER_FOUR_PI ;
610 Val->Val[1] *= ONE_OVER_FOUR_PI ;
611 Val->Val[2] *= ONE_OVER_FOUR_PI ;
612 break ;
613
614 default :
615 Message::Error("Unknown Element Type (%s) for 'GF_GradLaplacexForm'",
616 Get_StringForDefine(Element_Type, Element->ElementSource->Type));
617 }
618 break ;
619
620 default :
621 Message::Error("Unknown Dimension (%d) for 'GF_GradLaplacexForm'",
622 (int)Fct->Para[0]);
623
624 }
625 #endif
626 }
627
628 /* ------------------------------------------------------------------------ */
629 /* G F _ N P x G r a d L a p l a c e x F o r m */
630 /* ------------------------------------------------------------------------ */
631
GF_NPxGradLaplacexForm(GF_ARGX)632 void GF_NPxGradLaplacexForm(GF_ARGX)
633 {
634 #if !defined(HAVE_KERNEL)
635 Message::Error("GF_NPGradLaplacexForm requires Kernel");
636 #else
637 double xs[MAX_NODES], ys[MAX_NODES] ;
638 double xp[MAX_NODES], yp[MAX_NODES], N[3] ;
639 int Type_Int;
640 double a, b, c, d, m, n, Jp, i1, Is, I1=0 ;
641 struct Value ValGrad ;
642
643 Val->Type = SCALAR ;
644 Val->Val[MAX_DIM] = 0.0 ;
645
646 if (Element->Num == Element->ElementSource->Num) {
647 Val->Val[0] = 0.0 ;
648 return ;
649 }
650
651 switch ((int)Fct->Para[0]) {
652
653 case DIM_2D :
654
655 switch (Element->ElementSource->Type) {
656
657 case LINE :
658 if (Element->Type != LINE)
659 Message::Error("GF_NPxGradLaplacexForm not ready for mixed geometrical elements");
660
661 xs[0] = Element->ElementSource->x[0] ; ys[0] = Element->ElementSource->y[0] ;
662 xs[1] = Element->ElementSource->x[1] ; ys[1] = Element->ElementSource->y[1] ;
663
664 if(xFunctionBF == (void(*)())BF_Volume) {
665 if ((x == xs[0]) && (y == ys[0]))
666 Type_Int = 1 ;
667 else if ((x == xs[1]) && (y == ys[1]))
668 Type_Int = 2 ;
669 else
670 Type_Int = 3 ;
671
672 xp[0] = Element->x[0] ; yp[0] = Element->y[0] ;
673 xp[1] = Element->x[1] ; yp[1] = Element->y[1] ;
674
675 a = SQU(xs[0]-xs[1]) + SQU(ys[0]-ys[1]) ;
676 b = 2. * ((x-xs[0]) * (xs[0]-xs[1]) + (y-ys[0]) * (ys[0]-ys[1])) ;
677 c = SQU(x-xs[0]) + SQU(y-ys[0]) ;
678 d = 4.*a*c - b*b ;
679
680 switch (Type_Int) {
681 case 1 :
682 case 2 :
683 Message::Error("Degenerate case not done in 'GF_NPxGradLaplacexForm'");
684 break ;
685 case 3 :
686 if (fabs(d) < EPSILON2) {
687 I1 = 0.0 ;
688 }
689 else {
690 if(d<0) Message::Error("Unexpected value in 'GF_NPxGradLaplacexForm'");
691 i1 = sqrt(d) ;
692 Is = 2. / i1 * (atan((2.*a+b)/i1) - atan(b/i1)) ;
693 Jp = sqrt(SQU(xp[0]-xp[1])+SQU(yp[0]-yp[1])) ;
694 m = ((ys[0]-ys[1]) * (xp[0]-xp[1]) + (xs[0]-xs[1]) * (yp[1]-yp[0])) / Jp ;
695 n = ((yp[1]-yp[0]) * (x-xs[0]) + (xp[0]-xp[1]) * (y-ys[0])) / Jp ;
696 I1 = m /(2.*a) * log((a+b)/c+1.) + (n - m*b/(2.*a)) * Is ;
697 }
698 break ;
699 }
700 Val->Val[0] = - ONE_OVER_TWO_PI * I1 ;
701 }
702 else {
703 Message::Error("Unknown Basis Function Type for 'GF_NPxGradLaplacexForm'");
704 }
705 break ;
706
707 default :
708 Message::Error("Unknown Element Type (%s) for 'GF_NPxGradLaplacexForm'",
709 Get_StringForDefine(Element_Type, Element->ElementSource->Type));
710 }
711 break ;
712
713 case DIM_3D:
714 switch (Element->ElementSource->Type) {
715 case TRIANGLE :
716 Geo_CreateNormal(Element->Type,
717 Element->x,Element->y,Element->z, N);
718
719 GF_GradLaplacexForm(Element, Fct, xFunctionBF, EntityNum, x, y, z, &ValGrad) ;
720
721 Val->Val[0] = N[0]*ValGrad.Val[0] + N[1]*ValGrad.Val[1] + N[2]*ValGrad.Val[2] ;
722 break ;
723 default :
724 Message::Error("Unknown Element Type (%s) for 'GF_NPxGradLaplacexForm'",
725 Get_StringForDefine(Element_Type, Element->ElementSource->Type));
726 }
727 break ;
728
729 default :
730 Message::Error("Unknown Dimension (%d) for 'GF_NPxGradLaplacexForm'",
731 (int)Fct->Para[0]);
732 }
733 #endif
734 }
735
736
737
738 /* ------------------------------------------------------------------------ */
739 /* G F _ N S x G r a d L a p l a c e x F o r m */
740 /* ------------------------------------------------------------------------ */
741
GF_NSxGradLaplacexForm(GF_ARGX)742 void GF_NSxGradLaplacexForm(GF_ARGX)
743 {
744 Message::Error("Not done: 'GF_NSxGradLaplacexForm'");
745 }
746
747 /* ------------------------------------------------------------------------ */
748 /* G F _ A p p r o x i m a t e L a p l a c e x F o r m */
749 /* ------------------------------------------------------------------------ */
750
GF_ApproximateLaplacexForm(GF_ARGX)751 void GF_ApproximateLaplacexForm(GF_ARGX)
752 {
753 switch ((int)Fct->Para[1]) {
754
755 case 0 :
756 GF_LaplacexForm(Element, Fct, (void(*)())BF_Volume, 1, x, y, z, Val);
757 break ;
758
759 default :
760 Message::Error("Bad Parameter Value in 'GF_ApproximateLaplacexForm'");
761 break;
762
763 }
764 }
765