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/int_range.h"
20 : #include "libmesh/dof_map.h"
21 :
22 : #include "libmesh/elem.h"
23 : #include "libmesh/face_quad.h"
24 : #include "libmesh/face_quad4.h"
25 : #include "libmesh/fe_compute_data.h"
26 : #include "libmesh/fe_interface.h"
27 : #include "libmesh/id_types.h"
28 : #include "libmesh/int_range.h"
29 : #include "libmesh/libmesh_common.h"
30 : #include "libmesh/numeric_vector.h"
31 : #include "libmesh/explicit_system.h"
32 : #include "libmesh/plane.h"
33 : #include "libmesh/enum_to_string.h"
34 : #include <memory>
35 : #include "libmesh/quadrature_gauss.h"
36 : #include "libmesh/fe_base.h"
37 :
38 : using namespace libMesh;
39 :
40 546 : ParameterMesh::ParameterMesh(const FEType & param_type,
41 : const std::string & exodus_mesh,
42 : const bool find_closest,
43 546 : const unsigned int kdtree_candidates)
44 546 : : _communicator(MPI_COMM_SELF),
45 546 : _mesh(_communicator),
46 546 : _find_closest(find_closest),
47 546 : _kdtree_candidates(kdtree_candidates),
48 546 : _param_var_id(0),
49 546 : _dof_map(nullptr),
50 546 : _fe_type(param_type)
51 : {
52 : _mesh.allow_renumbering(false);
53 546 : _mesh.prepare_for_use();
54 1092 : _exodusII_io = std::make_unique<ExodusII_IO>(_mesh);
55 546 : _exodusII_io->read(exodus_mesh);
56 546 : _mesh.read(exodus_mesh);
57 : // Create system to store parameter values
58 1092 : _eq = std::make_unique<EquationSystems>(_mesh);
59 546 : _sys = &_eq->add_system<ExplicitSystem>("_parameter_mesh_sys");
60 546 : _sys->add_variable("_parameter_mesh_var", param_type);
61 :
62 : // Create point locator
63 1092 : _point_locator = PointLocatorBase::build(TREE_LOCAL_ELEMENTS, _mesh);
64 546 : _point_locator->enable_out_of_mesh_mode();
65 :
66 : // Initialize the equations systems
67 546 : _eq->init();
68 :
69 : // getting number of parameter dofs for size() function
70 546 : const unsigned short int var_id = _sys->variable_number("_parameter_mesh_var");
71 : std::set<dof_id_type> var_indices;
72 546 : _sys->local_dof_indices(var_id, var_indices);
73 546 : _param_dofs = var_indices.size();
74 :
75 546 : if (_find_closest)
76 : {
77 147264 : for (const auto & elem : _mesh.element_ptr_range())
78 73542 : if (elem->default_order() != FIRST)
79 92 : mooseError("Closet point projection currently does not support second order elements.");
80 : }
81 :
82 : // Initialize node-based KDTree optimization
83 544 : _mesh_nodes.clear();
84 : _node_to_elements.clear();
85 :
86 : // Extract all node coordinates
87 80952 : for (const auto & node : _mesh.node_ptr_range())
88 40476 : _mesh_nodes.push_back(*node);
89 :
90 : // Build node-to-elements connectivity map
91 163276 : for (const auto & elem : _mesh.element_ptr_range())
92 : {
93 492234 : for (const auto n : make_range(elem->n_nodes()))
94 : {
95 330046 : dof_id_type node_id = elem->node_id(n);
96 330046 : _node_to_elements[node_id].insert(elem);
97 : }
98 544 : }
99 :
100 : // Create KDTree from node coordinates
101 544 : if (!_mesh_nodes.empty())
102 1088 : _node_kdtree = std::make_unique<KDTree>(_mesh_nodes, 10);
103 : // Update cached values for gradient computations
104 544 : const_cast<unsigned short int &>(_param_var_id) = var_id;
105 544 : const_cast<const DofMap *&>(_dof_map) = &_sys->get_dof_map();
106 544 : const_cast<FEType &>(_fe_type) = _dof_map->variable_type(_param_var_id);
107 544 : }
108 :
109 : void
110 12974412 : ParameterMesh::getIndexAndWeight(const Point & pt,
111 : std::vector<dof_id_type> & dof_indices,
112 : std::vector<Real> & weights) const
113 : {
114 12974412 : Point test_point = (_find_closest ? projectToMesh(pt) : pt);
115 :
116 12974412 : const Elem * elem = (*_point_locator)(test_point);
117 12974412 : if (!elem)
118 0 : mooseError("No element was found to contain point ", test_point);
119 :
120 : // Get the dof_indices for our element
121 : // variable id is hard coded to _param_var_id
122 : // this is probably the only variable in the ParameterMesh system used by ParameterMeshFunction
123 12974412 : _dof_map->dof_indices(elem, dof_indices, _param_var_id);
124 :
125 : // Map the physical co-ordinates to the reference co-ordinates
126 12974412 : Point coor = FEMap::inverse_map(elem->dim(), elem, test_point);
127 : // get the shape function value via the FEInterface
128 12974412 : FEComputeData fe_data(*_eq, coor);
129 12974412 : FEInterface::compute_data(elem->dim(), _fe_type, elem, fe_data);
130 : // Set weights to the value of the shape functions
131 12974412 : weights = fe_data.shape;
132 :
133 12974412 : if (dof_indices.size() != weights.size())
134 0 : mooseError("Internal error: weights and DoF indices do not have the same size.");
135 12974412 : }
136 :
137 : void
138 16200 : ParameterMesh::getIndexAndWeight(const Point & pt,
139 : std::vector<dof_id_type> & dof_indices,
140 : std::vector<RealGradient> & weights) const
141 : {
142 16200 : if (!_sys->has_variable("_parameter_mesh_var"))
143 0 : mooseError("Internal error: System being read does not contain _parameter_mesh_var.");
144 16200 : Point test_point = (_find_closest ? projectToMesh(pt) : pt);
145 : // Locate the element the point is in
146 16200 : const Elem * elem = (*_point_locator)(test_point);
147 :
148 : // Get the dof_indices for our element
149 : // variable id is hard coded to _param_var_id
150 : // this is probably the only variable in the ParameterMesh system used by ParameterMeshFunction
151 16200 : _dof_map->dof_indices(elem, dof_indices, _param_var_id);
152 :
153 : // Map the physical co-ordinates to the reference co-ordinates
154 16200 : Point coor = FEMap::inverse_map(elem->dim(), elem, test_point);
155 : // get the shape function value via the FEInterface
156 16200 : FEComputeData fe_data(*_eq, coor);
157 16200 : fe_data.enable_derivative();
158 16200 : FEInterface::compute_data(elem->dim(), _fe_type, elem, fe_data);
159 : // Set weights to the value of the shape functions
160 16200 : weights = fe_data.dshape;
161 :
162 16200 : if (dof_indices.size() != weights.size())
163 0 : mooseError("Internal error: weights and DoF indices do not have the same size.");
164 16200 : }
165 :
166 : Point
167 41730 : ParameterMesh::projectToMesh(const Point & p) const
168 : {
169 : // quick path: p already inside an element
170 41730 : if ((*_point_locator)(p))
171 1800 : return p;
172 :
173 : // Lambda to find closest point from elements using squared distance for efficiency
174 39930 : auto findClosestElement = [&p, this](const auto & elements) -> Point
175 : {
176 : Real best_d2 = std::numeric_limits<Real>::max();
177 39930 : Point best_point = p;
178 :
179 793866 : for (const auto * elem : elements)
180 : {
181 753936 : Point trial = closestPoint(*elem, p);
182 753936 : Real d2 = (trial - p).norm_sq();
183 753936 : if (d2 < best_d2)
184 : {
185 : best_d2 = d2;
186 146059 : best_point = trial;
187 : }
188 : }
189 :
190 39930 : if (best_d2 == std::numeric_limits<Real>::max())
191 0 : mooseError("project_to_mesh failed - no candidate elements.");
192 39930 : return best_point;
193 39930 : };
194 :
195 : // Use KDTree optimization if available
196 39930 : if (_node_kdtree && !_mesh_nodes.empty())
197 : {
198 : // Find K nearest nodes using KDTree
199 : std::vector<std::size_t> nearest_node_indices;
200 39930 : _node_kdtree->neighborSearch(p, _kdtree_candidates, nearest_node_indices);
201 :
202 : // Collect all elements connected to these nodes
203 : std::set<const Elem *> candidate_elements;
204 239580 : for (auto node_idx : nearest_node_indices)
205 : {
206 : // Get the actual node from the mesh using the index
207 199650 : if (node_idx < _mesh.n_nodes())
208 : {
209 199650 : const Node * node = _mesh.node_ptr(node_idx);
210 199650 : dof_id_type node_id = node->id();
211 : auto it = _node_to_elements.find(node_id);
212 199650 : if (it != _node_to_elements.end())
213 : {
214 : const auto & connected_elems = it->second;
215 199650 : candidate_elements.insert(connected_elems.begin(), connected_elems.end());
216 : }
217 : }
218 : }
219 :
220 : // Convert set to vector for consistent type
221 : std::vector<const Elem *> candidate_vector(candidate_elements.begin(),
222 39930 : candidate_elements.end());
223 39930 : return findClosestElement(candidate_vector);
224 79860 : }
225 : else
226 : {
227 : // Fallback to original O(n) method if KDTree not available
228 : std::vector<const Elem *> all_elements;
229 0 : for (const auto & elem : _mesh.element_ptr_range())
230 0 : all_elements.push_back(elem);
231 :
232 0 : return findClosestElement(all_elements);
233 0 : }
234 : }
235 :
236 : Point
237 11197791 : ParameterMesh::closestPoint(const Elem & elem, const Point & p) const
238 : {
239 : mooseAssert(!elem.contains_point(p),
240 : "Points inside of elements shouldn't need to find closestPoint.");
241 :
242 : // Lambda to find closest point from range without storing temporary vectors
243 3101880 : auto findClosest = [&p](auto range, auto point_func) -> Point
244 : {
245 : Real min_distance = std::numeric_limits<Real>::max();
246 3101880 : Point min_point = p;
247 :
248 13545735 : for (const auto & item : range)
249 : {
250 10443855 : Point candidate = point_func(item);
251 10443855 : Real distance = (candidate - p).norm();
252 10443855 : if (distance < min_distance)
253 : {
254 : min_distance = distance;
255 4158198 : min_point = candidate;
256 : }
257 : }
258 3101880 : return min_point;
259 11197791 : };
260 :
261 11197791 : switch (elem.type())
262 : {
263 : case EDGE2:
264 : {
265 8056704 : LineSegment ls(*(elem.node_ptr(0)), *(elem.node_ptr(1)));
266 8056704 : return ls.closest_point(p);
267 : }
268 :
269 : case TRI3:
270 : {
271 2126376 : Point a = *(elem.node_ptr(0));
272 2126376 : Point b = *(elem.node_ptr(1));
273 2126376 : Point c = *(elem.node_ptr(2));
274 2126376 : Plane pl(a, b, c);
275 2126376 : Point trial = pl.closest_point(p);
276 2126376 : if (elem.contains_point(trial))
277 28899 : return trial;
278 :
279 2097477 : return findClosest(make_range(elem.n_edges()),
280 6292431 : [&](dof_id_type i) { return closestPoint(*elem.build_edge_ptr(i), p); });
281 2126376 : }
282 : case QUAD4:
283 : {
284 448800 : Point a = *(elem.node_ptr(0));
285 448800 : Point b = *(elem.node_ptr(1));
286 448800 : Point c = *(elem.node_ptr(2));
287 448800 : Point d = *(elem.node_ptr(3));
288 448800 : Plane pl1(a, b, c);
289 448800 : Plane pl2(b, c, d);
290 448800 : Point trial1 = pl1.closest_point(p);
291 448800 : Point trial2 = pl2.closest_point(p);
292 448800 : if (!trial1.absolute_fuzzy_equals(trial2, TOLERANCE * TOLERANCE))
293 0 : mooseError("Quad4 element is not coplanar");
294 :
295 448800 : if (elem.contains_point(trial1))
296 10308 : return trial1;
297 :
298 438492 : return findClosest(make_range(elem.n_edges()),
299 1753968 : [&](dof_id_type i) { return closestPoint(*elem.build_edge_ptr(i), p); });
300 448800 : }
301 :
302 565911 : default:
303 : {
304 565911 : if (elem.dim() == 3)
305 : {
306 565911 : return findClosest(make_range(elem.n_sides()),
307 2397456 : [&](dof_id_type i) { return closestPoint(*elem.build_side_ptr(i), p); });
308 : }
309 : else
310 : {
311 0 : mooseError("Unsupported element type ",
312 0 : Utility::enum_to_string(elem.type()),
313 : " for projection of parameter mesh.");
314 : }
315 : }
316 : }
317 : }
318 :
319 : template <typename T>
320 : T
321 3204 : ParameterMesh::computeRegularizationLoop(const std::vector<Real> & parameter_values,
322 : RegularizationType reg_type) const
323 : {
324 3204 : if (parameter_values.size() != _param_dofs)
325 0 : mooseError("Parameter values size (",
326 0 : parameter_values.size(),
327 : ") does not match mesh DOFs (",
328 0 : _param_dofs,
329 : ")");
330 :
331 : T result;
332 : if constexpr (std::is_same_v<T, Real>)
333 : result = 0.0;
334 : else if constexpr (std::is_same_v<T, std::vector<Real>>)
335 1602 : result.resize(_param_dofs, 0.0);
336 :
337 : // Iterate over all elements in the mesh
338 1355292 : for (const auto & elem : _mesh.element_ptr_range())
339 : {
340 : // Get DOF indices for this element
341 : std::vector<dof_id_type> dof_indices;
342 448560 : _dof_map->dof_indices(elem, dof_indices, _param_var_id);
343 :
344 : // Get quadrature rule for this element
345 448560 : const unsigned int dim = elem->dim();
346 448560 : QGauss qrule(dim, _fe_type.default_quadrature_order());
347 :
348 : // Create finite element objects
349 448560 : std::unique_ptr<FEBase> fe(FEBase::build(dim, _fe_type));
350 448560 : fe->attach_quadrature_rule(&qrule);
351 :
352 : // Request shape functions and derivatives before reinit
353 : const std::vector<Real> & JxW = fe->get_JxW();
354 : const std::vector<std::vector<Real>> & phi = fe->get_phi();
355 : const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
356 :
357 : // Reinitialize for current element
358 448560 : fe->reinit(elem);
359 :
360 2242800 : for (const auto qp : make_range(qrule.n_points()))
361 : {
362 : if constexpr (std::is_same_v<T, Real>)
363 897120 : result +=
364 897120 : computeRegularizationQp(parameter_values, phi, dphi, qp, dof_indices, JxW, reg_type);
365 : else if constexpr (std::is_same_v<T, std::vector<Real>>)
366 897120 : computeRegularizationGradientQp(
367 : parameter_values, phi, dphi, qp, dof_indices, JxW, reg_type, result);
368 : }
369 : }
370 :
371 3204 : return result;
372 0 : }
373 :
374 : Real
375 1602 : ParameterMesh::computeRegularizationObjective(const std::vector<Real> & parameter_values,
376 : RegularizationType reg_type) const
377 : {
378 1602 : return computeRegularizationLoop<Real>(parameter_values, reg_type);
379 : }
380 :
381 : std::vector<Real>
382 1602 : ParameterMesh::computeRegularizationGradient(const std::vector<Real> & parameter_values,
383 : RegularizationType reg_type) const
384 : {
385 1602 : return computeRegularizationLoop<std::vector<Real>>(parameter_values, reg_type);
386 : }
387 :
388 : Real
389 897120 : ParameterMesh::computeRegularizationQp(const std::vector<Real> & parameter_values,
390 : const std::vector<std::vector<Real>> & /*phi*/,
391 : const std::vector<std::vector<RealGradient>> & dphi,
392 : const unsigned int qp,
393 : const std::vector<dof_id_type> & dof_indices,
394 : const std::vector<Real> & JxW,
395 : RegularizationType reg_type) const
396 : {
397 : Real objective_contribution = 0.0;
398 :
399 : // Switch on regularization type
400 897120 : switch (reg_type)
401 : {
402 : case RegularizationType::L2_GRADIENT:
403 : {
404 : // Compute parameter gradient at this quadrature point
405 : RealGradient param_grad;
406 4485600 : for (const auto i : index_range(dof_indices))
407 3588480 : param_grad += parameter_values[dof_indices[i]] * dphi[i][qp];
408 :
409 : // Add L2 norm squared of gradient for regularization
410 897120 : objective_contribution = param_grad.norm_sq() * JxW[qp];
411 : break;
412 : }
413 0 : default:
414 0 : mooseError("Unknown Regularization Type");
415 : }
416 :
417 897120 : return objective_contribution;
418 : }
419 :
420 : void
421 897120 : ParameterMesh::computeRegularizationGradientQp(const std::vector<Real> & parameter_values,
422 : const std::vector<std::vector<Real>> & /*phi*/,
423 : const std::vector<std::vector<RealGradient>> & dphi,
424 : const unsigned int qp,
425 : const std::vector<dof_id_type> & dof_indices,
426 : const std::vector<Real> & JxW,
427 : RegularizationType reg_type,
428 : std::vector<Real> & gradient) const
429 : {
430 : // Switch on regularization type
431 897120 : switch (reg_type)
432 : {
433 : case RegularizationType::L2_GRADIENT:
434 : {
435 : // Compute parameter gradient at this quadrature point
436 : RealGradient param_grad;
437 4485600 : for (const auto i : index_range(dof_indices))
438 3588480 : param_grad += parameter_values[dof_indices[i]] * dphi[i][qp];
439 :
440 : // Compute gradient contribution: 2 * grad(p) * dphi_j
441 4485600 : for (const auto j : index_range(dof_indices))
442 3588480 : gradient[dof_indices[j]] += 2.0 * param_grad * dphi[j][qp] * JxW[qp];
443 : break;
444 : }
445 :
446 0 : default:
447 0 : mooseError("Unknown Regularization Type");
448 : }
449 897120 : }
|
