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