libMesh
generic_projector.h
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 #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 
58 template <typename T>
59 struct TypeToSend {
60  typedef T type;
61 };
62 
63 template <typename T>
64 const typename TypeToSend<T>::type convert_to_send (const T& in)
65 { return in; }
66 
67 template <typename SendT, typename T>
68 void convert_from_receive (SendT & received, T & converted)
69 { converted = received; }
70 
81 template <typename FFunctor, typename GFunctor,
82  typename FValue, typename ProjectionAction>
84 {
85 public:
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 
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 
118 
120 
121 public:
122  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  system(system_in),
129  master_f(f_in),
130  master_g(g_in),
131  master_action(act_in),
132  variables(variables_in),
133  nodes_to_elem(nodes_to_elem_in)
134  {
135  if (!nodes_to_elem_in)
136  {
139  }
140  }
141 
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),
148  variables(in.variables),
150  {}
151 
152  ~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 {
171 
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,
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<
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<
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<
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<
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.
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>>
248 
251 
252 
253  struct SubProjector : public SubFunctor {
255 
256  using SubFunctor::action;
257  using SubFunctor::f;
258  using SubFunctor::system;
259  using SubFunctor::context;
260  using SubFunctor::conts;
262  using SubFunctor::insert_id;
264 
265  protected:
266  // Projections of C1 elements require a gradient as well
267  std::unique_ptr<GFunctor> g;
268 
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,
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 {
282 
284 
285  using SubFunctor::action;
286  using SubFunctor::f;
287  using SubFunctor::system;
288  using SubFunctor::context;
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 
303  std::vector<const Elem *> interiors;
304  };
305 
306 
307  struct ProjectVertices : public SubProjector {
309 
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;
319 
322 
323  void operator() (const node_range & range);
324  };
325 
326 
327  struct ProjectEdges : public SubProjector {
329 
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;
339 
342 
343  void operator() (const node_range & range);
344  };
345 
346  struct ProjectSides : public SubProjector {
348 
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;
358 
361 
362  void operator() (const node_range & range);
363  };
364 
365 
368 
369  struct ProjectInteriors : public SubProjector {
371 
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;
381 
384 
385  void operator() (const interior_range & range);
386  };
387 
388  template <typename Value>
390  (std::unordered_map<dof_id_type, std::pair<Value, processor_id_type>> & ids_to_push,
391  ProjectionAction & action) const;
392 };
393 
394 
404 template <typename Val>
406 {
407 private:
409 
410 public:
411  typedef Val InsertInput;
412 
414  target_vector(target_vec) {}
415 
417  Val val)
418  {
419  // Lock the new vector since it is shared among threads.
420  {
421  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
422  target_vector.set(id, val);
423  }
424  }
425 
426 
427  void insert(const std::vector<dof_id_type> & dof_indices,
428  const DenseVector<Val> & Ue)
429  {
430  const numeric_index_type
433 
434  unsigned int size = Ue.size();
435 
436  libmesh_assert_equal_to(size, dof_indices.size());
437 
438  // Lock the new vector since it is shared among threads.
439  {
440  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
441 
442  for (unsigned int i = 0; i != size; ++i)
443  if ((dof_indices[i] >= first) && (dof_indices[i] < last))
444  target_vector.set(dof_indices[i], Ue(i));
445  }
446  }
447 
448 };
449 
450 
459 template <typename Output>
461 {
462 protected:
464 
465 public:
468  typedef Output FunctorValue;
469 
470  FEMFunctionWrapper(const FEMFunctionBase<Output> & f) : _f(f.clone()) {}
471 
473  _f(fw._f->clone()) {}
474 
475  void init_context (FEMContext & c) { _f->init_context(c); }
476 
477  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  { return _f->component(c, i, n, time); }
484 
485  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  { return _f->component(c, i, n, time); }
491 
492  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  { libmesh_error(); } // this is only for grid projections
498 
499  bool is_grid_projection() { return false; }
500 
501  void eval_old_dofs (const Elem &,
502  unsigned int,
503  unsigned int,
504  std::vector<dof_id_type> &,
505  std::vector<Output> &)
506  { libmesh_error(); }
507 
508  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  { libmesh_error(); }
515 
516 private:
517  std::unique_ptr<FEMFunctionBase<Output>> _f;
518 };
519 
520 
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>
538 {
539 protected:
541 public:
543 
545  const std::vector<unsigned int> * vars) :
546  last_elem(nullptr),
547  sys(sys_in),
548  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  auto make_lookup = [this](unsigned int v)
553  {
554  const unsigned int vcomp = sys.variable_scalar_number(v,0);
555  if (vcomp >= component_to_var.size())
556  component_to_var.resize(vcomp+1, static_cast<unsigned int>(-1));
557  component_to_var[vcomp] = v;
558  };
559 
560  if (vars)
561  for (auto v : *vars)
562  make_lookup(v);
563  else
564  for (auto v : make_range(sys.n_vars()))
565  make_lookup(v);
566  }
567 
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.
581  {
583 
584  const std::set<unsigned char> & elem_dims =
585  old_context.elem_dimensions();
586 
587  // Loop over variables and dimensions, to prerequest
588  for (const auto & dim : elem_dims)
589  {
590  FEAbstract * fe = nullptr;
591  if (old_context.active_vars())
592  for (auto var : *old_context.active_vars())
593  {
594  old_context.get_element_fe(var, fe, dim);
595  get_shape_outputs(*fe);
596  }
597  else
598  for (auto var : make_range(sys.n_vars()))
599  {
600  old_context.get_element_fe(var, fe, dim);
601  get_shape_outputs(*fe);
602  }
603  }
604  }
605 
606  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  bool check_old_context (const FEMContext & c, const Point & p)
643  {
644  LOG_SCOPE ("check_old_context(c,p)", "OldSolutionBase");
645  const Elem & elem = c.get_elem();
646  if (last_elem != &elem)
647  {
648  if (elem.refinement_flag() == Elem::JUST_REFINED)
649  {
650  old_context.pre_fe_reinit(sys, elem.parent());
651  }
652  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  const Real master_tol = out_of_elem_tol / elem.hmax() * 2;
662 
663  for (auto & child : elem.child_ref_range())
664  if (child.close_to_point(p, master_tol))
665  {
666  old_context.pre_fe_reinit(sys, &child);
667  break;
668  }
669 
671  (old_context.get_elem().close_to_point(p, master_tol));
672  }
673  else
674  {
675  if (!elem.get_old_dof_object())
676  return false;
677 
678  old_context.pre_fe_reinit(sys, &elem);
679  }
680 
681  last_elem = &elem;
682  }
683  else
684  {
685  libmesh_assert(old_context.has_elem());
686 
687  const Real master_tol = out_of_elem_tol / elem.hmax() * 2;
688 
689  if (!old_context.get_elem().close_to_point(p, master_tol))
690  {
691  libmesh_assert_equal_to
693 
694  for (auto & child : elem.child_ref_range())
695  if (child.close_to_point(p, master_tol))
696  {
697  old_context.pre_fe_reinit(sys, &child);
698  break;
699  }
700 
702  (old_context.get_elem().close_to_point(p, master_tol));
703  }
704  }
705 
706  return true;
707  }
708 
709 protected:
710  const Elem * last_elem;
711  const System & sys;
713  std::vector<unsigned int> component_to_var;
714 
715  static const Real out_of_elem_tol;
716 };
717 
718 
727 template <typename Output,
728  void (FEMContext::*point_output) (unsigned int,
729  const Point &,
730  Output &,
731  const Real) const>
733 {
734 public:
737  typedef Output FunctorValue;
739 
741  const NumericVector<Number> & old_sol,
742  const std::vector<unsigned int> * vars) :
743  OldSolutionBase<Output, point_output>(sys_in, vars),
744  old_solution(old_sol)
745  {
746  this->old_context.set_algebraic_type(FEMContext::OLD);
747  this->old_context.set_custom_solution(&old_solution);
748  }
749 
751  OldSolutionBase<Output, point_output>(in.sys, in.old_context.active_vars()),
752  old_solution(in.old_solution)
753  {
754  this->old_context.set_algebraic_type(FEMContext::OLD);
755  this->old_context.set_custom_solution(&old_solution);
756  }
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  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  LOG_SCOPE ("eval_at_point()", "OldSolutionValue");
773 
774  if (!skip_context_check)
775  if (!this->check_old_context(c, p))
776  return Output(0);
777 
778  // Handle offset from non-scalar components in previous variables
779  libmesh_assert_less(i, this->component_to_var.size());
780  unsigned int var = this->component_to_var[i];
781 
782  Output n;
783  (this->old_context.*point_output)(var, p, n, this->out_of_elem_tol);
784  return n;
785  }
786 
787  template <typename T = Output,
788  typename std::enable_if<std::is_same<T, Number>::value, int>::type = 0>
789  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  LOG_SCOPE ("eval_mixed_derivatives", "OldSolutionValue");
796 
797  // This should only be called on vertices
798  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  libmesh_assert_less(i, this->component_to_var.size());
803  unsigned int var = this->component_to_var[i];
804 
805  // We have 1 mixed derivative in 2D, 4 in 3D
806  const unsigned int n_mixed = (dim-1) * (dim-1);
807  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  const DofObject * old_dof_object = n.get_old_dof_object();
812  if (old_dof_object &&
813  old_dof_object->n_vars(this->sys.number()) &&
814  old_dof_object->n_comp(this->sys.number(), var))
815  {
816  const dof_id_type first_old_id =
817  old_dof_object->dof_number(this->sys.number(), var, dim);
818  std::vector<dof_id_type> old_ids(n_mixed);
819  std::iota(old_ids.begin(), old_ids.end(), first_old_id);
820  old_solution.get(old_ids, derivs);
821  }
822  else
823  {
824  std::fill(derivs.begin(), derivs.end(), 0);
825  }
826  }
827 
828  template <typename T = Output,
829  typename std::enable_if<!std::is_same<T, Number>::value, int>::type = 0>
831  const FEMContext &, unsigned int, unsigned int, const Node &, std::vector<Output> &)
832  {
833  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  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  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  indices.clear();
848 
849  this->sys.get_dof_map().dof_indices(elem, node_num, indices, var_num);
850 
851  std::vector<dof_id_type> old_indices;
852 
853  this->sys.get_dof_map().old_dof_indices(elem, node_num, old_indices, var_num);
854 
855  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  bool invalid_old_index = false;
860  for (const auto & di : old_indices)
861  if (di == DofObject::invalid_id)
862  invalid_old_index = true;
863 
864  values.resize(old_indices.size());
865  if (invalid_old_index)
866  {
867  for (auto i : index_range(old_indices))
868  {
869  const dof_id_type di = old_indices[i];
870  if (di == DofObject::invalid_id)
871  values[i] = 0;
872  else
873  values[i] = old_solution(di);
874  }
875  }
876  else
877  old_solution.get(old_indices, values);
878  }
879 
880 
881  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  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  const Elem & old_elem =
893  (elem.refinement_flag() == Elem::JUST_REFINED) ?
894  *elem.parent() : elem;
895 
896  // If there are any element-based DOF numbers, get them
897  const unsigned int nc = FEInterface::n_dofs_per_elem(fe_type, &elem);
898 
899  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  if (nc != 0)
907  {
908  const DofObject & old_dof_object = old_elem.get_old_dof_object_ref();
909  libmesh_assert_greater(elem.n_systems(), sys_num);
910 
911  const std::pair<unsigned int, unsigned int>
912  vg_and_offset = elem.var_to_vg_and_offset(sys_num,var_num);
913  const unsigned int vg = vg_and_offset.first;
914  const unsigned int vig = vg_and_offset.second;
915 
916  unsigned int n_comp = old_dof_object.n_comp_group(sys_num,vg);
917  n_comp = std::min(n_comp, nc);
918 
919  std::vector<dof_id_type> old_dof_indices(n_comp);
920 
921  for (unsigned int i=0; i != n_comp; ++i)
922  {
923  const dof_id_type d_old =
924  old_dof_object.dof_number(sys_num, vg, vig, i, n_comp);
925  const dof_id_type d_new =
926  elem.dof_number(sys_num, vg, vig, i, n_comp);
927  libmesh_assert_not_equal_to (d_old, DofObject::invalid_id);
928  libmesh_assert_not_equal_to (d_new, DofObject::invalid_id);
929 
930  old_dof_indices[i] = d_old;
931  indices[i] = d_new;
932  }
933 
934  values.resize(n_comp);
935  old_solution.get(old_dof_indices, values);
936 
937  for (unsigned int i=n_comp; i != nc; ++i)
938  {
939  const dof_id_type d_new =
940  elem.dof_number(sys_num, vg, vig, i, n_comp);
941  libmesh_assert_not_equal_to (d_new, DofObject::invalid_id);
942  indices[i] = d_new;
943  }
944 
945  values.resize(nc, 0);
946  }
947  else
948  values.clear();
949  }
950 
951 private:
953 };
954 
955 template<>
956 inline void
958 {
959  fe.request_phi();
960 }
961 
962 
963 template<>
964 inline void
966 {
967  fe.request_dphi();
968 }
969 
970 template<>
971 inline void
973 {
974  fe.request_phi();
975 }
976 
977 
978 template<>
979 inline void
981 {
982  fe.request_dphi();
983 }
984 
985 
986 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
987 template<>
988 inline void
990 {
991  fe.request_phi();
992 }
993 
994 
995 template<>
996 inline void
998 {
999  fe.request_dphi();
1000 }
1001 
1002 template<>
1003 inline void
1005 {
1006  fe.request_phi();
1007 }
1008 
1009 
1010 template<>
1011 inline void
1013 {
1014  fe.request_dphi();
1015 }
1016 #endif // LIBMESH_USE_COMPLEX_NUMBERS
1017 
1018 
1019 template<>
1020 inline
1021 Number
1024  unsigned int i,
1025  unsigned int /* elem_dim */,
1026  const Node & n,
1027  bool extra_hanging_dofs,
1028  Real /* time */)
1029 {
1030  LOG_SCOPE ("Number eval_at_node()", "OldSolutionValue");
1031 
1032  // This should only be called on vertices
1033  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  libmesh_assert_less(i, this->component_to_var.size());
1038  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  const Elem::RefinementState flag = c.get_elem().refinement_flag();
1051 
1052  const DofObject * old_dof_object = n.get_old_dof_object();
1053  if (old_dof_object &&
1054  (!extra_hanging_dofs ||
1055  flag == Elem::JUST_COARSENED ||
1056  flag == Elem::DO_NOTHING) &&
1057  old_dof_object->n_vars(sys.number()) &&
1058  old_dof_object->n_comp(sys.number(), var))
1059  {
1060  const dof_id_type old_id =
1061  old_dof_object->dof_number(sys.number(), var, 0);
1062  return old_solution(old_id);
1063  }
1064 
1065  return this->eval_at_point(c, i, n, 0, false);
1066 }
1067 
1068 template <>
1069 inline
1070 Gradient
1073  unsigned int i,
1074  unsigned int /* elem_dim */,
1075  const Node & n,
1076  bool extra_hanging_dofs,
1077  Real /* time */)
1078 {
1079  LOG_SCOPE ("Number eval_at_node()", "OldSolutionValue");
1080 
1081  // This should only be called on vertices
1082  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  libmesh_assert_less(i, this->component_to_var.size());
1087  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  const auto & elem = c.get_elem();
1100 
1101  const Elem::RefinementState flag = elem.refinement_flag();
1102 
1103  const DofObject * old_dof_object = n.get_old_dof_object();
1104  if (old_dof_object &&
1105  (!extra_hanging_dofs ||
1106  flag == Elem::JUST_COARSENED ||
1107  flag == Elem::DO_NOTHING) &&
1108  old_dof_object->n_vars(sys.number()) &&
1109  old_dof_object->n_comp(sys.number(), var))
1110  {
1111  Gradient return_val;
1112 
1113  for (auto dim : make_range(elem.dim()))
1114  {
1115  const dof_id_type old_id =
1116  old_dof_object->dof_number(sys.number(), var, dim);
1117  return_val(dim) = old_solution(old_id);
1118  }
1119 
1120  return return_val;
1121  }
1122 
1123  return this->eval_at_point(c, i, n, 0, false);
1124 }
1125 
1126 
1127 
1128 template<>
1129 inline
1130 Gradient
1133  unsigned int i,
1134  unsigned int elem_dim,
1135  const Node & n,
1136  bool extra_hanging_dofs,
1137  Real /* time */)
1138 {
1139  LOG_SCOPE ("Gradient eval_at_node()", "OldSolutionValue");
1140 
1141  // This should only be called on vertices
1142  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  libmesh_assert_less(i, this->component_to_var.size());
1147  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  const Elem::RefinementState flag = c.get_elem().refinement_flag();
1160 
1161  const DofObject * old_dof_object = n.get_old_dof_object();
1162  if (old_dof_object &&
1163  (!extra_hanging_dofs ||
1164  flag == Elem::JUST_COARSENED ||
1165  flag == Elem::DO_NOTHING) &&
1166  old_dof_object->n_vars(sys.number()) &&
1167  old_dof_object->n_comp(sys.number(), var))
1168  {
1169  Gradient g;
1170  for (unsigned int d = 0; d != elem_dim; ++d)
1171  {
1172  const dof_id_type old_id =
1173  old_dof_object->dof_number(sys.number(), var, d+1);
1174  g(d) = old_solution(old_id);
1175  }
1176  return g;
1177  }
1178 
1179  return this->eval_at_point(c, i, n, 0, false);
1180 }
1181 
1182 template<>
1183 inline
1184 Tensor
1187  unsigned int,
1188  unsigned int,
1189  const Node &,
1190  bool,
1191  Real)
1192 {
1193  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>
1204 
1205 #endif // LIBMESH_ENABLE_AMR
1206 
1211 template <typename FFunctor, typename GFunctor,
1212  typename FValue, typename ProjectionAction>
1214  (const ConstElemRange & range)
1215 {
1216  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  done_saving_ids = false;
1221 
1222  SortAndCopy sort_work(*this);
1223  Threads::parallel_reduce (range, sort_work);
1224  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  ids_to_push;
1229 
1230  ids_to_push.insert(sort_work.new_ids_to_push.begin(),
1231  sort_work.new_ids_to_push.end());
1232  ids_to_save.insert(sort_work.new_ids_to_save.begin(),
1233  sort_work.new_ids_to_save.end());
1234 
1235  std::vector<node_projection> vertices(sort_work.vertices.begin(),
1236  sort_work.vertices.end());
1237 
1238  done_saving_ids = sort_work.edges.empty() &&
1239  sort_work.sides.empty() && sort_work.interiors.empty();
1240 
1241  {
1242  ProjectVertices project_vertices(*this);
1243  Threads::parallel_reduce (node_range(&vertices), project_vertices);
1244  libmesh_merge_move(ids_to_push, project_vertices.new_ids_to_push);
1245  libmesh_merge_move(ids_to_save, project_vertices.new_ids_to_save);
1246  }
1247 
1248  done_saving_ids = sort_work.sides.empty() && sort_work.interiors.empty();
1249 
1250  this->send_and_insert_dof_values(ids_to_push, action);
1251 
1252  {
1253  std::vector<node_projection> edges(sort_work.edges.begin(), sort_work.edges.end());
1254  ProjectEdges project_edges(*this);
1255  Threads::parallel_reduce (node_range(&edges), project_edges);
1256  libmesh_merge_move(ids_to_push, project_edges.new_ids_to_push);
1257  libmesh_merge_move(ids_to_save, project_edges.new_ids_to_save);
1258  }
1259 
1260  done_saving_ids = sort_work.interiors.empty();
1261 
1262  this->send_and_insert_dof_values(ids_to_push, action);
1263 
1264  {
1265  std::vector<node_projection> sides(sort_work.sides.begin(), sort_work.sides.end());
1266  ProjectSides project_sides(*this);
1267  Threads::parallel_reduce (node_range(&sides), project_sides);
1268  libmesh_merge_move(ids_to_push, project_sides.new_ids_to_push);
1269  libmesh_merge_move(ids_to_save, project_sides.new_ids_to_save);
1270  }
1271 
1272  done_saving_ids = true;
1273  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  ProjectInteriors project_interiors(*this);
1279  project_interiors);
1280 }
1281 
1282 
1283 template <typename FFunctor, typename GFunctor,
1284  typename FValue, typename ProjectionAction>
1287  projector(p),
1288  action(p.master_action),
1289  f(p.master_f),
1290  context(p.system, &p.variables, /* allocate_local_matrices= */ false),
1291  conts(p.system.n_vars()),
1292  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  for (const auto & var : this->projector.variables)
1297  {
1298  // FIXME: Need to generalize this to vector-valued elements. [PB]
1299  FEAbstract * fe = nullptr;
1300  FEAbstract * side_fe = nullptr;
1301  FEAbstract * edge_fe = nullptr;
1302 
1303  const std::set<unsigned char> & elem_dims =
1304  context.elem_dimensions();
1305 
1306  for (const auto & dim : elem_dims)
1307  {
1308  context.get_element_fe( var, fe, dim );
1309  if (fe->get_fe_type().family == SCALAR)
1310  continue;
1311  if (dim > 1)
1312  context.get_side_fe( var, side_fe, dim );
1313  if (dim > 2)
1314  context.get_edge_fe( var, edge_fe );
1315 
1316  fe->get_JxW();
1317  fe->get_xyz();
1318  fe->get_JxW();
1319 
1320  fe->request_phi();
1321  if (dim > 1)
1322  {
1323  side_fe->get_JxW();
1324  side_fe->get_xyz();
1325  side_fe->request_phi();
1326  }
1327  if (dim > 2)
1328  {
1329  edge_fe->get_JxW();
1330  edge_fe->get_xyz();
1331  edge_fe->request_phi();
1332  }
1333 
1334  const FEContinuity cont = fe->get_continuity();
1335  this->conts[var] = cont;
1336  if (cont == C_ONE)
1337  {
1338  fe->request_dphi();
1339  if (dim > 1)
1340  side_fe->request_dphi();
1341  if (dim > 2)
1342  edge_fe->request_dphi();
1343  }
1344 
1345  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  f.init_context(context);
1351 }
1352 
1353 
1354 template <typename FFunctor, typename GFunctor,
1355  typename FValue, typename ProjectionAction>
1358  SubFunctor(p)
1359 {
1360 
1361  // Our C1 elements need gradient information
1362  for (const auto & var : this->projector.variables)
1363  if (this->conts[var] == C_ONE)
1364  {
1366  g = std::make_unique<GFunctor>(*p.master_g);
1367  g->init_context(context);
1368  return;
1369  }
1370 }
1371 
1372 template <typename FFunctor, typename GFunctor, typename FValue, typename ProjectionAction>
1373 template <typename InsertInput,
1374  typename std::enable_if<
1376  int>::type>
1377 void
1379  dof_id_type, const InsertInput & , processor_id_type)
1380 {
1381  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<
1389  int>::type>
1390 void
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  if (id == DofObject::invalid_id)
1397  return;
1398 
1399  auto iter = new_ids_to_push.find(id);
1400  if (iter == new_ids_to_push.end())
1401  action.insert(id, val);
1402  else
1403  {
1405  iter->second = std::make_pair(val, pid);
1406  }
1407  if (!this->projector.done_saving_ids)
1408  {
1409  libmesh_assert(!new_ids_to_save.count(id));
1410  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<
1418  int>::type>
1419 void
1421  const std::vector<dof_id_type> &,
1422  const std::vector<InsertInput> &,
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<
1433  int>::type>
1434 void
1436  const std::vector<dof_id_type> & ids,
1437  const std::vector<InsertInput> & vals,
1438  processor_id_type pid)
1439 {
1440  libmesh_assert_equal_to(ids.size(), vals.size());
1441  for (auto i : index_range(ids))
1442  {
1443  const dof_id_type id = ids[i];
1444 
1445  // We may see invalid ids when expanding a subdomain with a
1446  // restricted variable
1447  if (id == DofObject::invalid_id)
1448  continue;
1449 
1450  const InsertInput & val = vals[i];
1451 
1452  auto iter = new_ids_to_push.find(id);
1453  if (iter == new_ids_to_push.end())
1454  action.insert(id, val);
1455  else
1456  {
1458  iter->second = std::make_pair(val, pid);
1459  }
1460  if (!this->projector.done_saving_ids)
1461  {
1462  libmesh_assert(!new_ids_to_save.count(id));
1463  new_ids_to_save[id] = val;
1464  }
1465  }
1466 }
1467 
1468 template <typename FFunctor, typename GFunctor,
1469  typename FValue, typename ProjectionAction>
1472 {
1473  new_ids_to_push.insert(other.new_ids_to_push.begin(), other.new_ids_to_push.end());
1474  new_ids_to_save.insert(other.new_ids_to_save.begin(), other.new_ids_to_save.end());
1475 }
1476 
1477 
1478 template <typename FFunctor, typename GFunctor,
1479  typename FValue, typename ProjectionAction>
1481  (const ConstElemRange & range)
1482 {
1483  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  std::unordered_map<const Node *, std::pair<var_set, var_set>> copied_nodes;
1510 
1511  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  std::vector<unsigned short> extra_hanging_dofs;
1517  bool all_extra_hanging_dofs = true;
1518  for (auto v_num : this->projector.variables)
1519  {
1520  if (extra_hanging_dofs.size() <= v_num)
1521  extra_hanging_dofs.resize(v_num+1, false);
1522  extra_hanging_dofs[v_num] =
1524 
1525  if (!extra_hanging_dofs[v_num])
1526  all_extra_hanging_dofs = false;
1527  }
1528 
1529  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  bool copy_this_elem = false;
1534 
1535 #ifdef LIBMESH_ENABLE_AMR
1536  // If we're projecting from an old grid
1537  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  const DofObject * old_dof_object = elem->get_old_dof_object();
1544  const Elem::RefinementState h_flag = elem->refinement_flag();
1545  const Elem::RefinementState p_flag = elem->p_refinement_flag();
1546  if (!old_dof_object &&
1547  h_flag != Elem::JUST_REFINED &&
1548  h_flag != Elem::JUST_COARSENED)
1549  continue;
1550 
1551  // If this is an unchanged element, just copy everything
1552  if (h_flag != Elem::JUST_REFINED &&
1553  h_flag != Elem::JUST_COARSENED &&
1554  p_flag != Elem::JUST_REFINED &&
1555  p_flag != Elem::JUST_COARSENED)
1556  copy_this_elem = true;
1557  else
1558  {
1559  bool reinitted = false;
1560 
1561  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  const bool copy_possible =
1566  p_level == 0 &&
1567  h_flag != Elem::JUST_COARSENED &&
1568  p_flag != Elem::JUST_COARSENED;
1569 
1570  std::vector<typename FFunctor::ValuePushType> Ue(1);
1571  std::vector<dof_id_type> elem_dof_ids(1);
1572 
1573  for (auto v_num : this->projector.variables)
1574  {
1575  const Variable & var = system.variable(v_num);
1576  if (!var.active_on_subdomain(elem->subdomain_id()))
1577  continue;
1578  const FEType fe_type = var.type();
1579 
1580  if (fe_type.family == MONOMIAL &&
1581  fe_type.order == CONSTANT &&
1582  copy_possible)
1583  {
1584  if (!reinitted)
1585  {
1586  reinitted = true;
1587  context.pre_fe_reinit(system, elem);
1588  }
1589 
1590  f.eval_old_dofs(*elem, fe_type, sys_num, v_num,
1591  elem_dof_ids, Ue);
1592 
1593  action.insert(elem_dof_ids[0], Ue[0]);
1594  }
1595  }
1596  }
1597  }
1598 #endif // LIBMESH_ENABLE_AMR
1599 
1600  const int dim = elem->dim();
1601 
1602  const unsigned int n_vertices = elem->n_vertices();
1603  const unsigned int n_edges = elem->n_edges();
1604  const unsigned int n_nodes = elem->n_nodes();
1605 
1606  // In 1-D we already handle our sides as vertices
1607  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  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  const bool has_edge_nodes = (n_nodes > n_vertices && dim > 2);
1615 
1616  // If we have even more nodes, the next is a side node.
1617  const bool has_side_nodes =
1618  (n_nodes > n_vertices + ((dim > 2) * n_edges));
1619 
1620  // We may be out of nodes at this point or we have interior
1621  // nodes which may have DoFs to project too
1622  const bool has_interior_nodes =
1623  (n_nodes > n_vertices + ((dim > 2) * n_edges) + n_sides);
1624 
1625  for (auto v_num : this->projector.variables)
1626  {
1627  const Variable & var = this->projector.system.variable(v_num);
1628  if (!var.active_on_subdomain(elem->subdomain_id()))
1629  continue;
1630  const FEType fe_type = var.type();
1631  const auto & dof_map = this->projector.system.get_dof_map();
1632  const auto vg = dof_map.var_group_from_var_number(v_num);
1633  const bool add_p_level = dof_map.should_p_refine(vg);
1634 
1635  // If we're trying to do projections on an isogeometric
1636  // analysis mesh, only the finite element nodes constrained
1637  // by those spline nodes truly have delta function degrees
1638  // of freedom. We'll have to back out the spline degrees of
1639  // freedom indirectly later.
1640  if (fe_type.family == RATIONAL_BERNSTEIN &&
1641  elem->type() == NODEELEM)
1642  continue;
1643 
1644  if (FEInterface::n_dofs_at_node(fe_type, elem, 0, add_p_level))
1645  vertex_vars.insert(vertex_vars.end(), v_num);
1646 
1647  // The first non-vertex node is always an edge node if those
1648  // exist. All edge nodes have the same number of DoFs
1649  if (has_edge_nodes)
1650  if (FEInterface::n_dofs_at_node(fe_type, elem, n_vertices, add_p_level))
1651  edge_vars.insert(edge_vars.end(), v_num);
1652 
1653  if (has_side_nodes)
1654  {
1655  if (dim != 3)
1656  {
1657  if (FEInterface::n_dofs_at_node(fe_type, elem, n_vertices, add_p_level))
1658  side_vars.insert(side_vars.end(), v_num);
1659  }
1660  else
1661  // In 3D, not all face nodes always have the same number of
1662  // DoFs! We'll loop over all sides to be safe.
1663  for (unsigned int n = 0; n != n_nodes; ++n)
1664  if (elem->is_face(n))
1665  if (FEInterface::n_dofs_at_node(fe_type, elem, n, add_p_level))
1666  {
1667  side_vars.insert(side_vars.end(), v_num);
1668  break;
1669  }
1670  }
1671 
1672  if (FEInterface::n_dofs_per_elem(fe_type, elem, add_p_level) ||
1673  (has_interior_nodes &&
1674  FEInterface::n_dofs_at_node(fe_type, elem, n_nodes-1, add_p_level)))
1675  {
1676 #ifdef LIBMESH_ENABLE_AMR
1677  // We may have already just copied constant monomials,
1678  // or we may be about to copy the whole element
1679  if ((f.is_grid_projection() &&
1680  fe_type.family == MONOMIAL &&
1681  fe_type.order == CONSTANT &&
1682  elem->p_level() == 0 &&
1683  elem->refinement_flag() != Elem::JUST_COARSENED &&
1684  elem->p_refinement_flag() != Elem::JUST_COARSENED)
1685  || copy_this_elem
1686  )
1687  continue;
1688 #endif // LIBMESH_ENABLE_AMR
1689 
1690  // We need to project any other variables
1691  if (interiors.empty() || interiors.back() != elem)
1692  interiors.push_back(elem);
1693  }
1694  }
1695 
1696  // We'll use a greedy algorithm in most cases: if another
1697  // element has already claimed some of our DoFs, we'll let it do
1698  // the work.
1699  auto erase_covered_vars = []
1700  (const Node * node, var_set & remaining, const ves_multimap & covered)
1701  {
1702  auto covered_range = covered.equal_range(node);
1703  for (const auto & v_ent : as_range(covered_range))
1704  for (const unsigned int var_covered :
1705  std::get<2>(v_ent.second))
1706  remaining.erase(var_covered);
1707  };
1708 
1709  auto erase_nonhanging_vars = [&extra_hanging_dofs]
1710  (const Node * node, var_set & remaining, const ves_multimap & covered)
1711  {
1712  auto covered_range = covered.equal_range(node);
1713  for (const auto & v_ent : as_range(covered_range))
1714  for (const unsigned int var_covered :
1715  std::get<2>(v_ent.second))
1716  if (!extra_hanging_dofs[var_covered])
1717  remaining.erase(var_covered);
1718  };
1719 
1720  auto erase_copied_vars = [&copied_nodes, &extra_hanging_dofs]
1721  (const Node * node, bool is_vertex, var_set & remaining)
1722  {
1723  auto copying_range = copied_nodes.equal_range(node);
1724  for (const auto & v_ent : as_range(copying_range))
1725  {
1726  for (const unsigned int var_covered :
1727  v_ent.second.first)
1728  if (is_vertex || !extra_hanging_dofs[var_covered])
1729  remaining.erase(var_covered);
1730 
1731  for (const unsigned int var_covered :
1732  v_ent.second.second)
1733  if (!is_vertex || !extra_hanging_dofs[var_covered])
1734  remaining.erase(var_covered);
1735  }
1736  };
1737 
1738  for (unsigned int v=0; v != n_vertices; ++v)
1739  {
1740  const Node * node = elem->node_ptr(v);
1741 
1742  auto remaining_vars = vertex_vars;
1743 
1744  erase_covered_vars(node, remaining_vars, vertices);
1745 
1746  if (remaining_vars.empty())
1747  continue;
1748 
1749  if (!all_extra_hanging_dofs)
1750  {
1751  erase_nonhanging_vars(node, remaining_vars, edges);
1752  if (remaining_vars.empty())
1753  continue;
1754 
1755  erase_nonhanging_vars(node, remaining_vars, sides);
1756  if (remaining_vars.empty())
1757  continue;
1758  }
1759 
1760  erase_copied_vars(node, true, remaining_vars);
1761  if (remaining_vars.empty())
1762  continue;
1763 
1764  if (copy_this_elem)
1765  {
1766  std::vector<dof_id_type> node_dof_ids;
1767  std::vector<typename FFunctor::ValuePushType> values;
1768 
1769  for (auto var : remaining_vars)
1770  {
1771  f.eval_old_dofs(*elem, v, var, node_dof_ids, values);
1772  insert_ids(node_dof_ids, values, node->processor_id());
1773  }
1774  copied_nodes[node].first.insert(remaining_vars.begin(),
1775  remaining_vars.end());
1776  this->find_dofs_to_send(*node, *elem, v, remaining_vars);
1777  }
1778  else
1779  vertices.emplace
1780  (node, std::make_tuple(elem, v, std::move(remaining_vars)));
1781  }
1782 
1783  if (has_edge_nodes)
1784  {
1785  for (unsigned int e=0; e != n_edges; ++e)
1786  {
1787  const Node * node = elem->node_ptr(n_vertices+e);
1788 
1789  auto remaining_vars = edge_vars;
1790 
1791  erase_covered_vars(node, remaining_vars, edges);
1792  if (remaining_vars.empty())
1793  continue;
1794 
1795  erase_covered_vars(node, remaining_vars, sides);
1796  if (remaining_vars.empty())
1797  continue;
1798 
1799  if (!all_extra_hanging_dofs)
1800  {
1801  erase_nonhanging_vars(node, remaining_vars, vertices);
1802  if (remaining_vars.empty())
1803  continue;
1804  }
1805 
1806  erase_copied_vars(node, false, remaining_vars);
1807  if (remaining_vars.empty())
1808  continue;
1809 
1810  if (copy_this_elem)
1811  {
1812  std::vector<dof_id_type> edge_dof_ids;
1813  std::vector<typename FFunctor::ValuePushType> values;
1814 
1815  for (auto var : remaining_vars)
1816  {
1817  f.eval_old_dofs(*elem, n_vertices+e, var, edge_dof_ids, values);
1818  insert_ids(edge_dof_ids, values, node->processor_id());
1819  }
1820 
1821  copied_nodes[node].second.insert(remaining_vars.begin(),
1822  remaining_vars.end());
1823  this->find_dofs_to_send(*node, *elem, n_vertices+e, remaining_vars);
1824  }
1825  else
1826  edges.emplace
1827  (node, std::make_tuple(elem, e, std::move(remaining_vars)));
1828  }
1829  }
1830 
1831  if (has_side_nodes)
1832  {
1833  for (unsigned int side=0; side != n_sides; ++side)
1834  {
1835  const Node * node = nullptr;
1836  unsigned short node_num = n_vertices+(dim>2)*n_edges+side;
1837  if (dim != 3)
1838  node = elem->node_ptr(node_num);
1839  else
1840  {
1841  // In 3D only some sides may have nodes
1842  for (unsigned int n = 0; n != n_nodes; ++n)
1843  {
1844  if (!elem->is_face(n))
1845  continue;
1846 
1847  if (elem->is_node_on_side(n, side))
1848  {
1849  node_num = n;
1850  node = elem->node_ptr(node_num);
1851  break;
1852  }
1853  }
1854  }
1855 
1856  if (!node)
1857  continue;
1858 
1859  auto remaining_vars = side_vars;
1860 
1861  erase_covered_vars(node, remaining_vars, edges);
1862  if (remaining_vars.empty())
1863  continue;
1864 
1865  erase_covered_vars(node, remaining_vars, sides);
1866  if (remaining_vars.empty())
1867  continue;
1868 
1869  if (!all_extra_hanging_dofs)
1870  {
1871  erase_nonhanging_vars(node, remaining_vars, vertices);
1872  if (remaining_vars.empty())
1873  continue;
1874  }
1875 
1876  erase_copied_vars(node, false, remaining_vars);
1877  if (remaining_vars.empty())
1878  continue;
1879 
1880  if (copy_this_elem)
1881  {
1882  std::vector<dof_id_type> side_dof_ids;
1883  std::vector<typename FFunctor::ValuePushType> values;
1884 
1885  for (auto var : remaining_vars)
1886  {
1887  f.eval_old_dofs(*elem, node_num, var, side_dof_ids, values);
1888  insert_ids(side_dof_ids, values, node->processor_id());
1889  }
1890 
1891  copied_nodes[node].second.insert(remaining_vars.begin(),
1892  remaining_vars.end());
1893  this->find_dofs_to_send(*node, *elem, node_num, remaining_vars);
1894  }
1895  else
1896  sides.emplace
1897  (node, std::make_tuple(elem, side, std::move(remaining_vars)));
1898  }
1899  }
1900 
1901  // Elements with elemental dofs might need those copied too.
1902  if (copy_this_elem)
1903  {
1904  std::vector<typename FFunctor::ValuePushType> U;
1905  std::vector<dof_id_type> dof_ids;
1906 
1907  for (auto v_num : this->projector.variables)
1908  {
1909  const Variable & var = system.variable(v_num);
1910  if (!var.active_on_subdomain(elem->subdomain_id()))
1911  continue;
1912  FEType fe_type = var.type();
1913 
1914  f.eval_old_dofs(*elem, fe_type, sys_num, v_num,
1915  dof_ids, U);
1916  action.insert(dof_ids, U);
1917 
1918  if (has_interior_nodes)
1919  {
1920  f.eval_old_dofs(*elem, n_nodes-1, v_num, dof_ids, U);
1921  action.insert(dof_ids, U);
1922  }
1923  }
1924  }
1925  }
1926 }
1927 
1928 
1929 template <typename FFunctor, typename GFunctor,
1930  typename FValue, typename ProjectionAction>
1933 {
1934  auto merge_multimaps = [](ves_multimap & self, const ves_multimap & other_mm)
1935  {
1936  for (const auto & pair : other_mm)
1937  {
1938  const Node * node = pair.first;
1939  auto remaining_vars = std::get<2>(pair.second);
1940 
1941  auto my_range = self.equal_range(node);
1942  for (const auto & v_ent : as_range(my_range))
1943  for (const unsigned int var_covered :
1944  std::get<2>(v_ent.second))
1945  remaining_vars.erase(var_covered);
1946 
1947  if (!remaining_vars.empty())
1948  self.emplace
1949  (node, std::make_tuple(std::get<0>(pair.second),
1950  std::get<1>(pair.second),
1951  std::move(remaining_vars)));
1952  }
1953  };
1954 
1955  merge_multimaps(vertices, other.vertices);
1956  merge_multimaps(edges, other.edges);
1957  merge_multimaps(sides, other.sides);
1958 
1959  std::sort(interiors.begin(), interiors.end());
1960  std::vector<const Elem *> other_interiors = other.interiors;
1961  std::sort(other_interiors.begin(), other_interiors.end());
1962  std::vector<const Elem *> merged_interiors;
1963  std::set_union(interiors.begin(), interiors.end(),
1964  other_interiors.begin(), other_interiors.end(),
1965  std::back_inserter(merged_interiors));
1966  interiors.swap(merged_interiors);
1967 
1968  SubFunctor::join(other);
1969 }
1970 
1971 namespace
1972 {
1973 template <typename Output, typename Input>
1974 typename std::enable_if<ScalarTraits<Input>::value, Output>::type
1975 raw_value(const Input & input, unsigned int /*index*/)
1976 {
1977  return input;
1978 }
1979 
1980 template <typename Output, template <typename> class WrapperClass, typename T>
1983  Output>::type
1984 raw_value(const WrapperClass<T> & input, unsigned int index)
1985 {
1986  return input.slice(index);
1987 }
1988 
1989 template <typename T>
1990 typename T::value_type
1991 grad_component(const T &, unsigned int);
1992 
1993 template <typename T>
1995 grad_component(const VectorValue<T> & grad, unsigned int component)
1996 {
1997  return grad(component);
1998 }
1999 
2000 template <typename T>
2002 grad_component(const TensorValue<T> & grad, unsigned int component)
2003 {
2004  libmesh_error_msg("This function should only be called for Hermites. "
2005  "I don't know how you got here");
2006  return grad(component, component);
2007 }
2008 
2009 
2010 }
2011 
2012 template <typename FFunctor, typename GFunctor,
2013  typename FValue, typename ProjectionAction>
2014 void GenericProjector<FFunctor, GFunctor, FValue, ProjectionAction>::ProjectVertices::operator()
2015  (const node_range & range)
2016 {
2017  LOG_SCOPE ("project_vertices","GenericProjector");
2018 
2019  const unsigned int sys_num = system.number();
2020 
2021  // Variables with extra hanging dofs can't safely use eval_at_node
2022  // in as many places as variables without can.
2023  std::vector<unsigned short> extra_hanging_dofs;
2024  for (auto v_num : this->projector.variables)
2025  {
2026  if (extra_hanging_dofs.size() <= v_num)
2027  extra_hanging_dofs.resize(v_num+1, false);
2028  extra_hanging_dofs[v_num] =
2030  }
2031 
2032  for (const auto & v_pair : range)
2033  {
2034  const Node & vertex = *v_pair.first;
2035  const Elem & elem = *std::get<0>(v_pair.second);
2036  const unsigned int n = std::get<1>(v_pair.second);
2037  const var_set & vertex_vars = std::get<2>(v_pair.second);
2038 
2039  context.pre_fe_reinit(system, &elem);
2040 
2041  this->find_dofs_to_send(vertex, elem, n, vertex_vars);
2042 
2043  // Look at all the variables we're supposed to interpolate from
2044  // this element on this vertex
2045  for (const auto & var : vertex_vars)
2046  {
2047  const Variable & variable = system.variable(var);
2048  const FEType & base_fe_type = variable.type();
2049  const unsigned int var_component =
2051 
2052  if (base_fe_type.family == SCALAR)
2053  continue;
2054 
2055  const FEContinuity cont = this->conts[var];
2056  const FEFieldType field_type = this->field_types[var];
2057 
2058  if (cont == DISCONTINUOUS)
2059  {
2060  libmesh_assert_equal_to(vertex.n_comp(sys_num, var), 0);
2061  }
2062  else if (cont == C_ZERO ||
2063  cont == SIDE_DISCONTINUOUS)
2064  {
2065  if (cont == SIDE_DISCONTINUOUS &&
2066  elem.dim() != 1)
2067  {
2068  libmesh_assert_equal_to(vertex.n_comp(sys_num, var), 0);
2069  continue;
2070  }
2071 
2072  const FValue val = f.eval_at_node
2073  (context, var_component, /*dim=*/ 0, // Don't care w/C0
2074  vertex, extra_hanging_dofs[var], system.time);
2075 
2076  if (field_type == TYPE_VECTOR)
2077  {
2078  libmesh_assert_equal_to(vertex.n_comp(sys_num, var), elem.dim());
2079 
2080  // We will have a number of nodal value DoFs equal to the elem dim
2081  for (auto i : make_range(elem.dim()))
2082  {
2083  const dof_id_type id = vertex.dof_number(sys_num, var, i);
2084 
2085  // Need this conversion so that this method
2086  // will compile for TYPE_SCALAR instantiations
2087  const auto insert_val =
2088  raw_value<typename ProjectionAction::InsertInput>(val, i);
2089 
2090  insert_id(id, insert_val, vertex.processor_id());
2091  }
2092  }
2093  else
2094  {
2095  // C_ZERO elements have a single nodal value DoF at
2096  // vertices. We can't assert n_comp==1 here,
2097  // because if this is a hanging node then it may have
2098  // more face/edge DoFs, but we don't need to deal with
2099  // those here.
2100 
2101  const dof_id_type id = vertex.dof_number(sys_num, var, 0);
2102  insert_id(id, val, vertex.processor_id());
2103  }
2104  }
2105  else if (cont == C_ONE)
2106  {
2107  libmesh_assert(vertex.n_comp(sys_num, var));
2108  const dof_id_type first_id = vertex.dof_number(sys_num, var, 0);
2109 
2110  // C_ONE elements have a single nodal value and dim
2111  // gradient values at vertices, as well as cross
2112  // gradients for HERMITE. We need to have an element in
2113  // hand to figure out dim and to have in case this
2114  // vertex is a new vertex.
2115  const int dim = elem.dim();
2116 #ifndef NDEBUG
2117  // For now all C1 elements at a vertex had better have
2118  // the same dimension. If anyone hits these asserts let
2119  // me know; we could probably support a mixed-dimension
2120  // mesh IFF the 2D elements were all parallel to xy and
2121  // the 1D elements all parallel to x.
2122  for (const auto e_id : (*this->projector.nodes_to_elem)[vertex.id()])
2123  {
2124  const Elem & e = system.get_mesh().elem_ref(e_id);
2125  libmesh_assert_equal_to(dim, e.dim());
2126  }
2127 #endif
2128 #ifdef LIBMESH_ENABLE_AMR
2129  bool is_old_vertex = true;
2130  if (elem.refinement_flag() == Elem::JUST_REFINED)
2131  {
2132  const int i_am_child =
2133  elem.parent()->which_child_am_i(&elem);
2134  is_old_vertex =
2135  elem.parent()->is_vertex_on_parent(i_am_child, n);
2136  }
2137 #else
2138  const bool is_old_vertex = false;
2139 #endif
2140 
2141  // The hermite element vertex shape functions are weird
2142  if (base_fe_type.family == HERMITE)
2143  {
2144  const FValue val =
2145  f.eval_at_node(context,
2146  var_component,
2147  dim,
2148  vertex,
2149  extra_hanging_dofs[var],
2150  system.time);
2151  insert_id(first_id, val, vertex.processor_id());
2152 
2153  typename GFunctor::FunctorValue grad =
2154  is_old_vertex ?
2155  g->eval_at_node(context,
2156  var_component,
2157  dim,
2158  vertex,
2159  extra_hanging_dofs[var],
2160  system.time) :
2161  g->eval_at_point(context,
2162  var_component,
2163  vertex,
2164  system.time,
2165  false);
2166  // x derivative. Use slice because grad may be a tensor type
2167  insert_id(first_id+1, grad.slice(0),
2168  vertex.processor_id());
2169 #if LIBMESH_DIM > 1
2170  if (dim > 1 && is_old_vertex && f.is_grid_projection())
2171  {
2172  for (int i = 1; i < dim; ++i)
2173  insert_id(first_id+i+1, grad.slice(i),
2174  vertex.processor_id());
2175 
2176  // We can directly copy everything else too
2177  std::vector<FValue> derivs;
2178  f.eval_mixed_derivatives
2179  (context, var_component, dim, vertex, derivs);
2180  for (auto i : index_range(derivs))
2181  insert_id(first_id+dim+1+i, derivs[i],
2182  vertex.processor_id());
2183  }
2184  else if (dim > 1)
2185  {
2186  // We'll finite difference mixed derivatives.
2187  // This delta_x used to be TOLERANCE*hmin, but
2188  // the factor of 10 improved the accuracy in
2189  // some unit test projections
2190  Real delta_x = TOLERANCE * 10 * elem.hmin();
2191 
2192  Point nxminus = elem.point(n),
2193  nxplus = elem.point(n);
2194  nxminus(0) -= delta_x;
2195  nxplus(0) += delta_x;
2196  typename GFunctor::FunctorValue gxminus =
2197  g->eval_at_point(context,
2198  var_component,
2199  nxminus,
2200  system.time,
2201  true);
2202  typename GFunctor::FunctorValue gxplus =
2203  g->eval_at_point(context,
2204  var_component,
2205  nxplus,
2206  system.time,
2207  true);
2208  // y derivative
2209  insert_id(first_id+2, grad.slice(1),
2210  vertex.processor_id());
2211  // xy derivative
2212  insert_id(first_id+3,
2213  (grad_component(gxplus, 1) - grad_component(gxminus, 1)) / 2. / delta_x,
2214  vertex.processor_id());
2215 
2216 #if LIBMESH_DIM > 2
2217  if (dim > 2)
2218  {
2219  // z derivative
2220  insert_id(first_id+4, grad.slice(2),
2221  vertex.processor_id());
2222  // xz derivative
2223  insert_id(first_id+5,
2224  (grad_component(gxplus, 2) - grad_component(gxminus, 2)) / 2. / delta_x,
2225  vertex.processor_id());
2226 
2227  // We need new points for yz
2228  Point nyminus = elem.point(n),
2229  nyplus = elem.point(n);
2230  nyminus(1) -= delta_x;
2231  nyplus(1) += delta_x;
2232  typename GFunctor::FunctorValue gyminus =
2233  g->eval_at_point(context,
2234  var_component,
2235  nyminus,
2236  system.time,
2237  true);
2238  typename GFunctor::FunctorValue gyplus =
2239  g->eval_at_point(context,
2240  var_component,
2241  nyplus,
2242  system.time,
2243  true);
2244  // yz derivative
2245  insert_id(first_id+6,
2246  (grad_component(gyplus, 2) - grad_component(gyminus, 2)) / 2. / delta_x,
2247  vertex.processor_id());
2248  // Getting a 2nd order xyz is more tedious
2249  Point nxmym = elem.point(n),
2250  nxmyp = elem.point(n),
2251  nxpym = elem.point(n),
2252  nxpyp = elem.point(n);
2253  nxmym(0) -= delta_x;
2254  nxmym(1) -= delta_x;
2255  nxmyp(0) -= delta_x;
2256  nxmyp(1) += delta_x;
2257  nxpym(0) += delta_x;
2258  nxpym(1) -= delta_x;
2259  nxpyp(0) += delta_x;
2260  nxpyp(1) += delta_x;
2261  typename GFunctor::FunctorValue gxmym =
2262  g->eval_at_point(context,
2263  var_component,
2264  nxmym,
2265  system.time,
2266  true);
2267  typename GFunctor::FunctorValue gxmyp =
2268  g->eval_at_point(context,
2269  var_component,
2270  nxmyp,
2271  system.time,
2272  true);
2273  typename GFunctor::FunctorValue gxpym =
2274  g->eval_at_point(context,
2275  var_component,
2276  nxpym,
2277  system.time,
2278  true);
2279  typename GFunctor::FunctorValue gxpyp =
2280  g->eval_at_point(context,
2281  var_component,
2282  nxpyp,
2283  system.time,
2284  true);
2285  FValue gxzplus = (grad_component(gxpyp, 2) - grad_component(gxmyp, 2))
2286  / 2. / delta_x;
2287  FValue gxzminus = (grad_component(gxpym, 2) - grad_component(gxmym, 2))
2288  / 2. / delta_x;
2289  // xyz derivative
2290  insert_id(first_id+7,
2291  (gxzplus - gxzminus) / 2. / delta_x,
2292  vertex.processor_id());
2293  }
2294 #endif // LIBMESH_DIM > 2
2295  }
2296 #endif // LIBMESH_DIM > 1
2297  }
2298  else
2299  {
2300  // Currently other C_ONE elements have a single nodal
2301  // value shape function and nodal gradient component
2302  // shape functions
2303  libmesh_assert_equal_to(
2305  base_fe_type,
2306  &elem,
2307  elem.get_node_index(&vertex),
2308  this->projector.system.get_dof_map().should_p_refine(
2309  this->projector.system.get_dof_map().var_group_from_var_number(var))),
2310  (unsigned int)(1 + dim));
2311 
2312  const FValue val =
2313  f.eval_at_node(context, var_component, dim,
2314  vertex, extra_hanging_dofs[var],
2315  system.time);
2316  insert_id(first_id, val, vertex.processor_id());
2317  typename GFunctor::FunctorValue grad =
2318  is_old_vertex ?
2319  g->eval_at_node(context, var_component, dim,
2320  vertex, extra_hanging_dofs[var],
2321  system.time) :
2322  g->eval_at_point(context, var_component, vertex,
2323  system.time, false);
2324  for (int i=0; i!= dim; ++i)
2325  insert_id(first_id + i + 1, grad.slice(i),
2326  vertex.processor_id());
2327  }
2328  }
2329  else
2330  libmesh_error_msg("Unknown continuity " << cont);
2331  }
2332  }
2333 }
2334 
2335 
2336 template <typename FFunctor, typename GFunctor,
2337  typename FValue, typename ProjectionAction>
2339  (const node_range & range)
2340 {
2341  LOG_SCOPE ("project_edges","GenericProjector");
2342 
2343  const unsigned int sys_num = system.number();
2344 
2345  for (const auto & e_pair : range)
2346  {
2347  const Elem & elem = *std::get<0>(e_pair.second);
2348 
2349  // If this is an unchanged element then we already copied all
2350  // its dofs
2351 #ifdef LIBMESH_ENABLE_AMR
2352  if (f.is_grid_projection() &&
2353  (elem.refinement_flag() != Elem::JUST_REFINED &&
2357  continue;
2358 #endif // LIBMESH_ENABLE_AMR
2359 
2360  const Node & edge_node = *e_pair.first;
2361  const int dim = elem.dim();
2362  const var_set & edge_vars = std::get<2>(e_pair.second);
2363 
2364  const unsigned short edge_num = std::get<1>(e_pair.second);
2365  const unsigned short node_num = elem.n_vertices() + edge_num;
2366  context.edge = edge_num;
2367  context.pre_fe_reinit(system, &elem);
2368 
2369  this->find_dofs_to_send(edge_node, elem, node_num, edge_vars);
2370 
2371  // Look at all the variables we're supposed to interpolate from
2372  // this element on this edge
2373  for (const auto & var : edge_vars)
2374  {
2375  const Variable & variable = system.variable(var);
2376  const FEType & base_fe_type = variable.type();
2377  const unsigned int var_component =
2379 
2380  if (base_fe_type.family == SCALAR)
2381  continue;
2382 
2383  FEType fe_type = base_fe_type;
2384 
2385  // This may be a p refined element
2386  fe_type.order = fe_type.order + elem.p_level();
2387 
2388  // If this is a Lagrange element with DoFs on edges then by
2389  // convention we interpolate at the node rather than project
2390  // along the edge.
2391  if (fe_type.family == LAGRANGE || fe_type.family == LAGRANGE_VEC)
2392  {
2393  if (fe_type.order > 1)
2394  {
2395  const processor_id_type pid =
2396  edge_node.processor_id();
2397  FValue fval = f.eval_at_point
2398  (context, var_component, edge_node, system.time,
2399  false);
2400  if (fe_type.family == LAGRANGE_VEC)
2401  {
2402  // We will have a number of nodal value DoFs equal to the elem dim
2403  for (auto i : make_range(elem.dim()))
2404  {
2405  const dof_id_type dof_id =
2406  edge_node.dof_number(sys_num, var, i);
2407 
2408  // Need this conversion so that this method
2409  // will compile for TYPE_SCALAR instantiations
2410  const auto insert_val =
2411  raw_value<typename ProjectionAction::InsertInput>(fval, i);
2412 
2413  insert_id(dof_id, insert_val, pid);
2414  }
2415  }
2416  else // We are LAGRANGE
2417  {
2418  const dof_id_type dof_id =
2419  edge_node.dof_number(sys_num, var, 0);
2420  insert_id(dof_id, fval, pid);
2421  }
2422  }
2423  continue;
2424  }
2425 
2426 #ifdef LIBMESH_ENABLE_AMR
2427  // If this is a low order monomial element which has merely
2428  // been h refined then we already copied all its dofs
2429  if (fe_type.family == MONOMIAL &&
2430  fe_type.order == CONSTANT &&
2433  continue;
2434 #endif // LIBMESH_ENABLE_AMR
2435 
2436  // FIXME: Need to generalize this to vector-valued elements. [PB]
2438  context.get_element_fe( var, fe, dim );
2439  FEGenericBase<typename FFunctor::RealType> * edge_fe = nullptr;
2440  context.get_edge_fe( var, edge_fe );
2441 
2442  // If we're JUST_COARSENED we'll need a custom
2443  // evaluation, not just the standard edge FE
2445 #ifdef LIBMESH_ENABLE_AMR
2447  *fe :
2448 #endif
2449  *edge_fe;
2450 
2451 #ifdef LIBMESH_ENABLE_AMR
2452  if (elem.refinement_flag() == Elem::JUST_COARSENED)
2453  {
2454  std::vector<Point> fine_points;
2455 
2456  std::unique_ptr<FEGenericBase<typename FFunctor::RealType>> fine_fe
2458 
2459  std::unique_ptr<QBase> qrule
2460  (base_fe_type.default_quadrature_rule(1));
2461  fine_fe->attach_quadrature_rule(qrule.get());
2462 
2463  const std::vector<Point> & child_xyz =
2464  fine_fe->get_xyz();
2465 
2466  for (unsigned int c = 0, nc = elem.n_children();
2467  c != nc; ++c)
2468  {
2469  if (!elem.is_child_on_edge(c, context.edge))
2470  continue;
2471 
2472  fine_fe->edge_reinit(elem.child_ptr(c), context.edge);
2473  fine_points.insert(fine_points.end(),
2474  child_xyz.begin(),
2475  child_xyz.end());
2476  }
2477 
2478  std::vector<Point> fine_qp;
2479  FEMap::inverse_map (dim, &elem, fine_points, fine_qp);
2480 
2481  context.elem_fe_reinit(&fine_qp);
2482  }
2483  else
2484 #endif // LIBMESH_ENABLE_AMR
2485  context.edge_fe_reinit();
2486 
2487  const std::vector<dof_id_type> & dof_indices =
2488  context.get_dof_indices(var);
2489 
2490  std::vector<unsigned int> edge_dofs;
2491  FEInterface::dofs_on_edge(&elem, dim, base_fe_type,
2492  context.edge, edge_dofs);
2493 
2494  this->construct_projection
2495  (dof_indices, edge_dofs, var_component,
2496  &edge_node, proj_fe);
2497  }
2498  }
2499 }
2500 
2501 
2502 template <typename FFunctor, typename GFunctor,
2503  typename FValue, typename ProjectionAction>
2505  (const node_range & range)
2506 {
2507  LOG_SCOPE ("project_sides","GenericProjector");
2508 
2509  const unsigned int sys_num = system.number();
2510 
2511  for (const auto & s_pair : range)
2512  {
2513  const Elem & elem = *std::get<0>(s_pair.second);
2514 
2515  // If this is an unchanged element then we already copied all
2516  // its dofs
2517 #ifdef LIBMESH_ENABLE_AMR
2518  if (f.is_grid_projection() &&
2519  (elem.refinement_flag() != Elem::JUST_REFINED &&
2523  continue;
2524 #endif // LIBMESH_ENABLE_AMR
2525 
2526  const Node & side_node = *s_pair.first;
2527  const int dim = elem.dim();
2528  const var_set & side_vars = std::get<2>(s_pair.second);
2529 
2530  const unsigned int side_num = std::get<1>(s_pair.second);
2531  unsigned short node_num = elem.n_vertices()+side_num;
2532  // In 3D only some sides may have nodes
2533  if (dim == 3)
2534  for (auto n : make_range(elem.n_nodes()))
2535  {
2536  if (!elem.is_face(n))
2537  continue;
2538 
2539  if (elem.is_node_on_side(n, side_num))
2540  {
2541  node_num = n;
2542  break;
2543  }
2544  }
2545 
2546  context.side = side_num;
2547  context.pre_fe_reinit(system, &elem);
2548 
2549  this->find_dofs_to_send(side_node, elem, node_num, side_vars);
2550 
2551  // Look at all the variables we're supposed to interpolate from
2552  // this element on this side
2553  for (const auto & var : side_vars)
2554  {
2555  const Variable & variable = system.variable(var);
2556  const FEType & base_fe_type = variable.type();
2557  const unsigned int var_component =
2559 
2560  if (base_fe_type.family == SCALAR)
2561  continue;
2562 
2563  FEType fe_type = base_fe_type;
2564  const auto & dof_map = system.get_dof_map();
2565  const bool add_p_level = dof_map.should_p_refine_var(var);
2566 
2567  // This may be a p refined element
2568  fe_type.order = fe_type.order + add_p_level*elem.p_level();
2569 
2570  // If this is a Lagrange element with DoFs on sides then by
2571  // convention we interpolate at the node rather than project
2572  // along the side.
2573  if (fe_type.family == LAGRANGE || fe_type.family == LAGRANGE_VEC)
2574  {
2575  // If order==1 then we only have DoFs on vertices, so
2576  // skip all this.
2577  // If order>1 ... we might still have something tricky
2578  // like order==2 PRISM20, where some side nodes have
2579  // DoFs but others don't. Gotta check.
2580  if (fe_type.order > 1 &&
2581  side_node.n_comp(sys_num, var))
2582  {
2583  const processor_id_type pid =
2584  side_node.processor_id();
2585  FValue fval = f.eval_at_point
2586  (context, var_component, side_node, system.time,
2587  false);
2588 
2589  if (fe_type.family == LAGRANGE_VEC)
2590  {
2591  // We will have a number of nodal value DoFs equal to the elem dim
2592  for (auto i : make_range(elem.dim()))
2593  {
2594  const dof_id_type dof_id = side_node.dof_number(sys_num, var, i);
2595 
2596  // Need this conversion so that this method
2597  // will compile for TYPE_SCALAR instantiations
2598  const auto insert_val =
2599  raw_value<typename ProjectionAction::InsertInput>(fval, i);
2600 
2601  insert_id(dof_id, insert_val, pid);
2602  }
2603  }
2604  else // We are LAGRANGE
2605  {
2606  const dof_id_type dof_id =
2607  side_node.dof_number(sys_num, var, 0);
2608  insert_id(dof_id, fval, pid);
2609  }
2610  }
2611  continue;
2612  }
2613 
2614 #ifdef LIBMESH_ENABLE_AMR
2615  // If this is a low order monomial element which has merely
2616  // been h refined then we already copied all its dofs
2617  if (fe_type.family == MONOMIAL &&
2618  fe_type.order == CONSTANT &&
2621  continue;
2622 #endif // LIBMESH_ENABLE_AMR
2623 
2624  // FIXME: Need to generalize this to vector-valued elements. [PB]
2626  context.get_element_fe( var, fe, dim );
2627  FEGenericBase<typename FFunctor::RealType> * side_fe = nullptr;
2628  context.get_side_fe( var, side_fe );
2629 
2630  // If we're JUST_COARSENED we'll need a custom
2631  // evaluation, not just the standard side FE
2633 #ifdef LIBMESH_ENABLE_AMR
2635  *fe :
2636 #endif
2637  *side_fe;
2638 
2639 #ifdef LIBMESH_ENABLE_AMR
2640  if (elem.refinement_flag() == Elem::JUST_COARSENED)
2641  {
2642  std::vector<Point> fine_points;
2643 
2644  std::unique_ptr<FEGenericBase<typename FFunctor::RealType>> fine_fe
2646 
2647  std::unique_ptr<QBase> qrule
2648  (base_fe_type.default_quadrature_rule(1));
2649  fine_fe->attach_quadrature_rule(qrule.get());
2650 
2651  const std::vector<Point> & child_xyz =
2652  fine_fe->get_xyz();
2653 
2654  for (unsigned int c = 0, nc = elem.n_children();
2655  c != nc; ++c)
2656  {
2657  if (!elem.is_child_on_side(c, context.side))
2658  continue;
2659 
2660  fine_fe->reinit(elem.child_ptr(c), context.side);
2661  fine_points.insert(fine_points.end(),
2662  child_xyz.begin(),
2663  child_xyz.end());
2664  }
2665 
2666  std::vector<Point> fine_qp;
2667  FEMap::inverse_map (dim, &elem, fine_points, fine_qp);
2668 
2669  context.elem_fe_reinit(&fine_qp);
2670  }
2671  else
2672 #endif // LIBMESH_ENABLE_AMR
2673  context.side_fe_reinit();
2674 
2675  const std::vector<dof_id_type> & dof_indices =
2676  context.get_dof_indices(var);
2677 
2678  std::vector<unsigned int> side_dofs;
2679  FEInterface::dofs_on_side(&elem, dim, base_fe_type,
2680  context.side, side_dofs, add_p_level);
2681 
2682  this->construct_projection
2683  (dof_indices, side_dofs, var_component,
2684  &side_node, proj_fe);
2685  }
2686  }
2687 }
2688 
2689 
2690 template <typename FFunctor, typename GFunctor,
2691  typename FValue, typename ProjectionAction>
2693  (const interior_range & range)
2694 {
2695  LOG_SCOPE ("project_interiors","GenericProjector");
2696 
2697  const unsigned int sys_num = system.number();
2698 
2699  // Iterate over all dof-bearing element interiors in the range
2700  for (const auto & elem : range)
2701  {
2702  unsigned char dim = cast_int<unsigned char>(elem->dim());
2703 
2704  context.pre_fe_reinit(system, elem);
2705 
2706  // Loop over all the variables we've been requested to project, to
2707  // do the projection
2708  for (const auto & var : this->projector.variables)
2709  {
2710  const Variable & variable = system.variable(var);
2711 
2712  if (!variable.active_on_subdomain(elem->subdomain_id()))
2713  continue;
2714 
2715  const FEType & base_fe_type = variable.type();
2716 
2717  if (base_fe_type.family == SCALAR)
2718  continue;
2719 
2721  context.get_element_fe( var, fe, dim );
2722 
2723  FEType fe_type = base_fe_type;
2724  const auto & dof_map = system.get_dof_map();
2725  const auto vg = dof_map.var_group_from_var_number(var);
2726  const bool add_p_level = dof_map.should_p_refine(vg);
2727 
2728  // This may be a p refined element
2729  fe_type.order = fe_type.order + add_p_level * elem->p_level();
2730 
2731  const unsigned int var_component =
2733 
2734  // If this is a Lagrange element with interior DoFs then by
2735  // convention we interpolate at the interior node rather
2736  // than project along the interior.
2737  if (fe_type.family == LAGRANGE || fe_type.family == LAGRANGE_VEC)
2738  {
2739  if (fe_type.order > 1)
2740  {
2741  const unsigned int first_interior_node =
2742  (elem->n_vertices() +
2743  ((elem->dim() > 2) * elem->n_edges()) +
2744  ((elem->dim() > 1) * elem->n_sides()));
2745  const unsigned int n_nodes = elem->n_nodes();
2746 
2747  // < vs != is important here for HEX20, QUAD8!
2748  for (unsigned int n = first_interior_node; n < n_nodes; ++n)
2749  {
2750  const Node & interior_node = elem->node_ref(n);
2751 
2752  FValue fval = f.eval_at_point
2753  (context, var_component, interior_node,
2754  system.time, false);
2755  const processor_id_type pid =
2756  interior_node.processor_id();
2757 
2758  if (fe_type.family == LAGRANGE_VEC)
2759  {
2760  // We will have a number of nodal value DoFs equal to the elem dim
2761  for (unsigned int i = 0; i < elem->dim(); ++i)
2762  {
2763  const dof_id_type dof_id = interior_node.dof_number(sys_num, var, i);
2764 
2765  // Need this conversion so that this method
2766  // will compile for TYPE_SCALAR instantiations
2767  const auto insert_val =
2768  raw_value<typename ProjectionAction::InsertInput>(fval, i);
2769 
2770  insert_id(dof_id, insert_val, pid);
2771  }
2772  }
2773  else // We are LAGRANGE
2774  {
2775  const dof_id_type dof_id =
2776  interior_node.dof_number(sys_num, var, 0);
2777  insert_id(dof_id, fval, pid);
2778  }
2779  }
2780  }
2781  continue;
2782  }
2783 
2784 #ifdef LIBMESH_ENABLE_AMR
2785  if (elem->refinement_flag() == Elem::JUST_COARSENED)
2786  {
2787  std::vector<Point> fine_points;
2788 
2789  std::unique_ptr<FEGenericBase<typename FFunctor::RealType>> fine_fe
2791 
2792  std::unique_ptr<QBase> qrule
2793  (base_fe_type.default_quadrature_rule(dim));
2794  fine_fe->attach_quadrature_rule(qrule.get());
2795 
2796  const std::vector<Point> & child_xyz =
2797  fine_fe->get_xyz();
2798 
2799  for (auto & child : elem->child_ref_range())
2800  {
2801  fine_fe->reinit(&child);
2802  fine_points.insert(fine_points.end(),
2803  child_xyz.begin(),
2804  child_xyz.end());
2805  }
2806 
2807  std::vector<Point> fine_qp;
2808  FEMap::inverse_map (dim, elem, fine_points, fine_qp);
2809 
2810  context.elem_fe_reinit(&fine_qp);
2811  }
2812  else
2813 #endif // LIBMESH_ENABLE_AMR
2814  context.elem_fe_reinit();
2815 
2816  const std::vector<dof_id_type> & dof_indices =
2817  context.get_dof_indices(var);
2818 
2819  const unsigned int n_dofs =
2820  cast_int<unsigned int>(dof_indices.size());
2821 
2822  std::vector<unsigned int> all_dofs(n_dofs);
2823  std::iota(all_dofs.begin(), all_dofs.end(), 0);
2824 
2825  this->construct_projection
2826  (dof_indices, all_dofs, var_component,
2827  nullptr, *fe);
2828  } // end variables loop
2829  } // end elements loop
2830 }
2831 
2832 
2833 
2834 template <typename FFunctor, typename GFunctor,
2835  typename FValue, typename ProjectionAction>
2836 void
2837 GenericProjector<FFunctor, GFunctor, FValue,
2838  ProjectionAction>::SubFunctor::find_dofs_to_send
2839  (const Node & node, const Elem & elem, unsigned short node_num, const var_set & vars)
2840 {
2841  libmesh_assert (&node == elem.node_ptr(node_num));
2842 
2843  // Any ghosted node in our range might have an owner who needs our
2844  // data
2845  const processor_id_type owner = node.processor_id();
2846  if (owner != system.processor_id())
2847  {
2848  const MeshBase & mesh = system.get_mesh();
2849  const DofMap & dof_map = system.get_dof_map();
2850 
2851  // But let's check and see if we can be certain the owner can
2852  // compute any or all of its own dof coefficients on that node.
2853  std::vector<dof_id_type> node_dof_ids, patch_dof_ids;
2854  for (const auto & var : vars)
2855  {
2856  const Variable & variable = system.variable(var);
2857 
2858  if (!variable.active_on_subdomain(elem.subdomain_id()))
2859  continue;
2860 
2861  dof_map.dof_indices(elem, node_num, node_dof_ids, var);
2862  }
2863  libmesh_assert(std::is_sorted(node_dof_ids.begin(),
2864  node_dof_ids.end()));
2865  const std::vector<dof_id_type> & patch =
2866  (*this->projector.nodes_to_elem)[node.id()];
2867  for (const auto & elem_id : patch)
2868  {
2869  const Elem & patch_elem = mesh.elem_ref(elem_id);
2870  if (!patch_elem.active() || owner != patch_elem.processor_id())
2871  continue;
2872  dof_map.dof_indices(&patch_elem, patch_dof_ids);
2873  std::sort(patch_dof_ids.begin(), patch_dof_ids.end());
2874 
2875  // Remove any node_dof_ids that we see can be calculated on
2876  // this element
2877  std::vector<dof_id_type> diff_ids(node_dof_ids.size());
2878  auto it = std::set_difference(node_dof_ids.begin(), node_dof_ids.end(),
2879  patch_dof_ids.begin(), patch_dof_ids.end(), diff_ids.begin());
2880  diff_ids.resize(it-diff_ids.begin());
2881  node_dof_ids.swap(diff_ids);
2882  if (node_dof_ids.empty())
2883  break;
2884  }
2885 
2886  // Give ids_to_push default invalid pid: not yet computed
2887  for (auto id : node_dof_ids)
2888  new_ids_to_push[id].second = DofObject::invalid_processor_id;
2889  }
2890 }
2891 
2892 
2893 
2894 template <typename FFunctor, typename GFunctor,
2895  typename FValue, typename ProjectionAction>
2896 template <typename Value>
2897 void
2898 GenericProjector<FFunctor, GFunctor, FValue,
2899  ProjectionAction>::send_and_insert_dof_values
2900  (std::unordered_map<dof_id_type, std::pair<Value, processor_id_type>> & ids_to_push,
2901  ProjectionAction & action) const
2902 {
2903  LOG_SCOPE ("send_and_insert_dof_values", "GenericProjector");
2904 
2905  // See if we calculated any ids that need to be pushed; get them
2906  // ready to push.
2907  std::unordered_map<processor_id_type, std::vector<dof_id_type>>
2908  pushed_dof_ids, received_dof_ids;
2909  std::unordered_map<processor_id_type, std::vector<typename TypeToSend<Value>::type>>
2910  pushed_dof_values, received_dof_values;
2911  for (auto & id_val_pid : ids_to_push)
2912  {
2913  processor_id_type pid = id_val_pid.second.second;
2915  {
2916  pushed_dof_ids[pid].push_back(id_val_pid.first);
2917  pushed_dof_values[pid].push_back(convert_to_send(id_val_pid.second.first));
2918  }
2919  }
2920 
2921  // If and when we get ids pushed to us, act on them if we have
2922  // corresponding values or save them if not
2923  auto ids_action_functor =
2924  [&action, &received_dof_ids, &received_dof_values]
2925  (processor_id_type pid,
2926  const std::vector<dof_id_type> & data)
2927  {
2928  auto iter = received_dof_values.find(pid);
2929  if (iter == received_dof_values.end())
2930  {
2931  libmesh_assert(received_dof_ids.find(pid) ==
2932  received_dof_ids.end());
2933  received_dof_ids[pid] = data;
2934  }
2935  else
2936  {
2937  auto & vals = iter->second;
2938  std::size_t vals_size = vals.size();
2939  libmesh_assert_equal_to(vals_size, data.size());
2940  for (std::size_t i=0; i != vals_size; ++i)
2941  {
2942  Value val;
2943  convert_from_receive(vals[i], val);
2944  action.insert(data[i], val);
2945  }
2946  received_dof_values.erase(iter);
2947  }
2948  };
2949 
2950  // If and when we get values pushed to us, act on them if we have
2951  // corresponding ids or save them if not
2952  auto values_action_functor =
2953  [&action, &received_dof_ids, &received_dof_values]
2954  (processor_id_type pid,
2955  const std::vector<typename TypeToSend<Value>::type> & data)
2956  {
2957  auto iter = received_dof_ids.find(pid);
2958  if (iter == received_dof_ids.end())
2959  {
2960  // We get no more than 1 values vector from anywhere
2961  libmesh_assert(received_dof_values.find(pid) ==
2962  received_dof_values.end());
2963  received_dof_values[pid] = data;
2964  }
2965  else
2966  {
2967  auto & ids = iter->second;
2968  std::size_t ids_size = ids.size();
2969  libmesh_assert_equal_to(ids_size, data.size());
2970  for (std::size_t i=0; i != ids_size; ++i)
2971  {
2972  Value val;
2973  convert_from_receive(data[i], val);
2974  action.insert(ids[i], val);
2975  }
2976  received_dof_ids.erase(iter);
2977  }
2978  };
2979 
2980  Parallel::push_parallel_vector_data
2981  (system.comm(), pushed_dof_ids, ids_action_functor);
2982 
2983  Parallel::push_parallel_vector_data
2984  (system.comm(), pushed_dof_values, values_action_functor);
2985 
2986  // At this point we shouldn't have any unprocessed data left
2987  libmesh_assert(received_dof_ids.empty());
2988  libmesh_assert(received_dof_values.empty());
2989 
2990 }
2991 
2992 
2993 
2994 template <typename FFunctor, typename GFunctor,
2995  typename FValue, typename ProjectionAction>
2996 void
2998  (const std::vector<dof_id_type> & dof_indices_var,
2999  const std::vector<unsigned int> & involved_dofs,
3000  unsigned int var_component,
3001  const Node * node,
3003 {
3004  const auto & JxW = fe.get_JxW();
3005  const auto & phi = fe.get_phi();
3006  const std::vector<std::vector<typename FEGenericBase<typename FFunctor::RealType>::OutputGradient>> * dphi = nullptr;
3007  const std::vector<Point> & xyz_values = fe.get_xyz();
3008  const FEContinuity cont = fe.get_continuity();
3009  const std::unordered_map<dof_id_type, typename FFunctor::ValuePushType> & ids_to_save =
3010  this->projector.ids_to_save;
3011 
3012  if (cont == C_ONE)
3013  dphi = &(fe.get_dphi());
3014 
3015  const unsigned int n_involved_dofs =
3016  cast_int<unsigned int>(involved_dofs.size());
3017 
3018  std::vector<dof_id_type> free_dof_ids;
3019  DenseVector<typename FFunctor::ValuePushType> Uinvolved(n_involved_dofs);
3020  std::vector<char> dof_is_fixed(n_involved_dofs, false); // bools
3021 
3022  for (auto i : make_range(n_involved_dofs))
3023  {
3024  const dof_id_type id = dof_indices_var[involved_dofs[i]];
3025  auto iter = ids_to_save.find(id);
3026  if (iter == ids_to_save.end())
3027  free_dof_ids.push_back(id);
3028  else
3029  {
3030  dof_is_fixed[i] = true;
3031  Uinvolved(i) = iter->second;
3032  }
3033  }
3034 
3035  const unsigned int free_dofs = free_dof_ids.size();
3036 
3037  // There may be nothing to project
3038  if (!free_dofs)
3039  return;
3040 
3041  // The element matrix and RHS for projections.
3042  // Note that Ke is always real-valued, whereas
3043  // Fe may be complex valued if complex number
3044  // support is enabled
3045  DenseMatrix<Real> Ke(free_dofs, free_dofs);
3047  // The new degree of freedom coefficients to solve for
3049 
3050  const unsigned int n_qp =
3051  cast_int<unsigned int>(xyz_values.size());
3052 
3053  // Loop over the quadrature points
3054  for (unsigned int qp=0; qp<n_qp; qp++)
3055  {
3056  // solution at the quadrature point
3057  FValue fineval = f.eval_at_point(context,
3058  var_component,
3059  xyz_values[qp],
3060  system.time,
3061  false);
3062  // solution grad at the quadrature point
3063  typename GFunctor::FunctorValue finegrad;
3064  if (cont == C_ONE)
3065  finegrad = g->eval_at_point(context,
3066  var_component,
3067  xyz_values[qp],
3068  system.time,
3069  false);
3070 
3071  // Form edge projection matrix
3072  for (unsigned int sidei=0, freei=0;
3073  sidei != n_involved_dofs; ++sidei)
3074  {
3075  unsigned int i = involved_dofs[sidei];
3076  // fixed DoFs aren't test functions
3077  if (dof_is_fixed[sidei])
3078  continue;
3079  for (unsigned int sidej=0, freej=0;
3080  sidej != n_involved_dofs; ++sidej)
3081  {
3082  unsigned int j = involved_dofs[sidej];
3083  if (dof_is_fixed[sidej])
3084  Fe(freei) -= phi[i][qp] * phi[j][qp] *
3085  JxW[qp] * Uinvolved(sidej);
3086  else
3087  Ke(freei,freej) += phi[i][qp] *
3088  phi[j][qp] * JxW[qp];
3089  if (cont == C_ONE)
3090  {
3091  if (dof_is_fixed[sidej])
3092  Fe(freei) -= ( TensorTools::inner_product((*dphi)[i][qp],
3093  (*dphi)[j][qp]) ) *
3094  JxW[qp] * Uinvolved(sidej);
3095  else
3096  Ke(freei,freej) += ( TensorTools::inner_product((*dphi)[i][qp],
3097  (*dphi)[j][qp]) )
3098  * JxW[qp];
3099  }
3100  if (!dof_is_fixed[sidej])
3101  freej++;
3102  }
3103  Fe(freei) += phi[i][qp] * fineval * JxW[qp];
3104  if (cont == C_ONE)
3105  Fe(freei) += (TensorTools::inner_product(finegrad,
3106  (*dphi)[i][qp]) ) *
3107  JxW[qp];
3108  freei++;
3109  }
3110  }
3111 
3112  Ke.cholesky_solve(Fe, Ufree);
3113 
3114  // Transfer new edge solutions to element
3115  const processor_id_type pid = node ?
3117  insert_ids(free_dof_ids, Ufree.get_values(), pid);
3118 }
3119 
3120 
3121 } // namespace libMesh
3122 
3123 #endif // GENERIC_PROJECTOR_H
The GenericProjector class implements the core of other projection operations, using two input functo...
static unsigned int n_dofs_per_elem(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition: fe_interface.C:530
std::unordered_map< dof_id_type, std::pair< typename FFunctor::ValuePushType, processor_id_type > > new_ids_to_push
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:196
void send_and_insert_dof_values(std::unordered_map< dof_id_type, std::pair< Value, processor_id_type >> &ids_to_push, ProjectionAction &action) const
virtual bool is_vertex_on_parent(unsigned int c, unsigned int n) const
Definition: elem.C:3487
FEFamily family
The type of finite element.
Definition: fe_type.h:221
void check_old_context(const FEMContext &c)
Real time
For time-dependent problems, this is the time t at the beginning of the current timestep.
Definition: system.h:1615
OldSolutionBase(const libMesh::System &sys_in, const std::vector< unsigned int > *vars)
RefinementState refinement_flag() const
Definition: elem.h:3210
dof_id_type dof_number(const unsigned int s, const unsigned int var, const unsigned int comp) const
Definition: dof_object.h:1032
const Elem * parent() const
Definition: elem.h:3030
unsigned int variable_scalar_number(std::string_view var, unsigned int component) const
Definition: system.h:2489
void eval_old_dofs(const Elem &elem, const FEType &fe_type, unsigned int sys_num, unsigned int var_num, std::vector< dof_id_type > &indices, std::vector< DofValueType > &values)
OldSolutionValue(const OldSolutionValue &in)
A Node is like a Point, but with more information.
Definition: node.h:52
const Variable & variable(unsigned int var) const
Return a constant reference to Variable var.
Definition: system.h:2458
bool check_old_context(const FEMContext &c, const Point &p)
const NumericVector< Number > & old_solution
void eval_old_dofs(const Elem &elem, unsigned int node_num, unsigned int var_num, std::vector< dof_id_type > &indices, std::vector< DofValueType > &values)
unsigned int get_node_index(const Node *node_ptr) const
Definition: elem.h:2545
unsigned int n_comp(const unsigned int s, const unsigned int var) const
Definition: dof_object.h:1002
void operator()(const node_range &range)
void construct_projection(const std::vector< dof_id_type > &dof_indices_var, const std::vector< unsigned int > &involved_dofs, unsigned int var_component, const Node *node, const FEGenericBase< typename FFunctor::RealType > &fe)
std::unique_ptr< GFunctor > master_g_deepcopy
Needed for C1 type elements only.
virtual bool is_face(const unsigned int i) const =0
NodesToElemMap * nodes_to_elem
void init_context(FEMContext &c)
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Definition: dof_map.C:2164
static constexpr Real TOLERANCE
const Elem & get_elem() const
Accessor for current Elem object.
Definition: fem_context.h:908
unsigned int dim
static Point inverse_map(const unsigned int dim, const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true, const bool extra_checks=true)
Definition: fe_map.C:1628
ProjectSides(ProjectSides &p_s, Threads::split)
void libmesh_merge_move(T &target, T &source)
The OldSolutionBase input functor abstract base class is the root of the OldSolutionValue and OldSolu...
virtual bool is_child_on_side(const unsigned int c, const unsigned int s) const =0
void eval_mixed_derivatives(const FEMContext &, unsigned int, unsigned int, const Node &, std::vector< Output > &)
void eval_mixed_derivatives(const FEMContext &, unsigned int, unsigned int, const Node &, std::vector< Output > &)
Dummy "splitting object" used to distinguish splitting constructors from copy constructors.
Definition: threads_none.h:63
RefinementState p_refinement_flag() const
Definition: elem.h:3226
static void dofs_on_edge(const Elem *const elem, const unsigned int dim, const FEType &fe_t, unsigned int e, std::vector< unsigned int > &di, const bool add_p_level=true)
Fills the vector di with the local degree of freedom indices associated with edge e of element elem A...
Definition: fe_interface.C:611
std::vector< T > & get_values()
Definition: dense_vector.h:278
FEMFunctionWrapper(const FEMFunctionWrapper< Output > &fw)
std::unordered_multimap< const Node *, std::tuple< const Elem *, unsigned short, var_set > > ves_multimap
ProjectionAction & master_action
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
MeshBase & mesh
FEMFunctionWrapper(const FEMFunctionBase< Output > &f)
static FEFieldType field_type(const FEType &fe_type)
RefinementState
Enumeration of possible element refinement states.
Definition: elem.h:1452
virtual bool is_node_on_side(const unsigned int n, const unsigned int s) const =0
const Parallel::Communicator & comm() const
virtual unsigned int n_children() const =0
unsigned int p_level() const
Definition: elem.h:3108
OrderWrapper order
The approximation order of the element.
Definition: fe_type.h:215
void set_algebraic_type(const AlgebraicType atype)
Setting which determines whether to initialize algebraic structures (elem_*) on each element and set ...
Definition: fem_context.h:986
void build_nodes_to_elem_map(const MeshBase &mesh, std::vector< std::vector< dof_id_type >> &nodes_to_elem_map)
After calling this function the input vector nodes_to_elem_map will contain the node to element conne...
Definition: mesh_tools.C:449
ProjectVertices(ProjectVertices &p_v, Threads::split)
This class defines a vector in LIBMESH_DIM dimensional Real or Complex space.
The StoredRange class defines a contiguous, divisible set of objects.
Definition: stored_range.h:54
The libMesh namespace provides an interface to certain functionality in the library.
The OldSolutionValue input functor class can be used with GenericProjector to read values from a solu...
void eval_old_dofs(const Elem &, const FEType &, unsigned int, unsigned int, std::vector< dof_id_type > &, std::vector< Output > &)
virtual Real hmax() const
Definition: elem.C:722
NumericVector< Val > & target_vector
VectorSetAction(NumericVector< Val > &target_vec)
const MeshBase & get_mesh() const
Definition: system.h:2358
void insert_id(dof_id_type id, const InsertInput &val, processor_id_type pid)
void operator()(const interior_range &range)
void join(const SubFunctor &other)
OldSolutionBase(const OldSolutionBase &in)
void iota(ForwardIter first, ForwardIter last, T value)
Utility::iota was created back when std::iota was just an SGI STL extension.
Definition: utility.h:229
virtual void request_phi() const =0
request phi calculations
uint8_t processor_id_type
Definition: id_types.h:104
This is the MeshBase class.
Definition: mesh_base.h:75
SimpleRange< ChildRefIter > child_ref_range()
Returns a range with all children of a parent element, usable in range-based for loops.
Definition: elem.h:2346
virtual FEContinuity get_continuity() const =0
bool should_p_refine_var(unsigned int var) const
Whether the given variable should be p-refined.
Definition: dof_map.h:2437
unsigned int var_group_from_var_number(unsigned int var_num) const
Definition: dof_map.h:2430
FEType get_fe_type() const
Definition: fe_abstract.h:520
std::vector< const Elem * > interiors
static void get_shape_outputs(FEAbstract &fe)
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:179
std::unique_ptr< QBase > default_quadrature_rule(const unsigned int dim, const int extraorder=0) const
Definition: fe_type.C:34
const std::vector< std::vector< OutputGradient > > & get_dphi() const
Definition: fe_base.h:230
void eval_old_dofs(const Elem &, unsigned int, unsigned int, std::vector< dof_id_type > &, std::vector< Output > &)
const dof_id_type n_nodes
Definition: tecplot_io.C:67
unsigned int number() const
Definition: system.h:2350
std::unordered_map< dof_id_type, typename FFunctor::ValuePushType > ids_to_save
This class defines the notion of a variable in the system.
Definition: variable.h:50
dof_id_type id() const
Definition: dof_object.h:828
dof_id_type numeric_index_type
Definition: id_types.h:99
OldSolutionValue(const libMesh::System &sys_in, const NumericVector< Number > &old_sol, const std::vector< unsigned int > *vars)
virtual Real hmin() const
Definition: elem.C:702
void project(const ConstElemRange &range)
Function definitions.
static bool extra_hanging_dofs(const FEType &fe_t)
TensorTools::MakeBaseNumber< Output >::type DofValueType
virtual unsigned int n_nodes() const =0
void insert_ids(const std::vector< dof_id_type > &ids, const std::vector< InsertInput > &vals, processor_id_type pid)
bool is_sorted(const std::vector< KeyType > &v)
void operator()(const ConstElemRange &range)
static const processor_id_type invalid_processor_id
An invalid processor_id to distinguish DoFs that have not been assigned to a processor.
Definition: dof_object.h:493
unsigned int n_vars(const unsigned int s, const unsigned int vg) const
Definition: dof_object.h:967
StoredRange< std::vector< node_projection >::const_iterator, node_projection > node_range
Manages consistently variables, degrees of freedom, and coefficient vectors.
Definition: system.h:96
unsigned int n_systems() const
Definition: dof_object.h:937
std::unordered_map< dof_id_type, std::vector< dof_id_type > > NodesToElemMap
Convenience typedef for the Node-to-attached-Elem mapping that may be passed in to the constructor...
SimpleRange< IndexType > as_range(const std::pair< IndexType, IndexType > &p)
Helper function that allows us to treat a homogenous pair as a range.
Definition: simple_range.h:57
void insert(dof_id_type id, Val val)
ProjectEdges(ProjectEdges &p_e, Threads::split)
libmesh_assert(ctx)
std::vector< FEFieldType > field_types
static const dof_id_type invalid_id
An invalid id to distinguish an uninitialized DofObject.
Definition: dof_object.h:482
virtual_for_inffe const std::vector< Real > & get_JxW() const
Definition: fe_abstract.h:303
This class provides all data required for a physics package (e.g.
Definition: fem_context.h:62
void init_context(FEMContext &c)
std::unique_ptr< FEMFunctionBase< Output > > _f
virtual_for_inffe const std::vector< Point > & get_xyz() const
Definition: fe_abstract.h:280
DofObject * get_old_dof_object()
Pointer accessor for previously public old_dof_object.
Definition: dof_object.h:96
bool active_on_subdomain(subdomain_id_type sid) const
Definition: variable.h:167
void insert(const std::vector< dof_id_type > &dof_indices, const DenseVector< Val > &Ue)
const std::vector< unsigned int > & variables
void eval_mixed_derivatives(const FEMContext &libmesh_dbg_var(c), unsigned int i, unsigned int dim, const Node &n, std::vector< Output > &derivs)
unsigned int which_child_am_i(const Elem *e) const
Definition: elem.h:3192
DofObject & get_old_dof_object_ref()
As above, but do not use in situations where the old_dof_object may be nullptr, since this function a...
Definition: dof_object.h:104
Output eval_at_point(const FEMContext &c, unsigned int i, const Point &n, const Real time, bool)
virtual void request_dphi() const =0
request dphi calculations
static unsigned int n_dofs_at_node(const unsigned int dim, const FEType &fe_t, const ElemType t, const unsigned int n)
Definition: fe_interface.C:436
void operator()(const node_range &range)
TensorTools::MakeReal< Output >::type RealType
virtual numeric_index_type first_local_index() const =0
GenericProjector(const System &system_in, FFunctor &f_in, GFunctor *g_in, ProjectionAction &act_in, const std::vector< unsigned int > &variables_in, NodesToElemMap *nodes_to_elem_in=nullptr)
void find_dofs_to_send(const Node &node, const Elem &elem, unsigned short node_num, const var_set &vars)
virtual unsigned int n_vertices() const =0
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static const Real out_of_elem_tol
void convert_from_receive(SendT &received, T &converted)
subdomain_id_type subdomain_id() const
Definition: elem.h:2582
NumberTensorValue Tensor
virtual unsigned short dim() const =0
ProjectInteriors(ProjectInteriors &p_i, Threads::split)
FEContinuity
defines an enum for finite element types to libmesh_assert a certain level (or type? Hcurl?) of continuity.
System * system() const
Definition: variable.h:114
const Node * node_ptr(const unsigned int i) const
Definition: elem.h:2507
TensorTools::MakeReal< Output >::type RealType
For ease of communication, we allow users to translate their own value types to a more easily computa...
std::unordered_map< dof_id_type, typename FFunctor::ValuePushType > new_ids_to_save
static const bool value
Definition: xdr_io.C:54
virtual const Elem & elem_ref(const dof_id_type i) const
Definition: mesh_base.h:639
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:140
TensorTools::MakeBaseNumber< Output >::type DofValueType
The DofObject defines an abstract base class for objects that have degrees of freedom associated with...
Definition: dof_object.h:54
std::pair< const Node *, std::tuple< const Elem *, unsigned short, var_set > > node_projection
std::set< unsigned int > var_set
void parallel_reduce(const Range &range, Body &body)
Execute the provided reduction operation in parallel on the specified range.
Definition: threads_none.h:101
virtual unsigned int size() const override final
Definition: dense_vector.h:104
std::pair< unsigned int, unsigned int > var_to_vg_and_offset(const unsigned int s, const unsigned int var) const
Definition: dof_object.h:1200
This helper structure is used to determine whether a template class is one of our mathematical struct...
Definition: tensor_tools.h:345
SortAndCopy(SortAndCopy &other, Threads::split)
virtual void set(const numeric_index_type i, const T value)=0
Sets v(i) = value.
std::vector< unsigned int > component_to_var
Defines a dense vector for use in Finite Element-type computations.
Definition: dof_map.h:74
void cholesky_solve(const DenseVector< T2 > &b, DenseVector< T2 > &x)
For symmetric positive definite (SPD) matrices.
The VectorSetAction output functor class can be used with a GenericProjector to set projection values...
StoredRange< std::vector< const Elem * >::const_iterator, const Elem * > interior_range
unsigned int n_vars() const
Definition: system.h:2430
This class forms the foundation from which generic finite elements may be derived.
Definition: fe_abstract.h:99
const TypeToSend< T >::type convert_to_send(const T &in)
unsigned int n_comp_group(const unsigned int s, const unsigned int vg) const
Definition: dof_object.h:1015
processor_id_type processor_id() const
void operator()(const node_range &range)
Output eval_at_node(const FEMContext &c, unsigned int i, unsigned int elem_dim, const Node &n, bool, Real=0.)
FEMFunctionBase is a base class from which users can derive in order to define "function-like" object...
bool active() const
Definition: elem.h:2941
const DofMap & get_dof_map() const
Definition: system.h:2374
processor_id_type processor_id() const
Definition: dof_object.h:905
Output eval_at_node(const FEMContext &c, unsigned int i, unsigned int, const Node &n, bool, const Real time)
static void dofs_on_side(const Elem *const elem, const unsigned int dim, const FEType &fe_t, unsigned int s, std::vector< unsigned int > &di, const bool add_p_level=true)
Fills the vector di with the local degree of freedom indices associated with side s of element elem A...
Definition: fe_interface.C:597
The FEMFunctionWrapper input functor class can be used with a GenericProjector to read values from an...
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
const Point & point(const unsigned int i) const
Definition: elem.h:2453
void join(const SortAndCopy &other)
virtual numeric_index_type last_local_index() const =0
void ErrorVector unsigned int
Definition: adjoints_ex3.C:360
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition: int_range.h:117
virtual bool is_child_on_edge(const unsigned int c, const unsigned int e) const
Definition: elem.C:2306
const Elem & get(const ElemType type_in)
This class forms the foundation from which generic finite elements may be derived.
boostcopy::enable_if_c< ScalarTraits< T >::value &&ScalarTraits< T2 >::value, typename CompareTypes< T, T2 >::supertype >::type inner_product(const T &a, const T2 &b)
Definition: tensor_tools.h:51
std::vector< FEContinuity > conts
NodesToElemMap nodes_to_elem_ourcopy
nodes_to_elem is either a shallow copy of a map passed in to the constructor, or points to nodes_to_e...
GenericProjector(const GenericProjector &in)
const Elem * child_ptr(unsigned int i) const
Definition: elem.h:3163
This class defines a tensor in LIBMESH_DIM dimensional Real or Complex space.
uint8_t dof_id_type
Definition: id_types.h:67
Output eval_at_point(const FEMContext &c, unsigned int i, const Point &p, Real, bool skip_context_check)
const std::vector< std::vector< OutputShape > > & get_phi() const
Definition: fe_base.h:207
const FEType & type() const
Definition: variable.h:144
spin_mutex spin_mtx
A convenient spin mutex object which can be used for obtaining locks.
Definition: threads.C:30
FEFieldType
defines an enum for finite element field types - i.e.