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 17018204 : unsigned int idx(const ElemType type,
80 : const unsigned int nx,
81 : const unsigned int i,
82 : const unsigned int j)
83 : {
84 16589022 : switch(type)
85 : {
86 3431276 : case INVALID_ELEM:
87 : case QUAD4:
88 : case QUADSHELL4:
89 : case TRI3:
90 : case TRISHELL3:
91 : {
92 3431276 : return i + j*(nx+1);
93 : }
94 :
95 13586928 : case QUAD8:
96 : case QUADSHELL8:
97 : case QUAD9:
98 : case QUADSHELL9:
99 : case TRI6:
100 : case TRI7:
101 : {
102 13586928 : 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 45132646 : 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 44039910 : switch(type)
124 : {
125 8692612 : case INVALID_ELEM:
126 : case HEX8:
127 : case PRISM6:
128 : {
129 8692612 : return i + (nx+1)*(j + k*(ny+1));
130 : }
131 :
132 36440034 : 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 36440034 : 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 279288 : 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 279288 : mesh.clear();
337 :
338 7880 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
339 :
340 279288 : if (nz != 0)
341 : {
342 129633 : mesh.set_mesh_dimension(3);
343 129633 : mesh.set_spatial_dimension(3);
344 : }
345 149655 : else if (ny != 0)
346 : {
347 112532 : mesh.set_mesh_dimension(2);
348 112532 : mesh.set_spatial_dimension(2);
349 : }
350 37123 : else if (nx != 0)
351 : {
352 34851 : mesh.set_mesh_dimension(1);
353 34851 : mesh.set_spatial_dimension(1);
354 : }
355 : else
356 : {
357 : // Will we get here?
358 2272 : mesh.set_mesh_dimension(0);
359 2272 : mesh.set_spatial_dimension(0);
360 : }
361 :
362 279288 : switch (mesh.mesh_dimension())
363 : {
364 : //---------------------------------------------------------------------
365 : // Build a 0D point
366 2272 : 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 2336 : mesh.add_point (Point(0, 0, 0), 0);
376 2272 : Elem * elem = mesh.add_elem(Elem::build(NODEELEM));
377 2272 : elem->set_node(0, mesh.node_ptr(0));
378 :
379 2208 : break;
380 : }
381 :
382 :
383 :
384 : //---------------------------------------------------------------------
385 : // Build a 1D line
386 34851 : 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 34851 : case INVALID_ELEM:
397 : case EDGE2:
398 : case EDGE3:
399 : case EDGE4:
400 : {
401 34851 : 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 11425 : case INVALID_ELEM:
413 : case EDGE2:
414 : {
415 11425 : mesh.reserve_nodes(nx+1);
416 11103 : break;
417 : }
418 :
419 21367 : case EDGE3:
420 : {
421 21367 : mesh.reserve_nodes(2*nx+1);
422 20765 : break;
423 : }
424 :
425 2059 : case EDGE4:
426 : {
427 2059 : mesh.reserve_nodes(3*nx+1);
428 2001 : 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 332317 : for (unsigned int i=0; i<=nx; i++)
445 : {
446 330048 : const Node * const node = mesh.add_point (Point(static_cast<Real>(i)/nx, 0, 0), node_id++);
447 320892 : if (i == 0)
448 11425 : boundary_info.add_node(node, 0);
449 320892 : if (i == nx)
450 11425 : boundary_info.add_node(node, 1);
451 : }
452 :
453 322 : break;
454 : }
455 :
456 602 : case EDGE3:
457 : {
458 114222 : for (unsigned int i=0; i<=2*nx; i++)
459 : {
460 92855 : const Node * const node = mesh.add_point (Point(static_cast<Real>(i)/(2*nx), 0, 0), node_id++);
461 92855 : if (i == 0)
462 21367 : boundary_info.add_node(node, 0);
463 92855 : if (i == 2*nx)
464 21367 : boundary_info.add_node(node, 1);
465 : }
466 602 : break;
467 : }
468 :
469 58 : case EDGE4:
470 : {
471 20519 : for (unsigned int i=0; i<=3*nx; i++)
472 : {
473 18980 : const Node * const node = mesh.add_point (Point(static_cast<Real>(i)/(3*nx), 0, 0), node_id++);
474 18460 : if (i == 0)
475 2059 : boundary_info.add_node(node, 0);
476 18460 : if (i == 3*nx)
477 2059 : 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 320892 : for (unsigned int i=0; i<nx; i++)
495 : {
496 309467 : Elem * elem = mesh.add_elem(Elem::build_with_id(EDGE2, i));
497 309467 : elem->set_node(0, mesh.node_ptr(i));
498 309467 : elem->set_node(1, mesh.node_ptr(i+1));
499 :
500 309467 : if (i == 0)
501 11425 : boundary_info.add_side(elem, 0, 0);
502 :
503 309467 : if (i == (nx-1))
504 11425 : boundary_info.add_side(elem, 1, 1);
505 :
506 : }
507 322 : break;
508 : }
509 :
510 602 : case EDGE3:
511 : {
512 57111 : for (unsigned int i=0; i<nx; i++)
513 : {
514 35744 : Elem * elem = mesh.add_elem(Elem::build_with_id(EDGE3, i));
515 35744 : elem->set_node(0, mesh.node_ptr(2*i));
516 35744 : elem->set_node(2, mesh.node_ptr(2*i+1));
517 35744 : elem->set_node(1, mesh.node_ptr(2*i+2));
518 :
519 35744 : if (i == 0)
520 21367 : boundary_info.add_side(elem, 0, 0);
521 :
522 35744 : if (i == (nx-1))
523 21367 : boundary_info.add_side(elem, 1, 1);
524 : }
525 602 : break;
526 : }
527 :
528 58 : case EDGE4:
529 : {
530 7526 : for (unsigned int i=0; i<nx; i++)
531 : {
532 5467 : Elem * elem = mesh.add_elem(Elem::build_with_id(EDGE4, i));
533 5467 : elem->set_node(0, mesh.node_ptr(3*i));
534 5467 : elem->set_node(2, mesh.node_ptr(3*i+1));
535 5467 : elem->set_node(3, mesh.node_ptr(3*i+2));
536 5467 : elem->set_node(1, mesh.node_ptr(3*i+3));
537 :
538 5467 : if (i == 0)
539 2059 : boundary_info.add_side(elem, 0, 0);
540 :
541 5467 : if (i == (nx-1))
542 2059 : 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 34851 : 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 500927 : for (Node * node : mesh.node_ptr_range())
560 465094 : (*node)(0) = (*node)(0)*(xmax-xmin) + xmin;
561 : }
562 :
563 : // Add sideset names to boundary info
564 34851 : boundary_info.sideset_name(0) = "left";
565 34851 : boundary_info.sideset_name(1) = "right";
566 :
567 : // Add nodeset names to boundary info
568 34851 : boundary_info.nodeset_name(0) = "left";
569 34851 : 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 112532 : 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 57039 : case INVALID_ELEM:
598 : case QUAD4:
599 : case QUADSHELL4:
600 : case QUAD8:
601 : case QUADSHELL8:
602 : case QUAD9:
603 : case QUADSHELL9:
604 : {
605 57039 : mesh.reserve_elem (nx*ny);
606 55437 : break;
607 : }
608 :
609 54854 : case TRI3:
610 : case TRISHELL3:
611 : case TRI6:
612 : case TRI7:
613 : {
614 54854 : mesh.reserve_elem (2*nx*ny);
615 53316 : break;
616 : }
617 :
618 639 : case C0POLYGON:
619 : {
620 639 : mesh.reserve_elem ((nx + 1) * (ny + 1));
621 621 : 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 29870 : case INVALID_ELEM:
635 : case QUAD4:
636 : case QUADSHELL4:
637 : case TRI3:
638 : case TRISHELL3:
639 : {
640 29870 : mesh.reserve_nodes( (nx+1)*(ny+1) );
641 29026 : break;
642 : }
643 :
644 60733 : case QUAD8:
645 : case QUADSHELL8:
646 : case QUAD9:
647 : case QUADSHELL9:
648 : case TRI6:
649 : {
650 60733 : mesh.reserve_nodes( (2*nx+1)*(2*ny+1) );
651 59033 : break;
652 : }
653 :
654 21290 : case TRI7:
655 : {
656 21290 : mesh.reserve_nodes( (2*nx+1)*(2*ny+1) + 2*nx*ny );
657 20694 : break;
658 : }
659 639 : case C0POLYGON:
660 : {
661 639 : mesh.reserve_nodes (4 + 3*nx*ny + 2*nx + 2*ny);
662 621 : 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 139013 : for (unsigned int j=0; j<=ny; j++)
684 1131274 : for (unsigned int i=0; i<=nx; i++)
685 : {
686 : const Node * const node =
687 1988086 : mesh.add_point(Point(static_cast<Real>(i) / static_cast<Real>(nx),
688 1022131 : static_cast<Real>(j) / static_cast<Real>(ny),
689 : 0.),
690 1012350 : node_id++);
691 1022131 : if (j == 0)
692 103179 : boundary_info.add_node(node, 0);
693 1022131 : if (j == ny)
694 103179 : boundary_info.add_node(node, 2);
695 1022131 : if (i == 0)
696 109143 : boundary_info.add_node(node, 3);
697 1022131 : if (i == nx)
698 109143 : 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 515142 : for (unsigned int j=0; j<=(2*ny); j++)
712 6653534 : for (unsigned int i=0; i<=(2*nx); i++)
713 : {
714 : const Node * const node =
715 12129614 : mesh.add_point(Point(static_cast<Real>(i) / static_cast<Real>(2 * nx),
716 6220415 : static_cast<Real>(j) / static_cast<Real>(2 * ny),
717 : 0),
718 6144209 : node_id++);
719 6220415 : if (j == 0)
720 446719 : boundary_info.add_node(node, 0);
721 6220415 : if (j == 2*ny)
722 446719 : boundary_info.add_node(node, 2);
723 6220415 : if (i == 0)
724 433119 : boundary_info.add_node(node, 3);
725 6220415 : if (i == 2*nx)
726 433119 : 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 82023 : if (type == TRI7)
732 58461 : for (unsigned int j=0; j<(3*ny); j += 3)
733 270192 : for (unsigned int i=0; i<(3*nx); i += 3)
734 : {
735 : // The bottom-right triangle's center node
736 455238 : mesh.add_point(Point(static_cast<Real>(i+2) / static_cast<Real>(3 * nx),
737 233021 : static_cast<Real>(j+1) / static_cast<Real>(3 * ny),
738 : 0),
739 230320 : node_id++);
740 : // The top-left triangle's center node
741 233021 : mesh.add_point(Point(static_cast<Real>(i+1) / static_cast<Real>(3 * nx),
742 233021 : static_cast<Real>(j+2) / static_cast<Real>(3 * ny),
743 : 0),
744 230320 : 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 80603 : for (unsigned int j=0; j<ny; j++)
775 865562 : for (unsigned int i=0; i<nx; i++)
776 : {
777 825512 : Elem * elem = mesh.add_elem(Elem::build_with_id(type == INVALID_ELEM ? QUAD4 : type, elem_id++));
778 803399 : elem->set_node(0, mesh.node_ptr(idx(type,nx,i,j) ));
779 803399 : elem->set_node(1, mesh.node_ptr(idx(type,nx,i+1,j) ));
780 803399 : elem->set_node(2, mesh.node_ptr(idx(type,nx,i+1,j+1)));
781 803399 : elem->set_node(3, mesh.node_ptr(idx(type,nx,i,j+1) ));
782 :
783 803399 : if (j == 0)
784 56057 : boundary_info.add_side(elem, 0, 0);
785 :
786 803399 : if (j == (ny-1))
787 56057 : boundary_info.add_side(elem, 2, 2);
788 :
789 803399 : if (i == 0)
790 62163 : boundary_info.add_side(elem, 3, 3);
791 :
792 803399 : if (i == (nx-1))
793 62163 : boundary_info.add_side(elem, 1, 1);
794 : }
795 522 : break;
796 : }
797 :
798 :
799 322 : case TRI3:
800 : case TRISHELL3:
801 : {
802 28540 : for (unsigned int j=0; j<ny; j++)
803 53390 : for (unsigned int i=0; i<nx; i++)
804 : {
805 : // Add first Tri3
806 36280 : Elem * elem = mesh.add_elem(Elem::build_with_id(type, elem_id++));
807 36280 : elem->set_node(0, mesh.node_ptr(idx(type,nx,i,j) ));
808 36280 : elem->set_node(1, mesh.node_ptr(idx(type,nx,i+1,j) ));
809 36280 : elem->set_node(2, mesh.node_ptr(idx(type,nx,i+1,j+1)));
810 :
811 36280 : if (j == 0)
812 17252 : boundary_info.add_side(elem, 0, 0);
813 :
814 36280 : if (i == (nx-1))
815 17110 : boundary_info.add_side(elem, 1, 1);
816 :
817 : // Add second Tri3
818 36280 : elem = mesh.add_elem(Elem::build_with_id(type, elem_id++));
819 36280 : elem->set_node(0, mesh.node_ptr(idx(type,nx,i,j) ));
820 36280 : elem->set_node(1, mesh.node_ptr(idx(type,nx,i+1,j+1)));
821 36280 : elem->set_node(2, mesh.node_ptr(idx(type,nx,i,j+1) ));
822 :
823 36280 : if (j == (ny-1))
824 17252 : boundary_info.add_side(elem, 1, 2);
825 :
826 36280 : if (i == 0)
827 17110 : 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 133502 : for (unsigned int j=0; j<(2*ny); j += 2)
840 923089 : for (unsigned int i=0; i<(2*nx); i += 2)
841 : {
842 828186 : Elem * elem = mesh.add_elem(Elem::build_with_id(type, elem_id++));
843 828186 : elem->set_node(0, mesh.node_ptr(idx(type,nx,i,j) ));
844 828186 : elem->set_node(1, mesh.node_ptr(idx(type,nx,i+2,j) ));
845 828186 : elem->set_node(2, mesh.node_ptr(idx(type,nx,i+2,j+2)));
846 828186 : elem->set_node(3, mesh.node_ptr(idx(type,nx,i,j+2) ));
847 828186 : elem->set_node(4, mesh.node_ptr(idx(type,nx,i+1,j) ));
848 828186 : elem->set_node(5, mesh.node_ptr(idx(type,nx,i+2,j+1)));
849 828186 : elem->set_node(6, mesh.node_ptr(idx(type,nx,i+1,j+2)));
850 828186 : elem->set_node(7, mesh.node_ptr(idx(type,nx,i,j+1) ));
851 :
852 828186 : if (type == QUAD9 || type == QUADSHELL9)
853 631872 : elem->set_node(8, mesh.node_ptr(idx(type,nx,i+1,j+1)));
854 :
855 828186 : if (j == 0)
856 100512 : boundary_info.add_side(elem, 0, 0);
857 :
858 828186 : if (j == 2*(ny-1))
859 100512 : boundary_info.add_side(elem, 2, 2);
860 :
861 828186 : if (i == 0)
862 94903 : boundary_info.add_side(elem, 3, 3);
863 :
864 828186 : if (i == 2*(nx-1))
865 94903 : boundary_info.add_side(elem, 1, 1);
866 : }
867 1080 : break;
868 : }
869 :
870 :
871 1216 : case TRI6:
872 : case TRI7:
873 : {
874 124069 : for (unsigned int j=0; j<(2*ny); j += 2)
875 608109 : for (unsigned int i=0; i<(2*nx); i += 2)
876 : {
877 : // Add first Tri in the bottom-right of its quad
878 527464 : Elem * elem = mesh.add_elem(Elem::build_with_id(type, elem_id++));
879 527464 : elem->set_node(0, mesh.node_ptr(idx(type,nx,i,j) ));
880 527464 : elem->set_node(1, mesh.node_ptr(idx(type,nx,i+2,j) ));
881 527464 : elem->set_node(2, mesh.node_ptr(idx(type,nx,i+2,j+2)));
882 527464 : elem->set_node(3, mesh.node_ptr(idx(type,nx,i+1,j) ));
883 527464 : elem->set_node(4, mesh.node_ptr(idx(type,nx,i+2,j+1)));
884 527464 : elem->set_node(5, mesh.node_ptr(idx(type,nx,i+1,j+1)));
885 :
886 527464 : if (type == TRI7)
887 233021 : elem->set_node(6, mesh.node_ptr(elem->id()+(2*nx+1)*(2*ny+1)));
888 :
889 527464 : if (j == 0)
890 81836 : boundary_info.add_side(elem, 0, 0);
891 :
892 527464 : if (i == 2*(nx-1))
893 80645 : boundary_info.add_side(elem, 1, 1);
894 :
895 : // Add second Tri in the top left of its quad
896 527464 : elem = mesh.add_elem(Elem::build_with_id(type, elem_id++));
897 527464 : elem->set_node(0, mesh.node_ptr(idx(type,nx,i,j) ));
898 527464 : elem->set_node(1, mesh.node_ptr(idx(type,nx,i+2,j+2)));
899 527464 : elem->set_node(2, mesh.node_ptr(idx(type,nx,i,j+2) ));
900 527464 : elem->set_node(3, mesh.node_ptr(idx(type,nx,i+1,j+1)));
901 527464 : elem->set_node(4, mesh.node_ptr(idx(type,nx,i+1,j+2)));
902 527464 : elem->set_node(5, mesh.node_ptr(idx(type,nx,i,j+1) ));
903 :
904 527464 : if (type == TRI7)
905 233021 : elem->set_node(6, mesh.node_ptr(elem->id()+(2*nx+1)*(2*ny+1)));
906 :
907 527464 : if (j == 2*(ny-1))
908 81836 : boundary_info.add_side(elem, 1, 2);
909 :
910 527464 : if (i == 0)
911 80645 : 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 639 : const auto dx_tri = (xmax - xmin) / (nx);
925 639 : const auto dy_tri = (ymax - ymin) / (ny + 1);
926 621 : std::unique_ptr<Elem> new_elem;
927 4544 : for (const auto i : make_range(nx + 1))
928 : {
929 : // Make new nodes for bottom layer of triangles
930 : Node *node0, *node1, *node2;
931 3905 : if (i == 0)
932 : {
933 657 : node0 = mesh.add_point(Point(0., 0, 0.));
934 657 : node1 = mesh.add_point(Point(0., dy_tri / 2., 0.));
935 657 : node2 = mesh.add_point(Point(dx_tri / 2., 0., 0.));
936 639 : node_list.push_back(node0);
937 639 : node_list.push_back(node1);
938 639 : node_list.push_back(node2);
939 : }
940 3266 : else if (i < nx)
941 : {
942 2627 : node0 = node_list.back();
943 2701 : node1 = mesh.add_point(Point((i)*dx_tri, dy_tri / 2., 0.));
944 2701 : node2 = mesh.add_point(Point((i + 1. / 2.) * dx_tri, 0., 0.));
945 2627 : node_list.push_back(node1);
946 2627 : node_list.push_back(node2);
947 : }
948 : else
949 : {
950 639 : node0 = node_list.back();
951 657 : node1 = mesh.add_point(Point((i)*dx_tri, dy_tri / 2., 0.));
952 657 : node2 = mesh.add_point(Point((i)*dx_tri, 0., 0.));
953 639 : node_list.push_back(node1);
954 639 : node_list.push_back(node2);
955 : }
956 :
957 3905 : new_elem = std::make_unique<C0Polygon>(3);
958 : // Switch to Tri3 when exodus default output supports element type mixes
959 3905 : new_elem->set_node(0, node0);
960 3905 : new_elem->set_node(1, node1);
961 3905 : new_elem->set_node(2, node2);
962 4015 : auto * elem = mesh.add_elem(std::move(new_elem));
963 :
964 : // Set boundaries
965 3905 : if (i == 0)
966 639 : boundary_info.add_side(elem, 0, 3); // left
967 3266 : else if (i == nx)
968 639 : boundary_info.add_side(elem, 1, 1); // right
969 3905 : 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 639 : (ymax - ymin - (ny == 1 ? dy_tri : (1. + (ny - 1) / 2.) * dy_tri)) / ny;
977 3905 : for (const auto j : make_range(ny))
978 : {
979 22365 : for (const auto i : make_range(nx + (j % 2)))
980 : {
981 19099 : if ((j % 2 == 0) || ((i > 0) && (i < nx)))
982 : {
983 : Node *n0, *n1, *n2, *n3, *n4, *n5;
984 16117 : n0 = node_list[running_index++];
985 16117 : n1 = node_list[running_index++];
986 16117 : n2 = node_list[running_index];
987 :
988 16117 : if (i == 0)
989 : {
990 1825 : n3 = mesh.add_point(Point(*n0) + RealVectorValue(0, hex_side, 0));
991 1775 : node_list.push_back(n3);
992 : }
993 : else
994 14342 : n3 = node_list.back();
995 :
996 16571 : n4 = mesh.add_point(Point(*n1) + RealVectorValue(0, hex_side + dy_tri, 0));
997 16571 : n5 = mesh.add_point(Point(*n2) + RealVectorValue(0, hex_side, 0));
998 16117 : node_list.push_back(n4);
999 16117 : node_list.push_back(n5);
1000 :
1001 16117 : new_elem = std::make_unique<libMesh::C0Polygon>(6);
1002 16117 : new_elem->set_node(0, n0);
1003 16117 : new_elem->set_node(1, n1);
1004 16117 : new_elem->set_node(2, n2);
1005 16117 : new_elem->set_node(3, n5);
1006 16117 : new_elem->set_node(4, n4);
1007 16117 : new_elem->set_node(5, n3);
1008 16571 : auto * elem = mesh.add_elem(std::move(new_elem));
1009 :
1010 : // Set boundaries
1011 16117 : if (i == 0)
1012 1775 : boundary_info.add_side(elem, 5, 3); // left
1013 14342 : else if (i == nx)
1014 908 : boundary_info.add_side(elem, 2, 1); // right
1015 15209 : }
1016 : // The hexagons are offset, so we build on a quad on each external side to fill
1017 2982 : else if (i == 0 || i == nx)
1018 : {
1019 : Node *n0, *n1, *n2, *n3;
1020 2982 : n0 = node_list[running_index++];
1021 2982 : n1 = node_list[running_index];
1022 :
1023 2982 : if (i == 0)
1024 : {
1025 1533 : n2 = mesh.add_point(Point(*n0) + RealVectorValue(0, hex_side + dy_tri, 0));
1026 1491 : node_list.push_back(n2);
1027 1533 : n3 = mesh.add_point(Point(*n1) + RealVectorValue(0, hex_side, 0));
1028 : }
1029 : else
1030 : {
1031 1491 : n2 = node_list.back();
1032 1533 : n3 = mesh.add_point(Point(*n1) + RealVectorValue(0, hex_side + dy_tri, 0));
1033 : }
1034 2982 : node_list.push_back(n3);
1035 :
1036 2982 : new_elem = std::make_unique<C0Polygon>(4);
1037 : // Switch to Quad4 when exodus default output supports element type mixes
1038 2982 : new_elem->set_node(0, n0);
1039 2982 : new_elem->set_node(1, n1);
1040 2982 : new_elem->set_node(3, n2);
1041 2982 : new_elem->set_node(2, n3);
1042 3066 : auto * elem = mesh.add_elem(std::move(new_elem));
1043 :
1044 : // Set boundaries
1045 2982 : if (i == 0)
1046 1491 : boundary_info.add_side(elem, 3, 3); // left
1047 1491 : else if (i == nx)
1048 1575 : 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 3266 : running_index++;
1055 :
1056 : // Skip lower right corner node
1057 3266 : if (j == 0)
1058 639 : running_index++;
1059 : }
1060 :
1061 : // Build a final layer of triangles
1062 639 : const bool ny_odd = (ny % 2 == 1);
1063 4189 : for (const auto i : make_range(nx + ny_odd))
1064 : {
1065 : // Use existing nodes, except at the corners
1066 : Node *node0, *node1, *node2;
1067 3550 : if (i == 0 && ny_odd)
1068 : {
1069 292 : node0 = mesh.add_point(Point(0., ymax, 0.));
1070 284 : node1 = node_list[running_index++];
1071 292 : node2 = node_list[running_index];
1072 : }
1073 3266 : else if (i < nx)
1074 : {
1075 2982 : node0 = node_list[running_index++];
1076 2982 : node1 = node_list[running_index++];
1077 3066 : 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 284 : node0 = node_list[running_index++];
1083 284 : node1 = node_list[running_index];
1084 292 : node2 = mesh.add_point(Point(xmax, ymax, 0.));
1085 : }
1086 :
1087 3550 : new_elem = std::make_unique<C0Polygon>(3);
1088 : // Switch to Tri3 when exodus default output supports element type mixes
1089 3550 : new_elem->set_node(0, node0);
1090 3550 : new_elem->set_node(1, node1);
1091 3550 : new_elem->set_node(2, node2);
1092 3650 : auto * elem = mesh.add_elem(std::move(new_elem));
1093 :
1094 : // Set boundaries
1095 3550 : if (i == 0)
1096 639 : boundary_info.add_side(elem, 0, 3); // left
1097 2911 : else if (i == nx)
1098 284 : boundary_info.add_side(elem, 1, 1); // right
1099 3550 : boundary_info.add_side(elem, 2, 2); // top
1100 :
1101 : }
1102 18 : break;
1103 603 : }
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 112532 : 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 7977993 : for (Node * node : mesh.node_ptr_range())
1123 : {
1124 7756087 : (*node)(0) = ((*node)(0))*(xmax-xmin) + xmin;
1125 7756087 : (*node)(1) = ((*node)(1))*(ymax-ymin) + ymin;
1126 106216 : }
1127 : }
1128 :
1129 : // Add sideset names to boundary info
1130 112532 : boundary_info.sideset_name(0) = "bottom";
1131 112532 : boundary_info.sideset_name(1) = "right";
1132 112532 : boundary_info.sideset_name(2) = "top";
1133 112532 : boundary_info.sideset_name(3) = "left";
1134 :
1135 : // Add nodeset names to boundary info
1136 112532 : boundary_info.nodeset_name(0) = "bottom";
1137 112532 : boundary_info.nodeset_name(1) = "right";
1138 112532 : boundary_info.nodeset_name(2) = "top";
1139 112532 : 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 129633 : 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 81414 : 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 81414 : mesh.reserve_elem(nx*ny*nz);
1183 79100 : break;
1184 : }
1185 :
1186 48219 : case PRISM6:
1187 : case PRISM15:
1188 : case PRISM18:
1189 : case PRISM20:
1190 : case PRISM21:
1191 : {
1192 48219 : mesh.reserve_elem(2*nx*ny*nz);
1193 46857 : 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 21648 : case INVALID_ELEM:
1208 : case HEX8:
1209 : case PRISM6:
1210 : {
1211 21648 : mesh.reserve_nodes( (nx+1)*(ny+1)*(nz+1) );
1212 21012 : break;
1213 : }
1214 :
1215 73134 : 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 73134 : mesh.reserve_nodes( (2*nx+1)*(2*ny+1)*(2*nz+1) );
1231 71072 : break;
1232 : }
1233 :
1234 12202 : case TET14:
1235 : {
1236 13902 : mesh.reserve_nodes( (2*nx+1)*(2*ny+1)*(2*nz+1) +
1237 12202 : 24*nx*ny*nz +
1238 12202 : 4*(nx*ny + ny*nz + nx*nz) );
1239 11862 : break;
1240 : }
1241 :
1242 10295 : case PRISM20:
1243 : {
1244 11165 : mesh.reserve_nodes( (2*nx+1)*(2*ny+1)*(2*nz+1) +
1245 10295 : 2*nx*ny*(nz+1) );
1246 10005 : break;
1247 : }
1248 :
1249 12354 : case PRISM21:
1250 : {
1251 13398 : mesh.reserve_nodes( (2*nx+1)*(2*ny+1)*(2*nz+1) +
1252 12354 : 2*nx*ny*(2*nz+1) );
1253 12006 : 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 79978 : for (unsigned int k=0; k<=nz; k++)
1272 267616 : for (unsigned int j=0; j<=ny; j++)
1273 1874926 : for (unsigned int i=0; i<=nx; i++)
1274 : {
1275 : const Node * const node =
1276 3236784 : mesh.add_point(Point(static_cast<Real>(i) / static_cast<Real>(nx),
1277 1665640 : static_cast<Real>(j) / static_cast<Real>(ny),
1278 1665640 : static_cast<Real>(k) / static_cast<Real>(nz)),
1279 1646038 : node_id++);
1280 1665640 : if (k == 0)
1281 298888 : boundary_info.add_node(node, 0);
1282 1665640 : if (k == nz)
1283 298888 : boundary_info.add_node(node, 5);
1284 1665640 : if (j == 0)
1285 252028 : boundary_info.add_node(node, 1);
1286 1665640 : if (j == ny)
1287 252028 : boundary_info.add_node(node, 3);
1288 1665640 : if (i == 0)
1289 209286 : boundary_info.add_node(node, 4);
1290 1665640 : if (i == nx)
1291 209286 : 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 521076 : for (unsigned int k=0; k<=(2*nz); k++)
1312 2326450 : for (unsigned int j=0; j<=(2*ny); j++)
1313 17273344 : for (unsigned int i=0; i<=(2*nx); i++)
1314 : {
1315 : const Node * const node =
1316 29973906 : mesh.add_point(Point(static_cast<Real>(i) / static_cast<Real>(2 * nx),
1317 15359985 : static_cast<Real>(j) / static_cast<Real>(2 * ny),
1318 15359985 : static_cast<Real>(k) / static_cast<Real>(2 * nz)),
1319 15173994 : node_id++);
1320 15359985 : if (k == 0)
1321 2027851 : boundary_info.add_node(node, 0);
1322 15359985 : if (k == 2*nz)
1323 2027851 : boundary_info.add_node(node, 5);
1324 15359985 : if (j == 0)
1325 1972461 : boundary_info.add_node(node, 1);
1326 15359985 : if (j == 2*ny)
1327 1972461 : boundary_info.add_node(node, 3);
1328 15359985 : if (i == 0)
1329 1913359 : boundary_info.add_node(node, 4);
1330 15359985 : if (i == 2*nx)
1331 1913359 : boundary_info.add_node(node, 2);
1332 : }
1333 :
1334 107985 : if (type == PRISM20 ||
1335 : type == PRISM21)
1336 : {
1337 22649 : const unsigned int kmax = (type == PRISM20) ? nz : 2*nz;
1338 91377 : for (unsigned int k=0; k<=kmax; k++)
1339 166992 : for (unsigned int j=0; j<ny; j++)
1340 255600 : for (unsigned int i=0; i<nx; i++)
1341 : {
1342 : const Node * const node1 =
1343 305808 : mesh.add_point(Point((static_cast<Real>(i)+1/Real(3)) / static_cast<Real>(nx),
1344 157336 : (static_cast<Real>(j)+1/Real(3)) / static_cast<Real>(ny),
1345 157336 : static_cast<Real>(k) / static_cast<Real>(kmax)),
1346 155120 : node_id++);
1347 157336 : if (k == 0)
1348 44801 : boundary_info.add_node(node1, 0);
1349 157336 : if (k == kmax)
1350 44801 : boundary_info.add_node(node1, 5);
1351 :
1352 : const Node * const node2 =
1353 305808 : mesh.add_point(Point((static_cast<Real>(i)+2/Real(3)) / static_cast<Real>(nx),
1354 157336 : (static_cast<Real>(j)+2/Real(3)) / static_cast<Real>(ny),
1355 4432 : static_cast<Real>(k) / static_cast<Real>(kmax)),
1356 155120 : node_id++);
1357 157336 : if (k == 0)
1358 44801 : boundary_info.add_node(node2, 0);
1359 157336 : if (k == kmax)
1360 44801 : 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 39728 : for (unsigned int k=0; k<nz; k++)
1383 122272 : for (unsigned int j=0; j<ny; j++)
1384 1133649 : for (unsigned int i=0; i<nx; i++)
1385 : {
1386 1037480 : Elem * elem = mesh.add_elem(Elem::build_with_id(HEX8, elem_id++));
1387 1037480 : elem->set_node(0, mesh.node_ptr(idx(type,nx,ny,i,j,k) ));
1388 1037480 : elem->set_node(1, mesh.node_ptr(idx(type,nx,ny,i+1,j,k) ));
1389 1037480 : elem->set_node(2, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) ));
1390 1037480 : elem->set_node(3, mesh.node_ptr(idx(type,nx,ny,i,j+1,k) ));
1391 1037480 : elem->set_node(4, mesh.node_ptr(idx(type,nx,ny,i,j,k+1) ));
1392 1037480 : elem->set_node(5, mesh.node_ptr(idx(type,nx,ny,i+1,j,k+1) ));
1393 1037480 : elem->set_node(6, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1)));
1394 1037480 : elem->set_node(7, mesh.node_ptr(idx(type,nx,ny,i,j+1,k+1) ));
1395 :
1396 1037480 : if (k == 0)
1397 175760 : boundary_info.add_side(elem, 0, 0);
1398 :
1399 1037480 : if (k == (nz-1))
1400 175760 : boundary_info.add_side(elem, 5, 5);
1401 :
1402 1037480 : if (j == 0)
1403 130320 : boundary_info.add_side(elem, 1, 1);
1404 :
1405 1037480 : if (j == (ny-1))
1406 130320 : boundary_info.add_side(elem, 3, 3);
1407 :
1408 1037480 : if (i == 0)
1409 96169 : boundary_info.add_side(elem, 4, 4);
1410 :
1411 1037480 : if (i == (nx-1))
1412 96169 : boundary_info.add_side(elem, 2, 2);
1413 : }
1414 410 : break;
1415 : }
1416 :
1417 :
1418 :
1419 :
1420 226 : case PRISM6:
1421 : {
1422 18602 : for (unsigned int k=0; k<nz; k++)
1423 27264 : for (unsigned int j=0; j<ny; j++)
1424 49416 : for (unsigned int i=0; i<nx; i++)
1425 : {
1426 : // First Prism
1427 32731 : Elem * elem = mesh.add_elem(Elem::build_with_id(PRISM6, elem_id++));
1428 32731 : elem->set_node(0, mesh.node_ptr(idx(type,nx,ny,i,j,k) ));
1429 32731 : elem->set_node(1, mesh.node_ptr(idx(type,nx,ny,i+1,j,k) ));
1430 32731 : elem->set_node(2, mesh.node_ptr(idx(type,nx,ny,i,j+1,k) ));
1431 32731 : elem->set_node(3, mesh.node_ptr(idx(type,nx,ny,i,j,k+1) ));
1432 32731 : elem->set_node(4, mesh.node_ptr(idx(type,nx,ny,i+1,j,k+1) ));
1433 32731 : 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 32731 : if (i==0)
1437 16685 : boundary_info.add_side(elem, 3, 4);
1438 :
1439 32731 : if (j==0)
1440 16685 : boundary_info.add_side(elem, 1, 1);
1441 :
1442 32731 : if (k==0)
1443 16685 : boundary_info.add_side(elem, 0, 0);
1444 :
1445 32731 : if (k == (nz-1))
1446 16685 : boundary_info.add_side(elem, 4, 5);
1447 :
1448 : // Second Prism
1449 32731 : elem = mesh.add_elem(Elem::build_with_id(PRISM6, elem_id++));
1450 32731 : elem->set_node(0, mesh.node_ptr(idx(type,nx,ny,i+1,j,k) ));
1451 32731 : elem->set_node(1, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) ));
1452 32731 : elem->set_node(2, mesh.node_ptr(idx(type,nx,ny,i,j+1,k) ));
1453 32731 : elem->set_node(3, mesh.node_ptr(idx(type,nx,ny,i+1,j,k+1) ));
1454 32731 : elem->set_node(4, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1)));
1455 32731 : 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 32731 : if (i == (nx-1))
1459 16685 : boundary_info.add_side(elem, 1, 2);
1460 :
1461 32731 : if (j == (ny-1))
1462 16685 : boundary_info.add_side(elem, 2, 3);
1463 :
1464 32731 : if (k==0)
1465 16685 : boundary_info.add_side(elem, 0, 0);
1466 :
1467 32731 : if (k == (nz-1))
1468 16685 : 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 168700 : for (unsigned int k=0; k<(2*nz); k += 2)
1489 324917 : for (unsigned int j=0; j<(2*ny); j += 2)
1490 1426934 : for (unsigned int i=0; i<(2*nx); i += 2)
1491 : {
1492 1202928 : ElemType build_type = (type == HEX20) ? HEX20 : HEX27;
1493 1202928 : Elem * elem = mesh.add_elem(Elem::build_with_id(build_type, elem_id++));
1494 :
1495 1202928 : elem->set_node(0, mesh.node_ptr(idx(type,nx,ny,i, j, k) ));
1496 1202928 : elem->set_node(1, mesh.node_ptr(idx(type,nx,ny,i+2,j, k) ));
1497 1202928 : elem->set_node(2, mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k) ));
1498 1202928 : elem->set_node(3, mesh.node_ptr(idx(type,nx,ny,i, j+2,k) ));
1499 1202928 : elem->set_node(4, mesh.node_ptr(idx(type,nx,ny,i, j, k+2)));
1500 1202928 : elem->set_node(5, mesh.node_ptr(idx(type,nx,ny,i+2,j, k+2)));
1501 1202928 : elem->set_node(6, mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k+2)));
1502 1202928 : elem->set_node(7, mesh.node_ptr(idx(type,nx,ny,i, j+2,k+2)));
1503 1202928 : elem->set_node(8, mesh.node_ptr(idx(type,nx,ny,i+1,j, k) ));
1504 1202928 : elem->set_node(9, mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k) ));
1505 1202928 : elem->set_node(10, mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k) ));
1506 1202928 : elem->set_node(11, mesh.node_ptr(idx(type,nx,ny,i, j+1,k) ));
1507 1202928 : elem->set_node(12, mesh.node_ptr(idx(type,nx,ny,i, j, k+1)));
1508 1202928 : elem->set_node(13, mesh.node_ptr(idx(type,nx,ny,i+2,j, k+1)));
1509 1202928 : elem->set_node(14, mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k+1)));
1510 1202928 : elem->set_node(15, mesh.node_ptr(idx(type,nx,ny,i, j+2,k+1)));
1511 1202928 : elem->set_node(16, mesh.node_ptr(idx(type,nx,ny,i+1,j, k+2)));
1512 1202928 : elem->set_node(17, mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k+2)));
1513 1202928 : elem->set_node(18, mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k+2)));
1514 1202928 : elem->set_node(19, mesh.node_ptr(idx(type,nx,ny,i, j+1,k+2)));
1515 :
1516 1202928 : if ((type == HEX27) || (type == TET4) || (type == TET10) || (type == TET14) ||
1517 5426 : (type == PYRAMID5) || (type == PYRAMID13) || (type == PYRAMID14) ||
1518 1870 : (type == PYRAMID18))
1519 : {
1520 1167570 : elem->set_node(20, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) ));
1521 1167570 : elem->set_node(21, mesh.node_ptr(idx(type,nx,ny,i+1,j, k+1)));
1522 1167570 : elem->set_node(22, mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k+1)));
1523 1167570 : elem->set_node(23, mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k+1)));
1524 1167570 : elem->set_node(24, mesh.node_ptr(idx(type,nx,ny,i, j+1,k+1)));
1525 1167570 : elem->set_node(25, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+2)));
1526 1167570 : elem->set_node(26, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1)));
1527 : }
1528 :
1529 1202928 : if (k == 0)
1530 250098 : boundary_info.add_side(elem, 0, 0);
1531 :
1532 1202928 : if (k == 2*(nz-1))
1533 250098 : boundary_info.add_side(elem, 5, 5);
1534 :
1535 1202928 : if (j == 0)
1536 236511 : boundary_info.add_side(elem, 1, 1);
1537 :
1538 1202928 : if (j == 2*(ny-1))
1539 236511 : boundary_info.add_side(elem, 3, 3);
1540 :
1541 1202928 : if (i == 0)
1542 224006 : boundary_info.add_side(elem, 4, 4);
1543 :
1544 1202928 : if (i == 2*(nx-1))
1545 224006 : 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 91838 : for (unsigned int k=0; k<(2*nz); k += 2)
1559 126216 : for (unsigned int j=0; j<(2*ny); j += 2)
1560 195192 : for (unsigned int i=0; i<(2*nx); i += 2)
1561 : {
1562 : // First Prism
1563 120618 : Elem * elem = mesh.add_elem(Elem::build_with_id(type, elem_id++));
1564 120618 : elem->set_node(0, mesh.node_ptr(idx(type,nx,ny,i, j, k) ));
1565 120618 : elem->set_node(1, mesh.node_ptr(idx(type,nx,ny,i+2,j, k) ));
1566 120618 : elem->set_node(2, mesh.node_ptr(idx(type,nx,ny,i, j+2,k) ));
1567 120618 : elem->set_node(3, mesh.node_ptr(idx(type,nx,ny,i, j, k+2)));
1568 120618 : elem->set_node(4, mesh.node_ptr(idx(type,nx,ny,i+2,j, k+2)));
1569 120618 : elem->set_node(5, mesh.node_ptr(idx(type,nx,ny,i, j+2,k+2)));
1570 120618 : elem->set_node(6, mesh.node_ptr(idx(type,nx,ny,i+1,j, k) ));
1571 120618 : elem->set_node(7, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) ));
1572 120618 : elem->set_node(8, mesh.node_ptr(idx(type,nx,ny,i, j+1,k) ));
1573 120618 : elem->set_node(9, mesh.node_ptr(idx(type,nx,ny,i, j, k+1)));
1574 120618 : elem->set_node(10, mesh.node_ptr(idx(type,nx,ny,i+2,j, k+1)));
1575 120618 : elem->set_node(11, mesh.node_ptr(idx(type,nx,ny,i, j+2,k+1)));
1576 120618 : elem->set_node(12, mesh.node_ptr(idx(type,nx,ny,i+1,j, k+2)));
1577 120618 : elem->set_node(13, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+2)));
1578 120618 : elem->set_node(14, mesh.node_ptr(idx(type,nx,ny,i, j+1,k+2)));
1579 :
1580 124170 : if (type == PRISM18 ||
1581 118770 : type == PRISM20 ||
1582 : type == PRISM21)
1583 : {
1584 98324 : elem->set_node(15, mesh.node_ptr(idx(type,nx,ny,i+1,j, k+1)));
1585 98324 : elem->set_node(16, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1)));
1586 98324 : elem->set_node(17, mesh.node_ptr(idx(type,nx,ny,i, j+1,k+1)));
1587 : }
1588 :
1589 120618 : if (type == PRISM20)
1590 : {
1591 36139 : const dof_id_type base_idx = (2*nx+1)*(2*ny+1)*(2*nz+1);
1592 36139 : elem->set_node(18, mesh.node_ptr(base_idx+((k/2)*(nx*ny)+j/2*nx+i/2)*2));
1593 36139 : elem->set_node(19, mesh.node_ptr(base_idx+(((k/2)+1)*(nx*ny)+j/2*nx+i/2)*2));
1594 : }
1595 :
1596 120618 : if (type == PRISM21)
1597 : {
1598 38198 : const dof_id_type base_idx = (2*nx+1)*(2*ny+1)*(2*nz+1);
1599 38198 : elem->set_node(18, mesh.node_ptr(base_idx+(k*(nx*ny)+j/2*nx+i/2)*2));
1600 38198 : elem->set_node(19, mesh.node_ptr(base_idx+((k+2)*(nx*ny)+j/2*nx+i/2)*2));
1601 38198 : 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 120618 : if (i==0)
1606 74574 : boundary_info.add_side(elem, 3, 4);
1607 :
1608 120618 : if (j==0)
1609 74574 : boundary_info.add_side(elem, 1, 1);
1610 :
1611 120618 : if (k==0)
1612 74624 : boundary_info.add_side(elem, 0, 0);
1613 :
1614 120618 : if (k == 2*(nz-1))
1615 74624 : boundary_info.add_side(elem, 4, 5);
1616 :
1617 :
1618 : // Second Prism
1619 120618 : elem = mesh.add_elem(Elem::build_with_id(type, elem_id++));
1620 120618 : elem->set_node(0, mesh.node_ptr(idx(type,nx,ny,i+2,j,k) ));
1621 120618 : elem->set_node(1, mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k) ));
1622 120618 : elem->set_node(2, mesh.node_ptr(idx(type,nx,ny,i,j+2,k) ));
1623 120618 : elem->set_node(3, mesh.node_ptr(idx(type,nx,ny,i+2,j,k+2) ));
1624 120618 : elem->set_node(4, mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k+2) ));
1625 120618 : elem->set_node(5, mesh.node_ptr(idx(type,nx,ny,i,j+2,k+2) ));
1626 120618 : elem->set_node(6, mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k) ));
1627 120618 : elem->set_node(7, mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k) ));
1628 120618 : elem->set_node(8, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) ));
1629 120618 : elem->set_node(9, mesh.node_ptr(idx(type,nx,ny,i+2,j,k+1) ));
1630 120618 : elem->set_node(10, mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k+1)));
1631 120618 : elem->set_node(11, mesh.node_ptr(idx(type,nx,ny,i,j+2,k+1) ));
1632 120618 : elem->set_node(12, mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k+2)));
1633 120618 : elem->set_node(13, mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k+2)));
1634 120618 : elem->set_node(14, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+2)));
1635 :
1636 120618 : if (type == PRISM18 ||
1637 60492 : type == PRISM20 ||
1638 : type == PRISM21)
1639 : {
1640 98324 : elem->set_node(15, mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k+1)));
1641 98324 : elem->set_node(16, mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k+1)));
1642 98324 : elem->set_node(17, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1)));
1643 : }
1644 :
1645 120618 : if (type == PRISM20)
1646 : {
1647 36139 : const dof_id_type base_idx = (2*nx+1)*(2*ny+1)*(2*nz+1);
1648 36139 : elem->set_node(18, mesh.node_ptr(base_idx+((k/2)*(nx*ny)+j/2*nx+i/2)*2+1));
1649 36139 : elem->set_node(19, mesh.node_ptr(base_idx+(((k/2)+1)*(nx*ny)+j/2*nx+i/2)*2+1));
1650 : }
1651 :
1652 120618 : if (type == PRISM21)
1653 : {
1654 38198 : const dof_id_type base_idx = (2*nx+1)*(2*ny+1)*(2*nz+1);
1655 38198 : elem->set_node(18, mesh.node_ptr(base_idx+(k*(nx*ny)+j/2*nx+i/2)*2+1));
1656 38198 : elem->set_node(19, mesh.node_ptr(base_idx+((k+2)*(nx*ny)+j/2*nx+i/2)*2+1));
1657 38198 : 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 120618 : if (i == 2*(nx-1))
1662 74574 : boundary_info.add_side(elem, 1, 2);
1663 :
1664 120618 : if (j == 2*(ny-1))
1665 74574 : boundary_info.add_side(elem, 2, 3);
1666 :
1667 120618 : if (k==0)
1668 74624 : boundary_info.add_side(elem, 0, 0);
1669 :
1670 120618 : if (k == 2*(nz-1))
1671 74624 : 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 129633 : 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 17469930 : for (unsigned int p=0; p<mesh.n_nodes(); p++)
1700 : {
1701 17340297 : mesh.node_ref(p)(0) = (mesh.node_ref(p)(0))*(xmax-xmin) + xmin;
1702 17340297 : mesh.node_ref(p)(1) = (mesh.node_ref(p)(1))*(ymax-ymin) + ymin;
1703 17340297 : 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 133309 : if ((type == TET4) ||
1717 129133 : (type == TET10) ||
1718 128793 : (type == TET14) ||
1719 99669 : (type == PYRAMID5) ||
1720 2682 : (type == PYRAMID13) ||
1721 96886 : (type == PYRAMID14) ||
1722 91598 : (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 40536 : std::unique_ptr<Elem> side;
1729 :
1730 40536 : if ((type == TET4) || (type == TET10) || (type == TET14))
1731 29886 : new_elements.reserve(24*mesh.n_elem());
1732 : else
1733 10650 : new_elements.reserve(6*mesh.n_elem());
1734 :
1735 : // Create tetrahedra or pyramids
1736 548874 : for (auto & base_hex : mesh.element_ptr_range())
1737 : {
1738 : // Get a pointer to the node located at the HEX27 center
1739 240789 : 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 1685523 : for (auto s : base_hex->side_index_range())
1745 : {
1746 : // Get the boundary ID(s) for this side
1747 1444734 : 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 1444734 : boundary_id_type b_id = ids.empty() ? BoundaryInfo::invalid_id : ids[0];
1754 :
1755 : // Need to build the full-ordered side!
1756 1444734 : base_hex->build_side_ptr(side, s);
1757 :
1758 1444734 : if ((type == TET4) || (type == TET10) || (type == TET14))
1759 : {
1760 : // Build 4 sub-tets per side
1761 5010600 : for (unsigned int sub_tet=0; sub_tet<4; ++sub_tet)
1762 : {
1763 7915200 : new_elements.push_back( Elem::build(TET4) );
1764 101760 : auto & sub_elem = new_elements.back();
1765 4212000 : sub_elem->set_node(0, side->node_ptr(sub_tet));
1766 4212000 : sub_elem->set_node(1, side->node_ptr(8)); // center of the face
1767 4110240 : sub_elem->set_node(2, side->node_ptr(sub_tet==3 ? 0 : sub_tet+1 )); // wrap-around
1768 4008480 : 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 4008480 : if (b_id != BoundaryInfo::invalid_id)
1774 1726312 : 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 872760 : 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 467550 : sub_elem->set_node(0, side->node_ptr(0));
1788 467550 : sub_elem->set_node(1, side->node_ptr(3));
1789 467550 : sub_elem->set_node(2, side->node_ptr(2));
1790 467550 : sub_elem->set_node(3, side->node_ptr(1));
1791 :
1792 : // Set the apex
1793 442614 : 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 442614 : if (b_id != BoundaryInfo::invalid_id)
1798 222066 : boundary_info.add_side(sub_elem.get(), 4, b_id);
1799 : } // end else type==PYRAMID*
1800 : }
1801 38256 : }
1802 :
1803 :
1804 : // Delete the original HEX27 elements from the mesh, and the boundary info structure.
1805 548874 : for (auto & elem : mesh.element_ptr_range())
1806 : {
1807 240789 : boundary_info.remove(elem); // Safe even if elem has no boundary info.
1808 240789 : mesh.delete_elem(elem);
1809 38256 : }
1810 :
1811 : // Add the new elements
1812 4491630 : for (auto i : index_range(new_elements))
1813 : {
1814 399798 : new_elements[i]->set_id(i);
1815 4679550 : mesh.add_elem( std::move(new_elements[i]) );
1816 : }
1817 :
1818 38256 : } // 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 129633 : if ((type == TET10) || (type == PYRAMID14))
1823 11791 : mesh.all_second_order();
1824 :
1825 117842 : else if (type == PYRAMID13)
1826 2698 : mesh.all_second_order(/*full_ordered=*/false);
1827 :
1828 115144 : else if ((type == TET14) || (type == PYRAMID18))
1829 14687 : mesh.all_complete_order();
1830 :
1831 :
1832 : // Add sideset names to boundary info (Z axis out of the screen)
1833 129633 : boundary_info.sideset_name(0) = "back";
1834 129633 : boundary_info.sideset_name(1) = "bottom";
1835 129633 : boundary_info.sideset_name(2) = "right";
1836 129633 : boundary_info.sideset_name(3) = "top";
1837 129633 : boundary_info.sideset_name(4) = "left";
1838 129633 : boundary_info.sideset_name(5) = "front";
1839 :
1840 : // Add nodeset names to boundary info
1841 129633 : boundary_info.nodeset_name(0) = "back";
1842 129633 : boundary_info.nodeset_name(1) = "bottom";
1843 129633 : boundary_info.nodeset_name(2) = "right";
1844 129633 : boundary_info.nodeset_name(3) = "top";
1845 129633 : boundary_info.nodeset_name(4) = "left";
1846 129633 : 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 279288 : mesh.prepare_for_use ();
1857 279288 : }
1858 :
1859 :
1860 :
1861 994 : 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 994 : build_cube(mesh,
1870 : 0, 0, 0,
1871 : 0., 0.,
1872 : 0., 0.,
1873 : 0., 0.,
1874 : type,
1875 : gauss_lobatto_grid);
1876 994 : }
1877 :
1878 :
1879 4819 : 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 4819 : build_cube(mesh,
1890 : nx, 0, 0,
1891 : xmin, xmax,
1892 : 0., 0.,
1893 : 0., 0.,
1894 : type,
1895 : gauss_lobatto_grid);
1896 4819 : }
1897 :
1898 :
1899 :
1900 21227 : 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 21227 : build_cube (mesh,
1914 : nx, ny, 0,
1915 : xmin, xmax,
1916 : ymin, ymax,
1917 : 0., 0.,
1918 : type,
1919 : gauss_lobatto_grid);
1920 21227 : }
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 2281 : 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 2281 : cast_int<unsigned char>(mesh.mesh_dimension());
1961 2281 : mesh.clear();
1962 2281 : 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 2281 : if (mesh.mesh_dimension() == 1)
1971 : {
1972 2281 : switch (type)
1973 : {
1974 857 : case HEX8:
1975 : case HEX27:
1976 : case TET4:
1977 : case TET10:
1978 : case TET14:
1979 857 : mesh.set_mesh_dimension(3);
1980 26 : break;
1981 1140 : 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 1140 : mesh.set_mesh_dimension(2);
1991 34 : break;
1992 284 : case EDGE2:
1993 : case EDGE3:
1994 : case EDGE4:
1995 284 : 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 2281 : const bool is_replicated = mesh.is_replicated();
2009 :
2010 : // Sphere is centered at origin by default
2011 68 : const Point cent;
2012 :
2013 4562 : const Sphere sphere (cent, rad);
2014 :
2015 2281 : switch (mesh.mesh_dimension())
2016 : {
2017 : //-----------------------------------------------------------------
2018 : // Build a line in one dimension
2019 284 : case 1:
2020 : {
2021 284 : 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 1140 : 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 1140 : if (flat)
2040 : {
2041 34 : const Real sqrt_2 = std::sqrt(2.);
2042 1140 : const Real rad_2 = .25*rad;
2043 1140 : const Real rad_sqrt_2 = rad/sqrt_2;
2044 :
2045 : // (Temporary) convenient storage for node pointers
2046 1174 : std::vector<Node *> nodes(8);
2047 :
2048 : // Point 0
2049 1174 : nodes[0] = mesh.add_point (Point(-rad_2,-rad_2, 0.), node_id++);
2050 :
2051 : // Point 1
2052 1174 : nodes[1] = mesh.add_point (Point( rad_2,-rad_2, 0.), node_id++);
2053 :
2054 : // Point 2
2055 1174 : nodes[2] = mesh.add_point (Point( rad_2, rad_2, 0.), node_id++);
2056 :
2057 : // Point 3
2058 1174 : nodes[3] = mesh.add_point (Point(-rad_2, rad_2, 0.), node_id++);
2059 :
2060 : // Point 4
2061 1174 : nodes[4] = mesh.add_point (Point(-rad_sqrt_2,-rad_sqrt_2, 0.), node_id++);
2062 :
2063 : // Point 5
2064 1174 : nodes[5] = mesh.add_point (Point( rad_sqrt_2,-rad_sqrt_2, 0.), node_id++);
2065 :
2066 : // Point 6
2067 1174 : nodes[6] = mesh.add_point (Point( rad_sqrt_2, rad_sqrt_2, 0.), node_id++);
2068 :
2069 : // Point 7
2070 1174 : 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 1140 : Elem * elem0 = mesh.add_elem (Elem::build(QUAD4));
2077 1140 : elem0->set_node(0, nodes[0]);
2078 1140 : elem0->set_node(1, nodes[1]);
2079 1140 : elem0->set_node(2, nodes[2]);
2080 1140 : elem0->set_node(3, nodes[3]);
2081 : }
2082 :
2083 : // Element 1
2084 : {
2085 1140 : Elem * elem1 = mesh.add_elem (Elem::build(QUAD4));
2086 1140 : elem1->set_node(0, nodes[4]);
2087 1140 : elem1->set_node(1, nodes[0]);
2088 1140 : elem1->set_node(2, nodes[3]);
2089 1140 : elem1->set_node(3, nodes[7]);
2090 : }
2091 :
2092 : // Element 2
2093 : {
2094 1140 : Elem * elem2 = mesh.add_elem (Elem::build(QUAD4));
2095 1140 : elem2->set_node(0, nodes[4]);
2096 1140 : elem2->set_node(1, nodes[5]);
2097 1140 : elem2->set_node(2, nodes[1]);
2098 1140 : elem2->set_node(3, nodes[0]);
2099 : }
2100 :
2101 : // Element 3
2102 : {
2103 1140 : Elem * elem3 = mesh.add_elem (Elem::build(QUAD4));
2104 1140 : elem3->set_node(0, nodes[1]);
2105 1140 : elem3->set_node(1, nodes[5]);
2106 1140 : elem3->set_node(2, nodes[6]);
2107 1140 : elem3->set_node(3, nodes[2]);
2108 : }
2109 :
2110 : // Element 4
2111 : {
2112 1140 : Elem * elem4 = mesh.add_elem (Elem::build(QUAD4));
2113 1140 : elem4->set_node(0, nodes[3]);
2114 1140 : elem4->set_node(1, nodes[2]);
2115 1140 : elem4->set_node(2, nodes[6]);
2116 1140 : 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 857 : case 3:
2185 : {
2186 : // (Currently) supported types
2187 857 : 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 857 : r_small = 0.25*rad, // 0.25 *radius
2197 857 : r_med = (0.125*std::sqrt(2.)+0.5)*rad; // .67677*radius
2198 :
2199 : // (Temporary) convenient storage for node pointers
2200 883 : 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 883 : nodes[0] = mesh.add_point (Point(-r_small,-r_small, -r_small), node_id++);
2210 883 : nodes[1] = mesh.add_point (Point( r_small,-r_small, -r_small), node_id++);
2211 883 : nodes[2] = mesh.add_point (Point( r_small, r_small, -r_small), node_id++);
2212 883 : nodes[3] = mesh.add_point (Point(-r_small, r_small, -r_small), node_id++);
2213 883 : nodes[4] = mesh.add_point (Point(-r_small,-r_small, r_small), node_id++);
2214 883 : nodes[5] = mesh.add_point (Point( r_small,-r_small, r_small), node_id++);
2215 883 : nodes[6] = mesh.add_point (Point( r_small, r_small, r_small), node_id++);
2216 883 : 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 883 : nodes[8] = mesh.add_point (Point(-r_med,-r_med, -r_med), node_id++);
2220 883 : nodes[9] = mesh.add_point (Point( r_med,-r_med, -r_med), node_id++);
2221 883 : nodes[10] = mesh.add_point (Point( r_med, r_med, -r_med), node_id++);
2222 883 : nodes[11] = mesh.add_point (Point(-r_med, r_med, -r_med), node_id++);
2223 883 : nodes[12] = mesh.add_point (Point(-r_med,-r_med, r_med), node_id++);
2224 883 : nodes[13] = mesh.add_point (Point( r_med,-r_med, r_med), node_id++);
2225 883 : nodes[14] = mesh.add_point (Point( r_med, r_med, r_med), node_id++);
2226 883 : 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 857 : Elem * elem0 = mesh.add_elem(Elem::build(HEX8));
2232 857 : elem0->set_node(0, nodes[0]);
2233 857 : elem0->set_node(1, nodes[1]);
2234 857 : elem0->set_node(2, nodes[2]);
2235 857 : elem0->set_node(3, nodes[3]);
2236 857 : elem0->set_node(4, nodes[4]);
2237 857 : elem0->set_node(5, nodes[5]);
2238 857 : elem0->set_node(6, nodes[6]);
2239 857 : elem0->set_node(7, nodes[7]);
2240 : }
2241 :
2242 : // Element 1 - "bottom"
2243 : {
2244 857 : Elem * elem1 = mesh.add_elem(Elem::build(HEX8));
2245 857 : elem1->set_node(0, nodes[8]);
2246 857 : elem1->set_node(1, nodes[9]);
2247 857 : elem1->set_node(2, nodes[10]);
2248 857 : elem1->set_node(3, nodes[11]);
2249 857 : elem1->set_node(4, nodes[0]);
2250 857 : elem1->set_node(5, nodes[1]);
2251 857 : elem1->set_node(6, nodes[2]);
2252 857 : elem1->set_node(7, nodes[3]);
2253 : }
2254 :
2255 : // Element 2 - "front"
2256 : {
2257 857 : Elem * elem2 = mesh.add_elem(Elem::build(HEX8));
2258 857 : elem2->set_node(0, nodes[8]);
2259 857 : elem2->set_node(1, nodes[9]);
2260 857 : elem2->set_node(2, nodes[1]);
2261 857 : elem2->set_node(3, nodes[0]);
2262 857 : elem2->set_node(4, nodes[12]);
2263 857 : elem2->set_node(5, nodes[13]);
2264 857 : elem2->set_node(6, nodes[5]);
2265 857 : elem2->set_node(7, nodes[4]);
2266 : }
2267 :
2268 : // Element 3 - "right"
2269 : {
2270 857 : Elem * elem3 = mesh.add_elem(Elem::build(HEX8));
2271 857 : elem3->set_node(0, nodes[1]);
2272 857 : elem3->set_node(1, nodes[9]);
2273 857 : elem3->set_node(2, nodes[10]);
2274 857 : elem3->set_node(3, nodes[2]);
2275 857 : elem3->set_node(4, nodes[5]);
2276 857 : elem3->set_node(5, nodes[13]);
2277 857 : elem3->set_node(6, nodes[14]);
2278 857 : elem3->set_node(7, nodes[6]);
2279 : }
2280 :
2281 : // Element 4 - "back"
2282 : {
2283 857 : Elem * elem4 = mesh.add_elem(Elem::build(HEX8));
2284 857 : elem4->set_node(0, nodes[3]);
2285 857 : elem4->set_node(1, nodes[2]);
2286 857 : elem4->set_node(2, nodes[10]);
2287 857 : elem4->set_node(3, nodes[11]);
2288 857 : elem4->set_node(4, nodes[7]);
2289 857 : elem4->set_node(5, nodes[6]);
2290 857 : elem4->set_node(6, nodes[14]);
2291 857 : elem4->set_node(7, nodes[15]);
2292 : }
2293 :
2294 : // Element 5 - "left"
2295 : {
2296 857 : Elem * elem5 = mesh.add_elem(Elem::build(HEX8));
2297 857 : elem5->set_node(0, nodes[8]);
2298 857 : elem5->set_node(1, nodes[0]);
2299 857 : elem5->set_node(2, nodes[3]);
2300 857 : elem5->set_node(3, nodes[11]);
2301 857 : elem5->set_node(4, nodes[12]);
2302 857 : elem5->set_node(5, nodes[4]);
2303 857 : elem5->set_node(6, nodes[7]);
2304 857 : elem5->set_node(7, nodes[15]);
2305 : }
2306 :
2307 : // Element 6 - "top"
2308 : {
2309 857 : Elem * elem6 = mesh.add_elem(Elem::build(HEX8));
2310 857 : elem6->set_node(0, nodes[4]);
2311 857 : elem6->set_node(1, nodes[5]);
2312 857 : elem6->set_node(2, nodes[6]);
2313 857 : elem6->set_node(3, nodes[7]);
2314 857 : elem6->set_node(4, nodes[12]);
2315 857 : elem6->set_node(5, nodes[13]);
2316 857 : elem6->set_node(6, nodes[14]);
2317 857 : 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 4562 : MeshRefinement mesh_refinement (mesh);
2334 :
2335 : // For avoiding extraneous element side construction
2336 2281 : std::unique_ptr<Elem> side;
2337 :
2338 : // Loop over the elements, refine, pop nodes to boundary.
2339 5490 : 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 3209 : if (!is_replicated)
2345 2196 : mesh.prepare_for_use();
2346 :
2347 3209 : mesh_refinement.uniformly_refine(1);
2348 :
2349 436360 : for (const auto & elem : mesh.active_element_ptr_range())
2350 1510919 : for (auto s : elem->side_index_range())
2351 1333676 : if (elem->neighbor_ptr(s) == nullptr || (mesh.mesh_dimension() == 2 && !flat))
2352 : {
2353 53988 : 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 252675 : for (auto n : side->node_index_range())
2359 : {
2360 9272 : Node & side_node = side->node_ref(n);
2361 : side_node =
2362 198687 : sphere.closest_point(side->point(n));
2363 :
2364 200167 : if (!is_replicated &&
2365 104327 : side_node.processor_id() != mesh.processor_id())
2366 72871 : moved_ghost_nodes.insert(side_node.id());
2367 : }
2368 3017 : }
2369 :
2370 3209 : if (!is_replicated)
2371 : {
2372 24 : std::map<processor_id_type, std::vector<dof_id_type>> moved_nodes_map;
2373 28356 : for (auto id : moved_ghost_nodes)
2374 : {
2375 26160 : const Node & node = mesh.node_ref(id);
2376 26160 : moved_nodes_map[node.processor_id()].push_back(node.id());
2377 : }
2378 :
2379 : auto action_functor =
2380 6523 : [& mesh, & sphere]
2381 : (processor_id_type /* pid */,
2382 52856 : const std::vector<dof_id_type> & my_moved_nodes)
2383 : {
2384 32715 : for (auto id : my_moved_nodes)
2385 : {
2386 26160 : Node & node = mesh.node_ref(id);
2387 26160 : node = sphere.closest_point(node);
2388 : }
2389 2204 : };
2390 :
2391 : // First get new node positions to their owners
2392 : Parallel::push_parallel_vector_data
2393 2196 : (mesh.comm(), moved_nodes_map, action_functor);
2394 :
2395 : // Then get node positions to anyone else with them ghosted
2396 2196 : SyncNodalPositions sync_object(mesh);
2397 : Parallel::sync_dofobject_data_by_id
2398 4368 : (mesh.comm(), mesh.nodes_begin(), mesh.nodes_end(),
2399 : sync_object);
2400 : }
2401 : }
2402 :
2403 : // A DistributedMesh needs a little prep before flattening
2404 2281 : if (!is_replicated)
2405 1706 : 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 2281 : MeshTools::Modification::flatten(mesh);
2412 :
2413 : // Convert all the tensor product elements to simplices if requested
2414 2281 : 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 497 : if (is_replicated)
2419 106 : mesh.prepare_for_use();
2420 :
2421 497 : MeshTools::Modification::all_tri(mesh);
2422 : }
2423 :
2424 : // Convert to second-order elements if the user requested it.
2425 2349 : if (Elem::build(type)->default_order() != FIRST)
2426 : {
2427 1571 : 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 1571 : bool full_ordered = !((type==QUAD8) || (type==HEX20));
2436 1571 : mesh.all_second_order(full_ordered);
2437 : }
2438 :
2439 : // And pop to the boundary again...
2440 176670 : for (const auto & elem : mesh.active_element_ptr_range())
2441 612961 : for (auto s : elem->side_index_range())
2442 545469 : if (elem->neighbor_ptr(s) == nullptr)
2443 : {
2444 20319 : elem->build_side_ptr(side, s);
2445 :
2446 : // Pop each point to the sphere boundary
2447 183824 : for (auto n : side->node_index_range())
2448 173327 : side->point(n) =
2449 173327 : sphere.closest_point(side->point(n));
2450 1475 : }
2451 : }
2452 :
2453 :
2454 : // The meshes could probably use some smoothing.
2455 2281 : if (mesh.mesh_dimension() > 1)
2456 : {
2457 2057 : LaplaceMeshSmoother smoother(mesh, n_smooth);
2458 1997 : smoother.smooth();
2459 : }
2460 :
2461 : // We'll give the whole sphere surface a boundary id of 0
2462 810930 : for (const auto & elem : mesh.active_element_ptr_range())
2463 2343660 : for (auto s : elem->side_index_range())
2464 1988314 : if (!elem->neighbor_ptr(s))
2465 53135 : boundary_info.add_side(elem, s, 0);
2466 :
2467 : // Done building the mesh. Now prepare it for use.
2468 2281 : mesh.prepare_for_use();
2469 2281 : }
2470 :
2471 : #endif // #ifndef LIBMESH_ENABLE_AMR
2472 :
2473 :
2474 : // Meshes the tensor product of a 1D and a 1D-or-2D domain.
2475 497 : 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 497 : if (!cross_section.n_elem())
2484 0 : return;
2485 :
2486 497 : dof_id_type orig_elem = cross_section.n_elem();
2487 497 : dof_id_type orig_nodes = cross_section.n_nodes();
2488 :
2489 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
2490 497 : unique_id_type orig_unique_ids = cross_section.parallel_max_unique_id();
2491 : #endif
2492 :
2493 497 : 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 497 : if (!cross_section.is_serial())
2507 384 : mesh.delete_remote_elements();
2508 :
2509 : // We know a priori how many elements we'll need
2510 497 : mesh.reserve_elem(nz*orig_elem);
2511 :
2512 : // For straightforward meshes we need one or two additional layers per
2513 : // element.
2514 1843 : if (cross_section.elements_begin() != cross_section.elements_end() &&
2515 1225 : (*cross_section.elements_begin())->default_order() == SECOND)
2516 111 : order = 2;
2517 497 : mesh.comm().max(order);
2518 :
2519 497 : 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 15888 : for (const auto & node : cross_section.node_ptr_range())
2525 : {
2526 43092 : for (unsigned int k=0; k != order*nz+1; ++k)
2527 : {
2528 35296 : const dof_id_type new_node_id = node->id() + k * orig_nodes;
2529 35296 : Node * my_node = mesh.query_node_ptr(new_node_id);
2530 35296 : if (!my_node)
2531 : {
2532 : std::unique_ptr<Node> new_node = Node::build
2533 36814 : (*node + (extrusion_vector * k / nz / order),
2534 1518 : new_node_id);
2535 35296 : 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 35296 : const unique_id_type uid = (k == 0) ?
2542 684 : node->unique_id() :
2543 27500 : orig_unique_ids + (k-1)*(orig_nodes + orig_elem) + node->id();
2544 :
2545 1518 : new_node->set_unique_id(uid);
2546 : #endif
2547 :
2548 35296 : cross_section_boundary_info.boundary_ids(node, ids_to_copy);
2549 35296 : boundary_info.add_node(new_node.get(), ids_to_copy);
2550 :
2551 38332 : mesh.add_node(std::move(new_node));
2552 32260 : }
2553 : }
2554 469 : }
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 497 : 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 497 : cross_section.comm().max(next_side_id);
2566 :
2567 6566 : for (const auto & elem : cross_section.element_ptr_range())
2568 : {
2569 2923 : const ElemType etype = elem->type();
2570 :
2571 : // build_extrusion currently only works on coarse meshes
2572 130 : libmesh_assert (!elem->parent());
2573 :
2574 10181 : for (unsigned int k=0; k != nz; ++k)
2575 : {
2576 6982 : std::unique_ptr<Elem> new_elem;
2577 7258 : 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 860 : case TRI3:
2615 : {
2616 1624 : new_elem = Elem::build(PRISM6);
2617 908 : new_elem->set_node(0, mesh.node_ptr(elem->node_ptr(0)->id() + (k * orig_nodes)));
2618 908 : new_elem->set_node(1, mesh.node_ptr(elem->node_ptr(1)->id() + (k * orig_nodes)));
2619 908 : new_elem->set_node(2, mesh.node_ptr(elem->node_ptr(2)->id() + (k * orig_nodes)));
2620 908 : new_elem->set_node(3, mesh.node_ptr(elem->node_ptr(0)->id() + ((k+1) * orig_nodes)));
2621 908 : new_elem->set_node(4, mesh.node_ptr(elem->node_ptr(1)->id() + ((k+1) * orig_nodes)));
2622 908 : new_elem->set_node(5, mesh.node_ptr(elem->node_ptr(2)->id() + ((k+1) * orig_nodes)));
2623 :
2624 908 : if (elem->neighbor_ptr(0) == remote_elem)
2625 0 : new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
2626 1720 : if (elem->neighbor_ptr(1) == remote_elem)
2627 0 : new_elem->set_neighbor(2, const_cast<RemoteElem *>(remote_elem));
2628 908 : 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 4544 : case QUAD4:
2700 : {
2701 8832 : new_elem = Elem::build(HEX8);
2702 4672 : new_elem->set_node(0, mesh.node_ptr(elem->node_ptr(0)->id() + (k * orig_nodes)));
2703 4672 : new_elem->set_node(1, mesh.node_ptr(elem->node_ptr(1)->id() + (k * orig_nodes)));
2704 4672 : new_elem->set_node(2, mesh.node_ptr(elem->node_ptr(2)->id() + (k * orig_nodes)));
2705 4672 : new_elem->set_node(3, mesh.node_ptr(elem->node_ptr(3)->id() + (k * orig_nodes)));
2706 4672 : new_elem->set_node(4, mesh.node_ptr(elem->node_ptr(0)->id() + ((k+1) * orig_nodes)));
2707 4672 : new_elem->set_node(5, mesh.node_ptr(elem->node_ptr(1)->id() + ((k+1) * orig_nodes)));
2708 4672 : new_elem->set_node(6, mesh.node_ptr(elem->node_ptr(2)->id() + ((k+1) * orig_nodes)));
2709 4672 : new_elem->set_node(7, mesh.node_ptr(elem->node_ptr(3)->id() + ((k+1) * orig_nodes)));
2710 :
2711 4672 : if (elem->neighbor_ptr(0) == remote_elem)
2712 0 : new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
2713 4672 : if (elem->neighbor_ptr(1) == remote_elem)
2714 0 : new_elem->set_neighbor(2, const_cast<RemoteElem *>(remote_elem));
2715 4672 : if (elem->neighbor_ptr(2) == remote_elem)
2716 0 : new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
2717 4672 : 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 1854 : case QUAD9:
2723 : {
2724 3508 : new_elem = Elem::build(HEX27);
2725 1954 : new_elem->set_node(0, mesh.node_ptr(elem->node_ptr(0)->id() + (2*k * orig_nodes)));
2726 1954 : new_elem->set_node(1, mesh.node_ptr(elem->node_ptr(1)->id() + (2*k * orig_nodes)));
2727 1954 : new_elem->set_node(2, mesh.node_ptr(elem->node_ptr(2)->id() + (2*k * orig_nodes)));
2728 1954 : new_elem->set_node(3, mesh.node_ptr(elem->node_ptr(3)->id() + (2*k * orig_nodes)));
2729 1954 : new_elem->set_node(4, mesh.node_ptr(elem->node_ptr(0)->id() + ((2*k+2) * orig_nodes)));
2730 1954 : new_elem->set_node(5, mesh.node_ptr(elem->node_ptr(1)->id() + ((2*k+2) * orig_nodes)));
2731 1954 : new_elem->set_node(6, mesh.node_ptr(elem->node_ptr(2)->id() + ((2*k+2) * orig_nodes)));
2732 1954 : new_elem->set_node(7, mesh.node_ptr(elem->node_ptr(3)->id() + ((2*k+2) * orig_nodes)));
2733 1954 : new_elem->set_node(8, mesh.node_ptr(elem->node_ptr(4)->id() + (2*k * orig_nodes)));
2734 1954 : new_elem->set_node(9, mesh.node_ptr(elem->node_ptr(5)->id() + (2*k * orig_nodes)));
2735 1954 : new_elem->set_node(10, mesh.node_ptr(elem->node_ptr(6)->id() + (2*k * orig_nodes)));
2736 1954 : new_elem->set_node(11, mesh.node_ptr(elem->node_ptr(7)->id() + (2*k * orig_nodes)));
2737 1954 : new_elem->set_node(12, mesh.node_ptr(elem->node_ptr(0)->id() + ((2*k+1) * orig_nodes)));
2738 1954 : new_elem->set_node(13, mesh.node_ptr(elem->node_ptr(1)->id() + ((2*k+1) * orig_nodes)));
2739 1954 : new_elem->set_node(14, mesh.node_ptr(elem->node_ptr(2)->id() + ((2*k+1) * orig_nodes)));
2740 1954 : new_elem->set_node(15, mesh.node_ptr(elem->node_ptr(3)->id() + ((2*k+1) * orig_nodes)));
2741 1954 : new_elem->set_node(16, mesh.node_ptr(elem->node_ptr(4)->id() + ((2*k+2) * orig_nodes)));
2742 1954 : new_elem->set_node(17, mesh.node_ptr(elem->node_ptr(5)->id() + ((2*k+2) * orig_nodes)));
2743 1954 : new_elem->set_node(18, mesh.node_ptr(elem->node_ptr(6)->id() + ((2*k+2) * orig_nodes)));
2744 1954 : new_elem->set_node(19, mesh.node_ptr(elem->node_ptr(7)->id() + ((2*k+2) * orig_nodes)));
2745 1954 : new_elem->set_node(20, mesh.node_ptr(elem->node_ptr(8)->id() + (2*k * orig_nodes)));
2746 1954 : new_elem->set_node(21, mesh.node_ptr(elem->node_ptr(4)->id() + ((2*k+1) * orig_nodes)));
2747 1954 : new_elem->set_node(22, mesh.node_ptr(elem->node_ptr(5)->id() + ((2*k+1) * orig_nodes)));
2748 1954 : new_elem->set_node(23, mesh.node_ptr(elem->node_ptr(6)->id() + ((2*k+1) * orig_nodes)));
2749 1954 : new_elem->set_node(24, mesh.node_ptr(elem->node_ptr(7)->id() + ((2*k+1) * orig_nodes)));
2750 1954 : new_elem->set_node(25, mesh.node_ptr(elem->node_ptr(8)->id() + ((2*k+2) * orig_nodes)));
2751 1954 : new_elem->set_node(26, mesh.node_ptr(elem->node_ptr(8)->id() + ((2*k+1) * orig_nodes)));
2752 :
2753 1954 : if (elem->neighbor_ptr(0) == remote_elem)
2754 0 : new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
2755 1954 : if (elem->neighbor_ptr(1) == remote_elem)
2756 0 : new_elem->set_neighbor(2, const_cast<RemoteElem *>(remote_elem));
2757 1954 : if (elem->neighbor_ptr(2) == remote_elem)
2758 0 : new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
2759 1954 : 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 7258 : new_elem->set_id(elem->id() + (k * orig_elem));
2772 7258 : 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 7258 : const unique_id_type uid = (k == 0) ?
2779 260 : elem->unique_id() :
2780 4335 : 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 7258 : if (!elem_subdomain)
2786 : // maintain the subdomain_id
2787 2714 : new_elem->subdomain_id() = elem->subdomain_id();
2788 : else
2789 : // Allow the user to choose new subdomain_ids
2790 4544 : new_elem->subdomain_id() = elem_subdomain->get_subdomain_for_layer(elem, k);
2791 :
2792 7534 : Elem * added_elem = mesh.add_elem(std::move(new_elem));
2793 :
2794 : // Copy any old boundary ids on all sides
2795 35706 : for (auto s : elem->side_index_range())
2796 : {
2797 28172 : cross_section_boundary_info.boundary_ids(elem, s, ids_to_copy);
2798 :
2799 28172 : 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 28172 : boundary_info.add_side(added_elem,
2806 27116 : 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 7258 : if (k == 0)
2823 2923 : boundary_info.add_side(added_elem, 0, next_side_id);
2824 7258 : 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 2923 : const unsigned short top_id = added_elem->dim() == 3 ?
2830 2923 : cast_int<unsigned short>(elem->n_sides()+1) : 2;
2831 : boundary_info.add_side
2832 2923 : (added_elem, top_id,
2833 2923 : cast_int<boundary_id_type>(next_side_id+1));
2834 : }
2835 6706 : }
2836 469 : }
2837 :
2838 : // Done building the mesh. Now prepare it for use.
2839 497 : 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 378 : 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 378 : const Real xavg = (xmin + xmax)/2;
2969 378 : const Real yavg = (ymin + ymax)/2;
2970 378 : const Real zavg = (zmin + zmax)/2;
2971 390 : mesh.add_point(Point(xavg,yavg,zmin), 0);
2972 390 : mesh.add_point(Point(xmax,yavg,zavg), 1);
2973 390 : mesh.add_point(Point(xavg,ymax,zavg), 2);
2974 390 : mesh.add_point(Point(xmin,yavg,zavg), 3);
2975 390 : mesh.add_point(Point(xavg,ymin,zavg), 4);
2976 390 : mesh.add_point(Point(xavg,yavg,zmax), 5);
2977 :
2978 17024 : auto add_tri = [&mesh, flip_tris](std::array<dof_id_type,3> nodes)
2979 : {
2980 3024 : auto elem = mesh.add_elem(Elem::build(TRI3));
2981 3024 : elem->set_node(0, mesh.node_ptr(nodes[0]));
2982 3024 : elem->set_node(1, mesh.node_ptr(nodes[1]));
2983 3024 : elem->set_node(2, mesh.node_ptr(nodes[2]));
2984 3024 : if (flip_tris)
2985 1168 : elem->flip(&mesh.get_boundary_info());
2986 3036 : };
2987 :
2988 378 : add_tri({0,2,1});
2989 378 : add_tri({0,3,2});
2990 378 : add_tri({0,4,3});
2991 378 : add_tri({0,1,4});
2992 378 : add_tri({5,4,1});
2993 378 : add_tri({5,3,4});
2994 378 : add_tri({5,2,3});
2995 378 : add_tri({5,1,2});
2996 :
2997 378 : mesh.prepare_for_use();
2998 378 : }
2999 :
3000 :
3001 : } // namespace libMesh
|
