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 "ProDefine.h"
9 #include "GeoData.h"
10 #include "Get_DofOfElement.h"
11 #include "Cal_Quantity.h"
12 #include "Cal_Value.h"
13 #include "Get_Geometry.h"
14 #include "Get_FunctionValue.h"
15 #include "ExtendedGroup.h"
16 #include "Pos_Formulation.h"
17 #include "MallocUtils.h"
18 #include "Message.h"
19
20 extern struct Problem Problem_S ;
21 extern struct CurrentData Current ;
22
23 /* ------------------------------------------------------------------------ */
24 /* P o s _ L o c a l O r I n t e g r a l Q u a n t i t y */
25 /* ------------------------------------------------------------------------ */
26
27 static int Warning_NoJacobian = 0 ;
28
Pos_LocalOrIntegralQuantity(struct PostQuantity * PostQuantity_P,struct DefineQuantity * DefineQuantity_P0,struct QuantityStorage * QuantityStorage_P0,struct PostQuantityTerm * PostQuantityTerm_P,struct Element * Element,int Type_Quantity,double u,double v,double w,struct Value * Value)29 void Pos_LocalOrIntegralQuantity(struct PostQuantity *PostQuantity_P,
30 struct DefineQuantity *DefineQuantity_P0,
31 struct QuantityStorage *QuantityStorage_P0,
32 struct PostQuantityTerm *PostQuantityTerm_P,
33 struct Element *Element,
34 int Type_Quantity,
35 double u, double v, double w,
36 struct Value *Value)
37 {
38 struct FunctionSpace *FunctionSpace_P ;
39 struct DefineQuantity *DefineQuantity_P ;
40 struct QuantityStorage *QuantityStorage_P ;
41 struct Value TermValue, tmpValue;
42 struct JacobianCase *JacobianCase_P0 ;
43 struct IntegrationCase *IntegrationCase_P ;
44 struct Quadrature *Quadrature_P ;
45
46 List_T *IntegrationCase_L, *JacobianCase_L ;
47 double ui, vi, wi, weight, Factor ;
48 int Index_DefineQuantity ;
49 int i, j, Type_Dimension ;
50 int CriterionIndex, Nbr_IntPoints, i_IntPoint ;
51
52 double (*Get_Jacobian) (struct Element * Element, MATRIX3x3 * Jac) = 0;
53 void (*Get_IntPoint) (int Nbr_Points, int Num,
54 double * u, double * v, double * w, double * wght) ;
55
56 /* Get the functionspaces
57 Get the DoF for local quantities */
58
59 for (i = 0 ; i < PostQuantityTerm_P->NbrQuantityIndex ; i++) {
60 Index_DefineQuantity = PostQuantityTerm_P->QuantityIndexTable[i] ;
61 DefineQuantity_P = DefineQuantity_P0 + Index_DefineQuantity ;
62 QuantityStorage_P = QuantityStorage_P0 + Index_DefineQuantity ;
63
64 if (QuantityStorage_P->NumLastElementForFunctionSpace != Element->Num) {
65
66 QuantityStorage_P->NumLastElementForFunctionSpace = Element->Num ;
67
68 if (Type_Quantity != INTEGRALQUANTITY){
69 QuantityStorage_P->FunctionSpace = FunctionSpace_P =
70 (struct FunctionSpace*)
71 List_Pointer(Problem_S.FunctionSpace,
72 DefineQuantity_P->FunctionSpaceIndex) ;
73 if (DefineQuantity_P->Type == LOCALQUANTITY)
74 Get_DofOfElement(Element, FunctionSpace_P, QuantityStorage_P,
75 DefineQuantity_P->IndexInFunctionSpace) ;
76 }
77 else{ /* INTEGRALQUANTITY */
78 if(DefineQuantity_P->IntegralQuantity.DefineQuantityIndexDof >= 0)
79 QuantityStorage_P->FunctionSpace = (struct FunctionSpace*)
80 List_Pointer(Problem_S.FunctionSpace,
81 DefineQuantity_P->FunctionSpaceIndex) ;
82 /* Get the function space for the associated local quantities */
83 for (j = 0 ; j < DefineQuantity_P->IntegralQuantity.NbrQuantityIndex ; j++) {
84 Index_DefineQuantity = DefineQuantity_P->IntegralQuantity.QuantityIndexTable[j];
85 DefineQuantity_P = DefineQuantity_P0 + Index_DefineQuantity ;
86 QuantityStorage_P = QuantityStorage_P0 + Index_DefineQuantity ;
87 QuantityStorage_P->FunctionSpace = (struct FunctionSpace*)
88 List_Pointer(Problem_S.FunctionSpace,
89 DefineQuantity_P->FunctionSpaceIndex) ;
90 }
91 }
92
93 }
94 }
95
96 /* get the jacobian */
97
98 if (Element->Num != NO_ELEMENT) {
99 if (PostQuantityTerm_P->JacobianMethodIndex >= 0) {
100 JacobianCase_L =
101 ((struct JacobianMethod *)
102 List_Pointer(Problem_S.JacobianMethod,
103 PostQuantityTerm_P->JacobianMethodIndex))
104 ->JacobianCase ;
105 JacobianCase_P0 = (struct JacobianCase*)List_Pointer(JacobianCase_L, 0);
106
107 i = 0 ;
108 while ((i < List_Nbr(JacobianCase_L)) &&
109 ((j = (JacobianCase_P0 + i)->RegionIndex) >= 0) &&
110 !List_Search
111 (((struct Group *)List_Pointer(Problem_S.Group, j))
112 ->InitialList, &Element->Region, fcmp_int) ) i++ ;
113
114 if (i == List_Nbr(JacobianCase_L)){
115 Message::Error("Undefined Jacobian in Region %d", Element->Region) ;
116 return;
117 }
118
119 Element->JacobianCase = JacobianCase_P0 + i ;
120 Get_Jacobian = (double (*)(struct Element*, MATRIX3x3*))
121 Get_JacobianFunction(Element->JacobianCase->TypeJacobian,
122 Element->Type, &Type_Dimension) ;
123 }
124 else {
125 if(!Warning_NoJacobian){
126 Message::Warning("No Jacobian method specification in PostProcessing quantity: "
127 "using default Jacobian (Vol)");
128 Warning_NoJacobian = 1 ;
129 }
130 Get_Jacobian = (double (*)(struct Element*, MATRIX3x3*))
131 Get_JacobianFunction(JACOBIAN_VOL, Element->Type, &Type_Dimension) ;
132 }
133
134 Get_NodesCoordinatesOfElement(Element) ;
135 }
136
137 /* local evaluation at one point */
138
139 if(PostQuantityTerm_P->EvaluationType == LOCAL){
140
141 if (Element->Num != NO_ELEMENT) {
142 Get_BFGeoElement(Element, u, v, w) ;
143 Element->DetJac = Get_Jacobian(Element, &Element->Jac) ;
144 if (Element->DetJac != 0.)
145 Get_InverseMatrix(Type_Dimension, Element->Type, Element->DetJac,
146 &Element->Jac, &Element->InvJac) ;
147 }
148
149 Cal_WholeQuantity
150 (Current.Element = Element,
151 QuantityStorage_P0, PostQuantityTerm_P->WholeQuantity,
152 Current.u = u, Current.v = v, Current.w = w, -1, -1, &TermValue) ;
153
154 }
155
156 /* integral evaluation over the element */
157
158 else if(PostQuantityTerm_P->EvaluationType == INTEGRAL){
159
160 if(Element->Num == NO_ELEMENT){
161 Message::Error("No element in which to integrate");
162 return;
163 }
164
165 if(PostQuantityTerm_P->IntegrationMethodIndex < 0){
166 Message::Error("Missing Integration method in PostProcesssing Quantity '%s'",
167 PostQuantity_P->Name);
168 return;
169 }
170
171 IntegrationCase_L =
172 ((struct IntegrationMethod *)
173 List_Pointer(Problem_S.IntegrationMethod,
174 PostQuantityTerm_P->IntegrationMethodIndex))->IntegrationCase ;
175
176 CriterionIndex =
177 ((struct IntegrationMethod *)
178 List_Pointer(Problem_S.IntegrationMethod,
179 PostQuantityTerm_P->IntegrationMethodIndex))
180 ->CriterionIndex ;
181
182 IntegrationCase_P = Get_IntegrationCase(Element,
183 IntegrationCase_L,
184 CriterionIndex) ;
185
186 if(IntegrationCase_P->Type != GAUSS){
187 Message::Error("Only numerical integration is available "
188 "in Integral PostQuantities");
189 return;
190 }
191
192 Quadrature_P = (struct Quadrature*)
193 List_PQuery(IntegrationCase_P->Case, &Element->Type, fcmp_int);
194
195 if(!Quadrature_P){
196 Message::Error("Unknown type of Element (%s) for Integration method (%s) "
197 " in PostProcessing Quantity (%s)",
198 Get_StringForDefine(Element_Type, Element->Type),
199 ((struct IntegrationMethod *)
200 List_Pointer(Problem_S.IntegrationMethod,
201 PostQuantityTerm_P->IntegrationMethodIndex))->Name,
202 PostQuantity_P->Name);
203 return;
204 }
205
206 Cal_ZeroValue(&TermValue);
207
208 Nbr_IntPoints = Quadrature_P->NumberOfPoints ;
209 Get_IntPoint = (void (*) (int,int,double*,double*,double*,double*))
210 Quadrature_P->Function ;
211
212 for (i_IntPoint = 0 ; i_IntPoint < Nbr_IntPoints ; i_IntPoint++) {
213
214 Current.QuadraturePointIndex = i_IntPoint;
215
216 Get_IntPoint(Nbr_IntPoints, i_IntPoint, &ui, &vi, &wi, &weight) ;
217
218 Get_BFGeoElement (Element, ui, vi, wi) ;
219 Element->DetJac = Get_Jacobian(Element, &Element->Jac) ;
220 if (Element->DetJac != 0.){
221 Get_InverseMatrix(Type_Dimension, Element->Type, Element->DetJac,
222 &Element->Jac, &Element->InvJac) ;
223 }
224 else{
225 Message::Warning("Zero determinant in 'Cal_PostQuantity'");
226 }
227 Current.x = Current.y = Current.z = 0. ;
228 if (Type_Quantity == INTEGRALQUANTITY){
229 for (i = 0 ; i < Element->GeoElement->NbrNodes ; i++) {
230 Current.x += Element->x[i] * Element->n[i] ;
231 Current.y += Element->y[i] * Element->n[i] ;
232 Current.z += Element->z[i] * Element->n[i] ;
233 }
234 }
235
236 Cal_WholeQuantity
237 (Current.Element = Element,
238 QuantityStorage_P0, PostQuantityTerm_P->WholeQuantity,
239 Current.u = ui, Current.v = vi, Current.w = wi, -1, -1, &tmpValue) ;
240
241 Factor = weight * fabs(Element->DetJac) ;
242
243 TermValue.Type = tmpValue.Type ;
244 Cal_AddMultValue(&TermValue,&tmpValue,Factor,&TermValue);
245 }
246 }
247
248 Value->Type = TermValue.Type;
249 Cal_AddValue(Value,&TermValue,Value);
250 }
251
252 /* ------------------------------------------------------------------------ */
253 /* P o s _ G l o b a l Q u a n t i t y */
254 /* ------------------------------------------------------------------------ */
255
Pos_GlobalQuantity(struct PostQuantity * PostQuantity_P,struct DefineQuantity * DefineQuantity_P0,struct QuantityStorage * QuantityStorage_P0,struct PostQuantityTerm * PostQuantityTerm_P,struct Element * ElementEmpty,List_T * InRegion_L,List_T * Support_L,struct Value * Value,int Type_InRegion)256 void Pos_GlobalQuantity(struct PostQuantity *PostQuantity_P,
257 struct DefineQuantity *DefineQuantity_P0,
258 struct QuantityStorage *QuantityStorage_P0,
259 struct PostQuantityTerm *PostQuantityTerm_P,
260 struct Element *ElementEmpty,
261 List_T *InRegion_L, List_T *Support_L,
262 struct Value *Value, int Type_InRegion)
263 {
264 struct DefineQuantity *DefineQuantity_P ;
265 struct QuantityStorage *QuantityStorage_P ;
266 struct FunctionSpace *FunctionSpace_P ;
267 struct GlobalQuantity *GlobalQuantity_P ;
268 struct Value TermValue ;
269
270 int k, Index_DefineQuantity ;
271
272 int Nbr_Element, i_Element ;
273 struct Element Element ;
274 int Type_Quantity ;
275
276 if (PostQuantityTerm_P->EvaluationType == LOCAL &&
277 List_Search(InRegion_L, &Current.Region, fcmp_int)) {
278
279 for (k = 0 ; k < PostQuantityTerm_P->NbrQuantityIndex ; k++) {
280 Index_DefineQuantity = PostQuantityTerm_P->QuantityIndexTable[k] ;
281 DefineQuantity_P = DefineQuantity_P0 + Index_DefineQuantity ;
282 QuantityStorage_P = QuantityStorage_P0 + Index_DefineQuantity ;
283
284 if (QuantityStorage_P->NumLastElementForFunctionSpace != Current.Region) {
285 QuantityStorage_P->NumLastElementForFunctionSpace = Current.Region ;
286 QuantityStorage_P->FunctionSpace = FunctionSpace_P =
287 (struct FunctionSpace*)
288 List_Pointer(Problem_S.FunctionSpace,
289 DefineQuantity_P->FunctionSpaceIndex) ;
290 GlobalQuantity_P = (struct GlobalQuantity*)
291 List_Pointer
292 (QuantityStorage_P->FunctionSpace->GlobalQuantity,
293 *(int *)List_Pointer(DefineQuantity_P->IndexInFunctionSpace, 0)) ;
294
295 if (DefineQuantity_P->Type == GLOBALQUANTITY)
296 Get_DofOfRegion(Current.Region, GlobalQuantity_P,
297 FunctionSpace_P, QuantityStorage_P) ;
298 }
299 }
300
301 Cal_WholeQuantity
302 (Current.Element = ElementEmpty,
303 QuantityStorage_P0, PostQuantityTerm_P->WholeQuantity,
304 Current.u = 0., Current.v = 0., Current.w = 0., -1, -1, &TermValue) ;
305
306 Value->Type = TermValue.Type;
307 Cal_AddValue(Value,&TermValue,Value);
308
309 } /* if LOCAL && ... */
310
311 else if (PostQuantityTerm_P->EvaluationType == INTEGRAL) {
312
313 Nbr_Element = Geo_GetNbrGeoElements() ;
314 Get_InitDofOfElement(&Element) ;
315
316 Type_Quantity = LOCALQUANTITY ; /* Attention... il faut se comprendre: */
317 /* il s'agit de grandeurs locales qui seront integrees */
318 for (i_Element = 0 ; i_Element < Nbr_Element; i_Element++) {
319 Element.GeoElement = Geo_GetGeoElement(i_Element) ;
320 Element.Num = Element.GeoElement->Num ;
321 Element.Type = Element.GeoElement->Type ;
322 Current.Region = Element.Region = Element.GeoElement->Region ;
323
324 /* Filter: only elements in both InRegion_L and Support_L are considered */
325 if ((!InRegion_L ||
326 (List_Search(InRegion_L, (Type_InRegion==ELEMENTSOF ?
327 &Element.Num : &Element.Region), fcmp_int))) &&
328 (!Support_L || List_Search(Support_L, &Element.Region, fcmp_int))) {
329
330 Get_NodesCoordinatesOfElement(&Element) ;
331 Current.x = Element.x[0];
332 Current.y = Element.y[0];
333 Current.z = Element.z[0];
334 Pos_LocalOrIntegralQuantity(PostQuantity_P,
335 DefineQuantity_P0, QuantityStorage_P0,
336 PostQuantityTerm_P, &Element, Type_Quantity,
337 0., 0., 0., Value) ;
338 }
339 } /* for i_Element ... */
340
341 } /* if INTEGRAL ... */
342
343 }
344
345 /* ------------------------------------------------------------------------ */
346 /* C a l _ P o s t Q u a n t i t y */
347 /* ------------------------------------------------------------------------ */
348
Cal_PostQuantity(struct PostQuantity * PostQuantity_P,struct DefineQuantity * DefineQuantity_P0,struct QuantityStorage * QuantityStorage_P0,List_T * Support_L,struct Element * Element,double u,double v,double w,struct Value * Value)349 void Cal_PostQuantity(struct PostQuantity *PostQuantity_P,
350 struct DefineQuantity *DefineQuantity_P0,
351 struct QuantityStorage *QuantityStorage_P0,
352 List_T *Support_L,
353 struct Element *Element,
354 double u, double v, double w,
355 struct Value *Value)
356 {
357 struct PostQuantityTerm PostQuantityTerm ;
358
359 // InRegion_L is the region (group) on which the quantity in the
360 // PostProcessing section is defined (so it's actually a "support"...)
361
362
363 List_T *InRegion_L ;
364 int i_PQT, Type_Quantity, Type_InRegion ;
365 struct Group * Group_P ;/* For generating extended group */
366
367 /* mettre tout a zero: on ne connait pas a priori le type de retour */
368 /* (default type and value returned if Type_Quantity == -1) */
369
370 Cal_ZeroValue(Value);
371 Value->Type = SCALAR;
372
373 /* Loop on PostQuantity Terms */
374 /* ... with sum of results if common supports (In ...) */
375
376 for (i_PQT = 0 ; i_PQT < List_Nbr(PostQuantity_P->PostQuantityTerm) ; i_PQT++) {
377
378 List_Read(PostQuantity_P->PostQuantityTerm, i_PQT, &PostQuantityTerm) ;
379
380 /*
381 InRegion_L = (PostQuantityTerm.InIndex < 0)? NULL :
382 ((struct Group *)List_Pointer(Problem_S.Group,
383 PostQuantityTerm.InIndex))->InitialList ;
384 */
385
386 Group_P = (PostQuantityTerm.InIndex < 0)? NULL :
387 (struct Group *)List_Pointer(Problem_S.Group,
388 PostQuantityTerm.InIndex);
389 InRegion_L = Group_P ? Group_P->InitialList : NULL ;
390
391 if (PostQuantityTerm.SubRegion >=0) {
392 struct Group * GroupSubRegion_P = (struct Group *)
393 List_Pointer(Problem_S.Group,
394 PostQuantityTerm.SubRegion);
395 if (List_Nbr(GroupSubRegion_P->InitialList) == 1) {
396 List_Read(GroupSubRegion_P->InitialList, 0, &Current.SubRegion) ;
397 }
398 else {
399 Message::Error("One region allowed in SubRegion");
400 Current.SubRegion = -1;
401 }
402 }
403 //+++ Not in pos!!! else
404 // Current.SubRegion = -1;
405
406 Type_InRegion = Group_P ? Group_P->FunctionType : REGION;
407
408 /* Generating Extended Group if necessary */
409 if (Group_P && Group_P->FunctionType == ELEMENTSOF){
410 if (!Group_P->ExtendedList) Generate_ExtendedGroup(Group_P) ;
411 InRegion_L = Group_P->ExtendedList ;
412 }
413
414 if (!Support_L) Type_Quantity = PostQuantityTerm.Type ;
415 else Type_Quantity = GLOBALQUANTITY ; /* Always if Support */
416
417 if (InRegion_L) {
418 if (Element->Num != NO_ELEMENT) {
419 /* not correct for ElementsOf (i.e. ELEMENTLIST)
420 if (!List_Search(InRegion_L, &Element->Region, fcmp_int)) {
421 Type_Quantity = -1 ;
422 }
423 */
424 if (!((Group_P->Type != ELEMENTLIST &&
425 List_Search(Group_P->InitialList, &Element->Region, fcmp_int))
426 ||
427 (Group_P->Type == ELEMENTLIST &&
428 Check_IsEntityInExtendedGroup(Group_P, Element->Num, 0))
429 )) {
430 Type_Quantity = -1 ;
431 }
432 }
433 else {
434 if (Type_Quantity == GLOBALQUANTITY) {
435
436 /* Plus de test ici... vu que le OnRegion de la PostOperation n'a rien
437 a voir avec le support d'une integration ...
438 if (!List_Search(InRegion_L, &Current.Region, fcmp_int)) {
439 Type_Quantity = -1 ;
440 }
441 */
442 /* Il faut plutot voir si il existe au moins une region de InRegion_L
443 qui soit dans Support_L ... cela est fait apres, pour chaque element */
444 }
445 else if (Type_Quantity != INTEGRALQUANTITY) {
446 Type_Quantity = -1 ;
447 }
448 }
449 }
450 /* else if !InRegion_L -> No filter, i.e. globally defined quantity */
451
452 /* ---------------------------- */
453 /* Local or Integral quantities */
454 /* ---------------------------- */
455
456 if (Type_Quantity == LOCALQUANTITY || Type_Quantity == INTEGRALQUANTITY) {
457
458 Pos_LocalOrIntegralQuantity(PostQuantity_P,
459 DefineQuantity_P0, QuantityStorage_P0,
460 &PostQuantityTerm, Element, Type_Quantity,
461 u, v, w, Value) ;
462 }
463
464 /* ----------------- */
465 /* Global quantities */
466 /* ----------------- */
467
468 else if (Type_Quantity == GLOBALQUANTITY) {
469
470 Pos_GlobalQuantity(PostQuantity_P,
471 DefineQuantity_P0, QuantityStorage_P0,
472 &PostQuantityTerm, Element, InRegion_L, Support_L, Value, Type_InRegion) ;
473 }
474
475 } /* for i_PQT ... */
476
477 }
478
479 /* ------------------------------------------------------------------------ */
480 /* C a l _ P o s t C u m u l a t i v e Q u a n t i t y */
481 /* ------------------------------------------------------------------------ */
482
Cal_PostCumulativeQuantity(List_T * Region_L,int SupportIndex,List_T * TimeStep_L,struct PostQuantity * PostQuantity_P,struct DefineQuantity * DefineQuantity_P0,struct QuantityStorage * QuantityStorage_P0,struct Value ** Values)483 void Cal_PostCumulativeQuantity(List_T *Region_L,
484 int SupportIndex,
485 List_T *TimeStep_L,
486 struct PostQuantity *PostQuantity_P,
487 struct DefineQuantity *DefineQuantity_P0,
488 struct QuantityStorage *QuantityStorage_P0,
489 struct Value **Values)
490 {
491 struct Element Element ;
492 List_T *Support_L ;
493 int i, NbrTimeStep ;
494
495 Support_L = ((struct Group *)
496 List_Pointer(Problem_S.Group, SupportIndex))->InitialList ;
497
498 NbrTimeStep = List_Nbr(TimeStep_L) ;
499 *Values = (struct Value *)Malloc(NbrTimeStep*sizeof(struct Value)) ;
500
501 Element.Num = NO_ELEMENT ;
502
503 for(i = 0 ; i < NbrTimeStep ; i++) {
504 Pos_InitAllSolutions(TimeStep_L, i) ;
505
506 Cal_PostQuantity(PostQuantity_P, DefineQuantity_P0, QuantityStorage_P0,
507 Support_L, &Element, 0, 0, 0, &(*Values)[i]) ;
508 }
509 }
510
511 /* ------------------------------------------------------------------------ */
512 /* C o m b i n e _ P o s t Q u a n t i t y */
513 /* ------------------------------------------------------------------------ */
514
Combine_PostQuantity(int Type,int Order,struct Value * V1,struct Value * V2)515 void Combine_PostQuantity(int Type, int Order,
516 struct Value *V1, struct Value *V2)
517 {
518 switch(Type){
519 case MULTIPLICATION : Cal_ProductValue(V1, V2, V1) ; break ;
520 case ADDITION : Cal_AddValue(V1, V2, V1) ; break ;
521 case DIVISION : Cal_DivideValue(Order?V1:V2, Order?V2:V1, V1) ; break;
522 case SOUSTRACTION : Cal_SubstractValue(Order?V1:V2, Order?V2:V1, V1) ; break;
523 }
524 }
525