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 "LineSegment.h"
11 : #include "MooseError.h"
12 : #include "ParameterMesh.h"
13 :
14 : #include "libmesh/cell_hex8.h"
15 : #include "libmesh/cell_hex.h"
16 : #include "libmesh/edge_edge2.h"
17 : #include "libmesh/enum_elem_type.h"
18 : #include "libmesh/enum_point_locator_type.h"
19 : #include "libmesh/dof_map.h"
20 :
21 : #include "libmesh/elem.h"
22 : #include "libmesh/face_quad.h"
23 : #include "libmesh/face_quad4.h"
24 : #include "libmesh/fe_compute_data.h"
25 : #include "libmesh/fe_interface.h"
26 : #include "libmesh/id_types.h"
27 : #include "libmesh/int_range.h"
28 : #include "libmesh/libmesh_common.h"
29 : #include "libmesh/numeric_vector.h"
30 : #include "libmesh/explicit_system.h"
31 : #include "libmesh/plane.h"
32 : #include "libmesh/enum_to_string.h"
33 : #include <memory>
34 :
35 : using namespace libMesh;
36 :
37 508 : ParameterMesh::ParameterMesh(const FEType & param_type,
38 : const std::string & exodus_mesh,
39 : const std::vector<std::string> & var_names,
40 : const bool find_closest,
41 508 : const unsigned int kdtree_candidates)
42 508 : : _communicator(MPI_COMM_SELF),
43 508 : _mesh(_communicator),
44 508 : _find_closest(find_closest),
45 508 : _kdtree_candidates(kdtree_candidates)
46 : {
47 : _mesh.allow_renumbering(false);
48 508 : _mesh.prepare_for_use();
49 1016 : _exodusII_io = std::make_unique<ExodusII_IO>(_mesh);
50 508 : _exodusII_io->read(exodus_mesh);
51 508 : _mesh.read(exodus_mesh);
52 : // Create system to store parameter values
53 1016 : _eq = std::make_unique<EquationSystems>(_mesh);
54 508 : _sys = &_eq->add_system<ExplicitSystem>("_parameter_mesh_sys");
55 508 : _sys->add_variable("_parameter_mesh_var", param_type);
56 :
57 : // Create point locator
58 1016 : _point_locator = PointLocatorBase::build(TREE_LOCAL_ELEMENTS, _mesh);
59 508 : _point_locator->enable_out_of_mesh_mode();
60 :
61 508 : if (!var_names.empty())
62 : {
63 : // Make Exodus vars for equation system
64 87 : const std::vector<std::string> & all_nodal(_exodusII_io->get_nodal_var_names());
65 87 : const std::vector<std::string> & all_elemental(_exodusII_io->get_elem_var_names());
66 :
67 : std::vector<std::string> nodal_variables;
68 : std::vector<std::string> elemental_variables;
69 297 : for (const auto & var_name : var_names)
70 : {
71 210 : if (std::find(all_nodal.begin(), all_nodal.end(), var_name) != all_nodal.end())
72 161 : nodal_variables.push_back(var_name);
73 210 : if (std::find(all_elemental.begin(), all_elemental.end(), var_name) != all_elemental.end())
74 47 : elemental_variables.push_back(var_name);
75 : }
76 87 : if ((elemental_variables.size() + nodal_variables.size()) != var_names.size())
77 : {
78 2 : std::string out("\n Parameter Group Variables Requested: ");
79 4 : for (const auto & var : var_names)
80 4 : out += var + " ";
81 : out += "\n Exodus Nodal Variables: ";
82 14 : for (const auto & var : all_nodal)
83 24 : out += var + " ";
84 : out += "\n Exodus Elemental Variables: ";
85 8 : for (const auto & var : all_elemental)
86 12 : out += var + " ";
87 2 : mooseError("Exodus file did not contain all of the parameter names being intitialized.", out);
88 : }
89 :
90 85 : if (!nodal_variables.empty() && !elemental_variables.empty())
91 : {
92 2 : std::string out("\n Parameter Group Nodal Variables: ");
93 6 : for (const auto & var : nodal_variables)
94 8 : out += var + " ";
95 : out += "\n Parameter Group Elemental Variables: ";
96 4 : for (const auto & var : elemental_variables)
97 4 : out += var + " ";
98 2 : mooseError("Parameter group contains nodal and elemental variables, this is "
99 : "not allowed. ",
100 : out);
101 : }
102 : // Add the parameter group ics and bounds to the system
103 : // All parameters in a group will be the same type, only one of these loops will do anything
104 240 : for (const auto & var_name : nodal_variables)
105 157 : _sys->add_variable(var_name, param_type);
106 128 : for (const auto & var_name : elemental_variables)
107 45 : _sys->add_variable(var_name, param_type);
108 83 : }
109 : // Initialize the equations systems
110 504 : _eq->init();
111 :
112 504 : const unsigned short int var_id = _sys->variable_number("_parameter_mesh_var");
113 : std::set<dof_id_type> var_indices;
114 504 : _sys->local_dof_indices(var_id, var_indices);
115 504 : _param_dofs = var_indices.size();
116 :
117 504 : if (_find_closest)
118 : {
119 147264 : for (const auto & elem : _mesh.element_ptr_range())
120 73542 : if (elem->default_order() != FIRST)
121 92 : mooseError("Closet point projection currently does not support second order elements.");
122 : }
123 :
124 : // Initialize node-based KDTree optimization
125 502 : _mesh_nodes.clear();
126 : _node_to_elements.clear();
127 :
128 : // Extract all node coordinates
129 66888 : for (const auto & node : _mesh.node_ptr_range())
130 33444 : _mesh_nodes.push_back(*node);
131 :
132 : // Build node-to-elements connectivity map
133 151984 : for (const auto & elem : _mesh.element_ptr_range())
134 : {
135 383100 : for (const auto n : make_range(elem->n_nodes()))
136 : {
137 307610 : dof_id_type node_id = elem->node_id(n);
138 307610 : _node_to_elements[node_id].insert(elem);
139 : }
140 502 : }
141 :
142 : // Create KDTree from node coordinates
143 502 : if (!_mesh_nodes.empty())
144 1004 : _node_kdtree = std::make_unique<KDTree>(_mesh_nodes, 10);
145 502 : }
146 :
147 : void
148 6936316 : ParameterMesh::getIndexAndWeight(const Point & pt,
149 : std::vector<dof_id_type> & dof_indices,
150 : std::vector<Real> & weights) const
151 : {
152 6936316 : Point test_point = (_find_closest ? projectToMesh(pt) : pt);
153 :
154 6936316 : const Elem * elem = (*_point_locator)(test_point);
155 6936316 : if (!elem)
156 0 : mooseError("No element was found to contain point ", test_point);
157 :
158 : // Get the in the dof_indices for our element
159 : // variable name is hard coded to _parameter_mesh_var
160 : // this is probably the only variable in the ParameterMesh system used by ParameterMeshFunction
161 6936316 : const unsigned short int var = _sys->variable_number("_parameter_mesh_var");
162 6936316 : const DofMap & dof_map = _sys->get_dof_map();
163 6936316 : dof_map.dof_indices(elem, dof_indices, var);
164 :
165 : // Map the physical co-ordinates to the reference co-ordinates
166 6936316 : Point coor = FEMap::inverse_map(elem->dim(), elem, test_point);
167 : // get the shape function value via the FEInterface
168 6936316 : FEType fe_type = dof_map.variable_type(var);
169 6936316 : FEComputeData fe_data(*_eq, coor);
170 6936316 : FEInterface::compute_data(elem->dim(), fe_type, elem, fe_data);
171 : // Set weights to the value of the shape functions
172 6936316 : weights = fe_data.shape;
173 :
174 6936316 : if (dof_indices.size() != weights.size())
175 0 : mooseError("Internal error: weights and DoF indices do not have the same size.");
176 6936316 : }
177 :
178 : void
179 16200 : ParameterMesh::getIndexAndWeight(const Point & pt,
180 : std::vector<dof_id_type> & dof_indices,
181 : std::vector<RealGradient> & weights) const
182 : {
183 16200 : if (!_sys->has_variable("_parameter_mesh_var"))
184 0 : mooseError("Internal error: System being read does not contain _parameter_mesh_var.");
185 16200 : Point test_point = (_find_closest ? projectToMesh(pt) : pt);
186 : // Locate the element the point is in
187 16200 : const Elem * elem = (*_point_locator)(test_point);
188 :
189 : // Get the in the dof_indices for our element
190 : // variable name is hard coded to _parameter_mesh_var
191 : // this is probably the only variable in the ParameterMesh system used by ParameterMeshFunction
192 16200 : const unsigned short int var = _sys->variable_number("_parameter_mesh_var");
193 16200 : const DofMap & dof_map = _sys->get_dof_map();
194 16200 : dof_map.dof_indices(elem, dof_indices, var);
195 :
196 : // Map the physical co-ordinates to the reference co-ordinates
197 16200 : Point coor = FEMap::inverse_map(elem->dim(), elem, test_point);
198 : // get the shape function value via the FEInterface
199 32400 : FEType fe_type = dof_map.variable_type(var);
200 16200 : FEComputeData fe_data(*_eq, coor);
201 16200 : fe_data.enable_derivative();
202 16200 : FEInterface::compute_data(elem->dim(), fe_type, elem, fe_data);
203 : // Set weights to the value of the shape functions
204 16200 : weights = fe_data.dshape;
205 :
206 16200 : if (dof_indices.size() != weights.size())
207 0 : mooseError("Internal error: weights and DoF indices do not have the same size.");
208 16200 : }
209 :
210 : std::vector<Real>
211 366 : ParameterMesh::getParameterValues(std::string var_name, unsigned int time_step) const
212 : {
213 366 : if (!_sys->has_variable(var_name))
214 0 : mooseError("Exodus file being read does not contain ", var_name, ".");
215 : // get the exodus variable and put it into the equation system.
216 366 : unsigned int step_to_read = _exodusII_io->get_num_time_steps();
217 366 : if (time_step <= step_to_read)
218 : step_to_read = time_step;
219 110 : else if (time_step != std::numeric_limits<unsigned int>::max())
220 2 : mooseError("Invalid value passed as \"time_step\". Expected a valid integer "
221 : "less than ",
222 2 : _exodusII_io->get_num_time_steps(),
223 : ", received ",
224 : time_step);
225 :
226 : // determine what kind of variable you are trying to read from mesh
227 364 : const std::vector<std::string> & all_nodal(_exodusII_io->get_nodal_var_names());
228 364 : const std::vector<std::string> & all_elemental(_exodusII_io->get_elem_var_names());
229 364 : if (std::find(all_nodal.begin(), all_nodal.end(), var_name) != all_nodal.end())
230 530 : _exodusII_io->copy_nodal_solution(*_sys, var_name, var_name, step_to_read);
231 99 : else if (std::find(all_elemental.begin(), all_elemental.end(), var_name) != all_elemental.end())
232 198 : _exodusII_io->copy_elemental_solution(*_sys, var_name, var_name, step_to_read);
233 :
234 : // Update the equations systems
235 364 : _sys->update();
236 :
237 364 : const unsigned short int var_id = _sys->variable_number(var_name);
238 : std::set<dof_id_type> var_indices;
239 364 : _sys->local_dof_indices(var_id, var_indices); // Everything is local so this is fine
240 364 : std::vector<dof_id_type> var_indices_vector(var_indices.begin(), var_indices.end());
241 :
242 : std::vector<Real> parameter_values;
243 364 : _sys->solution->localize(parameter_values, var_indices_vector);
244 364 : return parameter_values;
245 364 : }
246 :
247 : Point
248 41730 : ParameterMesh::projectToMesh(const Point & p) const
249 : {
250 : // quick path: p already inside an element
251 41730 : if ((*_point_locator)(p))
252 1800 : return p;
253 :
254 : // Lambda to find closest point from elements using squared distance for efficiency
255 39930 : auto findClosestElement = [&p, this](const auto & elements) -> Point
256 : {
257 : Real best_d2 = std::numeric_limits<Real>::max();
258 39930 : Point best_point = p;
259 :
260 793866 : for (const auto * elem : elements)
261 : {
262 753936 : Point trial = closestPoint(*elem, p);
263 753936 : Real d2 = (trial - p).norm_sq();
264 753936 : if (d2 < best_d2)
265 : {
266 : best_d2 = d2;
267 146059 : best_point = trial;
268 : }
269 : }
270 :
271 39930 : if (best_d2 == std::numeric_limits<Real>::max())
272 0 : mooseError("project_to_mesh failed - no candidate elements.");
273 39930 : return best_point;
274 39930 : };
275 :
276 : // Use KDTree optimization if available
277 39930 : if (_node_kdtree && !_mesh_nodes.empty())
278 : {
279 : // Find K nearest nodes using KDTree
280 : std::vector<std::size_t> nearest_node_indices;
281 39930 : _node_kdtree->neighborSearch(p, _kdtree_candidates, nearest_node_indices);
282 :
283 : // Collect all elements connected to these nodes
284 : std::set<const Elem *> candidate_elements;
285 239580 : for (auto node_idx : nearest_node_indices)
286 : {
287 : // Get the actual node from the mesh using the index
288 199650 : if (node_idx < _mesh.n_nodes())
289 : {
290 199650 : const Node * node = _mesh.node_ptr(node_idx);
291 199650 : dof_id_type node_id = node->id();
292 : auto it = _node_to_elements.find(node_id);
293 199650 : if (it != _node_to_elements.end())
294 : {
295 : const auto & connected_elems = it->second;
296 199650 : candidate_elements.insert(connected_elems.begin(), connected_elems.end());
297 : }
298 : }
299 : }
300 :
301 : // Convert set to vector for consistent type
302 : std::vector<const Elem *> candidate_vector(candidate_elements.begin(),
303 39930 : candidate_elements.end());
304 39930 : return findClosestElement(candidate_vector);
305 79860 : }
306 : else
307 : {
308 : // Fallback to original O(n) method if KDTree not available
309 : std::vector<const Elem *> all_elements;
310 0 : for (const auto & elem : _mesh.element_ptr_range())
311 0 : all_elements.push_back(elem);
312 :
313 0 : return findClosestElement(all_elements);
314 0 : }
315 : }
316 :
317 : Point
318 11197791 : ParameterMesh::closestPoint(const Elem & elem, const Point & p) const
319 : {
320 : mooseAssert(!elem.contains_point(p),
321 : "Points inside of elements shouldn't need to find closestPoint.");
322 :
323 : // Lambda to find closest point from range without storing temporary vectors
324 3101880 : auto findClosest = [&p](auto range, auto point_func) -> Point
325 : {
326 : Real min_distance = std::numeric_limits<Real>::max();
327 3101880 : Point min_point = p;
328 :
329 13545735 : for (const auto & item : range)
330 : {
331 10443855 : Point candidate = point_func(item);
332 10443855 : Real distance = (candidate - p).norm();
333 10443855 : if (distance < min_distance)
334 : {
335 : min_distance = distance;
336 4158198 : min_point = candidate;
337 : }
338 : }
339 3101880 : return min_point;
340 11197791 : };
341 :
342 11197791 : switch (elem.type())
343 : {
344 : case EDGE2:
345 : {
346 8056704 : LineSegment ls(*(elem.node_ptr(0)), *(elem.node_ptr(1)));
347 8056704 : return ls.closest_point(p);
348 : }
349 :
350 : case TRI3:
351 : {
352 2126376 : Point a = *(elem.node_ptr(0));
353 2126376 : Point b = *(elem.node_ptr(1));
354 2126376 : Point c = *(elem.node_ptr(2));
355 2126376 : Plane pl(a, b, c);
356 2126376 : Point trial = pl.closest_point(p);
357 2126376 : if (elem.contains_point(trial))
358 28899 : return trial;
359 :
360 2097477 : return findClosest(make_range(elem.n_edges()),
361 6292431 : [&](dof_id_type i) { return closestPoint(*elem.build_edge_ptr(i), p); });
362 2126376 : }
363 : case QUAD4:
364 : {
365 448800 : Point a = *(elem.node_ptr(0));
366 448800 : Point b = *(elem.node_ptr(1));
367 448800 : Point c = *(elem.node_ptr(2));
368 448800 : Point d = *(elem.node_ptr(3));
369 448800 : Plane pl1(a, b, c);
370 448800 : Plane pl2(b, c, d);
371 448800 : Point trial1 = pl1.closest_point(p);
372 448800 : Point trial2 = pl2.closest_point(p);
373 448800 : if (!trial1.absolute_fuzzy_equals(trial2, TOLERANCE * TOLERANCE))
374 0 : mooseError("Quad4 element is not coplanar");
375 :
376 448800 : if (elem.contains_point(trial1))
377 10308 : return trial1;
378 :
379 438492 : return findClosest(make_range(elem.n_edges()),
380 1753968 : [&](dof_id_type i) { return closestPoint(*elem.build_edge_ptr(i), p); });
381 448800 : }
382 :
383 565911 : default:
384 : {
385 565911 : if (elem.dim() == 3)
386 : {
387 565911 : return findClosest(make_range(elem.n_sides()),
388 2397456 : [&](dof_id_type i) { return closestPoint(*elem.build_side_ptr(i), p); });
389 : }
390 : else
391 : {
392 0 : mooseError("Unsupported element type ",
393 0 : Utility::enum_to_string(elem.type()),
394 : " for projection of parameter mesh.");
395 : }
396 : }
397 : }
398 : }
|
