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