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