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 "PolyLineMeshFollowingNodeSetGenerator.h"
11 :
12 : #include "CastUniquePointer.h"
13 : #include "MooseMeshUtils.h"
14 : #include "MooseUtils.h"
15 :
16 : #include "libmesh/elem.h"
17 : #include "libmesh/int_range.h"
18 : #include "libmesh/unstructured_mesh.h"
19 :
20 : registerMooseObject("MooseApp", PolyLineMeshFollowingNodeSetGenerator);
21 :
22 : InputParameters
23 3077 : PolyLineMeshFollowingNodeSetGenerator::validParams()
24 : {
25 3077 : InputParameters params = MeshGenerator::validParams();
26 :
27 : // Path parameters
28 12308 : params.addRequiredParam<Point>("starting_point", "Starting point for the polyline");
29 12308 : params.addRequiredParam<Point>("starting_direction",
30 : "Starting value for the direction of the line");
31 12308 : params.addRequiredParam<MeshGeneratorName>("input", "The mesh we get the sideset from");
32 12308 : params.addRequiredParam<BoundaryName>("nodeset", "Nodeset to follow to form the polyline");
33 12308 : params.addRequiredParam<Real>("search_radius",
34 : "Radius of the sphere used to find points in the nodeset");
35 9231 : params.addParam<bool>(
36 : "ignore_nodes_behind",
37 6154 : false,
38 : "Ignore nodes in the nodeset that are behind the current point in the polyline");
39 12308 : params.addParam<bool>("loop", false, "Whether edges should form a closed loop");
40 :
41 : // Discretization parameters
42 : // NOTE: we could have another dx as a path search parameter, and decouple the two options
43 18462 : params.addRequiredRangeCheckedParam<Real>(
44 : "dx",
45 : "dx>0",
46 : "Approximate size of the edge elements (before any refinement with num_edges_between_points) "
47 : "and approximate step to advance by to search for the next point in the polyline");
48 9231 : params.addParam<unsigned int>(
49 6154 : "max_edges", 1000, "Maximum number of edges. Serves as a stopping criterion");
50 9231 : params.addParam<unsigned int>(
51 6154 : "num_edges_between_points", 1, "How many Edge elements to build between each point pair");
52 :
53 : // Naming of result parameters
54 12308 : params.addParam<SubdomainName>("line_subdomain", "line", "Subdomain name for the line");
55 12308 : params.addParam<BoundaryName>(
56 : "start_boundary", "start", "Boundary to assign to (non-looped) polyline start");
57 12308 : params.addParam<BoundaryName>(
58 : "end_boundary", "end", "Boundary to assign to (non-looped) polyline end");
59 :
60 9231 : params.addParam<bool>(
61 : "verbose",
62 6154 : false,
63 : "whether to output additional information to console during the line generation");
64 :
65 3077 : params.addClassDescription(
66 : "Generates a polyline (open ended or looped) of Edge elements by marching along a nodeeset "
67 : "and trying to be as close as possible to the nodes of the nodeset");
68 :
69 3077 : return params;
70 0 : }
71 :
72 8 : PolyLineMeshFollowingNodeSetGenerator::PolyLineMeshFollowingNodeSetGenerator(
73 8 : const InputParameters & parameters)
74 : : MeshGenerator(parameters),
75 8 : _input(getMesh("input")),
76 16 : _starting_point(getParam<Point>("starting_point")),
77 16 : _starting_direction(getParam<Point>("starting_direction")),
78 16 : _ignore_nodes_behind(getParam<bool>("ignore_nodes_behind")),
79 16 : _loop(getParam<bool>("loop")),
80 16 : _line_subdomain(getParam<SubdomainName>("line_subdomain")),
81 16 : _start_boundary(getParam<BoundaryName>("start_boundary")),
82 16 : _end_boundary(getParam<BoundaryName>("end_boundary")),
83 16 : _dx(getParam<Real>("dx")),
84 16 : _num_edges_between_points(getParam<unsigned int>("num_edges_between_points")),
85 24 : _verbose(getParam<bool>("verbose"))
86 : {
87 8 : if (_loop && (isParamSetByUser("start_boundary") || isParamSetByUser("end_boundary")))
88 0 : paramError("loop",
89 : "Loop does not have a start or end boundary. These parameters must not be passed.");
90 8 : }
91 :
92 : std::unique_ptr<MeshBase>
93 8 : PolyLineMeshFollowingNodeSetGenerator::generate()
94 : {
95 8 : auto uptr_mesh = buildMeshBaseObject();
96 8 : MeshBase & mesh = *uptr_mesh;
97 8 : std::unique_ptr<MeshBase> base_mesh = std::move(_input);
98 8 : if (!base_mesh->is_serial())
99 0 : paramError("input", "Input mesh must not be distributed");
100 :
101 : // We may rely on boundary info caches later
102 8 : if (!base_mesh->preparation().has_boundary_id_sets)
103 0 : base_mesh->get_boundary_info().regenerate_id_sets();
104 :
105 : // Get nodeset ID in input mesh
106 : const auto nodeset_id =
107 16 : MooseMeshUtils::getBoundaryID(getParam<BoundaryName>("nodeset"), *base_mesh);
108 :
109 16 : const auto search_radius_sq = libMesh::Utility::pow<2>(getParam<Real>("search_radius"));
110 16 : const auto n_points = getParam<unsigned int>("max_edges");
111 8 : Point current_point = _starting_point;
112 8 : Point previous_direction = _starting_direction;
113 8 : mesh.add_point(_starting_point, 0);
114 :
115 : // Pre-find all points in the nodeset to follow
116 : // TODO: Build a KNN tree to speed up the search later on
117 8 : const auto all_nodeset_tuples = base_mesh->get_boundary_info().build_node_list(
118 8 : libMesh::BoundaryInfo::NodeBCTupleSortBy::BOUNDARY_ID);
119 8 : std::vector<dof_id_type> nodeset_nodes;
120 16 : nodeset_nodes.reserve(all_nodeset_tuples.size() /
121 8 : base_mesh->get_boundary_info().n_boundary_ids());
122 680 : for (const auto & tup : all_nodeset_tuples)
123 672 : if (BoundaryID(std::get<1>(tup)) == nodeset_id)
124 304 : nodeset_nodes.push_back(std::get<0>(tup));
125 8 : if (_verbose)
126 24 : _console << "Total number of nodes in nodeset " << getParam<BoundaryName>("nodeset") << ": "
127 8 : << nodeset_nodes.size() << std::endl;
128 :
129 8 : unsigned int n_segments = 0;
130 328 : for (const auto i : make_range(n_points))
131 : {
132 : // Move the point forward in the search direction
133 320 : const auto previous_point = current_point;
134 320 : current_point += previous_direction * _dx;
135 :
136 : // Draw a sphere to find the next point as the barycenter of the nodes from the nodeset inside
137 : // the sphere.
138 : // NOTE: there are many heuristics we could use here. We could try to fit a cylinder
139 : // to the nodeset if we know the nodeset represents a cylinder for example.
140 320 : Point barycenter(0);
141 320 : unsigned int n_sum = 0;
142 12480 : for (const auto n_id : nodeset_nodes)
143 12160 : if ((current_point - base_mesh->node_ref(n_id)).norm_sq() < search_radius_sq)
144 : {
145 12336 : if (!_ignore_nodes_behind ||
146 12336 : ((base_mesh->node_ref(n_id) - current_point) * previous_direction >= 0))
147 : {
148 3496 : barycenter += base_mesh->node_ref(n_id);
149 3496 : n_sum++;
150 : }
151 : }
152 320 : if (n_sum > 0)
153 320 : barycenter /= n_sum;
154 0 : else if (!_ignore_nodes_behind)
155 0 : mooseError("Did not find any nodes in the nodeset near the current point at: ",
156 : current_point);
157 : else
158 0 : barycenter = previous_point;
159 :
160 320 : if (MooseUtils::absoluteFuzzyEqual((barycenter - previous_point).norm_sq(), 0))
161 : {
162 0 : mooseInfo("Barycenter did not move from ", barycenter, ". Returning!");
163 0 : goto done_drawing;
164 : }
165 :
166 : // Compute new direction
167 320 : const auto new_direction = (barycenter - previous_point).unit();
168 320 : previous_direction = new_direction;
169 :
170 : // Set the new point towards the barycenter
171 320 : n_segments++;
172 : // Note: this dx could be different than the dx used to search the barycenter
173 320 : current_point = previous_point + _dx * new_direction;
174 320 : mesh.add_point(current_point, (i + 1) * _num_edges_between_points);
175 320 : if (_verbose)
176 320 : _console << i << ": new point: " << current_point << " new direction " << previous_direction
177 320 : << std::endl;
178 :
179 : // Add the additional edges in between if requested
180 320 : if (_num_edges_between_points > 1)
181 : {
182 320 : auto p = previous_point;
183 320 : const Point pvec = (current_point - previous_point) / _num_edges_between_points;
184 640 : for (auto j : make_range(1u, _num_edges_between_points))
185 : {
186 320 : p += pvec;
187 320 : mesh.add_point(p, i * _num_edges_between_points + j);
188 : }
189 : }
190 : }
191 :
192 8 : done_drawing:
193 :
194 8 : const auto n_elem = n_segments * _num_edges_between_points - 1;
195 8 : const auto max_nodes = n_segments * _num_edges_between_points - 1;
196 640 : for (auto i : make_range(n_elem + _loop))
197 : {
198 632 : const auto ip1 = _loop ? (i + 1) % max_nodes : (i + 1);
199 632 : auto elem = Elem::build(EDGE2);
200 632 : elem->set_node(0, mesh.node_ptr(i));
201 632 : elem->set_node(1, mesh.node_ptr(ip1));
202 632 : elem->set_id() = i;
203 632 : mesh.add_elem(std::move(elem));
204 632 : }
205 :
206 : // Add the starting and end boundary
207 8 : if (!_loop)
208 : {
209 8 : BoundaryInfo & bi = mesh.get_boundary_info();
210 32 : std::vector<BoundaryName> bdy_names{_start_boundary, _end_boundary};
211 8 : std::vector<boundary_id_type> ids = MooseMeshUtils::getBoundaryIDs(mesh, bdy_names, true);
212 8 : bi.add_side(mesh.elem_ptr(0), 0, ids[0]);
213 8 : bi.add_side(mesh.elem_ptr(n_elem - 1), 1, ids[1]);
214 8 : }
215 :
216 8 : mesh.prepare_for_use();
217 :
218 16 : return uptr_mesh;
219 16 : }
|
