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 "FillBetweenPointVectorsTools.h"
11 : #include "MooseMeshUtils.h"
12 : #include "MooseMesh.h"
13 : #include "MeshGenerator.h"
14 : #include "MooseError.h"
15 :
16 : // libMesh includes
17 : #include "libmesh/int_range.h"
18 : #include "libmesh/mesh_base.h"
19 : #include "libmesh/mesh_generation.h"
20 : #include "libmesh/mesh_serializer.h"
21 : #include "libmesh/point.h"
22 : #include "libmesh/elem.h"
23 : #include "libmesh/node.h"
24 : #include "libmesh/face_tri3.h"
25 : #include "libmesh/face_quad4.h"
26 :
27 : using namespace libMesh;
28 :
29 : namespace FillBetweenPointVectorsTools
30 : {
31 : void
32 143 : fillBetweenPointVectorsGenerator(MeshBase & mesh, // an empty mesh is expected
33 : const std::vector<Point> & boundary_points_vec_1,
34 : const std::vector<Point> & boundary_points_vec_2,
35 : const unsigned int num_layers,
36 : const subdomain_id_type transition_layer_id,
37 : const boundary_id_type input_boundary_1_id,
38 : const boundary_id_type input_boundary_2_id,
39 : const boundary_id_type begin_side_boundary_id,
40 : const boundary_id_type end_side_boundary_id,
41 : const std::string type,
42 : const std::string name,
43 : const bool quad_elem,
44 : const Real bias_parameter,
45 : const Real sigma)
46 : {
47 286 : if (!MooseMeshUtils::isCoPlanar(
48 429 : boundary_points_vec_1, Point(0.0, 0.0, 1.0), Point(0.0, 0.0, 0.0)) ||
49 286 : !MooseMeshUtils::isCoPlanar(
50 : boundary_points_vec_2, Point(0.0, 0.0, 1.0), Point(0.0, 0.0, 0.0)))
51 4 : mooseError("In ",
52 : type,
53 : " ",
54 : name,
55 : ", the input vectors of points for "
56 : "FillBetweenPointVectorsTools::fillBetweenPointVectorsGenerator "
57 : "must be in XY plane.");
58 :
59 139 : const unsigned int vec_1_node_num = boundary_points_vec_1.size();
60 139 : const unsigned int vec_2_node_num = boundary_points_vec_2.size();
61 :
62 139 : if (vec_1_node_num < 2 || vec_2_node_num < 2)
63 4 : mooseError("In ",
64 : type,
65 : " ",
66 : name,
67 : ", the two input vectors of points for "
68 : "FillBetweenPointVectorsTools::fillBetweenPointVectorsGenerator "
69 : "must respectively contain at least two elements.");
70 :
71 135 : if (quad_elem && boundary_points_vec_1.size() != boundary_points_vec_2.size())
72 4 : mooseError("In ",
73 : type,
74 : " ",
75 : name,
76 : ", QUAD4 elements option can only be selected when the two input vectors of Points "
77 : "have the same length. In the current instance, the first vector has ",
78 4 : boundary_points_vec_1.size(),
79 : " points and the second ",
80 4 : boundary_points_vec_2.size(),
81 : " points");
82 :
83 131 : std::vector<Point> possibly_reoriented_boundary_points_vec_2;
84 131 : const std::vector<Point> * oriented_boundary_points_vec_2 = &boundary_points_vec_2;
85 :
86 131 : if (needFlip(boundary_points_vec_1, boundary_points_vec_2))
87 : {
88 0 : possibly_reoriented_boundary_points_vec_2.assign(boundary_points_vec_2.rbegin(),
89 0 : boundary_points_vec_2.rend());
90 0 : oriented_boundary_points_vec_2 = &possibly_reoriented_boundary_points_vec_2;
91 :
92 : // This isn't worth warning about. The way
93 : // MooseMeshUtils::makeOrderedNodeList works, we can end up
94 : // finding a flip necessary on one element numbering and
95 : // unnecessary on another.
96 : //
97 : // mooseWarning(
98 : // "In FillBetweenPointVectorsTools, one of the vector of Points must be flipped to ensure "
99 : // "correct transition layer shape.");
100 : }
101 :
102 131 : std::vector<Real> vec_1_index; // Unweighted index
103 131 : std::vector<Real> vec_2_index; // Unweighted index
104 :
105 131 : std::vector<Real> wt_1;
106 131 : std::vector<Real> index_1; // Weighted index
107 131 : std::vector<Real> wt_2;
108 131 : std::vector<Real> index_2; // Weighted index
109 :
110 : // create interpolations
111 131 : std::unique_ptr<LinearInterpolation> linear_vec_1_x;
112 131 : std::unique_ptr<LinearInterpolation> linear_vec_1_y;
113 131 : std::unique_ptr<SplineInterpolation> spline_vec_1_l;
114 131 : std::unique_ptr<LinearInterpolation> linear_vec_2_x;
115 131 : std::unique_ptr<LinearInterpolation> linear_vec_2_y;
116 131 : std::unique_ptr<SplineInterpolation> spline_vec_2_l;
117 :
118 131 : weightedInterpolator(vec_1_node_num,
119 : boundary_points_vec_1,
120 : vec_1_index,
121 : wt_1,
122 : index_1,
123 : sigma,
124 : linear_vec_1_x,
125 : linear_vec_1_y,
126 : spline_vec_1_l);
127 131 : weightedInterpolator(vec_2_node_num,
128 : *oriented_boundary_points_vec_2,
129 : vec_2_index,
130 : wt_2,
131 : index_2,
132 : sigma,
133 : linear_vec_2_x,
134 : linear_vec_2_y,
135 : spline_vec_2_l);
136 :
137 : // If the two input vectors have different sizes
138 : // The node numbers of the intermediate layers change linearly
139 131 : const Real increment = ((Real)vec_2_node_num - (Real)vec_1_node_num) / (Real)(num_layers);
140 : // Number of nodes in each sublayer
141 131 : std::vector<unsigned int> node_number_vec;
142 : // 2D vector of nodes
143 131 : std::vector<std::vector<Node *>> nodes(num_layers + 1);
144 : // Node counter
145 131 : unsigned int node_counter = 0;
146 :
147 720 : for (const auto i : make_range(num_layers + 1))
148 : {
149 : // calculate number of nodes in each sublayer
150 589 : node_number_vec.push_back(
151 589 : (unsigned int)(vec_1_node_num + (long)(increment * i + 0.5 - (increment < 0))));
152 : // Reserve memory for new nodes
153 589 : nodes[i] = std::vector<Node *>(node_number_vec[i]);
154 :
155 : // Calculate vectors of weighted surrogated index for side #1
156 589 : std::vector<Real> weighted_surrogate_index_1;
157 589 : std::vector<Real> unweighted_surrogate_index_1;
158 :
159 589 : surrogateGenerator(weighted_surrogate_index_1,
160 : unweighted_surrogate_index_1,
161 : node_number_vec,
162 : wt_1,
163 : vec_1_index,
164 : vec_1_node_num,
165 : i);
166 :
167 : // Calculate vectors of weighted surrogated index for side #2
168 589 : std::vector<Real> weighted_surrogate_index_2;
169 589 : std::vector<Real> unweighted_surrogate_index_2;
170 :
171 589 : surrogateGenerator(weighted_surrogate_index_2,
172 : unweighted_surrogate_index_2,
173 : node_number_vec,
174 : wt_2,
175 : vec_2_index,
176 : vec_2_node_num,
177 : i);
178 :
179 5312 : for (const auto j : make_range(node_number_vec[i]))
180 : {
181 : // Create surrogate Points on side #1 for Point #j on the sublayer
182 4723 : Point surrogate_pos_1 = Point(linear_vec_1_x->sample(weighted_surrogate_index_1[j]),
183 4723 : linear_vec_1_y->sample(weighted_surrogate_index_1[j]),
184 14169 : 0.0);
185 : // Create surrogate Points on side #2 for Point #j on the sublayer
186 4723 : Point surrogate_pos_2 = Point(linear_vec_2_x->sample(weighted_surrogate_index_2[j]),
187 4723 : linear_vec_2_y->sample(weighted_surrogate_index_2[j]),
188 14169 : 0.0);
189 : const Real l_ratio = bias_parameter <= 0.0
190 5871 : ? std::pow(spline_vec_2_l->sample(weighted_surrogate_index_2[j]) /
191 1148 : spline_vec_1_l->sample(weighted_surrogate_index_1[j]),
192 1148 : 1.0 / ((Real)num_layers - 1.0))
193 4723 : : bias_parameter;
194 : const Real index_factor =
195 4723 : MooseUtils::absoluteFuzzyEqual(l_ratio, 1.0)
196 4723 : ? (Real)i / (Real)num_layers
197 4723 : : (1.0 - std::pow(l_ratio, (Real)i)) / (1.0 - std::pow(l_ratio, (Real)num_layers));
198 4723 : Point tmp_point = surrogate_pos_2 * index_factor + surrogate_pos_1 * (1.0 - index_factor);
199 4723 : nodes[i][j] = mesh.add_point(tmp_point, j + node_counter);
200 : }
201 589 : node_counter += node_number_vec[i];
202 589 : }
203 : // Create triangular elements based on the 2D Node vector
204 131 : if (quad_elem)
205 49 : elementsCreationFromNodesVectorsQuad(mesh,
206 : nodes,
207 : num_layers,
208 : node_number_vec,
209 : transition_layer_id,
210 : input_boundary_1_id,
211 : input_boundary_2_id,
212 : begin_side_boundary_id,
213 : end_side_boundary_id);
214 : else
215 82 : elementsCreationFromNodesVectors(mesh,
216 : nodes,
217 : num_layers,
218 : node_number_vec,
219 : transition_layer_id,
220 : input_boundary_1_id,
221 : input_boundary_2_id,
222 : begin_side_boundary_id,
223 : end_side_boundary_id);
224 131 : }
225 :
226 : void
227 0 : fillBetweenPointVectorsGenerator(MeshBase & mesh,
228 : const std::vector<Point> & boundary_points_vec_1,
229 : const std::vector<Point> & boundary_points_vec_2,
230 : const unsigned int num_layers,
231 : const subdomain_id_type transition_layer_id,
232 : const boundary_id_type external_boundary_id,
233 : const std::string type,
234 : const std::string name,
235 : const bool quad_elem)
236 : {
237 0 : fillBetweenPointVectorsGenerator(mesh,
238 : boundary_points_vec_1,
239 : boundary_points_vec_2,
240 : num_layers,
241 : transition_layer_id,
242 : external_boundary_id,
243 : external_boundary_id,
244 : external_boundary_id,
245 : external_boundary_id,
246 : type,
247 : name,
248 : quad_elem);
249 0 : }
250 :
251 : void
252 49 : elementsCreationFromNodesVectorsQuad(MeshBase & mesh,
253 : const std::vector<std::vector<Node *>> & nodes,
254 : const unsigned int num_layers,
255 : const std::vector<unsigned int> & node_number_vec,
256 : const subdomain_id_type transition_layer_id,
257 : const boundary_id_type input_boundary_1_id,
258 : const boundary_id_type input_boundary_2_id,
259 : const boundary_id_type begin_side_boundary_id,
260 : const boundary_id_type end_side_boundary_id)
261 : {
262 49 : const unsigned int node_number = node_number_vec.front();
263 49 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
264 :
265 261 : for (const auto i : make_range(num_layers))
266 1736 : for (unsigned int j = 1; j < node_number; j++)
267 : {
268 1524 : Elem * elem = mesh.add_elem(new Quad4);
269 1524 : bool is_elem_flip = buildQuadElement(elem,
270 1524 : nodes[i][j - 1],
271 1524 : nodes[i + 1][j - 1],
272 1524 : nodes[i + 1][j],
273 1524 : nodes[i][j],
274 : transition_layer_id);
275 1524 : if (i == 0)
276 326 : boundary_info.add_side(elem, is_elem_flip ? 0 : 3, input_boundary_1_id);
277 1524 : if (i == num_layers - 1)
278 326 : boundary_info.add_side(elem, is_elem_flip ? 2 : 1, input_boundary_2_id);
279 1524 : if (j == 1)
280 212 : boundary_info.add_side(elem, is_elem_flip ? 3 : 0, begin_side_boundary_id);
281 1524 : if (j == node_number - 1)
282 212 : boundary_info.add_side(elem, is_elem_flip ? 1 : 2, end_side_boundary_id);
283 : }
284 49 : }
285 :
286 : void
287 82 : elementsCreationFromNodesVectors(MeshBase & mesh,
288 : const std::vector<std::vector<Node *>> & nodes,
289 : const unsigned int num_layers,
290 : const std::vector<unsigned int> & node_number_vec,
291 : const subdomain_id_type transition_layer_id,
292 : const boundary_id_type input_boundary_1_id,
293 : const boundary_id_type input_boundary_2_id,
294 : const boundary_id_type begin_side_boundary_id,
295 : const boundary_id_type end_side_boundary_id)
296 : {
297 82 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
298 :
299 328 : for (const auto i : make_range(num_layers))
300 : {
301 246 : unsigned int nodes_up_it = 0;
302 246 : unsigned int nodes_down_it = 0;
303 246 : const unsigned int node_number_up = node_number_vec[i + 1];
304 246 : const unsigned int node_number_down = node_number_vec[i];
305 :
306 3318 : while (nodes_up_it < node_number_up - 1 && nodes_down_it < node_number_down - 1 &&
307 3072 : nodes_up_it + nodes_down_it < node_number_up + node_number_down - 3)
308 : {
309 : // Define the two possible options and chose the one with shorter distance
310 3072 : Real dis1 = (*nodes[i + 1][nodes_up_it] - *nodes[i][nodes_down_it + 1]).norm();
311 3072 : Real dis2 = (*nodes[i + 1][nodes_up_it + 1] - *nodes[i][nodes_down_it]).norm();
312 3072 : if (MooseUtils::absoluteFuzzyGreaterThan(dis1, dis2))
313 : {
314 1552 : Elem * elem = mesh.add_elem(new Tri3);
315 1552 : bool is_elem_flip = buildTriElement(elem,
316 1552 : nodes[i + 1][nodes_up_it],
317 1552 : nodes[i][nodes_down_it],
318 1552 : nodes[i + 1][nodes_up_it + 1],
319 : transition_layer_id);
320 1552 : if (i == num_layers - 1)
321 452 : boundary_info.add_side(elem, is_elem_flip ? 0 : 2, input_boundary_2_id);
322 1552 : if (nodes_up_it == 0 && nodes_down_it == 0)
323 80 : boundary_info.add_side(elem, is_elem_flip ? 2 : 0, begin_side_boundary_id);
324 1552 : nodes_up_it++;
325 : }
326 : else
327 : {
328 1520 : Elem * elem = mesh.add_elem(new Tri3);
329 1520 : bool is_elem_flip = buildTriElement(elem,
330 1520 : nodes[i + 1][nodes_up_it],
331 1520 : nodes[i][nodes_down_it],
332 1520 : nodes[i][nodes_down_it + 1],
333 : transition_layer_id);
334 1520 : if (i == 0)
335 496 : boundary_info.add_side(elem, 1, input_boundary_1_id);
336 1520 : if (nodes_up_it == 0 && nodes_down_it == 0)
337 166 : boundary_info.add_side(elem, is_elem_flip ? 2 : 0, begin_side_boundary_id);
338 1520 : nodes_down_it++;
339 : }
340 : }
341 : // Handle the end
342 376 : while (nodes_up_it < node_number_up - 1)
343 : {
344 130 : Elem * elem = mesh.add_elem(new Tri3);
345 130 : bool is_elem_flip = buildTriElement(elem,
346 130 : nodes[i + 1][nodes_up_it],
347 130 : nodes[i][nodes_down_it],
348 130 : nodes[i + 1][nodes_up_it + 1],
349 : transition_layer_id);
350 130 : nodes_up_it++;
351 130 : if (i == num_layers - 1)
352 88 : boundary_info.add_side(elem, is_elem_flip ? 0 : 2, input_boundary_2_id);
353 130 : if (nodes_up_it == node_number_up - 1 && nodes_down_it == node_number_down - 1)
354 66 : boundary_info.add_side(elem, 1, end_side_boundary_id);
355 : }
356 470 : while (nodes_down_it < node_number_down - 1)
357 : {
358 224 : Elem * elem = mesh.add_elem(new Tri3);
359 224 : bool is_elem_flip = buildTriElement(elem,
360 224 : nodes[i + 1][nodes_up_it],
361 224 : nodes[i][nodes_down_it],
362 224 : nodes[i][nodes_down_it + 1],
363 : transition_layer_id);
364 224 : nodes_down_it++;
365 224 : if (i == 0)
366 106 : boundary_info.add_side(elem, 1, input_boundary_1_id);
367 224 : if (nodes_up_it == node_number_up - 1 && nodes_down_it == node_number_down - 1)
368 180 : boundary_info.add_side(elem, is_elem_flip ? 0 : 2, end_side_boundary_id);
369 : }
370 : }
371 82 : }
372 :
373 : void
374 262 : weightedInterpolator(const unsigned int vec_node_num,
375 : const std::vector<Point> & boundary_points_vec,
376 : std::vector<Real> & vec_index,
377 : std::vector<Real> & wt,
378 : std::vector<Real> & index,
379 : const Real sigma,
380 : std::unique_ptr<LinearInterpolation> & linear_vec_x,
381 : std::unique_ptr<LinearInterpolation> & linear_vec_y,
382 : std::unique_ptr<SplineInterpolation> & spline_vec_l)
383 : {
384 262 : std::vector<Real> pos_x;
385 262 : std::vector<Real> pos_y;
386 262 : std::vector<Real> dist_vec;
387 262 : std::vector<Real> pos_l;
388 :
389 2318 : for (const auto i : make_range(vec_node_num))
390 : {
391 : // Unweighted, the index interval is just uniform
392 : // Normalized range 0~1
393 2056 : vec_index.push_back((Real)i / ((Real)vec_node_num - 1.0));
394 : // X and Y coordinates cooresponding to the index
395 2056 : pos_x.push_back(boundary_points_vec[i](0));
396 2056 : pos_y.push_back(boundary_points_vec[i](1));
397 : // Use Point-to-Point distance as unnormalized weight
398 2056 : if (i > 0)
399 : {
400 1794 : wt.push_back((boundary_points_vec[i] - boundary_points_vec[i - 1]).norm());
401 1794 : dist_vec.push_back(wt.back());
402 : }
403 : // Accumulated unnormalized weights to get unnormalized weighted index
404 2056 : index.push_back(std::accumulate(wt.begin(), wt.end(), 0.0));
405 : }
406 262 : const Real dist_vec_total = index.back(); // Total accumulated distances
407 262 : const Real wt_norm_factor = dist_vec_total / ((Real)vec_node_num - 1.0); // Normalization factor
408 : // Normalization for both weights and weighted indices
409 262 : std::transform(
410 1794 : wt.begin(), wt.end(), wt.begin(), [wt_norm_factor](Real & c) { return c / wt_norm_factor; });
411 262 : std::transform(index.begin(),
412 : index.end(),
413 : index.begin(),
414 2056 : [dist_vec_total](Real & c) { return c / dist_vec_total; });
415 : // Use Gaussian blurring to smoothen local density
416 2318 : for (const auto i : make_range(vec_node_num))
417 : {
418 2056 : Real gaussian_factor(0.0);
419 2056 : Real sum_tmp(0.0);
420 : // Use interval as parameter now, consider distance in the future
421 17384 : for (const auto j : make_range(vec_node_num - 1))
422 : {
423 : // dis_vec and index are off by 0.5
424 : const Real tmp_factor =
425 15328 : exp(-((Real)(i - j) - 0.5) * ((Real)(i - j) - 0.5) / 2.0 / sigma / sigma);
426 15328 : gaussian_factor += tmp_factor;
427 15328 : sum_tmp += tmp_factor * dist_vec[j];
428 : }
429 2056 : pos_l.push_back(sum_tmp / gaussian_factor);
430 : }
431 : // Interpolate positions based on weighted indices
432 262 : linear_vec_x = std::make_unique<LinearInterpolation>(index, pos_x);
433 262 : linear_vec_y = std::make_unique<LinearInterpolation>(index, pos_y);
434 262 : spline_vec_l = std::make_unique<SplineInterpolation>(index, pos_l);
435 262 : }
436 :
437 : void
438 1178 : surrogateGenerator(std::vector<Real> & weighted_surrogate_index,
439 : std::vector<Real> & unweighted_surrogate_index,
440 : const std::vector<unsigned int> & node_number_vec,
441 : const std::vector<Real> & wt,
442 : const std::vector<Real> & index,
443 : const unsigned int boundary_node_num,
444 : const unsigned int i)
445 : {
446 : // First element is trivial
447 1178 : weighted_surrogate_index.push_back(0.0);
448 1178 : unweighted_surrogate_index.push_back(0.0);
449 9446 : for (unsigned int j = 1; j < node_number_vec[i]; j++)
450 : {
451 : // uniform interval for unweighted index
452 8268 : unweighted_surrogate_index.push_back((Real)j / ((Real)node_number_vec[i] - 1.0));
453 : // >
454 : const auto it_0 =
455 8268 : std::upper_bound(index.begin(), index.end(), unweighted_surrogate_index[j - 1]);
456 : // >=
457 8268 : const auto it_1 = std::lower_bound(index.begin(), index.end(), unweighted_surrogate_index[j]);
458 : //
459 8268 : const auto it_dist = std::distance(it_0, it_1);
460 : //
461 8268 : const auto it_start = std::distance(index.begin(), it_0);
462 :
463 8268 : if (it_0 == it_1)
464 11700 : weighted_surrogate_index.push_back(weighted_surrogate_index[j - 1] +
465 5850 : wt[it_start - 1] / ((Real)node_number_vec[i] - 1.0));
466 : else
467 : {
468 2418 : weighted_surrogate_index.push_back(weighted_surrogate_index[j - 1]);
469 2418 : weighted_surrogate_index[j] += (*it_0 - unweighted_surrogate_index[j - 1]) * wt[it_start - 1];
470 4836 : weighted_surrogate_index[j] +=
471 2418 : (unweighted_surrogate_index[j] - *(it_1 - 1)) * wt[it_start + it_dist - 1];
472 2716 : for (unsigned int k = 1; k < it_dist; k++)
473 298 : weighted_surrogate_index[j] += wt[it_start + k - 1] / ((Real)boundary_node_num - 1.0);
474 : }
475 : }
476 1178 : }
477 :
478 : bool
479 131 : needFlip(const std::vector<Point> & vec_pts_1, const std::vector<Point> & vec_pts_2)
480 : {
481 : const Real th1 =
482 131 : acos((vec_pts_1.back() - vec_pts_1.front()) * (vec_pts_2.front() - vec_pts_1.front()) /
483 131 : (vec_pts_1.back() - vec_pts_1.front()).norm() /
484 131 : (vec_pts_2.front() - vec_pts_1.front()).norm());
485 131 : const Real th2 = acos(
486 131 : (vec_pts_2.back() - vec_pts_1.back()) * (vec_pts_1.front() - vec_pts_1.back()) /
487 131 : (vec_pts_2.back() - vec_pts_1.back()).norm() / (vec_pts_1.front() - vec_pts_1.back()).norm());
488 131 : const Real th3 = acos(
489 131 : (vec_pts_2.front() - vec_pts_2.back()) * (vec_pts_1.back() - vec_pts_2.back()) /
490 131 : (vec_pts_2.front() - vec_pts_2.back()).norm() / (vec_pts_1.back() - vec_pts_2.back()).norm());
491 : const Real th4 =
492 131 : acos((vec_pts_1.front() - vec_pts_2.front()) * (vec_pts_2.back() - vec_pts_2.front()) /
493 131 : (vec_pts_1.front() - vec_pts_2.front()).norm() /
494 131 : (vec_pts_2.back() - vec_pts_2.front()).norm());
495 131 : if (MooseUtils::absoluteFuzzyEqual(th1 + th2 + th3 + th4, 2 * M_PI))
496 131 : return false;
497 0 : return true;
498 : }
499 :
500 : bool
501 150 : isBoundarySimpleClosedLoop(MeshBase & mesh,
502 : Real & max_node_radius,
503 : std::vector<dof_id_type> & boundary_ordered_node_list,
504 : const Point origin_pt,
505 : const boundary_id_type bid)
506 : {
507 : // This has no communication and expects elem_ptr to find any
508 : // element, so it only works on serialized meshes
509 150 : libMesh::MeshSerializer serial(mesh);
510 :
511 150 : max_node_radius = 0.0;
512 150 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
513 150 : auto side_list_tmp = boundary_info.build_side_list();
514 150 : std::vector<std::pair<dof_id_type, dof_id_type>> boundary_node_assm;
515 150 : std::vector<dof_id_type> boundary_midpoint_node_list;
516 3454 : for (const auto i : index_range(side_list_tmp))
517 : {
518 3304 : if (std::get<2>(side_list_tmp[i]) == bid)
519 : {
520 : // store two nodes of each side
521 1156 : const auto elem = mesh.elem_ptr(std::get<0>(side_list_tmp[i]));
522 1156 : const auto side = elem->side_ptr(std::get<1>(side_list_tmp[i]));
523 1156 : boundary_node_assm.push_back(std::make_pair(side->node_id(0), side->node_id(1)));
524 : // see if there is a midpoint
525 1156 : const auto & side_type = elem->side_type(std::get<1>(side_list_tmp[i]));
526 1156 : if (side_type == EDGE3)
527 0 : boundary_midpoint_node_list.push_back(
528 0 : elem->node_id(elem->n_vertices() + std::get<1>(side_list_tmp[i])));
529 : else
530 1156 : boundary_midpoint_node_list.push_back(DofObject::invalid_id);
531 1156 : }
532 : }
533 : bool is_closed_loop;
534 442 : isClosedLoop(mesh,
535 : max_node_radius,
536 : boundary_ordered_node_list,
537 : boundary_node_assm,
538 : boundary_midpoint_node_list,
539 : origin_pt,
540 : "external boundary",
541 : is_closed_loop);
542 4 : return is_closed_loop;
543 588 : }
544 :
545 : bool
546 0 : isBoundarySimpleClosedLoop(MeshBase & mesh,
547 : Real & max_node_radius,
548 : const Point origin_pt,
549 : const boundary_id_type bid)
550 : {
551 0 : std::vector<dof_id_type> dummy_boundary_ordered_node_list;
552 0 : return isBoundarySimpleClosedLoop(
553 0 : mesh, max_node_radius, dummy_boundary_ordered_node_list, origin_pt, bid);
554 0 : }
555 :
556 : bool
557 0 : isBoundarySimpleClosedLoop(MeshBase & mesh, const Point origin_pt, const boundary_id_type bid)
558 : {
559 : Real dummy_max_node_radius;
560 0 : return isBoundarySimpleClosedLoop(mesh, dummy_max_node_radius, origin_pt, bid);
561 : }
562 :
563 : bool
564 150 : isBoundaryOpenSingleSegment(MeshBase & mesh,
565 : Real & max_node_radius,
566 : std::vector<dof_id_type> & boundary_ordered_node_list,
567 : const Point origin_pt,
568 : const boundary_id_type bid)
569 : {
570 : try
571 : {
572 150 : isBoundarySimpleClosedLoop(mesh, max_node_radius, boundary_ordered_node_list, origin_pt, bid);
573 : }
574 146 : catch (MooseException & e)
575 : {
576 292 : if (((std::string)e.what())
577 146 : .compare("This mesh generator does not work for the provided external boundary as it "
578 146 : "is not a closed loop.") != 0)
579 4 : throw MooseException("The provided boundary is not an open single-segment boundary.");
580 : else
581 142 : return true;
582 146 : }
583 :
584 4 : throw MooseException("The provided boundary is closed loop, which is not supported.");
585 : }
586 :
587 : bool
588 158 : isExternalBoundary(MeshBase & mesh, const boundary_id_type bid)
589 : {
590 : // This has no communication and expects elem_ptr to find any
591 : // element, so it only works on serialized meshes
592 158 : libMesh::MeshSerializer serial(mesh);
593 :
594 158 : if (!mesh.is_prepared())
595 158 : mesh.find_neighbors();
596 158 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
597 158 : auto side_list = boundary_info.build_side_list();
598 3622 : for (const auto i : index_range(side_list))
599 : {
600 3464 : if (std::get<2>(side_list[i]) == bid)
601 1236 : if (mesh.elem_ptr(std::get<0>(side_list[i]))->neighbor_ptr(std::get<1>(side_list[i])) !=
602 : nullptr)
603 0 : return false;
604 : }
605 158 : return true;
606 158 : }
607 :
608 : bool
609 40 : isCurveSimpleClosedLoop(MeshBase & mesh,
610 : Real & max_node_radius,
611 : std::vector<dof_id_type> & ordered_node_list,
612 : const Point origin_pt)
613 : {
614 : // This has no communication and expects to loop over all elements
615 : // on every processor, so it only works on serialized meshes
616 40 : libMesh::MeshSerializer serial(mesh);
617 :
618 40 : max_node_radius = 0.0;
619 40 : std::vector<std::pair<dof_id_type, dof_id_type>> node_assm;
620 280 : for (auto it = mesh.active_elements_begin(); it != mesh.active_elements_end(); it++)
621 280 : node_assm.push_back(std::make_pair((*it)->node_id(0), (*it)->node_id(1)));
622 : bool is_closed_loop;
623 120 : isClosedLoop(
624 : mesh, max_node_radius, ordered_node_list, node_assm, origin_pt, "curve", is_closed_loop);
625 0 : return is_closed_loop;
626 80 : }
627 :
628 : bool
629 0 : isCurveSimpleClosedLoop(MeshBase & mesh, Real & max_node_radius, const Point origin_pt)
630 : {
631 0 : std::vector<dof_id_type> dummy_ordered_node_list;
632 0 : return isCurveSimpleClosedLoop(mesh, max_node_radius, dummy_ordered_node_list, origin_pt);
633 0 : }
634 :
635 : bool
636 0 : isCurveSimpleClosedLoop(MeshBase & mesh, const Point origin_pt)
637 : {
638 : Real dummy_max_node_radius;
639 0 : return isCurveSimpleClosedLoop(mesh, dummy_max_node_radius, origin_pt);
640 : }
641 :
642 : bool
643 40 : isCurveOpenSingleSegment(MeshBase & mesh,
644 : Real & max_node_radius,
645 : std::vector<dof_id_type> & ordered_node_list,
646 : const Point origin_pt)
647 : {
648 : try
649 : {
650 40 : isCurveSimpleClosedLoop(mesh, max_node_radius, ordered_node_list, origin_pt);
651 : }
652 40 : catch (MooseException & e)
653 : {
654 80 : if (((std::string)e.what())
655 40 : .compare("This mesh generator does not work for the provided curve as it is not a "
656 40 : "closed loop.") != 0)
657 0 : throw MooseException("The provided curve is not an open single-segment boundary.");
658 : else
659 40 : return true;
660 40 : }
661 0 : throw MooseException("The provided curve is closed loop, which is not supported.");
662 : return false;
663 : }
664 :
665 : void
666 246 : isClosedLoop(MeshBase & mesh,
667 : Real & max_node_radius,
668 : std::vector<dof_id_type> & ordered_node_list,
669 : std::vector<std::pair<dof_id_type, dof_id_type>> & node_assm,
670 : std::vector<dof_id_type> & midpoint_node_list,
671 : const Point origin_pt,
672 : const std::string input_type,
673 : bool & is_closed_loop,
674 : const bool suppress_exception)
675 : {
676 : // This has no communication and expects node_ptr to find any
677 : // node, so it only works on serialized meshes
678 246 : libMesh::MeshSerializer serial(mesh);
679 :
680 246 : std::vector<dof_id_type> dummy_elem_list = std::vector<dof_id_type>(node_assm.size(), 0);
681 246 : std::vector<dof_id_type> ordered_dummy_elem_list;
682 246 : is_closed_loop = false;
683 246 : MooseMeshUtils::makeOrderedNodeList(
684 : node_assm, dummy_elem_list, midpoint_node_list, ordered_node_list, ordered_dummy_elem_list);
685 : // If the code ever gets here, node_assm is empty.
686 : // If the ordered_node_list front and back are not the same, the boundary is not a loop.
687 : // This is not done inside the loop just for some potential applications in the future.
688 242 : if (ordered_node_list.front() != ordered_node_list.back())
689 : {
690 : // This is invalid type #2
691 202 : if (!suppress_exception)
692 182 : throw MooseException("This mesh generator does not work for the provided ",
693 : input_type,
694 364 : " as it is not a closed loop.");
695 : }
696 : // It the curve is a loop, check if azimuthal angles change monotonically
697 : else
698 : {
699 : // Utilize cross product here.
700 : // If azimuthal angles change monotonically,
701 : // the z components of the cross products are always negative or positive.
702 40 : std::vector<Real> ordered_node_azi_list;
703 538 : for (const auto i : make_range(ordered_node_list.size() - 1))
704 : {
705 498 : ordered_node_azi_list.push_back(
706 498 : (*mesh.node_ptr(ordered_node_list[i]) - origin_pt)
707 498 : .cross(*mesh.node_ptr(ordered_node_list[i + 1]) - origin_pt)(2));
708 : // Use this opportunity to calculate maximum radius
709 498 : max_node_radius =
710 498 : std::max((*mesh.node_ptr(ordered_node_list[i]) - origin_pt).norm(), max_node_radius);
711 : }
712 40 : std::sort(ordered_node_azi_list.begin(), ordered_node_azi_list.end());
713 40 : if (ordered_node_azi_list.front() * ordered_node_azi_list.back() < 0.0)
714 : {
715 : // This is invalid type #3
716 0 : if (!suppress_exception)
717 0 : throw MooseException(
718 : "This mesh generator does not work for the provided ",
719 : input_type,
720 0 : " as azimuthal angles of consecutive nodes do not change monotonically.");
721 : }
722 : else
723 40 : is_closed_loop = true;
724 40 : }
725 618 : }
726 :
727 : void
728 96 : isClosedLoop(MeshBase & mesh,
729 : Real & max_node_radius,
730 : std::vector<dof_id_type> & ordered_node_list,
731 : std::vector<std::pair<dof_id_type, dof_id_type>> & node_assm,
732 : const Point origin_pt,
733 : const std::string input_type,
734 : bool & is_closed_loop,
735 : const bool suppress_exception)
736 : {
737 96 : std::vector<dof_id_type> dummy_midpoint_node_list(node_assm.size(), DofObject::invalid_id);
738 136 : isClosedLoop(mesh,
739 : max_node_radius,
740 : ordered_node_list,
741 : node_assm,
742 : dummy_midpoint_node_list,
743 : origin_pt,
744 : input_type,
745 : is_closed_loop,
746 : suppress_exception);
747 96 : }
748 :
749 : bool
750 1524 : buildQuadElement(Elem * elem,
751 : Node * nd_0,
752 : Node * nd_1,
753 : Node * nd_2,
754 : Node * nd_3,
755 : const subdomain_id_type transition_layer_id)
756 : {
757 : // Adjust the order of nodes in an element so that the mesh can be extruded in (0 0 1)
758 : // direction.
759 1524 : elem->subdomain_id() = transition_layer_id;
760 1524 : if (((*nd_1 - *nd_0).cross((*nd_2 - *nd_0)).unit())(2) > 0)
761 : {
762 936 : elem->set_node(0, nd_0);
763 936 : elem->set_node(1, nd_1);
764 936 : elem->set_node(2, nd_2);
765 936 : elem->set_node(3, nd_3);
766 936 : return false;
767 : }
768 : else
769 : {
770 588 : elem->set_node(0, nd_0);
771 588 : elem->set_node(3, nd_1);
772 588 : elem->set_node(2, nd_2);
773 588 : elem->set_node(1, nd_3);
774 588 : return true;
775 : }
776 : }
777 :
778 : bool
779 3426 : buildTriElement(
780 : Elem * elem, Node * nd_0, Node * nd_1, Node * nd_2, const subdomain_id_type transition_layer_id)
781 : {
782 : // Adjust the order of nodes in an element so that the mesh can be extruded in (0 0 1)
783 : // direction.
784 3426 : elem->subdomain_id() = transition_layer_id;
785 3426 : if (((*nd_1 - *nd_0).cross((*nd_2 - *nd_0)).unit())(2) > 0)
786 : {
787 1170 : elem->set_node(0, nd_0);
788 1170 : elem->set_node(1, nd_1);
789 1170 : elem->set_node(2, nd_2);
790 1170 : return false;
791 : }
792 : else
793 : {
794 2256 : elem->set_node(0, nd_0);
795 2256 : elem->set_node(2, nd_1);
796 2256 : elem->set_node(1, nd_2);
797 2256 : return true;
798 : }
799 : }
800 : }
|
