Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
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 "SurfaceMeshGeneratorBase.h"
11 : #include "InputParameters.h"
12 : #include "MooseMeshUtils.h"
13 : #include "MeshTraversingUtils.h"
14 :
15 : #include "libmesh/mesh_generation.h"
16 : #include "libmesh/mesh.h"
17 : #include "libmesh/elem.h"
18 : #include "libmesh/remote_elem.h"
19 : #include "libmesh/string_to_enum.h"
20 :
21 : InputParameters
22 6354 : SurfaceMeshGeneratorBase::validParams()
23 : {
24 6354 : InputParameters params = MeshGenerator::validParams();
25 25416 : params.addRequiredParam<MeshGeneratorName>("input", "The mesh we want to modify");
26 : // NOTE: don't forget to suppress this if not using this base class to create subdomains
27 25416 : params.addRequiredParam<std::vector<SubdomainName>>(
28 : "new_subdomain", "The list of subdomain names to create on the supplied subdomain");
29 25416 : params.addParam<std::vector<SubdomainName>>(
30 : "included_subdomains",
31 : "A set of subdomain names or ids an element has to be previously a part of to be considered "
32 : "for modification.");
33 :
34 : // Painting/flooding parameters
35 : // Using normals of the surface element
36 19062 : params.addParam<Point>("normal",
37 12708 : Point(),
38 : "If supplied, only surface (2D) elements with normal equal to this, up to "
39 : "normal_tol, will be considered for modification");
40 31770 : params.addRangeCheckedParam<Real>("normal_tol",
41 12708 : 0.1,
42 : "normal_tol>=0 & normal_tol<=2",
43 : "Surface elements are "
44 : "only added if face_normal.normal_hat >= "
45 : "1 - normal_tol, where normal_hat = "
46 : "normal/|normal|. The normal can the normal specified by the "
47 : "parameter or by a specific mesh generator.");
48 19062 : params.addParam<bool>(
49 : "fixed_normal",
50 12708 : false,
51 : "Whether to move the normal vector as we paint/flood the geometry, or keep it "
52 : "fixed from the first element we started painting with");
53 19062 : params.addParam<bool>("check_painted_neighbor_normals",
54 12708 : false,
55 : "When examining the normal of a 2D element and comparing to the 'painting' "
56 : "normal, also check if the element neighbors in the 'painted' subdomain "
57 : "have a closer normal and accept the element into the new subdomain if so");
58 :
59 : // Parameters for handling surface elemets with flipped normals. These are unfortunately quite
60 : // common in STL files, and if we do not account for them the flooding algorithm often floods
61 : // elements all around the flipped normal element, creating a single-element patch 'hole' inside
62 : // the larger patch!
63 31770 : params.addRangeCheckedParam<Real>("flipped_normal_tol",
64 12708 : 0.1,
65 : "flipped_normal_tol>=0 & flipped_normal_tol<=2",
66 : "If 'consider_flipped_normals' is true, surface elements are "
67 : "also added if -1 * face_normal.normal_hat >= "
68 : "1 - normal_tol, where normal_hat = "
69 : "normal/|normal|");
70 19062 : params.addParam<bool>(
71 : "consider_flipped_normals",
72 12708 : false,
73 : "Whether to allow for surface elements to be considered "
74 : "when their normal is flipped with regard to the neighbor element we are painting from. A "
75 : "specific tolerance may be specified for these flipped normals using the "
76 : "'flipped_normal_tol'.");
77 19062 : params.addParam<bool>(
78 12708 : "flip_inverted_normals", false, "Whether to flip the elements with inverted normals");
79 :
80 : // Other painting / flooding parameters
81 25416 : params.addParam<std::vector<Real>>(
82 : "max_paint_size_centroids",
83 : "Maximum distance between element centroids (using a vertex average approximation) "
84 : "when painting/flooding the geometry. 0 means not applying a distance constraint, a single "
85 : "value in the vector is applied to all elements, multiple values can be specified for each "
86 : "'included_subdomains'");
87 19062 : params.addParam<bool>(
88 : "flood_elements_once",
89 12708 : false,
90 : "Whether to consider elements only once when painting/flooding the geometry");
91 19062 : params.addParam<unsigned int>(
92 : "flood_max_recursion",
93 12708 : 1e5,
94 : "Max number of recursions when flooding. Once this number is reached, the propagation is "
95 : "stopped until enough calls have returned");
96 :
97 25416 : params.addParamNamesToGroup("normal normal_tol fixed_normal check_painted_neighbor_normals",
98 : "Flooding using surface element normals");
99 25416 : params.addParamNamesToGroup("consider_flipped_normals flipped_normal_tol",
100 : "Flooding using surface element inverted normals");
101 19062 : params.addParamNamesToGroup("max_paint_size_centroids flood_elements_once", "Other flooding");
102 6354 : return params;
103 0 : }
104 :
105 116 : SurfaceMeshGeneratorBase::SurfaceMeshGeneratorBase(const InputParameters & parameters)
106 : : MeshGenerator(parameters),
107 116 : _input(getMesh("input")),
108 348 : _subdomain_names(isParamValid("new_subdomain")
109 116 : ? getParam<std::vector<SubdomainName>>("new_subdomain")
110 : : std::vector<SubdomainName>()),
111 232 : _check_subdomains(isParamValid("included_subdomains")),
112 116 : _included_subdomain_ids(std::vector<subdomain_id_type>()),
113 376 : _using_normal(isParamSetByUser("normal") || isParamValid("_using_normal")),
114 276 : _normal(isParamSetByUser("normal")
115 408 : ? Point(getParam<Point>("normal") / getParam<Point>("normal").norm())
116 260 : : getParam<Point>("normal")),
117 232 : _normal_tol(getParam<Real>("normal_tol")),
118 232 : _flipped_normal_tol(getParam<Real>("flipped_normal_tol")),
119 232 : _fixed_flooding_normal(getParam<bool>("fixed_normal")),
120 232 : _consider_flipped_normals(getParam<bool>("consider_flipped_normals")),
121 232 : _flip_inverted_normals(getParam<bool>("flip_inverted_normals")),
122 232 : _has_max_distance_criterion(isParamSetByUser("max_paint_size_centroids")),
123 232 : _flood_only_once(getParam<bool>("flood_elements_once")),
124 232 : _flood_max_recursion(getParam<unsigned int>("flood_max_recursion")),
125 580 : _check_painted_neighor_normals(getParam<bool>("check_painted_neighbor_normals"))
126 : {
127 116 : }
128 :
129 : void
130 116 : SurfaceMeshGeneratorBase::setup(MeshBase & mesh)
131 : {
132 : // To know the dimension of the mesh
133 116 : if (!mesh.preparation().has_cached_elem_data)
134 80 : mesh.cache_elem_data();
135 :
136 : // Get the subdomain ids from the names
137 348 : if (isParamValid("included_subdomains"))
138 : {
139 : // check that the subdomains exist in the mesh
140 16 : const auto & subdomains = getParam<std::vector<SubdomainName>>("included_subdomains");
141 16 : for (const auto & name : subdomains)
142 8 : if (!MooseMeshUtils::hasSubdomainName(mesh, name))
143 0 : paramError("included_subdomains", "The block '", name, "' was not found in the mesh");
144 :
145 8 : _included_subdomain_ids = MooseMeshUtils::getSubdomainIDs(mesh, subdomains);
146 : }
147 :
148 : // Set up the max elem-to-elem distance map
149 116 : if (_has_max_distance_criterion)
150 : {
151 48 : const auto & max_dists = getParam<std::vector<Real>>("max_paint_size_centroids");
152 24 : if (max_dists.size() != _included_subdomain_ids.size() && max_dists.size() != 1)
153 0 : paramError("max_paint_size_centroids",
154 : "Maximum distance should be specified uniformly for all subdomains (1 value) or a "
155 : "value for each 'included_subdomains'");
156 :
157 24 : if (_included_subdomain_ids.size())
158 0 : for (const auto i : index_range(_included_subdomain_ids))
159 : // Single value is applied for all subdomains
160 : // 0 is translated to a very big number which therefore won't impose the criterion
161 0 : _max_elem_distance[_included_subdomain_ids[i]] =
162 0 : (max_dists.size() == 1)
163 0 : ? max_dists[0]
164 0 : : (max_dists[i] > 0 ? max_dists[i]
165 0 : : std::pow(std::numeric_limits<Real>::max(), 0.3));
166 : else
167 24 : _max_elem_distance[static_cast<subdomain_id_type>(-1)] = max_dists[0];
168 : }
169 :
170 : // Clear the maps used to count visits
171 116 : _visited.clear();
172 116 : _acted_upon_once.clear();
173 116 : }
174 :
175 : void
176 16778 : SurfaceMeshGeneratorBase::flood(Elem * const elem,
177 : const Point & base_normal,
178 : const Elem & starting_elem,
179 : const subdomain_id_type & sub_id,
180 : MeshBase & mesh)
181 : {
182 : // Avoid re-considering the same elements
183 33476 : if (elem == nullptr || elem == remote_elem ||
184 33476 : (_visited[sub_id].find(elem) != _visited[sub_id].end()))
185 13082 : return;
186 7352 : if (_flood_only_once && _acted_upon_once.count(elem))
187 1584 : return;
188 :
189 : // Skip if element is not in specified subdomains
190 5768 : if (_check_subdomains &&
191 0 : !MeshTraversingUtils::elementSubdomainIdInList(elem, _included_subdomain_ids))
192 0 : return;
193 :
194 5768 : _visited[sub_id].insert(elem);
195 5768 : auto elem_normal = get2DElemNormal(elem);
196 :
197 5768 : bool criterion_met = false;
198 5768 : if (elementSatisfiesRequirements(elem, base_normal, starting_elem, elem_normal))
199 3696 : criterion_met = true;
200 : // Additional heuristic based on neighbor normal vs element normal
201 2072 : else if (_check_painted_neighor_normals)
202 : {
203 : // Try to flood from each side with the same subdomain
204 0 : for (const auto neighbor : make_range(elem->n_sides()))
205 0 : if (elem->neighbor_ptr(neighbor) &&
206 0 : (elem->neighbor_ptr(neighbor)->subdomain_id() == sub_id) &&
207 0 : elementSatisfiesRequirements(
208 0 : elem, get2DElemNormal(elem->neighbor_ptr(neighbor)), starting_elem, elem_normal))
209 0 : criterion_met = true;
210 : }
211 :
212 5768 : if (!criterion_met)
213 2072 : return;
214 :
215 : // Modify the subdomain
216 : // TODO: consider working with base class attributes rather than arguments
217 3696 : actOnElem(elem, base_normal, sub_id, mesh);
218 :
219 : // We don't want to remove the element from consideration too early
220 3696 : _acted_upon_once.insert(elem);
221 :
222 : // Flip the element if needed
223 : // NOTE: we have already passed "elementSatisfiesRequirements" here
224 3696 : if (_consider_flipped_normals && base_normal * elem_normal < 0)
225 : {
226 200 : elem_normal *= -1;
227 200 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
228 :
229 200 : if (_flip_inverted_normals)
230 200 : elem->flip(&boundary_info);
231 : }
232 :
233 18480 : for (const auto neighbor : make_range(elem->n_sides()))
234 : {
235 14784 : _flood_recursion_count++;
236 : // Flood to the neighboring elements
237 14784 : if (_flood_recursion_count < _flood_max_recursion)
238 14784 : flood(elem->neighbor_ptr(neighbor),
239 14784 : _fixed_flooding_normal ? base_normal : elem_normal,
240 : starting_elem,
241 : sub_id,
242 : mesh);
243 14784 : _flood_recursion_count--;
244 : }
245 : }
246 :
247 : bool
248 13288 : SurfaceMeshGeneratorBase::elementSatisfiesRequirements(const Elem * const elem,
249 : const Point & desired_normal,
250 : const Elem & base_elem,
251 : const Point & face_normal) const
252 : {
253 : // False if element is not in specified subdomains
254 13288 : if (_check_subdomains &&
255 0 : !MeshTraversingUtils::elementSubdomainIdInList(elem, _included_subdomain_ids))
256 0 : return false;
257 :
258 : // False if normal does not meet criteria
259 26576 : if (_using_normal &&
260 13288 : (!MeshTraversingUtils::normalsWithinTol(desired_normal, face_normal, _normal_tol) &&
261 7328 : (!_consider_flipped_normals ||
262 13656 : !MeshTraversingUtils::normalsWithinTol(desired_normal, -face_normal, _flipped_normal_tol))))
263 7080 : return false;
264 :
265 : // False if exceeding the patch size
266 6208 : if (_has_max_distance_criterion)
267 : {
268 : // The subdomain from which the element to paint over comes from is used to find the limitation
269 : // on the radius, which will effectively be applied onto the new subdomains (to which base_elem
270 : // already belongs)
271 : const auto max_dsq =
272 4080 : _check_subdomains
273 4080 : ? MathUtils::pow(libmesh_map_find(_max_elem_distance, elem->subdomain_id()), 2)
274 4080 : : MathUtils::pow(
275 4080 : libmesh_map_find(_max_elem_distance, static_cast<subdomain_id_type>(-1)), 2);
276 4080 : if ((elem->vertex_average() - base_elem.vertex_average()).norm_sq() > max_dsq)
277 768 : return false;
278 : }
279 :
280 5440 : return true;
281 : }
282 :
283 : Point
284 9104 : SurfaceMeshGeneratorBase::get2DElemNormal(const Elem * const elem) const
285 : {
286 : mooseAssert(elem->dim() == 2, "Should be a 2D element");
287 : mooseAssert(elem->default_order() == FIRST, "Should be a first order element");
288 : // TODO: add support for non-planar element
289 : // TODO: optimize to avoid the allocation
290 9104 : if (elem->n_nodes() > 3)
291 : {
292 9104 : std::vector<Point> vec_pts(elem->n_nodes());
293 45520 : for (const auto & node : elem->node_ref_range())
294 36416 : vec_pts.push_back(node);
295 9104 : if (!MooseMeshUtils::isCoPlanar(vec_pts))
296 9024 : mooseDoOnce(mooseWarning("A non planar element was detected. Normal is approximated with the "
297 : "first 3 nodes at time."););
298 9104 : }
299 :
300 9104 : const auto & p1 = elem->point(0);
301 9104 : const auto & p2 = elem->point(1);
302 9104 : const auto & p3 = elem->point(2);
303 9104 : auto elem_normal = (p1 - p2).cross(p1 - p3);
304 9104 : if (elem_normal.norm_sq() == 0)
305 : {
306 0 : mooseWarning("Colinear nodes on element " + std::to_string(elem->id()) + ", using 0 normal");
307 0 : return elem_normal;
308 : }
309 9104 : return elem_normal.unit();
310 : }
|
