libMesh
generic_projector.h
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2019 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 // Local 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/libmesh_logging.h"
34 #include "libmesh/mesh_tools.h"
35 #include "libmesh/numeric_vector.h"
36 #include "libmesh/parallel_sync.h"
37 #include "libmesh/quadrature.h"
38 #include "libmesh/system.h"
39 #include "libmesh/threads.h"
40 
41 namespace libMesh
42 {
43 
50 template <typename T>
51 struct TypeToSend {
52  typedef T type;
53 };
54 
55 template <typename T>
56 const typename TypeToSend<T>::type convert_to_send (const T& in)
57 { return in; }
58 
59 template <typename SendT, typename T>
60 void convert_from_receive (SendT & received, T & converted)
61 { converted = received; }
62 
73 template <typename FFunctor, typename GFunctor,
74  typename FValue, typename ProjectionAction>
76 {
77 private:
78  const System & system;
79 
80  // For TBB compatibility and thread safety we'll copy these in
81  // operator()
82  const FFunctor & master_f;
83  const GFunctor * master_g; // Needed for C1 type elements only
85  const ProjectionAction & master_action;
86  const std::vector<unsigned int> & variables;
87  std::unordered_map<dof_id_type, std::vector<dof_id_type>> * nodes_to_elem;
89 
90 public:
91  GenericProjector (const System & system_in,
92  const FFunctor & f_in,
93  const GFunctor * g_in,
94  const ProjectionAction & act_in,
95  const std::vector<unsigned int> & variables_in,
96  std::unordered_map<dof_id_type, std::vector<dof_id_type>> *
97  nodes_to_elem_in = nullptr) :
98  system(system_in),
99  master_f(f_in),
100  master_g(g_in),
101  g_was_copied(false),
102  map_was_created(!nodes_to_elem_in),
103  master_action(act_in),
104  variables(variables_in),
105  nodes_to_elem(nodes_to_elem_in)
106  {
107  if (map_was_created) // past tense misnomer here
108  {
109  nodes_to_elem = new
110  std::unordered_map<dof_id_type, std::vector<dof_id_type>>;
112  }
113  }
114 
116  system(in.system),
117  master_f(in.master_f),
118  master_g(in.master_g ? new GFunctor(*in.master_g) : nullptr),
121  variables(in.variables),
123  {}
124 
126  {
127  if (g_was_copied)
128  delete master_g;
129  if (map_was_created)
130  delete nodes_to_elem;
131  }
132 
133  // The single "project" function handles all the threaded projection
134  // calculations and intervening MPI communications
135  void project(const ConstElemRange & range);
136 
137  // We generally need to hang on to every value we've calculated
138  // until we're all done, because later projection calculations
139  // depend on boundary data from earlier calculations.
140  std::unordered_map<dof_id_type, FValue> ids_to_save;
141 
142  // This needs to be sorted so we can get sorted dof indices cheaply
143  // later
144  typedef std::set<unsigned int> var_set;
145 
146  // We now need to divide our threadable work into stages, to allow
147  // for the potential necessity of MPI communication in between them.
148  struct SubFunctor {
150 
152 
153  // While we're working on nodes, also figure out which ghosted
154  // nodes will have data we might need to send to their owners
155  // instead of being acted on by ourselves.
156  //
157  // We keep track of which dof ids we might need to send, and what
158  // values those ids should get (along with a pprocessor_id to
159  // leave invalid in case *we* can't compute those values either).
160  std::unordered_map<dof_id_type, std::pair<FValue, processor_id_type>> new_ids_to_push;
161 
162  // We'll hang on to any new ids to save on a per-thread basis
163  // because we won't need them until subsequent phases
164  std::unordered_map<dof_id_type, FValue> new_ids_to_save;
165 
166  // Helper function for filling new_ids_to_push
167  void find_dofs_to_send (const Node & node,
168  const Elem & elem,
169  unsigned short node_num,
170  const var_set & vars);
171 
172  // When we have new data to act on, we may also need to save it
173  // and get ready to push it.
174  void insert_id(dof_id_type id,
175  const FValue & val,
176  processor_id_type pid);
177 
178  void insert_ids(const std::vector<dof_id_type> & ids,
179  const std::vector<FValue> & vals,
180  processor_id_type pid);
181 
182  // Our subclasses will need to merge saved ids
183  void join (const SubFunctor & other);
184 
185  protected:
186  // For thread safety and efficiency we'll want thread-local copies
187  // of various functors
188  ProjectionAction action;
189  FFunctor f;
190 
191  // Each stage of the work we do will require us to prepare
192  // FEMContext objects to assist in the work.
194 
195  // These get used often enough to cache them
196  std::vector<FEContinuity> conts;
197 
198  const System & system;
199  };
200 
201 
202  typedef std::pair<const Node *, std::tuple<const Elem *, unsigned short, var_set>>
204 
207 
208 
209  struct SubProjector : public SubFunctor {
211 
212  using SubFunctor::action;
213  using SubFunctor::f;
214  using SubFunctor::system;
215  using SubFunctor::context;
216  using SubFunctor::conts;
217  using SubFunctor::insert_id;
219 
220  protected:
221  // Projections of C1 elements require a gradient as well
222  std::unique_ptr<GFunctor> g;
223 
225  (const std::vector<dof_id_type> & dof_indices_var,
226  const std::vector<unsigned int> & involved_dofs,
227  unsigned int var_component,
228  const Node * node,
229  const FEBase & fe);
230  };
231 
232 
233  // First we'll copy DoFs where we can, while sorting remaining DoFs
234  // for interpolation and projection later.
235  struct SortAndCopy : public SubFunctor {
237 
239 
240  using SubFunctor::action;
241  using SubFunctor::f;
242  using SubFunctor::system;
243  using SubFunctor::context;
245 
246  void operator() (const ConstElemRange & range);
247 
248  // We need to merge saved multimaps when working from multiple
249  // threads
250  void join (const SortAndCopy & other);
251 
252  // For vertices, edges and sides, we need to know what element,
253  // what local vertex/edge/side number, and what set of variables
254  // to evaluate.
255  typedef std::unordered_multimap<const Node *, std::tuple<const Elem *, unsigned short, var_set>> ves_multimap;
256 
258  std::vector<const Elem *> interiors;
259  };
260 
261 
262  struct ProjectVertices : public SubProjector {
264 
266 
267  using SubProjector::action;
268  using SubProjector::f;
269  using SubProjector::g;
270  using SubProjector::system;
271  using SubProjector::context;
272  using SubProjector::conts;
273 
277 
278  void operator() (const node_range & range);
279  };
280 
281 
282  struct ProjectEdges : public SubProjector {
284 
286 
287  using SubProjector::action;
288  using SubProjector::f;
289  using SubProjector::g;
290  using SubProjector::system;
291  using SubProjector::context;
292  using SubProjector::conts;
293 
297 
298  void operator() (const node_range & range);
299  };
300 
301  struct ProjectSides : public SubProjector {
303 
305 
306  using SubProjector::action;
307  using SubProjector::f;
308  using SubProjector::g;
309  using SubProjector::system;
310  using SubProjector::context;
311  using SubProjector::conts;
312 
316 
317  void operator() (const node_range & range);
318  };
319 
320 
323 
324  struct ProjectInteriors : public SubProjector {
326 
328 
329  using SubProjector::action;
330  using SubProjector::f;
331  using SubProjector::g;
332  using SubProjector::system;
333  using SubProjector::context;
334  using SubProjector::conts;
335 
339 
340  void operator() (const interior_range & range);
341  };
342 
344  (std::unordered_map<dof_id_type, std::pair<FValue, processor_id_type>> & ids_to_push,
345  ProjectionAction & action) const;
346 };
347 
348 
358 template <typename Val>
360 {
361 private:
363 
364 public:
366  target_vector(target_vec) {}
367 
369  Val val)
370  {
371  // Lock the new vector since it is shared among threads.
372  {
373  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
374  target_vector.set(id, val);
375  }
376  }
377 
378 
379  void insert(const std::vector<dof_id_type> & dof_indices,
380  const DenseVector<Val> & Ue)
381  {
382  const numeric_index_type
385 
386  unsigned int size = Ue.size();
387 
388  libmesh_assert_equal_to(size, dof_indices.size());
389 
390  // Lock the new vector since it is shared among threads.
391  {
392  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
393 
394  for (unsigned int i = 0; i != size; ++i)
395  if ((dof_indices[i] >= first) && (dof_indices[i] < last))
396  target_vector.set(dof_indices[i], Ue(i));
397  }
398  }
399 
400 };
401 
402 
411 template <typename Output>
413 {
414 public:
415  FEMFunctionWrapper(const FEMFunctionBase<Output> & f) : _f(f.clone()) {}
416 
418  _f(fw._f->clone()) {}
419 
420  void init_context (FEMContext & c) { _f->init_context(c); }
421 
422  Output eval_at_node (const FEMContext & c,
423  unsigned int i,
424  unsigned int /*elem_dim*/,
425  const Node & n,
426  bool /*extra_hanging_dofs*/,
427  const Real time)
428  { return _f->component(c, i, n, time); }
429 
430  Output eval_at_point (const FEMContext & c,
431  unsigned int i,
432  const Point & n,
433  const Real time)
434  { return _f->component(c, i, n, time); }
435 
436  bool is_grid_projection() { return false; }
437 
438  void eval_old_dofs (const Elem &,
439  unsigned int,
440  unsigned int,
441  std::vector<dof_id_type> &,
442  std::vector<Output> &)
443  { libmesh_error(); }
444 
445  void eval_old_dofs (const Elem &,
446  const FEType &,
447  unsigned int,
448  unsigned int,
449  std::vector<dof_id_type> &,
450  std::vector<Output> &)
451  { libmesh_error(); }
452 
453 private:
454  std::unique_ptr<FEMFunctionBase<Output>> _f;
455 };
456 
457 
468 #ifdef LIBMESH_ENABLE_AMR
469 template <typename Output,
470  void (FEMContext::*point_output) (unsigned int,
471  const Point &,
472  Output &,
473  const Real) const>
475 {
476 public:
478  last_elem(nullptr),
479  sys(sys_in),
480  old_context(sys_in)
481  {
482  // We'll be queried for components but we'll typically be looking
483  // up data by variables, and those indices don't always match
484  for (auto v : IntRange<unsigned int>(0, sys.n_vars()))
485  {
486  const unsigned int vcomp = sys.variable_scalar_number(v,0);
487  if (vcomp >= component_to_var.size())
488  component_to_var.resize(vcomp+1, static_cast<unsigned int>(-1));
489  component_to_var[vcomp] = v;
490  }
491  }
492 
494  last_elem(nullptr),
495  sys(in.sys),
496  old_context(sys),
497  component_to_var(in.component_to_var)
498  {
499  }
500 
501  static void get_shape_outputs(FEBase & fe);
502 
503  // Integrating on new mesh elements, we won't yet have an up to date
504  // current_local_solution.
506  {
508 
509  // Loop over variables, to prerequest
510  for (unsigned int var=0; var!=sys.n_vars(); ++var)
511  {
512  FEBase * fe = nullptr;
513  const std::set<unsigned char> & elem_dims =
514  old_context.elem_dimensions();
515 
516  for (const auto & dim : elem_dims)
517  {
518  old_context.get_element_fe(var, fe, dim);
519  get_shape_outputs(*fe);
520  }
521  }
522  }
523 
524  bool is_grid_projection() { return true; }
525 
526 protected:
527  void check_old_context (const FEMContext & c)
528  {
529  LOG_SCOPE ("check_old_context(c)", "OldSolutionBase");
530  const Elem & elem = c.get_elem();
531  if (last_elem != &elem)
532  {
533  if (elem.refinement_flag() == Elem::JUST_REFINED)
534  {
535  old_context.pre_fe_reinit(sys, elem.parent());
536  }
537  else if (elem.refinement_flag() == Elem::JUST_COARSENED)
538  {
539  libmesh_error();
540  }
541  else
542  {
543  if (!elem.old_dof_object)
544  {
545  libmesh_error();
546  }
547 
548  old_context.pre_fe_reinit(sys, &elem);
549  }
550 
551  last_elem = &elem;
552  }
553  else
554  {
555  libmesh_assert(old_context.has_elem());
556  }
557  }
558 
559 
560  bool check_old_context (const FEMContext & c, const Point & p)
561  {
562  LOG_SCOPE ("check_old_context(c,p)", "OldSolutionBase");
563  const Elem & elem = c.get_elem();
564  if (last_elem != &elem)
565  {
566  if (elem.refinement_flag() == Elem::JUST_REFINED)
567  {
568  old_context.pre_fe_reinit(sys, elem.parent());
569  }
570  else if (elem.refinement_flag() == Elem::JUST_COARSENED)
571  {
572  // Find the child with this point. Use out_of_elem_tol
573  // (in physical space, which may correspond to a large
574  // tolerance in master space!) to allow for out-of-element
575  // finite differencing of mixed gradient terms. Pray we
576  // have no quadrature locations which are within 1e-5 of
577  // the element subdivision boundary but are not exactly on
578  // that boundary.
579  const Real master_tol = out_of_elem_tol / elem.hmax() * 2;
580 
581  for (auto & child : elem.child_ref_range())
582  if (child.close_to_point(p, master_tol))
583  {
584  old_context.pre_fe_reinit(sys, &child);
585  break;
586  }
587 
588  libmesh_assert
589  (old_context.get_elem().close_to_point(p, master_tol));
590  }
591  else
592  {
593  if (!elem.old_dof_object)
594  return false;
595 
596  old_context.pre_fe_reinit(sys, &elem);
597  }
598 
599  last_elem = &elem;
600  }
601  else
602  {
603  libmesh_assert(old_context.has_elem());
604 
605  const Real master_tol = out_of_elem_tol / elem.hmax() * 2;
606 
607  if (!old_context.get_elem().close_to_point(p, master_tol))
608  {
609  libmesh_assert_equal_to
611 
612  for (auto & child : elem.child_ref_range())
613  if (child.close_to_point(p, master_tol))
614  {
615  old_context.pre_fe_reinit(sys, &child);
616  break;
617  }
618 
619  libmesh_assert
620  (old_context.get_elem().close_to_point(p, master_tol));
621  }
622  }
623 
624  return true;
625  }
626 
627 protected:
628  const Elem * last_elem;
629  const System & sys;
631  std::vector<unsigned int> component_to_var;
632 
633  static const Real out_of_elem_tol;
634 };
635 
636 
645 template <typename Output,
646  void (FEMContext::*point_output) (unsigned int,
647  const Point &,
648  Output &,
649  const Real) const>
651 {
652 public:
654  const NumericVector<Number> & old_sol) :
655  OldSolutionBase<Output, point_output>(sys_in),
656  old_solution(old_sol)
657  {
658  this->old_context.set_algebraic_type(FEMContext::OLD);
659  this->old_context.set_custom_solution(&old_solution);
660  }
661 
663  OldSolutionBase<Output, point_output>(in.sys),
664  old_solution(in.old_solution)
665  {
666  this->old_context.set_algebraic_type(FEMContext::OLD);
667  this->old_context.set_custom_solution(&old_solution);
668  }
669 
670 
671  Output eval_at_node (const FEMContext & c,
672  unsigned int i,
673  unsigned int elem_dim,
674  const Node & n,
675  bool /* extra_hanging_dofs */,
676  Real /* time */ =0.);
677 
678  Output eval_at_point(const FEMContext & c,
679  unsigned int i,
680  const Point & p,
681  Real /* time */ =0.)
682  {
683  LOG_SCOPE ("eval_at_point()", "OldSolutionValue");
684 
685  if (!this->check_old_context(c, p))
686  return 0;
687 
688  // Handle offset from non-scalar components in previous variables
689  libmesh_assert_less(i, this->component_to_var.size());
690  unsigned int var = this->component_to_var[i];
691 
692  Output n;
693  (this->old_context.*point_output)(var, p, n, this->out_of_elem_tol);
694  return n;
695  }
696 
697  void eval_old_dofs (const Elem & elem,
698  unsigned int node_num,
699  unsigned int var_num,
700  std::vector<dof_id_type> & indices,
701  std::vector<Output> & values)
702  {
703  LOG_SCOPE ("eval_old_dofs(node)", "OldSolutionValue");
704 
705  this->sys.get_dof_map().dof_indices(elem, node_num, indices, var_num);
706 
707  std::vector<dof_id_type> old_indices;
708 
709  this->sys.get_dof_map().old_dof_indices(elem, node_num, old_indices, var_num);
710 
711  libmesh_assert_equal_to (old_indices.size(), indices.size());
712 
713  values.resize(old_indices.size());
714 
715  old_solution.get(old_indices, values);
716  }
717 
718 
719  void eval_old_dofs (const Elem & elem,
720  const FEType & fe_type,
721  unsigned int sys_num,
722  unsigned int var_num,
723  std::vector<dof_id_type> & indices,
724  std::vector<Output> & values)
725  {
726  LOG_SCOPE ("eval_old_dofs(elem)", "OldSolutionValue");
727 
728  // We're only to be asked for old dofs on elements that can copy
729  // them through DO_NOTHING or through refinement.
730  const Elem & old_elem =
731  (elem.refinement_flag() == Elem::JUST_REFINED) ?
732  *elem.parent() : elem;
733 
734  // If there are any element-based DOF numbers, get them
735  const unsigned int nc = FEInterface::n_dofs_per_elem(elem.dim(),
736  fe_type,
737  elem.type());
738 
739  std::vector<dof_id_type> old_dof_indices(nc);
740  indices.resize(nc);
741 
742  // We should never have fewer dofs than necessary on an
743  // element unless we're getting indices on a parent element,
744  // and we should never need those indices
745  if (nc != 0)
746  {
747  libmesh_assert(old_elem.old_dof_object);
748 
749  const std::pair<unsigned int, unsigned int>
750  vg_and_offset = elem.var_to_vg_and_offset(sys_num,var_num);
751  const unsigned int vg = vg_and_offset.first;
752  const unsigned int vig = vg_and_offset.second;
753 
754  const unsigned int n_comp = elem.n_comp_group(sys_num,vg);
755  libmesh_assert_greater(elem.n_systems(), sys_num);
756  libmesh_assert_greater_equal(n_comp, nc);
757 
758  for (unsigned int i=0; i<nc; i++)
759  {
760  const dof_id_type d_old =
761  old_elem.old_dof_object->dof_number(sys_num, vg, vig, i, n_comp);
762  const dof_id_type d_new =
763  elem.dof_number(sys_num, vg, vig, i, n_comp);
764  libmesh_assert_not_equal_to (d_old, DofObject::invalid_id);
765  libmesh_assert_not_equal_to (d_new, DofObject::invalid_id);
766 
767  old_dof_indices[i] = d_old;
768  indices[i] = d_new;
769  }
770  }
771 
772  values.resize(nc);
773 
774  old_solution.get(old_dof_indices, values);
775  }
776 
777 private:
779 };
780 
781 
782 template<>
783 inline void
785 {
786  fe.get_phi();
787 }
788 
789 
790 template<>
791 inline void
793 {
794  fe.get_dphi();
795 }
796 
797 
798 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
799 template<>
800 inline void
802 {
803  fe.get_phi();
804 }
805 
806 
807 template<>
808 inline void
810 {
811  fe.get_dphi();
812 }
813 #endif // LIBMESH_USE_COMPLEX_NUMBERS
814 
815 
816 template<>
817 inline
818 Number
821  unsigned int i,
822  unsigned int /* elem_dim */,
823  const Node & n,
824  bool extra_hanging_dofs,
825  Real /* time */)
826 {
827  LOG_SCOPE ("Number eval_at_node()", "OldSolutionValue");
828 
829  // This should only be called on vertices
830  libmesh_assert_less(c.get_elem().get_node_index(&n),
831  c.get_elem().n_vertices());
832 
833  // Handle offset from non-scalar components in previous variables
834  libmesh_assert_less(i, this->component_to_var.size());
835  unsigned int var = this->component_to_var[i];
836 
837  // Optimize for the common case, where this node was part of the
838  // old solution.
839  //
840  // Be sure to handle cases where the variable wasn't defined on
841  // this node (due to changing subdomain support) or where the
842  // variable has no components on this node (due to Elem order
843  // exceeding FE order) or where the old_dof_object dofs might
844  // correspond to non-vertex dofs (due to extra_hanging_dofs and
845  // refinement)
846 
848 
849  if (n.old_dof_object &&
850  (!extra_hanging_dofs ||
851  flag == Elem::JUST_COARSENED ||
852  flag == Elem::DO_NOTHING) &&
853  n.old_dof_object->n_vars(sys.number()) &&
854  n.old_dof_object->n_comp(sys.number(), var))
855  {
856  const dof_id_type old_id =
857  n.old_dof_object->dof_number(sys.number(), var, 0);
858  return old_solution(old_id);
859  }
860 
861  return this->eval_at_point(c, i, n, 0);
862 }
863 
864 
865 
866 template<>
867 inline
868 Gradient
871  unsigned int i,
872  unsigned int elem_dim,
873  const Node & n,
874  bool extra_hanging_dofs,
875  Real /* time */)
876 {
877  LOG_SCOPE ("Gradient eval_at_node()", "OldSolutionValue");
878 
879  // This should only be called on vertices
880  libmesh_assert_less(c.get_elem().get_node_index(&n),
881  c.get_elem().n_vertices());
882 
883  // Handle offset from non-scalar components in previous variables
884  libmesh_assert_less(i, this->component_to_var.size());
885  unsigned int var = this->component_to_var[i];
886 
887  // Optimize for the common case, where this node was part of the
888  // old solution.
889  //
890  // Be sure to handle cases where the variable wasn't defined on
891  // this node (due to changing subdomain support) or where the
892  // variable has no components on this node (due to Elem order
893  // exceeding FE order) or where the old_dof_object dofs might
894  // correspond to non-vertex dofs (due to extra_hanging_dofs and
895  // refinement)
896 
898 
899  if (n.old_dof_object &&
900  (!extra_hanging_dofs ||
901  flag == Elem::JUST_COARSENED ||
902  flag == Elem::DO_NOTHING) &&
903  n.old_dof_object->n_vars(sys.number()) &&
904  n.old_dof_object->n_comp(sys.number(), var))
905  {
906  Gradient g;
907  for (unsigned int d = 0; d != elem_dim; ++d)
908  {
909  const dof_id_type old_id =
910  n.old_dof_object->dof_number(sys.number(), var, d+1);
911  g(d) = old_solution(old_id);
912  }
913  return g;
914  }
915 
916  return this->eval_at_point(c, i, n, 0);
917 }
918 
919 
920 
921 
922 
923 template <>
925 
926 template <>
928 
929 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
930 template <>
932 
933 template <>
935 #endif // LIBMESH_USE_COMPLEX_NUMBERS
936 
937 #endif // LIBMESH_ENABLE_AMR
938 
943 template <typename FFunctor, typename GFunctor,
944  typename FValue, typename ProjectionAction>
946  (const ConstElemRange & range)
947 {
948  LOG_SCOPE ("project", "GenericProjector");
949 
950  // Unless we split sort and copy into two passes we can't know for
951  // sure ahead of time whether we need to save the copied ids
952  done_saving_ids = false;
953 
954  SortAndCopy sort_work(*this);
955  Threads::parallel_reduce (range, sort_work);
956  ProjectionAction action(master_action);
957 
958  // Keep track of dof ids and values to send to other ranks
959  std::unordered_map<dof_id_type, std::pair<FValue, processor_id_type>> ids_to_push;
960 
961  ids_to_push.insert(sort_work.new_ids_to_push.begin(),
962  sort_work.new_ids_to_push.end());
963  ids_to_save.insert(sort_work.new_ids_to_save.begin(),
964  sort_work.new_ids_to_save.end());
965 
966  std::vector<node_projection> vertices(sort_work.vertices.begin(),
967  sort_work.vertices.end());
968 
969  done_saving_ids = sort_work.edges.empty() &&
970  sort_work.sides.empty() && sort_work.interiors.empty();
971  system.comm().max(done_saving_ids);
972 
973  {
974  ProjectVertices project_vertices(*this);
975  Threads::parallel_reduce (node_range(&vertices), project_vertices);
976  ids_to_push.insert(project_vertices.new_ids_to_push.begin(),
977  project_vertices.new_ids_to_push.end());
978  ids_to_save.insert(project_vertices.new_ids_to_save.begin(),
979  project_vertices.new_ids_to_save.end());
980  }
981 
982  done_saving_ids = sort_work.sides.empty() && sort_work.interiors.empty();
983  system.comm().max(done_saving_ids);
984 
985  this->send_and_insert_dof_values(ids_to_push, action);
986 
987  {
988  std::vector<node_projection> edges(sort_work.edges.begin(), sort_work.edges.end());
989  ProjectEdges project_edges(*this);
990  Threads::parallel_reduce (node_range(&edges), project_edges);
991  ids_to_push.insert(project_edges.new_ids_to_push.begin(),
992  project_edges.new_ids_to_push.end());
993  ids_to_save.insert(project_edges.new_ids_to_save.begin(),
994  project_edges.new_ids_to_save.end());
995  }
996 
997  done_saving_ids = sort_work.interiors.empty();
998  system.comm().max(done_saving_ids);
999 
1000  this->send_and_insert_dof_values(ids_to_push, action);
1001 
1002  {
1003  std::vector<node_projection> sides(sort_work.sides.begin(), sort_work.sides.end());
1004  ProjectSides project_sides(*this);
1005  Threads::parallel_reduce (node_range(&sides), project_sides);
1006  ids_to_push.insert(project_sides.new_ids_to_push.begin(),
1007  project_sides.new_ids_to_push.end());
1008  ids_to_save.insert(project_sides.new_ids_to_save.begin(),
1009  project_sides.new_ids_to_save.end());
1010  }
1011 
1012  done_saving_ids = true;
1013  this->send_and_insert_dof_values(ids_to_push, action);
1014 
1015  // No ids to save or push this time, but we still use a reduce since
1016  // nominally ProjectInteriors still has non-const operator()
1017  ProjectInteriors project_interiors(*this);
1019  project_interiors);
1020 }
1021 
1022 
1023 template <typename FFunctor, typename GFunctor,
1024  typename FValue, typename ProjectionAction>
1027  projector(p), action(p.master_action), f(p.master_f),
1028  context(p.system), conts(p.system.n_vars()), system(p.system)
1029 {
1030  // Loop over all the variables we've been requested to project, to
1031  // pre-request
1032  for (const auto & var : this->projector.variables)
1033  {
1034  // FIXME: Need to generalize this to vector-valued elements. [PB]
1035  FEBase * fe = nullptr;
1036  FEBase * side_fe = nullptr;
1037  FEBase * edge_fe = nullptr;
1038 
1039  const std::set<unsigned char> & elem_dims =
1040  context.elem_dimensions();
1041 
1042  for (const auto & dim : elem_dims)
1043  {
1044  context.get_element_fe( var, fe, dim );
1045  if (fe->get_fe_type().family == SCALAR)
1046  continue;
1047  if (dim > 1)
1048  context.get_side_fe( var, side_fe, dim );
1049  if (dim > 2)
1050  context.get_edge_fe( var, edge_fe );
1051 
1052  fe->get_JxW();
1053  fe->get_xyz();
1054  fe->get_JxW();
1055 
1056  fe->get_phi();
1057  if (dim > 1)
1058  {
1059  side_fe->get_JxW();
1060  side_fe->get_xyz();
1061  side_fe->get_phi();
1062  }
1063  if (dim > 2)
1064  {
1065  edge_fe->get_JxW();
1066  edge_fe->get_xyz();
1067  edge_fe->get_phi();
1068  }
1069 
1070  const FEContinuity cont = fe->get_continuity();
1071  this->conts[var] = cont;
1072  if (cont == C_ONE)
1073  {
1074  fe->get_dphi();
1075  if (dim > 1)
1076  side_fe->get_dphi();
1077  if (dim > 2)
1078  edge_fe->get_dphi();
1079  }
1080  }
1081  }
1082 
1083  // Now initialize any user requested shape functions, xyz vals, etc
1084  f.init_context(context);
1085 }
1086 
1087 
1088 template <typename FFunctor, typename GFunctor,
1089  typename FValue, typename ProjectionAction>
1092  SubFunctor(p)
1093 {
1094  if (p.master_g)
1095  g.reset(new GFunctor(*p.master_g));
1096 
1097 #ifndef NDEBUG
1098  // Our C1 elements need gradient information
1099  for (const auto & var : this->projector.variables)
1100  if (this->conts[var] == C_ONE)
1101  libmesh_assert(g);
1102 #endif
1103 
1104  if (g)
1105  g->init_context(context);
1106 }
1107 
1108 
1109 template <typename FFunctor, typename GFunctor,
1110  typename FValue, typename ProjectionAction>
1112  (dof_id_type id, const FValue & val, processor_id_type pid)
1113 {
1114  auto iter = new_ids_to_push.find(id);
1115  if (iter == new_ids_to_push.end())
1116  action.insert(id, val);
1117  else
1118  {
1119  libmesh_assert(pid != DofObject::invalid_processor_id);
1120  iter->second = std::make_pair(val, pid);
1121  }
1122  if (!this->projector.done_saving_ids)
1123  {
1124  libmesh_assert(!new_ids_to_save.count(id));
1125  new_ids_to_save[id] = val;
1126  }
1127 }
1128 
1129 
1130 template <typename FFunctor, typename GFunctor,
1131  typename FValue, typename ProjectionAction>
1133  (const std::vector<dof_id_type> & ids, const std::vector<FValue> & vals, processor_id_type pid)
1134 {
1135  libmesh_assert_equal_to(ids.size(), vals.size());
1136  for (auto i : index_range(ids))
1137  {
1138  const dof_id_type id = ids[i];
1139  const FValue & val = vals[i];
1140 
1141  auto iter = new_ids_to_push.find(id);
1142  if (iter == new_ids_to_push.end())
1143  action.insert(id, val);
1144  else
1145  {
1146  libmesh_assert(pid != DofObject::invalid_processor_id);
1147  iter->second = std::make_pair(val, pid);
1148  }
1149  if (!this->projector.done_saving_ids)
1150  {
1151  libmesh_assert(!new_ids_to_save.count(id));
1152  new_ids_to_save[id] = val;
1153  }
1154  }
1155 }
1156 
1157 
1158 template <typename FFunctor, typename GFunctor,
1159  typename FValue, typename ProjectionAction>
1162 {
1163  new_ids_to_push.insert(other.new_ids_to_push.begin(), other.new_ids_to_push.end());
1164  new_ids_to_save.insert(other.new_ids_to_save.begin(), other.new_ids_to_save.end());
1165 }
1166 
1167 
1168 template <typename FFunctor, typename GFunctor,
1169  typename FValue, typename ProjectionAction>
1171  (const ConstElemRange & range)
1172 {
1173  // Look at all the elements in the range. Directly copy values from
1174  // unchanged elements. For other elements, determine sets of
1175  // vertices, edge nodes, and side nodes to project.
1176  //
1177  // As per our other weird nomenclature, "sides" means faces in 3D
1178  // and edges in 2D, and "edges" gets skipped in 2D
1179  //
1180  // This gets tricky in the case of subdomain-restricted
1181  // variables, for which we might need to do the same projection
1182  // from different elements when evaluating different variables.
1183  // We'll keep track of which variables can be projected from which
1184  // elements.
1185 
1186  // If we copy DoFs on a node, keep track of it so we don't bother
1187  // with any redundant interpolations or projections later.
1188  //
1189  // We still need to keep track of *which* variables get copied,
1190  // since we may be able to copy elements which lack some
1191  // subdomain-restricted variables.
1192  //
1193  // For extra_hanging_dofs FE types, we'll keep track of which
1194  // variables were copied from vertices, and which from edges/sides,
1195  // because those types need copies made from *both* on hanging
1196  // nodes.
1197  std::unordered_map<const Node *, std::pair<var_set, var_set>> copied_nodes;
1198 
1199  const unsigned int sys_num = system.number();
1200 
1201  // At hanging nodes for variables with extra hanging dofs we'll need
1202  // to do projections *separately* from vertex elements and side/edge
1203  // elements.
1204  std::vector<unsigned short> extra_hanging_dofs;
1205  bool all_extra_hanging_dofs = true;
1206  for (auto v_num : this->projector.variables)
1207  {
1208  if (extra_hanging_dofs.size() <= v_num)
1209  extra_hanging_dofs.resize(v_num+1, false);
1210  extra_hanging_dofs[v_num] =
1211  FEInterface::extra_hanging_dofs(system.variable(v_num).type());
1212 
1213  if (!extra_hanging_dofs[v_num])
1214  all_extra_hanging_dofs = false;
1215  }
1216 
1217  for (const auto & elem : range)
1218  {
1219  // If we're doing AMR, we might be able to copy more DoFs than
1220  // we interpolate or project.
1221  bool copy_this_elem = false;
1222 
1223 #ifdef LIBMESH_ENABLE_AMR
1224  // If we're projecting from an old grid
1225  if (f.is_grid_projection())
1226  {
1227  // If this element doesn't have an old_dof_object, but it
1228  // wasn't just refined or just coarsened into activity, then
1229  // it must be newly added, so the user is responsible for
1230  // setting the new dofs on it during a grid projection.
1231  if (!elem->old_dof_object &&
1232  elem->refinement_flag() != Elem::JUST_REFINED &&
1233  elem->refinement_flag() != Elem::JUST_COARSENED)
1234  continue;
1235 
1236  // If this is an unchanged element, just copy everything
1237  if ((elem->refinement_flag() != Elem::JUST_REFINED &&
1238  elem->refinement_flag() != Elem::JUST_COARSENED &&
1239  elem->p_refinement_flag() != Elem::JUST_REFINED &&
1240  elem->p_refinement_flag() != Elem::JUST_COARSENED))
1241  copy_this_elem = true;
1242  else
1243  {
1244  bool reinitted = false;
1245 
1246  // If this element has a low order monomial which has
1247  // merely been h refined, copy it.
1248  for (auto v_num : this->projector.variables)
1249  {
1250  const Variable & var = system.variable(v_num);
1251  if (!var.active_on_subdomain(elem->subdomain_id()))
1252  continue;
1253  FEType fe_type = var.type();
1254 
1255  if (fe_type.family == MONOMIAL &&
1256  fe_type.order == CONSTANT &&
1257  elem->p_level() == 0 &&
1258  elem->refinement_flag() != Elem::JUST_COARSENED &&
1259  elem->p_refinement_flag() != Elem::JUST_COARSENED)
1260  {
1261  if (!reinitted)
1262  {
1263  reinitted = true;
1264  context.pre_fe_reinit(system, elem);
1265  }
1266 
1267  std::vector<FValue> Ue(1);
1268  std::vector<dof_id_type> elem_dof_ids(1);
1269 
1270  f.eval_old_dofs(*elem, fe_type, sys_num, v_num,
1271  elem_dof_ids, Ue);
1272 
1273  action.insert(elem_dof_ids[0], Ue[0]);
1274  }
1275  }
1276  }
1277  }
1278 #endif // LIBMESH_ENABLE_AMR
1279 
1280  const int dim = elem->dim();
1281 
1282  const unsigned int n_vertices = elem->n_vertices();
1283  const unsigned int n_edges = elem->n_edges();
1284  const unsigned int n_nodes = elem->n_nodes();
1285 
1286  // In 1-D we already handle our sides as vertices
1287  const unsigned int n_sides = (dim > 1) * elem->n_sides();
1288 
1289  // What variables are supported on each kind of node on this elem?
1290  var_set vertex_vars, edge_vars, side_vars;
1291 
1292  // If we have non-vertex nodes, the first is an edge node, but
1293  // if we're in 2D we'll call that a side node
1294  const bool has_edge_nodes = (n_nodes > n_vertices && dim > 2);
1295 
1296  // If we have even more nodes, the next is a side node.
1297  const bool has_side_nodes =
1298  (n_nodes > n_vertices + ((dim > 2) * n_edges));
1299 
1300  // We may be out of nodes at this point or we have interior
1301  // nodes which may have DoFs to project too
1302  const bool has_interior_nodes =
1303  (n_nodes > n_vertices + ((dim > 2) * n_edges) + n_sides);
1304 
1305  for (auto v_num : this->projector.variables)
1306  {
1307  const Variable & var = this->projector.system.variable(v_num);
1308  if (!var.active_on_subdomain(elem->subdomain_id()))
1309  continue;
1310  FEType fe_type = var.type();
1311  fe_type.order =
1312  libMesh::Order (fe_type.order + elem->p_level());
1313  const ElemType elem_type = elem->type();
1314 
1315  if (FEInterface::n_dofs_at_node(dim, fe_type, elem_type, 0))
1316  vertex_vars.insert(vertex_vars.end(), v_num);
1317 
1318  // The first non-vertex node is always an edge node if those
1319  // exist. All edge nodes have the same number of DoFs
1320  if (has_edge_nodes)
1321  if (FEInterface::n_dofs_at_node(dim, fe_type, elem_type, n_vertices))
1322  edge_vars.insert(edge_vars.end(), v_num);
1323 
1324  if (has_side_nodes)
1325  {
1326  if (dim != 3)
1327  {
1328  if (FEInterface::n_dofs_at_node(dim, fe_type, elem_type, n_vertices))
1329  side_vars.insert(side_vars.end(), v_num);
1330  }
1331  else
1332  // In 3D, not all face nodes always have the same number of
1333  // DoFs! We'll loop over all sides to be safe.
1334  for (unsigned int n = 0; n != n_nodes; ++n)
1335  if (elem->is_face(n))
1336  if (FEInterface::n_dofs_at_node(dim, fe_type,
1337  elem_type, n))
1338  {
1339  side_vars.insert(side_vars.end(), v_num);
1340  break;
1341  }
1342  }
1343 
1344  if (FEInterface::n_dofs_per_elem(dim, fe_type, elem_type) ||
1345  (has_interior_nodes &&
1346  FEInterface::n_dofs_at_node(dim, fe_type, elem_type, n_nodes-1)))
1347  {
1348  // We may have already just copied constant monomials,
1349  // or we may be about to copy the whole element
1350  if ((f.is_grid_projection() &&
1351  fe_type.family == MONOMIAL &&
1352  fe_type.order == CONSTANT &&
1353  elem->p_level() == 0 &&
1354  elem->refinement_flag() != Elem::JUST_COARSENED &&
1355  elem->p_refinement_flag() != Elem::JUST_COARSENED)
1356  || copy_this_elem
1357  )
1358  continue;
1359 
1360  // We need to project any other variables
1361  if (interiors.empty() || interiors.back() != elem)
1362  interiors.push_back(elem);
1363  }
1364  }
1365 
1366  // We'll use a greedy algorithm in most cases: if another
1367  // element has already claimed some of our DoFs, we'll let it do
1368  // the work.
1369  auto erase_covered_vars = []
1370  (const Node * node, var_set & remaining, const ves_multimap & covered)
1371  {
1372  auto covered_range = covered.equal_range(node);
1373  for (const auto & v_ent : as_range(covered_range))
1374  for (const unsigned int var_covered :
1375  std::get<2>(v_ent.second))
1376  remaining.erase(var_covered);
1377  };
1378 
1379  auto erase_nonhanging_vars = [&extra_hanging_dofs]
1380  (const Node * node, var_set & remaining, const ves_multimap & covered)
1381  {
1382  auto covered_range = covered.equal_range(node);
1383  for (const auto & v_ent : as_range(covered_range))
1384  for (const unsigned int var_covered :
1385  std::get<2>(v_ent.second))
1386  if (!extra_hanging_dofs[var_covered])
1387  remaining.erase(var_covered);
1388  };
1389 
1390  auto erase_copied_vars = [&copied_nodes, &extra_hanging_dofs]
1391  (const Node * node, bool is_vertex, var_set & remaining)
1392  {
1393  auto copying_range = copied_nodes.equal_range(node);
1394  for (const auto & v_ent : as_range(copying_range))
1395  {
1396  for (const unsigned int var_covered :
1397  v_ent.second.first)
1398  if (is_vertex || !extra_hanging_dofs[var_covered])
1399  remaining.erase(var_covered);
1400 
1401  for (const unsigned int var_covered :
1402  v_ent.second.second)
1403  if (!is_vertex || !extra_hanging_dofs[var_covered])
1404  remaining.erase(var_covered);
1405  }
1406  };
1407 
1408  for (unsigned int v=0; v != n_vertices; ++v)
1409  {
1410  const Node * node = elem->node_ptr(v);
1411 
1412  auto remaining_vars = vertex_vars;
1413 
1414  erase_covered_vars(node, remaining_vars, vertices);
1415 
1416  if (remaining_vars.empty())
1417  continue;
1418 
1419  if (!all_extra_hanging_dofs)
1420  {
1421  erase_nonhanging_vars(node, remaining_vars, edges);
1422  if (remaining_vars.empty())
1423  continue;
1424 
1425  erase_nonhanging_vars(node, remaining_vars, sides);
1426  if (remaining_vars.empty())
1427  continue;
1428  }
1429 
1430  erase_copied_vars(node, true, remaining_vars);
1431  if (remaining_vars.empty())
1432  continue;
1433 
1434  if (copy_this_elem)
1435  {
1436  for (auto var : remaining_vars)
1437  {
1438  std::vector<dof_id_type> node_dof_ids;
1439  std::vector<FValue> values;
1440 
1441  f.eval_old_dofs(*elem, v, var, node_dof_ids, values);
1442 
1443  insert_ids(node_dof_ids, values, node->processor_id());
1444  }
1445  copied_nodes[node].first.insert(remaining_vars.begin(),
1446  remaining_vars.end());
1447  this->find_dofs_to_send(*node, *elem, v, remaining_vars);
1448  }
1449  else
1450  vertices.emplace
1451  (node, std::make_tuple(elem, v, std::move(remaining_vars)));
1452  }
1453 
1454  if (has_edge_nodes)
1455  {
1456  for (unsigned int e=0; e != n_edges; ++e)
1457  {
1458  const Node * node = elem->node_ptr(n_vertices+e);
1459 
1460  auto remaining_vars = edge_vars;
1461 
1462  erase_covered_vars(node, remaining_vars, edges);
1463  if (remaining_vars.empty())
1464  continue;
1465 
1466  erase_covered_vars(node, remaining_vars, sides);
1467  if (remaining_vars.empty())
1468  continue;
1469 
1470  if (!all_extra_hanging_dofs)
1471  {
1472  erase_nonhanging_vars(node, remaining_vars, vertices);
1473  if (remaining_vars.empty())
1474  continue;
1475  }
1476 
1477  erase_copied_vars(node, false, remaining_vars);
1478  if (remaining_vars.empty())
1479  continue;
1480 
1481  if (copy_this_elem)
1482  {
1483  for (auto var : remaining_vars)
1484  {
1485  std::vector<dof_id_type> edge_dof_ids;
1486  std::vector<FValue> values;
1487 
1488  f.eval_old_dofs(*elem, n_vertices+e, var, edge_dof_ids, values);
1489 
1490  insert_ids(edge_dof_ids, values, node->processor_id());
1491  }
1492  copied_nodes[node].second.insert(remaining_vars.begin(),
1493  remaining_vars.end());
1494  this->find_dofs_to_send(*node, *elem, n_vertices+e, remaining_vars);
1495  }
1496  else
1497  edges.emplace
1498  (node, std::make_tuple(elem, e, std::move(remaining_vars)));
1499  }
1500  }
1501 
1502  if (has_side_nodes)
1503  {
1504  for (unsigned int side=0; side != n_sides; ++side)
1505  {
1506  const Node * node = nullptr;
1507  unsigned short node_num = n_vertices+(dim>2)*n_edges+side;
1508  if (dim != 3)
1509  node = elem->node_ptr(node_num);
1510  else
1511  {
1512  // In 3D only some sides may have nodes
1513  for (unsigned int n = 0; n != n_nodes; ++n)
1514  {
1515  if (!elem->is_face(n))
1516  continue;
1517 
1518  if (elem->is_node_on_side(n, side))
1519  {
1520  node_num = n;
1521  node = elem->node_ptr(node_num);
1522  break;
1523  }
1524  }
1525  }
1526 
1527  if (!node)
1528  continue;
1529 
1530  auto remaining_vars = side_vars;
1531 
1532  erase_covered_vars(node, remaining_vars, edges);
1533  if (remaining_vars.empty())
1534  continue;
1535 
1536  erase_covered_vars(node, remaining_vars, sides);
1537  if (remaining_vars.empty())
1538  continue;
1539 
1540  if (!all_extra_hanging_dofs)
1541  {
1542  erase_nonhanging_vars(node, remaining_vars, vertices);
1543  if (remaining_vars.empty())
1544  continue;
1545  }
1546 
1547  erase_copied_vars(node, false, remaining_vars);
1548  if (remaining_vars.empty())
1549  continue;
1550 
1551  if (copy_this_elem)
1552  {
1553  for (auto var : remaining_vars)
1554  {
1555  std::vector<dof_id_type> side_dof_ids;
1556  std::vector<FValue> values;
1557 
1558  f.eval_old_dofs(*elem, node_num, var, side_dof_ids, values);
1559 
1560  insert_ids(side_dof_ids, values, node->processor_id());
1561  }
1562  copied_nodes[node].second.insert(remaining_vars.begin(),
1563  remaining_vars.end());
1564  this->find_dofs_to_send(*node, *elem, node_num, remaining_vars);
1565  }
1566  else
1567  sides.emplace
1568  (node, std::make_tuple(elem, side, std::move(remaining_vars)));
1569  }
1570  }
1571 
1572  // Elements with elemental dofs might need those copied too.
1573  if (copy_this_elem)
1574  {
1575  for (auto v_num : this->projector.variables)
1576  {
1577  const Variable & var = system.variable(v_num);
1578  if (!var.active_on_subdomain(elem->subdomain_id()))
1579  continue;
1580  FEType fe_type = var.type();
1581 
1582  std::vector<FValue> Ue;
1583  std::vector<dof_id_type> elem_dof_ids;
1584  f.eval_old_dofs(*elem, fe_type, sys_num, v_num,
1585  elem_dof_ids, Ue);
1586  action.insert(elem_dof_ids, Ue);
1587 
1588  if (has_interior_nodes)
1589  {
1590  std::vector<FValue> Un;
1591  std::vector<dof_id_type> node_dof_ids;
1592 
1593  f.eval_old_dofs(*elem, n_nodes-1, v_num, node_dof_ids, Un);
1594  action.insert(node_dof_ids, Un);
1595  }
1596  }
1597  }
1598  }
1599 }
1600 
1601 
1602 template <typename FFunctor, typename GFunctor,
1603  typename FValue, typename ProjectionAction>
1606 {
1607  auto merge_multimaps = [](ves_multimap & self, const ves_multimap & other_mm)
1608  {
1609  for (const auto & pair : other_mm)
1610  {
1611  const Node * node = pair.first;
1612  auto remaining_vars = std::get<2>(pair.second);
1613 
1614  auto my_range = self.equal_range(node);
1615  for (const auto & v_ent : as_range(my_range))
1616  for (const unsigned int var_covered :
1617  std::get<2>(v_ent.second))
1618  remaining_vars.erase(var_covered);
1619 
1620  if (!remaining_vars.empty())
1621  self.emplace
1622  (node, std::make_tuple(std::get<0>(pair.second),
1623  std::get<1>(pair.second),
1624  std::move(remaining_vars)));
1625  }
1626  };
1627 
1628  merge_multimaps(vertices, other.vertices);
1629  merge_multimaps(edges, other.edges);
1630  merge_multimaps(sides, other.sides);
1631 }
1632 
1633 
1634 template <typename FFunctor, typename GFunctor,
1635  typename FValue, typename ProjectionAction>
1637  (const node_range & range)
1638 {
1639  LOG_SCOPE ("project_vertices","GenericProjector");
1640 
1641  const unsigned int sys_num = system.number();
1642 
1643  // Variables with extra hanging dofs can't safely use eval_at_node
1644  // in as many places as variables without can.
1645  std::vector<unsigned short> extra_hanging_dofs;
1646  for (auto v_num : this->projector.variables)
1647  {
1648  if (extra_hanging_dofs.size() <= v_num)
1649  extra_hanging_dofs.resize(v_num+1, false);
1650  extra_hanging_dofs[v_num] =
1651  FEInterface::extra_hanging_dofs(system.variable(v_num).type());
1652  }
1653 
1654  for (const auto & v_pair : range)
1655  {
1656  const Node & vertex = *v_pair.first;
1657  const Elem & elem = *std::get<0>(v_pair.second);
1658  const unsigned int n = std::get<1>(v_pair.second);
1659  const var_set & vertex_vars = std::get<2>(v_pair.second);
1660 
1661  context.pre_fe_reinit(system, &elem);
1662 
1663  this->find_dofs_to_send(vertex, elem, n, vertex_vars);
1664 
1665  // Look at all the variables we're supposed to interpolate from
1666  // this element on this vertex
1667  for (const auto & var : vertex_vars)
1668  {
1669  const Variable & variable = system.variable(var);
1670  const FEType & base_fe_type = variable.type();
1671  const unsigned int var_component =
1672  system.variable_scalar_number(var, 0);
1673 
1674  if (base_fe_type.family == SCALAR)
1675  continue;
1676 
1677  const FEContinuity & cont = this->conts[var];
1678  if (cont == DISCONTINUOUS)
1679  {
1680  libmesh_assert_equal_to(vertex.n_comp(sys_num, var), 0);
1681  }
1682  else if (cont == C_ZERO)
1683  {
1684  libmesh_assert(vertex.n_comp(sys_num, var));
1685  const dof_id_type id = vertex.dof_number(sys_num, var, 0);
1686  // C_ZERO elements have a single nodal value DoF at vertices
1687  const FValue val = f.eval_at_node
1688  (context, var_component, /*dim=*/ 0, // Don't care w/C0
1689  vertex, extra_hanging_dofs[var], system.time);
1690  insert_id(id, val, vertex.processor_id());
1691  }
1692  else if (cont == C_ONE)
1693  {
1694  libmesh_assert(vertex.n_comp(sys_num, var));
1695  const dof_id_type first_id = vertex.dof_number(sys_num, var, 0);
1696 
1697  // C_ONE elements have a single nodal value and dim
1698  // gradient values at vertices, as well as cross
1699  // gradients for HERMITE. We need to have an element in
1700  // hand to figure out dim and to have in case this
1701  // vertex is a new vertex.
1702  const int dim = elem.dim();
1703 #ifndef NDEBUG
1704  // For now all C1 elements at a vertex had better have
1705  // the same dimension. If anyone hits these asserts let
1706  // me know; we could probably support a mixed-dimension
1707  // mesh IFF the 2D elements were all parallel to xy and
1708  // the 1D elements all parallel to x.
1709  for (const auto e_id : (*this->projector.nodes_to_elem)[vertex.id()])
1710  {
1711  const Elem & e = system.get_mesh().elem_ref(e_id);
1712  libmesh_assert_equal_to(dim, e.dim());
1713  }
1714 #endif
1715 #ifdef LIBMESH_ENABLE_AMR
1716  bool is_parent_vertex = false;
1717  if (elem.parent())
1718  {
1719  const int i_am_child =
1720  elem.parent()->which_child_am_i(&elem);
1721  is_parent_vertex =
1722  elem.parent()->is_vertex_on_parent(i_am_child, n);
1723  }
1724 #else
1725  const bool is_parent_vertex = false;
1726 #endif
1727 
1728  // The hermite element vertex shape functions are weird
1729  if (base_fe_type.family == HERMITE)
1730  {
1731  const FValue val =
1732  f.eval_at_node(context,
1733  var_component,
1734  dim,
1735  vertex,
1736  extra_hanging_dofs[var],
1737  system.time);
1738  insert_id(first_id, val, vertex.processor_id());
1739 
1740  VectorValue<FValue> grad =
1741  is_parent_vertex ?
1742  g->eval_at_node(context,
1743  var_component,
1744  dim,
1745  vertex,
1746  extra_hanging_dofs[var],
1747  system.time) :
1748  g->eval_at_point(context,
1749  var_component,
1750  vertex,
1751  system.time);
1752  // x derivative
1753  insert_id(first_id+1, grad(0),
1754  vertex.processor_id());
1755  if (dim > 1)
1756  {
1757  // We'll finite difference mixed derivatives
1758  Real delta_x = TOLERANCE * elem.hmin();
1759 
1760  Point nxminus = elem.point(n),
1761  nxplus = elem.point(n);
1762  nxminus(0) -= delta_x;
1763  nxplus(0) += delta_x;
1764  VectorValue<FValue> gxminus =
1765  g->eval_at_point(context,
1766  var_component,
1767  nxminus,
1768  system.time);
1769  VectorValue<FValue> gxplus =
1770  g->eval_at_point(context,
1771  var_component,
1772  nxplus,
1773  system.time);
1774  // y derivative
1775  insert_id(first_id+2, grad(1),
1776  vertex.processor_id());
1777  // xy derivative
1778  insert_id(first_id+3,
1779  (gxplus(1) - gxminus(1)) / 2. / delta_x,
1780  vertex.processor_id());
1781 
1782  if (dim > 2)
1783  {
1784  // z derivative
1785  insert_id(first_id+4, grad(2),
1786  vertex.processor_id());
1787  // xz derivative
1788  insert_id(first_id+5,
1789  (gxplus(2) - gxminus(2)) / 2. / delta_x,
1790  vertex.processor_id());
1791 
1792  // We need new points for yz
1793  Point nyminus = elem.point(n),
1794  nyplus = elem.point(n);
1795  nyminus(1) -= delta_x;
1796  nyplus(1) += delta_x;
1797  VectorValue<FValue> gyminus =
1798  g->eval_at_point(context,
1799  var_component,
1800  nyminus,
1801  system.time);
1802  VectorValue<FValue> gyplus =
1803  g->eval_at_point(context,
1804  var_component,
1805  nyplus,
1806  system.time);
1807  // yz derivative
1808  insert_id(first_id+6,
1809  (gyplus(2) - gyminus(2)) / 2. / delta_x,
1810  vertex.processor_id());
1811  // Getting a 2nd order xyz is more tedious
1812  Point nxmym = elem.point(n),
1813  nxmyp = elem.point(n),
1814  nxpym = elem.point(n),
1815  nxpyp = elem.point(n);
1816  nxmym(0) -= delta_x;
1817  nxmym(1) -= delta_x;
1818  nxmyp(0) -= delta_x;
1819  nxmyp(1) += delta_x;
1820  nxpym(0) += delta_x;
1821  nxpym(1) -= delta_x;
1822  nxpyp(0) += delta_x;
1823  nxpyp(1) += delta_x;
1824  VectorValue<FValue> gxmym =
1825  g->eval_at_point(context,
1826  var_component,
1827  nxmym,
1828  system.time);
1829  VectorValue<FValue> gxmyp =
1830  g->eval_at_point(context,
1831  var_component,
1832  nxmyp,
1833  system.time);
1834  VectorValue<FValue> gxpym =
1835  g->eval_at_point(context,
1836  var_component,
1837  nxpym,
1838  system.time);
1839  VectorValue<FValue> gxpyp =
1840  g->eval_at_point(context,
1841  var_component,
1842  nxpyp,
1843  system.time);
1844  FValue gxzplus = (gxpyp(2) - gxmyp(2))
1845  / 2. / delta_x;
1846  FValue gxzminus = (gxpym(2) - gxmym(2))
1847  / 2. / delta_x;
1848  // xyz derivative
1849  insert_id(first_id+7,
1850  (gxzplus - gxzminus) / 2. / delta_x,
1851  vertex.processor_id());
1852  }
1853  }
1854  }
1855  else
1856  {
1857  // Currently other C_ONE elements have a single nodal
1858  // value shape function and nodal gradient component
1859  // shape functions
1860  libmesh_assert_equal_to
1862  (dim, base_fe_type, elem.type(),
1863  elem.get_node_index(&vertex)),
1864  (unsigned int)(1 + dim));
1865  const FValue val =
1866  f.eval_at_node(context, var_component, dim,
1867  vertex, extra_hanging_dofs[var],
1868  system.time);
1869  insert_id(first_id, val, vertex.processor_id());
1870  VectorValue<FValue> grad =
1871  is_parent_vertex ?
1872  g->eval_at_node(context, var_component, dim,
1873  vertex, extra_hanging_dofs[var],
1874  system.time) :
1875  g->eval_at_point(context, var_component, vertex,
1876  system.time);
1877  for (int i=0; i!= dim; ++i)
1878  insert_id(first_id + i + 1, grad(i),
1879  vertex.processor_id());
1880  }
1881  }
1882  else
1883  libmesh_error_msg("Unknown continuity " << cont);
1884  }
1885  }
1886 }
1887 
1888 
1889 template <typename FFunctor, typename GFunctor,
1890  typename FValue, typename ProjectionAction>
1892  (const node_range & range)
1893 {
1894  LOG_SCOPE ("project_edges","GenericProjector");
1895 
1896  const unsigned int sys_num = system.number();
1897 
1898  for (const auto & e_pair : range)
1899  {
1900  const Elem & elem = *std::get<0>(e_pair.second);
1901 
1902  // If this is an unchanged element then we already copied all
1903  // its dofs
1904 #ifdef LIBMESH_ENABLE_AMR
1905  if (f.is_grid_projection() &&
1906  (elem.refinement_flag() != Elem::JUST_REFINED &&
1910  continue;
1911 #endif // LIBMESH_ENABLE_AMR
1912 
1913  const Node & edge_node = *e_pair.first;
1914  const int dim = elem.dim();
1915  const var_set & edge_vars = std::get<2>(e_pair.second);
1916 
1917  const unsigned short edge_num = std::get<1>(e_pair.second);
1918  const unsigned short node_num = elem.n_vertices() + edge_num;
1919  context.edge = edge_num;
1920  context.pre_fe_reinit(system, &elem);
1921 
1922  this->find_dofs_to_send(edge_node, elem, node_num, edge_vars);
1923 
1924  // Look at all the variables we're supposed to interpolate from
1925  // this element on this edge
1926  for (const auto & var : edge_vars)
1927  {
1928  const Variable & variable = system.variable(var);
1929  const FEType & base_fe_type = variable.type();
1930  const unsigned int var_component =
1931  system.variable_scalar_number(var, 0);
1932 
1933  if (base_fe_type.family == SCALAR)
1934  continue;
1935 
1936  FEType fe_type = base_fe_type;
1937 
1938  // This may be a p refined element
1939  fe_type.order =
1940  libMesh::Order (fe_type.order + elem.p_level());
1941 
1942  // If this is a Lagrange element with DoFs on edges then by
1943  // convention we interpolate at the node rather than project
1944  // along the edge.
1945  if (fe_type.family == LAGRANGE)
1946  {
1947  if (fe_type.order > 1)
1948  {
1949  const dof_id_type dof_id =
1950  edge_node.dof_number(sys_num, var, 0);
1951  const processor_id_type pid =
1952  edge_node.processor_id();
1953  FValue fval = f.eval_at_point
1954  (context, var_component, edge_node, system.time);
1955  insert_id(dof_id, fval, pid);
1956  }
1957  continue;
1958  }
1959 
1960  // If this is a low order monomial element which has merely
1961  // been h refined then we already copied all its dofs
1962  if (fe_type.family == MONOMIAL &&
1963  fe_type.order == CONSTANT &&
1966  continue;
1967 
1968  // FIXME: Need to generalize this to vector-valued elements. [PB]
1969  FEBase * fe = nullptr;
1970  context.get_element_fe( var, fe, dim );
1971  FEBase * edge_fe = nullptr;
1972  context.get_edge_fe( var, edge_fe );
1973 
1974  // If we're JUST_COARSENED we'll need a custom
1975  // evaluation, not just the standard edge FE
1976  const FEBase & proj_fe =
1977 #ifdef LIBMESH_ENABLE_AMR
1979  *fe :
1980 #endif
1981  *edge_fe;
1982 
1983 #ifdef LIBMESH_ENABLE_AMR
1984  if (elem.refinement_flag() == Elem::JUST_COARSENED)
1985  {
1986  std::vector<Point> fine_points;
1987 
1988  std::unique_ptr<FEBase> fine_fe
1989  (FEBase::build (dim, base_fe_type));
1990 
1991  std::unique_ptr<QBase> qrule
1992  (base_fe_type.default_quadrature_rule(1));
1993  fine_fe->attach_quadrature_rule(qrule.get());
1994 
1995  const std::vector<Point> & child_xyz =
1996  fine_fe->get_xyz();
1997 
1998  for (unsigned int c = 0, nc = elem.n_children();
1999  c != nc; ++c)
2000  {
2001  if (!elem.is_child_on_edge(c, context.edge))
2002  continue;
2003 
2004  fine_fe->edge_reinit(elem.child_ptr(c), context.edge);
2005  fine_points.insert(fine_points.end(),
2006  child_xyz.begin(),
2007  child_xyz.end());
2008  }
2009 
2010  std::vector<Point> fine_qp;
2011  FEInterface::inverse_map (dim, base_fe_type, &elem,
2012  fine_points, fine_qp);
2013 
2014  context.elem_fe_reinit(&fine_qp);
2015  }
2016  else
2017 #endif // LIBMESH_ENABLE_AMR
2018  context.edge_fe_reinit();
2019 
2020  const std::vector<dof_id_type> & dof_indices =
2021  context.get_dof_indices(var);
2022 
2023  std::vector<unsigned int> edge_dofs;
2024  FEInterface::dofs_on_edge(&elem, dim, base_fe_type,
2025  context.edge, edge_dofs);
2026 
2027  this->construct_projection
2028  (dof_indices, edge_dofs, var_component,
2029  &edge_node, proj_fe);
2030  }
2031  }
2032 }
2033 
2034 
2035 template <typename FFunctor, typename GFunctor,
2036  typename FValue, typename ProjectionAction>
2038  (const node_range & range)
2039 {
2040  LOG_SCOPE ("project_sides","GenericProjector");
2041 
2042  const unsigned int sys_num = system.number();
2043 
2044  for (const auto & s_pair : range)
2045  {
2046  const Elem & elem = *std::get<0>(s_pair.second);
2047 
2048  // If this is an unchanged element then we already copied all
2049  // its dofs
2050 #ifdef LIBMESH_ENABLE_AMR
2051  if (f.is_grid_projection() &&
2052  (elem.refinement_flag() != Elem::JUST_REFINED &&
2056  continue;
2057 #endif // LIBMESH_ENABLE_AMR
2058 
2059  const Node & side_node = *s_pair.first;
2060  const int dim = elem.dim();
2061  const var_set & side_vars = std::get<2>(s_pair.second);
2062 
2063  const unsigned int side_num = std::get<1>(s_pair.second);
2064  unsigned short node_num = elem.n_vertices()+side_num;
2065  // In 3D only some sides may have nodes
2066  if (dim == 3)
2067  for (auto n : IntRange<unsigned int>(0, elem.n_nodes()))
2068  {
2069  if (!elem.is_face(n))
2070  continue;
2071 
2072  if (elem.is_node_on_side(n, side_num))
2073  {
2074  node_num = n;
2075  break;
2076  }
2077  }
2078 
2079  context.side = side_num;
2080  context.pre_fe_reinit(system, &elem);
2081 
2082  this->find_dofs_to_send(side_node, elem, node_num, side_vars);
2083 
2084  // Look at all the variables we're supposed to interpolate from
2085  // this element on this side
2086  for (const auto & var : side_vars)
2087  {
2088  const Variable & variable = system.variable(var);
2089  const FEType & base_fe_type = variable.type();
2090  const unsigned int var_component =
2091  system.variable_scalar_number(var, 0);
2092 
2093  if (base_fe_type.family == SCALAR)
2094  continue;
2095 
2096  FEType fe_type = base_fe_type;
2097 
2098  // This may be a p refined element
2099  fe_type.order =
2100  libMesh::Order (fe_type.order + elem.p_level());
2101 
2102  // If this is a Lagrange element with DoFs on sides then by
2103  // convention we interpolate at the node rather than project
2104  // along the side.
2105  if (fe_type.family == LAGRANGE)
2106  {
2107  if (fe_type.order > 1)
2108  {
2109  const dof_id_type dof_id =
2110  side_node.dof_number(sys_num, var, 0);
2111  const processor_id_type pid =
2112  side_node.processor_id();
2113  FValue fval = f.eval_at_point
2114  (context, var_component, side_node, system.time);
2115  insert_id(dof_id, fval, pid);
2116  }
2117  continue;
2118  }
2119 
2120  // If this is a low order monomial element which has merely
2121  // been h refined then we already copied all its dofs
2122  if (fe_type.family == MONOMIAL &&
2123  fe_type.order == CONSTANT &&
2126  continue;
2127 
2128  // FIXME: Need to generalize this to vector-valued elements. [PB]
2129  FEBase * fe = nullptr;
2130  context.get_element_fe( var, fe, dim );
2131  FEBase * side_fe = nullptr;
2132  context.get_side_fe( var, side_fe );
2133 
2134  // If we're JUST_COARSENED we'll need a custom
2135  // evaluation, not just the standard side FE
2136  const FEBase & proj_fe =
2137 #ifdef LIBMESH_ENABLE_AMR
2139  *fe :
2140 #endif
2141  *side_fe;
2142 
2143 #ifdef LIBMESH_ENABLE_AMR
2144  if (elem.refinement_flag() == Elem::JUST_COARSENED)
2145  {
2146  std::vector<Point> fine_points;
2147 
2148  std::unique_ptr<FEBase> fine_fe
2149  (FEBase::build (dim, base_fe_type));
2150 
2151  std::unique_ptr<QBase> qrule
2152  (base_fe_type.default_quadrature_rule(1));
2153  fine_fe->attach_quadrature_rule(qrule.get());
2154 
2155  const std::vector<Point> & child_xyz =
2156  fine_fe->get_xyz();
2157 
2158  for (unsigned int c = 0, nc = elem.n_children();
2159  c != nc; ++c)
2160  {
2161  if (!elem.is_child_on_side(c, context.side))
2162  continue;
2163 
2164  fine_fe->reinit(elem.child_ptr(c), context.side);
2165  fine_points.insert(fine_points.end(),
2166  child_xyz.begin(),
2167  child_xyz.end());
2168  }
2169 
2170  std::vector<Point> fine_qp;
2171  FEInterface::inverse_map (dim, base_fe_type, &elem,
2172  fine_points, fine_qp);
2173 
2174  context.elem_fe_reinit(&fine_qp);
2175  }
2176  else
2177 #endif // LIBMESH_ENABLE_AMR
2178  context.side_fe_reinit();
2179 
2180  const std::vector<dof_id_type> & dof_indices =
2181  context.get_dof_indices(var);
2182 
2183  std::vector<unsigned int> side_dofs;
2184  FEInterface::dofs_on_side(&elem, dim, base_fe_type,
2185  context.side, side_dofs);
2186 
2187  this->construct_projection
2188  (dof_indices, side_dofs, var_component,
2189  &side_node, proj_fe);
2190  }
2191  }
2192 }
2193 
2194 
2195 template <typename FFunctor, typename GFunctor,
2196  typename FValue, typename ProjectionAction>
2198  (const interior_range & range)
2199 {
2200  LOG_SCOPE ("project_interiors","GenericProjector");
2201 
2202  const unsigned int sys_num = system.number();
2203 
2204  // Iterate over all dof-bearing element interiors in the range
2205  for (const auto & elem : range)
2206  {
2207  unsigned char dim = cast_int<unsigned char>(elem->dim());
2208 
2209  context.pre_fe_reinit(system, elem);
2210 
2211  // Loop over all the variables we've been requested to project, to
2212  // do the projection
2213  for (const auto & var : this->projector.variables)
2214  {
2215  const Variable & variable = system.variable(var);
2216 
2217  if (!variable.active_on_subdomain(elem->subdomain_id()))
2218  continue;
2219 
2220  const FEType & base_fe_type = variable.type();
2221 
2222  if (base_fe_type.family == SCALAR)
2223  continue;
2224 
2225  FEBase * fe = nullptr;
2226  context.get_element_fe( var, fe, dim );
2227 
2228  FEType fe_type = base_fe_type;
2229 
2230  // This may be a p refined element
2231  fe_type.order =
2232  libMesh::Order (fe_type.order + elem->p_level());
2233 
2234  const unsigned int var_component =
2235  system.variable_scalar_number(var, 0);
2236 
2237  // If this is a Lagrange element with interior DoFs then by
2238  // convention we interpolate at the interior node rather
2239  // than project along the interior.
2240  if (fe_type.family == LAGRANGE)
2241  {
2242  if (fe_type.order > 1)
2243  {
2244  const unsigned int first_interior_node =
2245  (elem->n_vertices() +
2246  ((elem->dim() > 2) * elem->n_edges()) +
2247  ((elem->dim() > 1) * elem->n_sides()));
2248  const unsigned int n_nodes = elem->n_nodes();
2249 
2250  // < vs != is important here for HEX20, QUAD8!
2251  for (unsigned int n = first_interior_node; n < n_nodes; ++n)
2252  {
2253  const Node & interior_node = elem->node_ref(n);
2254  const dof_id_type dof_id =
2255  interior_node.dof_number(sys_num, var, 0);
2256  const processor_id_type pid =
2257  interior_node.processor_id();
2258  FValue fval = f.eval_at_point
2259  (context, var_component, interior_node, system.time);
2260  insert_id(dof_id, fval, pid);
2261  }
2262  }
2263  continue;
2264  }
2265 
2266 #ifdef LIBMESH_ENABLE_AMR
2267  if (elem->refinement_flag() == Elem::JUST_COARSENED)
2268  {
2269  std::vector<Point> fine_points;
2270 
2271  std::unique_ptr<FEBase> fine_fe
2272  (FEBase::build (dim, base_fe_type));
2273 
2274  std::unique_ptr<QBase> qrule
2275  (base_fe_type.default_quadrature_rule(dim));
2276  fine_fe->attach_quadrature_rule(qrule.get());
2277 
2278  const std::vector<Point> & child_xyz =
2279  fine_fe->get_xyz();
2280 
2281  for (auto & child : elem->child_ref_range())
2282  {
2283  fine_fe->reinit(&child);
2284  fine_points.insert(fine_points.end(),
2285  child_xyz.begin(),
2286  child_xyz.end());
2287  }
2288 
2289  std::vector<Point> fine_qp;
2290  FEInterface::inverse_map (dim, base_fe_type, elem,
2291  fine_points, fine_qp);
2292 
2293  context.elem_fe_reinit(&fine_qp);
2294  }
2295  else
2296 #endif // LIBMESH_ENABLE_AMR
2297  context.elem_fe_reinit();
2298 
2299  const std::vector<dof_id_type> & dof_indices =
2300  context.get_dof_indices(var);
2301 
2302  const unsigned int n_dofs =
2303  cast_int<unsigned int>(dof_indices.size());
2304 
2305  std::vector<unsigned int> all_dofs(n_dofs);
2306  std::iota(all_dofs.begin(), all_dofs.end(), 0);
2307 
2308  this->construct_projection
2309  (dof_indices, all_dofs, var_component,
2310  nullptr, *fe);
2311  } // end variables loop
2312  } // end elements loop
2313 }
2314 
2315 
2316 
2317 template <typename FFunctor, typename GFunctor,
2318  typename FValue, typename ProjectionAction>
2319 void
2320 GenericProjector<FFunctor, GFunctor, FValue,
2321  ProjectionAction>::SubFunctor::find_dofs_to_send
2322  (const Node & node, const Elem & elem, unsigned short node_num, const var_set & vars)
2323 {
2324  libmesh_assert (&node == elem.node_ptr(node_num));
2325 
2326  // Any ghosted node in our range might have an owner who needs our
2327  // data
2328  const processor_id_type owner = node.processor_id();
2329  if (owner != system.processor_id())
2330  {
2331  const MeshBase & mesh = system.get_mesh();
2332  const DofMap & dof_map = system.get_dof_map();
2333 
2334  // But let's check and see if we can be certain the owner can
2335  // compute any or all of its own dof coefficients on that node.
2336  std::vector<dof_id_type> node_dof_ids, patch_dof_ids;
2337  for (const auto & var : vars)
2338  {
2339  const Variable & variable = system.variable(var);
2340 
2341  if (!variable.active_on_subdomain(elem.subdomain_id()))
2342  continue;
2343 
2344  dof_map.dof_indices(elem, node_num, node_dof_ids, var);
2345  }
2346  libmesh_assert(std::is_sorted(node_dof_ids.begin(),
2347  node_dof_ids.end()));
2348  const std::vector<dof_id_type> & patch =
2349  (*this->projector.nodes_to_elem)[node.id()];
2350  for (const auto & elem_id : patch)
2351  {
2352  const Elem & patch_elem = mesh.elem_ref(elem_id);
2353  if (!patch_elem.active() || owner != patch_elem.processor_id())
2354  continue;
2355  dof_map.dof_indices(&patch_elem, patch_dof_ids);
2356  std::sort(patch_dof_ids.begin(), patch_dof_ids.end());
2357 
2358  // Remove any node_dof_ids that we see can be calculated on
2359  // this element
2360  std::vector<dof_id_type> diff_ids(node_dof_ids.size());
2361  auto it = std::set_difference(node_dof_ids.begin(), node_dof_ids.end(),
2362  patch_dof_ids.begin(), patch_dof_ids.end(), diff_ids.begin());
2363  diff_ids.resize(it-diff_ids.begin());
2364  node_dof_ids.swap(diff_ids);
2365  if (node_dof_ids.empty())
2366  break;
2367  }
2368 
2369  // Give ids_to_push default invalid pid: not yet computed
2370  for (auto id : node_dof_ids)
2371  new_ids_to_push[id].second = DofObject::invalid_processor_id;
2372  }
2373 }
2374 
2375 
2376 
2377 template <typename FFunctor, typename GFunctor,
2378  typename FValue, typename ProjectionAction>
2379 void
2380 GenericProjector<FFunctor, GFunctor, FValue,
2381  ProjectionAction>::send_and_insert_dof_values
2382  (std::unordered_map<dof_id_type, std::pair<FValue, processor_id_type>> & ids_to_push,
2383  ProjectionAction & action) const
2384 {
2385  // See if we calculated any ids that need to be pushed; get them
2386  // ready to push.
2387  std::unordered_map<processor_id_type, std::vector<dof_id_type>>
2388  pushed_dof_ids, received_dof_ids;
2389  std::unordered_map<processor_id_type, std::vector<typename TypeToSend<FValue>::type>>
2390  pushed_dof_values, received_dof_values;
2391  for (auto & id_val_pid : ids_to_push)
2392  {
2393  processor_id_type pid = id_val_pid.second.second;
2395  {
2396  pushed_dof_ids[pid].push_back(id_val_pid.first);
2397  pushed_dof_values[pid].push_back(convert_to_send(id_val_pid.second.first));
2398  }
2399  }
2400 
2401  // If and when we get ids pushed to us, act on them if we have
2402  // corresponding values or save them if not
2403  auto ids_action_functor =
2404  [&action, &received_dof_ids, &received_dof_values]
2405  (processor_id_type pid,
2406  const std::vector<dof_id_type> & data)
2407  {
2408  auto iter = received_dof_values.find(pid);
2409  if (iter == received_dof_values.end())
2410  {
2411  libmesh_assert(received_dof_ids.find(pid) ==
2412  received_dof_ids.end());
2413  received_dof_ids[pid] = data;
2414  }
2415  else
2416  {
2417  auto & vals = iter->second;
2418  std::size_t vals_size = vals.size();
2419  libmesh_assert_equal_to(vals_size, data.size());
2420  for (std::size_t i=0; i != vals_size; ++i)
2421  {
2422  FValue val;
2423  convert_from_receive(vals[i], val);
2424  action.insert(data[i], val);
2425  }
2426  received_dof_values.erase(iter);
2427  }
2428  };
2429 
2430  // If and when we get values pushed to us, act on them if we have
2431  // corresponding ids or save them if not
2432  auto values_action_functor =
2433  [&action, &received_dof_ids, &received_dof_values]
2434  (processor_id_type pid,
2435  const std::vector<typename TypeToSend<FValue>::type> & data)
2436  {
2437  auto iter = received_dof_ids.find(pid);
2438  if (iter == received_dof_ids.end())
2439  {
2440  // We get no more than 1 values vector from anywhere
2441  libmesh_assert(received_dof_values.find(pid) ==
2442  received_dof_values.end());
2443  received_dof_values[pid] = data;
2444  }
2445  else
2446  {
2447  auto & ids = iter->second;
2448  std::size_t ids_size = ids.size();
2449  libmesh_assert_equal_to(ids_size, data.size());
2450  for (std::size_t i=0; i != ids_size; ++i)
2451  {
2452  FValue val;
2453  convert_from_receive(data[i], val);
2454  action.insert(ids[i], val);
2455  }
2456  received_dof_ids.erase(iter);
2457  }
2458  };
2459 
2461  (system.comm(), pushed_dof_ids, ids_action_functor);
2462 
2464  (system.comm(), pushed_dof_values, values_action_functor);
2465 
2466  // At this point we shouldn't have any unprocessed data left
2467  libmesh_assert(received_dof_ids.empty());
2468  libmesh_assert(received_dof_values.empty());
2469 
2470 }
2471 
2472 
2473 
2474 template <typename FFunctor, typename GFunctor,
2475  typename FValue, typename ProjectionAction>
2476 void
2478  (const std::vector<dof_id_type> & dof_indices_var,
2479  const std::vector<unsigned int> & involved_dofs,
2480  unsigned int var_component,
2481  const Node * node,
2482  const FEBase & fe)
2483 {
2484  const std::vector<Real> & JxW = fe.get_JxW();
2485  const std::vector<std::vector<Real>> & phi = fe.get_phi();
2486  const std::vector<std::vector<RealGradient>> * dphi = nullptr;
2487  const std::vector<Point> & xyz_values = fe.get_xyz();
2488  const FEContinuity cont = fe.get_continuity();
2489  const std::unordered_map<dof_id_type, FValue> & ids_to_save =
2490  this->projector.ids_to_save;
2491 
2492  if (cont == C_ONE)
2493  dphi = &(fe.get_dphi());
2494 
2495  const unsigned int n_involved_dofs =
2496  cast_int<unsigned int>(involved_dofs.size());
2497 
2498  std::vector<dof_id_type> free_dof_ids;
2499  DenseVector<FValue> Uinvolved(n_involved_dofs);
2500  std::vector<char> dof_is_fixed(n_involved_dofs, false); // bools
2501 
2502  for (auto i : IntRange<unsigned int>(0, n_involved_dofs))
2503  {
2504  const dof_id_type id = dof_indices_var[involved_dofs[i]];
2505  auto iter = ids_to_save.find(id);
2506  if (iter == ids_to_save.end())
2507  free_dof_ids.push_back(id);
2508  else
2509  {
2510  dof_is_fixed[i] = true;
2511  Uinvolved(i) = iter->second;
2512  }
2513  }
2514 
2515  const unsigned int free_dofs = free_dof_ids.size();
2516 
2517  // There may be nothing to project
2518  if (!free_dofs)
2519  return;
2520 
2521  // The element matrix and RHS for projections.
2522  // Note that Ke is always real-valued, whereas
2523  // Fe may be complex valued if complex number
2524  // support is enabled
2525  DenseMatrix<Real> Ke(free_dofs, free_dofs);
2526  DenseVector<FValue> Fe(free_dofs);
2527  // The new degree of freedom coefficients to solve for
2528  DenseVector<FValue> Ufree(free_dofs);
2529 
2530  const unsigned int n_qp =
2531  cast_int<unsigned int>(xyz_values.size());
2532 
2533  // Loop over the quadrature points
2534  for (unsigned int qp=0; qp<n_qp; qp++)
2535  {
2536  // solution at the quadrature point
2537  FValue fineval = f.eval_at_point(context,
2538  var_component,
2539  xyz_values[qp],
2540  system.time);
2541  // solution grad at the quadrature point
2542  VectorValue<FValue> finegrad;
2543  if (cont == C_ONE)
2544  finegrad = g->eval_at_point(context,
2545  var_component,
2546  xyz_values[qp],
2547  system.time);
2548 
2549  // Form edge projection matrix
2550  for (unsigned int sidei=0, freei=0;
2551  sidei != n_involved_dofs; ++sidei)
2552  {
2553  unsigned int i = involved_dofs[sidei];
2554  // fixed DoFs aren't test functions
2555  if (dof_is_fixed[sidei])
2556  continue;
2557  for (unsigned int sidej=0, freej=0;
2558  sidej != n_involved_dofs; ++sidej)
2559  {
2560  unsigned int j = involved_dofs[sidej];
2561  if (dof_is_fixed[sidej])
2562  Fe(freei) -= phi[i][qp] * phi[j][qp] *
2563  JxW[qp] * Uinvolved(sidej);
2564  else
2565  Ke(freei,freej) += phi[i][qp] *
2566  phi[j][qp] * JxW[qp];
2567  if (cont == C_ONE)
2568  {
2569  if (dof_is_fixed[sidej])
2570  Fe(freei) -= ( (*dphi)[i][qp] *
2571  (*dphi)[j][qp] ) *
2572  JxW[qp] * Uinvolved(sidej);
2573  else
2574  Ke(freei,freej) += ( (*dphi)[i][qp] *
2575  (*dphi)[j][qp] )
2576  * JxW[qp];
2577  }
2578  if (!dof_is_fixed[sidej])
2579  freej++;
2580  }
2581  Fe(freei) += phi[i][qp] * fineval * JxW[qp];
2582  if (cont == C_ONE)
2583  Fe(freei) += (finegrad * (*dphi)[i][qp] ) *
2584  JxW[qp];
2585  freei++;
2586  }
2587  }
2588 
2589  Ke.cholesky_solve(Fe, Ufree);
2590 
2591  // Transfer new edge solutions to element
2592  const processor_id_type pid = node ?
2594  insert_ids(free_dof_ids, Ufree.get_values(), pid);
2595 }
2596 
2597 
2598 } // namespace libMesh
2599 
2600 #endif // GENERIC_PROJECTOR_H
The GenericProjector class implements the core of other projection operations, using two input functo...
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:179
virtual unsigned int size() const override
Definition: dense_vector.h:92
FEFamily family
The type of finite element.
Definition: fe_type.h:204
void check_old_context(const FEMContext &c)
RefinementState refinement_flag() const
Definition: elem.h:2599
ElemType
Defines an enum for geometric element types.
dof_id_type n_nodes(const MeshBase::const_node_iterator &begin, const MeshBase::const_node_iterator &end)
Count up the number of nodes of a specific type (as defined by an iterator range).
dof_id_type dof_number(const unsigned int s, const unsigned int var, const unsigned int comp) const
Definition: dof_object.h:927
const Elem * parent() const
Definition: elem.h:2453
Order
defines an enum for polynomial orders.
Definition: enum_order.h:40
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:2167
const std::vector< Point > & get_xyz() const
Definition: fe_abstract.h:238
bool check_old_context(const FEMContext &c, const Point &p)
const NumericVector< Number > & old_solution
unsigned int get_node_index(const Node *node_ptr) const
Definition: elem.h:2012
unsigned int n_comp(const unsigned int s, const unsigned int var) const
Definition: dof_object.h:897
std::unique_ptr< QBase > default_quadrature_rule(const unsigned int dim, const int extraorder=0) const
void operator()(const node_range &range)
virtual bool is_face(const unsigned int i) const =0
void init_context(FEMContext &c)
The IntRange templated class is intended to make it easy to loop over integers which are indices of a...
Definition: int_range.h:51
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Fills the vector di with the global degree of freedom indices for the element.
const Elem & get_elem() const
Accessor for current Elem object.
Definition: fem_context.h:883
ProjectSides(ProjectSides &p_s, Threads::split)
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
Dummy "splitting object" used to distinguish splitting constructors from copy constructors.
Definition: threads_none.h:63
virtual Real hmax() const
RefinementState p_refinement_flag() const
Definition: elem.h:2615
std::vector< T > & get_values()
Definition: dense_vector.h:256
FEMFunctionWrapper(const FEMFunctionWrapper< Output > &fw)
std::unordered_multimap< const Node *, std::tuple< const Elem *, unsigned short, var_set > > ves_multimap
GenericProjector(const System &system_in, const FFunctor &f_in, const GFunctor *g_in, const ProjectionAction &act_in, const std::vector< unsigned int > &variables_in, std::unordered_map< dof_id_type, std::vector< dof_id_type >> *nodes_to_elem_in=nullptr)
This is the base class from which all geometric element types are derived.
Definition: elem.h:101
IntRange< std::size_t > index_range(const std::vector< T > &vec)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition: int_range.h:104
FEMFunctionWrapper(const FEMFunctionBase< Output > &f)
RefinementState
Enumeration of possible element refinement states.
Definition: elem.h:1149
void insert_id(dof_id_type id, const FValue &val, processor_id_type pid)
virtual bool is_node_on_side(const unsigned int n, const unsigned int s) const =0
virtual unsigned int n_children() const =0
unsigned int p_level() const
Definition: elem.h:2529
void eval_old_dofs(const Elem &elem, unsigned int node_num, unsigned int var_num, std::vector< dof_id_type > &indices, std::vector< Output > &values)
OrderWrapper order
The approximation order of the element.
Definition: fe_type.h:198
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:960
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...
static const Real TOLERANCE
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:52
OldSolutionBase(const libMesh::System &sys_in)
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 > &)
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< Output > &values)
NumericVector< Val > & target_vector
VectorSetAction(NumericVector< Val > &target_vec)
const MeshBase & get_mesh() const
Definition: system.h:2067
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 is a duplication of the SGI STL extension std::iota.
Definition: utility.h:57
uint8_t processor_id_type
Definition: id_types.h:104
static void get_shape_outputs(FEBase &fe)
This is the MeshBase class.
Definition: mesh_base.h:77
SimpleRange< ChildRefIter > child_ref_range()
Returns a range with all children of a parent element, usable in range-based for loops.
Definition: elem.h:1807
virtual FEContinuity get_continuity() const =0
void push_parallel_vector_data(const Communicator &comm, const MapToVectors &data, const ActionFunctor &act_on_data)
Send and receive and act on vectors of data.
FEType get_fe_type() const
Definition: fe_abstract.h:428
std::vector< const Elem * > interiors
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:177
std::unordered_map< dof_id_type, FValue > new_ids_to_save
const std::vector< std::vector< OutputGradient > > & get_dphi() const
Definition: fe_base.h:215
virtual bool is_vertex_on_parent(unsigned int c, unsigned int n) const
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)
Fills the vector di with the local degree of freedom indices associated with edge e of element elem A...
void eval_old_dofs(const Elem &, unsigned int, unsigned int, std::vector< dof_id_type > &, std::vector< Output > &)
This class defines the notion of a variable in the system.
Definition: variable.h:49
dof_id_type id() const
Definition: dof_object.h:738
dof_id_type numeric_index_type
Definition: id_types.h:99
void project(const ConstElemRange &range)
Function definitions.
virtual unsigned int n_nodes() const =0
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:404
unsigned int n_vars(const unsigned int s, const unsigned int vg) const
Definition: dof_object.h:862
StoredRange< std::vector< node_projection >::const_iterator, node_projection > node_range
Manages consistently variables, degrees of freedom, and coefficient vectors.
Definition: system.h:92
unsigned int n_systems() const
Definition: dof_object.h:832
const std::vector< Real > & get_JxW() const
Definition: fe_abstract.h:245
void insert(dof_id_type id, Val val)
ProjectEdges(ProjectEdges &p_e, Threads::split)
static std::unique_ptr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.
SimpleRange< I > as_range(const std::pair< I, I > &p)
Helper function that allows us to treat a homogenous pair as a range.
Definition: simple_range.h:57
static const dof_id_type invalid_id
An invalid id to distinguish an uninitialized DofObject.
Definition: dof_object.h:393
This class provides all data required for a physics package (e.g.
Definition: fem_context.h:61
void init_context(FEMContext &c)
void send_and_insert_dof_values(std::unordered_map< dof_id_type, std::pair< FValue, processor_id_type >> &ids_to_push, ProjectionAction &action) const
std::unique_ptr< FEMFunctionBase< Output > > _f
bool active_on_subdomain(subdomain_id_type sid) const
Definition: variable.h:136
static Point inverse_map(const unsigned int dim, const FEType &fe_t, const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true)
const ProjectionAction & master_action
std::unordered_map< dof_id_type, FValue > ids_to_save
void insert(const std::vector< dof_id_type > &dof_indices, const DenseVector< Val > &Ue)
const std::vector< unsigned int > & variables
unsigned int which_child_am_i(const Elem *e) const
Definition: elem.h:2581
void operator()(const node_range &range)
virtual numeric_index_type first_local_index() const =0
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
static unsigned int n_dofs_per_elem(const unsigned int dim, const FEType &fe_t, const ElemType t)
std::unordered_map< dof_id_type, std::vector< dof_id_type > > * nodes_to_elem
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)
DofObject * old_dof_object
This object on the last mesh.
Definition: dof_object.h:79
virtual Real hmin() const
subdomain_id_type subdomain_id() const
Definition: elem.h:2034
Output eval_at_point(const FEMContext &c, unsigned int i, const Point &p, Real=0.)
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:92
const Node * node_ptr(const unsigned int i) const
Definition: elem.h:1974
Output eval_at_point(const FEMContext &c, unsigned int i, const Point &n, const Real time)
For ease of communication, we allow users to translate their own value types to a more easily computa...
virtual const Elem & elem_ref(const dof_id_type i) const
Definition: mesh_base.h:520
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
std::pair< unsigned int, unsigned int > var_to_vg_and_offset(const unsigned int s, const unsigned int var) const
Definition: dof_object.h:1042
OldSolutionValue(const libMesh::System &sys_in, const NumericVector< Number > &old_sol)
SortAndCopy(SortAndCopy &other, Threads::split)
void insert_ids(const std::vector< dof_id_type > &ids, const std::vector< FValue > &vals, processor_id_type pid)
virtual void set(const numeric_index_type i, const T value)=0
Sets v(i) = value.
std::vector< unsigned int > component_to_var
static unsigned int n_dofs_at_node(const unsigned int dim, const FEType &fe_t, const ElemType t, const unsigned int n)
void cholesky_solve(const DenseVector< T2 > &b, DenseVector< T2 > &x)
For symmetric positive definite (SPD) matrices.
Defines a dense vector for use in Finite Element-type computations.
std::unordered_map< dof_id_type, std::pair< FValue, processor_id_type > > new_ids_to_push
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
IterBase * data
Ideally this private member data should have protected access.
unsigned int n_vars() const
Definition: system.h:2139
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:910
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:2364
processor_id_type processor_id() const
Definition: dof_object.h:800
virtual ElemType type() const =0
Output eval_at_node(const FEMContext &c, unsigned int i, unsigned int, const Node &n, bool, const Real time)
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:38
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)
Fills the vector di with the local degree of freedom indices associated with side s of element elem A...
const Point & point(const unsigned int i) const
Definition: elem.h:1920
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 FEBase &fe)
void join(const SortAndCopy &other)
virtual bool is_child_on_edge(const unsigned int c, const unsigned int e) const
virtual numeric_index_type last_local_index() const =0
This class forms the foundation from which generic finite elements may be derived.
std::vector< FEContinuity > conts
GenericProjector(const GenericProjector &in)
const Elem * child_ptr(unsigned int i) const
Definition: elem.h:2552
uint8_t dof_id_type
Definition: id_types.h:67
static bool extra_hanging_dofs(const FEType &fe_t)
const std::vector< std::vector< OutputShape > > & get_phi() const
Definition: fe_base.h:207
const FEType & type() const
Definition: variable.h:119
spin_mutex spin_mtx
A convenient spin mutex object which can be used for obtaining locks.