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 //
9
10 #include "ProData.h"
11 #include "DofData.h"
12 #include "Message.h"
13 #include <math.h>
14
15 #define SQU(a) ((a)*(a))
16
17 extern struct CurrentData Current ;
18
19 static int Warning_Dt = 0, Warning_DtStatic = 0 ;
20 static int Warning_DtDt = 0, Warning_DtDtStatic = 0, Warning_DtDtFirstOrder = 0 ;
21
22 /* ------------------------------------------------------------------------ */
23 /* No Time Derivative */
24 /* ------------------------------------------------------------------------ */
25
Cal_AssembleTerm_NoDt(struct Dof * Equ,struct Dof * Dof,double Val[])26 void Cal_AssembleTerm_NoDt(struct Dof * Equ, struct Dof * Dof, double Val[])
27 {
28 int k ;
29 double tmp[2] ;
30
31 if(Current.TypeAssembly == ASSEMBLY_SEPARATE){
32 if (!Current.DofData->Flag_Init[1]) {
33 Current.DofData->Flag_Init[1] = 1 ;
34 LinAlg_CreateMatrix(&Current.DofData->M1, &Current.DofData->Solver,
35 Current.DofData->NbrDof, Current.DofData->NbrDof) ;
36 LinAlg_CreateVector(&Current.DofData->m1, &Current.DofData->Solver,
37 Current.DofData->NbrDof) ;
38 LinAlg_ZeroMatrix(&Current.DofData->M1);
39 LinAlg_ZeroVector(&Current.DofData->m1);
40 Current.DofData->m1s = List_Create(10, 10, sizeof(gVector));
41 for(int i = 0; i < List_Nbr(Current.DofData->TimeFunctionIndex); i++){
42 gVector m;
43 LinAlg_CreateVector(&m, &Current.DofData->Solver,
44 Current.DofData->NbrDof) ;
45 LinAlg_ZeroVector(&m);
46 List_Add(Current.DofData->m1s, &m);
47 }
48 }
49 for (k = 0 ; k < Current.NbrHar ; k += 2){
50 int incr = (gSCALAR_SIZE == 2) ? k / 2 : k;
51 Dof_AssembleInMat(Equ + incr, Dof + incr, Current.NbrHar, &Val[k],
52 &Current.DofData->M1, &Current.DofData->m1,
53 Current.DofData->m1s) ;
54 }
55 }
56 else {
57 if (Current.NbrHar == 1) {
58 switch (Current.TypeTime) {
59 case TIME_STATIC :
60 Dof_AssembleInMat(Equ, Dof, Current.NbrHar, &Val[0],
61 &Current.DofData->A, &Current.DofData->b) ;
62 break ;
63 case TIME_THETA :
64 tmp[0] = Val[0]*Current.Theta ;
65 Dof_AssembleInMat(Equ, Dof, Current.NbrHar, tmp,
66 &Current.DofData->A, &Current.DofData->b) ;
67 tmp[0] = Val[0]*(Current.Theta-1.) ;
68 Dof_AssembleInVec(Equ, Dof, Current.NbrHar, tmp,
69 Current.DofData->CurrentSolution-1,
70 &(Current.DofData->CurrentSolution-1)->x,
71 &Current.DofData->b) ;
72 break ;
73 case TIME_NEWMARK :
74 tmp[0] = Val[0]*Current.Beta ;
75 Dof_AssembleInMat(Equ, Dof, Current.NbrHar, tmp,
76 &Current.DofData->A, &Current.DofData->b) ;
77 tmp[0] = Val[0]*(2*Current.Beta-Current.Gamma-0.5) ;
78 Dof_AssembleInVec(Equ, Dof, Current.NbrHar, tmp,
79 Current.DofData->CurrentSolution-1,
80 &(Current.DofData->CurrentSolution-1)->x,
81 &Current.DofData->b) ;
82 tmp[0] = Val[0]*(Current.Gamma-Current.Beta-0.5) ;
83 Dof_AssembleInVec(Equ, Dof, Current.NbrHar, tmp,
84 Current.DofData->CurrentSolution-2,
85 &(Current.DofData->CurrentSolution-2)->x,
86 &Current.DofData->b) ;
87 break ;
88 case TIME_GEAR :
89 Dof_AssembleInMat(Equ, Dof, Current.NbrHar, Val,
90 &Current.DofData->A, &Current.DofData->b) ;
91 break ;
92 }
93 }
94 else {
95 for (k = 0 ; k < Current.NbrHar ; k += 2){
96 int incr = (gSCALAR_SIZE == 2) ? k / 2 : k;
97 Dof_AssembleInMat(Equ + incr, Dof + incr, Current.NbrHar, &Val[k],
98 &Current.DofData->A, &Current.DofData->b) ;
99 }
100 }
101 }
102 }
103
104 /* ------------------------------------------------------------------------ */
105 /* First order Time Derivative */
106 /* ------------------------------------------------------------------------ */
107
Cal_AssembleTerm_DtDof(struct Dof * Equ,struct Dof * Dof,double Val[])108 void Cal_AssembleTerm_DtDof(struct Dof * Equ, struct Dof * Dof, double Val[])
109 {
110 int k ;
111 double tmp[2] ;
112
113 if(Current.TypeAssembly == ASSEMBLY_SEPARATE){
114 if (!Current.DofData->Flag_Init[2]) {
115 Current.DofData->Flag_Init[2] = 1 ;
116 LinAlg_CreateMatrix(&Current.DofData->M2, &Current.DofData->Solver,
117 Current.DofData->NbrDof, Current.DofData->NbrDof) ;
118 LinAlg_CreateVector(&Current.DofData->m2, &Current.DofData->Solver,
119 Current.DofData->NbrDof) ;
120 LinAlg_ZeroMatrix(&Current.DofData->M2);
121 LinAlg_ZeroVector(&Current.DofData->m2);
122 Current.DofData->m2s = List_Create(10, 10, sizeof(gVector));
123 for(int i = 0; i < List_Nbr(Current.DofData->TimeFunctionIndex); i++){
124 gVector m;
125 LinAlg_CreateVector(&m, &Current.DofData->Solver,
126 Current.DofData->NbrDof) ;
127 LinAlg_ZeroVector(&m);
128 List_Add(Current.DofData->m2s, &m);
129 }
130 }
131 for (k = 0 ; k < Current.NbrHar ; k += 2){
132 int incr = (gSCALAR_SIZE == 2) ? k / 2 : k;
133 Dof_AssembleInMat(Equ + incr, Dof + incr, Current.NbrHar, &Val[k],
134 &Current.DofData->M2, &Current.DofData->m2,
135 Current.DofData->m2s) ;
136 }
137 }
138 else {
139 if (Current.NbrHar == 1) {
140 switch (Current.TypeTime) {
141 case TIME_STATIC :
142 if(!Warning_DtStatic){
143 Message::Info(3, "Discarded DtDof term in static analysis");
144 Warning_DtStatic = 1 ;
145 }
146 break;
147 case TIME_THETA :
148 tmp[0] = Val[0]/Current.DTime ;
149 Dof_AssembleInMat(Equ, Dof, Current.NbrHar, tmp,
150 &Current.DofData->A, &Current.DofData->b) ;
151 Dof_AssembleInVec(Equ, Dof, Current.NbrHar, tmp,
152 Current.DofData->CurrentSolution-1,
153 &(Current.DofData->CurrentSolution-1)->x,
154 &Current.DofData->b) ;
155 break ;
156 case TIME_NEWMARK :
157 tmp[0] = Val[0]*Current.Gamma/Current.DTime ;
158 Dof_AssembleInMat(Equ, Dof, Current.NbrHar, tmp,
159 &Current.DofData->A, &Current.DofData->b) ;
160 tmp[0] = Val[0]*(2.*Current.Gamma-1.)/Current.DTime ;
161 Dof_AssembleInVec(Equ, Dof, Current.NbrHar, tmp,
162 Current.DofData->CurrentSolution-1,
163 &(Current.DofData->CurrentSolution-1)->x,
164 &Current.DofData->b) ;
165 tmp[0] = Val[0]*(1.-Current.Gamma)/Current.DTime ;
166 Dof_AssembleInVec(Equ, Dof, Current.NbrHar, tmp,
167 Current.DofData->CurrentSolution-2,
168 &(Current.DofData->CurrentSolution-2)->x,
169 &Current.DofData->b) ;
170 break ;
171 case TIME_GEAR :
172 tmp[0] = Val[0]/(Current.bCorrCoeff*Current.DTime);
173 Dof_AssembleInMat(Equ, Dof, Current.NbrHar, tmp,
174 &Current.DofData->A, &Current.DofData->b) ;
175 for (int i=0; i<Current.CorrOrder; i++) {
176 tmp[0] = Val[0]*Current.aCorrCoeff[i]/(Current.bCorrCoeff*Current.DTime);
177 Dof_AssembleInVec(Equ, Dof, Current.NbrHar, tmp,
178 Current.DofData->CurrentSolution-1,
179 &(Current.DofData->CurrentSolution-1-i)->x,
180 &Current.DofData->b) ;
181 }
182 break ;
183 }
184 }
185 else {
186 for (k = 0 ; k < Current.NbrHar ; k += 2) {
187 //printf("====hola idiota=== k %d pul %g \n", k, Current.DofData->Val_Pulsation[k/2]);
188 tmp[0] = -Val[k+1] * Current.DofData->Val_Pulsation[k/2] ;
189 tmp[1] = Val[k] * Current.DofData->Val_Pulsation[k/2] ;
190 int incr = (gSCALAR_SIZE == 2) ? k / 2 : k;
191 Dof_AssembleInMat(Equ + incr, Dof + incr, Current.NbrHar, tmp,
192 &Current.DofData->A, &Current.DofData->b) ;
193 }
194 }
195 }
196 }
197
Cal_AssembleTerm_Dt(struct Dof * Equ,struct Dof * Dof,double Val[])198 void Cal_AssembleTerm_Dt(struct Dof * Equ, struct Dof * Dof, double Val[])
199 {
200 if(!Warning_Dt){
201 Message::Warning("Dt not implemented, using DtDof instead");
202 Warning_Dt = 1 ;
203 }
204 Cal_AssembleTerm_DtDof(Equ, Dof, Val);
205 }
206
207 /* En preparation ... */
Cal_AssembleTerm_DtNL(struct Dof * Equ,struct Dof * Dof,double Val[])208 void Cal_AssembleTerm_DtNL(struct Dof * Equ, struct Dof * Dof, double Val[])
209 {
210 double tmp[2] ;
211
212 if(Current.TypeAssembly == ASSEMBLY_SEPARATE){
213 Message::Error("DtNL not implemented for separate assembly");
214 }
215 else {
216 if (Current.NbrHar == 1) {
217 switch (Current.TypeTime) {
218 case TIME_STATIC :
219 if(!Warning_DtStatic){
220 Message::Info(3, "Discarded DtDof term in static analysis");
221 Warning_DtStatic = 1 ;
222 }
223 break;
224 case TIME_THETA :
225 tmp[0] = Val[0]/Current.DTime ;
226 Dof_AssembleInMat(Equ, Dof, Current.NbrHar, tmp,
227 &Current.DofData->A, &Current.DofData->b) ;
228 Dof_AssembleInVec(Equ, Dof, Current.NbrHar, tmp,
229 Current.DofData->CurrentSolution-1,
230 &(Current.DofData->CurrentSolution-1)->x,
231 &Current.DofData->b) ;
232 break ;
233 case TIME_NEWMARK :
234 Message::Error("DtNL not implemented for separate assembly with Newmark scheme");
235 return ;
236 case TIME_GEAR :
237 Message::Error("DtNL not implemented for Gear's method");
238 return ;
239 }
240 }
241 else {
242 Message::Error("DtNL not implemented for separate assembly in harmonic analysis");
243 }
244 }
245 }
246
247 /* ------------------------------------------------------------------------ */
248 /* Second order Time Derivative */
249 /* ------------------------------------------------------------------------ */
250
Cal_AssembleTerm_DtDtDof(struct Dof * Equ,struct Dof * Dof,double Val[])251 void Cal_AssembleTerm_DtDtDof(struct Dof * Equ, struct Dof * Dof, double Val[])
252 {
253 int k ;
254 double tmp[2] ;
255
256 if(Current.TypeAssembly == ASSEMBLY_SEPARATE){
257 if (!Current.DofData->Flag_Init[3]) {
258 Current.DofData->Flag_Init[3] = 1 ;
259 LinAlg_CreateMatrix(&Current.DofData->M3, &Current.DofData->Solver,
260 Current.DofData->NbrDof, Current.DofData->NbrDof) ;
261 LinAlg_CreateVector(&Current.DofData->m3, &Current.DofData->Solver,
262 Current.DofData->NbrDof) ;
263 LinAlg_ZeroMatrix(&Current.DofData->M3);
264 LinAlg_ZeroVector(&Current.DofData->m3);
265 Current.DofData->m3s = List_Create(10, 10, sizeof(gVector));
266 for(int i = 0; i < List_Nbr(Current.DofData->TimeFunctionIndex); i++){
267 gVector m;
268 LinAlg_CreateVector(&m, &Current.DofData->Solver,
269 Current.DofData->NbrDof) ;
270 LinAlg_ZeroVector(&m);
271 List_Add(Current.DofData->m3s, &m);
272 }
273 }
274 for (k = 0 ; k < Current.NbrHar ; k += 2) {
275 int incr = (gSCALAR_SIZE == 2) ? k / 2 : k;
276 Dof_AssembleInMat(Equ + incr, Dof + incr, Current.NbrHar, &Val[k],
277 &Current.DofData->M3, &Current.DofData->m3,
278 Current.DofData->m3s) ;
279 }
280 }
281 else {
282 if (Current.NbrHar == 1) {
283 switch (Current.TypeTime) {
284 case TIME_STATIC :
285 if(!Warning_DtDtStatic){
286 Message::Info(3, "Discarded DtDtDof term in static analysis");
287 Warning_DtDtStatic = 1 ;
288 }
289 break;
290 case TIME_THETA :
291 if(!Warning_DtDtFirstOrder){
292 Message::Info(3, "Discarded DtDtDof term in Theta time scheme");
293 Warning_DtDtFirstOrder = 1 ;
294 }
295 break;
296 case TIME_GEAR :
297 Message::Error("DtDtDof not implemented for Gear's method");
298 return ;
299 case TIME_NEWMARK :
300 tmp[0] = Val[0]/SQU(Current.DTime) ;
301 Dof_AssembleInMat(Equ, Dof, Current.NbrHar, tmp,
302 &Current.DofData->A, &Current.DofData->b) ;
303 tmp[0] = 2*Val[0]/SQU(Current.DTime) ;
304 Dof_AssembleInVec(Equ, Dof, Current.NbrHar, tmp,
305 Current.DofData->CurrentSolution-1,
306 &(Current.DofData->CurrentSolution-1)->x,
307 &Current.DofData->b) ;
308 tmp[0] = -Val[0]/SQU(Current.DTime) ;
309 Dof_AssembleInVec(Equ, Dof, Current.NbrHar, tmp,
310 Current.DofData->CurrentSolution-2,
311 &(Current.DofData->CurrentSolution-2)->x,
312 &Current.DofData->b) ;
313 break ;
314 }
315 }
316 else {
317 for (k = 0 ; k < Current.NbrHar ; k += 2) {
318 tmp[0] = - Val[k] * SQU(Current.DofData->Val_Pulsation[k/2]) ;
319 tmp[1] = - Val[k+1] * SQU(Current.DofData->Val_Pulsation[k/2]) ;
320 int incr = (gSCALAR_SIZE == 2) ? k / 2 : k;
321 Dof_AssembleInMat(Equ + incr, Dof + incr, Current.NbrHar, tmp,
322 &Current.DofData->A, &Current.DofData->b) ;
323 }
324 }
325 }
326 }
327
Cal_AssembleTerm_DtDt(struct Dof * Equ,struct Dof * Dof,double Val[])328 void Cal_AssembleTerm_DtDt(struct Dof * Equ, struct Dof * Dof, double Val[])
329 {
330 if(!Warning_DtDt){
331 Message::Warning("DtDt not implemented, using DtDtDof instead");
332 Warning_DtDt = 1 ;
333 }
334 Cal_AssembleTerm_DtDtDof(Equ, Dof, Val);
335 }
336
337 /* ------------------------------------------------- */
338 /* higher order Time Derivative for Polynomial EVP */
339 /* ------------------------------------------------- */
Cal_AssembleTerm_DtDtDtDof(struct Dof * Equ,struct Dof * Dof,double Val[])340 void Cal_AssembleTerm_DtDtDtDof(struct Dof * Equ, struct Dof * Dof, double Val[])
341 {
342 int k ;
343 if(Current.TypeAssembly == ASSEMBLY_SEPARATE){
344 if (!Current.DofData->Flag_Init[4]) {
345 Current.DofData->Flag_Init[4] = 1 ;
346 LinAlg_CreateMatrix(&Current.DofData->M4, &Current.DofData->Solver,
347 Current.DofData->NbrDof, Current.DofData->NbrDof) ;
348 LinAlg_CreateVector(&Current.DofData->m4, &Current.DofData->Solver,
349 Current.DofData->NbrDof) ;
350 LinAlg_ZeroMatrix(&Current.DofData->M4);
351 LinAlg_ZeroVector(&Current.DofData->m4);
352 Current.DofData->m4s = List_Create(10, 10, sizeof(gVector));
353 for(int i = 0; i < List_Nbr(Current.DofData->TimeFunctionIndex); i++){
354 gVector m;
355 LinAlg_CreateVector(&m, &Current.DofData->Solver,
356 Current.DofData->NbrDof) ;
357 LinAlg_ZeroVector(&m);
358 List_Add(Current.DofData->m4s, &m);
359 }
360 }
361 for (k = 0 ; k < Current.NbrHar ; k += 2) {
362 int incr = (gSCALAR_SIZE == 2) ? k / 2 : k;
363 Dof_AssembleInMat(Equ + incr, Dof + incr, Current.NbrHar, &Val[k],
364 &Current.DofData->M4, &Current.DofData->m4,
365 Current.DofData->m4s) ;
366 }
367 }
368 else {
369 Message::Error("DtDtDtDof only available with GenerateSeparate");
370 return ;
371 }
372 }
373
Cal_AssembleTerm_DtDtDtDtDof(struct Dof * Equ,struct Dof * Dof,double Val[])374 void Cal_AssembleTerm_DtDtDtDtDof(struct Dof * Equ, struct Dof * Dof, double Val[])
375 {
376 int k ;
377 if(Current.TypeAssembly == ASSEMBLY_SEPARATE){
378 if (!Current.DofData->Flag_Init[5]) {
379 Current.DofData->Flag_Init[5] = 1 ;
380 LinAlg_CreateMatrix(&Current.DofData->M5, &Current.DofData->Solver,
381 Current.DofData->NbrDof, Current.DofData->NbrDof) ;
382 LinAlg_CreateVector(&Current.DofData->m5, &Current.DofData->Solver,
383 Current.DofData->NbrDof) ;
384 LinAlg_ZeroMatrix(&Current.DofData->M5);
385 LinAlg_ZeroVector(&Current.DofData->m5);
386 Current.DofData->m5s = List_Create(10, 10, sizeof(gVector));
387 for(int i = 0; i < List_Nbr(Current.DofData->TimeFunctionIndex); i++){
388 gVector m;
389 LinAlg_CreateVector(&m, &Current.DofData->Solver,
390 Current.DofData->NbrDof) ;
391 LinAlg_ZeroVector(&m);
392 List_Add(Current.DofData->m5s, &m);
393 }
394 }
395 for (k = 0 ; k < Current.NbrHar ; k += 2) {
396 int incr = (gSCALAR_SIZE == 2) ? k / 2 : k;
397 Dof_AssembleInMat(Equ + incr, Dof + incr, Current.NbrHar, &Val[k],
398 &Current.DofData->M5, &Current.DofData->m5,
399 Current.DofData->m5s) ;
400 }
401 }
402 else {
403 Message::Error("DtDtDtDtDof only available with GenerateSeparate");
404 return ;
405 }
406 }
407
Cal_AssembleTerm_DtDtDtDtDtDof(struct Dof * Equ,struct Dof * Dof,double Val[])408 void Cal_AssembleTerm_DtDtDtDtDtDof(struct Dof * Equ, struct Dof * Dof, double Val[])
409 {
410 int k ;
411 if(Current.TypeAssembly == ASSEMBLY_SEPARATE){
412 if (!Current.DofData->Flag_Init[6]) {
413 Current.DofData->Flag_Init[6] = 1 ;
414 LinAlg_CreateMatrix(&Current.DofData->M6, &Current.DofData->Solver,
415 Current.DofData->NbrDof, Current.DofData->NbrDof) ;
416 LinAlg_CreateVector(&Current.DofData->m6, &Current.DofData->Solver,
417 Current.DofData->NbrDof) ;
418 LinAlg_ZeroMatrix(&Current.DofData->M6);
419 LinAlg_ZeroVector(&Current.DofData->m6);
420 Current.DofData->m6s = List_Create(10, 10, sizeof(gVector));
421 for(int i = 0; i < List_Nbr(Current.DofData->TimeFunctionIndex); i++){
422 gVector m;
423 LinAlg_CreateVector(&m, &Current.DofData->Solver,
424 Current.DofData->NbrDof) ;
425 LinAlg_ZeroVector(&m);
426 List_Add(Current.DofData->m6s, &m);
427 }
428 }
429 for (k = 0 ; k < Current.NbrHar ; k += 2) {
430 int incr = (gSCALAR_SIZE == 2) ? k / 2 : k;
431 Dof_AssembleInMat(Equ + incr, Dof + incr, Current.NbrHar, &Val[k],
432 &Current.DofData->M6, &Current.DofData->m6,
433 Current.DofData->m6s) ;
434 }
435 }
436 else {
437 Message::Error("DtDtDtDtDtDof only available with GenerateSeparate");
438 return ;
439 }
440 }
441
442 // nleigchange
Cal_AssembleTerm_NLEig1Dof(struct Dof * Equ,struct Dof * Dof,double Val[])443 void Cal_AssembleTerm_NLEig1Dof(struct Dof * Equ, struct Dof * Dof, double Val[])
444 {
445 int k ;
446 if(Current.TypeAssembly == ASSEMBLY_SEPARATE){
447 if (!Current.DofData->Flag_Init[7]) {
448 Current.DofData->Flag_Init[7] = 1 ;
449 LinAlg_CreateMatrix(&Current.DofData->M7, &Current.DofData->Solver,
450 Current.DofData->NbrDof, Current.DofData->NbrDof) ;
451 LinAlg_CreateVector(&Current.DofData->m7, &Current.DofData->Solver,
452 Current.DofData->NbrDof) ;
453 LinAlg_ZeroMatrix(&Current.DofData->M7);
454 LinAlg_ZeroVector(&Current.DofData->m7);
455 Current.DofData->m7s = List_Create(10, 10, sizeof(gVector));
456 for(int i = 0; i < List_Nbr(Current.DofData->TimeFunctionIndex); i++){
457 gVector m;
458 LinAlg_CreateVector(&m, &Current.DofData->Solver,
459 Current.DofData->NbrDof) ;
460 LinAlg_ZeroVector(&m);
461 List_Add(Current.DofData->m7s, &m);
462 }
463 }
464 for (k = 0 ; k < Current.NbrHar ; k += 2) {
465 int incr = (gSCALAR_SIZE == 2) ? k / 2 : k;
466 Dof_AssembleInMat(Equ + incr, Dof + incr, Current.NbrHar, &Val[k],
467 &Current.DofData->M7, &Current.DofData->m7,
468 Current.DofData->m7s) ;
469 }
470 }
471 else {
472 Message::Error("NLEig1Dof only available with GenerateSeparate");
473 return ;
474 }
475 }
476
Cal_AssembleTerm_NLEig2Dof(struct Dof * Equ,struct Dof * Dof,double Val[])477 void Cal_AssembleTerm_NLEig2Dof(struct Dof * Equ, struct Dof * Dof, double Val[])
478 {
479 int k ;
480 if(Current.TypeAssembly == ASSEMBLY_SEPARATE){
481 if (!Current.DofData->Flag_Init[6]) {
482 Current.DofData->Flag_Init[6] = 1 ;
483 LinAlg_CreateMatrix(&Current.DofData->M6, &Current.DofData->Solver,
484 Current.DofData->NbrDof, Current.DofData->NbrDof) ;
485 LinAlg_CreateVector(&Current.DofData->m6, &Current.DofData->Solver,
486 Current.DofData->NbrDof) ;
487 LinAlg_ZeroMatrix(&Current.DofData->M6);
488 LinAlg_ZeroVector(&Current.DofData->m6);
489 Current.DofData->m6s = List_Create(10, 10, sizeof(gVector));
490 for(int i = 0; i < List_Nbr(Current.DofData->TimeFunctionIndex); i++){
491 gVector m;
492 LinAlg_CreateVector(&m, &Current.DofData->Solver,
493 Current.DofData->NbrDof) ;
494 LinAlg_ZeroVector(&m);
495 List_Add(Current.DofData->m6s, &m);
496 }
497 }
498 for (k = 0 ; k < Current.NbrHar ; k += 2) {
499 int incr = (gSCALAR_SIZE == 2) ? k / 2 : k;
500 Dof_AssembleInMat(Equ + incr, Dof + incr, Current.NbrHar, &Val[k],
501 &Current.DofData->M6, &Current.DofData->m6,
502 Current.DofData->m6s) ;
503 }
504 }
505 else {
506 Message::Error("NLEig3Dof only available with GenerateSeparate");
507 return ;
508 }
509 }
510
Cal_AssembleTerm_NLEig3Dof(struct Dof * Equ,struct Dof * Dof,double Val[])511 void Cal_AssembleTerm_NLEig3Dof(struct Dof * Equ, struct Dof * Dof, double Val[])
512 {
513 int k ;
514 if(Current.TypeAssembly == ASSEMBLY_SEPARATE){
515 if (!Current.DofData->Flag_Init[5]) {
516 Current.DofData->Flag_Init[5] = 1 ;
517 LinAlg_CreateMatrix(&Current.DofData->M5, &Current.DofData->Solver,
518 Current.DofData->NbrDof, Current.DofData->NbrDof) ;
519 LinAlg_CreateVector(&Current.DofData->m5, &Current.DofData->Solver,
520 Current.DofData->NbrDof) ;
521 LinAlg_ZeroMatrix(&Current.DofData->M5);
522 LinAlg_ZeroVector(&Current.DofData->m5);
523 Current.DofData->m5s = List_Create(10, 10, sizeof(gVector));
524 for(int i = 0; i < List_Nbr(Current.DofData->TimeFunctionIndex); i++){
525 gVector m;
526 LinAlg_CreateVector(&m, &Current.DofData->Solver,
527 Current.DofData->NbrDof) ;
528 LinAlg_ZeroVector(&m);
529 List_Add(Current.DofData->m5s, &m);
530 }
531 }
532 for (k = 0 ; k < Current.NbrHar ; k += 2) {
533 int incr = (gSCALAR_SIZE == 2) ? k / 2 : k;
534 Dof_AssembleInMat(Equ + incr, Dof + incr, Current.NbrHar, &Val[k],
535 &Current.DofData->M5, &Current.DofData->m5,
536 Current.DofData->m5s) ;
537 }
538 }
539 else {
540 Message::Error("NLEig3Dof only available with GenerateSeparate");
541 return ;
542 }
543 }
544
Cal_AssembleTerm_NLEig4Dof(struct Dof * Equ,struct Dof * Dof,double Val[])545 void Cal_AssembleTerm_NLEig4Dof(struct Dof * Equ, struct Dof * Dof, double Val[])
546 {
547 int k ;
548 if(Current.TypeAssembly == ASSEMBLY_SEPARATE){
549 if (!Current.DofData->Flag_Init[4]) {
550 Current.DofData->Flag_Init[4] = 1 ;
551 LinAlg_CreateMatrix(&Current.DofData->M4, &Current.DofData->Solver,
552 Current.DofData->NbrDof, Current.DofData->NbrDof) ;
553 LinAlg_CreateVector(&Current.DofData->m4, &Current.DofData->Solver,
554 Current.DofData->NbrDof) ;
555 LinAlg_ZeroMatrix(&Current.DofData->M4);
556 LinAlg_ZeroVector(&Current.DofData->m4);
557 Current.DofData->m4s = List_Create(10, 10, sizeof(gVector));
558 for(int i = 0; i < List_Nbr(Current.DofData->TimeFunctionIndex); i++){
559 gVector m;
560 LinAlg_CreateVector(&m, &Current.DofData->Solver,
561 Current.DofData->NbrDof) ;
562 LinAlg_ZeroVector(&m);
563 List_Add(Current.DofData->m4s, &m);
564 }
565 }
566 for (k = 0 ; k < Current.NbrHar ; k += 2) {
567 int incr = (gSCALAR_SIZE == 2) ? k / 2 : k;
568 Dof_AssembleInMat(Equ + incr, Dof + incr, Current.NbrHar, &Val[k],
569 &Current.DofData->M4, &Current.DofData->m4,
570 Current.DofData->m4s) ;
571 }
572 }
573 else {
574 Message::Error("NLEig4Dof only available with GenerateSeparate");
575 return ;
576 }
577 }
578
Cal_AssembleTerm_NLEig5Dof(struct Dof * Equ,struct Dof * Dof,double Val[])579 void Cal_AssembleTerm_NLEig5Dof(struct Dof * Equ, struct Dof * Dof, double Val[])
580 {
581 int k ;
582 if(Current.TypeAssembly == ASSEMBLY_SEPARATE){
583 if (!Current.DofData->Flag_Init[3]) {
584 Current.DofData->Flag_Init[3] = 1 ;
585 LinAlg_CreateMatrix(&Current.DofData->M3, &Current.DofData->Solver,
586 Current.DofData->NbrDof, Current.DofData->NbrDof) ;
587 LinAlg_CreateVector(&Current.DofData->m3, &Current.DofData->Solver,
588 Current.DofData->NbrDof) ;
589 LinAlg_ZeroMatrix(&Current.DofData->M3);
590 LinAlg_ZeroVector(&Current.DofData->m3);
591 Current.DofData->m3s = List_Create(10, 10, sizeof(gVector));
592 for(int i = 0; i < List_Nbr(Current.DofData->TimeFunctionIndex); i++){
593 gVector m;
594 LinAlg_CreateVector(&m, &Current.DofData->Solver,
595 Current.DofData->NbrDof) ;
596 LinAlg_ZeroVector(&m);
597 List_Add(Current.DofData->m3s, &m);
598 }
599 }
600 for (k = 0 ; k < Current.NbrHar ; k += 2) {
601 int incr = (gSCALAR_SIZE == 2) ? k / 2 : k;
602 Dof_AssembleInMat(Equ + incr, Dof + incr, Current.NbrHar, &Val[k],
603 &Current.DofData->M3, &Current.DofData->m3,
604 Current.DofData->m3s) ;
605 }
606 }
607 else {
608 Message::Error("NLEig5Dof only available with GenerateSeparate");
609 return ;
610 }
611 }
612
Cal_AssembleTerm_NLEig6Dof(struct Dof * Equ,struct Dof * Dof,double Val[])613 void Cal_AssembleTerm_NLEig6Dof(struct Dof * Equ, struct Dof * Dof, double Val[])
614 {
615 int k ;
616 if(Current.TypeAssembly == ASSEMBLY_SEPARATE){
617 if (!Current.DofData->Flag_Init[2]) {
618 Current.DofData->Flag_Init[2] = 1 ;
619 LinAlg_CreateMatrix(&Current.DofData->M2, &Current.DofData->Solver,
620 Current.DofData->NbrDof, Current.DofData->NbrDof) ;
621 LinAlg_CreateVector(&Current.DofData->m2, &Current.DofData->Solver,
622 Current.DofData->NbrDof) ;
623 LinAlg_ZeroMatrix(&Current.DofData->M2);
624 LinAlg_ZeroVector(&Current.DofData->m2);
625 Current.DofData->m2s = List_Create(10, 10, sizeof(gVector));
626 for(int i = 0; i < List_Nbr(Current.DofData->TimeFunctionIndex); i++){
627 gVector m;
628 LinAlg_CreateVector(&m, &Current.DofData->Solver,
629 Current.DofData->NbrDof) ;
630 LinAlg_ZeroVector(&m);
631 List_Add(Current.DofData->m2s, &m);
632 }
633 }
634 for (k = 0 ; k < Current.NbrHar ; k += 2) {
635 int incr = (gSCALAR_SIZE == 2) ? k / 2 : k;
636 Dof_AssembleInMat(Equ + incr, Dof + incr, Current.NbrHar, &Val[k],
637 &Current.DofData->M2, &Current.DofData->m2,
638 Current.DofData->m2s) ;
639 }
640 }
641 else {
642 Message::Error("NLEig6Dof only available with GenerateSeparate");
643 return ;
644 }
645 }
646
647 /* ------------------------------------------------------------------------ */
648 /* Jacobian NonLinear */
649 /* ------------------------------------------------------------------------ */
650
Cal_AssembleTerm_JacNL(struct Dof * Equ,struct Dof * Dof,double Val[])651 void Cal_AssembleTerm_JacNL(struct Dof * Equ, struct Dof * Dof, double Val[])
652 {
653 int k ;
654
655 if(Current.TypeAssembly == ASSEMBLY_SEPARATE){
656 Message::Error("JacNL not implemented for separate assembly");
657 }
658 else{
659 if (Current.NbrHar == 1) {
660 switch (Current.TypeTime) {
661 case TIME_STATIC :
662 case TIME_THETA :
663 Dof_AssembleInMat(Equ, Dof, Current.NbrHar, &Val[0],
664 &Current.DofData->Jac, NULL) ;
665 break ;
666 case TIME_GEAR :
667 Dof_AssembleInMat(Equ, Dof, Current.NbrHar, &Val[0],
668 &Current.DofData->Jac, NULL) ;
669 break ;
670 case TIME_NEWMARK :
671 Message::Error("JacNL not implemented for Newmark's method");
672 return ;
673 }
674 }
675 else {
676 for (k = 0 ; k < Current.NbrHar ; k += 2){
677 int incr = (gSCALAR_SIZE == 2) ? k / 2 : k;
678 Dof_AssembleInMat(Equ + incr, Dof + incr, Current.NbrHar, &Val[k],
679 &Current.DofData->Jac, NULL) ;
680 }
681 }
682 }
683 }
684
Cal_AssembleTerm_DtDofJacNL(struct Dof * Equ,struct Dof * Dof,double Val[])685 void Cal_AssembleTerm_DtDofJacNL(struct Dof * Equ, struct Dof * Dof, double Val[])
686 {
687 double tmp[2] ;
688
689 if(Current.TypeAssembly == ASSEMBLY_SEPARATE)
690 Message::Error("DtDofJacNL not implemented for separate assembly");
691 else {
692 if (Current.NbrHar == 1) {
693 switch (Current.TypeTime) {
694 case TIME_STATIC :
695 if(!Warning_DtStatic){
696 Message::Info(3, "Discarded DtDofJacNL term in static analysis");
697 Warning_DtStatic = 1 ;
698 }
699 break;
700 case TIME_THETA :
701 if ( fabs(Current.Theta - 1.0) > 1e-3 ){
702 Message::Error("Theta method not implemented for nonlinear problems when "
703 "Theta != 1.0");
704 return;
705 }
706 tmp[0] = Val[0]/Current.DTime;
707 Dof_AssembleInMat(Equ, Dof, Current.NbrHar, tmp,
708 &Current.DofData->Jac, NULL);
709 break ;
710 case TIME_NEWMARK :
711 Message::Error("DtDofJacNL not implemented for Newmark scheme");
712 return ;
713 case TIME_GEAR :
714 tmp[0] = Val[0] / (Current.bCorrCoeff * Current.DTime);
715 Dof_AssembleInMat(Equ, Dof, Current.NbrHar, tmp,
716 &Current.DofData->Jac, NULL);
717 break ;
718 }
719 }
720 else
721 Message::Error("DtDofJacNL not implemented for harmonic analysis");
722 }
723 }
724
725 /* ------------------------------------------------------------------------ */
726 /* Never Time Derivative (provisoire mais tres important ... Patrick) */
727 /* ------------------------------------------------------------------------ */
728
Cal_AssembleTerm_NeverDt(struct Dof * Equ,struct Dof * Dof,double Val[])729 void Cal_AssembleTerm_NeverDt(struct Dof * Equ, struct Dof * Dof, double Val[])
730 {
731 int k ;
732
733 if(Current.TypeAssembly == ASSEMBLY_SEPARATE){
734 if (!Current.DofData->Flag_Init[1]) {
735 Current.DofData->Flag_Init[1] = 1 ;
736 LinAlg_CreateMatrix(&Current.DofData->M1, &Current.DofData->Solver,
737 Current.DofData->NbrDof, Current.DofData->NbrDof) ;
738 LinAlg_CreateVector(&Current.DofData->m1, &Current.DofData->Solver,
739 Current.DofData->NbrDof) ;
740 LinAlg_ZeroMatrix(&Current.DofData->M1);
741 LinAlg_ZeroVector(&Current.DofData->m1);
742 Current.DofData->m1s = List_Create(10, 10, sizeof(gVector));
743 for(int i = 0; i < List_Nbr(Current.DofData->TimeFunctionIndex); i++){
744 gVector m;
745 LinAlg_CreateVector(&m, &Current.DofData->Solver,
746 Current.DofData->NbrDof) ;
747 LinAlg_ZeroVector(&m);
748 List_Add(Current.DofData->m1s, &m);
749 }
750 }
751 for (k = 0 ; k < Current.NbrHar ; k += 2){
752 int incr = (gSCALAR_SIZE == 2) ? k / 2 : k;
753 Dof_AssembleInMat(Equ + incr, Dof + incr, Current.NbrHar, &Val[k],
754 &Current.DofData->M1, &Current.DofData->m1,
755 Current.DofData->m1s) ;
756 }
757 }
758 else{
759 if (Current.NbrHar == 1) {
760 switch (Current.TypeTime) {
761 case TIME_STATIC :
762 case TIME_THETA :
763 case TIME_NEWMARK :
764 case TIME_GEAR :
765 Dof_AssembleInMat(Equ, Dof, Current.NbrHar, &Val[0],
766 &Current.DofData->A, &Current.DofData->b) ;
767 break ;
768 }
769 }
770 else {
771 for (k = 0 ; k < Current.NbrHar ; k += 2) {
772 int incr = (gSCALAR_SIZE == 2) ? k / 2 : k;
773 Dof_AssembleInMat(Equ + incr, Dof + incr, Current.NbrHar, &Val[k],
774 &Current.DofData->A, &Current.DofData->b) ;
775 }
776 }
777 }
778 }
779
780 /* ------------------------------------------------------------------------ */
781 /* Multi-Harmonic with movement */
782 /* ------------------------------------------------------------------------ */
783
Cal_AssembleTerm_MHMoving(struct Dof * Equ,struct Dof * Dof,double Val[])784 void Cal_AssembleTerm_MHMoving(struct Dof * Equ, struct Dof * Dof, double Val[])
785 {
786 // MHMoving_assemblyType = 1 => Use current system A,b
787 // MHMoving_assemblyType = 2 => Use dedicated system A_MH_Moving, b_MH_Moving
788 // MHMoving_assemblyType = 3 => Look for unknowns and constraints in Moving Group
789
790 extern int MHMoving_assemblyType ;
791 extern double ** MH_Moving_Matrix ;
792 extern Tree_T * DofTree_MH_moving ;
793
794 // FIXME: this cannot work in complex arithmetic: AssembleInMat will
795 // need to assemble both real and imaginary parts at once -- See
796 // Cal_AssembleTerm for an example
797
798 if(MHMoving_assemblyType==1){
799 for (int k = 0 ; k < Current.NbrHar ; k++)
800 for (int l = 0 ; l < Current.NbrHar ; l++) {
801 double tmp = Val[0] * MH_Moving_Matrix[k][l] ;
802 // if (k==l)
803 Dof_AssembleInMat(Equ+k, Dof+l, 1, &tmp,
804 &Current.DofData->A, &Current.DofData->b) ;
805 }
806 }
807
808 if(MHMoving_assemblyType==2){
809 for (int k = 0 ; k < Current.NbrHar ; k++)
810 for (int l = 0 ; l < Current.NbrHar ; l++) {
811 double tmp = Val[0] * MH_Moving_Matrix[k][l] ;
812 // if (k==l)
813 Dof_AssembleInMat(Equ+k, Dof+l, 1, &tmp, &Current.DofData->A_MH_moving, NULL);
814 // &Current.DofData->A_MH_moving, &Current.DofData->b_MH_moving) ;
815 }
816 }
817
818 if(MHMoving_assemblyType==3){
819 if (Dof->Type == DOF_UNKNOWN && !Tree_PQuery(DofTree_MH_moving, Dof))
820 Tree_Add(DofTree_MH_moving,Dof) ;
821 else if (Dof->Type == DOF_LINK && !Tree_PQuery(DofTree_MH_moving, Dof->Case.Link.Dof))
822 Tree_Add(DofTree_MH_moving,Dof->Case.Link.Dof) ;
823
824 if (Equ->Type == DOF_UNKNOWN && !Tree_PQuery(DofTree_MH_moving, Equ))
825 Tree_Add(DofTree_MH_moving,Equ) ;
826 else if (Equ->Type == DOF_LINK && !Tree_PQuery(DofTree_MH_moving, Equ->Case.Link.Dof))
827 Tree_Add(DofTree_MH_moving,Equ->Case.Link.Dof) ;
828 }
829 }
830
831 /*
832 void Cal_AssembleTerm_MH_Moving_simple(struct Dof * Equ, struct Dof * Dof, double Val[])
833 {
834 int k, l ;
835 double tmp ;
836 extern double ** MH_Moving_Matrix ;
837
838 for (k = 0 ; k < Current.NbrHar ; k++)
839 for (l = 0 ; l < Current.NbrHar ; l++) {
840 tmp = Val[0] * MH_Moving_Matrix[k][l] ;
841 // if (k==l)
842 Dof_AssembleInMat(Equ+k, Dof+l, 1, &tmp,
843 &Current.DofData->A, &Current.DofData->b) ;
844 }
845 }
846
847 void Cal_AssembleTerm_MH_Moving_separate(struct Dof * Equ, struct Dof * Dof, double Val[])
848 {
849 int k, l ;
850 double tmp ;
851 extern double ** MH_Moving_Matrix ;
852
853 for (k = 0 ; k < Current.NbrHar ; k++)
854 for (l = 0 ; l < Current.NbrHar ; l++) {
855 tmp = Val[0] * MH_Moving_Matrix[k][l] ;
856 // if (k==l)
857 Dof_AssembleInMat(Equ+k, Dof+l, 1, &tmp,
858 &Current.DofData->A_MH_moving, &Current.DofData->b_MH_moving) ;
859 }
860 }
861
862 void Cal_AssembleTerm_MH_Moving_probe(struct Dof * Equ, struct Dof * Dof, double Val[])
863 {
864 extern Tree_T * DofTree_MH_moving ;
865
866 if (Dof->Type == DOF_UNKNOWN && !Tree_PQuery(DofTree_MH_moving, Dof))
867 Tree_Add(DofTree_MH_moving,Dof) ;
868 else if (Dof->Type == DOF_LINK && !Tree_PQuery(DofTree_MH_moving, Dof->Case.Link.Dof))
869 Tree_Add(DofTree_MH_moving,Dof->Case.Link.Dof) ;
870
871 if (Equ->Type == DOF_UNKNOWN && !Tree_PQuery(DofTree_MH_moving, Equ))
872 Tree_Add(DofTree_MH_moving,Equ) ;
873 else if (Equ->Type == DOF_LINK && !Tree_PQuery(DofTree_MH_moving, Equ->Case.Link.Dof))
874 Tree_Add(DofTree_MH_moving,Equ->Case.Link.Dof) ;
875
876 }
877 */
878