Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://www.mooseframework.org
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "Boundary2DDelaunayGenerator.h"
11 :
12 : #include "MooseMeshUtils.h"
13 : #include "CastUniquePointer.h"
14 :
15 : #include "libmesh/boundary_info.h"
16 : #include "libmesh/poly2tri_triangulator.h"
17 : #include "libmesh/mesh_triangle_holes.h"
18 : #include "libmesh/mesh_modification.h"
19 : #include "libmesh/parallel_algebra.h"
20 :
21 : registerMooseObject("MooseApp", Boundary2DDelaunayGenerator);
22 :
23 : InputParameters
24 3421 : Boundary2DDelaunayGenerator::validParams()
25 : {
26 3421 : InputParameters params = SurfaceDelaunayGeneratorBase::validParams();
27 3421 : params += LevelSetMeshingHelper::validParams();
28 :
29 6842 : params.addClassDescription(
30 : "Mesh generator that convert a 2D surface given as one or a few boundaries of a 3D mesh into "
31 : "a 2D mesh using Delaunay triangulation.");
32 13684 : params.addRequiredParam<MeshGeneratorName>("input", "The mesh we want to modify");
33 13684 : params.addRequiredParam<std::vector<BoundaryName>>("boundary_names", "The boundaries to be used");
34 10263 : params.addParam<std::vector<std::vector<BoundaryName>>>(
35 : "hole_boundary_names",
36 6842 : std::vector<std::vector<BoundaryName>>(),
37 : "The optional boundaries to be used as the holes in the mesh during triangulation. Note that "
38 : "this is a vector of vectors, which allows each hole to be defined as a combination of "
39 : "multiple boundaries.");
40 :
41 10263 : params.addParam<BoundaryName>(
42 : "output_external_boundary_name",
43 : "",
44 : "The optional name of the external boundary of the mesh to generate. If not provided, the "
45 : "external boundary will be have a trivial id of 0.");
46 :
47 3421 : return params;
48 0 : }
49 :
50 180 : Boundary2DDelaunayGenerator::Boundary2DDelaunayGenerator(const InputParameters & parameters)
51 : : SurfaceDelaunayGeneratorBase(parameters),
52 : LevelSetMeshingHelper(parameters),
53 180 : _input(getMesh("input")),
54 360 : _boundary_names(getParam<std::vector<BoundaryName>>("boundary_names")),
55 360 : _hole_boundary_names(getParam<std::vector<std::vector<BoundaryName>>>("hole_boundary_names")),
56 540 : _output_external_boundary_name(getParam<BoundaryName>("output_external_boundary_name"))
57 : {
58 180 : }
59 :
60 : std::unique_ptr<MeshBase>
61 100 : Boundary2DDelaunayGenerator::generate()
62 : {
63 100 : std::unique_ptr<MeshBase> mesh_3d = std::move(_input);
64 :
65 : // Generate a new 2D block based on the sidesets
66 100 : const auto new_block_id = MooseMeshUtils::getNextFreeSubdomainID(*mesh_3d);
67 : try
68 : {
69 209 : MooseMeshUtils::createSubdomainFromSidesets(
70 303 : *mesh_3d, _boundary_names, new_block_id, SubdomainName(), type());
71 : }
72 3 : catch (MooseException & e)
73 : {
74 6 : if (((std::string)e.what()).compare(0, 12, "The sideset ") == 0)
75 6 : paramError("boundary_names", e.what());
76 : else
77 0 : mooseError(e.what());
78 0 : }
79 :
80 : // If holes are provided, we need to create new blocks for them too
81 97 : std::vector<subdomain_id_type> hole_block_ids;
82 111 : for (const auto & hole : _hole_boundary_names)
83 : {
84 17 : hole_block_ids.push_back(MooseMeshUtils::getNextFreeSubdomainID(*mesh_3d));
85 : try
86 : {
87 57 : MooseMeshUtils::createSubdomainFromSidesets(
88 54 : *mesh_3d, hole, hole_block_ids.back(), SubdomainName(), type());
89 : }
90 3 : catch (MooseException & e)
91 : {
92 6 : if (((std::string)e.what()).compare(0, 12, "The sideset ") == 0)
93 6 : paramError("hole_boundary_names", e.what());
94 : else
95 0 : mooseError(e.what());
96 0 : }
97 : }
98 :
99 : // Create a 2D mesh form the 2D block
100 94 : auto mesh_2d = buildMeshBaseObject();
101 282 : MooseMeshUtils::convertBlockToMesh(*mesh_3d, *mesh_2d, {std::to_string(new_block_id)});
102 : // If holes are provided, we need to create a 2D mesh for each hole
103 94 : std::vector<std::unique_ptr<MeshBase>> hole_meshes_2d;
104 108 : for (const auto & hole_block_id : hole_block_ids)
105 : {
106 14 : hole_meshes_2d.push_back(buildMeshBaseObject());
107 56 : MooseMeshUtils::convertBlockToMesh(
108 28 : *mesh_3d, *hole_meshes_2d.back(), {std::to_string(hole_block_id)});
109 : }
110 :
111 : // We do not need the 3D mesh anymore
112 94 : mesh_3d->clear();
113 :
114 182 : return General2DDelaunay(mesh_2d, hole_meshes_2d);
115 196 : }
116 :
117 : std::unique_ptr<MeshBase>
118 94 : Boundary2DDelaunayGenerator::General2DDelaunay(
119 : std::unique_ptr<MeshBase> & mesh_2d, std::vector<std::unique_ptr<MeshBase>> & hole_meshes_2d)
120 : {
121 : // If a level set is provided, we need to check if the nodes in the original 2D mesh match the
122 : // level set
123 94 : if (_func_level_set)
124 : {
125 843 : for (const auto & node : mesh_2d->node_ptr_range())
126 : {
127 803 : if (std::abs(levelSetEvaluator(*node)) > libMesh::TOLERANCE)
128 : {
129 6 : paramError("level_set",
130 : "The level set function does not match the nodes in the given boundary of the "
131 : "input mesh.");
132 : }
133 40 : }
134 : }
135 :
136 : // Easier to work with a TRI3 mesh
137 : // all_tri() also prepares the mesh for use
138 91 : mesh_2d->prepare_for_use();
139 91 : MeshTools::Modification::all_tri(*mesh_2d);
140 91 : const auto mesh_2d_ext_bdry = MooseMeshUtils::getNextFreeBoundaryID(*mesh_2d);
141 3083 : for (const auto & elem : mesh_2d->active_element_ptr_range())
142 11968 : for (const auto & i_side : elem->side_index_range())
143 8976 : if (elem->neighbor_ptr(i_side) == nullptr)
144 1475 : mesh_2d->get_boundary_info().add_side(elem, i_side, mesh_2d_ext_bdry);
145 :
146 : // Create a clone of the 2D mesh to be used for the 1D mesh generation
147 91 : auto mesh_2d_dummy = dynamic_pointer_cast<MeshBase>(mesh_2d->clone());
148 : // Generate a new 1D block based on the external boundary
149 91 : const auto new_block_id_1d = MooseMeshUtils::getNextFreeSubdomainID(*mesh_2d_dummy);
150 :
151 455 : MooseMeshUtils::createSubdomainFromSidesets(
152 364 : *mesh_2d_dummy, {std::to_string(mesh_2d_ext_bdry)}, new_block_id_1d, SubdomainName(), type());
153 :
154 : // Create a 1D mesh form the 1D block
155 91 : auto mesh_1d = buildMeshBaseObject();
156 273 : MooseMeshUtils::convertBlockToMesh(*mesh_2d_dummy, *mesh_1d, {std::to_string(new_block_id_1d)});
157 91 : mesh_2d_dummy->clear();
158 :
159 : // If we have holes, we need to create a 1D mesh for each hole
160 91 : std::vector<std::unique_ptr<MeshBase>> hole_meshes_1d;
161 102 : for (auto & hole_mesh_2d : hole_meshes_2d)
162 : {
163 : // As we do not need these holes for reverse projection, we do not need to convert them to TRI3
164 : // meshes, but we still need to create a 1D mesh for each hole
165 11 : hole_mesh_2d->find_neighbors();
166 11 : const auto hole_mesh_2d_ext_bdry = MooseMeshUtils::getNextFreeBoundaryID(*hole_mesh_2d);
167 55 : for (const auto & elem : hole_mesh_2d->active_element_ptr_range())
168 220 : for (const auto & i_side : elem->side_index_range())
169 176 : if (elem->neighbor_ptr(i_side) == nullptr)
170 99 : hole_mesh_2d->get_boundary_info().add_side(elem, i_side, hole_mesh_2d_ext_bdry);
171 11 : const auto new_hole_block_id_1d = MooseMeshUtils::getNextFreeSubdomainID(*hole_mesh_2d);
172 44 : MooseMeshUtils::createSubdomainFromSidesets(*hole_mesh_2d,
173 11 : {std::to_string(hole_mesh_2d_ext_bdry)},
174 : new_hole_block_id_1d,
175 22 : SubdomainName(),
176 11 : type());
177 : // Create a 1D mesh form the 1D block
178 11 : hole_meshes_1d.push_back(buildMeshBaseObject());
179 44 : MooseMeshUtils::convertBlockToMesh(
180 22 : *hole_mesh_2d, *hole_meshes_1d.back(), {std::to_string(new_hole_block_id_1d)});
181 11 : hole_mesh_2d->clear();
182 : }
183 :
184 : // Find centroid of the 2D mesh
185 91 : const Point centroid = MooseMeshUtils::meshCentroidCalculator(*mesh_2d);
186 : // calculate an average normal vector of the 2D mesh
187 91 : const Point mesh_norm = meshNormal2D(*mesh_2d);
188 : // Check the deviation of the mesh normal vector from the global average normal vector
189 91 : if (meshNormalDeviation2D(*mesh_2d, mesh_norm) > _max_angle_deviation)
190 3 : paramError("boundary_names",
191 : "The normal vector of some elements in the 2D mesh deviates too much from the "
192 3 : "global average normal vector. The maximum deviation is " +
193 6 : std::to_string(_max_angle_deviation) +
194 : ". Consider dividing the boundary into several parts to "
195 : "reduce the angle deviation.");
196 :
197 : // Move both 2d and 1d meshes to the centroid of the 2D mesh
198 88 : MeshTools::Modification::translate(*mesh_1d, -centroid(0), -centroid(1), -centroid(2));
199 88 : MeshTools::Modification::translate(*mesh_2d, -centroid(0), -centroid(1), -centroid(2));
200 : // Alsop need to translate the 1D hole meshes if applicable
201 96 : for (auto & hole_mesh_1d : hole_meshes_1d)
202 8 : MeshTools::Modification::translate(*hole_mesh_1d, -centroid(0), -centroid(1), -centroid(2));
203 :
204 : // Calculate the Euler angles to rotate the meshes so that the 2D mesh is close to the XY plane
205 : // (i.e., the normal vector of the 2D mesh is aligned with the Z axis)
206 88 : const Real theta = std::acos(mesh_norm(2)) / M_PI * 180.0;
207 : const Real phi =
208 88 : (MooseUtils::absoluteFuzzyLessThan(mesh_norm(2), 1.0) ? std::atan2(mesh_norm(1), mesh_norm(0))
209 88 : : 0.0) /
210 88 : M_PI * 180.0;
211 88 : MeshTools::Modification::rotate(*mesh_1d, 90.0 - phi, theta, 0.0);
212 88 : MeshTools::Modification::rotate(*mesh_2d, 90.0 - phi, theta, 0.0);
213 : // Also rotate the 1D hole meshes if applicable
214 96 : for (auto & hole_mesh_1d : hole_meshes_1d)
215 8 : MeshTools::Modification::rotate(*hole_mesh_1d, 90.0 - phi, theta, 0.0);
216 :
217 : // Clone the 2D mesh to be used for reverse projection later
218 88 : auto mesh_2d_xyz = dynamic_pointer_cast<MeshBase>(mesh_2d->clone());
219 :
220 : // Project the 2D mesh to the XY plane so that XYDelaunay can be used
221 2256 : for (const auto & node : mesh_2d->node_ptr_range())
222 2256 : (*node)(2) = 0;
223 : // Project the 1D mesh to the XY plane as well
224 1400 : for (const auto & node : mesh_1d->node_ptr_range())
225 1400 : (*node)(2) = 0;
226 : // Also project the 1D hole meshes to the XY plane if applicable
227 96 : for (auto & hole_mesh_1d : hole_meshes_1d)
228 : {
229 72 : for (const auto & node : hole_mesh_1d->node_ptr_range())
230 72 : (*node)(2) = 0;
231 : }
232 :
233 88 : std::vector<libMesh::TriangulatorInterface::MeshedHole> meshed_holes;
234 88 : std::vector<libMesh::TriangulatorInterface::Hole *> triangulator_hole_ptrs(hole_meshes_1d.size());
235 88 : meshed_holes.reserve(hole_meshes_1d.size());
236 96 : for (auto hole_i : index_range(hole_meshes_1d))
237 : {
238 8 : hole_meshes_1d[hole_i]->prepare_for_use();
239 8 : meshed_holes.emplace_back(*hole_meshes_1d[hole_i]);
240 8 : triangulator_hole_ptrs[hole_i] = &meshed_holes.back();
241 : }
242 :
243 : // Finally, triangulation
244 : std::unique_ptr<UnstructuredMesh> mesh =
245 88 : dynamic_pointer_cast<UnstructuredMesh>(std::move(mesh_1d));
246 :
247 88 : Poly2TriTriangulator poly2tri(*mesh);
248 88 : poly2tri.triangulation_type() = TriangulatorInterface::PSLG;
249 :
250 88 : poly2tri.set_interpolate_boundary_points(0);
251 88 : poly2tri.set_refine_boundary_allowed(false);
252 88 : poly2tri.set_verify_hole_boundaries(false);
253 88 : poly2tri.desired_area() = 0;
254 88 : poly2tri.minimum_angle() = 0; // Not yet supported
255 88 : poly2tri.smooth_after_generating() = false;
256 88 : if (!triangulator_hole_ptrs.empty())
257 8 : poly2tri.attach_hole_list(&triangulator_hole_ptrs);
258 : // Future TODO: correct the area function based on the local normal vector
259 88 : if (_use_auto_area_func)
260 80 : poly2tri.set_auto_area_function(this->comm(),
261 80 : _auto_area_function_num_points,
262 80 : _auto_area_function_power,
263 80 : _auto_area_func_default_size,
264 80 : _auto_area_func_default_size_dist);
265 88 : poly2tri.triangulate();
266 :
267 : // Reverse the projection based on the original 2D mesh
268 2464 : for (const auto & node : mesh->node_ptr_range())
269 : {
270 2376 : bool node_mod = false;
271 : // Try to find the element in mesh_2d that contains the new node
272 41984 : for (const auto & elem : mesh_2d->active_element_ptr_range())
273 : {
274 41984 : if (elem->contains_point(Point((*node)(0), (*node)(1), 0.0)))
275 : {
276 : // Element id
277 2376 : const auto elem_id = elem->id();
278 : // element in xyz_in_xyz
279 2376 : const Elem & elem_xyz = *mesh_2d_xyz->elem_ptr(elem_id);
280 :
281 2376 : const Point elem_normal = elemNormal(elem_xyz);
282 2376 : const Point & elem_p = *mesh_2d_xyz->elem_ptr(elem_id)->node_ptr(0);
283 :
284 : // if the x and y values of the node is the same as the elem_p's first node, we can just
285 : // move it to that node's position
286 3176 : if (MooseUtils::absoluteFuzzyEqual((*node)(0), elem_p(0)) &&
287 800 : MooseUtils::absoluteFuzzyEqual((*node)(1), elem_p(1)))
288 : {
289 168 : (*node)(2) = elem_p(2);
290 168 : node_mod = true;
291 168 : break;
292 : }
293 : // Otherwise, we need to find a position inside the 2D element
294 : // It has the same x and y coordinates as the node in the projected mesh;
295 2208 : (*node)(2) = elem_p(2) - (((*node)(0) - elem_p(0)) * elem_normal(0) +
296 2208 : ((*node)(1) - elem_p(1)) * elem_normal(1)) /
297 2208 : elem_normal(2);
298 2208 : node_mod = true;
299 2208 : break;
300 : }
301 2376 : }
302 2376 : if (!node_mod)
303 0 : mooseError("Node not found in mesh_in_xy");
304 88 : }
305 :
306 : // Rotate the mesh back
307 88 : MeshTools::Modification::rotate(*mesh, 0.0, -theta, phi - 90.0);
308 : // Translate the mesh back
309 88 : MeshTools::Modification::translate(*mesh, centroid(0), centroid(1), centroid(2));
310 :
311 : // Correct the nodes based on the level set function
312 88 : if (_func_level_set)
313 : {
314 1040 : for (const auto & node : mesh->node_ptr_range())
315 : {
316 1000 : unsigned int iter_ct = 0;
317 4000 : while (iter_ct < _max_level_set_correction_iterations &&
318 1800 : std::abs(levelSetEvaluator(*node)) > libMesh::TOLERANCE * libMesh::TOLERANCE)
319 : {
320 1200 : levelSetCorrection(*node);
321 1200 : ++iter_ct;
322 : }
323 40 : }
324 : }
325 :
326 : // Assign the external boundary name if provided
327 88 : if (!_output_external_boundary_name.empty())
328 : {
329 : const std::vector<BoundaryID> output_boundary_id =
330 48 : MooseMeshUtils::getBoundaryIDs(*mesh, {_output_external_boundary_name}, true);
331 :
332 16 : libMesh::MeshTools::Modification::change_boundary_id(*mesh, 0, output_boundary_id[0]);
333 16 : mesh->get_boundary_info().sideset_name(output_boundary_id[0]) = _output_external_boundary_name;
334 16 : }
335 :
336 88 : mesh->unset_is_prepared();
337 :
338 176 : return mesh;
339 308 : }
|
