Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 :
20 : #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 1260 : const typename TypeToSend<T>::type convert_to_send (const T& in)
65 23340 : { return in; }
66 :
67 : template <typename SendT, typename T>
68 1260 : void convert_from_receive (SendT & received, T & converted)
69 23340 : { 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 256003 : 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 240835 : system(system_in),
129 240835 : master_f(f_in),
130 240835 : master_g(g_in),
131 240835 : master_action(act_in),
132 240835 : variables(variables_in),
133 271171 : nodes_to_elem(nodes_to_elem_in)
134 : {
135 256003 : if (!nodes_to_elem_in)
136 : {
137 256003 : MeshTools::build_nodes_to_elem_map (system.get_mesh(), nodes_to_elem_ourcopy);
138 256003 : nodes_to_elem = &nodes_to_elem_ourcopy;
139 : }
140 256003 : }
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 496838 : ~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 969401 : 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 256003 : SortAndCopy (GenericProjector & p) : SubFunctor(p) {}
282 :
283 7742 : 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 251219 : struct ProjectVertices : public SubProjector {
308 256003 : ProjectVertices (GenericProjector & p) : SubProjector(p) {}
309 :
310 5958 : 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 8162 : struct ProjectEdges : public SubProjector {
328 256003 : 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 9396 : struct ProjectSides : public SubProjector {
347 256003 : ProjectSides (GenericProjector & p) : SubProjector(p) {}
348 :
349 3684 : 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 871 : struct ProjectInteriors : public SubProjector {
370 256003 : ProjectInteriors (GenericProjector & p) : SubProjector(p) {}
371 :
372 1804 : 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 255373 : VectorSetAction(NumericVector<Val> & target_vec) :
414 255373 : target_vector(target_vec) {}
415 :
416 13992925 : void insert(dof_id_type id,
417 : Val val)
418 : {
419 : // Lock the new vector since it is shared among threads.
420 : {
421 2737484 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
422 13992925 : target_vector.set(id, val);
423 : }
424 13992925 : }
425 :
426 :
427 3065862 : void insert(const std::vector<dof_id_type> & dof_indices,
428 : const DenseVector<Val> & Ue)
429 : {
430 : const numeric_index_type
431 3065862 : first = target_vector.first_local_index(),
432 3065862 : last = target_vector.last_local_index();
433 :
434 331229 : unsigned int size = Ue.size();
435 :
436 331229 : libmesh_assert_equal_to(size, dof_indices.size());
437 :
438 : // Lock the new vector since it is shared among threads.
439 : {
440 662458 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
441 :
442 4358337 : for (unsigned int i = 0; i != size; ++i)
443 1405783 : if ((dof_indices[i] >= first) && (dof_indices[i] < last))
444 1405783 : target_vector.set(dof_indices[i], Ue(i));
445 : }
446 3065862 : }
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 402612 : 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 401484 : FEMFunctionWrapper(const FEMFunctionBase<Output> & f) : _f(f.clone()) {}
471 :
472 53293 : FEMFunctionWrapper(const FEMFunctionWrapper<Output> & fw) :
473 1043379 : _f(fw._f->clone()) {}
474 :
475 1025305 : void init_context (FEMContext & c) { _f->init_context(c); }
476 :
477 81086 : 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 903277 : { return _f->component(c, i, n, time); }
484 :
485 555104 : 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 6901700 : { 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 144337 : 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 119718 : 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 411560 : OldSolutionBase(const libMesh::System & sys_in,
545 : const std::vector<unsigned int> * vars) :
546 376081 : last_elem(nullptr),
547 376081 : sys(sys_in),
548 411560 : 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 1938488 : auto make_lookup = [this](unsigned int v)
553 : {
554 488997 : const unsigned int vcomp = sys.variable_scalar_number(v,0);
555 508982 : if (vcomp >= component_to_var.size())
556 488997 : component_to_var.resize(vcomp+1, static_cast<unsigned int>(-1));
557 508982 : component_to_var[vcomp] = v;
558 : };
559 :
560 411560 : if (vars)
561 900557 : for (auto v : *vars)
562 488997 : make_lookup(v);
563 : else
564 0 : for (auto v : make_range(sys.n_vars()))
565 0 : make_lookup(v);
566 411560 : }
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 302534 : void init_context (FEMContext & c)
581 : {
582 13923 : c.set_algebraic_type(FEMContext::DOFS_ONLY);
583 :
584 : const std::set<unsigned char> & elem_dims =
585 13923 : old_context.elem_dimensions();
586 :
587 : // Loop over variables and dimensions, to prerequest
588 606163 : for (const auto & dim : elem_dims)
589 : {
590 13963 : FEAbstract * fe = nullptr;
591 303629 : if (old_context.active_vars())
592 663099 : for (auto var : *old_context.active_vars())
593 : {
594 359470 : old_context.get_element_fe(var, fe, dim);
595 15740 : 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 302534 : }
605 :
606 722804 : 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 8070991 : bool check_old_context (const FEMContext & c, const Point & p)
643 : {
644 1363038 : LOG_SCOPE ("check_old_context(c,p)", "OldSolutionBase");
645 1363024 : const Elem & elem = c.get_elem();
646 8070991 : if (last_elem != &elem)
647 : {
648 3530520 : if (elem.refinement_flag() == Elem::JUST_REFINED)
649 : {
650 3273943 : old_context.pre_fe_reinit(sys, elem.parent());
651 : }
652 236346 : 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 228288 : const Real master_tol = out_of_elem_tol / elem.hmax() * 2;
662 :
663 396650 : for (auto & child : elem.child_ref_range())
664 396650 : if (child.close_to_point(p, master_tol))
665 : {
666 228288 : old_context.pre_fe_reinit(sys, &child);
667 19530 : break;
668 : }
669 :
670 19530 : libmesh_assert
671 : (old_context.get_elem().close_to_point(p, master_tol));
672 : }
673 : else
674 : {
675 8058 : if (!elem.get_old_dof_object())
676 0 : return false;
677 :
678 8058 : old_context.pre_fe_reinit(sys, &elem);
679 : }
680 :
681 3248082 : last_elem = &elem;
682 : }
683 : else
684 : {
685 398964 : libmesh_assert(old_context.has_elem());
686 :
687 4822909 : const Real master_tol = out_of_elem_tol / elem.hmax() * 2;
688 :
689 4822909 : if (!old_context.get_elem().close_to_point(p, master_tol))
690 : {
691 11638 : libmesh_assert_equal_to
692 : (elem.refinement_flag(), Elem::JUST_COARSENED);
693 :
694 423490 : for (auto & child : elem.child_ref_range())
695 423490 : if (child.close_to_point(p, master_tol))
696 : {
697 143497 : old_context.pre_fe_reinit(sys, &child);
698 11638 : break;
699 : }
700 :
701 11638 : libmesh_assert
702 : (old_context.get_elem().close_to_point(p, master_tol));
703 : }
704 : }
705 :
706 681519 : 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 120418 : 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 57659 : 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 55771 : old_solution(old_sol)
745 : {
746 3776 : this->old_context.set_algebraic_type(FEMContext::OLD);
747 3776 : this->old_context.set_custom_solution(&old_solution);
748 3776 : }
749 :
750 299228 : OldSolutionValue(const OldSolutionValue & in) :
751 : OldSolutionBase<Output, point_output>(in.sys, in.old_context.active_vars()),
752 299864 : old_solution(in.old_solution)
753 : {
754 13781 : this->old_context.set_algebraic_type(FEMContext::OLD);
755 13781 : this->old_context.set_custom_solution(&old_solution);
756 285438 : }
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 7970932 : 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 1345366 : LOG_SCOPE ("eval_at_point()", "OldSolutionValue");
773 :
774 7970932 : if (!skip_context_check)
775 7951552 : if (!this->check_old_context(c, p))
776 0 : return Output(0);
777 :
778 : // Handle offset from non-scalar components in previous variables
779 672683 : libmesh_assert_less(i, this->component_to_var.size());
780 7990857 : unsigned int var = this->component_to_var[i];
781 :
782 314725 : Output n;
783 7970932 : (this->old_context.*point_output)(var, p, n, this->out_of_elem_tol);
784 7970932 : 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 4606680 : 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 972762 : 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 486381 : indices.clear();
848 :
849 5092821 : this->sys.get_dof_map().dof_indices(elem, node_num, indices, var_num);
850 :
851 972762 : std::vector<dof_id_type> old_indices;
852 :
853 5092821 : this->sys.get_dof_map().old_dof_indices(elem, node_num, old_indices, var_num);
854 :
855 486381 : 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 486381 : bool invalid_old_index = false;
860 10053583 : for (const auto & di : old_indices)
861 5446903 : if (di == DofObject::invalid_id)
862 0 : invalid_old_index = true;
863 :
864 4606680 : values.resize(old_indices.size());
865 4606680 : 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 4606680 : old_solution.get(old_indices, values);
878 4606680 : }
879 :
880 :
881 2693470 : 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 592894 : 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 592396 : const Elem & old_elem =
893 2693470 : (elem.refinement_flag() == Elem::JUST_REFINED) ?
894 936 : *elem.parent() : elem;
895 :
896 : // If there are any element-based DOF numbers, get them
897 2693470 : const unsigned int nc = FEInterface::n_dofs_per_elem(fe_type, &elem);
898 :
899 2693470 : 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 2693470 : if (nc != 0)
907 : {
908 35768 : const DofObject & old_dof_object = old_elem.get_old_dof_object_ref();
909 35768 : libmesh_assert_greater(elem.n_systems(), sys_num);
910 :
911 : const std::pair<unsigned int, unsigned int>
912 392151 : vg_and_offset = elem.var_to_vg_and_offset(sys_num,var_num);
913 392151 : const unsigned int vg = vg_and_offset.first;
914 392151 : const unsigned int vig = vg_and_offset.second;
915 :
916 427919 : unsigned int n_comp = old_dof_object.n_comp_group(sys_num,vg);
917 427919 : n_comp = std::min(n_comp, nc);
918 :
919 427919 : std::vector<dof_id_type> old_dof_indices(n_comp);
920 :
921 1626926 : for (unsigned int i=0; i != n_comp; ++i)
922 : {
923 102725 : const dof_id_type d_old =
924 1096282 : old_dof_object.dof_number(sys_num, vg, vig, i, n_comp);
925 102725 : const dof_id_type d_new =
926 1096282 : elem.dof_number(sys_num, vg, vig, i, n_comp);
927 102725 : libmesh_assert_not_equal_to (d_old, DofObject::invalid_id);
928 102725 : libmesh_assert_not_equal_to (d_new, DofObject::invalid_id);
929 :
930 1199007 : old_dof_indices[i] = d_old;
931 1301732 : indices[i] = d_new;
932 : }
933 :
934 427919 : values.resize(n_comp);
935 427919 : old_solution.get(old_dof_indices, values);
936 :
937 427919 : 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 427919 : values.resize(nc, 0);
946 : }
947 : else
948 260679 : values.clear();
949 2693470 : }
950 :
951 : private:
952 : const NumericVector<Number> & old_solution;
953 : };
954 :
955 : template<>
956 : inline void
957 14134 : OldSolutionBase<Number, &FEMContext::point_value>::get_shape_outputs(FEAbstract & fe)
958 : {
959 313247 : fe.request_phi();
960 14134 : }
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 1876263 : 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 326856 : LOG_SCOPE ("Number eval_at_node()", "OldSolutionValue");
1031 :
1032 : // This should only be called on vertices
1033 163428 : 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 163428 : libmesh_assert_less(i, this->component_to_var.size());
1038 1876263 : 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 326813 : const Elem::RefinementState flag = c.get_elem().refinement_flag();
1051 :
1052 163428 : const DofObject * old_dof_object = n.get_old_dof_object();
1053 1399169 : if (old_dof_object &&
1054 1340616 : (!extra_hanging_dofs ||
1055 1140138 : flag == Elem::JUST_COARSENED ||
1056 1115205 : flag == Elem::DO_NOTHING) &&
1057 4115881 : old_dof_object->n_vars(sys.number()) &&
1058 1115205 : old_dof_object->n_comp(sys.number(), var))
1059 : {
1060 : const dof_id_type old_id =
1061 1092062 : old_dof_object->dof_number(sys.number(), var, 0);
1062 1092062 : return old_solution(old_id);
1063 : }
1064 :
1065 784201 : 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 12582 : 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 12582 : 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 2114 : const Elem::RefinementState flag = c.get_elem().refinement_flag();
1160 :
1161 1058 : const DofObject * old_dof_object = n.get_old_dof_object();
1162 13640 : if (old_dof_object &&
1163 13638 : (!extra_hanging_dofs ||
1164 12582 : flag == Elem::JUST_COARSENED ||
1165 464 : flag == Elem::DO_NOTHING) &&
1166 14529 : old_dof_object->n_vars(sys.number()) &&
1167 464 : old_dof_object->n_comp(sys.number(), var))
1168 : {
1169 39 : Gradient g;
1170 928 : for (unsigned int d = 0; d != elem_dim; ++d)
1171 : {
1172 : const dof_id_type old_id =
1173 464 : old_dof_object->dof_number(sys.number(), var, d+1);
1174 464 : g(d) = old_solution(old_id);
1175 : }
1176 464 : return g;
1177 : }
1178 :
1179 12118 : 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 256003 : void GenericProjector<FFunctor, GFunctor, FValue, ProjectionAction>::project
1214 : (const ConstElemRange & range)
1215 : {
1216 15168 : 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 256003 : done_saving_ids = false;
1221 :
1222 271171 : SortAndCopy sort_work(*this);
1223 256003 : Threads::parallel_reduce (range, sort_work);
1224 256003 : 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 15168 : ids_to_push;
1229 :
1230 248419 : ids_to_push.insert(sort_work.new_ids_to_push.begin(),
1231 : sort_work.new_ids_to_push.end());
1232 248419 : ids_to_save.insert(sort_work.new_ids_to_save.begin(),
1233 : sort_work.new_ids_to_save.end());
1234 :
1235 271171 : std::vector<node_projection> vertices(sort_work.vertices.begin(),
1236 : sort_work.vertices.end());
1237 :
1238 498615 : done_saving_ids = sort_work.edges.empty() &&
1239 464648 : sort_work.sides.empty() && sort_work.interiors.empty();
1240 :
1241 : {
1242 15168 : ProjectVertices project_vertices(*this);
1243 504422 : Threads::parallel_reduce (node_range(&vertices), project_vertices);
1244 7584 : libmesh_merge_move(ids_to_push, project_vertices.new_ids_to_push);
1245 7584 : libmesh_merge_move(ids_to_save, project_vertices.new_ids_to_save);
1246 : }
1247 :
1248 256003 : done_saving_ids = sort_work.sides.empty() && sort_work.interiors.empty();
1249 :
1250 256003 : this->send_and_insert_dof_values(ids_to_push, action);
1251 :
1252 : {
1253 271171 : std::vector<node_projection> edges(sort_work.edges.begin(), sort_work.edges.end());
1254 15168 : ProjectEdges project_edges(*this);
1255 504422 : Threads::parallel_reduce (node_range(&edges), project_edges);
1256 7584 : libmesh_merge_move(ids_to_push, project_edges.new_ids_to_push);
1257 7584 : libmesh_merge_move(ids_to_save, project_edges.new_ids_to_save);
1258 240835 : }
1259 :
1260 256003 : done_saving_ids = sort_work.interiors.empty();
1261 :
1262 256003 : this->send_and_insert_dof_values(ids_to_push, action);
1263 :
1264 : {
1265 271171 : std::vector<node_projection> sides(sort_work.sides.begin(), sort_work.sides.end());
1266 15168 : ProjectSides project_sides(*this);
1267 504422 : Threads::parallel_reduce (node_range(&sides), project_sides);
1268 7584 : libmesh_merge_move(ids_to_push, project_sides.new_ids_to_push);
1269 7584 : libmesh_merge_move(ids_to_save, project_sides.new_ids_to_save);
1270 240835 : }
1271 :
1272 256003 : done_saving_ids = true;
1273 256003 : 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 7584 : ProjectInteriors project_interiors(*this);
1278 504422 : Threads::parallel_reduce (interior_range(&sort_work.interiors),
1279 : project_interiors);
1280 496838 : }
1281 :
1282 :
1283 : template <typename FFunctor, typename GFunctor,
1284 : typename FValue, typename ProjectionAction>
1285 1306911 : GenericProjector<FFunctor, GFunctor, FValue, ProjectionAction>::SubFunctor::SubFunctor
1286 : (GenericProjector & p) :
1287 1212570 : projector(p),
1288 1306911 : action(p.master_action),
1289 1306911 : f(p.master_f),
1290 1306911 : context(p.system, &p.variables, /* allocate_local_matrices= */ false),
1291 1306911 : conts(p.system.n_vars()),
1292 1401252 : 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 2688108 : for (const auto & var : this->projector.variables)
1297 : {
1298 : // FIXME: Need to generalize this to vector-valued elements. [PB]
1299 49783 : FEAbstract * fe = nullptr;
1300 49783 : FEAbstract * side_fe = nullptr;
1301 49783 : FEAbstract * edge_fe = nullptr;
1302 :
1303 : const std::set<unsigned char> & elem_dims =
1304 49783 : context.elem_dimensions();
1305 :
1306 2765107 : for (const auto & dim : elem_dims)
1307 : {
1308 1383910 : context.get_element_fe( var, fe, dim );
1309 1383910 : if (fe->get_fe_type().family == SCALAR)
1310 0 : continue;
1311 1383910 : if (dim > 1)
1312 44636 : context.get_side_fe( var, side_fe, dim );
1313 1383910 : if (dim > 2)
1314 17742 : context.get_edge_fe( var, edge_fe );
1315 :
1316 118354 : fe->get_JxW();
1317 118354 : fe->get_xyz();
1318 118354 : fe->get_JxW();
1319 :
1320 1383910 : fe->request_phi();
1321 1383910 : if (dim > 1)
1322 : {
1323 105561 : side_fe->get_JxW();
1324 105561 : side_fe->get_xyz();
1325 1217256 : side_fe->request_phi();
1326 : }
1327 1383910 : if (dim > 2)
1328 : {
1329 42839 : edge_fe->get_JxW();
1330 42839 : edge_fe->get_xyz();
1331 530647 : edge_fe->request_phi();
1332 : }
1333 :
1334 1383910 : const FEContinuity cont = fe->get_continuity();
1335 1383910 : this->conts[var] = cont;
1336 1383910 : 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 1383910 : 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 316825 : f.init_context(context);
1351 1306911 : }
1352 :
1353 :
1354 : template <typename FFunctor, typename GFunctor,
1355 : typename FValue, typename ProjectionAction>
1356 1043166 : GenericProjector<FFunctor, GFunctor, FValue, ProjectionAction>::SubProjector::SubProjector
1357 : (GenericProjector & p) :
1358 1043166 : SubFunctor(p)
1359 : {
1360 :
1361 : // Our C1 elements need gradient information
1362 2106391 : for (const auto & var : this->projector.variables)
1363 1141085 : 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 5419099 : 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 5419099 : if (id == DofObject::invalid_id)
1397 0 : return;
1398 :
1399 456919 : auto iter = new_ids_to_push.find(id);
1400 5419099 : if (iter == new_ids_to_push.end())
1401 5413155 : action.insert(id, val);
1402 : else
1403 : {
1404 197 : libmesh_assert(pid != DofObject::invalid_processor_id);
1405 1095 : iter->second = std::make_pair(val, pid);
1406 : }
1407 5419099 : if (!this->projector.done_saving_ids)
1408 : {
1409 216841 : libmesh_assert(!new_ids_to_save.count(id));
1410 2581719 : 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 5231721 : 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 531609 : libmesh_assert_equal_to(ids.size(), vals.size());
1441 15321332 : for (auto i : index_range(ids))
1442 : {
1443 10089611 : const dof_id_type id = ids[i];
1444 :
1445 : // We may see invalid ids when expanding a subdomain with a
1446 : // restricted variable
1447 10089611 : if (id == DofObject::invalid_id)
1448 0 : continue;
1449 :
1450 1849726 : const InsertInput & val = vals[i];
1451 :
1452 924977 : auto iter = new_ids_to_push.find(id);
1453 10089611 : if (iter == new_ids_to_push.end())
1454 10084328 : action.insert(id, val);
1455 : else
1456 : {
1457 504 : libmesh_assert(pid != DofObject::invalid_processor_id);
1458 1008 : iter->second = std::make_pair(val, pid);
1459 : }
1460 10089611 : if (!this->projector.done_saving_ids)
1461 : {
1462 618069 : libmesh_assert(!new_ids_to_save.count(id));
1463 6280986 : new_ids_to_save[id] = val;
1464 : }
1465 : }
1466 5231721 : }
1467 :
1468 : template <typename FFunctor, typename GFunctor,
1469 : typename FValue, typename ProjectionAction>
1470 26896 : void GenericProjector<FFunctor, GFunctor, FValue, ProjectionAction>::SubFunctor::join
1471 : (const GenericProjector<FFunctor, GFunctor, FValue, ProjectionAction>::SubFunctor & other)
1472 : {
1473 17641 : new_ids_to_push.insert(other.new_ids_to_push.begin(), other.new_ids_to_push.end());
1474 17641 : new_ids_to_save.insert(other.new_ids_to_save.begin(), other.new_ids_to_save.end());
1475 26896 : }
1476 :
1477 :
1478 : template <typename FFunctor, typename GFunctor,
1479 : typename FValue, typename ProjectionAction>
1480 263745 : void GenericProjector<FFunctor, GFunctor, FValue, ProjectionAction>::SortAndCopy::operator()
1481 : (const ConstElemRange & range)
1482 : {
1483 20576 : 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 20576 : std::unordered_map<const Node *, std::pair<var_set, var_set>> copied_nodes;
1510 :
1511 65391 : 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 20576 : std::vector<unsigned short> extra_hanging_dofs;
1517 10288 : bool all_extra_hanging_dofs = true;
1518 542715 : for (auto v_num : this->projector.variables)
1519 : {
1520 289904 : if (extra_hanging_dofs.size() <= v_num)
1521 278970 : extra_hanging_dofs.resize(v_num+1, false);
1522 300838 : extra_hanging_dofs[v_num] =
1523 278970 : FEInterface::extra_hanging_dofs(system.variable(v_num).type());
1524 :
1525 278970 : if (!extra_hanging_dofs[v_num])
1526 8345 : all_extra_hanging_dofs = false;
1527 : }
1528 :
1529 5792492 : 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 554745 : bool copy_this_elem = false;
1534 :
1535 : #ifdef LIBMESH_ENABLE_AMR
1536 : // If we're projecting from an old grid
1537 554745 : 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 4446925 : const DofObject * old_dof_object = elem->get_old_dof_object();
1544 916804 : const Elem::RefinementState h_flag = elem->refinement_flag();
1545 916804 : const Elem::RefinementState p_flag = elem->p_refinement_flag();
1546 4905316 : if (!old_dof_object &&
1547 3988525 : 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 4446769 : if (h_flag != Elem::JUST_REFINED &&
1553 295380 : h_flag != Elem::JUST_COARSENED &&
1554 2657609 : p_flag != Elem::JUST_REFINED &&
1555 : p_flag != Elem::JUST_COARSENED)
1556 294753 : copy_this_elem = true;
1557 : else
1558 : {
1559 163647 : bool reinitted = false;
1560 :
1561 327302 : 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 163647 : const bool copy_possible =
1566 1958938 : p_level == 0 &&
1567 1958930 : h_flag != Elem::JUST_COARSENED &&
1568 : p_flag != Elem::JUST_COARSENED;
1569 :
1570 1961148 : std::vector<typename FFunctor::ValuePushType> Ue(1);
1571 1959660 : std::vector<dof_id_type> elem_dof_ids(1);
1572 :
1573 3644682 : for (auto v_num : this->projector.variables)
1574 : {
1575 1848669 : const Variable & var = system.variable(v_num);
1576 1848669 : if (!var.active_on_subdomain(elem->subdomain_id()))
1577 7368 : continue;
1578 1841301 : const FEType fe_type = var.type();
1579 :
1580 41682 : if (fe_type.family == MONOMIAL &&
1581 1849983 : 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 5528591 : const int dim = elem->dim();
1601 :
1602 5528591 : const unsigned int n_vertices = elem->n_vertices();
1603 5528591 : const unsigned int n_edges = elem->n_edges();
1604 5528591 : const unsigned int n_nodes = elem->n_nodes();
1605 :
1606 : // In 1-D we already handle our sides as vertices
1607 5528591 : 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 1109464 : 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 5528591 : const bool has_poly_midnode = (elem->type() == C0POLYHEDRON) && (n_nodes == n_vertices + 1);
1615 5528591 : const bool has_edge_nodes = (n_nodes > n_vertices + has_poly_midnode && dim > 2);
1616 :
1617 : // If we have even more nodes, the next is a side node.
1618 554732 : const bool has_side_nodes =
1619 5528591 : (n_nodes > n_vertices + ((dim > 2) * n_edges));
1620 :
1621 : // We may be out of nodes at this point or we have interior
1622 : // nodes which may have DoFs to project too
1623 554732 : const bool has_interior_nodes =
1624 5528591 : (n_nodes > n_vertices + ((dim > 2) * n_edges) + n_sides);
1625 :
1626 11159056 : for (auto v_num : this->projector.variables)
1627 : {
1628 5630465 : const Variable & var = this->projector.system.variable(v_num);
1629 5630465 : if (!var.active_on_subdomain(elem->subdomain_id()))
1630 548031 : continue;
1631 5615219 : const FEType fe_type = var.type();
1632 5615219 : const bool add_p_level = fe_type.p_refinement;
1633 :
1634 : // If we're trying to do projections on an isogeometric
1635 : // analysis mesh, only the finite element nodes constrained
1636 : // by those spline nodes truly have delta function degrees
1637 : // of freedom. We'll have to back out the spline degrees of
1638 : // freedom indirectly later.
1639 5627933 : if (fe_type.family == RATIONAL_BERNSTEIN &&
1640 12714 : elem->type() == NODEELEM)
1641 5435 : continue;
1642 :
1643 5609311 : if (FEInterface::n_dofs_at_node(fe_type, elem, 0, add_p_level))
1644 4843373 : vertex_vars.insert(vertex_vars.end(), v_num);
1645 :
1646 : // The first non-vertex node is always an edge node if those
1647 : // exist. All edge nodes have the same number of DoFs
1648 5609311 : if (has_edge_nodes)
1649 191472 : if (FEInterface::n_dofs_at_node(fe_type, elem, n_vertices, add_p_level))
1650 58013 : edge_vars.insert(edge_vars.end(), v_num);
1651 :
1652 5609311 : if (has_side_nodes)
1653 : {
1654 1587426 : if (dim != 3)
1655 : {
1656 1435746 : if (FEInterface::n_dofs_at_node(fe_type, elem, n_vertices, add_p_level))
1657 1218665 : side_vars.insert(side_vars.end(), v_num);
1658 : }
1659 : else
1660 : // In 3D, not all face nodes always have the same number of
1661 : // DoFs! We'll loop over all sides to be safe.
1662 2276558 : for (unsigned int n = 0; n != n_nodes; ++n)
1663 2235408 : if (elem->is_face(n))
1664 302116 : if (FEInterface::n_dofs_at_node(fe_type, elem, n, add_p_level))
1665 : {
1666 100722 : side_vars.insert(side_vars.end(), v_num);
1667 9808 : break;
1668 : }
1669 : }
1670 :
1671 5684418 : if (FEInterface::n_dofs_per_elem(fe_type, elem, add_p_level) ||
1672 855303 : (has_interior_nodes &&
1673 855303 : FEInterface::n_dofs_at_node(fe_type, elem, n_nodes-1, add_p_level)))
1674 : {
1675 : #ifdef LIBMESH_ENABLE_AMR
1676 : // We may have already just copied constant monomials,
1677 : // or we may be about to copy the whole element
1678 1783712 : if ((f.is_grid_projection() &&
1679 1344688 : fe_type.family == MONOMIAL &&
1680 12565 : fe_type.order == CONSTANT &&
1681 9213 : elem->p_level() == 0 &&
1682 9201 : elem->refinement_flag() != Elem::JUST_COARSENED &&
1683 1528 : elem->p_refinement_flag() != Elem::JUST_COARSENED)
1684 1457866 : || copy_this_elem
1685 : )
1686 480073 : continue;
1687 : #endif // LIBMESH_ENABLE_AMR
1688 :
1689 : // We need to project any other variables
1690 931916 : if (interiors.empty() || interiors.back() != elem)
1691 910596 : interiors.push_back(elem);
1692 : }
1693 : }
1694 :
1695 : // We'll use a greedy algorithm in most cases: if another
1696 : // element has already claimed some of our DoFs, we'll let it do
1697 : // the work.
1698 33693692 : auto erase_covered_vars = []
1699 6463860 : (const Node * node, var_set & remaining, const ves_multimap & covered)
1700 : {
1701 3232042 : auto covered_range = covered.equal_range(node);
1702 42468869 : for (const auto & v_ent : as_range(covered_range))
1703 17684614 : for (const unsigned int var_covered :
1704 794438 : std::get<2>(v_ent.second))
1705 807026 : remaining.erase(var_covered);
1706 : };
1707 :
1708 29216732 : auto erase_nonhanging_vars = [&extra_hanging_dofs]
1709 5258220 : (const Node * node, var_set & remaining, const ves_multimap & covered)
1710 : {
1711 2629170 : auto covered_range = covered.equal_range(node);
1712 23691836 : for (const auto & v_ent : as_range(covered_range))
1713 8406 : for (const unsigned int var_covered :
1714 257 : std::get<2>(v_ent.second))
1715 4711 : if (!extra_hanging_dofs[var_covered])
1716 273 : remaining.erase(var_covered);
1717 : };
1718 :
1719 24238060 : auto erase_copied_vars = [&copied_nodes, &extra_hanging_dofs]
1720 3498880 : (const Node * node, bool is_vertex, var_set & remaining)
1721 : {
1722 1749490 : auto copying_range = copied_nodes.equal_range(node);
1723 24725649 : for (const auto & v_ent : as_range(copying_range))
1724 : {
1725 14887076 : for (const unsigned int var_covered :
1726 : v_ent.second.first)
1727 7098571 : if (is_vertex || !extra_hanging_dofs[var_covered])
1728 783391 : remaining.erase(var_covered);
1729 :
1730 8717289 : for (const unsigned int var_covered :
1731 : v_ent.second.second)
1732 928784 : if (!is_vertex || !extra_hanging_dofs[var_covered])
1733 70262 : remaining.erase(var_covered);
1734 : }
1735 : };
1736 :
1737 : // Treat polyhedron midnode as a vertex
1738 : // NOTE: if we start having edge or side nodes on polyhedra, we need to use that +1 offset
1739 : // in the edge and side nodes code as well!
1740 : #ifndef NDEBUG
1741 1116338 : for (auto v_num : this->projector.variables)
1742 : {
1743 561606 : const auto & family = this->projector.system.variable(v_num).type().family;
1744 : // Add to the list once known to be correctly functioning with the midnode
1745 561606 : libmesh_assert(!has_poly_midnode || (family == LAGRANGE || family == MONOMIAL || family == XYZ));
1746 : }
1747 : #endif
1748 26405933 : for (const auto v : make_range(n_vertices + has_poly_midnode))
1749 : {
1750 22993536 : const Node * node = elem->node_ptr(v);
1751 :
1752 2116274 : auto remaining_vars = vertex_vars;
1753 :
1754 20877342 : erase_covered_vars(node, remaining_vars, vertices);
1755 :
1756 20877342 : if (remaining_vars.empty())
1757 7503562 : continue;
1758 :
1759 12629621 : if (!all_extra_hanging_dofs)
1760 : {
1761 10693835 : erase_nonhanging_vars(node, remaining_vars, edges);
1762 10693835 : if (remaining_vars.empty())
1763 0 : continue;
1764 :
1765 10693835 : erase_nonhanging_vars(node, remaining_vars, sides);
1766 10693835 : if (remaining_vars.empty())
1767 2169 : continue;
1768 : }
1769 :
1770 12627190 : erase_copied_vars(node, true, remaining_vars);
1771 12627190 : if (remaining_vars.empty())
1772 6223577 : continue;
1773 :
1774 4785037 : if (copy_this_elem)
1775 : {
1776 657832 : std::vector<dof_id_type> node_dof_ids;
1777 657832 : std::vector<typename FFunctor::ValuePushType> values;
1778 :
1779 5653062 : for (auto var : remaining_vars)
1780 : {
1781 2879056 : f.eval_old_dofs(*elem, v, var, node_dof_ids, values);
1782 2879056 : insert_ids(node_dof_ids, values, node->processor_id());
1783 : }
1784 2774006 : copied_nodes[node].first.insert(remaining_vars.begin(),
1785 : remaining_vars.end());
1786 2774006 : this->find_dofs_to_send(*node, *elem, v, remaining_vars);
1787 0 : }
1788 : else
1789 : vertices.emplace
1790 3103328 : (node, std::make_tuple(elem, v, std::move(remaining_vars)));
1791 : }
1792 :
1793 5528591 : if (has_edge_nodes)
1794 : {
1795 1491282 : for (unsigned int e=0; e != n_edges; ++e)
1796 : {
1797 1423830 : const Node * node = elem->node_ptr(n_vertices+e);
1798 :
1799 115212 : auto remaining_vars = edge_vars;
1800 :
1801 1308618 : erase_covered_vars(node, remaining_vars, edges);
1802 1308618 : if (remaining_vars.empty())
1803 921995 : continue;
1804 :
1805 300953 : erase_covered_vars(node, remaining_vars, sides);
1806 300953 : if (remaining_vars.empty())
1807 0 : continue;
1808 :
1809 300953 : if (!all_extra_hanging_dofs)
1810 : {
1811 205736 : erase_nonhanging_vars(node, remaining_vars, vertices);
1812 205736 : if (remaining_vars.empty())
1813 383 : continue;
1814 : }
1815 :
1816 300570 : erase_copied_vars(node, false, remaining_vars);
1817 300570 : if (remaining_vars.empty())
1818 16711 : continue;
1819 :
1820 127774 : if (copy_this_elem)
1821 : {
1822 7764 : std::vector<dof_id_type> edge_dof_ids;
1823 7764 : std::vector<typename FFunctor::ValuePushType> values;
1824 :
1825 52764 : for (auto var : remaining_vars)
1826 : {
1827 26382 : f.eval_old_dofs(*elem, n_vertices+e, var, edge_dof_ids, values);
1828 26382 : insert_ids(edge_dof_ids, values, node->processor_id());
1829 : }
1830 :
1831 26382 : copied_nodes[node].second.insert(remaining_vars.begin(),
1832 : remaining_vars.end());
1833 26382 : this->find_dofs_to_send(*node, *elem, n_vertices+e, remaining_vars);
1834 0 : }
1835 : else
1836 : edges.emplace
1837 273901 : (node, std::make_tuple(elem, e, std::move(remaining_vars)));
1838 : }
1839 : }
1840 :
1841 5528591 : if (has_side_nodes)
1842 : {
1843 7320750 : for (unsigned int side=0; side != n_sides; ++side)
1844 : {
1845 5795808 : const Node * node = nullptr;
1846 5795808 : unsigned short node_num = n_vertices+(dim>2)*n_edges+side;
1847 5795808 : if (dim != 3)
1848 5607419 : node = elem->node_ptr(node_num);
1849 : else
1850 : {
1851 : // In 3D only some sides may have nodes
1852 9840228 : for (unsigned int n = 0; n != n_nodes; ++n)
1853 : {
1854 9822816 : if (!elem->is_face(n))
1855 7371150 : continue;
1856 :
1857 1734426 : if (elem->is_node_on_side(n, side))
1858 : {
1859 152038 : node_num = n;
1860 617544 : node = elem->node_ptr(node_num);
1861 617544 : break;
1862 : }
1863 : }
1864 : }
1865 :
1866 5795808 : if (!node)
1867 2424957 : continue;
1868 :
1869 500057 : auto remaining_vars = side_vars;
1870 :
1871 5778396 : erase_covered_vars(node, remaining_vars, edges);
1872 5778396 : if (remaining_vars.empty())
1873 320913 : continue;
1874 :
1875 5428383 : erase_covered_vars(node, remaining_vars, sides);
1876 5428383 : if (remaining_vars.empty())
1877 1295374 : continue;
1878 :
1879 4010199 : if (!all_extra_hanging_dofs)
1880 : {
1881 2094453 : erase_nonhanging_vars(node, remaining_vars, vertices);
1882 2094453 : if (remaining_vars.empty())
1883 815 : continue;
1884 : }
1885 :
1886 4009384 : erase_copied_vars(node, false, remaining_vars);
1887 4009384 : if (remaining_vars.empty())
1888 576680 : continue;
1889 :
1890 3053817 : if (copy_this_elem)
1891 : {
1892 228174 : std::vector<dof_id_type> side_dof_ids;
1893 228174 : std::vector<typename FFunctor::ValuePushType> values;
1894 :
1895 2624612 : for (auto var : remaining_vars)
1896 : {
1897 1323180 : f.eval_old_dofs(*elem, node_num, var, side_dof_ids, values);
1898 1323180 : insert_ids(side_dof_ids, values, node->processor_id());
1899 : }
1900 :
1901 1301432 : copied_nodes[node].second.insert(remaining_vars.begin(),
1902 : remaining_vars.end());
1903 1301432 : this->find_dofs_to_send(*node, *elem, node_num, remaining_vars);
1904 0 : }
1905 : else
1906 : sides.emplace
1907 2241710 : (node, std::make_tuple(elem, side, std::move(remaining_vars)));
1908 : }
1909 : }
1910 :
1911 : // Elements with elemental dofs might need those copied too.
1912 4543101 : if (copy_this_elem)
1913 : {
1914 589506 : std::vector<typename FFunctor::ValuePushType> U;
1915 589506 : std::vector<dof_id_type> dof_ids;
1916 :
1917 5343374 : for (auto v_num : this->projector.variables)
1918 : {
1919 2692618 : const Variable & var = system.variable(v_num);
1920 2692618 : if (!var.active_on_subdomain(elem->subdomain_id()))
1921 4818 : continue;
1922 2687800 : FEType fe_type = var.type();
1923 :
1924 2687800 : f.eval_old_dofs(*elem, fe_type, sys_num, v_num,
1925 : dof_ids, U);
1926 4783672 : action.insert(dof_ids, U);
1927 :
1928 2687800 : if (has_interior_nodes)
1929 : {
1930 378062 : f.eval_old_dofs(*elem, n_nodes-1, v_num, dof_ids, U);
1931 720874 : action.insert(dof_ids, U);
1932 : }
1933 : }
1934 0 : }
1935 : }
1936 263745 : }
1937 :
1938 :
1939 : template <typename FFunctor, typename GFunctor,
1940 : typename FValue, typename ProjectionAction>
1941 7742 : void GenericProjector<FFunctor, GFunctor, FValue, ProjectionAction>::SortAndCopy::join
1942 : (const GenericProjector<FFunctor, GFunctor, FValue, ProjectionAction>::SortAndCopy & other)
1943 : {
1944 38340 : auto merge_multimaps = [](ves_multimap & self, const ves_multimap & other_mm)
1945 : {
1946 657710 : for (const auto & pair : other_mm)
1947 : {
1948 634484 : const Node * node = pair.first;
1949 433166 : auto remaining_vars = std::get<2>(pair.second);
1950 :
1951 216583 : auto my_range = self.equal_range(node);
1952 705579 : for (const auto & v_ent : as_range(my_range))
1953 143447 : for (const unsigned int var_covered :
1954 23937 : std::get<2>(v_ent.second))
1955 24356 : remaining_vars.erase(var_covered);
1956 :
1957 634484 : if (!remaining_vars.empty())
1958 : self.emplace
1959 756819 : (node, std::make_tuple(std::get<0>(pair.second),
1960 192844 : std::get<1>(pair.second),
1961 192844 : std::move(remaining_vars)));
1962 : }
1963 : };
1964 :
1965 7742 : merge_multimaps(vertices, other.vertices);
1966 7742 : merge_multimaps(edges, other.edges);
1967 7742 : merge_multimaps(sides, other.sides);
1968 :
1969 7742 : std::sort(interiors.begin(), interiors.end());
1970 10446 : std::vector<const Elem *> other_interiors = other.interiors;
1971 7742 : std::sort(other_interiors.begin(), other_interiors.end());
1972 5408 : std::vector<const Elem *> merged_interiors;
1973 7742 : std::set_union(interiors.begin(), interiors.end(),
1974 : other_interiors.begin(), other_interiors.end(),
1975 : std::back_inserter(merged_interiors));
1976 2704 : interiors.swap(merged_interiors);
1977 :
1978 7742 : SubFunctor::join(other);
1979 7742 : }
1980 :
1981 : namespace
1982 : {
1983 : template <typename Output, typename Input>
1984 : typename std::enable_if<ScalarTraits<Input>::value, Output>::type
1985 0 : raw_value(const Input & input, unsigned int /*index*/)
1986 : {
1987 0 : return input;
1988 : }
1989 :
1990 : template <typename Output, template <typename> class WrapperClass, typename T>
1991 : typename std::enable_if<ScalarTraits<T>::value &&
1992 : TensorTools::MathWrapperTraits<WrapperClass<T>>::value,
1993 : Output>::type
1994 16143 : raw_value(const WrapperClass<T> & input, unsigned int index)
1995 : {
1996 233178 : return input.slice(index);
1997 : }
1998 :
1999 : template <typename T>
2000 : typename T::value_type
2001 : grad_component(const T &, unsigned int);
2002 :
2003 : template <typename T>
2004 : typename VectorValue<T>::value_type
2005 2380 : grad_component(const VectorValue<T> & grad, unsigned int component)
2006 : {
2007 4709 : return grad(component);
2008 : }
2009 :
2010 : template <typename T>
2011 : typename TensorValue<T>::value_type
2012 0 : grad_component(const TensorValue<T> & grad, unsigned int component)
2013 : {
2014 0 : libmesh_error_msg("This function should only be called for Hermites. "
2015 : "I don't know how you got here");
2016 : return grad(component, component);
2017 : }
2018 :
2019 :
2020 : }
2021 :
2022 : template <typename FFunctor, typename GFunctor,
2023 : typename FValue, typename ProjectionAction>
2024 265122 : void GenericProjector<FFunctor, GFunctor, FValue, ProjectionAction>::ProjectVertices::operator()
2025 : (const node_range & range)
2026 : {
2027 21484 : LOG_SCOPE ("project_vertices","GenericProjector");
2028 :
2029 265122 : const unsigned int sys_num = system.number();
2030 :
2031 : // Variables with extra hanging dofs can't safely use eval_at_node
2032 : // in as many places as variables without can.
2033 21484 : std::vector<unsigned short> extra_hanging_dofs;
2034 545274 : for (auto v_num : this->projector.variables)
2035 : {
2036 291478 : if (extra_hanging_dofs.size() <= v_num)
2037 280152 : extra_hanging_dofs.resize(v_num+1, false);
2038 280152 : extra_hanging_dofs[v_num] =
2039 280152 : FEInterface::extra_hanging_dofs(system.variable(v_num).type());
2040 : }
2041 :
2042 3061038 : for (const auto & v_pair : range)
2043 : {
2044 2795916 : const Node & vertex = *v_pair.first;
2045 2795916 : const Elem & elem = *std::get<0>(v_pair.second);
2046 2795916 : const unsigned int n = std::get<1>(v_pair.second);
2047 241749 : const var_set & vertex_vars = std::get<2>(v_pair.second);
2048 :
2049 2795916 : context.pre_fe_reinit(system, &elem);
2050 :
2051 2795916 : this->find_dofs_to_send(vertex, elem, n, vertex_vars);
2052 :
2053 : // Look at all the variables we're supposed to interpolate from
2054 : // this element on this vertex
2055 5651175 : for (const auto & var : vertex_vars)
2056 : {
2057 2855259 : const Variable & variable = system.variable(var);
2058 246527 : const FEType & base_fe_type = variable.type();
2059 739538 : const unsigned int var_component =
2060 2855259 : system.variable_scalar_number(var, 0);
2061 :
2062 2855259 : if (base_fe_type.family == SCALAR)
2063 0 : continue;
2064 :
2065 2855259 : const FEContinuity cont = this->conts[var];
2066 2855259 : const FEFieldType field_type = this->field_types[var];
2067 :
2068 2855259 : if (cont == DISCONTINUOUS)
2069 : {
2070 0 : libmesh_assert_equal_to(vertex.n_comp(sys_num, var), 0);
2071 : }
2072 3101743 : else if (cont == C_ZERO ||
2073 2608732 : cont == SIDE_DISCONTINUOUS)
2074 : {
2075 2768372 : if (cont == SIDE_DISCONTINUOUS &&
2076 1518 : elem.dim() != 1)
2077 : {
2078 0 : libmesh_assert_equal_to(vertex.n_comp(sys_num, var), 0);
2079 0 : continue;
2080 : }
2081 :
2082 3554833 : const FValue val = f.eval_at_node
2083 1860979 : (context, var_component, /*dim=*/ 0, // Don't care w/C0
2084 2766854 : vertex, extra_hanging_dofs[var], system.time);
2085 :
2086 2766854 : if (field_type == TYPE_VECTOR)
2087 : {
2088 1302 : libmesh_assert_equal_to(vertex.n_comp(sys_num, var), elem.dim());
2089 :
2090 : // We will have a number of nodal value DoFs equal to the elem dim
2091 83664 : for (auto i : make_range(elem.dim()))
2092 : {
2093 61920 : const dof_id_type id = vertex.dof_number(sys_num, var, i);
2094 :
2095 : // Need this conversion so that this method
2096 : // will compile for TYPE_SCALAR instantiations
2097 59112 : const auto insert_val =
2098 6498 : raw_value<typename ProjectionAction::InsertInput>(val, i);
2099 :
2100 61920 : insert_id(id, insert_val, vertex.processor_id());
2101 : }
2102 : }
2103 : else
2104 : {
2105 : // C_ZERO elements have a single nodal value DoF at
2106 : // vertices. We can't assert n_comp==1 here,
2107 : // because if this is a hanging node then it may have
2108 : // more face/edge DoFs, but we don't need to deal with
2109 : // those here.
2110 :
2111 2745110 : const dof_id_type id = vertex.dof_number(sys_num, var, 0);
2112 2745110 : insert_id(id, val, vertex.processor_id());
2113 239040 : }
2114 : }
2115 88405 : else if (cont == C_ONE)
2116 : {
2117 7487 : libmesh_assert(vertex.n_comp(sys_num, var));
2118 88405 : const dof_id_type first_id = vertex.dof_number(sys_num, var, 0);
2119 :
2120 : // C_ONE elements have a single nodal value and dim
2121 : // gradient values at vertices, as well as cross
2122 : // gradients for HERMITE. We need to have an element in
2123 : // hand to figure out dim and to have in case this
2124 : // vertex is a new vertex.
2125 88405 : const int dim = elem.dim();
2126 : #ifndef NDEBUG
2127 : // For now all C1 elements at a vertex had better have
2128 : // the same dimension. If anyone hits these asserts let
2129 : // me know; we could probably support a mixed-dimension
2130 : // mesh IFF the 2D elements were all parallel to xy and
2131 : // the 1D elements all parallel to x.
2132 33290 : for (const auto e_id : (*this->projector.nodes_to_elem)[vertex.id()])
2133 : {
2134 25803 : const Elem & e = system.get_mesh().elem_ref(e_id);
2135 25803 : libmesh_assert_equal_to(dim, e.dim());
2136 : }
2137 : #endif
2138 : #ifdef LIBMESH_ENABLE_AMR
2139 7487 : bool is_old_vertex = true;
2140 88405 : if (elem.refinement_flag() == Elem::JUST_REFINED)
2141 : {
2142 3056 : const int i_am_child =
2143 42272 : elem.parent()->which_child_am_i(&elem);
2144 : is_old_vertex =
2145 42272 : elem.parent()->is_vertex_on_parent(i_am_child, n);
2146 : }
2147 : #else
2148 : const bool is_old_vertex = false;
2149 : #endif
2150 :
2151 : // The hermite element vertex shape functions are weird
2152 88405 : if (base_fe_type.family == HERMITE)
2153 : {
2154 98797 : const FValue val =
2155 29385 : f.eval_at_node(context,
2156 : var_component,
2157 : dim,
2158 : vertex,
2159 : extra_hanging_dofs[var],
2160 74698 : system.time);
2161 74698 : insert_id(first_id, val, vertex.processor_id());
2162 :
2163 13134 : typename GFunctor::FunctorValue grad =
2164 69103 : is_old_vertex ?
2165 13836 : g->eval_at_node(context,
2166 : var_component,
2167 : dim,
2168 : vertex,
2169 : extra_hanging_dofs[var],
2170 58018 : system.time) :
2171 15369 : g->eval_at_point(context,
2172 : var_component,
2173 : vertex,
2174 16680 : system.time,
2175 : false);
2176 : // x derivative. Use slice because grad may be a tensor type
2177 74698 : insert_id(first_id+1, grad.slice(0),
2178 : vertex.processor_id());
2179 : #if LIBMESH_DIM > 1
2180 70414 : if (dim > 1 && is_old_vertex && f.is_grid_projection())
2181 : {
2182 5632 : for (int i = 1; i < dim; ++i)
2183 3174 : insert_id(first_id+i+1, grad.slice(i),
2184 : vertex.processor_id());
2185 :
2186 : // We can directly copy everything else too
2187 716 : std::vector<FValue> derivs;
2188 1074 : f.eval_mixed_derivatives
2189 2458 : (context, var_component, dim, vertex, derivs);
2190 5632 : for (auto i : index_range(derivs))
2191 3174 : insert_id(first_id+dim+1+i, derivs[i],
2192 : vertex.processor_id());
2193 0 : }
2194 70442 : else if (dim > 1)
2195 : {
2196 : // We'll finite difference mixed derivatives.
2197 : // This delta_x used to be TOLERANCE*hmin, but
2198 : // the factor of 10 improved the accuracy in
2199 : // some unit test projections
2200 11418 : Real delta_x = TOLERANCE * 10 * elem.hmin();
2201 :
2202 11418 : Point nxminus = elem.point(n),
2203 11418 : nxplus = elem.point(n);
2204 11418 : nxminus(0) -= delta_x;
2205 11418 : nxplus(0) += delta_x;
2206 1660 : typename GFunctor::FunctorValue gxminus =
2207 9076 : g->eval_at_point(context,
2208 : var_component,
2209 : nxminus,
2210 11418 : system.time,
2211 : true);
2212 1660 : typename GFunctor::FunctorValue gxplus =
2213 9076 : g->eval_at_point(context,
2214 : var_component,
2215 : nxplus,
2216 11418 : system.time,
2217 : true);
2218 : // y derivative
2219 11418 : insert_id(first_id+2, grad.slice(1),
2220 : vertex.processor_id());
2221 : // xy derivative
2222 12320 : insert_id(first_id+3,
2223 11418 : (grad_component(gxplus, 1) - grad_component(gxminus, 1)) / 2. / delta_x,
2224 : vertex.processor_id());
2225 :
2226 : #if LIBMESH_DIM > 2
2227 11418 : if (dim > 2)
2228 : {
2229 : // z derivative
2230 864 : insert_id(first_id+4, grad.slice(2),
2231 : vertex.processor_id());
2232 : // xz derivative
2233 936 : insert_id(first_id+5,
2234 864 : (grad_component(gxplus, 2) - grad_component(gxminus, 2)) / 2. / delta_x,
2235 : vertex.processor_id());
2236 :
2237 : // We need new points for yz
2238 864 : Point nyminus = elem.point(n),
2239 864 : nyplus = elem.point(n);
2240 864 : nyminus(1) -= delta_x;
2241 864 : nyplus(1) += delta_x;
2242 72 : typename GFunctor::FunctorValue gyminus =
2243 72 : g->eval_at_point(context,
2244 : var_component,
2245 : nyminus,
2246 864 : system.time,
2247 : true);
2248 72 : typename GFunctor::FunctorValue gyplus =
2249 72 : g->eval_at_point(context,
2250 : var_component,
2251 : nyplus,
2252 864 : system.time,
2253 : true);
2254 : // yz derivative
2255 936 : insert_id(first_id+6,
2256 864 : (grad_component(gyplus, 2) - grad_component(gyminus, 2)) / 2. / delta_x,
2257 : vertex.processor_id());
2258 : // Getting a 2nd order xyz is more tedious
2259 864 : Point nxmym = elem.point(n),
2260 864 : nxmyp = elem.point(n),
2261 864 : nxpym = elem.point(n),
2262 864 : nxpyp = elem.point(n);
2263 864 : nxmym(0) -= delta_x;
2264 864 : nxmym(1) -= delta_x;
2265 864 : nxmyp(0) -= delta_x;
2266 864 : nxmyp(1) += delta_x;
2267 864 : nxpym(0) += delta_x;
2268 864 : nxpym(1) -= delta_x;
2269 864 : nxpyp(0) += delta_x;
2270 864 : nxpyp(1) += delta_x;
2271 72 : typename GFunctor::FunctorValue gxmym =
2272 72 : g->eval_at_point(context,
2273 : var_component,
2274 : nxmym,
2275 864 : system.time,
2276 : true);
2277 72 : typename GFunctor::FunctorValue gxmyp =
2278 72 : g->eval_at_point(context,
2279 : var_component,
2280 : nxmyp,
2281 864 : system.time,
2282 : true);
2283 72 : typename GFunctor::FunctorValue gxpym =
2284 72 : g->eval_at_point(context,
2285 : var_component,
2286 : nxpym,
2287 864 : system.time,
2288 : true);
2289 72 : typename GFunctor::FunctorValue gxpyp =
2290 72 : g->eval_at_point(context,
2291 : var_component,
2292 : nxpyp,
2293 864 : system.time,
2294 : true);
2295 864 : FValue gxzplus = (grad_component(gxpyp, 2) - grad_component(gxmyp, 2))
2296 792 : / 2. / delta_x;
2297 864 : FValue gxzminus = (grad_component(gxpym, 2) - grad_component(gxmym, 2))
2298 792 : / 2. / delta_x;
2299 : // xyz derivative
2300 936 : insert_id(first_id+7,
2301 864 : (gxzplus - gxzminus) / 2. / delta_x,
2302 : vertex.processor_id());
2303 : }
2304 : #endif // LIBMESH_DIM > 2
2305 : }
2306 : #endif // LIBMESH_DIM > 1
2307 : }
2308 : else
2309 : {
2310 : // Currently other C_ONE elements have a single nodal
2311 : // value shape function and nodal gradient component
2312 : // shape functions
2313 918 : libmesh_assert_equal_to(
2314 : FEInterface::n_dofs_at_node(
2315 : base_fe_type,
2316 : &elem,
2317 : elem.get_node_index(&vertex),
2318 : base_fe_type.p_refinement),
2319 : (unsigned int)(1 + dim));
2320 :
2321 16543 : const FValue val =
2322 11853 : f.eval_at_node(context, var_component, dim,
2323 : vertex, extra_hanging_dofs[var],
2324 13707 : system.time);
2325 13707 : insert_id(first_id, val, vertex.processor_id());
2326 1836 : typename GFunctor::FunctorValue grad =
2327 12875 : is_old_vertex ?
2328 2082 : g->eval_at_node(context, var_component, dim,
2329 : vertex, extra_hanging_dofs[var],
2330 3284 : system.time) :
2331 9699 : g->eval_at_point(context, var_component, vertex,
2332 10423 : system.time, false);
2333 41121 : for (int i=0; i!= dim; ++i)
2334 29250 : insert_id(first_id + i + 1, grad.slice(i),
2335 : vertex.processor_id());
2336 : }
2337 : }
2338 : else
2339 0 : libmesh_error_msg("Unknown continuity " << cont);
2340 : }
2341 : }
2342 265122 : }
2343 :
2344 :
2345 : template <typename FFunctor, typename GFunctor,
2346 : typename FValue, typename ProjectionAction>
2347 257739 : void GenericProjector<FFunctor, GFunctor, FValue, ProjectionAction>::ProjectEdges::operator()
2348 : (const node_range & range)
2349 : {
2350 16326 : LOG_SCOPE ("project_edges","GenericProjector");
2351 :
2352 257739 : const unsigned int sys_num = system.number();
2353 :
2354 502294 : for (const auto & e_pair : range)
2355 : {
2356 244555 : const Elem & elem = *std::get<0>(e_pair.second);
2357 :
2358 : // If this is an unchanged element then we already copied all
2359 : // its dofs
2360 : #ifdef LIBMESH_ENABLE_AMR
2361 103562 : if (f.is_grid_projection() &&
2362 13692 : (elem.refinement_flag() != Elem::JUST_REFINED &&
2363 0 : elem.refinement_flag() != Elem::JUST_COARSENED &&
2364 0 : elem.p_refinement_flag() != Elem::JUST_REFINED &&
2365 0 : elem.p_refinement_flag() != Elem::JUST_COARSENED))
2366 0 : continue;
2367 : #endif // LIBMESH_ENABLE_AMR
2368 :
2369 244555 : const Node & edge_node = *e_pair.first;
2370 244555 : const int dim = elem.dim();
2371 18274 : const var_set & edge_vars = std::get<2>(e_pair.second);
2372 :
2373 244555 : const unsigned short edge_num = std::get<1>(e_pair.second);
2374 244555 : const unsigned short node_num = elem.n_vertices() + edge_num;
2375 244555 : context.edge = edge_num;
2376 244555 : context.pre_fe_reinit(system, &elem);
2377 :
2378 244555 : this->find_dofs_to_send(edge_node, elem, node_num, edge_vars);
2379 :
2380 : // Look at all the variables we're supposed to interpolate from
2381 : // this element on this edge
2382 489110 : for (const auto & var : edge_vars)
2383 : {
2384 244555 : const Variable & variable = system.variable(var);
2385 18274 : const FEType & base_fe_type = variable.type();
2386 54822 : const unsigned int var_component =
2387 244555 : system.variable_scalar_number(var, 0);
2388 :
2389 244555 : if (base_fe_type.family == SCALAR)
2390 131941 : continue;
2391 :
2392 244555 : FEType fe_type = base_fe_type;
2393 :
2394 : // This may be a p refined element
2395 36548 : fe_type.order = fe_type.order + elem.p_level();
2396 :
2397 : // If this is a Lagrange element with DoFs on edges then by
2398 : // convention we interpolate at the node rather than project
2399 : // along the edge.
2400 244555 : if (fe_type.family == LAGRANGE || fe_type.family == LAGRANGE_VEC)
2401 : {
2402 131941 : if (fe_type.order > 1)
2403 : {
2404 : const processor_id_type pid =
2405 20288 : edge_node.processor_id();
2406 136192 : FValue fval = f.eval_at_point
2407 131941 : (context, var_component, edge_node, system.time,
2408 : false);
2409 131941 : if (fe_type.family == LAGRANGE_VEC)
2410 : {
2411 : // We will have a number of nodal value DoFs equal to the elem dim
2412 142392 : for (auto i : make_range(elem.dim()))
2413 : {
2414 : const dof_id_type dof_id =
2415 106794 : edge_node.dof_number(sys_num, var, i);
2416 :
2417 : // Need this conversion so that this method
2418 : // will compile for TYPE_SCALAR instantiations
2419 100350 : const auto insert_val =
2420 13896 : raw_value<typename ProjectionAction::InsertInput>(fval, i);
2421 :
2422 106794 : insert_id(dof_id, insert_val, pid);
2423 : }
2424 : }
2425 : else // We are LAGRANGE
2426 : {
2427 : const dof_id_type dof_id =
2428 96343 : edge_node.dof_number(sys_num, var, 0);
2429 96343 : insert_id(dof_id, fval, pid);
2430 : }
2431 : }
2432 131941 : continue;
2433 111653 : }
2434 :
2435 : #ifdef LIBMESH_ENABLE_AMR
2436 : // If this is a low order monomial element which has merely
2437 : // been h refined then we already copied all its dofs
2438 0 : if (fe_type.family == MONOMIAL &&
2439 0 : fe_type.order == CONSTANT &&
2440 112614 : elem.refinement_flag() != Elem::JUST_COARSENED &&
2441 0 : elem.p_refinement_flag() != Elem::JUST_COARSENED)
2442 0 : continue;
2443 : #endif // LIBMESH_ENABLE_AMR
2444 :
2445 : // FIXME: Need to generalize this to vector-valued elements. [PB]
2446 8130 : FEGenericBase<typename FFunctor::RealType> * fe = nullptr;
2447 112614 : context.get_element_fe( var, fe, dim );
2448 8130 : FEGenericBase<typename FFunctor::RealType> * edge_fe = nullptr;
2449 8130 : context.get_edge_fe( var, edge_fe );
2450 :
2451 : // If we're JUST_COARSENED we'll need a custom
2452 : // evaluation, not just the standard edge FE
2453 112614 : const FEGenericBase<typename FFunctor::RealType> & proj_fe =
2454 : #ifdef LIBMESH_ENABLE_AMR
2455 16260 : (elem.refinement_flag() == Elem::JUST_COARSENED) ?
2456 : *fe :
2457 : #endif
2458 : *edge_fe;
2459 :
2460 : #ifdef LIBMESH_ENABLE_AMR
2461 112614 : if (elem.refinement_flag() == Elem::JUST_COARSENED)
2462 : {
2463 0 : std::vector<Point> fine_points;
2464 :
2465 0 : std::unique_ptr<FEGenericBase<typename FFunctor::RealType>> fine_fe
2466 : (FEGenericBase<typename FFunctor::RealType>::build (dim, base_fe_type));
2467 :
2468 0 : std::unique_ptr<QBase> qrule
2469 : (base_fe_type.default_quadrature_rule(1));
2470 0 : fine_fe->attach_quadrature_rule(qrule.get());
2471 :
2472 0 : const std::vector<Point> & child_xyz =
2473 0 : fine_fe->get_xyz();
2474 :
2475 0 : for (unsigned int c = 0, nc = elem.n_children();
2476 0 : c != nc; ++c)
2477 : {
2478 0 : if (!elem.is_child_on_edge(c, context.edge))
2479 0 : continue;
2480 :
2481 0 : fine_fe->edge_reinit(elem.child_ptr(c), context.edge);
2482 0 : fine_points.insert(fine_points.end(),
2483 : child_xyz.begin(),
2484 : child_xyz.end());
2485 : }
2486 :
2487 0 : std::vector<Point> fine_qp;
2488 0 : FEMap::inverse_map (dim, &elem, fine_points, fine_qp);
2489 :
2490 0 : context.elem_fe_reinit(&fine_qp);
2491 0 : }
2492 : else
2493 : #endif // LIBMESH_ENABLE_AMR
2494 112614 : context.edge_fe_reinit();
2495 :
2496 16260 : const std::vector<dof_id_type> & dof_indices =
2497 96354 : context.get_dof_indices(var);
2498 :
2499 16260 : std::vector<unsigned int> edge_dofs;
2500 112614 : FEInterface::dofs_on_edge(&elem, dim, base_fe_type,
2501 112614 : context.edge, edge_dofs);
2502 :
2503 16260 : this->construct_projection
2504 96354 : (dof_indices, edge_dofs, var_component,
2505 : &edge_node, proj_fe);
2506 : }
2507 : }
2508 257739 : }
2509 :
2510 :
2511 : template <typename FFunctor, typename GFunctor,
2512 : typename FValue, typename ProjectionAction>
2513 261562 : void GenericProjector<FFunctor, GFunctor, FValue, ProjectionAction>::ProjectSides::operator()
2514 : (const node_range & range)
2515 : {
2516 18912 : LOG_SCOPE ("project_sides","GenericProjector");
2517 :
2518 261562 : const unsigned int sys_num = system.number();
2519 :
2520 2317880 : for (const auto & s_pair : range)
2521 : {
2522 2056318 : const Elem & elem = *std::get<0>(s_pair.second);
2523 :
2524 : // If this is an unchanged element then we already copied all
2525 : // its dofs
2526 : #ifdef LIBMESH_ENABLE_AMR
2527 1883835 : if (f.is_grid_projection() &&
2528 426982 : (elem.refinement_flag() != Elem::JUST_REFINED &&
2529 14467 : elem.refinement_flag() != Elem::JUST_COARSENED &&
2530 308 : elem.p_refinement_flag() != Elem::JUST_REFINED &&
2531 20 : elem.p_refinement_flag() != Elem::JUST_COARSENED))
2532 0 : continue;
2533 : #endif // LIBMESH_ENABLE_AMR
2534 :
2535 2056318 : const Node & side_node = *s_pair.first;
2536 2056318 : const int dim = elem.dim();
2537 167888 : const var_set & side_vars = std::get<2>(s_pair.second);
2538 :
2539 2056318 : const unsigned int side_num = std::get<1>(s_pair.second);
2540 2056318 : unsigned short node_num = elem.n_vertices()+side_num;
2541 : // In 3D only some sides may have nodes
2542 2056318 : if (dim == 3)
2543 5091294 : for (auto n : make_range(elem.n_nodes()))
2544 : {
2545 5091294 : if (!elem.is_face(n))
2546 3856492 : continue;
2547 :
2548 904677 : if (elem.is_node_on_side(n, side_num))
2549 : {
2550 322626 : node_num = n;
2551 322626 : break;
2552 : }
2553 : }
2554 :
2555 2056318 : context.side = side_num;
2556 2056318 : context.pre_fe_reinit(system, &elem);
2557 :
2558 2056318 : this->find_dofs_to_send(side_node, elem, node_num, side_vars);
2559 :
2560 : // Look at all the variables we're supposed to interpolate from
2561 : // this element on this side
2562 4162181 : for (const auto & var : side_vars)
2563 : {
2564 2105863 : const Variable & variable = system.variable(var);
2565 172091 : const FEType & base_fe_type = variable.type();
2566 516309 : const unsigned int var_component =
2567 2105863 : system.variable_scalar_number(var, 0);
2568 :
2569 2105863 : if (base_fe_type.family == SCALAR)
2570 1494554 : continue;
2571 :
2572 2105863 : FEType fe_type = base_fe_type;
2573 2105863 : const bool add_p_level = base_fe_type.p_refinement;
2574 :
2575 : // This may be a p refined element
2576 2277990 : fe_type.order = fe_type.order + add_p_level*elem.p_level();
2577 :
2578 : // If this is a Lagrange element with DoFs on sides then by
2579 : // convention we interpolate at the node rather than project
2580 : // along the side.
2581 2105863 : if (fe_type.family == LAGRANGE || fe_type.family == LAGRANGE_VEC)
2582 : {
2583 : // If order==1 then we only have DoFs on vertices, so
2584 : // skip all this.
2585 : // If order>1 ... we might still have something tricky
2586 : // like order==2 PRISM20, where some side nodes have
2587 : // DoFs but others don't. Gotta check.
2588 2989108 : if (fe_type.order > 1 &&
2589 1494554 : side_node.n_comp(sys_num, var))
2590 : {
2591 : const processor_id_type pid =
2592 245984 : side_node.processor_id();
2593 1946694 : FValue fval = f.eval_at_point
2594 1489775 : (context, var_component, side_node, system.time,
2595 : false);
2596 :
2597 1489775 : if (fe_type.family == LAGRANGE_VEC)
2598 : {
2599 : // We will have a number of nodal value DoFs equal to the elem dim
2600 84492 : for (auto i : make_range(elem.dim()))
2601 : {
2602 62736 : const dof_id_type dof_id = side_node.dof_number(sys_num, var, i);
2603 :
2604 : // Need this conversion so that this method
2605 : // will compile for TYPE_SCALAR instantiations
2606 58188 : const auto insert_val =
2607 9405 : raw_value<typename ProjectionAction::InsertInput>(fval, i);
2608 :
2609 62736 : insert_id(dof_id, insert_val, pid);
2610 : }
2611 : }
2612 : else // We are LAGRANGE
2613 : {
2614 : const dof_id_type dof_id =
2615 1468019 : side_node.dof_number(sys_num, var, 0);
2616 1468019 : insert_id(dof_id, fval, pid);
2617 : }
2618 : }
2619 1494554 : continue;
2620 1247838 : }
2621 :
2622 : #ifdef LIBMESH_ENABLE_AMR
2623 : // If this is a low order monomial element which has merely
2624 : // been h refined then we already copied all its dofs
2625 0 : if (fe_type.family == MONOMIAL &&
2626 0 : fe_type.order == CONSTANT &&
2627 611309 : elem.refinement_flag() != Elem::JUST_COARSENED &&
2628 0 : elem.p_refinement_flag() != Elem::JUST_COARSENED)
2629 0 : continue;
2630 : #endif // LIBMESH_ENABLE_AMR
2631 :
2632 : // FIXME: Need to generalize this to vector-valued elements. [PB]
2633 48751 : FEGenericBase<typename FFunctor::RealType> * fe = nullptr;
2634 611309 : context.get_element_fe( var, fe, dim );
2635 562558 : FEGenericBase<typename FFunctor::RealType> * side_fe = nullptr;
2636 562558 : context.get_side_fe( var, side_fe );
2637 :
2638 : // If we're JUST_COARSENED we'll need a custom
2639 : // evaluation, not just the standard side FE
2640 611309 : const FEGenericBase<typename FFunctor::RealType> & proj_fe =
2641 : #ifdef LIBMESH_ENABLE_AMR
2642 97502 : (elem.refinement_flag() == Elem::JUST_COARSENED) ?
2643 : *fe :
2644 : #endif
2645 : *side_fe;
2646 :
2647 : #ifdef LIBMESH_ENABLE_AMR
2648 611309 : if (elem.refinement_flag() == Elem::JUST_COARSENED)
2649 : {
2650 9252 : std::vector<Point> fine_points;
2651 :
2652 60105 : std::unique_ptr<FEGenericBase<typename FFunctor::RealType>> fine_fe
2653 : (FEGenericBase<typename FFunctor::RealType>::build (dim, base_fe_type));
2654 :
2655 60105 : std::unique_ptr<QBase> qrule
2656 : (base_fe_type.default_quadrature_rule(1));
2657 60105 : fine_fe->attach_quadrature_rule(qrule.get());
2658 :
2659 9252 : const std::vector<Point> & child_xyz =
2660 4626 : fine_fe->get_xyz();
2661 :
2662 277395 : for (unsigned int c = 0, nc = elem.n_children();
2663 277395 : c != nc; ++c)
2664 : {
2665 221916 : if (!elem.is_child_on_side(c, context.side))
2666 101706 : continue;
2667 :
2668 120210 : fine_fe->reinit(elem.child_ptr(c), context.side);
2669 120210 : fine_points.insert(fine_points.end(),
2670 : child_xyz.begin(),
2671 : child_xyz.end());
2672 : }
2673 :
2674 9252 : std::vector<Point> fine_qp;
2675 55479 : FEMap::inverse_map (dim, &elem, fine_points, fine_qp);
2676 :
2677 55479 : context.elem_fe_reinit(&fine_qp);
2678 46227 : }
2679 : else
2680 : #endif // LIBMESH_ENABLE_AMR
2681 555830 : context.side_fe_reinit();
2682 :
2683 97502 : const std::vector<dof_id_type> & dof_indices =
2684 513807 : context.get_dof_indices(var);
2685 :
2686 97502 : std::vector<unsigned int> side_dofs;
2687 708811 : FEInterface::dofs_on_side(&elem, dim, base_fe_type,
2688 611309 : context.side, side_dofs, add_p_level);
2689 :
2690 97502 : this->construct_projection
2691 513807 : (dof_indices, side_dofs, var_component,
2692 : &side_node, proj_fe);
2693 : }
2694 : }
2695 261562 : }
2696 :
2697 :
2698 : template <typename FFunctor, typename GFunctor,
2699 : typename FValue, typename ProjectionAction>
2700 258743 : void GenericProjector<FFunctor, GFunctor, FValue, ProjectionAction>::ProjectInteriors::operator()
2701 : (const interior_range & range)
2702 : {
2703 17034 : LOG_SCOPE ("project_interiors","GenericProjector");
2704 :
2705 258743 : const unsigned int sys_num = system.number();
2706 :
2707 : // Iterate over all dof-bearing element interiors in the range
2708 1169339 : for (const auto & elem : range)
2709 : {
2710 910596 : unsigned char dim = cast_int<unsigned char>(elem->dim());
2711 :
2712 910596 : context.pre_fe_reinit(system, elem);
2713 :
2714 : // Loop over all the variables we've been requested to project, to
2715 : // do the projection
2716 1871196 : for (const auto & var : this->projector.variables)
2717 : {
2718 960600 : const Variable & variable = system.variable(var);
2719 :
2720 960600 : if (!variable.active_on_subdomain(elem->subdomain_id()))
2721 680202 : continue;
2722 :
2723 81115 : const FEType & base_fe_type = variable.type();
2724 :
2725 957096 : if (base_fe_type.family == SCALAR)
2726 0 : continue;
2727 :
2728 81115 : FEGenericBase<typename FFunctor::RealType> * fe = nullptr;
2729 957096 : context.get_element_fe( var, fe, dim );
2730 :
2731 957096 : FEType fe_type = base_fe_type;
2732 957096 : const bool add_p_level = fe_type.p_refinement;
2733 :
2734 : // This may be a p refined element
2735 1038221 : fe_type.order = fe_type.order + add_p_level * elem->p_level();
2736 :
2737 243355 : const unsigned int var_component =
2738 957096 : system.variable_scalar_number(var, 0);
2739 :
2740 : // If this is a Lagrange element with interior DoFs then by
2741 : // convention we interpolate at the interior node rather
2742 : // than project along the interior.
2743 957096 : if (fe_type.family == LAGRANGE || fe_type.family == LAGRANGE_VEC)
2744 : {
2745 676698 : if (fe_type.order > 1)
2746 : {
2747 654728 : const unsigned int first_interior_node =
2748 654728 : (elem->n_vertices() +
2749 654728 : ((elem->dim() > 2) * elem->n_edges()) +
2750 654728 : ((elem->dim() > 1) * elem->n_sides()));
2751 654728 : const unsigned int n_nodes = elem->n_nodes();
2752 :
2753 : // < vs != is important here for HEX20, QUAD8!
2754 1309456 : for (unsigned int n = first_interior_node; n < n_nodes; ++n)
2755 : {
2756 654728 : const Node & interior_node = elem->node_ref(n);
2757 :
2758 868297 : FValue fval = f.eval_at_point
2759 573349 : (context, var_component, interior_node,
2760 654728 : system.time, false);
2761 : const processor_id_type pid =
2762 110906 : interior_node.processor_id();
2763 :
2764 654728 : if (fe_type.family == LAGRANGE_VEC)
2765 : {
2766 : // We will have a number of nodal value DoFs equal to the elem dim
2767 2448 : for (unsigned int i = 0; i < elem->dim(); ++i)
2768 : {
2769 1728 : const dof_id_type dof_id = interior_node.dof_number(sys_num, var, i);
2770 :
2771 : // Need this conversion so that this method
2772 : // will compile for TYPE_SCALAR instantiations
2773 1584 : const auto insert_val =
2774 288 : raw_value<typename ProjectionAction::InsertInput>(fval, i);
2775 :
2776 1728 : insert_id(dof_id, insert_val, pid);
2777 : }
2778 : }
2779 : else // We are LAGRANGE
2780 : {
2781 110786 : const dof_id_type dof_id =
2782 543222 : interior_node.dof_number(sys_num, var, 0);
2783 654008 : insert_id(dof_id, fval, pid);
2784 : }
2785 : }
2786 : }
2787 676698 : continue;
2788 561886 : }
2789 :
2790 : #ifdef LIBMESH_ENABLE_AMR
2791 280398 : if (elem->refinement_flag() == Elem::JUST_COARSENED)
2792 : {
2793 4242 : std::vector<Point> fine_points;
2794 :
2795 26106 : std::unique_ptr<FEGenericBase<typename FFunctor::RealType>> fine_fe
2796 : (FEGenericBase<typename FFunctor::RealType>::build (dim, base_fe_type));
2797 :
2798 26106 : std::unique_ptr<QBase> qrule
2799 : (base_fe_type.default_quadrature_rule(dim));
2800 26106 : fine_fe->attach_quadrature_rule(qrule.get());
2801 :
2802 4242 : const std::vector<Point> & child_xyz =
2803 2145 : fine_fe->get_xyz();
2804 :
2805 123198 : for (auto & child : elem->child_ref_range())
2806 : {
2807 97092 : fine_fe->reinit(&child);
2808 105672 : fine_points.insert(fine_points.end(),
2809 : child_xyz.begin(),
2810 : child_xyz.end());
2811 : }
2812 :
2813 4242 : std::vector<Point> fine_qp;
2814 23985 : FEMap::inverse_map (dim, elem, fine_points, fine_qp);
2815 :
2816 23985 : context.elem_fe_reinit(&fine_qp);
2817 19743 : }
2818 : else
2819 : #endif // LIBMESH_ENABLE_AMR
2820 256413 : context.elem_fe_reinit();
2821 :
2822 47428 : const std::vector<dof_id_type> & dof_indices =
2823 232970 : context.get_dof_indices(var);
2824 :
2825 : const unsigned int n_dofs =
2826 47428 : cast_int<unsigned int>(dof_indices.size());
2827 :
2828 304113 : std::vector<unsigned int> all_dofs(n_dofs);
2829 23715 : std::iota(all_dofs.begin(), all_dofs.end(), 0);
2830 :
2831 47428 : this->construct_projection
2832 232970 : (dof_indices, all_dofs, var_component,
2833 : nullptr, *fe);
2834 : } // end variables loop
2835 : } // end elements loop
2836 258743 : }
2837 :
2838 :
2839 :
2840 : template <typename FFunctor, typename GFunctor,
2841 : typename FValue, typename ProjectionAction>
2842 : void
2843 9198609 : GenericProjector<FFunctor, GFunctor, FValue,
2844 : ProjectionAction>::SubFunctor::find_dofs_to_send
2845 : (const Node & node, const Elem & elem, unsigned short node_num, const var_set & vars)
2846 : {
2847 874796 : libmesh_assert (&node == elem.node_ptr(node_num));
2848 :
2849 : // Any ghosted node in our range might have an owner who needs our
2850 : // data
2851 1749375 : const processor_id_type owner = node.processor_id();
2852 10073188 : if (owner != system.processor_id())
2853 : {
2854 77612 : const MeshBase & mesh = system.get_mesh();
2855 38839 : const DofMap & dof_map = system.get_dof_map();
2856 :
2857 : // But let's check and see if we can be certain the owner can
2858 : // compute any or all of its own dof coefficients on that node.
2859 77678 : std::vector<dof_id_type> node_dof_ids, patch_dof_ids;
2860 1925568 : for (const auto & var : vars)
2861 : {
2862 991573 : const Variable & variable = system.variable(var);
2863 :
2864 991573 : if (!variable.active_on_subdomain(elem.subdomain_id()))
2865 0 : continue;
2866 :
2867 991573 : dof_map.dof_indices(elem, node_num, node_dof_ids, var);
2868 : }
2869 38839 : libmesh_assert(std::is_sorted(node_dof_ids.begin(),
2870 : node_dof_ids.end()));
2871 77678 : const std::vector<dof_id_type> & patch =
2872 933995 : (*this->projector.nodes_to_elem)[node.id()];
2873 2621236 : for (const auto & elem_id : patch)
2874 : {
2875 2599447 : const Elem & patch_elem = mesh.elem_ref(elem_id);
2876 1451844 : if (!patch_elem.active() || owner != patch_elem.processor_id())
2877 1651007 : continue;
2878 948440 : dof_map.dof_indices(&patch_elem, patch_dof_ids);
2879 948440 : std::sort(patch_dof_ids.begin(), patch_dof_ids.end());
2880 :
2881 : // Remove any node_dof_ids that we see can be calculated on
2882 : // this element
2883 987909 : std::vector<dof_id_type> diff_ids(node_dof_ids.size());
2884 948440 : auto it = std::set_difference(node_dof_ids.begin(), node_dof_ids.end(),
2885 : patch_dof_ids.begin(), patch_dof_ids.end(), diff_ids.begin());
2886 948440 : diff_ids.resize(it-diff_ids.begin());
2887 39535 : node_dof_ids.swap(diff_ids);
2888 948440 : if (node_dof_ids.empty())
2889 37889 : break;
2890 : }
2891 :
2892 : // Give ids_to_push default invalid pid: not yet computed
2893 969362 : for (auto id : node_dof_ids)
2894 35367 : new_ids_to_push[id].second = DofObject::invalid_processor_id;
2895 : }
2896 9198609 : }
2897 :
2898 :
2899 :
2900 : template <typename FFunctor, typename GFunctor,
2901 : typename FValue, typename ProjectionAction>
2902 : template <typename Value>
2903 : void
2904 768009 : GenericProjector<FFunctor, GFunctor, FValue,
2905 : ProjectionAction>::send_and_insert_dof_values
2906 : (std::unordered_map<dof_id_type, std::pair<Value, processor_id_type>> & ids_to_push,
2907 : ProjectionAction & action) const
2908 : {
2909 45504 : LOG_SCOPE ("send_and_insert_dof_values", "GenericProjector");
2910 :
2911 : // See if we calculated any ids that need to be pushed; get them
2912 : // ready to push.
2913 : std::unordered_map<processor_id_type, std::vector<dof_id_type>>
2914 45504 : pushed_dof_ids, received_dof_ids;
2915 : std::unordered_map<processor_id_type, std::vector<typename TypeToSend<Value>::type>>
2916 45504 : pushed_dof_values, received_dof_values;
2917 864489 : for (auto & id_val_pid : ids_to_push)
2918 : {
2919 96480 : processor_id_type pid = id_val_pid.second.second;
2920 96480 : if (pid != DofObject::invalid_processor_id)
2921 : {
2922 24600 : pushed_dof_ids[pid].push_back(id_val_pid.first);
2923 24600 : pushed_dof_values[pid].push_back(convert_to_send(id_val_pid.second.first));
2924 : }
2925 : }
2926 :
2927 : // If and when we get ids pushed to us, act on them if we have
2928 : // corresponding values or save them if not
2929 768333 : auto ids_action_functor =
2930 8418 : [&action, &received_dof_ids, &received_dof_values]
2931 : (processor_id_type pid,
2932 324 : const std::vector<dof_id_type> & data)
2933 : {
2934 162 : auto iter = received_dof_values.find(pid);
2935 8742 : if (iter == received_dof_values.end())
2936 : {
2937 162 : libmesh_assert(received_dof_ids.find(pid) ==
2938 : received_dof_ids.end());
2939 8742 : received_dof_ids[pid] = data;
2940 : }
2941 : else
2942 : {
2943 0 : auto & vals = iter->second;
2944 0 : std::size_t vals_size = vals.size();
2945 0 : libmesh_assert_equal_to(vals_size, data.size());
2946 0 : for (std::size_t i=0; i != vals_size; ++i)
2947 : {
2948 0 : Value val;
2949 0 : convert_from_receive(vals[i], val);
2950 0 : action.insert(data[i], val);
2951 : }
2952 0 : received_dof_values.erase(iter);
2953 : }
2954 : };
2955 :
2956 : // If and when we get values pushed to us, act on them if we have
2957 : // corresponding ids or save them if not
2958 768333 : auto values_action_functor =
2959 38916 : [&action, &received_dof_ids, &received_dof_values]
2960 : (processor_id_type pid,
2961 324 : const std::vector<typename TypeToSend<Value>::type> & data)
2962 : {
2963 162 : auto iter = received_dof_ids.find(pid);
2964 8742 : if (iter == received_dof_ids.end())
2965 : {
2966 : // We get no more than 1 values vector from anywhere
2967 0 : libmesh_assert(received_dof_values.find(pid) ==
2968 : received_dof_values.end());
2969 0 : received_dof_values[pid] = data;
2970 : }
2971 : else
2972 : {
2973 162 : auto & ids = iter->second;
2974 324 : std::size_t ids_size = ids.size();
2975 162 : libmesh_assert_equal_to(ids_size, data.size());
2976 33342 : for (std::size_t i=0; i != ids_size; ++i)
2977 : {
2978 0 : Value val;
2979 2520 : convert_from_receive(data[i], val);
2980 25860 : action.insert(ids[i], val);
2981 : }
2982 162 : received_dof_ids.erase(iter);
2983 : }
2984 : };
2985 :
2986 : Parallel::push_parallel_vector_data
2987 768009 : (system.comm(), pushed_dof_ids, ids_action_functor);
2988 :
2989 : Parallel::push_parallel_vector_data
2990 768009 : (system.comm(), pushed_dof_values, values_action_functor);
2991 :
2992 : // At this point we shouldn't have any unprocessed data left
2993 22752 : libmesh_assert(received_dof_ids.empty());
2994 22752 : libmesh_assert(received_dof_values.empty());
2995 :
2996 768009 : }
2997 :
2998 :
2999 :
3000 : template <typename FFunctor, typename GFunctor,
3001 : typename FValue, typename ProjectionAction>
3002 : void
3003 1004321 : GenericProjector<FFunctor, GFunctor, FValue, ProjectionAction>::SubProjector::construct_projection
3004 : (const std::vector<dof_id_type> & dof_indices_var,
3005 : const std::vector<unsigned int> & involved_dofs,
3006 : unsigned int var_component,
3007 : const Node * node,
3008 : const FEGenericBase<typename FFunctor::RealType> & fe)
3009 : {
3010 1004321 : const auto & JxW = fe.get_JxW();
3011 80596 : const auto & phi = fe.get_phi();
3012 80596 : const std::vector<std::vector<typename FEGenericBase<typename FFunctor::RealType>::OutputGradient>> * dphi = nullptr;
3013 208786 : const std::vector<Point> & xyz_values = fe.get_xyz();
3014 1004321 : const FEContinuity cont = fe.get_continuity();
3015 80596 : const std::unordered_map<dof_id_type, typename FFunctor::ValuePushType> & ids_to_save =
3016 1004321 : this->projector.ids_to_save;
3017 :
3018 1004321 : if (cont == C_ONE)
3019 2806 : dphi = &(fe.get_dphi());
3020 :
3021 : const unsigned int n_involved_dofs =
3022 161190 : cast_int<unsigned int>(involved_dofs.size());
3023 :
3024 80596 : std::vector<dof_id_type> free_dof_ids;
3025 923727 : DenseVector<typename FFunctor::ValuePushType> Uinvolved(n_involved_dofs);
3026 1004321 : std::vector<char> dof_is_fixed(n_involved_dofs, false); // bools
3027 :
3028 9350574 : for (auto i : make_range(n_involved_dofs))
3029 : {
3030 9703017 : const dof_id_type id = dof_indices_var[involved_dofs[i]];
3031 678408 : auto iter = ids_to_save.find(id);
3032 8346253 : if (iter == ids_to_save.end())
3033 4741846 : free_dof_ids.push_back(id);
3034 : else
3035 : {
3036 3604407 : dof_is_fixed[i] = true;
3037 3901149 : Uinvolved(i) = iter->second;
3038 : }
3039 : }
3040 :
3041 1004321 : const unsigned int free_dofs = free_dof_ids.size();
3042 :
3043 : // There may be nothing to project
3044 1004321 : if (!free_dofs)
3045 118 : return;
3046 :
3047 : // The element matrix and RHS for projections.
3048 : // Note that Ke is always real-valued, whereas
3049 : // Fe may be complex valued if complex number
3050 : // support is enabled
3051 1164057 : DenseMatrix<Real> Ke(free_dofs, free_dofs);
3052 1003105 : DenseVector<typename FFunctor::ValuePushType> Fe(free_dofs);
3053 : // The new degree of freedom coefficients to solve for
3054 1003105 : DenseVector<typename FFunctor::ValuePushType> Ufree(free_dofs);
3055 :
3056 : const unsigned int n_qp =
3057 160954 : cast_int<unsigned int>(xyz_values.size());
3058 :
3059 : // Loop over the quadrature points
3060 12705333 : for (unsigned int qp=0; qp<n_qp; qp++)
3061 : {
3062 : // solution at the quadrature point
3063 11702230 : FValue fineval = f.eval_at_point(context,
3064 : var_component,
3065 : xyz_values[qp],
3066 11702230 : system.time,
3067 : false);
3068 : // solution grad at the quadrature point
3069 961601 : typename GFunctor::FunctorValue finegrad;
3070 11702230 : if (cont == C_ONE)
3071 158847 : finegrad = g->eval_at_point(context,
3072 : var_component,
3073 : xyz_values[qp],
3074 147174 : system.time,
3075 : false);
3076 :
3077 : // Form edge projection matrix
3078 229418392 : for (unsigned int sidei=0, freei=0;
3079 231003706 : sidei != n_involved_dofs; ++sidei)
3080 : {
3081 219301476 : unsigned int i = involved_dofs[sidei];
3082 : // fixed DoFs aren't test functions
3083 237456158 : if (dof_is_fixed[sidei])
3084 75457824 : continue;
3085 4997928672 : for (unsigned int sidej=0, freej=0;
3086 5123617488 : sidej != n_involved_dofs; ++sidej)
3087 : {
3088 4986664160 : unsigned int j = involved_dofs[sidej];
3089 5400324120 : if (dof_is_fixed[sidej])
3090 657767204 : Fe(freei) -= phi[i][qp] * phi[j][qp] *
3091 444996952 : JxW[qp] * Uinvolved(sidej);
3092 : else
3093 5679552392 : Ke(freei,freej) += phi[i][qp] *
3094 5679551674 : phi[j][qp] * JxW[qp];
3095 4986664160 : if (cont == C_ONE)
3096 : {
3097 1301847 : if (dof_is_fixed[sidej])
3098 1225224 : Fe(freei) -= ( TensorTools::inner_product((*dphi)[i][qp],
3099 1146224 : (*dphi)[j][qp]) ) *
3100 989772 : JxW[qp] * Uinvolved(sidej);
3101 : else
3102 268599 : Ke(freei,freej) += ( TensorTools::inner_product((*dphi)[i][qp],
3103 234007 : (*dphi)[j][qp]) )
3104 234007 : * JxW[qp];
3105 : }
3106 5400324120 : if (!dof_is_fixed[sidej])
3107 4548246016 : freej++;
3108 : }
3109 170746402 : Fe(freei) += phi[i][qp] * fineval * JxW[qp];
3110 136953328 : if (cont == C_ONE)
3111 174815 : Fe(freei) += (TensorTools::inner_product(finegrad,
3112 187539 : (*dphi)[i][qp]) ) *
3113 12921 : JxW[qp];
3114 136953328 : freei++;
3115 : }
3116 : }
3117 :
3118 1003103 : Ke.cholesky_solve(Fe, Ufree);
3119 :
3120 : // Transfer new edge solutions to element
3121 1003103 : const processor_id_type pid = node ?
3122 113574 : node->processor_id() : DofObject::invalid_processor_id;
3123 1003103 : insert_ids(free_dof_ids, Ufree.get_values(), pid);
3124 842149 : }
3125 :
3126 :
3127 : } // namespace libMesh
3128 :
3129 : #endif // GENERIC_PROJECTOR_H
|
