Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 :
20 : // Local includes
21 : #include "libmesh/node.h"
22 : #include "libmesh/elem.h"
23 : #include "libmesh/reference_elem.h"
24 : #include "libmesh/libmesh_singleton.h"
25 : #include "libmesh/threads.h"
26 : #include "libmesh/enum_to_string.h"
27 : #include "libmesh/enum_elem_type.h"
28 : #include "libmesh/fuzzy_equals.h"
29 :
30 : // C++ includes
31 : #include <map>
32 : #include <sstream>
33 : #include <memory> // std::unique_ptr
34 : #include <cmath> // std::sqrt, std::cbrt
35 : #include <utility> // std::pair
36 : #include <vector>
37 :
38 :
39 : //-----------------------------------------------
40 : // anonymous namespace for implementation details
41 : namespace
42 : {
43 : using namespace libMesh;
44 :
45 : namespace ElemDataStrings
46 : {
47 : // GCC 5.2.0 warns about overlength strings in the auto-generated
48 : // reference_elem.data file.
49 : #pragma GCC diagnostic ignored "-Woverlength-strings"
50 : #include "reference_elem.data"
51 : #pragma GCC diagnostic warning "-Woverlength-strings"
52 : }
53 :
54 : typedef Threads::spin_mutex InitMutex;
55 :
56 : // Mutex for thread safety.
57 : InitMutex init_mtx;
58 :
59 : // map from ElemType to reference element file system object name
60 : typedef std::map<ElemType, const char *> FileMapType;
61 : FileMapType ref_elem_file;
62 : Elem * ref_elem_map[INVALID_ELEM];
63 :
64 :
65 :
66 14 : class SingletonCache : public libMesh::Singleton
67 : {
68 : public:
69 712 : virtual ~SingletonCache() = default;
70 :
71 : std::vector<std::unique_ptr<Node>> node_list;
72 : std::vector<std::unique_ptr<Elem>> elem_list;
73 : };
74 :
75 : // From [0], regarding the lifetime of the singleton_cache variable:
76 : //
77 : // "All variables at namespace level, including the anonymous
78 : // namespace and function local static variable have static
79 : // storage duration unless they are declared thread_local."
80 : //
81 : // Variables with static storage duration are destroyed at the end of
82 : // the program execution. From [1],
83 : //
84 : // "If it is a pointer to the data which is static ... then like all
85 : // other dynamically allocated data, it will only be destructed when
86 : // you delete it. There are two frequent solutions:
87 : // * use a smart pointer, which has a destructor which deletes it, or
88 : // * don't delete it; in most cases, there's really no reason to call the
89 : // destructor, and if you happen to use the instance in the destructors
90 : // of other static objects, you'll run into an order of destruction
91 : // problem."
92 : //
93 : // I tried making the singleton_cache a std::unique_ptr, but this resulted
94 : // in a segfault during program shutdown which appeared to come from the
95 : // single_cache unique_ptr's destructor. I didn't investigate whether the
96 : // issue was caused by a double deletion or what, but it appears that the
97 : // first suggestion above may not be valid in general.
98 : //
99 : // [0]: https://stackoverflow.com/questions/24342393/how-anonymous-namespaces-avoids-making-global-static-variable
100 : // [1]: https://stackoverflow.com/questions/6850009/c-deleting-static-data
101 : SingletonCache * singleton_cache = nullptr;
102 :
103 :
104 :
105 9075 : void read_ref_elem (const ElemType type_in,
106 : std::istream & in)
107 : {
108 350 : libmesh_assert (singleton_cache != nullptr);
109 :
110 350 : std::string dummy;
111 : unsigned int n_elem, n_nodes, elem_type_read, nn;
112 : double x, y, z;
113 :
114 9075 : in >> dummy;
115 9075 : in >> n_elem; /**/ std::getline (in, dummy); libmesh_assert_equal_to (n_elem, 1);
116 9075 : in >> n_nodes; /**/ std::getline (in, dummy);
117 9075 : in >> dummy; /**/ std::getline (in, dummy);
118 9075 : in >> dummy; /**/ std::getline (in, dummy);
119 9075 : in >> dummy; /**/ std::getline (in, dummy);
120 9075 : in >> dummy; /**/ std::getline (in, dummy);
121 9075 : in >> n_elem; /**/ std::getline (in, dummy); libmesh_assert_equal_to (n_elem, 1);
122 :
123 350 : in >> elem_type_read;
124 :
125 350 : libmesh_assert_less (elem_type_read, INVALID_ELEM);
126 350 : libmesh_assert_equal_to (elem_type_read, static_cast<unsigned int>(type_in));
127 350 : libmesh_assert_equal_to (n_nodes, Elem::type_to_n_nodes_map[elem_type_read]);
128 :
129 : // Construct elem of appropriate type, store in the elem_list
130 9425 : auto & uelem = singleton_cache->elem_list.emplace_back(Elem::build(type_in));
131 :
132 : // We are expecting an identity map, so assert it!
133 103455 : for (unsigned int n=0; n<n_nodes; n++)
134 : {
135 3640 : in >> nn;
136 3640 : libmesh_assert_equal_to (n,nn);
137 : }
138 :
139 103455 : for (unsigned int n=0; n<n_nodes; n++)
140 : {
141 3640 : in >> x >> y >> z;
142 :
143 : auto & new_node =
144 98020 : singleton_cache->node_list.emplace_back(Node::build(x,y,z,n));
145 :
146 98020 : uelem->set_node(n, new_node.get());
147 : }
148 :
149 : // it is entirely possible we ran out of file or encountered
150 : // another error. If so, throw an error.
151 9425 : libmesh_error_msg_if(!in, "ERROR while creating element singleton!");
152 :
153 : // Also store a pointer to the newly created Elem in the ref_elem_map array.
154 9075 : ref_elem_map[type_in] = uelem.get();
155 9075 : }
156 :
157 :
158 :
159 425869 : void init_ref_elem_table()
160 : {
161 : // outside mutex - if this pointer is set, we can trust it.
162 425869 : if (singleton_cache != nullptr)
163 34580 : return;
164 :
165 : // playing with fire here - lock before touching shared
166 : // data structures
167 14 : InitMutex::scoped_lock lock(init_mtx);
168 :
169 : // inside mutex - pointer may have changed while waiting
170 : // for the lock to acquire, check it again.
171 363 : if (singleton_cache != nullptr)
172 0 : return;
173 :
174 : // OK, if we get here we have the lock and we are not
175 : // initialized. populate singleton. Note that we do not
176 : // use a smart pointer to manage the singleton_cache variable
177 : // since variables with static storage duration are destroyed
178 : // automatically at the end of program execution.
179 377 : singleton_cache = new SingletonCache;
180 :
181 : // initialize the reference file table
182 : {
183 14 : ref_elem_file.clear();
184 :
185 : // 0D elements
186 363 : ref_elem_file[NODEELEM] = ElemDataStrings::one_nodeelem;
187 :
188 : // 1D elements
189 363 : ref_elem_file[EDGE2] = ElemDataStrings::one_edge;
190 363 : ref_elem_file[EDGE3] = ElemDataStrings::one_edge3;
191 363 : ref_elem_file[EDGE4] = ElemDataStrings::one_edge4;
192 :
193 : // 2D elements
194 363 : ref_elem_file[TRI3] = ElemDataStrings::one_tri;
195 363 : ref_elem_file[TRI6] = ElemDataStrings::one_tri6;
196 363 : ref_elem_file[TRI7] = ElemDataStrings::one_tri7;
197 :
198 363 : ref_elem_file[QUAD4] = ElemDataStrings::one_quad;
199 363 : ref_elem_file[QUAD8] = ElemDataStrings::one_quad8;
200 363 : ref_elem_file[QUAD9] = ElemDataStrings::one_quad9;
201 :
202 : // 3D elements
203 363 : ref_elem_file[HEX8] = ElemDataStrings::one_hex;
204 363 : ref_elem_file[HEX20] = ElemDataStrings::one_hex20;
205 363 : ref_elem_file[HEX27] = ElemDataStrings::one_hex27;
206 :
207 363 : ref_elem_file[TET4] = ElemDataStrings::one_tet;
208 363 : ref_elem_file[TET10] = ElemDataStrings::one_tet10;
209 363 : ref_elem_file[TET14] = ElemDataStrings::one_tet14;
210 :
211 363 : ref_elem_file[PRISM6] = ElemDataStrings::one_prism;
212 363 : ref_elem_file[PRISM15] = ElemDataStrings::one_prism15;
213 363 : ref_elem_file[PRISM18] = ElemDataStrings::one_prism18;
214 363 : ref_elem_file[PRISM20] = ElemDataStrings::one_prism20;
215 363 : ref_elem_file[PRISM21] = ElemDataStrings::one_prism21;
216 :
217 363 : ref_elem_file[PYRAMID5] = ElemDataStrings::one_pyramid;
218 363 : ref_elem_file[PYRAMID13] = ElemDataStrings::one_pyramid13;
219 363 : ref_elem_file[PYRAMID14] = ElemDataStrings::one_pyramid14;
220 363 : ref_elem_file[PYRAMID18] = ElemDataStrings::one_pyramid18;
221 :
222 : // No entry here for C0POLYGON or C0POLYHEDRON - the reference
223 : // element there depends on the number of nodes and can't be
224 : // precomputed.
225 : }
226 :
227 : // Read'em
228 9438 : for (const auto & [elem_type, filename] : ref_elem_file)
229 : {
230 10125 : std::istringstream stream(filename);
231 9075 : read_ref_elem(elem_type, stream);
232 8375 : }
233 : }
234 :
235 :
236 : // no reason to do this at startup -
237 : // data structures will get initialized *if*
238 : // ReferenceElem::get() is ever called.
239 : // // Class to setup singleton data
240 : // class ReferenceElemSetup : public Singleton::Setup
241 : // {
242 : // void setup ()
243 : // {
244 : // init_ref_elem_table();
245 : // }
246 : // } reference_elem_setup;
247 :
248 : } // anonymous namespace
249 :
250 :
251 :
252 : //----------------------------------------------------------------------------
253 : // external API Implementation
254 : namespace libMesh
255 : {
256 : namespace ReferenceElem
257 : {
258 425869 : const Elem & get (const ElemType type_in)
259 : {
260 34594 : ElemType base_type = type_in;
261 :
262 : // For shell elements, use non shell type as the base type
263 425869 : if (type_in == TRISHELL3)
264 88 : base_type = TRI3;
265 :
266 424989 : if (type_in == QUADSHELL4)
267 3428 : base_type = QUAD4;
268 :
269 391589 : if (type_in == QUADSHELL8)
270 6756 : base_type = QUAD8;
271 :
272 358309 : if (type_in == QUADSHELL9)
273 100 : base_type = QUAD9;
274 :
275 425869 : init_ref_elem_table();
276 :
277 : // Throw an error if the user asked for an ElemType that we don't
278 : // have a reference element for.
279 425869 : libmesh_error_msg_if(ref_elem_map[base_type] == nullptr || type_in == INVALID_ELEM,
280 : "No reference elem data available for ElemType " << type_in
281 : << " = " << Utility::enum_to_string(type_in) << ".");
282 :
283 425869 : return *ref_elem_map[base_type];
284 : }
285 :
286 : std::pair<std::unique_ptr<Elem>, std::vector<std::unique_ptr<Node>>>
287 218896 : ideal_target (const ElemType type)
288 : {
289 : // Build target element
290 236986 : auto target_elem = Elem::build(type);
291 :
292 : // Volume of reference element
293 218896 : const auto ref_vol = target_elem->reference_elem()->volume();
294 :
295 : // Update the nodes of the target element, depending on type
296 18090 : const Real sqrt_2 = std::sqrt(Real(2));
297 18090 : const Real sqrt_3 = std::sqrt(Real(3));
298 54270 : std::vector<std::unique_ptr<Node>> owned_nodes;
299 :
300 236986 : const auto type_str = Utility::enum_to_string(type);
301 :
302 : // Elems deriving from Tri
303 218896 : if (type_str.compare(0, 3, "TRI") == 0)
304 : {
305 :
306 : // The target element will be an equilateral triangle with area equal to
307 : // the area of the reference element.
308 :
309 : // Equilateral triangle side length preserving area of the reference element
310 4213 : const auto side_length = std::sqrt(4. / sqrt_3 * ref_vol);
311 :
312 : // Define the nodal locations of the vertices
313 332 : const auto & s = side_length;
314 : // x y node_id
315 4545 : owned_nodes.emplace_back(Node::build(Point(0., 0.), 0));
316 4213 : owned_nodes.emplace_back(Node::build(Point(s, 0.), 1));
317 4545 : owned_nodes.emplace_back(Node::build(Point(0.5 * s, 0.5 * sqrt_3 * s), 2));
318 :
319 4213 : switch (type)
320 : {
321 170 : case TRI3:
322 : case TRISHELL3: {
323 : // Nothing to do here, vertices already added above
324 170 : break;
325 : }
326 :
327 1975 : case TRI6:
328 : case TRI7: {
329 : // Define the midpoint nodes of the equilateral triangle
330 : // x y node_id
331 1975 : owned_nodes.emplace_back(Node::build(Point(0.50 * s, 0.00), 3));
332 2137 : owned_nodes.emplace_back(Node::build(Point(0.75 * s, 0.25 * sqrt_3 * s), 4));
333 2137 : owned_nodes.emplace_back(Node::build(Point(0.25 * s, 0.25 * sqrt_3 * s), 5));
334 :
335 1975 : if (type == TRI7)
336 : {
337 : // Define the interior "bubble" node at the centroid
338 80 : const auto & on = owned_nodes;
339 1120 : owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[1] + *on[2]) / 3.), 6));
340 : }
341 :
342 162 : break;
343 : }
344 :
345 0 : default:
346 0 : libmesh_error_msg("Unsupported triangular element: " << type_str);
347 : break;
348 : }
349 : } // if Tri
350 :
351 : // Elems deriving from Prism
352 214683 : else if (type_str.compare(0, 5, "PRISM") == 0)
353 : {
354 :
355 : // The target element will be a prism with an equilateral triangular
356 : // base with volume equal to the volume of the reference element.
357 :
358 : // For an equilateral triangular base with side length s, the
359 : // base area is s^2 * sqrt(3) / 4.
360 : // The prism height that will result in equal face areas is
361 : // s * sqrt(3) / 4. We choose s such that the target element has
362 : // the same volume as the reference element:
363 : // v = (s^2 * sqrt(3) / 4) * (s * sqrt(3) / 4) = 3 * s^3 / 4
364 : // --> s = (16 * v / 3)^(1/3)
365 : // I have no particular motivation for imposing equal face areas,
366 : // so this can be updated if a more `optimal` target prism is
367 : // identified.
368 :
369 : // Side length that preserves the volume of the reference element
370 336 : const auto side_length = std::cbrt(16. * ref_vol / 3.);
371 : // Prism height with the property that all faces have equal area
372 336 : const auto target_height = 0.25 * side_length * sqrt_3;
373 :
374 14 : const auto & s = side_length;
375 14 : const auto & h = target_height;
376 : // x y z node_id
377 350 : owned_nodes.emplace_back(Node::build(Point(0., 0., 0.), 0));
378 336 : owned_nodes.emplace_back(Node::build(Point(s, 0., 0.), 1));
379 364 : owned_nodes.emplace_back(Node::build(Point(0.5 * s, 0.5 * sqrt_3 * s, 0.), 2));
380 350 : owned_nodes.emplace_back(Node::build(Point(0., 0., h), 3));
381 350 : owned_nodes.emplace_back(Node::build(Point(s, 0., h), 4));
382 336 : owned_nodes.emplace_back(Node::build(Point(0.5 * s, 0.5 * sqrt_3 * s, h), 5));
383 :
384 336 : if (type == PRISM15 || type == PRISM18 || type == PRISM20 || type == PRISM21)
385 : {
386 : // Define the edge midpoint nodes of the prism
387 10 : const auto & on = owned_nodes;
388 217 : owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[1]) / 2.), 6));
389 217 : owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[2]) / 2.), 7));
390 217 : owned_nodes.emplace_back(Node::build(Point((*on[2] + *on[0]) / 2.), 8));
391 217 : owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[3]) / 2.), 9));
392 217 : owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[4]) / 2.), 10));
393 217 : owned_nodes.emplace_back(Node::build(Point((*on[2] + *on[5]) / 2.), 11));
394 217 : owned_nodes.emplace_back(Node::build(Point((*on[3] + *on[4]) / 2.), 12));
395 217 : owned_nodes.emplace_back(Node::build(Point((*on[4] + *on[5]) / 2.), 13));
396 217 : owned_nodes.emplace_back(Node::build(Point((*on[5] + *on[3]) / 2.), 14));
397 :
398 207 : if (type == PRISM18 || type == PRISM20 || type == PRISM21)
399 : {
400 : // Define the rectangular face midpoint nodes of the prism
401 145 : owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[1] + *on[3] + *on[4]) / 4.), 15));
402 145 : owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[2] + *on[4] + *on[5]) / 4.), 16));
403 145 : owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[2] + *on[3] + *on[5]) / 4.), 17));
404 :
405 137 : if (type == PRISM20 || type == PRISM21)
406 : {
407 : // Define the triangular face midpoint nodes of the prism
408 72 : owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[1] + *on[2]) / 3.), 18));
409 72 : owned_nodes.emplace_back(Node::build(Point((*on[3] + *on[4] + *on[5]) / 3.), 19));
410 :
411 66 : if (type == PRISM21)
412 : // Define the interior point of the prism
413 52 : owned_nodes.emplace_back(Node::build(Point((*on[9] + *on[10] + *on[11]) / 3.), 20));
414 :
415 : }
416 10 : }
417 : }
418 :
419 129 : else if (type != PRISM6)
420 0 : libmesh_error_msg("Unsupported prism element: " << type_str);
421 :
422 : } // if Prism
423 :
424 : // Elems deriving from Pyramid
425 214347 : else if (type_str.compare(0, 7, "PYRAMID") == 0)
426 : {
427 :
428 : // The target element is a pyramid with an square base and
429 : // equilateral triangular sides with volume equal to the volume of the
430 : // reference element.
431 :
432 : // A pyramid with square base sidelength s and equilateral triangular
433 : // sides has height h = s / sqrt(2).
434 : // The volume is v = s^2 h / 3 = s^3 / ( 3 sqrt(2)).
435 : // Solving for s: s = (3 sqrt(2) v)^(1/3), where v is the volume of the
436 : // non-optimal reference element.
437 :
438 : // Side length that preserves the volume of the reference element
439 305 : const auto side_length = std::cbrt(3. * sqrt_2 * ref_vol);
440 : // Pyramid height with the property that all faces are equilateral triangles
441 305 : const auto target_height = side_length / sqrt_2;
442 :
443 10 : const auto & s = side_length;
444 10 : const auto & h = target_height;
445 :
446 : // x y z node_id
447 315 : owned_nodes.emplace_back(Node::build(Point(0., 0., 0.), 0));
448 315 : owned_nodes.emplace_back(Node::build(Point(s, 0., 0.), 1));
449 315 : owned_nodes.emplace_back(Node::build(Point(s, s, 0.), 2));
450 305 : owned_nodes.emplace_back(Node::build(Point(0., s, 0.), 3));
451 315 : owned_nodes.emplace_back(Node::build(Point(0.5 * s, 0.5 * s, h), 4));
452 :
453 305 : if (type == PYRAMID13 || type == PYRAMID14 || type == PYRAMID18)
454 : {
455 8 : const auto & on = owned_nodes;
456 : // Define the edge midpoint nodes of the pyramid
457 :
458 : // Base node to base node midpoint nodes
459 242 : owned_nodes.emplace_back(Node::build((*on[0] + *on[1]) / 2., 5));
460 242 : owned_nodes.emplace_back(Node::build((*on[1] + *on[2]) / 2., 6));
461 242 : owned_nodes.emplace_back(Node::build((*on[2] + *on[3]) / 2., 7));
462 242 : owned_nodes.emplace_back(Node::build((*on[3] + *on[0]) / 2., 8));
463 :
464 : // Base node to apex node midpoint nodes
465 242 : owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[4]) / 2.), 9));
466 242 : owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[4]) / 2.), 10));
467 242 : owned_nodes.emplace_back(Node::build(Point((*on[2] + *on[4]) / 2.), 11));
468 242 : owned_nodes.emplace_back(Node::build(Point((*on[3] + *on[4]) / 2.), 12));
469 :
470 234 : if (type == PYRAMID14 || type == PYRAMID18)
471 : {
472 : // Define the square face midpoint node of the pyramid
473 151 : owned_nodes.emplace_back(
474 169 : Node::build(Point((*on[0] + *on[1] + *on[2] + *on[3]) / 4.), 13));
475 :
476 163 : if (type == PYRAMID18)
477 : {
478 : // Define the triangular face nodes
479 96 : owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[1] + *on[4]) / 3.), 14));
480 96 : owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[2] + *on[4]) / 3.), 15));
481 96 : owned_nodes.emplace_back(Node::build(Point((*on[2] + *on[3] + *on[4]) / 3.), 16));
482 100 : owned_nodes.emplace_back(Node::build(Point((*on[3] + *on[0] + *on[4]) / 3.), 17));
483 : }
484 8 : }
485 : }
486 :
487 71 : else if (type != PYRAMID5)
488 0 : libmesh_error_msg("Unsupported pyramid element: " << type_str);
489 :
490 : } // if Pyramid
491 :
492 : // Elems deriving from Tet
493 214042 : else if (type_str.compare(0, 3, "TET") == 0)
494 : {
495 :
496 : // The ideal target element is a a regular tet with equilateral
497 : // triangles for all faces, with volume equal to the volume of the
498 : // reference element.
499 :
500 : // The volume of a tet is given by v = b * h / 3, where b is the area of
501 : // the base face and h is the height of the apex node. The area of an
502 : // equilateral triangle with side length s is b = sqrt(3) s^2 / 4.
503 : // For all faces to have side length s, the height of the apex node is
504 : // h = sqrt(2/3) * s. Then the volume is v = sqrt(2) * s^3 / 12.
505 : // Solving for s, the side length that will preserve the volume of the
506 : // reference element is s = (6 * sqrt(2) * v)^(1/3), where v is the volume
507 : // of the non-optimal reference element (i.e., a right tet).
508 :
509 : // Side length that preserves the volume of the reference element
510 193814 : const auto side_length = std::cbrt(6. * sqrt_2 * ref_vol);
511 : // tet height with the property that all faces are equilateral triangles
512 193814 : const auto target_height = sqrt_2 / sqrt_3 * side_length;
513 :
514 16136 : const auto & s = side_length;
515 16136 : const auto & h = target_height;
516 :
517 : // For regular tet
518 : // x y z node_id
519 209950 : owned_nodes.emplace_back(Node::build(Point(0., 0., 0.), 0));
520 193814 : owned_nodes.emplace_back(Node::build(Point(s, 0., 0.), 1));
521 209950 : owned_nodes.emplace_back(Node::build(Point(0.5 * s, 0.5 * sqrt_3 * s, 0.), 2));
522 209950 : owned_nodes.emplace_back(Node::build(Point(0.5 * s, sqrt_3 / 6. * s, h), 3));
523 :
524 193814 : if (type == TET10 || type == TET14)
525 : {
526 10758 : const auto & on = owned_nodes;
527 : // Define the edge midpoint nodes of the tet
528 :
529 : // Base node to base node midpoint nodes
530 139989 : owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[1]) / 2.), 4));
531 139989 : owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[2]) / 2.), 5));
532 139989 : owned_nodes.emplace_back(Node::build(Point((*on[2] + *on[0]) / 2.), 6));
533 : // Base node to apex node midpoint nodes
534 139989 : owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[3]) / 2.), 7));
535 139989 : owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[3]) / 2.), 8));
536 139989 : owned_nodes.emplace_back(Node::build(Point((*on[2] + *on[3]) / 2.), 9));
537 :
538 129231 : if (type == TET14)
539 : {
540 : // Define the face midpoint nodes of the tet
541 70028 : owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[1] + *on[2]) / 3.), 10));
542 70028 : owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[1] + *on[3]) / 3.), 11));
543 70028 : owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[2] + *on[3]) / 3.), 12));
544 75408 : owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[2] + *on[3]) / 3.), 13));
545 10758 : }
546 : }
547 :
548 64583 : else if (type != TET4)
549 0 : libmesh_error_msg("Unsupported tet element: " << type_str);
550 :
551 : } // if Tet
552 :
553 : // Set the target_elem equal to the reference elem
554 : else
555 334669 : for (const auto & node : target_elem->reference_elem()->node_ref_range())
556 363615 : owned_nodes.emplace_back(Node::build(node, node.id()));
557 :
558 : // Set nodes of target element
559 2368891 : for (const auto & node_ptr : owned_nodes)
560 2327905 : target_elem->set_node(node_ptr->id(), node_ptr.get());
561 :
562 18090 : libmesh_assert(relative_fuzzy_equals(target_elem->volume(), ref_vol, TOLERANCE));
563 :
564 236986 : return std::make_pair(std::move(target_elem), std::move(owned_nodes));
565 182716 : }
566 : } // namespace ReferenceElem
567 : } // namespace libMesh
|
