MFEM  v4.1.0
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-2020, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 #include "fem.hpp"
13 #include "../mesh/wedge.hpp"
14 
15 namespace mfem
16 {
17 
18 const char *Geometry::Name[NumGeom] =
19 { "Point", "Segment", "Triangle", "Square", "Tetrahedron", "Cube", "Prism" };
20 
21 const double Geometry::Volume[NumGeom] =
22 { 1.0, 1.0, 0.5, 1.0, 1./6, 1.0, 0.5 };
23 
24 Geometry::Geometry()
25 {
26  // Vertices for Geometry::POINT
27  GeomVert[0] = new IntegrationRule(1);
28  GeomVert[0]->IntPoint(0).x = 0.0;
29 
30  // Vertices for Geometry::SEGMENT
31  GeomVert[1] = new IntegrationRule(2);
32 
33  GeomVert[1]->IntPoint(0).x = 0.0;
34  GeomVert[1]->IntPoint(1).x = 1.0;
35 
36  // Vertices for Geometry::TRIANGLE
37  GeomVert[2] = new IntegrationRule(3);
38 
39  GeomVert[2]->IntPoint(0).x = 0.0;
40  GeomVert[2]->IntPoint(0).y = 0.0;
41 
42  GeomVert[2]->IntPoint(1).x = 1.0;
43  GeomVert[2]->IntPoint(1).y = 0.0;
44 
45  GeomVert[2]->IntPoint(2).x = 0.0;
46  GeomVert[2]->IntPoint(2).y = 1.0;
47 
48  // Vertices for Geometry::SQUARE
49  GeomVert[3] = new IntegrationRule(4);
50 
51  GeomVert[3]->IntPoint(0).x = 0.0;
52  GeomVert[3]->IntPoint(0).y = 0.0;
53 
54  GeomVert[3]->IntPoint(1).x = 1.0;
55  GeomVert[3]->IntPoint(1).y = 0.0;
56 
57  GeomVert[3]->IntPoint(2).x = 1.0;
58  GeomVert[3]->IntPoint(2).y = 1.0;
59 
60  GeomVert[3]->IntPoint(3).x = 0.0;
61  GeomVert[3]->IntPoint(3).y = 1.0;
62 
63  // Vertices for Geometry::TETRAHEDRON
64  GeomVert[4] = new IntegrationRule(4);
65  GeomVert[4]->IntPoint(0).x = 0.0;
66  GeomVert[4]->IntPoint(0).y = 0.0;
67  GeomVert[4]->IntPoint(0).z = 0.0;
68 
69  GeomVert[4]->IntPoint(1).x = 1.0;
70  GeomVert[4]->IntPoint(1).y = 0.0;
71  GeomVert[4]->IntPoint(1).z = 0.0;
72 
73  GeomVert[4]->IntPoint(2).x = 0.0;
74  GeomVert[4]->IntPoint(2).y = 1.0;
75  GeomVert[4]->IntPoint(2).z = 0.0;
76 
77  GeomVert[4]->IntPoint(3).x = 0.0;
78  GeomVert[4]->IntPoint(3).y = 0.0;
79  GeomVert[4]->IntPoint(3).z = 1.0;
80 
81  // Vertices for Geometry::CUBE
82  GeomVert[5] = new IntegrationRule(8);
83 
84  GeomVert[5]->IntPoint(0).x = 0.0;
85  GeomVert[5]->IntPoint(0).y = 0.0;
86  GeomVert[5]->IntPoint(0).z = 0.0;
87 
88  GeomVert[5]->IntPoint(1).x = 1.0;
89  GeomVert[5]->IntPoint(1).y = 0.0;
90  GeomVert[5]->IntPoint(1).z = 0.0;
91 
92  GeomVert[5]->IntPoint(2).x = 1.0;
93  GeomVert[5]->IntPoint(2).y = 1.0;
94  GeomVert[5]->IntPoint(2).z = 0.0;
95 
96  GeomVert[5]->IntPoint(3).x = 0.0;
97  GeomVert[5]->IntPoint(3).y = 1.0;
98  GeomVert[5]->IntPoint(3).z = 0.0;
99 
100  GeomVert[5]->IntPoint(4).x = 0.0;
101  GeomVert[5]->IntPoint(4).y = 0.0;
102  GeomVert[5]->IntPoint(4).z = 1.0;
103 
104  GeomVert[5]->IntPoint(5).x = 1.0;
105  GeomVert[5]->IntPoint(5).y = 0.0;
106  GeomVert[5]->IntPoint(5).z = 1.0;
107 
108  GeomVert[5]->IntPoint(6).x = 1.0;
109  GeomVert[5]->IntPoint(6).y = 1.0;
110  GeomVert[5]->IntPoint(6).z = 1.0;
111 
112  GeomVert[5]->IntPoint(7).x = 0.0;
113  GeomVert[5]->IntPoint(7).y = 1.0;
114  GeomVert[5]->IntPoint(7).z = 1.0;
115 
116  // Vertices for Geometry::PRISM
117  GeomVert[6] = new IntegrationRule(6);
118  GeomVert[6]->IntPoint(0).x = 0.0;
119  GeomVert[6]->IntPoint(0).y = 0.0;
120  GeomVert[6]->IntPoint(0).z = 0.0;
121 
122  GeomVert[6]->IntPoint(1).x = 1.0;
123  GeomVert[6]->IntPoint(1).y = 0.0;
124  GeomVert[6]->IntPoint(1).z = 0.0;
125 
126  GeomVert[6]->IntPoint(2).x = 0.0;
127  GeomVert[6]->IntPoint(2).y = 1.0;
128  GeomVert[6]->IntPoint(2).z = 0.0;
129 
130  GeomVert[6]->IntPoint(3).x = 0.0;
131  GeomVert[6]->IntPoint(3).y = 0.0;
132  GeomVert[6]->IntPoint(3).z = 1.0;
133 
134  GeomVert[6]->IntPoint(4).x = 1.0;
135  GeomVert[6]->IntPoint(4).y = 0.0;
136  GeomVert[6]->IntPoint(4).z = 1.0;
137 
138  GeomVert[6]->IntPoint(5).x = 0.0;
139  GeomVert[6]->IntPoint(5).y = 1.0;
140  GeomVert[6]->IntPoint(5).z = 1.0;
141 
142  GeomCenter[POINT].x = 0.0;
143  GeomCenter[POINT].y = 0.0;
144  GeomCenter[POINT].z = 0.0;
145 
146  GeomCenter[SEGMENT].x = 0.5;
147  GeomCenter[SEGMENT].y = 0.0;
148  GeomCenter[SEGMENT].z = 0.0;
149 
150  GeomCenter[TRIANGLE].x = 1.0 / 3.0;
151  GeomCenter[TRIANGLE].y = 1.0 / 3.0;
152  GeomCenter[TRIANGLE].z = 0.0;
153 
154  GeomCenter[SQUARE].x = 0.5;
155  GeomCenter[SQUARE].y = 0.5;
156  GeomCenter[SQUARE].z = 0.0;
157 
158  GeomCenter[TETRAHEDRON].x = 0.25;
159  GeomCenter[TETRAHEDRON].y = 0.25;
160  GeomCenter[TETRAHEDRON].z = 0.25;
161 
162  GeomCenter[CUBE].x = 0.5;
163  GeomCenter[CUBE].y = 0.5;
164  GeomCenter[CUBE].z = 0.5;
165 
166  GeomCenter[PRISM].x = 1.0 / 3.0;
167  GeomCenter[PRISM].y = 1.0 / 3.0;
168  GeomCenter[PRISM].z = 0.5;
169 
170  GeomToPerfGeomJac[POINT] = NULL;
171  GeomToPerfGeomJac[SEGMENT] = new DenseMatrix(1);
172  GeomToPerfGeomJac[TRIANGLE] = new DenseMatrix(2);
173  GeomToPerfGeomJac[SQUARE] = new DenseMatrix(2);
174  GeomToPerfGeomJac[TETRAHEDRON] = new DenseMatrix(3);
175  GeomToPerfGeomJac[CUBE] = new DenseMatrix(3);
176  GeomToPerfGeomJac[PRISM] = new DenseMatrix(3);
177 
178  PerfGeomToGeomJac[POINT] = NULL;
179  PerfGeomToGeomJac[SEGMENT] = NULL;
180  PerfGeomToGeomJac[TRIANGLE] = new DenseMatrix(2);
181  PerfGeomToGeomJac[SQUARE] = NULL;
182  PerfGeomToGeomJac[TETRAHEDRON] = new DenseMatrix(3);
183  PerfGeomToGeomJac[CUBE] = NULL;
184  PerfGeomToGeomJac[PRISM] = new DenseMatrix(3);
185 
186  GeomToPerfGeomJac[SEGMENT]->Diag(1.0, 1);
187  {
188  IsoparametricTransformation tri_T;
189  tri_T.SetFE(&TriangleFE);
190  GetPerfPointMat (TRIANGLE, tri_T.GetPointMat());
191  tri_T.FinalizeTransformation();
192  tri_T.SetIntPoint(&GeomCenter[TRIANGLE]);
193  *GeomToPerfGeomJac[TRIANGLE] = tri_T.Jacobian();
194  CalcInverse(tri_T.Jacobian(), *PerfGeomToGeomJac[TRIANGLE]);
195  }
196  GeomToPerfGeomJac[SQUARE]->Diag(1.0, 2);
197  {
198  IsoparametricTransformation tet_T;
199  tet_T.SetFE(&TetrahedronFE);
200  GetPerfPointMat (TETRAHEDRON, tet_T.GetPointMat());
201  tet_T.FinalizeTransformation();
202  tet_T.SetIntPoint(&GeomCenter[TETRAHEDRON]);
203  *GeomToPerfGeomJac[TETRAHEDRON] = tet_T.Jacobian();
204  CalcInverse(tet_T.Jacobian(), *PerfGeomToGeomJac[TETRAHEDRON]);
205  }
206  GeomToPerfGeomJac[CUBE]->Diag(1.0, 3);
207  {
208  IsoparametricTransformation pri_T;
209  pri_T.SetFE(&WedgeFE);
210  GetPerfPointMat (PRISM, pri_T.GetPointMat());
211  pri_T.FinalizeTransformation();
212  pri_T.SetIntPoint(&GeomCenter[PRISM]);
213  *GeomToPerfGeomJac[PRISM] = pri_T.Jacobian();
214  CalcInverse(pri_T.Jacobian(), *PerfGeomToGeomJac[PRISM]);
215  }
216 }
217 
218 Geometry::~Geometry()
219 {
220  for (int i = 0; i < NumGeom; i++)
221  {
222  delete PerfGeomToGeomJac[i];
223  delete GeomToPerfGeomJac[i];
224  delete GeomVert[i];
225  }
226 }
227 
228 const IntegrationRule * Geometry::GetVertices(int GeomType)
229 {
230  switch (GeomType)
231  {
232  case Geometry::POINT: return GeomVert[0];
233  case Geometry::SEGMENT: return GeomVert[1];
234  case Geometry::TRIANGLE: return GeomVert[2];
235  case Geometry::SQUARE: return GeomVert[3];
236  case Geometry::TETRAHEDRON: return GeomVert[4];
237  case Geometry::CUBE: return GeomVert[5];
238  case Geometry::PRISM: return GeomVert[6];
239  default:
240  mfem_error ("Geometry::GetVertices(...)");
241  }
242  // make some compilers happy.
243  return GeomVert[0];
244 }
245 
246 // static method
247 void Geometry::GetRandomPoint(int GeomType, IntegrationPoint &ip)
248 {
249  switch (GeomType)
250  {
251  case Geometry::POINT:
252  ip.x = 0.0;
253  break;
254  case Geometry::SEGMENT:
255  ip.x = double(rand()) / RAND_MAX;
256  break;
257  case Geometry::TRIANGLE:
258  ip.x = double(rand()) / RAND_MAX;
259  ip.y = double(rand()) / RAND_MAX;
260  if (ip.x + ip.y > 1.0)
261  {
262  ip.x = 1.0 - ip.x;
263  ip.y = 1.0 - ip.y;
264  }
265  break;
266  case Geometry::SQUARE:
267  ip.x = double(rand()) / RAND_MAX;
268  ip.y = double(rand()) / RAND_MAX;
269  break;
270  case Geometry::TETRAHEDRON:
271  ip.x = double(rand()) / RAND_MAX;
272  ip.y = double(rand()) / RAND_MAX;
273  ip.z = double(rand()) / RAND_MAX;
274  // map to the triangular prism obtained by extruding the reference
275  // triangle in z direction
276  if (ip.x + ip.y > 1.0)
277  {
278  ip.x = 1.0 - ip.x;
279  ip.y = 1.0 - ip.y;
280  }
281  // split the prism into 3 parts: 1 is the reference tet, and the
282  // other two tets (as given below) are mapped to the reference tet
283  if (ip.x + ip.z > 1.0)
284  {
285  // tet with vertices: (0,0,1),(1,0,1),(0,1,1),(1,0,0)
286  ip.x = ip.x + ip.z - 1.0;
287  // ip.y = ip.y;
288  ip.z = 1.0 - ip.z;
289  // mapped to: (0,0,0),(1,0,0),(0,1,0),(0,0,1)
290  }
291  else if (ip.x + ip.y + ip.z > 1.0)
292  {
293  // tet with vertices: (0,1,1),(0,1,0),(0,0,1),(1,0,0)
294  double x = ip.x;
295  ip.x = 1.0 - x - ip.z;
296  ip.y = 1.0 - x - ip.y;
297  ip.z = x;
298  // mapped to: (0,0,0),(1,0,0),(0,1,0),(0,0,1)
299  }
300  break;
301  case Geometry::CUBE:
302  ip.x = double(rand()) / RAND_MAX;
303  ip.y = double(rand()) / RAND_MAX;
304  ip.z = double(rand()) / RAND_MAX;
305  break;
306  case Geometry::PRISM:
307  ip.x = double(rand()) / RAND_MAX;
308  ip.y = double(rand()) / RAND_MAX;
309  ip.z = double(rand()) / RAND_MAX;
310  if (ip.x + ip.y > 1.0)
311  {
312  ip.x = 1.0 - ip.x;
313  ip.y = 1.0 - ip.y;
314  }
315  break;
316  default:
317  MFEM_ABORT("Unknown type of reference element!");
318  }
319 }
320 
321 
322 namespace internal
323 {
324 
325 // Fuzzy equality operator with absolute tolerance eps.
326 inline bool NearlyEqual(double x, double y, double eps)
327 {
328  return std::abs(x-y) <= eps;
329 }
330 
331 // Fuzzy greater than comparison operator with absolute tolerance eps.
332 // Returns true when x is greater than y by at least eps.
333 inline bool FuzzyGT(double x, double y, double eps)
334 {
335  return (x > y + eps);
336 }
337 
338 // Fuzzy less than comparison operator with absolute tolerance eps.
339 // Returns true when x is less than y by at least eps.
340 inline bool FuzzyLT(double x, double y, double eps)
341 {
342  return (x < y - eps);
343 }
344 
345 }
346 
347 // static method
348 bool Geometry::CheckPoint(int GeomType, const IntegrationPoint &ip)
349 {
350  switch (GeomType)
351  {
352  case Geometry::POINT:
353  if (ip.x != 0.0) { return false; }
354  break;
355  case Geometry::SEGMENT:
356  if (ip.x < 0.0 || ip.x > 1.0) { return false; }
357  break;
358  case Geometry::TRIANGLE:
359  if (ip.x < 0.0 || ip.y < 0.0 || ip.x+ip.y > 1.0) { return false; }
360  break;
361  case Geometry::SQUARE:
362  if (ip.x < 0.0 || ip.x > 1.0 || ip.y < 0.0 || ip.y > 1.0)
363  { return false; }
364  break;
365  case Geometry::TETRAHEDRON:
366  if (ip.x < 0.0 || ip.y < 0.0 || ip.z < 0.0 ||
367  ip.x+ip.y+ip.z > 1.0) { return false; }
368  break;
369  case Geometry::CUBE:
370  if (ip.x < 0.0 || ip.x > 1.0 || ip.y < 0.0 || ip.y > 1.0 ||
371  ip.z < 0.0 || ip.z > 1.0) { return false; }
372  break;
373  case Geometry::PRISM:
374  if (ip.x < 0.0 || ip.y < 0.0 || ip.x+ip.y > 1.0 ||
375  ip.z < 0.0 || ip.z > 1.0) { return false; }
376  break;
377  default:
378  MFEM_ABORT("Unknown type of reference element!");
379  }
380  return true;
381 }
382 
383 bool Geometry::CheckPoint(int GeomType, const IntegrationPoint &ip, double eps)
384 {
385  switch (GeomType)
386  {
387  case Geometry::POINT:
388  if (! internal::NearlyEqual(ip.x, 0.0, eps))
389  {
390  return false;
391  }
392  break;
393  case Geometry::SEGMENT:
394  if ( internal::FuzzyLT(ip.x, 0.0, eps)
395  || internal::FuzzyGT(ip.x, 1.0, eps) )
396  {
397  return false;
398  }
399  break;
400  case Geometry::TRIANGLE:
401  if ( internal::FuzzyLT(ip.x, 0.0, eps)
402  || internal::FuzzyLT(ip.y, 0.0, eps)
403  || internal::FuzzyGT(ip.x+ip.y, 1.0, eps) )
404  {
405  return false;
406  }
407  break;
408  case Geometry::SQUARE:
409  if ( internal::FuzzyLT(ip.x, 0.0, eps)
410  || internal::FuzzyGT(ip.x, 1.0, eps)
411  || internal::FuzzyLT(ip.y, 0.0, eps)
412  || internal::FuzzyGT(ip.y, 1.0, eps) )
413  {
414  return false;
415  }
416  break;
417  case Geometry::TETRAHEDRON:
418  if ( internal::FuzzyLT(ip.x, 0.0, eps)
419  || internal::FuzzyLT(ip.y, 0.0, eps)
420  || internal::FuzzyLT(ip.z, 0.0, eps)
421  || internal::FuzzyGT(ip.x+ip.y+ip.z, 1.0, eps) )
422  {
423  return false;
424  }
425  break;
426  case Geometry::CUBE:
427  if ( internal::FuzzyLT(ip.x, 0.0, eps)
428  || internal::FuzzyGT(ip.x, 1.0, eps)
429  || internal::FuzzyLT(ip.y, 0.0, eps)
430  || internal::FuzzyGT(ip.y, 1.0, eps)
431  || internal::FuzzyLT(ip.z, 0.0, eps)
432  || internal::FuzzyGT(ip.z, 1.0, eps) )
433  {
434  return false;
435  }
436  break;
437  case Geometry::PRISM:
438  if ( internal::FuzzyLT(ip.x, 0.0, eps)
439  || internal::FuzzyLT(ip.y, 0.0, eps)
440  || internal::FuzzyGT(ip.x+ip.y, 1.0, eps)
441  || internal::FuzzyLT(ip.z, 0.0, eps)
442  || internal::FuzzyGT(ip.z, 1.0, eps) )
443  {
444  return false;
445  }
446  break;
447  default:
448  MFEM_ABORT("Unknown type of reference element!");
449  }
450  return true;
451 }
452 
453 
454 namespace internal
455 {
456 
457 template <int N, int dim>
458 inline bool IntersectSegment(double lbeg[N], double lend[N],
459  IntegrationPoint &end)
460 {
461  double t = 1.0;
462  bool out = false;
463  for (int i = 0; i < N; i++)
464  {
465  lbeg[i] = std::max(lbeg[i], 0.0); // remove round-off
466  if (lend[i] < 0.0)
467  {
468  out = true;
469  t = std::min(t, lbeg[i]/(lbeg[i]-lend[i]));
470  }
471  }
472  if (out)
473  {
474  if (dim >= 1) { end.x = t*lend[0] + (1.0-t)*lbeg[0]; }
475  if (dim >= 2) { end.y = t*lend[1] + (1.0-t)*lbeg[1]; }
476  if (dim >= 3) { end.z = t*lend[2] + (1.0-t)*lbeg[2]; }
477  return false;
478  }
479  return true;
480 }
481 
482 inline bool ProjectTriangle(double &x, double &y)
483 {
484  if (x < 0.0)
485  {
486  x = 0.0;
487  if (y < 0.0) { y = 0.0; }
488  else if (y > 1.0) { y = 1.0; }
489  return false;
490  }
491  if (y < 0.0)
492  {
493  if (x > 1.0) { x = 1.0; }
494  y = 0.0;
495  return false;
496  }
497  const double l3 = 1.0-x-y;
498  if (l3 < 0.0)
499  {
500  if (y - x > 1.0) { x = 0.0; y = 1.0; }
501  else if (y - x < -1.0) { x = 1.0; y = 0.0; }
502  else { x += l3/2; y += l3/2; }
503  return false;
504  }
505  return true;
506 }
507 
508 }
509 
510 // static method
511 bool Geometry::ProjectPoint(int GeomType, const IntegrationPoint &beg,
512  IntegrationPoint &end)
513 {
514  switch (GeomType)
515  {
516  case Geometry::POINT:
517  {
518  if (end.x != 0.0) { end.x = 0.0; return false; }
519  break;
520  }
521  case Geometry::SEGMENT:
522  {
523  if (end.x < 0.0) { end.x = 0.0; return false; }
524  if (end.x > 1.0) { end.x = 1.0; return false; }
525  break;
526  }
527  case Geometry::TRIANGLE:
528  {
529  double lend[3] = { end.x, end.y, 1.0-end.x-end.y };
530  double lbeg[3] = { beg.x, beg.y, 1.0-beg.x-beg.y };
531  return internal::IntersectSegment<3,2>(lbeg, lend, end);
532  }
533  case Geometry::SQUARE:
534  {
535  double lend[4] = { end.x, end.y, 1.0-end.x, 1.0-end.y };
536  double lbeg[4] = { beg.x, beg.y, 1.0-beg.x, 1.0-beg.y };
537  return internal::IntersectSegment<4,2>(lbeg, lend, end);
538  }
539  case Geometry::TETRAHEDRON:
540  {
541  double lend[4] = { end.x, end.y, end.z, 1.0-end.x-end.y-end.z };
542  double lbeg[4] = { beg.x, beg.y, beg.z, 1.0-beg.x-beg.y-beg.z };
543  return internal::IntersectSegment<4,3>(lbeg, lend, end);
544  }
545  case Geometry::CUBE:
546  {
547  double lend[6] = { end.x, end.y, end.z,
548  1.0-end.x, 1.0-end.y, 1.0-end.z
549  };
550  double lbeg[6] = { beg.x, beg.y, beg.z,
551  1.0-beg.x, 1.0-beg.y, 1.0-beg.z
552  };
553  return internal::IntersectSegment<6,3>(lbeg, lend, end);
554  }
555  case Geometry::PRISM:
556  {
557  double lend[5] = { end.x, end.y, end.z, 1.0-end.x-end.y, 1.0-end.z };
558  double lbeg[5] = { beg.x, beg.y, beg.z, 1.0-beg.x-beg.y, 1.0-beg.z };
559  return internal::IntersectSegment<5,3>(lbeg, lend, end);
560  }
561  default:
562  MFEM_ABORT("Unknown type of reference element!");
563  }
564  return true;
565 }
566 
567 // static method
568 bool Geometry::ProjectPoint(int GeomType, IntegrationPoint &ip)
569 {
570  // If ip is outside the element, replace it with the point on the boundary
571  // that is closest to the original ip and return false; otherwise, return
572  // true without changing ip.
573 
574  switch (GeomType)
575  {
576  case SEGMENT:
577  {
578  if (ip.x < 0.0) { ip.x = 0.0; return false; }
579  else if (ip.x > 1.0) { ip.x = 1.0; return false; }
580  return true;
581  }
582 
583  case TRIANGLE:
584  {
585  return internal::ProjectTriangle(ip.x, ip.y);
586  }
587 
588  case SQUARE:
589  {
590  bool in_x, in_y;
591  if (ip.x < 0.0) { in_x = false; ip.x = 0.0; }
592  else if (ip.x > 1.0) { in_x = false; ip.x = 1.0; }
593  else { in_x = true; }
594  if (ip.y < 0.0) { in_y = false; ip.y = 0.0; }
595  else if (ip.y > 1.0) { in_y = false; ip.y = 1.0; }
596  else { in_y = true; }
597  return in_x && in_y;
598  }
599 
600  case TETRAHEDRON:
601  {
602  if (ip.z < 0.0)
603  {
604  ip.z = 0.0;
605  internal::ProjectTriangle(ip.x, ip.y);
606  return false;
607  }
608  if (ip.y < 0.0)
609  {
610  ip.y = 0.0;
611  internal::ProjectTriangle(ip.x, ip.z);
612  return false;
613  }
614  if (ip.x < 0.0)
615  {
616  ip.x = 0.0;
617  internal::ProjectTriangle(ip.y, ip.z);
618  return false;
619  }
620  const double l4 = 1.0-ip.x-ip.y-ip.z;
621  if (l4 < 0.0)
622  {
623  const double l4_3 = l4/3;
624  ip.x += l4_3;
625  ip.y += l4_3;
626  internal::ProjectTriangle(ip.x, ip.y);
627  ip.z = 1.0-ip.x-ip.y;
628  return false;
629  }
630  return true;
631  }
632 
633  case CUBE:
634  {
635  bool in_x, in_y, in_z;
636  if (ip.x < 0.0) { in_x = false; ip.x = 0.0; }
637  else if (ip.x > 1.0) { in_x = false; ip.x = 1.0; }
638  else { in_x = true; }
639  if (ip.y < 0.0) { in_y = false; ip.y = 0.0; }
640  else if (ip.y > 1.0) { in_y = false; ip.y = 1.0; }
641  else { in_y = true; }
642  if (ip.z < 0.0) { in_z = false; ip.z = 0.0; }
643  else if (ip.z > 1.0) { in_z = false; ip.z = 1.0; }
644  else { in_z = true; }
645  return in_x && in_y && in_z;
646  }
647 
648  case PRISM:
649  {
650  bool in_tri, in_z;
651  in_tri = internal::ProjectTriangle(ip.x, ip.y);
652  if (ip.z < 0.0) { in_z = false; ip.z = 0.0; }
653  else if (ip.z > 1.0) { in_z = false; ip.z = 1.0; }
654  else { in_z = true; }
655  return in_tri && in_z;
656  }
657 
658  default:
659  MFEM_ABORT("Reference element type is not supported!");
660  }
661  return true;
662 }
663 
664 void Geometry::GetPerfPointMat(int GeomType, DenseMatrix &pm)
665 {
666  switch (GeomType)
667  {
668  case Geometry::SEGMENT:
669  {
670  pm.SetSize (1, 2);
671  pm(0,0) = 0.0;
672  pm(0,1) = 1.0;
673  }
674  break;
675 
676  case Geometry::TRIANGLE:
677  {
678  pm.SetSize (2, 3);
679  pm(0,0) = 0.0; pm(1,0) = 0.0;
680  pm(0,1) = 1.0; pm(1,1) = 0.0;
681  pm(0,2) = 0.5; pm(1,2) = 0.86602540378443864676;
682  }
683  break;
684 
685  case Geometry::SQUARE:
686  {
687  pm.SetSize (2, 4);
688  pm(0,0) = 0.0; pm(1,0) = 0.0;
689  pm(0,1) = 1.0; pm(1,1) = 0.0;
690  pm(0,2) = 1.0; pm(1,2) = 1.0;
691  pm(0,3) = 0.0; pm(1,3) = 1.0;
692  }
693  break;
694 
695  case Geometry::TETRAHEDRON:
696  {
697  pm.SetSize (3, 4);
698  pm(0,0) = 0.0; pm(1,0) = 0.0; pm(2,0) = 0.0;
699  pm(0,1) = 1.0; pm(1,1) = 0.0; pm(2,1) = 0.0;
700  pm(0,2) = 0.5; pm(1,2) = 0.86602540378443864676; pm(2,2) = 0.0;
701  pm(0,3) = 0.5; pm(1,3) = 0.28867513459481288225;
702  pm(2,3) = 0.81649658092772603273;
703  }
704  break;
705 
706  case Geometry::CUBE:
707  {
708  pm.SetSize (3, 8);
709  pm(0,0) = 0.0; pm(1,0) = 0.0; pm(2,0) = 0.0;
710  pm(0,1) = 1.0; pm(1,1) = 0.0; pm(2,1) = 0.0;
711  pm(0,2) = 1.0; pm(1,2) = 1.0; pm(2,2) = 0.0;
712  pm(0,3) = 0.0; pm(1,3) = 1.0; pm(2,3) = 0.0;
713  pm(0,4) = 0.0; pm(1,4) = 0.0; pm(2,4) = 1.0;
714  pm(0,5) = 1.0; pm(1,5) = 0.0; pm(2,5) = 1.0;
715  pm(0,6) = 1.0; pm(1,6) = 1.0; pm(2,6) = 1.0;
716  pm(0,7) = 0.0; pm(1,7) = 1.0; pm(2,7) = 1.0;
717  }
718  break;
719 
720  case Geometry::PRISM:
721  {
722  pm.SetSize (3, 6);
723  pm(0,0) = 0.0; pm(1,0) = 0.0; pm(2,0) = 0.0;
724  pm(0,1) = 1.0; pm(1,1) = 0.0; pm(2,1) = 0.0;
725  pm(0,2) = 0.5; pm(1,2) = 0.86602540378443864676; pm(2,2) = 0.0;
726  pm(0,3) = 0.0; pm(1,3) = 0.0; pm(2,3) = 1.0;
727  pm(0,4) = 1.0; pm(1,4) = 0.0; pm(2,4) = 1.0;
728  pm(0,5) = 0.5; pm(1,5) = 0.86602540378443864676; pm(2,5) = 1.0;
729  }
730  break;
731 
732  default:
733  mfem_error ("Geometry::GetPerfPointMat (...)");
734  }
735 }
736 
737 void Geometry::JacToPerfJac(int GeomType, const DenseMatrix &J,
738  DenseMatrix &PJ) const
739 {
740  if (PerfGeomToGeomJac[GeomType])
741  {
742  Mult(J, *PerfGeomToGeomJac[GeomType], PJ);
743  }
744  else
745  {
746  PJ = J;
747  }
748 }
749 
750 const int Geometry::NumBdrArray[NumGeom] = { 0, 2, 3, 4, 4, 6, 5 };
751 const int Geometry::Dimension[NumGeom] = { 0, 1, 2, 2, 3, 3, 3 };
752 const int Geometry::DimStart[MaxDim+2] =
753 { POINT, SEGMENT, TRIANGLE, TETRAHEDRON, NUM_GEOMETRIES };
754 const int Geometry::NumVerts[NumGeom] = { 1, 2, 3, 4, 4, 8, 6 };
755 const int Geometry::NumEdges[NumGeom] = { 0, 1, 3, 4, 6, 12, 9 };
756 const int Geometry::NumFaces[NumGeom] = { 0, 0, 1, 1, 4, 6, 5 };
757 
758 const int Geometry::
759 Constants<Geometry::POINT>::Orient[1][1] = {{0}};
760 const int Geometry::
761 Constants<Geometry::POINT>::InvOrient[1] = {0};
762 
763 const int Geometry::
764 Constants<Geometry::SEGMENT>::Edges[1][2] = { {0, 1} };
765 const int Geometry::
766 Constants<Geometry::SEGMENT>::Orient[2][2] = { {0, 1}, {1, 0} };
767 const int Geometry::
768 Constants<Geometry::SEGMENT>::InvOrient[2] = { 0, 1 };
769 
770 const int Geometry::
771 Constants<Geometry::TRIANGLE>::Edges[3][2] = {{0, 1}, {1, 2}, {2, 0}};
772 const int Geometry::
773 Constants<Geometry::TRIANGLE>::VertToVert::I[3] = {0, 2, 3};
774 const int Geometry::
775 Constants<Geometry::TRIANGLE>::VertToVert::J[3][2] = {{1, 0}, {2, -3}, {2, 1}};
776 const int Geometry::
777 Constants<Geometry::TRIANGLE>::FaceVert[1][3] = {{0, 1, 2}};
778 const int Geometry::
779 Constants<Geometry::TRIANGLE>::Orient[6][3] =
780 {
781  {0, 1, 2}, {1, 0, 2}, {2, 0, 1},
782  {2, 1, 0}, {1, 2, 0}, {0, 2, 1}
783 };
784 const int Geometry::
785 Constants<Geometry::TRIANGLE>::InvOrient[6] = { 0, 1, 4, 3, 2, 5 };
786 
787 const int Geometry::
788 Constants<Geometry::SQUARE>::Edges[4][2] = {{0, 1}, {1, 2}, {2, 3}, {3, 0}};
789 const int Geometry::
790 Constants<Geometry::SQUARE>::VertToVert::I[4] = {0, 2, 3, 4};
791 const int Geometry::
792 Constants<Geometry::SQUARE>::VertToVert::J[4][2] =
793 {{1, 0}, {3, -4}, {2, 1}, {3, 2}};
794 const int Geometry::
795 Constants<Geometry::SQUARE>::FaceVert[1][4] = {{0, 1, 2, 3}};
796 const int Geometry::
797 Constants<Geometry::SQUARE>::Orient[8][4] =
798 {
799  {0, 1, 2, 3}, {0, 3, 2, 1}, {1, 2, 3, 0}, {1, 0, 3, 2},
800  {2, 3, 0, 1}, {2, 1, 0, 3}, {3, 0, 1, 2}, {3, 2, 1, 0}
801 };
802 const int Geometry::
803 Constants<Geometry::SQUARE>::InvOrient[8] = { 0, 1, 6, 3, 4, 5, 2, 7 };
804 
805 const int Geometry::
806 Constants<Geometry::TETRAHEDRON>::Edges[6][2] =
807 {{0, 1}, {0, 2}, {0, 3}, {1, 2}, {1, 3}, {2, 3}};
808 const int Geometry::
809 Constants<Geometry::TETRAHEDRON>::FaceTypes[4] =
810 {
811  Geometry::TRIANGLE, Geometry::TRIANGLE,
812  Geometry::TRIANGLE, Geometry::TRIANGLE
813 };
814 const int Geometry::
815 Constants<Geometry::TETRAHEDRON>::FaceVert[4][3] =
816 {{1, 2, 3}, {0, 3, 2}, {0, 1, 3}, {0, 2, 1}};
817 const int Geometry::
818 Constants<Geometry::TETRAHEDRON>::VertToVert::I[4] = {0, 3, 5, 6};
819 const int Geometry::
820 Constants<Geometry::TETRAHEDRON>::VertToVert::J[6][2] =
821 {
822  {1, 0}, {2, 1}, {3, 2}, // 0,1:0 0,2:1 0,3:2
823  {2, 3}, {3, 4}, // 1,2:3 1,3:4
824  {3, 5} // 2,3:5
825 };
826 
827 const int Geometry::
828 Constants<Geometry::CUBE>::Edges[12][2] =
829 {
830  {0, 1}, {1, 2}, {3, 2}, {0, 3}, {4, 5}, {5, 6},
831  {7, 6}, {4, 7}, {0, 4}, {1, 5}, {2, 6}, {3, 7}
832 };
833 const int Geometry::
834 Constants<Geometry::CUBE>::FaceTypes[6] =
835 {
836  Geometry::SQUARE, Geometry::SQUARE, Geometry::SQUARE,
837  Geometry::SQUARE, Geometry::SQUARE, Geometry::SQUARE
838 };
839 const int Geometry::
840 Constants<Geometry::CUBE>::FaceVert[6][4] =
841 {
842  {3, 2, 1, 0}, {0, 1, 5, 4}, {1, 2, 6, 5},
843  {2, 3, 7, 6}, {3, 0, 4, 7}, {4, 5, 6, 7}
844 };
845 const int Geometry::
846 Constants<Geometry::CUBE>::VertToVert::I[8] = {0, 3, 5, 7, 8, 10, 11, 12};
847 const int Geometry::
848 Constants<Geometry::CUBE>::VertToVert::J[12][2] =
849 {
850  {1, 0}, {3, 3}, {4, 8}, // 0,1:0 0,3:3 0,4:8
851  {2, 1}, {5, 9}, // 1,2:1 1,5:9
852  {3,-3}, {6,10}, // 2,3:-3 2,6:10
853  {7,11}, // 3,7:11
854  {5, 4}, {7, 7}, // 4,5:4 4,7:7
855  {6, 5}, // 5,6:5
856  {7,-7} // 6,7:-7
857 };
858 
859 const int Geometry::
860 Constants<Geometry::PRISM>::Edges[9][2] =
861 {{0, 1}, {1, 2}, {2, 0}, {3, 4}, {4, 5}, {5, 3}, {0, 3}, {1, 4}, {2, 5}};
862 const int Geometry::
863 Constants<Geometry::PRISM>::FaceTypes[5] =
864 {
865  Geometry::TRIANGLE, Geometry::TRIANGLE,
866  Geometry::SQUARE, Geometry::SQUARE, Geometry::SQUARE
867 };
868 const int Geometry::
869 Constants<Geometry::PRISM>::FaceVert[5][4] =
870 {{0, 2, 1, -1}, {3, 4, 5, -1}, {0, 1, 4, 3}, {1, 2, 5, 4}, {2, 0, 3, 5}};
871 const int Geometry::
872 Constants<Geometry::PRISM>::VertToVert::I[6] = {0, 3, 5, 6, 8, 9};
873 const int Geometry::
874 Constants<Geometry::PRISM>::VertToVert::J[9][2] =
875 {
876  {1, 0}, {2, -3}, {3, 6}, // 0,1:0 0,2:-3 0,3:6
877  {2, 1}, {4, 7}, // 1,2:1 1,4:7
878  {5, 8}, // 2,5:8
879  {4, 3}, {5, -6}, // 3,4:3 3,5:-6
880  {5, 4} // 4,5:4
881 };
882 
883 
884 GeometryRefiner::GeometryRefiner()
885 {
886  type = Quadrature1D::ClosedUniform;
887 }
888 
889 GeometryRefiner::~GeometryRefiner()
890 {
891  for (int i = 0; i < Geometry::NumGeom; i++)
892  {
893  for (int j = 0; j < RGeom[i].Size(); j++) { delete RGeom[i][j]; }
894  for (int j = 0; j < IntPts[i].Size(); j++) { delete IntPts[i][j]; }
895  }
896 }
897 
898 RefinedGeometry *GeometryRefiner::FindInRGeom(Geometry::Type Geom,
899  int Times, int ETimes,
900  int Type)
901 {
902  Array<RefinedGeometry *> &RGA = RGeom[Geom];
903  for (int i = 0; i < RGA.Size(); i++)
904  {
905  RefinedGeometry &RG = *RGA[i];
906  if (RG.Times == Times && RG.ETimes == ETimes && RG.Type == Type)
907  {
908  return &RG;
909  }
910  }
911  return NULL;
912 }
913 
914 IntegrationRule *GeometryRefiner::FindInIntPts(Geometry::Type Geom, int NPts)
915 {
916  Array<IntegrationRule *> &IPA = IntPts[Geom];
917  for (int i = 0; i < IPA.Size(); i++)
918  {
919  IntegrationRule &ir = *IPA[i];
920  if (ir.GetNPoints() == NPts) { return &ir; }
921  }
922  return NULL;
923 }
924 
925 RefinedGeometry * GeometryRefiner::Refine(Geometry::Type Geom,
926  int Times, int ETimes)
927 {
928  int i, j, k, l, m;
929 
930  Times = std::max(Times, 1);
931  ETimes = std::max(ETimes, 1);
932  const double *cp = poly1d.GetPoints(Times, BasisType::GetNodalBasis(type));
933 
934  RefinedGeometry *RG = FindInRGeom(Geom, Times, ETimes, type);
935  if (RG) { return RG; }
936 
937  switch (Geom)
938  {
939  case Geometry::POINT:
940  {
941  RG = new RefinedGeometry(1, 1, 0);
942  RG->Times = 1;
943  RG->ETimes = 0;
944  RG->Type = type;
945  RG->RefPts.IntPoint(0).x = cp[0];
946  RG->RefGeoms[0] = 0;
947 
948  RGeom[Geometry::POINT].Append(RG);
949  return RG;
950  }
951 
952  case Geometry::SEGMENT:
953  {
954  RG = new RefinedGeometry(Times+1, 2*Times, 0);
955  RG->Times = Times;
956  RG->ETimes = 0;
957  RG->Type = type;
958  for (i = 0; i <= Times; i++)
959  {
960  IntegrationPoint &ip = RG->RefPts.IntPoint(i);
961  ip.x = cp[i];
962  }
963  Array<int> &G = RG->RefGeoms;
964  for (i = 0; i < Times; i++)
965  {
966  G[2*i+0] = i;
967  G[2*i+1] = i+1;
968  }
969 
970  RGeom[Geometry::SEGMENT].Append(RG);
971  return RG;
972  }
973 
974  case Geometry::TRIANGLE:
975  {
976  RG = new RefinedGeometry((Times+1)*(Times+2)/2, 3*Times*Times,
977  3*Times*(ETimes+1), 3*Times);
978  RG->Times = Times;
979  RG->ETimes = ETimes;
980  RG->Type = type;
981  for (k = j = 0; j <= Times; j++)
982  for (i = 0; i <= Times-j; i++, k++)
983  {
984  IntegrationPoint &ip = RG->RefPts.IntPoint(k);
985  ip.x = cp[i]/(cp[i] + cp[j] + cp[Times-i-j]);
986  ip.y = cp[j]/(cp[i] + cp[j] + cp[Times-i-j]);
987  }
988  Array<int> &G = RG->RefGeoms;
989  for (l = k = j = 0; j < Times; j++, k++)
990  for (i = 0; i < Times-j; i++, k++)
991  {
992  G[l++] = k;
993  G[l++] = k+1;
994  G[l++] = k+Times-j+1;
995  if (i+j+1 < Times)
996  {
997  G[l++] = k+1;
998  G[l++] = k+Times-j+2;
999  G[l++] = k+Times-j+1;
1000  }
1001  }
1002  Array<int> &E = RG->RefEdges;
1003  int lb = 0, li = 2*RG->NumBdrEdges;
1004  // horizontal edges
1005  for (k = 0; k < Times; k += Times/ETimes)
1006  {
1007  int &lt = (k == 0) ? lb : li;
1008  j = k*(Times+1)-((k-1)*k)/2;
1009  for (i = 0; i < Times-k; i++)
1010  {
1011  E[lt++] = j; j++;
1012  E[lt++] = j;
1013  }
1014  }
1015  // diagonal edges
1016  for (k = Times; k > 0; k -= Times/ETimes)
1017  {
1018  int &lt = (k == Times) ? lb : li;
1019  j = k;
1020  for (i = 0; i < k; i++)
1021  {
1022  E[lt++] = j; j += Times-i;
1023  E[lt++] = j;
1024  }
1025  }
1026  // vertical edges
1027  for (k = 0; k < Times; k += Times/ETimes)
1028  {
1029  int &lt = (k == 0) ? lb : li;
1030  j = k;
1031  for (i = 0; i < Times-k; i++)
1032  {
1033  E[lt++] = j; j += Times-i+1;
1034  E[lt++] = j;
1035  }
1036  }
1037 
1038  RGeom[Geometry::TRIANGLE].Append(RG);
1039  return RG;
1040  }
1041 
1042  case Geometry::SQUARE:
1043  {
1044  RG = new RefinedGeometry((Times+1)*(Times+1), 4*Times*Times,
1045  4*(ETimes+1)*Times, 4*Times);
1046  RG->Times = Times;
1047  RG->ETimes = ETimes;
1048  RG->Type = type;
1049  for (k = j = 0; j <= Times; j++)
1050  for (i = 0; i <= Times; i++, k++)
1051  {
1052  IntegrationPoint &ip = RG->RefPts.IntPoint(k);
1053  ip.x = cp[i];
1054  ip.y = cp[j];
1055  }
1056  Array<int> &G = RG->RefGeoms;
1057  for (l = k = j = 0; j < Times; j++, k++)
1058  for (i = 0; i < Times; i++, k++)
1059  {
1060  G[l++] = k;
1061  G[l++] = k+1;
1062  G[l++] = k+Times+2;
1063  G[l++] = k+Times+1;
1064  }
1065  Array<int> &E = RG->RefEdges;
1066  int lb = 0, li = 2*RG->NumBdrEdges;
1067  // horizontal edges
1068  for (k = 0; k <= Times; k += Times/ETimes)
1069  {
1070  int &lt = (k == 0 || k == Times) ? lb : li;
1071  for (i = 0, j = k*(Times+1); i < Times; i++)
1072  {
1073  E[lt++] = j; j++;
1074  E[lt++] = j;
1075  }
1076  }
1077  // vertical edges (in right-to-left order)
1078  for (k = Times; k >= 0; k -= Times/ETimes)
1079  {
1080  int &lt = (k == Times || k == 0) ? lb : li;
1081  for (i = 0, j = k; i < Times; i++)
1082  {
1083  E[lt++] = j; j += Times+1;
1084  E[lt++] = j;
1085  }
1086  }
1087 
1088  RGeom[Geometry::SQUARE].Append(RG);
1089  return RG;
1090  }
1091 
1092  case Geometry::CUBE:
1093  {
1094  RG = new RefinedGeometry ((Times+1)*(Times+1)*(Times+1),
1095  8*Times*Times*Times, 0);
1096  RG->Times = Times;
1097  RG->ETimes = ETimes;
1098  RG->Type = type;
1099  for (l = k = 0; k <= Times; k++)
1100  for (j = 0; j <= Times; j++)
1101  for (i = 0; i <= Times; i++, l++)
1102  {
1103  IntegrationPoint &ip = RG->RefPts.IntPoint(l);
1104  ip.x = cp[i];
1105  ip.y = cp[j];
1106  ip.z = cp[k];
1107  }
1108  Array<int> &G = RG->RefGeoms;
1109  for (l = k = 0; k < Times; k++)
1110  for (j = 0; j < Times; j++)
1111  for (i = 0; i < Times; i++)
1112  {
1113  G[l++] = i+0 + (j+0 + (k+0) * (Times+1)) * (Times+1);
1114  G[l++] = i+1 + (j+0 + (k+0) * (Times+1)) * (Times+1);
1115  G[l++] = i+1 + (j+1 + (k+0) * (Times+1)) * (Times+1);
1116  G[l++] = i+0 + (j+1 + (k+0) * (Times+1)) * (Times+1);
1117  G[l++] = i+0 + (j+0 + (k+1) * (Times+1)) * (Times+1);
1118  G[l++] = i+1 + (j+0 + (k+1) * (Times+1)) * (Times+1);
1119  G[l++] = i+1 + (j+1 + (k+1) * (Times+1)) * (Times+1);
1120  G[l++] = i+0 + (j+1 + (k+1) * (Times+1)) * (Times+1);
1121  }
1122 
1123  RGeom[Geometry::CUBE].Append(RG);
1124  return RG;
1125  }
1126 
1127  case Geometry::TETRAHEDRON:
1128  {
1129  // subdivide the tetrahedron with vertices
1130  // (0,0,0), (0,0,1), (1,1,1), (0,1,1)
1131 
1132  // vertices: 0 <= i <= j <= k <= Times
1133  // (3-combination with repetitions)
1134  // number of vertices: (n+3)*(n+2)*(n+1)/6, n = Times
1135 
1136  // elements: the vertices are: v1=(i,j,k), v2=v1+u1, v3=v2+u2, v4=v3+u3
1137  // where 0 <= i <= j <= k <= n-1 and
1138  // u1,u2,u3 is a permutation of (1,0,0),(0,1,0),(0,0,1)
1139  // such that all v2,v3,v4 have non-decreasing components
1140  // number of elements: n^3
1141 
1142  const int n = Times;
1143  RG = new RefinedGeometry((n+3)*(n+2)*(n+1)/6, 4*n*n*n, 0);
1144  RG->Times = Times;
1145  RG->ETimes = ETimes;
1146  RG->Type = type;
1147  // enumerate and define the vertices
1148  Array<int> vi((n+1)*(n+1)*(n+1));
1149  vi = -1;
1150  m = 0;
1151  for (k = 0; k <= n; k++)
1152  for (j = 0; j <= k; j++)
1153  for (i = 0; i <= j; i++)
1154  {
1155  IntegrationPoint &ip = RG->RefPts.IntPoint(m);
1156  // map the coordinates to the reference tetrahedron
1157  // (0,0,0) -> (0,0,0)
1158  // (0,0,1) -> (1,0,0)
1159  // (1,1,1) -> (0,1,0)
1160  // (0,1,1) -> (0,0,1)
1161  double w = cp[k-j] + cp[i] + cp[j-i] + cp[Times-k];
1162  ip.x = cp[k-j]/w;
1163  ip.y = cp[i]/w;
1164  ip.z = cp[j-i]/w;
1165  l = i + (j + k * (n+1)) * (n+1);
1166  vi[l] = m;
1167  m++;
1168  }
1169  if (m != (n+3)*(n+2)*(n+1)/6)
1170  {
1171  mfem_error("GeometryRefiner::Refine() for TETRAHEDRON #1");
1172  }
1173  // elements
1174  Array<int> &G = RG->RefGeoms;
1175  m = 0;
1176  for (k = 0; k < n; k++)
1177  for (j = 0; j <= k; j++)
1178  for (i = 0; i <= j; i++)
1179  {
1180  // the ordering of the vertices is chosen to ensure:
1181  // 1) correct orientation
1182  // 2) the x,y,z edges are in the set of edges
1183  // {(0,1),(2,3), (0,2),(1,3)}
1184  // (goal is to ensure that subsequent refinement using
1185  // this procedure preserves the six tetrahedral shapes)
1186 
1187  // zyx: (i,j,k)-(i,j,k+1)-(i+1,j+1,k+1)-(i,j+1,k+1)
1188  G[m++] = vi[i+0 + (j+0 + (k+0) * (n+1)) * (n+1)];
1189  G[m++] = vi[i+0 + (j+0 + (k+1) * (n+1)) * (n+1)];
1190  G[m++] = vi[i+1 + (j+1 + (k+1) * (n+1)) * (n+1)];
1191  G[m++] = vi[i+0 + (j+1 + (k+1) * (n+1)) * (n+1)];
1192  if (j < k)
1193  {
1194  // yzx: (i,j,k)-(i+1,j+1,k+1)-(i,j+1,k)-(i,j+1,k+1)
1195  G[m++] = vi[i+0 + (j+0 + (k+0) * (n+1)) * (n+1)];
1196  G[m++] = vi[i+1 + (j+1 + (k+1) * (n+1)) * (n+1)];
1197  G[m++] = vi[i+0 + (j+1 + (k+0) * (n+1)) * (n+1)];
1198  G[m++] = vi[i+0 + (j+1 + (k+1) * (n+1)) * (n+1)];
1199  // yxz: (i,j,k)-(i,j+1,k)-(i+1,j+1,k+1)-(i+1,j+1,k)
1200  G[m++] = vi[i+0 + (j+0 + (k+0) * (n+1)) * (n+1)];
1201  G[m++] = vi[i+0 + (j+1 + (k+0) * (n+1)) * (n+1)];
1202  G[m++] = vi[i+1 + (j+1 + (k+1) * (n+1)) * (n+1)];
1203  G[m++] = vi[i+1 + (j+1 + (k+0) * (n+1)) * (n+1)];
1204  }
1205  if (i < j)
1206  {
1207  // xzy: (i,j,k)-(i+1,j,k)-(i+1,j+1,k+1)-(i+1,j,k+1)
1208  G[m++] = vi[i+0 + (j+0 + (k+0) * (n+1)) * (n+1)];
1209  G[m++] = vi[i+1 + (j+0 + (k+0) * (n+1)) * (n+1)];
1210  G[m++] = vi[i+1 + (j+1 + (k+1) * (n+1)) * (n+1)];
1211  G[m++] = vi[i+1 + (j+0 + (k+1) * (n+1)) * (n+1)];
1212  if (j < k)
1213  {
1214  // xyz: (i,j,k)-(i+1,j+1,k+1)-(i+1,j,k)-(i+1,j+1,k)
1215  G[m++] = vi[i+0 + (j+0 + (k+0) * (n+1)) * (n+1)];
1216  G[m++] = vi[i+1 + (j+1 + (k+1) * (n+1)) * (n+1)];
1217  G[m++] = vi[i+1 + (j+0 + (k+0) * (n+1)) * (n+1)];
1218  G[m++] = vi[i+1 + (j+1 + (k+0) * (n+1)) * (n+1)];
1219  }
1220  // zxy: (i,j,k)-(i+1,j+1,k+1)-(i,j,k+1)-(i+1,j,k+1)
1221  G[m++] = vi[i+0 + (j+0 + (k+0) * (n+1)) * (n+1)];
1222  G[m++] = vi[i+1 + (j+1 + (k+1) * (n+1)) * (n+1)];
1223  G[m++] = vi[i+0 + (j+0 + (k+1) * (n+1)) * (n+1)];
1224  G[m++] = vi[i+1 + (j+0 + (k+1) * (n+1)) * (n+1)];
1225  }
1226  }
1227  if (m != 4*n*n*n)
1228  {
1229  mfem_error("GeometryRefiner::Refine() for TETRAHEDRON #2");
1230  }
1231  for (i = 0; i < m; i++)
1232  if (G[i] < 0)
1233  {
1234  mfem_error("GeometryRefiner::Refine() for TETRAHEDRON #3");
1235  }
1236 
1237  RGeom[Geometry::TETRAHEDRON].Append(RG);
1238  return RG;
1239  }
1240 
1241  case Geometry::PRISM:
1242  {
1243  const int n = Times;
1244  RG = new RefinedGeometry ((n+1)*(n+1)*(n+2)/2, 6*n*n*n, 0);
1245  RG->Times = Times;
1246  RG->ETimes = ETimes;
1247  RG->Type = type;
1248  // enumerate and define the vertices
1249  m = 0;
1250  for (l = k = 0; k <= n; k++)
1251  for (j = 0; j <= n; j++)
1252  for (i = 0; i <= n-j; i++, l++)
1253  {
1254  IntegrationPoint &ip = RG->RefPts.IntPoint(l);
1255  if (type == 0)
1256  {
1257  ip.x = double(i) / n;
1258  ip.y = double(j) / n;
1259  ip.z = double(k) / n;
1260  }
1261  else
1262  {
1263  ip.x = cp[i]/(cp[i] + cp[j] + cp[n-i-j]);
1264  ip.y = cp[j]/(cp[i] + cp[j] + cp[n-i-j]);
1265  ip.z = cp[k];
1266  }
1267  m++;
1268  }
1269  if (m != (n+1)*(n+1)*(n+2)/2)
1270  {
1271  mfem_error("GeometryRefiner::Refine() for PRISM #1");
1272  }
1273  // elements
1274  Array<int> &G = RG->RefGeoms;
1275  m = 0;
1276  for (m = k = 0; k < n; k++)
1277  for (l = j = 0; j < n; j++, l++)
1278  for (i = 0; i < n-j; i++, l++)
1279  {
1280  G[m++] = l + (k+0) * (n+1) * (n+2) / 2;
1281  G[m++] = l + 1 + (k+0) * (n+1) * (n+2) / 2;
1282  G[m++] = l - j + (2 + (k+0) * (n+2)) * (n+1) / 2;
1283  G[m++] = l + (k+1) * (n+1) * (n+2) / 2;
1284  G[m++] = l + 1 + (k+1) * (n+1) * (Times+2) / 2;
1285  G[m++] = l - j + (2 + (k+1) * (n+2)) * (n+1) / 2;
1286  if (i+j+1 < n)
1287  {
1288  G[m++] = l + 1 + (k+0) * (n+1) * (n+2)/2;
1289  G[m++] = l - j + (2 + (k+0) * (n+1)) * (n+2) / 2;
1290  G[m++] = l - j + (2 + (k+0) * (n+2)) * (n+1) / 2;
1291  G[m++] = l + 1 + (k+1) * (n+1) * (n+2) / 2;
1292  G[m++] = l - j + (2 + (k+1) * (n+1)) * (n+2) / 2;
1293  G[m++] = l - j + (2 + (k+1) * (n+2)) * (n+1) / 2;
1294  }
1295  }
1296  if (m != 6*n*n*n)
1297  {
1298  mfem_error("GeometryRefiner::Refine() for PRISM #2");
1299  }
1300  for (i = 0; i < m; i++)
1301  if (G[i] < 0)
1302  {
1303  mfem_error("GeometryRefiner::Refine() for PRISM #3");
1304  }
1305 
1306  RGeom[Geometry::PRISM].Append(RG);
1307  return RG;
1308  }
1309 
1310  default:
1311 
1312  return NULL;
1313  }
1314 }
1315 
1316 const IntegrationRule *GeometryRefiner::RefineInterior(Geometry::Type Geom,
1317  int Times)
1318 {
1319  IntegrationRule *ir = NULL;
1320 
1321  switch (Geom)
1322  {
1323  case Geometry::SEGMENT:
1324  {
1325  if (Times < 2)
1326  {
1327  return NULL;
1328  }
1329  ir = FindInIntPts(Geom, Times-1);
1330  if (ir == NULL)
1331  {
1332  ir = new IntegrationRule(Times-1);
1333  for (int i = 1; i < Times; i++)
1334  {
1335  IntegrationPoint &ip = ir->IntPoint(i-1);
1336  ip.x = double(i) / Times;
1337  ip.y = ip.z = 0.0;
1338  }
1339  }
1340  }
1341  break;
1342 
1343  case Geometry::TRIANGLE:
1344  {
1345  if (Times < 3)
1346  {
1347  return NULL;
1348  }
1349  ir = FindInIntPts(Geom, ((Times-1)*(Times-2))/2);
1350  if (ir == NULL)
1351  {
1352  ir = new IntegrationRule(((Times-1)*(Times-2))/2);
1353  for (int k = 0, j = 1; j < Times-1; j++)
1354  for (int i = 1; i < Times-j; i++, k++)
1355  {
1356  IntegrationPoint &ip = ir->IntPoint(k);
1357  ip.x = double(i) / Times;
1358  ip.y = double(j) / Times;
1359  ip.z = 0.0;
1360  }
1361  }
1362  }
1363  break;
1364 
1365  case Geometry::SQUARE:
1366  {
1367  if (Times < 2)
1368  {
1369  return NULL;
1370  }
1371  ir = FindInIntPts(Geom, (Times-1)*(Times-1));
1372  if (ir == NULL)
1373  {
1374  ir = new IntegrationRule((Times-1)*(Times-1));
1375  for (int k = 0, j = 1; j < Times; j++)
1376  for (int i = 1; i < Times; i++, k++)
1377  {
1378  IntegrationPoint &ip = ir->IntPoint(k);
1379  ip.x = double(i) / Times;
1380  ip.y = double(j) / Times;
1381  ip.z = 0.0;
1382  }
1383  }
1384  }
1385  break;
1386 
1387  default:
1388  mfem_error("GeometryRefiner::RefineInterior(...)");
1389  }
1390 
1391  if (ir) { IntPts[Geom].Append(ir); }
1392  return ir;
1393 }
1394 
1395 GeometryRefiner GlobGeometryRefiner;
1396 
1397 }
mfem::Array::Size
int Size() const
Logical size of the array.
Definition: array.hpp:124
mfem::Geometry::Constants< Geometry::POINT >::InvOrient
static const int InvOrient[NumOrient]
Definition: geom.hpp:116
mfem::IntegrationRule
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
mfem::IsoparametricTransformation::GetPointMat
DenseMatrix & GetPointMat()
Read and write access to the underlying point matrix describing the transformation.
Definition: eltrans.hpp:321
mfem::BasisType::GetNodalBasis
static int GetNodalBasis(int qpt_type)
Return the nodal BasisType corresponding to the Quadrature1D type.
Definition: fe.hpp:76
mfem::Geometry::NumGeom
static const int NumGeom
Definition: geom.hpp:41
mfem::Geometry::JacToPerfJac
void JacToPerfJac(int GeomType, const DenseMatrix &J, DenseMatrix &PJ) const
Definition: geom.cpp:737
mfem::Geometry::GetRandomPoint
static void GetRandomPoint(int GeomType, IntegrationPoint &ip)
Get a random point in the reference element specified by GeomType.
Definition: geom.cpp:247
mfem::Mult
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:476
mfem::Geometry::Constants< Geometry::TETRAHEDRON >::FaceVert
static const int FaceVert[NumFaces][MaxFaceVert]
Definition: geom.hpp:186
mfem::ElementTransformation::SetIntPoint
void SetIntPoint(const IntegrationPoint *ip)
Definition: eltrans.hpp:56
mfem::Geometry::Constants< Geometry::PRISM >::FaceVert
static const int FaceVert[NumFaces][MaxFaceVert]
Definition: geom.hpp:222
mfem::Geometry::Constants< Geometry::PRISM >::Edges
static const int Edges[NumEdges][2]
Definition: geom.hpp:218
mfem::DenseMatrix
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
mfem::Geometry::Volume
static const double Volume[NumGeom]
Definition: geom.hpp:45
mfem::Geometry::NumEdges
static const int NumEdges[NumGeom]
Definition: geom.hpp:49
mfem::RefinedGeometry::RefEdges
Array< int > RefEdges
Definition: geom.hpp:240
mfem::Geometry::GetVertices
const IntegrationRule * GetVertices(int GeomType)
Return an IntegrationRule consisting of all vertices of the given Geometry::Type, GeomType...
Definition: geom.cpp:228
mfem::GeometryRefiner::RefineInterior
const IntegrationRule * RefineInterior(Geometry::Type Geom, int Times)
Definition: geom.cpp:1316
mfem::Poly_1D::GetPoints
const double * GetPoints(const int p, const int btype)
Get the coordinates of the points of the given BasisType, btype.
Definition: fe.cpp:7279
mfem::Geometry::NumFaces
static const int NumFaces[NumGeom]
Definition: geom.hpp:50
mfem::Geometry::Constants< Geometry::TRIANGLE >::InvOrient
static const int InvOrient[NumOrient]
Definition: geom.hpp:154
mfem::IntegrationRule::IntPoint
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:248
mfem::Geometry::Constants< Geometry::TRIANGLE >::Edges
static const int Edges[NumEdges][2]
Definition: geom.hpp:136
mfem::Geometry::Dimension
static const int Dimension[NumGeom]
Definition: geom.hpp:46
mfem::Geometry::Constants< Geometry::CUBE >::FaceTypes
static const int FaceTypes[NumFaces]
Definition: geom.hpp:202
mfem::mfem_error
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:153
mfem::Geometry::Constants< Geometry::PRISM >::FaceTypes
static const int FaceTypes[NumFaces]
Definition: geom.hpp:220
mfem::ElementTransformation::Jacobian
const DenseMatrix & Jacobian()
Return the Jacobian matrix of the transformation at the currently set IntegrationPoint, using the method SetIntPoint().
Definition: eltrans.hpp:71
mfem::Geometry::NumVerts
static const int NumVerts[NumGeom]
Definition: geom.hpp:48
mfem::Geometry::Constants< Geometry::SEGMENT >::Edges
static const int Edges[NumEdges][2]
Definition: geom.hpp:124
mfem::GlobGeometryRefiner
GeometryRefiner GlobGeometryRefiner
Definition: geom.cpp:1395
mfem::RefinedGeometry::RefPts
IntegrationRule RefPts
Definition: geom.hpp:239
mfem::Geometry::NumBdrArray
static const int NumBdrArray[NumGeom]
Definition: geom.hpp:43
mfem::GeometryRefiner::Refine
RefinedGeometry * Refine(Geometry::Type Geom, int Times, int ETimes=1)
Definition: geom.cpp:925
mfem::Geometry::Constants< Geometry::SQUARE >::InvOrient
static const int InvOrient[NumOrient]
Definition: geom.hpp:174
mfem::Geometry::Constants< Geometry::SEGMENT >::Orient
static const int Orient[NumOrient][NumVert]
Definition: geom.hpp:127
mfem::Geometry::Constants< Geometry::TETRAHEDRON >::Edges
static const int Edges[NumEdges][2]
Definition: geom.hpp:182
mfem::Geometry::Name
static const char * Name[NumGeom]
Definition: geom.hpp:44
mfem::WedgeFE
class H1_WedgeElement WedgeFE
mfem::CalcInverse
void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
Definition: densemat.cpp:2197
mfem::Geometry::Constants< Geometry::SEGMENT >::InvOrient
static const int InvOrient[NumOrient]
Definition: geom.hpp:128
mfem::Geometry::Constants< Geometry::SQUARE >::Edges
static const int Edges[NumEdges][2]
Definition: geom.hpp:162
mfem::Geometry::ProjectPoint
static bool ProjectPoint(int GeomType, const IntegrationPoint &beg, IntegrationPoint &end)
Project a point end, onto the given Geometry::Type, GeomType.
Definition: geom.cpp:511
mfem::Geometry::GetPerfPointMat
void GetPerfPointMat(int GeomType, DenseMatrix &pm)
Definition: geom.cpp:664
mfem::TetrahedronFE
class Linear3DFiniteElement TetrahedronFE
Definition: fe.cpp:12625
mfem::TriangleFE
Linear2DFiniteElement TriangleFE
Definition: fe.cpp:12621
mfem::Geometry::Constants< Geometry::CUBE >::FaceVert
static const int FaceVert[NumFaces][MaxFaceVert]
Definition: geom.hpp:204
mfem::Geometry::Constants< Geometry::TRIANGLE >::FaceVert
static const int FaceVert[NumFaces][NumVert]
Definition: geom.hpp:144
mfem::IntegrationPoint
Class for integration point with weight.
Definition: intrules.hpp:25
mfem::Geometry::DimStart
static const int DimStart[MaxDim+2]
Definition: geom.hpp:47
mfem::DenseMatrix::Diag
void Diag(double c, int n)
Creates n x n diagonal matrix with diagonal elements c.
Definition: densemat.cpp:1346
mfem::Geometry::Constants< Geometry::SQUARE >::FaceVert
static const int FaceVert[NumFaces][NumVert]
Definition: geom.hpp:170
mfem::Geometry::Constants< Geometry::SQUARE >::Orient
static const int Orient[NumOrient][NumVert]
Definition: geom.hpp:173
dim
int dim
Definition: ex24.cpp:43
mfem::IsoparametricTransformation::SetFE
void SetFE(const FiniteElement *FE)
Definition: eltrans.hpp:307
mfem::Geometry::Constants< Geometry::CUBE >::Edges
static const int Edges[NumEdges][2]
Definition: geom.hpp:200
mfem::Geometry::Constants< Geometry::POINT >::Orient
static const int Orient[NumOrient][NumVert]
Definition: geom.hpp:115
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:348
mfem::DenseMatrix::SetSize
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:90
mfem::out
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition: globals.hpp:66
mfem::RefinedGeometry::RefGeoms
Array< int > RefGeoms
Definition: geom.hpp:240
mfem::Geometry::Constants< Geometry::TRIANGLE >::Orient
static const int Orient[NumOrient][NumVert]
Definition: geom.hpp:152
mfem::Quadrature1D::ClosedUniform
aka closed Newton-Cotes
Definition: intrules.hpp:295
mfem::poly1d
Poly_1D poly1d
Definition: fe.cpp:7351
Generated on Fri Feb 18 2022 18:54:17 for MFEM by   doxygen 1.8.5