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 <stdlib.h>
7 #include <math.h>
8 #include "ProData.h"
9 #include "GeoData.h"
10 #include "DofData.h"
11 #include "Cal_Quantity.h"
12 #include "Get_Geometry.h"
13 #include "Message.h"
14
15 #define SQU(a) ((a)*(a))
16
17 extern struct Problem Problem_S ;
18 extern struct CurrentData Current ;
19
20 /* ------------------------------------------------------------------------ */
21 /* G e t _ V a l u e F r o m F o r m */
22 /* ------------------------------------------------------------------------ */
23
Get_ValueFromForm(int Form)24 int Get_ValueFromForm(int Form)
25 {
26 switch (Form) {
27 case FORM0 :
28 case FORM3 : case FORM3P :
29 case SCALAR :
30 return(SCALAR) ;
31
32 case FORM1 : case FORM1P : case FORM1S :
33 case FORM2 : case FORM2P : case FORM2S :
34 case VECTOR : case VECTORP :
35 return(VECTOR) ;
36
37 default :
38 Message::Error("Unknown Form type in 'Get_ValueFromForm'");
39 return(-1) ;
40 }
41 }
42
43 /* ------------------------------------------------------------------------ */
44 /* G e t _ I n t e g r a t i o n C a s e */
45 /* ------------------------------------------------------------------------ */
46
47 /*
48 Il faudrait reorganiser les 'Current.XXX'
49 Ca devient un peu le bordel.
50 */
51
Get_IntegrationCase(struct Element * Element,List_T * IntegrationCase_L,int CriterionIndex)52 struct IntegrationCase * Get_IntegrationCase (struct Element * Element,
53 List_T *IntegrationCase_L,
54 int CriterionIndex)
55 {
56 struct Value Criterion ;
57
58 if (CriterionIndex >= 0){
59 Current.Element = Element ;
60 Current.ElementSource = Element->ElementSource ;
61 Get_ValueOfExpression
62 ((struct Expression *)
63 List_Pointer(Problem_S.Expression, CriterionIndex),
64 NULL, 0., 0., 0., &Criterion) ;
65 if(Criterion.Val[0] < 0 || Criterion.Val[0] >= List_Nbr(IntegrationCase_L))
66 Message::Error("Integration criterion out of range");
67 }
68 else {
69 if(List_Nbr(IntegrationCase_L) > 1)
70 Message::Error("Missing integration criterion");
71 Criterion.Val[0] = 0;
72 }
73
74 return((struct IntegrationCase*) List_Pointer(IntegrationCase_L,
75 (int)Criterion.Val[0])) ;
76 }
77
78 /* ------------------------------------------------------------------------ */
79 /* G e t _ F u n c t i o n V a l u e */
80 /* ------------------------------------------------------------------------ */
81
Get_FunctionValue(int Nbr_Function,void (* xFunctionBF[])(),int Type_Operator,struct QuantityStorage * QuantityStorage_P,int * Type_Form)82 void Get_FunctionValue(int Nbr_Function,
83 void (*xFunctionBF[])(),
84 int Type_Operator,
85 struct QuantityStorage * QuantityStorage_P,
86 int * Type_Form)
87 {
88 int i ;
89
90 switch (Type_Operator) {
91
92 case NOOP :
93 *Type_Form = QuantityStorage_P->TypeQuantity ;
94 for (i = 0 ; i < Nbr_Function ; i++)
95 xFunctionBF[i] =
96 QuantityStorage_P->BasisFunction[i].BasisFunction->Function ;
97 break ;
98
99 case EXTDER :
100 *Type_Form = QuantityStorage_P->TypeQuantity + 1 ;
101 for (i = 0 ; i < Nbr_Function ; i++)
102 xFunctionBF[i] =
103 QuantityStorage_P->BasisFunction[i].BasisFunction->dFunction ;
104 break ;
105
106 case EXTDERINV :
107 *Type_Form = QuantityStorage_P->TypeQuantity - 1 ;
108 for (i = 0 ; i < Nbr_Function ; i++)
109 xFunctionBF[i] =
110 QuantityStorage_P->BasisFunction[i].BasisFunction->dInvFunction ;
111 break ;
112
113 case GRAD :
114 if (QuantityStorage_P->TypeQuantity == FORM0) {
115 *Type_Form = QuantityStorage_P->TypeQuantity + 1 ;
116 for (i = 0 ; i < Nbr_Function ; i++)
117 xFunctionBF[i] =
118 QuantityStorage_P->BasisFunction[i].BasisFunction->dFunction ;
119 }
120 else if (QuantityStorage_P->TypeQuantity == SCALAR) {
121 *Type_Form = VECTOR ;
122 }
123 else{
124 Message::Error("Cannot apply Grad operator to quantity type %d",
125 QuantityStorage_P->TypeQuantity);
126 }
127 break ;
128
129 case CURL :
130 if ((QuantityStorage_P->TypeQuantity == FORM1) ||
131 (QuantityStorage_P->TypeQuantity == FORM1P)) {
132 *Type_Form = QuantityStorage_P->TypeQuantity + 1 ;
133 for (i = 0 ; i < Nbr_Function ; i++)
134 xFunctionBF[i] =
135 QuantityStorage_P->BasisFunction[i].BasisFunction->dFunction ;
136 }
137 else if (QuantityStorage_P->TypeQuantity == VECTOR) {
138 *Type_Form = VECTOR ;
139 for (i = 0 ; i < Nbr_Function ; i++)
140 xFunctionBF[i] =
141 QuantityStorage_P->BasisFunction[i].BasisFunction->dFunction ;
142 }
143 else{
144 Message::Error("Cannot apply Curl operator to quantity type %d",
145 QuantityStorage_P->TypeQuantity);
146 }
147 break ;
148
149 case DIV :
150 if (QuantityStorage_P->TypeQuantity == FORM2) {
151 *Type_Form = QuantityStorage_P->TypeQuantity + 1 ;
152 for (i = 0 ; i < Nbr_Function ; i++)
153 xFunctionBF[i] =
154 QuantityStorage_P->BasisFunction[i].BasisFunction->dFunction ;
155 }
156 else if (QuantityStorage_P->TypeQuantity == VECTOR) {
157 *Type_Form = SCALAR ;
158 for (i = 0 ; i < Nbr_Function ; i++)
159 xFunctionBF[i] =
160 QuantityStorage_P->BasisFunction[i].BasisFunction->dInvFunction ;
161 }
162 else{
163 Message::Error("Cannot apply Div operator to quantity type %d",
164 QuantityStorage_P->TypeQuantity);
165 }
166 break ;
167
168 case GRADINV :
169 if (QuantityStorage_P->TypeQuantity == FORM1) {
170 *Type_Form = QuantityStorage_P->TypeQuantity - 1 ;
171 for (i = 0 ; i < Nbr_Function ; i++)
172 xFunctionBF[i] =
173 QuantityStorage_P->BasisFunction[i].BasisFunction->dInvFunction ;
174 }
175 else if (QuantityStorage_P->TypeQuantity == VECTOR) {
176 *Type_Form = SCALAR ;
177 }
178 else{
179 Message::Error("Cannot apply GradInv operator to quantity type %d",
180 QuantityStorage_P->TypeQuantity);
181 }
182 break ;
183
184 case CURLINV :
185 if (QuantityStorage_P->TypeQuantity == FORM2) {
186 *Type_Form = QuantityStorage_P->TypeQuantity - 1 ;
187 for (i = 0 ; i < Nbr_Function ; i++)
188 xFunctionBF[i] =
189 QuantityStorage_P->BasisFunction[i].BasisFunction->dInvFunction ;
190 }
191 else if (QuantityStorage_P->TypeQuantity == VECTOR) {
192 *Type_Form = VECTOR ;
193 }
194 else{
195 Message::Error("Cannot apply CurlInv operator to quantity type %d",
196 QuantityStorage_P->TypeQuantity);
197 }
198 break ;
199
200 case DIVINV :
201 if ((QuantityStorage_P->TypeQuantity == FORM3) ||
202 (QuantityStorage_P->TypeQuantity == FORM3P)) {
203 *Type_Form = QuantityStorage_P->TypeQuantity - 1 ;
204 for (i = 0 ; i < Nbr_Function ; i++)
205 xFunctionBF[i] =
206 QuantityStorage_P->BasisFunction[i].BasisFunction->dInvFunction ;
207 }
208 else if (QuantityStorage_P->TypeQuantity == SCALAR) {
209 *Type_Form = VECTOR ;
210 }
211 else{
212 Message::Error("Cannot apply DivInv operator to quantity type %d",
213 QuantityStorage_P->TypeQuantity);
214 }
215 break ;
216
217 case OP_D1 :
218 if (QuantityStorage_P->TypeQuantity == VECTOR) {
219 *Type_Form = VECTOR ;
220 for (i = 0 ; i < Nbr_Function ; i++)
221 xFunctionBF[i] =
222 QuantityStorage_P->BasisFunction[i].BasisFunction->dFunction ;
223 }
224 else{
225 Message::Error("Cannot apply D1 operator to quantity type %d",
226 QuantityStorage_P->TypeQuantity);
227 }
228 break ;
229
230 case OP_D2 :
231 if (QuantityStorage_P->TypeQuantity == VECTOR) {
232 *Type_Form = VECTOR ;
233 for (i = 0 ; i < Nbr_Function ; i++)
234 xFunctionBF[i] =
235 QuantityStorage_P->BasisFunction[i].BasisFunction->dInvFunction ;
236 }
237 else{
238 Message::Error("Cannot apply D2 operator to quantity type %d",
239 QuantityStorage_P->TypeQuantity);
240 }
241 break ;
242
243 case OP_D3 :
244 if (QuantityStorage_P->TypeQuantity == VECTOR) {
245 *Type_Form = VECTOR ;
246 for (i = 0 ; i < Nbr_Function ; i++)
247 xFunctionBF[i] =
248 QuantityStorage_P->BasisFunction[i].BasisFunction->dPlusFunction ;
249 }
250 else{
251 Message::Error("Cannot apply D3 operator to quantity type %d",
252 QuantityStorage_P->TypeQuantity);
253 }
254 break ;
255
256 default :
257 Message::Error("Unknown operator in 'Get_FunctionValue'");
258 break;
259 }
260 }
261
262
263 /* ------------------------------------------------------------------------ */
264 /* G e t _ I n i t F u n c t i o n V a l u e */
265 /* ------------------------------------------------------------------------ */
266
Get_InitFunctionValue(int Type_Operator,struct QuantityStorage * QuantityStorage_P,int * Type_Form)267 void Get_InitFunctionValue(int Type_Operator,
268 struct QuantityStorage * QuantityStorage_P,
269 int * Type_Form)
270 {
271 switch (Type_Operator) {
272
273 case NOOP :
274 *Type_Form = QuantityStorage_P->TypeQuantity ; break ;
275
276 case EXTDER :
277 *Type_Form = QuantityStorage_P->TypeQuantity + 1 ; break ;
278
279 case EXTDERINV :
280 *Type_Form = QuantityStorage_P->TypeQuantity - 1 ; break ;
281
282
283 case GRAD :
284 if (QuantityStorage_P->TypeQuantity == FORM0)
285 *Type_Form = QuantityStorage_P->TypeQuantity + 1 ;
286 else if (QuantityStorage_P->TypeQuantity == SCALAR) *Type_Form = VECTOR ;
287 break ;
288
289 case CURL :
290 if ((QuantityStorage_P->TypeQuantity == FORM1) ||
291 (QuantityStorage_P->TypeQuantity == FORM1P))
292 *Type_Form = QuantityStorage_P->TypeQuantity + 1 ;
293 else if (QuantityStorage_P->TypeQuantity == VECTOR) *Type_Form = VECTOR ;
294 break ;
295
296 case DIV :
297 if (QuantityStorage_P->TypeQuantity == FORM2)
298 *Type_Form = QuantityStorage_P->TypeQuantity + 1 ;
299 else if (QuantityStorage_P->TypeQuantity == VECTOR) *Type_Form = SCALAR ;
300 break ;
301
302 case GRADINV :
303 if (QuantityStorage_P->TypeQuantity == FORM1)
304 *Type_Form = QuantityStorage_P->TypeQuantity - 1 ;
305 else if (QuantityStorage_P->TypeQuantity == VECTOR) *Type_Form = SCALAR ;
306 break ;
307
308 case CURLINV :
309 if (QuantityStorage_P->TypeQuantity == FORM2)
310 *Type_Form = QuantityStorage_P->TypeQuantity - 1 ;
311 else if (QuantityStorage_P->TypeQuantity == VECTOR) *Type_Form = VECTOR ;
312 break ;
313
314 case DIVINV :
315 if ((QuantityStorage_P->TypeQuantity == FORM3) ||
316 (QuantityStorage_P->TypeQuantity == FORM3P))
317 *Type_Form = QuantityStorage_P->TypeQuantity - 1 ;
318 else if (QuantityStorage_P->TypeQuantity == SCALAR) *Type_Form = VECTOR ;
319 break ;
320
321 case OP_D1 :
322 case OP_D2 :
323 case OP_D3 :
324 if (QuantityStorage_P->TypeQuantity == VECTOR)
325 *Type_Form = VECTOR ;
326 else
327 *Type_Form = VECTOR ;
328 break ;
329
330 default :
331 Message::Error("Unknown operator in 'Get_InitFunctionValue'");
332 break;
333 }
334 }
335
336 /* ------------------------------------------------------------------------ */
337 /* C a l _ I n t e r p o l a t i o n O r d e r */
338 /* ------------------------------------------------------------------------ */
339
Cal_InterpolationOrder(struct Element * Element,struct QuantityStorage * QuantityStorage)340 double Cal_InterpolationOrder(struct Element * Element,
341 struct QuantityStorage * QuantityStorage)
342 {
343 int i ;
344 double Order = 0.0 ;
345
346 for(i = 0 ; i < QuantityStorage->NbrElementaryBasisFunction ; i++)
347 if(QuantityStorage->BasisFunction[i].Dof->Type == DOF_UNKNOWN)
348 Order = std::max(QuantityStorage->BasisFunction[i].BasisFunction->Order, Order) ;
349
350 return(Order) ;
351 }
352
353 /* ------------------------------------------------------------------------ */
354 /* C a l _ M a x E d g e L e n g t h */
355 /* ------------------------------------------------------------------------ */
356
Cal_MaxEdgeLength(struct Element * Element)357 double Cal_MaxEdgeLength(struct Element * Element)
358 {
359 int i, *IM, *N, NbrEdges ;
360 double l, lmax = 0.0 ;
361
362 IM = Geo_GetIM_Den(Element->Type, &NbrEdges) ;
363 for(i = 0 ; i < NbrEdges ; i++){
364 N = IM + i * NBR_MAX_SUBENTITIES_IN_ELEMENT ;
365 l = sqrt(SQU(Element->x[abs(N[1])-1]-Element->x[abs(N[0])-1]) +
366 SQU(Element->y[abs(N[1])-1]-Element->y[abs(N[0])-1]) +
367 SQU(Element->z[abs(N[1])-1]-Element->z[abs(N[0])-1])) ;
368 lmax = std::max(lmax, l) ;
369 }
370
371 return(lmax) ;
372 }
373