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 334 : ParameterMesh::ParameterMesh(const FEType & param_type,
41 : const std::string & exodus_mesh,
42 : const bool find_closest,
43 334 : const unsigned int kdtree_candidates)
44 334 : : _communicator(MPI_COMM_SELF),
45 334 : _mesh(_communicator),
46 334 : _find_closest(find_closest),
47 334 : _kdtree_candidates(kdtree_candidates),
48 334 : _param_var_id(0),
49 334 : _dof_map(nullptr),
50 334 : _fe_type(param_type)
51 : {
52 : _mesh.allow_renumbering(false);
53 334 : _mesh.prepare_for_use();
54 668 : _exodusII_io = std::make_unique<ExodusII_IO>(_mesh);
55 334 : _exodusII_io->read(exodus_mesh);
56 334 : _mesh.read(exodus_mesh);
57 : // Create system to store parameter values
58 668 : _eq = std::make_unique<EquationSystems>(_mesh);
59 334 : _sys = &_eq->add_system<ExplicitSystem>("_parameter_mesh_sys");
60 334 : _sys->add_variable("_parameter_mesh_var", param_type);
61 :
62 : // Create point locator
63 668 : _point_locator = PointLocatorBase::build(TREE_LOCAL_ELEMENTS, _mesh);
64 334 : _point_locator->enable_out_of_mesh_mode();
65 :
66 : // Initialize the equations systems
67 334 : _eq->init();
68 :
69 : // getting number of parameter dofs for size() function
70 334 : const unsigned short int var_id = _sys->variable_number("_parameter_mesh_var");
71 : std::set<dof_id_type> var_indices;
72 334 : _sys->local_dof_indices(var_id, var_indices);
73 334 : _param_dofs = var_indices.size();
74 :
75 334 : if (_find_closest)
76 : {
77 68096 : for (const auto & elem : _mesh.element_ptr_range())
78 34002 : if (elem->default_order() != FIRST)
79 48 : mooseError("Closet point projection currently does not support second order elements.");
80 : }
81 :
82 : // Initialize node-based KDTree optimization
83 332 : _mesh_nodes.clear();
84 : _node_to_elements.clear();
85 :
86 : // Extract all node coordinates
87 44228 : for (const auto & node : _mesh.node_ptr_range())
88 22114 : _mesh_nodes.push_back(*node);
89 :
90 : // Build node-to-elements connectivity map
91 78204 : for (const auto & elem : _mesh.element_ptr_range())
92 : {
93 237450 : for (const auto n : make_range(elem->n_nodes()))
94 : {
95 159910 : dof_id_type node_id = elem->node_id(n);
96 159910 : _node_to_elements[node_id].insert(elem);
97 : }
98 332 : }
99 :
100 : // Create KDTree from node coordinates
101 332 : if (!_mesh_nodes.empty())
102 664 : _node_kdtree = std::make_unique<KDTree>(_mesh_nodes, 10);
103 : // Update cached values for gradient computations
104 332 : const_cast<unsigned short int &>(_param_var_id) = var_id;
105 332 : const_cast<const DofMap *&>(_dof_map) = &_sys->get_dof_map();
106 332 : const_cast<FEType &>(_fe_type) = _dof_map->variable_type(_param_var_id);
107 332 : }
108 :
109 : void
110 8652643 : ParameterMesh::getIndexAndWeight(const Point & pt,
111 : std::vector<dof_id_type> & dof_indices,
112 : std::vector<Real> & weights) const
113 : {
114 8652643 : Point test_point = (_find_closest ? projectToMesh(pt) : pt);
115 :
116 8652643 : const Elem * elem = (*_point_locator)(test_point);
117 8652643 : 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 8652643 : _dof_map->dof_indices(elem, dof_indices, _param_var_id);
124 :
125 : // Map the physical co-ordinates to the reference co-ordinates
126 8652643 : Point coor = FEMap::inverse_map(elem->dim(), elem, test_point);
127 : // get the shape function value via the FEInterface
128 8652643 : FEComputeData fe_data(*_eq, coor);
129 8652643 : FEInterface::compute_data(elem->dim(), _fe_type, elem, fe_data);
130 : // Set weights to the value of the shape functions
131 8652643 : weights = fe_data.shape;
132 :
133 8652643 : if (dof_indices.size() != weights.size())
134 0 : mooseError("Internal error: weights and DoF indices do not have the same size.");
135 8652643 : }
136 :
137 : void
138 10800 : ParameterMesh::getIndexAndWeight(const Point & pt,
139 : std::vector<dof_id_type> & dof_indices,
140 : std::vector<RealGradient> & weights) const
141 : {
142 10800 : if (!_sys->has_variable("_parameter_mesh_var"))
143 0 : mooseError("Internal error: System being read does not contain _parameter_mesh_var.");
144 10800 : Point test_point = (_find_closest ? projectToMesh(pt) : pt);
145 : // Locate the element the point is in
146 10800 : 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 10800 : _dof_map->dof_indices(elem, dof_indices, _param_var_id);
152 :
153 : // Map the physical co-ordinates to the reference co-ordinates
154 10800 : Point coor = FEMap::inverse_map(elem->dim(), elem, test_point);
155 : // get the shape function value via the FEInterface
156 10800 : FEComputeData fe_data(*_eq, coor);
157 10800 : fe_data.enable_derivative();
158 10800 : FEInterface::compute_data(elem->dim(), _fe_type, elem, fe_data);
159 : // Set weights to the value of the shape functions
160 10800 : weights = fe_data.dshape;
161 :
162 10800 : if (dof_indices.size() != weights.size())
163 0 : mooseError("Internal error: weights and DoF indices do not have the same size.");
164 10800 : }
165 :
166 : Point
167 30482 : ParameterMesh::projectToMesh(const Point & p) const
168 : {
169 : // quick path: p already inside an element
170 30482 : if ((*_point_locator)(p))
171 1320 : return p;
172 :
173 : // Lambda to find closest point from elements using squared distance for efficiency
174 29162 : auto findClosestElement = [&p, this](const auto & elements) -> Point
175 : {
176 : Real best_d2 = std::numeric_limits<Real>::max();
177 29162 : Point best_point = p;
178 :
179 550831 : for (const auto * elem : elements)
180 : {
181 521669 : Point trial = closestPoint(*elem, p);
182 521669 : Real d2 = (trial - p).norm_sq();
183 521669 : if (d2 < best_d2)
184 : {
185 : best_d2 = d2;
186 106126 : best_point = trial;
187 : }
188 : }
189 :
190 29162 : if (best_d2 == std::numeric_limits<Real>::max())
191 0 : mooseError("project_to_mesh failed - no candidate elements.");
192 29162 : return best_point;
193 29162 : };
194 :
195 : // Use KDTree optimization if available
196 29162 : if (_node_kdtree && !_mesh_nodes.empty())
197 : {
198 : // Find K nearest nodes using KDTree
199 : std::vector<std::size_t> nearest_node_indices;
200 29162 : _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 174972 : for (auto node_idx : nearest_node_indices)
205 : {
206 : // Get the actual node from the mesh using the index
207 145810 : if (node_idx < _mesh.n_nodes())
208 : {
209 145810 : const Node * node = _mesh.node_ptr(node_idx);
210 145810 : dof_id_type node_id = node->id();
211 : auto it = _node_to_elements.find(node_id);
212 145810 : if (it != _node_to_elements.end())
213 : {
214 : const auto & connected_elems = it->second;
215 145810 : 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 29162 : candidate_elements.end());
223 29162 : return findClosestElement(candidate_vector);
224 58324 : }
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 7843473 : 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 2152153 : auto findClosest = [&p](auto range, auto point_func) -> Point
244 : {
245 : Real min_distance = std::numeric_limits<Real>::max();
246 2152153 : Point min_point = p;
247 :
248 9473957 : for (const auto & item : range)
249 : {
250 7321804 : Point candidate = point_func(item);
251 7321804 : Real distance = (candidate - p).norm();
252 7321804 : if (distance < min_distance)
253 : {
254 : min_distance = distance;
255 2932011 : min_point = candidate;
256 : }
257 : }
258 2152153 : return min_point;
259 7843473 : };
260 :
261 7843473 : switch (elem.type())
262 : {
263 : case EDGE2:
264 : {
265 5663464 : LineSegment ls(*(elem.node_ptr(0)), *(elem.node_ptr(1)));
266 5663464 : return ls.closest_point(p);
267 : }
268 :
269 : case TRI3:
270 : {
271 1417584 : Point a = *(elem.node_ptr(0));
272 1417584 : Point b = *(elem.node_ptr(1));
273 1417584 : Point c = *(elem.node_ptr(2));
274 1417584 : Plane pl(a, b, c);
275 1417584 : Point trial = pl.closest_point(p);
276 1417584 : if (elem.contains_point(trial))
277 19266 : return trial;
278 :
279 1398318 : return findClosest(make_range(elem.n_edges()),
280 4194954 : [&](dof_id_type i) { return closestPoint(*elem.build_edge_ptr(i), p); });
281 1417584 : }
282 : case QUAD4:
283 : {
284 374000 : Point a = *(elem.node_ptr(0));
285 374000 : Point b = *(elem.node_ptr(1));
286 374000 : Point c = *(elem.node_ptr(2));
287 374000 : Point d = *(elem.node_ptr(3));
288 374000 : Plane pl1(a, b, c);
289 374000 : Plane pl2(b, c, d);
290 374000 : Point trial1 = pl1.closest_point(p);
291 374000 : Point trial2 = pl2.closest_point(p);
292 374000 : if (!trial1.absolute_fuzzy_equals(trial2, TOLERANCE * TOLERANCE))
293 0 : mooseError("Quad4 element is not coplanar");
294 :
295 374000 : if (elem.contains_point(trial1))
296 8590 : return trial1;
297 :
298 365410 : return findClosest(make_range(elem.n_edges()),
299 1461640 : [&](dof_id_type i) { return closestPoint(*elem.build_edge_ptr(i), p); });
300 374000 : }
301 :
302 388425 : default:
303 : {
304 388425 : if (elem.dim() == 3)
305 : {
306 388425 : return findClosest(make_range(elem.n_sides()),
307 1665210 : [&](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 2136 : ParameterMesh::computeRegularizationLoop(const std::vector<Real> & parameter_values,
322 : RegularizationType reg_type) const
323 : {
324 2136 : 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 1068 : result.resize(_param_dofs, 0.0);
336 :
337 : // Iterate over all elements in the mesh
338 903528 : 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 299040 : _dof_map->dof_indices(elem, dof_indices, _param_var_id);
343 :
344 : // Get quadrature rule for this element
345 299040 : const unsigned int dim = elem->dim();
346 299040 : QGauss qrule(dim, _fe_type.default_quadrature_order());
347 :
348 : // Create finite element objects
349 299040 : std::unique_ptr<FEBase> fe(FEBase::build(dim, _fe_type));
350 299040 : 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 299040 : fe->reinit(elem);
359 :
360 1495200 : for (const auto qp : make_range(qrule.n_points()))
361 : {
362 : if constexpr (std::is_same_v<T, Real>)
363 598080 : result +=
364 598080 : computeRegularizationQp(parameter_values, phi, dphi, qp, dof_indices, JxW, reg_type);
365 : else if constexpr (std::is_same_v<T, std::vector<Real>>)
366 598080 : computeRegularizationGradientQp(
367 : parameter_values, phi, dphi, qp, dof_indices, JxW, reg_type, result);
368 : }
369 : }
370 :
371 2136 : return result;
372 0 : }
373 :
374 : Real
375 1068 : ParameterMesh::computeRegularizationObjective(const std::vector<Real> & parameter_values,
376 : RegularizationType reg_type) const
377 : {
378 1068 : return computeRegularizationLoop<Real>(parameter_values, reg_type);
379 : }
380 :
381 : std::vector<Real>
382 1068 : ParameterMesh::computeRegularizationGradient(const std::vector<Real> & parameter_values,
383 : RegularizationType reg_type) const
384 : {
385 1068 : return computeRegularizationLoop<std::vector<Real>>(parameter_values, reg_type);
386 : }
387 :
388 : Real
389 598080 : 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 598080 : switch (reg_type)
401 : {
402 : case RegularizationType::L2_GRADIENT:
403 : {
404 : // Compute parameter gradient at this quadrature point
405 : RealGradient param_grad;
406 2990400 : for (const auto i : index_range(dof_indices))
407 2392320 : param_grad += parameter_values[dof_indices[i]] * dphi[i][qp];
408 :
409 : // Add L2 norm squared of gradient for regularization
410 598080 : objective_contribution = param_grad.norm_sq() * JxW[qp];
411 : break;
412 : }
413 0 : default:
414 0 : mooseError("Unknown Regularization Type");
415 : }
416 :
417 598080 : return objective_contribution;
418 : }
419 :
420 : void
421 598080 : 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 598080 : switch (reg_type)
432 : {
433 : case RegularizationType::L2_GRADIENT:
434 : {
435 : // Compute parameter gradient at this quadrature point
436 : RealGradient param_grad;
437 2990400 : for (const auto i : index_range(dof_indices))
438 2392320 : param_grad += parameter_values[dof_indices[i]] * dphi[i][qp];
439 :
440 : // Compute gradient contribution: 2 * grad(p) * dphi_j
441 2990400 : for (const auto j : index_range(dof_indices))
442 2392320 : 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 598080 : }
|
