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