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 115637 : SideSetsGeneratorBase::validParams()
26 : {
27 115637 : InputParameters params = MeshGenerator::validParams();
28 462548 : params.addRequiredParam<MeshGeneratorName>("input", "The mesh we want to modify");
29 462548 : params.addRequiredParam<std::vector<BoundaryName>>(
30 : "new_boundary", "The list of boundary names to create on the supplied subdomain");
31 346911 : params.addParam<bool>("fixed_normal",
32 231274 : false,
33 : "This Boolean determines whether we fix our normal "
34 : "or allow it to vary to \"paint\" around curves");
35 :
36 346911 : params.addParam<bool>("replace",
37 231274 : false,
38 : "If true, replace the old sidesets. If false, the current sidesets (if "
39 : "any) will be preserved.");
40 :
41 462548 : 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 462548 : 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 462548 : 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 462548 : 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 346911 : params.addParam<bool>(
58 : "include_only_external_sides",
59 231274 : false,
60 : "Whether to only include external sides when considering sides to add to the sideset");
61 :
62 346911 : params.addParam<Point>("normal",
63 231274 : Point(),
64 : "If supplied, only faces with normal equal to this, up to "
65 : "normal_tol, will be added to the sidesets specified");
66 578185 : params.addRangeCheckedParam<Real>("normal_tol",
67 231274 : 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 :
74 : // Sideset restriction param group
75 346911 : params.addParamNamesToGroup(
76 : "included_boundaries excluded_boundaries included_subdomains included_neighbors "
77 : "include_only_external_sides normal normal_tol",
78 : "Sideset restrictions");
79 :
80 115637 : return params;
81 0 : }
82 :
83 5762 : SideSetsGeneratorBase::SideSetsGeneratorBase(const InputParameters & parameters)
84 : : MeshGenerator(parameters),
85 5762 : _input(getMesh("input")),
86 5762 : _boundary_names(std::vector<BoundaryName>()),
87 11524 : _fixed_normal(getParam<bool>("fixed_normal")),
88 17286 : _replace(getParam<bool>("replace")),
89 11524 : _check_included_boundaries(isParamValid("included_boundaries")),
90 11524 : _check_excluded_boundaries(isParamValid("excluded_boundaries")),
91 11524 : _check_subdomains(isParamValid("included_subdomains")),
92 11524 : _check_neighbor_subdomains(isParamValid("included_neighbors")),
93 5762 : _included_boundary_ids(std::vector<boundary_id_type>()),
94 5762 : _excluded_boundary_ids(std::vector<boundary_id_type>()),
95 5762 : _included_subdomain_ids(std::vector<subdomain_id_type>()),
96 5762 : _included_neighbor_subdomain_ids(std::vector<subdomain_id_type>()),
97 11524 : _include_only_external_sides(getParam<bool>("include_only_external_sides")),
98 11524 : _using_normal(isParamSetByUser("normal")),
99 15292 : _normal(_using_normal ? Point(getParam<Point>("normal") / getParam<Point>("normal").norm())
100 14774 : : getParam<Point>("normal")),
101 23048 : _normal_tol(getParam<Real>("normal_tol"))
102 : {
103 17286 : if (isParamValid("new_boundary"))
104 15060 : _boundary_names = getParam<std::vector<BoundaryName>>("new_boundary");
105 5762 : }
106 :
107 5494 : SideSetsGeneratorBase::~SideSetsGeneratorBase() {}
108 :
109 : void
110 5484 : SideSetsGeneratorBase::setup(MeshBase & mesh)
111 : {
112 : mooseAssert(_fe_face == nullptr, "FE Face has already been initialized");
113 :
114 : // To know the dimension of the mesh
115 5484 : if (!mesh.is_prepared())
116 4757 : mesh.prepare_for_use();
117 5484 : const auto dim = mesh.mesh_dimension();
118 :
119 : // Setup the FE Object so we can calculate normals
120 10968 : libMesh::FEType fe_type(Utility::string_to_enum<Order>("CONSTANT"),
121 5484 : Utility::string_to_enum<libMesh::FEFamily>("MONOMIAL"));
122 5484 : _fe_face = libMesh::FEBase::build(dim, fe_type);
123 5484 : _qface = std::make_unique<libMesh::QGauss>(dim - 1, FIRST);
124 5484 : _fe_face->attach_quadrature_rule(_qface.get());
125 : // Must always pre-request quantities you want to compute
126 5484 : _fe_face->get_normals();
127 :
128 : // Handle incompatible parameters
129 5484 : if (_include_only_external_sides && _check_neighbor_subdomains)
130 0 : paramError("include_only_external_sides", "External sides dont have neighbors");
131 :
132 5484 : if (_check_included_boundaries)
133 : {
134 702 : const auto & included_boundaries = getParam<std::vector<BoundaryName>>("included_boundaries");
135 698 : for (const auto & boundary_name : _boundary_names)
136 351 : if (std::find(included_boundaries.begin(), included_boundaries.end(), boundary_name) !=
137 702 : included_boundaries.end())
138 8 : paramError(
139 : "new_boundary",
140 : "A boundary cannot be both the new boundary and be included in the list of included "
141 : "boundaries. If you are trying to restrict an existing boundary, you must use a "
142 : "different name for 'new_boundary', delete the old boundary, and then rename the "
143 : "new boundary to the old boundary.");
144 :
145 347 : _included_boundary_ids = MooseMeshUtils::getBoundaryIDs(mesh, included_boundaries, false);
146 :
147 : // Check that the included boundary ids/names exist in the mesh
148 872 : for (const auto i : index_range(_included_boundary_ids))
149 529 : if (_included_boundary_ids[i] == Moose::INVALID_BOUNDARY_ID)
150 8 : paramError("included_boundaries",
151 : "The boundary '",
152 4 : included_boundaries[i],
153 : "' was not found within the mesh");
154 : }
155 :
156 5476 : if (_check_excluded_boundaries)
157 : {
158 42 : const auto & excluded_boundaries = getParam<std::vector<BoundaryName>>("excluded_boundaries");
159 38 : for (const auto & boundary_name : _boundary_names)
160 21 : if (std::find(excluded_boundaries.begin(), excluded_boundaries.end(), boundary_name) !=
161 42 : excluded_boundaries.end())
162 8 : paramError(
163 : "new_boundary",
164 : "A boundary cannot be both the new boundary and be excluded in the list of excluded "
165 : "boundaries.");
166 17 : _excluded_boundary_ids = MooseMeshUtils::getBoundaryIDs(mesh, excluded_boundaries, false);
167 :
168 : // Check that the excluded boundary ids/names exist in the mesh
169 30 : for (const auto i : index_range(_excluded_boundary_ids))
170 17 : if (_excluded_boundary_ids[i] == Moose::INVALID_BOUNDARY_ID)
171 8 : paramError("excluded_boundaries",
172 : "The boundary '",
173 4 : excluded_boundaries[i],
174 : "' was not found within the mesh");
175 :
176 13 : if (_check_included_boundaries)
177 : {
178 : // Check that included and excluded boundary lists do not overlap
179 4 : for (const auto & boundary_id : _included_boundary_ids)
180 4 : if (std::find(_excluded_boundary_ids.begin(), _excluded_boundary_ids.end(), boundary_id) !=
181 8 : _excluded_boundary_ids.end())
182 8 : paramError("excluded_boundaries",
183 : "'included_boundaries' and 'excluded_boundaries' lists should not overlap");
184 : }
185 : }
186 :
187 : // Get the boundary ids from the names
188 16392 : if (parameters().isParamValid("included_subdomains"))
189 : {
190 : // check that the subdomains exist in the mesh
191 9656 : const auto subdomains = getParam<std::vector<SubdomainName>>("included_subdomains");
192 9721 : for (const auto & name : subdomains)
193 4905 : if (!MooseMeshUtils::hasSubdomainName(mesh, name))
194 24 : paramError("included_subdomains", "The block '", name, "' was not found in the mesh");
195 :
196 4816 : _included_subdomain_ids = MooseMeshUtils::getSubdomainIDs(mesh, subdomains);
197 4816 : }
198 :
199 16356 : if (parameters().isParamValid("included_neighbors"))
200 : {
201 : // check that the subdomains exist in the mesh
202 7072 : const auto subdomains = getParam<std::vector<SubdomainName>>("included_neighbors");
203 7075 : for (const auto & name : subdomains)
204 3547 : if (!MooseMeshUtils::hasSubdomainName(mesh, name))
205 16 : paramError("included_neighbors", "The block '", name, "' was not found in the mesh");
206 :
207 3528 : _included_neighbor_subdomain_ids = MooseMeshUtils::getSubdomainIDs(mesh, subdomains);
208 3528 : }
209 :
210 : // We will want to Change the below code when we have more fine-grained control over advertising
211 : // what we need and how we satisfy those needs. For now we know we need to have neighbors per
212 : // #15823...and we do have an explicit `find_neighbors` call...but we don't have a
213 : // `neighbors_found` API and it seems off to do:
214 : //
215 : // if (!mesh.is_prepared())
216 : // mesh.find_neighbors()
217 5444 : }
218 :
219 : void
220 477 : SideSetsGeneratorBase::finalize()
221 : {
222 477 : _qface.reset();
223 477 : _fe_face.reset();
224 477 : }
225 :
226 : void
227 2224310 : SideSetsGeneratorBase::flood(const Elem * elem,
228 : const Point & normal,
229 : const boundary_id_type & side_id,
230 : MeshBase & mesh)
231 : {
232 4030773 : if (elem == nullptr || elem == remote_elem ||
233 4030773 : (_visited[side_id].find(elem) != _visited[side_id].end()))
234 1466080 : return;
235 :
236 : // Skip if element is not in specified subdomains
237 758230 : if (_check_subdomains && !elementSubdomainIdInList(elem, _included_subdomain_ids))
238 0 : return;
239 :
240 758230 : _visited[side_id].insert(elem);
241 :
242 : // Request to compute normal vectors
243 758230 : const std::vector<Point> & face_normals = _fe_face->get_normals();
244 :
245 4569020 : for (const auto side : make_range(elem->n_sides()))
246 : {
247 :
248 3810790 : _fe_face->reinit(elem, side);
249 : // We'll just use the normal of the first qp
250 3810790 : const Point face_normal = face_normals[0];
251 :
252 3810790 : if (!elemSideSatisfiesRequirements(elem, side, mesh, normal, face_normal))
253 3462179 : continue;
254 :
255 348611 : if (_replace)
256 1440 : mesh.get_boundary_info().remove_side(elem, side);
257 :
258 348611 : mesh.get_boundary_info().add_side(elem, side, side_id);
259 2105197 : for (const auto neighbor : make_range(elem->n_sides()))
260 : {
261 : // Flood to the neighboring elements using the current matching side normal from this
262 : // element.
263 : // This will allow us to tolerate small changes in the normals so we can "paint" around a
264 : // curve.
265 1756586 : flood(elem->neighbor_ptr(neighbor), _fixed_normal ? normal : face_normal, side_id, mesh);
266 : }
267 : }
268 : }
269 :
270 : bool
271 2355054 : SideSetsGeneratorBase::normalsWithinTol(const Point & normal_1,
272 : const Point & normal_2,
273 : const Real & tol) const
274 : {
275 2355054 : return (1.0 - normal_1 * normal_2) <= tol;
276 : }
277 :
278 : bool
279 1032182 : SideSetsGeneratorBase::elementSubdomainIdInList(
280 : const Elem * const elem, const std::vector<subdomain_id_type> & subdomain_id_list) const
281 : {
282 1032182 : subdomain_id_type curr_subdomain = elem->subdomain_id();
283 1032182 : return std::find(subdomain_id_list.begin(), subdomain_id_list.end(), curr_subdomain) !=
284 2064364 : subdomain_id_list.end();
285 : }
286 :
287 : bool
288 142218 : SideSetsGeneratorBase::elementSideInIncludedBoundaries(const Elem * const elem,
289 : const unsigned int side,
290 : const MeshBase & mesh) const
291 : {
292 282330 : for (const auto & bid : _included_boundary_ids)
293 146957 : if (mesh.get_boundary_info().has_boundary_id(elem, side, bid))
294 6845 : return true;
295 135373 : return false;
296 : }
297 :
298 : bool
299 1944 : SideSetsGeneratorBase::elementSideInExcludedBoundaries(const Elem * const elem,
300 : const unsigned int side,
301 : const MeshBase & mesh) const
302 : {
303 3870 : for (const auto bid : _excluded_boundary_ids)
304 1944 : if (mesh.get_boundary_info().has_boundary_id(elem, side, bid))
305 18 : return true;
306 1926 : return false;
307 : }
308 :
309 : bool
310 4646696 : SideSetsGeneratorBase::elemSideSatisfiesRequirements(const Elem * const elem,
311 : const unsigned int side,
312 : const MeshBase & mesh,
313 : const Point & desired_normal,
314 : const Point & face_normal)
315 : {
316 : // Skip if side has neighbor and we only want external sides
317 4646696 : if ((elem->neighbor_ptr(side) && _include_only_external_sides))
318 3456641 : return false;
319 :
320 : // Skip if side is not part of included boundaries
321 1190055 : if (_check_included_boundaries && !elementSideInIncludedBoundaries(elem, side, mesh))
322 135373 : return false;
323 : // Skip if side is part of excluded boundaries
324 1054682 : if (_check_excluded_boundaries && elementSideInExcludedBoundaries(elem, side, mesh))
325 18 : return false;
326 :
327 : // Skip if element does not have neighbor in specified subdomains
328 1054664 : if (_check_neighbor_subdomains)
329 : {
330 622420 : const Elem * const neighbor = elem->neighbor_ptr(side);
331 : // if the neighbor does not exist, then skip this face; we only add sidesets
332 : // between existing elems if _check_neighbor_subdomains is true
333 622420 : if (!(neighbor && elementSubdomainIdInList(neighbor, _included_neighbor_subdomain_ids)))
334 583913 : return false;
335 : }
336 :
337 470751 : if (_using_normal && !normalsWithinTol(desired_normal, face_normal, _normal_tol))
338 40654 : return false;
339 :
340 430097 : return true;
341 : }
|
