11 #include "axom/core.hpp"
12 #include "axom/fmt.hpp"
22 std::string msg = axom::fmt::format(
"Opening mesh file: '{0}'", mesh_file);
26 serac::logger::flush();
27 if (!axom::utilities::filesystem::pathExists(mesh_file)) {
28 msg = axom::fmt::format(
"Given mesh file does not exist: '{0}'", mesh_file);
34 mfem::named_ifgzstream imesh(mesh_file);
37 serac::logger::flush();
38 std::string err_msg = axom::fmt::format(
"Can not open mesh file: '{0}'", mesh_file);
39 SLIC_ERROR_ROOT(err_msg);
42 return mfem::Mesh{imesh, 1, 1,
true};
52 int num_vertices = mesh.GetNV();
53 int dim = mesh.SpaceDimension();
55 mfem::Vector vertices;
56 mesh.GetVertices(vertices);
57 mfem::Vector vertex(dim);
58 for (
int i = 0; i < num_vertices; i++) {
59 for (
int d = 0; d < dim; d++) {
60 vertex(d) = vertices[d * num_vertices + i];
63 double L1_norm = vertex.Norml1();
64 double L2_norm = vertex.Norml2();
65 vertex *= (L2_norm < 1.0e-6) ? 0.0 : (L1_norm / L2_norm);
67 for (
int d = 0; d < dim; d++) {
68 vertices[d * num_vertices + i] = vertex(d);
71 mesh.SetVertices(vertices);
76 static constexpr
int dim = 2;
77 static constexpr
int num_elems = 4;
78 static constexpr
int num_vertices = 5;
79 static constexpr
int num_boundary_elements = 4;
81 static constexpr
double vertices[num_vertices][dim] = {{0, 0}, {1, 0}, {0, 1}, {-1, 0}, {0, -1}};
82 static constexpr
int triangles[num_elems][3] = {{1, 2, 0}, {2, 3, 0}, {3, 4, 0}, {4, 1, 0}};
83 static constexpr
int segments[num_elems][2] = {{1, 2}, {2, 3}, {3, 4}, {4, 1}};
85 auto mesh = mfem::Mesh(dim, num_vertices, num_elems, num_boundary_elements);
87 for (
auto vertex : vertices) {
88 mesh.AddVertex(vertex);
90 for (
auto triangle : triangles) {
91 mesh.AddTriangle(triangle);
93 for (
auto segment : segments) {
94 mesh.AddBdrSegment(segment);
96 mesh.FinalizeTriMesh();
98 while (mesh.GetNE() < (0.5 * approx_number_of_elements)) {
99 mesh.UniformRefinement();
109 static constexpr
int dim = 3;
110 static constexpr
int num_elems = 8;
111 static constexpr
int num_vertices = 7;
112 static constexpr
int num_boundary_elements = 8;
114 static constexpr
double vertices[num_vertices][dim] = {{0, 0, 0}, {-1, 0, 0}, {0, 1, 0}, {0, 0, -1},
115 {0, 0, 1}, {0, -1, 0}, {1, 0, 0}};
116 static constexpr
int triangles[num_elems][3] = {{4, 5, 6}, {4, 6, 2}, {4, 2, 1}, {4, 1, 5},
117 {5, 1, 3}, {5, 3, 6}, {3, 1, 2}, {6, 3, 2}};
118 static constexpr
int tetrahedra[num_elems][4] = {{0, 4, 5, 6}, {0, 4, 6, 2}, {0, 4, 2, 1}, {0, 4, 1, 5},
119 {0, 5, 1, 3}, {0, 5, 3, 6}, {0, 3, 1, 2}, {0, 6, 3, 2}};
121 auto mesh = mfem::Mesh(dim, num_vertices, num_elems, num_boundary_elements);
123 for (
auto vertex : vertices) {
124 mesh.AddVertex(vertex);
126 for (
auto tetrahedron : tetrahedra) {
127 mesh.AddTet(tetrahedron);
129 for (
auto triangle : triangles) {
130 mesh.AddBdrTriangle(triangle);
133 mesh.FinalizeTetMesh();
135 while (mesh.GetNE() < (0.25 * approx_number_of_elements)) {
136 mesh.UniformRefinement();
146 return mfem::Mesh::MakeCartesian2D(elements_in_x, elements_in_y, mfem::Element::QUADRILATERAL,
true, size_x, size_y);
149 mfem::Mesh
buildCuboidMesh(
int elements_in_x,
int elements_in_y,
int elements_in_z,
double size_x,
double size_y,
152 return mfem::Mesh::MakeCartesian3D(elements_in_x, elements_in_y, elements_in_z, mfem::Element::HEXAHEDRON, size_x,
156 mfem::Mesh
buildCylinderMesh(
int radial_refinement,
int elements_lengthwise,
double radius,
double height)
158 static constexpr
int dim = 2;
159 static constexpr
int num_vertices = 17;
160 static constexpr
int num_elems = 12;
161 static constexpr
int num_boundary_elements = 8;
169 constexpr
double a = 1.3;
170 static constexpr
double vertices[num_vertices][dim] = {{0.0000000000000000, 0.0000000000000000},
171 {0.5773502691896258, 0.0000000000000000},
172 {0.4082482904638630 * a, 0.4082482904638630 * a},
173 {0.0000000000000000, 0.5773502691896258},
174 {-0.4082482904638630 * a, 0.4082482904638630 * a},
175 {-0.5773502691896258, 0.0000000000000000},
176 {-0.4082482904638630 * a, -0.4082482904638630 * a},
177 {0.0000000000000000, -0.5773502691896258},
178 {0.4082482904638630 * a, -0.4082482904638630 * a},
179 {1.000000000000000, 0.0000000000000000},
180 {0.7071067811865475, 0.7071067811865475},
181 {0.0000000000000000, 1.000000000000000},
182 {-0.707106781186548, 0.7071067811865475},
183 {-1.000000000000000, 0.0000000000000000},
184 {-0.707106781186548, -0.707106781186548},
185 {0.0000000000000000, -1.000000000000000},
186 {0.7071067811865475, -0.707106781186548}};
188 static constexpr
int elems[num_elems][4] = {{0, 1, 2, 3}, {0, 3, 4, 5}, {0, 5, 6, 7}, {0, 7, 8, 1},
189 {1, 9, 10, 2}, {2, 10, 11, 3}, {3, 11, 12, 4}, {4, 12, 13, 5},
190 {5, 13, 14, 6}, {6, 14, 15, 7}, {7, 15, 16, 8}, {8, 16, 9, 1}};
192 static constexpr
int boundary_elems[num_boundary_elements][2] = {{9, 10}, {10, 11}, {11, 12}, {12, 13},
193 {13, 14}, {14, 15}, {15, 16}, {16, 9}};
195 auto mesh = mfem::Mesh(dim, num_vertices, num_elems, num_boundary_elements);
197 for (
auto vertex : vertices) {
198 mesh.AddVertex(vertex);
200 for (
auto elem : elems) {
203 for (
auto boundary_elem : boundary_elems) {
204 mesh.AddBdrSegment(boundary_elem);
207 for (
int i = 0; i < radial_refinement; i++) {
208 mesh.UniformRefinement();
215 int n = mesh.GetNV();
217 mfem::Vector new_vertices;
218 mesh.GetVertices(new_vertices);
219 mfem::Vector vertex(dim);
220 for (
int i = 0; i < n; i++) {
221 for (
int d = 0; d < dim; d++) {
222 vertex(d) = new_vertices[d * n + i];
226 double theta =
atan2(vertex(1), vertex(0));
227 double phi = fmod(theta + M_PI, M_PI_4);
228 vertex *= radius * (
cos(phi) + (-1.0 +
sqrt(2.0)) *
sin(phi));
230 for (
int d = 0; d < dim; d++) {
231 new_vertices[d * n + i] = vertex(d);
234 mesh.SetVertices(new_vertices);
237 return mfem::Mesh(*mfem::Extrude2D(&mesh, elements_lengthwise, height));
241 mfem::Mesh
buildRing(
int radial_refinement,
double inner_radius,
double outer_radius,
double total_angle,
int sectors)
243 using index_type = int;
244 using size_type = std::vector<index_type>::size_type;
246 static constexpr
int dim = 2;
248 SLIC_ASSERT_MSG(total_angle > 0.,
"only positive angles supported");
251 total_angle =
std::min(total_angle, 2. * M_PI);
252 const double angle = total_angle / sectors;
254 auto num_elems =
static_cast<size_type
>(sectors);
255 auto num_vertices_ring =
static_cast<size_type
>((total_angle == 2. * M_PI) ? sectors : sectors + 1);
256 auto num_vertices = num_vertices_ring * 2;
257 auto num_boundary_elements = num_elems * 2;
259 SLIC_ERROR_ROOT_IF(outer_radius <= inner_radius,
260 "Outer radius is smaller than inner radius while building a cylinder mesh.");
262 std::vector<std::vector<double>> vertices(
static_cast<size_type
>(num_vertices), std::vector<double>(dim, 0.));
263 for (size_type i = 0; i < num_vertices_ring; i++) {
264 double s =
sin(angle *
static_cast<double>(i));
265 double c =
cos(angle *
static_cast<double>(i));
266 vertices[i][0] = inner_radius * c;
267 vertices[i][1] = inner_radius * s;
269 vertices[i + num_vertices_ring][0] = outer_radius * c;
270 vertices[i + num_vertices_ring][1] = outer_radius * s;
273 std::vector<std::vector<index_type>> elems(
static_cast<size_type
>(num_elems), std::vector<index_type>(4, 0));
274 std::vector<std::vector<index_type>> boundary_elems(
static_cast<size_type
>(num_boundary_elements),
275 std::vector<index_type>(2, 0));
276 for (size_type i = 0; i < num_elems; i++) {
277 elems[i][0] =
static_cast<index_type
>(i);
278 elems[i][1] =
static_cast<index_type
>(num_vertices_ring + i);
279 elems[i][2] =
static_cast<index_type
>(num_vertices_ring + (i + 1) % (num_vertices_ring));
280 elems[i][3] =
static_cast<index_type
>((i + 1) % num_vertices_ring);
283 boundary_elems[i][0] = elems[i][3];
284 boundary_elems[i][1] = elems[i][0];
287 boundary_elems[i + num_elems][0] = elems[i][1];
288 boundary_elems[i + num_elems][1] = elems[i][2];
291 mfem::Mesh mesh(dim,
static_cast<int>(num_vertices),
static_cast<int>(num_elems),
292 static_cast<int>(num_boundary_elements));
294 for (
auto vertex : vertices) {
295 mesh.AddVertex(vertex.data());
297 for (
auto elem : elems) {
298 mesh.AddQuad(elem[0], elem[1], elem[2], elem[3]);
300 for (
auto boundary_elem : boundary_elems) {
301 mesh.AddBdrSegment(boundary_elem[0], boundary_elem[1]);
304 for (
int i = 0; i < radial_refinement; i++) {
305 mesh.UniformRefinement();
312 int n = mesh.GetNV();
314 mfem::Vector new_vertices;
315 mesh.GetVertices(new_vertices);
316 mfem::Vector vertex(dim);
317 for (
int i = 0; i < n; i++) {
318 for (
int d = 0; d < dim; d++) {
319 vertex(d) = new_vertices[d * n + i];
324 double theta =
atan2(vertex(1), vertex(0));
325 double phi = fmod(theta + 2. * M_PI, angle);
331 double factor =
cos(fabs(0.5 * angle - phi)) /
cos(0.5 * angle);
334 for (
int d = 0; d < dim; d++) {
335 new_vertices[d * n + i] = vertex(d);
338 mesh.SetVertices(new_vertices);
344 mfem::Mesh
buildRingMesh(
int radial_refinement,
double inner_radius,
double outer_radius,
double total_angle,
347 return buildRing(radial_refinement, inner_radius, outer_radius, total_angle, sectors);
351 double outer_radius,
double height,
double total_angle,
int sectors)
353 auto mesh =
buildRing(radial_refinement, inner_radius, outer_radius, total_angle, sectors);
354 return mfem::Mesh(*mfem::Extrude2D(&mesh, elements_lengthwise, height));
359 std::size_t vertical_divisions,
double inner_radius,
double outer_radius,
362 constexpr
int dim = 3;
365 mfem::Mesh mesh = mfem::Mesh::MakeCartesian3D(
static_cast<int>(radial_divisions),
static_cast<int>(angular_divisions),
366 static_cast<int>(vertical_divisions), mfem::Element::HEXAHEDRON);
368 int num_vertices = mesh.GetNV();
369 mfem::Vector new_vertices;
370 mesh.GetVertices(new_vertices);
371 mfem::Vector vertex(dim);
372 for (
int i = 0; i < num_vertices; i++) {
373 for (
int d = 0; d < dim; d++) {
374 vertex(d) = new_vertices[d * num_vertices + i];
378 double r = inner_radius + (outer_radius - inner_radius) * vertex[0];
379 double theta = vertex[1] * M_PI_2;
380 vertex(0) = r *
cos(theta);
381 vertex(1) = r *
sin(theta);
382 vertex(2) = vertex(2) * height;
384 for (
int d = 0; d < dim; d++) {
385 new_vertices[d * num_vertices + i] = vertex(d);
388 mesh.SetVertices(new_vertices);
397 container.addInt(
"ser_ref_levels",
"Number of times to refine the mesh uniformly in serial.").defaultValue(0);
398 container.addInt(
"par_ref_levels",
"Number of times to refine the mesh uniformly in parallel.").defaultValue(0);
401 container.addString(
"type",
"Type of mesh").required().validValues({
"ball",
"box",
"disk",
"file"});
404 container.addString(
"mesh",
"Path to Mesh file");
407 auto& elements = container.addStruct(
"elements");
409 elements.addInt(
"x",
"x-dimension");
410 elements.addInt(
"y",
"y-dimension");
411 elements.addInt(
"z",
"z-dimension");
413 auto&
size = container.addStruct(
"size");
415 size.addDouble(
"x",
"Size in the x-dimension");
416 size.addDouble(
"y",
"Size in the y-dimension");
417 size.addDouble(
"z",
"Size in the z-dimension");
420 container.addInt(
"approx_elements",
"Approximate number of elements in an n-ball mesh");
425 std::optional<mfem::Mesh> serial_mesh;
427 if (
const auto file_opts = std::get_if<FileInputOptions>(&options.
extra_options)) {
428 SLIC_ERROR_ROOT_IF(file_opts->absolute_mesh_file_name.empty(),
429 "Absolute path to mesh file was not configured, did you forget to call findMeshFilePath?");
431 }
else if (
const auto box_opts = std::get_if<BoxInputOptions>(&options.
extra_options)) {
432 const auto& elems = box_opts->elements;
433 const auto& sizes = box_opts->overall_size;
434 if (elems.size() == 2) {
435 serial_mesh.emplace(
buildRectangleMesh(elems.at(0), elems.at(1), sizes.at(0), sizes.at(1)));
438 buildCuboidMesh(elems.at(0), elems.at(1), elems.at(2), sizes.at(0), sizes.at(1), sizes.at(2)));
440 }
else if (
const auto ball_opts = std::get_if<NBallInputOptions>(&options.
extra_options)) {
441 if (ball_opts->dimension == 2) {
442 serial_mesh.emplace(
buildDiskMesh(ball_opts->approx_elements));
444 serial_mesh.emplace(
buildBallMesh(ball_opts->approx_elements));
448 SLIC_ERROR_ROOT_IF(!serial_mesh,
"Mesh input options were invalid");
453 const int refine_parallel,
const MPI_Comm comm)
456 for (
int lev = 0; lev < refine_serial; lev++) {
457 serial_mesh.UniformRefinement();
461 auto parallel_mesh = std::make_unique<mfem::ParMesh>(comm, serial_mesh);
462 for (
int lev = 0; lev < refine_parallel; lev++) {
463 parallel_mesh->UniformRefinement();
466 parallel_mesh->EnsureNodes();
467 parallel_mesh->ExchangeFaceNbrData();
469 return parallel_mesh;
477 int ser_ref = base[
"ser_ref_levels"];
478 int par_ref = base[
"par_ref_levels"];
481 std::string mesh_type = base[
"type"];
482 if (mesh_type ==
"box") {
483 auto elements_input = base[
"elements"];
484 bool z_present = elements_input.contains(
"z");
486 std::vector<int> elements(z_present ? 3 : 2);
487 elements[0] = elements_input[
"x"];
488 elements[1] = elements_input[
"y"];
489 if (z_present) elements[2] = elements_input[
"z"];
491 std::vector<double> overall_size;
492 if (base.contains(
"size")) {
493 auto size_input = base[
"size"];
494 overall_size.push_back(size_input[
"x"]);
495 overall_size.push_back(size_input[
"y"]);
497 if (size_input.contains(
"z")) {
498 overall_size.push_back(size_input[
"z"]);
501 overall_size = std::vector<double>(elements.size(), 1.);
505 }
else if (mesh_type ==
"disk" || mesh_type ==
"ball") {
506 int approx_elements = base[
"approx_elements"];
508 if (mesh_type ==
"disk") {
512 }
else if (mesh_type ==
"file") {
513 std::string mesh_path = base[
"mesh"];
518 serac::logger::flush();
519 std::string err_msg = axom::fmt::format(
"Specified type not supported: '{0}'", mesh_type);
520 SLIC_ERROR_ROOT(err_msg);
This file contains the all the necessary functions and macros required for logging as well as a helpe...
This file contains helper functions for importing and managing various mesh objects.
std::unique_ptr< mfem::ParMesh > buildParallelMesh(const InputOptions &options, const MPI_Comm comm)
Constructs an MFEM parallel mesh from a set of input options.
std::unique_ptr< mfem::ParMesh > refineAndDistribute(mfem::Mesh &&serial_mesh, const int refine_serial, const int refine_parallel, const MPI_Comm comm)
Finalizes a serial mesh into a refined parallel mesh.
Accelerator functionality.
SERAC_HOST_DEVICE auto atan2(dual< gradient_type > y, dual< gradient_type > x)
implementation of atan2 for dual numbers
mfem::Mesh buildHollowCylinderMesh(int radial_refinement, int elements_lengthwise, double inner_radius, double outer_radius, double height, double total_angle, int sectors)
Constructs a 3D MFEM mesh of a hollow cylinder.
SERAC_HOST_DEVICE auto sin(dual< gradient_type > a)
implementation of sine for dual numbers
mfem::Mesh buildBallMesh(int approx_number_of_elements)
Constructs a 3D MFEM mesh of a unit ball, centered at the origin.
mfem::Mesh buildRing(int radial_refinement, double inner_radius, double outer_radius, double total_angle, int sectors)
Constructs a 2D MFEM mesh of a ring.
SERAC_HOST_DEVICE auto cos(dual< gradient_type > a)
implementation of cosine for dual numbers
SERAC_HOST_DEVICE auto min(dual< gradient_type > a, double b)
Implementation of min for dual numbers.
SERAC_HOST_DEVICE auto sqrt(dual< gradient_type > x)
implementation of square root for dual numbers
constexpr SERAC_HOST_DEVICE int size(const tensor< T, n... > &)
returns the total number of stored values in a tensor
mfem::Mesh buildMeshFromFile(const std::string &mesh_file)
Constructs an MFEM mesh from a file.
mfem::Mesh buildCylinderMesh(int radial_refinement, int elements_lengthwise, double radius, double height)
Constructs a 3D MFEM mesh of a cylinder.
mfem::Mesh buildRectangleMesh(int elements_in_x, int elements_in_y, double size_x, double size_y)
Constructs a 2D MFEM mesh of a rectangle.
void squish(mfem::Mesh &mesh)
a transformation from the unit disk/sphere (in L1 norm) to a unit disk/sphere (in L2 norm)
mfem::Mesh buildRingMesh(int radial_refinement, double inner_radius, double outer_radius, double total_angle, int sectors)
Constructs a 2D MFEM mesh of a ring.
mfem::Mesh buildCuboidMesh(int elements_in_x, int elements_in_y, int elements_in_z, double size_x, double size_y, double size_z)
Constructs a 3D MFEM mesh of a cuboid.
mfem::Mesh buildDiskMesh(int approx_number_of_elements)
Constructs a 2D MFEM mesh of a unit disk, centered at the origin.
mfem::Mesh build_hollow_quarter_cylinder(std::size_t radial_divisions, std::size_t angular_divisions, std::size_t vertical_divisions, double inner_radius, double outer_radius, double height)
Constructs an MFEM mesh of a hollow cylinder restricted to the first orthant.
Helper functions for exiting Serac cleanly.