libMesh
mesh_function.C
Go to the documentation of this file.
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
43  const NumericVector<Number> & vec,
44  const DofMap & dof_map,
45  std::vector<unsigned int> vars,
46  const FunctionBase<Number> * master) :
47  FunctionBase<Number> (master),
48  ParallelObject (eqn_systems),
49  _eqn_systems (eqn_systems),
50  _vector (vec),
51  _dof_map (dof_map),
52  _system_vars (std::move(vars)),
53  _out_of_mesh_mode (false)
54 {
55 }
56 
57 
58 
60  const NumericVector<Number> & vec,
61  const DofMap & dof_map,
62  const unsigned int var,
63  const FunctionBase<Number> * master) :
64  FunctionBase<Number> (master),
65  ParallelObject (eqn_systems),
66  _eqn_systems (eqn_systems),
67  _vector (vec),
68  _dof_map (dof_map),
69  _system_vars (1,var),
70  _out_of_mesh_mode (false)
71 {
72 }
73 
75  FunctionBase<Number> (mf._master),
76  ParallelObject (mf._eqn_systems),
77  _eqn_systems (mf._eqn_systems),
78  _vector (mf._vector),
79  _dof_map (mf._dof_map),
80  _system_vars (mf._system_vars),
81  _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  if(mf.initialized())
86  {
87  this->MeshFunction::init();
88 
91 
92  }
93 
94  if (mf._subdomain_ids)
96  std::make_unique<std::set<subdomain_id_type>>
97  (*mf._subdomain_ids);
98 }
99 
100 MeshFunction::~MeshFunction () = default;
101 
102 
103 
105 {
106  // are indices of the desired variable(s) provided?
107  libmesh_assert_greater (this->_system_vars.size(), 0);
108 
109  // Don't do twice...
110  if (this->_initialized)
111  {
113  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  const MeshBase & mesh = this->_eqn_systems.get_mesh();
120 
121  // ready for use
122  this->_initialized = true;
123 }
124 
125 
126 
127 #ifdef LIBMESH_ENABLE_DEPRECATED
128 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  this->init();
136 }
137 #endif // LIBMESH_ENABLE_DEPRECATED
138 
139 
140 
141 void
143 {
144  // only delete the point locator when we are the master
145  if (_point_locator && !_master)
146  _point_locator.reset();
147 
148  this->_initialized = false;
149 }
150 
151 
152 
153 std::unique_ptr<FunctionBase<Number>> MeshFunction::clone () const
154 {
155  return std::make_unique<MeshFunction>(*this);
156 }
157 
158 
159 
161  const Real time)
162 {
163  libmesh_assert (this->initialized());
164 
165  DenseVector<Number> buf (1);
166  this->operator() (p, time, buf);
167  return buf(0);
168 }
169 
170 std::map<const Elem *, Number> MeshFunction::discontinuous_value (const Point & p,
171  const Real time)
172 {
173  libmesh_assert (this->initialized());
174 
175  std::map<const Elem *, DenseVector<Number>> buffer;
176  this->discontinuous_value (p, time, buffer);
177  std::map<const Elem *, Number> return_value;
178  for (const auto & [elem, vec] : buffer)
179  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  return return_value;
183 }
184 
186  const Real time)
187 {
188  libmesh_assert (this->initialized());
189 
190  std::vector<Gradient> buf (1);
191  this->gradient(p, time, buf);
192  return buf[0];
193 }
194 
195 std::map<const Elem *, Gradient> MeshFunction::discontinuous_gradient (const Point & p,
196  const Real time)
197 {
198  libmesh_assert (this->initialized());
199 
200  std::map<const Elem *, std::vector<Gradient>> buffer;
201  this->discontinuous_gradient (p, time, buffer);
202  std::map<const Elem *, Gradient> return_value;
203  for (const auto & [elem, vec] : buffer)
204  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  return return_value;
208 }
209 
210 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
212  const Real time)
213 {
214  libmesh_assert (this->initialized());
215 
216  std::vector<Tensor> buf (1);
217  this->hessian(p, time, buf);
218  return buf[0];
219 }
220 #endif
221 
223  const Real time,
224  DenseVector<Number> & output)
225 {
226  this->operator() (p, time, output, this->_subdomain_ids.get());
227 }
228 
230  const Real,
231  DenseVector<Number> & output,
232  const std::set<subdomain_id_type> * subdomain_ids)
233 {
234  libmesh_assert (this->initialized());
235 
236  const Elem * element = this->find_element(p,subdomain_ids);
237 
238  if (!element)
239  {
240  // We'd better be in out_of_mesh_mode if we couldn't find an
241  // element in the mesh
243  output = _out_of_mesh_value;
244  }
245  else
246  {
247  {
248  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  const Point mapped_point (FEMap::inverse_map (dim, element,
255  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  unsigned int output_size = 0;
261  for (auto index : index_range(this->_system_vars))
262  {
263  const unsigned int var = _system_vars[index];
264 
265  if (var == libMesh::invalid_uint)
266  {
268  index < _out_of_mesh_value.size());
269  output(index) = _out_of_mesh_value(index);
270  continue;
271  }
272 
273  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  switch (fe_type.family)
278  {
279  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  output_size = output_size + dim;
289  break;
290  }
291  default:
292  {
293  output_size++;
294  break;
295  }
296  }
297 
298  }
299  output.resize(output_size);
300 
301  // Adding a counter to keep the output indexing correct for mix scalar-vector output sets
302  unsigned int vec_count = 0;
303 
304  // loop over all vars
305  for (auto index : index_range(this->_system_vars))
306  {
307  // the data for this variable
308  const unsigned int var = _system_vars[index];
309 
310  if (var == libMesh::invalid_uint)
311  {
313  index < _out_of_mesh_value.size());
314  output(index) = _out_of_mesh_value(index);
315  continue;
316  }
317 
318  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  switch (fe_type.family)
324  {
325  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  FEComputeData data (this->_eqn_systems, mapped_point);
336  const Point & pt = data.p;
337 
338  // where the solution values for the var-th variable are stored
339  std::vector<dof_id_type> dof_indices;
340  this->_dof_map.dof_indices (element, dof_indices, var);
341 
342  // Calling the properly-transformed shape function using get_phi()
343  std::unique_ptr<FEVectorBase> vec_fe = FEVectorBase::build(dim, fe_type);
344  std::vector<Point> vec_pt = {pt};
345  const std::vector<std::vector<RealGradient>> & vec_phi = vec_fe->get_phi();
346  vec_fe->reinit(element, &vec_pt);
347 
348  // interpolate the solution
349  {
350  for (unsigned int d = 0; d < dim; d++)
351  {
352  Number value = 0.;
353 
354  for (auto i : index_range(dof_indices))
355  value += this->_vector(dof_indices[i]) * (vec_phi[i][0](d));
356 
357  output(index+vec_count*(dim-1)+d) = value;
358  }
359  }
360 
361  vec_count++;
362  break;
363  }
364  default:
365  {
366  // Build an FEComputeData that contains both input and output data
367  // for the specific compute_data method.
368 
369  FEComputeData data (this->_eqn_systems, mapped_point);
370 
371  FEInterface::compute_data (dim, fe_type, element, data);
372 
373  // where the solution values for the var-th variable are stored
374  std::vector<dof_id_type> dof_indices;
375  this->_dof_map.dof_indices (element, dof_indices, var);
376 
377  // interpolate the solution
378  {
379  Number value = 0.;
380 
381  for (auto i : index_range(dof_indices))
382  value += this->_vector(dof_indices[i]) * data.shape[i];
383 
384  output(index+vec_count*(dim-1)) = value;
385  }
386  break;
387  }
388  }
389 
390  // next variable
391  }
392  }
393  }
394 }
395 
397  const Real time,
398  std::map<const Elem *, DenseVector<Number>> & output)
399 {
400  this->discontinuous_value (p, time, output, this->_subdomain_ids.get());
401 }
402 
403 
404 
406  const Real,
407  std::map<const Elem *, DenseVector<Number>> & output,
408  const std::set<subdomain_id_type> * subdomain_ids)
409 {
410  libmesh_assert (this->initialized());
411 
412  // clear the output map
413  output.clear();
414 
415  // get the candidate elements
416  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  for (const auto & element : candidate_element)
421  {
422  const unsigned int dim = element->dim();
423 
424  // define a temporary vector to store all values
425  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  const Point mapped_point (FEMap::inverse_map (dim, element, p));
431 
432  // loop over all vars
433  for (auto index : index_range(this->_system_vars))
434  {
435  // the data for this variable
436  const unsigned int var = _system_vars[index];
437 
438  if (var == libMesh::invalid_uint)
439  {
441  index < _out_of_mesh_value.size());
442  temp_output(index) = _out_of_mesh_value(index);
443  continue;
444  }
445 
446  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  FEComputeData data (this->_eqn_systems, mapped_point);
452 
453  FEInterface::compute_data (dim, fe_type, element, data);
454 
455  // where the solution values for the var-th variable are stored
456  std::vector<dof_id_type> dof_indices;
457  this->_dof_map.dof_indices (element, dof_indices, var);
458 
459  // interpolate the solution
460  {
461  Number value = 0.;
462 
463  for (auto i : index_range(dof_indices))
464  value += this->_vector(dof_indices[i]) * data.shape[i];
465 
466  temp_output(index) = value;
467  }
468 
469  }
470 
471  // next variable
472  }
473 
474  // Insert temp_output into output
475  output[element] = temp_output;
476  }
477 }
478 
479 
480 
482  const Real time,
483  std::vector<Gradient> & output)
484 {
485  this->gradient(p, time, output, this->_subdomain_ids.get());
486 }
487 
488 
489 
491  const Real,
492  std::vector<Gradient> & output,
493  const std::set<subdomain_id_type> * subdomain_ids)
494 {
495  libmesh_assert (this->initialized());
496 
497  const Elem * element = this->find_element(p,subdomain_ids);
498 
499  if (!element)
500  output.resize(0);
501  else
502  this->_gradient_on_elem(p, element, output);
503 }
504 
505 
507  const Real time,
508  std::map<const Elem *, std::vector<Gradient>> & output)
509 {
510  this->discontinuous_gradient (p, time, output, this->_subdomain_ids.get());
511 }
512 
513 
514 
516  const Real,
517  std::map<const Elem *, std::vector<Gradient>> & output,
518  const std::set<subdomain_id_type> * subdomain_ids)
519 {
520  libmesh_assert (this->initialized());
521 
522  // clear the output map
523  output.clear();
524 
525  // get the candidate elements
526  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  for (const auto & element : candidate_element)
531  {
532  // define a temporary vector to store all values
533  std::vector<Gradient> temp_output (cast_int<unsigned int>(this->_system_vars.size()));
534 
535  this->_gradient_on_elem(p, element, temp_output);
536 
537  // Insert temp_output into output
538  output.emplace(element, std::move(temp_output));
539  }
540 }
541 
542 
543 
545  const Elem * element,
546  std::vector<Gradient> & output)
547 {
548  libmesh_assert(element);
549 
550  // resize the output vector to the number of output values
551  // that the user told us
552  output.resize (this->_system_vars.size());
553 
554  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  const Point mapped_point (FEMap::inverse_map (dim, element,
560  p));
561 
562  std::vector<Point> point_list (1, mapped_point);
563 
564  // loop over all vars
565  for (auto index : index_range(this->_system_vars))
566  {
567  // the data for this variable
568  const unsigned int var = _system_vars[index];
569 
570  if (var == libMesh::invalid_uint)
571  {
573  index < _out_of_mesh_value.size());
574  output[index] = Gradient(_out_of_mesh_value(index));
575  continue;
576  }
577 
578  const FEType & fe_type = this->_dof_map.variable_type(var);
579 
580  // where the solution values for the var-th variable are stored
581  std::vector<dof_id_type> dof_indices;
582  this->_dof_map.dof_indices (element, dof_indices, var);
583 
584  // interpolate the solution
585  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  if (!element->infinite())
591  {
592  std::unique_ptr<FEBase> point_fe (FEBase::build(dim, fe_type));
593  const std::vector<std::vector<RealGradient>> & dphi = point_fe->get_dphi();
594  point_fe->reinit(element, &point_list);
595 
596  for (auto i : index_range(dof_indices))
597  grad.add_scaled(dphi[i][0], this->_vector(dof_indices[i]));
598 
599  }
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  FEComputeData data (this->_eqn_systems, mapped_point);
608  data.enable_derivative();
609  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  for (auto i : index_range(dof_indices))
614  {
615  // local coordinates
616  for (std::size_t v=0; v<dim; v++)
617  for (std::size_t xyz=0; xyz<LIBMESH_DIM; xyz++)
618  {
619  // FIXME: this needs better syntax: It is matrix-vector multiplication.
620  grad(xyz) += data.local_transform[v][xyz]
621  * data.dshape[i](v)
622  * this->_vector(dof_indices[i]);
623  }
624  }
625  }
626 #endif
627  output[index] = grad;
628  }
629 }
630 
631 
632 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
633 void MeshFunction::hessian (const Point & p,
634  const Real time,
635  std::vector<Tensor> & output)
636 {
637  this->hessian(p, time, output, this->_subdomain_ids.get());
638 }
639 
640 
641 
642 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  libmesh_assert (this->initialized());
648 
649  const Elem * element = this->find_element(p,subdomain_ids);
650 
651  if (!element)
652  {
653  output.resize(0);
654  }
655  else
656  {
657  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  output.resize (this->_system_vars.size());
666 
667 
668  {
669  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  const Point mapped_point (FEMap::inverse_map (dim, element,
676  p));
677 
678  std::vector<Point> point_list (1, mapped_point);
679 
680  // loop over all vars
681  for (auto index : index_range(this->_system_vars))
682  {
683  // the data for this variable
684  const unsigned int var = _system_vars[index];
685 
686  if (var == libMesh::invalid_uint)
687  {
689  index < _out_of_mesh_value.size());
690  output[index] = Tensor(_out_of_mesh_value(index));
691  continue;
692  }
693  const FEType & fe_type = this->_dof_map.variable_type(var);
694 
695  std::unique_ptr<FEBase> point_fe (FEBase::build(dim, fe_type));
696  const std::vector<std::vector<RealTensor>> & d2phi =
697  point_fe->get_d2phi();
698  point_fe->reinit(element, &point_list);
699 
700  // where the solution values for the var-th variable are stored
701  std::vector<dof_id_type> dof_indices;
702  this->_dof_map.dof_indices (element, dof_indices, var);
703 
704  // interpolate the solution
705  Tensor hess;
706 
707  for (auto i : index_range(dof_indices))
708  hess.add_scaled(d2phi[i][0], this->_vector(dof_indices[i]));
709 
710  output[index] = hess;
711  }
712  }
713  }
714 }
715 #endif
716 
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  if (this->_master != nullptr)
727  {
728  const MeshFunction * master =
729  cast_ptr<const MeshFunction *>(this->_master);
730  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  const Elem * element = (*_point_locator)(p, subdomain_ids);
738 
739  // Make sure that the element found is evaluable
740  return this->check_found_elem(element, p);
741 }
742 
743 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  if (this->_master != nullptr)
753  {
754  const MeshFunction * master =
755  cast_ptr<const MeshFunction *>(this->_master);
756  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  std::set<const Elem *> candidate_elements;
764  std::set<const Elem *> final_candidate_elements;
765  (*_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  for (const auto & element : candidate_elements)
770  final_candidate_elements.insert(this->check_found_elem(element, p));
771 
772  return final_candidate_elements;
773 }
774 
775 const Elem *
776 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  if (!element || (element->processor_id() == this->processor_id()))
782  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  if (_vector.type() == SERIAL)
790  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  if (_vector.type() == GHOSTED && _dof_map.is_evaluable(*element))
798  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  std::set<const Elem *> point_neighbors;
804  element->find_point_neighbors(p, point_neighbors);
805  element = nullptr;
806  for (const auto & elem : point_neighbors)
807  if (elem->processor_id() == this->processor_id())
808  {
809  element = elem;
810  break;
811  }
812 
813  return element;
814 }
815 
817 {
818  libmesh_assert (this->initialized());
819  return *_point_locator;
820 }
821 
823 {
824  libmesh_assert (this->initialized());
825  return *_point_locator;
826 }
827 
829 {
830  libmesh_assert (this->initialized());
831  _point_locator->enable_out_of_mesh_mode();
832  _out_of_mesh_mode = true;
834 }
835 
837 {
838  DenseVector<Number> v(1);
839  v(0) = value;
840  this->enable_out_of_mesh_mode(v);
841 }
842 
844 {
845  libmesh_assert (this->initialized());
846  _point_locator->disable_out_of_mesh_mode();
847  _out_of_mesh_mode = false;
848 }
849 
851 {
852  _point_locator->set_close_to_point_tol(tol);
853  _point_locator->set_contains_point_tol(tol);
854 }
855 
857 {
858  _point_locator->unset_close_to_point_tol();
859 }
860 
861 void MeshFunction::set_subdomain_ids(const std::set<subdomain_id_type> * subdomain_ids)
862 {
863  if (subdomain_ids)
864  _subdomain_ids = std::make_unique<std::set<subdomain_id_type>>(*subdomain_ids);
865  else
866  _subdomain_ids.reset();
867 }
868 
869 } // namespace libMesh
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:196
std::unique_ptr< std::set< subdomain_id_type > > _subdomain_ids
A default set of subdomain ids in which to search for points.
FEFamily family
The type of finite element.
Definition: fe_type.h:221
This is the EquationSystems class.
const EquationSystems & _eqn_systems
The equation systems handler, from which the data are gathered.
virtual std::unique_ptr< FunctionBase< Number > > clone() const override
void add_scaled(const TypeVector< T2 > &, const T &)
Add a scaled value to this vector without creating a temporary.
Definition: type_vector.h:629
BuildType
enum defining how to build the tree.
Definition: tree_base.h:58
const NumericVector< Number > & _vector
A reference to the vector that holds the data that is to be interpolated.
bool _out_of_mesh_mode
true if out-of-mesh mode is enabled.
const unsigned int invalid_uint
A number which is used quite often to represent an invalid or uninitialized value for an unsigned int...
Definition: libmesh.h:310
std::unique_ptr< PointLocatorBase > sub_point_locator() const
Definition: mesh_base.C:1675
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Definition: dof_map.C:2229
unsigned int dim
static Point inverse_map(const unsigned int dim, const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true, const bool extra_checks=true)
Definition: fe_map.C:1628
class FEComputeData hides arbitrary data to be passed to and from children of FEBase through the FEIn...
const Point & p
Holds the point where the data are to be computed.
const DofMap & _dof_map
Need access to the DofMap of the other system.
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:396
const Elem * find_element(const Point &p, const std::set< subdomain_id_type > *subdomain_ids=nullptr) const
Helper function to reduce code duplication.
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:2230
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
MeshBase & mesh
const FunctionBase * _master
Const pointer to our master, initialized to nullptr.
bool _initialized
When init() was called so that everything is ready for calls to operator() (...), then this bool is t...
Gradient gradient(const Point &p, const Real time=0.)
std::vector< Gradient > dshape
Storage for the computed shape derivative values.
This class defines a vector in LIBMESH_DIM dimensional Real or Complex space.
~MeshFunction()
Destructor.
std::vector< std::vector< Real > > local_transform
Storage for local to global mapping at p.
The libMesh namespace provides an interface to certain functionality in the library.
void disable_out_of_mesh_mode()
Disables out-of-mesh mode.
Real get_close_to_point_tol() const
Get the close-to-point tolerance.
This is the MeshBase class.
Definition: mesh_base.h:75
void _gradient_on_elem(const Point &p, const Elem *element, std::vector< Gradient > &output)
Helper function for finding a gradient as evaluated from a specific element.
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:179
DenseVector< Number > _out_of_mesh_value
Value to return outside the mesh if out-of-mesh mode is enabled.
void enable_out_of_mesh_mode(const DenseVector< Number > &value)
Enables out-of-mesh mode.
virtual void clear() override
Clears the function.
std::set< const Elem * > find_elements(const Point &p, const std::set< subdomain_id_type > *subdomain_ids=nullptr) const
std::map< const Elem *, Number > discontinuous_value(const Point &p, const Real time=0.)
void find_point_neighbors(const Point &p, std::set< const Elem *> &neighbor_set) const
This function finds all active elements (including this one) which are in the same manifold as this e...
Definition: elem.C:993
static void compute_data(const unsigned int dim, const FEType &fe_t, const Elem *elem, FEComputeData &data)
Lets the appropriate child of FEBase compute the requested data for the input specified in data...
NumberVectorValue Gradient
void unset_point_locator_tolerance()
Turn off the user-specified PointLocator tolerance.
virtual void init() override
Override the FunctionBase::init() member function.
static std::unique_ptr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.
libmesh_assert(ctx)
const PointLocatorBase & get_point_locator() const
This is the base class for point locators.
An object whose state is distributed along a set of processors.
void add_scaled(const TypeTensor< T2 > &, const T &)
Add a scaled tensor to this tensor without creating a temporary.
Definition: type_tensor.h:851
std::vector< Number > shape
Storage for the computed shape function values.
ParallelType type() const
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
NumberTensorValue Tensor
virtual unsigned short dim() const =0
virtual Number operator()(const Point &p, const Real time=0.) override
bool initialized() const
MeshFunction(const EquationSystems &eqn_systems, const NumericVector< Number > &vec, const DofMap &dof_map, std::vector< unsigned int > vars, const FunctionBase< Number > *master=nullptr)
Constructor for mesh based functions with vectors as return value.
Definition: mesh_function.C:42
const MeshBase & get_mesh() const
static const bool value
Definition: xdr_io.C:55
virtual unsigned int size() const override final
Definition: dense_vector.h:104
virtual bool infinite() const =0
This class provides function-like objects for data distributed over a mesh.
Definition: mesh_function.h:54
void set_point_locator_tolerance(Real tol)
We may want to specify a tolerance for the PointLocator to use, since in some cases the point we want...
const Elem * check_found_elem(const Elem *element, const Point &p) const
Helper function that is called by MeshFunction::find_element() and MeshFunction::find_elements() to e...
processor_id_type processor_id() const
void enable_derivative()
Enable the computation of shape gradients (dshape).
processor_id_type processor_id() const
Definition: dof_object.h:905
std::map< const Elem *, Gradient > discontinuous_gradient(const Point &p, const Real time=0.)
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
Tensor hessian(const Point &p, const Real time=0.)
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition: int_range.h:117
void set_subdomain_ids(const std::set< subdomain_id_type > *subdomain_ids)
Choose a default list of subdomain ids to be searched for points.
std::unique_ptr< PointLocatorBase > _point_locator
A point locator is needed to locate the points in the mesh.
bool is_evaluable(const DofObjectSubclass &obj, unsigned int var_num=libMesh::invalid_uint) const
Definition: dof_map.C:2677
This class defines a tensor in LIBMESH_DIM dimensional Real or Complex space.
const std::vector< unsigned int > _system_vars
The indices of the variables within the other system for which data are to be gathered.