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