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 "SideSetsGeneratorBase.h"
11 : #include "Parser.h"
12 : #include "InputParameters.h"
13 : #include "MooseMesh.h"
14 : #include "MooseMeshUtils.h"
15 :
16 : #include "libmesh/mesh_generation.h"
17 : #include "libmesh/mesh.h"
18 : #include "libmesh/string_to_enum.h"
19 : #include "libmesh/quadrature_gauss.h"
20 : #include "libmesh/point_locator_base.h"
21 : #include "libmesh/elem.h"
22 : #include "libmesh/remote_elem.h"
23 :
24 : InputParameters
25 110269 : SideSetsGeneratorBase::validParams()
26 : {
27 110269 : InputParameters params = MeshGenerator::validParams();
28 110269 : params.addRequiredParam<MeshGeneratorName>("input", "The mesh we want to modify");
29 110269 : params.addRequiredParam<std::vector<BoundaryName>>(
30 : "new_boundary", "The list of boundary names to create on the supplied subdomain");
31 330807 : params.addParam<bool>("fixed_normal",
32 220538 : false,
33 : "This Boolean determines whether we fix our normal "
34 : "or allow it to vary to \"paint\" around curves");
35 :
36 330807 : params.addParam<bool>("replace",
37 220538 : false,
38 : "If true, replace the old sidesets. If false, the current sidesets (if "
39 : "any) will be preserved.");
40 :
41 110269 : params.addParam<std::vector<BoundaryName>>(
42 : "included_boundaries",
43 : "A set of boundary names or ids whose sides will be included in the new sidesets. A side "
44 : "is only added if it also belongs to one of these boundaries.");
45 110269 : params.addParam<std::vector<BoundaryName>>(
46 : "excluded_boundaries",
47 : "A set of boundary names or ids whose sides will be excluded from the new sidesets. A side "
48 : "is only added if does not belong to any of these boundaries.");
49 110269 : params.addParam<std::vector<SubdomainName>>(
50 : "included_subdomains",
51 : "A set of subdomain names or ids whose sides will be included in the new sidesets. A side "
52 : "is only added if the subdomain id of the corresponding element is in this set.");
53 110269 : params.addParam<std::vector<SubdomainName>>("included_neighbors",
54 : "A set of neighboring subdomain names or ids. A face "
55 : "is only added if the subdomain id of the "
56 : "neighbor is in this set");
57 330807 : params.addParam<bool>(
58 : "include_only_external_sides",
59 220538 : false,
60 : "Whether to only include external sides when considering sides to add to the sideset");
61 :
62 330807 : params.addParam<Point>("normal",
63 220538 : Point(),
64 : "If supplied, only faces with normal equal to this, up to "
65 : "normal_tol, will be added to the sidesets specified");
66 330807 : params.addRangeCheckedParam<Real>("normal_tol",
67 220538 : 0.1,
68 : "normal_tol>=0 & normal_tol<=2",
69 : "If normal is supplied then faces are "
70 : "only added if face_normal.normal_hat >= "
71 : "1 - normal_tol, where normal_hat = "
72 : "normal/|normal|");
73 110269 : params.addParam<Real>("variance", "The variance allowed when comparing normals");
74 110269 : params.deprecateParam("variance", "normal_tol", "4/01/2025");
75 :
76 : // Sideset restriction param group
77 110269 : params.addParamNamesToGroup(
78 : "included_boundaries excluded_boundaries included_subdomains included_neighbors "
79 : "include_only_external_sides normal normal_tol",
80 : "Sideset restrictions");
81 :
82 110269 : return params;
83 0 : }
84 :
85 5185 : SideSetsGeneratorBase::SideSetsGeneratorBase(const InputParameters & parameters)
86 : : MeshGenerator(parameters),
87 5185 : _input(getMesh("input")),
88 5185 : _boundary_names(std::vector<BoundaryName>()),
89 5185 : _fixed_normal(getParam<bool>("fixed_normal")),
90 5185 : _replace(getParam<bool>("replace")),
91 5185 : _check_included_boundaries(isParamValid("included_boundaries")),
92 5185 : _check_excluded_boundaries(isParamValid("excluded_boundaries")),
93 5185 : _check_subdomains(isParamValid("included_subdomains")),
94 5185 : _check_neighbor_subdomains(isParamValid("included_neighbors")),
95 5185 : _included_boundary_ids(std::vector<boundary_id_type>()),
96 5185 : _excluded_boundary_ids(std::vector<boundary_id_type>()),
97 5185 : _included_subdomain_ids(std::vector<subdomain_id_type>()),
98 5185 : _included_neighbor_subdomain_ids(std::vector<subdomain_id_type>()),
99 5185 : _include_only_external_sides(getParam<bool>("include_only_external_sides")),
100 5185 : _using_normal(isParamSetByUser("normal")),
101 10370 : _normal(_using_normal ? Point(getParam<Point>("normal") / getParam<Point>("normal").norm())
102 5185 : : getParam<Point>("normal")),
103 10370 : _normal_tol(getParam<Real>("normal_tol"))
104 : {
105 5185 : if (isParamValid("new_boundary"))
106 4675 : _boundary_names = getParam<std::vector<BoundaryName>>("new_boundary");
107 5185 : }
108 :
109 4941 : SideSetsGeneratorBase::~SideSetsGeneratorBase() {}
110 :
111 : void
112 4962 : SideSetsGeneratorBase::setup(MeshBase & mesh)
113 : {
114 : mooseAssert(_fe_face == nullptr, "FE Face has already been initialized");
115 :
116 : // To know the dimension of the mesh
117 4962 : if (!mesh.is_prepared())
118 4344 : mesh.prepare_for_use();
119 4962 : const auto dim = mesh.mesh_dimension();
120 :
121 : // Setup the FE Object so we can calculate normals
122 9924 : libMesh::FEType fe_type(Utility::string_to_enum<Order>("CONSTANT"),
123 4962 : Utility::string_to_enum<libMesh::FEFamily>("MONOMIAL"));
124 4962 : _fe_face = libMesh::FEBase::build(dim, fe_type);
125 4962 : _qface = std::make_unique<libMesh::QGauss>(dim - 1, FIRST);
126 4962 : _fe_face->attach_quadrature_rule(_qface.get());
127 : // Must always pre-request quantities you want to compute
128 4962 : _fe_face->get_normals();
129 :
130 : // Handle incompatible parameters
131 4962 : if (_include_only_external_sides && _check_neighbor_subdomains)
132 0 : paramError("include_only_external_sides", "External sides dont have neighbors");
133 :
134 4962 : if (_check_included_boundaries)
135 : {
136 239 : const auto & included_boundaries = getParam<std::vector<BoundaryName>>("included_boundaries");
137 474 : for (const auto & boundary_name : _boundary_names)
138 239 : if (std::find(included_boundaries.begin(), included_boundaries.end(), boundary_name) !=
139 478 : included_boundaries.end())
140 4 : paramError(
141 : "new_boundary",
142 : "A boundary cannot be both the new boundary and be included in the list of included "
143 : "boundaries. If you are trying to restrict an existing boundary, you must use a "
144 : "different name for 'new_boundary', delete the old boundary, and then rename the "
145 : "new boundary to the old boundary.");
146 :
147 235 : _included_boundary_ids = MooseMeshUtils::getBoundaryIDs(mesh, included_boundaries, false);
148 :
149 : // Check that the included boundary ids/names exist in the mesh
150 630 : for (const auto i : index_range(_included_boundary_ids))
151 399 : if (_included_boundary_ids[i] == Moose::INVALID_BOUNDARY_ID)
152 4 : paramError("included_boundaries",
153 : "The boundary '",
154 4 : included_boundaries[i],
155 : "' was not found within the mesh");
156 : }
157 :
158 4954 : if (_check_excluded_boundaries)
159 : {
160 20 : const auto & excluded_boundaries = getParam<std::vector<BoundaryName>>("excluded_boundaries");
161 36 : for (const auto & boundary_name : _boundary_names)
162 20 : if (std::find(excluded_boundaries.begin(), excluded_boundaries.end(), boundary_name) !=
163 40 : excluded_boundaries.end())
164 4 : paramError(
165 : "new_boundary",
166 : "A boundary cannot be both the new boundary and be excluded in the list of excluded "
167 : "boundaries.");
168 16 : _excluded_boundary_ids = MooseMeshUtils::getBoundaryIDs(mesh, excluded_boundaries, false);
169 :
170 : // Check that the excluded boundary ids/names exist in the mesh
171 28 : for (const auto i : index_range(_excluded_boundary_ids))
172 16 : if (_excluded_boundary_ids[i] == Moose::INVALID_BOUNDARY_ID)
173 4 : paramError("excluded_boundaries",
174 : "The boundary '",
175 4 : excluded_boundaries[i],
176 : "' was not found within the mesh");
177 :
178 12 : if (_check_included_boundaries)
179 : {
180 : // Check that included and excluded boundary lists do not overlap
181 4 : for (const auto & boundary_id : _included_boundary_ids)
182 4 : if (std::find(_excluded_boundary_ids.begin(), _excluded_boundary_ids.end(), boundary_id) !=
183 8 : _excluded_boundary_ids.end())
184 4 : paramError("excluded_boundaries",
185 : "'included_boundaries' and 'excluded_boundaries' lists should not overlap");
186 : }
187 : }
188 :
189 : // Get the boundary ids from the names
190 4942 : if (parameters().isParamValid("included_subdomains"))
191 : {
192 : // check that the subdomains exist in the mesh
193 4479 : const auto subdomains = getParam<std::vector<SubdomainName>>("included_subdomains");
194 9016 : for (const auto & name : subdomains)
195 4549 : if (!MooseMeshUtils::hasSubdomainName(mesh, name))
196 12 : paramError("included_subdomains", "The block '", name, "' was not found in the mesh");
197 :
198 4467 : _included_subdomain_ids = MooseMeshUtils::getSubdomainIDs(mesh, subdomains);
199 4467 : }
200 :
201 4930 : if (parameters().isParamValid("included_neighbors"))
202 : {
203 : // check that the subdomains exist in the mesh
204 3288 : const auto subdomains = getParam<std::vector<SubdomainName>>("included_neighbors");
205 6578 : for (const auto & name : subdomains)
206 3298 : if (!MooseMeshUtils::hasSubdomainName(mesh, name))
207 8 : paramError("included_neighbors", "The block '", name, "' was not found in the mesh");
208 :
209 3280 : _included_neighbor_subdomain_ids = MooseMeshUtils::getSubdomainIDs(mesh, subdomains);
210 3280 : }
211 :
212 : // We will want to Change the below code when we have more fine-grained control over advertising
213 : // what we need and how we satisfy those needs. For now we know we need to have neighbors per
214 : // #15823...and we do have an explicit `find_neighbors` call...but we don't have a
215 : // `neighbors_found` API and it seems off to do:
216 : //
217 : // if (!mesh.is_prepared())
218 : // mesh.find_neighbors()
219 4922 : }
220 :
221 : void
222 319 : SideSetsGeneratorBase::finalize()
223 : {
224 319 : _qface.reset();
225 319 : _fe_face.reset();
226 319 : }
227 :
228 : void
229 1978518 : SideSetsGeneratorBase::flood(const Elem * elem,
230 : const Point & normal,
231 : const boundary_id_type & side_id,
232 : MeshBase & mesh)
233 : {
234 3585276 : if (elem == nullptr || elem == remote_elem ||
235 3585276 : (_visited[side_id].find(elem) != _visited[side_id].end()))
236 1303976 : return;
237 :
238 : // Skip if element is not in specified subdomains
239 674542 : if (_check_subdomains && !elementSubdomainIdInList(elem, _included_subdomain_ids))
240 0 : return;
241 :
242 674542 : _visited[side_id].insert(elem);
243 :
244 : // Request to compute normal vectors
245 674542 : const std::vector<Point> & face_normals = _fe_face->get_normals();
246 :
247 4064150 : for (const auto side : make_range(elem->n_sides()))
248 : {
249 :
250 3389608 : _fe_face->reinit(elem, side);
251 : // We'll just use the normal of the first qp
252 3389608 : const Point face_normal = face_normals[0];
253 :
254 3389608 : if (!elemSideSatisfiesRequirements(elem, side, mesh, normal, face_normal))
255 3079430 : continue;
256 :
257 310178 : if (_replace)
258 1280 : mesh.get_boundary_info().remove_side(elem, side);
259 :
260 310178 : mesh.get_boundary_info().add_side(elem, side, side_id);
261 1872794 : for (const auto neighbor : make_range(elem->n_sides()))
262 : {
263 : // Flood to the neighboring elements using the current matching side normal from this
264 : // element.
265 : // This will allow us to tolerate small changes in the normals so we can "paint" around a
266 : // curve.
267 1562616 : flood(elem->neighbor_ptr(neighbor), _fixed_normal ? normal : face_normal, side_id, mesh);
268 : }
269 : }
270 : }
271 :
272 : bool
273 2093118 : SideSetsGeneratorBase::normalsWithinTol(const Point & normal_1,
274 : const Point & normal_2,
275 : const Real & tol) const
276 : {
277 2093118 : return (1.0 - normal_1 * normal_2) <= tol;
278 : }
279 :
280 : bool
281 940856 : SideSetsGeneratorBase::elementSubdomainIdInList(
282 : const Elem * const elem, const std::vector<subdomain_id_type> & subdomain_id_list) const
283 : {
284 940856 : subdomain_id_type curr_subdomain = elem->subdomain_id();
285 940856 : return std::find(subdomain_id_list.begin(), subdomain_id_list.end(), curr_subdomain) !=
286 1881712 : subdomain_id_list.end();
287 : }
288 :
289 : bool
290 17792 : SideSetsGeneratorBase::elementSideInIncludedBoundaries(const Elem * const elem,
291 : const unsigned int side,
292 : const MeshBase & mesh) const
293 : {
294 38936 : for (const auto & bid : _included_boundary_ids)
295 22032 : if (mesh.get_boundary_info().has_boundary_id(elem, side, bid))
296 888 : return true;
297 16904 : return false;
298 : }
299 :
300 : bool
301 1728 : SideSetsGeneratorBase::elementSideInExcludedBoundaries(const Elem * const elem,
302 : const unsigned int side,
303 : const MeshBase & mesh) const
304 : {
305 3440 : for (const auto bid : _excluded_boundary_ids)
306 1728 : if (mesh.get_boundary_info().has_boundary_id(elem, side, bid))
307 16 : return true;
308 1712 : return false;
309 : }
310 :
311 : bool
312 4031226 : SideSetsGeneratorBase::elemSideSatisfiesRequirements(const Elem * const elem,
313 : const unsigned int side,
314 : const MeshBase & mesh,
315 : const Point & desired_normal,
316 : const Point & face_normal)
317 : {
318 : // Skip if side has neighbor and we only want external sides
319 4031226 : if ((elem->neighbor_ptr(side) && _include_only_external_sides))
320 3074426 : return false;
321 :
322 : // Skip if side is not part of included boundaries
323 956800 : if (_check_included_boundaries && !elementSideInIncludedBoundaries(elem, side, mesh))
324 16904 : return false;
325 : // Skip if side is part of excluded boundaries
326 939896 : if (_check_excluded_boundaries && elementSideInExcludedBoundaries(elem, side, mesh))
327 16 : return false;
328 :
329 : // Skip if element does not have neighbor in specified subdomains
330 939880 : if (_check_neighbor_subdomains)
331 : {
332 566464 : const Elem * const neighbor = elem->neighbor_ptr(side);
333 : // if the neighbor does not exist, then skip this face; we only add sidesets
334 : // between existing elems if _check_neighbor_subdomains is true
335 566464 : if (!(neighbor && elementSubdomainIdInList(neighbor, _included_neighbor_subdomain_ids)))
336 531326 : return false;
337 : }
338 :
339 408554 : if (_using_normal && !normalsWithinTol(desired_normal, face_normal, _normal_tol))
340 34959 : return false;
341 :
342 373595 : return true;
343 : }
|
