Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 :
20 : // Local Includes
21 : #include "libmesh/mesh_function.h"
22 : #include "libmesh/dense_vector.h"
23 : #include "libmesh/equation_systems.h"
24 : #include "libmesh/numeric_vector.h"
25 : #include "libmesh/dof_map.h"
26 : #include "libmesh/point_locator_tree.h"
27 : #include "libmesh/fe_base.h"
28 : #include "libmesh/fe_interface.h"
29 : #include "libmesh/fe_compute_data.h"
30 : #include "libmesh/mesh_base.h"
31 : #include "libmesh/point.h"
32 : #include "libmesh/elem.h"
33 : #include "libmesh/int_range.h"
34 : #include "libmesh/fe_map.h"
35 :
36 : namespace libMesh
37 : {
38 :
39 :
40 : //------------------------------------------------------------------
41 : // MeshFunction methods
42 66 : MeshFunction::MeshFunction (const EquationSystems & eqn_systems,
43 : const NumericVector<Number> & vec,
44 : const DofMap & dof_map,
45 : std::vector<unsigned int> vars,
46 66 : const FunctionBase<Number> * master) :
47 : FunctionBase<Number> (master),
48 : ParallelObject (eqn_systems),
49 26 : _eqn_systems (eqn_systems),
50 26 : _vector (vec),
51 26 : _dof_map (dof_map),
52 20 : _system_vars (std::move(vars)),
53 86 : _out_of_mesh_mode (false)
54 : {
55 66 : }
56 :
57 :
58 :
59 21 : MeshFunction::MeshFunction (const EquationSystems & eqn_systems,
60 : const NumericVector<Number> & vec,
61 : const DofMap & dof_map,
62 : const unsigned int var,
63 21 : const FunctionBase<Number> * master) :
64 : FunctionBase<Number> (master),
65 : ParallelObject (eqn_systems),
66 9 : _eqn_systems (eqn_systems),
67 9 : _vector (vec),
68 9 : _dof_map (dof_map),
69 18 : _system_vars (1,var),
70 21 : _out_of_mesh_mode (false)
71 : {
72 21 : }
73 :
74 56 : MeshFunction::MeshFunction (const MeshFunction & mf):
75 56 : FunctionBase<Number> (mf._master),
76 56 : ParallelObject (mf._eqn_systems),
77 40 : _eqn_systems (mf._eqn_systems),
78 56 : _vector (mf._vector),
79 56 : _dof_map (mf._dof_map),
80 56 : _system_vars (mf._system_vars),
81 56 : _out_of_mesh_mode (mf._out_of_mesh_mode)
82 : {
83 : // Initialize the mf and set the point locator if the
84 : // input mf had done so.
85 56 : if(mf.initialized())
86 : {
87 49 : this->MeshFunction::init();
88 :
89 49 : if(mf.get_point_locator().initialized())
90 49 : this->set_point_locator_tolerance(mf.get_point_locator().get_close_to_point_tol());
91 :
92 : }
93 :
94 56 : if (mf._subdomain_ids)
95 : _subdomain_ids =
96 : std::make_unique<std::set<subdomain_id_type>>
97 0 : (*mf._subdomain_ids);
98 56 : }
99 :
100 89 : MeshFunction::~MeshFunction () = default;
101 :
102 :
103 :
104 136 : void MeshFunction::init ()
105 : {
106 : // are indices of the desired variable(s) provided?
107 40 : libmesh_assert_greater (this->_system_vars.size(), 0);
108 :
109 : // Don't do twice...
110 136 : if (this->_initialized)
111 : {
112 0 : libmesh_assert(_point_locator);
113 0 : return;
114 : }
115 :
116 : // The Mesh owns the "master" PointLocator, while handing us a
117 : // PointLocator "proxy" that forwards all requests to the master.
118 136 : const MeshBase & mesh = this->_eqn_systems.get_mesh();
119 192 : _point_locator = mesh.sub_point_locator();
120 :
121 : // ready for use
122 136 : this->_initialized = true;
123 : }
124 :
125 :
126 :
127 : void
128 0 : MeshFunction::clear ()
129 : {
130 : // only delete the point locator when we are the master
131 0 : if (_point_locator && !_master)
132 0 : _point_locator.reset();
133 :
134 0 : this->_initialized = false;
135 0 : }
136 :
137 :
138 :
139 56 : std::unique_ptr<FunctionBase<Number>> MeshFunction::clone () const
140 : {
141 56 : return std::make_unique<MeshFunction>(*this);
142 : }
143 :
144 :
145 :
146 0 : Number MeshFunction::operator() (const Point & p,
147 : const Real time)
148 : {
149 0 : libmesh_assert (this->initialized());
150 :
151 0 : DenseVector<Number> buf (1);
152 0 : this->operator() (p, time, buf);
153 0 : return buf(0);
154 : }
155 :
156 21 : std::map<const Elem *, Number> MeshFunction::discontinuous_value (const Point & p,
157 : const Real time)
158 : {
159 6 : libmesh_assert (this->initialized());
160 :
161 12 : std::map<const Elem *, DenseVector<Number>> buffer;
162 21 : this->discontinuous_value (p, time, buffer);
163 6 : std::map<const Elem *, Number> return_value;
164 49 : for (const auto & [elem, vec] : buffer)
165 28 : return_value[elem] = vec(0);
166 : // NOTE: If no suitable element is found, then the map return_value is empty. This
167 : // puts burden on the user of this function but I don't really see a better way.
168 27 : return return_value;
169 : }
170 :
171 0 : Gradient MeshFunction::gradient (const Point & p,
172 : const Real time)
173 : {
174 0 : libmesh_assert (this->initialized());
175 :
176 0 : std::vector<Gradient> buf (1);
177 0 : this->gradient(p, time, buf);
178 0 : return buf[0];
179 : }
180 :
181 21 : std::map<const Elem *, Gradient> MeshFunction::discontinuous_gradient (const Point & p,
182 : const Real time)
183 : {
184 6 : libmesh_assert (this->initialized());
185 :
186 12 : std::map<const Elem *, std::vector<Gradient>> buffer;
187 21 : this->discontinuous_gradient (p, time, buffer);
188 6 : std::map<const Elem *, Gradient> return_value;
189 49 : for (const auto & [elem, vec] : buffer)
190 28 : return_value[elem] = vec[0];
191 : // NOTE: If no suitable element is found, then the map return_value is empty. This
192 : // puts burden on the user of this function but I don't really see a better way.
193 27 : return return_value;
194 : }
195 :
196 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
197 0 : Tensor MeshFunction::hessian (const Point & p,
198 : const Real time)
199 : {
200 0 : libmesh_assert (this->initialized());
201 :
202 0 : std::vector<Tensor> buf (1);
203 0 : this->hessian(p, time, buf);
204 0 : return buf[0];
205 : }
206 : #endif
207 :
208 75941 : void MeshFunction::operator() (const Point & p,
209 : const Real time,
210 : DenseVector<Number> & output)
211 : {
212 75941 : this->operator() (p, time, output, this->_subdomain_ids.get());
213 75941 : }
214 :
215 452829 : void MeshFunction::operator() (const Point & p,
216 : const Real,
217 : DenseVector<Number> & output,
218 : const std::set<subdomain_id_type> * subdomain_ids)
219 : {
220 113417 : libmesh_assert (this->initialized());
221 :
222 452829 : const Elem * element = this->find_element(p,subdomain_ids);
223 :
224 452829 : if (!element)
225 : {
226 : // We'd better be in out_of_mesh_mode if we couldn't find an
227 : // element in the mesh
228 80 : libmesh_assert (_out_of_mesh_mode);
229 80 : output = _out_of_mesh_value;
230 : }
231 : else
232 : {
233 : {
234 452509 : const unsigned int dim = element->dim();
235 :
236 :
237 : // Get local coordinates to feed these into compute_data().
238 : // Note that the fe_type can safely be used from the 0-variable,
239 : // since the inverse mapping is the same for all FEFamilies
240 225835 : const Point mapped_point (FEMap::inverse_map (dim, element,
241 226674 : p));
242 :
243 : // resize the output vector to the number of output values
244 : // that the user told us, this include multiple entries for
245 : // the spatial components of vector variable
246 113337 : unsigned int output_size = 0;
247 905042 : for (auto index : index_range(this->_system_vars))
248 : {
249 452533 : const unsigned int var = _system_vars[index];
250 :
251 452533 : if (var == libMesh::invalid_uint)
252 : {
253 0 : libmesh_assert (_out_of_mesh_mode &&
254 : index < _out_of_mesh_value.size());
255 0 : output(index) = _out_of_mesh_value(index);
256 0 : continue;
257 : }
258 :
259 452533 : const FEType & fe_type = this->_dof_map.variable_type(var);
260 :
261 : // Adding entries for each scalar variable and spatial components of
262 : // vector variables
263 452533 : switch (fe_type.family)
264 : {
265 24 : case HIERARCHIC_VEC:
266 : case L2_HIERARCHIC_VEC:
267 : case LAGRANGE_VEC:
268 : case L2_LAGRANGE_VEC:
269 : case MONOMIAL_VEC:
270 : case NEDELEC_ONE:
271 : case RAVIART_THOMAS:
272 : case L2_RAVIART_THOMAS:
273 : {
274 24 : output_size = output_size + dim;
275 24 : break;
276 : }
277 452509 : default:
278 : {
279 452509 : output_size++;
280 452509 : break;
281 : }
282 : }
283 :
284 : }
285 339172 : output.resize(output_size);
286 :
287 : // Adding a counter to keep the output indexing correct for mix scalar-vector output sets
288 113337 : unsigned int vec_count = 0;
289 :
290 : // loop over all vars
291 905042 : for (auto index : index_range(this->_system_vars))
292 : {
293 : // the data for this variable
294 452533 : const unsigned int var = _system_vars[index];
295 :
296 452533 : if (var == libMesh::invalid_uint)
297 : {
298 0 : libmesh_assert (_out_of_mesh_mode &&
299 : index < _out_of_mesh_value.size());
300 0 : output(index) = _out_of_mesh_value(index);
301 0 : continue;
302 : }
303 :
304 452533 : const FEType & fe_type = this->_dof_map.variable_type(var);
305 :
306 : // The method between vector and scalar variable are different:
307 : // For scalar variables, we use the previous established method of using FEInterface::compute_data
308 : // For vector variables, we call the phi() data directly to build the variable solution
309 452533 : switch (fe_type.family)
310 : {
311 24 : case HIERARCHIC_VEC:
312 : case L2_HIERARCHIC_VEC:
313 : case LAGRANGE_VEC:
314 : case L2_LAGRANGE_VEC:
315 : case MONOMIAL_VEC:
316 : case NEDELEC_ONE:
317 : case RAVIART_THOMAS:
318 : case L2_RAVIART_THOMAS:
319 : {
320 : // Calling the physical mesh point needed for the shape function
321 40 : FEComputeData data (this->_eqn_systems, mapped_point);
322 24 : const Point & pt = data.p;
323 :
324 : // where the solution values for the var-th variable are stored
325 16 : std::vector<dof_id_type> dof_indices;
326 24 : this->_dof_map.dof_indices (element, dof_indices, var);
327 :
328 : // Calling the properly-transformed shape function using get_phi()
329 32 : std::unique_ptr<FEVectorBase> vec_fe = FEVectorBase::build(dim, fe_type);
330 32 : std::vector<Point> vec_pt = {pt};
331 8 : const std::vector<std::vector<RealGradient>> & vec_phi = vec_fe->get_phi();
332 24 : vec_fe->reinit(element, &vec_pt);
333 :
334 : // interpolate the solution
335 : {
336 72 : for (unsigned int d = 0; d < dim; d++)
337 : {
338 16 : Number value = 0.;
339 :
340 288 : for (auto i : index_range(dof_indices))
341 320 : value += this->_vector(dof_indices[i]) * (vec_phi[i][0](d));
342 :
343 64 : output(index+vec_count*(dim-1)+d) = value;
344 : }
345 : }
346 :
347 24 : vec_count++;
348 8 : break;
349 16 : }
350 452509 : default:
351 : {
352 : // Build an FEComputeData that contains both input and output data
353 : // for the specific compute_data method.
354 :
355 679183 : FEComputeData data (this->_eqn_systems, mapped_point);
356 :
357 452509 : FEInterface::compute_data (dim, fe_type, element, data);
358 :
359 : // where the solution values for the var-th variable are stored
360 226674 : std::vector<dof_id_type> dof_indices;
361 452509 : this->_dof_map.dof_indices (element, dof_indices, var);
362 :
363 : // interpolate the solution
364 : {
365 113337 : Number value = 0.;
366 :
367 4976617 : for (auto i : index_range(dof_indices))
368 5659292 : value += this->_vector(dof_indices[i]) * data.shape[i];
369 :
370 565846 : output(index+vec_count*(dim-1)) = value;
371 : }
372 113337 : break;
373 225835 : }
374 : }
375 :
376 : // next variable
377 : }
378 : }
379 : }
380 452829 : }
381 :
382 21 : void MeshFunction::discontinuous_value (const Point & p,
383 : const Real time,
384 : std::map<const Elem *, DenseVector<Number>> & output)
385 : {
386 21 : this->discontinuous_value (p, time, output, this->_subdomain_ids.get());
387 21 : }
388 :
389 :
390 :
391 21 : void MeshFunction::discontinuous_value (const Point & p,
392 : const Real,
393 : std::map<const Elem *, DenseVector<Number>> & output,
394 : const std::set<subdomain_id_type> * subdomain_ids)
395 : {
396 6 : libmesh_assert (this->initialized());
397 :
398 : // clear the output map
399 6 : output.clear();
400 :
401 : // get the candidate elements
402 27 : std::set<const Elem *> candidate_element = this->find_elements(p,subdomain_ids);
403 :
404 : // loop through all candidates, if the set is empty this function will return an
405 : // empty map
406 49 : for (const auto & element : candidate_element)
407 : {
408 28 : if (!element)
409 : {
410 : // We'd better be in out_of_mesh_mode if we couldn't find an
411 : // element in the mesh
412 0 : libmesh_assert (_out_of_mesh_mode);
413 0 : output[element] = _out_of_mesh_value;
414 0 : continue;
415 : }
416 :
417 28 : const unsigned int dim = element->dim();
418 :
419 : // define a temporary vector to store all values
420 36 : DenseVector<Number> temp_output (cast_int<unsigned int>(this->_system_vars.size()));
421 :
422 : // Get local coordinates to feed these into compute_data().
423 : // Note that the fe_type can safely be used from the 0-variable,
424 : // since the inverse mapping is the same for all FEFamilies
425 28 : const Point mapped_point (FEMap::inverse_map (dim, element, p));
426 :
427 : // loop over all vars
428 56 : for (auto index : index_range(this->_system_vars))
429 : {
430 : // the data for this variable
431 28 : const unsigned int var = _system_vars[index];
432 :
433 28 : if (var == libMesh::invalid_uint)
434 : {
435 0 : libmesh_assert (_out_of_mesh_mode &&
436 : index < _out_of_mesh_value.size());
437 0 : temp_output(index) = _out_of_mesh_value(index);
438 0 : continue;
439 : }
440 :
441 28 : const FEType & fe_type = this->_dof_map.variable_type(var);
442 :
443 : // Build an FEComputeData that contains both input and output data
444 : // for the specific compute_data method.
445 : {
446 44 : FEComputeData data (this->_eqn_systems, mapped_point);
447 :
448 28 : FEInterface::compute_data (dim, fe_type, element, data);
449 :
450 : // where the solution values for the var-th variable are stored
451 8 : std::vector<dof_id_type> dof_indices;
452 28 : this->_dof_map.dof_indices (element, dof_indices, var);
453 :
454 : // interpolate the solution
455 : {
456 8 : Number value = 0.;
457 :
458 56 : for (auto i : index_range(dof_indices))
459 36 : value += this->_vector(dof_indices[i]) * data.shape[i];
460 :
461 28 : temp_output(index) = value;
462 : }
463 :
464 12 : }
465 :
466 : // next variable
467 : }
468 :
469 : // Insert temp_output into output
470 28 : output[element] = temp_output;
471 : }
472 21 : }
473 :
474 :
475 :
476 11 : void MeshFunction::gradient (const Point & p,
477 : const Real time,
478 : std::vector<Gradient> & output)
479 : {
480 11 : this->gradient(p, time, output, this->_subdomain_ids.get());
481 11 : }
482 :
483 :
484 :
485 375851 : void MeshFunction::gradient (const Point & p,
486 : const Real,
487 : std::vector<Gradient> & output,
488 : const std::set<subdomain_id_type> * subdomain_ids)
489 : {
490 93963 : libmesh_assert (this->initialized());
491 :
492 375851 : const Elem * element = this->find_element(p,subdomain_ids);
493 :
494 375851 : this->_gradient_on_elem(p, element, output);
495 375851 : }
496 :
497 :
498 21 : void MeshFunction::discontinuous_gradient (const Point & p,
499 : const Real time,
500 : std::map<const Elem *, std::vector<Gradient>> & output)
501 : {
502 21 : this->discontinuous_gradient (p, time, output, this->_subdomain_ids.get());
503 21 : }
504 :
505 :
506 :
507 21 : void MeshFunction::discontinuous_gradient (const Point & p,
508 : const Real,
509 : std::map<const Elem *, std::vector<Gradient>> & output,
510 : const std::set<subdomain_id_type> * subdomain_ids)
511 : {
512 6 : libmesh_assert (this->initialized());
513 :
514 : // clear the output map
515 6 : output.clear();
516 :
517 : // get the candidate elements
518 27 : std::set<const Elem *> candidate_element = this->find_elements(p,subdomain_ids);
519 :
520 : // loop through all candidates, if the set is empty this function will return an
521 : // empty map
522 49 : for (const auto & element : candidate_element)
523 : {
524 : // define a temporary vector to store all values
525 44 : std::vector<Gradient> temp_output (cast_int<unsigned int>(this->_system_vars.size()));
526 :
527 28 : this->_gradient_on_elem(p, element, temp_output);
528 :
529 : // Insert temp_output into output
530 20 : output.emplace(element, std::move(temp_output));
531 : }
532 21 : }
533 :
534 :
535 :
536 375879 : void MeshFunction::_gradient_on_elem (const Point & p,
537 : const Elem * element,
538 : std::vector<Gradient> & output)
539 : {
540 : // resize the output vector to the number of output values
541 : // that the user told us
542 469850 : output.resize (this->_system_vars.size());
543 :
544 375879 : if (!element)
545 : {
546 4 : libmesh_error_msg_if(!_out_of_mesh_mode,
547 : "MeshFunction couldn't find element at p=" <<
548 : p << ", but is not in out-of-mesh mode");
549 1 : libmesh_assert_equal_to (_out_of_mesh_value.size(),
550 : this->_system_vars.size());
551 16 : for (auto i : index_range(this->_system_vars))
552 15 : output[i] = Gradient(_out_of_mesh_value[i]);
553 4 : return;
554 : }
555 :
556 375875 : const unsigned int dim = element->dim();
557 :
558 : // Get local coordinates to feed these into compute_data().
559 : // Note that the fe_type can safely be used from the 0-variable,
560 : // since the inverse mapping is the same for all FEFamilies
561 187935 : const Point mapped_point (FEMap::inverse_map (dim, element,
562 187940 : p));
563 :
564 469845 : std::vector<Point> point_list (1, mapped_point);
565 :
566 : // loop over all vars
567 751764 : for (auto index : index_range(this->_system_vars))
568 : {
569 : // the data for this variable
570 375889 : const unsigned int var = _system_vars[index];
571 :
572 375889 : if (var == libMesh::invalid_uint)
573 : {
574 2 : libmesh_assert (_out_of_mesh_mode &&
575 : index < _out_of_mesh_value.size());
576 7 : output[index] = Gradient(_out_of_mesh_value(index));
577 7 : continue;
578 : }
579 :
580 375882 : const FEType & fe_type = this->_dof_map.variable_type(var);
581 :
582 : // where the solution values for the var-th variable are stored
583 93972 : std::vector<dof_id_type> dof_indices;
584 375882 : this->_dof_map.dof_indices (element, dof_indices, var);
585 :
586 : // interpolate the solution
587 93972 : Gradient grad(0.);
588 :
589 : // for performance-reasons, we use different algorithms now.
590 : // TODO: Check that both give the same result for finite elements.
591 : // Otherwive it is wrong...
592 281910 : if (!element->infinite())
593 : {
594 469854 : std::unique_ptr<FEBase> point_fe (FEBase::build(dim, fe_type));
595 93972 : const std::vector<std::vector<RealGradient>> & dphi = point_fe->get_dphi();
596 375882 : point_fe->reinit(element, &point_list);
597 :
598 3382749 : for (auto i : index_range(dof_indices))
599 3758589 : grad.add_scaled(dphi[i][0], this->_vector(dof_indices[i]));
600 :
601 187938 : }
602 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
603 : else
604 : {
605 : // Build an FEComputeData that contains both input and output data
606 : // for the specific compute_data method.
607 : //
608 : // TODO: enable this for a vector of points as well...
609 0 : FEComputeData data (this->_eqn_systems, mapped_point);
610 0 : data.enable_derivative();
611 0 : FEInterface::compute_data (dim, fe_type, element, data);
612 :
613 : // grad [x] = data.dshape[i](v) * dv/dx * dof_index [i]
614 : // sum over all indices
615 0 : for (auto i : index_range(dof_indices))
616 : {
617 : // local coordinates
618 0 : for (std::size_t v=0; v<dim; v++)
619 0 : for (std::size_t xyz=0; xyz<LIBMESH_DIM; xyz++)
620 : {
621 : // FIXME: this needs better syntax: It is matrix-vector multiplication.
622 0 : grad(xyz) += data.local_transform[v][xyz]
623 0 : * data.dshape[i](v)
624 0 : * this->_vector(dof_indices[i]);
625 : }
626 : }
627 0 : }
628 : #endif
629 469854 : output[index] = grad;
630 : }
631 : }
632 :
633 :
634 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
635 11 : void MeshFunction::hessian (const Point & p,
636 : const Real time,
637 : std::vector<Tensor> & output)
638 : {
639 11 : this->hessian(p, time, output, this->_subdomain_ids.get());
640 11 : }
641 :
642 :
643 :
644 375851 : void MeshFunction::hessian (const Point & p,
645 : const Real,
646 : std::vector<Tensor> & output,
647 : const std::set<subdomain_id_type> * subdomain_ids)
648 : {
649 93963 : libmesh_assert (this->initialized());
650 :
651 : // resize the output vector to the number of output values
652 : // that the user told us
653 469814 : output.resize (this->_system_vars.size());
654 :
655 375851 : const Elem * element = this->find_element(p,subdomain_ids);
656 :
657 375851 : if (!element)
658 : {
659 4 : libmesh_error_msg_if(!_out_of_mesh_mode,
660 : "MeshFunction couldn't find element at p=" <<
661 : p << ", but is not in out-of-mesh mode");
662 1 : libmesh_assert_equal_to (_out_of_mesh_value.size(),
663 : this->_system_vars.size());
664 16 : for (auto i : index_range(this->_system_vars))
665 15 : output[i] = Tensor(_out_of_mesh_value[i]);
666 1 : return;
667 : }
668 : else
669 : {
670 281885 : if(element->infinite())
671 : libmesh_warning("Warning: Requested the Hessian of an Infinite element."
672 : << "Second derivatives for Infinite elements"
673 : << " are not yet implemented!"
674 : << std::endl);
675 :
676 : {
677 375847 : const unsigned int dim = element->dim();
678 :
679 :
680 : // Get local coordinates to feed these into compute_data().
681 : // Note that the fe_type can safely be used from the 0-variable,
682 : // since the inverse mapping is the same for all FEFamilies
683 187923 : const Point mapped_point (FEMap::inverse_map (dim, element,
684 187924 : p));
685 :
686 469809 : std::vector<Point> point_list (1, mapped_point);
687 :
688 : // loop over all vars
689 751708 : for (auto index : index_range(this->_system_vars))
690 : {
691 : // the data for this variable
692 375861 : const unsigned int var = _system_vars[index];
693 :
694 375861 : if (var == libMesh::invalid_uint)
695 : {
696 2 : libmesh_assert (_out_of_mesh_mode &&
697 : index < _out_of_mesh_value.size());
698 7 : output[index] = Tensor(_out_of_mesh_value(index));
699 7 : continue;
700 : }
701 375854 : const FEType & fe_type = this->_dof_map.variable_type(var);
702 :
703 469818 : std::unique_ptr<FEBase> point_fe (FEBase::build(dim, fe_type));
704 : const std::vector<std::vector<RealTensor>> & d2phi =
705 93964 : point_fe->get_d2phi();
706 375854 : point_fe->reinit(element, &point_list);
707 :
708 : // where the solution values for the var-th variable are stored
709 187928 : std::vector<dof_id_type> dof_indices;
710 375854 : this->_dof_map.dof_indices (element, dof_indices, var);
711 :
712 : // interpolate the solution
713 93964 : Tensor hess;
714 :
715 3382644 : for (auto i : index_range(dof_indices))
716 3758490 : hess.add_scaled(d2phi[i][0], this->_vector(dof_indices[i]));
717 :
718 469818 : output[index] = hess;
719 187926 : }
720 : }
721 : }
722 : }
723 : #endif
724 :
725 1204531 : const Elem * MeshFunction::find_element(const Point & p,
726 : const std::set<subdomain_id_type> * subdomain_ids) const
727 : {
728 : // Ensure that in the case of a master mesh function, the
729 : // out-of-mesh mode is enabled either for both or for none. This is
730 : // important because the out-of-mesh mode is also communicated to
731 : // the point locator. Since this is time consuming, enable it only
732 : // in debug mode.
733 : #ifdef DEBUG
734 301343 : if (this->_master != nullptr)
735 : {
736 : const MeshFunction * master =
737 0 : cast_ptr<const MeshFunction *>(this->_master);
738 0 : libmesh_error_msg_if(_out_of_mesh_mode!=master->_out_of_mesh_mode,
739 : "ERROR: If you use out-of-mesh-mode in connection with master mesh "
740 : "functions, you must enable out-of-mesh mode for both the master and the slave mesh function.");
741 : }
742 : #endif
743 :
744 : // locate the point in the other mesh
745 1204531 : const Elem * element = (*_point_locator)(p, subdomain_ids);
746 :
747 : // Make sure that the element found is evaluable
748 1204531 : return this->check_found_elem(element, p);
749 : }
750 :
751 42 : std::set<const Elem *> MeshFunction::find_elements(const Point & p,
752 : const std::set<subdomain_id_type> * subdomain_ids) const
753 : {
754 : // Ensure that in the case of a master mesh function, the
755 : // out-of-mesh mode is enabled either for both or for none. This is
756 : // important because the out-of-mesh mode is also communicated to
757 : // the point locator. Since this is time consuming, enable it only
758 : // in debug mode.
759 : #ifdef DEBUG
760 12 : if (this->_master != nullptr)
761 : {
762 : const MeshFunction * master =
763 0 : cast_ptr<const MeshFunction *>(this->_master);
764 0 : libmesh_error_msg_if(_out_of_mesh_mode!=master->_out_of_mesh_mode,
765 : "ERROR: If you use out-of-mesh-mode in connection with master mesh "
766 : "functions, you must enable out-of-mesh mode for both the master and the slave mesh function.");
767 : }
768 : #endif
769 :
770 : // locate the point in the other mesh
771 24 : std::set<const Elem *> candidate_elements;
772 12 : std::set<const Elem *> final_candidate_elements;
773 42 : (*_point_locator)(p, candidate_elements, subdomain_ids);
774 :
775 : // For each candidate Elem, if it is evaluable, add it to the set of
776 : // final candidate Elems.
777 98 : for (const auto & element : candidate_elements)
778 56 : final_candidate_elements.insert(this->check_found_elem(element, p));
779 :
780 54 : return final_candidate_elements;
781 : }
782 :
783 : const Elem *
784 1204587 : MeshFunction::check_found_elem(const Elem * element, const Point & p) const
785 : {
786 : // If the PointLocator did not find an Element containing Point p,
787 : // OR if the PointLocator found a local element containting Point p,
788 : // we can simply return it.
789 1204587 : if (!element || (element->processor_id() == this->processor_id()))
790 292875 : return element;
791 :
792 : // Otherwise, the PointLocator returned a valid but non-local
793 : // element. Therefore, we have to do some further checks:
794 : //
795 : // 1.) If we have a SERIAL _vector, then we can just return the
796 : // non-local element, because all the DOFs are available.
797 25452 : if (_vector.type() == SERIAL)
798 8358 : return element;
799 :
800 : // 2.) If we have a GHOSTED _vector, then we can return a non-local
801 : // Element provided that all of its DOFs are ghosted on this
802 : // processor. That should be faster than the other option, which
803 : // is to search through all the Elem's point_neighbors() for a
804 : // suitable local Elem.
805 378 : if (_vector.type() == GHOSTED && _dof_map.is_evaluable(*element))
806 113 : return element;
807 :
808 : // 3.) If we have a PARALLEL vector, then we search the non-local
809 : // Elem's point neighbors for a local one to return instead. If
810 : // we don't eventually find one, we just return nullptr.
811 13 : std::set<const Elem *> point_neighbors;
812 39 : element->find_point_neighbors(p, point_neighbors);
813 13 : element = nullptr;
814 130 : for (const auto & elem : point_neighbors)
815 178 : if (elem->processor_id() == this->processor_id())
816 : {
817 13 : element = elem;
818 13 : break;
819 : }
820 :
821 13 : return element;
822 : }
823 :
824 98 : const PointLocatorBase & MeshFunction::get_point_locator () const
825 : {
826 28 : libmesh_assert (this->initialized());
827 98 : return *_point_locator;
828 : }
829 :
830 0 : PointLocatorBase & MeshFunction::get_point_locator ()
831 : {
832 0 : libmesh_assert (this->initialized());
833 0 : return *_point_locator;
834 : }
835 :
836 21 : void MeshFunction::enable_out_of_mesh_mode(const DenseVector<Number> & value)
837 : {
838 6 : libmesh_assert (this->initialized());
839 21 : _point_locator->enable_out_of_mesh_mode();
840 21 : _out_of_mesh_mode = true;
841 6 : _out_of_mesh_value = value;
842 21 : }
843 :
844 0 : void MeshFunction::enable_out_of_mesh_mode(const Number & value)
845 : {
846 0 : DenseVector<Number> v(1);
847 0 : v(0) = value;
848 0 : this->enable_out_of_mesh_mode(v);
849 0 : }
850 :
851 0 : void MeshFunction::disable_out_of_mesh_mode()
852 : {
853 0 : libmesh_assert (this->initialized());
854 0 : _point_locator->disable_out_of_mesh_mode();
855 0 : _out_of_mesh_mode = false;
856 0 : }
857 :
858 56 : void MeshFunction::set_point_locator_tolerance(Real tol)
859 : {
860 56 : _point_locator->set_close_to_point_tol(tol);
861 56 : _point_locator->set_contains_point_tol(tol);
862 56 : }
863 :
864 0 : void MeshFunction::unset_point_locator_tolerance()
865 : {
866 0 : _point_locator->unset_close_to_point_tol();
867 0 : }
868 :
869 7 : void MeshFunction::set_subdomain_ids(const std::set<subdomain_id_type> * subdomain_ids)
870 : {
871 7 : if (subdomain_ids)
872 12 : _subdomain_ids = std::make_unique<std::set<subdomain_id_type>>(*subdomain_ids);
873 : else
874 0 : _subdomain_ids.reset();
875 7 : }
876 :
877 : } // namespace libMesh
|
