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 : #ifndef GENERIC_PROJECTOR_H
21 : #define GENERIC_PROJECTOR_H
22 :
23 : // C++ includes
24 : #include <vector>
25 :
26 : // libMesh includes
27 : #include "libmesh/dense_matrix.h"
28 : #include "libmesh/dense_vector.h"
29 : #include "libmesh/dof_map.h"
30 : #include "libmesh/elem.h"
31 : #include "libmesh/fe_base.h"
32 : #include "libmesh/fe_interface.h"
33 : #include "libmesh/int_range.h"
34 : #include "libmesh/libmesh_logging.h"
35 : #include "libmesh/mesh_tools.h"
36 : #include "libmesh/numeric_vector.h"
37 : #include "libmesh/quadrature.h"
38 : #include "libmesh/system.h"
39 : #include "libmesh/threads.h"
40 : #include "libmesh/tensor_tools.h"
41 : #include "libmesh/libmesh_common.h"
42 :
43 : // TIMPI includes
44 : #include "timpi/parallel_sync.h"
45 :
46 : // C++ Includes
47 : #include <memory>
48 :
49 : namespace libMesh
50 : {
51 :
52 : /**
53 : * For ease of communication, we allow users to translate their own
54 : * value types to a more easily computable (typically a vector of some
55 : * fixed-size type) output, by specializing these calls using
56 : * different types.
57 : */
58 : template <typename T>
59 : struct TypeToSend {
60 : typedef T type;
61 : };
62 :
63 : template <typename T>
64 1065 : const typename TypeToSend<T>::type convert_to_send (const T& in)
65 18576 : { return in; }
66 :
67 : template <typename SendT, typename T>
68 1065 : void convert_from_receive (SendT & received, T & converted)
69 18576 : { converted = received; }
70 :
71 : /**
72 : * The GenericProjector class implements the core of other projection
73 : * operations, using two input functors to read values to be projected
74 : * and an output functor to set degrees of freedom in the result.
75 : *
76 : * This may be executed in parallel on multiple threads.
77 : *
78 : * \author Roy H. Stogner
79 : * \date 2016
80 : */
81 : template <typename FFunctor, typename GFunctor,
82 : typename FValue, typename ProjectionAction>
83 : class GenericProjector
84 : {
85 : public:
86 : /**
87 : * Convenience typedef for the Node-to-attached-Elem mapping that
88 : * may be passed in to the constructor.
89 : */
90 : typedef std::unordered_map<dof_id_type, std::vector<dof_id_type>> NodesToElemMap;
91 :
92 : private:
93 : const System & system;
94 :
95 : // For TBB compatibility and thread safety we'll copy these in
96 : // operator()
97 : FFunctor & master_f;
98 :
99 : /**
100 : * Needed for C1 type elements only. master_g is either a shallow
101 : * copy of a GFunctor pointer passed to the constructor, or to
102 : * master_g_deepcopy.get(), depending on which constructor was
103 : * called.
104 : */
105 : std::unique_ptr<GFunctor> master_g_deepcopy;
106 : GFunctor * master_g;
107 :
108 : ProjectionAction & master_action;
109 : const std::vector<unsigned int> & variables;
110 :
111 : /**
112 : * nodes_to_elem is either a shallow copy of a map passed in to
113 : * the constructor, or points to nodes_to_elem_ourcopy, if no
114 : * such map was provided.
115 : */
116 : NodesToElemMap nodes_to_elem_ourcopy;
117 : NodesToElemMap * nodes_to_elem;
118 :
119 : bool done_saving_ids;
120 :
121 : public:
122 243232 : GenericProjector (const System & system_in,
123 : FFunctor & f_in,
124 : GFunctor * g_in,
125 : ProjectionAction & act_in,
126 : const std::vector<unsigned int> & variables_in,
127 : NodesToElemMap * nodes_to_elem_in = nullptr) :
128 228792 : system(system_in),
129 228792 : master_f(f_in),
130 228792 : master_g(g_in),
131 228792 : master_action(act_in),
132 228792 : variables(variables_in),
133 257672 : nodes_to_elem(nodes_to_elem_in)
134 : {
135 243232 : if (!nodes_to_elem_in)
136 : {
137 243232 : MeshTools::build_nodes_to_elem_map (system.get_mesh(), nodes_to_elem_ourcopy);
138 243232 : nodes_to_elem = &nodes_to_elem_ourcopy;
139 : }
140 243232 : }
141 :
142 : GenericProjector (const GenericProjector & in) :
143 : system(in.system),
144 : master_f(in.master_f),
145 : master_g_deepcopy(in.master_g ? std::make_unique<GFunctor>(*in.master_g) : nullptr),
146 : master_g(in.master_g ? master_g_deepcopy.get() : nullptr),
147 : master_action(in.master_action),
148 : variables(in.variables),
149 : nodes_to_elem(in.nodes_to_elem)
150 : {}
151 :
152 472024 : ~GenericProjector() = default;
153 :
154 : // The single "project" function handles all the threaded projection
155 : // calculations and intervening MPI communications
156 : void project(const ConstElemRange & range);
157 :
158 : // We generally need to hang on to every value we've calculated
159 : // until we're all done, because later projection calculations
160 : // depend on boundary data from earlier calculations.
161 : std::unordered_map<dof_id_type, typename FFunctor::ValuePushType> ids_to_save;
162 :
163 : // This needs to be sorted so we can get sorted dof indices cheaply
164 : // later
165 : typedef std::set<unsigned int> var_set;
166 :
167 : // We now need to divide our threadable work into stages, to allow
168 : // for the potential necessity of MPI communication in between them.
169 : struct SubFunctor {
170 : GenericProjector & projector;
171 :
172 : SubFunctor (GenericProjector & p);
173 :
174 : // While we're working on nodes, also figure out which ghosted
175 : // nodes will have data we might need to send to their owners
176 : // instead of being acted on by ourselves.
177 : //
178 : // We keep track of which dof ids we might need to send, and what
179 : // values those ids should get (along with a pprocessor_id to
180 : // leave invalid in case *we* can't compute those values either).
181 : std::unordered_map<dof_id_type,
182 : std::pair<typename FFunctor::ValuePushType,
183 : processor_id_type>> new_ids_to_push;
184 :
185 : // We'll hang on to any new ids to save on a per-thread basis
186 : // because we won't need them until subsequent phases
187 : std::unordered_map<dof_id_type, typename FFunctor::ValuePushType> new_ids_to_save;
188 :
189 : // Helper function for filling new_ids_to_push
190 : void find_dofs_to_send (const Node & node,
191 : const Elem & elem,
192 : unsigned short node_num,
193 : const var_set & vars);
194 :
195 : // When we have new data to act on, we may also need to save it
196 : // and get ready to push it.
197 : template <typename InsertInput,
198 : typename std::enable_if<
199 : std::is_same<typename ProjectionAction::InsertInput, InsertInput>::value,
200 : int>::type = 0>
201 : void insert_id(dof_id_type id, const InsertInput & val, processor_id_type pid);
202 :
203 : template <typename InsertInput,
204 : typename std::enable_if<
205 : !std::is_same<typename ProjectionAction::InsertInput, InsertInput>::value,
206 : int>::type = 0>
207 : void insert_id(dof_id_type id, const InsertInput & val, processor_id_type pid);
208 :
209 : template <typename InsertInput,
210 : typename std::enable_if<
211 : std::is_same<typename ProjectionAction::InsertInput, InsertInput>::value,
212 : int>::type = 0>
213 : void insert_ids(const std::vector<dof_id_type> & ids,
214 : const std::vector<InsertInput> & vals,
215 : processor_id_type pid);
216 :
217 : template <typename InsertInput,
218 : typename std::enable_if<
219 : !std::is_same<typename ProjectionAction::InsertInput, InsertInput>::value,
220 : int>::type = 0>
221 : void insert_ids(const std::vector<dof_id_type> & ids,
222 : const std::vector<InsertInput> & vals,
223 : processor_id_type pid);
224 :
225 : // Our subclasses will need to merge saved ids
226 : void join (const SubFunctor & other);
227 :
228 : protected:
229 : // For thread safety and efficiency we'll want thread-local copies
230 : // of various functors
231 : ProjectionAction action;
232 : FFunctor f;
233 :
234 : // Each stage of the work we do will require us to prepare
235 : // FEMContext objects to assist in the work.
236 : FEMContext context;
237 :
238 : // These get used often enough to cache them
239 : std::vector<FEContinuity> conts;
240 : std::vector<FEFieldType> field_types;
241 :
242 : const System & system;
243 : };
244 :
245 :
246 : typedef std::pair<const Node *, std::tuple<const Elem *, unsigned short, var_set>>
247 : node_projection;
248 :
249 : typedef StoredRange<std::vector<node_projection>::const_iterator,
250 : node_projection> node_range;
251 :
252 :
253 921215 : struct SubProjector : public SubFunctor {
254 : SubProjector (GenericProjector & p);
255 :
256 : using SubFunctor::action;
257 : using SubFunctor::f;
258 : using SubFunctor::system;
259 : using SubFunctor::context;
260 : using SubFunctor::conts;
261 : using SubFunctor::field_types;
262 : using SubFunctor::insert_id;
263 : using SubFunctor::insert_ids;
264 :
265 : protected:
266 : // Projections of C1 elements require a gradient as well
267 : std::unique_ptr<GFunctor> g;
268 :
269 : void construct_projection
270 : (const std::vector<dof_id_type> & dof_indices_var,
271 : const std::vector<unsigned int> & involved_dofs,
272 : unsigned int var_component,
273 : const Node * node,
274 : const FEGenericBase<typename FFunctor::RealType> & fe);
275 : };
276 :
277 :
278 : // First we'll copy DoFs where we can, while sorting remaining DoFs
279 : // for interpolation and projection later.
280 : struct SortAndCopy : public SubFunctor {
281 243232 : SortAndCopy (GenericProjector & p) : SubFunctor(p) {}
282 :
283 7718 : SortAndCopy (SortAndCopy & other, Threads::split) : SubFunctor(other.projector) {}
284 :
285 : using SubFunctor::action;
286 : using SubFunctor::f;
287 : using SubFunctor::system;
288 : using SubFunctor::context;
289 : using SubFunctor::insert_ids;
290 :
291 : void operator() (const ConstElemRange & range);
292 :
293 : // We need to merge saved multimaps when working from multiple
294 : // threads
295 : void join (const SortAndCopy & other);
296 :
297 : // For vertices, edges and sides, we need to know what element,
298 : // what local vertex/edge/side number, and what set of variables
299 : // to evaluate.
300 : typedef std::unordered_multimap<const Node *, std::tuple<const Elem *, unsigned short, var_set>> ves_multimap;
301 :
302 : ves_multimap vertices, edges, sides;
303 : std::vector<const Elem *> interiors;
304 : };
305 :
306 :
307 238808 : struct ProjectVertices : public SubProjector {
308 243232 : ProjectVertices (GenericProjector & p) : SubProjector(p) {}
309 :
310 5950 : ProjectVertices (ProjectVertices & p_v, Threads::split) : SubProjector(p_v.projector) {}
311 :
312 : using SubProjector::action;
313 : using SubProjector::f;
314 : using SubProjector::g;
315 : using SubProjector::system;
316 : using SubProjector::context;
317 : using SubProjector::conts;
318 : using SubFunctor::field_types;
319 :
320 : using SubProjector::insert_id;
321 : using SubProjector::insert_ids;
322 :
323 : void operator() (const node_range & range);
324 : };
325 :
326 :
327 7798 : struct ProjectEdges : public SubProjector {
328 243232 : ProjectEdges (GenericProjector & p) : SubProjector(p) {}
329 :
330 1157 : ProjectEdges (ProjectEdges & p_e, Threads::split) : SubProjector(p_e.projector) {}
331 :
332 : using SubProjector::action;
333 : using SubProjector::f;
334 : using SubProjector::g;
335 : using SubProjector::system;
336 : using SubProjector::context;
337 : using SubProjector::conts;
338 : using SubFunctor::field_types;
339 :
340 : using SubProjector::insert_id;
341 : using SubProjector::insert_ids;
342 :
343 : void operator() (const node_range & range);
344 : };
345 :
346 9030 : struct ProjectSides : public SubProjector {
347 243232 : ProjectSides (GenericProjector & p) : SubProjector(p) {}
348 :
349 3680 : ProjectSides (ProjectSides & p_s, Threads::split) : SubProjector(p_s.projector) {}
350 :
351 : using SubProjector::action;
352 : using SubProjector::f;
353 : using SubProjector::g;
354 : using SubProjector::system;
355 : using SubProjector::context;
356 : using SubProjector::conts;
357 : using SubFunctor::field_types;
358 :
359 : using SubProjector::insert_id;
360 : using SubProjector::insert_ids;
361 :
362 : void operator() (const node_range & range);
363 : };
364 :
365 :
366 : typedef StoredRange<std::vector<const Elem *>::const_iterator,
367 : const Elem *> interior_range;
368 :
369 863 : struct ProjectInteriors : public SubProjector {
370 243232 : ProjectInteriors (GenericProjector & p) : SubProjector(p) {}
371 :
372 1788 : ProjectInteriors (ProjectInteriors & p_i, Threads::split) : SubProjector(p_i.projector) {}
373 :
374 : using SubProjector::action;
375 : using SubProjector::f;
376 : using SubProjector::g;
377 : using SubProjector::system;
378 : using SubProjector::context;
379 : using SubProjector::conts;
380 : using SubFunctor::field_types;
381 :
382 : using SubProjector::insert_id;
383 : using SubProjector::insert_ids;
384 :
385 : void operator() (const interior_range & range);
386 : };
387 :
388 : template <typename Value>
389 : void send_and_insert_dof_values
390 : (std::unordered_map<dof_id_type, std::pair<Value, processor_id_type>> & ids_to_push,
391 : ProjectionAction & action) const;
392 : };
393 :
394 :
395 : /**
396 : * The VectorSetAction output functor class can be used with a
397 : * GenericProjector to set projection values (which must be of type
398 : * Val) as coefficients of the given NumericVector.
399 : *
400 : * \author Roy H. Stogner
401 : * \date 2016
402 : */
403 :
404 : template <typename Val>
405 : class VectorSetAction
406 : {
407 : private:
408 : NumericVector<Val> & target_vector;
409 :
410 : public:
411 : typedef Val InsertInput;
412 :
413 242602 : VectorSetAction(NumericVector<Val> & target_vec) :
414 242602 : target_vector(target_vec) {}
415 :
416 10907403 : void insert(dof_id_type id,
417 : Val val)
418 : {
419 : // Lock the new vector since it is shared among threads.
420 : {
421 2197440 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
422 10907403 : target_vector.set(id, val);
423 : }
424 10907403 : }
425 :
426 :
427 2409837 : void insert(const std::vector<dof_id_type> & dof_indices,
428 : const DenseVector<Val> & Ue)
429 : {
430 : const numeric_index_type
431 2409837 : first = target_vector.first_local_index(),
432 2409837 : last = target_vector.last_local_index();
433 :
434 275534 : unsigned int size = Ue.size();
435 :
436 275534 : libmesh_assert_equal_to(size, dof_indices.size());
437 :
438 : // Lock the new vector since it is shared among threads.
439 : {
440 551068 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
441 :
442 2546322 : for (unsigned int i = 0; i != size; ++i)
443 150652 : if ((dof_indices[i] >= first) && (dof_indices[i] < last))
444 150652 : target_vector.set(dof_indices[i], Ue(i));
445 : }
446 2409837 : }
447 :
448 : };
449 :
450 :
451 : /**
452 : * The FEMFunctionWrapper input functor class can be used with a
453 : * GenericProjector to read values from an FEMFunction.
454 : *
455 : * \author Roy H. Stogner
456 : * \date 2016
457 : */
458 :
459 : template <typename Output>
460 398262 : class FEMFunctionWrapper
461 : {
462 : protected:
463 : typedef typename TensorTools::MakeBaseNumber<Output>::type DofValueType;
464 :
465 : public:
466 : typedef typename TensorTools::MakeReal<Output>::type RealType;
467 : typedef DofValueType ValuePushType;
468 : typedef Output FunctorValue;
469 :
470 396942 : FEMFunctionWrapper(const FEMFunctionBase<Output> & f) : _f(f.clone()) {}
471 :
472 52951 : FEMFunctionWrapper(const FEMFunctionWrapper<Output> & fw) :
473 1031958 : _f(fw._f->clone()) {}
474 :
475 1013884 : void init_context (FEMContext & c) { _f->init_context(c); }
476 :
477 80590 : Output eval_at_node (const FEMContext & c,
478 : unsigned int i,
479 : unsigned int /*elem_dim*/,
480 : const Node & n,
481 : bool /*extra_hanging_dofs*/,
482 : const Real time)
483 897035 : { return _f->component(c, i, n, time); }
484 :
485 539701 : Output eval_at_point (const FEMContext & c,
486 : unsigned int i,
487 : const Point & n,
488 : const Real time,
489 : bool /*skip_context_check*/)
490 6727421 : { return _f->component(c, i, n, time); }
491 :
492 0 : void eval_mixed_derivatives (const FEMContext & /*c*/,
493 : unsigned int /*i*/,
494 : unsigned int /*dim*/,
495 : const Node & /*n*/,
496 : std::vector<Output> & /*derivs*/)
497 0 : { libmesh_error(); } // this is only for grid projections
498 :
499 141720 : bool is_grid_projection() { return false; }
500 :
501 0 : void eval_old_dofs (const Elem &,
502 : unsigned int,
503 : unsigned int,
504 : std::vector<dof_id_type> &,
505 : std::vector<Output> &)
506 0 : { libmesh_error(); }
507 :
508 0 : void eval_old_dofs (const Elem &,
509 : const FEType &,
510 : unsigned int,
511 : unsigned int,
512 : std::vector<dof_id_type> &,
513 : std::vector<Output> &)
514 0 : { libmesh_error(); }
515 :
516 : private:
517 : std::unique_ptr<FEMFunctionBase<Output>> _f;
518 : };
519 :
520 :
521 : /**
522 : * The OldSolutionBase input functor abstract base class is the root
523 : * of the OldSolutionValue and OldSolutionCoefs classes which allow a
524 : * GenericProjector to read old solution values or solution
525 : * interpolation coefficients for a just-refined-and-coarsened mesh.
526 : *
527 : * \author Roy H. Stogner
528 : * \date 2016
529 : */
530 :
531 : #ifdef LIBMESH_ENABLE_AMR
532 : template <typename Output,
533 : void (FEMContext::*point_output) (unsigned int,
534 : const Point &,
535 : Output &,
536 : const Real) const>
537 99918 : class OldSolutionBase
538 : {
539 : protected:
540 : typedef typename TensorTools::MakeBaseNumber<Output>::type DofValueType;
541 : public:
542 : typedef typename TensorTools::MakeReal<Output>::type RealType;
543 :
544 338051 : OldSolutionBase(const libMesh::System & sys_in,
545 : const std::vector<unsigned int> * vars) :
546 306781 : last_elem(nullptr),
547 306781 : sys(sys_in),
548 338051 : old_context(sys_in, vars, /* allocate_local_matrices= */ false)
549 : {
550 : // We'll be queried for components but we'll typically be looking
551 : // up data by variables, and those indices don't always match
552 1240531 : auto make_lookup = [this](unsigned int v)
553 : {
554 415488 : const unsigned int vcomp = sys.variable_scalar_number(v,0);
555 433364 : if (vcomp >= component_to_var.size())
556 415488 : component_to_var.resize(vcomp+1, static_cast<unsigned int>(-1));
557 433364 : component_to_var[vcomp] = v;
558 : };
559 :
560 338051 : if (vars)
561 753539 : for (auto v : *vars)
562 415488 : make_lookup(v);
563 : else
564 0 : for (auto v : make_range(sys.n_vars()))
565 0 : make_lookup(v);
566 338051 : }
567 :
568 : OldSolutionBase(const OldSolutionBase & in) :
569 : last_elem(nullptr),
570 : sys(in.sys),
571 : old_context(sys, in.old_context.active_vars(), /* allocate_local_matrices= */ false),
572 : component_to_var(in.component_to_var)
573 : {
574 : }
575 :
576 : static void get_shape_outputs(FEAbstract & fe);
577 :
578 : // Integrating on new mesh elements, we won't yet have an up to date
579 : // current_local_solution.
580 250025 : void init_context (FEMContext & c)
581 : {
582 12423 : c.set_algebraic_type(FEMContext::DOFS_ONLY);
583 :
584 : const std::set<unsigned char> & elem_dims =
585 12423 : old_context.elem_dimensions();
586 :
587 : // Loop over variables and dimensions, to prerequest
588 501145 : for (const auto & dim : elem_dims)
589 : {
590 12463 : FEAbstract * fe = nullptr;
591 251120 : if (old_context.active_vars())
592 558081 : for (auto var : *old_context.active_vars())
593 : {
594 306961 : old_context.get_element_fe(var, fe, dim);
595 14240 : get_shape_outputs(*fe);
596 : }
597 : else
598 0 : for (auto var : make_range(sys.n_vars()))
599 : {
600 0 : old_context.get_element_fe(var, fe, dim);
601 0 : get_shape_outputs(*fe);
602 : }
603 : }
604 250025 : }
605 :
606 616957 : bool is_grid_projection() { return true; }
607 :
608 : protected:
609 : void check_old_context (const FEMContext & c)
610 : {
611 : LOG_SCOPE ("check_old_context(c)", "OldSolutionBase");
612 : const Elem & elem = c.get_elem();
613 : if (last_elem != &elem)
614 : {
615 : if (elem.refinement_flag() == Elem::JUST_REFINED)
616 : {
617 : old_context.pre_fe_reinit(sys, elem.parent());
618 : }
619 : else if (elem.refinement_flag() == Elem::JUST_COARSENED)
620 : {
621 : libmesh_error();
622 : }
623 : else
624 : {
625 : if (!elem.get_old_dof_object())
626 : {
627 : libmesh_error();
628 : }
629 :
630 : old_context.pre_fe_reinit(sys, &elem);
631 : }
632 :
633 : last_elem = &elem;
634 : }
635 : else
636 : {
637 : libmesh_assert(old_context.has_elem());
638 : }
639 : }
640 :
641 :
642 4409936 : bool check_old_context (const FEMContext & c, const Point & p)
643 : {
644 737282 : LOG_SCOPE ("check_old_context(c,p)", "OldSolutionBase");
645 737308 : const Elem & elem = c.get_elem();
646 4409936 : if (last_elem != &elem)
647 : {
648 3130220 : if (elem.refinement_flag() == Elem::JUST_REFINED)
649 : {
650 2951565 : old_context.pre_fe_reinit(sys, elem.parent());
651 : }
652 164841 : else if (elem.refinement_flag() == Elem::JUST_COARSENED)
653 : {
654 : // Find the child with this point. Use out_of_elem_tol
655 : // (in physical space, which may correspond to a large
656 : // tolerance in master space!) to allow for out-of-element
657 : // finite differencing of mixed gradient terms. Pray we
658 : // have no quadrature locations which are within 1e-5 of
659 : // the element subdivision boundary but are not exactly on
660 : // that boundary.
661 156809 : const Real master_tol = out_of_elem_tol / elem.hmax() * 2;
662 :
663 294154 : for (auto & child : elem.child_ref_range())
664 294154 : if (child.close_to_point(p, master_tol))
665 : {
666 156809 : old_context.pre_fe_reinit(sys, &child);
667 13151 : break;
668 : }
669 :
670 13151 : libmesh_assert
671 : (old_context.get_elem().close_to_point(p, master_tol));
672 : }
673 : else
674 : {
675 8032 : if (!elem.get_old_dof_object())
676 0 : return false;
677 :
678 8032 : old_context.pre_fe_reinit(sys, &elem);
679 : }
680 :
681 2879256 : last_elem = &elem;
682 : }
683 : else
684 : {
685 117691 : libmesh_assert(old_context.has_elem());
686 :
687 1530680 : const Real master_tol = out_of_elem_tol / elem.hmax() * 2;
688 :
689 1530680 : if (!old_context.get_elem().close_to_point(p, master_tol))
690 : {
691 428 : libmesh_assert_equal_to
692 : (elem.refinement_flag(), Elem::JUST_COARSENED);
693 :
694 30099 : for (auto & child : elem.child_ref_range())
695 30099 : if (child.close_to_point(p, master_tol))
696 : {
697 10731 : old_context.pre_fe_reinit(sys, &child);
698 428 : break;
699 : }
700 :
701 428 : libmesh_assert
702 : (old_context.get_elem().close_to_point(p, master_tol));
703 : }
704 : }
705 :
706 368641 : return true;
707 : }
708 :
709 : protected:
710 : const Elem * last_elem;
711 : const System & sys;
712 : FEMContext old_context;
713 : std::vector<unsigned int> component_to_var;
714 :
715 : static const Real out_of_elem_tol;
716 : };
717 :
718 :
719 : /**
720 : * The OldSolutionValue input functor class can be used with
721 : * GenericProjector to read values from a solution on a
722 : * just-refined-and-coarsened mesh.
723 : *
724 : * \author Roy H. Stogner
725 : * \date 2016
726 : */
727 : template <typename Output,
728 : void (FEMContext::*point_output) (unsigned int,
729 : const Point &,
730 : Output &,
731 : const Real) const>
732 100318 : class OldSolutionValue : public OldSolutionBase<Output, point_output>
733 : {
734 : public:
735 : using typename OldSolutionBase<Output, point_output>::RealType;
736 : using typename OldSolutionBase<Output, point_output>::DofValueType;
737 : typedef Output FunctorValue;
738 : typedef DofValueType ValuePushType;
739 :
740 46559 : OldSolutionValue(const libMesh::System & sys_in,
741 : const NumericVector<Number> & old_sol,
742 : const std::vector<unsigned int> * vars) :
743 : OldSolutionBase<Output, point_output>(sys_in, vars),
744 44971 : old_solution(old_sol)
745 : {
746 3176 : this->old_context.set_algebraic_type(FEMContext::OLD);
747 3176 : this->old_context.set_custom_solution(&old_solution);
748 3176 : }
749 :
750 246719 : OldSolutionValue(const OldSolutionValue & in) :
751 : OldSolutionBase<Output, point_output>(in.sys, in.old_context.active_vars()),
752 247355 : old_solution(in.old_solution)
753 : {
754 12281 : this->old_context.set_algebraic_type(FEMContext::OLD);
755 12281 : this->old_context.set_custom_solution(&old_solution);
756 234438 : }
757 :
758 :
759 : Output eval_at_node (const FEMContext & c,
760 : unsigned int i,
761 : unsigned int elem_dim,
762 : const Node & n,
763 : bool /* extra_hanging_dofs */,
764 : Real /* time */ =0.);
765 :
766 4309877 : Output eval_at_point(const FEMContext & c,
767 : unsigned int i,
768 : const Point & p,
769 : Real /* time */,
770 : bool skip_context_check)
771 : {
772 719610 : LOG_SCOPE ("eval_at_point()", "OldSolutionValue");
773 :
774 4309877 : if (!skip_context_check)
775 4290497 : if (!this->check_old_context(c, p))
776 0 : return Output(0);
777 :
778 : // Handle offset from non-scalar components in previous variables
779 359805 : libmesh_assert_less(i, this->component_to_var.size());
780 4329881 : unsigned int var = this->component_to_var[i];
781 :
782 314760 : Output n;
783 4309877 : (this->old_context.*point_output)(var, p, n, this->out_of_elem_tol);
784 4309877 : return n;
785 : }
786 :
787 : template <typename T = Output,
788 : typename std::enable_if<std::is_same<T, Number>::value, int>::type = 0>
789 2816 : void eval_mixed_derivatives(const FEMContext & libmesh_dbg_var(c),
790 : unsigned int i,
791 : unsigned int dim,
792 : const Node & n,
793 : std::vector<Output> & derivs)
794 : {
795 716 : LOG_SCOPE ("eval_mixed_derivatives", "OldSolutionValue");
796 :
797 : // This should only be called on vertices
798 358 : libmesh_assert_less(c.get_elem().get_node_index(&n),
799 : c.get_elem().n_vertices());
800 :
801 : // Handle offset from non-scalar components in previous variables
802 358 : libmesh_assert_less(i, this->component_to_var.size());
803 2816 : unsigned int var = this->component_to_var[i];
804 :
805 : // We have 1 mixed derivative in 2D, 4 in 3D
806 2816 : const unsigned int n_mixed = (dim-1) * (dim-1);
807 2816 : derivs.resize(n_mixed);
808 :
809 : // Be sure to handle cases where the variable wasn't defined on
810 : // this node (e.g. due to changing subdomain support)
811 358 : const DofObject * old_dof_object = n.get_old_dof_object();
812 2816 : if (old_dof_object &&
813 5274 : old_dof_object->n_vars(this->sys.number()) &&
814 2458 : old_dof_object->n_comp(this->sys.number(), var))
815 : {
816 716 : const dof_id_type first_old_id =
817 2816 : old_dof_object->dof_number(this->sys.number(), var, dim);
818 3174 : std::vector<dof_id_type> old_ids(n_mixed);
819 358 : std::iota(old_ids.begin(), old_ids.end(), first_old_id);
820 2816 : old_solution.get(old_ids, derivs);
821 : }
822 : else
823 : {
824 0 : std::fill(derivs.begin(), derivs.end(), 0);
825 : }
826 2816 : }
827 :
828 : template <typename T = Output,
829 : typename std::enable_if<!std::is_same<T, Number>::value, int>::type = 0>
830 0 : void eval_mixed_derivatives(
831 : const FEMContext &, unsigned int, unsigned int, const Node &, std::vector<Output> &)
832 : {
833 0 : libmesh_error_msg("eval_mixed_derivatives should only be applicable for Hermite finite element "
834 : "types. I don't know how you got here");
835 : }
836 :
837 2853864 : void eval_old_dofs (const Elem & elem,
838 : unsigned int node_num,
839 : unsigned int var_num,
840 : std::vector<dof_id_type> & indices,
841 : std::vector<DofValueType> & values)
842 : {
843 692394 : LOG_SCOPE ("eval_old_dofs(node)", "OldSolutionValue");
844 :
845 : // We may be reusing a std::vector here, but the following
846 : // dof_indices call appends without first clearing.
847 346197 : indices.clear();
848 :
849 3200061 : this->sys.get_dof_map().dof_indices(elem, node_num, indices, var_num);
850 :
851 692394 : std::vector<dof_id_type> old_indices;
852 :
853 3200061 : this->sys.get_dof_map().old_dof_indices(elem, node_num, old_indices, var_num);
854 :
855 346197 : libmesh_assert_equal_to (old_indices.size(), indices.size());
856 :
857 : // We may have invalid_id in cases where no old DoF existed, e.g.
858 : // due to expansion of a subdomain-restricted variable's subdomain
859 346197 : bool invalid_old_index = false;
860 5802109 : for (const auto & di : old_indices)
861 2948245 : if (di == DofObject::invalid_id)
862 0 : invalid_old_index = true;
863 :
864 2853864 : values.resize(old_indices.size());
865 2853864 : if (invalid_old_index)
866 : {
867 0 : for (auto i : index_range(old_indices))
868 : {
869 0 : const dof_id_type di = old_indices[i];
870 0 : if (di == DofObject::invalid_id)
871 0 : values[i] = 0;
872 : else
873 0 : values[i] = old_solution(di);
874 : }
875 : }
876 : else
877 2853864 : old_solution.get(old_indices, values);
878 2853864 : }
879 :
880 :
881 2287366 : void eval_old_dofs (const Elem & elem,
882 : const FEType & fe_type,
883 : unsigned int sys_num,
884 : unsigned int var_num,
885 : std::vector<dof_id_type> & indices,
886 : std::vector<DofValueType> & values)
887 : {
888 524980 : LOG_SCOPE ("eval_old_dofs(elem)", "OldSolutionValue");
889 :
890 : // We're only to be asked for old dofs on elements that can copy
891 : // them through DO_NOTHING or through refinement.
892 524512 : const Elem & old_elem =
893 2287366 : (elem.refinement_flag() == Elem::JUST_REFINED) ?
894 936 : *elem.parent() : elem;
895 :
896 : // If there are any element-based DOF numbers, get them
897 2287366 : const unsigned int nc = FEInterface::n_dofs_per_elem(fe_type, &elem);
898 :
899 2287366 : indices.resize(nc);
900 :
901 : // We should never have fewer dofs than necessary on an
902 : // element unless we're getting indices on a parent element,
903 : // which we should never need, or getting indices on a newly
904 : // expanded subdomain, in which case we initialize new values to
905 : // zero.
906 2287366 : if (nc != 0)
907 : {
908 1811 : const DofObject & old_dof_object = old_elem.get_old_dof_object_ref();
909 1811 : libmesh_assert_greater(elem.n_systems(), sys_num);
910 :
911 : const std::pair<unsigned int, unsigned int>
912 19965 : vg_and_offset = elem.var_to_vg_and_offset(sys_num,var_num);
913 19965 : const unsigned int vg = vg_and_offset.first;
914 19965 : const unsigned int vig = vg_and_offset.second;
915 :
916 21776 : unsigned int n_comp = old_dof_object.n_comp_group(sys_num,vg);
917 21776 : n_comp = std::min(n_comp, nc);
918 :
919 21776 : std::vector<dof_id_type> old_dof_indices(n_comp);
920 :
921 64751 : for (unsigned int i=0; i != n_comp; ++i)
922 : {
923 3554 : const dof_id_type d_old =
924 39421 : old_dof_object.dof_number(sys_num, vg, vig, i, n_comp);
925 3554 : const dof_id_type d_new =
926 39421 : elem.dof_number(sys_num, vg, vig, i, n_comp);
927 3554 : libmesh_assert_not_equal_to (d_old, DofObject::invalid_id);
928 3554 : libmesh_assert_not_equal_to (d_new, DofObject::invalid_id);
929 :
930 42975 : old_dof_indices[i] = d_old;
931 46529 : indices[i] = d_new;
932 : }
933 :
934 21776 : values.resize(n_comp);
935 21776 : old_solution.get(old_dof_indices, values);
936 :
937 21776 : for (unsigned int i=n_comp; i != nc; ++i)
938 : {
939 0 : const dof_id_type d_new =
940 0 : elem.dof_number(sys_num, vg, vig, i, n_comp);
941 0 : libmesh_assert_not_equal_to (d_new, DofObject::invalid_id);
942 0 : indices[i] = d_new;
943 : }
944 :
945 21776 : values.resize(nc, 0);
946 : }
947 : else
948 260679 : values.clear();
949 2287366 : }
950 :
951 : private:
952 : const NumericVector<Number> & old_solution;
953 : };
954 :
955 : template<>
956 : inline void
957 12634 : OldSolutionBase<Number, &FEMContext::point_value>::get_shape_outputs(FEAbstract & fe)
958 : {
959 260738 : fe.request_phi();
960 12634 : }
961 :
962 :
963 : template<>
964 : inline void
965 636 : OldSolutionBase<Gradient, &FEMContext::point_gradient>::get_shape_outputs(FEAbstract & fe)
966 : {
967 19588 : fe.request_dphi();
968 636 : }
969 :
970 : template<>
971 : inline void
972 970 : OldSolutionBase<Gradient, &FEMContext::point_value>::get_shape_outputs(FEAbstract & fe)
973 : {
974 26635 : fe.request_phi();
975 970 : }
976 :
977 :
978 : template<>
979 : inline void
980 0 : OldSolutionBase<Tensor, &FEMContext::point_gradient>::get_shape_outputs(FEAbstract & fe)
981 : {
982 0 : fe.request_dphi();
983 0 : }
984 :
985 :
986 : #ifdef LIBMESH_USE_COMPLEX_NUMBERS
987 : template<>
988 : inline void
989 : OldSolutionBase<Real, &FEMContext::point_value>::get_shape_outputs(FEAbstract & fe)
990 : {
991 0 : fe.request_phi();
992 : }
993 :
994 :
995 : template<>
996 : inline void
997 : OldSolutionBase<RealGradient, &FEMContext::point_gradient>::get_shape_outputs(FEAbstract & fe)
998 : {
999 0 : fe.request_dphi();
1000 : }
1001 :
1002 : template<>
1003 : inline void
1004 : OldSolutionBase<RealGradient, &FEMContext::point_value>::get_shape_outputs(FEAbstract & fe)
1005 : {
1006 : fe.request_phi();
1007 : }
1008 :
1009 :
1010 : template<>
1011 : inline void
1012 : OldSolutionBase<RealTensor, &FEMContext::point_gradient>::get_shape_outputs(FEAbstract & fe)
1013 : {
1014 : fe.request_dphi();
1015 : }
1016 : #endif // LIBMESH_USE_COMPLEX_NUMBERS
1017 :
1018 :
1019 : template<>
1020 : inline
1021 : Number
1022 1779507 : OldSolutionValue<Number, &FEMContext::point_value>::
1023 : eval_at_node(const FEMContext & c,
1024 : unsigned int i,
1025 : unsigned int /* elem_dim */,
1026 : const Node & n,
1027 : bool extra_hanging_dofs,
1028 : Real /* time */)
1029 : {
1030 312074 : LOG_SCOPE ("Number eval_at_node()", "OldSolutionValue");
1031 :
1032 : // This should only be called on vertices
1033 156037 : libmesh_assert_less(c.get_elem().get_node_index(&n),
1034 : c.get_elem().n_vertices());
1035 :
1036 : // Handle offset from non-scalar components in previous variables
1037 156037 : libmesh_assert_less(i, this->component_to_var.size());
1038 1779507 : unsigned int var = this->component_to_var[i];
1039 :
1040 : // Optimize for the common case, where this node was part of the
1041 : // old solution.
1042 : //
1043 : // Be sure to handle cases where the variable wasn't defined on
1044 : // this node (due to changing subdomain support) or where the
1045 : // variable has no components on this node (due to Elem order
1046 : // exceeding FE order) or where the old_dof_object dofs might
1047 : // correspond to non-vertex dofs (due to extra_hanging_dofs and
1048 : // refinement)
1049 :
1050 312075 : const Elem::RefinementState flag = c.get_elem().refinement_flag();
1051 :
1052 156037 : const DofObject * old_dof_object = n.get_old_dof_object();
1053 1295022 : if (old_dof_object &&
1054 1236513 : (!extra_hanging_dofs ||
1055 1045776 : flag == Elem::JUST_COARSENED ||
1056 1083563 : flag == Elem::DO_NOTHING) &&
1057 3950848 : old_dof_object->n_vars(sys.number()) &&
1058 1083563 : old_dof_object->n_comp(sys.number(), var))
1059 : {
1060 : const dof_id_type old_id =
1061 1060420 : old_dof_object->dof_number(sys.number(), var, 0);
1062 1060420 : return old_solution(old_id);
1063 : }
1064 :
1065 719087 : return this->eval_at_point(c, i, n, 0, false);
1066 : }
1067 :
1068 : template <>
1069 : inline
1070 : Gradient
1071 21744 : OldSolutionValue<Gradient, &FEMContext::point_value>::
1072 : eval_at_node(const FEMContext & c,
1073 : unsigned int i,
1074 : unsigned int /* elem_dim */,
1075 : const Node & n,
1076 : bool extra_hanging_dofs,
1077 : Real /* time */)
1078 : {
1079 2604 : LOG_SCOPE ("Number eval_at_node()", "OldSolutionValue");
1080 :
1081 : // This should only be called on vertices
1082 1302 : libmesh_assert_less(c.get_elem().get_node_index(&n),
1083 : c.get_elem().n_vertices());
1084 :
1085 : // Handle offset from non-scalar components in previous variables
1086 1302 : libmesh_assert_less(i, this->component_to_var.size());
1087 21744 : unsigned int var = this->component_to_var[i];
1088 :
1089 : // Optimize for the common case, where this node was part of the
1090 : // old solution.
1091 : //
1092 : // Be sure to handle cases where the variable wasn't defined on
1093 : // this node (due to changing subdomain support) or where the
1094 : // variable has no components on this node (due to Elem order
1095 : // exceeding FE order) or where the old_dof_object dofs might
1096 : // correspond to non-vertex dofs (due to extra_hanging_dofs and
1097 : // refinement)
1098 :
1099 2604 : const auto & elem = c.get_elem();
1100 :
1101 2604 : const Elem::RefinementState flag = elem.refinement_flag();
1102 :
1103 1302 : const DofObject * old_dof_object = n.get_old_dof_object();
1104 16518 : if (old_dof_object &&
1105 16152 : (!extra_hanging_dofs ||
1106 14280 : flag == Elem::JUST_COARSENED ||
1107 15216 : flag == Elem::DO_NOTHING) &&
1108 52176 : old_dof_object->n_vars(sys.number()) &&
1109 15216 : old_dof_object->n_comp(sys.number(), var))
1110 : {
1111 936 : Gradient return_val;
1112 :
1113 58437 : for (auto dim : make_range(elem.dim()))
1114 : {
1115 : const dof_id_type old_id =
1116 43221 : old_dof_object->dof_number(sys.number(), var, dim);
1117 43221 : return_val(dim) = old_solution(old_id);
1118 : }
1119 :
1120 15216 : return return_val;
1121 : }
1122 :
1123 6528 : return this->eval_at_point(c, i, n, 0, false);
1124 : }
1125 :
1126 :
1127 :
1128 : template<>
1129 : inline
1130 : Gradient
1131 12586 : OldSolutionValue<Gradient, &FEMContext::point_gradient>::
1132 : eval_at_node(const FEMContext & c,
1133 : unsigned int i,
1134 : unsigned int elem_dim,
1135 : const Node & n,
1136 : bool extra_hanging_dofs,
1137 : Real /* time */)
1138 : {
1139 2116 : LOG_SCOPE ("Gradient eval_at_node()", "OldSolutionValue");
1140 :
1141 : // This should only be called on vertices
1142 1058 : libmesh_assert_less(c.get_elem().get_node_index(&n),
1143 : c.get_elem().n_vertices());
1144 :
1145 : // Handle offset from non-scalar components in previous variables
1146 1058 : libmesh_assert_less(i, this->component_to_var.size());
1147 12586 : unsigned int var = this->component_to_var[i];
1148 :
1149 : // Optimize for the common case, where this node was part of the
1150 : // old solution.
1151 : //
1152 : // Be sure to handle cases where the variable wasn't defined on
1153 : // this node (due to changing subdomain support) or where the
1154 : // variable has no components on this node (due to Elem order
1155 : // exceeding FE order) or where the old_dof_object dofs might
1156 : // correspond to non-vertex dofs (due to extra_hanging_dofs and
1157 : // refinement)
1158 :
1159 2115 : const Elem::RefinementState flag = c.get_elem().refinement_flag();
1160 :
1161 1058 : const DofObject * old_dof_object = n.get_old_dof_object();
1162 13644 : if (old_dof_object &&
1163 13643 : (!extra_hanging_dofs ||
1164 12586 : flag == Elem::JUST_COARSENED ||
1165 451 : flag == Elem::DO_NOTHING) &&
1166 14511 : old_dof_object->n_vars(sys.number()) &&
1167 451 : old_dof_object->n_comp(sys.number(), var))
1168 : {
1169 35 : Gradient g;
1170 902 : for (unsigned int d = 0; d != elem_dim; ++d)
1171 : {
1172 : const dof_id_type old_id =
1173 451 : old_dof_object->dof_number(sys.number(), var, d+1);
1174 451 : g(d) = old_solution(old_id);
1175 : }
1176 451 : return g;
1177 : }
1178 :
1179 12135 : return this->eval_at_point(c, i, n, 0, false);
1180 : }
1181 :
1182 : template<>
1183 : inline
1184 : Tensor
1185 0 : OldSolutionValue<Tensor, &FEMContext::point_gradient>::
1186 : eval_at_node(const FEMContext &,
1187 : unsigned int,
1188 : unsigned int,
1189 : const Node &,
1190 : bool,
1191 : Real)
1192 : {
1193 0 : libmesh_error_msg("You shouldn't need to call eval_at_node for the gradient "
1194 : "functor for a vector-valued finite element type");
1195 : return Tensor();
1196 : }
1197 :
1198 : template <typename Output,
1199 : void (FEMContext::*point_output) (unsigned int,
1200 : const Point &,
1201 : Output &,
1202 : const Real) const>
1203 : const Real OldSolutionBase<Output, point_output>::out_of_elem_tol = 10 * TOLERANCE;
1204 :
1205 : #endif // LIBMESH_ENABLE_AMR
1206 :
1207 : /**
1208 : * Function definitions
1209 : */
1210 :
1211 : template <typename FFunctor, typename GFunctor,
1212 : typename FValue, typename ProjectionAction>
1213 243232 : void GenericProjector<FFunctor, GFunctor, FValue, ProjectionAction>::project
1214 : (const ConstElemRange & range)
1215 : {
1216 14440 : LOG_SCOPE ("project", "GenericProjector");
1217 :
1218 : // Unless we split sort and copy into two passes we can't know for
1219 : // sure ahead of time whether we need to save the copied ids
1220 243232 : done_saving_ids = false;
1221 :
1222 257672 : SortAndCopy sort_work(*this);
1223 243232 : Threads::parallel_reduce (range, sort_work);
1224 243232 : ProjectionAction action(master_action);
1225 :
1226 : // Keep track of dof ids and values to send to other ranks
1227 : std::unordered_map<dof_id_type, std::pair<typename FFunctor::ValuePushType, processor_id_type>>
1228 14440 : ids_to_push;
1229 :
1230 236012 : ids_to_push.insert(sort_work.new_ids_to_push.begin(),
1231 : sort_work.new_ids_to_push.end());
1232 236012 : ids_to_save.insert(sort_work.new_ids_to_save.begin(),
1233 : sort_work.new_ids_to_save.end());
1234 :
1235 257672 : std::vector<node_projection> vertices(sort_work.vertices.begin(),
1236 : sort_work.vertices.end());
1237 :
1238 473314 : done_saving_ids = sort_work.edges.empty() &&
1239 446004 : sort_work.sides.empty() && sort_work.interiors.empty();
1240 :
1241 : {
1242 14440 : ProjectVertices project_vertices(*this);
1243 479244 : Threads::parallel_reduce (node_range(&vertices), project_vertices);
1244 7220 : libmesh_merge_move(ids_to_push, project_vertices.new_ids_to_push);
1245 7220 : libmesh_merge_move(ids_to_save, project_vertices.new_ids_to_save);
1246 : }
1247 :
1248 243232 : done_saving_ids = sort_work.sides.empty() && sort_work.interiors.empty();
1249 :
1250 243232 : this->send_and_insert_dof_values(ids_to_push, action);
1251 :
1252 : {
1253 257672 : std::vector<node_projection> edges(sort_work.edges.begin(), sort_work.edges.end());
1254 14440 : ProjectEdges project_edges(*this);
1255 479244 : Threads::parallel_reduce (node_range(&edges), project_edges);
1256 7220 : libmesh_merge_move(ids_to_push, project_edges.new_ids_to_push);
1257 7220 : libmesh_merge_move(ids_to_save, project_edges.new_ids_to_save);
1258 228792 : }
1259 :
1260 243232 : done_saving_ids = sort_work.interiors.empty();
1261 :
1262 243232 : this->send_and_insert_dof_values(ids_to_push, action);
1263 :
1264 : {
1265 257672 : std::vector<node_projection> sides(sort_work.sides.begin(), sort_work.sides.end());
1266 14440 : ProjectSides project_sides(*this);
1267 479244 : Threads::parallel_reduce (node_range(&sides), project_sides);
1268 7220 : libmesh_merge_move(ids_to_push, project_sides.new_ids_to_push);
1269 7220 : libmesh_merge_move(ids_to_save, project_sides.new_ids_to_save);
1270 228792 : }
1271 :
1272 243232 : done_saving_ids = true;
1273 243232 : this->send_and_insert_dof_values(ids_to_push, action);
1274 :
1275 : // No ids to save or push this time, but we still use a reduce since
1276 : // nominally ProjectInteriors still has non-const operator()
1277 7220 : ProjectInteriors project_interiors(*this);
1278 479244 : Threads::parallel_reduce (interior_range(&sort_work.interiors),
1279 : project_interiors);
1280 472024 : }
1281 :
1282 :
1283 : template <typename FFunctor, typename GFunctor,
1284 : typename FValue, typename ProjectionAction>
1285 1242981 : GenericProjector<FFunctor, GFunctor, FValue, ProjectionAction>::SubFunctor::SubFunctor
1286 : (GenericProjector & p) :
1287 1152333 : projector(p),
1288 1242981 : action(p.master_action),
1289 1242981 : f(p.master_f),
1290 1242981 : context(p.system, &p.variables, /* allocate_local_matrices= */ false),
1291 1242981 : conts(p.system.n_vars()),
1292 1333629 : field_types(p.system.n_vars()), system(p.system)
1293 : {
1294 : // Loop over all the variables we've been requested to project, to
1295 : // pre-request
1296 2559129 : for (const auto & var : this->projector.variables)
1297 : {
1298 : // FIXME: Need to generalize this to vector-valued elements. [PB]
1299 47893 : FEAbstract * fe = nullptr;
1300 47893 : FEAbstract * side_fe = nullptr;
1301 47893 : FEAbstract * edge_fe = nullptr;
1302 :
1303 : const std::set<unsigned char> & elem_dims =
1304 47893 : context.elem_dimensions();
1305 :
1306 2635009 : for (const auto & dim : elem_dims)
1307 : {
1308 1318861 : context.get_element_fe( var, fe, dim );
1309 1318861 : if (fe->get_fe_type().family == SCALAR)
1310 0 : continue;
1311 1318861 : if (dim > 1)
1312 42746 : context.get_side_fe( var, side_fe, dim );
1313 1318861 : if (dim > 2)
1314 17472 : context.get_edge_fe( var, edge_fe );
1315 :
1316 114395 : fe->get_JxW();
1317 114395 : fe->get_xyz();
1318 114395 : fe->get_JxW();
1319 :
1320 1318861 : fe->request_phi();
1321 1318861 : if (dim > 1)
1322 : {
1323 101602 : side_fe->get_JxW();
1324 101602 : side_fe->get_xyz();
1325 1152207 : side_fe->request_phi();
1326 : }
1327 1318861 : if (dim > 2)
1328 : {
1329 42164 : edge_fe->get_JxW();
1330 42164 : edge_fe->get_xyz();
1331 521062 : edge_fe->request_phi();
1332 : }
1333 :
1334 1318861 : const FEContinuity cont = fe->get_continuity();
1335 1318861 : this->conts[var] = cont;
1336 1318861 : if (cont == C_ONE)
1337 : {
1338 48796 : fe->request_dphi();
1339 48796 : if (dim > 1)
1340 27153 : side_fe->request_dphi();
1341 48796 : if (dim > 2)
1342 3222 : edge_fe->request_dphi();
1343 : }
1344 :
1345 1318861 : this->field_types[var] = FEInterface::field_type(fe->get_fe_type());
1346 : }
1347 : }
1348 :
1349 : // Now initialize any user requested shape functions, xyz vals, etc
1350 263974 : f.init_context(context);
1351 1242981 : }
1352 :
1353 :
1354 : template <typename FFunctor, typename GFunctor,
1355 : typename FValue, typename ProjectionAction>
1356 992031 : GenericProjector<FFunctor, GFunctor, FValue, ProjectionAction>::SubProjector::SubProjector
1357 : (GenericProjector & p) :
1358 992031 : SubFunctor(p)
1359 : {
1360 :
1361 : // Our C1 elements need gradient information
1362 2003233 : for (const auto & var : this->projector.variables)
1363 1087547 : if (this->conts[var] == C_ONE)
1364 : {
1365 1306 : libmesh_assert(p.master_g);
1366 76698 : g = std::make_unique<GFunctor>(*p.master_g);
1367 39002 : g->init_context(context);
1368 1306 : return;
1369 : }
1370 0 : }
1371 :
1372 : template <typename FFunctor, typename GFunctor, typename FValue, typename ProjectionAction>
1373 : template <typename InsertInput,
1374 : typename std::enable_if<
1375 : !std::is_same<typename ProjectionAction::InsertInput, InsertInput>::value,
1376 : int>::type>
1377 : void
1378 0 : GenericProjector<FFunctor, GFunctor, FValue, ProjectionAction>::SubFunctor::insert_id(
1379 : dof_id_type, const InsertInput & , processor_id_type)
1380 : {
1381 0 : libmesh_error_msg("ID insertion should only occur when the provided input matches that "
1382 : "expected by the ProjectionAction");
1383 : }
1384 :
1385 : template <typename FFunctor, typename GFunctor, typename FValue, typename ProjectionAction>
1386 : template <typename InsertInput,
1387 : typename std::enable_if<
1388 : std::is_same<typename ProjectionAction::InsertInput, InsertInput>::value,
1389 : int>::type>
1390 : void
1391 5315035 : GenericProjector<FFunctor, GFunctor, FValue, ProjectionAction>::SubFunctor::insert_id(
1392 : dof_id_type id, const InsertInput & val, processor_id_type pid)
1393 : {
1394 : // We may see invalid ids when expanding a subdomain with a
1395 : // restricted variable
1396 5315035 : if (id == DofObject::invalid_id)
1397 0 : return;
1398 :
1399 448973 : auto iter = new_ids_to_push.find(id);
1400 5315035 : if (iter == new_ids_to_push.end())
1401 5309994 : action.insert(id, val);
1402 : else
1403 : {
1404 170 : libmesh_assert(pid != DofObject::invalid_processor_id);
1405 900 : iter->second = std::make_pair(val, pid);
1406 : }
1407 5315035 : if (!this->projector.done_saving_ids)
1408 : {
1409 208911 : libmesh_assert(!new_ids_to_save.count(id));
1410 2477907 : new_ids_to_save[id] = val;
1411 : }
1412 : }
1413 :
1414 : template <typename FFunctor, typename GFunctor, typename FValue, typename ProjectionAction>
1415 : template <typename InsertInput,
1416 : typename std::enable_if<
1417 : !std::is_same<typename ProjectionAction::InsertInput, InsertInput>::value,
1418 : int>::type>
1419 : void
1420 : GenericProjector<FFunctor, GFunctor, FValue, ProjectionAction>::SubFunctor::insert_ids(
1421 : const std::vector<dof_id_type> &,
1422 : const std::vector<InsertInput> &,
1423 : processor_id_type)
1424 : {
1425 : libmesh_error_msg("ID insertion should only occur when the provided input matches that "
1426 : "expected by the ProjectionAction");
1427 : }
1428 :
1429 : template <typename FFunctor, typename GFunctor, typename FValue, typename ProjectionAction>
1430 : template <typename InsertInput,
1431 : typename std::enable_if<
1432 : std::is_same<typename ProjectionAction::InsertInput, InsertInput>::value,
1433 : int>::type>
1434 : void
1435 3366765 : GenericProjector<FFunctor, GFunctor, FValue, ProjectionAction>::SubFunctor::insert_ids(
1436 : const std::vector<dof_id_type> & ids,
1437 : const std::vector<InsertInput> & vals,
1438 : processor_id_type pid)
1439 : {
1440 382722 : libmesh_assert_equal_to(ids.size(), vals.size());
1441 10206930 : for (auto i : index_range(ids))
1442 : {
1443 6840165 : const dof_id_type id = ids[i];
1444 :
1445 : // We may see invalid ids when expanding a subdomain with a
1446 : // restricted variable
1447 6840165 : if (id == DofObject::invalid_id)
1448 0 : continue;
1449 :
1450 1325909 : const InsertInput & val = vals[i];
1451 :
1452 662955 : auto iter = new_ids_to_push.find(id);
1453 6840165 : if (iter == new_ids_to_push.end())
1454 6837132 : action.insert(id, val);
1455 : else
1456 : {
1457 390 : libmesh_assert(pid != DofObject::invalid_processor_id);
1458 780 : iter->second = std::make_pair(val, pid);
1459 : }
1460 6840165 : if (!this->projector.done_saving_ids)
1461 : {
1462 378456 : libmesh_assert(!new_ids_to_save.count(id));
1463 3291284 : new_ids_to_save[id] = val;
1464 : }
1465 : }
1466 3366765 : }
1467 :
1468 : template <typename FFunctor, typename GFunctor,
1469 : typename FValue, typename ProjectionAction>
1470 26821 : void GenericProjector<FFunctor, GFunctor, FValue, ProjectionAction>::SubFunctor::join
1471 : (const GenericProjector<FFunctor, GFunctor, FValue, ProjectionAction>::SubFunctor & other)
1472 : {
1473 17597 : new_ids_to_push.insert(other.new_ids_to_push.begin(), other.new_ids_to_push.end());
1474 17597 : new_ids_to_save.insert(other.new_ids_to_save.begin(), other.new_ids_to_save.end());
1475 26821 : }
1476 :
1477 :
1478 : template <typename FFunctor, typename GFunctor,
1479 : typename FValue, typename ProjectionAction>
1480 250950 : void GenericProjector<FFunctor, GFunctor, FValue, ProjectionAction>::SortAndCopy::operator()
1481 : (const ConstElemRange & range)
1482 : {
1483 19832 : LOG_SCOPE ("SortAndCopy::operator()","GenericProjector");
1484 :
1485 : // Look at all the elements in the range. Directly copy values from
1486 : // unchanged elements. For other elements, determine sets of
1487 : // vertices, edge nodes, and side nodes to project.
1488 : //
1489 : // As per our other weird nomenclature, "sides" means faces in 3D
1490 : // and edges in 2D, and "edges" gets skipped in 2D
1491 : //
1492 : // This gets tricky in the case of subdomain-restricted
1493 : // variables, for which we might need to do the same projection
1494 : // from different elements when evaluating different variables.
1495 : // We'll keep track of which variables can be projected from which
1496 : // elements.
1497 :
1498 : // If we copy DoFs on a node, keep track of it so we don't bother
1499 : // with any redundant interpolations or projections later.
1500 : //
1501 : // We still need to keep track of *which* variables get copied,
1502 : // since we may be able to copy elements which lack some
1503 : // subdomain-restricted variables.
1504 : //
1505 : // For extra_hanging_dofs FE types, we'll keep track of which
1506 : // variables were copied from vertices, and which from edges/sides,
1507 : // because those types need copies made from *both* on hanging
1508 : // nodes.
1509 19832 : std::unordered_map<const Node *, std::pair<var_set, var_set>> copied_nodes;
1510 :
1511 54819 : const unsigned int sys_num = system.number();
1512 :
1513 : // At hanging nodes for variables with extra hanging dofs we'll need
1514 : // to do projections *separately* from vertex elements and side/edge
1515 : // elements.
1516 19832 : std::vector<unsigned short> extra_hanging_dofs;
1517 9916 : bool all_extra_hanging_dofs = true;
1518 516894 : for (auto v_num : this->projector.variables)
1519 : {
1520 276494 : if (extra_hanging_dofs.size() <= v_num)
1521 265944 : extra_hanging_dofs.resize(v_num+1, false);
1522 287044 : extra_hanging_dofs[v_num] =
1523 276494 : FEInterface::extra_hanging_dofs(system.variable(v_num).type());
1524 :
1525 265944 : if (!extra_hanging_dofs[v_num])
1526 8263 : all_extra_hanging_dofs = false;
1527 : }
1528 :
1529 5257258 : for (const auto & elem : range)
1530 : {
1531 : // If we're doing AMR, we might be able to copy more DoFs than
1532 : // we interpolate or project.
1533 510742 : bool copy_this_elem = false;
1534 :
1535 : #ifdef LIBMESH_ENABLE_AMR
1536 : // If we're projecting from an old grid
1537 510742 : if (f.is_grid_projection())
1538 : {
1539 : // If this element doesn't have an old_dof_object, but it
1540 : // wasn't just refined or just coarsened into activity, then
1541 : // it must be newly added, so the user is responsible for
1542 : // setting the new dofs on it during a grid projection.
1543 3932594 : const DofObject * old_dof_object = elem->get_old_dof_object();
1544 830279 : const Elem::RefinementState h_flag = elem->refinement_flag();
1545 830279 : const Elem::RefinementState p_flag = elem->p_refinement_flag();
1546 4347734 : if (!old_dof_object &&
1547 3517468 : h_flag != Elem::JUST_REFINED &&
1548 : h_flag != Elem::JUST_COARSENED)
1549 156 : continue;
1550 :
1551 : // If this is an unchanged element, just copy everything
1552 3932438 : if (h_flag != Elem::JUST_REFINED &&
1553 261419 : h_flag != Elem::JUST_COARSENED &&
1554 2251497 : p_flag != Elem::JUST_REFINED &&
1555 : p_flag != Elem::JUST_COARSENED)
1556 260796 : copy_this_elem = true;
1557 : else
1558 : {
1559 154330 : bool reinitted = false;
1560 :
1561 308661 : const unsigned int p_level = elem->p_level();
1562 :
1563 : // If this element has a low order monomial which has
1564 : // merely been h refined, copy it.
1565 154330 : const bool copy_possible =
1566 1841383 : p_level == 0 &&
1567 1841382 : h_flag != Elem::JUST_COARSENED &&
1568 : p_flag != Elem::JUST_COARSENED;
1569 :
1570 1843604 : std::vector<typename FFunctor::ValuePushType> Ue(1);
1571 1842116 : std::vector<dof_id_type> elem_dof_ids(1);
1572 :
1573 3428228 : for (auto v_num : this->projector.variables)
1574 : {
1575 1740442 : const Variable & var = system.variable(v_num);
1576 1740442 : if (!var.active_on_subdomain(elem->subdomain_id()))
1577 7368 : continue;
1578 1733074 : const FEType fe_type = var.type();
1579 :
1580 41682 : if (fe_type.family == MONOMIAL &&
1581 1741756 : fe_type.order == CONSTANT &&
1582 : copy_possible)
1583 : {
1584 5670 : if (!reinitted)
1585 : {
1586 468 : reinitted = true;
1587 5670 : context.pre_fe_reinit(system, elem);
1588 : }
1589 :
1590 5670 : f.eval_old_dofs(*elem, fe_type, sys_num, v_num,
1591 : elem_dof_ids, Ue);
1592 :
1593 5670 : action.insert(elem_dof_ids[0], Ue[0]);
1594 : }
1595 : }
1596 13392 : }
1597 : }
1598 : #endif // LIBMESH_ENABLE_AMR
1599 :
1600 5006152 : const int dim = elem->dim();
1601 :
1602 5006152 : const unsigned int n_vertices = elem->n_vertices();
1603 5006152 : const unsigned int n_edges = elem->n_edges();
1604 5006152 : const unsigned int n_nodes = elem->n_nodes();
1605 :
1606 : // In 1-D we already handle our sides as vertices
1607 5006152 : const unsigned int n_sides = (dim > 1) * elem->n_sides();
1608 :
1609 : // What variables are supported on each kind of node on this elem?
1610 1021458 : var_set vertex_vars, edge_vars, side_vars;
1611 :
1612 : // If we have non-vertex nodes, the first is an edge node, but
1613 : // if we're in 2D we'll call that a side node
1614 5006152 : const bool has_edge_nodes = (n_nodes > n_vertices && dim > 2);
1615 :
1616 : // If we have even more nodes, the next is a side node.
1617 510729 : const bool has_side_nodes =
1618 5006152 : (n_nodes > n_vertices + ((dim > 2) * n_edges));
1619 :
1620 : // We may be out of nodes at this point or we have interior
1621 : // nodes which may have DoFs to project too
1622 510729 : const bool has_interior_nodes =
1623 5006152 : (n_nodes > n_vertices + ((dim > 2) * n_edges) + n_sides);
1624 :
1625 10113494 : for (auto v_num : this->projector.variables)
1626 : {
1627 5107342 : const Variable & var = this->projector.system.variable(v_num);
1628 5107342 : if (!var.active_on_subdomain(elem->subdomain_id()))
1629 141930 : continue;
1630 5092096 : const FEType fe_type = var.type();
1631 516277 : const auto & dof_map = this->projector.system.get_dof_map();
1632 5092096 : const auto vg = dof_map.var_group_from_var_number(v_num);
1633 5092096 : const bool add_p_level = dof_map.should_p_refine(vg);
1634 :
1635 : // If we're trying to do projections on an isogeometric
1636 : // analysis mesh, only the finite element nodes constrained
1637 : // by those spline nodes truly have delta function degrees
1638 : // of freedom. We'll have to back out the spline degrees of
1639 : // freedom indirectly later.
1640 5104810 : if (fe_type.family == RATIONAL_BERNSTEIN &&
1641 12714 : elem->type() == NODEELEM)
1642 5435 : continue;
1643 :
1644 5086188 : if (FEInterface::n_dofs_at_node(fe_type, elem, 0, add_p_level))
1645 4365541 : vertex_vars.insert(vertex_vars.end(), v_num);
1646 :
1647 : // The first non-vertex node is always an edge node if those
1648 : // exist. All edge nodes have the same number of DoFs
1649 5086188 : if (has_edge_nodes)
1650 191472 : if (FEInterface::n_dofs_at_node(fe_type, elem, n_vertices, add_p_level))
1651 58013 : edge_vars.insert(edge_vars.end(), v_num);
1652 :
1653 5086188 : if (has_side_nodes)
1654 : {
1655 1065671 : if (dim != 3)
1656 : {
1657 913991 : if (FEInterface::n_dofs_at_node(fe_type, elem, n_vertices, add_p_level))
1658 740996 : side_vars.insert(side_vars.end(), v_num);
1659 : }
1660 : else
1661 : // In 3D, not all face nodes always have the same number of
1662 : // DoFs! We'll loop over all sides to be safe.
1663 2276558 : for (unsigned int n = 0; n != n_nodes; ++n)
1664 2235408 : if (elem->is_face(n))
1665 302116 : if (FEInterface::n_dofs_at_node(fe_type, elem, n, add_p_level))
1666 : {
1667 100722 : side_vars.insert(side_vars.end(), v_num);
1668 9808 : break;
1669 : }
1670 : }
1671 :
1672 5161279 : if (FEInterface::n_dofs_per_elem(fe_type, elem, add_p_level) ||
1673 855091 : (has_interior_nodes &&
1674 855091 : FEInterface::n_dofs_at_node(fe_type, elem, n_nodes-1, add_p_level)))
1675 : {
1676 : #ifdef LIBMESH_ENABLE_AMR
1677 : // We may have already just copied constant monomials,
1678 : // or we may be about to copy the whole element
1679 930387 : if ((f.is_grid_projection() &&
1680 830350 : fe_type.family == MONOMIAL &&
1681 12565 : fe_type.order == CONSTANT &&
1682 9213 : elem->p_level() == 0 &&
1683 9201 : elem->refinement_flag() != Elem::JUST_COARSENED &&
1684 1528 : elem->p_refinement_flag() != Elem::JUST_COARSENED)
1685 899484 : || copy_this_elem
1686 : )
1687 107899 : continue;
1688 : #endif // LIBMESH_ENABLE_AMR
1689 :
1690 : // We need to project any other variables
1691 815079 : if (interiors.empty() || interiors.back() != elem)
1692 794251 : interiors.push_back(elem);
1693 : }
1694 : }
1695 :
1696 : // We'll use a greedy algorithm in most cases: if another
1697 : // element has already claimed some of our DoFs, we'll let it do
1698 : // the work.
1699 28102228 : auto erase_covered_vars = []
1700 5517492 : (const Node * node, var_set & remaining, const ves_multimap & covered)
1701 : {
1702 2758745 : auto covered_range = covered.equal_range(node);
1703 36536873 : for (const auto & v_ent : as_range(covered_range))
1704 17003550 : for (const unsigned int var_covered :
1705 763644 : std::get<2>(v_ent.second))
1706 776232 : remaining.erase(var_covered);
1707 : };
1708 :
1709 28691333 : auto erase_nonhanging_vars = [&extra_hanging_dofs]
1710 5257968 : (const Node * node, var_set & remaining, const ves_multimap & covered)
1711 : {
1712 2628984 : auto covered_range = covered.equal_range(node);
1713 23688834 : for (const auto & v_ent : as_range(covered_range))
1714 8304 : for (const unsigned int var_covered :
1715 257 : std::get<2>(v_ent.second))
1716 4651 : if (!extra_hanging_dofs[var_covered])
1717 273 : remaining.erase(var_covered);
1718 : };
1719 :
1720 20027996 : auto erase_copied_vars = [&copied_nodes, &extra_hanging_dofs]
1721 2929954 : (const Node * node, bool is_vertex, var_set & remaining)
1722 : {
1723 1464976 : auto copying_range = copied_nodes.equal_range(node);
1724 19698856 : for (const auto & v_ent : as_range(copying_range))
1725 : {
1726 12209235 : for (const unsigned int var_covered :
1727 : v_ent.second.first)
1728 6063126 : if (is_vertex || !extra_hanging_dofs[var_covered])
1729 691150 : remaining.erase(var_covered);
1730 :
1731 6366686 : for (const unsigned int var_covered :
1732 : v_ent.second.second)
1733 220577 : if (!is_vertex || !extra_hanging_dofs[var_covered])
1734 24752 : remaining.erase(var_covered);
1735 : }
1736 : };
1737 :
1738 24016286 : for (unsigned int v=0; v != n_vertices; ++v)
1739 : {
1740 20968365 : const Node * node = elem->node_ptr(v);
1741 :
1742 1958229 : auto remaining_vars = vertex_vars;
1743 :
1744 19010134 : erase_covered_vars(node, remaining_vars, vertices);
1745 :
1746 19010134 : if (remaining_vars.empty())
1747 7301739 : continue;
1748 :
1749 10984073 : if (!all_extra_hanging_dofs)
1750 : {
1751 10692724 : erase_nonhanging_vars(node, remaining_vars, edges);
1752 10692724 : if (remaining_vars.empty())
1753 0 : continue;
1754 :
1755 10692724 : erase_nonhanging_vars(node, remaining_vars, sides);
1756 10692724 : if (remaining_vars.empty())
1757 2127 : continue;
1758 : }
1759 :
1760 10981693 : erase_copied_vars(node, true, remaining_vars);
1761 10981693 : if (remaining_vars.empty())
1762 5280538 : continue;
1763 :
1764 4180647 : if (copy_this_elem)
1765 : {
1766 581800 : std::vector<dof_id_type> node_dof_ids;
1767 581800 : std::vector<typename FFunctor::ValuePushType> values;
1768 :
1769 4638798 : for (auto var : remaining_vars)
1770 : {
1771 2371924 : f.eval_old_dofs(*elem, v, var, node_dof_ids, values);
1772 2371924 : insert_ids(node_dof_ids, values, node->processor_id());
1773 : }
1774 2266874 : copied_nodes[node].first.insert(remaining_vars.begin(),
1775 : remaining_vars.end());
1776 2266874 : this->find_dofs_to_send(*node, *elem, v, remaining_vars);
1777 0 : }
1778 : else
1779 : vertices.emplace
1780 2992443 : (node, std::make_tuple(elem, v, std::move(remaining_vars)));
1781 : }
1782 :
1783 5006152 : if (has_edge_nodes)
1784 : {
1785 1491282 : for (unsigned int e=0; e != n_edges; ++e)
1786 : {
1787 1423830 : const Node * node = elem->node_ptr(n_vertices+e);
1788 :
1789 115212 : auto remaining_vars = edge_vars;
1790 :
1791 1308618 : erase_covered_vars(node, remaining_vars, edges);
1792 1308618 : if (remaining_vars.empty())
1793 921995 : continue;
1794 :
1795 300953 : erase_covered_vars(node, remaining_vars, sides);
1796 300953 : if (remaining_vars.empty())
1797 0 : continue;
1798 :
1799 300953 : if (!all_extra_hanging_dofs)
1800 : {
1801 205736 : erase_nonhanging_vars(node, remaining_vars, vertices);
1802 205736 : if (remaining_vars.empty())
1803 383 : continue;
1804 : }
1805 :
1806 300570 : erase_copied_vars(node, false, remaining_vars);
1807 300570 : if (remaining_vars.empty())
1808 16711 : continue;
1809 :
1810 127774 : if (copy_this_elem)
1811 : {
1812 7764 : std::vector<dof_id_type> edge_dof_ids;
1813 7764 : std::vector<typename FFunctor::ValuePushType> values;
1814 :
1815 52764 : for (auto var : remaining_vars)
1816 : {
1817 26382 : f.eval_old_dofs(*elem, n_vertices+e, var, edge_dof_ids, values);
1818 26382 : insert_ids(edge_dof_ids, values, node->processor_id());
1819 : }
1820 :
1821 26382 : copied_nodes[node].second.insert(remaining_vars.begin(),
1822 : remaining_vars.end());
1823 26382 : this->find_dofs_to_send(*node, *elem, n_vertices+e, remaining_vars);
1824 0 : }
1825 : else
1826 : edges.emplace
1827 273901 : (node, std::make_tuple(elem, e, std::move(remaining_vars)));
1828 : }
1829 : }
1830 :
1831 5006152 : if (has_side_nodes)
1832 : {
1833 4937059 : for (unsigned int side=0; side != n_sides; ++side)
1834 : {
1835 3933680 : const Node * node = nullptr;
1836 3933680 : unsigned short node_num = n_vertices+(dim>2)*n_edges+side;
1837 3933680 : if (dim != 3)
1838 3587737 : node = elem->node_ptr(node_num);
1839 : else
1840 : {
1841 : // In 3D only some sides may have nodes
1842 9840228 : for (unsigned int n = 0; n != n_nodes; ++n)
1843 : {
1844 9822816 : if (!elem->is_face(n))
1845 7371150 : continue;
1846 :
1847 1734426 : if (elem->is_node_on_side(n, side))
1848 : {
1849 152038 : node_num = n;
1850 617544 : node = elem->node_ptr(node_num);
1851 617544 : break;
1852 : }
1853 : }
1854 : }
1855 :
1856 3933680 : if (!node)
1857 1804938 : continue;
1858 :
1859 342431 : auto remaining_vars = side_vars;
1860 :
1861 3916268 : erase_covered_vars(node, remaining_vars, edges);
1862 3916268 : if (remaining_vars.empty())
1863 320913 : continue;
1864 :
1865 3566255 : erase_covered_vars(node, remaining_vars, sides);
1866 3566255 : if (remaining_vars.empty())
1867 1183405 : continue;
1868 :
1869 2271299 : if (!all_extra_hanging_dofs)
1870 : {
1871 2093724 : erase_nonhanging_vars(node, remaining_vars, vertices);
1872 2093724 : if (remaining_vars.empty())
1873 815 : continue;
1874 : }
1875 :
1876 2270484 : erase_copied_vars(node, false, remaining_vars);
1877 2270484 : if (remaining_vars.empty())
1878 125453 : continue;
1879 :
1880 1823479 : if (copy_this_elem)
1881 : {
1882 67314 : std::vector<dof_id_type> side_dof_ids;
1883 67314 : std::vector<typename FFunctor::ValuePushType> values;
1884 :
1885 633086 : for (auto var : remaining_vars)
1886 : {
1887 327417 : f.eval_old_dofs(*elem, node_num, var, side_dof_ids, values);
1888 327417 : insert_ids(side_dof_ids, values, node->processor_id());
1889 : }
1890 :
1891 305669 : copied_nodes[node].second.insert(remaining_vars.begin(),
1892 : remaining_vars.end());
1893 305669 : this->find_dofs_to_send(*node, *elem, node_num, remaining_vars);
1894 0 : }
1895 : else
1896 : sides.emplace
1897 1974907 : (node, std::make_tuple(elem, side, std::move(remaining_vars)));
1898 : }
1899 : }
1900 :
1901 : // Elements with elemental dofs might need those copied too.
1902 4028041 : if (copy_this_elem)
1903 : {
1904 521592 : std::vector<typename FFunctor::ValuePushType> U;
1905 521592 : std::vector<dof_id_type> dof_ids;
1906 :
1907 4531166 : for (auto v_num : this->projector.variables)
1908 : {
1909 2286514 : const Variable & var = system.variable(v_num);
1910 2286514 : if (!var.active_on_subdomain(elem->subdomain_id()))
1911 4818 : continue;
1912 2281696 : FEType fe_type = var.type();
1913 :
1914 2281696 : f.eval_old_dofs(*elem, fe_type, sys_num, v_num,
1915 : dof_ids, U);
1916 4039348 : action.insert(dof_ids, U);
1917 :
1918 2281696 : if (has_interior_nodes)
1919 : {
1920 128141 : f.eval_old_dofs(*elem, n_nodes-1, v_num, dof_ids, U);
1921 242770 : action.insert(dof_ids, U);
1922 : }
1923 : }
1924 0 : }
1925 : }
1926 250950 : }
1927 :
1928 :
1929 : template <typename FFunctor, typename GFunctor,
1930 : typename FValue, typename ProjectionAction>
1931 7718 : void GenericProjector<FFunctor, GFunctor, FValue, ProjectionAction>::SortAndCopy::join
1932 : (const GenericProjector<FFunctor, GFunctor, FValue, ProjectionAction>::SortAndCopy & other)
1933 : {
1934 38220 : auto merge_multimaps = [](ves_multimap & self, const ves_multimap & other_mm)
1935 : {
1936 657167 : for (const auto & pair : other_mm)
1937 : {
1938 634013 : const Node * node = pair.first;
1939 433046 : auto remaining_vars = std::get<2>(pair.second);
1940 :
1941 216523 : auto my_range = self.equal_range(node);
1942 705066 : for (const auto & v_ent : as_range(my_range))
1943 143363 : for (const unsigned int var_covered :
1944 23924 : std::get<2>(v_ent.second))
1945 24343 : remaining_vars.erase(var_covered);
1946 :
1947 634013 : if (!remaining_vars.empty())
1948 : self.emplace
1949 756352 : (node, std::make_tuple(std::get<0>(pair.second),
1950 192797 : std::get<1>(pair.second),
1951 192797 : std::move(remaining_vars)));
1952 : }
1953 : };
1954 :
1955 7718 : merge_multimaps(vertices, other.vertices);
1956 7718 : merge_multimaps(edges, other.edges);
1957 7718 : merge_multimaps(sides, other.sides);
1958 :
1959 7718 : std::sort(interiors.begin(), interiors.end());
1960 10414 : std::vector<const Elem *> other_interiors = other.interiors;
1961 7718 : std::sort(other_interiors.begin(), other_interiors.end());
1962 5392 : std::vector<const Elem *> merged_interiors;
1963 7718 : std::set_union(interiors.begin(), interiors.end(),
1964 : other_interiors.begin(), other_interiors.end(),
1965 : std::back_inserter(merged_interiors));
1966 2696 : interiors.swap(merged_interiors);
1967 :
1968 7718 : SubFunctor::join(other);
1969 7718 : }
1970 :
1971 : namespace
1972 : {
1973 : template <typename Output, typename Input>
1974 : typename std::enable_if<ScalarTraits<Input>::value, Output>::type
1975 0 : raw_value(const Input & input, unsigned int /*index*/)
1976 : {
1977 0 : return input;
1978 : }
1979 :
1980 : template <typename Output, template <typename> class WrapperClass, typename T>
1981 : typename std::enable_if<ScalarTraits<T>::value &&
1982 : TensorTools::MathWrapperTraits<WrapperClass<T>>::value,
1983 : Output>::type
1984 16143 : raw_value(const WrapperClass<T> & input, unsigned int index)
1985 : {
1986 233178 : return input.slice(index);
1987 : }
1988 :
1989 : template <typename T>
1990 : typename T::value_type
1991 : grad_component(const T &, unsigned int);
1992 :
1993 : template <typename T>
1994 : typename VectorValue<T>::value_type
1995 2380 : grad_component(const VectorValue<T> & grad, unsigned int component)
1996 : {
1997 4709 : return grad(component);
1998 : }
1999 :
2000 : template <typename T>
2001 : typename TensorValue<T>::value_type
2002 0 : grad_component(const TensorValue<T> & grad, unsigned int component)
2003 : {
2004 0 : libmesh_error_msg("This function should only be called for Hermites. "
2005 : "I don't know how you got here");
2006 : return grad(component, component);
2007 : }
2008 :
2009 :
2010 : }
2011 :
2012 : template <typename FFunctor, typename GFunctor,
2013 : typename FValue, typename ProjectionAction>
2014 252336 : void GenericProjector<FFunctor, GFunctor, FValue, ProjectionAction>::ProjectVertices::operator()
2015 : (const node_range & range)
2016 : {
2017 20748 : LOG_SCOPE ("project_vertices","GenericProjector");
2018 :
2019 252336 : const unsigned int sys_num = system.number();
2020 :
2021 : // Variables with extra hanging dofs can't safely use eval_at_node
2022 : // in as many places as variables without can.
2023 20748 : std::vector<unsigned short> extra_hanging_dofs;
2024 519477 : for (auto v_num : this->projector.variables)
2025 : {
2026 278086 : if (extra_hanging_dofs.size() <= v_num)
2027 267141 : extra_hanging_dofs.resize(v_num+1, false);
2028 267141 : extra_hanging_dofs[v_num] =
2029 278086 : FEInterface::extra_hanging_dofs(system.variable(v_num).type());
2030 : }
2031 :
2032 2945254 : for (const auto & v_pair : range)
2033 : {
2034 2692918 : const Node & vertex = *v_pair.first;
2035 2692918 : const Elem & elem = *std::get<0>(v_pair.second);
2036 2692918 : const unsigned int n = std::get<1>(v_pair.second);
2037 233862 : const var_set & vertex_vars = std::get<2>(v_pair.second);
2038 :
2039 2692918 : context.pre_fe_reinit(system, &elem);
2040 :
2041 2692918 : this->find_dofs_to_send(vertex, elem, n, vertex_vars);
2042 :
2043 : // Look at all the variables we're supposed to interpolate from
2044 : // this element on this vertex
2045 5445179 : for (const auto & var : vertex_vars)
2046 : {
2047 2752261 : const Variable & variable = system.variable(var);
2048 238640 : const FEType & base_fe_type = variable.type();
2049 477280 : const unsigned int var_component =
2050 238640 : system.variable_scalar_number(var, 0);
2051 :
2052 2752261 : if (base_fe_type.family == SCALAR)
2053 0 : continue;
2054 :
2055 2752261 : const FEContinuity cont = this->conts[var];
2056 2752261 : const FEFieldType field_type = this->field_types[var];
2057 :
2058 2752261 : if (cont == DISCONTINUOUS)
2059 : {
2060 0 : libmesh_assert_equal_to(vertex.n_comp(sys_num, var), 0);
2061 : }
2062 2990902 : else if (cont == C_ZERO ||
2063 2513621 : cont == SIDE_DISCONTINUOUS)
2064 : {
2065 2665341 : if (cont == SIDE_DISCONTINUOUS &&
2066 1518 : elem.dim() != 1)
2067 : {
2068 0 : libmesh_assert_equal_to(vertex.n_comp(sys_num, var), 0);
2069 0 : continue;
2070 : }
2071 :
2072 3420288 : const FValue val = f.eval_at_node
2073 1771000 : (context, var_component, /*dim=*/ 0, // Don't care w/C0
2074 2663823 : vertex, extra_hanging_dofs[var], system.time);
2075 :
2076 2663823 : if (field_type == TYPE_VECTOR)
2077 : {
2078 1302 : libmesh_assert_equal_to(vertex.n_comp(sys_num, var), elem.dim());
2079 :
2080 : // We will have a number of nodal value DoFs equal to the elem dim
2081 83664 : for (auto i : make_range(elem.dim()))
2082 : {
2083 61920 : const dof_id_type id = vertex.dof_number(sys_num, var, i);
2084 :
2085 : // Need this conversion so that this method
2086 : // will compile for TYPE_SCALAR instantiations
2087 59112 : const auto insert_val =
2088 6498 : raw_value<typename ProjectionAction::InsertInput>(val, i);
2089 :
2090 61920 : insert_id(id, insert_val, vertex.processor_id());
2091 : }
2092 : }
2093 : else
2094 : {
2095 : // C_ZERO elements have a single nodal value DoF at
2096 : // vertices. We can't assert n_comp==1 here,
2097 : // because if this is a hanging node then it may have
2098 : // more face/edge DoFs, but we don't need to deal with
2099 : // those here.
2100 :
2101 2642079 : const dof_id_type id = vertex.dof_number(sys_num, var, 0);
2102 2642079 : insert_id(id, val, vertex.processor_id());
2103 231149 : }
2104 : }
2105 88438 : else if (cont == C_ONE)
2106 : {
2107 7491 : libmesh_assert(vertex.n_comp(sys_num, var));
2108 88438 : const dof_id_type first_id = vertex.dof_number(sys_num, var, 0);
2109 :
2110 : // C_ONE elements have a single nodal value and dim
2111 : // gradient values at vertices, as well as cross
2112 : // gradients for HERMITE. We need to have an element in
2113 : // hand to figure out dim and to have in case this
2114 : // vertex is a new vertex.
2115 88438 : const int dim = elem.dim();
2116 : #ifndef NDEBUG
2117 : // For now all C1 elements at a vertex had better have
2118 : // the same dimension. If anyone hits these asserts let
2119 : // me know; we could probably support a mixed-dimension
2120 : // mesh IFF the 2D elements were all parallel to xy and
2121 : // the 1D elements all parallel to x.
2122 33310 : for (const auto e_id : (*this->projector.nodes_to_elem)[vertex.id()])
2123 : {
2124 25819 : const Elem & e = system.get_mesh().elem_ref(e_id);
2125 25819 : libmesh_assert_equal_to(dim, e.dim());
2126 : }
2127 : #endif
2128 : #ifdef LIBMESH_ENABLE_AMR
2129 7491 : bool is_old_vertex = true;
2130 88438 : if (elem.refinement_flag() == Elem::JUST_REFINED)
2131 : {
2132 3064 : const int i_am_child =
2133 42333 : elem.parent()->which_child_am_i(&elem);
2134 : is_old_vertex =
2135 42333 : elem.parent()->is_vertex_on_parent(i_am_child, n);
2136 : }
2137 : #else
2138 : const bool is_old_vertex = false;
2139 : #endif
2140 :
2141 : // The hermite element vertex shape functions are weird
2142 88438 : if (base_fe_type.family == HERMITE)
2143 : {
2144 98835 : const FValue val =
2145 29398 : f.eval_at_node(context,
2146 : var_component,
2147 : dim,
2148 : vertex,
2149 : extra_hanging_dofs[var],
2150 74720 : system.time);
2151 74720 : insert_id(first_id, val, vertex.processor_id());
2152 :
2153 13147 : typename GFunctor::FunctorValue grad =
2154 69117 : is_old_vertex ?
2155 13839 : g->eval_at_node(context,
2156 : var_component,
2157 : dim,
2158 : vertex,
2159 : extra_hanging_dofs[var],
2160 58022 : system.time) :
2161 15379 : g->eval_at_point(context,
2162 : var_component,
2163 : vertex,
2164 16698 : system.time,
2165 : false);
2166 : // x derivative. Use slice because grad may be a tensor type
2167 74720 : insert_id(first_id+1, grad.slice(0),
2168 : vertex.processor_id());
2169 : #if LIBMESH_DIM > 1
2170 70436 : if (dim > 1 && is_old_vertex && f.is_grid_projection())
2171 : {
2172 5632 : for (int i = 1; i < dim; ++i)
2173 3174 : insert_id(first_id+i+1, grad.slice(i),
2174 : vertex.processor_id());
2175 :
2176 : // We can directly copy everything else too
2177 716 : std::vector<FValue> derivs;
2178 1074 : f.eval_mixed_derivatives
2179 2458 : (context, var_component, dim, vertex, derivs);
2180 5632 : for (auto i : index_range(derivs))
2181 3174 : insert_id(first_id+dim+1+i, derivs[i],
2182 : vertex.processor_id());
2183 0 : }
2184 70464 : else if (dim > 1)
2185 : {
2186 : // We'll finite difference mixed derivatives.
2187 : // This delta_x used to be TOLERANCE*hmin, but
2188 : // the factor of 10 improved the accuracy in
2189 : // some unit test projections
2190 11418 : Real delta_x = TOLERANCE * 10 * elem.hmin();
2191 :
2192 11418 : Point nxminus = elem.point(n),
2193 11418 : nxplus = elem.point(n);
2194 11418 : nxminus(0) -= delta_x;
2195 11418 : nxplus(0) += delta_x;
2196 1660 : typename GFunctor::FunctorValue gxminus =
2197 9076 : g->eval_at_point(context,
2198 : var_component,
2199 : nxminus,
2200 11418 : system.time,
2201 : true);
2202 1660 : typename GFunctor::FunctorValue gxplus =
2203 9076 : g->eval_at_point(context,
2204 : var_component,
2205 : nxplus,
2206 11418 : system.time,
2207 : true);
2208 : // y derivative
2209 11418 : insert_id(first_id+2, grad.slice(1),
2210 : vertex.processor_id());
2211 : // xy derivative
2212 12320 : insert_id(first_id+3,
2213 11418 : (grad_component(gxplus, 1) - grad_component(gxminus, 1)) / 2. / delta_x,
2214 : vertex.processor_id());
2215 :
2216 : #if LIBMESH_DIM > 2
2217 11418 : if (dim > 2)
2218 : {
2219 : // z derivative
2220 864 : insert_id(first_id+4, grad.slice(2),
2221 : vertex.processor_id());
2222 : // xz derivative
2223 936 : insert_id(first_id+5,
2224 864 : (grad_component(gxplus, 2) - grad_component(gxminus, 2)) / 2. / delta_x,
2225 : vertex.processor_id());
2226 :
2227 : // We need new points for yz
2228 864 : Point nyminus = elem.point(n),
2229 864 : nyplus = elem.point(n);
2230 864 : nyminus(1) -= delta_x;
2231 864 : nyplus(1) += delta_x;
2232 72 : typename GFunctor::FunctorValue gyminus =
2233 72 : g->eval_at_point(context,
2234 : var_component,
2235 : nyminus,
2236 864 : system.time,
2237 : true);
2238 72 : typename GFunctor::FunctorValue gyplus =
2239 72 : g->eval_at_point(context,
2240 : var_component,
2241 : nyplus,
2242 864 : system.time,
2243 : true);
2244 : // yz derivative
2245 936 : insert_id(first_id+6,
2246 864 : (grad_component(gyplus, 2) - grad_component(gyminus, 2)) / 2. / delta_x,
2247 : vertex.processor_id());
2248 : // Getting a 2nd order xyz is more tedious
2249 864 : Point nxmym = elem.point(n),
2250 864 : nxmyp = elem.point(n),
2251 864 : nxpym = elem.point(n),
2252 864 : nxpyp = elem.point(n);
2253 864 : nxmym(0) -= delta_x;
2254 864 : nxmym(1) -= delta_x;
2255 864 : nxmyp(0) -= delta_x;
2256 864 : nxmyp(1) += delta_x;
2257 864 : nxpym(0) += delta_x;
2258 864 : nxpym(1) -= delta_x;
2259 864 : nxpyp(0) += delta_x;
2260 864 : nxpyp(1) += delta_x;
2261 72 : typename GFunctor::FunctorValue gxmym =
2262 72 : g->eval_at_point(context,
2263 : var_component,
2264 : nxmym,
2265 864 : system.time,
2266 : true);
2267 72 : typename GFunctor::FunctorValue gxmyp =
2268 72 : g->eval_at_point(context,
2269 : var_component,
2270 : nxmyp,
2271 864 : system.time,
2272 : true);
2273 72 : typename GFunctor::FunctorValue gxpym =
2274 72 : g->eval_at_point(context,
2275 : var_component,
2276 : nxpym,
2277 864 : system.time,
2278 : true);
2279 72 : typename GFunctor::FunctorValue gxpyp =
2280 72 : g->eval_at_point(context,
2281 : var_component,
2282 : nxpyp,
2283 864 : system.time,
2284 : true);
2285 864 : FValue gxzplus = (grad_component(gxpyp, 2) - grad_component(gxmyp, 2))
2286 792 : / 2. / delta_x;
2287 864 : FValue gxzminus = (grad_component(gxpym, 2) - grad_component(gxmym, 2))
2288 792 : / 2. / delta_x;
2289 : // xyz derivative
2290 936 : insert_id(first_id+7,
2291 864 : (gxzplus - gxzminus) / 2. / delta_x,
2292 : vertex.processor_id());
2293 : }
2294 : #endif // LIBMESH_DIM > 2
2295 : }
2296 : #endif // LIBMESH_DIM > 1
2297 : }
2298 : else
2299 : {
2300 : // Currently other C_ONE elements have a single nodal
2301 : // value shape function and nodal gradient component
2302 : // shape functions
2303 918 : libmesh_assert_equal_to(
2304 : FEInterface::n_dofs_at_node(
2305 : base_fe_type,
2306 : &elem,
2307 : elem.get_node_index(&vertex),
2308 : this->projector.system.get_dof_map().should_p_refine(
2309 : this->projector.system.get_dof_map().var_group_from_var_number(var))),
2310 : (unsigned int)(1 + dim));
2311 :
2312 16554 : const FValue val =
2313 11864 : f.eval_at_node(context, var_component, dim,
2314 : vertex, extra_hanging_dofs[var],
2315 13718 : system.time);
2316 13718 : insert_id(first_id, val, vertex.processor_id());
2317 1836 : typename GFunctor::FunctorValue grad =
2318 12886 : is_old_vertex ?
2319 2082 : g->eval_at_node(context, var_component, dim,
2320 : vertex, extra_hanging_dofs[var],
2321 3284 : system.time) :
2322 9710 : g->eval_at_point(context, var_component, vertex,
2323 10434 : system.time, false);
2324 41154 : for (int i=0; i!= dim; ++i)
2325 29272 : insert_id(first_id + i + 1, grad.slice(i),
2326 : vertex.processor_id());
2327 : }
2328 : }
2329 : else
2330 0 : libmesh_error_msg("Unknown continuity " << cont);
2331 : }
2332 : }
2333 252336 : }
2334 :
2335 :
2336 : template <typename FFunctor, typename GFunctor,
2337 : typename FValue, typename ProjectionAction>
2338 244968 : void GenericProjector<FFunctor, GFunctor, FValue, ProjectionAction>::ProjectEdges::operator()
2339 : (const node_range & range)
2340 : {
2341 15598 : LOG_SCOPE ("project_edges","GenericProjector");
2342 :
2343 244968 : const unsigned int sys_num = system.number();
2344 :
2345 489523 : for (const auto & e_pair : range)
2346 : {
2347 244555 : const Elem & elem = *std::get<0>(e_pair.second);
2348 :
2349 : // If this is an unchanged element then we already copied all
2350 : // its dofs
2351 : #ifdef LIBMESH_ENABLE_AMR
2352 103562 : if (f.is_grid_projection() &&
2353 13692 : (elem.refinement_flag() != Elem::JUST_REFINED &&
2354 0 : elem.refinement_flag() != Elem::JUST_COARSENED &&
2355 0 : elem.p_refinement_flag() != Elem::JUST_REFINED &&
2356 0 : elem.p_refinement_flag() != Elem::JUST_COARSENED))
2357 0 : continue;
2358 : #endif // LIBMESH_ENABLE_AMR
2359 :
2360 244555 : const Node & edge_node = *e_pair.first;
2361 244555 : const int dim = elem.dim();
2362 18274 : const var_set & edge_vars = std::get<2>(e_pair.second);
2363 :
2364 244555 : const unsigned short edge_num = std::get<1>(e_pair.second);
2365 244555 : const unsigned short node_num = elem.n_vertices() + edge_num;
2366 244555 : context.edge = edge_num;
2367 244555 : context.pre_fe_reinit(system, &elem);
2368 :
2369 244555 : this->find_dofs_to_send(edge_node, elem, node_num, edge_vars);
2370 :
2371 : // Look at all the variables we're supposed to interpolate from
2372 : // this element on this edge
2373 489110 : for (const auto & var : edge_vars)
2374 : {
2375 244555 : const Variable & variable = system.variable(var);
2376 18274 : const FEType & base_fe_type = variable.type();
2377 36548 : const unsigned int var_component =
2378 18274 : system.variable_scalar_number(var, 0);
2379 :
2380 244555 : if (base_fe_type.family == SCALAR)
2381 131941 : continue;
2382 :
2383 244555 : FEType fe_type = base_fe_type;
2384 :
2385 : // This may be a p refined element
2386 36548 : fe_type.order = fe_type.order + elem.p_level();
2387 :
2388 : // If this is a Lagrange element with DoFs on edges then by
2389 : // convention we interpolate at the node rather than project
2390 : // along the edge.
2391 244555 : if (fe_type.family == LAGRANGE || fe_type.family == LAGRANGE_VEC)
2392 : {
2393 131941 : if (fe_type.order > 1)
2394 : {
2395 : const processor_id_type pid =
2396 20288 : edge_node.processor_id();
2397 136192 : FValue fval = f.eval_at_point
2398 131941 : (context, var_component, edge_node, system.time,
2399 : false);
2400 131941 : if (fe_type.family == LAGRANGE_VEC)
2401 : {
2402 : // We will have a number of nodal value DoFs equal to the elem dim
2403 142392 : for (auto i : make_range(elem.dim()))
2404 : {
2405 : const dof_id_type dof_id =
2406 106794 : edge_node.dof_number(sys_num, var, i);
2407 :
2408 : // Need this conversion so that this method
2409 : // will compile for TYPE_SCALAR instantiations
2410 100350 : const auto insert_val =
2411 13896 : raw_value<typename ProjectionAction::InsertInput>(fval, i);
2412 :
2413 106794 : insert_id(dof_id, insert_val, pid);
2414 : }
2415 : }
2416 : else // We are LAGRANGE
2417 : {
2418 : const dof_id_type dof_id =
2419 96343 : edge_node.dof_number(sys_num, var, 0);
2420 96343 : insert_id(dof_id, fval, pid);
2421 : }
2422 : }
2423 131941 : continue;
2424 111653 : }
2425 :
2426 : #ifdef LIBMESH_ENABLE_AMR
2427 : // If this is a low order monomial element which has merely
2428 : // been h refined then we already copied all its dofs
2429 0 : if (fe_type.family == MONOMIAL &&
2430 0 : fe_type.order == CONSTANT &&
2431 112614 : elem.refinement_flag() != Elem::JUST_COARSENED &&
2432 0 : elem.p_refinement_flag() != Elem::JUST_COARSENED)
2433 0 : continue;
2434 : #endif // LIBMESH_ENABLE_AMR
2435 :
2436 : // FIXME: Need to generalize this to vector-valued elements. [PB]
2437 8130 : FEGenericBase<typename FFunctor::RealType> * fe = nullptr;
2438 8130 : context.get_element_fe( var, fe, dim );
2439 8130 : FEGenericBase<typename FFunctor::RealType> * edge_fe = nullptr;
2440 8130 : context.get_edge_fe( var, edge_fe );
2441 :
2442 : // If we're JUST_COARSENED we'll need a custom
2443 : // evaluation, not just the standard edge FE
2444 112614 : const FEGenericBase<typename FFunctor::RealType> & proj_fe =
2445 : #ifdef LIBMESH_ENABLE_AMR
2446 16260 : (elem.refinement_flag() == Elem::JUST_COARSENED) ?
2447 : *fe :
2448 : #endif
2449 : *edge_fe;
2450 :
2451 : #ifdef LIBMESH_ENABLE_AMR
2452 112614 : if (elem.refinement_flag() == Elem::JUST_COARSENED)
2453 : {
2454 0 : std::vector<Point> fine_points;
2455 :
2456 0 : std::unique_ptr<FEGenericBase<typename FFunctor::RealType>> fine_fe
2457 : (FEGenericBase<typename FFunctor::RealType>::build (dim, base_fe_type));
2458 :
2459 0 : std::unique_ptr<QBase> qrule
2460 : (base_fe_type.default_quadrature_rule(1));
2461 0 : fine_fe->attach_quadrature_rule(qrule.get());
2462 :
2463 0 : const std::vector<Point> & child_xyz =
2464 0 : fine_fe->get_xyz();
2465 :
2466 0 : for (unsigned int c = 0, nc = elem.n_children();
2467 0 : c != nc; ++c)
2468 : {
2469 0 : if (!elem.is_child_on_edge(c, context.edge))
2470 0 : continue;
2471 :
2472 0 : fine_fe->edge_reinit(elem.child_ptr(c), context.edge);
2473 0 : fine_points.insert(fine_points.end(),
2474 : child_xyz.begin(),
2475 : child_xyz.end());
2476 : }
2477 :
2478 0 : std::vector<Point> fine_qp;
2479 0 : FEMap::inverse_map (dim, &elem, fine_points, fine_qp);
2480 :
2481 0 : context.elem_fe_reinit(&fine_qp);
2482 0 : }
2483 : else
2484 : #endif // LIBMESH_ENABLE_AMR
2485 112614 : context.edge_fe_reinit();
2486 :
2487 16260 : const std::vector<dof_id_type> & dof_indices =
2488 96354 : context.get_dof_indices(var);
2489 :
2490 16260 : std::vector<unsigned int> edge_dofs;
2491 112614 : FEInterface::dofs_on_edge(&elem, dim, base_fe_type,
2492 112614 : context.edge, edge_dofs);
2493 :
2494 16260 : this->construct_projection
2495 96354 : (dof_indices, edge_dofs, var_component,
2496 : &edge_node, proj_fe);
2497 : }
2498 : }
2499 244968 : }
2500 :
2501 :
2502 : template <typename FFunctor, typename GFunctor,
2503 : typename FValue, typename ProjectionAction>
2504 248782 : void GenericProjector<FFunctor, GFunctor, FValue, ProjectionAction>::ProjectSides::operator()
2505 : (const node_range & range)
2506 : {
2507 18180 : LOG_SCOPE ("project_sides","GenericProjector");
2508 :
2509 248782 : const unsigned int sys_num = system.number();
2510 :
2511 2058766 : for (const auto & s_pair : range)
2512 : {
2513 1809984 : const Elem & elem = *std::get<0>(s_pair.second);
2514 :
2515 : // If this is an unchanged element then we already copied all
2516 : // its dofs
2517 : #ifdef LIBMESH_ENABLE_AMR
2518 1629965 : if (f.is_grid_projection() &&
2519 342121 : (elem.refinement_flag() != Elem::JUST_REFINED &&
2520 9841 : elem.refinement_flag() != Elem::JUST_COARSENED &&
2521 308 : elem.p_refinement_flag() != Elem::JUST_REFINED &&
2522 20 : elem.p_refinement_flag() != Elem::JUST_COARSENED))
2523 0 : continue;
2524 : #endif // LIBMESH_ENABLE_AMR
2525 :
2526 1809984 : const Node & side_node = *s_pair.first;
2527 1809984 : const int dim = elem.dim();
2528 147471 : const var_set & side_vars = std::get<2>(s_pair.second);
2529 :
2530 1809984 : const unsigned int side_num = std::get<1>(s_pair.second);
2531 1809984 : unsigned short node_num = elem.n_vertices()+side_num;
2532 : // In 3D only some sides may have nodes
2533 1809984 : if (dim == 3)
2534 5091294 : for (auto n : make_range(elem.n_nodes()))
2535 : {
2536 5091294 : if (!elem.is_face(n))
2537 3856492 : continue;
2538 :
2539 904677 : if (elem.is_node_on_side(n, side_num))
2540 : {
2541 322626 : node_num = n;
2542 322626 : break;
2543 : }
2544 : }
2545 :
2546 1809984 : context.side = side_num;
2547 1809984 : context.pre_fe_reinit(system, &elem);
2548 :
2549 1809984 : this->find_dofs_to_send(side_node, elem, node_num, side_vars);
2550 :
2551 : // Look at all the variables we're supposed to interpolate from
2552 : // this element on this side
2553 3669513 : for (const auto & var : side_vars)
2554 : {
2555 1859529 : const Variable & variable = system.variable(var);
2556 151674 : const FEType & base_fe_type = variable.type();
2557 303348 : const unsigned int var_component =
2558 151674 : system.variable_scalar_number(var, 0);
2559 :
2560 1859529 : if (base_fe_type.family == SCALAR)
2561 1493696 : continue;
2562 :
2563 1859529 : FEType fe_type = base_fe_type;
2564 151674 : const auto & dof_map = system.get_dof_map();
2565 1859529 : const bool add_p_level = dof_map.should_p_refine_var(var);
2566 :
2567 : // This may be a p refined element
2568 2011203 : fe_type.order = fe_type.order + add_p_level*elem.p_level();
2569 :
2570 : // If this is a Lagrange element with DoFs on sides then by
2571 : // convention we interpolate at the node rather than project
2572 : // along the side.
2573 1859529 : if (fe_type.family == LAGRANGE || fe_type.family == LAGRANGE_VEC)
2574 : {
2575 : // If order==1 then we only have DoFs on vertices, so
2576 : // skip all this.
2577 : // If order>1 ... we might still have something tricky
2578 : // like order==2 PRISM20, where some side nodes have
2579 : // DoFs but others don't. Gotta check.
2580 2987392 : if (fe_type.order > 1 &&
2581 1493696 : side_node.n_comp(sys_num, var))
2582 : {
2583 : const processor_id_type pid =
2584 245854 : side_node.processor_id();
2585 1945613 : FValue fval = f.eval_at_point
2586 1488917 : (context, var_component, side_node, system.time,
2587 : false);
2588 :
2589 1488917 : if (fe_type.family == LAGRANGE_VEC)
2590 : {
2591 : // We will have a number of nodal value DoFs equal to the elem dim
2592 84492 : for (auto i : make_range(elem.dim()))
2593 : {
2594 62736 : const dof_id_type dof_id = side_node.dof_number(sys_num, var, i);
2595 :
2596 : // Need this conversion so that this method
2597 : // will compile for TYPE_SCALAR instantiations
2598 58188 : const auto insert_val =
2599 9405 : raw_value<typename ProjectionAction::InsertInput>(fval, i);
2600 :
2601 62736 : insert_id(dof_id, insert_val, pid);
2602 : }
2603 : }
2604 : else // We are LAGRANGE
2605 : {
2606 : const dof_id_type dof_id =
2607 1467161 : side_node.dof_number(sys_num, var, 0);
2608 1467161 : insert_id(dof_id, fval, pid);
2609 : }
2610 : }
2611 1493696 : continue;
2612 1247110 : }
2613 :
2614 : #ifdef LIBMESH_ENABLE_AMR
2615 : // If this is a low order monomial element which has merely
2616 : // been h refined then we already copied all its dofs
2617 0 : if (fe_type.family == MONOMIAL &&
2618 0 : fe_type.order == CONSTANT &&
2619 365833 : elem.refinement_flag() != Elem::JUST_COARSENED &&
2620 0 : elem.p_refinement_flag() != Elem::JUST_COARSENED)
2621 0 : continue;
2622 : #endif // LIBMESH_ENABLE_AMR
2623 :
2624 : // FIXME: Need to generalize this to vector-valued elements. [PB]
2625 28381 : FEGenericBase<typename FFunctor::RealType> * fe = nullptr;
2626 365833 : context.get_element_fe( var, fe, dim );
2627 337452 : FEGenericBase<typename FFunctor::RealType> * side_fe = nullptr;
2628 337452 : context.get_side_fe( var, side_fe );
2629 :
2630 : // If we're JUST_COARSENED we'll need a custom
2631 : // evaluation, not just the standard side FE
2632 365833 : const FEGenericBase<typename FFunctor::RealType> & proj_fe =
2633 : #ifdef LIBMESH_ENABLE_AMR
2634 56762 : (elem.refinement_flag() == Elem::JUST_COARSENED) ?
2635 : *fe :
2636 : #endif
2637 : *side_fe;
2638 :
2639 : #ifdef LIBMESH_ENABLE_AMR
2640 365833 : if (elem.refinement_flag() == Elem::JUST_COARSENED)
2641 : {
2642 0 : std::vector<Point> fine_points;
2643 :
2644 0 : std::unique_ptr<FEGenericBase<typename FFunctor::RealType>> fine_fe
2645 : (FEGenericBase<typename FFunctor::RealType>::build (dim, base_fe_type));
2646 :
2647 0 : std::unique_ptr<QBase> qrule
2648 : (base_fe_type.default_quadrature_rule(1));
2649 0 : fine_fe->attach_quadrature_rule(qrule.get());
2650 :
2651 0 : const std::vector<Point> & child_xyz =
2652 0 : fine_fe->get_xyz();
2653 :
2654 0 : for (unsigned int c = 0, nc = elem.n_children();
2655 0 : c != nc; ++c)
2656 : {
2657 0 : if (!elem.is_child_on_side(c, context.side))
2658 0 : continue;
2659 :
2660 0 : fine_fe->reinit(elem.child_ptr(c), context.side);
2661 0 : fine_points.insert(fine_points.end(),
2662 : child_xyz.begin(),
2663 : child_xyz.end());
2664 : }
2665 :
2666 0 : std::vector<Point> fine_qp;
2667 0 : FEMap::inverse_map (dim, &elem, fine_points, fine_qp);
2668 :
2669 0 : context.elem_fe_reinit(&fine_qp);
2670 0 : }
2671 : else
2672 : #endif // LIBMESH_ENABLE_AMR
2673 365833 : context.side_fe_reinit();
2674 :
2675 56762 : const std::vector<dof_id_type> & dof_indices =
2676 309071 : context.get_dof_indices(var);
2677 :
2678 56762 : std::vector<unsigned int> side_dofs;
2679 422595 : FEInterface::dofs_on_side(&elem, dim, base_fe_type,
2680 365833 : context.side, side_dofs, add_p_level);
2681 :
2682 56762 : this->construct_projection
2683 309071 : (dof_indices, side_dofs, var_component,
2684 : &side_node, proj_fe);
2685 : }
2686 : }
2687 248782 : }
2688 :
2689 :
2690 : template <typename FFunctor, typename GFunctor,
2691 : typename FValue, typename ProjectionAction>
2692 245945 : void GenericProjector<FFunctor, GFunctor, FValue, ProjectionAction>::ProjectInteriors::operator()
2693 : (const interior_range & range)
2694 : {
2695 16290 : LOG_SCOPE ("project_interiors","GenericProjector");
2696 :
2697 245945 : const unsigned int sys_num = system.number();
2698 :
2699 : // Iterate over all dof-bearing element interiors in the range
2700 1040196 : for (const auto & elem : range)
2701 : {
2702 794251 : unsigned char dim = cast_int<unsigned char>(elem->dim());
2703 :
2704 794251 : context.pre_fe_reinit(system, elem);
2705 :
2706 : // Loop over all the variables we've been requested to project, to
2707 : // do the projection
2708 1637822 : for (const auto & var : this->projector.variables)
2709 : {
2710 843571 : const Variable & variable = system.variable(var);
2711 :
2712 843571 : if (!variable.active_on_subdomain(elem->subdomain_id()))
2713 679758 : continue;
2714 :
2715 71012 : const FEType & base_fe_type = variable.type();
2716 :
2717 840067 : if (base_fe_type.family == SCALAR)
2718 0 : continue;
2719 :
2720 71012 : FEGenericBase<typename FFunctor::RealType> * fe = nullptr;
2721 71012 : context.get_element_fe( var, fe, dim );
2722 :
2723 840067 : FEType fe_type = base_fe_type;
2724 71012 : const auto & dof_map = system.get_dof_map();
2725 840067 : const auto vg = dof_map.var_group_from_var_number(var);
2726 840067 : const bool add_p_level = dof_map.should_p_refine(vg);
2727 :
2728 : // This may be a p refined element
2729 911082 : fe_type.order = fe_type.order + add_p_level * elem->p_level();
2730 :
2731 213039 : const unsigned int var_component =
2732 840067 : system.variable_scalar_number(var, 0);
2733 :
2734 : // If this is a Lagrange element with interior DoFs then by
2735 : // convention we interpolate at the interior node rather
2736 : // than project along the interior.
2737 840067 : if (fe_type.family == LAGRANGE || fe_type.family == LAGRANGE_VEC)
2738 : {
2739 676254 : if (fe_type.order > 1)
2740 : {
2741 654476 : const unsigned int first_interior_node =
2742 654476 : (elem->n_vertices() +
2743 654476 : ((elem->dim() > 2) * elem->n_edges()) +
2744 654476 : ((elem->dim() > 1) * elem->n_sides()));
2745 654476 : const unsigned int n_nodes = elem->n_nodes();
2746 :
2747 : // < vs != is important here for HEX20, QUAD8!
2748 1308952 : for (unsigned int n = first_interior_node; n < n_nodes; ++n)
2749 : {
2750 654476 : const Node & interior_node = elem->node_ref(n);
2751 :
2752 867973 : FValue fval = f.eval_at_point
2753 573125 : (context, var_component, interior_node,
2754 654476 : system.time, false);
2755 : const processor_id_type pid =
2756 110862 : interior_node.processor_id();
2757 :
2758 654476 : if (fe_type.family == LAGRANGE_VEC)
2759 : {
2760 : // We will have a number of nodal value DoFs equal to the elem dim
2761 2448 : for (unsigned int i = 0; i < elem->dim(); ++i)
2762 : {
2763 1728 : const dof_id_type dof_id = interior_node.dof_number(sys_num, var, i);
2764 :
2765 : // Need this conversion so that this method
2766 : // will compile for TYPE_SCALAR instantiations
2767 1584 : const auto insert_val =
2768 288 : raw_value<typename ProjectionAction::InsertInput>(fval, i);
2769 :
2770 1728 : insert_id(dof_id, insert_val, pid);
2771 : }
2772 : }
2773 : else // We are LAGRANGE
2774 : {
2775 110742 : const dof_id_type dof_id =
2776 543014 : interior_node.dof_number(sys_num, var, 0);
2777 653756 : insert_id(dof_id, fval, pid);
2778 : }
2779 : }
2780 : }
2781 676254 : continue;
2782 561518 : }
2783 :
2784 : #ifdef LIBMESH_ENABLE_AMR
2785 163813 : if (elem->refinement_flag() == Elem::JUST_COARSENED)
2786 : {
2787 48 : std::vector<Point> fine_points;
2788 :
2789 312 : std::unique_ptr<FEGenericBase<typename FFunctor::RealType>> fine_fe
2790 : (FEGenericBase<typename FFunctor::RealType>::build (dim, base_fe_type));
2791 :
2792 312 : std::unique_ptr<QBase> qrule
2793 : (base_fe_type.default_quadrature_rule(dim));
2794 312 : fine_fe->attach_quadrature_rule(qrule.get());
2795 :
2796 48 : const std::vector<Point> & child_xyz =
2797 48 : fine_fe->get_xyz();
2798 :
2799 2616 : for (auto & child : elem->child_ref_range())
2800 : {
2801 2304 : fine_fe->reinit(&child);
2802 2496 : fine_points.insert(fine_points.end(),
2803 : child_xyz.begin(),
2804 : child_xyz.end());
2805 : }
2806 :
2807 48 : std::vector<Point> fine_qp;
2808 288 : FEMap::inverse_map (dim, elem, fine_points, fine_qp);
2809 :
2810 288 : context.elem_fe_reinit(&fine_qp);
2811 240 : }
2812 : else
2813 : #endif // LIBMESH_ENABLE_AMR
2814 163525 : context.elem_fe_reinit();
2815 :
2816 27291 : const std::vector<dof_id_type> & dof_indices =
2817 136522 : context.get_dof_indices(var);
2818 :
2819 : const unsigned int n_dofs =
2820 27291 : cast_int<unsigned int>(dof_indices.size());
2821 :
2822 177457 : std::vector<unsigned int> all_dofs(n_dofs);
2823 13644 : std::iota(all_dofs.begin(), all_dofs.end(), 0);
2824 :
2825 27291 : this->construct_projection
2826 136522 : (dof_indices, all_dofs, var_component,
2827 : nullptr, *fe);
2828 : } // end variables loop
2829 : } // end elements loop
2830 245945 : }
2831 :
2832 :
2833 :
2834 : template <typename FFunctor, typename GFunctor,
2835 : typename FValue, typename ProjectionAction>
2836 : void
2837 7346382 : GenericProjector<FFunctor, GFunctor, FValue,
2838 : ProjectionAction>::SubFunctor::find_dofs_to_send
2839 : (const Node & node, const Elem & elem, unsigned short node_num, const var_set & vars)
2840 : {
2841 728046 : libmesh_assert (&node == elem.node_ptr(node_num));
2842 :
2843 : // Any ghosted node in our range might have an owner who needs our
2844 : // data
2845 1456093 : const processor_id_type owner = node.processor_id();
2846 8074429 : if (owner != system.processor_id())
2847 : {
2848 66080 : const MeshBase & mesh = system.get_mesh();
2849 33034 : const DofMap & dof_map = system.get_dof_map();
2850 :
2851 : // But let's check and see if we can be certain the owner can
2852 : // compute any or all of its own dof coefficients on that node.
2853 66068 : std::vector<dof_id_type> node_dof_ids, patch_dof_ids;
2854 1521170 : for (const auto & var : vars)
2855 : {
2856 789374 : const Variable & variable = system.variable(var);
2857 :
2858 789374 : if (!variable.active_on_subdomain(elem.subdomain_id()))
2859 0 : continue;
2860 :
2861 789374 : dof_map.dof_indices(elem, node_num, node_dof_ids, var);
2862 : }
2863 33034 : libmesh_assert(std::is_sorted(node_dof_ids.begin(),
2864 : node_dof_ids.end()));
2865 66068 : const std::vector<dof_id_type> & patch =
2866 731796 : (*this->projector.nodes_to_elem)[node.id()];
2867 1887309 : for (const auto & elem_id : patch)
2868 : {
2869 1877631 : const Elem & patch_elem = mesh.elem_ref(elem_id);
2870 1074436 : if (!patch_elem.active() || owner != patch_elem.processor_id())
2871 1139980 : continue;
2872 737651 : dof_map.dof_indices(&patch_elem, patch_dof_ids);
2873 737651 : std::sort(patch_dof_ids.begin(), patch_dof_ids.end());
2874 :
2875 : // Remove any node_dof_ids that we see can be calculated on
2876 : // this element
2877 771165 : std::vector<dof_id_type> diff_ids(node_dof_ids.size());
2878 737651 : auto it = std::set_difference(node_dof_ids.begin(), node_dof_ids.end(),
2879 : patch_dof_ids.begin(), patch_dof_ids.end(), diff_ids.begin());
2880 737651 : diff_ids.resize(it-diff_ids.begin());
2881 33502 : node_dof_ids.swap(diff_ids);
2882 737651 : if (node_dof_ids.empty())
2883 32402 : break;
2884 : }
2885 :
2886 : // Give ids_to_push default invalid pid: not yet computed
2887 747608 : for (auto id : node_dof_ids)
2888 15812 : new_ids_to_push[id].second = DofObject::invalid_processor_id;
2889 : }
2890 7346382 : }
2891 :
2892 :
2893 :
2894 : template <typename FFunctor, typename GFunctor,
2895 : typename FValue, typename ProjectionAction>
2896 : template <typename Value>
2897 : void
2898 729696 : GenericProjector<FFunctor, GFunctor, FValue,
2899 : ProjectionAction>::send_and_insert_dof_values
2900 : (std::unordered_map<dof_id_type, std::pair<Value, processor_id_type>> & ids_to_push,
2901 : ProjectionAction & action) const
2902 : {
2903 43320 : LOG_SCOPE ("send_and_insert_dof_values", "GenericProjector");
2904 :
2905 : // See if we calculated any ids that need to be pushed; get them
2906 : // ready to push.
2907 : std::unordered_map<processor_id_type, std::vector<dof_id_type>>
2908 43320 : pushed_dof_ids, received_dof_ids;
2909 : std::unordered_map<processor_id_type, std::vector<typename TypeToSend<Value>::type>>
2910 43320 : pushed_dof_values, received_dof_values;
2911 772011 : for (auto & id_val_pid : ids_to_push)
2912 : {
2913 42315 : processor_id_type pid = id_val_pid.second.second;
2914 42315 : if (pid != DofObject::invalid_processor_id)
2915 : {
2916 19641 : pushed_dof_ids[pid].push_back(id_val_pid.first);
2917 19641 : pushed_dof_values[pid].push_back(convert_to_send(id_val_pid.second.first));
2918 : }
2919 : }
2920 :
2921 : // If and when we get ids pushed to us, act on them if we have
2922 : // corresponding values or save them if not
2923 729828 : auto ids_action_functor =
2924 5166 : [&action, &received_dof_ids, &received_dof_values]
2925 : (processor_id_type pid,
2926 132 : const std::vector<dof_id_type> & data)
2927 : {
2928 66 : auto iter = received_dof_values.find(pid);
2929 5298 : if (iter == received_dof_values.end())
2930 : {
2931 66 : libmesh_assert(received_dof_ids.find(pid) ==
2932 : received_dof_ids.end());
2933 5298 : received_dof_ids[pid] = data;
2934 : }
2935 : else
2936 : {
2937 0 : auto & vals = iter->second;
2938 0 : std::size_t vals_size = vals.size();
2939 0 : libmesh_assert_equal_to(vals_size, data.size());
2940 0 : for (std::size_t i=0; i != vals_size; ++i)
2941 : {
2942 0 : Value val;
2943 0 : convert_from_receive(vals[i], val);
2944 0 : action.insert(data[i], val);
2945 : }
2946 0 : received_dof_values.erase(iter);
2947 : }
2948 : };
2949 :
2950 : // If and when we get values pushed to us, act on them if we have
2951 : // corresponding ids or save them if not
2952 729828 : auto values_action_functor =
2953 27843 : [&action, &received_dof_ids, &received_dof_values]
2954 : (processor_id_type pid,
2955 132 : const std::vector<typename TypeToSend<Value>::type> & data)
2956 : {
2957 66 : auto iter = received_dof_ids.find(pid);
2958 5298 : if (iter == received_dof_ids.end())
2959 : {
2960 : // We get no more than 1 values vector from anywhere
2961 0 : libmesh_assert(received_dof_values.find(pid) ==
2962 : received_dof_values.end());
2963 0 : received_dof_values[pid] = data;
2964 : }
2965 : else
2966 : {
2967 66 : auto & ids = iter->second;
2968 132 : std::size_t ids_size = ids.size();
2969 66 : libmesh_assert_equal_to(ids_size, data.size());
2970 24939 : for (std::size_t i=0; i != ids_size; ++i)
2971 : {
2972 0 : Value val;
2973 2130 : convert_from_receive(data[i], val);
2974 20706 : action.insert(ids[i], val);
2975 : }
2976 66 : received_dof_ids.erase(iter);
2977 : }
2978 : };
2979 :
2980 : Parallel::push_parallel_vector_data
2981 729696 : (system.comm(), pushed_dof_ids, ids_action_functor);
2982 :
2983 : Parallel::push_parallel_vector_data
2984 729696 : (system.comm(), pushed_dof_values, values_action_functor);
2985 :
2986 : // At this point we shouldn't have any unprocessed data left
2987 21660 : libmesh_assert(received_dof_ids.empty());
2988 21660 : libmesh_assert(received_dof_values.empty());
2989 :
2990 729696 : }
2991 :
2992 :
2993 :
2994 : template <typename FFunctor, typename GFunctor,
2995 : typename FValue, typename ProjectionAction>
2996 : void
2997 642260 : GenericProjector<FFunctor, GFunctor, FValue, ProjectionAction>::SubProjector::construct_projection
2998 : (const std::vector<dof_id_type> & dof_indices_var,
2999 : const std::vector<unsigned int> & involved_dofs,
3000 : unsigned int var_component,
3001 : const Node * node,
3002 : const FEGenericBase<typename FFunctor::RealType> & fe)
3003 : {
3004 642260 : const auto & JxW = fe.get_JxW();
3005 50155 : const auto & phi = fe.get_phi();
3006 50155 : const std::vector<std::vector<typename FEGenericBase<typename FFunctor::RealType>::OutputGradient>> * dphi = nullptr;
3007 147795 : const std::vector<Point> & xyz_values = fe.get_xyz();
3008 642260 : const FEContinuity cont = fe.get_continuity();
3009 50155 : const std::unordered_map<dof_id_type, typename FFunctor::ValuePushType> & ids_to_save =
3010 642260 : this->projector.ids_to_save;
3011 :
3012 642260 : if (cont == C_ONE)
3013 2810 : dphi = &(fe.get_dphi());
3014 :
3015 : const unsigned int n_involved_dofs =
3016 100313 : cast_int<unsigned int>(involved_dofs.size());
3017 :
3018 50155 : std::vector<dof_id_type> free_dof_ids;
3019 592102 : DenseVector<typename FFunctor::ValuePushType> Uinvolved(n_involved_dofs);
3020 642260 : std::vector<char> dof_is_fixed(n_involved_dofs, false); // bools
3021 :
3022 6568242 : for (auto i : make_range(n_involved_dofs))
3023 : {
3024 6871356 : const dof_id_type id = dof_indices_var[involved_dofs[i]];
3025 472676 : auto iter = ids_to_save.find(id);
3026 5925982 : if (iter == ids_to_save.end())
3027 3991100 : free_dof_ids.push_back(id);
3028 : else
3029 : {
3030 1934882 : dof_is_fixed[i] = true;
3031 2089061 : Uinvolved(i) = iter->second;
3032 : }
3033 : }
3034 :
3035 642260 : const unsigned int free_dofs = free_dof_ids.size();
3036 :
3037 : // There may be nothing to project
3038 642260 : if (!free_dofs)
3039 118 : return;
3040 :
3041 : // The element matrix and RHS for projections.
3042 : // Note that Ke is always real-valued, whereas
3043 : // Fe may be complex valued if complex number
3044 : // support is enabled
3045 741119 : DenseMatrix<Real> Ke(free_dofs, free_dofs);
3046 641039 : DenseVector<typename FFunctor::ValuePushType> Fe(free_dofs);
3047 : // The new degree of freedom coefficients to solve for
3048 641039 : DenseVector<typename FFunctor::ValuePushType> Ufree(free_dofs);
3049 :
3050 : const unsigned int n_qp =
3051 100077 : cast_int<unsigned int>(xyz_values.size());
3052 :
3053 : // Loop over the quadrature points
3054 8573860 : for (unsigned int qp=0; qp<n_qp; qp++)
3055 : {
3056 : // solution at the quadrature point
3057 7932818 : FValue fineval = f.eval_at_point(context,
3058 : var_component,
3059 : xyz_values[qp],
3060 7932818 : system.time,
3061 : false);
3062 : // solution grad at the quadrature point
3063 638341 : typename GFunctor::FunctorValue finegrad;
3064 7932818 : if (cont == C_ONE)
3065 159130 : finegrad = g->eval_at_point(context,
3066 : var_component,
3067 : xyz_values[qp],
3068 147430 : system.time,
3069 : false);
3070 :
3071 : // Form edge projection matrix
3072 188517696 : for (unsigned int sidei=0, freei=0;
3073 189776932 : sidei != n_involved_dofs; ++sidei)
3074 : {
3075 181844114 : unsigned int i = involved_dofs[sidei];
3076 : // fixed DoFs aren't test functions
3077 196764971 : if (dof_is_fixed[sidei])
3078 49210600 : continue;
3079 4896386950 : for (unsigned int sidej=0, freej=0;
3080 5014099588 : sidej != n_involved_dofs; ++sidej)
3081 : {
3082 4885867566 : unsigned int j = involved_dofs[sidej];
3083 5290864547 : if (dof_is_fixed[sidej])
3084 545975656 : Fe(freei) -= phi[i][qp] * phi[j][qp] *
3085 364952224 : JxW[qp] * Uinvolved(sidej);
3086 : else
3087 5645509480 : Ke(freei,freej) += phi[i][qp] *
3088 5645509383 : phi[j][qp] * JxW[qp];
3089 4885867566 : if (cont == C_ONE)
3090 : {
3091 1306973 : if (dof_is_fixed[sidej])
3092 1229476 : Fe(freei) -= ( TensorTools::inner_product((*dphi)[i][qp],
3093 1150160 : (*dphi)[j][qp]) ) *
3094 992852 : JxW[qp] * Uinvolved(sidej);
3095 : else
3096 271712 : Ke(freei,freej) += ( TensorTools::inner_product((*dphi)[i][qp],
3097 236053 : (*dphi)[j][qp]) )
3098 236053 : * JxW[qp];
3099 : }
3100 5290864547 : if (!dof_is_fixed[sidej])
3101 4521144894 : freej++;
3102 : }
3103 159790117 : Fe(freei) += phi[i][qp] * fineval * JxW[qp];
3104 128232022 : if (cont == C_ONE)
3105 175400 : Fe(freei) += (TensorTools::inner_product(finegrad,
3106 188552 : (*dphi)[i][qp]) ) *
3107 13000 : JxW[qp];
3108 128232022 : freei++;
3109 : }
3110 : }
3111 :
3112 641042 : Ke.cholesky_solve(Fe, Ufree);
3113 :
3114 : // Transfer new edge solutions to element
3115 641042 : const processor_id_type pid = node ?
3116 72834 : node->processor_id() : DofObject::invalid_processor_id;
3117 641042 : insert_ids(free_dof_ids, Ufree.get_values(), pid);
3118 540965 : }
3119 :
3120 :
3121 : } // namespace libMesh
3122 :
3123 : #endif // GENERIC_PROJECTOR_H
|
