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 <math.h>
11 #include "GetDPConfig.h"
12 #include "ProData.h"
13 #include "ProDefine.h"
14 #include "F.h"
15 #include "MallocUtils.h"
16 #include "Message.h"
17 #include "Cal_Quantity.h"
18
19 #if defined(HAVE_KERNEL)
20 #include "DofData.h"
21 #include "Get_Geometry.h"
22 #include "Get_FunctionValue.h"
23 #endif
24
25 #define TWO_PI 6.2831853071795865
26
27 extern struct Problem Problem_S ;
28 extern struct CurrentData Current ;
29
30 #if defined(HAVE_KERNEL)
31
32 struct MH_InitData{
33 int Case ;
34 int NbrPoints, NbrPointsX; /* number of samples per smallest and fundamental period resp. */
35 struct DofData *DofData ;
36 double **H, ***HH ;
37 double *t, *w;
38 };
39
40 List_T * MH_InitData_L = NULL;
41
fcmp_MH_InitData(const void * a,const void * b)42 int fcmp_MH_InitData(const void * a, const void * b)
43 {
44 int Result ;
45 if ((Result = ((struct MH_InitData *)a)->DofData - ((struct MH_InitData *)b)->DofData) != 0)
46 return Result ;
47 if ((Result = ((struct MH_InitData *)a)->Case - ((struct MH_InitData *)b)->Case) != 0)
48 return Result ;
49 if (((struct MH_InitData *)a)->Case != 3)
50 return ((struct MH_InitData *)a)->NbrPoints - ((struct MH_InitData *)b)->NbrPoints ;
51 else
52 return ((struct MH_InitData *)a)->NbrPointsX - ((struct MH_InitData *)b)->NbrPointsX ;
53 }
54
55 #endif
56
NbrValues_Type(int Type)57 int NbrValues_Type (int Type)
58 {
59 switch (Type){
60 case SCALAR : return 1 ;
61 case VECTOR : case TENSOR_DIAG : return 3 ;
62 case TENSOR_SYM : return 6 ;
63 case TENSOR : return 9 ;
64 default :
65 Message::Error("Unknown type in NbrValues_Type");
66 return 0;
67 }
68 }
69
Product_SCALARxSCALARxSCALAR(double * V1,double * V2,double * V3)70 double Product_SCALARxSCALARxSCALAR (double *V1, double *V2, double *V3)
71 {
72 return V1[0] * V2[0] * V3[0] ;
73 }
74
Product_VECTORxTENSOR_SYMxVECTOR(double * V1,double * V2,double * V3)75 double Product_VECTORxTENSOR_SYMxVECTOR (double *V1, double *V2, double *V3)
76 {
77 return V3[0] * (V1[0] * V2[0] + V1[1] * V2[1] + V1[2] * V2[2]) +
78 V3[1] * (V1[0] * V2[1] + V1[1] * V2[3] + V1[2] * V2[4]) +
79 V3[2] * (V1[0] * V2[2] + V1[1] * V2[4] + V1[2] * V2[5]) ;
80 }
81
Product_VECTORxTENSOR_DIAGxVECTOR(double * V1,double * V2,double * V3)82 double Product_VECTORxTENSOR_DIAGxVECTOR (double *V1, double *V2, double *V3)
83 {
84 return V1[0] * V2[0] * V3[0] + V1[1] * V2[1] * V3[1] + V1[2] * V2[2] * V3[2] ;
85 }
86
Product_VECTORxSCALARxVECTOR(double * V1,double * V2,double * V3)87 double Product_VECTORxSCALARxVECTOR (double *V1, double *V2, double *V3)
88 {
89 return V2[0] * (V1[0] * V3[0] + V1[1] * V3[1] + V1[2] * V3[2]) ;
90 }
91
Product_VECTORxTENSORxVECTOR(double * V1,double * V2,double * V3)92 double Product_VECTORxTENSORxVECTOR (double *V1, double *V2, double *V3)
93 {
94 return V3[0] * (V1[0] * V2[0] + V1[1] * V2[3] + V1[2] * V2[6]) +
95 V3[1] * (V1[0] * V2[1] + V1[1] * V2[4] + V1[2] * V2[7]) +
96 V3[2] * (V1[0] * V2[2] + V1[1] * V2[5] + V1[2] * V2[8]) ;
97 }
98
Get_RealProductFunction_Type1xType2xType1(int Type1,int Type2)99 void *Get_RealProductFunction_Type1xType2xType1 (int Type1, int Type2)
100 {
101 if (Type1 == SCALAR && Type2 == SCALAR) {
102 return (void *)Product_SCALARxSCALARxSCALAR;
103 }
104 else if (Type1 == VECTOR && Type2 == TENSOR_SYM) {
105 return (void *)Product_VECTORxTENSOR_SYMxVECTOR;
106 }
107 else if (Type1 == VECTOR && Type2 == TENSOR_DIAG) {
108 return (void *)Product_VECTORxTENSOR_DIAGxVECTOR;
109 }
110 else if (Type1 == VECTOR && Type2 == SCALAR) {
111 return (void *)Product_VECTORxSCALARxVECTOR;
112 }
113 else if (Type1 == VECTOR && Type2 == TENSOR) {
114 return (void *)Product_VECTORxTENSORxVECTOR;
115 }
116 else {
117 Message::Error("Not allowed types in Get_RealProductFunction_Type1xType2xType1");
118 return 0;
119 }
120 }
121
122 /* ------------------------------------------------------------------------ */
123 /* MH_Get_InitData */
124 /* ------------------------------------------------------------------------ */
125
126 /*
127 Case = 1 : MHTransform NbrPoints (samples per smallest period) is given,
128 NbrPointsX (samples per fundamental period) is derived
129 Case = 2 : MHBilinear NbrPoints given, NbrPointsX derived
130 Case = 3 : HarmonicToTime NbrPointsX given, NbrPoints derived
131 */
132
MH_Get_InitData(int Case,int NbrPoints,int * NbrPointsX_P,double *** H_P,double **** HH_P,double ** t_P,double ** w_P)133 void MH_Get_InitData(int Case, int NbrPoints, int *NbrPointsX_P,
134 double ***H_P, double ****HH_P, double **t_P, double **w_P)
135 {
136 #if !defined(HAVE_KERNEL)
137 Message::Error("MH_Get_InitData requires Kernal");
138 #else
139 int NbrHar, iPul, iTime, iHar, jHar, NbrPointsX ;
140 double *Val_Pulsation, MaxPuls, MinPuls ;
141 double **H, ***HH = 0, *t, *w ;
142 struct MH_InitData MH_InitData_S, *MH_InitData_P ;
143
144 MH_InitData_S.Case = Case;
145 MH_InitData_S.DofData = Current.DofData;
146 MH_InitData_S.NbrPoints = NbrPoints;
147 MH_InitData_S.NbrPointsX = NbrPointsX = *NbrPointsX_P;
148
149 if (MH_InitData_L == NULL)
150 MH_InitData_L = List_Create(1, 1, sizeof(struct MH_InitData)) ;
151
152 if ((MH_InitData_P = (struct MH_InitData *)
153 List_PQuery(MH_InitData_L, &MH_InitData_S, fcmp_MH_InitData))){
154 *H_P = MH_InitData_P->H;
155 *HH_P = MH_InitData_P->HH;
156 *t_P = MH_InitData_P->t;
157 *w_P = MH_InitData_P->w;
158 *NbrPointsX_P = MH_InitData_P->NbrPointsX;
159 return;
160 }
161
162 NbrHar = Current.NbrHar;
163
164 Val_Pulsation = Current.DofData->Val_Pulsation;
165 MaxPuls = 0. ; MinPuls = 1.e99 ;
166
167 for (iPul = 0 ; iPul < NbrHar/2 ; iPul++) {
168 if (Val_Pulsation[iPul]){
169 MinPuls = std::min(Val_Pulsation[iPul], MinPuls);
170 MaxPuls = std::max(Val_Pulsation[iPul], MaxPuls);
171 }
172 /*
173 // Ruth: Previous implementation
174 if (Val_Pulsation[iPul] && Val_Pulsation[iPul] < MinPuls)
175 MinPuls = Val_Pulsation[iPul] ;
176 if (Val_Pulsation[iPul] && Val_Pulsation[iPul] > MaxPuls)
177 MaxPuls = Val_Pulsation[iPul] ;
178 */
179 }
180
181 if (Case != 3)
182 NbrPointsX = (int)((MaxPuls/MinPuls*(double)NbrPoints));
183 else
184 NbrPoints = (int)((MinPuls/MaxPuls*(double)NbrPointsX));
185
186 if(Case==1)
187 Message::Info("MH_Get_InitData (MHTransform) => NbrHar = %d NbrPoints = %d|%d",
188 NbrHar, NbrPoints, NbrPointsX);
189 if(Case==2)
190 Message::Info("MH_Get_InitData (MHBilinear) => NbrHar = %d NbrPoints = %d|%d",
191 NbrHar, NbrPoints, NbrPointsX);
192 if(Case==3)
193 Message::Info("MH_Get_InitData (HarmonicToTime) => NbrHar = %d NbrPoints = %d|%d",
194 NbrHar, NbrPoints, NbrPointsX);
195
196
197 t = (double *)Malloc(sizeof(double)*NbrPointsX) ;
198
199 if (Case != 3)
200 for (iTime = 0 ; iTime < NbrPointsX ; iTime++)
201 t[iTime] = (double)iTime/(double)NbrPointsX/(MinPuls/TWO_PI) ;
202 else
203 for (iTime = 0 ; iTime < NbrPointsX ; iTime++)
204 t[iTime] = (double)iTime/((double)NbrPointsX-1.)/(MinPuls/TWO_PI) ;
205
206 w = (double *)Malloc(sizeof(double)*NbrHar) ;
207 for (iPul = 0 ; iPul < NbrHar/2 ; iPul++)
208 if (Val_Pulsation[iPul]){
209 w[2*iPul ] = 2. / (double)NbrPointsX ;
210 w[2*iPul+1] = 2. / (double)NbrPointsX ;
211 }
212 else{
213 w[2*iPul ] = 1. / (double)NbrPointsX ;
214 w[2*iPul+1] = 1. / (double)NbrPointsX ;
215 }
216
217 H = (double **)Malloc(sizeof(double *)*NbrPointsX) ;
218 for (iTime = 0 ; iTime < NbrPointsX ; iTime++){
219 H[iTime] = (double *)Malloc(sizeof(double)*NbrHar) ;
220 for (iPul = 0 ; iPul < NbrHar/2 ; iPul++) {
221 H[iTime][2*iPul ] = cos(Val_Pulsation[iPul] * t[iTime]) ;
222 H[iTime][2*iPul+1] = - sin(Val_Pulsation[iPul] * t[iTime]) ;
223 }
224 }
225
226 /*
227 for (iHar = 0 ; iHar < NbrHar ; iHar++)
228 for (jHar = iHar ; jHar < NbrHar ; jHar++){
229 sum = 0.;
230 for (iTime = 0 ; iTime < NbrPointsX ; iTime++)
231 sum += w[iTime] * H[iTime][iHar] * H[iTime][jHar] ;
232 sum -= (iHar==jHar)? 1. : 0. ;
233 printf("iHar %d jHar %d sum %e\n", iHar, jHar, sum);
234 }
235 */
236
237
238 if (Case == 2) {
239 //if(Current.DofData->Flag_Init[0] < 2)
240 // Message::Error("Jacobian system not initialized (missing GenerateJac?)");
241
242 HH = (double ***)Malloc(sizeof(double **)*NbrPointsX) ;
243 for (iTime = 0 ; iTime < NbrPointsX ; iTime++){
244 HH[iTime] = (double **)Malloc(sizeof(double *)*NbrHar) ;
245 for (iHar = 0 ; iHar < NbrHar ; iHar++){
246 HH[iTime][iHar] = (double *)Malloc(sizeof(double)*NbrHar) ;
247 for (jHar = 0 ; jHar < NbrHar ; jHar++){
248 if (Val_Pulsation [iHar/2] && Val_Pulsation [jHar/2] )
249 HH[iTime][iHar][jHar] = 2. / (double)NbrPointsX * H[iTime][iHar] * H[iTime][jHar] ;
250 else
251 HH[iTime][iHar][jHar] = 1. / (double)NbrPointsX * H[iTime][iHar] * H[iTime][jHar] ;
252 }
253 }
254 }
255 }
256
257 *H_P = MH_InitData_S.H = H;
258 *t_P = MH_InitData_S.t = t;
259 *w_P = MH_InitData_S.w = w;
260 *HH_P = MH_InitData_S.HH = HH;
261 *NbrPointsX_P = MH_InitData_S.NbrPointsX = NbrPointsX;
262 List_Add (MH_InitData_L, &MH_InitData_S);
263 #endif
264 }
265
266 /* ------------------------------------------------------------------------ */
267 /* F_MHToTime0 (HarmonicToTime in PostOperation) */
268 /* ------------------------------------------------------------------------ */
269
F_MHToTime0(int init,struct Value * A,struct Value * V,int iTime,int NbrPointsX,double * TimeMH)270 void F_MHToTime0(int init, struct Value * A, struct Value * V,
271 int iTime, int NbrPointsX, double * TimeMH)
272 {
273 static double **H, ***HH, *t, *weight;
274 int iVal, nVal, iHar;
275
276 if (Current.NbrHar == 1) return;
277
278 if (!init) MH_Get_InitData(3, 0, &NbrPointsX, &H, &HH, &t, &weight);
279
280 *TimeMH = t[iTime] ;
281 V->Type = A->Type ;
282 nVal = NbrValues_Type (A->Type) ;
283
284 for (iVal = 0 ; iVal < nVal ; iVal++){
285 V->Val[iVal] = 0;
286 for (iHar = 0 ; iHar < Current.NbrHar ; iHar++)
287 V->Val[iVal] += H[iTime][iHar] * A->Val[iHar*MAX_DIM+iVal] ;
288 }
289
290 }
291
292
293
294 /* ---------------------------------------------------------------------- */
295 /* F_MHToTime */
296 /* ---------------------------------------------------------------------- */
297
F_MHToTime(struct Function * Fct,struct Value * A,struct Value * V)298 void F_MHToTime (struct Function * Fct, struct Value * A, struct Value * V)
299 {
300 #if !defined(HAVE_KERNEL)
301 Message::Error("F_MHToTime requires Kernel");
302 #else
303 int iHar, iVal, nVal ;
304 double time, H[NBR_MAX_HARMONIC];
305 struct Value Vtemp;
306
307 if (Current.NbrHar == 1) Message::Error("'F_MHToTime' only for Multi-Harmonic stuff") ;
308 if((A+1)->Type != SCALAR)
309 Message::Error("'F_MHToTime' requires second scalar argument (time)");
310 time = (A+1)->Val[0] ;
311
312 for (iHar = 0 ; iHar < Current.NbrHar/2 ; iHar++) {
313 /* if (Current.DofData->Val_Pulsation [iHar]){ */
314 H[2*iHar ] = cos(Current.DofData->Val_Pulsation[iHar] * time) ;
315 H[2*iHar+1] = - sin(Current.DofData->Val_Pulsation[iHar] * time) ;
316 /*
317 }
318 else {
319 H[2*iHar ] = 0.5 ;
320 H[2*iHar+1] = 0 ;
321 }
322 */
323 }
324
325 nVal = NbrValues_Type (A->Type) ;
326
327 for (iVal = 0 ; iVal < MAX_DIM ; iVal++)
328 for (iHar = 0 ; iHar < Current.NbrHar ; iHar++)
329 Vtemp.Val[iHar*MAX_DIM+iVal] = 0.;
330
331 for (iVal = 0 ; iVal < nVal ; iVal++)
332 for (iHar = 0 ; iHar < Current.NbrHar ; iHar++)
333 Vtemp.Val[iVal] += H[iHar] * A->Val[iHar*MAX_DIM+iVal] ;
334
335 V->Type = A->Type ;
336 for (iVal = 0 ; iVal < MAX_DIM ; iVal++)
337 for (iHar = 0 ; iHar < Current.NbrHar ; iHar++)
338 V->Val[iHar*MAX_DIM+iVal] = Vtemp.Val[iHar*MAX_DIM+iVal] ;
339 #endif
340 }
341
342
343 /* ------------------------------------------------------------------------ */
344 /* MHTransform */
345 /* ------------------------------------------------------------------------ */
346
MHTransform(struct Element * Element,struct QuantityStorage * QuantityStorage_P0,double u,double v,double w,std::vector<struct Value> & MH_Inputs,struct Expression * Expression_P,int NbrPoints,struct Value & MH_Output)347 void MHTransform(struct Element * Element, struct QuantityStorage * QuantityStorage_P0,
348 double u, double v, double w, std::vector<struct Value> &MH_Inputs,
349 struct Expression * Expression_P, int NbrPoints, struct Value &MH_Output)
350 {
351 double **H, ***HH, *t, *weight ;
352 int NbrPointsX;
353 MH_Get_InitData(1, NbrPoints, &NbrPointsX, &H, &HH, &t, &weight);
354
355 int NbrHar = Current.NbrHar; // save NbrHar
356 Current.NbrHar = 1; // evaluation in time domain!
357
358 for (int iVal = 0 ; iVal < MAX_DIM ; iVal++)
359 for (int iHar = 0 ; iHar < NbrHar ; iHar++)
360 MH_Output.Val[iHar*MAX_DIM+iVal] = 0. ;
361
362 int N = MH_Inputs.size(), nVal2 = 0;
363 std::vector<struct Value> t_Values(N + 1); // in case N==0
364
365 for (int iTime = 0 ; iTime < NbrPointsX ; iTime++) {
366 // evaluate MH_Inputs at iTime-th time point
367 for(int j = 0; j < N; j++){
368 int nVal1 = NbrValues_Type(MH_Inputs[j].Type);
369 t_Values[j].Type = MH_Inputs[j].Type;
370 for (int iVal = 0 ; iVal < nVal1 ; iVal++){
371 t_Values[j].Val[iVal] = 0.;
372 for (int iHar = 0 ; iHar < NbrHar ; iHar++)
373 t_Values[j].Val[iVal] += H[iTime][iHar] * MH_Inputs[j].Val[iHar*MAX_DIM+iVal] ;
374 }
375 }
376
377 // evaluate the function, passing the N time-domain values as arguments
378 Get_ValueOfExpression(Expression_P, QuantityStorage_P0, u, v, w, &t_Values[0], N);
379
380 if(!iTime){
381 nVal2 = NbrValues_Type(t_Values[0].Type) ;
382 MH_Output.Type = t_Values[0].Type;
383 }
384 for (int iVal = 0 ; iVal < nVal2 ; iVal++)
385 for (int iHar = 0 ; iHar < NbrHar ; iHar++)
386 MH_Output.Val[iHar*MAX_DIM+iVal] +=
387 weight[iHar] * H[iTime][iHar] * t_Values[0].Val[iVal] ;
388 }
389
390 Current.NbrHar = NbrHar ; // reset NbrHar
391 }
392
393 #if defined(HAVE_KERNEL)
394
395 /* ----------------------------------------------------------------------------------- */
396 /* C a l _ I n i t G a l e r k i n T e r m O f F e m E q u a t i o n _ M H J a c N L */
397 /* ----------------------------------------------------------------------------------- */
398
Cal_InitGalerkinTermOfFemEquation_MHBilinear(struct EquationTerm * EquationTerm_P)399 void Cal_InitGalerkinTermOfFemEquation_MHBilinear(struct EquationTerm * EquationTerm_P)
400 {
401 struct FemLocalTermActive * FI ;
402 List_T * WholeQuantity_L;
403 struct WholeQuantity *WholeQuantity_P0 ;
404 int i_WQ ;
405
406 FI = EquationTerm_P->Case.LocalTerm.Active ;
407 FI->MHBilinear = 0 ;
408
409 /* search for MHBilinear-term(s) */
410 WholeQuantity_L = EquationTerm_P->Case.LocalTerm.Term.WholeQuantity ;
411 WholeQuantity_P0 = (struct WholeQuantity*)List_Pointer(WholeQuantity_L, 0) ;
412 i_WQ = 0 ; while ( i_WQ < List_Nbr(WholeQuantity_L) &&
413 (WholeQuantity_P0 + i_WQ)->Type != WQ_MHBILINEAR) i_WQ++ ;
414
415 if (i_WQ == List_Nbr(WholeQuantity_L) )
416 return; /* no MHBilinear stuff, let's get the hell out of here ! */
417
418 /* check if Galerkin term produces symmetrical contribution to system matrix */
419 if (!FI->SymmetricalMatrix)
420 Message::Error("Galerkin term with MHBilinear must be symmetrical");
421
422 if(EquationTerm_P->Case.LocalTerm.Term.CanonicalWholeQuantity_Equ != CWQ_NONE)
423 Message::Error("Not allowed expression in Galerkin term with MHBilinear");
424
425 if (List_Nbr(WholeQuantity_L) == 3){
426 if (i_WQ != 0 ||
427 EquationTerm_P->Case.LocalTerm.Term.DofIndexInWholeQuantity != 1 ||
428 (WholeQuantity_P0 + 2)->Type != WQ_BINARYOPERATOR ||
429 (WholeQuantity_P0 + 2)->Case.Operator.TypeOperator != OP_TIME)
430 Message::Error("Not allowed expression in Galerkin term with MHBilinear");
431 }
432 else {
433 Message::Error("Not allowed expression in Galerkin term with MHBilinear (%d terms) ",
434 List_Nbr(WholeQuantity_L));
435 }
436
437 FI->MHBilinear = 1 ;
438 // index of function, e.g. dhdb[{d a}]
439 FI->MHBilinear_Index = (WholeQuantity_P0 + i_WQ)->Case.MHBilinear.Index ;
440 // list of quantities to transform (arguments of the function)
441 FI->MHBilinear_WholeQuantity_L = (WholeQuantity_P0 + i_WQ)->Case.MHBilinear.WholeQuantity_L ;
442 if(Message::GetVerbosity() == 10)
443 Message::Info("FreqOffSet in 'MHBilinear' == %d ", (WholeQuantity_P0 + i_WQ)->Case.MHBilinear.FreqOffSet) ;
444 FI->MHBilinear_HarOffSet = 2 * (WholeQuantity_P0 + i_WQ)->Case.MHBilinear.FreqOffSet ;
445 if (FI->MHBilinear_HarOffSet > Current.NbrHar-2){
446 Message::Warning("FreqOffSet in 'MHBilinear' cannot exceed %d => set to %d ",
447 Current.NbrHar/2-1, Current.NbrHar/2-1) ;
448 FI->MHBilinear_HarOffSet = Current.NbrHar-2 ;
449 }
450 FI->MHBilinear_JacNL = (EquationTerm_P->Case.LocalTerm.Term.TypeTimeDerivative == JACNL_);
451
452 MH_Get_InitData(2, (WholeQuantity_P0 + i_WQ)->Case.MHBilinear.NbrPoints,
453 &FI->MHBilinear_NbrPointsX, &FI->MHBilinear_H, &FI->MHBilinear_HH,
454 &FI->MHBilinear_t, &FI->MHBilinear_w) ;
455 }
456
457 #define OFFSET (iHar < NbrHar-OffSet)? 0 : iHar-NbrHar+OffSet+2-iHar%2
458
459 /* --------------------------------------------------------------------------- */
460 /* C a l _ G a l e r k i n T e r m O f F e m E q u a t i o n _ M H J a c N L */
461 /* --------------------------------------------------------------------------- */
462
Cal_GalerkinTermOfFemEquation_MHBilinear(struct Element * Element,struct EquationTerm * EquationTerm_P,struct QuantityStorage * QuantityStorage_P0)463 void Cal_GalerkinTermOfFemEquation_MHBilinear(struct Element * Element,
464 struct EquationTerm * EquationTerm_P,
465 struct QuantityStorage * QuantityStorage_P0)
466 {
467 struct FemLocalTermActive * FI ;
468 struct QuantityStorage * QuantityStorage_P;
469 struct IntegrationCase * IntegrationCase_P ;
470 struct Quadrature * Quadrature_P ;
471
472 int Nbr_Dof, NbrHar ;
473 double vBFuDof [NBR_MAX_BASISFUNCTIONS][MAX_DIM] ;
474 double vBFxDof [NBR_MAX_BASISFUNCTIONS][MAX_DIM] ;
475 double Val_Dof [NBR_MAX_BASISFUNCTIONS][NBR_MAX_HARMONIC] ;
476 double Val [3*NBR_MAX_BASISFUNCTIONS];
477
478 int i, j, k, Type_Dimension, Nbr_IntPoints, i_IntPoint ;
479 int iTime, iDof, jDof, iHar, jHar, nVal1, nVal2 = 0, iVal1, iVal2, Type1;
480 double **H, ***HH, plus, plus0, weightIntPoint;
481 int NbrPointsX, OffSet;
482 struct Expression * Expression_P;
483 struct Dof * Dofi, *Dofj;
484
485 // double one = 1.0 ;
486 int iPul, ZeroHarmonic, DcHarmonic;
487 double E_D[NBR_MAX_HARMONIC][NBR_MAX_HARMONIC][MAX_DIM];
488
489 void (*xFunctionBFDof[NBR_MAX_BASISFUNCTIONS])
490 (struct Element * Element, int NumEntity,
491 double u, double v, double w, double Value[] ) ;
492 double (*Get_Jacobian)(struct Element*, MATRIX3x3*) ;
493 void (*Get_IntPoint)(int,int,double*,double*,double*,double*);
494 double (*Get_Product)(double*,double*,double*) = 0;
495
496 FI = EquationTerm_P->Case.LocalTerm.Active ;
497 QuantityStorage_P = FI->QuantityStorageDof_P ;
498
499 /* ---------------------------------------------------------------------- */
500 /* G e t F u n c t i o n V a l u e f o r t e s t f u n c t i o n s */
501 /* ---------------------------------------------------------------------- */
502
503 if (!(Nbr_Dof = QuantityStorage_P->NbrElementaryBasisFunction)){
504 return;
505 }
506
507 // test!
508 std::vector<std::vector<std::vector<std::vector<double> > > > E_MH(Nbr_Dof);
509 for(unsigned int i = 0; i < E_MH.size(); i++){
510 E_MH[i].resize(Nbr_Dof);
511 for(unsigned int j = 0; j < E_MH[i].size(); j++){
512 E_MH[i][j].resize(Current.NbrHar);
513 for(unsigned int k = 0; k < E_MH[i][j].size(); k++){
514 E_MH[i][j][k].resize(Current.NbrHar, 0.);
515 }
516 }
517 }
518
519 Get_FunctionValue(Nbr_Dof, (void (**)())xFunctionBFDof,
520 EquationTerm_P->Case.LocalTerm.Term.TypeOperatorDof,
521 QuantityStorage_P, &FI->Type_FormDof) ;
522
523
524 for (j = 0 ; j < Nbr_Dof ; j++)
525 for (k = 0 ; k < Current.NbrHar ; k+=2)
526 Dof_GetComplexDofValue
527 (QuantityStorage_P->FunctionSpace->DofData,
528 QuantityStorage_P->BasisFunction[j].Dof + k/2*gCOMPLEX_INCREMENT,
529 &Val_Dof[j][k], &Val_Dof[j][k+1]) ;
530
531 nVal1 = NbrValues_Type (Type1 = Get_ValueFromForm(FI->Type_FormDof)) ;
532
533 /* ------------------------------------------------------------------- */
534 /* G e t J a c o b i a n M e t h o d */
535 /* ------------------------------------------------------------------- */
536
537 i = 0 ; while ((i < FI->NbrJacobianCase) &&
538 ((j = (FI->JacobianCase_P0 + i)->RegionIndex) >= 0) &&
539 !List_Search
540 (((struct Group *)List_Pointer(Problem_S.Group, j))
541 ->InitialList, &Element->Region, fcmp_int) ) i++ ;
542
543 if (i == FI->NbrJacobianCase)
544 Message::Error("Undefined Jacobian in Region %d", Element->Region);
545
546 Element->JacobianCase = FI->JacobianCase_P0 + i ;
547
548 Get_Jacobian = (double (*)(struct Element*, MATRIX3x3*))
549 Get_JacobianFunction(Element->JacobianCase->TypeJacobian,
550 Element->Type, &Type_Dimension) ;
551
552 if (FI->Flag_ChangeCoord) Get_NodesCoordinatesOfElement(Element) ;
553
554 /* integration in space */
555
556 IntegrationCase_P =
557 Get_IntegrationCase(Element, FI->IntegrationCase_L, FI->CriterionIndex);
558
559 switch (IntegrationCase_P->Type) {
560 case ANALYTIC :
561 Message::Error("Analytical integration not implemented for MHBilinear");
562 }
563
564 Quadrature_P = (struct Quadrature*)
565 List_PQuery(IntegrationCase_P->Case, &Element->Type, fcmp_int);
566
567 if(!Quadrature_P)
568 Message::Error("Unknown type of Element (%s) for Integration method (%s)",
569 Get_StringForDefine(Element_Type, Element->Type),
570 ((struct IntegrationMethod *)
571 List_Pointer(Problem_S.IntegrationMethod,
572 EquationTerm_P->Case.LocalTerm.IntegrationMethodIndex))->Name);
573
574 Nbr_IntPoints = Quadrature_P->NumberOfPoints ;
575 Get_IntPoint = (void(*)(int,int,double*,double*,double*,double*))
576 Quadrature_P->Function ;
577
578 /* integration in fundamental time period */
579
580 NbrPointsX = FI->MHBilinear_NbrPointsX;
581 HH = FI->MHBilinear_HH;
582 H = FI->MHBilinear_H ;
583 Expression_P = (struct Expression*)List_Pointer(Problem_S.Expression, FI->MHBilinear_Index);
584 OffSet = FI->MHBilinear_HarOffSet;
585
586 /* ------------------------------------------------------------------------ */
587 /* ------------------------------------------------------------------------ */
588 /* C o m p u t a t i o n o f E l e m e n t a r y m a t r i x */
589 /* ------------------------------------------------------------------------ */
590 /* ------------------------------------------------------------------------ */
591
592 NbrHar = Current.NbrHar ;
593
594 /* special treatment of DC-term and associated dummy sinus-term */
595 DcHarmonic = NbrHar;
596 ZeroHarmonic = 0;
597 for (iPul = 0 ; iPul < NbrHar/2 ; iPul++)
598 if (!Current.DofData->Val_Pulsation[iPul]){
599 DcHarmonic = 2*iPul ;
600 ZeroHarmonic = 2*iPul+1 ;
601 break;
602 }
603
604 /* volume integration over element */
605 for (i_IntPoint = 0 ; i_IntPoint < Nbr_IntPoints ; i_IntPoint++) {
606
607 Get_IntPoint(Nbr_IntPoints, i_IntPoint,
608 &Current.u, &Current.v, &Current.w, &weightIntPoint) ;
609
610 if (FI->Flag_ChangeCoord) {
611 Get_BFGeoElement(Element, Current.u, Current.v, Current.w) ;
612
613 Element->DetJac = Get_Jacobian(Element, &Element->Jac) ;
614 weightIntPoint *= fabs(Element->DetJac) ;
615
616 if (FI->Flag_InvJac)
617 Get_InverseMatrix(Type_Dimension, Element->Type, Element->DetJac,
618 &Element->Jac, &Element->InvJac) ;
619 }
620
621 // Test and shape Functions (are the same)
622 for (i = 0 ; i < Nbr_Dof ; i++) {
623 xFunctionBFDof[i]
624 (Element,
625 QuantityStorage_P->BasisFunction[i].NumEntityInElement+1,
626 Current.u, Current.v, Current.w, vBFuDof[i]) ;
627 ((void (*)(struct Element*, double*, double*))
628 FI->xChangeOfCoordinatesEqu) (Element, vBFuDof[i], vBFxDof[i]) ;
629 }
630
631 switch (Type1) {
632 case SCALAR :
633 for (k = 0 ; k < NbrHar ; k++){
634 Val[k] = 0.;
635 for (j = 0 ; j < Nbr_Dof ; j++)
636 Val[k] += vBFxDof[j][0] * Val_Dof[j][k] ;
637 }
638 break ;
639 case VECTOR :
640 for (k = 0 ; k < NbrHar ; k++){
641 Val[3*k] = Val[3*k+1] = Val[3*k+2] = 0.;
642 for (j = 0 ; j < Nbr_Dof ; j++){
643 Val[3*k ] += vBFxDof[j][0] * Val_Dof[j][k] ;
644 Val[3*k+1] += vBFxDof[j][1] * Val_Dof[j][k] ;
645 Val[3*k+2] += vBFxDof[j][2] * Val_Dof[j][k] ;
646 }
647 }
648 break ;
649 }
650
651 // evaluate additional (multi-harmonic) arguments
652 int N = List_Nbr(FI->MHBilinear_WholeQuantity_L);
653
654 std::vector<struct Value> MH_Values(N);
655 for(int j = 1; j < N; j++){
656 List_T *WQ; List_Read(FI->MHBilinear_WholeQuantity_L, j, &WQ);
657 Cal_WholeQuantity(Element, QuantityStorage_P0, WQ, Current.u, Current.v, Current.w, -1, 0,
658 &MH_Values[j]) ;
659 }
660
661 Current.NbrHar = 1; /* evaluation in time domain */
662
663 int saveTime = Current.Time;
664 int saveTimeStep = Current.TimeStep;
665
666 std::vector<struct Value> t_Values(N + 1); // in case N==0
667
668 // time integration over fundamental period
669 for (iTime = 0 ; iTime < NbrPointsX ; iTime++) {
670
671 Current.TimeStep = iTime;
672 Current.Time = iTime * 2. * M_PI / (double)NbrPointsX;
673
674 t_Values[0].Type = Type1 ;
675 for (iVal1 = 0 ; iVal1 < nVal1 ; iVal1++){
676 t_Values[0].Val[iVal1] = 0;
677 for (iHar = 0 ; iHar < NbrHar ; iHar++)
678 t_Values[0].Val[iVal1] += H[iTime][iHar] * Val[iHar*nVal1+iVal1] ;
679 }
680 for(int j = 1; j < N; j++){
681 int nVal1 = NbrValues_Type(MH_Values[j].Type);
682 t_Values[j].Type = MH_Values[j].Type;
683 for (int iVal = 0 ; iVal < nVal1 ; iVal++){
684 t_Values[j].Val[iVal] = 0.;
685 for (int iHar = 0 ; iHar < NbrHar ; iHar++)
686 t_Values[j].Val[iVal] += H[iTime][iHar] * MH_Values[j].Val[iHar*MAX_DIM+iVal] ;
687 }
688 }
689
690 Get_ValueOfExpression(Expression_P, QuantityStorage_P0,
691 Current.u, Current.v, Current.w, &t_Values[0], N);
692
693 if (!iTime){
694 if (!i_IntPoint){
695 nVal2 = NbrValues_Type (t_Values[0].Type) ;
696 Get_Product = (double(*)(double*,double*,double*))
697 Get_RealProductFunction_Type1xType2xType1 (Type1, t_Values[0].Type) ;
698 }
699 for (iHar = 0 ; iHar < NbrHar ; iHar++)
700 for (jHar = OFFSET ; jHar <= iHar ; jHar++)
701 for (iVal2 = 0 ; iVal2 < nVal2 ; iVal2++)
702 E_D[iHar][jHar][iVal2] = 0. ;
703 }
704
705 for (iHar = 0 ; iHar < NbrHar ; iHar++)
706 for (jHar = OFFSET ; jHar <= iHar ; jHar++){
707 for (iVal2 = 0 ; iVal2 < nVal2 ; iVal2++)
708 E_D[iHar][jHar][iVal2] += HH[iTime][iHar][jHar] * t_Values[0].Val[iVal2] ;
709 }
710
711 } /* for iTime ... */
712
713 Current.TimeStep = saveTimeStep;
714 Current.Time = saveTime;
715
716 for (iDof = 0 ; iDof < Nbr_Dof ; iDof++)
717 for (jDof = 0 ; jDof <= iDof ; jDof++)
718 for (iHar = 0 ; iHar < NbrHar ; iHar++)
719 for (jHar = OFFSET ; jHar <= iHar ; jHar++){
720 E_MH[iDof][jDof][iHar][jHar] += weightIntPoint *
721 Get_Product(vBFxDof[iDof], E_D[iHar][jHar], vBFxDof[jDof]) ;
722 Message::Debug("E_MH[%d][%d][%d][%d] = %e",
723 iDof, jDof, iHar, jHar, E_MH[iDof][jDof][iHar][jHar]) ;
724 }
725
726 Current.NbrHar = NbrHar ;
727
728 } /* for i_IntPoint ... */
729
730 /* set imaginary part = to real part for 0th harmonic; this replaces the dummy
731 regularization that we did before (see below, "dummy 1's on the diagonal...") */
732 if(ZeroHarmonic){
733 for (iDof = 0 ; iDof < Nbr_Dof ; iDof++){
734 for (jDof = 0 ; jDof < Nbr_Dof ; jDof++){
735 E_MH[iDof][jDof][1][1] = E_MH[iDof][jDof][0][0];
736 }
737 }
738 }
739
740 /* -------------------------------------------------------------------- */
741 /* A d d c o n t r i b u t i o n t o J a c o b i a n M a t r i x */
742 /* -------------------------------------------------------------------- */
743
744 gMatrix *matrix;
745 gVector *rhs = NULL;
746 if(FI->MHBilinear_JacNL){
747 matrix = &Current.DofData->Jac;
748 }
749 else{
750 matrix = &Current.DofData->A;
751 rhs = &Current.DofData->b;
752 }
753
754 for (iDof = 0 ; iDof < Nbr_Dof ; iDof++){
755 Dofi = QuantityStorage_P->BasisFunction[iDof].Dof ;
756 for (jDof = 0 ; jDof <= iDof ; jDof++){
757 Dofj = QuantityStorage_P->BasisFunction[jDof].Dof ;
758
759 for (iHar = 0 ; iHar < NbrHar ; iHar++)
760 for (jHar = OFFSET ; jHar <= iHar ; jHar++){
761 plus = plus0 = E_MH[iDof][jDof][iHar][jHar] ;
762
763 if(jHar==DcHarmonic && iHar!=DcHarmonic) { plus0 *= 1. ; plus *= 2. ;}
764 // FIXME: this cannot work in complex arithmetic: AssembleInMat will
765 // need to assemble both real and imaginary parts at once -- See
766 // Cal_AssembleTerm for an example
767 Dof_AssembleInMat(Dofi+iHar, Dofj+jHar, 1, &plus, matrix, rhs) ;
768 if(iHar != jHar)
769 Dof_AssembleInMat(Dofi+jHar, Dofj+iHar, 1, &plus0, matrix, rhs) ;
770 if(iDof != jDof){
771 Dof_AssembleInMat(Dofj+iHar, Dofi+jHar, 1, &plus, matrix, rhs) ;
772 if(iHar != jHar)
773 Dof_AssembleInMat(Dofj+jHar, Dofi+iHar, 1, &plus0, matrix, rhs) ;
774 }
775 }
776 }
777 }
778
779 /* dummy 1's on the diagonal for sinus-term of dc-component */
780 /*
781 if (ZeroHarmonic) {
782 for (iDof = 0 ; iDof < Nbr_Dof ; iDof++){
783 Dofi = QuantityStorage_P->BasisFunction[iDof].Dof + ZeroHarmonic ;
784 // FIXME: this cannot work in complex arithmetic: AssembleInMat will
785 // need to assemble both real and imaginary parts at once -- See
786 // Cal_AssembleTerm for an example
787 Dof_AssembleInMat(Dofi, Dofi, 1, &one, matrix, rhs) ;
788 }
789 }
790 */
791
792 }
793
794 #undef OFFSET
795
796 #endif
797