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 : #include "MeshTraversingUtils.h"
16 :
17 : #include "libmesh/mesh_generation.h"
18 : #include "libmesh/mesh.h"
19 : #include "libmesh/string_to_enum.h"
20 : #include "libmesh/quadrature_gauss.h"
21 : #include "libmesh/point_locator_base.h"
22 : #include "libmesh/elem.h"
23 : #include "libmesh/remote_elem.h"
24 :
25 : InputParameters
26 34649 : SideSetsGeneratorBase::validParams()
27 : {
28 34649 : InputParameters params = MeshGenerator::validParams();
29 138596 : params.addRequiredParam<MeshGeneratorName>("input", "The mesh we want to modify");
30 138596 : params.addRequiredParam<std::vector<BoundaryName>>(
31 : "new_boundary", "The list of boundary names to create on the supplied subdomain");
32 103947 : params.addParam<bool>("fixed_normal",
33 69298 : false,
34 : "This Boolean determines whether we fix our normal "
35 : "or allow it to vary to \"paint\" around curves");
36 :
37 103947 : params.addParam<bool>("replace",
38 69298 : false,
39 : "If true, replace the old sidesets. If false, the current sidesets (if "
40 : "any) will be preserved.");
41 :
42 138596 : params.addParam<std::vector<BoundaryName>>(
43 : "included_boundaries",
44 : "A set of boundary names or ids whose sides will be included in the new sidesets. A side "
45 : "is only added if it also belongs to one of these boundaries.");
46 138596 : params.addParam<std::vector<BoundaryName>>(
47 : "excluded_boundaries",
48 : "A set of boundary names or ids whose sides will be excluded from the new sidesets. A side "
49 : "is only added if does not belong to any of these boundaries.");
50 138596 : params.addParam<std::vector<SubdomainName>>(
51 : "included_subdomains",
52 : "A set of subdomain names or ids whose sides will be included in the new sidesets. A side "
53 : "is only added if the subdomain id of the corresponding element is in this set.");
54 138596 : params.addParam<std::vector<SubdomainName>>("included_neighbors",
55 : "A set of neighboring subdomain names or ids. A face "
56 : "is only added if the subdomain id of the "
57 : "neighbor is in this set");
58 103947 : params.addParam<bool>(
59 : "include_only_external_sides",
60 69298 : false,
61 : "Whether to only include external sides when considering sides to add to the sideset");
62 :
63 103947 : params.addParam<Point>("normal",
64 69298 : Point(),
65 : "If supplied, only faces with normal equal to this, up to "
66 : "normal_tol, will be added to the sidesets specified");
67 173245 : params.addRangeCheckedParam<Real>("normal_tol",
68 69298 : 0.1,
69 : "normal_tol>=0 & normal_tol<=2",
70 : "If normal is supplied then faces are "
71 : "only added if face_normal.normal_hat >= "
72 : "1 - normal_tol, where normal_hat = "
73 : "normal/|normal|");
74 :
75 : // Sideset restriction param group
76 103947 : params.addParamNamesToGroup(
77 : "included_boundaries excluded_boundaries included_subdomains included_neighbors "
78 : "include_only_external_sides normal normal_tol",
79 : "Sideset restrictions");
80 :
81 34649 : return params;
82 0 : }
83 :
84 5064 : SideSetsGeneratorBase::SideSetsGeneratorBase(const InputParameters & parameters)
85 : : MeshGenerator(parameters),
86 5064 : _input(getMesh("input")),
87 5064 : _boundary_names(std::vector<BoundaryName>()),
88 15192 : _fixed_normal(getParam<bool>("fixed_normal")),
89 10128 : _replace(getParam<bool>("replace")),
90 10128 : _check_included_boundaries(isParamValid("included_boundaries")),
91 10128 : _check_excluded_boundaries(isParamValid("excluded_boundaries")),
92 10128 : _check_subdomains(isParamValid("included_subdomains")),
93 10128 : _check_neighbor_subdomains(isParamValid("included_neighbors")),
94 5064 : _included_boundary_ids(std::vector<boundary_id_type>()),
95 5064 : _excluded_boundary_ids(std::vector<boundary_id_type>()),
96 5064 : _included_subdomain_ids(std::vector<subdomain_id_type>()),
97 5064 : _included_neighbor_subdomain_ids(std::vector<subdomain_id_type>()),
98 10128 : _include_only_external_sides(getParam<bool>("include_only_external_sides")),
99 10128 : _using_normal(isParamSetByUser("normal")),
100 13143 : _normal(_using_normal ? Point(getParam<Point>("normal") / getParam<Point>("normal").norm())
101 13182 : : getParam<Point>("normal")),
102 20256 : _normal_tol(getParam<Real>("normal_tol"))
103 : {
104 15192 : if (isParamValid("new_boundary"))
105 13209 : _boundary_names = getParam<std::vector<BoundaryName>>("new_boundary");
106 5064 : }
107 :
108 4863 : SideSetsGeneratorBase::~SideSetsGeneratorBase() {}
109 :
110 : void
111 4790 : SideSetsGeneratorBase::setup(MeshBase & mesh)
112 : {
113 : mooseAssert(_fe_face == nullptr, "FE Face has already been initialized");
114 :
115 : // To know the dimension of the mesh
116 4790 : if (!mesh.is_prepared())
117 4129 : mesh.prepare_for_use();
118 4790 : const auto dim = mesh.mesh_dimension();
119 :
120 : // Setup the FE Object so we can calculate normals
121 9580 : libMesh::FEType fe_type(Utility::string_to_enum<Order>("CONSTANT"),
122 4790 : Utility::string_to_enum<libMesh::FEFamily>("MONOMIAL"));
123 4790 : _fe_face = libMesh::FEBase::build(dim, fe_type);
124 4790 : _qface = std::make_unique<libMesh::QGauss>(dim - 1, FIRST);
125 4790 : _fe_face->attach_quadrature_rule(_qface.get());
126 : // Must always pre-request quantities you want to compute
127 4790 : _fe_face->get_normals();
128 :
129 : // Handle incompatible parameters
130 4790 : if (_include_only_external_sides && _check_neighbor_subdomains)
131 0 : paramError("include_only_external_sides", "External sides dont have neighbors");
132 :
133 4790 : if (_check_included_boundaries)
134 : {
135 640 : const auto & included_boundaries = getParam<std::vector<BoundaryName>>("included_boundaries");
136 637 : for (const auto & boundary_name : _boundary_names)
137 320 : if (std::find(included_boundaries.begin(), included_boundaries.end(), boundary_name) !=
138 640 : included_boundaries.end())
139 6 : paramError(
140 : "new_boundary",
141 : "A boundary cannot be both the new boundary and be included in the list of included "
142 : "boundaries. If you are trying to restrict an existing boundary, you must use a "
143 : "different name for 'new_boundary', delete the old boundary, and then rename the "
144 : "new boundary to the old boundary.");
145 :
146 317 : _included_boundary_ids = MooseMeshUtils::getBoundaryIDs(mesh, included_boundaries, false);
147 :
148 : // Check that the included boundary ids/names exist in the mesh
149 791 : for (const auto i : index_range(_included_boundary_ids))
150 477 : if (_included_boundary_ids[i] == Moose::INVALID_BOUNDARY_ID)
151 6 : paramError("included_boundaries",
152 : "The boundary '",
153 3 : included_boundaries[i],
154 : "' was not found within the mesh");
155 : }
156 :
157 4784 : if (_check_excluded_boundaries)
158 : {
159 34 : const auto & excluded_boundaries = getParam<std::vector<BoundaryName>>("excluded_boundaries");
160 31 : for (const auto & boundary_name : _boundary_names)
161 17 : if (std::find(excluded_boundaries.begin(), excluded_boundaries.end(), boundary_name) !=
162 34 : excluded_boundaries.end())
163 6 : paramError(
164 : "new_boundary",
165 : "A boundary cannot be both the new boundary and be excluded in the list of excluded "
166 : "boundaries.");
167 14 : _excluded_boundary_ids = MooseMeshUtils::getBoundaryIDs(mesh, excluded_boundaries, false);
168 :
169 : // Check that the excluded boundary ids/names exist in the mesh
170 25 : for (const auto i : index_range(_excluded_boundary_ids))
171 14 : if (_excluded_boundary_ids[i] == Moose::INVALID_BOUNDARY_ID)
172 6 : paramError("excluded_boundaries",
173 : "The boundary '",
174 3 : excluded_boundaries[i],
175 : "' was not found within the mesh");
176 :
177 11 : if (_check_included_boundaries)
178 : {
179 : // Check that included and excluded boundary lists do not overlap
180 3 : for (const auto & boundary_id : _included_boundary_ids)
181 3 : if (std::find(_excluded_boundary_ids.begin(), _excluded_boundary_ids.end(), boundary_id) !=
182 6 : _excluded_boundary_ids.end())
183 6 : paramError("excluded_boundaries",
184 : "'included_boundaries' and 'excluded_boundaries' lists should not overlap");
185 : }
186 : }
187 :
188 : // Get the boundary ids from the names
189 14325 : if (parameters().isParamValid("included_subdomains"))
190 : {
191 : // check that the subdomains exist in the mesh
192 8404 : const auto subdomains = getParam<std::vector<SubdomainName>>("included_subdomains");
193 8464 : for (const auto & name : subdomains)
194 4271 : if (!MooseMeshUtils::hasSubdomainName(mesh, name))
195 18 : paramError("included_subdomains", "The block '", name, "' was not found in the mesh");
196 :
197 4193 : _included_subdomain_ids = MooseMeshUtils::getSubdomainIDs(mesh, subdomains);
198 4193 : }
199 :
200 14298 : if (parameters().isParamValid("included_neighbors"))
201 : {
202 : // check that the subdomains exist in the mesh
203 6336 : const auto subdomains = getParam<std::vector<SubdomainName>>("included_neighbors");
204 6340 : for (const auto & name : subdomains)
205 3178 : if (!MooseMeshUtils::hasSubdomainName(mesh, name))
206 12 : paramError("included_neighbors", "The block '", name, "' was not found in the mesh");
207 :
208 3162 : _included_neighbor_subdomain_ids = MooseMeshUtils::getSubdomainIDs(mesh, subdomains);
209 3162 : }
210 :
211 : // We will want to Change the below code when we have more fine-grained control over advertising
212 : // what we need and how we satisfy those needs. For now we know we need to have neighbors per
213 : // #15823...and we do have an explicit `find_neighbors` call...but we don't have a
214 : // `neighbors_found` API and it seems off to do:
215 : //
216 : // if (!mesh.is_prepared())
217 : // mesh.find_neighbors()
218 4760 : }
219 :
220 : void
221 420 : SideSetsGeneratorBase::finalize()
222 : {
223 420 : _qface.reset();
224 420 : _fe_face.reset();
225 420 : }
226 :
227 : void
228 1977202 : SideSetsGeneratorBase::flood(const Elem * elem,
229 : const Point & normal,
230 : const boundary_id_type & side_id,
231 : MeshBase & mesh)
232 : {
233 3582976 : if (elem == nullptr || elem == remote_elem ||
234 3582976 : (_visited[side_id].find(elem) != _visited[side_id].end()))
235 1303184 : return;
236 :
237 : // Skip if element is not in specified subdomains
238 674018 : if (_check_subdomains &&
239 0 : !MeshTraversingUtils::elementSubdomainIdInList(elem, _included_subdomain_ids))
240 0 : return;
241 :
242 674018 : _visited[side_id].insert(elem);
243 :
244 : // Request to compute normal vectors
245 674018 : const std::vector<Point> & face_normals = _fe_face->get_normals();
246 :
247 4061494 : for (const auto side : make_range(elem->n_sides()))
248 : {
249 :
250 3387476 : _fe_face->reinit(elem, side);
251 : // We'll just use the normal of the first qp
252 3387476 : const Point face_normal = face_normals[0];
253 :
254 3387476 : if (!elemSideSatisfiesRequirements(elem, side, mesh, normal, face_normal))
255 3077590 : continue;
256 :
257 309886 : if (_replace)
258 1280 : mesh.get_boundary_info().remove_side(elem, side);
259 :
260 309886 : mesh.get_boundary_info().add_side(elem, side, side_id);
261 1871322 : 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 1561436 : flood(elem->neighbor_ptr(neighbor), _fixed_normal ? normal : face_normal, side_id, mesh);
268 : }
269 : }
270 : }
271 :
272 : bool
273 121908 : SideSetsGeneratorBase::elementSideInIncludedBoundaries(const Elem * const elem,
274 : const unsigned int side,
275 : const MeshBase & mesh) const
276 : {
277 242054 : for (const auto & bid : _included_boundary_ids)
278 126086 : if (mesh.get_boundary_info().has_boundary_id(elem, side, bid))
279 5940 : return true;
280 115968 : return false;
281 : }
282 :
283 : bool
284 1728 : SideSetsGeneratorBase::elementSideInExcludedBoundaries(const Elem * const elem,
285 : const unsigned int side,
286 : const MeshBase & mesh) const
287 : {
288 3440 : for (const auto bid : _excluded_boundary_ids)
289 1728 : if (mesh.get_boundary_info().has_boundary_id(elem, side, bid))
290 16 : return true;
291 1712 : return false;
292 : }
293 :
294 : bool
295 4231483 : SideSetsGeneratorBase::elemSideSatisfiesRequirements(const Elem * const elem,
296 : const unsigned int side,
297 : const MeshBase & mesh,
298 : const Point & desired_normal,
299 : const Point & face_normal)
300 : {
301 : // Skip if side has neighbor and we only want external sides
302 4231483 : if ((elem->neighbor_ptr(side) && _include_only_external_sides))
303 3072642 : return false;
304 :
305 : // Skip if side is not part of included boundaries
306 1158841 : if (_check_included_boundaries && !elementSideInIncludedBoundaries(elem, side, mesh))
307 115968 : return false;
308 : // Skip if side is part of excluded boundaries
309 1042873 : if (_check_excluded_boundaries && elementSideInExcludedBoundaries(elem, side, mesh))
310 16 : return false;
311 :
312 : // Skip if element does not have neighbor in specified subdomains
313 1042857 : if (_check_neighbor_subdomains)
314 : {
315 656859 : const Elem * const neighbor = elem->neighbor_ptr(side);
316 : // if the neighbor does not exist, then skip this face; we only add sidesets
317 : // between existing elems if _check_neighbor_subdomains is true
318 1313510 : if (!(neighbor && MeshTraversingUtils::elementSubdomainIdInList(
319 656651 : neighbor, _included_neighbor_subdomain_ids)))
320 617775 : return false;
321 : }
322 :
323 484212 : if (_using_normal &&
324 59130 : !MeshTraversingUtils::normalsWithinTol(desired_normal, face_normal, _normal_tol))
325 37695 : return false;
326 :
327 387387 : return true;
328 : }
|
