MFEM  v4.2.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.SetIntPoint(&GeomCenter[TRIANGLE]);
192  *GeomToPerfGeomJac[TRIANGLE] = tri_T.Jacobian();
193  CalcInverse(tri_T.Jacobian(), *PerfGeomToGeomJac[TRIANGLE]);
194  }
195  GeomToPerfGeomJac[SQUARE]->Diag(1.0, 2);
196  {
197  IsoparametricTransformation tet_T;
198  tet_T.SetFE(&TetrahedronFE);
199  GetPerfPointMat (TETRAHEDRON, tet_T.GetPointMat());
200  tet_T.SetIntPoint(&GeomCenter[TETRAHEDRON]);
201  *GeomToPerfGeomJac[TETRAHEDRON] = tet_T.Jacobian();
202  CalcInverse(tet_T.Jacobian(), *PerfGeomToGeomJac[TETRAHEDRON]);
203  }
204  GeomToPerfGeomJac[CUBE]->Diag(1.0, 3);
205  {
206  IsoparametricTransformation pri_T;
207  pri_T.SetFE(&WedgeFE);
208  GetPerfPointMat (PRISM, pri_T.GetPointMat());
209  pri_T.SetIntPoint(&GeomCenter[PRISM]);
210  *GeomToPerfGeomJac[PRISM] = pri_T.Jacobian();
211  CalcInverse(pri_T.Jacobian(), *PerfGeomToGeomJac[PRISM]);
212  }
213 }
214 
215 Geometry::~Geometry()
216 {
217  for (int i = 0; i < NumGeom; i++)
218  {
219  delete PerfGeomToGeomJac[i];
220  delete GeomToPerfGeomJac[i];
221  delete GeomVert[i];
222  }
223 }
224 
225 const IntegrationRule * Geometry::GetVertices(int GeomType)
226 {
227  switch (GeomType)
228  {
229  case Geometry::POINT: return GeomVert[0];
230  case Geometry::SEGMENT: return GeomVert[1];
231  case Geometry::TRIANGLE: return GeomVert[2];
232  case Geometry::SQUARE: return GeomVert[3];
233  case Geometry::TETRAHEDRON: return GeomVert[4];
234  case Geometry::CUBE: return GeomVert[5];
235  case Geometry::PRISM: return GeomVert[6];
236  default:
237  mfem_error ("Geometry::GetVertices(...)");
238  }
239  // make some compilers happy.
240  return GeomVert[0];
241 }
242 
243 // static method
244 void Geometry::GetRandomPoint(int GeomType, IntegrationPoint &ip)
245 {
246  switch (GeomType)
247  {
248  case Geometry::POINT:
249  ip.x = 0.0;
250  break;
251  case Geometry::SEGMENT:
252  ip.x = double(rand()) / RAND_MAX;
253  break;
254  case Geometry::TRIANGLE:
255  ip.x = double(rand()) / RAND_MAX;
256  ip.y = double(rand()) / RAND_MAX;
257  if (ip.x + ip.y > 1.0)
258  {
259  ip.x = 1.0 - ip.x;
260  ip.y = 1.0 - ip.y;
261  }
262  break;
263  case Geometry::SQUARE:
264  ip.x = double(rand()) / RAND_MAX;
265  ip.y = double(rand()) / RAND_MAX;
266  break;
267  case Geometry::TETRAHEDRON:
268  ip.x = double(rand()) / RAND_MAX;
269  ip.y = double(rand()) / RAND_MAX;
270  ip.z = double(rand()) / RAND_MAX;
271  // map to the triangular prism obtained by extruding the reference
272  // triangle in z direction
273  if (ip.x + ip.y > 1.0)
274  {
275  ip.x = 1.0 - ip.x;
276  ip.y = 1.0 - ip.y;
277  }
278  // split the prism into 3 parts: 1 is the reference tet, and the
279  // other two tets (as given below) are mapped to the reference tet
280  if (ip.x + ip.z > 1.0)
281  {
282  // tet with vertices: (0,0,1),(1,0,1),(0,1,1),(1,0,0)
283  ip.x = ip.x + ip.z - 1.0;
284  // ip.y = ip.y;
285  ip.z = 1.0 - ip.z;
286  // mapped to: (0,0,0),(1,0,0),(0,1,0),(0,0,1)
287  }
288  else if (ip.x + ip.y + ip.z > 1.0)
289  {
290  // tet with vertices: (0,1,1),(0,1,0),(0,0,1),(1,0,0)
291  double x = ip.x;
292  ip.x = 1.0 - x - ip.z;
293  ip.y = 1.0 - x - ip.y;
294  ip.z = x;
295  // mapped to: (0,0,0),(1,0,0),(0,1,0),(0,0,1)
296  }
297  break;
298  case Geometry::CUBE:
299  ip.x = double(rand()) / RAND_MAX;
300  ip.y = double(rand()) / RAND_MAX;
301  ip.z = double(rand()) / RAND_MAX;
302  break;
303  case Geometry::PRISM:
304  ip.x = double(rand()) / RAND_MAX;
305  ip.y = double(rand()) / RAND_MAX;
306  ip.z = double(rand()) / RAND_MAX;
307  if (ip.x + ip.y > 1.0)
308  {
309  ip.x = 1.0 - ip.x;
310  ip.y = 1.0 - ip.y;
311  }
312  break;
313  default:
314  MFEM_ABORT("Unknown type of reference element!");
315  }
316 }
317 
318 
319 namespace internal
320 {
321 
322 // Fuzzy equality operator with absolute tolerance eps.
323 inline bool NearlyEqual(double x, double y, double eps)
324 {
325  return std::abs(x-y) <= eps;
326 }
327 
328 // Fuzzy greater than comparison operator with absolute tolerance eps.
329 // Returns true when x is greater than y by at least eps.
330 inline bool FuzzyGT(double x, double y, double eps)
331 {
332  return (x > y + eps);
333 }
334 
335 // Fuzzy less than comparison operator with absolute tolerance eps.
336 // Returns true when x is less than y by at least eps.
337 inline bool FuzzyLT(double x, double y, double eps)
338 {
339  return (x < y - eps);
340 }
341 
342 }
343 
344 // static method
345 bool Geometry::CheckPoint(int GeomType, const IntegrationPoint &ip)
346 {
347  switch (GeomType)
348  {
349  case Geometry::POINT:
350  if (ip.x != 0.0) { return false; }
351  break;
352  case Geometry::SEGMENT:
353  if (ip.x < 0.0 || ip.x > 1.0) { return false; }
354  break;
355  case Geometry::TRIANGLE:
356  if (ip.x < 0.0 || ip.y < 0.0 || ip.x+ip.y > 1.0) { return false; }
357  break;
358  case Geometry::SQUARE:
359  if (ip.x < 0.0 || ip.x > 1.0 || ip.y < 0.0 || ip.y > 1.0)
360  { return false; }
361  break;
362  case Geometry::TETRAHEDRON:
363  if (ip.x < 0.0 || ip.y < 0.0 || ip.z < 0.0 ||
364  ip.x+ip.y+ip.z > 1.0) { return false; }
365  break;
366  case Geometry::CUBE:
367  if (ip.x < 0.0 || ip.x > 1.0 || ip.y < 0.0 || ip.y > 1.0 ||
368  ip.z < 0.0 || ip.z > 1.0) { return false; }
369  break;
370  case Geometry::PRISM:
371  if (ip.x < 0.0 || ip.y < 0.0 || ip.x+ip.y > 1.0 ||
372  ip.z < 0.0 || ip.z > 1.0) { return false; }
373  break;
374  default:
375  MFEM_ABORT("Unknown type of reference element!");
376  }
377  return true;
378 }
379 
380 bool Geometry::CheckPoint(int GeomType, const IntegrationPoint &ip, double eps)
381 {
382  switch (GeomType)
383  {
384  case Geometry::POINT:
385  if (! internal::NearlyEqual(ip.x, 0.0, eps))
386  {
387  return false;
388  }
389  break;
390  case Geometry::SEGMENT:
391  if ( internal::FuzzyLT(ip.x, 0.0, eps)
392  || internal::FuzzyGT(ip.x, 1.0, eps) )
393  {
394  return false;
395  }
396  break;
397  case Geometry::TRIANGLE:
398  if ( internal::FuzzyLT(ip.x, 0.0, eps)
399  || internal::FuzzyLT(ip.y, 0.0, eps)
400  || internal::FuzzyGT(ip.x+ip.y, 1.0, eps) )
401  {
402  return false;
403  }
404  break;
405  case Geometry::SQUARE:
406  if ( internal::FuzzyLT(ip.x, 0.0, eps)
407  || internal::FuzzyGT(ip.x, 1.0, eps)
408  || internal::FuzzyLT(ip.y, 0.0, eps)
409  || internal::FuzzyGT(ip.y, 1.0, eps) )
410  {
411  return false;
412  }
413  break;
414  case Geometry::TETRAHEDRON:
415  if ( internal::FuzzyLT(ip.x, 0.0, eps)
416  || internal::FuzzyLT(ip.y, 0.0, eps)
417  || internal::FuzzyLT(ip.z, 0.0, eps)
418  || internal::FuzzyGT(ip.x+ip.y+ip.z, 1.0, eps) )
419  {
420  return false;
421  }
422  break;
423  case Geometry::CUBE:
424  if ( internal::FuzzyLT(ip.x, 0.0, eps)
425  || internal::FuzzyGT(ip.x, 1.0, eps)
426  || internal::FuzzyLT(ip.y, 0.0, eps)
427  || internal::FuzzyGT(ip.y, 1.0, eps)
428  || internal::FuzzyLT(ip.z, 0.0, eps)
429  || internal::FuzzyGT(ip.z, 1.0, eps) )
430  {
431  return false;
432  }
433  break;
434  case Geometry::PRISM:
435  if ( internal::FuzzyLT(ip.x, 0.0, eps)
436  || internal::FuzzyLT(ip.y, 0.0, eps)
437  || internal::FuzzyGT(ip.x+ip.y, 1.0, eps)
438  || internal::FuzzyLT(ip.z, 0.0, eps)
439  || internal::FuzzyGT(ip.z, 1.0, eps) )
440  {
441  return false;
442  }
443  break;
444  default:
445  MFEM_ABORT("Unknown type of reference element!");
446  }
447  return true;
448 }
449 
450 
451 namespace internal
452 {
453 
454 template <int N, int dim>
455 inline bool IntersectSegment(double lbeg[N], double lend[N],
456  IntegrationPoint &end)
457 {
458  double t = 1.0;
459  bool out = false;
460  for (int i = 0; i < N; i++)
461  {
462  lbeg[i] = std::max(lbeg[i], 0.0); // remove round-off
463  if (lend[i] < 0.0)
464  {
465  out = true;
466  t = std::min(t, lbeg[i]/(lbeg[i]-lend[i]));
467  }
468  }
469  if (out)
470  {
471  if (dim >= 1) { end.x = t*lend[0] + (1.0-t)*lbeg[0]; }
472  if (dim >= 2) { end.y = t*lend[1] + (1.0-t)*lbeg[1]; }
473  if (dim >= 3) { end.z = t*lend[2] + (1.0-t)*lbeg[2]; }
474  return false;
475  }
476  return true;
477 }
478 
479 inline bool ProjectTriangle(double &x, double &y)
480 {
481  if (x < 0.0)
482  {
483  x = 0.0;
484  if (y < 0.0) { y = 0.0; }
485  else if (y > 1.0) { y = 1.0; }
486  return false;
487  }
488  if (y < 0.0)
489  {
490  if (x > 1.0) { x = 1.0; }
491  y = 0.0;
492  return false;
493  }
494  const double l3 = 1.0-x-y;
495  if (l3 < 0.0)
496  {
497  if (y - x > 1.0) { x = 0.0; y = 1.0; }
498  else if (y - x < -1.0) { x = 1.0; y = 0.0; }
499  else { x += l3/2; y += l3/2; }
500  return false;
501  }
502  return true;
503 }
504 
505 }
506 
507 // static method
508 bool Geometry::ProjectPoint(int GeomType, const IntegrationPoint &beg,
509  IntegrationPoint &end)
510 {
511  switch (GeomType)
512  {
513  case Geometry::POINT:
514  {
515  if (end.x != 0.0) { end.x = 0.0; return false; }
516  break;
517  }
518  case Geometry::SEGMENT:
519  {
520  if (end.x < 0.0) { end.x = 0.0; return false; }
521  if (end.x > 1.0) { end.x = 1.0; return false; }
522  break;
523  }
524  case Geometry::TRIANGLE:
525  {
526  double lend[3] = { end.x, end.y, 1.0-end.x-end.y };
527  double lbeg[3] = { beg.x, beg.y, 1.0-beg.x-beg.y };
528  return internal::IntersectSegment<3,2>(lbeg, lend, end);
529  }
530  case Geometry::SQUARE:
531  {
532  double lend[4] = { end.x, end.y, 1.0-end.x, 1.0-end.y };
533  double lbeg[4] = { beg.x, beg.y, 1.0-beg.x, 1.0-beg.y };
534  return internal::IntersectSegment<4,2>(lbeg, lend, end);
535  }
536  case Geometry::TETRAHEDRON:
537  {
538  double lend[4] = { end.x, end.y, end.z, 1.0-end.x-end.y-end.z };
539  double lbeg[4] = { beg.x, beg.y, beg.z, 1.0-beg.x-beg.y-beg.z };
540  return internal::IntersectSegment<4,3>(lbeg, lend, end);
541  }
542  case Geometry::CUBE:
543  {
544  double lend[6] = { end.x, end.y, end.z,
545  1.0-end.x, 1.0-end.y, 1.0-end.z
546  };
547  double lbeg[6] = { beg.x, beg.y, beg.z,
548  1.0-beg.x, 1.0-beg.y, 1.0-beg.z
549  };
550  return internal::IntersectSegment<6,3>(lbeg, lend, end);
551  }
552  case Geometry::PRISM:
553  {
554  double lend[5] = { end.x, end.y, end.z, 1.0-end.x-end.y, 1.0-end.z };
555  double lbeg[5] = { beg.x, beg.y, beg.z, 1.0-beg.x-beg.y, 1.0-beg.z };
556  return internal::IntersectSegment<5,3>(lbeg, lend, end);
557  }
558  default:
559  MFEM_ABORT("Unknown type of reference element!");
560  }
561  return true;
562 }
563 
564 // static method
565 bool Geometry::ProjectPoint(int GeomType, IntegrationPoint &ip)
566 {
567  // If ip is outside the element, replace it with the point on the boundary
568  // that is closest to the original ip and return false; otherwise, return
569  // true without changing ip.
570 
571  switch (GeomType)
572  {
573  case SEGMENT:
574  {
575  if (ip.x < 0.0) { ip.x = 0.0; return false; }
576  else if (ip.x > 1.0) { ip.x = 1.0; return false; }
577  return true;
578  }
579 
580  case TRIANGLE:
581  {
582  return internal::ProjectTriangle(ip.x, ip.y);
583  }
584 
585  case SQUARE:
586  {
587  bool in_x, in_y;
588  if (ip.x < 0.0) { in_x = false; ip.x = 0.0; }
589  else if (ip.x > 1.0) { in_x = false; ip.x = 1.0; }
590  else { in_x = true; }
591  if (ip.y < 0.0) { in_y = false; ip.y = 0.0; }
592  else if (ip.y > 1.0) { in_y = false; ip.y = 1.0; }
593  else { in_y = true; }
594  return in_x && in_y;
595  }
596 
597  case TETRAHEDRON:
598  {
599  if (ip.z < 0.0)
600  {
601  ip.z = 0.0;
602  internal::ProjectTriangle(ip.x, ip.y);
603  return false;
604  }
605  if (ip.y < 0.0)
606  {
607  ip.y = 0.0;
608  internal::ProjectTriangle(ip.x, ip.z);
609  return false;
610  }
611  if (ip.x < 0.0)
612  {
613  ip.x = 0.0;
614  internal::ProjectTriangle(ip.y, ip.z);
615  return false;
616  }
617  const double l4 = 1.0-ip.x-ip.y-ip.z;
618  if (l4 < 0.0)
619  {
620  const double l4_3 = l4/3;
621  ip.x += l4_3;
622  ip.y += l4_3;
623  internal::ProjectTriangle(ip.x, ip.y);
624  ip.z = 1.0-ip.x-ip.y;
625  return false;
626  }
627  return true;
628  }
629 
630  case CUBE:
631  {
632  bool in_x, in_y, in_z;
633  if (ip.x < 0.0) { in_x = false; ip.x = 0.0; }
634  else if (ip.x > 1.0) { in_x = false; ip.x = 1.0; }
635  else { in_x = true; }
636  if (ip.y < 0.0) { in_y = false; ip.y = 0.0; }
637  else if (ip.y > 1.0) { in_y = false; ip.y = 1.0; }
638  else { in_y = true; }
639  if (ip.z < 0.0) { in_z = false; ip.z = 0.0; }
640  else if (ip.z > 1.0) { in_z = false; ip.z = 1.0; }
641  else { in_z = true; }
642  return in_x && in_y && in_z;
643  }
644 
645  case PRISM:
646  {
647  bool in_tri, in_z;
648  in_tri = internal::ProjectTriangle(ip.x, ip.y);
649  if (ip.z < 0.0) { in_z = false; ip.z = 0.0; }
650  else if (ip.z > 1.0) { in_z = false; ip.z = 1.0; }
651  else { in_z = true; }
652  return in_tri && in_z;
653  }
654 
655  default:
656  MFEM_ABORT("Reference element type is not supported!");
657  }
658  return true;
659 }
660 
661 void Geometry::GetPerfPointMat(int GeomType, DenseMatrix &pm)
662 {
663  switch (GeomType)
664  {
665  case Geometry::SEGMENT:
666  {
667  pm.SetSize (1, 2);
668  pm(0,0) = 0.0;
669  pm(0,1) = 1.0;
670  }
671  break;
672 
673  case Geometry::TRIANGLE:
674  {
675  pm.SetSize (2, 3);
676  pm(0,0) = 0.0; pm(1,0) = 0.0;
677  pm(0,1) = 1.0; pm(1,1) = 0.0;
678  pm(0,2) = 0.5; pm(1,2) = 0.86602540378443864676;
679  }
680  break;
681 
682  case Geometry::SQUARE:
683  {
684  pm.SetSize (2, 4);
685  pm(0,0) = 0.0; pm(1,0) = 0.0;
686  pm(0,1) = 1.0; pm(1,1) = 0.0;
687  pm(0,2) = 1.0; pm(1,2) = 1.0;
688  pm(0,3) = 0.0; pm(1,3) = 1.0;
689  }
690  break;
691 
692  case Geometry::TETRAHEDRON:
693  {
694  pm.SetSize (3, 4);
695  pm(0,0) = 0.0; pm(1,0) = 0.0; pm(2,0) = 0.0;
696  pm(0,1) = 1.0; pm(1,1) = 0.0; pm(2,1) = 0.0;
697  pm(0,2) = 0.5; pm(1,2) = 0.86602540378443864676; pm(2,2) = 0.0;
698  pm(0,3) = 0.5; pm(1,3) = 0.28867513459481288225;
699  pm(2,3) = 0.81649658092772603273;
700  }
701  break;
702 
703  case Geometry::CUBE:
704  {
705  pm.SetSize (3, 8);
706  pm(0,0) = 0.0; pm(1,0) = 0.0; pm(2,0) = 0.0;
707  pm(0,1) = 1.0; pm(1,1) = 0.0; pm(2,1) = 0.0;
708  pm(0,2) = 1.0; pm(1,2) = 1.0; pm(2,2) = 0.0;
709  pm(0,3) = 0.0; pm(1,3) = 1.0; pm(2,3) = 0.0;
710  pm(0,4) = 0.0; pm(1,4) = 0.0; pm(2,4) = 1.0;
711  pm(0,5) = 1.0; pm(1,5) = 0.0; pm(2,5) = 1.0;
712  pm(0,6) = 1.0; pm(1,6) = 1.0; pm(2,6) = 1.0;
713  pm(0,7) = 0.0; pm(1,7) = 1.0; pm(2,7) = 1.0;
714  }
715  break;
716 
717  case Geometry::PRISM:
718  {
719  pm.SetSize (3, 6);
720  pm(0,0) = 0.0; pm(1,0) = 0.0; pm(2,0) = 0.0;
721  pm(0,1) = 1.0; pm(1,1) = 0.0; pm(2,1) = 0.0;
722  pm(0,2) = 0.5; pm(1,2) = 0.86602540378443864676; pm(2,2) = 0.0;
723  pm(0,3) = 0.0; pm(1,3) = 0.0; pm(2,3) = 1.0;
724  pm(0,4) = 1.0; pm(1,4) = 0.0; pm(2,4) = 1.0;
725  pm(0,5) = 0.5; pm(1,5) = 0.86602540378443864676; pm(2,5) = 1.0;
726  }
727  break;
728 
729  default:
730  mfem_error ("Geometry::GetPerfPointMat (...)");
731  }
732 }
733 
734 void Geometry::JacToPerfJac(int GeomType, const DenseMatrix &J,
735  DenseMatrix &PJ) const
736 {
737  if (PerfGeomToGeomJac[GeomType])
738  {
739  Mult(J, *PerfGeomToGeomJac[GeomType], PJ);
740  }
741  else
742  {
743  PJ = J;
744  }
745 }
746 
747 const int Geometry::NumBdrArray[NumGeom] = { 0, 2, 3, 4, 4, 6, 5 };
748 const int Geometry::Dimension[NumGeom] = { 0, 1, 2, 2, 3, 3, 3 };
749 const int Geometry::DimStart[MaxDim+2] =
750 { POINT, SEGMENT, TRIANGLE, TETRAHEDRON, NUM_GEOMETRIES };
751 const int Geometry::NumVerts[NumGeom] = { 1, 2, 3, 4, 4, 8, 6 };
752 const int Geometry::NumEdges[NumGeom] = { 0, 1, 3, 4, 6, 12, 9 };
753 const int Geometry::NumFaces[NumGeom] = { 0, 0, 1, 1, 4, 6, 5 };
754 
755 const int Geometry::
756 Constants<Geometry::POINT>::Orient[1][1] = {{0}};
757 const int Geometry::
758 Constants<Geometry::POINT>::InvOrient[1] = {0};
759 
760 const int Geometry::
761 Constants<Geometry::SEGMENT>::Edges[1][2] = { {0, 1} };
762 const int Geometry::
763 Constants<Geometry::SEGMENT>::Orient[2][2] = { {0, 1}, {1, 0} };
764 const int Geometry::
765 Constants<Geometry::SEGMENT>::InvOrient[2] = { 0, 1 };
766 
767 const int Geometry::
768 Constants<Geometry::TRIANGLE>::Edges[3][2] = {{0, 1}, {1, 2}, {2, 0}};
769 const int Geometry::
770 Constants<Geometry::TRIANGLE>::VertToVert::I[3] = {0, 2, 3};
771 const int Geometry::
772 Constants<Geometry::TRIANGLE>::VertToVert::J[3][2] = {{1, 0}, {2, -3}, {2, 1}};
773 const int Geometry::
774 Constants<Geometry::TRIANGLE>::FaceVert[1][3] = {{0, 1, 2}};
775 const int Geometry::
776 Constants<Geometry::TRIANGLE>::Orient[6][3] =
777 {
778  {0, 1, 2}, {1, 0, 2}, {2, 0, 1},
779  {2, 1, 0}, {1, 2, 0}, {0, 2, 1}
780 };
781 const int Geometry::
782 Constants<Geometry::TRIANGLE>::InvOrient[6] = { 0, 1, 4, 3, 2, 5 };
783 
784 const int Geometry::
785 Constants<Geometry::SQUARE>::Edges[4][2] = {{0, 1}, {1, 2}, {2, 3}, {3, 0}};
786 const int Geometry::
787 Constants<Geometry::SQUARE>::VertToVert::I[4] = {0, 2, 3, 4};
788 const int Geometry::
789 Constants<Geometry::SQUARE>::VertToVert::J[4][2] =
790 {{1, 0}, {3, -4}, {2, 1}, {3, 2}};
791 const int Geometry::
792 Constants<Geometry::SQUARE>::FaceVert[1][4] = {{0, 1, 2, 3}};
793 const int Geometry::
794 Constants<Geometry::SQUARE>::Orient[8][4] =
795 {
796  {0, 1, 2, 3}, {0, 3, 2, 1}, {1, 2, 3, 0}, {1, 0, 3, 2},
797  {2, 3, 0, 1}, {2, 1, 0, 3}, {3, 0, 1, 2}, {3, 2, 1, 0}
798 };
799 const int Geometry::
800 Constants<Geometry::SQUARE>::InvOrient[8] = { 0, 1, 6, 3, 4, 5, 2, 7 };
801 
802 const int Geometry::
803 Constants<Geometry::TETRAHEDRON>::Edges[6][2] =
804 {{0, 1}, {0, 2}, {0, 3}, {1, 2}, {1, 3}, {2, 3}};
805 const int Geometry::
806 Constants<Geometry::TETRAHEDRON>::FaceTypes[4] =
807 {
808  Geometry::TRIANGLE, Geometry::TRIANGLE,
809  Geometry::TRIANGLE, Geometry::TRIANGLE
810 };
811 const int Geometry::
812 Constants<Geometry::TETRAHEDRON>::FaceVert[4][3] =
813 {{1, 2, 3}, {0, 3, 2}, {0, 1, 3}, {0, 2, 1}};
814 const int Geometry::
815 Constants<Geometry::TETRAHEDRON>::VertToVert::I[4] = {0, 3, 5, 6};
816 const int Geometry::
817 Constants<Geometry::TETRAHEDRON>::VertToVert::J[6][2] =
818 {
819  {1, 0}, {2, 1}, {3, 2}, // 0,1:0 0,2:1 0,3:2
820  {2, 3}, {3, 4}, // 1,2:3 1,3:4
821  {3, 5} // 2,3:5
822 };
823 const int Geometry::
824 Constants<Geometry::TETRAHEDRON>::Orient[24][4] =
825 {
826  {0, 1, 2, 3}, {0, 1, 3, 2}, {0, 2, 3, 1}, {0, 2, 1, 3},
827  {0, 3, 1, 2}, {0, 3, 2, 1},
828  {1, 2, 0, 3}, {1, 2, 3, 0}, {1, 3, 2, 0}, {1, 3, 0, 2},
829  {1, 0, 3, 2}, {1, 0, 2, 3},
830  {2, 3, 0, 1}, {2, 3, 1, 0}, {2, 0, 1, 3}, {2, 0, 3, 1},
831  {2, 1, 3, 0}, {2, 1, 0, 3},
832  {3, 0, 2, 1}, {3, 0, 1, 2}, {3, 1, 0, 2}, {3, 1, 2, 0},
833  {3, 2, 1, 0}, {3, 2, 0, 1}
834 };
835 const int Geometry::
836 Constants<Geometry::TETRAHEDRON>::InvOrient[24] =
837 {
838  0, 1, 4, 3, 2, 5,
839  14, 19, 18, 15, 10, 11,
840  12, 23, 6, 9, 20, 17,
841  8, 7, 16, 21, 22, 13
842 };
843 
844 const int Geometry::
845 Constants<Geometry::CUBE>::Edges[12][2] =
846 {
847  {0, 1}, {1, 2}, {3, 2}, {0, 3}, {4, 5}, {5, 6},
848  {7, 6}, {4, 7}, {0, 4}, {1, 5}, {2, 6}, {3, 7}
849 };
850 const int Geometry::
851 Constants<Geometry::CUBE>::FaceTypes[6] =
852 {
853  Geometry::SQUARE, Geometry::SQUARE, Geometry::SQUARE,
854  Geometry::SQUARE, Geometry::SQUARE, Geometry::SQUARE
855 };
856 const int Geometry::
857 Constants<Geometry::CUBE>::FaceVert[6][4] =
858 {
859  {3, 2, 1, 0}, {0, 1, 5, 4}, {1, 2, 6, 5},
860  {2, 3, 7, 6}, {3, 0, 4, 7}, {4, 5, 6, 7}
861 };
862 const int Geometry::
863 Constants<Geometry::CUBE>::VertToVert::I[8] = {0, 3, 5, 7, 8, 10, 11, 12};
864 const int Geometry::
865 Constants<Geometry::CUBE>::VertToVert::J[12][2] =
866 {
867  {1, 0}, {3, 3}, {4, 8}, // 0,1:0 0,3:3 0,4:8
868  {2, 1}, {5, 9}, // 1,2:1 1,5:9
869  {3,-3}, {6,10}, // 2,3:-3 2,6:10
870  {7,11}, // 3,7:11
871  {5, 4}, {7, 7}, // 4,5:4 4,7:7
872  {6, 5}, // 5,6:5
873  {7,-7} // 6,7:-7
874 };
875 
876 const int Geometry::
877 Constants<Geometry::PRISM>::Edges[9][2] =
878 {{0, 1}, {1, 2}, {2, 0}, {3, 4}, {4, 5}, {5, 3}, {0, 3}, {1, 4}, {2, 5}};
879 const int Geometry::
880 Constants<Geometry::PRISM>::FaceTypes[5] =
881 {
882  Geometry::TRIANGLE, Geometry::TRIANGLE,
883  Geometry::SQUARE, Geometry::SQUARE, Geometry::SQUARE
884 };
885 const int Geometry::
886 Constants<Geometry::PRISM>::FaceVert[5][4] =
887 {{0, 2, 1, -1}, {3, 4, 5, -1}, {0, 1, 4, 3}, {1, 2, 5, 4}, {2, 0, 3, 5}};
888 const int Geometry::
889 Constants<Geometry::PRISM>::VertToVert::I[6] = {0, 3, 5, 6, 8, 9};
890 const int Geometry::
891 Constants<Geometry::PRISM>::VertToVert::J[9][2] =
892 {
893  {1, 0}, {2, -3}, {3, 6}, // 0,1:0 0,2:-3 0,3:6
894  {2, 1}, {4, 7}, // 1,2:1 1,4:7
895  {5, 8}, // 2,5:8
896  {4, 3}, {5, -6}, // 3,4:3 3,5:-6
897  {5, 4} // 4,5:4
898 };
899 
900 
901 GeometryRefiner::GeometryRefiner()
902 {
903  type = Quadrature1D::ClosedUniform;
904 }
905 
906 GeometryRefiner::~GeometryRefiner()
907 {
908  for (int i = 0; i < Geometry::NumGeom; i++)
909  {
910  for (int j = 0; j < RGeom[i].Size(); j++) { delete RGeom[i][j]; }
911  for (int j = 0; j < IntPts[i].Size(); j++) { delete IntPts[i][j]; }
912  }
913 }
914 
915 RefinedGeometry *GeometryRefiner::FindInRGeom(Geometry::Type Geom,
916  int Times, int ETimes,
917  int Type)
918 {
919  Array<RefinedGeometry *> &RGA = RGeom[Geom];
920  for (int i = 0; i < RGA.Size(); i++)
921  {
922  RefinedGeometry &RG = *RGA[i];
923  if (RG.Times == Times && RG.ETimes == ETimes && RG.Type == Type)
924  {
925  return &RG;
926  }
927  }
928  return NULL;
929 }
930 
931 IntegrationRule *GeometryRefiner::FindInIntPts(Geometry::Type Geom, int NPts)
932 {
933  Array<IntegrationRule *> &IPA = IntPts[Geom];
934  for (int i = 0; i < IPA.Size(); i++)
935  {
936  IntegrationRule &ir = *IPA[i];
937  if (ir.GetNPoints() == NPts) { return &ir; }
938  }
939  return NULL;
940 }
941 
942 RefinedGeometry * GeometryRefiner::Refine(Geometry::Type Geom,
943  int Times, int ETimes)
944 {
945  int i, j, k, l, m;
946 
947  Times = std::max(Times, 1);
948  ETimes = std::max(ETimes, 1);
949  const double *cp = poly1d.GetPoints(Times, BasisType::GetNodalBasis(type));
950 
951  RefinedGeometry *RG = FindInRGeom(Geom, Times, ETimes, type);
952  if (RG) { return RG; }
953 
954  switch (Geom)
955  {
956  case Geometry::POINT:
957  {
958  RG = new RefinedGeometry(1, 1, 0);
959  RG->Times = 1;
960  RG->ETimes = 0;
961  RG->Type = type;
962  RG->RefPts.IntPoint(0).x = cp[0];
963  RG->RefGeoms[0] = 0;
964 
965  RGeom[Geometry::POINT].Append(RG);
966  return RG;
967  }
968 
969  case Geometry::SEGMENT:
970  {
971  RG = new RefinedGeometry(Times+1, 2*Times, 0);
972  RG->Times = Times;
973  RG->ETimes = 0;
974  RG->Type = type;
975  for (i = 0; i <= Times; i++)
976  {
977  IntegrationPoint &ip = RG->RefPts.IntPoint(i);
978  ip.x = cp[i];
979  }
980  Array<int> &G = RG->RefGeoms;
981  for (i = 0; i < Times; i++)
982  {
983  G[2*i+0] = i;
984  G[2*i+1] = i+1;
985  }
986 
987  RGeom[Geometry::SEGMENT].Append(RG);
988  return RG;
989  }
990 
991  case Geometry::TRIANGLE:
992  {
993  RG = new RefinedGeometry((Times+1)*(Times+2)/2, 3*Times*Times,
994  3*Times*(ETimes+1), 3*Times);
995  RG->Times = Times;
996  RG->ETimes = ETimes;
997  RG->Type = type;
998  for (k = j = 0; j <= Times; j++)
999  for (i = 0; i <= Times-j; i++, k++)
1000  {
1001  IntegrationPoint &ip = RG->RefPts.IntPoint(k);
1002  ip.x = cp[i]/(cp[i] + cp[j] + cp[Times-i-j]);
1003  ip.y = cp[j]/(cp[i] + cp[j] + cp[Times-i-j]);
1004  }
1005  Array<int> &G = RG->RefGeoms;
1006  for (l = k = j = 0; j < Times; j++, k++)
1007  for (i = 0; i < Times-j; i++, k++)
1008  {
1009  G[l++] = k;
1010  G[l++] = k+1;
1011  G[l++] = k+Times-j+1;
1012  if (i+j+1 < Times)
1013  {
1014  G[l++] = k+1;
1015  G[l++] = k+Times-j+2;
1016  G[l++] = k+Times-j+1;
1017  }
1018  }
1019  Array<int> &E = RG->RefEdges;
1020  int lb = 0, li = 2*RG->NumBdrEdges;
1021  // horizontal edges
1022  for (k = 0; k < Times; k += Times/ETimes)
1023  {
1024  int &lt = (k == 0) ? lb : li;
1025  j = k*(Times+1)-((k-1)*k)/2;
1026  for (i = 0; i < Times-k; i++)
1027  {
1028  E[lt++] = j; j++;
1029  E[lt++] = j;
1030  }
1031  }
1032  // diagonal edges
1033  for (k = Times; k > 0; k -= Times/ETimes)
1034  {
1035  int &lt = (k == Times) ? lb : li;
1036  j = k;
1037  for (i = 0; i < k; i++)
1038  {
1039  E[lt++] = j; j += Times-i;
1040  E[lt++] = j;
1041  }
1042  }
1043  // vertical edges
1044  for (k = 0; k < Times; k += Times/ETimes)
1045  {
1046  int &lt = (k == 0) ? lb : li;
1047  j = k;
1048  for (i = 0; i < Times-k; i++)
1049  {
1050  E[lt++] = j; j += Times-i+1;
1051  E[lt++] = j;
1052  }
1053  }
1054 
1055  RGeom[Geometry::TRIANGLE].Append(RG);
1056  return RG;
1057  }
1058 
1059  case Geometry::SQUARE:
1060  {
1061  RG = new RefinedGeometry((Times+1)*(Times+1), 4*Times*Times,
1062  4*(ETimes+1)*Times, 4*Times);
1063  RG->Times = Times;
1064  RG->ETimes = ETimes;
1065  RG->Type = type;
1066  for (k = j = 0; j <= Times; j++)
1067  for (i = 0; i <= Times; i++, k++)
1068  {
1069  IntegrationPoint &ip = RG->RefPts.IntPoint(k);
1070  ip.x = cp[i];
1071  ip.y = cp[j];
1072  }
1073  Array<int> &G = RG->RefGeoms;
1074  for (l = k = j = 0; j < Times; j++, k++)
1075  for (i = 0; i < Times; i++, k++)
1076  {
1077  G[l++] = k;
1078  G[l++] = k+1;
1079  G[l++] = k+Times+2;
1080  G[l++] = k+Times+1;
1081  }
1082  Array<int> &E = RG->RefEdges;
1083  int lb = 0, li = 2*RG->NumBdrEdges;
1084  // horizontal edges
1085  for (k = 0; k <= Times; k += Times/ETimes)
1086  {
1087  int &lt = (k == 0 || k == Times) ? lb : li;
1088  for (i = 0, j = k*(Times+1); i < Times; i++)
1089  {
1090  E[lt++] = j; j++;
1091  E[lt++] = j;
1092  }
1093  }
1094  // vertical edges (in right-to-left order)
1095  for (k = Times; k >= 0; k -= Times/ETimes)
1096  {
1097  int &lt = (k == Times || k == 0) ? lb : li;
1098  for (i = 0, j = k; i < Times; i++)
1099  {
1100  E[lt++] = j; j += Times+1;
1101  E[lt++] = j;
1102  }
1103  }
1104 
1105  RGeom[Geometry::SQUARE].Append(RG);
1106  return RG;
1107  }
1108 
1109  case Geometry::CUBE:
1110  {
1111  RG = new RefinedGeometry ((Times+1)*(Times+1)*(Times+1),
1112  8*Times*Times*Times, 0);
1113  RG->Times = Times;
1114  RG->ETimes = ETimes;
1115  RG->Type = type;
1116  for (l = k = 0; k <= Times; k++)
1117  for (j = 0; j <= Times; j++)
1118  for (i = 0; i <= Times; i++, l++)
1119  {
1120  IntegrationPoint &ip = RG->RefPts.IntPoint(l);
1121  ip.x = cp[i];
1122  ip.y = cp[j];
1123  ip.z = cp[k];
1124  }
1125  Array<int> &G = RG->RefGeoms;
1126  for (l = k = 0; k < Times; k++)
1127  for (j = 0; j < Times; j++)
1128  for (i = 0; i < Times; i++)
1129  {
1130  G[l++] = i+0 + (j+0 + (k+0) * (Times+1)) * (Times+1);
1131  G[l++] = i+1 + (j+0 + (k+0) * (Times+1)) * (Times+1);
1132  G[l++] = i+1 + (j+1 + (k+0) * (Times+1)) * (Times+1);
1133  G[l++] = i+0 + (j+1 + (k+0) * (Times+1)) * (Times+1);
1134  G[l++] = i+0 + (j+0 + (k+1) * (Times+1)) * (Times+1);
1135  G[l++] = i+1 + (j+0 + (k+1) * (Times+1)) * (Times+1);
1136  G[l++] = i+1 + (j+1 + (k+1) * (Times+1)) * (Times+1);
1137  G[l++] = i+0 + (j+1 + (k+1) * (Times+1)) * (Times+1);
1138  }
1139 
1140  RGeom[Geometry::CUBE].Append(RG);
1141  return RG;
1142  }
1143 
1144  case Geometry::TETRAHEDRON:
1145  {
1146  // subdivide the tetrahedron with vertices
1147  // (0,0,0), (0,0,1), (1,1,1), (0,1,1)
1148 
1149  // vertices: 0 <= i <= j <= k <= Times
1150  // (3-combination with repetitions)
1151  // number of vertices: (n+3)*(n+2)*(n+1)/6, n = Times
1152 
1153  // elements: the vertices are: v1=(i,j,k), v2=v1+u1, v3=v2+u2, v4=v3+u3
1154  // where 0 <= i <= j <= k <= n-1 and
1155  // u1,u2,u3 is a permutation of (1,0,0),(0,1,0),(0,0,1)
1156  // such that all v2,v3,v4 have non-decreasing components
1157  // number of elements: n^3
1158 
1159  const int n = Times;
1160  RG = new RefinedGeometry((n+3)*(n+2)*(n+1)/6, 4*n*n*n, 0);
1161  RG->Times = Times;
1162  RG->ETimes = ETimes;
1163  RG->Type = type;
1164  // enumerate and define the vertices
1165  Array<int> vi((n+1)*(n+1)*(n+1));
1166  vi = -1;
1167  m = 0;
1168  for (k = 0; k <= n; k++)
1169  for (j = 0; j <= k; j++)
1170  for (i = 0; i <= j; i++)
1171  {
1172  IntegrationPoint &ip = RG->RefPts.IntPoint(m);
1173  // map the coordinates to the reference tetrahedron
1174  // (0,0,0) -> (0,0,0)
1175  // (0,0,1) -> (1,0,0)
1176  // (1,1,1) -> (0,1,0)
1177  // (0,1,1) -> (0,0,1)
1178  double w = cp[k-j] + cp[i] + cp[j-i] + cp[Times-k];
1179  ip.x = cp[k-j]/w;
1180  ip.y = cp[i]/w;
1181  ip.z = cp[j-i]/w;
1182  l = i + (j + k * (n+1)) * (n+1);
1183  vi[l] = m;
1184  m++;
1185  }
1186  if (m != (n+3)*(n+2)*(n+1)/6)
1187  {
1188  mfem_error("GeometryRefiner::Refine() for TETRAHEDRON #1");
1189  }
1190  // elements
1191  Array<int> &G = RG->RefGeoms;
1192  m = 0;
1193  for (k = 0; k < n; k++)
1194  for (j = 0; j <= k; j++)
1195  for (i = 0; i <= j; i++)
1196  {
1197  // the ordering of the vertices is chosen to ensure:
1198  // 1) correct orientation
1199  // 2) the x,y,z edges are in the set of edges
1200  // {(0,1),(2,3), (0,2),(1,3)}
1201  // (goal is to ensure that subsequent refinement using
1202  // this procedure preserves the six tetrahedral shapes)
1203 
1204  // zyx: (i,j,k)-(i,j,k+1)-(i+1,j+1,k+1)-(i,j+1,k+1)
1205  G[m++] = vi[i+0 + (j+0 + (k+0) * (n+1)) * (n+1)];
1206  G[m++] = vi[i+0 + (j+0 + (k+1) * (n+1)) * (n+1)];
1207  G[m++] = vi[i+1 + (j+1 + (k+1) * (n+1)) * (n+1)];
1208  G[m++] = vi[i+0 + (j+1 + (k+1) * (n+1)) * (n+1)];
1209  if (j < k)
1210  {
1211  // yzx: (i,j,k)-(i+1,j+1,k+1)-(i,j+1,k)-(i,j+1,k+1)
1212  G[m++] = vi[i+0 + (j+0 + (k+0) * (n+1)) * (n+1)];
1213  G[m++] = vi[i+1 + (j+1 + (k+1) * (n+1)) * (n+1)];
1214  G[m++] = vi[i+0 + (j+1 + (k+0) * (n+1)) * (n+1)];
1215  G[m++] = vi[i+0 + (j+1 + (k+1) * (n+1)) * (n+1)];
1216  // yxz: (i,j,k)-(i,j+1,k)-(i+1,j+1,k+1)-(i+1,j+1,k)
1217  G[m++] = vi[i+0 + (j+0 + (k+0) * (n+1)) * (n+1)];
1218  G[m++] = vi[i+0 + (j+1 + (k+0) * (n+1)) * (n+1)];
1219  G[m++] = vi[i+1 + (j+1 + (k+1) * (n+1)) * (n+1)];
1220  G[m++] = vi[i+1 + (j+1 + (k+0) * (n+1)) * (n+1)];
1221  }
1222  if (i < j)
1223  {
1224  // xzy: (i,j,k)-(i+1,j,k)-(i+1,j+1,k+1)-(i+1,j,k+1)
1225  G[m++] = vi[i+0 + (j+0 + (k+0) * (n+1)) * (n+1)];
1226  G[m++] = vi[i+1 + (j+0 + (k+0) * (n+1)) * (n+1)];
1227  G[m++] = vi[i+1 + (j+1 + (k+1) * (n+1)) * (n+1)];
1228  G[m++] = vi[i+1 + (j+0 + (k+1) * (n+1)) * (n+1)];
1229  if (j < k)
1230  {
1231  // xyz: (i,j,k)-(i+1,j+1,k+1)-(i+1,j,k)-(i+1,j+1,k)
1232  G[m++] = vi[i+0 + (j+0 + (k+0) * (n+1)) * (n+1)];
1233  G[m++] = vi[i+1 + (j+1 + (k+1) * (n+1)) * (n+1)];
1234  G[m++] = vi[i+1 + (j+0 + (k+0) * (n+1)) * (n+1)];
1235  G[m++] = vi[i+1 + (j+1 + (k+0) * (n+1)) * (n+1)];
1236  }
1237  // zxy: (i,j,k)-(i+1,j+1,k+1)-(i,j,k+1)-(i+1,j,k+1)
1238  G[m++] = vi[i+0 + (j+0 + (k+0) * (n+1)) * (n+1)];
1239  G[m++] = vi[i+1 + (j+1 + (k+1) * (n+1)) * (n+1)];
1240  G[m++] = vi[i+0 + (j+0 + (k+1) * (n+1)) * (n+1)];
1241  G[m++] = vi[i+1 + (j+0 + (k+1) * (n+1)) * (n+1)];
1242  }
1243  }
1244  if (m != 4*n*n*n)
1245  {
1246  mfem_error("GeometryRefiner::Refine() for TETRAHEDRON #2");
1247  }
1248  for (i = 0; i < m; i++)
1249  if (G[i] < 0)
1250  {
1251  mfem_error("GeometryRefiner::Refine() for TETRAHEDRON #3");
1252  }
1253 
1254  RGeom[Geometry::TETRAHEDRON].Append(RG);
1255  return RG;
1256  }
1257 
1258  case Geometry::PRISM:
1259  {
1260  const int n = Times;
1261  RG = new RefinedGeometry ((n+1)*(n+1)*(n+2)/2, 6*n*n*n, 0);
1262  RG->Times = Times;
1263  RG->ETimes = ETimes;
1264  RG->Type = type;
1265  // enumerate and define the vertices
1266  m = 0;
1267  for (l = k = 0; k <= n; k++)
1268  for (j = 0; j <= n; j++)
1269  for (i = 0; i <= n-j; i++, l++)
1270  {
1271  IntegrationPoint &ip = RG->RefPts.IntPoint(l);
1272  if (type == 0)
1273  {
1274  ip.x = double(i) / n;
1275  ip.y = double(j) / n;
1276  ip.z = double(k) / n;
1277  }
1278  else
1279  {
1280  ip.x = cp[i]/(cp[i] + cp[j] + cp[n-i-j]);
1281  ip.y = cp[j]/(cp[i] + cp[j] + cp[n-i-j]);
1282  ip.z = cp[k];
1283  }
1284  m++;
1285  }
1286  if (m != (n+1)*(n+1)*(n+2)/2)
1287  {
1288  mfem_error("GeometryRefiner::Refine() for PRISM #1");
1289  }
1290  // elements
1291  Array<int> &G = RG->RefGeoms;
1292  m = 0;
1293  for (m = k = 0; k < n; k++)
1294  for (l = j = 0; j < n; j++, l++)
1295  for (i = 0; i < n-j; i++, l++)
1296  {
1297  G[m++] = l + (k+0) * (n+1) * (n+2) / 2;
1298  G[m++] = l + 1 + (k+0) * (n+1) * (n+2) / 2;
1299  G[m++] = l - j + (2 + (k+0) * (n+2)) * (n+1) / 2;
1300  G[m++] = l + (k+1) * (n+1) * (n+2) / 2;
1301  G[m++] = l + 1 + (k+1) * (n+1) * (Times+2) / 2;
1302  G[m++] = l - j + (2 + (k+1) * (n+2)) * (n+1) / 2;
1303  if (i+j+1 < n)
1304  {
1305  G[m++] = l + 1 + (k+0) * (n+1) * (n+2)/2;
1306  G[m++] = l - j + (2 + (k+0) * (n+1)) * (n+2) / 2;
1307  G[m++] = l - j + (2 + (k+0) * (n+2)) * (n+1) / 2;
1308  G[m++] = l + 1 + (k+1) * (n+1) * (n+2) / 2;
1309  G[m++] = l - j + (2 + (k+1) * (n+1)) * (n+2) / 2;
1310  G[m++] = l - j + (2 + (k+1) * (n+2)) * (n+1) / 2;
1311  }
1312  }
1313  if (m != 6*n*n*n)
1314  {
1315  mfem_error("GeometryRefiner::Refine() for PRISM #2");
1316  }
1317  for (i = 0; i < m; i++)
1318  if (G[i] < 0)
1319  {
1320  mfem_error("GeometryRefiner::Refine() for PRISM #3");
1321  }
1322 
1323  RGeom[Geometry::PRISM].Append(RG);
1324  return RG;
1325  }
1326 
1327  default:
1328 
1329  return NULL;
1330  }
1331 }
1332 
1333 const IntegrationRule *GeometryRefiner::RefineInterior(Geometry::Type Geom,
1334  int Times)
1335 {
1336  IntegrationRule *ir = NULL;
1337 
1338  switch (Geom)
1339  {
1340  case Geometry::SEGMENT:
1341  {
1342  if (Times < 2)
1343  {
1344  return NULL;
1345  }
1346  ir = FindInIntPts(Geom, Times-1);
1347  if (ir) { return ir; }
1348 
1349  ir = new IntegrationRule(Times-1);
1350  for (int i = 1; i < Times; i++)
1351  {
1352  IntegrationPoint &ip = ir->IntPoint(i-1);
1353  ip.x = double(i) / Times;
1354  ip.y = ip.z = 0.0;
1355  }
1356  }
1357  break;
1358 
1359  case Geometry::TRIANGLE:
1360  {
1361  if (Times < 3)
1362  {
1363  return NULL;
1364  }
1365  ir = FindInIntPts(Geom, ((Times-1)*(Times-2))/2);
1366  if (ir) { return ir; }
1367 
1368  ir = new IntegrationRule(((Times-1)*(Times-2))/2);
1369  for (int k = 0, j = 1; j < Times-1; j++)
1370  for (int i = 1; i < Times-j; i++, k++)
1371  {
1372  IntegrationPoint &ip = ir->IntPoint(k);
1373  ip.x = double(i) / Times;
1374  ip.y = double(j) / Times;
1375  ip.z = 0.0;
1376  }
1377  }
1378  break;
1379 
1380  case Geometry::SQUARE:
1381  {
1382  if (Times < 2)
1383  {
1384  return NULL;
1385  }
1386  ir = FindInIntPts(Geom, (Times-1)*(Times-1));
1387  if (ir) { return ir; }
1388 
1389  ir = new IntegrationRule((Times-1)*(Times-1));
1390  for (int k = 0, j = 1; j < Times; j++)
1391  for (int i = 1; i < Times; i++, k++)
1392  {
1393  IntegrationPoint &ip = ir->IntPoint(k);
1394  ip.x = double(i) / Times;
1395  ip.y = double(j) / Times;
1396  ip.z = 0.0;
1397  }
1398  }
1399  break;
1400 
1401  default:
1402  mfem_error("GeometryRefiner::RefineInterior(...)");
1403  }
1404 
1405  MFEM_ASSERT(ir != NULL, "Failed to construct the refined IntegrationRule.");
1406  IntPts[Geom].Append(ir);
1407 
1408  return ir;
1409 }
1410 
1411 
1412 int GeometryRefiner::GetRefinementLevelFromPoints(Geometry::Type geom, int Npts)
1413 {
1414  switch (geom)
1415  {
1416  case Geometry::POINT:
1417  {
1418  return -1;
1419  }
1420  case Geometry::SEGMENT:
1421  {
1422  return Npts -1;
1423  }
1424  case Geometry::TRIANGLE:
1425  {
1426  for (int n = 0, np = 0; (n < 15) && (np < Npts) ; n++)
1427  {
1428  np = (n+1)*(n+2)/2;
1429  if (np == Npts) { return n; }
1430  }
1431  return -1;
1432  }
1433  case Geometry::SQUARE:
1434  {
1435  for (int n = 0, np = 0; (n < 15) && (np < Npts) ; n++)
1436  {
1437  np = (n+1)*(n+1);
1438  if (np == Npts) { return n; }
1439  }
1440  return -1;
1441  }
1442  case Geometry::CUBE:
1443  {
1444  for (int n = 0, np = 0; (n < 15) && (np < Npts) ; n++)
1445  {
1446  np = (n+1)*(n+1)*(n+1);
1447  if (np == Npts) { return n; }
1448  }
1449  return -1;
1450  }
1451  case Geometry::TETRAHEDRON:
1452  {
1453  for (int n = 0, np = 0; (n < 15) && (np < Npts) ; n++)
1454  {
1455  np = (n+3)*(n+2)*(n+1)/6;
1456  if (np == Npts) { return n; }
1457  }
1458  return -1;
1459  }
1460  case Geometry::PRISM:
1461  {
1462  for (int n = 0, np = 0; (n < 15) && (np < Npts) ; n++)
1463  {
1464  np = (n+1)*(n+1)*(n+2)/2;
1465  if (np == Npts) { return n; }
1466  }
1467  return -1;
1468  }
1469  default:
1470  {
1471  mfem_error("Non existing Geometry.");
1472  }
1473  }
1474 
1475  return -1;
1476 }
1477 
1478 
1479 int GeometryRefiner::GetRefinementLevelFromElems(Geometry::Type geom, int Nels)
1480 {
1481  switch (geom)
1482  {
1483  case Geometry::POINT:
1484  {
1485  return -1;
1486  }
1487  case Geometry::SEGMENT:
1488  {
1489  return Nels;
1490  }
1491  case Geometry::TRIANGLE:
1492  case Geometry::SQUARE:
1493  {
1494  for (int n = 0; (n < 15) && (n*n < Nels+1) ; n++)
1495  {
1496  if (n*n == Nels) { return n-1; }
1497  }
1498  return -1;
1499  }
1500  case Geometry::CUBE:
1501  case Geometry::TETRAHEDRON:
1502  case Geometry::PRISM:
1503  {
1504  for (int n = 0; (n < 15) && (n*n*n < Nels+1) ; n++)
1505  {
1506  if (n*n*n == Nels) { return n-1; }
1507  }
1508  return -1;
1509  }
1510  default:
1511  {
1512  mfem_error("Non existing Geometry.");
1513  }
1514  }
1515 
1516  return -1;
1517 }
1518 
1519 
1520 GeometryRefiner GlobGeometryRefiner;
1521 
1522 }
mfem::Array::Size
int Size() const
Return the 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::BasisType::GetNodalBasis
static int GetNodalBasis(int qpt_type)
Return the nodal BasisType corresponding to the Quadrature1D type.
Definition: fe.hpp:78
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:734
geom
const Geometry::Type geom
Definition: ex1.cpp:40
mfem::Geometry::GetRandomPoint
static void GetRandomPoint(int GeomType, IntegrationPoint &ip)
Get a random point in the reference element specified by GeomType.
Definition: geom.cpp:244
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)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition: eltrans.hpp:85
mfem::Quadrature1D::ClosedUniform
aka closed Newton-Cotes
Definition: intrules.hpp:296
mfem::Geometry::Constants< Geometry::PRISM >::FaceVert
static const int FaceVert[NumFaces][MaxFaceVert]
Definition: geom.hpp:226
mfem::Geometry::Constants< Geometry::PRISM >::Edges
static const int Edges[NumEdges][2]
Definition: geom.hpp:222
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:244
mfem::Geometry::GetVertices
const IntegrationRule * GetVertices(int GeomType)
Return an IntegrationRule consisting of all vertices of the given Geometry::Type, GeomType...
Definition: geom.cpp:225
mfem::GeometryRefiner::RefineInterior
const IntegrationRule * RefineInterior(Geometry::Type Geom, int Times)
Definition: geom.cpp:1333
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:7404
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:206
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:224
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:111
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:1520
mfem::RefinedGeometry::RefPts
IntegrationRule RefPts
Definition: geom.hpp:243
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:942
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::GeometryRefiner::GetRefinementLevelFromElems
virtual int GetRefinementLevelFromElems(Geometry::Type geom, int Npts)
Get the Refinement level based on number of elements.
Definition: geom.cpp:1479
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::Geometry::Constants< Geometry::TETRAHEDRON >::InvOrient
static const int InvOrient[NumOrient]
Definition: geom.hpp:196
mfem::WedgeFE
class H1_WedgeElement WedgeFE
mfem::CalcInverse
void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
Definition: densemat.cpp:2199
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:508
mfem::Geometry::GetPerfPointMat
void GetPerfPointMat(int GeomType, DenseMatrix &pm)
Definition: geom.cpp:661
mfem::TetrahedronFE
class Linear3DFiniteElement TetrahedronFE
Definition: fe.cpp:12837
mfem::TriangleFE
Linear2DFiniteElement TriangleFE
Definition: fe.cpp:12833
mfem::Geometry::Constants< Geometry::CUBE >::FaceVert
static const int FaceVert[NumFaces][MaxFaceVert]
Definition: geom.hpp:208
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::IsoparametricTransformation
A standard isoparametric element transformation.
Definition: eltrans.hpp:348
mfem::GeometryRefiner::GetRefinementLevelFromPoints
virtual int GetRefinementLevelFromPoints(Geometry::Type Geom, int Npts)
Get the Refinement level based on number of points.
Definition: geom.cpp:1412
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:1348
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:53
mfem::IsoparametricTransformation::SetFE
void SetFE(const FiniteElement *FE)
Set the element that will be used to compute the transformations.
Definition: eltrans.hpp:368
mfem::Geometry::Constants< Geometry::CUBE >::Edges
static const int Edges[NumEdges][2]
Definition: geom.hpp:204
mfem::Geometry::Constants< Geometry::TETRAHEDRON >::Orient
static const int Orient[NumOrient][NumVert]
Definition: geom.hpp:195
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:345
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::IsoparametricTransformation::GetPointMat
const DenseMatrix & GetPointMat() const
Return the stored point matrix.
Definition: eltrans.hpp:391
mfem::RefinedGeometry::RefGeoms
Array< int > RefGeoms
Definition: geom.hpp:244
mfem::Geometry::Constants< Geometry::TRIANGLE >::Orient
static const int Orient[NumOrient][NumVert]
Definition: geom.hpp:152
mfem::poly1d
Poly_1D poly1d
Definition: fe.cpp:7476
Generated on Fri Oct 30 2020 15:47:59 for MFEM by   doxygen 1.8.5