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 "PatternedPolygonPeripheralModifierBase.h"
11 :
12 : #include "MooseMeshUtils.h"
13 : #include "FillBetweenPointVectorsTools.h"
14 : #include "KDTree.h"
15 :
16 : InputParameters
17 360 : PatternedPolygonPeripheralModifierBase::validParams()
18 : {
19 360 : InputParameters params = PolygonMeshGeneratorBase::validParams();
20 720 : params.addRequiredParam<MeshGeneratorName>(
21 : "input",
22 : "The input mesh to be modified. Note that this generator only works with "
23 : "PatternedHex/CartesianMeshGenerator and its derived classes such as "
24 : "HexIDPatternedMeshGenerator.");
25 720 : params.addRequiredParam<BoundaryName>("input_mesh_external_boundary",
26 : "The external boundary of the input mesh.");
27 720 : params.addRequiredParam<unsigned int>("new_num_sector",
28 : "Number of sectors of each side for the new mesh.");
29 720 : params.addParam<subdomain_id_type>(
30 : "transition_layer_id",
31 720 : TRANSITION_LAYER_DEFAULT,
32 : "Optional customized block id for the transition layer block.");
33 720 : params.addParam<SubdomainName>("transition_layer_name",
34 : "transition_layer",
35 : "Optional customized block name for the transition layer block.");
36 1080 : params.addRangeCheckedParam<unsigned int>(
37 720 : "num_layers", 1, "num_layers>0", "Layers of elements for transition.");
38 720 : params.addParam<std::vector<std::string>>(
39 : "extra_id_names_to_modify",
40 : "Names of the element extra ids in the peripheral region that should be modified");
41 360 : params.addParam<std::vector<dof_id_type>>(
42 : "new_extra_id_values_to_assign",
43 360 : std::vector<dof_id_type>(),
44 : "Values of the modified extra ids in the peripheral region.");
45 360 : params.addClassDescription("PatternedPolygonPeripheralModifierBase is the base class for "
46 : "PatternedCartPeripheralModifier and PatternedHexPeripheralModifier.");
47 :
48 360 : return params;
49 0 : }
50 :
51 176 : PatternedPolygonPeripheralModifierBase::PatternedPolygonPeripheralModifierBase(
52 176 : const InputParameters & parameters)
53 : : PolygonMeshGeneratorBase(parameters),
54 0 : _input_name(getParam<MeshGeneratorName>("input")),
55 176 : _input_mesh_external_boundary(getParam<BoundaryName>("input_mesh_external_boundary")),
56 352 : _new_num_sector(getParam<unsigned int>("new_num_sector")),
57 352 : _transition_layer_id(getParam<subdomain_id_type>("transition_layer_id")),
58 176 : _transition_layer_name(getParam<SubdomainName>("transition_layer_name")),
59 352 : _num_layers(getParam<unsigned int>("num_layers")),
60 374 : _extra_id_names_to_modify(isParamValid("extra_id_names_to_modify")
61 176 : ? getParam<std::vector<std::string>>("extra_id_names_to_modify")
62 : : std::vector<std::string>()),
63 352 : _new_extra_id_values_to_assign(
64 : getParam<std::vector<dof_id_type>>("new_extra_id_values_to_assign")),
65 : // Use CartesianConcentricCircleAdaptiveBoundaryMeshGenerator for cartesian control drum meshing
66 : // Use HexagonConcentricCircleAdaptiveBoundaryMeshGenerator for hexagonal control drum meshing
67 352 : _mesh(getMeshByName(_input_name))
68 : {
69 176 : declareMeshProperty<Real>("pattern_pitch_meta", 0.0);
70 352 : declareMeshProperty<bool>("is_control_drum_meta", false);
71 176 : if (_extra_id_names_to_modify.size() != _new_extra_id_values_to_assign.size())
72 2 : paramError("new_extra_id_values_to_assign",
73 : "This parameter must have the same size as extra_id_names_to_modify.");
74 174 : declareMeshProperty<bool>("peripheral_modifier_compatible", false);
75 174 : }
76 :
77 : std::unique_ptr<MeshBase>
78 166 : PatternedPolygonPeripheralModifierBase::generate()
79 : {
80 : // Transmit the Mesh Metadata
81 332 : auto pattern_pitch_meta = setMeshProperty(
82 166 : "pattern_pitch_meta", getMeshProperty<Real>("pattern_pitch_meta", _input_name));
83 : // Load the input mesh as ReplicatedMesh
84 166 : auto input_mesh = dynamic_pointer_cast<ReplicatedMesh>(_mesh);
85 : // Load boundary information
86 166 : _input_mesh_external_bid =
87 166 : MooseMeshUtils::getBoundaryID(_input_mesh_external_boundary, *input_mesh);
88 166 : if (_input_mesh_external_bid == Moose::INVALID_BOUNDARY_ID)
89 0 : paramError("input_mesh_external_boundary",
90 : "External boundary does not exist in the input mesh");
91 : std::vector<std::tuple<dof_id_type, unsigned short int, boundary_id_type>> side_list =
92 166 : input_mesh->get_boundary_info().build_side_list();
93 166 : input_mesh->get_boundary_info().build_node_list_from_side_list();
94 : std::vector<std::tuple<dof_id_type, boundary_id_type>> node_list =
95 166 : input_mesh->get_boundary_info().build_node_list();
96 : // Find neighbors
97 166 : if (!input_mesh->is_prepared())
98 166 : input_mesh->find_neighbors();
99 : // Load all the extra integer indexing information from the input mesh
100 : const unsigned int n_extra_integer_input = input_mesh->n_elem_integers();
101 : std::vector<std::string> extra_integer_names;
102 186 : for (unsigned int i = 0; i < n_extra_integer_input; i++)
103 20 : extra_integer_names.push_back(input_mesh->get_elem_integer_name(i));
104 :
105 : // Check elem extra integers names and convert them into indices
106 166 : std::vector<bool> extra_id_modify_flags = std::vector<bool>(n_extra_integer_input, false);
107 166 : std::vector<dof_id_type> extra_modify_values = std::vector<dof_id_type>(n_extra_integer_input, 0);
108 175 : for (unsigned int i = 0; i < _extra_id_names_to_modify.size(); i++)
109 : {
110 11 : if (!input_mesh->has_elem_integer(_extra_id_names_to_modify[i]))
111 2 : paramError(
112 : "extra_id_names_to_modify",
113 : "The parameter contains an extra element integer that does not exist in the input mesh.");
114 : else
115 : {
116 9 : extra_id_modify_flags[input_mesh->get_elem_integer_index(_extra_id_names_to_modify[i])] =
117 : true;
118 9 : extra_modify_values[input_mesh->get_elem_integer_index(_extra_id_names_to_modify[i])] =
119 : _new_extra_id_values_to_assign[i];
120 : }
121 : }
122 :
123 : std::vector<dof_id_type> bid_elem_list;
124 : std::vector<dof_id_type> new_bid_elem_list;
125 : std::vector<dof_id_type> ext_node_list;
126 : // We check the type of the QUAD elements on the boundary
127 : unsigned short quad_type = 0;
128 39956 : for (unsigned int i = 0; i < side_list.size(); i++)
129 : {
130 39792 : if (std::get<2>(side_list[i]) == _input_mesh_external_bid)
131 : {
132 : // List of elements on the boundary
133 6744 : bid_elem_list.push_back(std::get<0>(side_list[i]));
134 : // Check if the element is QUAD4 or QUAD8/9
135 6744 : const auto elem_type = input_mesh->elem_ptr(std::get<0>(side_list[i]))->type();
136 6744 : if (elem_type != QUAD4 && elem_type != QUAD8 && elem_type != QUAD9)
137 0 : paramError("input",
138 : "The input mesh has non-QUAD4/QUAD8/QUAD9 elements in its peripheral area, "
139 : "which is not supported.");
140 : // Generally, we need to check if the sideset involves mixed order elements
141 : // But this is not possible for a mesh generated by PatternedHex/CartesianMeshGenerator
142 : // So we only need to check the first element
143 6744 : else if (quad_type == 0)
144 164 : quad_type = elem_type == QUAD4 ? 1 : (elem_type == QUAD8 ? 8 : 9);
145 :
146 : // Note this object only works with quad elements
147 : // (side_id + 2)%4 gives the opposite side's neighboring element
148 : // Thus, list of elements on the new boundary is made
149 : new_bid_elem_list.push_back(input_mesh->elem_ptr(std::get<0>(side_list[i]))
150 6744 : ->neighbor_ptr((std::get<1>(side_list[i]) + 2) % 4)
151 6744 : ->id());
152 : }
153 : }
154 : // This eliminates the duplicate elements, which happen at the vertices.
155 164 : auto bid_elem_list_it = std::unique(bid_elem_list.begin(), bid_elem_list.end());
156 164 : auto new_bid_elem_list_it = std::unique(new_bid_elem_list.begin(), new_bid_elem_list.end());
157 164 : bid_elem_list.resize(std::distance(bid_elem_list.begin(), bid_elem_list_it));
158 164 : new_bid_elem_list.resize(std::distance(new_bid_elem_list.begin(), new_bid_elem_list_it));
159 : // Create a data structure to record centroid positions and extra element IDS of the
160 : // elements to be deleted
161 : std::vector<std::pair<Point, std::vector<dof_id_type>>> del_elem_extra_ids;
162 : // Delete the outer layer QUAD elements so that they can be replaced by TRI elements later.
163 6908 : for (unsigned int i = 0; i < bid_elem_list.size(); i++)
164 : {
165 6744 : const auto & tmp_elem_ptr = input_mesh->elem_ptr(bid_elem_list[i]);
166 : // Retain all the extra integer information first if any of them need to be retained
167 6744 : if (n_extra_integer_input > _extra_id_names_to_modify.size())
168 : {
169 : del_elem_extra_ids.push_back(
170 648 : std::make_pair(tmp_elem_ptr->true_centroid(), std::vector<dof_id_type>()));
171 648 : for (unsigned int j = 0; j < n_extra_integer_input; j++)
172 324 : del_elem_extra_ids[i].second.push_back(
173 324 : extra_id_modify_flags[j] ? extra_modify_values[j] : tmp_elem_ptr->get_extra_integer(j));
174 : }
175 6744 : input_mesh->delete_elem(tmp_elem_ptr);
176 : }
177 : // if all the extra id will need to be re-assigned, just a dummy reference vector is needed.
178 164 : if (n_extra_integer_input == _extra_id_names_to_modify.size())
179 : del_elem_extra_ids.push_back(
180 310 : std::make_pair(Point(0.0, 0.0, 0.0), _new_extra_id_values_to_assign));
181 : // As some elements are deleted, update the neighbor list
182 164 : input_mesh->find_neighbors();
183 : // Identify the new external boundary
184 : BoundaryInfo & boundary_info = input_mesh->get_boundary_info();
185 6908 : for (unsigned int i = 0; i < new_bid_elem_list.size(); i++)
186 : {
187 : // Assign default external Sideset ID to the new boundary
188 33720 : for (unsigned int j = 0; j < input_mesh->elem_ptr(new_bid_elem_list[i])->n_sides(); j++)
189 : {
190 26976 : if (input_mesh->elem_ptr(new_bid_elem_list[i])->neighbor_ptr(j) == nullptr)
191 6744 : boundary_info.add_side(new_bid_elem_list[i], j, OUTER_SIDESET_ID);
192 : }
193 : }
194 :
195 : // Refresh the mesh as it lost some elements
196 164 : input_mesh->contract();
197 164 : input_mesh->prepare_for_use();
198 :
199 : // Clone the original mesh without outer layer for information extraction
200 164 : auto input_mesh_origin = dynamic_pointer_cast<ReplicatedMesh>(input_mesh->clone());
201 :
202 : // Make four (cartesian) or six (hexagonal) meshes of new outer layer for four or six sides.
203 1000 : for (unsigned i_side = 0; i_side < _num_sides; i_side++)
204 : {
205 : // Use azimuthalAnglesCollector to assign boundary Points to a reference vector
206 : std::vector<Point> inner_pts;
207 : std::vector<Real> new_bdry_azi = azimuthalAnglesCollector(*input_mesh_origin,
208 : inner_pts,
209 : -180.0 / (Real)_num_sides,
210 836 : 180.0 / (Real)_num_sides,
211 836 : ANGLE_DEGREE);
212 836 : MeshTools::Modification::rotate(*input_mesh_origin, -360.0 / (Real)_num_sides, 0, 0);
213 :
214 : // Setting up of the outer boundary is done by defining the two ends
215 : // And then use libMesh interpolation tool to make the straight line
216 : std::vector<Point> outer_pts;
217 : const Point outer_end_1 = Point(pattern_pitch_meta / 2.0,
218 836 : -pattern_pitch_meta / 2.0 * std::tan(M_PI / (Real)_num_sides),
219 836 : 0.0);
220 : const Point outer_end_2 = Point(pattern_pitch_meta / 2.0,
221 : pattern_pitch_meta / 2.0 * std::tan(M_PI / (Real)_num_sides),
222 836 : 0.0);
223 :
224 836 : const std::vector<Real> input_arg{0.0, 1.0};
225 836 : const std::vector<Real> input_x{outer_end_1(0), outer_end_2(0)};
226 836 : const std::vector<Real> input_y{outer_end_1(1), outer_end_2(1)};
227 : std::unique_ptr<LinearInterpolation> linear_outer_x =
228 836 : std::make_unique<LinearInterpolation>(input_arg, input_x);
229 : std::unique_ptr<LinearInterpolation> linear_outer_y =
230 836 : std::make_unique<LinearInterpolation>(input_arg, input_y);
231 :
232 : // Uniformly spaced nodes on the outer boundary to facilitate stitching
233 9492 : for (unsigned int i = 0; i < _new_num_sector + 1; i++)
234 : {
235 8656 : outer_pts.push_back(Point(linear_outer_x->sample((Real)i / (Real)_new_num_sector),
236 17312 : linear_outer_y->sample((Real)i / (Real)_new_num_sector),
237 : 0.0));
238 : }
239 : // Now we have the inner and outer boundary points as input for the transition layer
240 : // For linear elements, that's all we need;
241 : // For quadratic elements, the inner point vector contains both vertices and midpoints.
242 : // fillBetweenPointVectorsGenerator only supports linear elements
243 : // So we need to
244 : // (1) remove the midpoints from the inner point vector
245 : // (2) generate linear mesh using fillBetweenPointVectorsGenerator
246 : // (3) convert the linear mesh to quadratic mesh
247 : // Generally, in the input mesh, the midpoint might not be located in the middle of the two
248 : // vertices. (4) In that case, the midpoint needs to be moved to match the input mesh In the
249 : // case of PatternedHex/CartesianMeshGenerator, the midpoints of the inner boundary shoule be
250 : // the average of the vertices unless the non-circular regions are deformed. So this step still
251 : // needs to be done.
252 :
253 : // Use fillBetweenPointVectorsGenerator to generate the transition layer
254 :
255 : // Remove the midpoints from the inner point vector
256 : // But record the midpoints for later use
257 : // In the meantime, calculate the average of the two vertices
258 : std::vector<std::pair<Point, Point>> inner_midpts_pairs;
259 836 : if (quad_type > 4)
260 : {
261 90 : auto inner_pts_tmp(inner_pts);
262 90 : inner_pts.clear();
263 1260 : for (const auto i : index_range(inner_pts_tmp))
264 1170 : if (i % 2 == 0)
265 630 : inner_pts.push_back(inner_pts_tmp[i]);
266 : else
267 : {
268 : // We check if the midpoint is located in the middle of the two vertices
269 : // If not, we record the them for later use
270 540 : const Point avgpt = 0.5 * (inner_pts_tmp[i - 1] + inner_pts_tmp[i + 1]);
271 540 : if (!MooseUtils::absoluteFuzzyEqual((avgpt - inner_pts_tmp[i]).norm(), 0.0))
272 108 : inner_midpts_pairs.push_back(std::make_pair(avgpt, inner_pts_tmp[i]));
273 : }
274 90 : inner_pts_tmp.clear();
275 90 : }
276 :
277 836 : auto mesh = buildReplicatedMesh(2);
278 836 : FillBetweenPointVectorsTools::fillBetweenPointVectorsGenerator(*mesh,
279 : inner_pts,
280 : outer_pts,
281 836 : _num_layers,
282 836 : _transition_layer_id,
283 : OUTER_SIDESET_ID,
284 836 : _type,
285 836 : _name);
286 :
287 : // Convert the linear mesh to quadratic mesh
288 836 : if (quad_type == 8)
289 54 : mesh->all_second_order();
290 782 : else if (quad_type == 9)
291 36 : mesh->all_complete_order();
292 :
293 : // Move the midpoints to match the input mesh if applicable
294 836 : if (inner_midpts_pairs.size())
295 : {
296 54 : mesh->get_boundary_info().build_node_list_from_side_list();
297 54 : auto tl_node_list = mesh->get_boundary_info().build_node_list();
298 : std::vector<dof_id_type> bdry_node_ids;
299 : std::vector<Point> bdry_node_pts;
300 2214 : for (const auto & tl_node_tuple : tl_node_list)
301 : {
302 2160 : if (std::get<1>(tl_node_tuple) == OUTER_SIDESET_ID)
303 : {
304 2160 : bdry_node_ids.push_back(std::get<0>(tl_node_tuple));
305 2160 : bdry_node_pts.push_back(*(mesh->node_ptr(std::get<0>(tl_node_tuple))));
306 : }
307 : }
308 54 : KDTree ref_kd_tree(bdry_node_pts, 4);
309 162 : for (const auto & midpt_pair : inner_midpts_pairs)
310 : {
311 : std::vector<std::size_t> nn_id;
312 108 : ref_kd_tree.neighborSearch(midpt_pair.first, 1, nn_id);
313 108 : *(mesh->node_ptr(bdry_node_ids[nn_id.front()])) = midpt_pair.second;
314 108 : }
315 54 : }
316 :
317 836 : MeshTools::Modification::rotate(*mesh, (Real)i_side * 360.0 / (Real)_num_sides, 0, 0);
318 : // Assign extra element id based on the nearest deleted element
319 836 : mesh->add_elem_integers(extra_integer_names, true);
320 836 : transferExtraElemIntegers(*mesh, del_elem_extra_ids);
321 836 : mesh->prepare_for_use();
322 836 : input_mesh->stitch_meshes(*mesh, OUTER_SIDESET_ID, OUTER_SIDESET_ID, TOLERANCE, true);
323 836 : mesh->clear();
324 836 : }
325 164 : input_mesh->subdomain_name(_transition_layer_id) = _transition_layer_name;
326 164 : return input_mesh;
327 328 : }
328 :
329 : void
330 836 : PatternedPolygonPeripheralModifierBase::transferExtraElemIntegers(
331 : ReplicatedMesh & mesh,
332 : const std::vector<std::pair<Point, std::vector<dof_id_type>>> ref_extra_ids)
333 : {
334 : // Build points vector for k-d tree constructor
335 : std::vector<Point> ref_pts;
336 3562 : for (auto & pt_extra_id : ref_extra_ids)
337 2726 : ref_pts.push_back(pt_extra_id.first);
338 : // K-d tree construction
339 836 : KDTree ref_kd_tree(ref_pts, 4);
340 : // Use the k-d tree for nearest neighbor searching
341 58604 : for (const auto & elem : as_range(mesh.active_elements_begin(), mesh.active_elements_end()))
342 : {
343 : std::vector<std::size_t> nn_id;
344 28048 : ref_kd_tree.neighborSearch(elem->true_centroid(), 1, nn_id);
345 31504 : for (unsigned int i = 0; i < ref_extra_ids[nn_id.front()].second.size(); i++)
346 3456 : elem->set_extra_integer(i, ref_extra_ids[nn_id.front()].second[i]);
347 28884 : }
348 836 : }
|
