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 : // libmesh includes
21 : #include "libmesh/mesh_generation.h"
22 : #include "libmesh/unstructured_mesh.h"
23 : #include "libmesh/mesh_refinement.h"
24 : #include "libmesh/edge_edge2.h"
25 : #include "libmesh/edge_edge3.h"
26 : #include "libmesh/edge_edge4.h"
27 : #include "libmesh/face_tri3.h"
28 : #include "libmesh/face_tri6.h"
29 : #include "libmesh/face_tri7.h"
30 : #include "libmesh/face_quad4.h"
31 : #include "libmesh/face_quad8.h"
32 : #include "libmesh/face_quad9.h"
33 : #include "libmesh/face_c0polygon.h"
34 : #include "libmesh/cell_c0polyhedron.h"
35 : #include "libmesh/cell_hex8.h"
36 : #include "libmesh/cell_hex20.h"
37 : #include "libmesh/cell_hex27.h"
38 : #include "libmesh/cell_prism6.h"
39 : #include "libmesh/cell_prism15.h"
40 : #include "libmesh/cell_prism18.h"
41 : #include "libmesh/cell_prism20.h"
42 : #include "libmesh/cell_prism21.h"
43 : #include "libmesh/cell_tet4.h"
44 : #include "libmesh/cell_pyramid5.h"
45 : #include "libmesh/libmesh_logging.h"
46 : #include "libmesh/boundary_info.h"
47 : #include "libmesh/remote_elem.h"
48 : #include "libmesh/sphere.h"
49 : #include "libmesh/mesh_modification.h"
50 : #include "libmesh/mesh_smoother_laplace.h"
51 : #include "libmesh/node_elem.h"
52 : #include "libmesh/vector_value.h"
53 : #include "libmesh/function_base.h"
54 : #include "libmesh/enum_order.h"
55 : #include "libmesh/int_range.h"
56 : #include "libmesh/parallel.h"
57 : #include "libmesh/parallel_ghost_sync.h"
58 : #include "libmesh/enum_to_string.h"
59 :
60 : // C++ includes
61 : #include <array>
62 : #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
63 : #include <cmath> // for std::sqrt
64 : #include <unordered_set>
65 :
66 :
67 : namespace libMesh
68 : {
69 :
70 : namespace MeshTools {
71 : namespace Generation {
72 : namespace Private {
73 : /**
74 : * A useful inline function which replaces the macros
75 : * used previously. Not private since this is a namespace,
76 : * but would be if this were a class. The first one returns
77 : * the proper node number for 2D elements while the second
78 : * one returns the node number for 3D elements.
79 : */
80 : inline
81 1538884 : unsigned int idx(const ElemType type,
82 : const unsigned int nx,
83 : const unsigned int i,
84 : const unsigned int j)
85 : {
86 1109486 : switch(type)
87 : {
88 317784 : case INVALID_ELEM:
89 : case QUAD4:
90 : case QUADSHELL4:
91 : case TRI3:
92 : case TRISHELL3:
93 : {
94 317784 : return i + j*(nx+1);
95 : }
96 :
97 1221100 : case QUAD8:
98 : case QUADSHELL8:
99 : case QUAD9:
100 : case QUADSHELL9:
101 : case TRI6:
102 : case TRI7:
103 : {
104 1221100 : return i + j*(2*nx+1);
105 : }
106 :
107 0 : default:
108 0 : libmesh_error_msg("ERROR: Unrecognized 2D element type == " << Utility::enum_to_string(type));
109 : }
110 :
111 : return libMesh::invalid_uint;
112 : }
113 :
114 :
115 :
116 : // Same as the function above, but for 3D elements
117 : inline
118 3991142 : unsigned int idx(const ElemType type,
119 : const unsigned int nx,
120 : const unsigned int ny,
121 : const unsigned int i,
122 : const unsigned int j,
123 : const unsigned int k)
124 : {
125 2897766 : switch(type)
126 : {
127 834884 : case INVALID_ELEM:
128 : case HEX8:
129 : case PRISM6:
130 : case C0POLYHEDRON:
131 : {
132 834884 : return i + (nx+1)*(j + k*(ny+1));
133 : }
134 :
135 3156258 : case HEX20:
136 : case HEX27:
137 : case TET4: // TET4's are created from an initial HEX27 discretization
138 : case TET10: // TET10's are created from an initial HEX27 discretization
139 : case TET14: // TET14's are created from an initial HEX27 discretization
140 : case PYRAMID5: // PYRAMID5's are created from an initial HEX27 discretization
141 : case PYRAMID13:
142 : case PYRAMID14:
143 : case PYRAMID18:
144 : case PRISM15:
145 : case PRISM18:
146 : case PRISM20:
147 : case PRISM21:
148 : {
149 3156258 : return i + (2*nx+1)*(j + k*(2*ny+1));
150 : }
151 :
152 0 : default:
153 0 : libmesh_error_msg("ERROR: Unrecognized element type == " << Utility::enum_to_string(type));
154 : }
155 :
156 : return libMesh::invalid_uint;
157 : }
158 :
159 :
160 : /**
161 : * This object is passed to MeshTools::Modification::redistribute() to
162 : * redistribute the points on a uniform grid into the Gauss-Lobatto
163 : * points on the actual grid.
164 : */
165 : class GaussLobattoRedistributionFunction : public FunctionBase<Real>
166 : {
167 : public:
168 : /**
169 : * Constructor.
170 : */
171 0 : GaussLobattoRedistributionFunction(unsigned int nx,
172 : Real xmin,
173 : Real xmax,
174 : unsigned int ny=0,
175 : Real ymin=0,
176 : Real ymax=0,
177 : unsigned int nz=0,
178 : Real zmin=0,
179 0 : Real zmax=0) :
180 0 : FunctionBase<Real>(nullptr)
181 : {
182 0 : _nelem.resize(3);
183 0 : _nelem[0] = nx;
184 0 : _nelem[1] = ny;
185 0 : _nelem[2] = nz;
186 :
187 0 : _mins.resize(3);
188 0 : _mins[0] = xmin;
189 0 : _mins[1] = ymin;
190 0 : _mins[2] = zmin;
191 :
192 0 : _widths.resize(3);
193 0 : _widths[0] = xmax - xmin;
194 0 : _widths[1] = ymax - ymin;
195 0 : _widths[2] = zmax - zmin;
196 :
197 : // Precompute the cosine values.
198 0 : _cosines.resize(3);
199 0 : for (unsigned dir=0; dir<3; ++dir)
200 0 : if (_nelem[dir] != 0)
201 : {
202 0 : _cosines[dir].resize(_nelem[dir]+1);
203 0 : for (auto i : index_range(_cosines[dir]))
204 0 : _cosines[dir][i] = std::cos(libMesh::pi * Real(i) / _nelem[dir]);
205 : }
206 0 : }
207 :
208 : /**
209 : * The 5 special functions can be defaulted for this class.
210 : */
211 : GaussLobattoRedistributionFunction (GaussLobattoRedistributionFunction &&) = default;
212 0 : GaussLobattoRedistributionFunction (const GaussLobattoRedistributionFunction &) = default;
213 : GaussLobattoRedistributionFunction & operator= (const GaussLobattoRedistributionFunction &) = default;
214 : GaussLobattoRedistributionFunction & operator= (GaussLobattoRedistributionFunction &&) = default;
215 0 : virtual ~GaussLobattoRedistributionFunction () = default;
216 :
217 : /**
218 : * We must provide a way to clone ourselves to satisfy the pure
219 : * virtual interface. We use the autogenerated copy constructor.
220 : */
221 0 : virtual std::unique_ptr<FunctionBase<Real>> clone () const override
222 : {
223 0 : return std::make_unique<GaussLobattoRedistributionFunction>(*this);
224 : }
225 :
226 : /**
227 : * This is the actual function that
228 : * MeshTools::Modification::redistribute() calls. Moves the points
229 : * of the grid to the Gauss-Lobatto points.
230 : */
231 0 : virtual void operator() (const Point & p,
232 : const Real /*time*/,
233 : DenseVector<Real> & output) override
234 : {
235 0 : output.resize(3);
236 :
237 0 : for (unsigned dir=0; dir<3; ++dir)
238 0 : if (_nelem[dir] != 0)
239 : {
240 : // Figure out the index of the current point.
241 0 : Real float_index = (p(dir) - _mins[dir]) * _nelem[dir] / _widths[dir];
242 :
243 : // std::modf separates the fractional and integer parts of the index.
244 0 : Real integer_part_f = 0;
245 0 : const Real fractional_part = std::modf(float_index, &integer_part_f);
246 :
247 0 : const int integer_part = int(integer_part_f);
248 :
249 : // Vertex node?
250 0 : if (std::abs(fractional_part) < TOLERANCE || std::abs(fractional_part - 1.0) < TOLERANCE)
251 : {
252 0 : int index = int(round(float_index));
253 :
254 : // Move node to the Gauss-Lobatto position.
255 0 : output(dir) = _mins[dir] + _widths[dir] * 0.5 * (1.0 - _cosines[dir][index]);
256 : }
257 :
258 : // Mid-edge (quadratic) node?
259 0 : else if (std::abs(fractional_part - 0.5) < TOLERANCE)
260 : {
261 : // Move node to the Gauss-Lobatto position, which is the average of
262 : // the node to the left and the node to the right.
263 0 : output(dir) = _mins[dir] + _widths[dir] * 0.5 *
264 0 : (1.0 - 0.5*(_cosines[dir][integer_part] + _cosines[dir][integer_part+1]));
265 : }
266 :
267 : // 1D only: Left interior (cubic) node?
268 0 : else if (std::abs(fractional_part - 1./3.) < TOLERANCE)
269 : {
270 : // Move node to the Gauss-Lobatto position, which is
271 : // 2/3*left_vertex + 1/3*right_vertex.
272 0 : output(dir) = _mins[dir] + _widths[dir] * 0.5 *
273 0 : (1.0 - 2./3.*_cosines[dir][integer_part] - 1./3.*_cosines[dir][integer_part+1]);
274 : }
275 :
276 : // 1D only: Right interior (cubic) node?
277 0 : else if (std::abs(fractional_part - 2./3.) < TOLERANCE)
278 : {
279 : // Move node to the Gauss-Lobatto position, which is
280 : // 1/3*left_vertex + 2/3*right_vertex.
281 0 : output(dir) = _mins[dir] + _widths[dir] * 0.5 *
282 0 : (1.0 - 1./3.*_cosines[dir][integer_part] - 2./3.*_cosines[dir][integer_part+1]);
283 : }
284 :
285 : else
286 0 : libmesh_error_msg("Cannot redistribute node: " << p);
287 : }
288 0 : }
289 :
290 : /**
291 : * We must also override operator() which returns a Real, but this function
292 : * should never be called, so it's left unimplemented.
293 : */
294 0 : virtual Real operator() (const Point & /*p*/,
295 : const Real /*time*/) override
296 : {
297 0 : libmesh_not_implemented();
298 : }
299 :
300 : protected:
301 : // Stored data
302 : std::vector<Real> _mins;
303 : std::vector<unsigned int> _nelem;
304 : std::vector<Real> _widths;
305 :
306 : // Precomputed values
307 : std::vector<std::vector<Real>> _cosines;
308 : };
309 :
310 :
311 : } // namespace Private
312 : } // namespace Generation
313 : } // namespace MeshTools
314 :
315 : // ------------------------------------------------------------
316 : // MeshTools::Generation function for mesh generation
317 27567 : void MeshTools::Generation::build_cube(UnstructuredMesh & mesh,
318 : const unsigned int nx,
319 : const unsigned int ny,
320 : const unsigned int nz,
321 : const Real xmin, const Real xmax,
322 : const Real ymin, const Real ymax,
323 : const Real zmin, const Real zmax,
324 : const ElemType type,
325 : const bool gauss_lobatto_grid)
326 : {
327 15792 : LOG_SCOPE("build_cube()", "MeshTools::Generation");
328 :
329 : // Declare that we are using the indexing utility routine
330 : // in the "Private" part of our current namespace. If this doesn't
331 : // work in GCC 2.95.3 we can either remove it or stop supporting
332 : // 2.95.3 altogether.
333 : // Changing this to import the whole namespace... just importing idx
334 : // causes an internal compiler error for Intel Compiler 11.0 on Linux
335 : // in debug mode.
336 : using namespace MeshTools::Generation::Private;
337 :
338 : // Clear the mesh and start from scratch
339 27567 : mesh.clear();
340 :
341 7896 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
342 :
343 27567 : if (nz != 0)
344 : {
345 12868 : mesh.set_mesh_dimension(3);
346 12868 : mesh.set_spatial_dimension(3);
347 : }
348 14699 : else if (ny != 0)
349 : {
350 11046 : mesh.set_mesh_dimension(2);
351 11046 : mesh.set_spatial_dimension(2);
352 : }
353 3653 : else if (nx != 0)
354 : {
355 3429 : mesh.set_mesh_dimension(1);
356 3429 : mesh.set_spatial_dimension(1);
357 : }
358 : else
359 : {
360 : // Will we get here?
361 224 : mesh.set_mesh_dimension(0);
362 224 : mesh.set_spatial_dimension(0);
363 : }
364 :
365 27567 : switch (mesh.mesh_dimension())
366 : {
367 : //---------------------------------------------------------------------
368 : // Build a 0D point
369 224 : case 0:
370 : {
371 64 : libmesh_assert_equal_to (nx, 0);
372 64 : libmesh_assert_equal_to (ny, 0);
373 64 : libmesh_assert_equal_to (nz, 0);
374 :
375 64 : libmesh_assert (type == INVALID_ELEM || type == NODEELEM);
376 :
377 : // Build one nodal element for the mesh
378 288 : mesh.add_point (Point(0, 0, 0), 0);
379 224 : Elem * elem = mesh.add_elem(Elem::build(NODEELEM));
380 224 : elem->set_node(0, mesh.node_ptr(0));
381 :
382 160 : break;
383 : }
384 :
385 :
386 :
387 : //---------------------------------------------------------------------
388 : // Build a 1D line
389 3429 : case 1:
390 : {
391 982 : libmesh_assert_not_equal_to (nx, 0);
392 982 : libmesh_assert_equal_to (ny, 0);
393 982 : libmesh_assert_equal_to (nz, 0);
394 982 : libmesh_assert_less (xmin, xmax);
395 :
396 : // Reserve elements
397 : switch (type)
398 : {
399 3429 : case INVALID_ELEM:
400 : case EDGE2:
401 : case EDGE3:
402 : case EDGE4:
403 : {
404 3429 : mesh.reserve_elem (nx);
405 982 : break;
406 : }
407 :
408 0 : default:
409 0 : libmesh_error_msg("ERROR: Unrecognized 1D element type == " << Utility::enum_to_string(type));
410 : }
411 :
412 : // Reserve nodes
413 : switch (type)
414 : {
415 1121 : case INVALID_ELEM:
416 : case EDGE2:
417 : {
418 1121 : mesh.reserve_nodes(nx+1);
419 799 : break;
420 : }
421 :
422 2105 : case EDGE3:
423 : {
424 2105 : mesh.reserve_nodes(2*nx+1);
425 1503 : break;
426 : }
427 :
428 203 : case EDGE4:
429 : {
430 203 : mesh.reserve_nodes(3*nx+1);
431 145 : break;
432 : }
433 :
434 0 : default:
435 0 : libmesh_error_msg("ERROR: Unrecognized 1D element type == " << Utility::enum_to_string(type));
436 : }
437 :
438 :
439 : // Build the nodes, depends on whether we're using linears,
440 : // quadratics or cubics and whether using uniform grid or Gauss-Lobatto
441 982 : unsigned int node_id = 0;
442 : switch(type)
443 : {
444 322 : case INVALID_ELEM:
445 : case EDGE2:
446 : {
447 29021 : for (unsigned int i=0; i<=nx; i++)
448 : {
449 37056 : const Node * const node = mesh.add_point (Point(static_cast<Real>(i)/nx, 0, 0), node_id++);
450 27900 : if (i == 0)
451 1121 : boundary_info.add_node(node, 0);
452 27900 : if (i == nx)
453 1121 : boundary_info.add_node(node, 1);
454 : }
455 :
456 322 : break;
457 : }
458 :
459 602 : case EDGE3:
460 : {
461 11226 : for (unsigned int i=0; i<=2*nx; i++)
462 : {
463 11739 : const Node * const node = mesh.add_point (Point(static_cast<Real>(i)/(2*nx), 0, 0), node_id++);
464 9121 : if (i == 0)
465 2105 : boundary_info.add_node(node, 0);
466 9121 : if (i == 2*nx)
467 2105 : boundary_info.add_node(node, 1);
468 : }
469 602 : break;
470 : }
471 :
472 58 : case EDGE4:
473 : {
474 2023 : for (unsigned int i=0; i<=3*nx; i++)
475 : {
476 2340 : const Node * const node = mesh.add_point (Point(static_cast<Real>(i)/(3*nx), 0, 0), node_id++);
477 1820 : if (i == 0)
478 203 : boundary_info.add_node(node, 0);
479 1820 : if (i == 3*nx)
480 203 : boundary_info.add_node(node, 1);
481 : }
482 :
483 58 : break;
484 : }
485 :
486 0 : default:
487 0 : libmesh_error_msg("ERROR: Unrecognized 1D element type == " << Utility::enum_to_string(type));
488 :
489 : }
490 :
491 : // Build the elements of the mesh
492 : switch(type)
493 : {
494 322 : case INVALID_ELEM:
495 : case EDGE2:
496 : {
497 27900 : for (unsigned int i=0; i<nx; i++)
498 : {
499 26779 : Elem * elem = mesh.add_elem(Elem::build_with_id(EDGE2, i));
500 26779 : elem->set_node(0, mesh.node_ptr(i));
501 26779 : elem->set_node(1, mesh.node_ptr(i+1));
502 :
503 26779 : if (i == 0)
504 1121 : boundary_info.add_side(elem, 0, 0);
505 :
506 26779 : if (i == (nx-1))
507 1121 : boundary_info.add_side(elem, 1, 1);
508 :
509 : }
510 322 : break;
511 : }
512 :
513 602 : case EDGE3:
514 : {
515 5613 : for (unsigned int i=0; i<nx; i++)
516 : {
517 3508 : Elem * elem = mesh.add_elem(Elem::build_with_id(EDGE3, i));
518 3508 : elem->set_node(0, mesh.node_ptr(2*i));
519 3508 : elem->set_node(2, mesh.node_ptr(2*i+1));
520 3508 : elem->set_node(1, mesh.node_ptr(2*i+2));
521 :
522 3508 : if (i == 0)
523 2105 : boundary_info.add_side(elem, 0, 0);
524 :
525 3508 : if (i == (nx-1))
526 2105 : boundary_info.add_side(elem, 1, 1);
527 : }
528 602 : break;
529 : }
530 :
531 58 : case EDGE4:
532 : {
533 742 : for (unsigned int i=0; i<nx; i++)
534 : {
535 539 : Elem * elem = mesh.add_elem(Elem::build_with_id(EDGE4, i));
536 539 : elem->set_node(0, mesh.node_ptr(3*i));
537 539 : elem->set_node(2, mesh.node_ptr(3*i+1));
538 539 : elem->set_node(3, mesh.node_ptr(3*i+2));
539 539 : elem->set_node(1, mesh.node_ptr(3*i+3));
540 :
541 539 : if (i == 0)
542 203 : boundary_info.add_side(elem, 0, 0);
543 :
544 539 : if (i == (nx-1))
545 203 : boundary_info.add_side(elem, 1, 1);
546 : }
547 58 : break;
548 : }
549 :
550 0 : default:
551 0 : libmesh_error_msg("ERROR: Unrecognized 1D element type == " << Utility::enum_to_string(type));
552 : }
553 :
554 : // Move the nodes to their final locations.
555 3429 : if (gauss_lobatto_grid)
556 : {
557 0 : GaussLobattoRedistributionFunction func(nx, xmin, xmax);
558 0 : MeshTools::Modification::redistribute(mesh, func);
559 0 : }
560 : else // !gauss_lobatto_grid
561 : {
562 44717 : for (Node * node : mesh.node_ptr_range())
563 40306 : (*node)(0) = (*node)(0)*(xmax-xmin) + xmin;
564 : }
565 :
566 : // Add sideset names to boundary info
567 3429 : boundary_info.sideset_name(0) = "left";
568 3429 : boundary_info.sideset_name(1) = "right";
569 :
570 : // Add nodeset names to boundary info
571 3429 : boundary_info.nodeset_name(0) = "left";
572 3429 : boundary_info.nodeset_name(1) = "right";
573 :
574 982 : break;
575 : }
576 :
577 :
578 :
579 :
580 :
581 :
582 :
583 :
584 :
585 :
586 : //---------------------------------------------------------------------
587 : // Build a 2D quadrilateral
588 11046 : case 2:
589 : {
590 3164 : libmesh_assert_not_equal_to (nx, 0);
591 3164 : libmesh_assert_not_equal_to (ny, 0);
592 3164 : libmesh_assert_equal_to (nz, 0);
593 3164 : libmesh_assert_less (xmin, xmax);
594 3164 : libmesh_assert_less (ymin, ymax);
595 :
596 : // Reserve elements. The TRI3 and TRI6 meshes
597 : // have twice as many elements...
598 : switch (type)
599 : {
600 5601 : case INVALID_ELEM:
601 : case QUAD4:
602 : case QUADSHELL4:
603 : case QUAD8:
604 : case QUADSHELL8:
605 : case QUAD9:
606 : case QUADSHELL9:
607 : {
608 5601 : mesh.reserve_elem (nx*ny);
609 3993 : break;
610 : }
611 :
612 5382 : case TRI3:
613 : case TRISHELL3:
614 : case TRI6:
615 : case TRI7:
616 : {
617 5382 : mesh.reserve_elem (2*nx*ny);
618 3844 : break;
619 : }
620 :
621 63 : case C0POLYGON:
622 : {
623 63 : mesh.reserve_elem ((nx + 1) * (ny + 1));
624 45 : break;
625 : }
626 :
627 0 : default:
628 0 : libmesh_error_msg("ERROR: Unrecognized 2D element type == " << Utility::enum_to_string(type));
629 : }
630 :
631 :
632 :
633 : // Reserve nodes. The quadratic element types
634 : // need to reserve more nodes than the linear types.
635 : switch (type)
636 : {
637 2945 : case INVALID_ELEM:
638 : case QUAD4:
639 : case QUADSHELL4:
640 : case TRI3:
641 : case TRISHELL3:
642 : {
643 2945 : mesh.reserve_nodes( (nx+1)*(ny+1) );
644 2095 : break;
645 : }
646 :
647 5948 : case QUAD8:
648 : case QUADSHELL8:
649 : case QUAD9:
650 : case QUADSHELL9:
651 : case TRI6:
652 : {
653 5948 : mesh.reserve_nodes( (2*nx+1)*(2*ny+1) );
654 4248 : break;
655 : }
656 :
657 2090 : case TRI7:
658 : {
659 2090 : mesh.reserve_nodes( (2*nx+1)*(2*ny+1) + 2*nx*ny );
660 1494 : break;
661 : }
662 63 : case C0POLYGON:
663 : {
664 63 : mesh.reserve_nodes (4 + 3*nx*ny + 2*nx + 2*ny);
665 45 : break;
666 : }
667 :
668 0 : default:
669 0 : libmesh_error_msg("ERROR: Unrecognized 2D element type == " << Utility::enum_to_string(type));
670 : }
671 :
672 :
673 :
674 : // Build the nodes. Depends on whether you are using a linear
675 : // or quadratic element, and whether you are using a uniform
676 : // grid or the Gauss-Lobatto grid points.
677 3164 : unsigned int node_id = 0;
678 : switch (type)
679 : {
680 850 : case INVALID_ELEM:
681 : case QUAD4:
682 : case QUADSHELL4:
683 : case TRI3:
684 : case TRISHELL3:
685 : {
686 13496 : for (unsigned int j=0; j<=ny; j++)
687 105778 : for (unsigned int i=0; i<=nx; i++)
688 : {
689 : const Node * const node =
690 134086 : mesh.add_point(Point(static_cast<Real>(i) / static_cast<Real>(nx),
691 95227 : static_cast<Real>(j) / static_cast<Real>(ny),
692 : 0.),
693 85398 : node_id++);
694 95227 : if (j == 0)
695 9963 : boundary_info.add_node(node, 0);
696 95227 : if (j == ny)
697 9963 : boundary_info.add_node(node, 2);
698 95227 : if (i == 0)
699 10551 : boundary_info.add_node(node, 3);
700 95227 : if (i == nx)
701 10551 : boundary_info.add_node(node, 1);
702 : }
703 :
704 850 : break;
705 : }
706 :
707 2296 : case QUAD8:
708 : case QUADSHELL8:
709 : case QUAD9:
710 : case QUADSHELL9:
711 : case TRI6:
712 : case TRI7:
713 : {
714 49328 : for (unsigned int j=0; j<=(2*ny); j++)
715 605456 : for (unsigned int i=0; i<=(2*nx); i++)
716 : {
717 : const Node * const node =
718 817116 : mesh.add_point(Point(static_cast<Real>(i) / static_cast<Real>(2 * nx),
719 564166 : static_cast<Real>(j) / static_cast<Real>(2 * ny),
720 : 0),
721 487960 : node_id++);
722 564166 : if (j == 0)
723 42602 : boundary_info.add_node(node, 0);
724 564166 : if (j == 2*ny)
725 42602 : boundary_info.add_node(node, 2);
726 564166 : if (i == 0)
727 41290 : boundary_info.add_node(node, 3);
728 564166 : if (i == 2*nx)
729 41290 : boundary_info.add_node(node, 1);
730 : }
731 :
732 : // We'll add any interior Tri7 nodes last, to keep from
733 : // messing with our idx function
734 8038 : if (type == TRI7)
735 5597 : for (unsigned int j=0; j<(3*ny); j += 3)
736 23664 : for (unsigned int i=0; i<(3*nx); i += 3)
737 : {
738 : // The bottom-right triangle's center node
739 29510 : mesh.add_point(Point(static_cast<Real>(i+2) / static_cast<Real>(3 * nx),
740 20157 : static_cast<Real>(j+1) / static_cast<Real>(3 * ny),
741 : 0),
742 17456 : node_id++);
743 : // The top-left triangle's center node
744 20157 : mesh.add_point(Point(static_cast<Real>(i+1) / static_cast<Real>(3 * nx),
745 20157 : static_cast<Real>(j+2) / static_cast<Real>(3 * ny),
746 : 0),
747 17456 : node_id++);
748 : }
749 :
750 2296 : break;
751 : }
752 :
753 18 : case C0POLYGON:
754 : {
755 : // we create the nodes at the same time as the elements
756 18 : break;
757 : }
758 :
759 0 : default:
760 0 : libmesh_error_msg("ERROR: Unrecognized 2D element type == " << Utility::enum_to_string(type));
761 : }
762 :
763 :
764 :
765 :
766 :
767 :
768 : // Build the elements. Each one is a bit different.
769 3164 : unsigned int elem_id = 0;
770 : switch (type)
771 : {
772 :
773 528 : case INVALID_ELEM:
774 : case QUAD4:
775 : case QUADSHELL4:
776 : {
777 7739 : for (unsigned int j=0; j<ny; j++)
778 80002 : for (unsigned int i=0; i<nx; i++)
779 : {
780 75651 : Elem * elem = mesh.add_elem(Elem::build_with_id(type == INVALID_ELEM ? QUAD4 : type, elem_id++));
781 74082 : elem->set_node(0, mesh.node_ptr(idx(type,nx,i,j) ));
782 74082 : elem->set_node(1, mesh.node_ptr(idx(type,nx,i+1,j) ));
783 74082 : elem->set_node(2, mesh.node_ptr(idx(type,nx,i+1,j+1)));
784 74082 : elem->set_node(3, mesh.node_ptr(idx(type,nx,i,j+1) ));
785 :
786 74082 : if (j == 0)
787 5318 : boundary_info.add_side(elem, 0, 0);
788 :
789 74082 : if (j == (ny-1))
790 5318 : boundary_info.add_side(elem, 2, 2);
791 :
792 74082 : if (i == 0)
793 5920 : boundary_info.add_side(elem, 3, 3);
794 :
795 74082 : if (i == (nx-1))
796 5920 : boundary_info.add_side(elem, 1, 1);
797 : }
798 528 : break;
799 : }
800 :
801 :
802 322 : case TRI3:
803 : case TRISHELL3:
804 : {
805 2812 : for (unsigned int j=0; j<ny; j++)
806 5262 : for (unsigned int i=0; i<nx; i++)
807 : {
808 : // Add first Tri3
809 3576 : Elem * elem = mesh.add_elem(Elem::build_with_id(type, elem_id++));
810 3576 : elem->set_node(0, mesh.node_ptr(idx(type,nx,i,j) ));
811 3576 : elem->set_node(1, mesh.node_ptr(idx(type,nx,i+1,j) ));
812 3576 : elem->set_node(2, mesh.node_ptr(idx(type,nx,i+1,j+1)));
813 :
814 3576 : if (j == 0)
815 1700 : boundary_info.add_side(elem, 0, 0);
816 :
817 3576 : if (i == (nx-1))
818 1686 : boundary_info.add_side(elem, 1, 1);
819 :
820 : // Add second Tri3
821 3576 : elem = mesh.add_elem(Elem::build_with_id(type, elem_id++));
822 3576 : elem->set_node(0, mesh.node_ptr(idx(type,nx,i,j) ));
823 3576 : elem->set_node(1, mesh.node_ptr(idx(type,nx,i+1,j+1)));
824 3576 : elem->set_node(2, mesh.node_ptr(idx(type,nx,i,j+1) ));
825 :
826 3576 : if (j == (ny-1))
827 1700 : boundary_info.add_side(elem, 1, 2);
828 :
829 3576 : if (i == 0)
830 1686 : boundary_info.add_side(elem, 2, 3);
831 : }
832 322 : break;
833 : }
834 :
835 :
836 :
837 1080 : case QUAD8:
838 : case QUADSHELL8:
839 : case QUAD9:
840 : case QUADSHELL9:
841 : {
842 12787 : for (unsigned int j=0; j<(2*ny); j += 2)
843 84771 : for (unsigned int i=0; i<(2*nx); i += 2)
844 : {
845 75766 : Elem * elem = mesh.add_elem(Elem::build_with_id(type, elem_id++));
846 75766 : elem->set_node(0, mesh.node_ptr(idx(type,nx,i,j) ));
847 75766 : elem->set_node(1, mesh.node_ptr(idx(type,nx,i+2,j) ));
848 75766 : elem->set_node(2, mesh.node_ptr(idx(type,nx,i+2,j+2)));
849 75766 : elem->set_node(3, mesh.node_ptr(idx(type,nx,i,j+2) ));
850 75766 : elem->set_node(4, mesh.node_ptr(idx(type,nx,i+1,j) ));
851 75766 : elem->set_node(5, mesh.node_ptr(idx(type,nx,i+2,j+1)));
852 75766 : elem->set_node(6, mesh.node_ptr(idx(type,nx,i+1,j+2)));
853 75766 : elem->set_node(7, mesh.node_ptr(idx(type,nx,i,j+1) ));
854 :
855 75766 : if (type == QUAD9 || type == QUADSHELL9)
856 59228 : elem->set_node(8, mesh.node_ptr(idx(type,nx,i+1,j+1)));
857 :
858 75766 : if (j == 0)
859 9558 : boundary_info.add_side(elem, 0, 0);
860 :
861 75766 : if (j == 2*(ny-1))
862 9558 : boundary_info.add_side(elem, 2, 2);
863 :
864 75766 : if (i == 0)
865 9005 : boundary_info.add_side(elem, 3, 3);
866 :
867 75766 : if (i == 2*(nx-1))
868 9005 : boundary_info.add_side(elem, 1, 1);
869 : }
870 1080 : break;
871 : }
872 :
873 :
874 1216 : case TRI6:
875 : case TRI7:
876 : {
877 11877 : for (unsigned int j=0; j<(2*ny); j += 2)
878 53933 : for (unsigned int i=0; i<(2*nx); i += 2)
879 : {
880 : // Add first Tri in the bottom-right of its quad
881 46312 : Elem * elem = mesh.add_elem(Elem::build_with_id(type, elem_id++));
882 46312 : elem->set_node(0, mesh.node_ptr(idx(type,nx,i,j) ));
883 46312 : elem->set_node(1, mesh.node_ptr(idx(type,nx,i+2,j) ));
884 46312 : elem->set_node(2, mesh.node_ptr(idx(type,nx,i+2,j+2)));
885 46312 : elem->set_node(3, mesh.node_ptr(idx(type,nx,i+1,j) ));
886 46312 : elem->set_node(4, mesh.node_ptr(idx(type,nx,i+2,j+1)));
887 46312 : elem->set_node(5, mesh.node_ptr(idx(type,nx,i+1,j+1)));
888 :
889 46312 : if (type == TRI7)
890 20157 : elem->set_node(6, mesh.node_ptr(elem->id()+(2*nx+1)*(2*ny+1)));
891 :
892 46312 : if (j == 0)
893 7724 : boundary_info.add_side(elem, 0, 0);
894 :
895 46312 : if (i == 2*(nx-1))
896 7621 : boundary_info.add_side(elem, 1, 1);
897 :
898 : // Add second Tri in the top left of its quad
899 46312 : elem = mesh.add_elem(Elem::build_with_id(type, elem_id++));
900 46312 : elem->set_node(0, mesh.node_ptr(idx(type,nx,i,j) ));
901 46312 : elem->set_node(1, mesh.node_ptr(idx(type,nx,i+2,j+2)));
902 46312 : elem->set_node(2, mesh.node_ptr(idx(type,nx,i,j+2) ));
903 46312 : elem->set_node(3, mesh.node_ptr(idx(type,nx,i+1,j+1)));
904 46312 : elem->set_node(4, mesh.node_ptr(idx(type,nx,i+1,j+2)));
905 46312 : elem->set_node(5, mesh.node_ptr(idx(type,nx,i,j+1) ));
906 :
907 46312 : if (type == TRI7)
908 20157 : elem->set_node(6, mesh.node_ptr(elem->id()+(2*nx+1)*(2*ny+1)));
909 :
910 46312 : if (j == 2*(ny-1))
911 7724 : boundary_info.add_side(elem, 1, 2);
912 :
913 46312 : if (i == 0)
914 7621 : boundary_info.add_side(elem, 2, 3);
915 : }
916 1216 : break;
917 : };
918 :
919 18 : case C0POLYGON:
920 : {
921 : // Build a 2D paving using hexagons (center), quads (part of y-boundary)
922 : // and triangles (x-boundaries).
923 : // Vector to re-use previously created nodes
924 36 : std::vector<Node *> node_list;
925 :
926 : // Start with a layer of triangles on the boundary
927 63 : const auto dx_tri = Real(1) / nx;
928 63 : const auto dy_tri = Real(1) / (ny + 1);
929 45 : std::unique_ptr<Elem> new_elem;
930 448 : for (const auto i : make_range(nx + 1))
931 : {
932 : // Make new nodes for bottom layer of triangles
933 : Node *node0, *node1, *node2;
934 385 : if (i == 0)
935 : {
936 81 : node0 = mesh.add_point(Point(0., 0, 0.));
937 81 : node1 = mesh.add_point(Point(0., dy_tri / 2., 0.));
938 81 : node2 = mesh.add_point(Point(dx_tri / 2., 0., 0.));
939 63 : node_list.push_back(node0);
940 63 : node_list.push_back(node1);
941 63 : node_list.push_back(node2);
942 : }
943 322 : else if (i < nx)
944 : {
945 259 : node0 = node_list.back();
946 333 : node1 = mesh.add_point(Point((i)*dx_tri, dy_tri / 2., 0.));
947 333 : node2 = mesh.add_point(Point((i + 1. / 2.) * dx_tri, 0., 0.));
948 259 : node_list.push_back(node1);
949 259 : node_list.push_back(node2);
950 : }
951 : else
952 : {
953 63 : node0 = node_list.back();
954 81 : node1 = mesh.add_point(Point((i)*dx_tri, dy_tri / 2., 0.));
955 81 : node2 = mesh.add_point(Point((i)*dx_tri, 0., 0.));
956 63 : node_list.push_back(node1);
957 63 : node_list.push_back(node2);
958 : }
959 :
960 385 : new_elem = std::make_unique<C0Polygon>(3);
961 : // Switch to Tri3 when exodus default output supports element type mixes
962 385 : new_elem->set_node(0, node0);
963 385 : new_elem->set_node(1, node1);
964 385 : new_elem->set_node(2, node2);
965 495 : auto * elem = mesh.add_elem(std::move(new_elem));
966 :
967 : // Set boundaries
968 385 : if (i == 0)
969 63 : boundary_info.add_side(elem, 0, 3); // left
970 322 : else if (i == nx)
971 63 : boundary_info.add_side(elem, 1, 1); // right
972 385 : boundary_info.add_side(elem, 2, 0); // bottom
973 : }
974 : // Start with the second node to build hexagons
975 18 : unsigned int running_index = 1;
976 :
977 : // Build layers of hexagons
978 18 : const auto hex_side =
979 63 : (Real(1) - (ny == 1 ?
980 : dy_tri :
981 63 : (Real(1) + (ny - 1) / 2.) * dy_tri)) / ny;
982 385 : for (const auto j : make_range(ny))
983 : {
984 2205 : for (const auto i : make_range(nx + (j % 2)))
985 : {
986 1883 : if ((j % 2 == 0) || ((i > 0) && (i < nx)))
987 : {
988 : Node *n0, *n1, *n2, *n3, *n4, *n5;
989 1589 : n0 = node_list[running_index++];
990 1589 : n1 = node_list[running_index++];
991 1589 : n2 = node_list[running_index];
992 :
993 1589 : if (i == 0)
994 : {
995 225 : n3 = mesh.add_point(Point(*n0) + RealVectorValue(0, hex_side, 0));
996 175 : node_list.push_back(n3);
997 : }
998 : else
999 1414 : n3 = node_list.back();
1000 :
1001 2043 : n4 = mesh.add_point(Point(*n1) + RealVectorValue(0, hex_side + dy_tri, 0));
1002 2043 : n5 = mesh.add_point(Point(*n2) + RealVectorValue(0, hex_side, 0));
1003 1589 : node_list.push_back(n4);
1004 1589 : node_list.push_back(n5);
1005 :
1006 1589 : new_elem = std::make_unique<libMesh::C0Polygon>(6);
1007 1589 : new_elem->set_node(0, n0);
1008 1589 : new_elem->set_node(1, n1);
1009 1589 : new_elem->set_node(2, n2);
1010 1589 : new_elem->set_node(3, n5);
1011 1589 : new_elem->set_node(4, n4);
1012 1589 : new_elem->set_node(5, n3);
1013 2043 : auto * elem = mesh.add_elem(std::move(new_elem));
1014 :
1015 : // Set boundaries
1016 1589 : if (i == 0)
1017 175 : boundary_info.add_side(elem, 5, 3); // left
1018 1414 : else if (i == nx)
1019 908 : boundary_info.add_side(elem, 2, 1); // right
1020 681 : }
1021 : // The hexagons are offset, so we build on a quad on each external side to fill
1022 294 : else if (i == 0 || i == nx)
1023 : {
1024 : Node *n0, *n1, *n2, *n3;
1025 294 : n0 = node_list[running_index++];
1026 294 : n1 = node_list[running_index];
1027 :
1028 294 : if (i == 0)
1029 : {
1030 189 : n2 = mesh.add_point(Point(*n0) + RealVectorValue(0, hex_side + dy_tri, 0));
1031 147 : node_list.push_back(n2);
1032 189 : n3 = mesh.add_point(Point(*n1) + RealVectorValue(0, hex_side, 0));
1033 : }
1034 : else
1035 : {
1036 147 : n2 = node_list.back();
1037 189 : n3 = mesh.add_point(Point(*n1) + RealVectorValue(0, hex_side + dy_tri, 0));
1038 : }
1039 294 : node_list.push_back(n3);
1040 :
1041 294 : new_elem = std::make_unique<C0Polygon>(4);
1042 : // Switch to Quad4 when exodus default output supports element type mixes
1043 294 : new_elem->set_node(0, n0);
1044 294 : new_elem->set_node(1, n1);
1045 294 : new_elem->set_node(3, n2);
1046 294 : new_elem->set_node(2, n3);
1047 378 : auto * elem = mesh.add_elem(std::move(new_elem));
1048 :
1049 : // Set boundaries
1050 294 : if (i == 0)
1051 147 : boundary_info.add_side(elem, 3, 3); // left
1052 147 : else if (i == nx)
1053 231 : boundary_info.add_side(elem, 1, 1); // right
1054 : }
1055 : else
1056 0 : libmesh_assert(false);
1057 : }
1058 : // Increment once to switch to next 'row' of nodes
1059 322 : running_index++;
1060 :
1061 : // Skip lower right corner node
1062 322 : if (j == 0)
1063 63 : running_index++;
1064 : }
1065 :
1066 : // Build a final layer of triangles
1067 63 : const bool ny_odd = (ny % 2 == 1);
1068 413 : for (const auto i : make_range(nx + ny_odd))
1069 : {
1070 : // Use existing nodes, except at the corners
1071 : Node *node0, *node1, *node2;
1072 350 : if (i == 0 && ny_odd)
1073 : {
1074 36 : node0 = mesh.add_point(Point(0., 1., 0.));
1075 28 : node1 = node_list[running_index++];
1076 36 : node2 = node_list[running_index];
1077 : }
1078 322 : else if (i < nx)
1079 : {
1080 294 : node0 = node_list[running_index++];
1081 294 : node1 = node_list[running_index++];
1082 378 : node2 = node_list[running_index];
1083 : }
1084 : // This case only reached if ny is odd and we are using a triangle in top right corner
1085 : else
1086 : {
1087 28 : node0 = node_list[running_index++];
1088 28 : node1 = node_list[running_index];
1089 36 : node2 = mesh.add_point(Point(1., 1., 0.));
1090 : }
1091 :
1092 350 : new_elem = std::make_unique<C0Polygon>(3);
1093 : // Switch to Tri3 when exodus default output supports element type mixes
1094 350 : new_elem->set_node(0, node0);
1095 350 : new_elem->set_node(1, node1);
1096 350 : new_elem->set_node(2, node2);
1097 450 : auto * elem = mesh.add_elem(std::move(new_elem));
1098 :
1099 : // Set boundaries
1100 350 : if (i == 0)
1101 63 : boundary_info.add_side(elem, 0, 3); // left
1102 287 : else if (i == nx)
1103 28 : boundary_info.add_side(elem, 1, 1); // right
1104 350 : boundary_info.add_side(elem, 2, 2); // top
1105 :
1106 : }
1107 18 : break;
1108 27 : }
1109 :
1110 :
1111 0 : default:
1112 0 : libmesh_error_msg("ERROR: Unrecognized 2D element type == " << Utility::enum_to_string(type));
1113 : }
1114 :
1115 :
1116 :
1117 :
1118 : // Scale the nodal positions
1119 11046 : if (gauss_lobatto_grid)
1120 : {
1121 : GaussLobattoRedistributionFunction func(nx, xmin, xmax,
1122 0 : ny, ymin, ymax);
1123 0 : MeshTools::Modification::redistribute(mesh, func);
1124 0 : }
1125 : else // !gauss_lobatto_grid
1126 : {
1127 723318 : for (Node * node : mesh.node_ptr_range())
1128 : {
1129 704390 : (*node)(0) = ((*node)(0))*(xmax-xmin) + xmin;
1130 704390 : (*node)(1) = ((*node)(1))*(ymax-ymin) + ymin;
1131 4718 : }
1132 : }
1133 :
1134 : // Add sideset names to boundary info
1135 11046 : boundary_info.sideset_name(0) = "bottom";
1136 11046 : boundary_info.sideset_name(1) = "right";
1137 11046 : boundary_info.sideset_name(2) = "top";
1138 11046 : boundary_info.sideset_name(3) = "left";
1139 :
1140 : // Add nodeset names to boundary info
1141 11046 : boundary_info.nodeset_name(0) = "bottom";
1142 11046 : boundary_info.nodeset_name(1) = "right";
1143 11046 : boundary_info.nodeset_name(2) = "top";
1144 11046 : boundary_info.nodeset_name(3) = "left";
1145 :
1146 3164 : break;
1147 : }
1148 :
1149 :
1150 :
1151 :
1152 :
1153 :
1154 :
1155 :
1156 :
1157 :
1158 :
1159 : //---------------------------------------------------------------------
1160 : // Build a 3D mesh using hexes, tets, prisms, or pyramids.
1161 12868 : case 3:
1162 : {
1163 3686 : libmesh_assert_not_equal_to (nx, 0);
1164 3686 : libmesh_assert_not_equal_to (ny, 0);
1165 3686 : libmesh_assert_not_equal_to (nz, 0);
1166 3686 : libmesh_assert_less (xmin, xmax);
1167 3686 : libmesh_assert_less (ymin, ymax);
1168 3686 : libmesh_assert_less (zmin, zmax);
1169 :
1170 :
1171 : // Reserve elements. Meshes with prismatic elements require
1172 : // twice as many elements.
1173 : switch (type)
1174 : {
1175 8105 : case INVALID_ELEM:
1176 : case HEX8:
1177 : case HEX20:
1178 : case HEX27:
1179 : case C0POLYHEDRON:
1180 : case TET4: // TET4's are created from an initial HEX27 discretization
1181 : case TET10: // TET10's are created from an initial HEX27 discretization
1182 : case TET14: // TET14's are created from an initial HEX27 discretization
1183 : case PYRAMID5: // PYRAMIDs are created from an initial HEX27 discretization
1184 : case PYRAMID13:
1185 : case PYRAMID14:
1186 : case PYRAMID18:
1187 : {
1188 8105 : mesh.reserve_elem(nx*ny*nz);
1189 5781 : break;
1190 : }
1191 :
1192 4763 : case PRISM6:
1193 : case PRISM15:
1194 : case PRISM18:
1195 : case PRISM20:
1196 : case PRISM21:
1197 : {
1198 4763 : mesh.reserve_elem(2*nx*ny*nz);
1199 3401 : break;
1200 : }
1201 :
1202 0 : default:
1203 0 : libmesh_error_msg("ERROR: Unrecognized 3D element type == " << Utility::enum_to_string(type));
1204 : }
1205 :
1206 :
1207 :
1208 :
1209 :
1210 : // Reserve nodes. Quadratic elements need twice as many nodes as linear elements.
1211 : switch (type)
1212 : {
1213 2227 : case INVALID_ELEM:
1214 : case HEX8:
1215 : case PRISM6:
1216 : case C0POLYHEDRON:
1217 : {
1218 : const dof_id_type grid_nodes =
1219 2227 : cast_int<dof_id_type>((nx+1)*(ny+1)*(nz+1));
1220 :
1221 : // Reserve one interior node per polyhedron for the robust
1222 : // fallback tetrahedralization used when the preferred
1223 : // tetrahedralization cannot be constructed.
1224 : const dof_id_type mid_polyhedron_nodes =
1225 2227 : (type == C0POLYHEDRON) ?
1226 35 : cast_int<dof_id_type>(nx*ny*nz) : 0;
1227 :
1228 2227 : mesh.reserve_nodes(grid_nodes + mid_polyhedron_nodes);
1229 1581 : break;
1230 : }
1231 :
1232 7214 : case HEX20:
1233 : case HEX27:
1234 : case TET4: // TET4's are created from an initial HEX27 discretization
1235 : case TET10: // TET10's are created from an initial HEX27 discretization
1236 : case PYRAMID5: // PYRAMIDs are created from an initial HEX27 discretization
1237 : case PYRAMID13:
1238 : case PYRAMID14:
1239 : case PYRAMID18:
1240 : case PRISM15:
1241 : case PRISM18:
1242 : {
1243 : // FYI: The resulting TET4 mesh will have exactly
1244 : // 5*(nx*ny*nz) + 2*(nx*ny + nx*nz + ny*nz) + (nx+ny+nz) + 1
1245 : // nodes once the additional mid-edge nodes for the HEX27 discretization
1246 : // have been deleted.
1247 7214 : mesh.reserve_nodes( (2*nx+1)*(2*ny+1)*(2*nz+1) );
1248 5152 : break;
1249 : }
1250 :
1251 1194 : case TET14:
1252 : {
1253 2894 : mesh.reserve_nodes( (2*nx+1)*(2*ny+1)*(2*nz+1) +
1254 1194 : 24*nx*ny*nz +
1255 1194 : 4*(nx*ny + ny*nz + nx*nz) );
1256 854 : break;
1257 : }
1258 :
1259 1015 : case PRISM20:
1260 : {
1261 1885 : mesh.reserve_nodes( (2*nx+1)*(2*ny+1)*(2*nz+1) +
1262 1015 : 2*nx*ny*(nz+1) );
1263 725 : break;
1264 : }
1265 :
1266 1218 : case PRISM21:
1267 : {
1268 2262 : mesh.reserve_nodes( (2*nx+1)*(2*ny+1)*(2*nz+1) +
1269 1218 : 2*nx*ny*(2*nz+1) );
1270 870 : break;
1271 : }
1272 :
1273 0 : default:
1274 0 : libmesh_error_msg("ERROR: Unrecognized 3D element type == " << Utility::enum_to_string(type));
1275 : }
1276 :
1277 :
1278 :
1279 :
1280 : // Build the nodes.
1281 3686 : unsigned int node_id = 0;
1282 : switch (type)
1283 : {
1284 646 : case INVALID_ELEM:
1285 : case HEX8:
1286 : case PRISM6:
1287 : case C0POLYHEDRON:
1288 : {
1289 8182 : for (unsigned int k=0; k<=nz; k++)
1290 26884 : for (unsigned int j=0; j<=ny; j++)
1291 182746 : for (unsigned int i=0; i<=nx; i++)
1292 : {
1293 : const Node * const node =
1294 228598 : mesh.add_point(Point(static_cast<Real>(i) / static_cast<Real>(nx),
1295 161817 : static_cast<Real>(j) / static_cast<Real>(ny),
1296 161817 : static_cast<Real>(k) / static_cast<Real>(nz)),
1297 142080 : node_id++);
1298 161817 : if (k == 0)
1299 29763 : boundary_info.add_node(node, 0);
1300 161817 : if (k == nz)
1301 29763 : boundary_info.add_node(node, 5);
1302 161817 : if (j == 0)
1303 25143 : boundary_info.add_node(node, 1);
1304 161817 : if (j == ny)
1305 25143 : boundary_info.add_node(node, 3);
1306 161817 : if (i == 0)
1307 20929 : boundary_info.add_node(node, 4);
1308 161817 : if (i == nx)
1309 20929 : boundary_info.add_node(node, 2);
1310 : }
1311 :
1312 646 : break;
1313 : }
1314 :
1315 3040 : case HEX20:
1316 : case HEX27:
1317 : case TET4: // TET4's are created from an initial HEX27 discretization
1318 : case TET10: // TET10's are created from an initial HEX27 discretization
1319 : case TET14: // TET14's are created from an initial HEX27 discretization
1320 : case PYRAMID5: // PYRAMIDs are created from an initial HEX27 discretization
1321 : case PYRAMID13:
1322 : case PYRAMID14:
1323 : case PYRAMID18:
1324 : case PRISM15:
1325 : case PRISM18:
1326 : case PRISM20:
1327 : case PRISM21:
1328 : {
1329 51188 : for (unsigned int k=0; k<=(2*nz); k++)
1330 224050 : for (unsigned int j=0; j<=(2*ny); j++)
1331 1552768 : for (unsigned int i=0; i<=(2*nx); i++)
1332 : {
1333 : const Node * const node =
1334 1992466 : mesh.add_point(Point(static_cast<Real>(i) / static_cast<Real>(2 * nx),
1335 1369265 : static_cast<Real>(j) / static_cast<Real>(2 * ny),
1336 1369265 : static_cast<Real>(k) / static_cast<Real>(2 * nz)),
1337 1183274 : node_id++);
1338 1369265 : if (k == 0)
1339 194827 : boundary_info.add_node(node, 0);
1340 1369265 : if (k == 2*nz)
1341 194827 : boundary_info.add_node(node, 5);
1342 1369265 : if (j == 0)
1343 189357 : boundary_info.add_node(node, 1);
1344 1369265 : if (j == 2*ny)
1345 189357 : boundary_info.add_node(node, 3);
1346 1369265 : if (i == 0)
1347 183503 : boundary_info.add_node(node, 4);
1348 1369265 : if (i == 2*nx)
1349 183503 : boundary_info.add_node(node, 2);
1350 : }
1351 :
1352 10641 : if (type == PRISM20 ||
1353 : type == PRISM21)
1354 : {
1355 2233 : const unsigned int kmax = (type == PRISM20) ? nz : 2*nz;
1356 9009 : for (unsigned int k=0; k<=kmax; k++)
1357 16464 : for (unsigned int j=0; j<ny; j++)
1358 25200 : for (unsigned int i=0; i<nx; i++)
1359 : {
1360 : const Node * const node1 =
1361 22160 : mesh.add_point(Point((static_cast<Real>(i)+1/Real(3)) / static_cast<Real>(nx),
1362 15512 : (static_cast<Real>(j)+1/Real(3)) / static_cast<Real>(ny),
1363 15512 : static_cast<Real>(k) / static_cast<Real>(kmax)),
1364 13296 : node_id++);
1365 15512 : if (k == 0)
1366 4417 : boundary_info.add_node(node1, 0);
1367 15512 : if (k == kmax)
1368 4417 : boundary_info.add_node(node1, 5);
1369 :
1370 : const Node * const node2 =
1371 22160 : mesh.add_point(Point((static_cast<Real>(i)+2/Real(3)) / static_cast<Real>(nx),
1372 15512 : (static_cast<Real>(j)+2/Real(3)) / static_cast<Real>(ny),
1373 4432 : static_cast<Real>(k) / static_cast<Real>(kmax)),
1374 13296 : node_id++);
1375 15512 : if (k == 0)
1376 4417 : boundary_info.add_node(node2, 0);
1377 15512 : if (k == kmax)
1378 4417 : boundary_info.add_node(node2, 5);
1379 : }
1380 : }
1381 :
1382 3040 : break;
1383 : }
1384 :
1385 :
1386 0 : default:
1387 0 : libmesh_error_msg("ERROR: Unrecognized 3D element type == " << Utility::enum_to_string(type));
1388 : }
1389 :
1390 :
1391 :
1392 :
1393 : // Build the elements.
1394 3686 : unsigned int elem_id = 0;
1395 : switch (type)
1396 : {
1397 410 : case INVALID_ELEM:
1398 : case HEX8:
1399 : {
1400 4016 : for (unsigned int k=0; k<nz; k++)
1401 11936 : for (unsigned int j=0; j<ny; j++)
1402 108561 : for (unsigned int i=0; i<nx; i++)
1403 : {
1404 99240 : Elem * elem = mesh.add_elem(Elem::build_with_id(HEX8, elem_id++));
1405 99240 : elem->set_node(0, mesh.node_ptr(idx(type,nx,ny,i,j,k) ));
1406 99240 : elem->set_node(1, mesh.node_ptr(idx(type,nx,ny,i+1,j,k) ));
1407 99240 : elem->set_node(2, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) ));
1408 99240 : elem->set_node(3, mesh.node_ptr(idx(type,nx,ny,i,j+1,k) ));
1409 99240 : elem->set_node(4, mesh.node_ptr(idx(type,nx,ny,i,j,k+1) ));
1410 99240 : elem->set_node(5, mesh.node_ptr(idx(type,nx,ny,i+1,j,k+1) ));
1411 99240 : elem->set_node(6, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1)));
1412 99240 : elem->set_node(7, mesh.node_ptr(idx(type,nx,ny,i,j+1,k+1) ));
1413 :
1414 99240 : if (k == 0)
1415 17168 : boundary_info.add_side(elem, 0, 0);
1416 :
1417 99240 : if (k == (nz-1))
1418 17168 : boundary_info.add_side(elem, 5, 5);
1419 :
1420 99240 : if (j == 0)
1421 12688 : boundary_info.add_side(elem, 1, 1);
1422 :
1423 99240 : if (j == (ny-1))
1424 12688 : boundary_info.add_side(elem, 3, 3);
1425 :
1426 99240 : if (i == 0)
1427 9321 : boundary_info.add_side(elem, 4, 4);
1428 :
1429 99240 : if (i == (nx-1))
1430 9321 : boundary_info.add_side(elem, 2, 2);
1431 : }
1432 410 : break;
1433 : }
1434 :
1435 :
1436 35 : case C0POLYHEDRON:
1437 : {
1438 35 : const std::array<std::array<unsigned int, 4>, 6> side_nodes =
1439 : {{{0, 1, 2, 3}, // min z
1440 : {0, 1, 5, 4}, // min y
1441 : {2, 6, 5, 1}, // max x
1442 : {2, 3, 7, 6}, // max y
1443 : {0, 4, 7, 3}, // min x
1444 : {5, 6, 7, 4}}}; // max z
1445 :
1446 105 : for (unsigned int k=0; k<nz; k++)
1447 210 : for (unsigned int j=0; j<ny; j++)
1448 420 : for (unsigned int i=0; i<nx; i++)
1449 : {
1450 : std::array<Node *, 8> elem_nodes =
1451 280 : {{mesh.node_ptr(idx(type,nx,ny,i,j,k) ),
1452 280 : mesh.node_ptr(idx(type,nx,ny,i+1,j,k) ),
1453 280 : mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) ),
1454 280 : mesh.node_ptr(idx(type,nx,ny,i,j+1,k) ),
1455 280 : mesh.node_ptr(idx(type,nx,ny,i,j,k+1) ),
1456 280 : mesh.node_ptr(idx(type,nx,ny,i+1,j,k+1) ),
1457 280 : mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1)),
1458 1960 : mesh.node_ptr(idx(type,nx,ny,i,j+1,k+1) )}};
1459 :
1460 440 : std::vector<std::shared_ptr<Polygon>> sides(side_nodes.size());
1461 1960 : for (auto s : index_range(side_nodes))
1462 : {
1463 2160 : sides[s] = std::make_shared<C0Polygon>(side_nodes[s].size());
1464 8400 : for (auto n : index_range(side_nodes[s]))
1465 10560 : sides[s]->set_node(n, elem_nodes[side_nodes[s][n]]);
1466 : }
1467 :
1468 280 : std::unique_ptr<Node> mid_elem_node;
1469 : std::unique_ptr<Elem> new_elem =
1470 360 : std::make_unique<C0Polyhedron>(sides, mid_elem_node);
1471 280 : if (mid_elem_node)
1472 0 : mesh.add_node(std::move(mid_elem_node));
1473 :
1474 280 : new_elem->set_id() = elem_id++;
1475 360 : Elem * elem = mesh.add_elem(std::move(new_elem));
1476 :
1477 280 : if (k == 0)
1478 140 : boundary_info.add_side(elem, 0, 0);
1479 :
1480 280 : if (k == (nz-1))
1481 140 : boundary_info.add_side(elem, 5, 5);
1482 :
1483 280 : if (j == 0)
1484 140 : boundary_info.add_side(elem, 1, 1);
1485 :
1486 280 : if (j == (ny-1))
1487 140 : boundary_info.add_side(elem, 3, 3);
1488 :
1489 280 : if (i == 0)
1490 140 : boundary_info.add_side(elem, 4, 4);
1491 :
1492 280 : if (i == (nx-1))
1493 140 : boundary_info.add_side(elem, 2, 2);
1494 120 : }
1495 10 : break;
1496 : }
1497 :
1498 :
1499 :
1500 :
1501 226 : case PRISM6:
1502 : {
1503 1834 : for (unsigned int k=0; k<nz; k++)
1504 2688 : for (unsigned int j=0; j<ny; j++)
1505 4872 : for (unsigned int i=0; i<nx; i++)
1506 : {
1507 : // First Prism
1508 3227 : Elem * elem = mesh.add_elem(Elem::build_with_id(PRISM6, elem_id++));
1509 3227 : elem->set_node(0, mesh.node_ptr(idx(type,nx,ny,i,j,k) ));
1510 3227 : elem->set_node(1, mesh.node_ptr(idx(type,nx,ny,i+1,j,k) ));
1511 3227 : elem->set_node(2, mesh.node_ptr(idx(type,nx,ny,i,j+1,k) ));
1512 3227 : elem->set_node(3, mesh.node_ptr(idx(type,nx,ny,i,j,k+1) ));
1513 3227 : elem->set_node(4, mesh.node_ptr(idx(type,nx,ny,i+1,j,k+1) ));
1514 3227 : elem->set_node(5, mesh.node_ptr(idx(type,nx,ny,i,j+1,k+1) ));
1515 :
1516 : // Add sides for first prism to boundary info object
1517 3227 : if (i==0)
1518 1645 : boundary_info.add_side(elem, 3, 4);
1519 :
1520 3227 : if (j==0)
1521 1645 : boundary_info.add_side(elem, 1, 1);
1522 :
1523 3227 : if (k==0)
1524 1645 : boundary_info.add_side(elem, 0, 0);
1525 :
1526 3227 : if (k == (nz-1))
1527 1645 : boundary_info.add_side(elem, 4, 5);
1528 :
1529 : // Second Prism
1530 3227 : elem = mesh.add_elem(Elem::build_with_id(PRISM6, elem_id++));
1531 3227 : elem->set_node(0, mesh.node_ptr(idx(type,nx,ny,i+1,j,k) ));
1532 3227 : elem->set_node(1, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) ));
1533 3227 : elem->set_node(2, mesh.node_ptr(idx(type,nx,ny,i,j+1,k) ));
1534 3227 : elem->set_node(3, mesh.node_ptr(idx(type,nx,ny,i+1,j,k+1) ));
1535 3227 : elem->set_node(4, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1)));
1536 3227 : elem->set_node(5, mesh.node_ptr(idx(type,nx,ny,i,j+1,k+1) ));
1537 :
1538 : // Add sides for second prism to boundary info object
1539 3227 : if (i == (nx-1))
1540 1645 : boundary_info.add_side(elem, 1, 2);
1541 :
1542 3227 : if (j == (ny-1))
1543 1645 : boundary_info.add_side(elem, 2, 3);
1544 :
1545 3227 : if (k==0)
1546 1645 : boundary_info.add_side(elem, 0, 0);
1547 :
1548 3227 : if (k == (nz-1))
1549 1645 : boundary_info.add_side(elem, 4, 5);
1550 : }
1551 226 : break;
1552 : }
1553 :
1554 :
1555 :
1556 :
1557 :
1558 :
1559 1904 : case HEX20:
1560 : case HEX27:
1561 : case TET4: // TET4's are created from an initial HEX27 discretization
1562 : case TET10: // TET10's are created from an initial HEX27 discretization
1563 : case TET14: // TET14's are created from an initial HEX27 discretization
1564 : case PYRAMID5: // PYRAMIDs are created from an initial HEX27 discretization
1565 : case PYRAMID13:
1566 : case PYRAMID14:
1567 : case PYRAMID18:
1568 : {
1569 16508 : for (unsigned int k=0; k<(2*nz); k += 2)
1570 30645 : for (unsigned int j=0; j<(2*ny); j += 2)
1571 122742 : for (unsigned int i=0; i<(2*nx); i += 2)
1572 : {
1573 101936 : ElemType build_type = (type == HEX20) ? HEX20 : HEX27;
1574 101936 : Elem * elem = mesh.add_elem(Elem::build_with_id(build_type, elem_id++));
1575 :
1576 101936 : elem->set_node(0, mesh.node_ptr(idx(type,nx,ny,i, j, k) ));
1577 101936 : elem->set_node(1, mesh.node_ptr(idx(type,nx,ny,i+2,j, k) ));
1578 101936 : elem->set_node(2, mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k) ));
1579 101936 : elem->set_node(3, mesh.node_ptr(idx(type,nx,ny,i, j+2,k) ));
1580 101936 : elem->set_node(4, mesh.node_ptr(idx(type,nx,ny,i, j, k+2)));
1581 101936 : elem->set_node(5, mesh.node_ptr(idx(type,nx,ny,i+2,j, k+2)));
1582 101936 : elem->set_node(6, mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k+2)));
1583 101936 : elem->set_node(7, mesh.node_ptr(idx(type,nx,ny,i, j+2,k+2)));
1584 101936 : elem->set_node(8, mesh.node_ptr(idx(type,nx,ny,i+1,j, k) ));
1585 101936 : elem->set_node(9, mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k) ));
1586 101936 : elem->set_node(10, mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k) ));
1587 101936 : elem->set_node(11, mesh.node_ptr(idx(type,nx,ny,i, j+1,k) ));
1588 101936 : elem->set_node(12, mesh.node_ptr(idx(type,nx,ny,i, j, k+1)));
1589 101936 : elem->set_node(13, mesh.node_ptr(idx(type,nx,ny,i+2,j, k+1)));
1590 101936 : elem->set_node(14, mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k+1)));
1591 101936 : elem->set_node(15, mesh.node_ptr(idx(type,nx,ny,i, j+2,k+1)));
1592 101936 : elem->set_node(16, mesh.node_ptr(idx(type,nx,ny,i+1,j, k+2)));
1593 101936 : elem->set_node(17, mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k+2)));
1594 101936 : elem->set_node(18, mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k+2)));
1595 101936 : elem->set_node(19, mesh.node_ptr(idx(type,nx,ny,i, j+1,k+2)));
1596 :
1597 101936 : if ((type == HEX27) || (type == TET4) || (type == TET10) || (type == TET14) ||
1598 5426 : (type == PYRAMID5) || (type == PYRAMID13) || (type == PYRAMID14) ||
1599 1870 : (type == PYRAMID18))
1600 : {
1601 98450 : elem->set_node(20, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) ));
1602 98450 : elem->set_node(21, mesh.node_ptr(idx(type,nx,ny,i+1,j, k+1)));
1603 98450 : elem->set_node(22, mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k+1)));
1604 98450 : elem->set_node(23, mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k+1)));
1605 98450 : elem->set_node(24, mesh.node_ptr(idx(type,nx,ny,i, j+1,k+1)));
1606 98450 : elem->set_node(25, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+2)));
1607 98450 : elem->set_node(26, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1)));
1608 : }
1609 :
1610 101936 : if (k == 0)
1611 23346 : boundary_info.add_side(elem, 0, 0);
1612 :
1613 101936 : if (k == 2*(nz-1))
1614 23346 : boundary_info.add_side(elem, 5, 5);
1615 :
1616 101936 : if (j == 0)
1617 22047 : boundary_info.add_side(elem, 1, 1);
1618 :
1619 101936 : if (j == 2*(ny-1))
1620 22047 : boundary_info.add_side(elem, 3, 3);
1621 :
1622 101936 : if (i == 0)
1623 20806 : boundary_info.add_side(elem, 4, 4);
1624 :
1625 101936 : if (i == 2*(nx-1))
1626 20806 : boundary_info.add_side(elem, 2, 2);
1627 : }
1628 1904 : break;
1629 : }
1630 :
1631 :
1632 :
1633 :
1634 1136 : case PRISM15:
1635 : case PRISM18:
1636 : case PRISM20:
1637 : case PRISM21:
1638 : {
1639 9086 : for (unsigned int k=0; k<(2*nz); k += 2)
1640 12552 : for (unsigned int j=0; j<(2*ny); j += 2)
1641 19704 : for (unsigned int i=0; i<(2*nx); i += 2)
1642 : {
1643 : // First Prism
1644 12266 : Elem * elem = mesh.add_elem(Elem::build_with_id(type, elem_id++));
1645 12266 : elem->set_node(0, mesh.node_ptr(idx(type,nx,ny,i, j, k) ));
1646 12266 : elem->set_node(1, mesh.node_ptr(idx(type,nx,ny,i+2,j, k) ));
1647 12266 : elem->set_node(2, mesh.node_ptr(idx(type,nx,ny,i, j+2,k) ));
1648 12266 : elem->set_node(3, mesh.node_ptr(idx(type,nx,ny,i, j, k+2)));
1649 12266 : elem->set_node(4, mesh.node_ptr(idx(type,nx,ny,i+2,j, k+2)));
1650 12266 : elem->set_node(5, mesh.node_ptr(idx(type,nx,ny,i, j+2,k+2)));
1651 12266 : elem->set_node(6, mesh.node_ptr(idx(type,nx,ny,i+1,j, k) ));
1652 12266 : elem->set_node(7, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) ));
1653 12266 : elem->set_node(8, mesh.node_ptr(idx(type,nx,ny,i, j+1,k) ));
1654 12266 : elem->set_node(9, mesh.node_ptr(idx(type,nx,ny,i, j, k+1)));
1655 12266 : elem->set_node(10, mesh.node_ptr(idx(type,nx,ny,i+2,j, k+1)));
1656 12266 : elem->set_node(11, mesh.node_ptr(idx(type,nx,ny,i, j+2,k+1)));
1657 12266 : elem->set_node(12, mesh.node_ptr(idx(type,nx,ny,i+1,j, k+2)));
1658 12266 : elem->set_node(13, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+2)));
1659 12266 : elem->set_node(14, mesh.node_ptr(idx(type,nx,ny,i, j+1,k+2)));
1660 :
1661 15818 : if (type == PRISM18 ||
1662 10418 : type == PRISM20 ||
1663 : type == PRISM21)
1664 : {
1665 10068 : elem->set_node(15, mesh.node_ptr(idx(type,nx,ny,i+1,j, k+1)));
1666 10068 : elem->set_node(16, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1)));
1667 10068 : elem->set_node(17, mesh.node_ptr(idx(type,nx,ny,i, j+1,k+1)));
1668 : }
1669 :
1670 12266 : if (type == PRISM20)
1671 : {
1672 3563 : const dof_id_type base_idx = (2*nx+1)*(2*ny+1)*(2*nz+1);
1673 3563 : elem->set_node(18, mesh.node_ptr(base_idx+((k/2)*(nx*ny)+j/2*nx+i/2)*2));
1674 3563 : elem->set_node(19, mesh.node_ptr(base_idx+(((k/2)+1)*(nx*ny)+j/2*nx+i/2)*2));
1675 : }
1676 :
1677 12266 : if (type == PRISM21)
1678 : {
1679 3766 : const dof_id_type base_idx = (2*nx+1)*(2*ny+1)*(2*nz+1);
1680 3766 : elem->set_node(18, mesh.node_ptr(base_idx+(k*(nx*ny)+j/2*nx+i/2)*2));
1681 3766 : elem->set_node(19, mesh.node_ptr(base_idx+((k+2)*(nx*ny)+j/2*nx+i/2)*2));
1682 3766 : elem->set_node(20, mesh.node_ptr(base_idx+((k+1)*(nx*ny)+j/2*nx+i/2)*2));
1683 : }
1684 :
1685 : // Add sides for first prism to boundary info object
1686 12266 : if (i==0)
1687 7438 : boundary_info.add_side(elem, 3, 4);
1688 :
1689 12266 : if (j==0)
1690 7438 : boundary_info.add_side(elem, 1, 1);
1691 :
1692 12266 : if (k==0)
1693 7488 : boundary_info.add_side(elem, 0, 0);
1694 :
1695 12266 : if (k == 2*(nz-1))
1696 7488 : boundary_info.add_side(elem, 4, 5);
1697 :
1698 :
1699 : // Second Prism
1700 12266 : elem = mesh.add_elem(Elem::build_with_id(type, elem_id++));
1701 12266 : elem->set_node(0, mesh.node_ptr(idx(type,nx,ny,i+2,j,k) ));
1702 12266 : elem->set_node(1, mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k) ));
1703 12266 : elem->set_node(2, mesh.node_ptr(idx(type,nx,ny,i,j+2,k) ));
1704 12266 : elem->set_node(3, mesh.node_ptr(idx(type,nx,ny,i+2,j,k+2) ));
1705 12266 : elem->set_node(4, mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k+2) ));
1706 12266 : elem->set_node(5, mesh.node_ptr(idx(type,nx,ny,i,j+2,k+2) ));
1707 12266 : elem->set_node(6, mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k) ));
1708 12266 : elem->set_node(7, mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k) ));
1709 12266 : elem->set_node(8, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) ));
1710 12266 : elem->set_node(9, mesh.node_ptr(idx(type,nx,ny,i+2,j,k+1) ));
1711 12266 : elem->set_node(10, mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k+1)));
1712 12266 : elem->set_node(11, mesh.node_ptr(idx(type,nx,ny,i,j+2,k+1) ));
1713 12266 : elem->set_node(12, mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k+2)));
1714 12266 : elem->set_node(13, mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k+2)));
1715 12266 : elem->set_node(14, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+2)));
1716 :
1717 12266 : if (type == PRISM18 ||
1718 5964 : type == PRISM20 ||
1719 : type == PRISM21)
1720 : {
1721 10068 : elem->set_node(15, mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k+1)));
1722 10068 : elem->set_node(16, mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k+1)));
1723 10068 : elem->set_node(17, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1)));
1724 : }
1725 :
1726 12266 : if (type == PRISM20)
1727 : {
1728 3563 : const dof_id_type base_idx = (2*nx+1)*(2*ny+1)*(2*nz+1);
1729 3563 : elem->set_node(18, mesh.node_ptr(base_idx+((k/2)*(nx*ny)+j/2*nx+i/2)*2+1));
1730 3563 : elem->set_node(19, mesh.node_ptr(base_idx+(((k/2)+1)*(nx*ny)+j/2*nx+i/2)*2+1));
1731 : }
1732 :
1733 12266 : if (type == PRISM21)
1734 : {
1735 3766 : const dof_id_type base_idx = (2*nx+1)*(2*ny+1)*(2*nz+1);
1736 3766 : elem->set_node(18, mesh.node_ptr(base_idx+(k*(nx*ny)+j/2*nx+i/2)*2+1));
1737 3766 : elem->set_node(19, mesh.node_ptr(base_idx+((k+2)*(nx*ny)+j/2*nx+i/2)*2+1));
1738 3766 : elem->set_node(20, mesh.node_ptr(base_idx+((k+1)*(nx*ny)+j/2*nx+i/2)*2+1));
1739 : }
1740 :
1741 : // Add sides for second prism to boundary info object
1742 12266 : if (i == 2*(nx-1))
1743 7438 : boundary_info.add_side(elem, 1, 2);
1744 :
1745 12266 : if (j == 2*(ny-1))
1746 7438 : boundary_info.add_side(elem, 2, 3);
1747 :
1748 12266 : if (k==0)
1749 7488 : boundary_info.add_side(elem, 0, 0);
1750 :
1751 12266 : if (k == 2*(nz-1))
1752 7488 : boundary_info.add_side(elem, 4, 5);
1753 :
1754 : }
1755 1136 : break;
1756 : }
1757 :
1758 :
1759 :
1760 :
1761 :
1762 0 : default:
1763 0 : libmesh_error_msg("ERROR: Unrecognized 3D element type == " << Utility::enum_to_string(type));
1764 : }
1765 :
1766 :
1767 :
1768 :
1769 : //.......................................
1770 : // Scale the nodal positions
1771 12868 : if (gauss_lobatto_grid)
1772 : {
1773 : GaussLobattoRedistributionFunction func(nx, xmin, xmax,
1774 : ny, ymin, ymax,
1775 0 : nz, zmin, zmax);
1776 0 : MeshTools::Modification::redistribute(mesh, func);
1777 0 : }
1778 : else // !gauss_lobatto_grid
1779 : {
1780 1574974 : for (unsigned int p=0; p<mesh.n_nodes(); p++)
1781 : {
1782 1562106 : mesh.node_ref(p)(0) = (mesh.node_ref(p)(0))*(xmax-xmin) + xmin;
1783 1562106 : mesh.node_ref(p)(1) = (mesh.node_ref(p)(1))*(ymax-ymin) + ymin;
1784 1562106 : mesh.node_ref(p)(2) = (mesh.node_ref(p)(2))*(zmax-zmin) + zmin;
1785 : }
1786 : }
1787 :
1788 :
1789 :
1790 : // Additional work for tets and pyramids: we take the existing
1791 : // HEX27 discretization and split each element into 24
1792 : // sub-tets or 6 sub-pyramids.
1793 : //
1794 : // 24 isn't the minimum-possible number of tets, but it
1795 : // obviates any concerns about the edge orientations between
1796 : // the various elements.
1797 16554 : if ((type == TET4) ||
1798 12368 : (type == TET10) ||
1799 12028 : (type == TET14) ||
1800 9848 : (type == PYRAMID5) ||
1801 2692 : (type == PYRAMID13) ||
1802 12003 : (type == PYRAMID14) ||
1803 6695 : (type == PYRAMID18))
1804 : {
1805 : // Temporary storage for new elements. (24 tets per hex, 6 pyramids)
1806 3420 : std::vector<std::unique_ptr<Elem>> new_elements;
1807 :
1808 : // For avoiding extraneous construction of element sides
1809 3992 : std::unique_ptr<Elem> side;
1810 :
1811 3992 : if ((type == TET4) || (type == TET10) || (type == TET14))
1812 2942 : new_elements.reserve(24*mesh.n_elem());
1813 : else
1814 1050 : new_elements.reserve(6*mesh.n_elem());
1815 :
1816 : // Create tetrahedra or pyramids
1817 39434 : for (auto & base_hex : mesh.element_ptr_range())
1818 : {
1819 : // Get a pointer to the node located at the HEX27 center
1820 22613 : Node * apex_node = base_hex->node_ptr(26);
1821 :
1822 : // Container to catch ids handed back from BoundaryInfo
1823 12636 : std::vector<boundary_id_type> ids;
1824 :
1825 158291 : for (auto s : base_hex->side_index_range())
1826 : {
1827 : // Get the boundary ID(s) for this side
1828 135678 : boundary_info.boundary_ids(base_hex, s, ids);
1829 :
1830 : // We're creating this Mesh, so there should be 0 or 1 boundary IDs.
1831 37908 : libmesh_assert(ids.size() <= 1);
1832 :
1833 : // A convenient name for the side's ID.
1834 135678 : boundary_id_type b_id = ids.empty() ? BoundaryInfo::invalid_id : ids[0];
1835 :
1836 : // Need to build the full-ordered side!
1837 135678 : base_hex->build_side_ptr(side, s);
1838 :
1839 135678 : if ((type == TET4) || (type == TET10) || (type == TET14))
1840 : {
1841 : // Build 4 sub-tets per side
1842 460200 : for (unsigned int sub_tet=0; sub_tet<4; ++sub_tet)
1843 : {
1844 634560 : new_elements.push_back( Elem::build(TET4) );
1845 101760 : auto & sub_elem = new_elements.back();
1846 571680 : sub_elem->set_node(0, side->node_ptr(sub_tet));
1847 571680 : sub_elem->set_node(1, side->node_ptr(8)); // center of the face
1848 469920 : sub_elem->set_node(2, side->node_ptr(sub_tet==3 ? 0 : sub_tet+1 )); // wrap-around
1849 368160 : sub_elem->set_node(3, apex_node); // apex node always used!
1850 :
1851 : // If the original hex was a boundary hex, add the new sub_tet's side
1852 : // 0 with the same b_id. Note: the tets are all aligned so that their
1853 : // side 0 is on the boundary.
1854 368160 : if (b_id != BoundaryInfo::invalid_id)
1855 206696 : boundary_info.add_side(sub_elem.get(), 0, b_id);
1856 25440 : }
1857 : } // end if ((type == TET4) || (type == TET10) || (type == TET14))
1858 :
1859 : else // type==PYRAMID*
1860 : {
1861 : // Build 1 sub-pyramid per side.
1862 74808 : new_elements.push_back( Elem::build(PYRAMID5) );
1863 12468 : auto & sub_elem = new_elements.back();
1864 :
1865 : // Set the base. Note that since the apex is *inside* the base_hex,
1866 : // and the pyramid uses a counter-clockwise base numbering, we need to
1867 : // reverse the [1] and [3] node indices.
1868 68574 : sub_elem->set_node(0, side->node_ptr(0));
1869 68574 : sub_elem->set_node(1, side->node_ptr(3));
1870 68574 : sub_elem->set_node(2, side->node_ptr(2));
1871 68574 : sub_elem->set_node(3, side->node_ptr(1));
1872 :
1873 : // Set the apex
1874 43638 : sub_elem->set_node(4, apex_node);
1875 :
1876 : // If the original hex was a boundary hex, add the new sub_pyr's side
1877 : // 4 (the square base) with the same b_id.
1878 43638 : if (b_id != BoundaryInfo::invalid_id)
1879 27378 : boundary_info.add_side(sub_elem.get(), 4, b_id);
1880 : } // end else type==PYRAMID*
1881 : }
1882 1712 : }
1883 :
1884 :
1885 : // Delete the original HEX27 elements from the mesh, and the boundary info structure.
1886 39434 : for (auto & elem : mesh.element_ptr_range())
1887 : {
1888 22613 : boundary_info.remove(elem); // Safe even if elem has no boundary info.
1889 22613 : mesh.delete_elem(elem);
1890 1712 : }
1891 :
1892 : // Add the new elements
1893 415790 : for (auto i : index_range(new_elements))
1894 : {
1895 399798 : new_elements[i]->set_id(i);
1896 640254 : mesh.add_elem( std::move(new_elements[i]) );
1897 : }
1898 :
1899 1712 : } // end if (type == TET*,PYRAMID*)
1900 :
1901 :
1902 : // Use all_second_order to convert the TET4's to TET10's or PYRAMID5's to PYRAMID14's
1903 12868 : if ((type == TET10) || (type == PYRAMID14))
1904 1167 : mesh.all_second_order();
1905 :
1906 11701 : else if (type == PYRAMID13)
1907 266 : mesh.all_second_order(/*full_ordered=*/false);
1908 :
1909 11435 : else if ((type == TET14) || (type == PYRAMID18))
1910 1439 : mesh.all_complete_order();
1911 :
1912 :
1913 : // Add sideset names to boundary info (Z axis out of the screen)
1914 12868 : boundary_info.sideset_name(0) = "back";
1915 12868 : boundary_info.sideset_name(1) = "bottom";
1916 12868 : boundary_info.sideset_name(2) = "right";
1917 12868 : boundary_info.sideset_name(3) = "top";
1918 12868 : boundary_info.sideset_name(4) = "left";
1919 12868 : boundary_info.sideset_name(5) = "front";
1920 :
1921 : // Add nodeset names to boundary info
1922 12868 : boundary_info.nodeset_name(0) = "back";
1923 12868 : boundary_info.nodeset_name(1) = "bottom";
1924 12868 : boundary_info.nodeset_name(2) = "right";
1925 12868 : boundary_info.nodeset_name(3) = "top";
1926 12868 : boundary_info.nodeset_name(4) = "left";
1927 12868 : boundary_info.nodeset_name(5) = "front";
1928 :
1929 3686 : break;
1930 : } // end case dim==3
1931 :
1932 0 : default:
1933 0 : libmesh_error_msg("Unknown dimension " << mesh.mesh_dimension());
1934 : }
1935 :
1936 : // Done building the mesh. Now prepare it for use.
1937 27567 : mesh.prepare_for_use ();
1938 27567 : }
1939 :
1940 :
1941 :
1942 98 : void MeshTools::Generation::build_point (UnstructuredMesh & mesh,
1943 : const ElemType type,
1944 : const bool gauss_lobatto_grid)
1945 : {
1946 : // This method only makes sense in 0D!
1947 : // But we now just turn a non-0D mesh into a 0D mesh
1948 : //libmesh_assert_equal_to (mesh.mesh_dimension(), 1);
1949 :
1950 98 : build_cube(mesh,
1951 : 0, 0, 0,
1952 : 0., 0.,
1953 : 0., 0.,
1954 : 0., 0.,
1955 : type,
1956 : gauss_lobatto_grid);
1957 98 : }
1958 :
1959 :
1960 469 : void MeshTools::Generation::build_line (UnstructuredMesh & mesh,
1961 : const unsigned int nx,
1962 : const Real xmin, const Real xmax,
1963 : const ElemType type,
1964 : const bool gauss_lobatto_grid)
1965 : {
1966 : // This method only makes sense in 1D!
1967 : // But we now just turn a non-1D mesh into a 1D mesh
1968 : //libmesh_assert_equal_to (mesh.mesh_dimension(), 1);
1969 :
1970 469 : build_cube(mesh,
1971 : nx, 0, 0,
1972 : xmin, xmax,
1973 : 0., 0.,
1974 : 0., 0.,
1975 : type,
1976 : gauss_lobatto_grid);
1977 469 : }
1978 :
1979 :
1980 :
1981 2045 : void MeshTools::Generation::build_square (UnstructuredMesh & mesh,
1982 : const unsigned int nx,
1983 : const unsigned int ny,
1984 : const Real xmin, const Real xmax,
1985 : const Real ymin, const Real ymax,
1986 : const ElemType type,
1987 : const bool gauss_lobatto_grid)
1988 : {
1989 : // This method only makes sense in 2D!
1990 : // But we now just turn a non-2D mesh into a 2D mesh
1991 : //libmesh_assert_equal_to (mesh.mesh_dimension(), 2);
1992 :
1993 : // Call the build_cube() member to actually do the work for us.
1994 2045 : build_cube (mesh,
1995 : nx, ny, 0,
1996 : xmin, xmax,
1997 : ymin, ymax,
1998 : 0., 0.,
1999 : type,
2000 : gauss_lobatto_grid);
2001 2045 : }
2002 :
2003 :
2004 :
2005 :
2006 :
2007 :
2008 :
2009 :
2010 :
2011 : #ifndef LIBMESH_ENABLE_AMR
2012 : void MeshTools::Generation::build_sphere (UnstructuredMesh &,
2013 : const Real,
2014 : const unsigned int,
2015 : const ElemType,
2016 : const unsigned int,
2017 : const bool)
2018 : {
2019 : libmesh_error_msg("Building a circle/sphere only works with AMR.");
2020 : }
2021 :
2022 : #else
2023 :
2024 233 : void MeshTools::Generation::build_sphere (UnstructuredMesh & mesh,
2025 : const Real rad,
2026 : const unsigned int nr,
2027 : const ElemType type,
2028 : const unsigned int n_smooth,
2029 : const bool flat)
2030 : {
2031 68 : libmesh_assert_greater (rad, 0.);
2032 : //libmesh_assert_greater (nr, 0); // must refine at least once otherwise will end up with a square/cube
2033 :
2034 136 : LOG_SCOPE("build_sphere()", "MeshTools::Generation");
2035 :
2036 : // Clear the mesh and start from scratch, but save the original
2037 : // mesh_dimension, since the original intent of this function was to
2038 : // allow the geometric entity (line, circle, ball, sphere)
2039 : // constructed to be determined by the mesh's dimension.
2040 : unsigned char orig_mesh_dimension =
2041 233 : cast_int<unsigned char>(mesh.mesh_dimension());
2042 233 : mesh.clear();
2043 233 : mesh.set_mesh_dimension(orig_mesh_dimension);
2044 :
2045 : // If mesh.mesh_dimension()==1, it *could* be because the user
2046 : // constructed a Mesh without specifying a dimension (since this is
2047 : // allowed now) and hence it got the default dimension of 1. In
2048 : // this case, we will try to infer the dimension they *really*
2049 : // wanted from the requested ElemType, and if they don't match, go
2050 : // with the ElemType.
2051 233 : if (mesh.mesh_dimension() == 1)
2052 : {
2053 233 : switch (type)
2054 : {
2055 89 : case HEX8:
2056 : case HEX27:
2057 : case TET4:
2058 : case TET10:
2059 : case TET14:
2060 89 : mesh.set_mesh_dimension(3);
2061 26 : break;
2062 116 : case TRI3:
2063 : case TRI6:
2064 : case TRI7:
2065 : case QUAD4:
2066 : case QUADSHELL4:
2067 : case QUAD8:
2068 : case QUADSHELL8:
2069 : case QUAD9:
2070 : case QUADSHELL9:
2071 116 : mesh.set_mesh_dimension(2);
2072 34 : break;
2073 28 : case EDGE2:
2074 : case EDGE3:
2075 : case EDGE4:
2076 28 : mesh.set_mesh_dimension(1);
2077 8 : break;
2078 0 : case INVALID_ELEM:
2079 : // Just keep the existing dimension
2080 0 : break;
2081 0 : default:
2082 0 : libmesh_error_msg("build_sphere(): Please specify a mesh dimension or a valid ElemType (EDGE{2,3,4}, TRI{3,6,7}, QUAD{4,8,9}, HEX{8,27}, TET{4,10,14})");
2083 : }
2084 : }
2085 :
2086 68 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
2087 :
2088 : // Building while distributed is a little more complicated
2089 233 : const bool is_replicated = mesh.is_replicated();
2090 :
2091 : // Sphere is centered at origin by default
2092 68 : const Point cent;
2093 :
2094 466 : const Sphere sphere (cent, rad);
2095 :
2096 233 : switch (mesh.mesh_dimension())
2097 : {
2098 : //-----------------------------------------------------------------
2099 : // Build a line in one dimension
2100 28 : case 1:
2101 : {
2102 28 : build_line (mesh, 3, -rad, rad, type);
2103 :
2104 8 : break;
2105 : }
2106 :
2107 :
2108 :
2109 :
2110 : //-----------------------------------------------------------------
2111 : // Build a circle or hollow sphere in two dimensions
2112 116 : case 2:
2113 : {
2114 : // For DistributedMesh, if we don't specify node IDs the Mesh
2115 : // will try to pick an appropriate (unique) one for us. But
2116 : // since we are adding these nodes on all processors, we want
2117 : // to be sure they have consistent IDs across all processors.
2118 34 : unsigned node_id = 0;
2119 :
2120 116 : if (flat)
2121 : {
2122 34 : const Real sqrt_2 = std::sqrt(2.);
2123 116 : const Real rad_2 = .25*rad;
2124 116 : const Real rad_sqrt_2 = rad/sqrt_2;
2125 :
2126 : // (Temporary) convenient storage for node pointers
2127 150 : std::vector<Node *> nodes(8);
2128 :
2129 : // Point 0
2130 150 : nodes[0] = mesh.add_point (Point(-rad_2,-rad_2, 0.), node_id++);
2131 :
2132 : // Point 1
2133 150 : nodes[1] = mesh.add_point (Point( rad_2,-rad_2, 0.), node_id++);
2134 :
2135 : // Point 2
2136 150 : nodes[2] = mesh.add_point (Point( rad_2, rad_2, 0.), node_id++);
2137 :
2138 : // Point 3
2139 150 : nodes[3] = mesh.add_point (Point(-rad_2, rad_2, 0.), node_id++);
2140 :
2141 : // Point 4
2142 150 : nodes[4] = mesh.add_point (Point(-rad_sqrt_2,-rad_sqrt_2, 0.), node_id++);
2143 :
2144 : // Point 5
2145 150 : nodes[5] = mesh.add_point (Point( rad_sqrt_2,-rad_sqrt_2, 0.), node_id++);
2146 :
2147 : // Point 6
2148 150 : nodes[6] = mesh.add_point (Point( rad_sqrt_2, rad_sqrt_2, 0.), node_id++);
2149 :
2150 : // Point 7
2151 150 : nodes[7] = mesh.add_point (Point(-rad_sqrt_2, rad_sqrt_2, 0.), node_id++);
2152 :
2153 : // Build the elements & set node pointers
2154 :
2155 : // Element 0
2156 : {
2157 116 : Elem * elem0 = mesh.add_elem (Elem::build(QUAD4));
2158 116 : elem0->set_node(0, nodes[0]);
2159 116 : elem0->set_node(1, nodes[1]);
2160 116 : elem0->set_node(2, nodes[2]);
2161 116 : elem0->set_node(3, nodes[3]);
2162 : }
2163 :
2164 : // Element 1
2165 : {
2166 116 : Elem * elem1 = mesh.add_elem (Elem::build(QUAD4));
2167 116 : elem1->set_node(0, nodes[4]);
2168 116 : elem1->set_node(1, nodes[0]);
2169 116 : elem1->set_node(2, nodes[3]);
2170 116 : elem1->set_node(3, nodes[7]);
2171 : }
2172 :
2173 : // Element 2
2174 : {
2175 116 : Elem * elem2 = mesh.add_elem (Elem::build(QUAD4));
2176 116 : elem2->set_node(0, nodes[4]);
2177 116 : elem2->set_node(1, nodes[5]);
2178 116 : elem2->set_node(2, nodes[1]);
2179 116 : elem2->set_node(3, nodes[0]);
2180 : }
2181 :
2182 : // Element 3
2183 : {
2184 116 : Elem * elem3 = mesh.add_elem (Elem::build(QUAD4));
2185 116 : elem3->set_node(0, nodes[1]);
2186 116 : elem3->set_node(1, nodes[5]);
2187 116 : elem3->set_node(2, nodes[6]);
2188 116 : elem3->set_node(3, nodes[2]);
2189 : }
2190 :
2191 : // Element 4
2192 : {
2193 116 : Elem * elem4 = mesh.add_elem (Elem::build(QUAD4));
2194 116 : elem4->set_node(0, nodes[3]);
2195 116 : elem4->set_node(1, nodes[2]);
2196 116 : elem4->set_node(2, nodes[6]);
2197 116 : elem4->set_node(3, nodes[7]);
2198 : }
2199 :
2200 : }
2201 : else
2202 : {
2203 : // Create the 12 vertices of a regular unit icosahedron
2204 0 : Real t = 0.5 * (1 + std::sqrt(5.0));
2205 0 : Real s = rad / std::sqrt(1 + t*t);
2206 0 : t *= s;
2207 :
2208 0 : mesh.add_point (Point(-s, t, 0), node_id++);
2209 0 : mesh.add_point (Point( s, t, 0), node_id++);
2210 0 : mesh.add_point (Point(-s, -t, 0), node_id++);
2211 0 : mesh.add_point (Point( s, -t, 0), node_id++);
2212 :
2213 0 : mesh.add_point (Point( 0, -s, t), node_id++);
2214 0 : mesh.add_point (Point( 0, s, t), node_id++);
2215 0 : mesh.add_point (Point( 0, -s, -t), node_id++);
2216 0 : mesh.add_point (Point( 0, s, -t), node_id++);
2217 :
2218 0 : mesh.add_point (Point( t, 0, -s), node_id++);
2219 0 : mesh.add_point (Point( t, 0, s), node_id++);
2220 0 : mesh.add_point (Point(-t, 0, -s), node_id++);
2221 0 : mesh.add_point (Point(-t, 0, s), node_id++);
2222 :
2223 : // Create the 20 triangles of the icosahedron
2224 : static const unsigned int idx1 [6] = {11, 5, 1, 7, 10, 11};
2225 : static const unsigned int idx2 [6] = {9, 4, 2, 6, 8, 9};
2226 : static const unsigned int idx3 [6] = {1, 5, 11, 10, 7, 1};
2227 :
2228 0 : for (unsigned int i = 0; i < 5; ++i)
2229 : {
2230 : // 5 elems around point 0
2231 0 : Elem * new_elem = mesh.add_elem(Elem::build(TRI3));
2232 0 : new_elem->set_node(0, mesh.node_ptr(0));
2233 0 : new_elem->set_node(1, mesh.node_ptr(idx1[i]));
2234 0 : new_elem->set_node(2, mesh.node_ptr(idx1[i+1]));
2235 :
2236 : // 5 adjacent elems
2237 0 : new_elem = mesh.add_elem(Elem::build(TRI3));
2238 0 : new_elem->set_node(0, mesh.node_ptr(idx3[i]));
2239 0 : new_elem->set_node(1, mesh.node_ptr(idx3[i+1]));
2240 0 : new_elem->set_node(2, mesh.node_ptr(idx2[i]));
2241 :
2242 : // 5 elems around point 3
2243 0 : new_elem = mesh.add_elem(Elem::build(TRI3));
2244 0 : new_elem->set_node(0, mesh.node_ptr(3));
2245 0 : new_elem->set_node(1, mesh.node_ptr(idx2[i]));
2246 0 : new_elem->set_node(2, mesh.node_ptr(idx2[i+1]));
2247 :
2248 : // 5 adjacent elems
2249 0 : new_elem = mesh.add_elem(Elem::build(TRI3));
2250 0 : new_elem->set_node(0, mesh.node_ptr(idx2[i+1]));
2251 0 : new_elem->set_node(1, mesh.node_ptr(idx2[i]));
2252 0 : new_elem->set_node(2, mesh.node_ptr(idx3[i+1]));
2253 : }
2254 : }
2255 :
2256 34 : break;
2257 : } // end case 2
2258 :
2259 :
2260 :
2261 :
2262 :
2263 : //-----------------------------------------------------------------
2264 : // Build a sphere in three dimensions
2265 89 : case 3:
2266 : {
2267 : // (Currently) supported types
2268 89 : if (!((type == HEX8) || (type == HEX27) || (type == TET4) ||
2269 : (type == TET10) || (type == TET14)))
2270 : {
2271 0 : libmesh_error_msg("Error: Only HEX8/27 and TET4/10/14 are currently supported in 3D.");
2272 : }
2273 :
2274 :
2275 : // 3D analog of 2D initial grid:
2276 : const Real
2277 89 : r_small = 0.25*rad, // 0.25 *radius
2278 89 : r_med = (0.125*std::sqrt(2.)+0.5)*rad; // .67677*radius
2279 :
2280 : // (Temporary) convenient storage for node pointers
2281 115 : std::vector<Node *> nodes(16);
2282 :
2283 : // For DistributedMesh, if we don't specify node IDs the Mesh
2284 : // will try to pick an appropriate (unique) one for us. But
2285 : // since we are adding these nodes on all processors, we want
2286 : // to be sure they have consistent IDs across all processors.
2287 26 : unsigned node_id = 0;
2288 :
2289 : // Points 0-7 are the initial HEX8
2290 115 : nodes[0] = mesh.add_point (Point(-r_small,-r_small, -r_small), node_id++);
2291 115 : nodes[1] = mesh.add_point (Point( r_small,-r_small, -r_small), node_id++);
2292 115 : nodes[2] = mesh.add_point (Point( r_small, r_small, -r_small), node_id++);
2293 115 : nodes[3] = mesh.add_point (Point(-r_small, r_small, -r_small), node_id++);
2294 115 : nodes[4] = mesh.add_point (Point(-r_small,-r_small, r_small), node_id++);
2295 115 : nodes[5] = mesh.add_point (Point( r_small,-r_small, r_small), node_id++);
2296 115 : nodes[6] = mesh.add_point (Point( r_small, r_small, r_small), node_id++);
2297 115 : nodes[7] = mesh.add_point (Point(-r_small, r_small, r_small), node_id++);
2298 :
2299 : // Points 8-15 are for the outer hexes, we number them in the same way
2300 115 : nodes[8] = mesh.add_point (Point(-r_med,-r_med, -r_med), node_id++);
2301 115 : nodes[9] = mesh.add_point (Point( r_med,-r_med, -r_med), node_id++);
2302 115 : nodes[10] = mesh.add_point (Point( r_med, r_med, -r_med), node_id++);
2303 115 : nodes[11] = mesh.add_point (Point(-r_med, r_med, -r_med), node_id++);
2304 115 : nodes[12] = mesh.add_point (Point(-r_med,-r_med, r_med), node_id++);
2305 115 : nodes[13] = mesh.add_point (Point( r_med,-r_med, r_med), node_id++);
2306 115 : nodes[14] = mesh.add_point (Point( r_med, r_med, r_med), node_id++);
2307 115 : nodes[15] = mesh.add_point (Point(-r_med, r_med, r_med), node_id++);
2308 :
2309 : // Now create the elements and add them to the mesh
2310 : // Element 0 - center element
2311 : {
2312 89 : Elem * elem0 = mesh.add_elem(Elem::build(HEX8));
2313 89 : elem0->set_node(0, nodes[0]);
2314 89 : elem0->set_node(1, nodes[1]);
2315 89 : elem0->set_node(2, nodes[2]);
2316 89 : elem0->set_node(3, nodes[3]);
2317 89 : elem0->set_node(4, nodes[4]);
2318 89 : elem0->set_node(5, nodes[5]);
2319 89 : elem0->set_node(6, nodes[6]);
2320 89 : elem0->set_node(7, nodes[7]);
2321 : }
2322 :
2323 : // Element 1 - "bottom"
2324 : {
2325 89 : Elem * elem1 = mesh.add_elem(Elem::build(HEX8));
2326 89 : elem1->set_node(0, nodes[8]);
2327 89 : elem1->set_node(1, nodes[9]);
2328 89 : elem1->set_node(2, nodes[10]);
2329 89 : elem1->set_node(3, nodes[11]);
2330 89 : elem1->set_node(4, nodes[0]);
2331 89 : elem1->set_node(5, nodes[1]);
2332 89 : elem1->set_node(6, nodes[2]);
2333 89 : elem1->set_node(7, nodes[3]);
2334 : }
2335 :
2336 : // Element 2 - "front"
2337 : {
2338 89 : Elem * elem2 = mesh.add_elem(Elem::build(HEX8));
2339 89 : elem2->set_node(0, nodes[8]);
2340 89 : elem2->set_node(1, nodes[9]);
2341 89 : elem2->set_node(2, nodes[1]);
2342 89 : elem2->set_node(3, nodes[0]);
2343 89 : elem2->set_node(4, nodes[12]);
2344 89 : elem2->set_node(5, nodes[13]);
2345 89 : elem2->set_node(6, nodes[5]);
2346 89 : elem2->set_node(7, nodes[4]);
2347 : }
2348 :
2349 : // Element 3 - "right"
2350 : {
2351 89 : Elem * elem3 = mesh.add_elem(Elem::build(HEX8));
2352 89 : elem3->set_node(0, nodes[1]);
2353 89 : elem3->set_node(1, nodes[9]);
2354 89 : elem3->set_node(2, nodes[10]);
2355 89 : elem3->set_node(3, nodes[2]);
2356 89 : elem3->set_node(4, nodes[5]);
2357 89 : elem3->set_node(5, nodes[13]);
2358 89 : elem3->set_node(6, nodes[14]);
2359 89 : elem3->set_node(7, nodes[6]);
2360 : }
2361 :
2362 : // Element 4 - "back"
2363 : {
2364 89 : Elem * elem4 = mesh.add_elem(Elem::build(HEX8));
2365 89 : elem4->set_node(0, nodes[3]);
2366 89 : elem4->set_node(1, nodes[2]);
2367 89 : elem4->set_node(2, nodes[10]);
2368 89 : elem4->set_node(3, nodes[11]);
2369 89 : elem4->set_node(4, nodes[7]);
2370 89 : elem4->set_node(5, nodes[6]);
2371 89 : elem4->set_node(6, nodes[14]);
2372 89 : elem4->set_node(7, nodes[15]);
2373 : }
2374 :
2375 : // Element 5 - "left"
2376 : {
2377 89 : Elem * elem5 = mesh.add_elem(Elem::build(HEX8));
2378 89 : elem5->set_node(0, nodes[8]);
2379 89 : elem5->set_node(1, nodes[0]);
2380 89 : elem5->set_node(2, nodes[3]);
2381 89 : elem5->set_node(3, nodes[11]);
2382 89 : elem5->set_node(4, nodes[12]);
2383 89 : elem5->set_node(5, nodes[4]);
2384 89 : elem5->set_node(6, nodes[7]);
2385 89 : elem5->set_node(7, nodes[15]);
2386 : }
2387 :
2388 : // Element 6 - "top"
2389 : {
2390 89 : Elem * elem6 = mesh.add_elem(Elem::build(HEX8));
2391 89 : elem6->set_node(0, nodes[4]);
2392 89 : elem6->set_node(1, nodes[5]);
2393 89 : elem6->set_node(2, nodes[6]);
2394 89 : elem6->set_node(3, nodes[7]);
2395 89 : elem6->set_node(4, nodes[12]);
2396 89 : elem6->set_node(5, nodes[13]);
2397 89 : elem6->set_node(6, nodes[14]);
2398 89 : elem6->set_node(7, nodes[15]);
2399 : }
2400 :
2401 26 : break;
2402 : } // end case 3
2403 :
2404 0 : default:
2405 0 : libmesh_error_msg("Unknown dimension " << mesh.mesh_dimension());
2406 :
2407 :
2408 :
2409 : } // end switch (dim)
2410 :
2411 : // Now we have the beginnings of a sphere.
2412 : // Add some more elements by doing uniform refinements and
2413 : // popping nodes to the boundary.
2414 466 : MeshRefinement mesh_refinement (mesh);
2415 :
2416 : // For avoiding extraneous element side construction
2417 233 : std::unique_ptr<Elem> side;
2418 :
2419 : // Loop over the elements, refine, pop nodes to boundary.
2420 562 : for (unsigned int r=0; r<nr; r++)
2421 : {
2422 : // A DistributedMesh needs a little prep before refinement, and
2423 : // may need us to keep track of ghost node movement.
2424 96 : std::unordered_set<dof_id_type> moved_ghost_nodes;
2425 329 : if (!is_replicated)
2426 84 : mesh.prepare_for_use();
2427 :
2428 329 : mesh_refinement.uniformly_refine(1);
2429 :
2430 : const bool move_only_boundary_nodes =
2431 329 : mesh.mesh_dimension() != 2 || flat;
2432 : MeshTools::Modification::interpolate_surface
2433 658 : (mesh, sphere, /*ids=*/{}, move_only_boundary_nodes);
2434 : }
2435 :
2436 : // A DistributedMesh needs a little prep before flattening
2437 233 : if (!is_replicated)
2438 42 : mesh.prepare_for_use();
2439 :
2440 : // The mesh now contains a refinement hierarchy due to the refinements
2441 : // used to generate the grid. In order to call other support functions
2442 : // like all_tri() and all_second_order, you need a "flat" mesh file (with no
2443 : // refinement trees) so
2444 233 : MeshTools::Modification::flatten(mesh);
2445 :
2446 : // Convert all the tensor product elements to simplices if requested
2447 233 : if ((type == TRI7) || (type == TRI6) || (type == TRI3) ||
2448 116 : (type == TET4) || (type == TET10) || (type == TET14))
2449 : {
2450 : // A DistributedMesh needs a little prep before all_tri()
2451 49 : if (is_replicated)
2452 42 : mesh.prepare_for_use();
2453 :
2454 49 : MeshTools::Modification::all_tri(mesh);
2455 : }
2456 :
2457 : // Convert to second-order elements if the user requested it.
2458 301 : if (Elem::build(type)->default_order() != FIRST)
2459 : {
2460 163 : if (type == TET14)
2461 0 : mesh.all_complete_order();
2462 : else
2463 : {
2464 : // type is second-order, determine if it is the
2465 : // "full-ordered" second-order element, or the "serendipity"
2466 : // second order element. Note also that all_second_order
2467 : // can't be called once the mesh has been refined.
2468 163 : bool full_ordered = !((type==QUAD8) || (type==HEX20));
2469 163 : mesh.all_second_order(full_ordered);
2470 : }
2471 :
2472 : // And pop to the boundary again...
2473 25924 : for (const auto & elem : mesh.active_element_ptr_range())
2474 124685 : for (auto s : elem->side_index_range())
2475 131158 : if (elem->neighbor_ptr(s) == nullptr)
2476 : {
2477 4178 : elem->build_side_ptr(side, s);
2478 :
2479 : // Pop each point to the sphere boundary
2480 37044 : for (auto n : side->node_index_range())
2481 42688 : side->point(n) =
2482 42688 : sphere.closest_point(side->point(n));
2483 67 : }
2484 : }
2485 :
2486 :
2487 : // The meshes could probably use some smoothing.
2488 233 : if (mesh.mesh_dimension() > 1)
2489 : {
2490 265 : LaplaceMeshSmoother smoother(mesh, n_smooth);
2491 205 : smoother.smooth();
2492 : }
2493 :
2494 : // We'll give the whole sphere surface a boundary id of 0
2495 91940 : for (const auto & elem : mesh.active_element_ptr_range())
2496 376577 : for (auto s : elem->side_index_range())
2497 378678 : if (!elem->neighbor_ptr(s))
2498 8513 : boundary_info.add_side(elem, s, 0);
2499 :
2500 : // Done building the mesh. Now prepare it for use.
2501 233 : mesh.prepare_for_use();
2502 233 : }
2503 :
2504 : #endif // #ifndef LIBMESH_ENABLE_AMR
2505 :
2506 :
2507 : // Meshes the tensor product of a 1D and a 1D-or-2D domain.
2508 49 : void MeshTools::Generation::build_extrusion (UnstructuredMesh & mesh,
2509 : const MeshBase & cross_section,
2510 : const unsigned int nz,
2511 : RealVectorValue extrusion_vector,
2512 : QueryElemSubdomainIDBase * elem_subdomain)
2513 : {
2514 14 : LOG_SCOPE("build_extrusion()", "MeshTools::Generation");
2515 :
2516 49 : if (!cross_section.n_elem())
2517 0 : return;
2518 :
2519 49 : dof_id_type orig_elem = cross_section.n_elem();
2520 49 : dof_id_type orig_nodes = cross_section.n_nodes();
2521 :
2522 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
2523 49 : unique_id_type orig_unique_ids = cross_section.parallel_max_unique_id();
2524 : #endif
2525 :
2526 49 : unsigned int order = 1;
2527 :
2528 14 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
2529 14 : const BoundaryInfo & cross_section_boundary_info = cross_section.get_boundary_info();
2530 :
2531 : // Copy name maps from old to new boundary. We won't copy the whole
2532 : // BoundaryInfo because that copies bc ids too, and we need to set
2533 : // those more carefully.
2534 14 : boundary_info.set_sideset_name_map() = cross_section_boundary_info.get_sideset_name_map();
2535 14 : boundary_info.set_nodeset_name_map() = cross_section_boundary_info.get_nodeset_name_map();
2536 14 : boundary_info.set_edgeset_name_map() = cross_section_boundary_info.get_edgeset_name_map();
2537 :
2538 : // If cross_section is distributed, so is its extrusion
2539 49 : if (!cross_section.is_serial())
2540 0 : mesh.delete_remote_elements();
2541 :
2542 : // We know a priori how many elements we'll need
2543 49 : mesh.reserve_elem(nz*orig_elem);
2544 :
2545 : // For straightforward meshes we need one or two additional layers per
2546 : // element.
2547 182 : if (cross_section.elements_begin() != cross_section.elements_end() &&
2548 143 : (*cross_section.elements_begin())->default_order() == SECOND)
2549 14 : order = 2;
2550 49 : mesh.comm().max(order);
2551 :
2552 49 : mesh.reserve_nodes((order*nz+1)*orig_nodes);
2553 :
2554 : // Container to catch the boundary IDs handed back by the BoundaryInfo object
2555 28 : std::vector<boundary_id_type> ids_to_copy;
2556 :
2557 1794 : for (const auto & node : cross_section.node_ptr_range())
2558 : {
2559 6510 : for (unsigned int k=0; k != order*nz+1; ++k)
2560 : {
2561 5313 : const dof_id_type new_node_id = node->id() + k * orig_nodes;
2562 5313 : Node * my_node = mesh.query_node_ptr(new_node_id);
2563 5313 : if (!my_node)
2564 : {
2565 : std::unique_ptr<Node> new_node = Node::build
2566 6831 : (*node + (extrusion_vector * k / nz / order),
2567 1518 : new_node_id);
2568 5313 : new_node->processor_id() = node->processor_id();
2569 :
2570 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
2571 : // Let's give the base of the extruded mesh the same
2572 : // unique_ids as the source mesh, in case anyone finds that
2573 : // a useful map to preserve.
2574 5313 : const unique_id_type uid = (k == 0) ?
2575 684 : node->unique_id() :
2576 4116 : orig_unique_ids + (k-1)*(orig_nodes + orig_elem) + node->id();
2577 :
2578 1518 : new_node->set_unique_id(uid);
2579 : #endif
2580 :
2581 5313 : cross_section_boundary_info.boundary_ids(node, ids_to_copy);
2582 5313 : boundary_info.add_node(new_node.get(), ids_to_copy);
2583 :
2584 8349 : mesh.add_node(std::move(new_node));
2585 2277 : }
2586 : }
2587 21 : }
2588 :
2589 : const std::set<boundary_id_type> & side_ids =
2590 14 : cross_section_boundary_info.get_side_boundary_ids();
2591 :
2592 14 : boundary_id_type next_side_id = side_ids.empty() ?
2593 49 : 0 : cast_int<boundary_id_type>(*side_ids.rbegin() + 1);
2594 :
2595 : // side_ids may not include ids from remote elements, in which case
2596 : // some processors may have underestimated the next_side_id; let's
2597 : // fix that.
2598 49 : cross_section.comm().max(next_side_id);
2599 :
2600 734 : for (const auto & elem : cross_section.element_ptr_range())
2601 : {
2602 455 : const ElemType etype = elem->type();
2603 :
2604 : // build_extrusion currently only works on coarse meshes
2605 130 : libmesh_assert (!elem->parent());
2606 :
2607 1421 : for (unsigned int k=0; k != nz; ++k)
2608 : {
2609 690 : std::unique_ptr<Elem> new_elem;
2610 966 : switch (etype)
2611 : {
2612 0 : case EDGE2:
2613 : {
2614 0 : new_elem = Elem::build(QUAD4);
2615 0 : new_elem->set_node(0, mesh.node_ptr(elem->node_ptr(0)->id() + (k * orig_nodes)));
2616 0 : new_elem->set_node(1, mesh.node_ptr(elem->node_ptr(1)->id() + (k * orig_nodes)));
2617 0 : new_elem->set_node(2, mesh.node_ptr(elem->node_ptr(1)->id() + ((k+1) * orig_nodes)));
2618 0 : new_elem->set_node(3, mesh.node_ptr(elem->node_ptr(0)->id() + ((k+1) * orig_nodes)));
2619 :
2620 0 : if (elem->neighbor_ptr(0) == remote_elem)
2621 0 : new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
2622 0 : if (elem->neighbor_ptr(1) == remote_elem)
2623 0 : new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
2624 :
2625 0 : break;
2626 : }
2627 0 : case EDGE3:
2628 : {
2629 0 : new_elem = Elem::build(QUAD9);
2630 0 : new_elem->set_node(0, mesh.node_ptr(elem->node_ptr(0)->id() + (2*k * orig_nodes)));
2631 0 : new_elem->set_node(1, mesh.node_ptr(elem->node_ptr(1)->id() + (2*k * orig_nodes)));
2632 0 : new_elem->set_node(2, mesh.node_ptr(elem->node_ptr(1)->id() + ((2*k+2) * orig_nodes)));
2633 0 : new_elem->set_node(3, mesh.node_ptr(elem->node_ptr(0)->id() + ((2*k+2) * orig_nodes)));
2634 0 : new_elem->set_node(4, mesh.node_ptr(elem->node_ptr(2)->id() + (2*k * orig_nodes)));
2635 0 : new_elem->set_node(5, mesh.node_ptr(elem->node_ptr(1)->id() + ((2*k+1) * orig_nodes)));
2636 0 : new_elem->set_node(6, mesh.node_ptr(elem->node_ptr(2)->id() + ((2*k+2) * orig_nodes)));
2637 0 : new_elem->set_node(7, mesh.node_ptr(elem->node_ptr(0)->id() + ((2*k+1) * orig_nodes)));
2638 0 : new_elem->set_node(8, mesh.node_ptr(elem->node_ptr(2)->id() + ((2*k+1) * orig_nodes)));
2639 :
2640 0 : if (elem->neighbor_ptr(0) == remote_elem)
2641 0 : new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
2642 0 : if (elem->neighbor_ptr(1) == remote_elem)
2643 0 : new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
2644 :
2645 0 : break;
2646 : }
2647 168 : case TRI3:
2648 : {
2649 240 : new_elem = Elem::build(PRISM6);
2650 216 : new_elem->set_node(0, mesh.node_ptr(elem->node_ptr(0)->id() + (k * orig_nodes)));
2651 216 : new_elem->set_node(1, mesh.node_ptr(elem->node_ptr(1)->id() + (k * orig_nodes)));
2652 216 : new_elem->set_node(2, mesh.node_ptr(elem->node_ptr(2)->id() + (k * orig_nodes)));
2653 216 : new_elem->set_node(3, mesh.node_ptr(elem->node_ptr(0)->id() + ((k+1) * orig_nodes)));
2654 216 : new_elem->set_node(4, mesh.node_ptr(elem->node_ptr(1)->id() + ((k+1) * orig_nodes)));
2655 216 : new_elem->set_node(5, mesh.node_ptr(elem->node_ptr(2)->id() + ((k+1) * orig_nodes)));
2656 :
2657 216 : if (elem->neighbor_ptr(0) == remote_elem)
2658 0 : new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
2659 216 : if (elem->neighbor_ptr(1) == remote_elem)
2660 0 : new_elem->set_neighbor(2, const_cast<RemoteElem *>(remote_elem));
2661 216 : if (elem->neighbor_ptr(2) == remote_elem)
2662 0 : new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
2663 :
2664 48 : break;
2665 : }
2666 0 : case TRI6:
2667 : {
2668 0 : new_elem = Elem::build(PRISM18);
2669 0 : new_elem->set_node(0, mesh.node_ptr(elem->node_ptr(0)->id() + (2*k * orig_nodes)));
2670 0 : new_elem->set_node(1, mesh.node_ptr(elem->node_ptr(1)->id() + (2*k * orig_nodes)));
2671 0 : new_elem->set_node(2, mesh.node_ptr(elem->node_ptr(2)->id() + (2*k * orig_nodes)));
2672 0 : new_elem->set_node(3, mesh.node_ptr(elem->node_ptr(0)->id() + ((2*k+2) * orig_nodes)));
2673 0 : new_elem->set_node(4, mesh.node_ptr(elem->node_ptr(1)->id() + ((2*k+2) * orig_nodes)));
2674 0 : new_elem->set_node(5, mesh.node_ptr(elem->node_ptr(2)->id() + ((2*k+2) * orig_nodes)));
2675 0 : new_elem->set_node(6, mesh.node_ptr(elem->node_ptr(3)->id() + (2*k * orig_nodes)));
2676 0 : new_elem->set_node(7, mesh.node_ptr(elem->node_ptr(4)->id() + (2*k * orig_nodes)));
2677 0 : new_elem->set_node(8, mesh.node_ptr(elem->node_ptr(5)->id() + (2*k * orig_nodes)));
2678 0 : new_elem->set_node(9, mesh.node_ptr(elem->node_ptr(0)->id() + ((2*k+1) * orig_nodes)));
2679 0 : new_elem->set_node(10, mesh.node_ptr(elem->node_ptr(1)->id() + ((2*k+1) * orig_nodes)));
2680 0 : new_elem->set_node(11, mesh.node_ptr(elem->node_ptr(2)->id() + ((2*k+1) * orig_nodes)));
2681 0 : new_elem->set_node(12, mesh.node_ptr(elem->node_ptr(3)->id() + ((2*k+2) * orig_nodes)));
2682 0 : new_elem->set_node(13, mesh.node_ptr(elem->node_ptr(4)->id() + ((2*k+2) * orig_nodes)));
2683 0 : new_elem->set_node(14, mesh.node_ptr(elem->node_ptr(5)->id() + ((2*k+2) * orig_nodes)));
2684 0 : new_elem->set_node(15, mesh.node_ptr(elem->node_ptr(3)->id() + ((2*k+1) * orig_nodes)));
2685 0 : new_elem->set_node(16, mesh.node_ptr(elem->node_ptr(4)->id() + ((2*k+1) * orig_nodes)));
2686 0 : new_elem->set_node(17, mesh.node_ptr(elem->node_ptr(5)->id() + ((2*k+1) * orig_nodes)));
2687 :
2688 0 : if (elem->neighbor_ptr(0) == remote_elem)
2689 0 : new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
2690 0 : if (elem->neighbor_ptr(1) == remote_elem)
2691 0 : new_elem->set_neighbor(2, const_cast<RemoteElem *>(remote_elem));
2692 0 : if (elem->neighbor_ptr(2) == remote_elem)
2693 0 : new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
2694 :
2695 0 : break;
2696 : }
2697 0 : case TRI7:
2698 : {
2699 0 : new_elem = Elem::build(PRISM21);
2700 0 : new_elem->set_node(0, mesh.node_ptr(elem->node_ptr(0)->id() + (2*k * orig_nodes)));
2701 0 : new_elem->set_node(1, mesh.node_ptr(elem->node_ptr(1)->id() + (2*k * orig_nodes)));
2702 0 : new_elem->set_node(2, mesh.node_ptr(elem->node_ptr(2)->id() + (2*k * orig_nodes)));
2703 0 : new_elem->set_node(3, mesh.node_ptr(elem->node_ptr(0)->id() + ((2*k+2) * orig_nodes)));
2704 0 : new_elem->set_node(4, mesh.node_ptr(elem->node_ptr(1)->id() + ((2*k+2) * orig_nodes)));
2705 0 : new_elem->set_node(5, mesh.node_ptr(elem->node_ptr(2)->id() + ((2*k+2) * orig_nodes)));
2706 0 : new_elem->set_node(6, mesh.node_ptr(elem->node_ptr(3)->id() + (2*k * orig_nodes)));
2707 0 : new_elem->set_node(7, mesh.node_ptr(elem->node_ptr(4)->id() + (2*k * orig_nodes)));
2708 0 : new_elem->set_node(8, mesh.node_ptr(elem->node_ptr(5)->id() + (2*k * orig_nodes)));
2709 0 : new_elem->set_node(9, mesh.node_ptr(elem->node_ptr(0)->id() + ((2*k+1) * orig_nodes)));
2710 0 : new_elem->set_node(10, mesh.node_ptr(elem->node_ptr(1)->id() + ((2*k+1) * orig_nodes)));
2711 0 : new_elem->set_node(11, mesh.node_ptr(elem->node_ptr(2)->id() + ((2*k+1) * orig_nodes)));
2712 0 : new_elem->set_node(12, mesh.node_ptr(elem->node_ptr(3)->id() + ((2*k+2) * orig_nodes)));
2713 0 : new_elem->set_node(13, mesh.node_ptr(elem->node_ptr(4)->id() + ((2*k+2) * orig_nodes)));
2714 0 : new_elem->set_node(14, mesh.node_ptr(elem->node_ptr(5)->id() + ((2*k+2) * orig_nodes)));
2715 0 : new_elem->set_node(15, mesh.node_ptr(elem->node_ptr(3)->id() + ((2*k+1) * orig_nodes)));
2716 0 : new_elem->set_node(16, mesh.node_ptr(elem->node_ptr(4)->id() + ((2*k+1) * orig_nodes)));
2717 0 : new_elem->set_node(17, mesh.node_ptr(elem->node_ptr(5)->id() + ((2*k+1) * orig_nodes)));
2718 :
2719 0 : new_elem->set_node(18, mesh.node_ptr(elem->node_ptr(6)->id() + (2*k * orig_nodes)));
2720 0 : new_elem->set_node(19, mesh.node_ptr(elem->node_ptr(6)->id() + ((2*k+2) * orig_nodes)));
2721 0 : new_elem->set_node(20, mesh.node_ptr(elem->node_ptr(6)->id() + ((2*k+1) * orig_nodes)));
2722 :
2723 0 : if (elem->neighbor_ptr(0) == remote_elem)
2724 0 : new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
2725 0 : if (elem->neighbor_ptr(1) == remote_elem)
2726 0 : new_elem->set_neighbor(2, const_cast<RemoteElem *>(remote_elem));
2727 0 : if (elem->neighbor_ptr(2) == remote_elem)
2728 0 : new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
2729 :
2730 0 : break;
2731 : }
2732 448 : case QUAD4:
2733 : {
2734 640 : new_elem = Elem::build(HEX8);
2735 576 : new_elem->set_node(0, mesh.node_ptr(elem->node_ptr(0)->id() + (k * orig_nodes)));
2736 576 : new_elem->set_node(1, mesh.node_ptr(elem->node_ptr(1)->id() + (k * orig_nodes)));
2737 576 : new_elem->set_node(2, mesh.node_ptr(elem->node_ptr(2)->id() + (k * orig_nodes)));
2738 576 : new_elem->set_node(3, mesh.node_ptr(elem->node_ptr(3)->id() + (k * orig_nodes)));
2739 576 : new_elem->set_node(4, mesh.node_ptr(elem->node_ptr(0)->id() + ((k+1) * orig_nodes)));
2740 576 : new_elem->set_node(5, mesh.node_ptr(elem->node_ptr(1)->id() + ((k+1) * orig_nodes)));
2741 576 : new_elem->set_node(6, mesh.node_ptr(elem->node_ptr(2)->id() + ((k+1) * orig_nodes)));
2742 576 : new_elem->set_node(7, mesh.node_ptr(elem->node_ptr(3)->id() + ((k+1) * orig_nodes)));
2743 :
2744 576 : if (elem->neighbor_ptr(0) == remote_elem)
2745 0 : new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
2746 576 : if (elem->neighbor_ptr(1) == remote_elem)
2747 0 : new_elem->set_neighbor(2, const_cast<RemoteElem *>(remote_elem));
2748 576 : if (elem->neighbor_ptr(2) == remote_elem)
2749 0 : new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
2750 576 : if (elem->neighbor_ptr(3) == remote_elem)
2751 0 : new_elem->set_neighbor(4, const_cast<RemoteElem *>(remote_elem));
2752 :
2753 128 : break;
2754 : }
2755 350 : case QUAD9:
2756 : {
2757 500 : new_elem = Elem::build(HEX27);
2758 450 : new_elem->set_node(0, mesh.node_ptr(elem->node_ptr(0)->id() + (2*k * orig_nodes)));
2759 450 : new_elem->set_node(1, mesh.node_ptr(elem->node_ptr(1)->id() + (2*k * orig_nodes)));
2760 450 : new_elem->set_node(2, mesh.node_ptr(elem->node_ptr(2)->id() + (2*k * orig_nodes)));
2761 450 : new_elem->set_node(3, mesh.node_ptr(elem->node_ptr(3)->id() + (2*k * orig_nodes)));
2762 450 : new_elem->set_node(4, mesh.node_ptr(elem->node_ptr(0)->id() + ((2*k+2) * orig_nodes)));
2763 450 : new_elem->set_node(5, mesh.node_ptr(elem->node_ptr(1)->id() + ((2*k+2) * orig_nodes)));
2764 450 : new_elem->set_node(6, mesh.node_ptr(elem->node_ptr(2)->id() + ((2*k+2) * orig_nodes)));
2765 450 : new_elem->set_node(7, mesh.node_ptr(elem->node_ptr(3)->id() + ((2*k+2) * orig_nodes)));
2766 450 : new_elem->set_node(8, mesh.node_ptr(elem->node_ptr(4)->id() + (2*k * orig_nodes)));
2767 450 : new_elem->set_node(9, mesh.node_ptr(elem->node_ptr(5)->id() + (2*k * orig_nodes)));
2768 450 : new_elem->set_node(10, mesh.node_ptr(elem->node_ptr(6)->id() + (2*k * orig_nodes)));
2769 450 : new_elem->set_node(11, mesh.node_ptr(elem->node_ptr(7)->id() + (2*k * orig_nodes)));
2770 450 : new_elem->set_node(12, mesh.node_ptr(elem->node_ptr(0)->id() + ((2*k+1) * orig_nodes)));
2771 450 : new_elem->set_node(13, mesh.node_ptr(elem->node_ptr(1)->id() + ((2*k+1) * orig_nodes)));
2772 450 : new_elem->set_node(14, mesh.node_ptr(elem->node_ptr(2)->id() + ((2*k+1) * orig_nodes)));
2773 450 : new_elem->set_node(15, mesh.node_ptr(elem->node_ptr(3)->id() + ((2*k+1) * orig_nodes)));
2774 450 : new_elem->set_node(16, mesh.node_ptr(elem->node_ptr(4)->id() + ((2*k+2) * orig_nodes)));
2775 450 : new_elem->set_node(17, mesh.node_ptr(elem->node_ptr(5)->id() + ((2*k+2) * orig_nodes)));
2776 450 : new_elem->set_node(18, mesh.node_ptr(elem->node_ptr(6)->id() + ((2*k+2) * orig_nodes)));
2777 450 : new_elem->set_node(19, mesh.node_ptr(elem->node_ptr(7)->id() + ((2*k+2) * orig_nodes)));
2778 450 : new_elem->set_node(20, mesh.node_ptr(elem->node_ptr(8)->id() + (2*k * orig_nodes)));
2779 450 : new_elem->set_node(21, mesh.node_ptr(elem->node_ptr(4)->id() + ((2*k+1) * orig_nodes)));
2780 450 : new_elem->set_node(22, mesh.node_ptr(elem->node_ptr(5)->id() + ((2*k+1) * orig_nodes)));
2781 450 : new_elem->set_node(23, mesh.node_ptr(elem->node_ptr(6)->id() + ((2*k+1) * orig_nodes)));
2782 450 : new_elem->set_node(24, mesh.node_ptr(elem->node_ptr(7)->id() + ((2*k+1) * orig_nodes)));
2783 450 : new_elem->set_node(25, mesh.node_ptr(elem->node_ptr(8)->id() + ((2*k+2) * orig_nodes)));
2784 450 : new_elem->set_node(26, mesh.node_ptr(elem->node_ptr(8)->id() + ((2*k+1) * orig_nodes)));
2785 :
2786 450 : if (elem->neighbor_ptr(0) == remote_elem)
2787 0 : new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
2788 450 : if (elem->neighbor_ptr(1) == remote_elem)
2789 0 : new_elem->set_neighbor(2, const_cast<RemoteElem *>(remote_elem));
2790 450 : if (elem->neighbor_ptr(2) == remote_elem)
2791 0 : new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
2792 450 : if (elem->neighbor_ptr(3) == remote_elem)
2793 0 : new_elem->set_neighbor(4, const_cast<RemoteElem *>(remote_elem));
2794 :
2795 100 : break;
2796 : }
2797 0 : default:
2798 : {
2799 0 : libmesh_not_implemented();
2800 : break;
2801 : }
2802 : }
2803 :
2804 966 : new_elem->set_id(elem->id() + (k * orig_elem));
2805 966 : new_elem->processor_id() = elem->processor_id();
2806 :
2807 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
2808 : // Let's give the base of the extruded mesh the same
2809 : // unique_ids as the source mesh, in case anyone finds that
2810 : // a useful map to preserve.
2811 966 : const unique_id_type uid = (k == 0) ?
2812 260 : elem->unique_id() :
2813 511 : orig_unique_ids + (k-1)*(orig_nodes + orig_elem) + orig_nodes + elem->id();
2814 :
2815 276 : new_elem->set_unique_id(uid);
2816 : #endif
2817 :
2818 966 : if (!elem_subdomain)
2819 : // maintain the subdomain_id
2820 518 : new_elem->subdomain_id() = elem->subdomain_id();
2821 : else
2822 : // Allow the user to choose new subdomain_ids
2823 448 : new_elem->subdomain_id() = elem_subdomain->get_subdomain_for_layer(elem, k);
2824 :
2825 1242 : Elem * added_elem = mesh.add_elem(std::move(new_elem));
2826 :
2827 : // Copy any old boundary ids on all sides
2828 4938 : for (auto s : elem->side_index_range())
2829 : {
2830 3696 : cross_section_boundary_info.boundary_ids(elem, s, ids_to_copy);
2831 :
2832 3696 : if (added_elem->dim() == 3)
2833 : {
2834 : // For 2D->3D extrusion, we give the boundary IDs
2835 : // for side s on the old element to side s+1 on the
2836 : // new element. This is just a happy coincidence as
2837 : // far as I can tell...
2838 3696 : boundary_info.add_side(added_elem,
2839 2640 : cast_int<unsigned short>(s+1),
2840 : ids_to_copy);
2841 : }
2842 : else
2843 : {
2844 : // For 1D->2D extrusion, the boundary IDs map as:
2845 : // Old elem -> New elem
2846 : // 0 -> 3
2847 : // 1 -> 1
2848 0 : libmesh_assert_less(s, 2);
2849 0 : const unsigned short sidemap[2] = {3, 1};
2850 0 : boundary_info.add_side(added_elem, sidemap[s], ids_to_copy);
2851 : }
2852 : }
2853 :
2854 : // Give new boundary ids to bottom and top
2855 966 : if (k == 0)
2856 455 : boundary_info.add_side(added_elem, 0, next_side_id);
2857 966 : if (k == nz-1)
2858 : {
2859 : // For 2D->3D extrusion, the "top" ID is 1+the original
2860 : // element's number of sides. For 1D->2D extrusion, the
2861 : // "top" ID is side 2.
2862 455 : const unsigned short top_id = added_elem->dim() == 3 ?
2863 455 : cast_int<unsigned short>(elem->n_sides()+1) : 2;
2864 : boundary_info.add_side
2865 455 : (added_elem, top_id,
2866 455 : cast_int<boundary_id_type>(next_side_id+1));
2867 : }
2868 414 : }
2869 21 : }
2870 :
2871 : // Done building the mesh. Now prepare it for use.
2872 49 : mesh.prepare_for_use();
2873 : }
2874 :
2875 :
2876 :
2877 :
2878 : #if defined(LIBMESH_HAVE_TRIANGLE) && LIBMESH_DIM > 1
2879 :
2880 : // Triangulates a 2D rectangular region with or without holes
2881 0 : void MeshTools::Generation::build_delaunay_square(UnstructuredMesh & mesh,
2882 : const unsigned int nx, // num. of elements in x-dir
2883 : const unsigned int ny, // num. of elements in y-dir
2884 : const Real xmin, const Real xmax,
2885 : const Real ymin, const Real ymax,
2886 : const ElemType type,
2887 : const std::vector<TriangleInterface::Hole*> * holes)
2888 : {
2889 : // Check for reasonable size
2890 0 : libmesh_assert_greater_equal (nx, 1); // need at least 1 element in x-direction
2891 0 : libmesh_assert_greater_equal (ny, 1); // need at least 1 element in y-direction
2892 0 : libmesh_assert_less (xmin, xmax);
2893 0 : libmesh_assert_less (ymin, ymax);
2894 :
2895 : // Clear out any data which may have been in the Mesh
2896 0 : mesh.clear();
2897 :
2898 0 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
2899 :
2900 : // Make sure the new Mesh will be 2D
2901 0 : mesh.set_mesh_dimension(2);
2902 :
2903 : // The x and y spacing between boundary points
2904 0 : const Real delta_x = (xmax-xmin) / static_cast<Real>(nx);
2905 0 : const Real delta_y = (ymax-ymin) / static_cast<Real>(ny);
2906 :
2907 : // Bottom
2908 0 : for (unsigned int p=0; p<=nx; ++p)
2909 0 : mesh.add_point(Point(xmin + p*delta_x, ymin));
2910 :
2911 : // Right side
2912 0 : for (unsigned int p=1; p<ny; ++p)
2913 0 : mesh.add_point(Point(xmax, ymin + p*delta_y));
2914 :
2915 : // Top
2916 0 : for (unsigned int p=0; p<=nx; ++p)
2917 0 : mesh.add_point(Point(xmax - p*delta_x, ymax));
2918 :
2919 : // Left side
2920 0 : for (unsigned int p=1; p<ny; ++p)
2921 0 : mesh.add_point(Point(xmin, ymax - p*delta_y));
2922 :
2923 : // Be sure we added as many points as we thought we did
2924 0 : libmesh_assert_equal_to (mesh.n_nodes(), 2*(nx+ny));
2925 :
2926 : // Construct the Triangle Interface object
2927 0 : TriangleInterface t(mesh);
2928 :
2929 : // Set custom variables for the triangulation
2930 0 : t.desired_area() = 0.5 * (xmax-xmin)*(ymax-ymin) / static_cast<Real>(nx*ny);
2931 0 : t.triangulation_type() = TriangleInterface::PSLG;
2932 0 : t.elem_type() = type;
2933 :
2934 0 : if (holes != nullptr)
2935 0 : t.attach_hole_list(holes);
2936 :
2937 : // Triangulate!
2938 0 : t.triangulate();
2939 :
2940 : // For avoiding extraneous side element construction
2941 0 : std::unique_ptr<const Elem> side;
2942 :
2943 : // The mesh is now generated, but we still need to mark the boundaries
2944 : // to be consistent with the other build_square routines. Note that all
2945 : // hole boundary elements get the same ID, 4.
2946 0 : for (auto & elem : mesh.element_ptr_range())
2947 0 : for (auto s : elem->side_index_range())
2948 0 : if (elem->neighbor_ptr(s) == nullptr)
2949 : {
2950 0 : elem->build_side_ptr(side, s);
2951 :
2952 : // Check the location of the side's midpoint. Since
2953 : // the square has straight sides, the midpoint is not
2954 : // on the corner and thus it is uniquely on one of the
2955 : // sides.
2956 0 : Point side_midpoint= 0.5f*( side->point(0) + side->point(1) );
2957 :
2958 : // The boundary ids are set following the same convention as Quad4 sides
2959 : // bottom = 0
2960 : // right = 1
2961 : // top = 2
2962 : // left = 3
2963 : // hole = 4
2964 0 : boundary_id_type bc_id=4;
2965 :
2966 : // bottom
2967 0 : if (std::fabs(side_midpoint(1) - ymin) < TOLERANCE)
2968 0 : bc_id=0;
2969 :
2970 : // right
2971 0 : else if (std::fabs(side_midpoint(0) - xmax) < TOLERANCE)
2972 0 : bc_id=1;
2973 :
2974 : // top
2975 0 : else if (std::fabs(side_midpoint(1) - ymax) < TOLERANCE)
2976 0 : bc_id=2;
2977 :
2978 : // left
2979 0 : else if (std::fabs(side_midpoint(0) - xmin) < TOLERANCE)
2980 0 : bc_id=3;
2981 :
2982 : // If the point is not on any of the external boundaries, it
2983 : // is on one of the holes....
2984 :
2985 : // Finally, add this element's information to the boundary info object.
2986 0 : boundary_info.add_side(elem->id(), s, bc_id);
2987 : }
2988 :
2989 0 : } // end build_delaunay_square
2990 :
2991 : #endif // LIBMESH_HAVE_TRIANGLE && LIBMESH_DIM > 1
2992 :
2993 :
2994 42 : void MeshTools::Generation::surface_octahedron
2995 : (UnstructuredMesh & mesh,
2996 : Real xmin, Real xmax,
2997 : Real ymin, Real ymax,
2998 : Real zmin, Real zmax,
2999 : bool flip_tris)
3000 : {
3001 42 : const Real xavg = (xmin + xmax)/2;
3002 42 : const Real yavg = (ymin + ymax)/2;
3003 42 : const Real zavg = (zmin + zmax)/2;
3004 54 : mesh.add_point(Point(xavg,yavg,zmin), 0);
3005 54 : mesh.add_point(Point(xmax,yavg,zavg), 1);
3006 54 : mesh.add_point(Point(xavg,ymax,zavg), 2);
3007 54 : mesh.add_point(Point(xmin,yavg,zavg), 3);
3008 54 : mesh.add_point(Point(xavg,ymin,zavg), 4);
3009 54 : mesh.add_point(Point(xavg,yavg,zmax), 5);
3010 :
3011 2560 : auto add_tri = [&mesh, flip_tris](std::array<dof_id_type,3> nodes)
3012 : {
3013 336 : auto elem = mesh.add_elem(Elem::build(TRI3));
3014 336 : elem->set_node(0, mesh.node_ptr(nodes[0]));
3015 336 : elem->set_node(1, mesh.node_ptr(nodes[1]));
3016 336 : elem->set_node(2, mesh.node_ptr(nodes[2]));
3017 336 : if (flip_tris)
3018 144 : elem->flip(&mesh.get_boundary_info());
3019 348 : };
3020 :
3021 42 : add_tri({0,2,1});
3022 42 : add_tri({0,3,2});
3023 42 : add_tri({0,4,3});
3024 42 : add_tri({0,1,4});
3025 42 : add_tri({5,4,1});
3026 42 : add_tri({5,3,4});
3027 42 : add_tri({5,2,3});
3028 42 : add_tri({5,1,2});
3029 :
3030 42 : mesh.prepare_for_use();
3031 42 : }
3032 :
3033 :
3034 : } // namespace libMesh
|
