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 #include <math.h>
7 #include "ProData.h"
8 #include "GF.h"
9 #include "Message.h"
10
11 #define SQU(a) ((a)*(a))
12 #define CUB(a) ((a)*(a)*(a))
13 #define ONE_OVER_TWO_PI 1.5915494309189534E-01
14 #define ONE_OVER_FOUR_PI 7.9577471545947668E-02
15
16 extern struct CurrentData Current ;
17
18 /* ------------------------------------------------------------------------ */
19 /* G F _ L a p l a c e */
20 /* ------------------------------------------------------------------------ */
21
GF_Laplace(GF_ARG)22 void GF_Laplace(GF_ARG)
23 {
24 double d;
25
26 switch((int)Fct->Para[0]){
27
28 case DIM_1D : /* r/2 */
29 V->Val[0] = 0.5 * sqrt(SQU(Current.x - Current.xs)) ;
30 break;
31
32 case DIM_2D : /* 1/(2*Pi) * ln(1/r) = -1/(4*Pi)*ln(r^2) */
33 d = SQU(Current.x - Current.xs) + SQU(Current.y - Current.ys) ;
34 if(!d) Message::Error("Log(0) in 'GF_Laplace'") ;
35 V->Val[0] = - ONE_OVER_FOUR_PI * log(d) ;
36 V->Val[MAX_DIM] = 0. ;
37 break;
38
39 case DIM_3D : /* 1/(4*Pi*r) */
40 d = SQU(Current.x - Current.xs) + SQU(Current.y - Current.ys) +
41 SQU(Current.z - Current.zs) ;
42 if(!d) Message::Error("1/0 in 'GF_Laplace'") ;
43 V->Val[0] = ONE_OVER_FOUR_PI / sqrt(d) ;
44 V->Val[MAX_DIM] = 0. ;
45 break;
46
47 default :
48 Message::Error("Bad Parameter for 'GF_Laplace' (%d)", (int)Fct->Para[0]);
49 break;
50 }
51
52 V->Type = SCALAR ;
53 V->Val[MAX_DIM] = 0. ;
54 }
55
56
57 /* ------------------------------------------------------------------------ */
58 /* G F _ G r a d L a p l a c e */
59 /* ------------------------------------------------------------------------ */
60
61 /* the gradient is taken relative to the destination point (x,y,z) */
62
GF_GradLaplace(GF_ARG)63 void GF_GradLaplace(GF_ARG)
64 {
65 double xxs, yys, zzs, r;
66
67 V->Type = VECTOR ;
68
69 switch((int)Fct->Para[0]){
70 case DIM_2D :
71 xxs = Current.x-Current.xs ;
72 yys = Current.y-Current.ys ;
73 r = SQU(xxs) + SQU(yys) ;
74 if(!r) Message::Error("1/0 in 'GF_GradLaplace'") ;
75 V->Val[0] = - ONE_OVER_TWO_PI * xxs / r ;
76 V->Val[1] = - ONE_OVER_TWO_PI * yys / r ;
77 V->Val[2] = 0.0 ;
78 V->Val[MAX_DIM ] = V->Val[MAX_DIM + 1] = V->Val[MAX_DIM + 2] = 0. ;
79 break ;
80
81 case DIM_3D :
82 xxs = Current.x-Current.xs ;
83 yys = Current.y-Current.ys ;
84 zzs = Current.z-Current.zs ;
85 r = CUB(sqrt(SQU(xxs) + SQU(yys) + SQU(zzs))) ;
86 if(!r) Message::Error("1/0 in 'GF_GradLaplace'") ;
87 V->Val[0] = - ONE_OVER_FOUR_PI * xxs / r ;
88 V->Val[1] = - ONE_OVER_FOUR_PI * yys / r ;
89 V->Val[2] = - ONE_OVER_FOUR_PI * zzs / r ;
90 V->Val[MAX_DIM ] = V->Val[MAX_DIM + 1] = V->Val[MAX_DIM + 2] = 0. ;
91 break ;
92
93 default :
94 Message::Error("Bad Parameter for 'GF_GradLaplace' (%d)", (int)Fct->Para[0]);
95 break;
96 }
97
98 V->Type = VECTOR ;
99
100 V->Val[MAX_DIM+0] = 0. ;
101 V->Val[MAX_DIM+1] = 0. ;
102 V->Val[MAX_DIM+2] = 0. ;
103 }
104
105 /* ------------------------------------------------------------------------ */
106 /* G F _ N P x G r a d L a p l a c e */
107 /* ------------------------------------------------------------------------ */
108
GF_NPxGradLaplace(GF_ARG)109 void GF_NPxGradLaplace(GF_ARG)
110 {
111 double x1x0, x2x0, y1y0, y2y0, z1z0, z2z0, xxs, yys, zzs, a, b, c, r ;
112
113 V->Type = SCALAR ;
114 V->Val[MAX_DIM] = 0. ;
115
116 if (Current.Element->Num == Current.ElementSource->Num) {
117 V->Val[0 ] = 0. ;
118 V->Val[MAX_DIM] = 0. ;
119 return ;
120 }
121
122 switch((int)Fct->Para[0]){
123 case DIM_2D :
124 x1x0 = Current.Element->x[1] - Current.Element->x[0] ;
125 y1y0 = Current.Element->y[1] - Current.Element->y[0] ;
126 xxs = Current.x-Current.xs ;
127 yys = Current.y-Current.ys ;
128 r = SQU(xxs)+SQU(yys) ;
129 if(!r) V->Val[0] = 0 ;
130 else
131 V->Val[0] = - ONE_OVER_TWO_PI * (y1y0 * xxs - x1x0 * yys)
132 / (sqrt(SQU(x1x0)+SQU(y1y0)) * r) ;
133
134 V->Val[MAX_DIM] = 0. ;
135 break;
136
137 case DIM_3D :
138 x1x0 = Current.Element->x[1] - Current.Element->x[0] ;
139 y1y0 = Current.Element->y[1] - Current.Element->y[0] ;
140 z1z0 = Current.Element->z[1] - Current.Element->z[0] ;
141 x2x0 = Current.Element->x[2] - Current.Element->x[0] ;
142 y2y0 = Current.Element->y[2] - Current.Element->y[0] ;
143 z2z0 = Current.Element->z[2] - Current.Element->z[0] ;
144 a = y1y0 * z2z0 - z1z0 * y2y0 ;
145 b = z1z0 * x2x0 - x1x0 * z2z0 ;
146 c = x1x0 * y2y0 - y1y0 * x2x0 ;
147 xxs = Current.x-Current.xs ;
148 yys = Current.y-Current.ys ;
149 zzs = Current.z-Current.zs ;
150 V->Val[0] = - ONE_OVER_FOUR_PI * (a * xxs + b * yys + c * zzs)
151 / ( (sqrt(SQU(a) + SQU(b) + SQU(c)) *
152 CUB(sqrt(SQU(xxs) + SQU(yys) + SQU(zzs)))) ) ;
153 V->Val[MAX_DIM] = 0. ;
154 break ;
155
156 default :
157 Message::Error("Bad Parameter for 'GF_NPxGradLaplace' (%d)", (int)Fct->Para[0]);
158 break;
159 }
160 }
161
162 /* ------------------------------------------------------------------------ */
163 /* G F _ N S x G r a d L a p l a c e */
164 /* ------------------------------------------------------------------------ */
165
GF_NSxGradLaplace(GF_ARG)166 void GF_NSxGradLaplace(GF_ARG)
167 {
168 double x1x0, x2x0, y1y0, y2y0, z1z0, z2z0, xxs, yys, zzs, a, b, c;
169
170 V->Type = SCALAR ;
171 V->Val[MAX_DIM] = 0. ;
172
173 if (Current.Element->Num == Current.ElementSource->Num) {
174 V->Val[0] = 0. ;
175 V->Val[MAX_DIM] = 0. ;
176 return ;
177 }
178
179 switch((int)Fct->Para[0]){
180 case DIM_2D :
181 x1x0 = Current.ElementSource->x[1] - Current.ElementSource->x[0] ;
182 y1y0 = Current.ElementSource->y[1] - Current.ElementSource->y[0] ;
183 xxs = Current.x-Current.xs ;
184 yys = Current.y-Current.ys ;
185 V->Val[0] = ONE_OVER_TWO_PI * (y1y0 * xxs - x1x0 * yys)
186 / (sqrt(SQU(x1x0)+SQU(y1y0)) * (SQU(xxs)+SQU(yys)) ) ;
187 V->Val[MAX_DIM] = 0. ;
188 break ;
189 case DIM_3D :
190 x1x0 = Current.ElementSource->x[1] - Current.ElementSource->x[0] ;
191 y1y0 = Current.ElementSource->y[1] - Current.ElementSource->y[0] ;
192 z1z0 = Current.ElementSource->z[1] - Current.ElementSource->z[0] ;
193 x2x0 = Current.ElementSource->x[2] - Current.ElementSource->x[0] ;
194 y2y0 = Current.ElementSource->y[2] - Current.ElementSource->y[0] ;
195 z2z0 = Current.ElementSource->z[2] - Current.ElementSource->z[0] ;
196 a = y1y0 * z2z0 - z1z0 * y2y0 ;
197 b = z1z0 * x2x0 - x1x0 * z2z0 ;
198 c = x1x0 * y2y0 - y1y0 * x2x0 ;
199 xxs = Current.x-Current.xs ;
200 yys = Current.y-Current.ys ;
201 zzs = Current.z-Current.zs ;
202 V->Val[0] = ONE_OVER_FOUR_PI * (a * xxs + b * yys + c * zzs)
203 / ( (sqrt(SQU(a)+SQU(b)+SQU(c)) *
204 CUB(sqrt(SQU(xxs)+SQU(yys)+SQU(zzs)))) ) ;
205 V->Val[MAX_DIM] = 0. ;
206 break ;
207 default :
208 Message::Error("Bad Parameter for 'GF_NSxGradLaplace' (%d)", (int)Fct->Para[0]);
209 break;
210 }
211 }
212
213 /* ------------------------------------------------------------------------ */
214 /* G F _ A p p r o x i m a t e L a p l a c e */
215 /* ------------------------------------------------------------------------ */
216
GF_ApproximateLaplace(GF_ARG)217 void GF_ApproximateLaplace(GF_ARG)
218 {
219 Message::Error("The Approximate Integral Kernels can only be Integrated Analytically");
220 }
221