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 1617406 : unsigned int idx(const ElemType type,
83 : const unsigned int nx,
84 : const unsigned int i,
85 : const unsigned int j)
86 : {
87 1165292 : switch(type)
88 : {
89 328248 : case INVALID_ELEM:
90 : case QUAD4:
91 : case QUADSHELL4:
92 : case TRI3:
93 : case TRISHELL3:
94 : {
95 328248 : return i + j*(nx+1);
96 : }
97 :
98 1289158 : case QUAD8:
99 : case QUADSHELL8:
100 : case QUAD9:
101 : case QUADSHELL9:
102 : case TRI6:
103 : case TRI7:
104 : {
105 1289158 : 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 4297350 : 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 3116486 : switch(type)
127 : {
128 851852 : case INVALID_ELEM:
129 : case HEX8:
130 : case PRISM6:
131 : case C0POLYHEDRON:
132 : {
133 851852 : return i + (nx+1)*(j + k*(ny+1));
134 : }
135 :
136 3445498 : 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 3445498 : 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 40530 : 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 34809 : 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 40530 : mesh.clear();
341 :
342 11603 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
343 :
344 40530 : if (nz != 0)
345 : {
346 18552 : mesh.set_mesh_dimension(3);
347 18552 : mesh.set_spatial_dimension(3);
348 : }
349 21978 : else if (ny != 0)
350 : {
351 15838 : mesh.set_mesh_dimension(2);
352 15838 : mesh.set_spatial_dimension(2);
353 : }
354 6140 : else if (nx != 0)
355 : {
356 5849 : mesh.set_mesh_dimension(1);
357 5849 : mesh.set_spatial_dimension(1);
358 : }
359 : else
360 : {
361 : // Will we get here?
362 291 : mesh.set_mesh_dimension(0);
363 291 : mesh.set_spatial_dimension(0);
364 : }
365 :
366 40530 : switch (mesh.mesh_dimension())
367 : {
368 : //---------------------------------------------------------------------
369 : // Build a 0D point
370 291 : 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 375 : mesh.add_point (Point(0, 0, 0), 0);
380 291 : Elem * elem = mesh.add_elem(Elem::build(NODEELEM));
381 291 : elem->set_node(0, mesh.node_ptr(0));
382 :
383 207 : break;
384 : }
385 :
386 :
387 :
388 : //---------------------------------------------------------------------
389 : // Build a 1D line
390 5849 : 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 5849 : case INVALID_ELEM:
401 : case EDGE2:
402 : case EDGE3:
403 : case EDGE4:
404 : {
405 5849 : 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 1917 : case INVALID_ELEM:
417 : case EDGE2:
418 : {
419 1917 : mesh.reserve_nodes(nx+1);
420 1367 : break;
421 : }
422 :
423 3225 : case EDGE3:
424 : {
425 3225 : mesh.reserve_nodes(2*nx+1);
426 2303 : break;
427 : }
428 :
429 707 : case EDGE4:
430 : {
431 707 : mesh.reserve_nodes(3*nx+1);
432 505 : 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 31868 : for (unsigned int i=0; i<=nx; i++)
449 : {
450 39697 : const Node * const node = mesh.add_point (Point(static_cast<Real>(i)/nx, 0, 0), node_id++);
451 29951 : if (i == 0)
452 1917 : boundary_info.add_node(node, 0);
453 29951 : if (i == nx)
454 1917 : boundary_info.add_node(node, 1);
455 : }
456 :
457 550 : break;
458 : }
459 :
460 922 : case EDGE3:
461 : {
462 15860 : for (unsigned int i=0; i<=2*nx; i++)
463 : {
464 16257 : const Node * const node = mesh.add_point (Point(static_cast<Real>(i)/(2*nx), 0, 0), node_id++);
465 12635 : if (i == 0)
466 3225 : boundary_info.add_node(node, 0);
467 12635 : if (i == 2*nx)
468 3225 : boundary_info.add_node(node, 1);
469 : }
470 922 : break;
471 : }
472 :
473 202 : case EDGE4:
474 : {
475 4690 : for (unsigned int i=0; i<=3*nx; i++)
476 : {
477 5121 : const Node * const node = mesh.add_point (Point(static_cast<Real>(i)/(3*nx), 0, 0), node_id++);
478 3983 : if (i == 0)
479 707 : boundary_info.add_node(node, 0);
480 3983 : if (i == 3*nx)
481 707 : 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 29951 : for (unsigned int i=0; i<nx; i++)
499 : {
500 28034 : Elem * elem = mesh.add_elem(Elem::build_with_id(EDGE2, i));
501 28034 : elem->set_node(0, mesh.node_ptr(i));
502 28034 : elem->set_node(1, mesh.node_ptr(i+1));
503 :
504 28034 : if (i == 0)
505 1917 : boundary_info.add_side(elem, 0, 0);
506 :
507 28034 : if (i == (nx-1))
508 1917 : boundary_info.add_side(elem, 1, 1);
509 :
510 : }
511 550 : break;
512 : }
513 :
514 922 : case EDGE3:
515 : {
516 7930 : for (unsigned int i=0; i<nx; i++)
517 : {
518 4705 : Elem * elem = mesh.add_elem(Elem::build_with_id(EDGE3, i));
519 4705 : elem->set_node(0, mesh.node_ptr(2*i));
520 4705 : elem->set_node(2, mesh.node_ptr(2*i+1));
521 4705 : elem->set_node(1, mesh.node_ptr(2*i+2));
522 :
523 4705 : if (i == 0)
524 3225 : boundary_info.add_side(elem, 0, 0);
525 :
526 4705 : if (i == (nx-1))
527 3225 : boundary_info.add_side(elem, 1, 1);
528 : }
529 922 : break;
530 : }
531 :
532 202 : case EDGE4:
533 : {
534 1799 : for (unsigned int i=0; i<nx; i++)
535 : {
536 1092 : Elem * elem = mesh.add_elem(Elem::build_with_id(EDGE4, i));
537 1092 : elem->set_node(0, mesh.node_ptr(3*i));
538 1092 : elem->set_node(2, mesh.node_ptr(3*i+1));
539 1092 : elem->set_node(3, mesh.node_ptr(3*i+2));
540 1092 : elem->set_node(1, mesh.node_ptr(3*i+3));
541 :
542 1092 : if (i == 0)
543 707 : boundary_info.add_side(elem, 0, 0);
544 :
545 1092 : if (i == (nx-1))
546 707 : 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 5849 : 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 5849 : (mesh.node_stored_range(),
565 61075 : [xmin, xmax](Node * node)
566 : {
567 46569 : (*node)(0) = (*node)(0)*(xmax-xmin) + xmin;
568 14506 : });
569 : }
570 :
571 : // Add sideset names to boundary info
572 5849 : boundary_info.sideset_name(0) = "left";
573 5849 : boundary_info.sideset_name(1) = "right";
574 :
575 : // Add nodeset names to boundary info
576 5849 : boundary_info.nodeset_name(0) = "left";
577 5849 : 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 15838 : case 2:
594 : {
595 4535 : libmesh_assert_not_equal_to (nx, 0);
596 4535 : libmesh_assert_not_equal_to (ny, 0);
597 4535 : libmesh_assert_equal_to (nz, 0);
598 4535 : libmesh_assert_less (xmin, xmax);
599 4535 : 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 8405 : case INVALID_ELEM:
606 : case QUAD4:
607 : case QUADSHELL4:
608 : case QUAD8:
609 : case QUADSHELL8:
610 : case QUAD9:
611 : case QUADSHELL9:
612 : {
613 8405 : mesh.reserve_elem (nx*ny);
614 5994 : break;
615 : }
616 :
617 7349 : case TRI3:
618 : case TRISHELL3:
619 : case TRI6:
620 : case TRI7:
621 : {
622 7349 : mesh.reserve_elem (2*nx*ny);
623 5249 : break;
624 : }
625 :
626 84 : case C0POLYGON:
627 : {
628 84 : mesh.reserve_elem ((nx + 1) * (ny + 1));
629 60 : 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 4194 : case INVALID_ELEM:
643 : case QUAD4:
644 : case QUADSHELL4:
645 : case TRI3:
646 : case TRISHELL3:
647 : {
648 4194 : mesh.reserve_nodes( (nx+1)*(ny+1) );
649 2987 : break;
650 : }
651 :
652 8679 : case QUAD8:
653 : case QUADSHELL8:
654 : case QUAD9:
655 : case QUADSHELL9:
656 : case TRI6:
657 : {
658 8679 : mesh.reserve_nodes( (2*nx+1)*(2*ny+1) );
659 6197 : break;
660 : }
661 :
662 2881 : case TRI7:
663 : {
664 2881 : mesh.reserve_nodes( (2*nx+1)*(2*ny+1) + 2*nx*ny );
665 2059 : break;
666 : }
667 84 : case C0POLYGON:
668 : {
669 84 : mesh.reserve_nodes (4 + 3*nx*ny + 2*nx + 2*ny);
670 60 : 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 4535 : unsigned int node_id = 0;
683 : switch (type)
684 : {
685 1207 : case INVALID_ELEM:
686 : case QUAD4:
687 : case QUADSHELL4:
688 : case TRI3:
689 : case TRISHELL3:
690 : {
691 17534 : for (unsigned int j=0; j<=ny; j++)
692 115180 : for (unsigned int i=0; i<=nx; i++)
693 : {
694 : const Node * const node =
695 143510 : mesh.add_point(Point(static_cast<Real>(i) / static_cast<Real>(nx),
696 101840 : static_cast<Real>(j) / static_cast<Real>(ny),
697 : 0.),
698 91101 : node_id++);
699 101840 : if (j == 0)
700 12784 : boundary_info.add_node(node, 0);
701 101840 : if (j == ny)
702 12784 : boundary_info.add_node(node, 2);
703 101840 : if (i == 0)
704 13340 : boundary_info.add_node(node, 3);
705 101840 : if (i == nx)
706 13340 : boundary_info.add_node(node, 1);
707 : }
708 :
709 1207 : break;
710 : }
711 :
712 3304 : case QUAD8:
713 : case QUADSHELL8:
714 : case QUAD9:
715 : case QUADSHELL9:
716 : case TRI6:
717 : case TRI7:
718 : {
719 64612 : for (unsigned int j=0; j<=(2*ny); j++)
720 665396 : for (unsigned int i=0; i<=(2*nx); i++)
721 : {
722 : const Node * const node =
723 885664 : mesh.add_point(Point(static_cast<Real>(i) / static_cast<Real>(2 * nx),
724 612344 : static_cast<Real>(j) / static_cast<Real>(2 * ny),
725 : 0),
726 529672 : node_id++);
727 612344 : if (j == 0)
728 54364 : boundary_info.add_node(node, 0);
729 612344 : if (j == 2*ny)
730 54364 : boundary_info.add_node(node, 2);
731 612344 : if (i == 0)
732 53052 : boundary_info.add_node(node, 3);
733 612344 : if (i == 2*nx)
734 53052 : 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 11560 : if (type == TRI7)
740 7228 : for (unsigned int j=0; j<(3*ny); j += 3)
741 25442 : for (unsigned int i=0; i<(3*nx); i += 3)
742 : {
743 : // The bottom-right triangle's center node
744 30850 : mesh.add_point(Point(static_cast<Real>(i+2) / static_cast<Real>(3 * nx),
745 21095 : static_cast<Real>(j+1) / static_cast<Real>(3 * ny),
746 : 0),
747 18260 : node_id++);
748 : // The top-left triangle's center node
749 21095 : mesh.add_point(Point(static_cast<Real>(i+1) / static_cast<Real>(3 * nx),
750 21095 : static_cast<Real>(j+2) / static_cast<Real>(3 * ny),
751 : 0),
752 18260 : node_id++);
753 : }
754 :
755 3304 : 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 4535 : unsigned int elem_id = 0;
775 : switch (type)
776 : {
777 :
778 767 : case INVALID_ELEM:
779 : case QUAD4:
780 : case QUADSHELL4:
781 : {
782 9597 : for (unsigned int j=0; j<ny; j++)
783 82548 : for (unsigned int i=0; i<nx; i++)
784 : {
785 77240 : Elem * elem = mesh.add_elem(Elem::build_with_id(type == INVALID_ELEM ? QUAD4 : type, elem_id++));
786 75606 : elem->set_node(0, mesh.node_ptr(idx(type,nx,i,j) ));
787 75606 : elem->set_node(1, mesh.node_ptr(idx(type,nx,i+1,j) ));
788 75606 : elem->set_node(2, mesh.node_ptr(idx(type,nx,i+1,j+1)));
789 75606 : elem->set_node(3, mesh.node_ptr(idx(type,nx,i,j+1) ));
790 :
791 75606 : if (j == 0)
792 6372 : boundary_info.add_side(elem, 0, 0);
793 :
794 75606 : if (j == (ny-1))
795 6372 : boundary_info.add_side(elem, 2, 2);
796 :
797 75606 : if (i == 0)
798 6942 : boundary_info.add_side(elem, 3, 3);
799 :
800 75606 : if (i == (nx-1))
801 6942 : boundary_info.add_side(elem, 1, 1);
802 : }
803 767 : break;
804 : }
805 :
806 :
807 440 : case TRI3:
808 : case TRISHELL3:
809 : {
810 3743 : for (unsigned int j=0; j<ny; j++)
811 6508 : for (unsigned int i=0; i<nx; i++)
812 : {
813 : // Add first Tri3
814 4304 : Elem * elem = mesh.add_elem(Elem::build_with_id(type, elem_id++));
815 4304 : elem->set_node(0, mesh.node_ptr(idx(type,nx,i,j) ));
816 4304 : elem->set_node(1, mesh.node_ptr(idx(type,nx,i+1,j) ));
817 4304 : elem->set_node(2, mesh.node_ptr(idx(type,nx,i+1,j+1)));
818 :
819 4304 : if (j == 0)
820 2218 : boundary_info.add_side(elem, 0, 0);
821 :
822 4304 : if (i == (nx-1))
823 2204 : boundary_info.add_side(elem, 1, 1);
824 :
825 : // Add second Tri3
826 4304 : elem = mesh.add_elem(Elem::build_with_id(type, elem_id++));
827 4304 : elem->set_node(0, mesh.node_ptr(idx(type,nx,i,j) ));
828 4304 : elem->set_node(1, mesh.node_ptr(idx(type,nx,i+1,j+1)));
829 4304 : elem->set_node(2, mesh.node_ptr(idx(type,nx,i,j+1) ));
830 :
831 4304 : if (j == (ny-1))
832 2218 : boundary_info.add_side(elem, 1, 2);
833 :
834 4304 : if (i == 0)
835 2204 : boundary_info.add_side(elem, 2, 3);
836 : }
837 440 : break;
838 : }
839 :
840 :
841 :
842 1644 : case QUAD8:
843 : case QUADSHELL8:
844 : case QUAD9:
845 : case QUADSHELL9:
846 : {
847 17209 : for (unsigned int j=0; j<(2*ny); j += 2)
848 92379 : for (unsigned int i=0; i<(2*nx); i += 2)
849 : {
850 80920 : Elem * elem = mesh.add_elem(Elem::build_with_id(type, elem_id++));
851 80920 : elem->set_node(0, mesh.node_ptr(idx(type,nx,i,j) ));
852 80920 : elem->set_node(1, mesh.node_ptr(idx(type,nx,i+2,j) ));
853 80920 : elem->set_node(2, mesh.node_ptr(idx(type,nx,i+2,j+2)));
854 80920 : elem->set_node(3, mesh.node_ptr(idx(type,nx,i,j+2) ));
855 80920 : elem->set_node(4, mesh.node_ptr(idx(type,nx,i+1,j) ));
856 80920 : elem->set_node(5, mesh.node_ptr(idx(type,nx,i+2,j+1)));
857 80920 : elem->set_node(6, mesh.node_ptr(idx(type,nx,i+1,j+2)));
858 80920 : elem->set_node(7, mesh.node_ptr(idx(type,nx,i,j+1) ));
859 :
860 80920 : if (type == QUAD9 || type == QUADSHELL9)
861 63374 : elem->set_node(8, mesh.node_ptr(idx(type,nx,i+1,j+1)));
862 :
863 80920 : if (j == 0)
864 12012 : boundary_info.add_side(elem, 0, 0);
865 :
866 80920 : if (j == 2*(ny-1))
867 12012 : boundary_info.add_side(elem, 2, 2);
868 :
869 80920 : if (i == 0)
870 11459 : boundary_info.add_side(elem, 3, 3);
871 :
872 80920 : if (i == 2*(nx-1))
873 11459 : boundary_info.add_side(elem, 1, 1);
874 : }
875 1644 : break;
876 : }
877 :
878 :
879 1660 : case TRI6:
880 : case TRI7:
881 : {
882 15097 : for (unsigned int j=0; j<(2*ny); j += 2)
883 57489 : for (unsigned int i=0; i<(2*nx); i += 2)
884 : {
885 : // Add first Tri in the bottom-right of its quad
886 48202 : Elem * elem = mesh.add_elem(Elem::build_with_id(type, elem_id++));
887 48202 : elem->set_node(0, mesh.node_ptr(idx(type,nx,i,j) ));
888 48202 : elem->set_node(1, mesh.node_ptr(idx(type,nx,i+2,j) ));
889 48202 : elem->set_node(2, mesh.node_ptr(idx(type,nx,i+2,j+2)));
890 48202 : elem->set_node(3, mesh.node_ptr(idx(type,nx,i+1,j) ));
891 48202 : elem->set_node(4, mesh.node_ptr(idx(type,nx,i+2,j+1)));
892 48202 : elem->set_node(5, mesh.node_ptr(idx(type,nx,i+1,j+1)));
893 :
894 48202 : if (type == TRI7)
895 21095 : elem->set_node(6, mesh.node_ptr(elem->id()+(2*nx+1)*(2*ny+1)));
896 :
897 48202 : if (j == 0)
898 9390 : boundary_info.add_side(elem, 0, 0);
899 :
900 48202 : if (i == 2*(nx-1))
901 9287 : boundary_info.add_side(elem, 1, 1);
902 :
903 : // Add second Tri in the top left of its quad
904 48202 : elem = mesh.add_elem(Elem::build_with_id(type, elem_id++));
905 48202 : elem->set_node(0, mesh.node_ptr(idx(type,nx,i,j) ));
906 48202 : elem->set_node(1, mesh.node_ptr(idx(type,nx,i+2,j+2)));
907 48202 : elem->set_node(2, mesh.node_ptr(idx(type,nx,i,j+2) ));
908 48202 : elem->set_node(3, mesh.node_ptr(idx(type,nx,i+1,j+1)));
909 48202 : elem->set_node(4, mesh.node_ptr(idx(type,nx,i+1,j+2)));
910 48202 : elem->set_node(5, mesh.node_ptr(idx(type,nx,i,j+1) ));
911 :
912 48202 : if (type == TRI7)
913 21095 : elem->set_node(6, mesh.node_ptr(elem->id()+(2*nx+1)*(2*ny+1)));
914 :
915 48202 : if (j == 2*(ny-1))
916 9390 : boundary_info.add_side(elem, 1, 2);
917 :
918 48202 : if (i == 0)
919 9287 : 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 84 : const auto dx_tri = Real(1) / nx;
933 84 : const auto dy_tri = Real(1) / (ny + 1);
934 60 : std::unique_ptr<Elem> new_elem;
935 686 : for (const auto i : make_range(nx + 1))
936 : {
937 : // Make new nodes for bottom layer of triangles
938 : Node *node0, *node1, *node2;
939 602 : if (i == 0)
940 : {
941 108 : node0 = mesh.add_point(Point(0., 0, 0.));
942 108 : node1 = mesh.add_point(Point(0., dy_tri / 2., 0.));
943 108 : node2 = mesh.add_point(Point(dx_tri / 2., 0., 0.));
944 84 : node_list.push_back(node0);
945 84 : node_list.push_back(node1);
946 84 : node_list.push_back(node2);
947 : }
948 518 : else if (i < nx)
949 : {
950 434 : node0 = node_list.back();
951 558 : node1 = mesh.add_point(Point((i)*dx_tri, dy_tri / 2., 0.));
952 558 : node2 = mesh.add_point(Point((i + 1. / 2.) * dx_tri, 0., 0.));
953 434 : node_list.push_back(node1);
954 434 : node_list.push_back(node2);
955 : }
956 : else
957 : {
958 84 : node0 = node_list.back();
959 108 : node1 = mesh.add_point(Point((i)*dx_tri, dy_tri / 2., 0.));
960 108 : node2 = mesh.add_point(Point((i)*dx_tri, 0., 0.));
961 84 : node_list.push_back(node1);
962 84 : node_list.push_back(node2);
963 : }
964 :
965 602 : new_elem = std::make_unique<C0Polygon>(3);
966 : // Switch to Tri3 when exodus default output supports element type mixes
967 602 : new_elem->set_node(0, node0);
968 602 : new_elem->set_node(1, node1);
969 602 : new_elem->set_node(2, node2);
970 774 : auto * elem = mesh.add_elem(std::move(new_elem));
971 :
972 : // Set boundaries
973 602 : if (i == 0)
974 84 : boundary_info.add_side(elem, 0, 3); // left
975 518 : else if (i == nx)
976 84 : boundary_info.add_side(elem, 1, 1); // right
977 602 : 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 84 : (Real(1) - (ny == 1 ?
985 : dy_tri :
986 84 : (Real(1) + (ny - 1) / 2.) * dy_tri)) / ny;
987 602 : for (const auto j : make_range(ny))
988 : {
989 4851 : for (const auto i : make_range(nx + (j % 2)))
990 : {
991 4333 : if ((j % 2 == 0) || ((i > 0) && (i < nx)))
992 : {
993 : Node *n0, *n1, *n2, *n3, *n4, *n5;
994 3843 : n0 = node_list[running_index++];
995 3843 : n1 = node_list[running_index++];
996 3843 : n2 = node_list[running_index];
997 :
998 3843 : if (i == 0)
999 : {
1000 351 : n3 = mesh.add_point(Point(*n0) + RealVectorValue(0, hex_side, 0));
1001 273 : node_list.push_back(n3);
1002 : }
1003 : else
1004 3570 : n3 = node_list.back();
1005 :
1006 4941 : n4 = mesh.add_point(Point(*n1) + RealVectorValue(0, hex_side + dy_tri, 0));
1007 4941 : n5 = mesh.add_point(Point(*n2) + RealVectorValue(0, hex_side, 0));
1008 3843 : node_list.push_back(n4);
1009 3843 : node_list.push_back(n5);
1010 :
1011 3843 : new_elem = std::make_unique<libMesh::C0Polygon>(6);
1012 3843 : new_elem->set_node(0, n0);
1013 3843 : new_elem->set_node(1, n1);
1014 3843 : new_elem->set_node(2, n2);
1015 3843 : new_elem->set_node(3, n5);
1016 3843 : new_elem->set_node(4, n4);
1017 3843 : new_elem->set_node(5, n3);
1018 4941 : auto * elem = mesh.add_elem(std::move(new_elem));
1019 :
1020 : // Set boundaries
1021 3843 : if (i == 0)
1022 273 : boundary_info.add_side(elem, 5, 3); // left
1023 3570 : else if (i == nx)
1024 2196 : boundary_info.add_side(elem, 2, 1); // right
1025 1647 : }
1026 : // The hexagons are offset, so we build on a quad on each external side to fill
1027 490 : else if (i == 0 || i == nx)
1028 : {
1029 : Node *n0, *n1, *n2, *n3;
1030 490 : n0 = node_list[running_index++];
1031 490 : n1 = node_list[running_index];
1032 :
1033 490 : if (i == 0)
1034 : {
1035 315 : n2 = mesh.add_point(Point(*n0) + RealVectorValue(0, hex_side + dy_tri, 0));
1036 245 : node_list.push_back(n2);
1037 315 : n3 = mesh.add_point(Point(*n1) + RealVectorValue(0, hex_side, 0));
1038 : }
1039 : else
1040 : {
1041 245 : n2 = node_list.back();
1042 315 : n3 = mesh.add_point(Point(*n1) + RealVectorValue(0, hex_side + dy_tri, 0));
1043 : }
1044 490 : node_list.push_back(n3);
1045 :
1046 490 : new_elem = std::make_unique<C0Polygon>(4);
1047 : // Switch to Quad4 when exodus default output supports element type mixes
1048 490 : new_elem->set_node(0, n0);
1049 490 : new_elem->set_node(1, n1);
1050 490 : new_elem->set_node(3, n2);
1051 490 : new_elem->set_node(2, n3);
1052 630 : auto * elem = mesh.add_elem(std::move(new_elem));
1053 :
1054 : // Set boundaries
1055 490 : if (i == 0)
1056 245 : boundary_info.add_side(elem, 3, 3); // left
1057 245 : else if (i == nx)
1058 385 : 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 518 : running_index++;
1065 :
1066 : // Skip lower right corner node
1067 518 : if (j == 0)
1068 84 : running_index++;
1069 : }
1070 :
1071 : // Build a final layer of triangles
1072 84 : const bool ny_odd = (ny % 2 == 1);
1073 630 : for (const auto i : make_range(nx + ny_odd))
1074 : {
1075 : // Use existing nodes, except at the corners
1076 : Node *node0, *node1, *node2;
1077 546 : if (i == 0 && ny_odd)
1078 : {
1079 36 : node0 = mesh.add_point(Point(0., 1., 0.));
1080 28 : node1 = node_list[running_index++];
1081 36 : node2 = node_list[running_index];
1082 : }
1083 518 : else if (i < nx)
1084 : {
1085 490 : node0 = node_list[running_index++];
1086 490 : node1 = node_list[running_index++];
1087 630 : 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 28 : node0 = node_list[running_index++];
1093 28 : node1 = node_list[running_index];
1094 36 : node2 = mesh.add_point(Point(1., 1., 0.));
1095 : }
1096 :
1097 546 : new_elem = std::make_unique<C0Polygon>(3);
1098 : // Switch to Tri3 when exodus default output supports element type mixes
1099 546 : new_elem->set_node(0, node0);
1100 546 : new_elem->set_node(1, node1);
1101 546 : new_elem->set_node(2, node2);
1102 702 : auto * elem = mesh.add_elem(std::move(new_elem));
1103 :
1104 : // Set boundaries
1105 546 : if (i == 0)
1106 84 : boundary_info.add_side(elem, 0, 3); // left
1107 462 : else if (i == nx)
1108 28 : boundary_info.add_side(elem, 1, 1); // right
1109 546 : boundary_info.add_side(elem, 2, 2); // top
1110 :
1111 : }
1112 24 : break;
1113 36 : }
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 15838 : 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 15838 : (mesh.node_stored_range(),
1134 1510171 : [xmin, xmax, ymin, ymax](Node * node)
1135 : {
1136 980217 : (*node)(0) = (*node)(0)*(xmax-xmin) + xmin;
1137 766412 : (*node)(1) = ((*node)(1))*(ymax-ymin) + ymin;
1138 552607 : });
1139 : }
1140 :
1141 : // Add sideset names to boundary info
1142 15838 : boundary_info.sideset_name(0) = "bottom";
1143 15838 : boundary_info.sideset_name(1) = "right";
1144 15838 : boundary_info.sideset_name(2) = "top";
1145 15838 : boundary_info.sideset_name(3) = "left";
1146 :
1147 : // Add nodeset names to boundary info
1148 15838 : boundary_info.nodeset_name(0) = "bottom";
1149 15838 : boundary_info.nodeset_name(1) = "right";
1150 15838 : boundary_info.nodeset_name(2) = "top";
1151 15838 : boundary_info.nodeset_name(3) = "left";
1152 :
1153 4535 : 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 18552 : 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 12032 : 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 12032 : mesh.reserve_elem(nx*ny*nz);
1196 8586 : break;
1197 : }
1198 :
1199 6520 : case PRISM6:
1200 : case PRISM15:
1201 : case PRISM18:
1202 : case PRISM20:
1203 : case PRISM21:
1204 : {
1205 6520 : mesh.reserve_elem(2*nx*ny*nz);
1206 4656 : 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 3242 : case INVALID_ELEM:
1221 : case HEX8:
1222 : case PRISM6:
1223 : case C0POLYHEDRON:
1224 : {
1225 : const dof_id_type grid_nodes =
1226 3242 : 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 3242 : (type == C0POLYHEDRON) ?
1233 35 : cast_int<dof_id_type>(nx*ny*nz) : 0;
1234 :
1235 3242 : mesh.reserve_nodes(grid_nodes + mid_polyhedron_nodes);
1236 2306 : break;
1237 : }
1238 :
1239 10686 : 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 10686 : mesh.reserve_nodes( (2*nx+1)*(2*ny+1)*(2*nz+1) );
1255 7632 : break;
1256 : }
1257 :
1258 1621 : case TET14:
1259 : {
1260 3931 : mesh.reserve_nodes( (2*nx+1)*(2*ny+1)*(2*nz+1) +
1261 1621 : 24*nx*ny*nz +
1262 1621 : 4*(nx*ny + ny*nz + nx*nz) );
1263 1159 : break;
1264 : }
1265 :
1266 1344 : case PRISM20:
1267 : {
1268 2496 : mesh.reserve_nodes( (2*nx+1)*(2*ny+1)*(2*nz+1) +
1269 1344 : 2*nx*ny*(nz+1) );
1270 960 : break;
1271 : }
1272 :
1273 1659 : case PRISM21:
1274 : {
1275 3081 : mesh.reserve_nodes( (2*nx+1)*(2*ny+1)*(2*nz+1) +
1276 1659 : 2*nx*ny*(2*nz+1) );
1277 1185 : 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 11339 : for (unsigned int k=0; k<=nz; k++)
1297 33646 : for (unsigned int j=0; j<=ny; j++)
1298 197614 : for (unsigned int i=0; i<=nx; i++)
1299 : {
1300 : const Node * const node =
1301 243238 : mesh.add_point(Point(static_cast<Real>(i) / static_cast<Real>(nx),
1302 172065 : static_cast<Real>(j) / static_cast<Real>(ny),
1303 172065 : static_cast<Real>(k) / static_cast<Real>(nz)),
1304 150864 : node_id++);
1305 172065 : if (k == 0)
1306 34383 : boundary_info.add_node(node, 0);
1307 172065 : if (k == nz)
1308 34383 : boundary_info.add_node(node, 5);
1309 172065 : if (j == 0)
1310 29763 : boundary_info.add_node(node, 1);
1311 172065 : if (j == ny)
1312 29763 : boundary_info.add_node(node, 3);
1313 172065 : if (i == 0)
1314 25549 : boundary_info.add_node(node, 4);
1315 172065 : if (i == nx)
1316 25549 : 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 71418 : for (unsigned int k=0; k<=(2*nz); k++)
1337 294064 : for (unsigned int j=0; j<=(2*ny); j++)
1338 1809430 : for (unsigned int i=0; i<=(2*nx); i++)
1339 : {
1340 : const Node * const node =
1341 2281336 : mesh.add_point(Point(static_cast<Real>(i) / static_cast<Real>(2 * nx),
1342 1571474 : static_cast<Real>(j) / static_cast<Real>(2 * ny),
1343 1571474 : static_cast<Real>(k) / static_cast<Real>(2 * nz)),
1344 1356596 : node_id++);
1345 1571474 : if (k == 0)
1346 249280 : boundary_info.add_node(node, 0);
1347 1571474 : if (k == 2*nz)
1348 249280 : boundary_info.add_node(node, 5);
1349 1571474 : if (j == 0)
1350 243810 : boundary_info.add_node(node, 1);
1351 1571474 : if (j == 2*ny)
1352 243810 : boundary_info.add_node(node, 3);
1353 1571474 : if (i == 0)
1354 237956 : boundary_info.add_node(node, 4);
1355 1571474 : if (i == 2*nx)
1356 237956 : boundary_info.add_node(node, 2);
1357 : }
1358 :
1359 15310 : if (type == PRISM20 ||
1360 : type == PRISM21)
1361 : {
1362 3003 : const unsigned int kmax = (type == PRISM20) ? nz : 2*nz;
1363 11963 : for (unsigned int k=0; k<=kmax; k++)
1364 21364 : for (unsigned int j=0; j<ny; j++)
1365 31696 : for (unsigned int i=0; i<nx; i++)
1366 : {
1367 : const Node * const node1 =
1368 27560 : mesh.add_point(Point((static_cast<Real>(i)+1/Real(3)) / static_cast<Real>(nx),
1369 19292 : (static_cast<Real>(j)+1/Real(3)) / static_cast<Real>(ny),
1370 19292 : static_cast<Real>(k) / static_cast<Real>(kmax)),
1371 16536 : node_id++);
1372 19292 : if (k == 0)
1373 5565 : boundary_info.add_node(node1, 0);
1374 19292 : if (k == kmax)
1375 5565 : boundary_info.add_node(node1, 5);
1376 :
1377 : const Node * const node2 =
1378 27560 : mesh.add_point(Point((static_cast<Real>(i)+2/Real(3)) / static_cast<Real>(nx),
1379 19292 : (static_cast<Real>(j)+2/Real(3)) / static_cast<Real>(ny),
1380 5512 : static_cast<Real>(k) / static_cast<Real>(kmax)),
1381 16536 : node_id++);
1382 19292 : if (k == 0)
1383 5565 : boundary_info.add_node(node2, 0);
1384 19292 : if (k == kmax)
1385 5565 : 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 5507 : for (unsigned int k=0; k<nz; k++)
1408 13616 : for (unsigned int j=0; j<ny; j++)
1409 110619 : for (unsigned int i=0; i<nx; i++)
1410 : {
1411 100395 : Elem * elem = mesh.add_elem(Elem::build_with_id(HEX8, elem_id++));
1412 100395 : elem->set_node(0, mesh.node_ptr(idx(type,nx,ny,i,j,k) ));
1413 100395 : elem->set_node(1, mesh.node_ptr(idx(type,nx,ny,i+1,j,k) ));
1414 100395 : elem->set_node(2, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) ));
1415 100395 : elem->set_node(3, mesh.node_ptr(idx(type,nx,ny,i,j+1,k) ));
1416 100395 : elem->set_node(4, mesh.node_ptr(idx(type,nx,ny,i,j,k+1) ));
1417 100395 : elem->set_node(5, mesh.node_ptr(idx(type,nx,ny,i+1,j,k+1) ));
1418 100395 : elem->set_node(6, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1)));
1419 100395 : elem->set_node(7, mesh.node_ptr(idx(type,nx,ny,i,j+1,k+1) ));
1420 :
1421 100395 : if (k == 0)
1422 18071 : boundary_info.add_side(elem, 0, 0);
1423 :
1424 100395 : if (k == (nz-1))
1425 18071 : boundary_info.add_side(elem, 5, 5);
1426 :
1427 100395 : if (j == 0)
1428 13591 : boundary_info.add_side(elem, 1, 1);
1429 :
1430 100395 : if (j == (ny-1))
1431 13591 : boundary_info.add_side(elem, 3, 3);
1432 :
1433 100395 : if (i == 0)
1434 10224 : boundary_info.add_side(elem, 4, 4);
1435 :
1436 100395 : if (i == (nx-1))
1437 10224 : boundary_info.add_side(elem, 2, 2);
1438 : }
1439 614 : break;
1440 : }
1441 :
1442 :
1443 35 : case C0POLYHEDRON:
1444 : {
1445 35 : 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 105 : for (unsigned int k=0; k<nz; k++)
1454 210 : for (unsigned int j=0; j<ny; j++)
1455 420 : for (unsigned int i=0; i<nx; i++)
1456 : {
1457 : std::array<Node *, 8> elem_nodes =
1458 280 : {{mesh.node_ptr(idx(type,nx,ny,i,j,k) ),
1459 280 : mesh.node_ptr(idx(type,nx,ny,i+1,j,k) ),
1460 280 : mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) ),
1461 280 : mesh.node_ptr(idx(type,nx,ny,i,j+1,k) ),
1462 280 : mesh.node_ptr(idx(type,nx,ny,i,j,k+1) ),
1463 280 : mesh.node_ptr(idx(type,nx,ny,i+1,j,k+1) ),
1464 280 : mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1)),
1465 1960 : mesh.node_ptr(idx(type,nx,ny,i,j+1,k+1) )}};
1466 :
1467 440 : std::vector<std::shared_ptr<Polygon>> sides(side_nodes.size());
1468 1960 : for (auto s : index_range(side_nodes))
1469 : {
1470 2160 : sides[s] = std::make_shared<C0Polygon>(side_nodes[s].size());
1471 8400 : for (auto n : index_range(side_nodes[s]))
1472 10560 : sides[s]->set_node(n, elem_nodes[side_nodes[s][n]]);
1473 : }
1474 :
1475 280 : std::unique_ptr<Node> mid_elem_node;
1476 : std::unique_ptr<Elem> new_elem =
1477 360 : std::make_unique<C0Polyhedron>(sides, mid_elem_node);
1478 280 : if (mid_elem_node)
1479 0 : mesh.add_node(std::move(mid_elem_node));
1480 :
1481 280 : new_elem->set_id() = elem_id++;
1482 360 : Elem * elem = mesh.add_elem(std::move(new_elem));
1483 :
1484 280 : if (k == 0)
1485 140 : boundary_info.add_side(elem, 0, 0);
1486 :
1487 280 : if (k == (nz-1))
1488 140 : boundary_info.add_side(elem, 5, 5);
1489 :
1490 280 : if (j == 0)
1491 140 : boundary_info.add_side(elem, 1, 1);
1492 :
1493 280 : if (j == (ny-1))
1494 140 : boundary_info.add_side(elem, 3, 3);
1495 :
1496 280 : if (i == 0)
1497 140 : boundary_info.add_side(elem, 4, 4);
1498 :
1499 280 : if (i == (nx-1))
1500 140 : boundary_info.add_side(elem, 2, 2);
1501 120 : }
1502 10 : break;
1503 : }
1504 :
1505 :
1506 :
1507 :
1508 312 : case PRISM6:
1509 : {
1510 2485 : for (unsigned int k=0; k<nz; k++)
1511 3486 : for (unsigned int j=0; j<ny; j++)
1512 5964 : for (unsigned int i=0; i<nx; i++)
1513 : {
1514 : // First Prism
1515 3871 : Elem * elem = mesh.add_elem(Elem::build_with_id(PRISM6, elem_id++));
1516 3871 : elem->set_node(0, mesh.node_ptr(idx(type,nx,ny,i,j,k) ));
1517 3871 : elem->set_node(1, mesh.node_ptr(idx(type,nx,ny,i+1,j,k) ));
1518 3871 : elem->set_node(2, mesh.node_ptr(idx(type,nx,ny,i,j+1,k) ));
1519 3871 : elem->set_node(3, mesh.node_ptr(idx(type,nx,ny,i,j,k+1) ));
1520 3871 : elem->set_node(4, mesh.node_ptr(idx(type,nx,ny,i+1,j,k+1) ));
1521 3871 : 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 3871 : if (i==0)
1525 2093 : boundary_info.add_side(elem, 3, 4);
1526 :
1527 3871 : if (j==0)
1528 2093 : boundary_info.add_side(elem, 1, 1);
1529 :
1530 3871 : if (k==0)
1531 2093 : boundary_info.add_side(elem, 0, 0);
1532 :
1533 3871 : if (k == (nz-1))
1534 2093 : boundary_info.add_side(elem, 4, 5);
1535 :
1536 : // Second Prism
1537 3871 : elem = mesh.add_elem(Elem::build_with_id(PRISM6, elem_id++));
1538 3871 : elem->set_node(0, mesh.node_ptr(idx(type,nx,ny,i+1,j,k) ));
1539 3871 : elem->set_node(1, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) ));
1540 3871 : elem->set_node(2, mesh.node_ptr(idx(type,nx,ny,i,j+1,k) ));
1541 3871 : elem->set_node(3, mesh.node_ptr(idx(type,nx,ny,i+1,j,k+1) ));
1542 3871 : elem->set_node(4, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1)));
1543 3871 : 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 3871 : if (i == (nx-1))
1547 2093 : boundary_info.add_side(elem, 1, 2);
1548 :
1549 3871 : if (j == (ny-1))
1550 2093 : boundary_info.add_side(elem, 2, 3);
1551 :
1552 3871 : if (k==0)
1553 2093 : boundary_info.add_side(elem, 0, 0);
1554 :
1555 3871 : if (k == (nz-1))
1556 2093 : 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 23487 : for (unsigned int k=0; k<(2*nz); k += 2)
1577 39283 : for (unsigned int j=0; j<(2*ny); j += 2)
1578 134698 : for (unsigned int i=0; i<(2*nx); i += 2)
1579 : {
1580 109020 : ElemType build_type = (type == HEX20) ? HEX20 : HEX27;
1581 109020 : Elem * elem = mesh.add_elem(Elem::build_with_id(build_type, elem_id++));
1582 :
1583 109020 : elem->set_node(0, mesh.node_ptr(idx(type,nx,ny,i, j, k) ));
1584 109020 : elem->set_node(1, mesh.node_ptr(idx(type,nx,ny,i+2,j, k) ));
1585 109020 : elem->set_node(2, mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k) ));
1586 109020 : elem->set_node(3, mesh.node_ptr(idx(type,nx,ny,i, j+2,k) ));
1587 109020 : elem->set_node(4, mesh.node_ptr(idx(type,nx,ny,i, j, k+2)));
1588 109020 : elem->set_node(5, mesh.node_ptr(idx(type,nx,ny,i+2,j, k+2)));
1589 109020 : elem->set_node(6, mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k+2)));
1590 109020 : elem->set_node(7, mesh.node_ptr(idx(type,nx,ny,i, j+2,k+2)));
1591 109020 : elem->set_node(8, mesh.node_ptr(idx(type,nx,ny,i+1,j, k) ));
1592 109020 : elem->set_node(9, mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k) ));
1593 109020 : elem->set_node(10, mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k) ));
1594 109020 : elem->set_node(11, mesh.node_ptr(idx(type,nx,ny,i, j+1,k) ));
1595 109020 : elem->set_node(12, mesh.node_ptr(idx(type,nx,ny,i, j, k+1)));
1596 109020 : elem->set_node(13, mesh.node_ptr(idx(type,nx,ny,i+2,j, k+1)));
1597 109020 : elem->set_node(14, mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k+1)));
1598 109020 : elem->set_node(15, mesh.node_ptr(idx(type,nx,ny,i, j+2,k+1)));
1599 109020 : elem->set_node(16, mesh.node_ptr(idx(type,nx,ny,i+1,j, k+2)));
1600 109020 : elem->set_node(17, mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k+2)));
1601 109020 : elem->set_node(18, mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k+2)));
1602 109020 : elem->set_node(19, mesh.node_ptr(idx(type,nx,ny,i, j+1,k+2)));
1603 :
1604 109020 : if ((type == HEX27) || (type == TET4) || (type == TET10) || (type == TET14) ||
1605 6884 : (type == PYRAMID5) || (type == PYRAMID13) || (type == PYRAMID14) ||
1606 2418 : (type == PYRAMID18))
1607 : {
1608 104554 : elem->set_node(20, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) ));
1609 104554 : elem->set_node(21, mesh.node_ptr(idx(type,nx,ny,i+1,j, k+1)));
1610 104554 : elem->set_node(22, mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k+1)));
1611 104554 : elem->set_node(23, mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k+1)));
1612 104554 : elem->set_node(24, mesh.node_ptr(idx(type,nx,ny,i, j+1,k+1)));
1613 104554 : elem->set_node(25, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+2)));
1614 104554 : elem->set_node(26, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1)));
1615 : }
1616 :
1617 109020 : if (k == 0)
1618 28218 : boundary_info.add_side(elem, 0, 0);
1619 :
1620 109020 : if (k == 2*(nz-1))
1621 28218 : boundary_info.add_side(elem, 5, 5);
1622 :
1623 109020 : if (j == 0)
1624 26919 : boundary_info.add_side(elem, 1, 1);
1625 :
1626 109020 : if (j == 2*(ny-1))
1627 26919 : boundary_info.add_side(elem, 3, 3);
1628 :
1629 109020 : if (i == 0)
1630 25678 : boundary_info.add_side(elem, 4, 4);
1631 :
1632 109020 : if (i == 2*(nx-1))
1633 25678 : 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 12222 : for (unsigned int k=0; k<(2*nz); k += 2)
1647 16360 : for (unsigned int j=0; j<(2*ny); j += 2)
1648 24856 : for (unsigned int i=0; i<(2*nx); i += 2)
1649 : {
1650 : // First Prism
1651 15290 : Elem * elem = mesh.add_elem(Elem::build_with_id(type, elem_id++));
1652 15290 : elem->set_node(0, mesh.node_ptr(idx(type,nx,ny,i, j, k) ));
1653 15290 : elem->set_node(1, mesh.node_ptr(idx(type,nx,ny,i+2,j, k) ));
1654 15290 : elem->set_node(2, mesh.node_ptr(idx(type,nx,ny,i, j+2,k) ));
1655 15290 : elem->set_node(3, mesh.node_ptr(idx(type,nx,ny,i, j, k+2)));
1656 15290 : elem->set_node(4, mesh.node_ptr(idx(type,nx,ny,i+2,j, k+2)));
1657 15290 : elem->set_node(5, mesh.node_ptr(idx(type,nx,ny,i, j+2,k+2)));
1658 15290 : elem->set_node(6, mesh.node_ptr(idx(type,nx,ny,i+1,j, k) ));
1659 15290 : elem->set_node(7, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) ));
1660 15290 : elem->set_node(8, mesh.node_ptr(idx(type,nx,ny,i, j+1,k) ));
1661 15290 : elem->set_node(9, mesh.node_ptr(idx(type,nx,ny,i, j, k+1)));
1662 15290 : elem->set_node(10, mesh.node_ptr(idx(type,nx,ny,i+2,j, k+1)));
1663 15290 : elem->set_node(11, mesh.node_ptr(idx(type,nx,ny,i, j+2,k+1)));
1664 15290 : elem->set_node(12, mesh.node_ptr(idx(type,nx,ny,i+1,j, k+2)));
1665 15290 : elem->set_node(13, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+2)));
1666 15290 : elem->set_node(14, mesh.node_ptr(idx(type,nx,ny,i, j+1,k+2)));
1667 :
1668 19706 : if (type == PRISM18 ||
1669 13050 : type == PRISM20 ||
1670 : type == PRISM21)
1671 : {
1672 12420 : elem->set_node(15, mesh.node_ptr(idx(type,nx,ny,i+1,j, k+1)));
1673 12420 : elem->set_node(16, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1)));
1674 12420 : elem->set_node(17, mesh.node_ptr(idx(type,nx,ny,i, j+1,k+1)));
1675 : }
1676 :
1677 15290 : if (type == PRISM20)
1678 : {
1679 4235 : const dof_id_type base_idx = (2*nx+1)*(2*ny+1)*(2*nz+1);
1680 4235 : elem->set_node(18, mesh.node_ptr(base_idx+((k/2)*(nx*ny)+j/2*nx+i/2)*2));
1681 4235 : elem->set_node(19, mesh.node_ptr(base_idx+(((k/2)+1)*(nx*ny)+j/2*nx+i/2)*2));
1682 : }
1683 :
1684 15290 : if (type == PRISM21)
1685 : {
1686 4746 : const dof_id_type base_idx = (2*nx+1)*(2*ny+1)*(2*nz+1);
1687 4746 : elem->set_node(18, mesh.node_ptr(base_idx+(k*(nx*ny)+j/2*nx+i/2)*2));
1688 4746 : elem->set_node(19, mesh.node_ptr(base_idx+((k+2)*(nx*ny)+j/2*nx+i/2)*2));
1689 4746 : 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 15290 : if (i==0)
1694 9566 : boundary_info.add_side(elem, 3, 4);
1695 :
1696 15290 : if (j==0)
1697 9566 : boundary_info.add_side(elem, 1, 1);
1698 :
1699 15290 : if (k==0)
1700 9616 : boundary_info.add_side(elem, 0, 0);
1701 :
1702 15290 : if (k == 2*(nz-1))
1703 9616 : boundary_info.add_side(elem, 4, 5);
1704 :
1705 :
1706 : // Second Prism
1707 15290 : elem = mesh.add_elem(Elem::build_with_id(type, elem_id++));
1708 15290 : elem->set_node(0, mesh.node_ptr(idx(type,nx,ny,i+2,j,k) ));
1709 15290 : elem->set_node(1, mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k) ));
1710 15290 : elem->set_node(2, mesh.node_ptr(idx(type,nx,ny,i,j+2,k) ));
1711 15290 : elem->set_node(3, mesh.node_ptr(idx(type,nx,ny,i+2,j,k+2) ));
1712 15290 : elem->set_node(4, mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k+2) ));
1713 15290 : elem->set_node(5, mesh.node_ptr(idx(type,nx,ny,i,j+2,k+2) ));
1714 15290 : elem->set_node(6, mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k) ));
1715 15290 : elem->set_node(7, mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k) ));
1716 15290 : elem->set_node(8, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) ));
1717 15290 : elem->set_node(9, mesh.node_ptr(idx(type,nx,ny,i+2,j,k+1) ));
1718 15290 : elem->set_node(10, mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k+1)));
1719 15290 : elem->set_node(11, mesh.node_ptr(idx(type,nx,ny,i,j+2,k+1) ));
1720 15290 : elem->set_node(12, mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k+2)));
1721 15290 : elem->set_node(13, mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k+2)));
1722 15290 : elem->set_node(14, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+2)));
1723 :
1724 15290 : if (type == PRISM18 ||
1725 7616 : type == PRISM20 ||
1726 : type == PRISM21)
1727 : {
1728 12420 : elem->set_node(15, mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k+1)));
1729 12420 : elem->set_node(16, mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k+1)));
1730 12420 : elem->set_node(17, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1)));
1731 : }
1732 :
1733 15290 : if (type == PRISM20)
1734 : {
1735 4235 : const dof_id_type base_idx = (2*nx+1)*(2*ny+1)*(2*nz+1);
1736 4235 : elem->set_node(18, mesh.node_ptr(base_idx+((k/2)*(nx*ny)+j/2*nx+i/2)*2+1));
1737 4235 : elem->set_node(19, mesh.node_ptr(base_idx+(((k/2)+1)*(nx*ny)+j/2*nx+i/2)*2+1));
1738 : }
1739 :
1740 15290 : if (type == PRISM21)
1741 : {
1742 4746 : const dof_id_type base_idx = (2*nx+1)*(2*ny+1)*(2*nz+1);
1743 4746 : elem->set_node(18, mesh.node_ptr(base_idx+(k*(nx*ny)+j/2*nx+i/2)*2+1));
1744 4746 : elem->set_node(19, mesh.node_ptr(base_idx+((k+2)*(nx*ny)+j/2*nx+i/2)*2+1));
1745 4746 : 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 15290 : if (i == 2*(nx-1))
1750 9566 : boundary_info.add_side(elem, 1, 2);
1751 :
1752 15290 : if (j == 2*(ny-1))
1753 9566 : boundary_info.add_side(elem, 2, 3);
1754 :
1755 15290 : if (k==0)
1756 9616 : boundary_info.add_side(elem, 0, 0);
1757 :
1758 15290 : if (k == 2*(nz-1))
1759 9616 : 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 18552 : 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 18552 : (mesh.node_stored_range(),
1789 4938624 : [xmin, xmax, ymin, ymax, zmin, zmax](Node * node)
1790 : {
1791 2274399 : (*node)(0) = ((*node)(0))*(xmax-xmin) + xmin;
1792 2274399 : (*node)(1) = ((*node)(1))*(ymax-ymin) + ymin;
1793 1782123 : (*node)(2) = ((*node)(2))*(zmax-zmin) + zmin;
1794 1289847 : });
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 23862 : if ((type == TET4) ||
1807 17864 : (type == TET10) ||
1808 17402 : (type == TET14) ||
1809 14423 : (type == PYRAMID5) ||
1810 3960 : (type == PYRAMID13) ||
1811 17685 : (type == PYRAMID14) ||
1812 9865 : (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 5490 : std::unique_ptr<Elem> side;
1819 :
1820 5490 : if ((type == TET4) || (type == TET10) || (type == TET14))
1821 4027 : new_elements.reserve(24*mesh.n_elem());
1822 : else
1823 1463 : new_elements.reserve(6*mesh.n_elem());
1824 :
1825 : // Create tetrahedra or pyramids
1826 47852 : for (auto & base_hex : mesh.element_ptr_range())
1827 : {
1828 : // Get a pointer to the node located at the HEX27 center
1829 26708 : 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 186956 : for (auto s : base_hex->side_index_range())
1835 : {
1836 : // Get the boundary ID(s) for this side
1837 160248 : 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 160248 : boundary_id_type b_id = ids.empty() ? BoundaryInfo::invalid_id : ids[0];
1844 :
1845 : // Need to build the full-ordered side!
1846 160248 : base_hex->build_side_ptr(side, s);
1847 :
1848 160248 : if ((type == TET4) || (type == TET10) || (type == TET14))
1849 : {
1850 : // Build 4 sub-tets per side
1851 529500 : for (unsigned int sub_tet=0; sub_tet<4; ++sub_tet)
1852 : {
1853 729600 : new_elements.push_back( Elem::build(TET4) );
1854 117600 : auto & sub_elem = new_elements.back();
1855 658800 : sub_elem->set_node(0, side->node_ptr(sub_tet));
1856 658800 : sub_elem->set_node(1, side->node_ptr(8)); // center of the face
1857 541200 : sub_elem->set_node(2, side->node_ptr(sub_tet==3 ? 0 : sub_tet+1 )); // wrap-around
1858 423600 : 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 423600 : if (b_id != BoundaryInfo::invalid_id)
1864 256376 : 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 93168 : 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 85404 : sub_elem->set_node(0, side->node_ptr(0));
1878 85404 : sub_elem->set_node(1, side->node_ptr(3));
1879 85404 : sub_elem->set_node(2, side->node_ptr(2));
1880 85404 : sub_elem->set_node(3, side->node_ptr(1));
1881 :
1882 : // Set the apex
1883 54348 : 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 54348 : if (b_id != BoundaryInfo::invalid_id)
1888 35100 : boundary_info.add_side(sub_elem.get(), 4, b_id);
1889 : } // end else type==PYRAMID*
1890 : }
1891 2354 : }
1892 :
1893 :
1894 : // Delete the original HEX27 elements from the mesh, and the boundary info structure.
1895 47852 : for (auto & elem : mesh.element_ptr_range())
1896 : {
1897 26708 : boundary_info.remove(elem); // Safe even if elem has no boundary info.
1898 26708 : mesh.delete_elem(elem);
1899 2354 : }
1900 :
1901 : // Add the new elements
1902 483438 : for (auto i : index_range(new_elements))
1903 : {
1904 465948 : new_elements[i]->set_id(i);
1905 744204 : mesh.add_elem( std::move(new_elements[i]) );
1906 : }
1907 :
1908 2354 : } // 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 18552 : if ((type == TET10) || (type == PYRAMID14))
1913 1580 : mesh.all_second_order();
1914 :
1915 16972 : else if (type == PYRAMID13)
1916 343 : mesh.all_second_order(/*full_ordered=*/false);
1917 :
1918 16629 : else if ((type == TET14) || (type == PYRAMID18))
1919 2034 : mesh.all_complete_order();
1920 :
1921 :
1922 : // Add sideset names to boundary info (Z axis out of the screen)
1923 18552 : boundary_info.sideset_name(0) = "back";
1924 18552 : boundary_info.sideset_name(1) = "bottom";
1925 18552 : boundary_info.sideset_name(2) = "right";
1926 18552 : boundary_info.sideset_name(3) = "top";
1927 18552 : boundary_info.sideset_name(4) = "left";
1928 18552 : boundary_info.sideset_name(5) = "front";
1929 :
1930 : // Add nodeset names to boundary info
1931 18552 : boundary_info.nodeset_name(0) = "back";
1932 18552 : boundary_info.nodeset_name(1) = "bottom";
1933 18552 : boundary_info.nodeset_name(2) = "right";
1934 18552 : boundary_info.nodeset_name(3) = "top";
1935 18552 : boundary_info.nodeset_name(4) = "left";
1936 18552 : 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 40530 : mesh.prepare_for_use ();
1947 40530 : }
1948 :
1949 :
1950 :
1951 116 : 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 116 : build_cube(mesh,
1960 : 0, 0, 0,
1961 : 0., 0.,
1962 : 0., 0.,
1963 : 0., 0.,
1964 : type,
1965 : gauss_lobatto_grid);
1966 116 : }
1967 :
1968 :
1969 579 : 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 579 : build_cube(mesh,
1980 : nx, 0, 0,
1981 : xmin, xmax,
1982 : 0., 0.,
1983 : 0., 0.,
1984 : type,
1985 : gauss_lobatto_grid);
1986 579 : }
1987 :
1988 :
1989 :
1990 2280 : 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 2280 : build_cube (mesh,
2004 : nx, ny, 0,
2005 : xmin, xmax,
2006 : ymin, ymax,
2007 : 0., 0.,
2008 : type,
2009 : gauss_lobatto_grid);
2010 2280 : }
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 233 : 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 233 : cast_int<unsigned char>(mesh.mesh_dimension());
2051 233 : mesh.clear();
2052 233 : 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 233 : if (mesh.mesh_dimension() == 1)
2061 : {
2062 233 : switch (type)
2063 : {
2064 89 : case HEX8:
2065 : case HEX27:
2066 : case TET4:
2067 : case TET10:
2068 : case TET14:
2069 89 : mesh.set_mesh_dimension(3);
2070 26 : break;
2071 116 : 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 116 : mesh.set_mesh_dimension(2);
2081 34 : break;
2082 28 : case EDGE2:
2083 : case EDGE3:
2084 : case EDGE4:
2085 28 : 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 233 : const bool is_replicated = mesh.is_replicated();
2099 :
2100 : // Sphere is centered at origin by default
2101 68 : const Point cent;
2102 :
2103 466 : const Sphere sphere (cent, rad);
2104 :
2105 233 : switch (mesh.mesh_dimension())
2106 : {
2107 : //-----------------------------------------------------------------
2108 : // Build a line in one dimension
2109 28 : case 1:
2110 : {
2111 28 : 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 116 : 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 116 : if (flat)
2130 : {
2131 34 : const Real sqrt_2 = std::sqrt(2.);
2132 116 : const Real rad_2 = .25*rad;
2133 116 : const Real rad_sqrt_2 = rad/sqrt_2;
2134 :
2135 : // (Temporary) convenient storage for node pointers
2136 150 : std::vector<Node *> nodes(8);
2137 :
2138 : // Point 0
2139 150 : nodes[0] = mesh.add_point (Point(-rad_2,-rad_2, 0.), node_id++);
2140 :
2141 : // Point 1
2142 150 : nodes[1] = mesh.add_point (Point( rad_2,-rad_2, 0.), node_id++);
2143 :
2144 : // Point 2
2145 150 : nodes[2] = mesh.add_point (Point( rad_2, rad_2, 0.), node_id++);
2146 :
2147 : // Point 3
2148 150 : nodes[3] = mesh.add_point (Point(-rad_2, rad_2, 0.), node_id++);
2149 :
2150 : // Point 4
2151 150 : nodes[4] = mesh.add_point (Point(-rad_sqrt_2,-rad_sqrt_2, 0.), node_id++);
2152 :
2153 : // Point 5
2154 150 : nodes[5] = mesh.add_point (Point( rad_sqrt_2,-rad_sqrt_2, 0.), node_id++);
2155 :
2156 : // Point 6
2157 150 : nodes[6] = mesh.add_point (Point( rad_sqrt_2, rad_sqrt_2, 0.), node_id++);
2158 :
2159 : // Point 7
2160 150 : 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 116 : Elem * elem0 = mesh.add_elem (Elem::build(QUAD4));
2167 116 : elem0->set_node(0, nodes[0]);
2168 116 : elem0->set_node(1, nodes[1]);
2169 116 : elem0->set_node(2, nodes[2]);
2170 116 : elem0->set_node(3, nodes[3]);
2171 : }
2172 :
2173 : // Element 1
2174 : {
2175 116 : Elem * elem1 = mesh.add_elem (Elem::build(QUAD4));
2176 116 : elem1->set_node(0, nodes[4]);
2177 116 : elem1->set_node(1, nodes[0]);
2178 116 : elem1->set_node(2, nodes[3]);
2179 116 : elem1->set_node(3, nodes[7]);
2180 : }
2181 :
2182 : // Element 2
2183 : {
2184 116 : Elem * elem2 = mesh.add_elem (Elem::build(QUAD4));
2185 116 : elem2->set_node(0, nodes[4]);
2186 116 : elem2->set_node(1, nodes[5]);
2187 116 : elem2->set_node(2, nodes[1]);
2188 116 : elem2->set_node(3, nodes[0]);
2189 : }
2190 :
2191 : // Element 3
2192 : {
2193 116 : Elem * elem3 = mesh.add_elem (Elem::build(QUAD4));
2194 116 : elem3->set_node(0, nodes[1]);
2195 116 : elem3->set_node(1, nodes[5]);
2196 116 : elem3->set_node(2, nodes[6]);
2197 116 : elem3->set_node(3, nodes[2]);
2198 : }
2199 :
2200 : // Element 4
2201 : {
2202 116 : Elem * elem4 = mesh.add_elem (Elem::build(QUAD4));
2203 116 : elem4->set_node(0, nodes[3]);
2204 116 : elem4->set_node(1, nodes[2]);
2205 116 : elem4->set_node(2, nodes[6]);
2206 116 : 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 89 : case 3:
2275 : {
2276 : // (Currently) supported types
2277 89 : 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 89 : r_small = 0.25*rad, // 0.25 *radius
2287 89 : r_med = (0.125*std::sqrt(2.)+0.5)*rad; // .67677*radius
2288 :
2289 : // (Temporary) convenient storage for node pointers
2290 115 : 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 115 : nodes[0] = mesh.add_point (Point(-r_small,-r_small, -r_small), node_id++);
2300 115 : nodes[1] = mesh.add_point (Point( r_small,-r_small, -r_small), node_id++);
2301 115 : nodes[2] = mesh.add_point (Point( r_small, r_small, -r_small), node_id++);
2302 115 : nodes[3] = mesh.add_point (Point(-r_small, r_small, -r_small), node_id++);
2303 115 : nodes[4] = mesh.add_point (Point(-r_small,-r_small, r_small), node_id++);
2304 115 : nodes[5] = mesh.add_point (Point( r_small,-r_small, r_small), node_id++);
2305 115 : nodes[6] = mesh.add_point (Point( r_small, r_small, r_small), node_id++);
2306 115 : 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 115 : nodes[8] = mesh.add_point (Point(-r_med,-r_med, -r_med), node_id++);
2310 115 : nodes[9] = mesh.add_point (Point( r_med,-r_med, -r_med), node_id++);
2311 115 : nodes[10] = mesh.add_point (Point( r_med, r_med, -r_med), node_id++);
2312 115 : nodes[11] = mesh.add_point (Point(-r_med, r_med, -r_med), node_id++);
2313 115 : nodes[12] = mesh.add_point (Point(-r_med,-r_med, r_med), node_id++);
2314 115 : nodes[13] = mesh.add_point (Point( r_med,-r_med, r_med), node_id++);
2315 115 : nodes[14] = mesh.add_point (Point( r_med, r_med, r_med), node_id++);
2316 115 : 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 89 : Elem * elem0 = mesh.add_elem(Elem::build(HEX8));
2322 89 : elem0->set_node(0, nodes[0]);
2323 89 : elem0->set_node(1, nodes[1]);
2324 89 : elem0->set_node(2, nodes[2]);
2325 89 : elem0->set_node(3, nodes[3]);
2326 89 : elem0->set_node(4, nodes[4]);
2327 89 : elem0->set_node(5, nodes[5]);
2328 89 : elem0->set_node(6, nodes[6]);
2329 89 : elem0->set_node(7, nodes[7]);
2330 : }
2331 :
2332 : // Element 1 - "bottom"
2333 : {
2334 89 : Elem * elem1 = mesh.add_elem(Elem::build(HEX8));
2335 89 : elem1->set_node(0, nodes[8]);
2336 89 : elem1->set_node(1, nodes[9]);
2337 89 : elem1->set_node(2, nodes[10]);
2338 89 : elem1->set_node(3, nodes[11]);
2339 89 : elem1->set_node(4, nodes[0]);
2340 89 : elem1->set_node(5, nodes[1]);
2341 89 : elem1->set_node(6, nodes[2]);
2342 89 : elem1->set_node(7, nodes[3]);
2343 : }
2344 :
2345 : // Element 2 - "front"
2346 : {
2347 89 : Elem * elem2 = mesh.add_elem(Elem::build(HEX8));
2348 89 : elem2->set_node(0, nodes[8]);
2349 89 : elem2->set_node(1, nodes[9]);
2350 89 : elem2->set_node(2, nodes[1]);
2351 89 : elem2->set_node(3, nodes[0]);
2352 89 : elem2->set_node(4, nodes[12]);
2353 89 : elem2->set_node(5, nodes[13]);
2354 89 : elem2->set_node(6, nodes[5]);
2355 89 : elem2->set_node(7, nodes[4]);
2356 : }
2357 :
2358 : // Element 3 - "right"
2359 : {
2360 89 : Elem * elem3 = mesh.add_elem(Elem::build(HEX8));
2361 89 : elem3->set_node(0, nodes[1]);
2362 89 : elem3->set_node(1, nodes[9]);
2363 89 : elem3->set_node(2, nodes[10]);
2364 89 : elem3->set_node(3, nodes[2]);
2365 89 : elem3->set_node(4, nodes[5]);
2366 89 : elem3->set_node(5, nodes[13]);
2367 89 : elem3->set_node(6, nodes[14]);
2368 89 : elem3->set_node(7, nodes[6]);
2369 : }
2370 :
2371 : // Element 4 - "back"
2372 : {
2373 89 : Elem * elem4 = mesh.add_elem(Elem::build(HEX8));
2374 89 : elem4->set_node(0, nodes[3]);
2375 89 : elem4->set_node(1, nodes[2]);
2376 89 : elem4->set_node(2, nodes[10]);
2377 89 : elem4->set_node(3, nodes[11]);
2378 89 : elem4->set_node(4, nodes[7]);
2379 89 : elem4->set_node(5, nodes[6]);
2380 89 : elem4->set_node(6, nodes[14]);
2381 89 : elem4->set_node(7, nodes[15]);
2382 : }
2383 :
2384 : // Element 5 - "left"
2385 : {
2386 89 : Elem * elem5 = mesh.add_elem(Elem::build(HEX8));
2387 89 : elem5->set_node(0, nodes[8]);
2388 89 : elem5->set_node(1, nodes[0]);
2389 89 : elem5->set_node(2, nodes[3]);
2390 89 : elem5->set_node(3, nodes[11]);
2391 89 : elem5->set_node(4, nodes[12]);
2392 89 : elem5->set_node(5, nodes[4]);
2393 89 : elem5->set_node(6, nodes[7]);
2394 89 : elem5->set_node(7, nodes[15]);
2395 : }
2396 :
2397 : // Element 6 - "top"
2398 : {
2399 89 : Elem * elem6 = mesh.add_elem(Elem::build(HEX8));
2400 89 : elem6->set_node(0, nodes[4]);
2401 89 : elem6->set_node(1, nodes[5]);
2402 89 : elem6->set_node(2, nodes[6]);
2403 89 : elem6->set_node(3, nodes[7]);
2404 89 : elem6->set_node(4, nodes[12]);
2405 89 : elem6->set_node(5, nodes[13]);
2406 89 : elem6->set_node(6, nodes[14]);
2407 89 : 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 466 : MeshRefinement mesh_refinement (mesh);
2424 :
2425 : // For avoiding extraneous element side construction
2426 233 : std::unique_ptr<Elem> side;
2427 :
2428 : // Loop over the elements, refine, pop nodes to boundary.
2429 562 : 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 329 : if (!is_replicated)
2435 84 : mesh.prepare_for_use();
2436 :
2437 329 : mesh_refinement.uniformly_refine(1);
2438 :
2439 : const bool move_only_boundary_nodes =
2440 329 : mesh.mesh_dimension() != 2 || flat;
2441 : MeshTools::Modification::interpolate_surface
2442 658 : (mesh, sphere, /*ids=*/{}, move_only_boundary_nodes);
2443 : }
2444 :
2445 : // A DistributedMesh needs a little prep before flattening
2446 233 : if (!is_replicated)
2447 42 : 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 233 : MeshTools::Modification::flatten(mesh);
2454 :
2455 : // Convert all the tensor product elements to simplices if requested
2456 233 : 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 49 : if (is_replicated)
2461 42 : mesh.prepare_for_use();
2462 :
2463 49 : MeshTools::Modification::all_tri(mesh);
2464 : }
2465 :
2466 : // Convert to second-order elements if the user requested it.
2467 301 : if (Elem::build(type)->default_order() != FIRST)
2468 : {
2469 163 : 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 163 : bool full_ordered = !((type==QUAD8) || (type==HEX20));
2478 163 : mesh.all_second_order(full_ordered);
2479 : }
2480 :
2481 : // And pop to the boundary again...
2482 25924 : for (const auto & elem : mesh.active_element_ptr_range())
2483 124685 : for (auto s : elem->side_index_range())
2484 131158 : if (elem->neighbor_ptr(s) == nullptr)
2485 : {
2486 4178 : elem->build_side_ptr(side, s);
2487 :
2488 : // Pop each point to the sphere boundary
2489 37044 : for (auto n : side->node_index_range())
2490 42688 : side->point(n) =
2491 42688 : sphere.closest_point(side->point(n));
2492 67 : }
2493 : }
2494 :
2495 :
2496 : // The meshes could probably use some smoothing.
2497 233 : if (mesh.mesh_dimension() > 1)
2498 : {
2499 265 : LaplaceMeshSmoother smoother(mesh, n_smooth);
2500 205 : smoother.smooth();
2501 : }
2502 :
2503 : // We'll give the whole sphere surface a boundary id of 0
2504 91940 : for (const auto & elem : mesh.active_element_ptr_range())
2505 376577 : for (auto s : elem->side_index_range())
2506 378678 : if (!elem->neighbor_ptr(s))
2507 8513 : boundary_info.add_side(elem, s, 0);
2508 :
2509 : // Done building the mesh. Now prepare it for use.
2510 233 : mesh.prepare_for_use();
2511 233 : }
2512 :
2513 : #endif // #ifndef LIBMESH_ENABLE_AMR
2514 :
2515 :
2516 : // Meshes the tensor product of a 1D and a 1D-or-2D domain.
2517 49 : 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 49 : if (!cross_section.n_elem())
2526 0 : return;
2527 :
2528 49 : dof_id_type orig_elem = cross_section.n_elem();
2529 49 : dof_id_type orig_nodes = cross_section.n_nodes();
2530 :
2531 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
2532 49 : unique_id_type orig_unique_ids = cross_section.parallel_max_unique_id();
2533 : #endif
2534 :
2535 49 : 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 49 : if (!cross_section.is_serial())
2549 0 : mesh.delete_remote_elements();
2550 :
2551 : // We know a priori how many elements we'll need
2552 49 : mesh.reserve_elem(nz*orig_elem);
2553 :
2554 : // For straightforward meshes we need one or two additional layers per
2555 : // element.
2556 182 : if (cross_section.elements_begin() != cross_section.elements_end() &&
2557 143 : (*cross_section.elements_begin())->default_order() == SECOND)
2558 14 : order = 2;
2559 49 : mesh.comm().max(order);
2560 :
2561 49 : 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 1794 : for (const auto & node : cross_section.node_ptr_range())
2567 : {
2568 6510 : for (unsigned int k=0; k != order*nz+1; ++k)
2569 : {
2570 5313 : const dof_id_type new_node_id = node->id() + k * orig_nodes;
2571 5313 : Node * my_node = mesh.query_node_ptr(new_node_id);
2572 5313 : if (!my_node)
2573 : {
2574 : std::unique_ptr<Node> new_node = Node::build
2575 6831 : (*node + (extrusion_vector * k / nz / order),
2576 1518 : new_node_id);
2577 5313 : 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 5313 : const unique_id_type uid = (k == 0) ?
2584 684 : node->unique_id() :
2585 4116 : orig_unique_ids + (k-1)*(orig_nodes + orig_elem) + node->id();
2586 :
2587 1518 : new_node->set_unique_id(uid);
2588 : #endif
2589 :
2590 5313 : cross_section_boundary_info.boundary_ids(node, ids_to_copy);
2591 5313 : boundary_info.add_node(new_node.get(), ids_to_copy);
2592 :
2593 8349 : mesh.add_node(std::move(new_node));
2594 2277 : }
2595 : }
2596 21 : }
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 49 : 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 49 : cross_section.comm().max(next_side_id);
2608 :
2609 734 : for (const auto & elem : cross_section.element_ptr_range())
2610 : {
2611 455 : const ElemType etype = elem->type();
2612 :
2613 : // build_extrusion currently only works on coarse meshes
2614 130 : libmesh_assert (!elem->parent());
2615 :
2616 1421 : for (unsigned int k=0; k != nz; ++k)
2617 : {
2618 690 : std::unique_ptr<Elem> new_elem;
2619 966 : 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 168 : case TRI3:
2657 : {
2658 240 : new_elem = Elem::build(PRISM6);
2659 216 : new_elem->set_node(0, mesh.node_ptr(elem->node_ptr(0)->id() + (k * orig_nodes)));
2660 216 : new_elem->set_node(1, mesh.node_ptr(elem->node_ptr(1)->id() + (k * orig_nodes)));
2661 216 : new_elem->set_node(2, mesh.node_ptr(elem->node_ptr(2)->id() + (k * orig_nodes)));
2662 216 : new_elem->set_node(3, mesh.node_ptr(elem->node_ptr(0)->id() + ((k+1) * orig_nodes)));
2663 216 : new_elem->set_node(4, mesh.node_ptr(elem->node_ptr(1)->id() + ((k+1) * orig_nodes)));
2664 216 : new_elem->set_node(5, mesh.node_ptr(elem->node_ptr(2)->id() + ((k+1) * orig_nodes)));
2665 :
2666 216 : if (elem->neighbor_ptr(0) == remote_elem)
2667 0 : new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
2668 216 : if (elem->neighbor_ptr(1) == remote_elem)
2669 0 : new_elem->set_neighbor(2, const_cast<RemoteElem *>(remote_elem));
2670 216 : 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 448 : case QUAD4:
2742 : {
2743 640 : new_elem = Elem::build(HEX8);
2744 576 : new_elem->set_node(0, mesh.node_ptr(elem->node_ptr(0)->id() + (k * orig_nodes)));
2745 576 : new_elem->set_node(1, mesh.node_ptr(elem->node_ptr(1)->id() + (k * orig_nodes)));
2746 576 : new_elem->set_node(2, mesh.node_ptr(elem->node_ptr(2)->id() + (k * orig_nodes)));
2747 576 : new_elem->set_node(3, mesh.node_ptr(elem->node_ptr(3)->id() + (k * orig_nodes)));
2748 576 : new_elem->set_node(4, mesh.node_ptr(elem->node_ptr(0)->id() + ((k+1) * orig_nodes)));
2749 576 : new_elem->set_node(5, mesh.node_ptr(elem->node_ptr(1)->id() + ((k+1) * orig_nodes)));
2750 576 : new_elem->set_node(6, mesh.node_ptr(elem->node_ptr(2)->id() + ((k+1) * orig_nodes)));
2751 576 : new_elem->set_node(7, mesh.node_ptr(elem->node_ptr(3)->id() + ((k+1) * orig_nodes)));
2752 :
2753 576 : if (elem->neighbor_ptr(0) == remote_elem)
2754 0 : new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
2755 576 : if (elem->neighbor_ptr(1) == remote_elem)
2756 0 : new_elem->set_neighbor(2, const_cast<RemoteElem *>(remote_elem));
2757 576 : if (elem->neighbor_ptr(2) == remote_elem)
2758 0 : new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
2759 576 : 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 350 : case QUAD9:
2765 : {
2766 500 : new_elem = Elem::build(HEX27);
2767 450 : new_elem->set_node(0, mesh.node_ptr(elem->node_ptr(0)->id() + (2*k * orig_nodes)));
2768 450 : new_elem->set_node(1, mesh.node_ptr(elem->node_ptr(1)->id() + (2*k * orig_nodes)));
2769 450 : new_elem->set_node(2, mesh.node_ptr(elem->node_ptr(2)->id() + (2*k * orig_nodes)));
2770 450 : new_elem->set_node(3, mesh.node_ptr(elem->node_ptr(3)->id() + (2*k * orig_nodes)));
2771 450 : new_elem->set_node(4, mesh.node_ptr(elem->node_ptr(0)->id() + ((2*k+2) * orig_nodes)));
2772 450 : new_elem->set_node(5, mesh.node_ptr(elem->node_ptr(1)->id() + ((2*k+2) * orig_nodes)));
2773 450 : new_elem->set_node(6, mesh.node_ptr(elem->node_ptr(2)->id() + ((2*k+2) * orig_nodes)));
2774 450 : new_elem->set_node(7, mesh.node_ptr(elem->node_ptr(3)->id() + ((2*k+2) * orig_nodes)));
2775 450 : new_elem->set_node(8, mesh.node_ptr(elem->node_ptr(4)->id() + (2*k * orig_nodes)));
2776 450 : new_elem->set_node(9, mesh.node_ptr(elem->node_ptr(5)->id() + (2*k * orig_nodes)));
2777 450 : new_elem->set_node(10, mesh.node_ptr(elem->node_ptr(6)->id() + (2*k * orig_nodes)));
2778 450 : new_elem->set_node(11, mesh.node_ptr(elem->node_ptr(7)->id() + (2*k * orig_nodes)));
2779 450 : new_elem->set_node(12, mesh.node_ptr(elem->node_ptr(0)->id() + ((2*k+1) * orig_nodes)));
2780 450 : new_elem->set_node(13, mesh.node_ptr(elem->node_ptr(1)->id() + ((2*k+1) * orig_nodes)));
2781 450 : new_elem->set_node(14, mesh.node_ptr(elem->node_ptr(2)->id() + ((2*k+1) * orig_nodes)));
2782 450 : new_elem->set_node(15, mesh.node_ptr(elem->node_ptr(3)->id() + ((2*k+1) * orig_nodes)));
2783 450 : new_elem->set_node(16, mesh.node_ptr(elem->node_ptr(4)->id() + ((2*k+2) * orig_nodes)));
2784 450 : new_elem->set_node(17, mesh.node_ptr(elem->node_ptr(5)->id() + ((2*k+2) * orig_nodes)));
2785 450 : new_elem->set_node(18, mesh.node_ptr(elem->node_ptr(6)->id() + ((2*k+2) * orig_nodes)));
2786 450 : new_elem->set_node(19, mesh.node_ptr(elem->node_ptr(7)->id() + ((2*k+2) * orig_nodes)));
2787 450 : new_elem->set_node(20, mesh.node_ptr(elem->node_ptr(8)->id() + (2*k * orig_nodes)));
2788 450 : new_elem->set_node(21, mesh.node_ptr(elem->node_ptr(4)->id() + ((2*k+1) * orig_nodes)));
2789 450 : new_elem->set_node(22, mesh.node_ptr(elem->node_ptr(5)->id() + ((2*k+1) * orig_nodes)));
2790 450 : new_elem->set_node(23, mesh.node_ptr(elem->node_ptr(6)->id() + ((2*k+1) * orig_nodes)));
2791 450 : new_elem->set_node(24, mesh.node_ptr(elem->node_ptr(7)->id() + ((2*k+1) * orig_nodes)));
2792 450 : new_elem->set_node(25, mesh.node_ptr(elem->node_ptr(8)->id() + ((2*k+2) * orig_nodes)));
2793 450 : new_elem->set_node(26, mesh.node_ptr(elem->node_ptr(8)->id() + ((2*k+1) * orig_nodes)));
2794 :
2795 450 : if (elem->neighbor_ptr(0) == remote_elem)
2796 0 : new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
2797 450 : if (elem->neighbor_ptr(1) == remote_elem)
2798 0 : new_elem->set_neighbor(2, const_cast<RemoteElem *>(remote_elem));
2799 450 : if (elem->neighbor_ptr(2) == remote_elem)
2800 0 : new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
2801 450 : 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 966 : new_elem->set_id(elem->id() + (k * orig_elem));
2814 966 : 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 966 : const unique_id_type uid = (k == 0) ?
2821 260 : elem->unique_id() :
2822 511 : 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 966 : if (!elem_subdomain)
2828 : // maintain the subdomain_id
2829 518 : new_elem->subdomain_id() = elem->subdomain_id();
2830 : else
2831 : // Allow the user to choose new subdomain_ids
2832 448 : new_elem->subdomain_id() = elem_subdomain->get_subdomain_for_layer(elem, k);
2833 :
2834 1242 : Elem * added_elem = mesh.add_elem(std::move(new_elem));
2835 :
2836 : // Copy any old boundary ids on all sides
2837 4938 : for (auto s : elem->side_index_range())
2838 : {
2839 3696 : cross_section_boundary_info.boundary_ids(elem, s, ids_to_copy);
2840 :
2841 3696 : 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 3696 : boundary_info.add_side(added_elem,
2848 2640 : 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 966 : if (k == 0)
2865 455 : boundary_info.add_side(added_elem, 0, next_side_id);
2866 966 : 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 455 : const unsigned short top_id = added_elem->dim() == 3 ?
2872 455 : cast_int<unsigned short>(elem->n_sides()+1) : 2;
2873 : boundary_info.add_side
2874 455 : (added_elem, top_id,
2875 455 : cast_int<boundary_id_type>(next_side_id+1));
2876 : }
2877 414 : }
2878 21 : }
2879 :
2880 : // Done building the mesh. Now prepare it for use.
2881 49 : 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 70 : 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 70 : const Real xavg = (xmin + xmax)/2;
3011 70 : const Real yavg = (ymin + ymax)/2;
3012 70 : const Real zavg = (zmin + zmax)/2;
3013 90 : mesh.add_point(Point(xavg,yavg,zmin), 0);
3014 90 : mesh.add_point(Point(xmax,yavg,zavg), 1);
3015 90 : mesh.add_point(Point(xavg,ymax,zavg), 2);
3016 90 : mesh.add_point(Point(xmin,yavg,zavg), 3);
3017 90 : mesh.add_point(Point(xavg,ymin,zavg), 4);
3018 90 : mesh.add_point(Point(xavg,yavg,zmax), 5);
3019 :
3020 4192 : auto add_tri = [&mesh, flip_tris](std::array<dof_id_type,3> nodes)
3021 : {
3022 560 : auto elem = mesh.add_elem(Elem::build(TRI3));
3023 560 : elem->set_node(0, mesh.node_ptr(nodes[0]));
3024 560 : elem->set_node(1, mesh.node_ptr(nodes[1]));
3025 560 : elem->set_node(2, mesh.node_ptr(nodes[2]));
3026 560 : if (flip_tris)
3027 144 : elem->flip(&mesh.get_boundary_info());
3028 580 : };
3029 :
3030 70 : add_tri({0,2,1});
3031 70 : add_tri({0,3,2});
3032 70 : add_tri({0,4,3});
3033 70 : add_tri({0,1,4});
3034 70 : add_tri({5,4,1});
3035 70 : add_tri({5,3,4});
3036 70 : add_tri({5,2,3});
3037 70 : add_tri({5,1,2});
3038 :
3039 70 : mesh.prepare_for_use();
3040 70 : }
3041 :
3042 :
3043 : } // namespace libMesh
|
