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 // Johan Gyselinck
8 // Ruth Sabariego
9 //
10
11 #include <map>
12 #include <complex>
13 #include <string.h>
14 #include <math.h>
15 #include "GetDPVersion.h"
16 #include "ProData.h"
17 #include "DofData.h"
18 #include "GeoData.h"
19 #include "ExtendedGroup.h"
20 #include "ListUtils.h"
21 #include "TreeUtils.h"
22 #include "MallocUtils.h"
23 #include "Message.h"
24 #include "ProParser.h"
25 #include "OS.h"
26
27 #define TWO_PI 6.2831853071795865
28
29 extern struct Problem Problem_S ;
30 extern struct CurrentData Current ;
31
32 FILE * File_PRE = 0, * File_RES = 0, * File_TMP = 0 ;
33
34 struct DofData * CurrentDofData ;
35
fcmp_Dof(const void * a,const void * b)36 int fcmp_Dof(const void * a, const void * b)
37 {
38 int Result ;
39
40 if((Result = ((struct Dof *)a)->NumType - ((struct Dof *)b)->NumType) != 0)
41 return Result ;
42 if((Result = ((struct Dof *)a)->Entity - ((struct Dof *)b)->Entity) != 0)
43 return Result ;
44 return ((struct Dof *)a)->Harmonic - ((struct Dof *)b)->Harmonic ;
45 }
46
47 /* ------------------------------------------------------------------------ */
48 /* D o f _ I n i t D o f D a t a */
49 /* ------------------------------------------------------------------------ */
50
Dof_InitDofData(struct DofData * DofData_P,int Num,int ResolutionIndex,int SystemIndex,char * Name_SolverDataFile)51 void Dof_InitDofData(struct DofData * DofData_P, int Num,
52 int ResolutionIndex, int SystemIndex,
53 char * Name_SolverDataFile)
54 {
55 int Index ;
56
57 DofData_P->Num = Num ;
58
59 DofData_P->ResolutionIndex = ResolutionIndex ;
60 DofData_P->SystemIndex = SystemIndex ;
61
62 DofData_P->FunctionSpaceIndex = NULL ;
63 DofData_P->TimeFunctionIndex = List_Create(10, 5, sizeof(int)) ;
64 Index = 0 ; List_Add(DofData_P->TimeFunctionIndex, &Index) ;
65 DofData_P->Pulsation = NULL ; DofData_P->Val_Pulsation = NULL ;
66 DofData_P->NbrHar = 1 ;
67
68 DofData_P->NbrAnyDof = 0 ; DofData_P->NbrDof = 0 ;
69 DofData_P->DofTree = Tree_Create(sizeof(struct Dof), fcmp_Dof) ;
70 DofData_P->DofList = NULL ;
71
72 DofData_P->SolverDataFileName = Name_SolverDataFile ;
73 DofData_P->Flag_RHS = 0 ;
74 DofData_P->Flag_ListOfRHS = 0 ;
75 DofData_P->CounterOfRHS = 0 ;
76 DofData_P->TotalNumberOfRHS = 0 ;
77
78 DofData_P->Flag_Init[0] = 0 ;
79 DofData_P->Flag_Init[1] = 0 ;
80 DofData_P->Flag_Init[2] = 0 ;
81 DofData_P->Flag_Init[3] = 0 ;
82 DofData_P->Flag_Init[4] = 0 ;
83 DofData_P->Flag_Init[5] = 0 ;
84 DofData_P->Flag_Init[6] = 0 ;
85 DofData_P->Flag_Init[7] = 0 ;
86
87 DofData_P->Flag_Only = 0 ;
88 DofData_P->Flag_InitOnly[0] = 0 ;
89 DofData_P->Flag_InitOnly[1] = 0 ;
90 DofData_P->Flag_InitOnly[2] = 0 ;
91
92 DofData_P->OnlyTheseMatrices = NULL;
93
94 DofData_P->Solutions = NULL ;
95 DofData_P->CurrentSolution = NULL ;
96
97 DofData_P->CorrectionSolutions.Flag = 0;
98 DofData_P->CorrectionSolutions.AllSolutions = NULL;
99
100 DofData_P->DummyDof = NULL ;
101 }
102
103 /* ------------------------------------------------------------------------ */
104 /* D o f _ F r e e D o f D a t a */
105 /* ------------------------------------------------------------------------ */
106
Dof_FreeDofData(struct DofData * DofData_P)107 void Dof_FreeDofData(struct DofData * DofData_P)
108 {
109 Message::Debug("Freeing DofData %d", DofData_P->Num);
110
111 List_Delete(DofData_P->FunctionSpaceIndex);
112 List_Delete(DofData_P->TimeFunctionIndex);
113 List_Delete(DofData_P->Pulsation);
114 Tree_Delete(DofData_P->DofTree);
115 List_Delete(DofData_P->DofList);
116 Free(DofData_P->DummyDof);
117
118 if(DofData_P->Solutions){
119 for(int i = 0; i < List_Nbr(DofData_P->Solutions); i++){
120 Solution *Solution_P = (struct Solution*)List_Pointer(DofData_P->Solutions, i);
121 if(Solution_P->SolutionExist){
122 LinAlg_DestroyVector(&Solution_P->x);
123 Free(Solution_P->TimeFunctionValues) ;
124 Solution_P->TimeFunctionValues = NULL;
125 Solution_P->SolutionExist = 0;
126 }
127 }
128 List_Delete(DofData_P->Solutions);
129 }
130
131 List_Delete(DofData_P->OnlyTheseMatrices);
132
133 if(DofData_P->Flag_Init[0] == 1 || DofData_P->Flag_Init[0] == 2){
134 LinAlg_DestroyMatrix(&DofData_P->A);
135 LinAlg_DestroyVector(&DofData_P->b);
136 LinAlg_DestroyVector(&DofData_P->res);
137 LinAlg_DestroyVector(&DofData_P->dx);
138 LinAlg_DestroySolver(&DofData_P->Solver);
139 }
140
141 if(DofData_P->Flag_Init[0] == 2){
142 LinAlg_DestroyMatrix(&DofData_P->Jac);
143 }
144
145 if(DofData_P->Flag_Init[1] == 1){
146 LinAlg_DestroyMatrix(&DofData_P->M1);
147 LinAlg_DestroyVector(&DofData_P->m1);
148 for(int i = 0; i < List_Nbr(DofData_P->m1s); i++)
149 LinAlg_DestroyVector((gVector*)List_Pointer(DofData_P->m1s, i));
150 List_Delete(DofData_P->m1s);
151 }
152
153 if(DofData_P->Flag_Init[2] == 1){
154 LinAlg_DestroyMatrix(&DofData_P->M2);
155 LinAlg_DestroyVector(&DofData_P->m2);
156 for(int i = 0; i < List_Nbr(DofData_P->m2s); i++)
157 LinAlg_DestroyVector((gVector*)List_Pointer(DofData_P->m2s, i));
158 List_Delete(DofData_P->m2s);
159 }
160
161 if(DofData_P->Flag_Init[3] == 1){
162 LinAlg_DestroyMatrix(&DofData_P->M3);
163 LinAlg_DestroyVector(&DofData_P->m3);
164 for(int i = 0; i < List_Nbr(DofData_P->m3s); i++)
165 LinAlg_DestroyVector((gVector*)List_Pointer(DofData_P->m3s, i));
166 List_Delete(DofData_P->m3s);
167 }
168
169 if(DofData_P->Flag_Init[4] == 1){
170 LinAlg_DestroyMatrix(&DofData_P->M4);
171 LinAlg_DestroyVector(&DofData_P->m4);
172 for(int i = 0; i < List_Nbr(DofData_P->m4s); i++)
173 LinAlg_DestroyVector((gVector*)List_Pointer(DofData_P->m4s, i));
174 List_Delete(DofData_P->m4s);
175 }
176
177 if(DofData_P->Flag_Init[5] == 1){
178 LinAlg_DestroyMatrix(&DofData_P->M5);
179 LinAlg_DestroyVector(&DofData_P->m5);
180 for(int i = 0; i < List_Nbr(DofData_P->m5s); i++)
181 LinAlg_DestroyVector((gVector*)List_Pointer(DofData_P->m5s, i));
182 List_Delete(DofData_P->m5s);
183 }
184
185 if(DofData_P->Flag_Init[6] == 1){
186 LinAlg_DestroyMatrix(&DofData_P->M6);
187 LinAlg_DestroyVector(&DofData_P->m6);
188 for(int i = 0; i < List_Nbr(DofData_P->m6s); i++)
189 LinAlg_DestroyVector((gVector*)List_Pointer(DofData_P->m6s, i));
190 List_Delete(DofData_P->m6s);
191 }
192 //nleigchange
193 if(DofData_P->Flag_Init[7] == 1){
194 LinAlg_DestroyMatrix(&DofData_P->M7);
195 LinAlg_DestroyVector(&DofData_P->m7);
196 for(int i = 0; i < List_Nbr(DofData_P->m7s); i++)
197 LinAlg_DestroyVector((gVector*)List_Pointer(DofData_P->m7s, i));
198 List_Delete(DofData_P->m7s);
199 }
200
201 if(DofData_P->Flag_Only){
202 if(DofData_P->Flag_InitOnly[0] == 1){
203 LinAlg_DestroyMatrix(&DofData_P->A1);
204 LinAlg_DestroyVector(&DofData_P->b1);
205 }
206 if(DofData_P->Flag_InitOnly[1] == 1){
207 LinAlg_DestroyMatrix(&DofData_P->A2);
208 LinAlg_DestroyVector(&DofData_P->b2);
209 }
210 if(DofData_P->Flag_InitOnly[2] == 1){
211 LinAlg_DestroyMatrix(&DofData_P->A3);
212 LinAlg_DestroyVector(&DofData_P->b3);
213 }
214 }
215
216 // TODO: handle MH data and CorrectionSolutions
217 }
218
219 /* ------------------------------------------------------------------------ */
220 /* D o f _ S e t C u r r e n t D o f D a t a */
221 /* ------------------------------------------------------------------------ */
222
Dof_SetCurrentDofData(struct DofData * DofData_P)223 void Dof_SetCurrentDofData(struct DofData * DofData_P)
224 {
225 CurrentDofData = DofData_P ;
226 }
227
228 /* ------------------------------------------------------------------------ */
229 /* F i l e s . . . */
230 /* ------------------------------------------------------------------------ */
231
232 /* ------------------------------------------------------------------------ */
233 /* D o f _ O p e n F i l e */
234 /* ------------------------------------------------------------------------ */
235
Dof_OpenFile(int Type,char * Name,const char * Mode)236 void Dof_OpenFile(int Type, char * Name, const char * Mode)
237 {
238 if((Message::GetIsCommWorld() && Message::GetCommRank()) &&
239 (Mode[0] == 'w' || Mode[0] == 'a')){
240 switch (Type) {
241 case DOF_PRE : File_PRE = 0 ; break ;
242 case DOF_RES : File_RES = 0 ; break ;
243 case DOF_TMP : File_RES = 0 ; break ;
244 default : break ;
245 }
246 return;
247 }
248
249 const char * Extension;
250 char FileName[256] ;
251 FILE * File_X ;
252
253 switch (Type) {
254 case DOF_PRE : Extension = ".pre" ; break ;
255 case DOF_RES : Extension = "" ; break ;
256 case DOF_TMP : Extension = "" ; break ;
257 default : Extension = ".pre" ; break ;
258 }
259
260 strcpy(FileName, Name) ; strcat(FileName, Extension) ;
261
262 if(!(File_X = FOpen(FileName, Mode)))
263 Message::Error("Unable to open file '%s'", FileName) ;
264
265 switch (Type) {
266 case DOF_PRE : File_PRE = File_X ; break ;
267 case DOF_RES : File_RES = File_X ; break ;
268 case DOF_TMP : File_TMP = File_RES ; File_RES = File_X ; break ;
269 default : break ;
270 }
271 }
272
273 /* ------------------------------------------------------------------------ */
274 /* D o f _ C l o s e F i l e */
275 /* ------------------------------------------------------------------------ */
276
Dof_CloseFile(int Type)277 void Dof_CloseFile(int Type)
278 {
279 switch (Type) {
280 case DOF_PRE : if(File_PRE) fclose(File_PRE) ; break ;
281 case DOF_RES : if(File_RES) fclose(File_RES) ; break ;
282 case DOF_TMP : if(File_RES) fclose(File_RES) ; File_RES = File_TMP ; break ;
283 }
284 }
285
286 /* ------------------------------------------------------------------------ */
287 /* D o f _ F l u s h F i l e */
288 /* ------------------------------------------------------------------------ */
289
Dof_FlushFile(int Type)290 void Dof_FlushFile(int Type)
291 {
292 switch (Type) {
293 case DOF_PRE : fflush(File_PRE) ; break ;
294 case DOF_RES : fflush(File_RES) ; break ;
295 }
296 }
297
298 /* ------------------------------------------------------------------------ */
299 /* D o f _ W r i t e F i l e P R E 0 */
300 /* ------------------------------------------------------------------------ */
301
Dof_WriteFilePRE0(int Num_Resolution,char * Name_Resolution,int Nbr_DofData)302 void Dof_WriteFilePRE0(int Num_Resolution, char * Name_Resolution,
303 int Nbr_DofData)
304 {
305 if(Message::GetIsCommWorld() && Message::GetCommRank()) return;
306
307 fprintf(File_PRE, "$Resolution /* '%s' */\n", Name_Resolution) ;
308 fprintf(File_PRE, "%d %d\n", Num_Resolution, Nbr_DofData) ;
309 fprintf(File_PRE, "$EndResolution\n") ;
310 }
311
312 /* ------------------------------------------------------------------------ */
313 /* D o f _ R e a d F i l e P R E 0 */
314 /* ------------------------------------------------------------------------ */
315
Dof_ReadFilePRE0(int * Num_Resolution,int * Nbr_DofData)316 void Dof_ReadFilePRE0(int * Num_Resolution, int * Nbr_DofData)
317 {
318 Message::Barrier();
319
320 char String[256] ;
321
322 do {
323 if(!fgets(String, sizeof(String), File_PRE)) break;
324 if(feof(File_PRE)) break;
325 } while (String[0] != '$');
326
327 if(feof(File_PRE)){
328 Message::Error("$Resolution field not found in file");
329 return;
330 }
331
332 if(!strncmp(&String[1], "Resolution", 10)) {
333 if(fscanf(File_PRE, "%d %d", Num_Resolution, Nbr_DofData) != 2){
334 Message::Error("Could not read resolution data");
335 return;
336 }
337 }
338
339 do {
340 if(!fgets(String, sizeof(String), File_PRE) || feof(File_PRE)){
341 Message::Error("Premature end of file");
342 break;
343 }
344 } while (String[0] != '$') ;
345 }
346
347 /* ------------------------------------------------------------------------ */
348 /* D o f _ W r i t e F i l e P R E */
349 /* ------------------------------------------------------------------------ */
350
Dof_WriteFilePRE(struct DofData * DofData_P)351 void Dof_WriteFilePRE(struct DofData * DofData_P)
352 {
353 if(Message::GetIsCommWorld() && Message::GetCommRank()) return;
354
355 int i, Nbr_Index ;
356 struct Dof * Dof_P0 ;
357
358 fprintf(File_PRE, "$DofData /* #%d */\n", DofData_P->Num) ;
359
360 fprintf(File_PRE, "%d %d\n",
361 DofData_P->ResolutionIndex, DofData_P->SystemIndex) ;
362
363 Nbr_Index = List_Nbr(DofData_P->FunctionSpaceIndex) ;
364 fprintf(File_PRE, "%d", Nbr_Index) ;
365 for(i = 0 ; i < Nbr_Index ; i++)
366 fprintf(File_PRE, " %d",
367 *((int *)List_Pointer(DofData_P->FunctionSpaceIndex, i))) ;
368 fprintf(File_PRE, "\n") ;
369
370 Nbr_Index = List_Nbr(DofData_P->TimeFunctionIndex) ;
371 fprintf(File_PRE, "%d", Nbr_Index) ;
372 for(i = 0 ; i < Nbr_Index ; i++)
373 fprintf(File_PRE, " %d",
374 *((int *)List_Pointer(DofData_P->TimeFunctionIndex, i))) ;
375 fprintf(File_PRE, "\n") ;
376
377 fprintf(File_PRE, "%d", 1) ;
378 for(i = 0 ; i < 1 ; i++)
379 fprintf(File_PRE, " %d", 0) ;
380 fprintf(File_PRE, "\n") ;
381
382 fprintf(File_PRE, "%d %d\n",
383 (DofData_P->DofTree)? Tree_Nbr(DofData_P->DofTree) : DofData_P->NbrAnyDof,
384 DofData_P->NbrDof) ;
385
386 if(DofData_P->DofTree)
387 Tree_Action(DofData_P->DofTree, Dof_WriteDofPRE) ;
388 else {
389 if(DofData_P->NbrAnyDof){
390 Dof_P0 = (struct Dof *)List_Pointer(DofData_P->DofList, 0) ;
391 for(i = 0 ; i < DofData_P->NbrAnyDof ; i++)
392 Dof_WriteDofPRE(Dof_P0 + i, NULL) ;
393 }
394 }
395 fprintf(File_PRE, "$EndDofData\n") ;
396 fflush(File_PRE) ;
397 }
398
399 /* ------------------------------- */
400 /* D o f _ W r i t e D o f P R E */
401 /* ------------------------------- */
402
Dof_WriteDofPRE(void * a,void * b)403 void Dof_WriteDofPRE(void * a, void * b)
404 {
405 struct Dof * Dof_P ;
406
407 Dof_P = (struct Dof *) a ;
408
409 fprintf(File_PRE, "%d %d %d %d ",
410 Dof_P->NumType, Dof_P->Entity, Dof_P->Harmonic, Dof_P->Type) ;
411
412 switch (Dof_P->Type) {
413 case DOF_UNKNOWN :
414 fprintf(File_PRE, "%d %d\n", Dof_P->Case.Unknown.NumDof,
415 Dof_P->Case.Unknown.NonLocal ? -1 : 1) ;
416 break ;
417 case DOF_FIXEDWITHASSOCIATE :
418 fprintf(File_PRE, "%d ", Dof_P->Case.FixedAssociate.NumDof) ;
419 LinAlg_PrintScalar(File_PRE, &Dof_P->Val);
420 fprintf(File_PRE, " %d\n", Dof_P->Case.FixedAssociate.TimeFunctionIndex) ;
421 break ;
422 case DOF_FIXED :
423 LinAlg_PrintScalar(File_PRE, &Dof_P->Val);
424 fprintf(File_PRE, " %d\n", Dof_P->Case.FixedAssociate.TimeFunctionIndex) ;
425 break ;
426 case DOF_FIXED_SOLVE :
427 fprintf(File_PRE, "%d\n", Dof_P->Case.FixedAssociate.TimeFunctionIndex) ;
428 break ;
429 case DOF_UNKNOWN_INIT :
430 fprintf(File_PRE, "%d ", Dof_P->Case.Unknown.NumDof) ;
431 LinAlg_PrintScalar(File_PRE, &Dof_P->Val);
432 fprintf(File_PRE, " ") ;
433 LinAlg_PrintScalar(File_PRE, &Dof_P->Val2);
434 fprintf(File_PRE, " %d\n", Dof_P->Case.Unknown.NonLocal ? -1 : 1) ;
435 break ;
436 case DOF_LINK :
437 fprintf(File_PRE, "%.16g %d\n",
438 Dof_P->Case.Link.Coef, Dof_P->Case.Link.EntityRef) ;
439 break ;
440 case DOF_LINKCPLX :
441 fprintf(File_PRE, "%.16g %.16g %d\n",
442 Dof_P->Case.Link.Coef, Dof_P->Case.Link.Coef2,
443 Dof_P->Case.Link.EntityRef) ;
444 break ;
445 }
446 }
447
448 /* ------------------------------------------------------------------------ */
449 /* D o f _ R e a d F i l e P R E */
450 /* ------------------------------------------------------------------------ */
451
Dof_ReadFilePRE(struct DofData * DofData_P)452 void Dof_ReadFilePRE(struct DofData * DofData_P)
453 {
454 Message::Barrier();
455
456 int i, Nbr_Index, Int, Dummy ;
457 struct Dof Dof ;
458 char String[256] ;
459
460 do {
461 if(!fgets(String, sizeof(String), File_PRE)) break;
462 if(feof(File_PRE)) break;
463 } while (String[0] != '$') ;
464
465 if(feof(File_PRE)){
466 Message::Error("$DofData field not found in file");
467 return;
468 }
469
470 if(!strncmp(&String[1], "DofData", 7)) {
471
472 if(fscanf(File_PRE, "%d %d",
473 &DofData_P->ResolutionIndex, &DofData_P->SystemIndex) != 2){
474 Message::Error("Could not read Resolution and DofData indices");
475 return;
476 }
477
478 if(fscanf(File_PRE, "%d", &Nbr_Index) != 1){
479 Message::Error("Could not read number of function spaces");
480 return;
481 }
482 DofData_P->FunctionSpaceIndex = List_Create(Nbr_Index, 1, sizeof(int)) ;
483 for(i = 0 ; i < Nbr_Index ; i++) {
484 if(fscanf(File_PRE, "%d", &Int) != 1){
485 Message::Error("Could not read FunctionSpace index");
486 return;
487 }
488 List_Add(DofData_P->FunctionSpaceIndex, &Int) ;
489 }
490
491 if(fscanf(File_PRE, "%d", &Nbr_Index) != 1){
492 Message::Error("Could not read number of time functions");
493 return;
494 }
495 DofData_P->TimeFunctionIndex = List_Create(Nbr_Index, 1, sizeof(int)) ;
496 for(i = 0 ; i < Nbr_Index ; i++) {
497 if(fscanf(File_PRE, "%d", &Int) != 1){
498 Message::Error("Could not read TimeFunction index");
499 return;
500 }
501 List_Add(DofData_P->TimeFunctionIndex, &Int) ;
502 }
503
504 if(fscanf(File_PRE, "%d", &Dummy) != 1){
505 Message::Error("Format error");
506 return;
507 }
508 for(i = 0 ; i < 1 ; i++){
509 if(fscanf(File_PRE, "%d", &Dummy) != 1){
510 Message::Error("Format error");
511 return;
512 }
513 }
514
515 if(fscanf(File_PRE, "%d %d", &DofData_P->NbrAnyDof, &DofData_P->NbrDof) != 2){
516 Message::Error("Could not read number of dofs");
517 return;
518 }
519
520 DofData_P->DofList = List_Create(DofData_P->NbrAnyDof, 1, sizeof(struct Dof)) ;
521 Tree_Delete(DofData_P->DofTree);
522 DofData_P->DofTree = NULL;
523
524 for(i = 0 ; i < DofData_P->NbrAnyDof ; i++) {
525
526 if(fscanf(File_PRE, "%d %d %d %d",
527 &Dof.NumType, &Dof.Entity, &Dof.Harmonic, &Dof.Type) != 4){
528 Message::Error("Could not dof");
529 return;
530 }
531
532 switch (Dof.Type) {
533 case DOF_UNKNOWN :
534 if(fscanf(File_PRE, "%d %d", &Dof.Case.Unknown.NumDof, &Dummy) != 2){
535 Message::Error("Could not read DOF_UNKNOWN");
536 return;
537 }
538 Dof.Case.Unknown.NonLocal = (Dummy < 0) ? true : false;
539 if(Dummy < 0)
540 DofData_P->NonLocalEquations.push_back(Dof.Case.Unknown.NumDof);
541 break ;
542 case DOF_FIXEDWITHASSOCIATE :
543 if(fscanf(File_PRE, "%d", &Dof.Case.FixedAssociate.NumDof) != 1){
544 Message::Error("Could not read DOF_FIXEDWITHASSOCIATE");
545 return;
546 }
547 LinAlg_ScanScalar(File_PRE, &Dof.Val) ;
548 if(fscanf(File_PRE, "%d", &Dof.Case.FixedAssociate.TimeFunctionIndex) != 1){
549 Message::Error("Could not read DOF_FIXEDWITHASSOCIATE");
550 return;
551 }
552 break ;
553 case DOF_FIXED :
554 LinAlg_ScanScalar(File_PRE, &Dof.Val) ;
555 if(fscanf(File_PRE, "%d", &Dof.Case.FixedAssociate.TimeFunctionIndex) != 1){
556 Message::Error("Could not read DOF_FIXED");
557 return;
558 }
559 break ;
560 case DOF_FIXED_SOLVE :
561 if(fscanf(File_PRE, "%d", &Dof.Case.FixedAssociate.TimeFunctionIndex) != 1){
562 Message::Error("Could not read DOF_FIXED_SOLVED");
563 return;
564 }
565 break ;
566 case DOF_UNKNOWN_INIT :
567 if(fscanf(File_PRE, "%d", &Dof.Case.Unknown.NumDof) != 1){
568 Message::Error("Could not read DOF_UNKNOWN_INIT");
569 return;
570 }
571 LinAlg_ScanScalar(File_PRE, &Dof.Val) ;
572 LinAlg_ScanScalar(File_PRE, &Dof.Val2) ;
573 if(fscanf(File_PRE, "%d", &Dummy) != 1){
574 Message::Error("Could not read DOF_UNKNOWN_INIT");
575 return;
576 }
577 Dof.Case.Unknown.NonLocal = (Dummy < 0) ? true : false;
578 if(Dummy < 0)
579 DofData_P->NonLocalEquations.push_back(Dof.Case.Unknown.NumDof);
580 break ;
581 case DOF_LINK :
582 if(fscanf(File_PRE, "%lf %d",
583 &Dof.Case.Link.Coef, &Dof.Case.Link.EntityRef) != 2){
584 Message::Error("Could not read DOF_LINK");
585 return;
586 }
587 Dof.Case.Link.Dof = NULL ;
588 break ;
589 case DOF_LINKCPLX :
590 if(fscanf(File_PRE, "%lf %lf %d",
591 &Dof.Case.Link.Coef, &Dof.Case.Link.Coef2,
592 &Dof.Case.Link.EntityRef) != 3){
593 Message::Error("Could not read DOF_LINKCPLX");
594 return;
595 }
596 Dof.Case.Link.Dof = NULL ;
597 break ;
598 }
599
600 List_Add(DofData_P->DofList, &Dof) ;
601 }
602 }
603
604 do {
605 if(!fgets(String, sizeof(String), File_PRE) || feof(File_PRE)){
606 Message::Error("Premature end of file");
607 break;
608 }
609 } while (String[0] != '$') ;
610
611 Dof_InitDofType(DofData_P) ;
612 }
613
614 /* ------------------------------------------------------------------------ */
615 /* D o f _ W r i t e F i l e R E S 0 */
616 /* ------------------------------------------------------------------------ */
617
Dof_WriteFileRES0(char * Name_File,int Format)618 void Dof_WriteFileRES0(char * Name_File, int Format)
619 {
620 if(Message::GetIsCommWorld() && Message::GetCommRank()) return;
621
622 Dof_OpenFile(DOF_RES, Name_File, (char*)(Format ? "wb" : "w")) ;
623 fprintf(File_RES, "$ResFormat /* GetDP %s, %s */\n", GETDP_VERSION,
624 Format ? "binary" : "ascii") ;
625 fprintf(File_RES, "1.1 %d\n", Format) ;
626 fprintf(File_RES, "$EndResFormat\n") ;
627 Dof_CloseFile(DOF_RES) ;
628 }
629
630 /* ------------------------------------------------------------------------ */
631 /* D o f _ W r i t e F i l e R E S _ E x t e n d M H */
632 /* ------------------------------------------------------------------------ */
633
Dof_WriteFileRES_ExtendMH(char * Name_File,struct DofData * DofData_P,int Format,int NbrH)634 void Dof_WriteFileRES_ExtendMH(char * Name_File, struct DofData * DofData_P,
635 int Format, int NbrH)
636 {
637 if(Message::GetIsCommWorld() && Message::GetCommRank()) return;
638
639 if(!DofData_P->CurrentSolution){
640 Message::Warning("No ExtendMH solution to save");
641 return;
642 }
643
644 gVector x;
645 double d;
646 int i, inew;
647
648 Dof_OpenFile(DOF_RES, Name_File, (char*)(Format ? "ab" : "a")) ;
649
650 fprintf(File_RES, "$Solution /* DofData #%d */\n", DofData_P->Num) ;
651 fprintf(File_RES, "%d 0 0 0 \n", DofData_P->Num) ;
652
653 LinAlg_CreateVector(&x, &DofData_P->Solver, (DofData_P->NbrDof / Current.NbrHar) * NbrH) ;
654
655 LinAlg_ZeroVector(&x) ;
656
657 for(i=0 ; i<DofData_P->NbrDof ; i++){
658 LinAlg_GetDoubleInVector(&d, &DofData_P->CurrentSolution->x, i);
659 inew = (i / Current.NbrHar) * NbrH + i % Current.NbrHar;
660 LinAlg_SetDoubleInVector(d, &x, inew);
661 }
662
663 Format ?
664 LinAlg_WriteVector(File_RES,&x) :
665 LinAlg_PrintVector(File_RES,&x) ;
666
667 fprintf(File_RES, "$EndSolution\n") ;
668
669 Dof_CloseFile(DOF_RES) ;
670
671 LinAlg_DestroyVector(&x) ;
672 }
673
Dof_WriteFileRES_MHtoTime(char * Name_File,struct DofData * DofData_P,int Format,List_T * Time_L)674 void Dof_WriteFileRES_MHtoTime(char * Name_File, struct DofData * DofData_P,
675 int Format, List_T * Time_L)
676 {
677 if(Message::GetIsCommWorld() && Message::GetCommRank()) return;
678
679 if(!DofData_P->CurrentSolution){
680 Message::Warning("No solution to save");
681 return;
682 }
683
684 gVector x;
685 double Time, d1, d2, d, *Pulsation;
686 int iT, i, j, k;
687
688 Dof_OpenFile(DOF_RES, Name_File, (char*)(Format ? "ab" : "a")) ;
689
690
691 for(iT=0 ; iT<List_Nbr(Time_L) ; iT++){
692
693 List_Read(Time_L, iT, &Time);
694
695 fprintf(File_RES, "$Solution /* DofData #%d */\n", DofData_P->Num) ;
696 fprintf(File_RES, "%d %e 0 %d \n", DofData_P->Num,Time,iT) ;
697
698 Pulsation = DofData_P->Val_Pulsation ;
699
700 LinAlg_CreateVector(&x, &DofData_P->Solver, DofData_P->NbrDof/Current.NbrHar) ;
701
702 LinAlg_ZeroVector(&x) ;
703
704 for(i=0 ; i<DofData_P->NbrDof/Current.NbrHar ; i++) {
705 d = 0;
706 for(k=0 ; k<Current.NbrHar/2 ; k++) {
707 j = i * Current.NbrHar + 2*k ;
708 LinAlg_GetDoubleInVector(&d1, &DofData_P->CurrentSolution->x, j) ;
709 LinAlg_GetDoubleInVector(&d2, &DofData_P->CurrentSolution->x, j+1) ;
710 d += d1 * cos(Pulsation[k]*Time) - d2 * sin(Pulsation[k]*Time) ;
711 }
712 LinAlg_SetDoubleInVector(d, &x, i) ;
713 }
714
715 Format ?
716 LinAlg_WriteVector(File_RES,&x) :
717 LinAlg_PrintVector(File_RES,&x) ;
718
719 fprintf(File_RES, "$EndSolution\n") ;
720 }
721
722
723 Dof_CloseFile(DOF_RES) ;
724
725 LinAlg_DestroyVector(&x) ;
726 }
727
728 /* ------------------------------------------------------------------------ */
729 /* D o f _ W r i t e F i l e R E S */
730 /* ------------------------------------------------------------------------ */
731
Dof_WriteFileRES(char * Name_File,struct DofData * DofData_P,int Format,double Val_Time,double Val_TimeImag,int Val_TimeStep)732 void Dof_WriteFileRES(char * Name_File, struct DofData * DofData_P, int Format,
733 double Val_Time, double Val_TimeImag, int Val_TimeStep)
734 {
735 if(Message::GetIsCommWorld() && Message::GetCommRank()) return;
736
737 if(!DofData_P->CurrentSolution){
738 Message::Warning("No solution to save");
739 return;
740 }
741
742 Dof_OpenFile(DOF_RES, Name_File, (char*)(Format ? "ab" : "a")) ;
743
744 fprintf(File_RES, "$Solution /* DofData #%d */\n", DofData_P->Num) ;
745 fprintf(File_RES, "%d %.16g %.16g %d\n", DofData_P->Num, Val_Time,
746 Val_TimeImag, Val_TimeStep) ;
747
748 Format ?
749 LinAlg_WriteVector(File_RES, &DofData_P->CurrentSolution->x) :
750 LinAlg_PrintVector(File_RES, &DofData_P->CurrentSolution->x) ;
751
752 fprintf(File_RES, "$EndSolution\n") ;
753
754 Dof_CloseFile(DOF_RES) ;
755 }
756
757 /* ------------------------------------------------------------------------ */
758 /* D o f _ W r i t e F i l e R E S _ W i t h E n t i t y N u m */
759 /* ------------------------------------------------------------------------ */
760
Dof_WriteFileRES_WithEntityNum(char * Name_File,struct DofData * DofData_P,struct GeoData * GeoData_P0,struct Group * Group_P,bool saveFixed)761 void Dof_WriteFileRES_WithEntityNum(char * Name_File, struct DofData * DofData_P,
762 struct GeoData * GeoData_P0, struct Group *Group_P,
763 bool saveFixed)
764 {
765 if(Message::GetIsCommWorld() && Message::GetCommRank()) return;
766
767 if(!DofData_P->CurrentSolution){
768 Message::Warning("No solution to save");
769 return;
770 }
771
772 std::map<int, std::map<int, std::complex<double> > > unknowns;
773
774 List_T *l = !DofData_P->DofList ? Tree2List(DofData_P->DofTree) : 0;
775 int N = l ? List_Nbr(l) : List_Nbr(DofData_P->DofList);
776 for(int i = 0; i < N; i++){
777 struct Dof *dof;
778 if(l)
779 dof = (Dof*)List_Pointer(l, i);
780 else
781 dof = (Dof*)List_Pointer(DofData_P->DofList, i);
782 if(dof->Type == DOF_UNKNOWN || dof->Type == DOF_UNKNOWN_INIT){
783 gScalar s;
784 LinAlg_GetScalarInVector(&s, &DofData_P->CurrentSolution->x,
785 dof->Case.Unknown.NumDof - 1);
786 unknowns[dof->NumType][dof->Entity] = s.s;
787 }
788 if(saveFixed && dof->Type == DOF_FIXED){
789 unknowns[dof->NumType][dof->Entity] = dof->Val.s;
790 }
791 }
792
793 for(std::map<int, std::map<int, std::complex<double> > >::iterator it =
794 unknowns.begin(); it != unknowns.end(); it++){
795
796 // create files that can be interpreted by ListFromFile and Value/VectorFromIndex
797 char FileRe[256], FileIm[256] ;
798 if(unknowns.size() > 1){
799 sprintf(FileRe, "%s_%d", Name_File, it->first);
800 sprintf(FileIm, "%s_%d", Name_File, it->first);
801 }
802 else{
803 strcpy(FileRe, Name_File);
804 strcpy(FileIm, Name_File);
805 }
806 if(Current.NbrHar > 1){
807 strcat(FileRe, "_Re.txt");
808 strcat(FileIm, "_Im.txt");
809 }
810 else{
811 strcat(FileRe, ".txt");
812 strcat(FileIm, ".txt");
813 }
814 FILE *fpRe = FOpen(FileRe, "w");
815 if(!fpRe){
816 Message::Error("Unable to open file '%s'", FileRe) ;
817 return;
818 }
819 FILE *fpIm = 0;
820 if(Current.NbrHar > 1){
821 fpIm = FOpen(FileIm, "w");
822 if(!fpIm){
823 Message::Error("Unable to open file '%s'", FileIm) ;
824 return;
825 }
826 }
827
828 // create vectors that can be shared as lists
829 std::vector<double> exportRe, exportIm;
830
831 if(!Group_P){
832 int n = (int)it->second.size();
833 fprintf(fpRe, "%d\n", n);
834 exportRe.push_back(n);
835 if(fpIm){
836 fprintf(fpIm, "%d\n", n);
837 exportIm.push_back(n);
838 }
839 for(std::map<int, std::complex<double> >::iterator it2 = it->second.begin();
840 it2 != it->second.end(); it2++){
841 fprintf(fpRe, "%d %.16g\n", it2->first, it2->second.real());
842 exportRe.push_back(it2->first);
843 exportRe.push_back(it2->second.real());
844 if(fpIm){
845 fprintf(fpIm, "%d %.16g\n", it2->first, it2->second.imag());
846 exportIm.push_back(it2->first);
847 exportIm.push_back(it2->second.imag());
848 }
849 }
850 }
851 else{
852 Message::Info("Writing solution for all entities in group '%s'", Group_P->Name) ;
853
854 // force generation of extended list (necessary when using multiple
855 // meshes)
856 List_Delete(Group_P->ExtendedList);
857 Generate_ExtendedGroup(Group_P) ;
858 int n = List_Nbr(Group_P->ExtendedList);
859 fprintf(fpRe, "%d\n", n);
860 exportRe.push_back(n);
861 if(fpIm){
862 fprintf(fpIm, "%d\n", n);
863 exportIm.push_back(n);
864 }
865 for(int i = 0; i < List_Nbr(Group_P->ExtendedList); i++){
866 int num;
867 List_Read(Group_P->ExtendedList, i, &num);
868 if(!Group_P->InitialSuppList ||
869 (!List_Search(Group_P->ExtendedSuppList, &num, fcmp_int))){
870 // SuppList assumed to be "Not"!
871 if(it->second.count(num)){
872 std::complex<double> s = it->second[num];
873 fprintf(fpRe, "%d %.16g\n", num, s.real());
874 exportRe.push_back(num);
875 exportRe.push_back(s.real());
876 if(fpIm){
877 fprintf(fpIm, "%d %.16g\n", num, s.imag());
878 exportIm.push_back(num);
879 exportIm.push_back(s.imag());
880 }
881 }
882 else{
883 // yes, write zero: that's on purpose for the iterative schemes
884 fprintf(fpRe, "%d 0\n", num);
885 exportRe.push_back(num);
886 exportRe.push_back(0);
887 if(fpIm){
888 fprintf(fpIm, "%d 0\n", num);
889 exportIm.push_back(num);
890 exportIm.push_back(0);
891 }
892 }
893 }
894 }
895 }
896 fclose(fpRe);
897 GetDPNumbers[FileRe] = exportRe;
898 if(fpIm){
899 fclose(fpIm);
900 GetDPNumbers[FileIm] = exportIm;
901 }
902 }
903
904 List_Delete(l);
905
906 }
907
908 /* ------------------------------------------------------------------------ */
909 /* D o f _ R e a d F i l e R E S */
910 /* ------------------------------------------------------------------------ */
911
Dof_ReadFileRES(List_T * DofData_L,struct DofData * Read_DofData_P,int Read_DofData,double * Time,double * TimeImag,double * TimeStep)912 void Dof_ReadFileRES(List_T * DofData_L, struct DofData * Read_DofData_P,
913 int Read_DofData, double *Time, double *TimeImag,
914 double *TimeStep)
915 {
916 Message::Barrier();
917
918 int Num_DofData, Val_TimeStep, Format = 0, Read ;
919 double Val_Time, Val_TimeImag = 0., Version = 0.;
920 struct DofData * DofData_P = NULL ;
921 struct Solution Solution_S ;
922 char String[256] ;
923
924 while (1) {
925
926 do {
927 if(!fgets(String, sizeof(String), File_RES)) break;
928 if(feof(File_RES)) break ;
929 } while (String[0] != '$') ;
930
931 if(feof(File_RES)) break ;
932
933 /* F o r m a t */
934
935 if(!strncmp(&String[1], "ResFormat", 9)) {
936 if(fscanf(File_RES, "%lf %d\n", &Version, &Format) != 2){
937 Message::Error("Could not read ResFormat");
938 return;
939 }
940 }
941
942 /* S o l u t i o n */
943
944 if(!strncmp(&String[1], "Solution", 8)) {
945
946 /* don't use fscanf directly on the stream here: the data that
947 follows can be binary, and the first character could be
948 e.g. 0x0d, which would cause fscanf to eat-up one character
949 too much, leading to an offset in fread */
950 if(!fgets(String, sizeof(String), File_RES)){
951 Message::Error("Could not read Solution");
952 return;
953 }
954 if(Version <= 1.0)
955 sscanf(String, "%d %lf %d", &Num_DofData, &Val_Time, &Val_TimeStep) ;
956 else
957 sscanf(String, "%d %lf %lf %d", &Num_DofData, &Val_Time, &Val_TimeImag,
958 &Val_TimeStep) ;
959
960 if(Read_DofData < 0){
961 Read = 1 ; DofData_P = (struct DofData*)List_Pointer(DofData_L, Num_DofData) ;
962 }
963 else if(Num_DofData == Read_DofData) {
964 Read = 1 ; DofData_P = Read_DofData_P ;
965 }
966 else {
967 Read = 0 ;
968 }
969
970 if(Read){
971 Solution_S.Time = Val_Time ;
972 Solution_S.TimeImag = Val_TimeImag ;
973 Solution_S.TimeStep = Val_TimeStep ;
974 Solution_S.SolutionExist = 1 ;
975 Solution_S.TimeFunctionValues = NULL ;
976
977 LinAlg_CreateVector(&Solution_S.x, &DofData_P->Solver, DofData_P->NbrDof) ;
978 Format ?
979 LinAlg_ReadVector(File_RES, &Solution_S.x) :
980 LinAlg_ScanVector(File_RES, &Solution_S.x) ;
981
982 if(DofData_P->Solutions == NULL)
983 DofData_P->Solutions = List_Create( 20, 20, sizeof(struct Solution)) ;
984 List_Add(DofData_P->Solutions, &Solution_S) ;
985 }
986 }
987
988 do {
989 if(!fgets(String, sizeof(String), File_RES) || feof(File_RES)){
990 Message::Warning("Premature end of file (Time Step %d)", Val_TimeStep);
991 break;
992 }
993 } while (String[0] != '$') ;
994
995 } /* while 1 ... */
996
997 *Time = Val_Time ;
998 *TimeImag = Val_TimeImag ;
999 *TimeStep = (double)Val_TimeStep ;
1000 }
1001
1002 /* ------------------------------------------------------------------------ */
1003 /* D o f _ T r a n s f e r D o f T r e e T o L i s t */
1004 /* ------------------------------------------------------------------------ */
1005
Dof_TransferDofTreeToList(struct DofData * DofData_P)1006 void Dof_TransferDofTreeToList(struct DofData * DofData_P)
1007 {
1008 if(DofData_P->DofTree) {
1009 DofData_P->DofList = Tree2List(DofData_P->DofTree) ;
1010 Tree_Delete(DofData_P->DofTree) ;
1011 DofData_P->DofTree = NULL ;
1012 DofData_P->NbrAnyDof = List_Nbr(DofData_P->DofList) ;
1013 }
1014
1015 Dof_InitDofType(DofData_P) ;
1016 }
1017
1018 /* ------------------------------------------------------------------------ */
1019 /* D o f _ I n i t D o f T y p e */
1020 /* ------------------------------------------------------------------------ */
1021
Dof_InitDofType(struct DofData * DofData_P)1022 void Dof_InitDofType(struct DofData * DofData_P)
1023 {
1024 struct Dof * Dof_P, * Dof_P0 ;
1025 int i ;
1026
1027 if(!DofData_P->NbrAnyDof){
1028 return;
1029 }
1030
1031 Dof_P0 = (struct Dof *)List_Pointer(DofData_P->DofList, 0) ;
1032
1033 for(i = 0 ; i < DofData_P->NbrAnyDof ; i++) {
1034 Dof_P = Dof_P0 + i ;
1035
1036 switch (Dof_P->Type) {
1037 case DOF_LINK :
1038 case DOF_LINKCPLX :
1039 Dof_P->Case.Link.Dof =
1040 Dof_GetDofStruct(DofData_P, Dof_P->NumType,
1041 Dof_P->Case.Link.EntityRef, Dof_P->Harmonic) ;
1042 if(Dof_P->Case.Link.Dof == NULL
1043 ||
1044 Dof_P->Case.Link.Dof ==
1045 Dof_GetDofStruct(DofData_P, Dof_P->NumType,
1046 Dof_P->Entity, Dof_P->Harmonic)
1047 ) {
1048 Dof_P->Case.Link.Dof = /* Attention: bricolage ... */
1049 Dof_GetDofStruct(DofData_P, Dof_P->NumType-1,
1050 Dof_P->Case.Link.EntityRef, Dof_P->Harmonic) ;
1051 if(Dof_P->Case.Link.Dof == NULL)
1052 Message::Error("Wrong Link Constraint: reference Dof (%d %d %d) does not exist",
1053 Dof_P->NumType, Dof_P->Case.Link.EntityRef, Dof_P->Harmonic);
1054 }
1055 /*
1056 if(Dof_P->Case.Link.Dof == NULL)
1057 Message::Error("Wrong Link Constraint: reference Dof (%d %d %d) does not exist",
1058 Dof_P->NumType, Dof_P->Case.Link.EntityRef, Dof_P->Harmonic);
1059 */
1060 break ;
1061 default :
1062 break ;
1063 }
1064
1065 }
1066 }
1067
1068 /* ------------------------------------------------------------------------ */
1069 /* P R E P R O C E S S I N G ( C o d e s i n T r e e ) */
1070 /* ------------------------------------------------------------------------ */
1071
1072 /* ------------------------------------------------------------------------ */
1073 /* D o f _ A d d F u n c t i o n S p a c e I n d e x */
1074 /* ------------------------------------------------------------------------ */
1075
Dof_AddFunctionSpaceIndex(int Index_FunctionSpace)1076 void Dof_AddFunctionSpaceIndex(int Index_FunctionSpace)
1077 {
1078 if(CurrentDofData->FunctionSpaceIndex == NULL)
1079 CurrentDofData->FunctionSpaceIndex = List_Create(10, 5, sizeof(int)) ;
1080
1081 if(List_PQuery
1082 (CurrentDofData->FunctionSpaceIndex, &Index_FunctionSpace, fcmp_int) == NULL) {
1083 List_Add(CurrentDofData->FunctionSpaceIndex, &Index_FunctionSpace) ;
1084 List_Sort(CurrentDofData->FunctionSpaceIndex, fcmp_int) ;
1085 }
1086 }
1087
1088 /* ------------------------------------------------------------------------ */
1089 /* D o f _ A d d T i m e F u n c t i o n I n d e x */
1090 /* ------------------------------------------------------------------------ */
1091
Dof_AddTimeFunctionIndex(int Index_TimeFunction)1092 void Dof_AddTimeFunctionIndex(int Index_TimeFunction)
1093 {
1094 if(List_PQuery
1095 (CurrentDofData->TimeFunctionIndex, &Index_TimeFunction, fcmp_int) == NULL) {
1096 List_Add(CurrentDofData->TimeFunctionIndex, &Index_TimeFunction) ;
1097 List_Sort(CurrentDofData->TimeFunctionIndex, fcmp_int) ;
1098 }
1099 }
1100
1101 /* ------------------------------------------------------------------------ */
1102 /* D o f _ A d d P u l s a t i o n */
1103 /* ------------------------------------------------------------------------ */
1104
Dof_AddPulsation(struct DofData * DofData_P,double Val_Pulsation)1105 void Dof_AddPulsation(struct DofData * DofData_P, double Val_Pulsation)
1106 {
1107 if(DofData_P->Pulsation == NULL)
1108 DofData_P->Pulsation = List_Create(1, 2, sizeof(double)) ;
1109 List_Add(DofData_P->Pulsation, &Val_Pulsation) ;
1110 /*
1111 if(List_PQuery
1112 (DofData_P->Pulsation, &Val_Pulsation, fcmp_double) == NULL) {
1113 List_Add(DofData_P->Pulsation, &Val_Pulsation) ;
1114 List_Sort(DofData_P->Pulsation, fcmp_double) ;
1115 }
1116 */
1117 }
1118
1119 /* ------------------------------------------------------------------------ */
1120 /* D o f _ D e f i n e A s s i g n F i x e d D o f */
1121 /* ------------------------------------------------------------------------ */
1122
Dof_DefineAssignFixedDof(int D1,int D2,int NbrHar,double * Val,int Index_TimeFunction)1123 void Dof_DefineAssignFixedDof(int D1, int D2, int NbrHar, double *Val,
1124 int Index_TimeFunction)
1125 {
1126 struct Dof Dof, * Dof_P ;
1127 int k ;
1128
1129 Dof.NumType = D1 ; Dof.Entity = D2 ;
1130
1131 for(k=0 ; k<NbrHar ; k+=gSCALAR_SIZE){
1132 Dof.Harmonic = k ;
1133 if(!(Dof_P = (struct Dof *)Tree_PQuery(CurrentDofData->DofTree, &Dof))) {
1134 Dof.Type = DOF_FIXED ;
1135 LinAlg_SetScalar(&Dof.Val, &Val[k]) ;
1136 Dof.Case.FixedAssociate.TimeFunctionIndex = Index_TimeFunction + 1 ;
1137 Dof_AddTimeFunctionIndex(Index_TimeFunction + 1) ;
1138 Tree_Add(CurrentDofData->DofTree, &Dof) ;
1139 }
1140 else if(Dof_P->Type == DOF_UNKNOWN) {
1141 if(Message::GetVerbosity() == 10)
1142 Message::Info("Overriding unknown Dof with fixed Dof");
1143 Dof_P->Type = DOF_FIXED ;
1144 LinAlg_SetScalar(&Dof_P->Val, &Val[k]) ;
1145 Dof_P->Case.FixedAssociate.TimeFunctionIndex = Index_TimeFunction + 1 ;
1146 Dof_AddTimeFunctionIndex(Index_TimeFunction + 1) ;
1147 }
1148 }
1149 }
1150
1151 /* ------------------------------------------------------------------------ */
1152 /* D o f _ D e f i n e A s s i g n S o l v e D o f */
1153 /* ------------------------------------------------------------------------ */
1154
Dof_DefineAssignSolveDof(int D1,int D2,int NbrHar,int Index_TimeFunction)1155 void Dof_DefineAssignSolveDof(int D1, int D2, int NbrHar, int Index_TimeFunction)
1156 {
1157 struct Dof Dof ;
1158 int k ;
1159
1160 Dof.NumType = D1 ; Dof.Entity = D2 ;
1161
1162 for(k=0 ; k<NbrHar ; k+=gSCALAR_SIZE){
1163 Dof.Harmonic = k ;
1164 if(!Tree_PQuery(CurrentDofData->DofTree, &Dof)) {
1165 Dof.Type = DOF_FIXED_SOLVE ;
1166 Dof.Case.FixedAssociate.TimeFunctionIndex = Index_TimeFunction + 1 ;
1167 Dof_AddTimeFunctionIndex(Index_TimeFunction + 1) ;
1168 Tree_Add(CurrentDofData->DofTree, &Dof) ;
1169 }
1170 }
1171 }
1172
1173 /* ------------------------------------------------------------------------ */
1174 /* D o f _ D e f i n e I n i t F i x e d D o f */
1175 /* ------------------------------------------------------------------------ */
1176
Dof_DefineInitFixedDof(int D1,int D2,int NbrHar,double * Val,double * Val2,bool NonLocal)1177 void Dof_DefineInitFixedDof(int D1, int D2, int NbrHar, double *Val,
1178 double *Val2, bool NonLocal)
1179 {
1180 struct Dof Dof ;
1181 int k ;
1182
1183 Dof.NumType = D1 ; Dof.Entity = D2 ;
1184
1185 for(k=0 ; k<NbrHar ; k+=gSCALAR_SIZE){
1186 Dof.Harmonic = k ;
1187 if(!Tree_PQuery(CurrentDofData->DofTree, &Dof)) {
1188 Dof.Type = DOF_UNKNOWN_INIT ;
1189 LinAlg_SetScalar(&Dof.Val, &Val[k]) ;
1190 LinAlg_SetScalar(&Dof.Val2, &Val2[k]) ;
1191 Dof.Case.Unknown.NumDof = ++(CurrentDofData->NbrDof) ;
1192 Dof.Case.Unknown.NonLocal = NonLocal;
1193 Tree_Add(CurrentDofData->DofTree, &Dof) ;
1194 }
1195 }
1196 }
1197
1198 /* ------------------------------------------------------------------------ */
1199 /* D o f _ D e f i n e I n i t S o l v e D o f */
1200 /* ------------------------------------------------------------------------ */
1201
Dof_DefineInitSolveDof(int D1,int D2,int NbrHar)1202 void Dof_DefineInitSolveDof(int D1, int D2, int NbrHar)
1203 {
1204 struct Dof Dof ;
1205 int k ;
1206
1207 Dof.NumType = D1 ; Dof.Entity = D2 ;
1208
1209 for(k=0 ; k<NbrHar ; k+=gSCALAR_SIZE){
1210 Dof.Harmonic = k ;
1211 if(!Tree_PQuery(CurrentDofData->DofTree, &Dof)) {
1212 Dof.Type = DOF_UNKNOWN_INIT ;
1213 Dof.Case.Unknown.NumDof = ++(CurrentDofData->NbrDof) ;
1214 Dof.Case.Unknown.NonLocal = false ;
1215 Tree_Add(CurrentDofData->DofTree, &Dof) ;
1216 }
1217 }
1218 }
1219
1220 /* ------------------------------------------------------------------------ */
1221 /* D o f _ D e f i n e L i n k D o f */
1222 /* ------------------------------------------------------------------------ */
1223
Dof_DefineLinkDof(int D1,int D2,int NbrHar,double Value[],int D2_Link)1224 void Dof_DefineLinkDof(int D1, int D2, int NbrHar, double Value[], int D2_Link)
1225 {
1226 struct Dof Dof ;
1227 int k ;
1228
1229 Dof.NumType = D1 ; Dof.Entity = D2 ;
1230
1231 for(k=0 ; k<NbrHar ; k+=gSCALAR_SIZE){
1232 Dof.Harmonic = k ;
1233 if(!Tree_PQuery(CurrentDofData->DofTree, &Dof)) {
1234 Dof.Type = DOF_LINK ;
1235 Dof.Case.Link.Coef = Value[0] ;
1236 Dof.Case.Link.EntityRef = D2_Link ;
1237 Dof.Case.Link.Dof = NULL ;
1238 Tree_Add(CurrentDofData->DofTree, &Dof) ;
1239 }
1240 }
1241 }
1242
Dof_DefineLinkCplxDof(int D1,int D2,int NbrHar,double Value[],int D2_Link)1243 void Dof_DefineLinkCplxDof(int D1, int D2, int NbrHar, double Value[], int D2_Link)
1244 {
1245 struct Dof Dof ;
1246 int k ;
1247
1248 Dof.NumType = D1 ; Dof.Entity = D2 ;
1249
1250 for(k=0 ; k<NbrHar ; k+=gSCALAR_SIZE){
1251 Dof.Harmonic = k ;
1252 if(!Tree_PQuery(CurrentDofData->DofTree, &Dof)) {
1253 Dof.Type = DOF_LINKCPLX ;
1254 Dof.Case.Link.Coef = Value[0] ;
1255 Dof.Case.Link.Coef2 = Value[1] ;
1256 Dof.Case.Link.EntityRef = D2_Link ;
1257 Dof.Case.Link.Dof = NULL ;
1258 Tree_Add(CurrentDofData->DofTree, &Dof) ;
1259 }
1260 }
1261 }
1262
1263 /* ------------------------------------------------------------------------ */
1264 /* D o f _ D e f i n e U n k n o w n D o f */
1265 /* ------------------------------------------------------------------------ */
1266
Dof_DefineUnknownDof(int D1,int D2,int NbrHar,bool NonLocal)1267 void Dof_DefineUnknownDof(int D1, int D2, int NbrHar, bool NonLocal)
1268 {
1269 struct Dof Dof ;
1270 int k ;
1271
1272 Dof.NumType = D1 ; Dof.Entity = D2 ;
1273
1274 for(k=0 ; k<NbrHar ; k+=gSCALAR_SIZE){
1275 Dof.Harmonic = k ;
1276 if(!Tree_PQuery(CurrentDofData->DofTree, &Dof)) {
1277 Dof.Type = DOF_UNKNOWN ;
1278 /* Dof.Case.Unknown.NumDof = ++(CurrentDofData->NbrDof) ; */
1279 Dof.Case.Unknown.NumDof = -1 ;
1280 Dof.Case.Unknown.NonLocal = NonLocal ;
1281 Tree_Add(CurrentDofData->DofTree, &Dof) ;
1282 }
1283 }
1284 }
1285
NumberUnknownDof(void * a,void * b)1286 static void NumberUnknownDof (void *a, void *b)
1287 {
1288 struct Dof * Dof_P ;
1289
1290 Dof_P = (struct Dof *)a ;
1291
1292 if(Dof_P->Type == DOF_UNKNOWN){
1293 if(Dof_P->Case.Unknown.NumDof == -1)
1294 Dof_P->Case.Unknown.NumDof = ++(CurrentDofData->NbrDof) ;
1295 if(Dof_P->Case.Unknown.NonLocal)
1296 CurrentDofData->NonLocalEquations.push_back(Dof_P->Case.Unknown.NumDof);
1297 }
1298 }
1299
Dof_NumberUnknownDof(void)1300 void Dof_NumberUnknownDof(void)
1301 {
1302 if(CurrentDofData->DofTree)
1303 Tree_Action(CurrentDofData->DofTree, NumberUnknownDof) ;
1304 else
1305 List_Action(CurrentDofData->DofList, NumberUnknownDof) ;
1306 }
1307
1308 /* ------------------------------------------------------------------------ */
1309 /* D o f _ D e f i n e A s s o c i a t e D o f */
1310 /* ------------------------------------------------------------------------ */
1311
Dof_DefineAssociateDof(int E1,int E2,int D1,int D2,int NbrHar,int init,double * Val)1312 void Dof_DefineAssociateDof(int E1, int E2, int D1, int D2, int NbrHar,
1313 int init, double *Val)
1314 {
1315 struct Dof Dof, Equ, * Equ_P ;
1316 int k ;
1317
1318 Equ.NumType = E1 ; Equ.Entity = E2 ;
1319
1320 for(k=0 ; k<NbrHar ; k+=gSCALAR_SIZE){
1321 Equ.Harmonic = k ;
1322 if((Equ_P = (struct Dof*)Tree_PQuery(CurrentDofData->DofTree, &Equ))) {
1323 switch (Equ_P->Type) {
1324 case DOF_FIXED :
1325 Equ_P->Type = DOF_FIXEDWITHASSOCIATE ;
1326 Equ_P->Case.FixedAssociate.NumDof = ++(CurrentDofData->NbrDof) ;
1327 /* To be modified (Patrick): strange to define a new NumDof for Equ if
1328 associate-Dof already exists */
1329 Dof.NumType = D1 ; Dof.Entity = D2 ; Dof.Harmonic = k ;
1330 if(!Tree_PQuery(CurrentDofData->DofTree, &Dof)) {
1331 if(!init) {
1332 Dof.Type = DOF_UNKNOWN ;
1333 }
1334 else {
1335 Dof.Type = DOF_UNKNOWN_INIT ;
1336 LinAlg_SetScalar(&Dof.Val, &Val[k]) ;
1337 LinAlg_ZeroScalar(&Dof.Val2) ;
1338 }
1339 Dof.Case.Unknown.NumDof = CurrentDofData->NbrDof ;
1340 Dof.Case.Unknown.NonLocal = true ;
1341 Tree_Add(CurrentDofData->DofTree, &Dof) ;
1342 }
1343 break ;
1344 case DOF_FIXED_SOLVE :
1345 Equ_P->Type = DOF_FIXEDWITHASSOCIATE_SOLVE ;
1346 Equ_P->Case.FixedAssociate.NumDof = ++(CurrentDofData->NbrDof) ;
1347 Dof.NumType = D1 ; Dof.Entity = D2 ; Dof.Harmonic = k ;
1348 if(!Tree_PQuery(CurrentDofData->DofTree, &Dof)) {
1349 if(!init) {
1350 Dof.Type = DOF_UNKNOWN ;
1351 }
1352 else {
1353 Dof.Type = DOF_UNKNOWN_INIT ;
1354 LinAlg_SetScalar(&Dof.Val, &Val[k]) ;
1355 LinAlg_ZeroScalar(&Dof.Val2) ;
1356 }
1357 Dof.Case.Unknown.NumDof = CurrentDofData->NbrDof ;
1358 Dof.Case.Unknown.NonLocal = true ;
1359 Tree_Add(CurrentDofData->DofTree, &Dof) ;
1360 }
1361 break ;
1362 case DOF_UNKNOWN : case DOF_UNKNOWN_INIT :
1363 Dof_DefineUnknownDof(D1, D2, NbrHar) ;
1364 break ;
1365 }
1366 }
1367 }
1368 }
1369
1370 /* ------------------------------------------------------------------------ */
1371 /* P R O C E S S I N G ( C o d e s i n L i s t ) */
1372 /* ------------------------------------------------------------------------ */
1373
1374 /* ------------------------------------------------------------------------ */
1375 /* D o f _ G e t D o f S t r u c t */
1376 /* ------------------------------------------------------------------------ */
1377
Dof_GetDofStruct(struct DofData * DofData_P,int D1,int D2,int D3)1378 struct Dof *Dof_GetDofStruct(struct DofData * DofData_P, int D1, int D2, int D3)
1379 {
1380 struct Dof Dof ;
1381
1382 Dof.NumType = D1 ; Dof.Entity = D2 ; Dof.Harmonic = D3 ;
1383
1384 return (struct Dof *)List_PQuery(DofData_P->DofList, &Dof, fcmp_Dof);
1385 }
1386
1387 /* ------------------------------------------------------------------------ */
1388 /* D o f _ U p d a t e A s s i g n F i x e d D o f */
1389 /* ------------------------------------------------------------------------ */
1390
Dof_UpdateAssignFixedDof(int D1,int D2,int NbrHar,double * Val,double * Val2)1391 void Dof_UpdateAssignFixedDof(int D1, int D2, int NbrHar, double *Val, double *Val2)
1392 {
1393 struct Dof Dof, * Dof_P ;
1394 int k ;
1395
1396 Dof.NumType = D1 ; Dof.Entity = D2 ;
1397
1398 for(k=0 ; k<NbrHar ; k+=gSCALAR_SIZE){
1399 Dof.Harmonic = k ;
1400 if(CurrentDofData->DofTree)
1401 Dof_P = (struct Dof *)Tree_PQuery(CurrentDofData->DofTree, &Dof);
1402 else
1403 Dof_P = (struct Dof *)List_PQuery(CurrentDofData->DofList, &Dof, fcmp_Dof);
1404 LinAlg_SetScalar(&Dof_P->Val, &Val[Dof_P->Harmonic]) ;
1405 LinAlg_SetScalar(&Dof_P->Val2, &Val2[Dof_P->Harmonic]) ;
1406 }
1407 }
1408
1409 /* ------------------------------------------------------------------------ */
1410 /* D o f _ U p d a t e L i n k D o f */
1411 /* ------------------------------------------------------------------------ */
1412
Dof_UpdateLinkDof(int D1,int D2,int NbrHar,double Value[],int D2_Link)1413 void Dof_UpdateLinkDof(int D1, int D2, int NbrHar, double Value[], int D2_Link)
1414 {
1415 struct Dof Dof, * Dof_P ;
1416 int k ;
1417
1418 Dof.NumType = D1 ; Dof.Entity = D2 ;
1419
1420 for(k=0 ; k<NbrHar ; k+=gSCALAR_SIZE){
1421 Dof.Harmonic = k ;
1422 if(CurrentDofData->DofTree)
1423 Dof_P = (struct Dof *)Tree_PQuery(CurrentDofData->DofTree, &Dof);
1424 else
1425 Dof_P = (struct Dof *)List_PQuery(CurrentDofData->DofList, &Dof, fcmp_Dof);
1426
1427 if(Dof_P->Type == DOF_LINK || Dof_P->Type == DOF_LINKCPLX) {
1428 /*
1429 fprintf(stderr,"===> %d %d %.16g\n", Dof_P->NumType, Dof_P->Entity, Value[0]) ;
1430 */
1431 Dof_P->Case.Link.Coef = Value[0] ;
1432 if(Dof_P->Type == DOF_LINKCPLX)
1433 Dof_P->Case.Link.Coef2 = Value[1] ;
1434 Dof_P->Case.Link.EntityRef = D2_Link ;
1435 Dof_P->Case.Link.Dof = NULL ;
1436 }
1437 }
1438 }
1439
1440 /* ------------------------------------------------------------------------ */
1441 /* D o f _ A s s e m b l e I n M a t */
1442 /* ------------------------------------------------------------------------ */
1443
Dof_AssembleInMat(struct Dof * Equ_P,struct Dof * Dof_P,int NbrHar,double * Val,gMatrix * Mat,gVector * Vec,List_T * Vecs)1444 void Dof_AssembleInMat(struct Dof * Equ_P, struct Dof * Dof_P, int NbrHar,
1445 double * Val, gMatrix * Mat, gVector * Vec,
1446 List_T * Vecs)
1447 {
1448 gScalar tmp, tmp2 ;
1449 double valtmp[2], d1, d2 ;
1450
1451 switch (Equ_P->Type) {
1452 case DOF_UNKNOWN :
1453 case DOF_FIXEDWITHASSOCIATE :
1454
1455 switch (Dof_P->Type) {
1456
1457 case DOF_UNKNOWN :
1458 if(Current.DofData->Flag_RHS) break;
1459 if(NbrHar==1){
1460 LinAlg_AddDoubleInMatrix
1461 (Val[0], Mat,
1462 Equ_P->Case.Unknown.NumDof-1, Dof_P->Case.Unknown.NumDof-1) ;
1463 }
1464 else
1465 LinAlg_AddComplexInMatrix
1466 (Val[0], Val[1], Mat,
1467 Equ_P->Case.Unknown.NumDof-1, Dof_P->Case.Unknown.NumDof-1,
1468 (gSCALAR_SIZE==1)?((Equ_P+1)->Case.Unknown.NumDof-1):-1,
1469 (gSCALAR_SIZE==1)?((Dof_P+1)->Case.Unknown.NumDof-1):-1) ;
1470 break ;
1471
1472 case DOF_FIXED :
1473 case DOF_FIXEDWITHASSOCIATE :
1474 if(Vec){
1475 if(NbrHar==1){
1476 if(Val[0]){
1477 LinAlg_ProdScalarDouble
1478 (&Dof_P->Val,
1479 CurrentDofData->CurrentSolution->
1480 TimeFunctionValues[Dof_P->Case.FixedAssociate.TimeFunctionIndex],
1481 &tmp);
1482 LinAlg_ProdScalarDouble(&tmp, -Val[0], &tmp) ;
1483 LinAlg_AddScalarInVector(&tmp, Vec, Equ_P->Case.Unknown.NumDof-1) ;
1484 if(Vecs){ // experimental
1485 int index = List_ISearchSeq(Current.DofData->TimeFunctionIndex,
1486 &Dof_P->Case.FixedAssociate.TimeFunctionIndex,
1487 fcmp_int);
1488 if(index >= 0 && index < List_Nbr(Vecs)){
1489 gVector *v = (gVector*)List_Pointer(Vecs, index);
1490 LinAlg_AddScalarInVector(&tmp, v, Equ_P->Case.Unknown.NumDof-1) ;
1491 }
1492 else{
1493 Message::Error("Something wrong in multi-vec assembly");
1494 }
1495 }
1496 }
1497 }
1498 else{
1499 LinAlg_ProdScalarDouble
1500 (&Dof_P->Val,
1501 CurrentDofData->CurrentSolution->
1502 TimeFunctionValues[Dof_P->Case.FixedAssociate.TimeFunctionIndex],
1503 &tmp);
1504 if(gSCALAR_SIZE == 2){
1505 LinAlg_ProdScalarComplex(&tmp, -Val[0], -Val[1], &valtmp[0], &valtmp[1]) ;
1506 }
1507 else{
1508 LinAlg_GetDoubleInScalar(&d1, &tmp);
1509 LinAlg_ProdScalarDouble
1510 (&(Dof_P+1)->Val,
1511 CurrentDofData->CurrentSolution->
1512 TimeFunctionValues[Dof_P->Case.FixedAssociate.TimeFunctionIndex],
1513 &tmp2);
1514 LinAlg_GetDoubleInScalar(&d2, &tmp2);
1515 valtmp[0] = -d1*Val[0] + d2*Val[1] ;
1516 valtmp[1] = -d1*Val[1] - d2*Val[0] ;
1517 }
1518 LinAlg_AddComplexInVector
1519 (valtmp[0], valtmp[1], Vec,
1520 Equ_P->Case.Unknown.NumDof-1,
1521 (gSCALAR_SIZE==1)?((Equ_P+1)->Case.Unknown.NumDof-1):-1) ;
1522 }
1523 }
1524 break ;
1525
1526 case DOF_LINK :
1527 if(NbrHar==1)
1528 valtmp[0] = Val[0] * Dof_P->Case.Link.Coef ;
1529 else{
1530 valtmp[0] = Val[0] * Dof_P->Case.Link.Coef ;
1531 valtmp[1] = Val[1] * Dof_P->Case.Link.Coef ;
1532 }
1533 Dof_AssembleInMat(Equ_P, Dof_P->Case.Link.Dof, NbrHar, valtmp, Mat, Vec, Vecs) ;
1534 break ;
1535
1536 case DOF_LINKCPLX :
1537 if(NbrHar==1)
1538 Message::Error("LinkCplx only valid for Complex systems") ;
1539 else{
1540 valtmp[0] = Val[0] * Dof_P->Case.Link.Coef - Val[1] * Dof_P->Case.Link.Coef2 ;
1541 valtmp[1] = Val[1] * Dof_P->Case.Link.Coef + Val[0] * Dof_P->Case.Link.Coef2 ;
1542 }
1543 Dof_AssembleInMat(Equ_P, Dof_P->Case.Link.Dof, NbrHar, valtmp, Mat, Vec, Vecs) ;
1544 break ;
1545
1546 case DOF_FIXED_SOLVE : case DOF_FIXEDWITHASSOCIATE_SOLVE :
1547 Message::Error("Wrong Constraints: "
1548 "remaining Dof(s) waiting to be fixed by a Resolution");
1549 break;
1550
1551 case DOF_UNKNOWN_INIT :
1552 Message::Error("Wrong Initial Constraints: "
1553 "remaining Dof(s) with non-fixed initial conditions");
1554 break;
1555 }
1556
1557 break ;
1558
1559 case DOF_LINK :
1560 if(NbrHar==1)
1561 valtmp[0] = Val[0] * Equ_P->Case.Link.Coef ;
1562 else{
1563 valtmp[0] = Val[0] * Equ_P->Case.Link.Coef ;
1564 valtmp[1] = Val[1] * Equ_P->Case.Link.Coef ;
1565 }
1566 Dof_AssembleInMat(Equ_P->Case.Link.Dof, Dof_P, NbrHar, valtmp, Mat, Vec, Vecs) ;
1567 break ;
1568
1569 case DOF_LINKCPLX :
1570 if(NbrHar==1)
1571 Message::Error("LinkCplx only valid for Complex systems") ;
1572 else{ /* Warning: conjugate! */
1573 valtmp[0] = Val[0] * Equ_P->Case.Link.Coef + Val[1] * Equ_P->Case.Link.Coef2 ;
1574 valtmp[1] = Val[1] * Equ_P->Case.Link.Coef - Val[0] * Equ_P->Case.Link.Coef2 ;
1575 }
1576 Dof_AssembleInMat(Equ_P->Case.Link.Dof, Dof_P, NbrHar, valtmp, Mat, Vec, Vecs) ;
1577 break ;
1578
1579 }
1580 }
1581
1582 /* ------------------------------------------------------------------------ */
1583 /* D o f _ A s s e m b l e I n V e c */
1584 /* ------------------------------------------------------------------------ */
1585
Dof_AssembleInVec(struct Dof * Equ_P,struct Dof * Dof_P,int NbrHar,double * Val,struct Solution * OtherSolution,gVector * Vec0,gVector * Vec)1586 void Dof_AssembleInVec(struct Dof * Equ_P, struct Dof * Dof_P, int NbrHar,
1587 double * Val, struct Solution * OtherSolution,
1588 gVector * Vec0, gVector * Vec)
1589 {
1590 gScalar tmp ;
1591 double valtmp[2] ;
1592 double a, b, c, d ;
1593
1594 switch (Equ_P->Type) {
1595 case DOF_UNKNOWN :
1596 case DOF_FIXEDWITHASSOCIATE :
1597
1598 switch (Dof_P->Type) {
1599
1600 case DOF_UNKNOWN :
1601 if(NbrHar==1){
1602 if(Val[0]){
1603 LinAlg_GetDoubleInVector(&a, Vec0, Dof_P->Case.Unknown.NumDof-1) ;
1604 a *= Val[0] ;
1605 LinAlg_AddDoubleInVector(a, Vec, Equ_P->Case.Unknown.NumDof-1) ;
1606 }
1607 }
1608 else{
1609 LinAlg_GetComplexInVector(&a, &b, Vec0,
1610 Dof_P->Case.Unknown.NumDof-1,
1611 (gSCALAR_SIZE==1)?((Dof_P+1)->Case.Unknown.NumDof-1):-1) ;
1612 c = a * Val[0] - b * Val[1] ;
1613 d = a * Val[1] + b * Val[0] ;
1614 LinAlg_AddComplexInVector(c, d, Vec,
1615 Equ_P->Case.Unknown.NumDof-1,
1616 (gSCALAR_SIZE==1)?((Equ_P+1)->Case.Unknown.NumDof-1):-1) ;
1617 }
1618 break ;
1619
1620 case DOF_FIXED :
1621 case DOF_FIXEDWITHASSOCIATE :
1622 if(NbrHar==1){
1623 if(Val[0]){
1624 LinAlg_ProdScalarDouble
1625 (&Dof_P->Val,
1626 Val[0] * OtherSolution->
1627 TimeFunctionValues[Dof_P->Case.FixedAssociate.TimeFunctionIndex],
1628 &tmp) ;
1629 LinAlg_AddScalarInVector(&tmp, Vec, Equ_P->Case.Unknown.NumDof-1) ;
1630 }
1631 }
1632 else{
1633 if(gSCALAR_SIZE == 2){
1634 LinAlg_ProdScalarComplex
1635 (&Dof_P->Val,
1636 Val[0] * OtherSolution->
1637 TimeFunctionValues[Dof_P->Case.FixedAssociate.TimeFunctionIndex],
1638 Val[1] * OtherSolution->
1639 TimeFunctionValues[Dof_P->Case.FixedAssociate.TimeFunctionIndex],
1640 &a, &b) ;
1641 LinAlg_AddComplexInVector(a, b, Vec,
1642 Equ_P->Case.Unknown.NumDof-1,
1643 (gSCALAR_SIZE==1)?((Equ_P+1)->Case.Unknown.NumDof-1):-1) ;
1644 }
1645 else{
1646 Message::Error("Assemby in vectors with more than one harmonic not yet implemented") ;
1647 }
1648 }
1649 break ;
1650
1651 case DOF_LINK :
1652 if(NbrHar==1)
1653 valtmp[0] = Val[0] * Dof_P->Case.Link.Coef ;
1654 else{
1655 valtmp[0] = Val[0] * Dof_P->Case.Link.Coef ;
1656 valtmp[1] = Val[1] * Dof_P->Case.Link.Coef ;
1657 }
1658 Dof_AssembleInVec(Equ_P, Dof_P->Case.Link.Dof, NbrHar,
1659 valtmp, OtherSolution, Vec0, Vec) ;
1660 break ;
1661
1662 case DOF_LINKCPLX :
1663 if(NbrHar==1)
1664 Message::Error("LinkCplx only valid for Complex systems") ;
1665 else{
1666 valtmp[0] = Val[0] * Dof_P->Case.Link.Coef - Val[1] * Dof_P->Case.Link.Coef2 ;
1667 valtmp[1] = Val[1] * Dof_P->Case.Link.Coef + Val[0] * Dof_P->Case.Link.Coef2 ;
1668 }
1669 Dof_AssembleInVec(Equ_P, Dof_P->Case.Link.Dof, NbrHar,
1670 valtmp, OtherSolution, Vec0, Vec) ;
1671 break ;
1672
1673 case DOF_FIXED_SOLVE : case DOF_FIXEDWITHASSOCIATE_SOLVE :
1674 Message::Error("Wrong Constraints: "
1675 "remaining Dof(s) waiting to be fixed by a Resolution");
1676 break;
1677
1678 case DOF_UNKNOWN_INIT :
1679 Message::Error("Wrong Initial Constraints: "
1680 "remaining Dof(s) with non-fixed initial conditions");
1681 break;
1682 }
1683 break ;
1684
1685 case DOF_LINK :
1686 if(NbrHar==1)
1687 valtmp[0] = Val[0] * Equ_P->Case.Link.Coef ;
1688 else{
1689 valtmp[0] = Val[0] * Equ_P->Case.Link.Coef ;
1690 valtmp[1] = Val[1] * Equ_P->Case.Link.Coef ;
1691 }
1692 Dof_AssembleInVec(Equ_P->Case.Link.Dof, Dof_P, NbrHar,
1693 valtmp, OtherSolution, Vec0, Vec) ;
1694 break ;
1695
1696 case DOF_LINKCPLX :
1697 if(NbrHar==1)
1698 Message::Error("LinkCplx only valid for Complex systems") ;
1699 else{ /* Warning: conjugate! */
1700 valtmp[0] = Val[0] * Equ_P->Case.Link.Coef + Val[1] * Equ_P->Case.Link.Coef2 ;
1701 valtmp[1] = Val[1] * Equ_P->Case.Link.Coef - Val[0] * Equ_P->Case.Link.Coef2 ;
1702 }
1703 Dof_AssembleInVec(Equ_P->Case.Link.Dof, Dof_P, NbrHar,
1704 valtmp, OtherSolution, Vec0, Vec) ;
1705 break ;
1706
1707 }
1708 }
1709
1710 /* ------------------------------------------------------------------------ */
1711 /* D o f _ T r a n s f e r S o l u t i o n T o C o n s t r a i n t */
1712 /* ------------------------------------------------------------------------ */
1713
Dof_TransferSolutionToConstraint(struct DofData * DofData_P)1714 void Dof_TransferSolutionToConstraint(struct DofData * DofData_P)
1715 {
1716 struct Dof * Dof_P, * Dof_P0 ;
1717 int i ;
1718
1719 if(!DofData_P->NbrAnyDof){
1720 return;
1721 }
1722
1723 Dof_P0 = (struct Dof *)List_Pointer(DofData_P->DofList, 0) ;
1724
1725 for(i = 0 ; i < DofData_P->NbrAnyDof ; i++) {
1726 Dof_P = Dof_P0 + i ;
1727
1728 switch (Dof_P->Type) {
1729
1730 case DOF_UNKNOWN :
1731 Dof_P->Type = DOF_FIXED ;
1732 LinAlg_GetScalarInVector(&Dof_P->Val,
1733 &DofData_P->CurrentSolution->x,
1734 Dof_P->Case.Unknown.NumDof-1) ;
1735 Dof_P->Case.FixedAssociate.TimeFunctionIndex = 0 ;
1736 break ;
1737
1738 case DOF_FIXED :
1739 case DOF_FIXEDWITHASSOCIATE :
1740 case DOF_LINK :
1741 case DOF_LINKCPLX :
1742 break ;
1743
1744 default : break ;
1745 }
1746 }
1747
1748 DofData_P->NbrDof = 0 ;
1749 }
1750
1751 /* ------------------------------------------------------------------------ */
1752 /* D o f _ G e t D o f V a l u e */
1753 /* ------------------------------------------------------------------------ */
1754
Dof_GetDofValue(struct DofData * DofData_P,struct Dof * Dof_P)1755 gScalar Dof_GetDofValue(struct DofData * DofData_P, struct Dof * Dof_P)
1756 {
1757 gScalar tmp ;
1758
1759 switch (Dof_P->Type) {
1760
1761 case DOF_UNKNOWN :
1762 if(!DofData_P->CurrentSolution->SolutionExist)
1763 Message::Error("Empty solution in DofData %d", DofData_P->Num);
1764 else
1765 LinAlg_GetScalarInVector(&tmp, &DofData_P->CurrentSolution->x,
1766 Dof_P->Case.Unknown.NumDof-1) ;
1767 break ;
1768
1769 case DOF_FIXED :
1770 case DOF_FIXEDWITHASSOCIATE :
1771 LinAlg_ProdScalarDouble(&Dof_P->Val,
1772 ((Dof_P->Case.FixedAssociate.TimeFunctionIndex)?
1773 DofData_P->CurrentSolution->TimeFunctionValues
1774 [Dof_P->Case.FixedAssociate.TimeFunctionIndex] :
1775 1.),
1776 &tmp);
1777 break ;
1778
1779 case DOF_LINK :
1780 tmp = Dof_GetDofValue(DofData_P, Dof_P->Case.Link.Dof) ;
1781 LinAlg_ProdScalarDouble(&tmp, Dof_P->Case.Link.Coef, &tmp) ;
1782 break ;
1783
1784 case DOF_LINKCPLX :
1785 /* Too soon to treat LinkCplx: we need the real and imaginary parts */
1786 Message::Error("Cannot call Dof_GetDofValue for LinkCplx");
1787 break ;
1788
1789 default :
1790 LinAlg_ZeroScalar(&tmp) ;
1791 break ;
1792 }
1793
1794 return tmp ;
1795 }
1796
Dof_GetRealDofValue(struct DofData * DofData_P,struct Dof * Dof_P,double * d)1797 void Dof_GetRealDofValue(struct DofData * DofData_P, struct Dof * Dof_P, double *d)
1798 {
1799 gScalar tmp ;
1800
1801 if(Dof_P->Type == DOF_LINKCPLX) {
1802 Message::Error("Cannot call Dof_GetRealDofValue for LinkCplx");
1803 return;
1804 }
1805
1806 tmp = Dof_GetDofValue(DofData_P, Dof_P) ;
1807 LinAlg_GetDoubleInScalar(d, &tmp) ;
1808 }
1809
Dof_GetComplexDofValue(struct DofData * DofData_P,struct Dof * Dof_P,double * d1,double * d2)1810 void Dof_GetComplexDofValue(struct DofData * DofData_P, struct Dof * Dof_P,
1811 double *d1, double *d2)
1812 {
1813 gScalar tmp1, tmp2 ;
1814 double valtmp[2] ;
1815
1816 if(gSCALAR_SIZE == 1){
1817 if(Dof_P->Type == DOF_LINKCPLX) { /* Can only be done here */
1818 if(Dof_P->Case.Link.Dof->Type == DOF_LINKCPLX) { /* recurse */
1819 Dof_GetComplexDofValue(DofData_P, Dof_P->Case.Link.Dof, d1, d2);
1820 }
1821 else{
1822 tmp1 = Dof_GetDofValue(DofData_P, Dof_P->Case.Link.Dof) ;
1823 tmp2 = Dof_GetDofValue(DofData_P, (Dof_P+1)->Case.Link.Dof) ;
1824 LinAlg_GetDoubleInScalar(d1, &tmp1) ;
1825 LinAlg_GetDoubleInScalar(d2, &tmp2) ;
1826 }
1827 }
1828 else{
1829 tmp1 = Dof_GetDofValue(DofData_P, Dof_P) ;
1830 tmp2 = Dof_GetDofValue(DofData_P, Dof_P+1) ;
1831 LinAlg_GetDoubleInScalar(d1, &tmp1) ;
1832 LinAlg_GetDoubleInScalar(d2, &tmp2) ;
1833 }
1834 }
1835 else{
1836 if(Dof_P->Type == DOF_LINKCPLX) { /* Can only be done here */
1837 if(Dof_P->Case.Link.Dof->Type == DOF_LINKCPLX) { /* recurse */
1838 Dof_GetComplexDofValue(DofData_P, Dof_P->Case.Link.Dof, d1, d2);
1839 }
1840 else{
1841 tmp1 = Dof_GetDofValue(DofData_P, Dof_P->Case.Link.Dof) ;
1842 LinAlg_GetComplexInScalar(d1, d2, &tmp1) ;
1843 }
1844 }
1845 else{
1846 tmp1 = Dof_GetDofValue(DofData_P, Dof_P) ;
1847 LinAlg_GetComplexInScalar(d1, d2, &tmp1) ;
1848 }
1849 }
1850
1851 if(Dof_P->Type == DOF_LINKCPLX){
1852 valtmp[0] = Dof_P->Case.Link.Coef*(*d1) - Dof_P->Case.Link.Coef2*(*d2) ;
1853 valtmp[1] = Dof_P->Case.Link.Coef*(*d2) + Dof_P->Case.Link.Coef2*(*d1) ;
1854 *d1 = valtmp[0] ;
1855 *d2 = valtmp[1] ;
1856 }
1857 }
1858
1859 /* ------------------------------------------------------------------------- */
1860 /* D o f _ D e f i n e Unknown D o f F r o m Solve o r Init D o f */
1861 /* ------------------------------------------------------------------------- */
1862
Dof_DefineUnknownDofFromSolveOrInitDof(struct DofData ** DofData_P)1863 void Dof_DefineUnknownDofFromSolveOrInitDof(struct DofData ** DofData_P)
1864 {
1865 int i, Nbr_AnyDof ;
1866 struct Dof * Dof_P ;
1867
1868 Nbr_AnyDof = List_Nbr((*DofData_P)->DofList) ;
1869
1870 for(i = 0 ; i < Nbr_AnyDof ; i++) {
1871
1872 Dof_P = (struct Dof*)List_Pointer((*DofData_P)->DofList, i) ;
1873
1874 switch (Dof_P->Type) {
1875 case DOF_FIXED_SOLVE :
1876 case DOF_FIXEDWITHASSOCIATE_SOLVE :
1877 Dof_P->Type = DOF_UNKNOWN ;
1878 Dof_P->Case.Unknown.NumDof = ++((*DofData_P)->NbrDof) ;
1879 break ;
1880 case DOF_UNKNOWN_INIT :
1881 Dof_P->Type = DOF_UNKNOWN ;
1882 break ;
1883 }
1884 }
1885 }
1886
1887 /* ------------------------------------------------------------------------ */
1888 /* D o f _ T r a n s f e r D o f */
1889 /* ------------------------------------------------------------------------ */
1890
Dof_TransferDof(struct DofData * DofData_P1,struct DofData ** DofData_P2)1891 void Dof_TransferDof(struct DofData * DofData_P1,
1892 struct DofData ** DofData_P2)
1893 {
1894 int i, Nbr_AnyDof ;
1895 struct Dof Dof, * Dof_P ;
1896 struct Solution * Solutions_P0 ;
1897
1898 Nbr_AnyDof = List_Nbr(DofData_P1->DofList) ;
1899 Solutions_P0 = (struct Solution*)List_Pointer(DofData_P1->Solutions, 0) ;
1900 DofData_P1->CurrentSolution = Solutions_P0 ;
1901
1902 for(i = 0; i < Nbr_AnyDof; i++) {
1903 Dof = *(struct Dof *)List_Pointer(DofData_P1->DofList, i) ;
1904 if((Dof_P = (struct Dof*)Tree_PQuery((*DofData_P2)->DofTree, &Dof))){
1905 switch (Dof_P->Type) {
1906 case DOF_FIXED_SOLVE :
1907 Dof_P->Type = DOF_FIXED ;
1908 Dof_P->Val = Dof_GetDofValue(DofData_P1, &Dof) ;
1909 break ;
1910 case DOF_FIXEDWITHASSOCIATE_SOLVE :
1911 Dof_P->Type = DOF_FIXEDWITHASSOCIATE ;
1912 Dof_P->Val = Dof_GetDofValue(DofData_P1, &Dof) ;
1913 break ;
1914 case DOF_UNKNOWN_INIT :
1915 /* A DOF_UNKNOWN_INIT will always use the value obtained by
1916 pre-resolution even if a simple Init contraint is given; we should
1917 introduce DOF_UNKNOWN_INIT_SOLVE */
1918 Dof_P->Val = Dof_GetDofValue(DofData_P1, &Dof) ;
1919 if((DofData_P1->CurrentSolution - Solutions_P0) > 0){
1920 DofData_P1->CurrentSolution -= 1 ;
1921 Dof_P->Val2 = Dof_GetDofValue(DofData_P1, &Dof) ;
1922 DofData_P1->CurrentSolution += 1 ;
1923 }
1924 else{
1925 LinAlg_ZeroScalar(&Dof_P->Val2);
1926 }
1927 break ;
1928 }
1929 }
1930 }
1931 }
1932
1933 /* ------------------------------------------------------------------------ */
1934 /* D o f _ I n i t D o f F o r N o D o f */
1935 /* ------------------------------------------------------------------------ */
1936
Dof_InitDofForNoDof(struct Dof * DofForNoDof,int NbrHar)1937 void Dof_InitDofForNoDof(struct Dof * DofForNoDof, int NbrHar)
1938 {
1939 int k ;
1940 double Val[2] = {1.,0.} ;
1941
1942 for(k=0 ; k<NbrHar ; k+=gSCALAR_SIZE) {
1943 int incr = (gSCALAR_SIZE == 2) ? k / 2 : k;
1944 struct Dof * D = DofForNoDof + incr;
1945 D->Type = DOF_FIXED ;
1946 LinAlg_SetScalar(&D->Val, &Val[k%2]) ;
1947 D->Case.FixedAssociate.TimeFunctionIndex = 0 ;
1948 }
1949 }
1950
1951 /* ------------------------------------------------------- */
1952 /* P r i n t _ D o f N u m b e r */
1953 /* ------------------------------------------------------- */
1954
Print_DofNumber(struct Dof * Dof_P)1955 void Print_DofNumber(struct Dof *Dof_P)
1956 {
1957 switch(Dof_P->Type){
1958 case DOF_UNKNOWN :
1959 printf("%d(%d) ", Dof_P->Case.Unknown.NumDof, Dof_P->Entity) ;
1960 break ;
1961 case DOF_FIXED :
1962 printf("Fixed(%d) ", Dof_P->Entity) ;
1963 break ;
1964 case DOF_FIXEDWITHASSOCIATE :
1965 printf("Assoc-%d ", Dof_P->Case.FixedAssociate.NumDof) ;
1966 break ;
1967 case DOF_LINK :
1968 printf("Link-");
1969 Print_DofNumber(Dof_P->Case.Link.Dof);
1970 break ;
1971 case DOF_LINKCPLX :
1972 printf("LinkCplx-");
1973 Print_DofNumber(Dof_P->Case.Link.Dof);
1974 break ;
1975 default :
1976 printf(" ? ") ;
1977 break ;
1978 }
1979 }
1980
1981 /* ------------------------------------------------------- */
1982 /* D u m m y D o f s */
1983 /* ------------------------------------------------------- */
1984
Dof_GetDummies(struct DefineSystem * DefineSystem_P,struct DofData * DofData_P)1985 void Dof_GetDummies(struct DefineSystem * DefineSystem_P, struct DofData * DofData_P)
1986 {
1987 struct Formulation * Formulation_P ;
1988 struct DefineQuantity * DefineQuantity_P ;
1989 struct FunctionSpace * FunctionSpace_P ;
1990 struct BasisFunction * BasisFunction_P ;
1991 struct GlobalQuantity * GlobalQuantity_P ;
1992 struct Dof * Dof_P ;
1993
1994 int i, j, k, l, iDof, ii, iit, iNum, iHar;
1995 int Nbr_Formulation, Index_Formulation ;
1996 int *DummyDof;
1997 double FrequencySpectrum, *Val_Pulsation;
1998
1999 if(!(Val_Pulsation = Current.DofData->Val_Pulsation)){
2000 Message::Error("Dof_GetDummies can only be used for harmonic problems");
2001 return;
2002 }
2003
2004 DummyDof = DofData_P->DummyDof = (int *)Malloc(DofData_P->NbrDof*sizeof(int));
2005 for(iDof = 0 ; iDof < DofData_P->NbrDof ; iDof++) DummyDof[iDof]=0;
2006
2007 Nbr_Formulation = List_Nbr(DefineSystem_P->FormulationIndex) ;
2008
2009 for(i = 0 ; i < Nbr_Formulation ; i++) {
2010 List_Read(DefineSystem_P->FormulationIndex, i, &Index_Formulation) ;
2011 Formulation_P = (struct Formulation*)
2012 List_Pointer(Problem_S.Formulation, Index_Formulation) ;
2013 for(j = 0 ; j < List_Nbr(Formulation_P->DefineQuantity) ; j++) {
2014 DefineQuantity_P = (struct DefineQuantity*)
2015 List_Pointer(Formulation_P->DefineQuantity, j) ;
2016 for(l = 0 ; l < List_Nbr(DefineQuantity_P->FrequencySpectrum) ; l++) {
2017 FrequencySpectrum = *(double *)List_Pointer(DefineQuantity_P->FrequencySpectrum, l) ;
2018
2019 iHar=-1;
2020 for(k = 0 ; k < Current.NbrHar/2 ; k++)
2021 if(fabs(Val_Pulsation[k]-TWO_PI*FrequencySpectrum) <= 1e-10*Val_Pulsation[k]) {
2022 iHar = 2*k; break;
2023 }
2024 if(iHar>=0) {
2025 FunctionSpace_P = (struct FunctionSpace*)
2026 List_Pointer(Problem_S.FunctionSpace, DefineQuantity_P->FunctionSpaceIndex) ;
2027
2028 for(k = 0 ; k < List_Nbr(FunctionSpace_P->BasisFunction) ; k++) {
2029 BasisFunction_P = (struct BasisFunction *)
2030 List_Pointer(FunctionSpace_P->BasisFunction, k) ;
2031 iNum = ((struct BasisFunction *)BasisFunction_P)->Num;
2032 ii=iit=0;
2033 for(iDof = 0 ; iDof < List_Nbr(DofData_P->DofList) ; iDof++) {
2034 Dof_P = (struct Dof *)List_Pointer(DofData_P->DofList, iDof) ;
2035 if(Dof_P->Type == DOF_UNKNOWN && Dof_P->NumType == iNum) {
2036 iit++;
2037 if(Dof_P->Harmonic == iHar || Dof_P->Harmonic == iHar+1) {
2038 DummyDof[Dof_P->Case.Unknown.NumDof-1]=1; ii++;
2039 }
2040 }
2041 }
2042 if(ii)
2043 Message::Info("Freq %4lg (%d/%d) Formulation %s Quantity %s "
2044 "(BF %d) #DofsNotInFreqSpectrum %d/%d",
2045 Val_Pulsation[iHar/2]/TWO_PI, iHar/2, Current.NbrHar/2,
2046 Formulation_P->Name, DefineQuantity_P->Name,
2047 ((struct BasisFunction *)BasisFunction_P)->Num, ii, iit) ;
2048 }
2049
2050 for(k = 0 ; k < List_Nbr(FunctionSpace_P->GlobalQuantity) ; k++) {
2051 GlobalQuantity_P = (struct GlobalQuantity *)
2052 List_Pointer(FunctionSpace_P->GlobalQuantity, k) ;
2053 iNum = ((struct GlobalQuantity *)GlobalQuantity_P)->Num;
2054 ii=iit=0;
2055 for(iDof = 0 ; iDof < List_Nbr(DofData_P->DofList) ; iDof++) {
2056 Dof_P = (struct Dof *)List_Pointer(DofData_P->DofList, iDof) ;
2057 if(Dof_P->Type == DOF_UNKNOWN && Dof_P->NumType == iNum) {
2058 iit++;
2059 if(Dof_P->Harmonic == iHar || Dof_P->Harmonic == iHar+1) {
2060 DummyDof[Dof_P->Case.Unknown.NumDof-1]=1; ii++;
2061 }
2062 }
2063 }
2064 if(ii)
2065 Message::Info("Freq %4lg (%d/%d) Formulation %s GlobalQuantity %s "
2066 "(BF %d) #DofsNotInFrequencySpectrum %d/%d",
2067 Val_Pulsation[iHar/2]/TWO_PI, iHar/2, Current.NbrHar/2,
2068 Formulation_P->Name, GlobalQuantity_P->Name,
2069 ((struct GlobalQuantity *)GlobalQuantity_P)->Num, ii, iit) ;
2070 }
2071
2072 } /* end FrequencySpectrum in DofData */
2073 } /* end FrequencySpectrum in Quantity */
2074 } /* end Quantity */
2075 } /* end Formulation */
2076
2077 i=0;
2078 for(iDof = 0 ; iDof < DofData_P->NbrDof ; iDof++) {
2079 if(DummyDof[iDof]) i++;
2080
2081 if(Message::GetVerbosity() == 99){
2082 Dof_P = (struct Dof *)List_Pointer(DofData_P->DofList, iDof) ;
2083 Message::Debug("Dof Num iHar, Entity %d %d %d",
2084 iDof, Dof_P->NumType, Dof_P->Harmonic, Dof_P->Entity);
2085 }
2086 }
2087
2088 Message::Info("N: %d - N with FrequencySpectrum: %d", DofData_P->NbrDof, i) ;
2089 }
2090