MFEM  v3.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
geom.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
11 
12 #include "fem.hpp"
13 
14 namespace mfem
15 {
16 
17 const char *Geometry::Name[NumGeom] =
18 { "Point", "Segment", "Triangle", "Square", "Tetrahedron", "Cube" };
19 
20 const double Geometry::Volume[NumGeom] =
21 { 1.0, 1.0, 0.5, 1.0, 1./6, 1.0 };
22 
23 Geometry::Geometry()
24 {
25  // Vertices for Geometry::POINT
26  GeomVert[0] = NULL; // No vertices, dimension is 0
27 
28  // Vertices for Geometry::SEGMENT
29  GeomVert[1] = new IntegrationRule(2);
30 
31  GeomVert[1]->IntPoint(0).x = 0.0;
32  GeomVert[1]->IntPoint(1).x = 1.0;
33 
34  // Vertices for Geometry::TRIANGLE
35  GeomVert[2] = new IntegrationRule(3);
36 
37  GeomVert[2]->IntPoint(0).x = 0.0;
38  GeomVert[2]->IntPoint(0).y = 0.0;
39 
40  GeomVert[2]->IntPoint(1).x = 1.0;
41  GeomVert[2]->IntPoint(1).y = 0.0;
42 
43  GeomVert[2]->IntPoint(2).x = 0.0;
44  GeomVert[2]->IntPoint(2).y = 1.0;
45 
46  // Vertices for Geometry::SQUARE
47  GeomVert[3] = new IntegrationRule(4);
48 
49  GeomVert[3]->IntPoint(0).x = 0.0;
50  GeomVert[3]->IntPoint(0).y = 0.0;
51 
52  GeomVert[3]->IntPoint(1).x = 1.0;
53  GeomVert[3]->IntPoint(1).y = 0.0;
54 
55  GeomVert[3]->IntPoint(2).x = 1.0;
56  GeomVert[3]->IntPoint(2).y = 1.0;
57 
58  GeomVert[3]->IntPoint(3).x = 0.0;
59  GeomVert[3]->IntPoint(3).y = 1.0;
60 
61  // Vertices for Geometry::TETRAHEDRON
62  GeomVert[4] = new IntegrationRule(4);
63  GeomVert[4]->IntPoint(0).x = 0.0;
64  GeomVert[4]->IntPoint(0).y = 0.0;
65  GeomVert[4]->IntPoint(0).z = 0.0;
66 
67  GeomVert[4]->IntPoint(1).x = 1.0;
68  GeomVert[4]->IntPoint(1).y = 0.0;
69  GeomVert[4]->IntPoint(1).z = 0.0;
70 
71  GeomVert[4]->IntPoint(2).x = 0.0;
72  GeomVert[4]->IntPoint(2).y = 1.0;
73  GeomVert[4]->IntPoint(2).z = 0.0;
74 
75  GeomVert[4]->IntPoint(3).x = 0.0;
76  GeomVert[4]->IntPoint(3).y = 0.0;
77  GeomVert[4]->IntPoint(3).z = 1.0;
78 
79  // Vertices for Geometry::CUBE
80  GeomVert[5] = new IntegrationRule(8);
81 
82  GeomVert[5]->IntPoint(0).x = 0.0;
83  GeomVert[5]->IntPoint(0).y = 0.0;
84  GeomVert[5]->IntPoint(0).z = 0.0;
85 
86  GeomVert[5]->IntPoint(1).x = 1.0;
87  GeomVert[5]->IntPoint(1).y = 0.0;
88  GeomVert[5]->IntPoint(1).z = 0.0;
89 
90  GeomVert[5]->IntPoint(2).x = 1.0;
91  GeomVert[5]->IntPoint(2).y = 1.0;
92  GeomVert[5]->IntPoint(2).z = 0.0;
93 
94  GeomVert[5]->IntPoint(3).x = 0.0;
95  GeomVert[5]->IntPoint(3).y = 1.0;
96  GeomVert[5]->IntPoint(3).z = 0.0;
97 
98  GeomVert[5]->IntPoint(4).x = 0.0;
99  GeomVert[5]->IntPoint(4).y = 0.0;
100  GeomVert[5]->IntPoint(4).z = 1.0;
101 
102  GeomVert[5]->IntPoint(5).x = 1.0;
103  GeomVert[5]->IntPoint(5).y = 0.0;
104  GeomVert[5]->IntPoint(5).z = 1.0;
105 
106  GeomVert[5]->IntPoint(6).x = 1.0;
107  GeomVert[5]->IntPoint(6).y = 1.0;
108  GeomVert[5]->IntPoint(6).z = 1.0;
109 
110  GeomVert[5]->IntPoint(7).x = 0.0;
111  GeomVert[5]->IntPoint(7).y = 1.0;
112  GeomVert[5]->IntPoint(7).z = 1.0;
113 
114  GeomCenter[POINT].x = 0.0;
115  GeomCenter[POINT].y = 0.0;
116  GeomCenter[POINT].z = 0.0;
117 
118  GeomCenter[SEGMENT].x = 0.5;
119  GeomCenter[SEGMENT].y = 0.0;
120  GeomCenter[SEGMENT].z = 0.0;
121 
122  GeomCenter[TRIANGLE].x = 1.0 / 3.0;
123  GeomCenter[TRIANGLE].y = 1.0 / 3.0;
124  GeomCenter[TRIANGLE].z = 0.0;
125 
126  GeomCenter[SQUARE].x = 0.5;
127  GeomCenter[SQUARE].y = 0.5;
128  GeomCenter[SQUARE].z = 0.0;
129 
130  GeomCenter[TETRAHEDRON].x = 0.25;
131  GeomCenter[TETRAHEDRON].y = 0.25;
132  GeomCenter[TETRAHEDRON].z = 0.25;
133 
134  GeomCenter[CUBE].x = 0.5;
135  GeomCenter[CUBE].y = 0.5;
136  GeomCenter[CUBE].z = 0.5;
137 
138  PerfGeomToGeomJac[POINT] = NULL;
139  PerfGeomToGeomJac[SEGMENT] = NULL;
140  PerfGeomToGeomJac[TRIANGLE] = new DenseMatrix(2);
141  PerfGeomToGeomJac[SQUARE] = NULL;
142  PerfGeomToGeomJac[TETRAHEDRON] = new DenseMatrix(3);
143  PerfGeomToGeomJac[CUBE] = NULL;
144 
145  {
146  Linear2DFiniteElement TriFE;
147  IsoparametricTransformation tri_T;
148  tri_T.SetFE(&TriFE);
149  GetPerfPointMat (TRIANGLE, tri_T.GetPointMat());
150  tri_T.SetIntPoint(&GeomCenter[TRIANGLE]);
151  CalcInverse(tri_T.Jacobian(), *PerfGeomToGeomJac[TRIANGLE]);
152  }
153  {
154  Linear3DFiniteElement TetFE;
155  IsoparametricTransformation tet_T;
156  tet_T.SetFE(&TetFE);
157  GetPerfPointMat (TETRAHEDRON, tet_T.GetPointMat());
158  tet_T.SetIntPoint(&GeomCenter[TETRAHEDRON]);
159  CalcInverse(tet_T.Jacobian(), *PerfGeomToGeomJac[TETRAHEDRON]);
160  }
161 }
162 
163 Geometry::~Geometry()
164 {
165  for (int i = 0; i < NumGeom; i++)
166  {
167  delete PerfGeomToGeomJac[i];
168  delete GeomVert[i];
169  }
170 }
171 
172 const IntegrationRule * Geometry::GetVertices(int GeomType)
173 {
174  switch (GeomType)
175  {
176  case Geometry::POINT: return GeomVert[0];
177  case Geometry::SEGMENT: return GeomVert[1];
178  case Geometry::TRIANGLE: return GeomVert[2];
179  case Geometry::SQUARE: return GeomVert[3];
180  case Geometry::TETRAHEDRON: return GeomVert[4];
181  case Geometry::CUBE: return GeomVert[5];
182  default:
183  mfem_error ("Geometry::GetVertices(...)");
184  }
185  // make some compilers happy.
186  return GeomVert[0];
187 }
188 
189 // static method
190 void Geometry::GetRandomPoint(int GeomType, IntegrationPoint &ip)
191 {
192  switch (GeomType)
193  {
194  case Geometry::POINT:
195  ip.x = 0.0;
196  break;
197  case Geometry::SEGMENT:
198  ip.x = double(rand()) / RAND_MAX;
199  break;
200  case Geometry::TRIANGLE:
201  ip.x = double(rand()) / RAND_MAX;
202  ip.y = double(rand()) / RAND_MAX;
203  if (ip.x + ip.y > 1.0)
204  {
205  ip.x = 1.0 - ip.x;
206  ip.y = 1.0 - ip.y;
207  }
208  break;
209  case Geometry::SQUARE:
210  ip.x = double(rand()) / RAND_MAX;
211  ip.y = double(rand()) / RAND_MAX;
212  break;
213  case Geometry::TETRAHEDRON:
214  ip.x = double(rand()) / RAND_MAX;
215  ip.y = double(rand()) / RAND_MAX;
216  ip.z = double(rand()) / RAND_MAX;
217  // map to the triangular prism obtained by extruding the reference
218  // triangle in z direction
219  if (ip.x + ip.y > 1.0)
220  {
221  ip.x = 1.0 - ip.x;
222  ip.y = 1.0 - ip.y;
223  }
224  // split the prism into 3 parts: 1 is the reference tet, and the
225  // other two tets (as given below) are mapped to the reference tet
226  if (ip.x + ip.z > 1.0)
227  {
228  // tet with vertices: (0,0,1),(1,0,1),(0,1,1),(1,0,0)
229  ip.x = ip.x + ip.z - 1.0;
230  // ip.y = ip.y;
231  ip.z = 1.0 - ip.z;
232  // mapped to: (0,0,0),(1,0,0),(0,1,0),(0,0,1)
233  }
234  else if (ip.x + ip.y + ip.z > 1.0)
235  {
236  // tet with vertices: (0,1,1),(0,1,0),(0,0,1),(1,0,0)
237  double x = ip.x;
238  ip.x = 1.0 - x - ip.z;
239  ip.y = 1.0 - x - ip.y;
240  ip.z = x;
241  // mapped to: (0,0,0),(1,0,0),(0,1,0),(0,0,1)
242  }
243  break;
244  case Geometry::CUBE:
245  ip.x = double(rand()) / RAND_MAX;
246  ip.y = double(rand()) / RAND_MAX;
247  ip.z = double(rand()) / RAND_MAX;
248  break;
249  default:
250  MFEM_ABORT("Unknown type of reference element!");
251  }
252 }
253 
254 // static method
255 bool Geometry::CheckPoint(int GeomType, const IntegrationPoint &ip)
256 {
257  switch (GeomType)
258  {
259  case Geometry::POINT:
260  if (ip.x != 0.0) { return false; }
261  break;
262  case Geometry::SEGMENT:
263  if (ip.x < 0.0 || ip.x > 1.0) { return false; }
264  break;
265  case Geometry::TRIANGLE:
266  if (ip.x < 0.0 || ip.y < 0.0 || ip.x+ip.y > 1.0) { return false; }
267  break;
268  case Geometry::SQUARE:
269  if (ip.x < 0.0 || ip.x > 1.0 || ip.y < 0.0 || ip.y > 1.0)
270  { return false; }
271  break;
272  case Geometry::TETRAHEDRON:
273  if (ip.x < 0.0 || ip.y < 0.0 || ip.z < 0.0 ||
274  ip.x+ip.y+ip.z > 1.0) { return false; }
275  break;
276  case Geometry::CUBE:
277  if (ip.x < 0.0 || ip.x > 1.0 || ip.y < 0.0 || ip.y > 1.0 ||
278  ip.z < 0.0 || ip.z > 1.0) { return false; }
279  break;
280  default:
281  MFEM_ABORT("Unknown type of reference element!");
282  }
283  return true;
284 }
285 
286 namespace internal
287 {
288 
289 template <int N, int dim>
290 inline bool IntersectSegment(double lbeg[N], double lend[N],
291  IntegrationPoint &end)
292 {
293  double t = 1.0;
294  bool out = false;
295  for (int i = 0; i < N; i++)
296  {
297  lbeg[i] = std::max(lbeg[i], 0.0); // remove round-off
298  if (lend[i] < 0.0)
299  {
300  out = true;
301  t = std::min(t, lbeg[i]/(lbeg[i]-lend[i]));
302  }
303  }
304  if (out)
305  {
306  if (dim >= 1) { end.x = t*lend[0] + (1.0-t)*lbeg[0]; }
307  if (dim >= 2) { end.y = t*lend[1] + (1.0-t)*lbeg[1]; }
308  if (dim >= 3) { end.z = t*lend[2] + (1.0-t)*lbeg[2]; }
309  return false;
310  }
311  return true;
312 }
313 
314 }
315 
316 // static method
317 bool Geometry::ProjectPoint(int GeomType, const IntegrationPoint &beg,
318  IntegrationPoint &end)
319 {
320  switch (GeomType)
321  {
322  case Geometry::POINT:
323  {
324  if (end.x != 0.0) { end.x = 0.0; return false; }
325  break;
326  }
327  case Geometry::SEGMENT:
328  {
329  if (end.x < 0.0) { end.x = 0.0; return false; }
330  if (end.x > 1.0) { end.x = 1.0; return false; }
331  break;
332  }
333  case Geometry::TRIANGLE:
334  {
335  double lend[3] = { end.x, end.y, 1-end.x-end.y };
336  double lbeg[3] = { beg.x, beg.y, 1-beg.x-beg.y };
337  return internal::IntersectSegment<3,2>(lbeg, lend, end);
338  }
339  case Geometry::SQUARE:
340  {
341  double lend[4] = { end.x, end.y, 1-end.x, 1.0-end.y };
342  double lbeg[4] = { beg.x, beg.y, 1-beg.x, 1.0-beg.y };
343  return internal::IntersectSegment<4,2>(lbeg, lend, end);
344  }
345  case Geometry::TETRAHEDRON:
346  {
347  double lend[4] = { end.x, end.y, end.z, 1.0-end.x-end.y-end.z };
348  double lbeg[4] = { beg.x, beg.y, beg.z, 1.0-beg.x-beg.y-beg.z };
349  return internal::IntersectSegment<4,3>(lbeg, lend, end);
350  }
351  case Geometry::CUBE:
352  {
353  double lend[6] = { end.x, end.y, end.z,
354  1.0-end.x, 1.0-end.y, 1.0-end.z
355  };
356  double lbeg[6] = { beg.x, beg.y, beg.z,
357  1.0-beg.x, 1.0-beg.y, 1.0-beg.z
358  };
359  return internal::IntersectSegment<6,3>(lbeg, lend, end);
360  }
361  default:
362  MFEM_ABORT("Unknown type of reference element!");
363  }
364  return true;
365 }
366 
367 void Geometry::GetPerfPointMat(int GeomType, DenseMatrix &pm)
368 {
369  switch (GeomType)
370  {
371  case Geometry::SEGMENT:
372  {
373  pm.SetSize (1, 2);
374  pm(0,0) = 0.0;
375  pm(0,1) = 1.0;
376  }
377  break;
378 
379  case Geometry::TRIANGLE:
380  {
381  pm.SetSize (2, 3);
382  pm(0,0) = 0.0; pm(1,0) = 0.0;
383  pm(0,1) = 1.0; pm(1,1) = 0.0;
384  pm(0,2) = 0.5; pm(1,2) = 0.86602540378443864676;
385  }
386  break;
387 
388  case Geometry::SQUARE:
389  {
390  pm.SetSize (2, 4);
391  pm(0,0) = 0.0; pm(1,0) = 0.0;
392  pm(0,1) = 1.0; pm(1,1) = 0.0;
393  pm(0,2) = 1.0; pm(1,2) = 1.0;
394  pm(0,3) = 0.0; pm(1,3) = 1.0;
395  }
396  break;
397 
398  case Geometry::TETRAHEDRON:
399  {
400  pm.SetSize (3, 4);
401  pm(0,0) = 0.0; pm(1,0) = 0.0; pm(2,0) = 0.0;
402  pm(0,1) = 1.0; pm(1,1) = 0.0; pm(2,1) = 0.0;
403  pm(0,2) = 0.5; pm(1,2) = 0.86602540378443864676; pm(2,2) = 0.0;
404  pm(0,3) = 0.5; pm(1,3) = 0.28867513459481288225;
405  pm(2,3) = 0.81649658092772603273;
406  }
407  break;
408 
409  case Geometry::CUBE:
410  {
411  pm.SetSize (3, 8);
412  pm(0,0) = 0.0; pm(1,0) = 0.0; pm(2,0) = 0.0;
413  pm(0,1) = 1.0; pm(1,1) = 0.0; pm(2,1) = 0.0;
414  pm(0,2) = 1.0; pm(1,2) = 1.0; pm(2,2) = 0.0;
415  pm(0,3) = 0.0; pm(1,3) = 1.0; pm(2,3) = 0.0;
416  pm(0,4) = 0.0; pm(1,4) = 0.0; pm(2,4) = 1.0;
417  pm(0,5) = 1.0; pm(1,5) = 0.0; pm(2,5) = 1.0;
418  pm(0,6) = 1.0; pm(1,6) = 1.0; pm(2,6) = 1.0;
419  pm(0,7) = 0.0; pm(1,7) = 1.0; pm(2,7) = 1.0;
420  }
421  break;
422 
423  default:
424  mfem_error ("Geometry::GetPerfPointMat (...)");
425  }
426 }
427 
428 void Geometry::JacToPerfJac(int GeomType, const DenseMatrix &J,
429  DenseMatrix &PJ) const
430 {
431  if (PerfGeomToGeomJac[GeomType])
432  {
433  Mult(J, *PerfGeomToGeomJac[GeomType], PJ);
434  }
435  else
436  {
437  PJ = J;
438  }
439 }
440 
441 const int Geometry::NumBdrArray[] = { 0, 2, 3, 4, 4, 6 };
442 
443 Geometry Geometries;
444 
445 
446 GeometryRefiner::GeometryRefiner()
447 {
448  type = 0;
449  for (int i = 0; i < Geometry::NumGeom; i++)
450  {
451  RGeom[i] = NULL;
452  IntPts[i] = NULL;
453  }
454 }
455 
456 GeometryRefiner::~GeometryRefiner()
457 {
458  for (int i = 0; i < Geometry::NumGeom; i++)
459  {
460  delete RGeom[i];
461  delete IntPts[i];
462  }
463 }
464 
465 RefinedGeometry * GeometryRefiner::Refine (int Geom, int Times, int ETimes)
466 {
467  int i, j, k, l;
468 
469  const double *cp = NULL;
470  if (type)
471  {
472  cp = poly1d.ClosedPoints(Times);
473  }
474 
475  switch (Geom)
476  {
477  case Geometry::SEGMENT:
478  {
479  const int g = Geometry::SEGMENT;
480  if (RGeom[g] != NULL && RGeom[g]->Times == Times)
481  {
482  return RGeom[g];
483  }
484  delete RGeom[g];
485  RGeom[g] = new RefinedGeometry(Times+1, 2*Times, 0);
486  RGeom[g]->Times = Times;
487  RGeom[g]->ETimes = 0;
488  for (i = 0; i <= Times; i++)
489  {
490  IntegrationPoint &ip = RGeom[g]->RefPts.IntPoint(i);
491  ip.x = (type == 0) ? double(i) / Times : cp[i];
492  }
493  Array<int> &G = RGeom[g]->RefGeoms;
494  for (i = 0; i < Times; i++)
495  {
496  G[2*i+0] = i;
497  G[2*i+1] = i+1;
498  }
499 
500  return RGeom[g];
501  }
502 
503  case Geometry::TRIANGLE:
504  {
505  if (RGeom[2] != NULL && RGeom[2]->Times == Times &&
506  RGeom[2]->ETimes == ETimes)
507  {
508  return RGeom[2];
509  }
510 
511  if (RGeom[2] != NULL)
512  {
513  delete RGeom[2];
514  }
515  RGeom[2] = new RefinedGeometry ((Times+1)*(Times+2)/2, 3*Times*Times,
516  3*Times*(ETimes+1));
517  RGeom[2]->Times = Times;
518  RGeom[2]->ETimes = ETimes;
519  for (k = j = 0; j <= Times; j++)
520  for (i = 0; i <= Times-j; i++, k++)
521  {
522  IntegrationPoint &ip = RGeom[2]->RefPts.IntPoint(k);
523  if (type == 0)
524  {
525  ip.x = double(i) / Times;
526  ip.y = double(j) / Times;
527  }
528  else
529  {
530  ip.x = cp[i]/(cp[i] + cp[j] + cp[Times-i-j]);
531  ip.y = cp[j]/(cp[i] + cp[j] + cp[Times-i-j]);
532  }
533  }
534  Array<int> &G = RGeom[2]->RefGeoms;
535  for (l = k = j = 0; j < Times; j++, k++)
536  for (i = 0; i < Times-j; i++, k++)
537  {
538  G[l++] = k;
539  G[l++] = k+1;
540  G[l++] = k+Times-j+1;
541  if (i+j+1 < Times)
542  {
543  G[l++] = k+1;
544  G[l++] = k+Times-j+2;
545  G[l++] = k+Times-j+1;
546  }
547  }
548  Array<int> &E = RGeom[2]->RefEdges;
549  // horizontal edges
550  for (l = k = 0; k < Times; k += Times/ETimes)
551  {
552  j = k*(Times+1)-((k-1)*k)/2;
553  for (i = 0; i < Times-k; i++)
554  {
555  E[l++] = j; j++;
556  E[l++] = j;
557  }
558  }
559  // diagonal edges
560  for (k = Times; k > 0; k -= Times/ETimes)
561  {
562  j = k;
563  for (i = 0; i < k; i++)
564  {
565  E[l++] = j; j += Times-i;
566  E[l++] = j;
567  }
568  }
569  // vertical edges
570  for (k = 0; k < Times; k += Times/ETimes)
571  {
572  j = k;
573  for (i = 0; i < Times-k; i++)
574  {
575  E[l++] = j; j += Times-i+1;
576  E[l++] = j;
577  }
578  }
579 
580  return RGeom[2];
581  }
582 
583  case Geometry::SQUARE:
584  {
585  if (RGeom[3] != NULL && RGeom[3]->Times == Times &&
586  RGeom[3]->ETimes == ETimes)
587  {
588  return RGeom[3];
589  }
590 
591  if (RGeom[3] != NULL)
592  {
593  delete RGeom[3];
594  }
595  RGeom[3] = new RefinedGeometry ((Times+1)*(Times+1), 4*Times*Times,
596  4*(ETimes+1)*Times);
597  RGeom[3]->Times = Times;
598  RGeom[3]->ETimes = ETimes;
599  for (k = j = 0; j <= Times; j++)
600  for (i = 0; i <= Times; i++, k++)
601  {
602  IntegrationPoint &ip = RGeom[3]->RefPts.IntPoint(k);
603  if (type == 0)
604  {
605  ip.x = double(i) / Times;
606  ip.y = double(j) / Times;
607  }
608  else
609  {
610  ip.x = cp[i];
611  ip.y = cp[j];
612  }
613  }
614  Array<int> &G = RGeom[3]->RefGeoms;
615  for (l = k = j = 0; j < Times; j++, k++)
616  for (i = 0; i < Times; i++, k++)
617  {
618  G[l++] = k;
619  G[l++] = k+1;
620  G[l++] = k+Times+2;
621  G[l++] = k+Times+1;
622  }
623  Array<int> &E = RGeom[3]->RefEdges;
624  // horizontal edges
625  for (l = k = 0; k <= Times; k += Times/ETimes)
626  {
627  for (i = 0, j = k*(Times+1); i < Times; i++)
628  {
629  E[l++] = j; j++;
630  E[l++] = j;
631  }
632  }
633  // vertical edges (in right-to-left order)
634  for (k = Times; k >= 0; k -= Times/ETimes)
635  {
636  for (i = 0, j = k; i < Times; i++)
637  {
638  E[l++] = j; j += Times+1;
639  E[l++] = j;
640  }
641  }
642 
643  return RGeom[3];
644  }
645 
646  case Geometry::CUBE:
647  {
648  const int g = Geometry::CUBE;
649  if (RGeom[g] != NULL && RGeom[g]->Times == Times &&
650  RGeom[g]->ETimes == ETimes)
651  {
652  return RGeom[g];
653  }
654 
655  if (RGeom[g] != NULL)
656  {
657  delete RGeom[g];
658  }
659  RGeom[g] = new RefinedGeometry ((Times+1)*(Times+1)*(Times+1),
660  8*Times*Times*Times, 0);
661  RGeom[g]->Times = Times;
662  RGeom[g]->ETimes = ETimes;
663  for (l = k = 0; k <= Times; k++)
664  for (j = 0; j <= Times; j++)
665  for (i = 0; i <= Times; i++, l++)
666  {
667  IntegrationPoint &ip = RGeom[g]->RefPts.IntPoint(l);
668  if (type == 0)
669  {
670  ip.x = double(i) / Times;
671  ip.y = double(j) / Times;
672  ip.z = double(k) / Times;
673  }
674  else
675  {
676  ip.x = cp[i];
677  ip.y = cp[j];
678  ip.z = cp[k];
679  }
680  }
681  Array<int> &G = RGeom[g]->RefGeoms;
682  for (l = k = 0; k < Times; k++)
683  for (j = 0; j < Times; j++)
684  for (i = 0; i < Times; i++)
685  {
686  G[l++] = i+0 + (j+0 + (k+0) * (Times+1)) * (Times+1);
687  G[l++] = i+1 + (j+0 + (k+0) * (Times+1)) * (Times+1);
688  G[l++] = i+1 + (j+1 + (k+0) * (Times+1)) * (Times+1);
689  G[l++] = i+0 + (j+1 + (k+0) * (Times+1)) * (Times+1);
690  G[l++] = i+0 + (j+0 + (k+1) * (Times+1)) * (Times+1);
691  G[l++] = i+1 + (j+0 + (k+1) * (Times+1)) * (Times+1);
692  G[l++] = i+1 + (j+1 + (k+1) * (Times+1)) * (Times+1);
693  G[l++] = i+0 + (j+1 + (k+1) * (Times+1)) * (Times+1);
694  }
695 
696  return RGeom[g];
697  }
698 
699  case Geometry::TETRAHEDRON:
700  {
701  const int g = Geometry::TETRAHEDRON;
702  if (RGeom[g] != NULL && RGeom[g]->Times == Times &&
703  RGeom[g]->ETimes == ETimes)
704  {
705  return RGeom[g];
706  }
707 
708  if (RGeom[g] != NULL)
709  {
710  delete RGeom[g];
711  }
712 
713  // subdivide the tetrahedron with vertices
714  // (0,0,0), (0,0,1), (1,1,1), (0,1,1)
715 
716  // vertices: 0 <= i <= j <= k <= Times
717  // (3-combination with repetitions)
718  // number of vertices: (n+3)*(n+2)*(n+1)/6, n = Times
719 
720  // elements: the vertices are: v1=(i,j,k), v2=v1+u1, v3=v2+u2, v4=v3+u3
721  // where 0 <= i <= j <= k <= n-1 and
722  // u1,u2,u3 is a permutation of (1,0,0),(0,1,0),(0,0,1)
723  // such that all v2,v3,v4 have non-decreasing components
724  // number of elements: n^3
725 
726  const int n = Times;
727  RGeom[g] = new RefinedGeometry((n+3)*(n+2)*(n+1)/6, 4*n*n*n, 0);
728  // enumerate and define the vertices
729  Array<int> vi((n+1)*(n+1)*(n+1));
730  vi = -1;
731  int m = 0;
732  for (k = 0; k <= n; k++)
733  for (j = 0; j <= k; j++)
734  for (i = 0; i <= j; i++)
735  {
736  IntegrationPoint &ip = RGeom[g]->RefPts.IntPoint(m);
737  // map the coordinates to the reference tetrahedron
738  // (0,0,0) -> (0,0,0)
739  // (0,0,1) -> (1,0,0)
740  // (1,1,1) -> (0,1,0)
741  // (0,1,1) -> (0,0,1)
742  if (type == 0)
743  {
744  ip.x = double(k - j) / n;
745  ip.y = double(i) / n;
746  ip.z = double(j - i) / n;
747  }
748  else
749  {
750  double w = cp[k-j] + cp[i] + cp[j-i] + cp[Times-k];
751  ip.x = cp[k-j]/w;
752  ip.y = cp[i]/w;
753  ip.z = cp[j-i]/w;
754  }
755  l = i + (j + k * (n+1)) * (n+1);
756  vi[l] = m;
757  m++;
758  }
759  if (m != (n+3)*(n+2)*(n+1)/6)
760  {
761  mfem_error("GeometryRefiner::Refine() for TETRAHEDRON #1");
762  }
763  // elements
764  Array<int> &G = RGeom[g]->RefGeoms;
765  m = 0;
766  for (k = 0; k < n; k++)
767  for (j = 0; j <= k; j++)
768  for (i = 0; i <= j; i++)
769  {
770  // the ordering of the vertices is chosen to ensure:
771  // 1) correct orientation
772  // 2) the x,y,z edges are in the set of edges
773  // {(0,1),(2,3), (0,2),(1,3)}
774  // (goal is to ensure that subsequent refinement using
775  // this procedure preserves the six tetrahedral shapes)
776 
777  // zyx: (i,j,k)-(i,j,k+1)-(i+1,j+1,k+1)-(i,j+1,k+1)
778  G[m++] = vi[i+0 + (j+0 + (k+0) * (n+1)) * (n+1)];
779  G[m++] = vi[i+0 + (j+0 + (k+1) * (n+1)) * (n+1)];
780  G[m++] = vi[i+1 + (j+1 + (k+1) * (n+1)) * (n+1)];
781  G[m++] = vi[i+0 + (j+1 + (k+1) * (n+1)) * (n+1)];
782  if (j < k)
783  {
784  // yzx: (i,j,k)-(i+1,j+1,k+1)-(i,j+1,k)-(i,j+1,k+1)
785  G[m++] = vi[i+0 + (j+0 + (k+0) * (n+1)) * (n+1)];
786  G[m++] = vi[i+1 + (j+1 + (k+1) * (n+1)) * (n+1)];
787  G[m++] = vi[i+0 + (j+1 + (k+0) * (n+1)) * (n+1)];
788  G[m++] = vi[i+0 + (j+1 + (k+1) * (n+1)) * (n+1)];
789  // yxz: (i,j,k)-(i,j+1,k)-(i+1,j+1,k+1)-(i+1,j+1,k)
790  G[m++] = vi[i+0 + (j+0 + (k+0) * (n+1)) * (n+1)];
791  G[m++] = vi[i+0 + (j+1 + (k+0) * (n+1)) * (n+1)];
792  G[m++] = vi[i+1 + (j+1 + (k+1) * (n+1)) * (n+1)];
793  G[m++] = vi[i+1 + (j+1 + (k+0) * (n+1)) * (n+1)];
794  }
795  if (i < j)
796  {
797  // xzy: (i,j,k)-(i+1,j,k)-(i+1,j+1,k+1)-(i+1,j,k+1)
798  G[m++] = vi[i+0 + (j+0 + (k+0) * (n+1)) * (n+1)];
799  G[m++] = vi[i+1 + (j+0 + (k+0) * (n+1)) * (n+1)];
800  G[m++] = vi[i+1 + (j+1 + (k+1) * (n+1)) * (n+1)];
801  G[m++] = vi[i+1 + (j+0 + (k+1) * (n+1)) * (n+1)];
802  if (j < k)
803  {
804  // xyz: (i,j,k)-(i+1,j+1,k+1)-(i+1,j,k)-(i+1,j+1,k)
805  G[m++] = vi[i+0 + (j+0 + (k+0) * (n+1)) * (n+1)];
806  G[m++] = vi[i+1 + (j+1 + (k+1) * (n+1)) * (n+1)];
807  G[m++] = vi[i+1 + (j+0 + (k+0) * (n+1)) * (n+1)];
808  G[m++] = vi[i+1 + (j+1 + (k+0) * (n+1)) * (n+1)];
809  }
810  // zxy: (i,j,k)-(i+1,j+1,k+1)-(i,j,k+1)-(i+1,j,k+1)
811  G[m++] = vi[i+0 + (j+0 + (k+0) * (n+1)) * (n+1)];
812  G[m++] = vi[i+1 + (j+1 + (k+1) * (n+1)) * (n+1)];
813  G[m++] = vi[i+0 + (j+0 + (k+1) * (n+1)) * (n+1)];
814  G[m++] = vi[i+1 + (j+0 + (k+1) * (n+1)) * (n+1)];
815  }
816  }
817  if (m != 4*n*n*n)
818  {
819  mfem_error("GeometryRefiner::Refine() for TETRAHEDRON #2");
820  }
821  for (i = 0; i < m; i++)
822  if (G[i] < 0)
823  {
824  mfem_error("GeometryRefiner::Refine() for TETRAHEDRON #3");
825  }
826 
827  return RGeom[g];
828  }
829 
830  default:
831 
832  return RGeom[0];
833  }
834 }
835 
836 const IntegrationRule *GeometryRefiner::RefineInterior(int Geom, int Times)
837 {
838  int g = Geom;
839 
840  switch (g)
841  {
842  case Geometry::SEGMENT:
843  {
844  if (Times < 2)
845  {
846  return NULL;
847  }
848  if (IntPts[g] == NULL || IntPts[g]->GetNPoints() != Times-1)
849  {
850  delete IntPts[g];
851  IntPts[g] = new IntegrationRule(Times-1);
852  for (int i = 1; i < Times; i++)
853  {
854  IntegrationPoint &ip = IntPts[g]->IntPoint(i-1);
855  ip.x = double(i) / Times;
856  ip.y = ip.z = 0.0;
857  }
858  }
859  }
860  break;
861 
862  case Geometry::TRIANGLE:
863  {
864  if (Times < 3)
865  {
866  return NULL;
867  }
868  if (IntPts[g] == NULL ||
869  IntPts[g]->GetNPoints() != ((Times-1)*(Times-2))/2)
870  {
871  delete IntPts[g];
872  IntPts[g] = new IntegrationRule(((Times-1)*(Times-2))/2);
873  for (int k = 0, j = 1; j < Times-1; j++)
874  for (int i = 1; i < Times-j; i++, k++)
875  {
876  IntegrationPoint &ip = IntPts[g]->IntPoint(k);
877  ip.x = double(i) / Times;
878  ip.y = double(j) / Times;
879  ip.z = 0.0;
880  }
881  }
882  }
883  break;
884 
885  case Geometry::SQUARE:
886  {
887  if (Times < 2)
888  {
889  return NULL;
890  }
891  if (IntPts[g] == NULL || IntPts[g]->GetNPoints() != (Times-1)*(Times-1))
892  {
893  delete IntPts[g];
894  IntPts[g] = new IntegrationRule((Times-1)*(Times-1));
895  for (int k = 0, j = 1; j < Times; j++)
896  for (int i = 1; i < Times; i++, k++)
897  {
898  IntegrationPoint &ip = IntPts[g]->IntPoint(k);
899  ip.x = double(i) / Times;
900  ip.y = double(j) / Times;
901  ip.z = 0.0;
902  }
903  }
904  }
905  break;
906 
907  default:
908  mfem_error("GeometryRefiner::RefineInterior(...)");
909  }
910 
911  return IntPts[g];
912 }
913 
914 GeometryRefiner GlobGeometryRefiner;
915 
916 }
mfem::IntegrationRule
Class for integration rule.
Definition: intrules.hpp:83
mfem::Geometry::NumGeom
static const int NumGeom
Definition: geom.hpp:34
mfem::Geometry::JacToPerfJac
void JacToPerfJac(int GeomType, const DenseMatrix &J, DenseMatrix &PJ) const
Definition: geom.cpp:428
mfem::Geometry::GetRandomPoint
static void GetRandomPoint(int GeomType, IntegrationPoint &ip)
Definition: geom.cpp:190
mfem::Mult
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:445
mfem::ElementTransformation::SetIntPoint
void SetIntPoint(const IntegrationPoint *ip)
Definition: eltrans.hpp:33
mfem::GeometryRefiner::Refine
RefinedGeometry * Refine(int Geom, int Times, int ETimes=1)
Definition: geom.cpp:465
mfem::DenseMatrix
Data type dense matrix using column-major storage.
Definition: densemat.hpp:22
mfem::Geometry::Volume
static const double Volume[NumGeom]
Definition: geom.hpp:37
mfem::Poly_1D::ClosedPoints
const double * ClosedPoints(const int p)
Definition: fe.cpp:6523
mfem::Geometry::NumBdrArray
static const int NumBdrArray[]
Definition: geom.hpp:35
mfem::RefinedGeometry::RefEdges
Array< int > RefEdges
Definition: geom.hpp:81
mfem::Geometry::GetVertices
const IntegrationRule * GetVertices(int GeomType)
Definition: geom.cpp:172
mfem::Geometries
Geometry Geometries
Definition: geom.cpp:443
mfem::IntegrationRule::IntPoint
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:231
dim
int dim
Definition: ex3.cpp:48
mfem::GeometryRefiner::RefineInterior
const IntegrationRule * RefineInterior(int Geom, int Times)
Definition: geom.cpp:836
mfem::GlobGeometryRefiner
GeometryRefiner GlobGeometryRefiner
Definition: geom.cpp:914
mfem::Linear3DFiniteElement
Class for linear FE on tetrahedron.
Definition: fe.hpp:655
mfem::RefinedGeometry::RefPts
IntegrationRule RefPts
Definition: geom.hpp:80
mfem::Linear2DFiniteElement
Class for linear FE on triangle.
Definition: fe.hpp:374
mfem::Geometry::Name
static const char * Name[NumGeom]
Definition: geom.hpp:36
mfem::IsoparametricTransformation::Jacobian
virtual const DenseMatrix & Jacobian()
Definition: eltrans.cpp:52
mfem::CalcInverse
void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
Definition: densemat.cpp:2972
mfem::Geometry::ProjectPoint
static bool ProjectPoint(int GeomType, const IntegrationPoint &beg, IntegrationPoint &end)
Definition: geom.cpp:317
mfem::Geometry::GetPerfPointMat
void GetPerfPointMat(int GeomType, DenseMatrix &pm)
Definition: geom.cpp:367
mfem::mfem_error
void mfem_error(const char *msg)
Definition: error.cpp:23
mfem::IntegrationPoint
Class for integration point with weight.
Definition: intrules.hpp:25
mfem::IsoparametricTransformation::SetFE
void SetFE(const FiniteElement *FE)
Definition: eltrans.hpp:79
mfem::Geometry::CheckPoint
static bool CheckPoint(int GeomType, const IntegrationPoint &ip)
Check if the given point is inside the given reference element.
Definition: geom.cpp:255
mfem::DenseMatrix::SetSize
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:71
mfem::RefinedGeometry::RefGeoms
Array< int > RefGeoms
Definition: geom.hpp:81
mfem::poly1d
Poly_1D poly1d
Definition: fe.cpp:6614
Generated on Fri Feb 18 2022 18:59:17 for MFEM by   doxygen 1.8.5