ParMooN
 All Classes Functions Variables Friends Pages
Benchmark_DEVSS.h
1 // Navier-Stokes coupled with conformation stress tensor problem, Benchmark problem
2 
3 
4 void ExampleFile()
5 {
6  OutPut("Example: Benchmark_FlowPastCylinder_DEVSS" << endl) ;
7 }
8 
9 #define __BENCH__
10 
11 #include <Joint.h>
12 #include <BoundEdge.h>
13 #include <BoundComp.h>
14 #include <FE2D.h>
15 #include <FEDesc2D.h>
16 // ========================================================================
17 // exact solution
18 // ========================================================================
19 void ExactU1(double x, double y, double *values)
20 {
21  values[0] = 0;
22  values[1] = 0;
23  values[2] = 0;
24  values[3] = 0;
25 }
26 
27 void ExactU2(double x, double y, double *values)
28 {
29  values[0] = 0;
30  values[1] = 0;
31  values[2] = 0;
32  values[3] = 0;
33 }
34 
35 void ExactP(double x, double y, double *values)
36 {
37  values[0] = 0;
38  values[1] = 0;
39  values[2] = 0;
40  values[3] = 0;
41 }
42 
43 // ========================================================================
44 // boundary conditions
45 // ========================================================================
46 void BoundCondition(int i, double t, BoundCond &cond)
47 {
48  cond = DIRICHLET;
49 }
50 
51 void U1BoundValue(int BdComp, double Param, double &value)
52 {
53  switch(BdComp)
54  {
55  case 0: value = 0;
56  break;
57  case 1: value=1.2*Param*(1.0-Param);
58  break;
59  case 2: value = 0;
60  break;
61  case 3: value=1.2*Param*(1.0-Param) ;
62  break;
63  case 4: value=0;
64  break;
65  default: cout << "wrong boundary part number: " << BdComp << endl;
66  }
67 }
68 
69 void U2BoundValue(int BdComp, double Param, double &value)
70 {
71  value = 0;
72  if(BdComp>4) cout << "wrong boundary part number: " << BdComp << endl;
73 }
74 
75 // ========================================================================
76 // coefficients for Stokes form: A, B1, B2, f1, f2
77 // ========================================================================
78 void LinCoeffs(int n_points, double *x, double *y,
79  double **parameters, double **coeffs)
80 {
81  double eps = 1/TDatabase::ParamDB->RE_NR, beta = TDatabase::ParamDB->P3;
82  int i;
83  double *coeff, nondim;
84 
85  if (TDatabase::ParamDB->TENSOR_TYPE == 1)
86  nondim = beta*eps;
87  else if (TDatabase::ParamDB->TENSOR_TYPE == 2)
88  nondim = eps;
89 
90  for(i=0;i<n_points;i++)
91  {
92  coeff = coeffs[i];
93 
94  coeff[0] = eps;
95  coeff[1] = 0; // f1
96  coeff[2] = 0; // f2
97  }
98 }
99 
100 // ========================================================================
101 // exact solution
102 // ========================================================================
103 void ExactS1(double x, double y, double *values)
104 {
105  values[0] = 0;
106  values[1] = 0;
107  values[2] = 0;
108  values[3] = 0;
109 }
110 
111 void ExactS2(double x, double y, double *values)
112 {
113  values[0] = 0;
114  values[1] = 0;
115  values[2] = 0;
116  values[3] = 0;
117 }
118 
119 void ExactS3(double x, double y, double *values)
120 {
121  values[0] = 0;
122  values[1] = 0;
123  values[2] = 0;
124  values[3] = 0;
125 }
126 
127 // ========================================================================
128 // boundary conditions
129 // ========================================================================
130 void BoundCondition_CST(int i, double t, BoundCond &cond)
131 {
132  cond = DIRICHLET;
133 
134 }
135 
136 void S1BoundValue(int BdComp, double Param, double &value)
137 {
138  double Wei = TDatabase::ParamDB->WEI_NR;
139  switch(BdComp)
140  {
141  case 0: value=1.0;
142  break;
143  case 1: value= (pow(Wei*1.2*(1.0-(2.0*Param))/0.41,2) * 2.0 )+1.0;
144  break;
145  case 2: value=1.0;
146  break;
147  case 3: value=(pow(Wei*1.2*(1.0-(2.0*(1.0-Param)))/0.41,2)* 2.0 )+1.0;
148  break;
149  case 4: value=1.0;
150  break;
151  default: cout << "wrong boundary part number" << endl;
152  break;
153  }
154 }
155 
156 void S2BoundValue(int BdComp, double Param, double &value)
157 {
158  double Wei = TDatabase::ParamDB->WEI_NR;
159  switch(BdComp)
160  {
161  case 0: value=0;
162  break;
163  case 1: value=Wei*1.2*(1.0-(2.0*Param))/0.41;
164  break;
165  case 2: value=0;
166  break;
167  case 3: value=Wei*1.2*(1.0-(2.0*(1.0-Param)))/0.41;
168  break;
169  case 4: value=0;
170  break;
171  default: cout << "wrong boundary part number" << endl;
172  break;
173  }
174  return;
175 }
176 
177 
178 
179 void S3BoundValue(int BdComp, double Param, double &value)
180 {
181  value=1.0;
182  if(BdComp>4) cout << "wrong boundary part number: " << BdComp << endl;
183 }
184 
185 void LinCoeffs_CST(int n_points, double *X, double *Y,
186  double **parameters, double **coeffs)
187 {
188  double nu=1./TDatabase::ParamDB->WEI_NR;
189 
190  int i;
191  double *coeff;
192 
193  for(i=0;i<n_points;i++)
194  {
195  coeff = coeffs[i];
196 
197  coeff[0] = nu;
198  coeff[1] = nu; //f1
199  coeff[2] = 0.0; //f2
200  coeff[3] = nu; //f3
201 
202 
203  }
204 }
205 
206 
207 // ========================================================================
208 // exact solution
209 // ========================================================================
210 void ExactD1(double x, double y, double *values)
211 {
212  values[0] = 0;
213  values[1] = 0;
214  values[2] = 0;
215  values[3] = 0;
216 }
217 
218 void ExactD2(double x, double y, double *values)
219 {
220  values[0] = 0;
221  values[1] = 0;
222  values[2] = 0;
223  values[3] = 0;
224 }
225 
226 void ExactD3(double x, double y, double *values)
227 {
228  values[0] = 0;
229  values[1] = 0;
230  values[2] = 0;
231  values[3] = 0;
232 }
233 
234 // ========================================================================
235 // boundary conditions
236 // ========================================================================
237 void BoundCondition_DFT(int i, double t, BoundCond &cond)
238 {
239  cond = DIRICHLET;
240 
241 }
242 
243 void D1BoundValue(int BdComp, double Param, double &value)
244 {
245  switch(BdComp)
246  {
247  case 0: value=0;
248  break;
249  case 1: value=0;
250  break;
251  case 2: value=0;
252  break;
253  case 3: value=0;
254  break;
255  case 4: value=0;
256  break;
257  default: cout << "wrong boundary part number" << endl;
258  break;
259  }
260  return;
261 }
262 
263 void D2BoundValue(int BdComp, double Param, double &value)
264 {
265  switch(BdComp)
266  {
267  case 0: value=0;
268  break;
269  case 1: value=0.6*(1-2*Param)/0.41;
270  break;
271  case 2: value=0;
272  break;
273  case 3: value=0.6*(1-2*(1-Param))/0.41;
274  break;
275  case 4: value=0;
276  break;
277  default: cout << "wrong boundary part number" << endl;
278  break;
279  }
280  return;
281 }
282 
283 
284 
285 void D3BoundValue(int BdComp, double Param, double &value)
286 {
287 switch(BdComp)
288  {
289  case 0: value=0;
290  break;
291  case 1: value=0;
292  break;
293  case 2: value=0;
294  break;
295  case 3: value=0;
296  break;
297  case 4: value=0;
298  break;
299  default: cout << "wrong boundary part number" << endl;
300  break;
301  }
302  return;
303 }
304 
305 
306 void LinCoeffs_DFT(int n_points, double *X, double *Y,
307  double **parameters, double **coeffs)
308 {
309  int i;
310  double *coeff, x, y, *param;
311  double u1x, u1y, u2x, u2y;
312 
313  for(i=0;i<n_points;i++)
314  {
315  coeff = coeffs[i];
316 
317  x = X[i];
318  y = Y[i];
319 
320 // param = parameters[i];
321 //
322 // u1x = param[2];
323 // u1y = param[4];
324 // u2x = param[3];
325 // u2y = param[5];
326 
327  coeff[0] = 1;
328 // coeff[1] = u1x;
329 // coeff[2] = (u1y+u2x)*0.5;
330 // coeff[3] = u2y;
331 
332  }
333 }
334 
335 
337 void GetCdCl(TFEFunction2D *u1fct, TFEFunction2D *u2fct,
338  TFEFunction2D *pfct, TFEFunction2D *tau1fct,
339  TFEFunction2D *tau2fct, TFEFunction2D *tau3fct,
340  double &cd, double &cl)
341 {
342  int i,j,k,l, N_;
343  int N_Points,N_Edges,comp,ed_nr;
344  double *weights, *xi, *eta;
345  double X[MaxN_QuadPoints_2D];
346  double Y[MaxN_QuadPoints_2D];
347  double AbsDetjk[MaxN_QuadPoints_2D];
348  int N_LocalUsedElements;
349  FE2D LocalUsedElements[3], CurrentElement;
350  int *DOF;
351  double **OrigFEValues, *Orig;
352  bool SecondDer[3] = { FALSE, FALSE, FALSE };
353  double *u1, *u2, *p, *tau1, *tau2, *tau3;
354  TFESpace2D *USpace, *PSpace, *TauSpace;
355  int *UGlobalNumbers, *UBeginIndex;
356  int *PGlobalNumbers, *PBeginIndex;
357  int *TauGlobalNumbers, *TauBeginIndex;
358  int *N_BaseFunct, N_Cells;
359  BaseFunct2D BaseFunct, *BaseFuncts;
360  TCollection *Coll;
361  TBaseCell *cell;
362  double value, value1, value2, value3, value4, value5, value6;
363  double FEFunctValues[MaxN_BaseFunctions2D];
364  double FEFunctValues1[MaxN_BaseFunctions2D];
365  double FEFunctValues2[MaxN_BaseFunctions2D];
366  double FEFunctValues3[MaxN_BaseFunctions2D];
367  double FEFunctValues4[MaxN_BaseFunctions2D];
368  double FEFunctValues5[MaxN_BaseFunctions2D];
369  double FEFunctValues6[MaxN_BaseFunctions2D];
370 
371 
372  int N_DerivativesU = 3;
373  double *Derivatives[MaxN_BaseFunctions2D];
374  MultiIndex2D NeededDerivatives[3] = { D00, D10, D01 };
375  TFEFunction2D *vfct;
376  double *v, nu = 1/TDatabase::ParamDB->RE_NR;
377  double wei = 1/TDatabase::ParamDB->WEI_NR, beta = TDatabase::ParamDB->P3;
378  double *Der, *aux;
379  TJoint *joint;
380  TBoundEdge *boundedge;
381  TBoundComp *BoundComp;
382  TFE2D *eleCell;
383  FE2D FEEle;
384  TFEDesc2D *FEDesc;
385  int N_DOF_Circ, *DOF_Circ;
386  char VString[] = "v";
387 
388  u1 = u1fct->GetValues();
389  u2 = u2fct->GetValues();
390  p = pfct->GetValues();
391  tau1 = tau1fct->GetValues();
392  tau2 = tau2fct->GetValues();
393  tau3 = tau3fct->GetValues();
394 
395  USpace = u1fct->GetFESpace2D();
396  PSpace = pfct->GetFESpace2D();
397  TauSpace = tau1fct->GetFESpace2D();
398 
399  UGlobalNumbers = USpace->GetGlobalNumbers();
400  UBeginIndex = USpace->GetBeginIndex();
401 
402  PGlobalNumbers = PSpace->GetGlobalNumbers();
403  PBeginIndex = PSpace->GetBeginIndex();
404 
405  TauGlobalNumbers = TauSpace->GetGlobalNumbers();
406  TauBeginIndex = TauSpace->GetBeginIndex();
407 
408  BaseFuncts = TFEDatabase2D::GetBaseFunct2D_IDFromFE2D();
409  N_BaseFunct = TFEDatabase2D::GetN_BaseFunctFromFE2D();
410 
411  aux = new double [MaxN_QuadPoints_2D*13];
412  for(j=0;j<MaxN_QuadPoints_2D;j++)
413  Derivatives[j] = aux + j*13;
414 
415  N_ = u1fct->GetLength();
416  v = new double[N_];
417  memset(v,0,N_*SizeOfDouble);
418  vfct = new TFEFunction2D(USpace, VString, VString, v, N_);
419 
420 // ########################################################################
421 // loop over all cells
422 // ########################################################################
423  Coll = USpace->GetCollection(); // all spaces use same Coll
424  N_Cells = Coll->GetN_Cells();
425 
426  for(i=0;i<N_Cells;i++)
427  {
428  cell = Coll->GetCell(i);
429  N_Edges=cell->GetN_Edges();
430  for(j=0;j<N_Edges;j++) // loop over all edges of cell
431  {
432  joint=cell->GetJoint(j);
433  if ((joint->GetType() == BoundaryEdge)||
434  (joint->GetType() == IsoBoundEdge)) // boundary edge
435  {
436 
437  boundedge=(TBoundEdge *)joint;
438  BoundComp=boundedge->GetBoundComp(); // get boundary component
439  comp=BoundComp->GetID(); // boundary id
440  if (comp==4)
441  {
442  FEEle = USpace->GetFE2D(i,cell); // finite element of cell
443  eleCell = TFEDatabase2D::GetFE2D(FEEle);
444  FEDesc = eleCell->GetFEDesc2D(); // fe descriptor
445  N_DOF_Circ = FEDesc->GetN_JointDOF(); // # local dofs on joints
446  DOF_Circ = FEDesc->GetJointDOF(j); // local dofs on joint j
447  DOF = UGlobalNumbers + UBeginIndex[i]; // pointer to global dofs
448  for (k=0;k<N_DOF_Circ;k++) // set fe on circle to 1
449  v[DOF[DOF_Circ[k]]] = 1;
450  }
451  }
452  }
453  }
454 
455  cd = 0;
456  cl = 0;
457 
458 // ########################################################################
459 // loop over all cells
460 // ########################################################################
461  Coll = USpace->GetCollection(); // all spaces use same Coll
462  N_Cells = Coll->GetN_Cells();
463  for(i=0;i<N_Cells;i++)
464  {
465  cell = Coll->GetCell(i);
466 
467  // ####################################################################
468  // find local used elements on this cell
469  // ####################################################################
470  N_LocalUsedElements = 3;
471  LocalUsedElements[0] = USpace->GetFE2D(i, cell);
472  LocalUsedElements[1] = PSpace->GetFE2D(i, cell);
473  LocalUsedElements[2] = TauSpace->GetFE2D(i, cell);
474 
475  // ####################################################################
476  // calculate values on original element
477  // ####################################################################
478  TFEDatabase2D::GetOrig(N_LocalUsedElements, LocalUsedElements, Coll, cell, SecondDer, N_Points, xi, eta, weights, X, Y, AbsDetjk);
479 
480  // calculate all needed values of p
481  CurrentElement = LocalUsedElements[1];
482  BaseFunct = BaseFuncts[CurrentElement];
483  N_ = N_BaseFunct[CurrentElement];
484 
485  DOF = PGlobalNumbers + PBeginIndex[i];
486  for(l=0;l<N_;l++)
487  FEFunctValues[l] = p[DOF[l]];
488 
489  OrigFEValues = TFEDatabase2D::GetOrigElementValues(BaseFunct, D00);
490 
491  for(j=0;j<N_Points;j++)
492  {
493  Orig = OrigFEValues[j];
494  value = 0;
495  for(l=0;l<N_;l++)
496  value += FEFunctValues[l] * Orig[l];
497 
498  Derivatives[j][0] = value;
499  }
500 
501  // calculate all needed values of u1, u2
502  CurrentElement = LocalUsedElements[0];
503  BaseFunct = BaseFuncts[CurrentElement];
504  N_ = N_BaseFunct[CurrentElement];
505 
506  DOF = UGlobalNumbers + UBeginIndex[i];
507  for(l=0;l<N_;l++)
508  {
509  FEFunctValues1[l] = u1[DOF[l]];
510  FEFunctValues2[l] = u2[DOF[l]];
511  FEFunctValues3[l] = v[DOF[l]];
512  }
513 
514  for(k=0;k<N_DerivativesU;k++)
515  {
516  OrigFEValues = TFEDatabase2D::GetOrigElementValues(BaseFunct,
517  NeededDerivatives[k]);
518  for(j=0;j<N_Points;j++)
519  {
520  Orig = OrigFEValues[j];
521  value1 = 0;
522  value2 = 0;
523  value3 = 0;
524  for(l=0;l<N_;l++)
525  {
526  value1 += FEFunctValues1[l] * Orig[l];
527  value2 += FEFunctValues2[l] * Orig[l];
528  value3 += FEFunctValues3[l] * Orig[l];
529  } // endfor l
530  Derivatives[j][k+1] = value1;
531  Derivatives[j][k+4] = value2;
532  Derivatives[j][k+7] = value3;
533  } // endfor j
534  } // endfor k
535 
536 
537  // calculate all needed values of tau1, tau2, tau3
538  CurrentElement = LocalUsedElements[2];
539  BaseFunct = BaseFuncts[CurrentElement];
540  N_ = N_BaseFunct[CurrentElement];
541 
542  DOF = TauGlobalNumbers + TauBeginIndex[i];
543  for(l=0;l<N_;l++)
544  {
545  FEFunctValues4[l] = tau1[DOF[l]];
546  FEFunctValues5[l] = tau2[DOF[l]];
547  FEFunctValues6[l] = tau3[DOF[l]];
548 
549  }
550  OrigFEValues = TFEDatabase2D::GetOrigElementValues(BaseFunct, D00);
551 
552  for(j=0;j<N_Points;j++)
553  {
554  Orig = OrigFEValues[j];
555  value4 = 0;
556  value5 = 0;
557  value6 = 0;
558  for(l=0;l<N_;l++)
559  {
560  value4 += FEFunctValues4[l] * Orig[l];
561  value5 += FEFunctValues5[l] * Orig[l];
562  value6 += FEFunctValues6[l] * Orig[l];
563  } // endfor l
564  Derivatives[j][10] = value4;
565  Derivatives[j][11] = value5;
566  Derivatives[j][12] = value6;
567  }
568 
569 
570 
571  // calculation
572  for(j=0;j<N_Points;j++)
573  {
574  Der = Derivatives[j];
575 
576  // beta*nu * (u1_x*v_x, u1_y*v_y), v= (v,0)
577  value1 = beta*nu*(Der[2]*Der[8]+Der[3]*Der[9]);
578  // (u1 * u1_x + u2* u1_y) * (1,0)
579  value1 += (Der[1]*Der[2]+Der[4]*Der[3])*Der[7];
580  // pressure times divergence of test function (1,0)
581  value1 -= Der[0]*Der[8];
582  // (1-beta)*nu*wei* (tau1*v_x + tau2*v_y)
583  value1 += (1.0-beta)*nu*wei*(Der[10]*Der[8] + Der[11]*Der[9]);
584 
585  value2 = beta*nu*(Der[5]*Der[8]+Der[6]*Der[9]);
586  value2 += (Der[1]*Der[5]+Der[4]*Der[6])*Der[7];
587  value2 -= Der[0]*Der[9];
588  value2 += (1.0-beta)*nu*wei*(Der[11]*Der[8] + Der[12]*Der[9]);
589 
590  cd += AbsDetjk[j]*weights[j] * value1;
591  cl += AbsDetjk[j]*weights[j] * value2;
592  }
593 
594  } // endfor i
595 
596  cd *= -500;
597  cl *= -500;
598 
599  delete Derivatives[0];
600  delete vfct;
601  delete v;
602 }
TFEFunction2D::GetValues
double * GetValues()
Definition: FEFunction2D.h:67
TFEDatabase2D::GetFE2D
static TFE2D * GetFE2D(FE2D FE)
Definition: FEDatabase2D.h:353
TFEFunction2D::GetLength
int GetLength()
Definition: FEFunction2D.h:63
TFESpace::GetCollection
TCollection * GetCollection() const
Definition: FESpace.h:131
TJoint::GetType
JointType GetType()
Definition: Joint.h:75
TCollection::GetCell
TBaseCell * GetCell(int i) const
return Cell with index i in Cells-array
Definition: Collection.h:50
TBoundEdge::GetBoundComp
TBoundComp2D * GetBoundComp() const
Definition: BoundEdge.h:77
TFE2D
Definition: FE2D.h:20
TFESpace2D
Definition: FESpace2D.h:28
TParaDB::RE_NR
double RE_NR
Definition: Database.h:313
TFEDatabase2D::GetOrigElementValues
static double ** GetOrigElementValues(BaseFunct1D BaseFunct, MultiIndex1D MultiIndex)
Definition: FEDatabase2D.h:300
TCollection
store cells in an array, used by cell iterators
Definition: Collection.h:18
TCollection::GetN_Cells
int GetN_Cells() const
return number of cells
Definition: Collection.h:46
TBaseCell::GetJoint
TJoint * GetJoint(int J_i)
return the pointer to face with number i
Definition: BaseCell.h:175
TJoint
Definition: Joint.h:48
TFEFunction2D::GetFESpace2D
TFESpace2D * GetFESpace2D()
Definition: FEFunction2D.h:59
TFESpace2D::GetFE2D
FE2D GetFE2D(int i, TBaseCell *cell)
Definition: FESpace2D.C:1184
TFE2D::GetFEDesc2D
TFEDesc2D * GetFEDesc2D() const
Definition: FE2D.h:97
TFEDesc2D::GetJointDOF
int ** GetJointDOF() const
Definition: FEDesc2D.h:86
TFEDatabase2D::GetBaseFunct2D_IDFromFE2D
static BaseFunct2D * GetBaseFunct2D_IDFromFE2D()
Definition: FEDatabase2D.h:417
TFEDatabase2D::GetN_BaseFunctFromFE2D
static int * GetN_BaseFunctFromFE2D()
Definition: FEDatabase2D.h:421
TBoundComp::GetID
int GetID() const
Definition: BoundComp.h:49
TBaseCell
information for finite element data structure
Definition: BaseCell.h:25
TFESpace::GetGlobalNumbers
int * GetGlobalNumbers() const
Definition: FESpace.h:135
TBaseCell::GetN_Edges
int GetN_Edges()
return the number of edges of the cell
Definition: BaseCell.h:182
TBoundEdge
Definition: BoundEdge.h:19
TBoundComp
Definition: BoundComp.h:27
TFEDatabase2D::GetOrig
static RefTrans2D GetOrig(int N_LocalUsedElements, FE2D *LocalUsedElements, TCollection *Coll, TBaseCell *cell, bool *Needs2ndDer, int &N_Points, double *&xi, double *&eta, double *&weights, double *X, double *Y, double *absdetjk)
Definition: FEDatabase2D.C:1765
TFEDesc2D
Definition: FEDesc2D.h:15
TFESpace::GetBeginIndex
int * GetBeginIndex() const
Definition: FESpace.h:142
TFEDesc2D::GetN_JointDOF
int GetN_JointDOF() const
Definition: FEDesc2D.h:61
TDatabase::ParamDB
static TParamDB * ParamDB
Definition: Database.h:1134
TFEFunction2D
Definition: FEFunction2D.h:24
Generated on Tue May 10 2016 09:07:56 for ParMooN by   doxygen 1.8.5