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 706 : MeshFunction::MeshFunction (const EquationSystems & eqn_systems,
43 : const NumericVector<Number> & vec,
44 : const DofMap & dof_map,
45 : std::vector<unsigned int> vars,
46 706 : const FunctionBase<Number> * master) :
47 : FunctionBase<Number> (master),
48 : ParallelObject (eqn_systems),
49 666 : _eqn_systems (eqn_systems),
50 666 : _vector (vec),
51 666 : _dof_map (dof_map),
52 20 : _system_vars (std::move(vars)),
53 726 : _out_of_mesh_mode (false)
54 : {
55 706 : }
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 568 : MeshFunction::MeshFunction (const MeshFunction & mf):
75 568 : FunctionBase<Number> (mf._master),
76 568 : ParallelObject (mf._eqn_systems),
77 552 : _eqn_systems (mf._eqn_systems),
78 568 : _vector (mf._vector),
79 568 : _dof_map (mf._dof_map),
80 568 : _system_vars (mf._system_vars),
81 568 : _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 568 : if(mf.initialized())
86 : {
87 497 : this->MeshFunction::init();
88 :
89 497 : if(mf.get_point_locator().initialized())
90 497 : this->set_point_locator_tolerance(mf.get_point_locator().get_close_to_point_tol());
91 :
92 : }
93 :
94 568 : if (mf._subdomain_ids)
95 : _subdomain_ids =
96 : std::make_unique<std::set<subdomain_id_type>>
97 0 : (*mf._subdomain_ids);
98 568 : }
99 :
100 2073 : MeshFunction::~MeshFunction () = default;
101 :
102 :
103 :
104 1416 : 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 1416 : 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 1416 : const MeshBase & mesh = this->_eqn_systems.get_mesh();
119 2752 : _point_locator = mesh.sub_point_locator();
120 :
121 : // ready for use
122 1416 : 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 568 : std::unique_ptr<FunctionBase<Number>> MeshFunction::clone () const
140 : {
141 568 : 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 213 : 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 213 : this->discontinuous_value (p, time, buffer);
163 6 : std::map<const Elem *, Number> return_value;
164 497 : for (const auto & [elem, vec] : buffer)
165 284 : 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 219 : 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 213 : 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 213 : this->discontinuous_gradient (p, time, buffer);
188 6 : std::map<const Elem *, Gradient> return_value;
189 497 : for (const auto & [elem, vec] : buffer)
190 284 : 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 219 : 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 243382 : void MeshFunction::operator() (const Point & p,
209 : const Real time,
210 : DenseVector<Number> & output)
211 : {
212 243382 : this->operator() (p, time, output, this->_subdomain_ids.get());
213 243382 : }
214 :
215 1374254 : 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 1374254 : const Elem * element = this->find_element(p,subdomain_ids);
223 :
224 1374254 : 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 1373294 : 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 1146620 : 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 2746868 : for (auto index : index_range(this->_system_vars))
248 : {
249 1373574 : const unsigned int var = _system_vars[index];
250 :
251 1373574 : 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 1373574 : 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 1373574 : switch (fe_type.family)
264 : {
265 280 : 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 280 : output_size = output_size + dim;
275 280 : break;
276 : }
277 1373294 : default:
278 : {
279 1373294 : output_size++;
280 1373294 : break;
281 : }
282 : }
283 :
284 : }
285 1259957 : 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 2746868 : for (auto index : index_range(this->_system_vars))
292 : {
293 : // the data for this variable
294 1373574 : const unsigned int var = _system_vars[index];
295 :
296 1373574 : 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 1373574 : 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 1373574 : switch (fe_type.family)
310 : {
311 280 : 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 296 : FEComputeData data (this->_eqn_systems, mapped_point);
322 280 : 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 280 : this->_dof_map.dof_indices (element, dof_indices, var);
327 :
328 : // Calling the properly-transformed shape function using get_phi()
329 288 : std::unique_ptr<FEVectorBase> vec_fe = FEVectorBase::build(dim, fe_type);
330 288 : std::vector<Point> vec_pt = {pt};
331 8 : const std::vector<std::vector<RealGradient>> & vec_phi = vec_fe->get_phi();
332 280 : vec_fe->reinit(element, &vec_pt);
333 :
334 : // interpolate the solution
335 : {
336 840 : for (unsigned int d = 0; d < dim; d++)
337 : {
338 16 : Number value = 0.;
339 :
340 3360 : for (auto i : index_range(dof_indices))
341 2880 : value += this->_vector(dof_indices[i]) * (vec_phi[i][0](d));
342 :
343 576 : output(index+vec_count*(dim-1)+d) = value;
344 : }
345 : }
346 :
347 280 : vec_count++;
348 8 : break;
349 528 : }
350 1373294 : default:
351 : {
352 : // Build an FEComputeData that contains both input and output data
353 : // for the specific compute_data method.
354 :
355 1599968 : FEComputeData data (this->_eqn_systems, mapped_point);
356 :
357 1373294 : 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 1373294 : this->_dof_map.dof_indices (element, dof_indices, var);
362 :
363 : // interpolate the solution
364 : {
365 113337 : Number value = 0.;
366 :
367 15257006 : for (auto i : index_range(dof_indices))
368 15018896 : value += this->_vector(dof_indices[i]) * data.shape[i];
369 :
370 1486631 : output(index+vec_count*(dim-1)) = value;
371 : }
372 113337 : break;
373 1146620 : }
374 : }
375 :
376 : // next variable
377 : }
378 : }
379 : }
380 1374254 : }
381 :
382 213 : void MeshFunction::discontinuous_value (const Point & p,
383 : const Real time,
384 : std::map<const Elem *, DenseVector<Number>> & output)
385 : {
386 213 : this->discontinuous_value (p, time, output, this->_subdomain_ids.get());
387 213 : }
388 :
389 :
390 :
391 213 : 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 219 : 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 497 : for (const auto & element : candidate_element)
407 : {
408 284 : 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 284 : const unsigned int dim = element->dim();
418 :
419 : // define a temporary vector to store all values
420 292 : 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 284 : const Point mapped_point (FEMap::inverse_map (dim, element, p));
426 :
427 : // loop over all vars
428 568 : for (auto index : index_range(this->_system_vars))
429 : {
430 : // the data for this variable
431 284 : const unsigned int var = _system_vars[index];
432 :
433 284 : 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 284 : 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 300 : FEComputeData data (this->_eqn_systems, mapped_point);
447 :
448 284 : 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 284 : this->_dof_map.dof_indices (element, dof_indices, var);
453 :
454 : // interpolate the solution
455 : {
456 8 : Number value = 0.;
457 :
458 568 : for (auto i : index_range(dof_indices))
459 292 : value += this->_vector(dof_indices[i]) * data.shape[i];
460 :
461 284 : temp_output(index) = value;
462 : }
463 :
464 268 : }
465 :
466 : // next variable
467 : }
468 :
469 : // Insert temp_output into output
470 284 : output[element] = temp_output;
471 : }
472 213 : }
473 :
474 :
475 :
476 83 : void MeshFunction::gradient (const Point & p,
477 : const Real time,
478 : std::vector<Gradient> & output)
479 : {
480 83 : this->gradient(p, time, output, this->_subdomain_ids.get());
481 83 : }
482 :
483 :
484 :
485 1127603 : 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 1127603 : const Elem * element = this->find_element(p,subdomain_ids);
493 :
494 1127603 : this->_gradient_on_elem(p, element, output);
495 1127603 : }
496 :
497 :
498 213 : void MeshFunction::discontinuous_gradient (const Point & p,
499 : const Real time,
500 : std::map<const Elem *, std::vector<Gradient>> & output)
501 : {
502 213 : this->discontinuous_gradient (p, time, output, this->_subdomain_ids.get());
503 213 : }
504 :
505 :
506 :
507 213 : 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 219 : 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 497 : for (const auto & element : candidate_element)
523 : {
524 : // define a temporary vector to store all values
525 300 : std::vector<Gradient> temp_output (cast_int<unsigned int>(this->_system_vars.size()));
526 :
527 284 : this->_gradient_on_elem(p, element, temp_output);
528 :
529 : // Insert temp_output into output
530 276 : output.emplace(element, std::move(temp_output));
531 : }
532 213 : }
533 :
534 :
535 :
536 1127887 : 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 1221858 : output.resize (this->_system_vars.size());
543 :
544 1127887 : if (!element)
545 : {
546 50 : 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 200 : for (auto i : index_range(this->_system_vars))
552 153 : output[i] = Gradient(_out_of_mesh_value[i]);
553 50 : return;
554 : }
555 :
556 1127837 : 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 939897 : const Point mapped_point (FEMap::inverse_map (dim, element,
562 187940 : p));
563 :
564 1221807 : std::vector<Point> point_list (1, mapped_point);
565 :
566 : // loop over all vars
567 2255740 : for (auto index : index_range(this->_system_vars))
568 : {
569 : // the data for this variable
570 1127903 : const unsigned int var = _system_vars[index];
571 :
572 1127903 : if (var == libMesh::invalid_uint)
573 : {
574 2 : libmesh_assert (_out_of_mesh_mode &&
575 : index < _out_of_mesh_value.size());
576 33 : output[index] = Gradient(_out_of_mesh_value(index));
577 33 : continue;
578 : }
579 :
580 1127870 : 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 1127870 : 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 1221842 : 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 1127870 : point_fe->reinit(element, &point_list);
597 :
598 10149331 : for (auto i : index_range(dof_indices))
599 9773183 : grad.add_scaled(dphi[i][0], this->_vector(dof_indices[i]));
600 :
601 939926 : }
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 1221842 : output[index] = grad;
630 : }
631 : }
632 :
633 :
634 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
635 83 : void MeshFunction::hessian (const Point & p,
636 : const Real time,
637 : std::vector<Tensor> & output)
638 : {
639 83 : this->hessian(p, time, output, this->_subdomain_ids.get());
640 83 : }
641 :
642 :
643 :
644 1127603 : 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 1221566 : output.resize (this->_system_vars.size());
654 :
655 1127603 : const Elem * element = this->find_element(p,subdomain_ids);
656 :
657 1127603 : if (!element)
658 : {
659 48 : 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 192 : for (auto i : index_range(this->_system_vars))
665 147 : 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 1127555 : 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 939631 : const Point mapped_point (FEMap::inverse_map (dim, element,
684 187924 : p));
685 :
686 1221517 : std::vector<Point> point_list (1, mapped_point);
687 :
688 : // loop over all vars
689 2255180 : for (auto index : index_range(this->_system_vars))
690 : {
691 : // the data for this variable
692 1127625 : const unsigned int var = _system_vars[index];
693 :
694 1127625 : if (var == libMesh::invalid_uint)
695 : {
696 2 : libmesh_assert (_out_of_mesh_mode &&
697 : index < _out_of_mesh_value.size());
698 35 : output[index] = Tensor(_out_of_mesh_value(index));
699 35 : continue;
700 : }
701 1127590 : const FEType & fe_type = this->_dof_map.variable_type(var);
702 :
703 1221554 : 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 1127590 : 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 1127590 : this->_dof_map.dof_indices (element, dof_indices, var);
711 :
712 : // interpolate the solution
713 93964 : Tensor hess;
714 :
715 10148100 : for (auto i : index_range(dof_indices))
716 9772210 : hess.add_scaled(d2phi[i][0], this->_vector(dof_indices[i]));
717 :
718 1221554 : output[index] = hess;
719 939662 : }
720 : }
721 : }
722 : }
723 : #endif
724 :
725 3629460 : 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 3629460 : const Elem * element = (*_point_locator)(p, subdomain_ids);
746 :
747 : // Make sure that the element found is evaluable
748 3629460 : return this->check_found_elem(element, p);
749 : }
750 :
751 426 : 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 426 : (*_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 994 : for (const auto & element : candidate_elements)
778 568 : final_candidate_elements.insert(this->check_found_elem(element, p));
779 :
780 438 : return final_candidate_elements;
781 : }
782 :
783 : const Elem *
784 3630028 : 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 3630028 : 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 1903138 : 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 1944 : 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 278 : element->find_point_neighbors(p, point_neighbors);
813 13 : element = nullptr;
814 891 : for (const auto & elem : point_neighbors)
815 879 : if (elem->processor_id() == this->processor_id())
816 : {
817 13 : element = elem;
818 13 : break;
819 : }
820 :
821 13 : return element;
822 : }
823 :
824 994 : const PointLocatorBase & MeshFunction::get_point_locator () const
825 : {
826 28 : libmesh_assert (this->initialized());
827 994 : 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 213 : void MeshFunction::enable_out_of_mesh_mode(const DenseVector<Number> & value)
837 : {
838 6 : libmesh_assert (this->initialized());
839 213 : _point_locator->enable_out_of_mesh_mode();
840 213 : _out_of_mesh_mode = true;
841 6 : _out_of_mesh_value = value;
842 213 : }
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 568 : void MeshFunction::set_point_locator_tolerance(Real tol)
859 : {
860 568 : _point_locator->set_close_to_point_tol(tol);
861 568 : _point_locator->set_contains_point_tol(tol);
862 568 : }
863 :
864 0 : void MeshFunction::unset_point_locator_tolerance()
865 : {
866 0 : _point_locator->unset_close_to_point_tol();
867 0 : }
868 :
869 71 : void MeshFunction::set_subdomain_ids(const std::set<subdomain_id_type> * subdomain_ids)
870 : {
871 71 : if (subdomain_ids)
872 140 : _subdomain_ids = std::make_unique<std::set<subdomain_id_type>>(*subdomain_ids);
873 : else
874 0 : _subdomain_ids.reset();
875 71 : }
876 :
877 : } // namespace libMesh
|
