ParMooN
 All Classes Functions Variables Friends Pages
Hemker1996.h
1 // ======================================================================
2 // Example from P.W. Hemker
3 // ======================================================================
4 #define __HEMKER1996__
5 
6 void ExampleFile()
7 {
8  OutPut("Example: Hemker1996.h" << endl) ;
9 }
10 
11 // exact solution
12 void Exact(double x, double y, double *values)
13 {
14  values[0] = 0;
15  values[1] = 0;
16  values[2] = 0;
17  values[3] = 0;
18 }
19 
20 void BoundCondition(int BdComp, double t, BoundCond &cond)
21 {
22  switch(BdComp)
23  {
24  case 0:
25  case 1:
26  case 2:
27  cond = NEUMANN;
28  break;
29  default:
30  cond = DIRICHLET;
31  }
32 }
33 
34 
35 // value of boundary condition
36 void BoundValue(int BdComp, double Param, double &value)
37 {
38  switch(BdComp)
39  {
40  case 1:
41  value = 0;
42  break;
43  case 4:
44  value = 1;
45  break;
46  default:
47  value = 0;
48  }
49 }
50 
51 // initial conditon
52 void InitialCondition(double x, double y, double *values)
53 {
54  values[0] = 0;
55 }
56 
57 void BoundConditionAdjoint(int BdComp, double t, BoundCond &cond)
58 {
59  switch(BdComp)
60  {
61  case 0:
62  case 2:
63  case 3:
64  cond = NEUMANN;
65  break;
66  default:
67  cond = DIRICHLET;
68  }
69 }
70 
71 
72 // value of boundary condition
73 void BoundValueAdjoint(int BdComp, double Param, double &value)
74 {
75  value = 0;
76 }
77 
78 void BilinearCoeffs(int n_points, double *x, double *y,
79  double **parameters, double **coeffs)
80 {
81  double eps=1/TDatabase::ParamDB->RE_NR;
82  double angle = 0, v1, v2;
83  int i;
84  double *coeff;
85 
86  v1 = cos(angle);
87  v2 = sin(angle);
88 
89  for(i=0;i<n_points;i++)
90  {
91  coeff = coeffs[i];
92 
93  coeff[0] = eps;
94  coeff[1] = v1;
95  coeff[2] = v2;
96  coeff[3] = 0;
97 
98  coeff[4] = 0;
99  }
100 }
101 
102 void ComputeExtremalValues(int N, double *sol, double *values)
103 {
104  int i;
105  double max, min;
106 
107  min = 1e10;
108  max = -1e10;
109 
110  for(i=0;i<N;i++)
111  {
112  if(sol[i]-1 > max)
113  max = sol[i]-1;
114  if(sol[i] < min)
115  min = sol[i];
116  }
117 
118  values[0] = min;
119  values[1] = max;
120 }
121 
123 void ComputeOutflowBoundary(int level, TFEFunction2D *ufct)
124 {
125  double h, x=4,values[3],y;
126  int i, bound_points = 401;
127  h = 6.0/(bound_points-1);
128  for (i=0;i<bound_points; i++)
129  {
130  y = -3+i*h;
131  ufct->FindGradient(x,y,values);
132  OutPut("cutline " << x << " " << y <<
133  " " << values[0] << endl);
134  }
135 }
136 
137 // computation of some global errors, only for P1 or Q1 !!!
138 //
139 // values[0] : absolute value of largest negative undershoot in
140 // a circle around the cylinder
141 // values[1] : difference of largest positive value in a circle
142 // around the cylinder and 1
143 // values[2] : absolute value of largest negative undershoot
144 // for x > 2
145 // values[3] : difference of largest positive value for x>2 and 1
146 //
147 void ComputeLocalExtrema(TFEFunction2D *ufct, double *values)
148 {
149  TBaseCell *cell;
150  TCollection *Coll;
151  TFESpace2D *FESpace2D;
152  RefTrans2D RefTrans;
153  TBaseFunct2D *bf;
154  FE2D FE_ID;
155  TFE2D *FE_Obj;
156  int N_BaseFunct;
157  double xi, eta;
158  double *uorig, *uxorig, *uyorig, *uref, *uxiref, *uetaref, u;
159  double *Values, x, y, val;
160  int *GlobalNumbers, *BeginIndex, N_Cells, N_Edges;
161  int i, j, k, *Numbers;
162  double extr[4];
163 
164  extr[0] = -1;
165  extr[1] = -1;
166  extr[2] = -1;
167  extr[3] = 0;
168 
169  FESpace2D = ufct->GetFESpace2D();
170  BeginIndex = FESpace2D->GetBeginIndex();
171  GlobalNumbers = FESpace2D->GetGlobalNumbers();
172  Values = ufct->GetValues();
173 
174  Coll = FESpace2D->GetCollection();
175  N_Cells = Coll->GetN_Cells();
176 
177  // loop over all edges
178  for(i=0;i<N_Cells;i++)
179  {
180  cell = Coll->GetCell(i);
181  N_Edges=cell->GetN_Edges();
182  //double diam = cell->GetDiameter();
183 
184  FE_ID = FESpace2D->GetFE2D(i, cell);
185  FE_Obj = TFEDatabase2D::GetFE2D(FE_ID);
186  RefTrans = FE_Obj->GetRefTransID();
187 
188  // get base function object
189  bf = FE_Obj->GetBaseFunct2D();
190  N_BaseFunct = bf->GetDimension();
191 
192  uorig = new double[N_BaseFunct];
193  uxorig = new double[N_BaseFunct];
194  uyorig = new double[N_BaseFunct];
195 
196  uref = new double[N_BaseFunct];
197  uxiref = new double[N_BaseFunct];
198  uetaref = new double[N_BaseFunct];
199 
200  // set cell for reference transformation
201  TFEDatabase2D::SetCellForRefTrans(cell, RefTrans);
202  for (j=0;j<N_Edges;j++)
203  {
204  // compute coordinates
205  x = cell->GetVertex(j)->GetX();
206  y = cell->GetVertex(j)->GetY();
207  if (x<-1.5)
208  continue;
209  // find local coordinates of the given point
210  //cout << " x: " << x << endl;
211  //cout << " y: " << y << endl;
212  TFEDatabase2D::GetRefFromOrig(RefTrans, x, y, xi, eta);
213  //cout << " xi: " << xi << endl;
214  //cout << "eta: " << eta << endl;
215 
216  bf->GetDerivatives(D00, xi, eta, uref);
217  bf->GetDerivatives(D10, xi, eta, uxiref);
218  bf->GetDerivatives(D01, xi, eta, uetaref);
219 
220  TFEDatabase2D::GetOrigValues(RefTrans, xi, eta, bf, Coll, (TGridCell *)cell,
221  uref, uxiref, uetaref, uorig, uxorig, uyorig);
222  // compute value of fe function at (x,y)
223  u = 0;
224  Numbers = GlobalNumbers + BeginIndex[i];
225  for(k=0;k<N_BaseFunct;k++)
226  {
227  val = Values[Numbers[k]];
228  u += uorig[k]*val;
229  }
230  //OutPut(x << " " << y << " " << u << endl);
231  // strip with the circle
232  if ((x>-1.5)&&(x<1.5))
233  {
234  if ((u<=0)&&(fabs(u)>extr[0]))
235  extr[0] = fabs(u);
236  if ((u>=1)&&(fabs(u-1)>extr[1]))
237  extr[1] = fabs(u-1);
238  }
239  if (x>2)
240  {
241  if ((u<=0)&&(fabs(u)>extr[2]))
242  extr[2] = fabs(u);
243  if ((u>=1)&&(fabs(u-1)>extr[3]))
244  extr[3] = fabs(u-1);
245  }
246  } // endfor (j) N_Edges
247 
248  delete uorig;
249  delete uxorig;
250  delete uyorig;
251  delete uref;
252  delete uxiref;
253  delete uetaref;
254  } // endfor
255 
256  values[0] = extr[0];
257  values[1] = extr[1];
258  values[2] = extr[2];
259  values[3] = extr[3];
260  }
261 
262 /****************************************************************/
263 //
264 // for FEM_TVD
265 //
266 /****************************************************************/
267 
268 void CheckWrongNeumannNodes(TCollection *Coll, TFESpace2D *fespace,
269  int &N_neum_to_diri, int* &neum_to_diri,
270  int* &neum_to_diri_bdry,
271  double* &neum_to_diri_param)
272 {
273  const int max_entries = 50000;
274  int i, j, min_val, type;
275  int N_Cells, N_V, diri_counter = 0, found, diri_counter_1 = 0;
276  int *global_numbers, *begin_index, *dof;
277  int boundary_vertices[4], tmp_diri[max_entries], tmp_bdry[max_entries];
278  double x[4], y[4], eps = 1e-6, tmp_param[max_entries];
279  TBaseCell *cell;
280  TVertex *vertex;
281  FE2D CurrentElement;
282 
283  // number of mesh cells
284  N_Cells = Coll->GetN_Cells();
285  // array with global numbers of d.o.f.
286  global_numbers = fespace->GetGlobalNumbers();
287  // array which points to the beginning of the global numbers in
288  // global_numbers for each mesh cell
289  begin_index = fespace->GetBeginIndex();
290 
291  diri_counter = 0;
292  for(i=0;i<N_Cells;i++)
293  {
294  cell = Coll->GetCell(i);
295  N_V = cell->GetN_Vertices();
296  found = 0;
297  for (j=0;j<N_V;j++)
298  {
299  // read coordinates of the mesh cell
300  boundary_vertices[j] = 0;
301  vertex = cell->GetVertex(j);
302  vertex->GetCoords(x[j], y[j]);
303  if ((fabs(x[j]+3)<eps)//||(fabs(y[j]+3)<eps)||(fabs(y[j]-3)<eps)
304  || (fabs(sqrt(x[j]*x[j]+y[j]*y[j])-1)<eps))
305  {
306  boundary_vertices[j] = 1;
307  found++;
308  }
309  }
310  // no cell with edge with vertex on the boundary
311  if (found<2)
312  continue;
313  // finite element on the mesh cell
314  CurrentElement = fespace->GetFE2D(i, cell);
315  // number of basis functions (= number of d.o.f.)
316  //int N_ = TFEDatabase2D::GetN_BaseFunctFromFE2D(CurrentElement);
317  // the array which gives the mapping of the local to the global d.o.f.
318  dof = global_numbers+begin_index[i];
319  switch(CurrentElement)
320  {
321  // P_1, Q_1
322  case C_P1_2D_T_A:
323  case C_Q1_2D_Q_A:
324  case C_Q1_2D_Q_M:
325  for (j=0;j<N_V;j++)
326  {
327  // vertex on the boundary
328  if (boundary_vertices[j])
329  {
330  if (CurrentElement==C_P1_2D_T_A)
331  tmp_diri[diri_counter] = dof[j];
332  else
333  {
334  if (j<2){
335  tmp_diri[diri_counter] = dof[j];
336  }
337  else
338  {
339  if (j==2)
340  tmp_diri[diri_counter] = dof[3];
341  else
342  tmp_diri[diri_counter] = dof[2];
343  }
344  }
345  if (diri_counter > max_entries)
346  {
347  OutPut("tmp_diri too short !!!"<<endl);
348  exit(4711);
349  }
350  // inflow x = -3
351  if (fabs(x[j]+3)<eps)
352  {
353  tmp_bdry[diri_counter] = 3;
354  tmp_param[diri_counter] = (-y[j]+3)/6.0;
355  }
356  // circle
357  if (fabs(sqrt(x[j]*x[j]+y[j]*y[j])-1)<eps)
358  {
359  tmp_bdry[diri_counter] = 4;
360  // parameter does not matter, since b.c. equal to 1
361  tmp_param[diri_counter] = 0;
362  }
363  diri_counter++;
364  }
365  }
366  break;
367  // P_2, Q_2
368  case C_P2_2D_T_A:
369  case C_Q2_2D_Q_A:
370  case C_Q2_2D_Q_M:
371  // loop over the edges
372  for (j=0;j<N_V;j++)
373  {
374  // check of edge j is on boundary
375  if (boundary_vertices[j] && boundary_vertices[(j+1)%N_V])
376  {
377  // check if this is a boundary edge
378  type = cell->GetJoint(j)->GetType();
379  if (!((type == BoundaryEdge)||(type == IsoBoundEdge)))
380  continue;
381  switch(j)
382  {
383  case 0:
384  tmp_diri[diri_counter] = dof[0];
385  tmp_diri[diri_counter+1] = dof[1];
386  tmp_diri[diri_counter+2] = dof[2];
387  break;
388  case 1:
389  if (N_V==3)
390  {
391  tmp_diri[diri_counter] = dof[2];
392  tmp_diri[diri_counter+1] = dof[4];
393  tmp_diri[diri_counter+2] = dof[5];
394  }
395  else
396  {
397  tmp_diri[diri_counter] = dof[2];
398  tmp_diri[diri_counter+1] = dof[5];
399  tmp_diri[diri_counter+2] = dof[8];
400  }
401  break;
402  case 2:
403  if (N_V==3)
404  {
405  tmp_diri[diri_counter] = dof[5];
406  tmp_diri[diri_counter+1] = dof[3];
407  tmp_diri[diri_counter+2] = dof[0];
408  }
409  else
410  {
411  tmp_diri[diri_counter] = dof[8];
412  tmp_diri[diri_counter+1] = dof[7];
413  tmp_diri[diri_counter+2] = dof[6];
414  }
415  break;
416  case 3:
417  tmp_diri[diri_counter] = dof[6];
418  tmp_diri[diri_counter+1] = dof[3];
419  tmp_diri[diri_counter+2] = dof[0];
420  break;
421 
422  }
423 
424  if (diri_counter+2 > max_entries)
425  {
426  OutPut("tmp_diri too short !!!"<<endl);
427  exit(4711);
428  }
429 
430  // inflow x = -3
431  if ((fabs(x[j]+3)<eps)&&(fabs(x[(j+1)%N_V]+3)<eps))
432  {
433  tmp_bdry[diri_counter] = 3;
434  tmp_bdry[diri_counter+1] = 3;
435  tmp_bdry[diri_counter+2] = 3;
436  tmp_param[diri_counter] = (-y[j]+3)/6.0;
437  tmp_param[diri_counter+2] = (-y[(j+1)%N_V]+3)/6.0;
438  tmp_param[diri_counter+1] = (tmp_param[diri_counter] + tmp_param[diri_counter+2])/2.0;
439  }
440  // circle
441  if ((fabs(sqrt(x[j]*x[j]+y[j]*y[j])-1)<eps) &&
442  (fabs(sqrt(x[(j+1)%N_V]*x[(j+1)%N_V]+y[(j+1)%N_V]*y[(j+1)%N_V])-1)<eps))
443  {
444  tmp_bdry[diri_counter] = 4;
445  tmp_bdry[diri_counter+1] = 4;
446  tmp_bdry[diri_counter+2] = 4;
447  // parameter does not matter, since b.c. equal to 1
448  tmp_param[diri_counter] = 0;
449  tmp_param[diri_counter+1] = 0;
450  tmp_param[diri_counter+2] = 0;
451  }
452  diri_counter +=3;
453  }
454  }
455  break;
456  // P_3, Q_3
457  case C_P3_2D_T_A:
458  case C_Q3_2D_Q_A:
459  case C_Q3_2D_Q_M:
460  // loop over the edges
461  for (j=0;j<N_V;j++)
462  {
463  // check of edge j is on boundary
464  if (boundary_vertices[j] && boundary_vertices[(j+1)%N_V])
465  {
466  // check if this is a boundary edge
467  type = cell->GetJoint(j)->GetType();
468  if (!((type == BoundaryEdge)||(type == IsoBoundEdge)))
469  continue;
470 
471  // P3: local dof 0, 1, 2, 3 are on the boundary
472  // Q3: local dof 0, 1, 2, 3 are on the boundary
473  switch(j)
474  {
475  case 0:
476  tmp_diri[diri_counter] = dof[0];
477  tmp_diri[diri_counter+1] = dof[1];
478  tmp_diri[diri_counter+2] = dof[2];
479  tmp_diri[diri_counter+3] = dof[3];
480  break;
481  case 1:
482  if (N_V==3)
483  {
484  tmp_diri[diri_counter] = dof[3];
485  tmp_diri[diri_counter+1] = dof[6];
486  tmp_diri[diri_counter+2] = dof[8];
487  tmp_diri[diri_counter+3] = dof[9];
488  }
489  else
490  {
491  tmp_diri[diri_counter] = dof[3];
492  tmp_diri[diri_counter+1] = dof[7];
493  tmp_diri[diri_counter+2] = dof[11];
494  tmp_diri[diri_counter+3] = dof[15];
495  }
496  break;
497  case 2:
498  if (N_V==3)
499  {
500  tmp_diri[diri_counter] = dof[9];
501  tmp_diri[diri_counter+1] = dof[7];
502  tmp_diri[diri_counter+2] = dof[4];
503  tmp_diri[diri_counter+3] = dof[0];
504  }
505  else
506  {
507  tmp_diri[diri_counter] = dof[15];
508  tmp_diri[diri_counter+1] = dof[14];
509  tmp_diri[diri_counter+2] = dof[13];
510  tmp_diri[diri_counter+3] = dof[12];
511  }
512  break;
513  case 3:
514  tmp_diri[diri_counter] = dof[12];
515  tmp_diri[diri_counter+1] = dof[8];
516  tmp_diri[diri_counter+2] = dof[4];
517  tmp_diri[diri_counter+3] = dof[0];
518  break;
519  }
520 
521  if (diri_counter+3 > max_entries)
522  {
523  OutPut("tmp_diri too short !!!"<<endl);
524  exit(4711);
525  }
526 
527  // inflow x = -3
528  if ((fabs(x[j]+3)<eps)&&(fabs(x[(j+1)%N_V]+3)<eps))
529  {
530  tmp_bdry[diri_counter] = 3;
531  tmp_bdry[diri_counter+1] = 3;
532  tmp_bdry[diri_counter+2] = 3;
533  tmp_bdry[diri_counter+3] = 3;
534  tmp_param[diri_counter] = (-y[j]+3)/6.0;
535  tmp_param[diri_counter+3] = (-y[(j+1)%N_V]+3)/6.0;
536  tmp_param[diri_counter+1] = (2*tmp_param[diri_counter] + tmp_param[diri_counter+3])/3.0;
537  tmp_param[diri_counter+2] = (tmp_param[diri_counter] + 2*tmp_param[diri_counter+3])/2.0;
538  }
539  // circle
540  if ((fabs(sqrt(x[j]*x[j]+y[j]*y[j])-1)<eps) &&
541  (fabs(sqrt(x[(j+1)%N_V]*x[(j+1)%N_V]+y[(j+1)%N_V]*y[(j+1)%N_V])-1)<eps))
542  {
543  tmp_bdry[diri_counter] = 4;
544  tmp_bdry[diri_counter+1] = 4;
545  tmp_bdry[diri_counter+2] = 4;
546  tmp_bdry[diri_counter+3] = 4;
547  // parameter does not matter, since b.c. equal to 1
548  tmp_param[diri_counter] = 0;
549  tmp_param[diri_counter+1] = 0;
550  tmp_param[diri_counter+2] = 0;
551  tmp_param[diri_counter+3] = 0;
552  }
553  diri_counter +=4;
554  }
555  }
556  break;
557  default:
558  OutPut("CheckNeumannNodesForVelocity not implemented for element "
559  << CurrentElement << endl);
560  OutPut("code can be run without CheckNeumannNodesForVelocity, just delete the exit" << endl);
561  exit(4711);
562  }
563  }
564 
565  // condense
566  for (i=0;i<diri_counter;i++)
567  {
568  if (tmp_diri[i] == -1)
569  continue;
570  diri_counter_1++;
571  for (j=i+1;j<diri_counter;j++)
572  {
573  if (tmp_diri[i] == tmp_diri[j])
574  {
575  tmp_diri[j] = -1;
576  }
577  }
578  }
579 
580  //OutPut("CheckNeumannNodesForVelocity: N_neum_to_diri " << diri_counter_1 << endl);
581  N_neum_to_diri = diri_counter_1;
582  // allocate array for the indices
583  neum_to_diri = new int[diri_counter_1];
584  // allocate array for the corresponding boundary numbers
585  neum_to_diri_bdry = new int[diri_counter_1];
586  // allocate array for the corresponding boundary parameters
587  neum_to_diri_param = new double[diri_counter_1];
588  // fill array and sort
589  for (i=0;i<diri_counter_1;i++)
590  {
591  min_val = tmp_diri[0];
592  found = 0;
593  for (j=1;j<diri_counter;j++)
594  {
595  if ((tmp_diri[j]>0) && ((tmp_diri[j] < min_val) ||
596  (min_val == -1)))
597  {
598  min_val = tmp_diri[j];
599  found = j;
600  }
601  }
602  neum_to_diri[i] = tmp_diri[found];
603  neum_to_diri_bdry[i] = tmp_bdry[found];
604  neum_to_diri_param[i] = tmp_param[found];
605  tmp_diri[found] = -1;
606  }
607 /*
608  for (i=0;i<diri_counter_1;i++)
609  {
610  OutPut(i << " " << neum_to_diri[i] << " " << neum_to_diri_bdry[i] <<
611  " " << neum_to_diri_param[i] << endl);
612  }
613 */
614 }
615 
616 void SetDirichletNodesFromNeumannNodes(TSquareMatrix2D **SQMATRICES,
617  double *rhs,
618  int N_neum_to_diri,
619  int *neum_to_diri,
620  int *neum_to_diri_bdry,
621  double *neum_to_diri_param)
622 {
623  TSquareMatrix2D *MatrixA;
624  double *Entries_A;
625  int i, l, l0, l1, index, *RowPtr_A, *KCol_A;
626 
627  MatrixA = SQMATRICES[0];
628  RowPtr_A = MatrixA->GetRowPtr();
629  KCol_A = MatrixA->GetKCol();
630  Entries_A = MatrixA->GetEntries();
631  // loop over dof to change
632  for (i=0;i<N_neum_to_diri;i++)
633  {
634  index = neum_to_diri[i];
635  l0 = RowPtr_A[index];
636  l1 = RowPtr_A[index+1];
637  for (l=l0;l<l1;l++)
638  {
639  // diagonal entry
640  if (KCol_A[l]==index)
641  Entries_A[l] = 1;
642  else
643  Entries_A[l] = 0;
644  }
645  // set boundary condition
646  BoundValue(neum_to_diri_bdry[i], neum_to_diri_param[i], rhs[index]);
647  }
648 }
649 
650 void ComputeDifferenceToCoarseLevel(TCollection *Coll_fine,
651  TCollection *Coll_coarse,
652  TFEFunction2D *u_fine,
653  TFEFunction2D *u_coarse)
654 {
655  int i, j, k, N_Cells, N_Edges, coarse_no;
656  double x, y, x_c, y_c, val_fine[4], val_coarse[4], c1err = -1, c1err_coarse = -1;
657  double x_err, y_err, x_err_c, y_err_c;
658  TBaseCell *cell, *parent;
659 
660  // number of cells
661  N_Cells = Coll_fine->GetN_Cells();
662 
663  // loop over all edges
664  for(i=0;i<N_Cells;i++)
665  {
666  // cell
667  cell = Coll_fine->GetCell(i);
668  // parent cell
669  parent = cell->GetParent();
670  coarse_no = Coll_coarse->GetIndex(parent);
671  //OutPut(coarse_no << " ");
672  // number of edges
673  N_Edges=cell->GetN_Edges();
674  for (j=0;j<N_Edges;j++)
675  {
676  cell->GetVertex(j)->GetCoords(x, y);
677  u_fine->FindGradientLocal(cell, i, x, y, val_fine);
678  u_coarse->FindGradientLocal(parent, coarse_no, x, y, val_coarse);
679  if (fabs(val_fine[0] - val_coarse[0]) > c1err)
680  {
681  c1err = fabs(val_fine[0] - val_coarse[0]);
682  x_err = x;
683  y_err = y;
684  }
685  for (k=0;k<N_Edges;k++)
686  {
687  parent->GetVertex(k)->GetCoords(x_c, y_c);
688  if ((fabs(x_c -x ) < 1e-6) && (fabs(y_c -y ) < 1e-6))
689  {
690  if (fabs(val_fine[0] - val_coarse[0]) > c1err_coarse)
691  {
692  c1err_coarse = fabs(val_fine[0] - val_coarse[0]);
693  x_err_c = x;
694  y_err_c = y;
695  }
696  }
697  }
698  }
699  }
700 
701  //OutPut("C1 error f " << c1err << " \\& " << x_err <<"," << y_err << endl);
702  //OutPut(" C1 error c " << c1err_coarse << " \\& " << x_err_c <<"," << y_err_c << endl);
703 
704  OutPut("C1 error f "<<" & " << c1err << " & ( " << x_err <<"," << y_err << ")"<< "\\\\\\hline" << endl);
705  OutPut("C1 error c "<< " & "<< c1err_coarse << " & ( " << x_err_c <<"," << y_err_c << ")" << "\\\\\\hline"<< endl);
706 }
707 void ComputeCutLines_X(TCollection *Coll, TFEFunction2D *ufct, int level)
708 {
709  double h, val[3], x, y, *cutvalues, y0=0;
710  int i, j, N_Cells, bound_points = 10001;
711  TBaseCell *cell;
712 
713  h = 3.0/(bound_points-1);
714 
715  cutvalues = new double[6*bound_points];
716  memset(cutvalues , 0 , 6*bound_points*SizeOfDouble);
717 
718  N_Cells = Coll->GetN_Cells();
719  for(i=0;i<N_Cells;i++)
720  {
721  cell = Coll->GetCell(i);
722  x = -1;
723  for (j=0;j<bound_points;j++)
724  {
725  y = y0 + h * j;
726  cutvalues[j] = y;
727  if(cell->PointInCell(x,y))
728  {
729  ufct->FindGradientLocal(cell, i, x, y, val);
730  cutvalues[j+bound_points] = val[0];
731  }
732  }
733  x = 0;
734  for (j=0;j<bound_points;j++)
735  {
736  y = y0 + h * j;
737  if(cell->PointInCell(x,y))
738  {
739  ufct->FindGradientLocal(cell, i, x, y, val);
740  cutvalues[j+2*bound_points] = val[0];
741  }
742  }
743  x = 1;
744  for (j=0;j<bound_points;j++)
745  {
746  y = y0 + h * j;
747  if(cell->PointInCell(x,y))
748  {
749  ufct->FindGradientLocal(cell, i, x, y, val);
750  cutvalues[j+3*bound_points] = val[0];
751  }
752  }
753  x = 2;
754  for (j=0;j<bound_points;j++)
755  {
756  y = y0 + h * j;
757  if(cell->PointInCell(x,y))
758  {
759  ufct->FindGradientLocal(cell, i, x, y, val);
760  cutvalues[j+4*bound_points] = val[0];
761  }
762  }
763  x = 4;
764  for (j=0;j<bound_points;j++)
765  {
766  y = y0 + h * j;
767  if(cell->PointInCell(x,y))
768  {
769  ufct->FindGradientLocal(cell, i, x, y, val);
770  cutvalues[j+5*bound_points] = val[0];
771  }
772  }
773  }
774 
775  for (j=0;j<bound_points;j++)
776  {
777  OutPut("cutx"<< level << " " << cutvalues[j] << " " << cutvalues[j+bound_points]
778  << " " << cutvalues[j+2*bound_points] << " " << cutvalues[j+3*bound_points] << " "
779  << cutvalues[j+4*bound_points] << " " << cutvalues[j+5*bound_points] << endl);
780  }
781 }
782 
783 void ComputeCutLines_Y(TCollection *Coll, TFEFunction2D *ufct, int level)
784 {
785  double h, val[3], x, y, *cutvalues;
786  int i, j, N_Cells, bound_points = 20001;
787  TBaseCell *cell;
788 
789  h = 10.0/(bound_points-1);
790 
791  cutvalues = new double[3*bound_points];
792  memset(cutvalues , 0 , 3*bound_points*SizeOfDouble);
793 
794  N_Cells = Coll->GetN_Cells();
795  for(i=0;i<N_Cells;i++)
796  {
797  cell = Coll->GetCell(i);
798  y = 0;
799  for (j=0;j<bound_points;j++)
800  {
801  x = -2 + h * j;
802  cutvalues[j] = x;
803  if(cell->PointInCell(x,y))
804  {
805  ufct->FindGradientLocal(cell, i, x, y, val);
806  cutvalues[j+bound_points] = val[0];
807  }
808  }
809  y = 1;
810  for (j=0;j<bound_points;j++)
811  {
812  x = -2 + h * j;
813  if(cell->PointInCell(x,y))
814  {
815  ufct->FindGradientLocal(cell, i, x, y, val);
816  /*if (cutvalues[j+2*bound_points]!=4711)
817  {
818  OutPut("belegt"<<endl);
819  exit(1);
820  }
821  else*/
822  cutvalues[j+2*bound_points] = val[0];
823  }
824  }
825  }
826 
827  for (j=0;j<bound_points;j++)
828  {
829  OutPut("cuty"<< level << " " << cutvalues[j] << " " << cutvalues[j+bound_points]
830  << " " << cutvalues[j+2*bound_points] << endl);
831  }
832 }
833 
834 void ComputeCutLines_epsY(TCollection *Coll, TFEFunction2D *ufct, int level)
835 {
836  double h, val[3], x, y, *cutvalues;
837  int i, j, N_Cells, bound_points = 20001;
838  TBaseCell *cell;
839  double eps = 1.0/TDatabase::ParamDB->RE_NR;
840 
841  h = 10.0/(bound_points-1);
842 
843  cutvalues = new double[3*bound_points];
844  memset(cutvalues , 0 , 3*bound_points*SizeOfDouble);
845 
846  N_Cells = Coll->GetN_Cells();
847  for(i=0;i<N_Cells;i++)
848  {
849  cell = Coll->GetCell(i);
850  y = 1-sqrt(eps);
851  for (j=0;j<bound_points;j++)
852  {
853  x = -2 + h * j;
854  cutvalues[j] = x;
855  if(cell->PointInCell(x,y))
856  {
857  ufct->FindGradientLocal(cell, i, x, y, val);
858  cutvalues[j+bound_points] = val[0];
859  }
860  }
861  y = 1+sqrt(eps);
862  for (j=0;j<bound_points;j++)
863  {
864  x = -2 + h * j;
865  if(cell->PointInCell(x,y))
866  {
867  ufct->FindGradientLocal(cell, i, x, y, val);
868  cutvalues[j+2*bound_points] = val[0];
869  }
870  }
871  }
872 
873  for (j=0;j<bound_points;j++)
874  {
875  OutPut("cutyeps"<< level << " " << cutvalues[j] << " " << cutvalues[j+bound_points]
876  << " " << cutvalues[j+2*bound_points] << endl);
877  }
878 }
879 void ComputeCutLines_eps_radial(TCollection *Coll, TFEFunction2D *ufct, int level)
880 {
881  double h, val[3], x, y, *cutvalues, tmp, r;
882  int i, j, N_Cells, bound_points = 10001;
883  TBaseCell *cell;
884  double eps = 1.0/TDatabase::ParamDB->RE_NR;
885 
886  h = 2*Pi/(bound_points-1);
887 
888  cutvalues = new double[10*bound_points];
889  memset(cutvalues , 0 , 10*bound_points*SizeOfDouble);
890 
891  N_Cells = Coll->GetN_Cells();
892  for(i=0;i<N_Cells;i++)
893  {
894  cell = Coll->GetCell(i);
895  r = (1+eps)*(1+eps);
896  for (j=0;j<bound_points;j++)
897  {
898  x = r*cos(h * j-Pi);
899  y = r*sin(h*j-Pi);
900  cutvalues[j] = h*j -Pi;
901  cutvalues[j+bound_points] = x;
902  cutvalues[j+2*bound_points] = y;
903  if(cell->PointInCell(x,y))
904  {
905  ufct->FindGradientLocal(cell, i, x, y, val);
906  cutvalues[j+3*bound_points] = val[0];
907  }
908  }
909  tmp=pow(eps,2.0/3.0);
910  r = (1+tmp)*(1+tmp);
911  for (j=0;j<bound_points;j++)
912  {
913  x = r*cos(h * j-Pi);
914  y = r*sin(h*j-Pi);
915  cutvalues[j+4*bound_points] = x;
916  cutvalues[j+5*bound_points] = y;
917  if(cell->PointInCell(x,y))
918  {
919  ufct->FindGradientLocal(cell, i, x, y, val);
920  cutvalues[j+6*bound_points] = val[0];
921  }
922  }
923  tmp=sqrt(eps);
924  r = (1+tmp)*(1+tmp);
925  for (j=0;j<bound_points;j++)
926  {
927  x = r*cos(h * j-Pi);
928  y = r*sin(h*j-Pi);
929  cutvalues[j+7*bound_points] = x;
930  cutvalues[j+8*bound_points] = y;
931  if(cell->PointInCell(x,y))
932  {
933  ufct->FindGradientLocal(cell, i, x, y, val);
934  cutvalues[j+9*bound_points] = val[0];
935  }
936  }
937 
938  }
939 
940  for (j=0;j<bound_points;j++)
941  {
942  OutPut("cutradeps"<< level << " " << cutvalues[j] << " " << cutvalues[j+bound_points]
943  << " " << cutvalues[j+2*bound_points] << " " << cutvalues[j+3*bound_points]
944  << " " << cutvalues[j+4*bound_points] << " " << cutvalues[j+5*bound_points]
945  << " " << cutvalues[j+6*bound_points] << " " << cutvalues[j+7*bound_points]
946  << " " << cutvalues[j+8*bound_points] << " " << cutvalues[j+9*bound_points] << endl);
947  }
948 }
TFEFunction2D::GetValues
double * GetValues()
Definition: FEFunction2D.h:67
TFEFunction2D::FindGradientLocal
void FindGradientLocal(TBaseCell *cell, int cell_no, double x, double y, double *values)
Definition: FEFunction2D.C:1188
TFE2D::GetBaseFunct2D
TBaseFunct2D * GetBaseFunct2D() const
Definition: FE2D.h:67
TFEDatabase2D::GetFE2D
static TFE2D * GetFE2D(FE2D FE)
Definition: FEDatabase2D.h:353
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
TSquareMatrix2D
Definition: SquareMatrix2D.h:20
TFE2D
Definition: FE2D.h:20
TFESpace2D
Definition: FESpace2D.h:28
TParaDB::RE_NR
double RE_NR
Definition: Database.h:313
TFEDatabase2D::SetCellForRefTrans
static void SetCellForRefTrans(TBaseCell *cell, RefTrans2D reftrans)
Definition: FEDatabase2D.C:2300
TFEDatabase2D::GetOrigValues
static void GetOrigValues(RefTrans2D RefTrans, double xi, double eta, TBaseFunct2D *bf, TCollection *Coll, TGridCell *cell, double *uref, double *uxiref, double *uetaref, double *uorig, double *uxorig, double *uyorig)
Definition: FEDatabase2D.C:2150
TMatrix::GetEntries
double * GetEntries() const
Definition: Matrix.h:87
TBaseFunct2D::GetDimension
int GetDimension() const
Definition: BaseFunct2D.h:119
TCollection
store cells in an array, used by cell iterators
Definition: Collection.h:18
TBaseCell::GetVertex
virtual TVertex * GetVertex(int Vert_i)=0
return the pointer to vertex with number i
TCollection::GetIndex
int GetIndex(TBaseCell *cell)
return Index of cell in Cells-array
Definition: Collection.C:115
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
TMatrix::GetRowPtr
int * GetRowPtr() const
Definition: Matrix.h:79
TFE2D::GetRefTransID
RefTrans2D GetRefTransID() const
Definition: FE2D.h:89
TFEFunction2D::GetFESpace2D
TFESpace2D * GetFESpace2D()
Definition: FEFunction2D.h:59
TVertex::GetCoords
void GetCoords(double &x, double &y, double &z) const
Definition: Vertex.h:106
TFESpace2D::GetFE2D
FE2D GetFE2D(int i, TBaseCell *cell)
Definition: FESpace2D.C:1184
TVertex::GetY
double GetY() const
Definition: Vertex.h:99
TBaseCell::GetN_Vertices
int GetN_Vertices()
return the number of vertices of the cell
Definition: BaseCell.h:179
TBaseCell::GetParent
virtual TBaseCell * GetParent()=0
return the parent cell
TBaseCell
information for finite element data structure
Definition: BaseCell.h:25
TBaseCell::PointInCell
virtual bool PointInCell(double X, double Y)=0
return whether a point is inside a cell
TFESpace::GetGlobalNumbers
int * GetGlobalNumbers() const
Definition: FESpace.h:135
TVertex
Definition: Vertex.h:19
TBaseCell::GetN_Edges
int GetN_Edges()
return the number of edges of the cell
Definition: BaseCell.h:182
TMatrix::GetKCol
int * GetKCol() const
Definition: Matrix.h:67
TBaseFunct2D
Definition: BaseFunct2D.h:27
TFEFunction2D::FindGradient
void FindGradient(double x, double y, double *values)
Definition: FEFunction2D.C:1078
TFEDatabase2D::GetRefFromOrig
static void GetRefFromOrig(RefTrans2D RefTrans, double X, double Y, double &xi, double &eta)
Definition: FEDatabase2D.C:2269
TVertex::GetX
double GetX() const
Definition: Vertex.h:96
TGridCell
represent geometric information of the cell
Definition: GridCell.h:15
TFESpace::GetBeginIndex
int * GetBeginIndex() const
Definition: FESpace.h:142
TDatabase::ParamDB
static TParamDB * ParamDB
Definition: Database.h:1134
TBaseFunct2D::GetDerivatives
void GetDerivatives(MultiIndex2D MultiIndex, double xi, double eta, double *values)
Definition: BaseFunct2D.h:127
TFEFunction2D
Definition: FEFunction2D.h:24
Generated on Tue May 10 2016 09:07:56 for ParMooN by   doxygen 1.8.5