MFEM  v4.4.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-2022, 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 #include "../mesh/pyramid.hpp"
15 
16 namespace mfem
17 {
18 
19 const char *Geometry::Name[NumGeom] =
20 {
21  "Point", "Segment", "Triangle", "Square", "Tetrahedron", "Cube", "Prism",
22  "Pyramid"
23 };
24 
25 const double Geometry::Volume[NumGeom] =
26 { 1.0, 1.0, 0.5, 1.0, 1./6, 1.0, 0.5, 1./3 };
27 
28 Geometry::Geometry()
29 {
30  // Vertices for Geometry::POINT
31  GeomVert[0] = new IntegrationRule(1);
32  GeomVert[0]->IntPoint(0).x = 0.0;
33 
34  // Vertices for Geometry::SEGMENT
35  GeomVert[1] = new IntegrationRule(2);
36 
37  GeomVert[1]->IntPoint(0).x = 0.0;
38  GeomVert[1]->IntPoint(1).x = 1.0;
39 
40  // Vertices for Geometry::TRIANGLE
41  GeomVert[2] = new IntegrationRule(3);
42 
43  GeomVert[2]->IntPoint(0).x = 0.0;
44  GeomVert[2]->IntPoint(0).y = 0.0;
45 
46  GeomVert[2]->IntPoint(1).x = 1.0;
47  GeomVert[2]->IntPoint(1).y = 0.0;
48 
49  GeomVert[2]->IntPoint(2).x = 0.0;
50  GeomVert[2]->IntPoint(2).y = 1.0;
51 
52  // Vertices for Geometry::SQUARE
53  GeomVert[3] = new IntegrationRule(4);
54 
55  GeomVert[3]->IntPoint(0).x = 0.0;
56  GeomVert[3]->IntPoint(0).y = 0.0;
57 
58  GeomVert[3]->IntPoint(1).x = 1.0;
59  GeomVert[3]->IntPoint(1).y = 0.0;
60 
61  GeomVert[3]->IntPoint(2).x = 1.0;
62  GeomVert[3]->IntPoint(2).y = 1.0;
63 
64  GeomVert[3]->IntPoint(3).x = 0.0;
65  GeomVert[3]->IntPoint(3).y = 1.0;
66 
67  // Vertices for Geometry::TETRAHEDRON
68  GeomVert[4] = new IntegrationRule(4);
69  GeomVert[4]->IntPoint(0).x = 0.0;
70  GeomVert[4]->IntPoint(0).y = 0.0;
71  GeomVert[4]->IntPoint(0).z = 0.0;
72 
73  GeomVert[4]->IntPoint(1).x = 1.0;
74  GeomVert[4]->IntPoint(1).y = 0.0;
75  GeomVert[4]->IntPoint(1).z = 0.0;
76 
77  GeomVert[4]->IntPoint(2).x = 0.0;
78  GeomVert[4]->IntPoint(2).y = 1.0;
79  GeomVert[4]->IntPoint(2).z = 0.0;
80 
81  GeomVert[4]->IntPoint(3).x = 0.0;
82  GeomVert[4]->IntPoint(3).y = 0.0;
83  GeomVert[4]->IntPoint(3).z = 1.0;
84 
85  // Vertices for Geometry::CUBE
86  GeomVert[5] = new IntegrationRule(8);
87 
88  GeomVert[5]->IntPoint(0).x = 0.0;
89  GeomVert[5]->IntPoint(0).y = 0.0;
90  GeomVert[5]->IntPoint(0).z = 0.0;
91 
92  GeomVert[5]->IntPoint(1).x = 1.0;
93  GeomVert[5]->IntPoint(1).y = 0.0;
94  GeomVert[5]->IntPoint(1).z = 0.0;
95 
96  GeomVert[5]->IntPoint(2).x = 1.0;
97  GeomVert[5]->IntPoint(2).y = 1.0;
98  GeomVert[5]->IntPoint(2).z = 0.0;
99 
100  GeomVert[5]->IntPoint(3).x = 0.0;
101  GeomVert[5]->IntPoint(3).y = 1.0;
102  GeomVert[5]->IntPoint(3).z = 0.0;
103 
104  GeomVert[5]->IntPoint(4).x = 0.0;
105  GeomVert[5]->IntPoint(4).y = 0.0;
106  GeomVert[5]->IntPoint(4).z = 1.0;
107 
108  GeomVert[5]->IntPoint(5).x = 1.0;
109  GeomVert[5]->IntPoint(5).y = 0.0;
110  GeomVert[5]->IntPoint(5).z = 1.0;
111 
112  GeomVert[5]->IntPoint(6).x = 1.0;
113  GeomVert[5]->IntPoint(6).y = 1.0;
114  GeomVert[5]->IntPoint(6).z = 1.0;
115 
116  GeomVert[5]->IntPoint(7).x = 0.0;
117  GeomVert[5]->IntPoint(7).y = 1.0;
118  GeomVert[5]->IntPoint(7).z = 1.0;
119 
120  // Vertices for Geometry::PRISM
121  GeomVert[6] = new IntegrationRule(6);
122  GeomVert[6]->IntPoint(0).x = 0.0;
123  GeomVert[6]->IntPoint(0).y = 0.0;
124  GeomVert[6]->IntPoint(0).z = 0.0;
125 
126  GeomVert[6]->IntPoint(1).x = 1.0;
127  GeomVert[6]->IntPoint(1).y = 0.0;
128  GeomVert[6]->IntPoint(1).z = 0.0;
129 
130  GeomVert[6]->IntPoint(2).x = 0.0;
131  GeomVert[6]->IntPoint(2).y = 1.0;
132  GeomVert[6]->IntPoint(2).z = 0.0;
133 
134  GeomVert[6]->IntPoint(3).x = 0.0;
135  GeomVert[6]->IntPoint(3).y = 0.0;
136  GeomVert[6]->IntPoint(3).z = 1.0;
137 
138  GeomVert[6]->IntPoint(4).x = 1.0;
139  GeomVert[6]->IntPoint(4).y = 0.0;
140  GeomVert[6]->IntPoint(4).z = 1.0;
141 
142  GeomVert[6]->IntPoint(5).x = 0.0;
143  GeomVert[6]->IntPoint(5).y = 1.0;
144  GeomVert[6]->IntPoint(5).z = 1.0;
145 
146  // Vertices for Geometry::PYRAMID
147  GeomVert[7] = new IntegrationRule(5);
148  GeomVert[7]->IntPoint(0).x = 0.0;
149  GeomVert[7]->IntPoint(0).y = 0.0;
150  GeomVert[7]->IntPoint(0).z = 0.0;
151 
152  GeomVert[7]->IntPoint(1).x = 1.0;
153  GeomVert[7]->IntPoint(1).y = 0.0;
154  GeomVert[7]->IntPoint(1).z = 0.0;
155 
156  GeomVert[7]->IntPoint(2).x = 1.0;
157  GeomVert[7]->IntPoint(2).y = 1.0;
158  GeomVert[7]->IntPoint(2).z = 0.0;
159 
160  GeomVert[7]->IntPoint(3).x = 0.0;
161  GeomVert[7]->IntPoint(3).y = 1.0;
162  GeomVert[7]->IntPoint(3).z = 0.0;
163 
164  GeomVert[7]->IntPoint(4).x = 0.0;
165  GeomVert[7]->IntPoint(4).y = 0.0;
166  GeomVert[7]->IntPoint(4).z = 1.0;
167 
168  GeomCenter[POINT].x = 0.0;
169  GeomCenter[POINT].y = 0.0;
170  GeomCenter[POINT].z = 0.0;
171 
172  GeomCenter[SEGMENT].x = 0.5;
173  GeomCenter[SEGMENT].y = 0.0;
174  GeomCenter[SEGMENT].z = 0.0;
175 
176  GeomCenter[TRIANGLE].x = 1.0 / 3.0;
177  GeomCenter[TRIANGLE].y = 1.0 / 3.0;
178  GeomCenter[TRIANGLE].z = 0.0;
179 
180  GeomCenter[SQUARE].x = 0.5;
181  GeomCenter[SQUARE].y = 0.5;
182  GeomCenter[SQUARE].z = 0.0;
183 
184  GeomCenter[TETRAHEDRON].x = 0.25;
185  GeomCenter[TETRAHEDRON].y = 0.25;
186  GeomCenter[TETRAHEDRON].z = 0.25;
187 
188  GeomCenter[CUBE].x = 0.5;
189  GeomCenter[CUBE].y = 0.5;
190  GeomCenter[CUBE].z = 0.5;
191 
192  GeomCenter[PRISM].x = 1.0 / 3.0;
193  GeomCenter[PRISM].y = 1.0 / 3.0;
194  GeomCenter[PRISM].z = 0.5;
195 
196  GeomCenter[PYRAMID].x = 0.375;
197  GeomCenter[PYRAMID].y = 0.375;
198  GeomCenter[PYRAMID].z = 0.25;
199 
200  GeomToPerfGeomJac[POINT] = NULL;
201  GeomToPerfGeomJac[SEGMENT] = new DenseMatrix(1);
202  GeomToPerfGeomJac[TRIANGLE] = new DenseMatrix(2);
203  GeomToPerfGeomJac[SQUARE] = new DenseMatrix(2);
204  GeomToPerfGeomJac[TETRAHEDRON] = new DenseMatrix(3);
205  GeomToPerfGeomJac[CUBE] = new DenseMatrix(3);
206  GeomToPerfGeomJac[PRISM] = new DenseMatrix(3);
207  GeomToPerfGeomJac[PYRAMID] = new DenseMatrix(3);
208 
209  PerfGeomToGeomJac[POINT] = NULL;
210  PerfGeomToGeomJac[SEGMENT] = NULL;
211  PerfGeomToGeomJac[TRIANGLE] = new DenseMatrix(2);
212  PerfGeomToGeomJac[SQUARE] = NULL;
213  PerfGeomToGeomJac[TETRAHEDRON] = new DenseMatrix(3);
214  PerfGeomToGeomJac[CUBE] = NULL;
215  PerfGeomToGeomJac[PRISM] = new DenseMatrix(3);
216  PerfGeomToGeomJac[PYRAMID] = new DenseMatrix(3);
217 
218  GeomToPerfGeomJac[SEGMENT]->Diag(1.0, 1);
219  {
220  IsoparametricTransformation tri_T;
221  tri_T.SetFE(&TriangleFE);
222  GetPerfPointMat (TRIANGLE, tri_T.GetPointMat());
223  tri_T.SetIntPoint(&GeomCenter[TRIANGLE]);
224  *GeomToPerfGeomJac[TRIANGLE] = tri_T.Jacobian();
225  CalcInverse(tri_T.Jacobian(), *PerfGeomToGeomJac[TRIANGLE]);
226  }
227  GeomToPerfGeomJac[SQUARE]->Diag(1.0, 2);
228  {
229  IsoparametricTransformation tet_T;
230  tet_T.SetFE(&TetrahedronFE);
231  GetPerfPointMat (TETRAHEDRON, tet_T.GetPointMat());
232  tet_T.SetIntPoint(&GeomCenter[TETRAHEDRON]);
233  *GeomToPerfGeomJac[TETRAHEDRON] = tet_T.Jacobian();
234  CalcInverse(tet_T.Jacobian(), *PerfGeomToGeomJac[TETRAHEDRON]);
235  }
236  GeomToPerfGeomJac[CUBE]->Diag(1.0, 3);
237  {
238  IsoparametricTransformation pri_T;
239  pri_T.SetFE(&WedgeFE);
240  GetPerfPointMat (PRISM, pri_T.GetPointMat());
241  pri_T.SetIntPoint(&GeomCenter[PRISM]);
242  *GeomToPerfGeomJac[PRISM] = pri_T.Jacobian();
243  CalcInverse(pri_T.Jacobian(), *PerfGeomToGeomJac[PRISM]);
244  }
245  {
246  IsoparametricTransformation pyr_T;
247  pyr_T.SetFE(&PyramidFE);
248  GetPerfPointMat (PYRAMID, pyr_T.GetPointMat());
249  pyr_T.SetIntPoint(&GeomCenter[PYRAMID]);
250  *GeomToPerfGeomJac[PYRAMID] = pyr_T.Jacobian();
251  CalcInverse(pyr_T.Jacobian(), *PerfGeomToGeomJac[PYRAMID]);
252  }
253 }
254 
255 Geometry::~Geometry()
256 {
257  for (int i = 0; i < NumGeom; i++)
258  {
259  delete PerfGeomToGeomJac[i];
260  delete GeomToPerfGeomJac[i];
261  delete GeomVert[i];
262  }
263 }
264 
265 const IntegrationRule * Geometry::GetVertices(int GeomType)
266 {
267  switch (GeomType)
268  {
269  case Geometry::POINT: return GeomVert[0];
270  case Geometry::SEGMENT: return GeomVert[1];
271  case Geometry::TRIANGLE: return GeomVert[2];
272  case Geometry::SQUARE: return GeomVert[3];
273  case Geometry::TETRAHEDRON: return GeomVert[4];
274  case Geometry::CUBE: return GeomVert[5];
275  case Geometry::PRISM: return GeomVert[6];
276  case Geometry::PYRAMID: return GeomVert[7];
277  default:
278  mfem_error ("Geometry::GetVertices(...)");
279  }
280  // make some compilers happy.
281  return GeomVert[0];
282 }
283 
284 // static method
285 void Geometry::GetRandomPoint(int GeomType, IntegrationPoint &ip)
286 {
287  switch (GeomType)
288  {
289  case Geometry::POINT:
290  ip.x = 0.0;
291  break;
292  case Geometry::SEGMENT:
293  ip.x = double(rand()) / RAND_MAX;
294  break;
295  case Geometry::TRIANGLE:
296  ip.x = double(rand()) / RAND_MAX;
297  ip.y = double(rand()) / RAND_MAX;
298  if (ip.x + ip.y > 1.0)
299  {
300  ip.x = 1.0 - ip.x;
301  ip.y = 1.0 - ip.y;
302  }
303  break;
304  case Geometry::SQUARE:
305  ip.x = double(rand()) / RAND_MAX;
306  ip.y = double(rand()) / RAND_MAX;
307  break;
308  case Geometry::TETRAHEDRON:
309  ip.x = double(rand()) / RAND_MAX;
310  ip.y = double(rand()) / RAND_MAX;
311  ip.z = double(rand()) / RAND_MAX;
312  // map to the triangular prism obtained by extruding the reference
313  // triangle in z direction
314  if (ip.x + ip.y > 1.0)
315  {
316  ip.x = 1.0 - ip.x;
317  ip.y = 1.0 - ip.y;
318  }
319  // split the prism into 3 parts: 1 is the reference tet, and the
320  // other two tets (as given below) are mapped to the reference tet
321  if (ip.x + ip.z > 1.0)
322  {
323  // tet with vertices: (0,0,1),(1,0,1),(0,1,1),(1,0,0)
324  ip.x = ip.x + ip.z - 1.0;
325  // ip.y = ip.y;
326  ip.z = 1.0 - ip.z;
327  // mapped to: (0,0,0),(1,0,0),(0,1,0),(0,0,1)
328  }
329  else if (ip.x + ip.y + ip.z > 1.0)
330  {
331  // tet with vertices: (0,1,1),(0,1,0),(0,0,1),(1,0,0)
332  double x = ip.x;
333  ip.x = 1.0 - x - ip.z;
334  ip.y = 1.0 - x - ip.y;
335  ip.z = x;
336  // mapped to: (0,0,0),(1,0,0),(0,1,0),(0,0,1)
337  }
338  break;
339  case Geometry::CUBE:
340  ip.x = double(rand()) / RAND_MAX;
341  ip.y = double(rand()) / RAND_MAX;
342  ip.z = double(rand()) / RAND_MAX;
343  break;
344  case Geometry::PRISM:
345  ip.x = double(rand()) / RAND_MAX;
346  ip.y = double(rand()) / RAND_MAX;
347  ip.z = double(rand()) / RAND_MAX;
348  if (ip.x + ip.y > 1.0)
349  {
350  ip.x = 1.0 - ip.x;
351  ip.y = 1.0 - ip.y;
352  }
353  break;
354  case Geometry::PYRAMID:
355  ip.x = double(rand()) / RAND_MAX;
356  ip.y = double(rand()) / RAND_MAX;
357  ip.z = double(rand()) / RAND_MAX;
358  if (ip.x + ip.z > 1.0 && ip.y < ip.x)
359  {
360  double x = ip.x;
361  ip.x = ip.y;
362  ip.y = 1.0 - ip.z;
363  ip.z = 1.0 - x;
364  }
365  else if (ip.y + ip.z > 1.0)
366  {
367  double z = ip.z;
368  ip.z = 1.0 - ip.y;
369  ip.y = ip.x;
370  ip.x = 1.0 - z;
371  }
372  break;
373  default:
374  MFEM_ABORT("Unknown type of reference element!");
375  }
376 }
377 
378 
379 namespace internal
380 {
381 
382 // Fuzzy equality operator with absolute tolerance eps.
383 inline bool NearlyEqual(double x, double y, double eps)
384 {
385  return std::abs(x-y) <= eps;
386 }
387 
388 // Fuzzy greater than comparison operator with absolute tolerance eps.
389 // Returns true when x is greater than y by at least eps.
390 inline bool FuzzyGT(double x, double y, double eps)
391 {
392  return (x > y + eps);
393 }
394 
395 // Fuzzy less than comparison operator with absolute tolerance eps.
396 // Returns true when x is less than y by at least eps.
397 inline bool FuzzyLT(double x, double y, double eps)
398 {
399  return (x < y - eps);
400 }
401 
402 }
403 
404 // static method
405 bool Geometry::CheckPoint(int GeomType, const IntegrationPoint &ip)
406 {
407  switch (GeomType)
408  {
409  case Geometry::POINT:
410  if (ip.x != 0.0) { return false; }
411  break;
412  case Geometry::SEGMENT:
413  if (ip.x < 0.0 || ip.x > 1.0) { return false; }
414  break;
415  case Geometry::TRIANGLE:
416  if (ip.x < 0.0 || ip.y < 0.0 || ip.x+ip.y > 1.0) { return false; }
417  break;
418  case Geometry::SQUARE:
419  if (ip.x < 0.0 || ip.x > 1.0 || ip.y < 0.0 || ip.y > 1.0)
420  { return false; }
421  break;
422  case Geometry::TETRAHEDRON:
423  if (ip.x < 0.0 || ip.y < 0.0 || ip.z < 0.0 ||
424  ip.x+ip.y+ip.z > 1.0) { return false; }
425  break;
426  case Geometry::CUBE:
427  if (ip.x < 0.0 || ip.x > 1.0 || ip.y < 0.0 || ip.y > 1.0 ||
428  ip.z < 0.0 || ip.z > 1.0) { return false; }
429  break;
430  case Geometry::PRISM:
431  if (ip.x < 0.0 || ip.y < 0.0 || ip.x+ip.y > 1.0 ||
432  ip.z < 0.0 || ip.z > 1.0) { return false; }
433  break;
434  case Geometry::PYRAMID:
435  if (ip.x < 0.0 || ip.y < 0.0 || ip.x+ip.z > 1.0 || ip.y+ip.z > 1.0 ||
436  ip.z < 0.0 || ip.z > 1.0) { return false; }
437  break;
438  default:
439  MFEM_ABORT("Unknown type of reference element!");
440  }
441  return true;
442 }
443 
444 bool Geometry::CheckPoint(int GeomType, const IntegrationPoint &ip, double eps)
445 {
446  switch (GeomType)
447  {
448  case Geometry::POINT:
449  if (! internal::NearlyEqual(ip.x, 0.0, eps))
450  {
451  return false;
452  }
453  break;
454  case Geometry::SEGMENT:
455  if ( internal::FuzzyLT(ip.x, 0.0, eps)
456  || internal::FuzzyGT(ip.x, 1.0, eps) )
457  {
458  return false;
459  }
460  break;
461  case Geometry::TRIANGLE:
462  if ( internal::FuzzyLT(ip.x, 0.0, eps)
463  || internal::FuzzyLT(ip.y, 0.0, eps)
464  || internal::FuzzyGT(ip.x+ip.y, 1.0, eps) )
465  {
466  return false;
467  }
468  break;
469  case Geometry::SQUARE:
470  if ( internal::FuzzyLT(ip.x, 0.0, eps)
471  || internal::FuzzyGT(ip.x, 1.0, eps)
472  || internal::FuzzyLT(ip.y, 0.0, eps)
473  || internal::FuzzyGT(ip.y, 1.0, eps) )
474  {
475  return false;
476  }
477  break;
478  case Geometry::TETRAHEDRON:
479  if ( internal::FuzzyLT(ip.x, 0.0, eps)
480  || internal::FuzzyLT(ip.y, 0.0, eps)
481  || internal::FuzzyLT(ip.z, 0.0, eps)
482  || internal::FuzzyGT(ip.x+ip.y+ip.z, 1.0, eps) )
483  {
484  return false;
485  }
486  break;
487  case Geometry::CUBE:
488  if ( internal::FuzzyLT(ip.x, 0.0, eps)
489  || internal::FuzzyGT(ip.x, 1.0, eps)
490  || internal::FuzzyLT(ip.y, 0.0, eps)
491  || internal::FuzzyGT(ip.y, 1.0, eps)
492  || internal::FuzzyLT(ip.z, 0.0, eps)
493  || internal::FuzzyGT(ip.z, 1.0, eps) )
494  {
495  return false;
496  }
497  break;
498  case Geometry::PRISM:
499  if ( internal::FuzzyLT(ip.x, 0.0, eps)
500  || internal::FuzzyLT(ip.y, 0.0, eps)
501  || internal::FuzzyGT(ip.x+ip.y, 1.0, eps)
502  || internal::FuzzyLT(ip.z, 0.0, eps)
503  || internal::FuzzyGT(ip.z, 1.0, eps) )
504  {
505  return false;
506  }
507  break;
508  case Geometry::PYRAMID:
509  if (internal::FuzzyLT(ip.x, 0.0, eps)
510  || internal::FuzzyLT(ip.y, 0.0, eps)
511  || internal::FuzzyGT(ip.x+ip.z, 1.0, eps)
512  || internal::FuzzyGT(ip.y+ip.z, 1.0, eps)
513  || internal::FuzzyLT(ip.z, 0.0, eps)
514  || internal::FuzzyGT(ip.z, 1.0, eps) )
515  {
516  return false;
517  }
518  break;
519  default:
520  MFEM_ABORT("Unknown type of reference element!");
521  }
522  return true;
523 }
524 
525 
526 namespace internal
527 {
528 
529 template <int N, int dim>
530 inline bool IntersectSegment(double lbeg[N], double lend[N],
531  IntegrationPoint &end)
532 {
533  double t = 1.0;
534  bool out = false;
535  for (int i = 0; i < N; i++)
536  {
537  lbeg[i] = std::max(lbeg[i], 0.0); // remove round-off
538  if (lend[i] < 0.0)
539  {
540  out = true;
541  t = std::min(t, lbeg[i]/(lbeg[i]-lend[i]));
542  }
543  }
544  if (out)
545  {
546  if (dim >= 1) { end.x = t*lend[0] + (1.0-t)*lbeg[0]; }
547  if (dim >= 2) { end.y = t*lend[1] + (1.0-t)*lbeg[1]; }
548  if (dim >= 3) { end.z = t*lend[2] + (1.0-t)*lbeg[2]; }
549  return false;
550  }
551  return true;
552 }
553 
554 inline bool ProjectTriangle(double &x, double &y)
555 {
556  if (x < 0.0)
557  {
558  x = 0.0;
559  if (y < 0.0) { y = 0.0; }
560  else if (y > 1.0) { y = 1.0; }
561  return false;
562  }
563  if (y < 0.0)
564  {
565  if (x > 1.0) { x = 1.0; }
566  y = 0.0;
567  return false;
568  }
569  const double l3 = 1.0-x-y;
570  if (l3 < 0.0)
571  {
572  if (y - x > 1.0) { x = 0.0; y = 1.0; }
573  else if (y - x < -1.0) { x = 1.0; y = 0.0; }
574  else { x += l3/2; y += l3/2; }
575  return false;
576  }
577  return true;
578 }
579 
580 }
581 
582 // static method
583 bool Geometry::ProjectPoint(int GeomType, const IntegrationPoint &beg,
584  IntegrationPoint &end)
585 {
586  switch (GeomType)
587  {
588  case Geometry::POINT:
589  {
590  if (end.x != 0.0) { end.x = 0.0; return false; }
591  break;
592  }
593  case Geometry::SEGMENT:
594  {
595  if (end.x < 0.0) { end.x = 0.0; return false; }
596  if (end.x > 1.0) { end.x = 1.0; return false; }
597  break;
598  }
599  case Geometry::TRIANGLE:
600  {
601  double lend[3] = { end.x, end.y, 1.0-end.x-end.y };
602  double lbeg[3] = { beg.x, beg.y, 1.0-beg.x-beg.y };
603  return internal::IntersectSegment<3,2>(lbeg, lend, end);
604  }
605  case Geometry::SQUARE:
606  {
607  double lend[4] = { end.x, end.y, 1.0-end.x, 1.0-end.y };
608  double lbeg[4] = { beg.x, beg.y, 1.0-beg.x, 1.0-beg.y };
609  return internal::IntersectSegment<4,2>(lbeg, lend, end);
610  }
611  case Geometry::TETRAHEDRON:
612  {
613  double lend[4] = { end.x, end.y, end.z, 1.0-end.x-end.y-end.z };
614  double lbeg[4] = { beg.x, beg.y, beg.z, 1.0-beg.x-beg.y-beg.z };
615  return internal::IntersectSegment<4,3>(lbeg, lend, end);
616  }
617  case Geometry::CUBE:
618  {
619  double lend[6] = { end.x, end.y, end.z,
620  1.0-end.x, 1.0-end.y, 1.0-end.z
621  };
622  double lbeg[6] = { beg.x, beg.y, beg.z,
623  1.0-beg.x, 1.0-beg.y, 1.0-beg.z
624  };
625  return internal::IntersectSegment<6,3>(lbeg, lend, end);
626  }
627  case Geometry::PRISM:
628  {
629  double lend[5] = { end.x, end.y, end.z, 1.0-end.x-end.y, 1.0-end.z };
630  double lbeg[5] = { beg.x, beg.y, beg.z, 1.0-beg.x-beg.y, 1.0-beg.z };
631  return internal::IntersectSegment<5,3>(lbeg, lend, end);
632  }
633  case Geometry::PYRAMID:
634  {
635  double lend[6] = { end.x, end.y, end.z,
636  1.0-end.x-end.z, 1.0-end.y-end.z, 1.0-end.z
637  };
638  double lbeg[6] = { beg.x, beg.y, beg.z,
639  1.0-beg.x-beg.z, 1.0-beg.y-beg.z, 1.0-beg.z
640  };
641  return internal::IntersectSegment<6,3>(lbeg, lend, end);
642  }
643  default:
644  MFEM_ABORT("Unknown type of reference element!");
645  }
646  return true;
647 }
648 
649 // static method
650 bool Geometry::ProjectPoint(int GeomType, IntegrationPoint &ip)
651 {
652  // If ip is outside the element, replace it with the point on the boundary
653  // that is closest to the original ip and return false; otherwise, return
654  // true without changing ip.
655 
656  switch (GeomType)
657  {
658  case SEGMENT:
659  {
660  if (ip.x < 0.0) { ip.x = 0.0; return false; }
661  else if (ip.x > 1.0) { ip.x = 1.0; return false; }
662  return true;
663  }
664 
665  case TRIANGLE:
666  {
667  return internal::ProjectTriangle(ip.x, ip.y);
668  }
669 
670  case SQUARE:
671  {
672  bool in_x, in_y;
673  if (ip.x < 0.0) { in_x = false; ip.x = 0.0; }
674  else if (ip.x > 1.0) { in_x = false; ip.x = 1.0; }
675  else { in_x = true; }
676  if (ip.y < 0.0) { in_y = false; ip.y = 0.0; }
677  else if (ip.y > 1.0) { in_y = false; ip.y = 1.0; }
678  else { in_y = true; }
679  return in_x && in_y;
680  }
681 
682  case TETRAHEDRON:
683  {
684  if (ip.z < 0.0)
685  {
686  ip.z = 0.0;
687  internal::ProjectTriangle(ip.x, ip.y);
688  return false;
689  }
690  if (ip.y < 0.0)
691  {
692  ip.y = 0.0;
693  internal::ProjectTriangle(ip.x, ip.z);
694  return false;
695  }
696  if (ip.x < 0.0)
697  {
698  ip.x = 0.0;
699  internal::ProjectTriangle(ip.y, ip.z);
700  return false;
701  }
702  const double l4 = 1.0-ip.x-ip.y-ip.z;
703  if (l4 < 0.0)
704  {
705  const double l4_3 = l4/3;
706  ip.x += l4_3;
707  ip.y += l4_3;
708  internal::ProjectTriangle(ip.x, ip.y);
709  ip.z = 1.0-ip.x-ip.y;
710  return false;
711  }
712  return true;
713  }
714 
715  case CUBE:
716  {
717  bool in_x, in_y, in_z;
718  if (ip.x < 0.0) { in_x = false; ip.x = 0.0; }
719  else if (ip.x > 1.0) { in_x = false; ip.x = 1.0; }
720  else { in_x = true; }
721  if (ip.y < 0.0) { in_y = false; ip.y = 0.0; }
722  else if (ip.y > 1.0) { in_y = false; ip.y = 1.0; }
723  else { in_y = true; }
724  if (ip.z < 0.0) { in_z = false; ip.z = 0.0; }
725  else if (ip.z > 1.0) { in_z = false; ip.z = 1.0; }
726  else { in_z = true; }
727  return in_x && in_y && in_z;
728  }
729 
730  case PRISM:
731  {
732  bool in_tri, in_z;
733  in_tri = internal::ProjectTriangle(ip.x, ip.y);
734  if (ip.z < 0.0) { in_z = false; ip.z = 0.0; }
735  else if (ip.z > 1.0) { in_z = false; ip.z = 1.0; }
736  else { in_z = true; }
737  return in_tri && in_z;
738  }
739 
740  case PYRAMID:
741  {
742  if (ip.x < 0.0)
743  {
744  ip.x = 0.0;
745  internal::ProjectTriangle(ip.y, ip.z);
746  return false;
747  }
748  if (ip.y < 0.0)
749  {
750  ip.y = 0.0;
751  internal::ProjectTriangle(ip.x, ip.z);
752  return false;
753  }
754  if (ip.z < 0.0)
755  {
756  ip.z = 0.0;
757  if (ip.x > 1.0) { ip.x = 1.0; }
758  if (ip.y > 1.0) { ip.y = 1.0; }
759  return false;
760  }
761  if (ip.x >= ip.y)
762  {
763  bool in_y = true;
764  bool in_tri = internal::ProjectTriangle(ip.x, ip.z);
765  if (ip.y > ip.z) { in_y = false; ip.y = ip.z; }
766  return in_tri && in_y;
767  }
768  else
769  {
770  bool in_x = true;
771  bool in_tri = internal::ProjectTriangle(ip.y, ip.z);
772  if (ip.x > ip.z) { in_x = false; ip.x = ip.z; }
773  return in_tri && in_x;
774  }
775  }
776 
777  default:
778  MFEM_ABORT("Reference element type is not supported!");
779  }
780  return true;
781 }
782 
783 void Geometry::GetPerfPointMat(int GeomType, DenseMatrix &pm)
784 {
785  switch (GeomType)
786  {
787  case Geometry::SEGMENT:
788  {
789  pm.SetSize (1, 2);
790  pm(0,0) = 0.0;
791  pm(0,1) = 1.0;
792  }
793  break;
794 
795  case Geometry::TRIANGLE:
796  {
797  pm.SetSize (2, 3);
798  pm(0,0) = 0.0; pm(1,0) = 0.0;
799  pm(0,1) = 1.0; pm(1,1) = 0.0;
800  pm(0,2) = 0.5; pm(1,2) = 0.86602540378443864676;
801  }
802  break;
803 
804  case Geometry::SQUARE:
805  {
806  pm.SetSize (2, 4);
807  pm(0,0) = 0.0; pm(1,0) = 0.0;
808  pm(0,1) = 1.0; pm(1,1) = 0.0;
809  pm(0,2) = 1.0; pm(1,2) = 1.0;
810  pm(0,3) = 0.0; pm(1,3) = 1.0;
811  }
812  break;
813 
814  case Geometry::TETRAHEDRON:
815  {
816  pm.SetSize (3, 4);
817  pm(0,0) = 0.0; pm(1,0) = 0.0; pm(2,0) = 0.0;
818  pm(0,1) = 1.0; pm(1,1) = 0.0; pm(2,1) = 0.0;
819  pm(0,2) = 0.5; pm(1,2) = 0.86602540378443864676; pm(2,2) = 0.0;
820  pm(0,3) = 0.5; pm(1,3) = 0.28867513459481288225;
821  pm(2,3) = 0.81649658092772603273;
822  }
823  break;
824 
825  case Geometry::CUBE:
826  {
827  pm.SetSize (3, 8);
828  pm(0,0) = 0.0; pm(1,0) = 0.0; pm(2,0) = 0.0;
829  pm(0,1) = 1.0; pm(1,1) = 0.0; pm(2,1) = 0.0;
830  pm(0,2) = 1.0; pm(1,2) = 1.0; pm(2,2) = 0.0;
831  pm(0,3) = 0.0; pm(1,3) = 1.0; pm(2,3) = 0.0;
832  pm(0,4) = 0.0; pm(1,4) = 0.0; pm(2,4) = 1.0;
833  pm(0,5) = 1.0; pm(1,5) = 0.0; pm(2,5) = 1.0;
834  pm(0,6) = 1.0; pm(1,6) = 1.0; pm(2,6) = 1.0;
835  pm(0,7) = 0.0; pm(1,7) = 1.0; pm(2,7) = 1.0;
836  }
837  break;
838 
839  case Geometry::PRISM:
840  {
841  pm.SetSize (3, 6);
842  pm(0,0) = 0.0; pm(1,0) = 0.0; pm(2,0) = 0.0;
843  pm(0,1) = 1.0; pm(1,1) = 0.0; pm(2,1) = 0.0;
844  pm(0,2) = 0.5; pm(1,2) = 0.86602540378443864676; pm(2,2) = 0.0;
845  pm(0,3) = 0.0; pm(1,3) = 0.0; pm(2,3) = 1.0;
846  pm(0,4) = 1.0; pm(1,4) = 0.0; pm(2,4) = 1.0;
847  pm(0,5) = 0.5; pm(1,5) = 0.86602540378443864676; pm(2,5) = 1.0;
848  }
849  break;
850 
851  case Geometry::PYRAMID:
852  {
853  pm.SetSize (3, 5);
854  pm(0,0) = 0.0; pm(1,0) = 0.0; pm(2,0) = 0.0;
855  pm(0,1) = 1.0; pm(1,1) = 0.0; pm(2,1) = 0.0;
856  pm(0,2) = 1.0; pm(1,2) = 1.0; pm(2,2) = 0.0;
857  pm(0,3) = 0.0; pm(1,3) = 1.0; pm(2,3) = 0.0;
858  pm(0,4) = 0.5; pm(1,4) = 0.5; pm(2,4) = 0.7071067811865475;
859  }
860  break;
861 
862  default:
863  mfem_error ("Geometry::GetPerfPointMat (...)");
864  }
865 }
866 
867 void Geometry::JacToPerfJac(int GeomType, const DenseMatrix &J,
868  DenseMatrix &PJ) const
869 {
870  if (PerfGeomToGeomJac[GeomType])
871  {
872  Mult(J, *PerfGeomToGeomJac[GeomType], PJ);
873  }
874  else
875  {
876  PJ = J;
877  }
878 }
879 
880 const int Geometry::NumBdrArray[NumGeom] = { 0, 2, 3, 4, 4, 6, 5, 5 };
881 const int Geometry::Dimension[NumGeom] = { 0, 1, 2, 2, 3, 3, 3, 3 };
882 const int Geometry::DimStart[MaxDim+2] =
883 { POINT, SEGMENT, TRIANGLE, TETRAHEDRON, NUM_GEOMETRIES };
884 const int Geometry::NumVerts[NumGeom] = { 1, 2, 3, 4, 4, 8, 6, 5 };
885 const int Geometry::NumEdges[NumGeom] = { 0, 1, 3, 4, 6, 12, 9, 8 };
886 const int Geometry::NumFaces[NumGeom] = { 0, 0, 1, 1, 4, 6, 5, 5 };
887 
888 const int Geometry::
889 Constants<Geometry::POINT>::Orient[1][1] = {{0}};
890 const int Geometry::
891 Constants<Geometry::POINT>::InvOrient[1] = {0};
892 
893 const int Geometry::
894 Constants<Geometry::SEGMENT>::Edges[1][2] = { {0, 1} };
895 const int Geometry::
896 Constants<Geometry::SEGMENT>::Orient[2][2] = { {0, 1}, {1, 0} };
897 const int Geometry::
898 Constants<Geometry::SEGMENT>::InvOrient[2] = { 0, 1 };
899 
900 const int Geometry::
901 Constants<Geometry::TRIANGLE>::Edges[3][2] = {{0, 1}, {1, 2}, {2, 0}};
902 const int Geometry::
903 Constants<Geometry::TRIANGLE>::VertToVert::I[3] = {0, 2, 3};
904 const int Geometry::
905 Constants<Geometry::TRIANGLE>::VertToVert::J[3][2] = {{1, 0}, {2, -3}, {2, 1}};
906 const int Geometry::
907 Constants<Geometry::TRIANGLE>::FaceVert[1][3] = {{0, 1, 2}};
908 const int Geometry::
909 Constants<Geometry::TRIANGLE>::Orient[6][3] =
910 {
911  {0, 1, 2}, {1, 0, 2}, {2, 0, 1},
912  {2, 1, 0}, {1, 2, 0}, {0, 2, 1}
913 };
914 const int Geometry::
915 Constants<Geometry::TRIANGLE>::InvOrient[6] = { 0, 1, 4, 3, 2, 5 };
916 
917 const int Geometry::
918 Constants<Geometry::SQUARE>::Edges[4][2] = {{0, 1}, {1, 2}, {2, 3}, {3, 0}};
919 const int Geometry::
920 Constants<Geometry::SQUARE>::VertToVert::I[4] = {0, 2, 3, 4};
921 const int Geometry::
922 Constants<Geometry::SQUARE>::VertToVert::J[4][2] =
923 {{1, 0}, {3, -4}, {2, 1}, {3, 2}};
924 const int Geometry::
925 Constants<Geometry::SQUARE>::FaceVert[1][4] = {{0, 1, 2, 3}};
926 const int Geometry::
927 Constants<Geometry::SQUARE>::Orient[8][4] =
928 {
929  {0, 1, 2, 3}, {0, 3, 2, 1}, {1, 2, 3, 0}, {1, 0, 3, 2},
930  {2, 3, 0, 1}, {2, 1, 0, 3}, {3, 0, 1, 2}, {3, 2, 1, 0}
931 };
932 const int Geometry::
933 Constants<Geometry::SQUARE>::InvOrient[8] = { 0, 1, 6, 3, 4, 5, 2, 7 };
934 
935 const int Geometry::
936 Constants<Geometry::TETRAHEDRON>::Edges[6][2] =
937 {{0, 1}, {0, 2}, {0, 3}, {1, 2}, {1, 3}, {2, 3}};
938 const int Geometry::
939 Constants<Geometry::TETRAHEDRON>::FaceTypes[4] =
940 {
941  Geometry::TRIANGLE, Geometry::TRIANGLE,
942  Geometry::TRIANGLE, Geometry::TRIANGLE
943 };
944 const int Geometry::
945 Constants<Geometry::TETRAHEDRON>::FaceVert[4][3] =
946 {{1, 2, 3}, {0, 3, 2}, {0, 1, 3}, {0, 2, 1}};
947 const int Geometry::
948 Constants<Geometry::TETRAHEDRON>::VertToVert::I[4] = {0, 3, 5, 6};
949 const int Geometry::
950 Constants<Geometry::TETRAHEDRON>::VertToVert::J[6][2] =
951 {
952  {1, 0}, {2, 1}, {3, 2}, // 0,1:0 0,2:1 0,3:2
953  {2, 3}, {3, 4}, // 1,2:3 1,3:4
954  {3, 5} // 2,3:5
955 };
956 const int Geometry::
957 Constants<Geometry::TETRAHEDRON>::Orient[24][4] =
958 {
959  {0, 1, 2, 3}, {0, 1, 3, 2}, {0, 2, 3, 1}, {0, 2, 1, 3},
960  {0, 3, 1, 2}, {0, 3, 2, 1},
961  {1, 2, 0, 3}, {1, 2, 3, 0}, {1, 3, 2, 0}, {1, 3, 0, 2},
962  {1, 0, 3, 2}, {1, 0, 2, 3},
963  {2, 3, 0, 1}, {2, 3, 1, 0}, {2, 0, 1, 3}, {2, 0, 3, 1},
964  {2, 1, 3, 0}, {2, 1, 0, 3},
965  {3, 0, 2, 1}, {3, 0, 1, 2}, {3, 1, 0, 2}, {3, 1, 2, 0},
966  {3, 2, 1, 0}, {3, 2, 0, 1}
967 };
968 const int Geometry::
969 Constants<Geometry::TETRAHEDRON>::InvOrient[24] =
970 {
971  0, 1, 4, 3, 2, 5,
972  14, 19, 18, 15, 10, 11,
973  12, 23, 6, 9, 20, 17,
974  8, 7, 16, 21, 22, 13
975 };
976 
977 const int Geometry::
978 Constants<Geometry::CUBE>::Edges[12][2] =
979 {
980  {0, 1}, {1, 2}, {3, 2}, {0, 3}, {4, 5}, {5, 6},
981  {7, 6}, {4, 7}, {0, 4}, {1, 5}, {2, 6}, {3, 7}
982 };
983 const int Geometry::
984 Constants<Geometry::CUBE>::FaceTypes[6] =
985 {
986  Geometry::SQUARE, Geometry::SQUARE, Geometry::SQUARE,
987  Geometry::SQUARE, Geometry::SQUARE, Geometry::SQUARE
988 };
989 const int Geometry::
990 Constants<Geometry::CUBE>::FaceVert[6][4] =
991 {
992  {3, 2, 1, 0}, {0, 1, 5, 4}, {1, 2, 6, 5},
993  {2, 3, 7, 6}, {3, 0, 4, 7}, {4, 5, 6, 7}
994 };
995 const int Geometry::
996 Constants<Geometry::CUBE>::VertToVert::I[8] = {0, 3, 5, 7, 8, 10, 11, 12};
997 const int Geometry::
998 Constants<Geometry::CUBE>::VertToVert::J[12][2] =
999 {
1000  {1, 0}, {3, 3}, {4, 8}, // 0,1:0 0,3:3 0,4:8
1001  {2, 1}, {5, 9}, // 1,2:1 1,5:9
1002  {3,-3}, {6,10}, // 2,3:-3 2,6:10
1003  {7,11}, // 3,7:11
1004  {5, 4}, {7, 7}, // 4,5:4 4,7:7
1005  {6, 5}, // 5,6:5
1006  {7,-7} // 6,7:-7
1007 };
1008 
1009 const int Geometry::
1010 Constants<Geometry::PRISM>::Edges[9][2] =
1011 {{0, 1}, {1, 2}, {2, 0}, {3, 4}, {4, 5}, {5, 3}, {0, 3}, {1, 4}, {2, 5}};
1012 const int Geometry::
1013 Constants<Geometry::PRISM>::FaceTypes[5] =
1014 {
1015  Geometry::TRIANGLE, Geometry::TRIANGLE,
1016  Geometry::SQUARE, Geometry::SQUARE, Geometry::SQUARE
1017 };
1018 const int Geometry::
1019 Constants<Geometry::PRISM>::FaceVert[5][4] =
1020 {{0, 2, 1, -1}, {3, 4, 5, -1}, {0, 1, 4, 3}, {1, 2, 5, 4}, {2, 0, 3, 5}};
1021 const int Geometry::
1022 Constants<Geometry::PRISM>::VertToVert::I[6] = {0, 3, 5, 6, 8, 9};
1023 const int Geometry::
1024 Constants<Geometry::PRISM>::VertToVert::J[9][2] =
1025 {
1026  {1, 0}, {2, -3}, {3, 6}, // 0,1:0 0,2:-3 0,3:6
1027  {2, 1}, {4, 7}, // 1,2:1 1,4:7
1028  {5, 8}, // 2,5:8
1029  {4, 3}, {5, -6}, // 3,4:3 3,5:-6
1030  {5, 4} // 4,5:4
1031 };
1032 
1033 const int Geometry::
1034 Constants<Geometry::PYRAMID>::Edges[8][2] =
1035 {{0, 1}, {1, 2}, {3, 2}, {0, 3}, {0, 4}, {1, 4}, {2, 4}, {3, 4}};
1036 const int Geometry::
1037 Constants<Geometry::PYRAMID>::FaceTypes[5] =
1038 {
1039  Geometry::SQUARE,
1040  Geometry::TRIANGLE, Geometry::TRIANGLE,
1041  Geometry::TRIANGLE, Geometry::TRIANGLE
1042 };
1043 const int Geometry::
1044 Constants<Geometry::PYRAMID>::FaceVert[5][4] =
1045 {{3, 2, 1, 0}, {0, 1, 4, -1}, {1, 2, 4, -1}, {2, 3, 4, -1}, {3, 0, 4, -1}};
1046 const int Geometry::
1047 Constants<Geometry::PYRAMID>::VertToVert::I[5] = {0, 3, 5, 7, 8};
1048 const int Geometry::
1049 Constants<Geometry::PYRAMID>::VertToVert::J[8][2] =
1050 {
1051  {1, 0}, {3, 3}, {4, 4}, // 0,1:0 0,3:3 0,4:4
1052  {2, 1}, {4, 5}, // 1,2:1 1,4:5
1053  {3,-3}, {4, 6}, // 2,3:-3 2,4:6
1054  {4, 7} // 3,4:7
1055 };
1056 
1057 
1058 GeometryRefiner::GeometryRefiner()
1059 {
1060  type = Quadrature1D::ClosedUniform;
1061 }
1062 
1063 GeometryRefiner::~GeometryRefiner()
1064 {
1065  for (int i = 0; i < Geometry::NumGeom; i++)
1066  {
1067  for (int j = 0; j < RGeom[i].Size(); j++) { delete RGeom[i][j]; }
1068  for (int j = 0; j < IntPts[i].Size(); j++) { delete IntPts[i][j]; }
1069  }
1070 }
1071 
1072 RefinedGeometry *GeometryRefiner::FindInRGeom(Geometry::Type Geom,
1073  int Times, int ETimes,
1074  int Type)
1075 {
1076  Array<RefinedGeometry *> &RGA = RGeom[Geom];
1077  for (int i = 0; i < RGA.Size(); i++)
1078  {
1079  RefinedGeometry &RG = *RGA[i];
1080  if (RG.Times == Times && RG.ETimes == ETimes && RG.Type == Type)
1081  {
1082  return &RG;
1083  }
1084  }
1085  return NULL;
1086 }
1087 
1088 IntegrationRule *GeometryRefiner::FindInIntPts(Geometry::Type Geom, int NPts)
1089 {
1090  Array<IntegrationRule *> &IPA = IntPts[Geom];
1091  for (int i = 0; i < IPA.Size(); i++)
1092  {
1093  IntegrationRule &ir = *IPA[i];
1094  if (ir.GetNPoints() == NPts) { return &ir; }
1095  }
1096  return NULL;
1097 }
1098 
1099 RefinedGeometry * GeometryRefiner::Refine(Geometry::Type Geom,
1100  int Times, int ETimes)
1101 {
1102  int i, j, k, l, m;
1103 
1104  Times = std::max(Times, 1);
1105  ETimes = Geometry::Dimension[Geom] <= 1 ? 0 : std::max(ETimes, 1);
1106  const double *cp = poly1d.GetPoints(Times, BasisType::GetNodalBasis(type));
1107 
1108  RefinedGeometry *RG = FindInRGeom(Geom, Times, ETimes, type);
1109  if (RG) { return RG; }
1110 
1111  switch (Geom)
1112  {
1113  case Geometry::POINT:
1114  {
1115  RG = new RefinedGeometry(1, 1, 0);
1116  RG->Times = 1;
1117  RG->ETimes = 0;
1118  RG->Type = type;
1119  RG->RefPts.IntPoint(0).x = cp[0];
1120  RG->RefGeoms[0] = 0;
1121 
1122  RGeom[Geometry::POINT].Append(RG);
1123  return RG;
1124  }
1125 
1126  case Geometry::SEGMENT:
1127  {
1128  RG = new RefinedGeometry(Times+1, 2*Times, 0);
1129  RG->Times = Times;
1130  RG->ETimes = 0;
1131  RG->Type = type;
1132  for (i = 0; i <= Times; i++)
1133  {
1134  IntegrationPoint &ip = RG->RefPts.IntPoint(i);
1135  ip.x = cp[i];
1136  }
1137  Array<int> &G = RG->RefGeoms;
1138  for (i = 0; i < Times; i++)
1139  {
1140  G[2*i+0] = i;
1141  G[2*i+1] = i+1;
1142  }
1143 
1144  RGeom[Geometry::SEGMENT].Append(RG);
1145  return RG;
1146  }
1147 
1148  case Geometry::TRIANGLE:
1149  {
1150  RG = new RefinedGeometry((Times+1)*(Times+2)/2, 3*Times*Times,
1151  3*Times*(ETimes+1), 3*Times);
1152  RG->Times = Times;
1153  RG->ETimes = ETimes;
1154  RG->Type = type;
1155  for (k = j = 0; j <= Times; j++)
1156  for (i = 0; i <= Times-j; i++, k++)
1157  {
1158  IntegrationPoint &ip = RG->RefPts.IntPoint(k);
1159  ip.x = cp[i]/(cp[i] + cp[j] + cp[Times-i-j]);
1160  ip.y = cp[j]/(cp[i] + cp[j] + cp[Times-i-j]);
1161  }
1162  Array<int> &G = RG->RefGeoms;
1163  for (l = k = j = 0; j < Times; j++, k++)
1164  for (i = 0; i < Times-j; i++, k++)
1165  {
1166  G[l++] = k;
1167  G[l++] = k+1;
1168  G[l++] = k+Times-j+1;
1169  if (i+j+1 < Times)
1170  {
1171  G[l++] = k+1;
1172  G[l++] = k+Times-j+2;
1173  G[l++] = k+Times-j+1;
1174  }
1175  }
1176  Array<int> &E = RG->RefEdges;
1177  int lb = 0, li = 2*RG->NumBdrEdges;
1178  // horizontal edges
1179  for (k = 0; k < Times; k += Times/ETimes)
1180  {
1181  int &lt = (k == 0) ? lb : li;
1182  j = k*(Times+1)-((k-1)*k)/2;
1183  for (i = 0; i < Times-k; i++)
1184  {
1185  E[lt++] = j; j++;
1186  E[lt++] = j;
1187  }
1188  }
1189  // diagonal edges
1190  for (k = Times; k > 0; k -= Times/ETimes)
1191  {
1192  int &lt = (k == Times) ? lb : li;
1193  j = k;
1194  for (i = 0; i < k; i++)
1195  {
1196  E[lt++] = j; j += Times-i;
1197  E[lt++] = j;
1198  }
1199  }
1200  // vertical edges
1201  for (k = 0; k < Times; k += Times/ETimes)
1202  {
1203  int &lt = (k == 0) ? lb : li;
1204  j = k;
1205  for (i = 0; i < Times-k; i++)
1206  {
1207  E[lt++] = j; j += Times-i+1;
1208  E[lt++] = j;
1209  }
1210  }
1211 
1212  RGeom[Geometry::TRIANGLE].Append(RG);
1213  return RG;
1214  }
1215 
1216  case Geometry::SQUARE:
1217  {
1218  RG = new RefinedGeometry((Times+1)*(Times+1), 4*Times*Times,
1219  4*(ETimes+1)*Times, 4*Times);
1220  RG->Times = Times;
1221  RG->ETimes = ETimes;
1222  RG->Type = type;
1223  for (k = j = 0; j <= Times; j++)
1224  for (i = 0; i <= Times; i++, k++)
1225  {
1226  IntegrationPoint &ip = RG->RefPts.IntPoint(k);
1227  ip.x = cp[i];
1228  ip.y = cp[j];
1229  }
1230  Array<int> &G = RG->RefGeoms;
1231  for (l = k = j = 0; j < Times; j++, k++)
1232  for (i = 0; i < Times; i++, k++)
1233  {
1234  G[l++] = k;
1235  G[l++] = k+1;
1236  G[l++] = k+Times+2;
1237  G[l++] = k+Times+1;
1238  }
1239  Array<int> &E = RG->RefEdges;
1240  int lb = 0, li = 2*RG->NumBdrEdges;
1241  // horizontal edges
1242  for (k = 0; k <= Times; k += Times/ETimes)
1243  {
1244  int &lt = (k == 0 || k == Times) ? lb : li;
1245  for (i = 0, j = k*(Times+1); i < Times; i++)
1246  {
1247  E[lt++] = j; j++;
1248  E[lt++] = j;
1249  }
1250  }
1251  // vertical edges (in right-to-left order)
1252  for (k = Times; k >= 0; k -= Times/ETimes)
1253  {
1254  int &lt = (k == Times || k == 0) ? lb : li;
1255  for (i = 0, j = k; i < Times; i++)
1256  {
1257  E[lt++] = j; j += Times+1;
1258  E[lt++] = j;
1259  }
1260  }
1261 
1262  RGeom[Geometry::SQUARE].Append(RG);
1263  return RG;
1264  }
1265 
1266  case Geometry::CUBE:
1267  {
1268  RG = new RefinedGeometry ((Times+1)*(Times+1)*(Times+1),
1269  8*Times*Times*Times, 0);
1270  RG->Times = Times;
1271  RG->ETimes = ETimes;
1272  RG->Type = type;
1273  for (l = k = 0; k <= Times; k++)
1274  for (j = 0; j <= Times; j++)
1275  for (i = 0; i <= Times; i++, l++)
1276  {
1277  IntegrationPoint &ip = RG->RefPts.IntPoint(l);
1278  ip.x = cp[i];
1279  ip.y = cp[j];
1280  ip.z = cp[k];
1281  }
1282  Array<int> &G = RG->RefGeoms;
1283  for (l = k = 0; k < Times; k++)
1284  for (j = 0; j < Times; j++)
1285  for (i = 0; i < Times; i++)
1286  {
1287  G[l++] = i+0 + (j+0 + (k+0) * (Times+1)) * (Times+1);
1288  G[l++] = i+1 + (j+0 + (k+0) * (Times+1)) * (Times+1);
1289  G[l++] = i+1 + (j+1 + (k+0) * (Times+1)) * (Times+1);
1290  G[l++] = i+0 + (j+1 + (k+0) * (Times+1)) * (Times+1);
1291  G[l++] = i+0 + (j+0 + (k+1) * (Times+1)) * (Times+1);
1292  G[l++] = i+1 + (j+0 + (k+1) * (Times+1)) * (Times+1);
1293  G[l++] = i+1 + (j+1 + (k+1) * (Times+1)) * (Times+1);
1294  G[l++] = i+0 + (j+1 + (k+1) * (Times+1)) * (Times+1);
1295  }
1296 
1297  RGeom[Geometry::CUBE].Append(RG);
1298  return RG;
1299  }
1300 
1301  case Geometry::TETRAHEDRON:
1302  {
1303  // subdivide the tetrahedron with vertices
1304  // (0,0,0), (0,0,1), (1,1,1), (0,1,1)
1305 
1306  // vertices: 0 <= i <= j <= k <= Times
1307  // (3-combination with repetitions)
1308  // number of vertices: (n+3)*(n+2)*(n+1)/6, n = Times
1309 
1310  // elements: the vertices are: v1=(i,j,k), v2=v1+u1, v3=v2+u2, v4=v3+u3
1311  // where 0 <= i <= j <= k <= n-1 and
1312  // u1,u2,u3 is a permutation of (1,0,0),(0,1,0),(0,0,1)
1313  // such that all v2,v3,v4 have non-decreasing components
1314  // number of elements: n^3
1315 
1316  const int n = Times;
1317  RG = new RefinedGeometry((n+3)*(n+2)*(n+1)/6, 4*n*n*n, 0);
1318  RG->Times = Times;
1319  RG->ETimes = ETimes;
1320  RG->Type = type;
1321  // enumerate and define the vertices
1322  Array<int> vi((n+1)*(n+1)*(n+1));
1323  vi = -1;
1324  m = 0;
1325 
1326  // vertices are given in lexicographic ordering on the reference
1327  // element
1328  for (int kk = 0; kk <= n; kk++)
1329  for (int jj = 0; jj <= n-kk; jj++)
1330  for (int ii = 0; ii <= n-jj-kk; ii++)
1331  {
1332  IntegrationPoint &ip = RG->RefPts.IntPoint(m);
1333  double w = cp[ii] + cp[jj] + cp[kk] + cp[Times-ii-jj-kk];
1334  ip.x = cp[ii]/w;
1335  ip.y = cp[jj]/w;
1336  ip.z = cp[kk]/w;
1337  // (ii,jj,kk) are coordinates in the reference tetrahedron,
1338  // transform to coordinates (i,j,k) in the auxiliary
1339  // tetrahedron defined by (0,0,0), (0,0,1), (1,1,1), (0,1,1)
1340  i = jj;
1341  j = jj+kk;
1342  k = ii+jj+kk;
1343  l = i + (j + k * (n+1)) * (n+1);
1344  // map from linear Cartesian hex index in the auxiliary tet
1345  // to lexicographic in the reference tet
1346  vi[l] = m;
1347  m++;
1348  }
1349 
1350  if (m != (n+3)*(n+2)*(n+1)/6)
1351  {
1352  mfem_error("GeometryRefiner::Refine() for TETRAHEDRON #1");
1353  }
1354  // elements
1355  Array<int> &G = RG->RefGeoms;
1356  m = 0;
1357  for (k = 0; k < n; k++)
1358  for (j = 0; j <= k; j++)
1359  for (i = 0; i <= j; i++)
1360  {
1361  // the ordering of the vertices is chosen to ensure:
1362  // 1) correct orientation
1363  // 2) the x,y,z edges are in the set of edges
1364  // {(0,1),(2,3), (0,2),(1,3)}
1365  // (goal is to ensure that subsequent refinement using
1366  // this procedure preserves the six tetrahedral shapes)
1367 
1368  // zyx: (i,j,k)-(i,j,k+1)-(i+1,j+1,k+1)-(i,j+1,k+1)
1369  G[m++] = vi[i+0 + (j+0 + (k+0) * (n+1)) * (n+1)];
1370  G[m++] = vi[i+0 + (j+0 + (k+1) * (n+1)) * (n+1)];
1371  G[m++] = vi[i+1 + (j+1 + (k+1) * (n+1)) * (n+1)];
1372  G[m++] = vi[i+0 + (j+1 + (k+1) * (n+1)) * (n+1)];
1373  if (j < k)
1374  {
1375  // yzx: (i,j,k)-(i+1,j+1,k+1)-(i,j+1,k)-(i,j+1,k+1)
1376  G[m++] = vi[i+0 + (j+0 + (k+0) * (n+1)) * (n+1)];
1377  G[m++] = vi[i+1 + (j+1 + (k+1) * (n+1)) * (n+1)];
1378  G[m++] = vi[i+0 + (j+1 + (k+0) * (n+1)) * (n+1)];
1379  G[m++] = vi[i+0 + (j+1 + (k+1) * (n+1)) * (n+1)];
1380  // yxz: (i,j,k)-(i,j+1,k)-(i+1,j+1,k+1)-(i+1,j+1,k)
1381  G[m++] = vi[i+0 + (j+0 + (k+0) * (n+1)) * (n+1)];
1382  G[m++] = vi[i+0 + (j+1 + (k+0) * (n+1)) * (n+1)];
1383  G[m++] = vi[i+1 + (j+1 + (k+1) * (n+1)) * (n+1)];
1384  G[m++] = vi[i+1 + (j+1 + (k+0) * (n+1)) * (n+1)];
1385  }
1386  if (i < j)
1387  {
1388  // xzy: (i,j,k)-(i+1,j,k)-(i+1,j+1,k+1)-(i+1,j,k+1)
1389  G[m++] = vi[i+0 + (j+0 + (k+0) * (n+1)) * (n+1)];
1390  G[m++] = vi[i+1 + (j+0 + (k+0) * (n+1)) * (n+1)];
1391  G[m++] = vi[i+1 + (j+1 + (k+1) * (n+1)) * (n+1)];
1392  G[m++] = vi[i+1 + (j+0 + (k+1) * (n+1)) * (n+1)];
1393  if (j < k)
1394  {
1395  // xyz: (i,j,k)-(i+1,j+1,k+1)-(i+1,j,k)-(i+1,j+1,k)
1396  G[m++] = vi[i+0 + (j+0 + (k+0) * (n+1)) * (n+1)];
1397  G[m++] = vi[i+1 + (j+1 + (k+1) * (n+1)) * (n+1)];
1398  G[m++] = vi[i+1 + (j+0 + (k+0) * (n+1)) * (n+1)];
1399  G[m++] = vi[i+1 + (j+1 + (k+0) * (n+1)) * (n+1)];
1400  }
1401  // zxy: (i,j,k)-(i+1,j+1,k+1)-(i,j,k+1)-(i+1,j,k+1)
1402  G[m++] = vi[i+0 + (j+0 + (k+0) * (n+1)) * (n+1)];
1403  G[m++] = vi[i+1 + (j+1 + (k+1) * (n+1)) * (n+1)];
1404  G[m++] = vi[i+0 + (j+0 + (k+1) * (n+1)) * (n+1)];
1405  G[m++] = vi[i+1 + (j+0 + (k+1) * (n+1)) * (n+1)];
1406  }
1407  }
1408  if (m != 4*n*n*n)
1409  {
1410  mfem_error("GeometryRefiner::Refine() for TETRAHEDRON #2");
1411  }
1412  for (i = 0; i < m; i++)
1413  if (G[i] < 0)
1414  {
1415  mfem_error("GeometryRefiner::Refine() for TETRAHEDRON #3");
1416  }
1417 
1418  RGeom[Geometry::TETRAHEDRON].Append(RG);
1419  return RG;
1420  }
1421 
1422  case Geometry::PYRAMID:
1423  {
1424  const int n = Times;
1425  RG = new RefinedGeometry ((n+1)*(n+2)*(2*n+3)/6,
1426  5*n*(2*n-1)*(2*n+1)/3, 0);
1427  RG->Times = Times;
1428  RG->ETimes = ETimes;
1429  RG->Type = type;
1430  // enumerate and define the vertices
1431  m = 0;
1432  for (k = 0; k <= n; k++)
1433  {
1434  const double *cpij =
1435  poly1d.GetPoints(Times - k, BasisType::GetNodalBasis(type));
1436  for (j = 0; j <= n - k; j++)
1437  for (i = 0; i <= n - k; i++)
1438  {
1439  IntegrationPoint &ip = RG->RefPts.IntPoint(m);
1440  if (type == 0)
1441  {
1442  ip.x = (n > k) ? (double(i) / (n - k)) : 0.0;
1443  ip.y = (n > k) ? (double(j) / (n - k)) : 0.0;
1444  ip.z = double(k) / n;
1445  }
1446  else
1447  {
1448  ip.x = cpij[i] * (1.0 - cp[k]);
1449  ip.y = cpij[j] * (1.0 - cp[k]);
1450  ip.z = cp[k];
1451  }
1452  m++;
1453  }
1454  }
1455  if (m != (n+1)*(n+2)*(2*n+3)/6)
1456  {
1457  mfem_error("GeometryRefiner::Refine() for PYRAMID #1");
1458  }
1459  // elements
1460  Array<int> &G = RG->RefGeoms;
1461  m = 0;
1462  for (k = 0; k < n; k++)
1463  {
1464  int lk = k * (k * (2 * k - 6 * n - 9) + 6 * n * (n + 3) + 13) / 6;
1465  int lkp1 = (k + 1) *
1466  (k * (2 * k - 6 * n -5) + 6 * n * (n + 2) + 6) / 6;
1467  for (j = 0; j < n - k; j++)
1468  {
1469  for (i = 0; i < n - k; i++)
1470  {
1471  G[m++] = lk + j * (n - k + 1) + i;
1472  G[m++] = lk + j * (n - k + 1) + i + 1;
1473  G[m++] = lk + (j + 1) * (n - k + 1) + i + 1;
1474  G[m++] = lk + (j + 1) * (n - k + 1) + i;
1475  G[m++] = lkp1 + j * (n - k) + i;
1476  }
1477  }
1478  for (j = 0; j < n - k - 1; j++)
1479  {
1480  for (i = 0; i < n - k - 1; i++)
1481  {
1482  G[m++] = lkp1 + j * (n - k) + i;
1483  G[m++] = lkp1 + (j + 1) * (n - k) + i;
1484  G[m++] = lkp1 + (j + 1) * (n - k) + i + 1;
1485  G[m++] = lkp1 + j * (n - k) + i + 1;
1486  G[m++] = lk + (j + 1) * (n - k + 1) + i + 1;
1487  }
1488  }
1489  for (j = 0; j < n - k; j++)
1490  {
1491  for (i = 0; i < n - k - 1; i++)
1492  {
1493  G[m++] = lk + j * (n - k + 1) + i + 1;
1494  G[m++] = lk + (j + 1) * (n - k + 1) + i + 1;
1495  G[m++] = lkp1 + j * (n - k) + i;
1496  G[m++] = lkp1 + j * (n - k) + i + 1;
1497  G[m++] = -1;
1498  }
1499  }
1500  for (j = 0; j < n - k - 1; j++)
1501  {
1502  for (i = 0; i < n - k; i++)
1503  {
1504  G[m++] = lk + (j + 1) * (n - k + 1) + i;
1505  G[m++] = lk + (j + 1) * (n - k + 1) + i + 1;
1506  G[m++] = lkp1 + (j + 1) * (n - k) + i;
1507  G[m++] = lkp1 + j * (n - k) + i;
1508  G[m++] = -1;
1509  }
1510  }
1511  }
1512  if (m != 5*n*(2*n-1)*(2*n+1)/3)
1513  {
1514  mfem_error("GeometryRefiner::Refine() for PYRAMID #2");
1515  }
1516  RGeom[Geometry::PYRAMID].Append(RG);
1517  return RG;
1518  }
1519 
1520  case Geometry::PRISM:
1521  {
1522  const int n = Times;
1523  RG = new RefinedGeometry ((n+1)*(n+1)*(n+2)/2, 6*n*n*n, 0);
1524  RG->Times = Times;
1525  RG->ETimes = ETimes;
1526  RG->Type = type;
1527  // enumerate and define the vertices
1528  m = 0;
1529  for (l = k = 0; k <= n; k++)
1530  for (j = 0; j <= n; j++)
1531  for (i = 0; i <= n-j; i++, l++)
1532  {
1533  IntegrationPoint &ip = RG->RefPts.IntPoint(l);
1534  ip.x = cp[i]/(cp[i] + cp[j] + cp[n-i-j]);
1535  ip.y = cp[j]/(cp[i] + cp[j] + cp[n-i-j]);
1536  ip.z = cp[k];
1537  m++;
1538  }
1539  if (m != (n+1)*(n+1)*(n+2)/2)
1540  {
1541  mfem_error("GeometryRefiner::Refine() for PRISM #1");
1542  }
1543  // elements
1544  Array<int> &G = RG->RefGeoms;
1545  m = 0;
1546  for (m = k = 0; k < n; k++)
1547  for (l = j = 0; j < n; j++, l++)
1548  for (i = 0; i < n-j; i++, l++)
1549  {
1550  G[m++] = l + (k+0) * (n+1) * (n+2) / 2;
1551  G[m++] = l + 1 + (k+0) * (n+1) * (n+2) / 2;
1552  G[m++] = l - j + (2 + (k+0) * (n+2)) * (n+1) / 2;
1553  G[m++] = l + (k+1) * (n+1) * (n+2) / 2;
1554  G[m++] = l + 1 + (k+1) * (n+1) * (Times+2) / 2;
1555  G[m++] = l - j + (2 + (k+1) * (n+2)) * (n+1) / 2;
1556  if (i+j+1 < n)
1557  {
1558  G[m++] = l + 1 + (k+0) * (n+1) * (n+2)/2;
1559  G[m++] = l - j + (2 + (k+0) * (n+1)) * (n+2) / 2;
1560  G[m++] = l - j + (2 + (k+0) * (n+2)) * (n+1) / 2;
1561  G[m++] = l + 1 + (k+1) * (n+1) * (n+2) / 2;
1562  G[m++] = l - j + (2 + (k+1) * (n+1)) * (n+2) / 2;
1563  G[m++] = l - j + (2 + (k+1) * (n+2)) * (n+1) / 2;
1564  }
1565  }
1566  if (m != 6*n*n*n)
1567  {
1568  mfem_error("GeometryRefiner::Refine() for PRISM #2");
1569  }
1570  for (i = 0; i < m; i++)
1571  if (G[i] < 0)
1572  {
1573  mfem_error("GeometryRefiner::Refine() for PRISM #3");
1574  }
1575 
1576  RGeom[Geometry::PRISM].Append(RG);
1577  return RG;
1578  }
1579 
1580  default:
1581 
1582  return NULL;
1583  }
1584 }
1585 
1586 const IntegrationRule *GeometryRefiner::RefineInterior(Geometry::Type Geom,
1587  int Times)
1588 {
1589  IntegrationRule *ir = NULL;
1590 
1591  switch (Geom)
1592  {
1593  case Geometry::SEGMENT:
1594  {
1595  if (Times < 2)
1596  {
1597  return NULL;
1598  }
1599  ir = FindInIntPts(Geom, Times-1);
1600  if (ir) { return ir; }
1601 
1602  ir = new IntegrationRule(Times-1);
1603  for (int i = 1; i < Times; i++)
1604  {
1605  IntegrationPoint &ip = ir->IntPoint(i-1);
1606  ip.x = double(i) / Times;
1607  ip.y = ip.z = 0.0;
1608  }
1609  }
1610  break;
1611 
1612  case Geometry::TRIANGLE:
1613  {
1614  if (Times < 3)
1615  {
1616  return NULL;
1617  }
1618  ir = FindInIntPts(Geom, ((Times-1)*(Times-2))/2);
1619  if (ir) { return ir; }
1620 
1621  ir = new IntegrationRule(((Times-1)*(Times-2))/2);
1622  for (int k = 0, j = 1; j < Times-1; j++)
1623  for (int i = 1; i < Times-j; i++, k++)
1624  {
1625  IntegrationPoint &ip = ir->IntPoint(k);
1626  ip.x = double(i) / Times;
1627  ip.y = double(j) / Times;
1628  ip.z = 0.0;
1629  }
1630  }
1631  break;
1632 
1633  case Geometry::SQUARE:
1634  {
1635  if (Times < 2)
1636  {
1637  return NULL;
1638  }
1639  ir = FindInIntPts(Geom, (Times-1)*(Times-1));
1640  if (ir) { return ir; }
1641 
1642  ir = new IntegrationRule((Times-1)*(Times-1));
1643  for (int k = 0, j = 1; j < Times; j++)
1644  for (int i = 1; i < Times; i++, k++)
1645  {
1646  IntegrationPoint &ip = ir->IntPoint(k);
1647  ip.x = double(i) / Times;
1648  ip.y = double(j) / Times;
1649  ip.z = 0.0;
1650  }
1651  }
1652  break;
1653 
1654  default:
1655  mfem_error("GeometryRefiner::RefineInterior(...)");
1656  }
1657 
1658  MFEM_ASSERT(ir != NULL, "Failed to construct the refined IntegrationRule.");
1659  IntPts[Geom].Append(ir);
1660 
1661  return ir;
1662 }
1663 
1664 
1665 int GeometryRefiner::GetRefinementLevelFromPoints(Geometry::Type geom, int Npts)
1666 {
1667  switch (geom)
1668  {
1669  case Geometry::POINT:
1670  {
1671  return -1;
1672  }
1673  case Geometry::SEGMENT:
1674  {
1675  return Npts -1;
1676  }
1677  case Geometry::TRIANGLE:
1678  {
1679  for (int n = 0, np = 0; (n < 15) && (np < Npts) ; n++)
1680  {
1681  np = (n+1)*(n+2)/2;
1682  if (np == Npts) { return n; }
1683  }
1684  return -1;
1685  }
1686  case Geometry::SQUARE:
1687  {
1688  for (int n = 0, np = 0; (n < 15) && (np < Npts) ; n++)
1689  {
1690  np = (n+1)*(n+1);
1691  if (np == Npts) { return n; }
1692  }
1693  return -1;
1694  }
1695  case Geometry::CUBE:
1696  {
1697  for (int n = 0, np = 0; (n < 15) && (np < Npts) ; n++)
1698  {
1699  np = (n+1)*(n+1)*(n+1);
1700  if (np == Npts) { return n; }
1701  }
1702  return -1;
1703  }
1704  case Geometry::TETRAHEDRON:
1705  {
1706  for (int n = 0, np = 0; (n < 15) && (np < Npts) ; n++)
1707  {
1708  np = (n+3)*(n+2)*(n+1)/6;
1709  if (np == Npts) { return n; }
1710  }
1711  return -1;
1712  }
1713  case Geometry::PRISM:
1714  {
1715  for (int n = 0, np = 0; (n < 15) && (np < Npts) ; n++)
1716  {
1717  np = (n+1)*(n+1)*(n+2)/2;
1718  if (np == Npts) { return n; }
1719  }
1720  return -1;
1721  }
1722  default:
1723  {
1724  mfem_error("Non existing Geometry.");
1725  }
1726  }
1727 
1728  return -1;
1729 }
1730 
1731 
1732 int GeometryRefiner::GetRefinementLevelFromElems(Geometry::Type geom, int Nels)
1733 {
1734  switch (geom)
1735  {
1736  case Geometry::POINT:
1737  {
1738  return -1;
1739  }
1740  case Geometry::SEGMENT:
1741  {
1742  return Nels;
1743  }
1744  case Geometry::TRIANGLE:
1745  case Geometry::SQUARE:
1746  {
1747  for (int n = 0; (n < 15) && (n*n < Nels+1) ; n++)
1748  {
1749  if (n*n == Nels) { return n-1; }
1750  }
1751  return -1;
1752  }
1753  case Geometry::CUBE:
1754  case Geometry::TETRAHEDRON:
1755  case Geometry::PRISM:
1756  {
1757  for (int n = 0; (n < 15) && (n*n*n < Nels+1) ; n++)
1758  {
1759  if (n*n*n == Nels) { return n-1; }
1760  }
1761  return -1;
1762  }
1763  default:
1764  {
1765  mfem_error("Non existing Geometry.");
1766  }
1767  }
1768 
1769  return -1;
1770 }
1771 
1772 
1773 GeometryRefiner GlobGeometryRefiner;
1774 
1775 }
mfem::Array::Size
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
mfem::Geometry::Constants< Geometry::POINT >::InvOrient
static const int InvOrient[NumOrient]
Definition: geom.hpp:136
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_base.hpp:77
mfem::PyramidFE
class LinearPyramidFiniteElement PyramidFE
Definition: fe.cpp:44
mfem::Geometry::NumGeom
static const int NumGeom
Definition: geom.hpp:42
mfem::Geometry::JacToPerfJac
void JacToPerfJac(int GeomType, const DenseMatrix &J, DenseMatrix &PJ) const
Definition: geom.cpp:867
mfem::Geometry::GetRandomPoint
static void GetRandomPoint(int GeomType, IntegrationPoint &ip)
Get a random point in the reference element specified by GeomType.
Definition: geom.cpp:285
mfem::Mult
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:472
mfem::Geometry::Constants< Geometry::TETRAHEDRON >::FaceVert
static const int FaceVert[NumFaces][MaxFaceVert]
Definition: geom.hpp:206
mfem::ElementTransformation::SetIntPoint
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition: eltrans.hpp:93
mfem::Quadrature1D::ClosedUniform
aka closed Newton-Cotes
Definition: intrules.hpp:298
mfem::Geometry::Constants< Geometry::PRISM >::FaceVert
static const int FaceVert[NumFaces][MaxFaceVert]
Definition: geom.hpp:246
mfem::Geometry::Constants< Geometry::PRISM >::Edges
static const int Edges[NumEdges][2]
Definition: geom.hpp:242
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:46
mfem::Geometry::NumEdges
static const int NumEdges[NumGeom]
Definition: geom.hpp:50
mfem::RefinedGeometry::RefEdges
Array< int > RefEdges
Definition: geom.hpp:283
mfem::Geometry::GetVertices
const IntegrationRule * GetVertices(int GeomType)
Return an IntegrationRule consisting of all vertices of the given Geometry::Type, GeomType...
Definition: geom.cpp:265
mfem::GeometryRefiner::RefineInterior
const IntegrationRule * RefineInterior(Geometry::Type Geom, int Times)
Definition: geom.cpp:1586
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_base.cpp:2131
mfem::Geometry::NumFaces
static const int NumFaces[NumGeom]
Definition: geom.hpp:51
mfem::Geometry::Constants< Geometry::TRIANGLE >::InvOrient
static const int InvOrient[NumOrient]
Definition: geom.hpp:174
mfem::IntegrationRule::IntPoint
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:250
mfem::Geometry::Constants< Geometry::TRIANGLE >::Edges
static const int Edges[NumEdges][2]
Definition: geom.hpp:156
mfem::Geometry::Dimension
static const int Dimension[NumGeom]
Definition: geom.hpp:47
mfem::Geometry::Constants< Geometry::CUBE >::FaceTypes
static const int FaceTypes[NumFaces]
Definition: geom.hpp:226
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:154
mfem::Geometry::Constants< Geometry::PRISM >::FaceTypes
static const int FaceTypes[NumFaces]
Definition: geom.hpp:244
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:119
mfem::Geometry::NumVerts
static const int NumVerts[NumGeom]
Definition: geom.hpp:49
mfem::Geometry::Constants< Geometry::PYRAMID >::FaceTypes
static const int FaceTypes[NumFaces]
Definition: geom.hpp:262
mfem::Geometry::Constants< Geometry::SEGMENT >::Edges
static const int Edges[NumEdges][2]
Definition: geom.hpp:144
mfem::GlobGeometryRefiner
GeometryRefiner GlobGeometryRefiner
Definition: geom.cpp:1773
mfem::RefinedGeometry::RefPts
IntegrationRule RefPts
Definition: geom.hpp:282
mfem::Geometry::NumBdrArray
static const int NumBdrArray[NumGeom]
Definition: geom.hpp:44
mfem::GeometryRefiner::Refine
RefinedGeometry * Refine(Geometry::Type Geom, int Times, int ETimes=1)
Definition: geom.cpp:1099
mfem::Geometry::Constants< Geometry::SQUARE >::InvOrient
static const int InvOrient[NumOrient]
Definition: geom.hpp:194
mfem::Geometry::Constants< Geometry::SEGMENT >::Orient
static const int Orient[NumOrient][NumVert]
Definition: geom.hpp:147
mfem::GeometryRefiner::GetRefinementLevelFromElems
virtual int GetRefinementLevelFromElems(Geometry::Type geom, int Npts)
Get the Refinement level based on number of elements.
Definition: geom.cpp:1732
mfem::Geometry::Constants< Geometry::TETRAHEDRON >::Edges
static const int Edges[NumEdges][2]
Definition: geom.hpp:202
mfem::Geometry::Name
static const char * Name[NumGeom]
Definition: geom.hpp:45
mfem::Geometry::Constants< Geometry::TETRAHEDRON >::InvOrient
static const int InvOrient[NumOrient]
Definition: geom.hpp:216
mfem::CalcInverse
void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
Definition: densemat.cpp:2182
mfem::Geometry::Constants< Geometry::SEGMENT >::InvOrient
static const int InvOrient[NumOrient]
Definition: geom.hpp:148
mfem::Geometry::Constants< Geometry::SQUARE >::Edges
static const int Edges[NumEdges][2]
Definition: geom.hpp:182
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:583
mfem::Geometry::GetPerfPointMat
void GetPerfPointMat(int GeomType, DenseMatrix &pm)
Definition: geom.cpp:783
mfem::TetrahedronFE
class Linear3DFiniteElement TetrahedronFE
Definition: fe.cpp:36
mfem::TriangleFE
Linear2DFiniteElement TriangleFE
Definition: fe.cpp:32
mfem::Geometry::Constants< Geometry::CUBE >::FaceVert
static const int FaceVert[NumFaces][MaxFaceVert]
Definition: geom.hpp:228
mfem::Geometry::Constants< Geometry::TRIANGLE >::FaceVert
static const int FaceVert[NumFaces][NumVert]
Definition: geom.hpp:164
mfem::Geometry::Constants< Geometry::PYRAMID >::Edges
static const int Edges[NumEdges][2]
Definition: geom.hpp:260
mfem::IntegrationPoint
Class for integration point with weight.
Definition: intrules.hpp:25
mfem::IsoparametricTransformation
A standard isoparametric element transformation.
Definition: eltrans.hpp:361
mfem::GeometryRefiner::GetRefinementLevelFromPoints
virtual int GetRefinementLevelFromPoints(Geometry::Type Geom, int Npts)
Get the Refinement level based on number of points.
Definition: geom.cpp:1665
mfem::Geometry::DimStart
static const int DimStart[MaxDim+2]
Definition: geom.hpp:48
mfem::DenseMatrix::Diag
void Diag(double c, int n)
Creates n x n diagonal matrix with diagonal elements c.
Definition: densemat.cpp:1329
mfem::WedgeFE
class LinearWedgeFiniteElement WedgeFE
Definition: fe.cpp:40
mfem::Geometry::Constants< Geometry::SQUARE >::FaceVert
static const int FaceVert[NumFaces][NumVert]
Definition: geom.hpp:190
mfem::Geometry::Constants< Geometry::SQUARE >::Orient
static const int Orient[NumOrient][NumVert]
Definition: geom.hpp:193
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:381
mfem::Geometry::Constants< Geometry::CUBE >::Edges
static const int Edges[NumEdges][2]
Definition: geom.hpp:224
mfem::Geometry::Constants< Geometry::TETRAHEDRON >::Orient
static const int Orient[NumOrient][NumVert]
Definition: geom.hpp:215
t
RefCoord t[3]
Definition: ncmesh_tables.hpp:158
mfem::Geometry::Constants< Geometry::POINT >::Orient
static const int Orient[NumOrient][NumVert]
Definition: geom.hpp:135
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:405
mfem::DenseMatrix::SetSize
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:105
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:404
mfem::RefinedGeometry::RefGeoms
Array< int > RefGeoms
Definition: geom.hpp:283
mfem::Geometry::Constants< Geometry::TRIANGLE >::Orient
static const int Orient[NumOrient][NumVert]
Definition: geom.hpp:172
mfem::poly1d
Poly_1D poly1d
Definition: fe.cpp:28
mfem::Geometry::Constants< Geometry::PYRAMID >::FaceVert
static const int FaceVert[NumFaces][MaxFaceVert]
Definition: geom.hpp:264
Generated on Mon Mar 21 2022 16:49:49 for MFEM by   doxygen 1.8.5