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 153 : 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 306 : if (!MooseMeshUtils::isCoPlanar(
48 459 : boundary_points_vec_1, Point(0.0, 0.0, 1.0), Point(0.0, 0.0, 0.0)) ||
49 306 : !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 149 : const unsigned int vec_1_node_num = boundary_points_vec_1.size();
60 149 : const unsigned int vec_2_node_num = boundary_points_vec_2.size();
61 :
62 149 : 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 145 : 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 141 : std::vector<Point> possibly_reoriented_boundary_points_vec_2;
84 141 : const std::vector<Point> * oriented_boundary_points_vec_2 = &boundary_points_vec_2;
85 :
86 141 : 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 141 : std::vector<Real> vec_1_index; // Unweighted index
103 141 : std::vector<Real> vec_2_index; // Unweighted index
104 :
105 141 : std::vector<Real> wt_1;
106 141 : std::vector<Real> index_1; // Weighted index
107 141 : std::vector<Real> wt_2;
108 141 : std::vector<Real> index_2; // Weighted index
109 :
110 : // create interpolations
111 141 : std::unique_ptr<LinearInterpolation> linear_vec_1_x;
112 141 : std::unique_ptr<LinearInterpolation> linear_vec_1_y;
113 141 : std::unique_ptr<SplineInterpolation> spline_vec_1_l;
114 141 : std::unique_ptr<LinearInterpolation> linear_vec_2_x;
115 141 : std::unique_ptr<LinearInterpolation> linear_vec_2_y;
116 141 : std::unique_ptr<SplineInterpolation> spline_vec_2_l;
117 :
118 141 : 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 141 : 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 141 : const Real increment = ((Real)vec_2_node_num - (Real)vec_1_node_num) / (Real)(num_layers);
140 : // Number of nodes in each sublayer
141 141 : std::vector<unsigned int> node_number_vec;
142 : // 2D vector of nodes
143 141 : std::vector<std::vector<Node *>> nodes(num_layers + 1);
144 : // Node counter
145 141 : unsigned int node_counter = 0;
146 :
147 770 : for (const auto i : make_range(num_layers + 1))
148 : {
149 : // calculate number of nodes in each sublayer
150 629 : node_number_vec.push_back(
151 629 : (unsigned int)(vec_1_node_num + (long)(increment * i + 0.5 - (increment < 0))));
152 : // Reserve memory for new nodes
153 629 : nodes[i] = std::vector<Node *>(node_number_vec[i]);
154 :
155 : // Calculate vectors of weighted surrogated index for side #1
156 629 : std::vector<Real> weighted_surrogate_index_1;
157 629 : std::vector<Real> unweighted_surrogate_index_1;
158 :
159 629 : 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 629 : std::vector<Real> weighted_surrogate_index_2;
169 629 : std::vector<Real> unweighted_surrogate_index_2;
170 :
171 629 : 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 5662 : 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 5033 : Point surrogate_pos_1 = Point(linear_vec_1_x->sample(weighted_surrogate_index_1[j]),
183 5033 : linear_vec_1_y->sample(weighted_surrogate_index_1[j]),
184 15099 : 0.0);
185 : // Create surrogate Points on side #2 for Point #j on the sublayer
186 5033 : Point surrogate_pos_2 = Point(linear_vec_2_x->sample(weighted_surrogate_index_2[j]),
187 5033 : linear_vec_2_y->sample(weighted_surrogate_index_2[j]),
188 15099 : 0.0);
189 : const Real l_ratio = bias_parameter <= 0.0
190 6303 : ? std::pow(spline_vec_2_l->sample(weighted_surrogate_index_2[j]) /
191 1270 : spline_vec_1_l->sample(weighted_surrogate_index_1[j]),
192 1270 : 1.0 / ((Real)num_layers - 1.0))
193 5033 : : bias_parameter;
194 : const Real index_factor =
195 5033 : MooseUtils::absoluteFuzzyEqual(l_ratio, 1.0)
196 5033 : ? (Real)i / (Real)num_layers
197 5033 : : (1.0 - std::pow(l_ratio, (Real)i)) / (1.0 - std::pow(l_ratio, (Real)num_layers));
198 5033 : Point tmp_point = surrogate_pos_2 * index_factor + surrogate_pos_1 * (1.0 - index_factor);
199 5033 : nodes[i][j] = mesh.add_point(tmp_point, j + node_counter);
200 : }
201 629 : node_counter += node_number_vec[i];
202 629 : }
203 : // Create triangular elements based on the 2D Node vector
204 141 : if (quad_elem)
205 50 : 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 91 : 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 141 : }
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 50 : 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 50 : const unsigned int node_number = node_number_vec.front();
263 50 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
264 :
265 265 : for (const auto i : make_range(num_layers))
266 1751 : for (unsigned int j = 1; j < node_number; j++)
267 : {
268 1536 : Elem * elem = mesh.add_elem(new Quad4);
269 1536 : bool is_elem_flip = buildQuadElement(elem,
270 1536 : nodes[i][j - 1],
271 1536 : nodes[i + 1][j - 1],
272 1536 : nodes[i + 1][j],
273 1536 : nodes[i][j],
274 : transition_layer_id);
275 1536 : if (i == 0)
276 330 : boundary_info.add_side(elem, is_elem_flip ? 0 : 3, input_boundary_1_id);
277 1536 : if (i == num_layers - 1)
278 330 : boundary_info.add_side(elem, is_elem_flip ? 2 : 1, input_boundary_2_id);
279 1536 : if (j == 1)
280 215 : boundary_info.add_side(elem, is_elem_flip ? 3 : 0, begin_side_boundary_id);
281 1536 : if (j == node_number - 1)
282 215 : boundary_info.add_side(elem, is_elem_flip ? 1 : 2, end_side_boundary_id);
283 : }
284 50 : }
285 :
286 : void
287 91 : 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 91 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
298 :
299 364 : for (const auto i : make_range(num_layers))
300 : {
301 273 : unsigned int nodes_up_it = 0;
302 273 : unsigned int nodes_down_it = 0;
303 273 : const unsigned int node_number_up = node_number_vec[i + 1];
304 273 : const unsigned int node_number_down = node_number_vec[i];
305 :
306 3686 : while (nodes_up_it < node_number_up - 1 && nodes_down_it < node_number_down - 1 &&
307 3413 : 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 3413 : Real dis1 = (*nodes[i + 1][nodes_up_it] - *nodes[i][nodes_down_it + 1]).norm();
311 3413 : Real dis2 = (*nodes[i + 1][nodes_up_it + 1] - *nodes[i][nodes_down_it]).norm();
312 3413 : if (MooseUtils::absoluteFuzzyGreaterThan(dis1, dis2))
313 : {
314 1725 : Elem * elem = mesh.add_elem(new Tri3);
315 1725 : bool is_elem_flip = buildTriElement(elem,
316 1725 : nodes[i + 1][nodes_up_it],
317 1725 : nodes[i][nodes_down_it],
318 1725 : nodes[i + 1][nodes_up_it + 1],
319 : transition_layer_id);
320 1725 : if (i == num_layers - 1)
321 503 : boundary_info.add_side(elem, is_elem_flip ? 0 : 2, input_boundary_2_id);
322 1725 : if (nodes_up_it == 0 && nodes_down_it == 0)
323 89 : boundary_info.add_side(elem, is_elem_flip ? 2 : 0, begin_side_boundary_id);
324 1725 : nodes_up_it++;
325 : }
326 : else
327 : {
328 1688 : Elem * elem = mesh.add_elem(new Tri3);
329 1688 : bool is_elem_flip = buildTriElement(elem,
330 1688 : nodes[i + 1][nodes_up_it],
331 1688 : nodes[i][nodes_down_it],
332 1688 : nodes[i][nodes_down_it + 1],
333 : transition_layer_id);
334 1688 : if (i == 0)
335 549 : boundary_info.add_side(elem, 1, input_boundary_1_id);
336 1688 : if (nodes_up_it == 0 && nodes_down_it == 0)
337 184 : boundary_info.add_side(elem, is_elem_flip ? 2 : 0, begin_side_boundary_id);
338 1688 : nodes_down_it++;
339 : }
340 : }
341 : // Handle the end
342 419 : while (nodes_up_it < node_number_up - 1)
343 : {
344 146 : Elem * elem = mesh.add_elem(new Tri3);
345 146 : bool is_elem_flip = buildTriElement(elem,
346 146 : nodes[i + 1][nodes_up_it],
347 146 : nodes[i][nodes_down_it],
348 146 : nodes[i + 1][nodes_up_it + 1],
349 : transition_layer_id);
350 146 : nodes_up_it++;
351 146 : if (i == num_layers - 1)
352 99 : boundary_info.add_side(elem, is_elem_flip ? 0 : 2, input_boundary_2_id);
353 146 : if (nodes_up_it == node_number_up - 1 && nodes_down_it == node_number_down - 1)
354 74 : boundary_info.add_side(elem, 1, end_side_boundary_id);
355 : }
356 521 : while (nodes_down_it < node_number_down - 1)
357 : {
358 248 : Elem * elem = mesh.add_elem(new Tri3);
359 248 : bool is_elem_flip = buildTriElement(elem,
360 248 : nodes[i + 1][nodes_up_it],
361 248 : nodes[i][nodes_down_it],
362 248 : nodes[i][nodes_down_it + 1],
363 : transition_layer_id);
364 248 : nodes_down_it++;
365 248 : if (i == 0)
366 118 : boundary_info.add_side(elem, 1, input_boundary_1_id);
367 248 : if (nodes_up_it == node_number_up - 1 && nodes_down_it == node_number_down - 1)
368 199 : boundary_info.add_side(elem, is_elem_flip ? 0 : 2, end_side_boundary_id);
369 : }
370 : }
371 91 : }
372 :
373 : void
374 282 : 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 282 : std::vector<Real> pos_x;
385 282 : std::vector<Real> pos_y;
386 282 : std::vector<Real> dist_vec;
387 282 : std::vector<Real> pos_l;
388 :
389 2493 : for (const auto i : make_range(vec_node_num))
390 : {
391 : // Unweighted, the index interval is just uniform
392 : // Normalized range 0~1
393 2211 : vec_index.push_back((Real)i / ((Real)vec_node_num - 1.0));
394 : // X and Y coordinates cooresponding to the index
395 2211 : pos_x.push_back(boundary_points_vec[i](0));
396 2211 : pos_y.push_back(boundary_points_vec[i](1));
397 : // Use Point-to-Point distance as unnormalized weight
398 2211 : if (i > 0)
399 : {
400 1929 : wt.push_back((boundary_points_vec[i] - boundary_points_vec[i - 1]).norm());
401 1929 : dist_vec.push_back(wt.back());
402 : }
403 : // Accumulated unnormalized weights to get unnormalized weighted index
404 2211 : index.push_back(std::accumulate(wt.begin(), wt.end(), 0.0));
405 : }
406 282 : const Real dist_vec_total = index.back(); // Total accumulated distances
407 282 : 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 282 : std::transform(
410 1929 : wt.begin(), wt.end(), wt.begin(), [wt_norm_factor](Real & c) { return c / wt_norm_factor; });
411 282 : std::transform(index.begin(),
412 : index.end(),
413 : index.begin(),
414 2211 : [dist_vec_total](Real & c) { return c / dist_vec_total; });
415 : // Use Gaussian blurring to smoothen local density
416 2493 : for (const auto i : make_range(vec_node_num))
417 : {
418 2211 : Real gaussian_factor(0.0);
419 2211 : Real sum_tmp(0.0);
420 : // Use interval as parameter now, consider distance in the future
421 18689 : 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 16478 : exp(-((Real)(i - j) - 0.5) * ((Real)(i - j) - 0.5) / 2.0 / sigma / sigma);
426 16478 : gaussian_factor += tmp_factor;
427 16478 : sum_tmp += tmp_factor * dist_vec[j];
428 : }
429 2211 : pos_l.push_back(sum_tmp / gaussian_factor);
430 : }
431 : // Interpolate positions based on weighted indices
432 282 : linear_vec_x = std::make_unique<LinearInterpolation>(index, pos_x);
433 282 : linear_vec_y = std::make_unique<LinearInterpolation>(index, pos_y);
434 282 : spline_vec_l = std::make_unique<SplineInterpolation>(index, pos_l);
435 282 : }
436 :
437 : void
438 1258 : 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 1258 : weighted_surrogate_index.push_back(0.0);
448 1258 : unweighted_surrogate_index.push_back(0.0);
449 10066 : for (unsigned int j = 1; j < node_number_vec[i]; j++)
450 : {
451 : // uniform interval for unweighted index
452 8808 : unweighted_surrogate_index.push_back((Real)j / ((Real)node_number_vec[i] - 1.0));
453 : // >
454 : const auto it_0 =
455 8808 : std::upper_bound(index.begin(), index.end(), unweighted_surrogate_index[j - 1]);
456 : // >=
457 8808 : const auto it_1 = std::lower_bound(index.begin(), index.end(), unweighted_surrogate_index[j]);
458 : //
459 8808 : const auto it_dist = std::distance(it_0, it_1);
460 : //
461 8808 : const auto it_start = std::distance(index.begin(), it_0);
462 :
463 8808 : if (it_0 == it_1)
464 12242 : weighted_surrogate_index.push_back(weighted_surrogate_index[j - 1] +
465 6121 : wt[it_start - 1] / ((Real)node_number_vec[i] - 1.0));
466 : else
467 : {
468 2687 : weighted_surrogate_index.push_back(weighted_surrogate_index[j - 1]);
469 2687 : weighted_surrogate_index[j] += (*it_0 - unweighted_surrogate_index[j - 1]) * wt[it_start - 1];
470 5374 : weighted_surrogate_index[j] +=
471 2687 : (unweighted_surrogate_index[j] - *(it_1 - 1)) * wt[it_start + it_dist - 1];
472 3018 : for (unsigned int k = 1; k < it_dist; k++)
473 331 : weighted_surrogate_index[j] += wt[it_start + k - 1] / ((Real)boundary_node_num - 1.0);
474 : }
475 : }
476 1258 : }
477 :
478 : bool
479 141 : needFlip(const std::vector<Point> & vec_pts_1, const std::vector<Point> & vec_pts_2)
480 : {
481 : const Real th1 =
482 141 : acos((vec_pts_1.back() - vec_pts_1.front()) * (vec_pts_2.front() - vec_pts_1.front()) /
483 141 : (vec_pts_1.back() - vec_pts_1.front()).norm() /
484 141 : (vec_pts_2.front() - vec_pts_1.front()).norm());
485 141 : const Real th2 = acos(
486 141 : (vec_pts_2.back() - vec_pts_1.back()) * (vec_pts_1.front() - vec_pts_1.back()) /
487 141 : (vec_pts_2.back() - vec_pts_1.back()).norm() / (vec_pts_1.front() - vec_pts_1.back()).norm());
488 141 : const Real th3 = acos(
489 141 : (vec_pts_2.front() - vec_pts_2.back()) * (vec_pts_1.back() - vec_pts_2.back()) /
490 141 : (vec_pts_2.front() - vec_pts_2.back()).norm() / (vec_pts_1.back() - vec_pts_2.back()).norm());
491 : const Real th4 =
492 141 : acos((vec_pts_1.front() - vec_pts_2.front()) * (vec_pts_2.back() - vec_pts_2.front()) /
493 141 : (vec_pts_1.front() - vec_pts_2.front()).norm() /
494 141 : (vec_pts_2.back() - vec_pts_2.front()).norm());
495 141 : if (MooseUtils::absoluteFuzzyEqual(th1 + th2 + th3 + th4, 2 * M_PI))
496 141 : return false;
497 0 : return true;
498 : }
499 :
500 : bool
501 158 : 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 158 : libMesh::MeshSerializer serial(mesh);
510 :
511 158 : max_node_radius = 0.0;
512 158 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
513 158 : auto side_list_tmp = boundary_info.build_side_list();
514 158 : std::vector<std::pair<dof_id_type, dof_id_type>> boundary_node_assm;
515 158 : std::vector<dof_id_type> boundary_midpoint_node_list;
516 3590 : for (const auto i : index_range(side_list_tmp))
517 : {
518 3432 : if (std::get<2>(side_list_tmp[i]) == bid)
519 : {
520 : // store two nodes of each side
521 1220 : const auto elem = mesh.elem_ptr(std::get<0>(side_list_tmp[i]));
522 1220 : const auto side = elem->side_ptr(std::get<1>(side_list_tmp[i]));
523 1220 : boundary_node_assm.push_back(std::make_pair(side->node_id(0), side->node_id(1)));
524 : // see if there is a midpoint
525 1220 : const auto & side_type = elem->side_type(std::get<1>(side_list_tmp[i]));
526 1220 : 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 1220 : boundary_midpoint_node_list.push_back(DofObject::invalid_id);
531 1220 : }
532 : }
533 : bool is_closed_loop;
534 466 : 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 620 : }
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 158 : 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 158 : isBoundarySimpleClosedLoop(mesh, max_node_radius, boundary_ordered_node_list, origin_pt, bid);
573 : }
574 154 : catch (MooseException & e)
575 : {
576 308 : if (((std::string)e.what())
577 154 : .compare("This mesh generator does not work for the provided external boundary as it "
578 154 : "is not a closed loop.") != 0)
579 4 : throw MooseException("The provided boundary is not an open single-segment boundary.");
580 : else
581 150 : return true;
582 154 : }
583 :
584 4 : throw MooseException("The provided boundary is closed loop, which is not supported.");
585 : }
586 :
587 : bool
588 166 : 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 166 : libMesh::MeshSerializer serial(mesh);
593 :
594 166 : if (!mesh.is_prepared())
595 166 : mesh.find_neighbors();
596 166 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
597 166 : auto side_list = boundary_info.build_side_list();
598 3758 : for (const auto i : index_range(side_list))
599 : {
600 3592 : if (std::get<2>(side_list[i]) == bid)
601 1300 : 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 166 : return true;
606 166 : }
607 :
608 : bool
609 44 : 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 44 : libMesh::MeshSerializer serial(mesh);
617 :
618 44 : max_node_radius = 0.0;
619 44 : std::vector<std::pair<dof_id_type, dof_id_type>> node_assm;
620 308 : for (auto it = mesh.active_elements_begin(); it != mesh.active_elements_end(); it++)
621 308 : node_assm.push_back(std::make_pair((*it)->node_id(0), (*it)->node_id(1)));
622 : bool is_closed_loop;
623 132 : isClosedLoop(
624 : mesh, max_node_radius, ordered_node_list, node_assm, origin_pt, "curve", is_closed_loop);
625 0 : return is_closed_loop;
626 88 : }
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 44 : 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 44 : isCurveSimpleClosedLoop(mesh, max_node_radius, ordered_node_list, origin_pt);
651 : }
652 44 : catch (MooseException & e)
653 : {
654 88 : if (((std::string)e.what())
655 44 : .compare("This mesh generator does not work for the provided curve as it is not a "
656 44 : "closed loop.") != 0)
657 0 : throw MooseException("The provided curve is not an open single-segment boundary.");
658 : else
659 44 : return true;
660 44 : }
661 0 : throw MooseException("The provided curve is closed loop, which is not supported.");
662 : return false;
663 : }
664 :
665 : void
666 262 : 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 262 : libMesh::MeshSerializer serial(mesh);
679 :
680 262 : std::vector<dof_id_type> dummy_elem_list = std::vector<dof_id_type>(node_assm.size(), 0);
681 262 : std::vector<dof_id_type> ordered_dummy_elem_list;
682 262 : is_closed_loop = false;
683 262 : 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 258 : if (ordered_node_list.front() != ordered_node_list.back())
689 : {
690 : // This is invalid type #2
691 216 : if (!suppress_exception)
692 194 : throw MooseException("This mesh generator does not work for the provided ",
693 : input_type,
694 388 : " 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 42 : std::vector<Real> ordered_node_azi_list;
703 565 : for (const auto i : make_range(ordered_node_list.size() - 1))
704 : {
705 523 : ordered_node_azi_list.push_back(
706 523 : (*mesh.node_ptr(ordered_node_list[i]) - origin_pt)
707 523 : .cross(*mesh.node_ptr(ordered_node_list[i + 1]) - origin_pt)(2));
708 : // Use this opportunity to calculate maximum radius
709 523 : max_node_radius =
710 523 : std::max((*mesh.node_ptr(ordered_node_list[i]) - origin_pt).norm(), max_node_radius);
711 : }
712 42 : std::sort(ordered_node_azi_list.begin(), ordered_node_azi_list.end());
713 42 : 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 42 : is_closed_loop = true;
724 42 : }
725 658 : }
726 :
727 : void
728 104 : 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 104 : std::vector<dof_id_type> dummy_midpoint_node_list(node_assm.size(), DofObject::invalid_id);
738 148 : 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 104 : }
748 :
749 : bool
750 1536 : 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 1536 : elem->subdomain_id() = transition_layer_id;
760 1536 : 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 600 : elem->set_node(0, nd_0);
771 600 : elem->set_node(3, nd_1);
772 600 : elem->set_node(2, nd_2);
773 600 : elem->set_node(1, nd_3);
774 600 : return true;
775 : }
776 : }
777 :
778 : bool
779 3807 : 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 3807 : elem->subdomain_id() = transition_layer_id;
785 3807 : if (((*nd_1 - *nd_0).cross((*nd_2 - *nd_0)).unit())(2) > 0)
786 : {
787 1287 : elem->set_node(0, nd_0);
788 1287 : elem->set_node(1, nd_1);
789 1287 : elem->set_node(2, nd_2);
790 1287 : return false;
791 : }
792 : else
793 : {
794 2520 : elem->set_node(0, nd_0);
795 2520 : elem->set_node(2, nd_1);
796 2520 : elem->set_node(1, nd_2);
797 2520 : return true;
798 : }
799 : }
800 : }
|
