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