1 // GetDP - Copyright (C) 1997-2013 P. Dular, C. Geuzaine
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 "ProData.h"
7 #include "GeoData.h"
8 #include "DofData.h"
9 #include "Get_DofOfElement.h"
10 #include "Get_ElementSource.h"
11 #include "Pre_TermOfFemEquation.h"
12 #include "Cal_GalerkinTermOfFemEquation.h"
13 #include "Cal_SmallFemTermOfFemEquation.h"
14 #include "Cal_GlobalTermOfFemEquation.h"
15 #include "Cal_AssembleTerm.h"
16 #include "Generate_Network.h"
17 #include "ExtendedGroup.h"
18 #include "Message.h"
19
20 extern struct Problem Problem_S ;
21 extern struct CurrentData Current ;
22
23 extern int TreatmentStatus ;
24
25 extern List_T * GeoData_L ;
26
27 /* ------------------------------------------------------------------------ */
28 /* C a l _ F e m G l o b a l E q u a t i o n */
29 /* ------------------------------------------------------------------------ */
30
Cal_FemGlobalEquation2(int Index_DefineQuantity,int Num_Region,struct DefineQuantity * DefineQuantity_P0,struct QuantityStorage * QuantityStorage_P0)31 struct Dof * Cal_FemGlobalEquation2(int Index_DefineQuantity, int Num_Region,
32 struct DefineQuantity * DefineQuantity_P0,
33 struct QuantityStorage * QuantityStorage_P0)
34 {
35 struct DefineQuantity * DefineQuantity_P ;
36 struct QuantityStorage * QuantityStorage_P ;
37 struct GlobalQuantity * GlobalQuantity_P ;
38 struct QuantityStorage QuaSto_S ;
39
40 DefineQuantity_P = DefineQuantity_P0 + Index_DefineQuantity ;
41 QuantityStorage_P = QuantityStorage_P0 + Index_DefineQuantity ;
42 GlobalQuantity_P = (struct GlobalQuantity*)
43 List_Pointer(QuantityStorage_P->FunctionSpace->GlobalQuantity,
44 *(int *)List_Pointer(DefineQuantity_P->IndexInFunctionSpace, 0)) ;
45 Get_DofOfRegion(Num_Region, GlobalQuantity_P,
46 QuantityStorage_P->FunctionSpace, &QuaSto_S) ;
47
48 if (QuaSto_S.NbrElementaryBasisFunction == 1){
49 return QuaSto_S.BasisFunction[0].Dof ;
50 }
51 else {
52 Message::Error( "Not 1 Dof associated with GlobalQuantity (Region #%d)", Num_Region) ;
53 return NULL ;
54 }
55 }
56
Cal_FemGlobalEquation(struct EquationTerm * EquationTerm_P,struct DefineQuantity * DefineQuantity_P0,struct QuantityStorage * QuantityStorage_P0)57 void Cal_FemGlobalEquation(struct EquationTerm * EquationTerm_P,
58 struct DefineQuantity * DefineQuantity_P0,
59 struct QuantityStorage * QuantityStorage_P0)
60 {
61 int Nbr_GlobalEquationTerm, i_GlobalEquationTerm ;
62
63 struct Constraint * Constraint_P ;
64 struct GlobalEquationTerm * GlobalEquationTerm_P ;
65
66 int Nbr_EquAndDof ;
67
68 List_T * InitialListInIndex_L, * RegionIndex_L ;
69 int Nbr_Region, i_Region, Num_Region, k ;
70
71
72 int Nbr_MCPR, i_MCPR, Nbr_CPR, i_CPR, i_Node, i_Loop, j_Branch, Num_Equ ;
73 struct MultiConstraintPerRegion * MCPR_P ;
74 struct ConstraintPerRegion * CPR_P ;
75 struct Group * Group_P ;
76 double Val[NBR_MAX_HARMONIC] ;
77
78 struct DofGlobal { int NumRegion ; struct Dof * Dof ; } ;
79 List_T * DofGlobal_Equ_L, * DofGlobal_DofNode_L, * DofGlobal_DofLoop_L ;
80 struct DofGlobal * DofGlobal_Equ, * DofGlobal_DofNode, * DofGlobal_DofLoop ;
81 struct DofGlobal DofGlobal_S, * DofGlobal_P ;
82
83 /* Liste des Regions auxquelles on associe des Equations de Type 'Network' */
84
85 RegionIndex_L = List_Create(50,50, sizeof(int)) ;
86
87 Constraint_P = (struct Constraint*)
88 List_Pointer(Problem_S.Constraint,
89 EquationTerm_P->Case.GlobalEquation.ConstraintIndex) ;
90 Nbr_MCPR = List_Nbr(Constraint_P->MultiConstraintPerRegion) ;
91 for (i_MCPR = 0 ; i_MCPR < Nbr_MCPR ; i_MCPR++) {
92 MCPR_P = (struct MultiConstraintPerRegion*)
93 List_Pointer(Constraint_P->MultiConstraintPerRegion, i_MCPR) ;
94 Nbr_CPR = List_Nbr(MCPR_P->ConstraintPerRegion) ;
95 for (i_CPR = 0 ; i_CPR < Nbr_CPR ; i_CPR++) {
96 CPR_P = (struct ConstraintPerRegion*)
97 List_Pointer(MCPR_P->ConstraintPerRegion, i_CPR) ;
98 Group_P = (struct Group *)
99 List_Pointer(Problem_S.Group, CPR_P->RegionIndex) ;
100 List_Read(Group_P->InitialList, 0, &Num_Region) ;
101 if (!List_Search(RegionIndex_L, &Num_Region, fcmp_int))
102 List_Add(RegionIndex_L, &Num_Region) ;
103 else {
104 Message::Error("2 occurences of Elementary Region #%d in Contraint '%s'",
105 Num_Region, Constraint_P->Name);
106 return;
107 }
108 }
109 }
110 Nbr_EquAndDof = List_Nbr(RegionIndex_L) ;
111 if (!Nbr_EquAndDof){
112 return ;
113 }
114
115 DofGlobal_Equ_L = List_Create(Nbr_EquAndDof, 1, sizeof(struct DofGlobal)) ;
116 DofGlobal_DofNode_L = List_Create(Nbr_EquAndDof, 1, sizeof(struct DofGlobal)) ;
117 DofGlobal_DofLoop_L = List_Create(Nbr_EquAndDof, 1, sizeof(struct DofGlobal)) ;
118
119 /* Construction des listes de Dof globaux pour Equ, DofNode, DofLoop */
120
121 Nbr_GlobalEquationTerm =
122 List_Nbr(EquationTerm_P->Case.GlobalEquation.GlobalEquationTerm) ;
123 for (i_GlobalEquationTerm = 0 ;
124 i_GlobalEquationTerm < Nbr_GlobalEquationTerm ; i_GlobalEquationTerm++) {
125 GlobalEquationTerm_P = (struct GlobalEquationTerm*)
126 List_Pointer(EquationTerm_P->Case.GlobalEquation.GlobalEquationTerm,
127 i_GlobalEquationTerm) ;
128 InitialListInIndex_L =
129 ((struct Group *)List_Pointer(Problem_S.Group,
130 GlobalEquationTerm_P->InIndex))->InitialList ;
131 Nbr_Region = List_Nbr(InitialListInIndex_L) ;
132 List_Sort(InitialListInIndex_L, fcmp_int) ;
133
134 for (i_Region = 0 ; i_Region < Nbr_Region ; i_Region++) {
135 List_Read(InitialListInIndex_L, i_Region, &Num_Region) ;
136 if (List_Search(RegionIndex_L, &Num_Region, fcmp_int)) {
137 DofGlobal_S.NumRegion = Num_Region ;
138 DofGlobal_S.Dof = Cal_FemGlobalEquation2
139 (GlobalEquationTerm_P->DefineQuantityIndexEqu, Num_Region,
140 DefineQuantity_P0, QuantityStorage_P0) ;
141 List_Add(DofGlobal_Equ_L, &DofGlobal_S) ;
142 DofGlobal_S.Dof = Cal_FemGlobalEquation2
143 (GlobalEquationTerm_P->DefineQuantityIndexNode, Num_Region,
144 DefineQuantity_P0, QuantityStorage_P0) ;
145 List_Add(DofGlobal_DofNode_L, &DofGlobal_S) ;
146 DofGlobal_S.Dof = Cal_FemGlobalEquation2
147 (GlobalEquationTerm_P->DefineQuantityIndexLoop, Num_Region,
148 DefineQuantity_P0, QuantityStorage_P0) ;
149 List_Add(DofGlobal_DofLoop_L, &DofGlobal_S) ;
150 }
151 }
152 }
153 if (List_Nbr(DofGlobal_Equ_L) != Nbr_EquAndDof) {
154 Message::Error("Incompatible number of equations with Contraint '%s' "
155 "(%d equations obtained while %d branches are defined)",
156 Constraint_P->Name, List_Nbr(DofGlobal_Equ_L), Nbr_EquAndDof);
157 return;
158 }
159
160 DofGlobal_Equ = (struct DofGlobal*)List_Pointer(DofGlobal_Equ_L , 0) ;
161 DofGlobal_DofNode = (struct DofGlobal*)List_Pointer(DofGlobal_DofNode_L, 0) ;
162 DofGlobal_DofLoop = (struct DofGlobal*)List_Pointer(DofGlobal_DofLoop_L, 0) ;
163
164 for (k = 0 ; k < List_Nbr(DofGlobal_Equ_L) ; k++) {
165 if (DofGlobal_Equ[k].Dof->Type == DOF_FIXED ||
166 DofGlobal_Equ[k].Dof->Type == DOF_LINK) {
167 if (DofGlobal_Equ[k].Dof == DofGlobal_DofNode[k].Dof)
168 DofGlobal_Equ[k].Dof = DofGlobal_DofLoop[k].Dof ;
169 else
170 DofGlobal_Equ[k].Dof = DofGlobal_DofNode[k].Dof ;
171 }
172 }
173
174 /* Construction des equations (assemblage) */
175
176 Num_Equ = 0 ;
177
178 Nbr_MCPR = List_Nbr(Constraint_P->MultiConstraintPerRegion) ;
179 for (i_MCPR = 0 ; i_MCPR < Nbr_MCPR ; i_MCPR++) {
180 MCPR_P = (struct MultiConstraintPerRegion*)
181 List_Pointer(Constraint_P->MultiConstraintPerRegion, i_MCPR) ;
182
183 if (!MCPR_P->Active)
184 MCPR_P->Active = Generate_Network(MCPR_P->Name, MCPR_P->ConstraintPerRegion) ;
185
186 for (i_Node = 0 ; i_Node < MCPR_P->Active->Case.Network.NbrNode ; i_Node++) {
187 for (j_Branch = 0 ;
188 j_Branch < MCPR_P->Active->Case.Network.NbrBranch ; j_Branch++) {
189
190 if (MCPR_P->Active->Case.Network.MatNode[i_Node][j_Branch]) {
191
192 Group_P = (struct Group *)
193 List_Pointer
194 (Problem_S.Group,
195 ((struct ConstraintPerRegion *)
196 List_Pointer(MCPR_P->ConstraintPerRegion, j_Branch))->RegionIndex) ;
197 List_Read(Group_P->InitialList, 0, &Num_Region) ;
198
199 DofGlobal_P = (struct DofGlobal*)
200 List_PQuery(DofGlobal_DofNode_L, &Num_Region, fcmp_int) ;
201
202 Val[0] =
203 (double)(MCPR_P->Active->Case.Network.MatNode[i_Node][j_Branch]) ;
204 if (Current.NbrHar > 1) {
205 Val[1] = 0. ;
206 for (k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2) {
207 Val[k] = Val[0] ; Val[k+1] = 0. ;
208 }
209 }
210 /*
211 Message::Info("Node: eq.(%d) [%d, %d], dof [%d, %d] : %.16g\n",
212 Num_Equ,
213 DofGlobal_Equ[Num_Equ].Dof->NumType, DofGlobal_Equ[Num_Equ].Dof->Entity,
214 DofGlobal_P->Dof->NumType, DofGlobal_P->Dof->Entity,
215 Val[0]
216 ) ;
217 */
218 Cal_AssembleTerm_NeverDt(DofGlobal_Equ[Num_Equ].Dof, DofGlobal_P->Dof, Val) ;
219 }
220 }
221
222 Num_Equ++ ;
223 } /* for i_Node ... */
224
225 for (i_Loop = 0 ; i_Loop < MCPR_P->Active->Case.Network.NbrLoop ; i_Loop++) {
226 for (j_Branch = 0 ;
227 j_Branch < MCPR_P->Active->Case.Network.NbrBranch ; j_Branch++) {
228
229 if (MCPR_P->Active->Case.Network.MatLoop[i_Loop][j_Branch]) {
230
231 Group_P = (struct Group *)
232 List_Pointer
233 (Problem_S.Group,
234 ((struct ConstraintPerRegion *)
235 List_Pointer(MCPR_P->ConstraintPerRegion, j_Branch))->RegionIndex) ;
236 List_Read(Group_P->InitialList, 0, &Num_Region) ;
237
238 DofGlobal_P = (struct DofGlobal*)
239 List_PQuery(DofGlobal_DofLoop_L, &Num_Region, fcmp_int) ;
240
241 Val[0] =
242 (double)(MCPR_P->Active->Case.Network.MatLoop[i_Loop][j_Branch]) ;
243 if (Current.NbrHar > 1) {
244 Val[1] = 0. ;
245 for (k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2) {
246 Val[k] = Val[0] ; Val[k+1] = 0. ;
247 }
248 }
249 /*
250 Message::Info("Loop: eq.(%d) [%d, %d], dof [%d, %d] : %.16g\n",
251 Num_Equ,
252 DofGlobal_Equ[Num_Equ].Dof->NumType, DofGlobal_Equ[Num_Equ].Dof->Entity,
253 DofGlobal_P->Dof->NumType, DofGlobal_P->Dof->Entity,
254 Val[0]
255 ) ;
256 */
257 Cal_AssembleTerm_NeverDt(DofGlobal_Equ[Num_Equ].Dof, DofGlobal_P->Dof, Val) ;
258
259 }
260 }
261
262 Num_Equ++ ;
263 } /* for i_Loop ... */
264
265 } /* for i_MCPR ... */
266
267 List_Delete(DofGlobal_Equ_L) ;
268 List_Delete(DofGlobal_DofNode_L) ; List_Delete(DofGlobal_DofLoop_L) ;
269 List_Delete(RegionIndex_L) ;
270 }
271
272 /* ------------------------------------------------------------------------ */
273 /* T r e a t m e n t _ F e m F o r m u l a t i o n */
274 /* ------------------------------------------------------------------------ */
275
Treatment_FemFormulation(struct Formulation * Formulation_P)276 void Treatment_FemFormulation(struct Formulation * Formulation_P)
277 {
278 struct Element Element ;
279
280 struct QuantityStorage * QuantityStorage_P0, * QuantityStorage_P ;
281 struct QuantityStorage QuantityStorage_S, QuantityStorage_GlobalEqu_S ;
282 struct Dof DofForNoDof_P [NBR_MAX_HARMONIC] ;
283 struct EquationTerm * EquationTerm_P0 , * EquationTerm_P ;
284 struct GlobalQuantity * GlobalQuantity_P ;
285
286 int Nbr_DefineQuantity ;
287 struct DefineQuantity * DefineQuantity_P0 , * DefineQuantity_P ;
288
289 List_T * FemLocalTermActive_L ;
290 struct FemLocalTermActive FemLocalTermActive_S ;
291 List_T * QuantityStorage_L ;
292
293 struct FemGlobalTermActive FemGlobalTermActive_S;
294
295 struct Group * GroupIn_P ;
296
297 int i, j, Nbr_Element, i_Element, Nbr_EquationTerm, i_EquTerm ;
298 int Index_DefineQuantity, TraceGroupIndex_DefineQuantity ;
299
300 List_T * InitialListInIndex_L ;
301 int Nbr_Region, i_Region, Num_Region ;
302 int i_Region_Dof = 0, Num_Region_Dof = 0, i_Region_Dof_ini = 0;
303 int i_Region_Dof_end = 0, i_Region_Dof_skip = 0;
304
305 extern struct Group * Generate_Group ;
306 extern double ** MH_Moving_Matrix ;
307
308 gMatrix A = Current.DofData->A ;
309 gVector b = Current.DofData->b ;
310
311 int Flag_Only ;
312
313 /* --------------------------------------------------------------- */
314 /* 0. Initialization of an active zone (QuantityStorage) for each */
315 /* DefineQuantity */
316 /* --------------------------------------------------------------- */
317
318 if (!(Nbr_EquationTerm = List_Nbr(Formulation_P->Equation))){
319 Message::Error("No equation in Formulation '%s'", Formulation_P->Name);
320 return;
321 }
322
323 if (!(Nbr_DefineQuantity = List_Nbr(Formulation_P->DefineQuantity))){
324 Message::Error("No Quantity in Formulation '%s'", Formulation_P->Name);
325 return;
326 }
327
328 DefineQuantity_P0 = (struct DefineQuantity*)
329 List_Pointer(Formulation_P->DefineQuantity, 0) ;
330
331 QuantityStorage_L = List_Create(Nbr_DefineQuantity, 1,
332 sizeof (struct QuantityStorage) ) ;
333
334 QuantityStorage_S.NumLastElementForFunctionSpace =
335 QuantityStorage_S.NumLastElementForDofDefinition =
336 QuantityStorage_S.NumLastElementForEquDefinition = 0 ;
337
338 for (i = 0 ; i < Nbr_DefineQuantity ; i++) {
339 QuantityStorage_S.DefineQuantity = DefineQuantity_P0 + i ;
340
341 if(QuantityStorage_S.DefineQuantity->Type == INTEGRALQUANTITY &&
342 QuantityStorage_S.DefineQuantity->IntegralQuantity.DefineQuantityIndexDof < 0){
343 QuantityStorage_S.FunctionSpace = NULL ;
344 QuantityStorage_S.TypeQuantity = VECTOR ; /* to change */
345 }
346 else{
347 QuantityStorage_S.FunctionSpace = (struct FunctionSpace*)
348 List_Pointer(Problem_S.FunctionSpace,
349 QuantityStorage_S.DefineQuantity->FunctionSpaceIndex) ;
350 QuantityStorage_S.TypeQuantity = QuantityStorage_S.FunctionSpace->Type ;
351 }
352
353 List_Add(QuantityStorage_L, &QuantityStorage_S) ;
354 }
355
356 QuantityStorage_P0 = (struct QuantityStorage*)List_Pointer(QuantityStorage_L, 0) ;
357
358 Get_InitDofOfElement(&Element) ;
359
360
361 /* --------------------------------------------------------------- */
362 /* 1. Initialization of equation terms */
363 /* --------------------------------------------------------------- */
364
365 EquationTerm_P0 = (struct EquationTerm*)List_Pointer(Formulation_P->Equation, 0) ;
366 FemLocalTermActive_L = List_Create(Nbr_EquationTerm, 1,
367 sizeof (struct FemLocalTermActive) ) ;
368
369 for (i_EquTerm = 0 ; i_EquTerm < Nbr_EquationTerm ; i_EquTerm++) {
370 List_Add(FemLocalTermActive_L, &FemLocalTermActive_S) ;
371 EquationTerm_P = EquationTerm_P0 + i_EquTerm ;
372
373 switch(EquationTerm_P->Type){
374
375 case GALERKIN :
376 EquationTerm_P->Case.LocalTerm.Active = (struct FemLocalTermActive*)
377 List_Pointer(FemLocalTermActive_L, i_EquTerm) ;
378 switch (TreatmentStatus) {
379 case STATUS_PRE :
380 Pre_InitTermOfFemEquation(EquationTerm_P, QuantityStorage_P0) ;
381 break ;
382 case STATUS_CAL :
383 Cal_InitGalerkinTermOfFemEquation(EquationTerm_P, QuantityStorage_P0,
384 &QuantityStorage_S, DofForNoDof_P) ;
385 break ;
386 }
387 break;
388
389 case GLOBALTERM :
390 switch (TreatmentStatus) {
391 case STATUS_PRE :
392 Pre_InitGlobalTermOfFemEquation(EquationTerm_P, QuantityStorage_P0) ;
393 break ;
394 }
395 break;
396
397 case GLOBALEQUATION :
398 break ;
399
400 default :
401 Message::Error("Unknown type of equation term") ;
402 break ;
403 }
404
405 }
406
407 /* ---------------------------------------------------------- */
408 /* 2. Loop on geometrical elements : */
409 /* Treatment of eventual GALERKIN terms */
410 /* --------------------------------------------------------- */
411
412 Nbr_Element = Geo_GetNbrGeoElements() ;
413
414 Message::ResetProgressMeter();
415
416 for (i_Element = 0 ; i_Element < Nbr_Element; i_Element++) {
417
418 if (Generate_Group) {
419 Element.Region = Geo_GetGeoElement(i_Element)->Region ;
420 while (i_Element < Nbr_Element &&
421 !List_Search(Generate_Group->InitialList,
422 &Element.Region, fcmp_int) ) {
423 i_Element++ ;
424 if (i_Element < Nbr_Element) Element.Region = Geo_GetGeoElement(i_Element)->Region ;
425 }
426 if (i_Element == Nbr_Element) break ;
427 }
428
429 Element.GeoElement = Geo_GetGeoElement(i_Element) ;
430 Element.Num = Element.GeoElement->Num ;
431 Element.Type = Element.GeoElement->Type ;
432 Current.Region = Element.Region = Element.GeoElement->Region ;
433
434 /* ---------------------------- */
435 /* 2.1. Loop on equation terms */
436 /* ---------------------------- */
437
438 for (i_EquTerm = 0 ; i_EquTerm < Nbr_EquationTerm ; i_EquTerm++) {
439 EquationTerm_P = EquationTerm_P0 + i_EquTerm ;
440
441 if (EquationTerm_P->Type == GALERKIN) {
442
443 /* if the element is in the support of integration of the term */
444 /* ----------------------------------------------------------- */
445
446 GroupIn_P = (struct Group *)
447 List_Pointer(Problem_S.Group,
448 EquationTerm_P->Case.LocalTerm.InIndex);
449
450
451 if ((GroupIn_P->Type != ELEMENTLIST &&
452 List_Search(GroupIn_P->InitialList, &Element.Region, fcmp_int))
453 ||
454 (GroupIn_P->Type == ELEMENTLIST &&
455 Check_IsEntityInExtendedGroup(GroupIn_P, Element.Num, 0))
456 ) {
457
458 /*
459 if (List_Search(((struct Group *)
460 List_Pointer(Problem_S.Group,
461 EquationTerm_P->Case.
462 LocalTerm.InIndex))->InitialList,
463 &Element.Region, fcmp_int ) ) {
464 */
465
466 if(Message::GetVerbosity() == 10)
467 printf("==> Element #%d, EquationTerm #%d/%d\n",
468 Element.Num, i_EquTerm+1, Nbr_EquationTerm) ;
469
470 Current.IntegrationSupportIndex = EquationTerm_P->Case.LocalTerm.InIndex ;
471
472 if (EquationTerm_P->Case.LocalTerm.SubRegion >=0) {
473 struct Group * GroupSubRegion_P = (struct Group *)
474 List_Pointer(Problem_S.Group,
475 EquationTerm_P->Case.LocalTerm.SubRegion);
476 if (List_Nbr(GroupSubRegion_P->InitialList) == 1) {
477 List_Read(GroupSubRegion_P->InitialList, 0, &Current.SubRegion) ;
478 }
479 else {
480 Message::Error("One region allowed in SubRegion");
481 Current.SubRegion = -1;
482 }
483 }
484 else
485 Current.SubRegion = -1;
486
487
488 /* ---------------------------------------------------------- */
489 /* 2.1.1. Loop on quantities (test fcts and shape functions) */
490 /* ---------------------------------------------------------- */
491
492 for (i = 0 ; i < EquationTerm_P->Case.LocalTerm.Term.NbrQuantityIndex ; i++) {
493
494 Index_DefineQuantity =
495 EquationTerm_P->Case.LocalTerm.Term.QuantityIndexTable[i] ;
496 DefineQuantity_P = DefineQuantity_P0 + Index_DefineQuantity ;
497 QuantityStorage_P = QuantityStorage_P0 + Index_DefineQuantity ;
498
499 TraceGroupIndex_DefineQuantity =
500 EquationTerm_P->Case.LocalTerm.Term.QuantityTraceGroupIndexTable[i] ;
501
502 /* Only one analysis for each function space */
503
504 /*
505 * Attention : l'operateur de trace ne fonctionne que si le champ
506 * dont on prend la trace n'intervient qu'une seule fois dans le terme.
507 * du a - manque de generalite du code au niveau de la gestion des
508 * espaces fonctionnels des fcts tests pour 'Trace de Dof'
509 * - et Christophe fatigue
510 */
511
512 if (QuantityStorage_P->NumLastElementForFunctionSpace != Element.Num ||
513 TraceGroupIndex_DefineQuantity >= 0 ) {
514
515 QuantityStorage_P->NumLastElementForFunctionSpace = Element.Num ;
516
517 switch (DefineQuantity_P->Type) {
518 case LOCALQUANTITY :
519 if(TraceGroupIndex_DefineQuantity >= 0){
520 Get_ElementTrace(&Element, TraceGroupIndex_DefineQuantity) ;
521 QuantityStorage_P->NumLastElementForFunctionSpace = Element.ElementTrace->Num ;
522 Get_DofOfElement
523 (Element.ElementTrace, QuantityStorage_P->FunctionSpace, QuantityStorage_P,
524 DefineQuantity_P->IndexInFunctionSpace) ;
525 }
526 else{
527 Get_DofOfElement
528 (&Element, QuantityStorage_P->FunctionSpace, QuantityStorage_P,
529 DefineQuantity_P->IndexInFunctionSpace) ;
530 }
531 break ;
532 case INTEGRALQUANTITY :
533 QuantityStorage_P->NbrElementaryBasisFunction = 0 ;
534 break ;
535 default :
536 Message::Error("Bad kind of Quantity in Formulation '%s'",
537 Formulation_P->Name);
538 break;
539 }
540 }
541 } /* for i = 0, 1 ... */
542
543 /* -------------------------------------- */
544 /* 2.1.2. Treatment of the equation term */
545 /* -------------------------------------- */
546
547 switch (TreatmentStatus) {
548 case STATUS_PRE :
549 Pre_TermOfFemEquation(&Element, EquationTerm_P, QuantityStorage_P0) ;
550 break ;
551 case STATUS_CAL :
552 Flag_Only = 0 ;
553
554 if (Current.DofData->Flag_Only){
555 A = Current.DofData->A ;
556 b = Current.DofData->b ;
557
558 if (EquationTerm_P->Case.LocalTerm.MatrixIndex == -1)
559 EquationTerm_P->Case.LocalTerm.MatrixIndex = 0 ;
560
561 j = List_ISearch(Current.DofData->OnlyTheseMatrices,
562 &EquationTerm_P->Case.LocalTerm.MatrixIndex, fcmp_int) ;
563 if(j!=-1){
564 Flag_Only = 1 ;
565 switch(EquationTerm_P->Case.LocalTerm.MatrixIndex){
566 case 1 :
567 Current.DofData->A = Current.DofData->A1 ;
568 Current.DofData->b = Current.DofData->b1 ;
569 break;
570 case 2 :
571 Current.DofData->A = Current.DofData->A2 ;
572 Current.DofData->b = Current.DofData->b2 ;
573 break;
574 case 3 :
575 Current.DofData->A = Current.DofData->A3 ;
576 Current.DofData->b = Current.DofData->b3 ;
577 break;
578 }
579 }
580 }/* Only the matrices that vary are recalculated */
581
582 if (!Current.DofData->Flag_Only || (Current.DofData->Flag_Only && Flag_Only) ){
583 QuantityStorage_P = QuantityStorage_P0 +
584 EquationTerm_P->Case.LocalTerm.Term.DefineQuantityIndexEqu ;
585
586 if(EquationTerm_P->Type == GALERKIN){
587 #if defined(HAVE_SMALLFEM)
588 Cal_SmallFemTermOfFemEquation(&Element, EquationTerm_P, QuantityStorage_P0) ;
589 #else
590 Cal_GalerkinTermOfFemEquation(&Element, EquationTerm_P, QuantityStorage_P0) ;
591 #endif
592 }
593 if (Current.DofData->Flag_Only && Flag_Only){
594 Current.DofData->A = A ;
595 Current.DofData->b = b ;
596 }
597
598 }/* Flag_Only */
599 break ;
600
601 case STATUS_CST :
602 Cst_TermOfFemEquation(&Element, EquationTerm_P, QuantityStorage_P0) ;
603 break ;
604
605 }
606
607 }/* if Support ... */
608
609 } /* if GALERKIN ... */
610
611 } /* for i_EquTerm ... */
612
613 Message::ProgressMeter(i_Element + 1, Nbr_Element, (TreatmentStatus == STATUS_PRE) ?
614 "Pre-processing" : "Processing (Generate)");
615 if(Message::GetErrorCount()) break;
616 } /* for i_Element ... */
617
618
619 if (MH_Moving_Matrix) {
620 List_Delete(FemLocalTermActive_L) ; List_Delete(QuantityStorage_L) ;
621 return;
622 }
623
624 /* ------------------------------------------------------ */
625 /* 3. Loop on equation terms : */
626 /* Treatment of eventual GLOBAL terms */
627 /* ------------------------------------------------------ */
628
629 for (i_EquTerm = 0 ; i_EquTerm < Nbr_EquationTerm ; i_EquTerm++) {
630
631 EquationTerm_P = EquationTerm_P0 + i_EquTerm ;
632
633 if (EquationTerm_P->Type == GLOBALTERM) {
634
635 EquationTerm_P->Case.GlobalTerm.Active = &FemGlobalTermActive_S;
636
637 InitialListInIndex_L =
638 ((struct Group *)List_Pointer(Problem_S.Group,
639 EquationTerm_P->Case.GlobalTerm.InIndex))
640 ->InitialList ;
641 List_Sort(InitialListInIndex_L, fcmp_int) ;
642 Nbr_Region = List_Nbr(InitialListInIndex_L) ;
643
644 /* ---------------------------------------------- */
645 /* 3.1. Loop on Regions belonging to the support */
646 /* ---------------------------------------------- */
647
648 for (i_Region = 0 ; i_Region < Nbr_Region ; i_Region++) {
649 List_Read(InitialListInIndex_L, i_Region, &Num_Region) ;
650 Current.Region = Num_Region ;
651
652 switch (EquationTerm_P->Case.GlobalTerm.SubType) {
653 case EQ_ST_SELF:
654 i_Region_Dof_ini = i_Region; i_Region_Dof_end = i_Region+1;
655 i_Region_Dof_skip = -1;
656 break;
657 case EQ_ST_MUTUAL:
658 i_Region_Dof_ini = 0; i_Region_Dof_end = Nbr_Region;
659 i_Region_Dof_skip = i_Region;
660 break;
661 case EQ_ST_SELFANDMUTUAL:
662 i_Region_Dof_ini = 0; i_Region_Dof_end = Nbr_Region;
663 i_Region_Dof_skip = -1;
664 break;
665 }
666
667 // Possible mutual terms (need of double-piece-wise Function fct[r1,r2])
668 for (i_Region_Dof = i_Region_Dof_ini; i_Region_Dof < i_Region_Dof_end;
669 i_Region_Dof++) {
670
671 if (i_Region_Dof != i_Region_Dof_skip) {
672
673 List_Read(InitialListInIndex_L, i_Region_Dof, &Num_Region_Dof) ;
674 Current.SubRegion = Num_Region_Dof; // used in double-piece-wise Function
675
676 /* ---------------------------------------------------------------- */
677 /* 3.1.1. Loop on Quantities (test functions and shape functions) */
678 /* ---------------------------------------------------------------- */
679
680 for (i = 0 ; i < EquationTerm_P->Case.GlobalTerm.Term.NbrQuantityIndex ; i++) {
681
682 Index_DefineQuantity =
683 EquationTerm_P->Case.GlobalTerm.Term.QuantityIndexTable[i] ;
684 DefineQuantity_P = DefineQuantity_P0 + Index_DefineQuantity ;
685 QuantityStorage_P = QuantityStorage_P0 + Index_DefineQuantity ;
686 GlobalQuantity_P = (struct GlobalQuantity*)
687 List_Pointer(QuantityStorage_P->FunctionSpace->GlobalQuantity,
688 *(int*)List_Pointer(DefineQuantity_P->IndexInFunctionSpace, 0)) ;
689
690 /* Only one Function space analysis */
691 /* -------------------------------- */
692 if (QuantityStorage_P->NumLastElementForFunctionSpace != Num_Region_Dof) {
693 QuantityStorage_P->NumLastElementForFunctionSpace = Num_Region_Dof ;
694
695 switch (DefineQuantity_P->Type) {
696 case GLOBALQUANTITY :
697 Get_DofOfRegion
698 (Num_Region_Dof, GlobalQuantity_P,
699 QuantityStorage_P->FunctionSpace, QuantityStorage_P) ;
700 break ;
701 default :
702 Message::Error("Bad kind of Quantity in Formulation '%s'",
703 Formulation_P->Name);
704 break;
705 }
706 }
707
708 // Particular QuantityStorage for Equ
709 if (Num_Region != Num_Region_Dof &&
710 Index_DefineQuantity ==
711 EquationTerm_P->Case.GlobalTerm.Term.DefineQuantityIndexEqu) {
712 switch (DefineQuantity_P->Type) {
713 case GLOBALQUANTITY :
714 Get_DofOfRegion
715 (Num_Region, GlobalQuantity_P,
716 QuantityStorage_P->FunctionSpace, &QuantityStorage_GlobalEqu_S) ;
717 break ;
718 default :
719 Message::Error("Bad kind of Quantity in Formulation '%s'",
720 Formulation_P->Name);
721 break;
722 }
723 }
724
725 } /* for i = 0, 1 ... */
726
727 // QuantityStorage for Equ and Dof (can differ for SubType Mutual)
728 EquationTerm_P->Case.GlobalTerm.Active->QuantityStorageEqu_P =
729 (Num_Region == Num_Region_Dof)?
730 QuantityStorage_P0 +
731 EquationTerm_P->Case.GlobalTerm.Term.DefineQuantityIndexEqu
732 :
733 &QuantityStorage_GlobalEqu_S
734 ;
735
736 EquationTerm_P->Case.GlobalTerm.Active->flag_Dof =
737 (EquationTerm_P->Case.GlobalTerm.Term.DefineQuantityIndexDof >= 0);
738
739 EquationTerm_P->Case.GlobalTerm.Active->QuantityStorageDof_P =
740 (EquationTerm_P->Case.GlobalTerm.Active->flag_Dof)?
741 QuantityStorage_P0 +
742 EquationTerm_P->Case.GlobalTerm.Term.DefineQuantityIndexDof
743 : NULL;
744
745 /* ------------------------------ */
746 /* 3.1.2. Treatment of the term */
747 /* ------------------------------ */
748
749 switch (TreatmentStatus) {
750 case STATUS_PRE :
751 Pre_GlobalTermOfFemEquation(Num_Region, EquationTerm_P,
752 QuantityStorage_P0) ;
753 break ;
754 case STATUS_CAL :
755 Cal_GlobalTermOfFemEquation(Num_Region, EquationTerm_P,
756 QuantityStorage_P0,
757 &QuantityStorage_S, DofForNoDof_P) ;
758 break ;
759 case STATUS_CST :
760 Cst_GlobalTermOfFemEquation(Num_Region, EquationTerm_P,
761 QuantityStorage_P0) ;
762 break ;
763 }
764
765 }
766 } // for i_Region_Dof
767 } // for i_Region
768
769 } /* if GLOBALTERM ... */
770 } /* for i_EquTerm ... */
771
772
773 /* --------------------------------------------------------- */
774 /* 4. Loop on equation terms : */
775 /* Treatment of eventual GLOBAL EQUATION terms */
776 /* --------------------------------------------------------- */
777
778 for (i_EquTerm = 0 ; i_EquTerm < Nbr_EquationTerm ; i_EquTerm++) {
779 EquationTerm_P = EquationTerm_P0 + i_EquTerm ;
780
781 if (EquationTerm_P->Type == GLOBALEQUATION) {
782
783 if (EquationTerm_P->Case.GlobalEquation.ConstraintIndex >= 0)
784
785 switch (TreatmentStatus) {
786 case STATUS_PRE :
787 Pre_FemGlobalEquation(EquationTerm_P,
788 DefineQuantity_P0, QuantityStorage_P0) ;
789 break ;
790 case STATUS_CAL :
791 Cal_FemGlobalEquation(EquationTerm_P,
792 DefineQuantity_P0, QuantityStorage_P0) ;
793 break ;
794 }
795
796 } /* if GLOBALEQUATION ... */
797 } /* for i_EquTerm ... */
798
799
800 /* -------------------------- */
801 /* 5. End of FEM treatment */
802 /* -------------------------- */
803
804 List_Delete(FemLocalTermActive_L) ; List_Delete(QuantityStorage_L) ;
805 Cal_EndGalerkinTermOfFemEquation();
806 }
807
808 /* ------------------------------------------------------------------------ */
809 /* T r e a t m e n t _ G l o b a l F o r m u l a t i o n */
810 /* ------------------------------------------------------------------------ */
811
Treatment_GlobalFormulation(struct Formulation * Formulation_P)812 void Treatment_GlobalFormulation(struct Formulation * Formulation_P)
813 {
814 Message::Error("You should not be here!") ;
815 }
816
817 /* ------------------------------------------------------------------------ */
818 /* T r e a t m e n t _ F o r m u l a t i o n */
819 /* ------------------------------------------------------------------------ */
820
Treatment_Formulation(struct Formulation * Formulation_P)821 void Treatment_Formulation(struct Formulation * Formulation_P)
822 {
823 switch (Formulation_P->Type) {
824
825 case FEMEQUATION :
826 Treatment_FemFormulation(Formulation_P) ;
827 break ;
828
829 case GLOBALEQUATION :
830 Treatment_GlobalFormulation(Formulation_P) ;
831 break ;
832
833 default :
834 Message::Error("Unknown type for Formulation '%s'", Formulation_P->Name) ;
835 break ;
836 }
837 }
838