libMesh
system_projection.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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 // C++ includes
21 #include <vector>
22 
23 // Local includes
24 #include "libmesh/boundary_info.h"
25 #include "libmesh/dense_matrix.h"
26 #include "libmesh/dense_vector.h"
27 #include "libmesh/dirichlet_boundaries.h"
28 #include "libmesh/dof_map.h"
29 #include "libmesh/elem.h"
30 #include "libmesh/equation_systems.h"
31 #include "libmesh/fe_base.h"
32 #include "libmesh/fe_interface.h"
33 #include "libmesh/libmesh_logging.h"
34 #include "libmesh/mesh_base.h"
35 #include "libmesh/numeric_vector.h"
36 #include "libmesh/quadrature_gauss.h"
37 #include "libmesh/system.h"
38 #include "libmesh/threads.h"
39 #include "libmesh/wrapped_function.h"
40 #include "libmesh/wrapped_functor.h"
41 
42 namespace libMesh
43 {
44 
45 // ------------------------------------------------------------
46 // Helper class definitions
47 
52 template <typename FFunctor, typename GFunctor,
53  typename FValue, typename ProjectionAction>
55 {
56 private:
57  const System & system;
58 
59  // For TBB compatibility and thread safety we'll copy these in
60  // operator()
61  const FFunctor & master_f;
62  const GFunctor * master_g; // Needed for C1 type elements only
64  const ProjectionAction & master_action;
65  const std::vector<unsigned int> & variables;
66 
67 public:
68  GenericProjector (const System & system_in,
69  const FFunctor & f_in,
70  const GFunctor * g_in,
71  const ProjectionAction & act_in,
72  const std::vector<unsigned int> & variables_in) :
73  system(system_in),
74  master_f(f_in),
75  master_g(g_in),
76  g_was_copied(false),
77  master_action(act_in),
78  variables(variables_in)
79  {}
80 
82  system(in.system),
83  master_f(in.master_f),
84  master_g(in.master_g ? new GFunctor(*in.master_g) : libmesh_nullptr),
85  g_was_copied(in.master_g),
86  master_action(in.master_action),
87  variables(in.variables)
88  {}
89 
91  {
92  if (g_was_copied)
93  delete master_g;
94  }
95 
96  void operator() (const ConstElemRange & range) const;
97 };
98 
99 
106 template <typename Val>
108 {
109 private:
111 
112 public:
114  target_vector(target_vec) {}
115 
116  void insert(const FEMContext & c,
117  unsigned int var_num,
118  const DenseVector<Val> & Ue)
119  {
120  const numeric_index_type
121  first = target_vector.first_local_index(),
122  last = target_vector.last_local_index();
123 
124  const std::vector<dof_id_type> & dof_indices =
125  c.get_dof_indices(var_num);
126 
127  unsigned int size = Ue.size();
128 
129  libmesh_assert_equal_to(size, dof_indices.size());
130 
131  // Lock the new vector since it is shared among threads.
132  {
133  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
134 
135  for (unsigned int i = 0; i != size; ++i)
136  if ((dof_indices[i] >= first) && (dof_indices[i] < last))
137  target_vector.set(dof_indices[i], Ue(i));
138  }
139  }
140 };
141 
142 
143 template <typename Output>
145 {
146 public:
147  FEMFunctionWrapper(const FEMFunctionBase<Output> & f) : _f(f.clone()) {}
148 
150  _f(fw._f->clone()) {}
151 
152  void init_context (FEMContext & c) { _f->init_context(c); }
153 
154  Output eval_at_node (const FEMContext & c,
155  unsigned int i,
156  unsigned int /*elem_dim*/,
157  const Node & n,
158  const Real time)
159  { return _f->component(c, i, n, time); }
160 
161  Output eval_at_point (const FEMContext & c,
162  unsigned int i,
163  const Point & n,
164  const Real time)
165  { return _f->component(c, i, n, time); }
166 
167  bool is_grid_projection() { return false; }
168 
169  void eval_old_dofs (const FEMContext & /* c */,
170  unsigned int /* var_component */,
171  std::vector<Output> /* values */)
172  { libmesh_error(); }
173 
174 private:
176 };
177 
178 
179 #ifdef LIBMESH_ENABLE_AMR
180 template <typename Output,
181  void (FEMContext::*point_output) (unsigned int,
182  const Point &,
183  Output &,
184  const Real) const>
186 {
187 public:
189  const NumericVector<Number> & old_sol) :
190  last_elem(libmesh_nullptr),
191  sys(sys_in),
192  old_context(sys_in),
193  old_solution(old_sol)
194  {
195  old_context.set_algebraic_type(FEMContext::OLD);
196  old_context.set_custom_solution(&old_solution);
197  }
198 
200  last_elem(libmesh_nullptr),
201  sys(in.sys),
202  old_context(sys),
203  old_solution(in.old_solution)
204  {
205  old_context.set_algebraic_type(FEMContext::OLD);
206  old_context.set_custom_solution(&old_solution);
207  }
208 
209  static void get_shape_outputs(FEBase & fe);
210 
211  // Integrating on new mesh elements, we won't yet have an up to date
212  // current_local_solution.
214  {
216 
217  // Loop over variables, to prerequest
218  for (unsigned int var=0; var!=sys.n_vars(); ++var)
219  {
220  FEBase * fe = libmesh_nullptr;
221  const std::set<unsigned char> & elem_dims =
222  old_context.elem_dimensions();
223 
224  for (std::set<unsigned char>::const_iterator dim_it =
225  elem_dims.begin(); dim_it != elem_dims.end(); ++dim_it)
226  {
227  const unsigned char dim = *dim_it;
228  old_context.get_element_fe( var, fe, dim );
229  get_shape_outputs(*fe);
230  }
231  }
232  }
233 
234 
235 
236  Output eval_at_node (const FEMContext & c,
237  unsigned int i,
238  unsigned int elem_dim,
239  const Node & n,
240  Real /* time */ =0.);
241 
242  Output eval_at_point(const FEMContext & c,
243  unsigned int i,
244  const Point & p,
245  Real /* time */ =0.)
246  {
247  LOG_SCOPE ("eval_at_point()", "OldSolutionValue");
248 
249  if (!this->check_old_context(c, p))
250  return 0;
251 
252  Output n;
253  (old_context.*point_output)(i, p, n, out_of_elem_tol);
254  return n;
255  }
256 
257  bool is_grid_projection() { return true; }
258 
259  void eval_old_dofs (const FEMContext & c,
260  unsigned int var,
261  std::vector<Output> & values)
262  {
263  LOG_SCOPE ("eval_old_dofs()", "OldSolutionValue");
264 
265  this->check_old_context(c);
266 
267  const std::vector<dof_id_type> & old_dof_indices =
268  old_context.get_dof_indices(var);
269 
270  libmesh_assert_equal_to (old_dof_indices.size(), values.size());
271 
272  old_solution.get(old_dof_indices, values);
273  }
274 
275 protected:
276  void check_old_context (const FEMContext & c)
277  {
278  LOG_SCOPE ("check_old_context(c)", "OldSolutionValue");
279  const Elem & elem = c.get_elem();
280  if (last_elem != &elem)
281  {
282  if (elem.refinement_flag() == Elem::JUST_REFINED)
283  {
284  old_context.pre_fe_reinit(sys, elem.parent());
285  }
286  else if (elem.refinement_flag() == Elem::JUST_COARSENED)
287  {
288  libmesh_error();
289  }
290  else
291  {
292  if (!elem.old_dof_object)
293  {
294  libmesh_error();
295  }
296 
297  old_context.pre_fe_reinit(sys, &elem);
298  }
299 
300  last_elem = &elem;
301  }
302  else
303  {
304  libmesh_assert(old_context.has_elem());
305  }
306  }
307 
308 
309  bool check_old_context (const FEMContext & c, const Point & p)
310  {
311  LOG_SCOPE ("check_old_context(c,p)", "OldSolutionValue");
312  const Elem & elem = c.get_elem();
313  if (last_elem != &elem)
314  {
315  if (elem.refinement_flag() == Elem::JUST_REFINED)
316  {
317  old_context.pre_fe_reinit(sys, elem.parent());
318  }
319  else if (elem.refinement_flag() == Elem::JUST_COARSENED)
320  {
321  // Find the child with this point. Use out_of_elem_tol
322  // (in physical space, which may correspond to a large
323  // tolerance in master space!) to allow for out-of-element
324  // finite differencing of mixed gradient terms. Pray we
325  // have no quadrature locations which are within 1e-5 of
326  // the element subdivision boundary but are not exactly on
327  // that boundary.
328  const Real master_tol = out_of_elem_tol / elem.hmax() * 2;
329 
330  for (auto & child : elem.child_ref_range())
331  if (child.close_to_point(p, master_tol))
332  {
333  old_context.pre_fe_reinit(sys, &child);
334  break;
335  }
336 
338  (old_context.get_elem().close_to_point(p, master_tol));
339  }
340  else
341  {
342  if (!elem.old_dof_object)
343  return false;
344 
345  old_context.pre_fe_reinit(sys, &elem);
346  }
347 
348  last_elem = &elem;
349  }
350  else
351  {
352  libmesh_assert(old_context.has_elem());
353 
354  const Real master_tol = out_of_elem_tol / elem.hmax() * 2;
355 
356  if (!old_context.get_elem().close_to_point(p, master_tol))
357  {
358  libmesh_assert_equal_to
360 
361  for (auto & child : elem.child_ref_range())
362  if (child.close_to_point(p, master_tol))
363  {
364  old_context.pre_fe_reinit(sys, &child);
365  break;
366  }
367 
369  (old_context.get_elem().close_to_point(p, master_tol));
370  }
371  }
372 
373  return true;
374  }
375 
376 private:
377  const Elem * last_elem;
378  const System & sys;
381 
382  static const Real out_of_elem_tol;
383 };
384 
385 
386 template<>
387 inline void
389 {
390  fe.get_phi();
391 }
392 
393 
394 template<>
395 inline void
397 {
398  fe.get_dphi();
399 }
400 
401 
402 template<>
403 inline
404 Number
407  unsigned int i,
408  unsigned int /* elem_dim */,
409  const Node & n,
410  Real /* time */)
411 {
412  LOG_SCOPE ("Number eval_at_node()", "OldSolutionValue");
413 
414  // Optimize for the common case, where this node was part of the
415  // old solution.
416  //
417  // Be sure to handle cases where the variable wasn't defined on
418  // this node (due to changing subdomain support) or where the
419  // variable has no components on this node (due to Elem order
420  // exceeding FE order)
421  if (n.old_dof_object &&
423  n.old_dof_object->n_comp(sys.number(), i))
424  {
425  const dof_id_type old_id =
426  n.old_dof_object->dof_number(sys.number(), i, 0);
427  return old_solution(old_id);
428  }
429 
430  return this->eval_at_point(c, i, n, 0);
431 }
432 
433 
434 
435 template<>
436 inline
437 Gradient
440  unsigned int i,
441  unsigned int elem_dim,
442  const Node & n,
443  Real /* time */)
444 {
445  LOG_SCOPE ("Gradient eval_at_node()", "OldSolutionValue");
446 
447  // Optimize for the common case, where this node was part of the
448  // old solution.
449  //
450  // Be sure to handle cases where the variable wasn't defined on
451  // this node (due to changing subdomain support) or where the
452  // variable has no components on this node (due to Elem order
453  // exceeding FE order)
454  if (n.old_dof_object &&
456  n.old_dof_object->n_comp(sys.number(), i))
457  {
458  Gradient g;
459  for (unsigned int d = 0; d != elem_dim; ++d)
460  {
461  const dof_id_type old_id =
462  n.old_dof_object->dof_number(sys.number(), i, d+1);
463  g(d) = old_solution(old_id);
464  }
465  return g;
466  }
467 
468  return this->eval_at_point(c, i, n, 0);
469 }
470 
471 
472 
473 
474 
475 template <>
477 
478 template <>
480 
481 
492 {
493 private:
494  const System & system;
495 
496 public:
497  BuildProjectionList (const System & system_in) :
498  system(system_in),
499  send_list()
500  {}
501 
503  system(other.system),
504  send_list()
505  {}
506 
507  void unique();
508  void operator()(const ConstElemRange & range);
509  void join (const BuildProjectionList & other);
510  std::vector<dof_id_type> send_list;
511 };
512 
513 #endif // LIBMESH_ENABLE_AMR
514 
515 
522 {
523 private:
524  const std::set<boundary_id_type> & b;
525  const std::vector<unsigned int> & variables;
526  const System & system;
531 
532 public:
533  BoundaryProjectSolution (const std::set<boundary_id_type> & b_in,
534  const std::vector<unsigned int> & variables_in,
535  const System & system_in,
536  FunctionBase<Number> * f_in,
537  FunctionBase<Gradient> * g_in,
538  const Parameters & parameters_in,
539  NumericVector<Number> & new_v_in) :
540  b(b_in),
541  variables(variables_in),
542  system(system_in),
543  f(f_in ? f_in->clone() : UniquePtr<FunctionBase<Number>>()),
544  g(g_in ? g_in->clone() : UniquePtr<FunctionBase<Gradient>>()),
545  parameters(parameters_in),
546  new_vector(new_v_in)
547  {
548  libmesh_assert(f.get());
549  f->init();
550  if (g.get())
551  g->init();
552  }
553 
555  b(in.b),
556  variables(in.variables),
557  system(in.system),
558  f(in.f.get() ? in.f->clone() : UniquePtr<FunctionBase<Number>>()),
559  g(in.g.get() ? in.g->clone() : UniquePtr<FunctionBase<Gradient>>()),
560  parameters(in.parameters),
561  new_vector(in.new_vector)
562  {
563  libmesh_assert(f.get());
564  f->init();
565  if (g.get())
566  g->init();
567  }
568 
569  void operator()(const ConstElemRange & range) const;
570 };
571 
572 
573 
574 // ------------------------------------------------------------
575 // System implementation
577  int is_adjoint) const
578 {
579  // Create a copy of the vector, which currently
580  // contains the old data.
582  old_vector (vector.clone());
583 
584  // Project the old vector to the new vector
585  this->project_vector (*old_vector, vector, is_adjoint);
586 }
587 
588 
595  NumericVector<Number> & new_v,
596  int
597 #ifdef LIBMESH_ENABLE_AMR
598  is_adjoint
599 #endif
600  ) const
601 {
602  LOG_SCOPE ("project_vector(old,new)", "System");
603 
610  new_v.clear();
611 
612 #ifdef LIBMESH_ENABLE_AMR
613 
614  // Resize the new vector and get a serial version.
615  NumericVector<Number> * new_vector_ptr = libmesh_nullptr;
616  UniquePtr<NumericVector<Number>> new_vector_built;
617  NumericVector<Number> * local_old_vector;
618  UniquePtr<NumericVector<Number>> local_old_vector_built;
619  const NumericVector<Number> * old_vector_ptr = libmesh_nullptr;
620 
621  ConstElemRange active_local_elem_range
622  (this->get_mesh().active_local_elements_begin(),
623  this->get_mesh().active_local_elements_end());
624 
625  // If the old vector was uniprocessor, make the new
626  // vector uniprocessor
627  if (old_v.type() == SERIAL)
628  {
629  new_v.init (this->n_dofs(), false, SERIAL);
630  new_vector_ptr = &new_v;
631  old_vector_ptr = &old_v;
632  }
633 
634  // Otherwise it is a parallel, distributed vector, which
635  // we need to localize.
636  else if (old_v.type() == PARALLEL)
637  {
638  // Build a send list for efficient localization
639  BuildProjectionList projection_list(*this);
640  Threads::parallel_reduce (active_local_elem_range,
641  projection_list);
642 
643  // Create a sorted, unique send_list
644  projection_list.unique();
645 
646  new_v.init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
647  new_vector_built = NumericVector<Number>::build(this->comm());
648  local_old_vector_built = NumericVector<Number>::build(this->comm());
649  new_vector_ptr = new_vector_built.get();
650  local_old_vector = local_old_vector_built.get();
651  new_vector_ptr->init(this->n_dofs(), false, SERIAL);
652  local_old_vector->init(old_v.size(), false, SERIAL);
653  old_v.localize(*local_old_vector, projection_list.send_list);
654  local_old_vector->close();
655  old_vector_ptr = local_old_vector;
656  }
657  else if (old_v.type() == GHOSTED)
658  {
659  // Build a send list for efficient localization
660  BuildProjectionList projection_list(*this);
661  Threads::parallel_reduce (active_local_elem_range,
662  projection_list);
663 
664  // Create a sorted, unique send_list
665  projection_list.unique();
666 
667  new_v.init (this->n_dofs(), this->n_local_dofs(),
668  this->get_dof_map().get_send_list(), false, GHOSTED);
669 
670  local_old_vector_built = NumericVector<Number>::build(this->comm());
671  new_vector_ptr = &new_v;
672  local_old_vector = local_old_vector_built.get();
673  local_old_vector->init(old_v.size(), old_v.local_size(),
674  projection_list.send_list, false, GHOSTED);
675  old_v.localize(*local_old_vector, projection_list.send_list);
676  local_old_vector->close();
677  old_vector_ptr = local_old_vector;
678  }
679  else // unknown old_v.type()
680  libmesh_error_msg("ERROR: Unknown old_v.type() == " << old_v.type());
681 
682  // Note that the above will have zeroed the new_vector.
683  // Just to be sure, assert that new_vector_ptr and old_vector_ptr
684  // were successfully set before trying to deref them.
685  libmesh_assert(new_vector_ptr);
686  libmesh_assert(old_vector_ptr);
687 
688  NumericVector<Number> & new_vector = *new_vector_ptr;
689  const NumericVector<Number> & old_vector = *old_vector_ptr;
690 
691  const unsigned int n_variables = this->n_vars();
692 
693  if (n_variables)
694  {
695  std::vector<unsigned int> vars(n_variables);
696  for (unsigned int i=0; i != n_variables; ++i)
697  vars[i] = i;
698 
699  // Use a typedef to make the calling sequence for parallel_for() a bit more readable
700  typedef
703  Number, VectorSetAction<Number>> FEMProjector;
704 
706  OldSolutionValue<Gradient, &FEMContext::point_gradient> g(*this, old_vector);
707  VectorSetAction<Number> setter(new_vector);
708 
709  Threads::parallel_for (active_local_elem_range,
710  FEMProjector(*this, f, &g, setter, vars));
711 
712  // Copy the SCALAR dofs from old_vector to new_vector
713  // Note: We assume that all SCALAR dofs are on the
714  // processor with highest ID
715  if (this->processor_id() == (this->n_processors()-1))
716  {
717  const DofMap & dof_map = this->get_dof_map();
718  for (unsigned int var=0; var<this->n_vars(); var++)
719  if (this->variable(var).type().family == SCALAR)
720  {
721  // We can just map SCALAR dofs directly across
722  std::vector<dof_id_type> new_SCALAR_indices, old_SCALAR_indices;
723  dof_map.SCALAR_dof_indices (new_SCALAR_indices, var, false);
724  dof_map.SCALAR_dof_indices (old_SCALAR_indices, var, true);
725  const unsigned int new_n_dofs =
726  cast_int<unsigned int>(new_SCALAR_indices.size());
727 
728  for (unsigned int i=0; i<new_n_dofs; i++)
729  {
730  new_vector.set( new_SCALAR_indices[i], old_vector(old_SCALAR_indices[i]) );
731  }
732  }
733  }
734  }
735 
736  new_vector.close();
737 
738  // If the old vector was serial, we probably need to send our values
739  // to other processors
740  //
741  // FIXME: I'm not sure how to make a NumericVector do that without
742  // creating a temporary parallel vector to use localize! - RHS
743  if (old_v.type() == SERIAL)
744  {
746  dist_v->init(this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
747  dist_v->close();
748 
749  for (dof_id_type i=0; i!=dist_v->size(); i++)
750  if (new_vector(i) != 0.0)
751  dist_v->set(i, new_vector(i));
752 
753  dist_v->close();
754 
755  dist_v->localize (new_v, this->get_dof_map().get_send_list());
756  new_v.close();
757  }
758  // If the old vector was parallel, we need to update it
759  // and free the localized copies
760  else if (old_v.type() == PARALLEL)
761  {
762  // We may have to set dof values that this processor doesn't
763  // own in certain special cases, like LAGRANGE FIRST or
764  // HERMITE THIRD elements on second-order meshes
765  for (dof_id_type i=0; i!=new_v.size(); i++)
766  if (new_vector(i) != 0.0)
767  new_v.set(i, new_vector(i));
768  new_v.close();
769  }
770 
771  if (is_adjoint == -1)
772  this->get_dof_map().enforce_constraints_exactly(*this, &new_v);
773  else if (is_adjoint >= 0)
774  this->get_dof_map().enforce_adjoint_constraints_exactly(new_v,
775  is_adjoint);
776 
777 #else
778 
779  // AMR is disabled: simply copy the vector
780  new_v = old_v;
781 
782 #endif // #ifdef LIBMESH_ENABLE_AMR
783 }
784 
785 
786 
791 void System::project_solution (Number fptr(const Point & p,
792  const Parameters & parameters,
793  const std::string & sys_name,
794  const std::string & unknown_name),
795  Gradient gptr(const Point & p,
796  const Parameters & parameters,
797  const std::string & sys_name,
798  const std::string & unknown_name),
799  const Parameters & parameters) const
800 {
801  WrappedFunction<Number> f(*this, fptr, &parameters);
802  WrappedFunction<Gradient> g(*this, gptr, &parameters);
803  this->project_solution(&f, &g);
804 }
805 
806 
812  FunctionBase<Gradient> * g) const
813 {
814  this->project_vector(*solution, f, g);
815 
816  solution->localize(*current_local_solution, _dof_map->get_send_list());
817 }
818 
819 
825  FEMFunctionBase<Gradient> * g) const
826 {
827  this->project_vector(*solution, f, g);
828 
829  solution->localize(*current_local_solution, _dof_map->get_send_list());
830 }
831 
832 
837 void System::project_vector (Number fptr(const Point & p,
838  const Parameters & parameters,
839  const std::string & sys_name,
840  const std::string & unknown_name),
841  Gradient gptr(const Point & p,
842  const Parameters & parameters,
843  const std::string & sys_name,
844  const std::string & unknown_name),
845  const Parameters & parameters,
846  NumericVector<Number> & new_vector,
847  int is_adjoint) const
848 {
849  WrappedFunction<Number> f(*this, fptr, &parameters);
850  WrappedFunction<Gradient> g(*this, gptr, &parameters);
851  this->project_vector(new_vector, &f, &g, is_adjoint);
852 }
853 
861  int is_adjoint) const
862 {
863  LOG_SCOPE ("project_vector(FunctionBase)", "System");
864 
865  libmesh_assert(f);
866 
867  WrappedFunctor<Number> f_fem(*f);
868 
869  if (g)
870  {
871  WrappedFunctor<Gradient> g_fem(*g);
872 
873  this->project_vector(new_vector, &f_fem, &g_fem, is_adjoint);
874  }
875  else
876  this->project_vector(new_vector, &f_fem, libmesh_nullptr, is_adjoint);
877 }
878 
879 
887  int is_adjoint) const
888 {
889  LOG_SCOPE ("project_fem_vector()", "System");
890 
891  libmesh_assert (f);
892 
893  ConstElemRange active_local_range
894  (this->get_mesh().active_local_elements_begin(),
895  this->get_mesh().active_local_elements_end() );
896 
897  VectorSetAction<Number> setter(new_vector);
898 
899  const unsigned int n_variables = this->n_vars();
900 
901  std::vector<unsigned int> vars(n_variables);
902  for (unsigned int i=0; i != n_variables; ++i)
903  vars[i] = i;
904 
905  // Use a typedef to make the calling sequence for parallel_for() a bit more readable
906  typedef
908  Number, VectorSetAction<Number>> FEMProjector;
909 
911 
912  if (g)
913  {
914  FEMFunctionWrapper<Gradient> gw(*g);
915 
917  (active_local_range,
918  FEMProjector(*this, fw, &gw, setter, vars));
919  }
920  else
922  (active_local_range,
923  FEMProjector(*this, fw, libmesh_nullptr, setter, vars));
924 
925  // Also, load values into the SCALAR dofs
926  // Note: We assume that all SCALAR dofs are on the
927  // processor with highest ID
928  if (this->processor_id() == (this->n_processors()-1))
929  {
930  // FIXME: Do we want to first check for SCALAR vars before building this? [PB]
931  FEMContext context( *this );
932 
933  const DofMap & dof_map = this->get_dof_map();
934  for (unsigned int var=0; var<this->n_vars(); var++)
935  if (this->variable(var).type().family == SCALAR)
936  {
937  // FIXME: We reinit with an arbitrary element in case the user
938  // doesn't override FEMFunctionBase::component. Is there
939  // any use case we're missing? [PB]
940  Elem * el = const_cast<Elem *>(*(this->get_mesh().active_local_elements_begin()));
941  context.pre_fe_reinit(*this, el);
942 
943  std::vector<dof_id_type> SCALAR_indices;
944  dof_map.SCALAR_dof_indices (SCALAR_indices, var);
945  const unsigned int n_SCALAR_dofs =
946  cast_int<unsigned int>(SCALAR_indices.size());
947 
948  for (unsigned int i=0; i<n_SCALAR_dofs; i++)
949  {
950  const dof_id_type global_index = SCALAR_indices[i];
951  const unsigned int component_index =
952  this->variable_scalar_number(var,i);
953 
954  new_vector.set(global_index, f->component(context, component_index, Point(), this->time));
955  }
956  }
957  }
958 
959  new_vector.close();
960 
961 #ifdef LIBMESH_ENABLE_CONSTRAINTS
962  if (is_adjoint == -1)
963  this->get_dof_map().enforce_constraints_exactly(*this, &new_vector);
964  else if (is_adjoint >= 0)
965  this->get_dof_map().enforce_adjoint_constraints_exactly(new_vector,
966  is_adjoint);
967 #endif
968 }
969 
970 
976 void System::boundary_project_solution (const std::set<boundary_id_type> & b,
977  const std::vector<unsigned int> & variables,
978  Number fptr(const Point & p,
979  const Parameters & parameters,
980  const std::string & sys_name,
981  const std::string & unknown_name),
982  Gradient gptr(const Point & p,
983  const Parameters & parameters,
984  const std::string & sys_name,
985  const std::string & unknown_name),
986  const Parameters & parameters)
987 {
988  WrappedFunction<Number> f(*this, fptr, &parameters);
989  WrappedFunction<Gradient> g(*this, gptr, &parameters);
990  this->boundary_project_solution(b, variables, &f, &g);
991 }
992 
993 
999 void System::boundary_project_solution (const std::set<boundary_id_type> & b,
1000  const std::vector<unsigned int> & variables,
1003 {
1004  this->boundary_project_vector(b, variables, *solution, f, g);
1005 
1006  solution->localize(*current_local_solution);
1007 }
1008 
1009 
1010 
1011 
1012 
1017 void System::boundary_project_vector (const std::set<boundary_id_type> & b,
1018  const std::vector<unsigned int> & variables,
1019  Number fptr(const Point & p,
1020  const Parameters & parameters,
1021  const std::string & sys_name,
1022  const std::string & unknown_name),
1023  Gradient gptr(const Point & p,
1024  const Parameters & parameters,
1025  const std::string & sys_name,
1026  const std::string & unknown_name),
1027  const Parameters & parameters,
1028  NumericVector<Number> & new_vector,
1029  int is_adjoint) const
1030 {
1031  WrappedFunction<Number> f(*this, fptr, &parameters);
1032  WrappedFunction<Gradient> g(*this, gptr, &parameters);
1033  this->boundary_project_vector(b, variables, new_vector, &f, &g,
1034  is_adjoint);
1035 }
1036 
1041 void System::boundary_project_vector (const std::set<boundary_id_type> & b,
1042  const std::vector<unsigned int> & variables,
1043  NumericVector<Number> & new_vector,
1046  int is_adjoint) const
1047 {
1048  LOG_SCOPE ("boundary_project_vector()", "System");
1049 
1051  (ConstElemRange (this->get_mesh().active_local_elements_begin(),
1052  this->get_mesh().active_local_elements_end() ),
1053  BoundaryProjectSolution(b, variables, *this, f, g,
1054  this->get_equation_systems().parameters,
1055  new_vector)
1056  );
1057 
1058  // We don't do SCALAR dofs when just projecting the boundary, so
1059  // we're done here.
1060 
1061  new_vector.close();
1062 
1063 #ifdef LIBMESH_ENABLE_CONSTRAINTS
1064  if (is_adjoint == -1)
1065  this->get_dof_map().enforce_constraints_exactly(*this, &new_vector);
1066  else if (is_adjoint >= 0)
1067  this->get_dof_map().enforce_adjoint_constraints_exactly(new_vector,
1068  is_adjoint);
1069 #endif
1070 }
1071 
1072 
1073 
1074 template <typename FFunctor, typename GFunctor,
1075  typename FValue, typename ProjectionAction>
1077  (const ConstElemRange & range) const
1078 {
1079  LOG_SCOPE ("operator()", "GenericProjector");
1080 
1081  ProjectionAction action(master_action);
1082  FFunctor f(master_f);
1084  if (master_g)
1085  g.reset(new GFunctor(*master_g));
1086 
1087  // The DofMap for this system
1088  const DofMap & dof_map = system.get_dof_map();
1089 
1090  // The element matrix and RHS for projections.
1091  // Note that Ke is always real-valued, whereas
1092  // Fe may be complex valued if complex number
1093  // support is enabled
1094  DenseMatrix<Real> Ke;
1096  // The new element degree of freedom coefficients
1098 
1099  // Context objects to contain all our required FE objects
1100  FEMContext context( system );
1101 
1102  // Loop over all the variables we've been requested to project, to
1103  // pre-request
1104  for (std::size_t v=0; v!=variables.size(); v++)
1105  {
1106  const unsigned int var = variables[v];
1107 
1108  // FIXME: Need to generalize this to vector-valued elements. [PB]
1109  FEBase * fe = libmesh_nullptr;
1110  FEBase * side_fe = libmesh_nullptr;
1111  FEBase * edge_fe = libmesh_nullptr;
1112 
1113  const std::set<unsigned char> & elem_dims =
1114  context.elem_dimensions();
1115 
1116  for (std::set<unsigned char>::const_iterator dim_it =
1117  elem_dims.begin(); dim_it != elem_dims.end(); ++dim_it)
1118  {
1119  const unsigned char dim = *dim_it;
1120 
1121  context.get_element_fe( var, fe, dim );
1122  if (fe->get_fe_type().family == SCALAR)
1123  continue;
1124  if (dim > 1)
1125  context.get_side_fe( var, side_fe, dim );
1126  if (dim > 2)
1127  context.get_edge_fe( var, edge_fe );
1128 
1129  fe->get_xyz();
1130 
1131  fe->get_phi();
1132  if (dim > 1)
1133  side_fe->get_phi();
1134  if (dim > 2)
1135  edge_fe->get_phi();
1136 
1137  const FEContinuity cont = fe->get_continuity();
1138  if (cont == C_ONE)
1139  {
1140  // Our C1 elements need gradient information
1141  libmesh_assert(g);
1142 
1143  fe->get_dphi();
1144  if (dim > 1)
1145  side_fe->get_dphi();
1146  if (dim > 2)
1147  edge_fe->get_dphi();
1148  }
1149  }
1150  }
1151 
1152  // Now initialize any user requested shape functions, xyz vals, etc
1153  f.init_context(context);
1154  if (g.get())
1155  g->init_context(context);
1156 
1157  // this->init_context(context);
1158 
1159  // Iterate over all the elements in the range
1160  for (ConstElemRange::const_iterator elem_it=range.begin(); elem_it != range.end();
1161  ++elem_it)
1162  {
1163  const Elem * elem = *elem_it;
1164 
1165  unsigned int dim = elem->dim();
1166 
1167  context.pre_fe_reinit(system, elem);
1168 
1169  // If we're doing AMR, this might be a grid projection with a cheap
1170  // early exit.
1171 #ifdef LIBMESH_ENABLE_AMR
1172  // If this element doesn't have an old_dof_object, but it
1173  // wasn't just refined or just coarsened into activity, then
1174  // it must be newly added, so the user is responsible for
1175  // setting the new dofs on it during a grid projection.
1176  if (!elem->old_dof_object &&
1177  elem->refinement_flag() != Elem::JUST_REFINED &&
1179  f.is_grid_projection())
1180  continue;
1181 #endif // LIBMESH_ENABLE_AMR
1182 
1183  // Loop over all the variables we've been requested to project, to
1184  // do the projection
1185  for (std::size_t v=0; v!=variables.size(); v++)
1186  {
1187  const unsigned int var = variables[v];
1188 
1189  const Variable & variable = dof_map.variable(var);
1190 
1191  const FEType & base_fe_type = variable.type();
1192 
1193  FEType fe_type = base_fe_type;
1194 
1195  // This may be a p refined element
1196  fe_type.order =
1197  libMesh::Order (fe_type.order + elem->p_level());
1198 
1199  if (fe_type.family == SCALAR)
1200  continue;
1201 
1202  // Per-subdomain variables don't need to be projected on
1203  // elements where they're not active
1204  if (!variable.active_on_subdomain(elem->subdomain_id()))
1205  continue;
1206 
1207  const std::vector<dof_id_type> & dof_indices =
1208  context.get_dof_indices(var);
1209 
1210  // The number of DOFs on the element
1211  const unsigned int n_dofs =
1212  cast_int<unsigned int>(dof_indices.size());
1213 
1214  const unsigned int var_component =
1216 
1217  // Zero the interpolated values
1218  Ue.resize (n_dofs); Ue.zero();
1219 
1220  // If we're projecting from an old grid
1221 #ifdef LIBMESH_ENABLE_AMR
1222  if (f.is_grid_projection() &&
1223  // And either this is an unchanged element
1224  ((elem->refinement_flag() != Elem::JUST_REFINED &&
1228  // Or this is a low order monomial element which has merely
1229  // been h refined
1230  (fe_type.family == MONOMIAL &&
1231  fe_type.order == CONSTANT &&
1232  elem->p_level() == 0 &&
1235  )
1236  // then we can simply copy its old dof
1237  // values to new indices.
1238  {
1239  LOG_SCOPE ("copy_dofs", "GenericProjector");
1240 
1241  f.eval_old_dofs(context, var_component, Ue.get_values());
1242 
1243  action.insert(context, var, Ue);
1244 
1245  continue;
1246  }
1247 #endif // LIBMESH_ENABLE_AMR
1248 
1249  FEBase * fe = libmesh_nullptr;
1250  FEBase * side_fe = libmesh_nullptr;
1251  FEBase * edge_fe = libmesh_nullptr;
1252 
1253  context.get_element_fe( var, fe, dim );
1254  if (fe->get_fe_type().family == SCALAR)
1255  continue;
1256  if (dim > 1)
1257  context.get_side_fe( var, side_fe, dim );
1258  if (dim > 2)
1259  context.get_edge_fe( var, edge_fe );
1260 
1261  const FEContinuity cont = fe->get_continuity();
1262 
1263  std::vector<unsigned int> side_dofs;
1264 
1265  // Fixed vs. free DoFs on edge/face projections
1266  std::vector<char> dof_is_fixed(n_dofs, false); // bools
1267  std::vector<int> free_dof(n_dofs, 0);
1268 
1269  // The element type
1270  const ElemType elem_type = elem->type();
1271 
1272  // The number of nodes on the new element
1273  const unsigned int n_nodes = elem->n_nodes();
1274 
1275  START_LOG ("project_nodes","GenericProjector");
1276 
1277  // When interpolating C1 elements we need to know which
1278  // vertices were also parent vertices; we'll cache an
1279  // intermediate fact outside the nodes loop.
1280  int i_am_child = -1;
1281 #ifdef LIBMESH_ENABLE_AMR
1282  if (elem->parent())
1283  i_am_child = elem->parent()->which_child_am_i(elem);
1284 #endif
1285  // In general, we need a series of
1286  // projections to ensure a unique and continuous
1287  // solution. We start by interpolating nodes, then
1288  // hold those fixed and project edges, then
1289  // hold those fixed and project faces, then
1290  // hold those fixed and project interiors
1291  //
1292  // In the LAGRANGE case, we will save a lot of solution
1293  // evaluations (at a slight cost in accuracy) by simply
1294  // interpolating all nodes rather than projecting.
1295 
1296  // Interpolate vertex (or for LAGRANGE, all node) values first.
1297  unsigned int current_dof = 0;
1298  for (unsigned int n=0; n!= n_nodes; ++n)
1299  {
1300  // FIXME: this should go through the DofMap,
1301  // not duplicate dof_indices code badly!
1302  const unsigned int nc =
1303  FEInterface::n_dofs_at_node (dim, fe_type, elem_type, n);
1304 
1305  if (!elem->is_vertex(n) &&
1306  fe_type.family != LAGRANGE)
1307  {
1308  current_dof += nc;
1309  continue;
1310  }
1311 
1312  if (cont == DISCONTINUOUS)
1313  {
1314  libmesh_assert_equal_to (nc, 0);
1315  }
1316  else if (!nc)
1317  {
1318  // This should only occur for first-order LAGRANGE
1319  // FE on non-vertices of higher-order elements
1320  libmesh_assert (!elem->is_vertex(n));
1321  libmesh_assert_equal_to(fe_type.family, LAGRANGE);
1322  }
1323  // Assume that C_ZERO elements have a single nodal
1324  // value shape function at vertices
1325  else if (cont == C_ZERO)
1326  {
1327  Ue(current_dof) = f.eval_at_node(context,
1328  var_component,
1329  dim,
1330  elem->node_ref(n),
1331  system.time);
1332  dof_is_fixed[current_dof] = true;
1333  current_dof++;
1334  }
1335  else if (cont == C_ONE)
1336  {
1337  const bool is_parent_vertex = (i_am_child == -1) ||
1338  elem->parent()->is_vertex_on_parent(i_am_child, n);
1339 
1340  // The hermite element vertex shape functions are weird
1341  if (fe_type.family == HERMITE)
1342  {
1343  Ue(current_dof) =
1344  f.eval_at_node(context,
1345  var_component,
1346  dim,
1347  elem->node_ref(n),
1348  system.time);
1349  dof_is_fixed[current_dof] = true;
1350  current_dof++;
1351  VectorValue<FValue> grad =
1352  is_parent_vertex ?
1353  g->eval_at_node(context,
1354  var_component,
1355  dim,
1356  elem->node_ref(n),
1357  system.time) :
1358  g->eval_at_point(context,
1359  var_component,
1360  elem->node_ref(n),
1361  system.time);
1362  // x derivative
1363  Ue(current_dof) = grad(0);
1364  dof_is_fixed[current_dof] = true;
1365  current_dof++;
1366  if (dim > 1)
1367  {
1368  // We'll finite difference mixed derivatives
1369  Real delta_x = TOLERANCE * elem->hmin();
1370 
1371  Point nxminus = elem->point(n),
1372  nxplus = elem->point(n);
1373  nxminus(0) -= delta_x;
1374  nxplus(0) += delta_x;
1375  VectorValue<FValue> gxminus =
1376  g->eval_at_point(context,
1377  var_component,
1378  nxminus,
1379  system.time);
1380  VectorValue<FValue> gxplus =
1381  g->eval_at_point(context,
1382  var_component,
1383  nxplus,
1384  system.time);
1385  // y derivative
1386  Ue(current_dof) = grad(1);
1387  dof_is_fixed[current_dof] = true;
1388  current_dof++;
1389  // xy derivative
1390  Ue(current_dof) = (gxplus(1) - gxminus(1))
1391  / 2. / delta_x;
1392  dof_is_fixed[current_dof] = true;
1393  current_dof++;
1394 
1395  if (dim > 2)
1396  {
1397  // z derivative
1398  Ue(current_dof) = grad(2);
1399  dof_is_fixed[current_dof] = true;
1400  current_dof++;
1401  // xz derivative
1402  Ue(current_dof) = (gxplus(2) - gxminus(2))
1403  / 2. / delta_x;
1404  dof_is_fixed[current_dof] = true;
1405  current_dof++;
1406  // We need new points for yz
1407  Point nyminus = elem->point(n),
1408  nyplus = elem->point(n);
1409  nyminus(1) -= delta_x;
1410  nyplus(1) += delta_x;
1411  VectorValue<FValue> gyminus =
1412  g->eval_at_point(context,
1413  var_component,
1414  nyminus,
1415  system.time);
1416  VectorValue<FValue> gyplus =
1417  g->eval_at_point(context,
1418  var_component,
1419  nyplus,
1420  system.time);
1421  // xz derivative
1422  Ue(current_dof) = (gyplus(2) - gyminus(2))
1423  / 2. / delta_x;
1424  dof_is_fixed[current_dof] = true;
1425  current_dof++;
1426  // Getting a 2nd order xyz is more tedious
1427  Point nxmym = elem->point(n),
1428  nxmyp = elem->point(n),
1429  nxpym = elem->point(n),
1430  nxpyp = elem->point(n);
1431  nxmym(0) -= delta_x;
1432  nxmym(1) -= delta_x;
1433  nxmyp(0) -= delta_x;
1434  nxmyp(1) += delta_x;
1435  nxpym(0) += delta_x;
1436  nxpym(1) -= delta_x;
1437  nxpyp(0) += delta_x;
1438  nxpyp(1) += delta_x;
1439  VectorValue<FValue> gxmym =
1440  g->eval_at_point(context,
1441  var_component,
1442  nxmym,
1443  system.time);
1444  VectorValue<FValue> gxmyp =
1445  g->eval_at_point(context,
1446  var_component,
1447  nxmyp,
1448  system.time);
1449  VectorValue<FValue> gxpym =
1450  g->eval_at_point(context,
1451  var_component,
1452  nxpym,
1453  system.time);
1454  VectorValue<FValue> gxpyp =
1455  g->eval_at_point(context,
1456  var_component,
1457  nxpyp,
1458  system.time);
1459  FValue gxzplus = (gxpyp(2) - gxmyp(2))
1460  / 2. / delta_x;
1461  FValue gxzminus = (gxpym(2) - gxmym(2))
1462  / 2. / delta_x;
1463  // xyz derivative
1464  Ue(current_dof) = (gxzplus - gxzminus)
1465  / 2. / delta_x;
1466  dof_is_fixed[current_dof] = true;
1467  current_dof++;
1468  }
1469  }
1470  }
1471  else
1472  {
1473  // Assume that other C_ONE elements have a single nodal
1474  // value shape function and nodal gradient component
1475  // shape functions
1476  libmesh_assert_equal_to (nc, 1 + dim);
1477  Ue(current_dof) = f.eval_at_node(context,
1478  var_component,
1479  dim,
1480  elem->node_ref(n),
1481  system.time);
1482  dof_is_fixed[current_dof] = true;
1483  current_dof++;
1484  VectorValue<FValue> grad =
1485  is_parent_vertex ?
1486  g->eval_at_node(context,
1487  var_component,
1488  dim,
1489  elem->node_ref(n),
1490  system.time) :
1491  g->eval_at_point(context,
1492  var_component,
1493  elem->node_ref(n),
1494  system.time);
1495  for (unsigned int i=0; i!= dim; ++i)
1496  {
1497  Ue(current_dof) = grad(i);
1498  dof_is_fixed[current_dof] = true;
1499  current_dof++;
1500  }
1501  }
1502  }
1503  else
1504  libmesh_error_msg("Unknown continuity " << cont);
1505  }
1506 
1507  STOP_LOG ("project_nodes","GenericProjector");
1508 
1509  START_LOG ("project_edges","GenericProjector");
1510 
1511  // In 3D with non-LAGRANGE, project any edge values next
1512  if (dim > 2 &&
1513  cont != DISCONTINUOUS &&
1514  (fe_type.family != LAGRANGE
1515 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1516  || elem->infinite()
1517 #endif
1518  ))
1519  {
1520  // If we're JUST_COARSENED we'll need a custom
1521  // evaluation, not just the standard edge FE
1522  const std::vector<Point> & xyz_values =
1523 #ifdef LIBMESH_ENABLE_AMR
1524  (elem->refinement_flag() == Elem::JUST_COARSENED) ?
1525  fe->get_xyz() :
1526 #endif
1527  edge_fe->get_xyz();
1528  const std::vector<Real> & JxW =
1529 #ifdef LIBMESH_ENABLE_AMR
1530  (elem->refinement_flag() == Elem::JUST_COARSENED) ?
1531  fe->get_JxW() :
1532 #endif
1533  edge_fe->get_JxW();
1534 
1535  const std::vector<std::vector<Real>> & phi =
1536 #ifdef LIBMESH_ENABLE_AMR
1537  (elem->refinement_flag() == Elem::JUST_COARSENED) ?
1538  fe->get_phi() :
1539 #endif
1540  edge_fe->get_phi();
1541  const std::vector<std::vector<RealGradient>> * dphi = libmesh_nullptr;
1542  if (cont == C_ONE)
1543  dphi =
1544 #ifdef LIBMESH_ENABLE_AMR
1545  (elem->refinement_flag() == Elem::JUST_COARSENED) ?
1546  &(fe->get_dphi()) :
1547 #endif
1548  &(edge_fe->get_dphi());
1549 
1550  for (auto e : elem->edge_index_range())
1551  {
1552  context.edge = e;
1553 
1554 #ifdef LIBMESH_ENABLE_AMR
1555  if (elem->refinement_flag() == Elem::JUST_COARSENED)
1556  {
1557  std::vector<Point> fine_points;
1558 
1559  UniquePtr<FEBase> fine_fe
1560  (FEBase::build (dim, base_fe_type));
1561 
1562  UniquePtr<QBase> qrule
1563  (base_fe_type.default_quadrature_rule(1));
1564  fine_fe->attach_quadrature_rule(qrule.get());
1565 
1566  const std::vector<Point> & child_xyz =
1567  fine_fe->get_xyz();
1568 
1569  for (unsigned int c = 0, nc = elem->n_children();
1570  c != nc; ++c)
1571  {
1572  if (!elem->is_child_on_edge(c, e))
1573  continue;
1574 
1575  fine_fe->edge_reinit(elem->child_ptr(c), e);
1576  fine_points.insert(fine_points.end(),
1577  child_xyz.begin(),
1578  child_xyz.end());
1579  }
1580 
1581  std::vector<Point> fine_qp;
1582  FEInterface::inverse_map (dim, base_fe_type, elem,
1583  fine_points, fine_qp);
1584 
1585  context.elem_fe_reinit(&fine_qp);
1586  }
1587  else
1588 #endif // LIBMESH_ENABLE_AMR
1589  context.edge_fe_reinit();
1590 
1591  const unsigned int n_qp = xyz_values.size();
1592 
1593  FEInterface::dofs_on_edge(elem, dim, base_fe_type,
1594  e, side_dofs);
1595 
1596  // Some edge dofs are on nodes and already
1597  // fixed, others are free to calculate
1598  unsigned int free_dofs = 0;
1599  for (std::size_t i=0; i != side_dofs.size(); ++i)
1600  if (!dof_is_fixed[side_dofs[i]])
1601  free_dof[free_dofs++] = i;
1602 
1603  // There may be nothing to project
1604  if (!free_dofs)
1605  continue;
1606 
1607  Ke.resize (free_dofs, free_dofs); Ke.zero();
1608  Fe.resize (free_dofs); Fe.zero();
1609  // The new edge coefficients
1610  DenseVector<FValue> Uedge(free_dofs);
1611 
1612  // Loop over the quadrature points
1613  for (unsigned int qp=0; qp<n_qp; qp++)
1614  {
1615  // solution at the quadrature point
1616  FValue fineval = f.eval_at_point(context,
1617  var_component,
1618  xyz_values[qp],
1619  system.time);
1620  // solution grad at the quadrature point
1621  VectorValue<FValue> finegrad;
1622  if (cont == C_ONE)
1623  finegrad = g->eval_at_point(context,
1624  var_component,
1625  xyz_values[qp],
1626  system.time);
1627 
1628  // Form edge projection matrix
1629  for (std::size_t sidei=0, freei=0;
1630  sidei != side_dofs.size(); ++sidei)
1631  {
1632  unsigned int i = side_dofs[sidei];
1633  // fixed DoFs aren't test functions
1634  if (dof_is_fixed[i])
1635  continue;
1636  for (std::size_t sidej=0, freej=0;
1637  sidej != side_dofs.size(); ++sidej)
1638  {
1639  unsigned int j = side_dofs[sidej];
1640  if (dof_is_fixed[j])
1641  Fe(freei) -= phi[i][qp] * phi[j][qp] *
1642  JxW[qp] * Ue(j);
1643  else
1644  Ke(freei,freej) += phi[i][qp] *
1645  phi[j][qp] * JxW[qp];
1646  if (cont == C_ONE)
1647  {
1648  if (dof_is_fixed[j])
1649  Fe(freei) -= ( (*dphi)[i][qp] *
1650  (*dphi)[j][qp] ) *
1651  JxW[qp] * Ue(j);
1652  else
1653  Ke(freei,freej) += ( (*dphi)[i][qp] *
1654  (*dphi)[j][qp] )
1655  * JxW[qp];
1656  }
1657  if (!dof_is_fixed[j])
1658  freej++;
1659  }
1660  Fe(freei) += phi[i][qp] * fineval * JxW[qp];
1661  if (cont == C_ONE)
1662  Fe(freei) += (finegrad * (*dphi)[i][qp] ) *
1663  JxW[qp];
1664  freei++;
1665  }
1666  }
1667 
1668  Ke.cholesky_solve(Fe, Uedge);
1669 
1670  // Transfer new edge solutions to element
1671  for (unsigned int i=0; i != free_dofs; ++i)
1672  {
1673  FValue & ui = Ue(side_dofs[free_dof[i]]);
1675  std::abs(ui - Uedge(i)) < TOLERANCE);
1676  ui = Uedge(i);
1677  dof_is_fixed[side_dofs[free_dof[i]]] = true;
1678  }
1679  }
1680  } // end if (dim > 2, !DISCONTINUOUS, !LAGRANGE)
1681 
1682  STOP_LOG ("project_edges","GenericProjector");
1683 
1684  START_LOG ("project_sides","GenericProjector");
1685 
1686  // With non-LAGRANGE, project any side values (edges in 2D,
1687  // faces in 3D) next.
1688  if (dim > 1 &&
1689  cont != DISCONTINUOUS &&
1690  (fe_type.family != LAGRANGE
1691 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1692  || elem->infinite()
1693 #endif
1694  ))
1695  {
1696  // If we're JUST_COARSENED we'll need a custom
1697  // evaluation, not just the standard side FE
1698  const std::vector<Point> & xyz_values =
1699 #ifdef LIBMESH_ENABLE_AMR
1700  (elem->refinement_flag() == Elem::JUST_COARSENED) ?
1701  fe->get_xyz() :
1702 #endif // LIBMESH_ENABLE_AMR
1703  side_fe->get_xyz();
1704  const std::vector<Real> & JxW =
1705 #ifdef LIBMESH_ENABLE_AMR
1706  (elem->refinement_flag() == Elem::JUST_COARSENED) ?
1707  fe->get_JxW() :
1708 #endif // LIBMESH_ENABLE_AMR
1709  side_fe->get_JxW();
1710  const std::vector<std::vector<Real>> & phi =
1711 #ifdef LIBMESH_ENABLE_AMR
1712  (elem->refinement_flag() == Elem::JUST_COARSENED) ?
1713  fe->get_phi() :
1714 #endif // LIBMESH_ENABLE_AMR
1715  side_fe->get_phi();
1716  const std::vector<std::vector<RealGradient>> * dphi = libmesh_nullptr;
1717  if (cont == C_ONE)
1718  dphi =
1719 #ifdef LIBMESH_ENABLE_AMR
1720  (elem->refinement_flag() == Elem::JUST_COARSENED) ? &(fe->get_dphi()) :
1721 #endif // LIBMESH_ENABLE_AMR
1722  &(side_fe->get_dphi());
1723 
1724  for (auto s : elem->side_index_range())
1725  {
1726  FEInterface::dofs_on_side(elem, dim, base_fe_type,
1727  s, side_dofs);
1728 
1729  // Some side dofs are on nodes/edges and already
1730  // fixed, others are free to calculate
1731  unsigned int free_dofs = 0;
1732  for (std::size_t i=0; i != side_dofs.size(); ++i)
1733  if (!dof_is_fixed[side_dofs[i]])
1734  free_dof[free_dofs++] = i;
1735 
1736  // There may be nothing to project
1737  if (!free_dofs)
1738  continue;
1739 
1740  Ke.resize (free_dofs, free_dofs); Ke.zero();
1741  Fe.resize (free_dofs); Fe.zero();
1742  // The new side coefficients
1743  DenseVector<FValue> Uside(free_dofs);
1744 
1745  context.side = s;
1746 
1747 #ifdef LIBMESH_ENABLE_AMR
1748  if (elem->refinement_flag() == Elem::JUST_COARSENED)
1749  {
1750  std::vector<Point> fine_points;
1751 
1752  UniquePtr<FEBase> fine_fe
1753  (FEBase::build (dim, base_fe_type));
1754 
1755  UniquePtr<QBase> qrule
1756  (base_fe_type.default_quadrature_rule(dim-1));
1757  fine_fe->attach_quadrature_rule(qrule.get());
1758 
1759  const std::vector<Point> & child_xyz =
1760  fine_fe->get_xyz();
1761 
1762  for (unsigned int c = 0, nc = elem->n_children();
1763  c != nc; ++c)
1764  {
1765  if (!elem->is_child_on_side(c, s))
1766  continue;
1767 
1768  fine_fe->reinit(elem->child_ptr(c), s);
1769  fine_points.insert(fine_points.end(),
1770  child_xyz.begin(),
1771  child_xyz.end());
1772  }
1773 
1774  std::vector<Point> fine_qp;
1775  FEInterface::inverse_map (dim, base_fe_type, elem,
1776  fine_points, fine_qp);
1777 
1778  context.elem_fe_reinit(&fine_qp);
1779  }
1780  else
1781 #endif // LIBMESH_ENABLE_AMR
1782  context.side_fe_reinit();
1783 
1784  const unsigned int n_qp = xyz_values.size();
1785 
1786  // Loop over the quadrature points
1787  for (unsigned int qp=0; qp<n_qp; qp++)
1788  {
1789  // solution at the quadrature point
1790  FValue fineval = f.eval_at_point(context,
1791  var_component,
1792  xyz_values[qp],
1793  system.time);
1794  // solution grad at the quadrature point
1795  VectorValue<FValue> finegrad;
1796  if (cont == C_ONE)
1797  finegrad = g->eval_at_point(context,
1798  var_component,
1799  xyz_values[qp],
1800  system.time);
1801 
1802  // Form side projection matrix
1803  for (std::size_t sidei=0, freei=0;
1804  sidei != side_dofs.size(); ++sidei)
1805  {
1806  unsigned int i = side_dofs[sidei];
1807  // fixed DoFs aren't test functions
1808  if (dof_is_fixed[i])
1809  continue;
1810  for (std::size_t sidej=0, freej=0;
1811  sidej != side_dofs.size(); ++sidej)
1812  {
1813  unsigned int j = side_dofs[sidej];
1814  if (dof_is_fixed[j])
1815  Fe(freei) -= phi[i][qp] * phi[j][qp] *
1816  JxW[qp] * Ue(j);
1817  else
1818  Ke(freei,freej) += phi[i][qp] *
1819  phi[j][qp] * JxW[qp];
1820  if (cont == C_ONE)
1821  {
1822  if (dof_is_fixed[j])
1823  Fe(freei) -= ( (*dphi)[i][qp] *
1824  (*dphi)[j][qp] ) *
1825  JxW[qp] * Ue(j);
1826  else
1827  Ke(freei,freej) += ( (*dphi)[i][qp] *
1828  (*dphi)[j][qp] )
1829  * JxW[qp];
1830  }
1831  if (!dof_is_fixed[j])
1832  freej++;
1833  }
1834  Fe(freei) += (fineval * phi[i][qp]) * JxW[qp];
1835  if (cont == C_ONE)
1836  Fe(freei) += (finegrad * (*dphi)[i][qp] ) *
1837  JxW[qp];
1838  freei++;
1839  }
1840  }
1841 
1842  Ke.cholesky_solve(Fe, Uside);
1843 
1844  // Transfer new side solutions to element
1845  for (unsigned int i=0; i != free_dofs; ++i)
1846  {
1847  FValue & ui = Ue(side_dofs[free_dof[i]]);
1849  std::abs(ui - Uside(i)) < TOLERANCE);
1850  ui = Uside(i);
1851  dof_is_fixed[side_dofs[free_dof[i]]] = true;
1852  }
1853  }
1854  } // end if (dim > 1, !DISCONTINUOUS, !LAGRANGE)
1855 
1856  STOP_LOG ("project_sides","GenericProjector");
1857 
1858  START_LOG ("project_interior","GenericProjector");
1859 
1860  // Project the interior values, finally
1861 
1862  // Some interior dofs are on nodes/edges/sides and
1863  // already fixed, others are free to calculate
1864  unsigned int free_dofs = 0;
1865  for (unsigned int i=0; i != n_dofs; ++i)
1866  if (!dof_is_fixed[i])
1867  free_dof[free_dofs++] = i;
1868 
1869  // Project any remaining (interior) dofs in the non-LAGRANGE
1870  // case.
1871  if (free_dofs &&
1872  (fe_type.family != LAGRANGE
1873 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1874  || elem->infinite()
1875 #endif
1876  ))
1877  {
1878  const std::vector<Point> & xyz_values = fe->get_xyz();
1879  const std::vector<Real> & JxW = fe->get_JxW();
1880 
1881  const std::vector<std::vector<Real>> & phi = fe->get_phi();
1882  const std::vector<std::vector<RealGradient>> * dphi = libmesh_nullptr;
1883  if (cont == C_ONE)
1884  dphi = &(fe->get_dphi());
1885 
1886 #ifdef LIBMESH_ENABLE_AMR
1887  if (elem->refinement_flag() == Elem::JUST_COARSENED)
1888  {
1889  std::vector<Point> fine_points;
1890 
1891  UniquePtr<FEBase> fine_fe
1892  (FEBase::build (dim, base_fe_type));
1893 
1894  UniquePtr<QBase> qrule
1895  (base_fe_type.default_quadrature_rule(dim));
1896  fine_fe->attach_quadrature_rule(qrule.get());
1897 
1898  const std::vector<Point> & child_xyz =
1899  fine_fe->get_xyz();
1900 
1901  for (auto & child : elem->child_ref_range())
1902  {
1903  fine_fe->reinit(&child);
1904  fine_points.insert(fine_points.end(),
1905  child_xyz.begin(),
1906  child_xyz.end());
1907  }
1908 
1909  std::vector<Point> fine_qp;
1910  FEInterface::inverse_map (dim, base_fe_type, elem,
1911  fine_points, fine_qp);
1912 
1913  context.elem_fe_reinit(&fine_qp);
1914  }
1915  else
1916 #endif // LIBMESH_ENABLE_AMR
1917  context.elem_fe_reinit();
1918 
1919  const unsigned int n_qp = xyz_values.size();
1920 
1921  Ke.resize (free_dofs, free_dofs); Ke.zero();
1922  Fe.resize (free_dofs); Fe.zero();
1923  // The new interior coefficients
1924  DenseVector<FValue> Uint(free_dofs);
1925 
1926  // Loop over the quadrature points
1927  for (unsigned int qp=0; qp<n_qp; qp++)
1928  {
1929  // solution at the quadrature point
1930  FValue fineval = f.eval_at_point(context,
1931  var_component,
1932  xyz_values[qp],
1933  system.time);
1934  // solution grad at the quadrature point
1935  VectorValue<FValue> finegrad;
1936  if (cont == C_ONE)
1937  finegrad = g->eval_at_point(context,
1938  var_component,
1939  xyz_values[qp],
1940  system.time);
1941 
1942  // Form interior projection matrix
1943  for (unsigned int i=0, freei=0; i != n_dofs; ++i)
1944  {
1945  // fixed DoFs aren't test functions
1946  if (dof_is_fixed[i])
1947  continue;
1948  for (unsigned int j=0, freej=0; j != n_dofs; ++j)
1949  {
1950  if (dof_is_fixed[j])
1951  Fe(freei) -= phi[i][qp] * phi[j][qp] * JxW[qp]
1952  * Ue(j);
1953  else
1954  Ke(freei,freej) += phi[i][qp] * phi[j][qp] *
1955  JxW[qp];
1956  if (cont == C_ONE)
1957  {
1958  if (dof_is_fixed[j])
1959  Fe(freei) -= ( (*dphi)[i][qp] *
1960  (*dphi)[j][qp] ) * JxW[qp] *
1961  Ue(j);
1962  else
1963  Ke(freei,freej) += ( (*dphi)[i][qp] *
1964  (*dphi)[j][qp] ) *
1965  JxW[qp];
1966  }
1967  if (!dof_is_fixed[j])
1968  freej++;
1969  }
1970  Fe(freei) += phi[i][qp] * fineval * JxW[qp];
1971  if (cont == C_ONE)
1972  Fe(freei) += (finegrad * (*dphi)[i][qp] ) * JxW[qp];
1973  freei++;
1974  }
1975  }
1976  Ke.cholesky_solve(Fe, Uint);
1977 
1978  // Transfer new interior solutions to element
1979  for (unsigned int i=0; i != free_dofs; ++i)
1980  {
1981  FValue & ui = Ue(free_dof[i]);
1983  std::abs(ui - Uint(i)) < TOLERANCE);
1984  ui = Uint(i);
1985  dof_is_fixed[free_dof[i]] = true;
1986  }
1987 
1988  } // if there are free interior dofs
1989 
1990  STOP_LOG ("project_interior","GenericProjector");
1991 
1992  // Make sure every DoF got reached!
1993  for (unsigned int i=0; i != n_dofs; ++i)
1994  libmesh_assert(dof_is_fixed[i]);
1995 
1996  action.insert(context, var, Ue);
1997  } // end variables loop
1998  } // end elements loop
1999 }
2000 
2001 
2002 
2003 #ifdef LIBMESH_ENABLE_AMR
2005 {
2006  // Sort the send list. After this duplicated
2007  // elements will be adjacent in the vector
2008  std::sort(this->send_list.begin(),
2009  this->send_list.end());
2010 
2011  // Now use std::unique to remove duplicate entries
2012  std::vector<dof_id_type>::iterator new_end =
2013  std::unique (this->send_list.begin(),
2014  this->send_list.end());
2015 
2016  // Remove the end of the send_list. Use the "swap trick"
2017  // from Effective STL
2018  std::vector<dof_id_type>
2019  (this->send_list.begin(), new_end).swap (this->send_list);
2020 }
2021 
2022 
2023 
2025 {
2026  // The DofMap for this system
2027  const DofMap & dof_map = system.get_dof_map();
2028 
2029  const dof_id_type first_old_dof = dof_map.first_old_dof();
2030  const dof_id_type end_old_dof = dof_map.end_old_dof();
2031 
2032  // We can handle all the variables at once.
2033  // The old global DOF indices
2034  std::vector<dof_id_type> di;
2035 
2036  // Iterate over the elements in the range
2037  for (ConstElemRange::const_iterator elem_it=range.begin(); elem_it != range.end(); ++elem_it)
2038  {
2039  const Elem * elem = *elem_it;
2040  // If this element doesn't have an old_dof_object with dofs for the
2041  // current system, then it must be newly added, so the user
2042  // is responsible for setting the new dofs.
2043 
2044  // ... but we need a better way to test for that; the code
2045  // below breaks on any FE type for which the elem stores no
2046  // dofs.
2047  // if (!elem->old_dof_object || !elem->old_dof_object->has_dofs(system.number()))
2048  // continue;
2049 
2050  // Examining refinement flags instead should distinguish
2051  // between refinement-added and user-added elements lacking
2052  // old_dof_object
2053  if (!elem->old_dof_object &&
2054  elem->refinement_flag() != Elem::JUST_REFINED &&
2056  continue;
2057 
2058  const Elem * parent = elem->parent();
2059 
2060  if (elem->refinement_flag() == Elem::JUST_REFINED)
2061  {
2062  libmesh_assert(parent);
2063 
2064  // We used to hack_p_level here, but that wasn't thread-safe
2065  // so now we take p refinement flags into account in
2066  // old_dof_indices
2067 
2068  dof_map.old_dof_indices (parent, di);
2069 
2070  for (auto & node : elem->node_ref_range())
2071  {
2072  const DofObject * old_dofs = node.old_dof_object;
2073 
2074  if (old_dofs)
2075  {
2076  const unsigned int sysnum = system.number();
2077  const unsigned int nv = old_dofs->n_vars(sysnum);
2078  for (unsigned int v=0; v != nv; ++v)
2079  {
2080  const unsigned int nc =
2081  old_dofs->n_comp(sysnum, v);
2082  for (unsigned int c=0; c != nc; ++c)
2083  di.push_back
2084  (old_dofs->dof_number(sysnum, v, c));
2085  }
2086  }
2087  }
2088 
2089  std::sort(di.begin(), di.end());
2090  std::vector<dof_id_type>::iterator new_end =
2091  std::unique(di.begin(), di.end());
2092  std::vector<dof_id_type>(di.begin(), new_end).swap(di);
2093  }
2094  else if (elem->refinement_flag() == Elem::JUST_COARSENED)
2095  {
2096  std::vector<dof_id_type> di_child;
2097  di.clear();
2098  for (auto & child : elem->child_ref_range())
2099  {
2100  dof_map.old_dof_indices (&child, di_child);
2101  di.insert(di.end(), di_child.begin(), di_child.end());
2102  }
2103  }
2104  else
2105  dof_map.old_dof_indices (elem, di);
2106 
2107  for (std::size_t i=0; i != di.size(); ++i)
2108  if (di[i] < first_old_dof || di[i] >= end_old_dof)
2109  this->send_list.push_back(di[i]);
2110  } // end elem loop
2111 }
2112 
2113 
2114 
2116 {
2117  // Joining simply requires I add the dof indices from the other object
2118  this->send_list.insert(this->send_list.end(),
2119  other.send_list.begin(),
2120  other.send_list.end());
2121 }
2122 #endif // LIBMESH_ENABLE_AMR
2123 
2124 
2125 
2127 {
2128  // We need data to project
2129  libmesh_assert(f.get());
2130 
2138  // The dimensionality of the current mesh
2139  const unsigned int dim = system.get_mesh().mesh_dimension();
2140 
2141  // The DofMap for this system
2142  const DofMap & dof_map = system.get_dof_map();
2143 
2144  // Boundary info for the current mesh
2145  const BoundaryInfo & boundary_info =
2147 
2148  // The element matrix and RHS for projections.
2149  // Note that Ke is always real-valued, whereas
2150  // Fe may be complex valued if complex number
2151  // support is enabled
2152  DenseMatrix<Real> Ke;
2154  // The new element coefficients
2156 
2157 
2158  // Loop over all the variables we've been requested to project
2159  for (std::size_t v=0; v!=variables.size(); v++)
2160  {
2161  const unsigned int var = variables[v];
2162 
2163  const Variable & variable = dof_map.variable(var);
2164 
2165  const FEType & fe_type = variable.type();
2166 
2167  if (fe_type.family == SCALAR)
2168  continue;
2169 
2170  const unsigned int var_component =
2172 
2173  // Get FE objects of the appropriate type
2174  UniquePtr<FEBase> fe (FEBase::build(dim, fe_type));
2175 
2176  // Prepare variables for projection
2177  UniquePtr<QBase> qedgerule (fe_type.default_quadrature_rule(1));
2178  UniquePtr<QBase> qsiderule (fe_type.default_quadrature_rule(dim-1));
2179 
2180  // The values of the shape functions at the quadrature
2181  // points
2182  const std::vector<std::vector<Real>> & phi = fe->get_phi();
2183 
2184  // The gradients of the shape functions at the quadrature
2185  // points on the child element.
2186  const std::vector<std::vector<RealGradient>> * dphi = libmesh_nullptr;
2187 
2188  const FEContinuity cont = fe->get_continuity();
2189 
2190  if (cont == C_ONE)
2191  {
2192  // We'll need gradient data for a C1 projection
2193  libmesh_assert(g.get());
2194 
2195  const std::vector<std::vector<RealGradient>> &
2196  ref_dphi = fe->get_dphi();
2197  dphi = &ref_dphi;
2198  }
2199 
2200  // The Jacobian * quadrature weight at the quadrature points
2201  const std::vector<Real> & JxW =
2202  fe->get_JxW();
2203 
2204  // The XYZ locations of the quadrature points
2205  const std::vector<Point> & xyz_values =
2206  fe->get_xyz();
2207 
2208  // The global DOF indices
2209  std::vector<dof_id_type> dof_indices;
2210  // Side/edge DOF indices
2211  std::vector<unsigned int> side_dofs;
2212 
2213  // Container to catch IDs passed back from BoundaryInfo.
2214  std::vector<boundary_id_type> bc_ids;
2215 
2216  // Iterate over all the elements in the range
2217  for (ConstElemRange::const_iterator elem_it=range.begin(); elem_it != range.end(); ++elem_it)
2218  {
2219  const Elem * elem = *elem_it;
2220 
2221  // Per-subdomain variables don't need to be projected on
2222  // elements where they're not active
2223  if (!variable.active_on_subdomain(elem->subdomain_id()))
2224  continue;
2225 
2226  const unsigned short n_nodes = elem->n_nodes();
2227  const unsigned short n_edges = elem->n_edges();
2228  const unsigned short n_sides = elem->n_sides();
2229 
2230  // Find out which nodes, edges and sides are on a requested
2231  // boundary:
2232  std::vector<bool> is_boundary_node(n_nodes, false),
2233  is_boundary_edge(n_edges, false),
2234  is_boundary_side(n_sides, false);
2235  for (unsigned char s=0; s != n_sides; ++s)
2236  {
2237  // First see if this side has been requested
2238  boundary_info.boundary_ids (elem, s, bc_ids);
2239  bool do_this_side = false;
2240  for (std::vector<boundary_id_type>::iterator i=bc_ids.begin();
2241  i!=bc_ids.end(); ++i)
2242  if (b.count(*i))
2243  {
2244  do_this_side = true;
2245  break;
2246  }
2247  if (!do_this_side)
2248  continue;
2249 
2250  is_boundary_side[s] = true;
2251 
2252  // Then see what nodes and what edges are on it
2253  for (unsigned int n=0; n != n_nodes; ++n)
2254  if (elem->is_node_on_side(n,s))
2255  is_boundary_node[n] = true;
2256  for (unsigned int e=0; e != n_edges; ++e)
2257  if (elem->is_edge_on_side(e,s))
2258  is_boundary_edge[e] = true;
2259  }
2260 
2261  // Update the DOF indices for this element based on
2262  // the current mesh
2263  dof_map.dof_indices (elem, dof_indices, var);
2264 
2265  // The number of DOFs on the element
2266  const unsigned int n_dofs =
2267  cast_int<unsigned int>(dof_indices.size());
2268 
2269  // Fixed vs. free DoFs on edge/face projections
2270  std::vector<char> dof_is_fixed(n_dofs, false); // bools
2271  std::vector<int> free_dof(n_dofs, 0);
2272 
2273  // The element type
2274  const ElemType elem_type = elem->type();
2275 
2276  // Zero the interpolated values
2277  Ue.resize (n_dofs); Ue.zero();
2278 
2279  // In general, we need a series of
2280  // projections to ensure a unique and continuous
2281  // solution. We start by interpolating boundary nodes, then
2282  // hold those fixed and project boundary edges, then hold
2283  // those fixed and project boundary faces,
2284 
2285  // Interpolate node values first
2286  unsigned int current_dof = 0;
2287  for (unsigned short n = 0; n != n_nodes; ++n)
2288  {
2289  // FIXME: this should go through the DofMap,
2290  // not duplicate dof_indices code badly!
2291  const unsigned int nc =
2292  FEInterface::n_dofs_at_node (dim, fe_type, elem_type,
2293  n);
2294  if (!elem->is_vertex(n) || !is_boundary_node[n])
2295  {
2296  current_dof += nc;
2297  continue;
2298  }
2299  if (cont == DISCONTINUOUS)
2300  {
2301  libmesh_assert_equal_to (nc, 0);
2302  }
2303  // Assume that C_ZERO elements have a single nodal
2304  // value shape function
2305  else if (cont == C_ZERO)
2306  {
2307  libmesh_assert_equal_to (nc, 1);
2308  Ue(current_dof) = f->component(var_component,
2309  elem->point(n),
2310  system.time);
2311  dof_is_fixed[current_dof] = true;
2312  current_dof++;
2313  }
2314  // The hermite element vertex shape functions are weird
2315  else if (fe_type.family == HERMITE)
2316  {
2317  Ue(current_dof) = f->component(var_component,
2318  elem->point(n),
2319  system.time);
2320  dof_is_fixed[current_dof] = true;
2321  current_dof++;
2322  Gradient grad = g->component(var_component,
2323  elem->point(n),
2324  system.time);
2325  // x derivative
2326  Ue(current_dof) = grad(0);
2327  dof_is_fixed[current_dof] = true;
2328  current_dof++;
2329  if (dim > 1)
2330  {
2331  // We'll finite difference mixed derivatives
2332  Point nxminus = elem->point(n),
2333  nxplus = elem->point(n);
2334  nxminus(0) -= TOLERANCE;
2335  nxplus(0) += TOLERANCE;
2336  Gradient gxminus = g->component(var_component,
2337  nxminus,
2338  system.time);
2339  Gradient gxplus = g->component(var_component,
2340  nxplus,
2341  system.time);
2342  // y derivative
2343  Ue(current_dof) = grad(1);
2344  dof_is_fixed[current_dof] = true;
2345  current_dof++;
2346  // xy derivative
2347  Ue(current_dof) = (gxplus(1) - gxminus(1))
2348  / 2. / TOLERANCE;
2349  dof_is_fixed[current_dof] = true;
2350  current_dof++;
2351 
2352  if (dim > 2)
2353  {
2354  // z derivative
2355  Ue(current_dof) = grad(2);
2356  dof_is_fixed[current_dof] = true;
2357  current_dof++;
2358  // xz derivative
2359  Ue(current_dof) = (gxplus(2) - gxminus(2))
2360  / 2. / TOLERANCE;
2361  dof_is_fixed[current_dof] = true;
2362  current_dof++;
2363  // We need new points for yz
2364  Point nyminus = elem->point(n),
2365  nyplus = elem->point(n);
2366  nyminus(1) -= TOLERANCE;
2367  nyplus(1) += TOLERANCE;
2368  Gradient gyminus = g->component(var_component,
2369  nyminus,
2370  system.time);
2371  Gradient gyplus = g->component(var_component,
2372  nyplus,
2373  system.time);
2374  // xz derivative
2375  Ue(current_dof) = (gyplus(2) - gyminus(2))
2376  / 2. / TOLERANCE;
2377  dof_is_fixed[current_dof] = true;
2378  current_dof++;
2379  // Getting a 2nd order xyz is more tedious
2380  Point nxmym = elem->point(n),
2381  nxmyp = elem->point(n),
2382  nxpym = elem->point(n),
2383  nxpyp = elem->point(n);
2384  nxmym(0) -= TOLERANCE;
2385  nxmym(1) -= TOLERANCE;
2386  nxmyp(0) -= TOLERANCE;
2387  nxmyp(1) += TOLERANCE;
2388  nxpym(0) += TOLERANCE;
2389  nxpym(1) -= TOLERANCE;
2390  nxpyp(0) += TOLERANCE;
2391  nxpyp(1) += TOLERANCE;
2392  Gradient gxmym = g->component(var_component,
2393  nxmym,
2394  system.time);
2395  Gradient gxmyp = g->component(var_component,
2396  nxmyp,
2397  system.time);
2398  Gradient gxpym = g->component(var_component,
2399  nxpym,
2400  system.time);
2401  Gradient gxpyp = g->component(var_component,
2402  nxpyp,
2403  system.time);
2404  Number gxzplus = (gxpyp(2) - gxmyp(2))
2405  / 2. / TOLERANCE;
2406  Number gxzminus = (gxpym(2) - gxmym(2))
2407  / 2. / TOLERANCE;
2408  // xyz derivative
2409  Ue(current_dof) = (gxzplus - gxzminus)
2410  / 2. / TOLERANCE;
2411  dof_is_fixed[current_dof] = true;
2412  current_dof++;
2413  }
2414  }
2415  }
2416  // Assume that other C_ONE elements have a single nodal
2417  // value shape function and nodal gradient component
2418  // shape functions
2419  else if (cont == C_ONE)
2420  {
2421  libmesh_assert_equal_to (nc, 1 + dim);
2422  Ue(current_dof) = f->component(var_component,
2423  elem->point(n),
2424  system.time);
2425  dof_is_fixed[current_dof] = true;
2426  current_dof++;
2427  Gradient grad = g->component(var_component,
2428  elem->point(n),
2429  system.time);
2430  for (unsigned int i=0; i!= dim; ++i)
2431  {
2432  Ue(current_dof) = grad(i);
2433  dof_is_fixed[current_dof] = true;
2434  current_dof++;
2435  }
2436  }
2437  else
2438  libmesh_error_msg("Unknown continuity " << cont);
2439  }
2440 
2441  // In 3D, project any edge values next
2442  if (dim > 2 && cont != DISCONTINUOUS)
2443  for (unsigned short e = 0; e != n_edges; ++e)
2444  {
2445  if (!is_boundary_edge[e])
2446  continue;
2447 
2448  FEInterface::dofs_on_edge(elem, dim, fe_type, e,
2449  side_dofs);
2450 
2451  // Some edge dofs are on nodes and already
2452  // fixed, others are free to calculate
2453  unsigned int free_dofs = 0;
2454  for (std::size_t i=0; i != side_dofs.size(); ++i)
2455  if (!dof_is_fixed[side_dofs[i]])
2456  free_dof[free_dofs++] = i;
2457 
2458  // There may be nothing to project
2459  if (!free_dofs)
2460  continue;
2461 
2462  Ke.resize (free_dofs, free_dofs); Ke.zero();
2463  Fe.resize (free_dofs); Fe.zero();
2464  // The new edge coefficients
2465  DenseVector<Number> Uedge(free_dofs);
2466 
2467  // Initialize FE data on the edge
2468  fe->attach_quadrature_rule (qedgerule.get());
2469  fe->edge_reinit (elem, e);
2470  const unsigned int n_qp = qedgerule->n_points();
2471 
2472  // Loop over the quadrature points
2473  for (unsigned int qp=0; qp<n_qp; qp++)
2474  {
2475  // solution at the quadrature point
2476  Number fineval = f->component(var_component,
2477  xyz_values[qp],
2478  system.time);
2479  // solution grad at the quadrature point
2480  Gradient finegrad;
2481  if (cont == C_ONE)
2482  finegrad = g->component(var_component,
2483  xyz_values[qp],
2484  system.time);
2485 
2486  // Form edge projection matrix
2487  for (std::size_t sidei=0, freei=0;
2488  sidei != side_dofs.size(); ++sidei)
2489  {
2490  unsigned int i = side_dofs[sidei];
2491  // fixed DoFs aren't test functions
2492  if (dof_is_fixed[i])
2493  continue;
2494  for (std::size_t sidej=0, freej=0;
2495  sidej != side_dofs.size(); ++sidej)
2496  {
2497  unsigned int j = side_dofs[sidej];
2498  if (dof_is_fixed[j])
2499  Fe(freei) -= phi[i][qp] * phi[j][qp] *
2500  JxW[qp] * Ue(j);
2501  else
2502  Ke(freei,freej) += phi[i][qp] *
2503  phi[j][qp] * JxW[qp];
2504  if (cont == C_ONE)
2505  {
2506  if (dof_is_fixed[j])
2507  Fe(freei) -= ((*dphi)[i][qp] *
2508  (*dphi)[j][qp]) *
2509  JxW[qp] * Ue(j);
2510  else
2511  Ke(freei,freej) += ((*dphi)[i][qp] *
2512  (*dphi)[j][qp])
2513  * JxW[qp];
2514  }
2515  if (!dof_is_fixed[j])
2516  freej++;
2517  }
2518  Fe(freei) += phi[i][qp] * fineval * JxW[qp];
2519  if (cont == C_ONE)
2520  Fe(freei) += (finegrad * (*dphi)[i][qp]) *
2521  JxW[qp];
2522  freei++;
2523  }
2524  }
2525 
2526  Ke.cholesky_solve(Fe, Uedge);
2527 
2528  // Transfer new edge solutions to element
2529  for (unsigned int i=0; i != free_dofs; ++i)
2530  {
2531  Number & ui = Ue(side_dofs[free_dof[i]]);
2533  std::abs(ui - Uedge(i)) < TOLERANCE);
2534  ui = Uedge(i);
2535  dof_is_fixed[side_dofs[free_dof[i]]] = true;
2536  }
2537  }
2538 
2539  // Project any side values (edges in 2D, faces in 3D)
2540  if (dim > 1 && cont != DISCONTINUOUS)
2541  for (unsigned short s = 0; s != n_sides; ++s)
2542  {
2543  if (!is_boundary_side[s])
2544  continue;
2545 
2546  FEInterface::dofs_on_side(elem, dim, fe_type, s,
2547  side_dofs);
2548 
2549  // Some side dofs are on nodes/edges and already
2550  // fixed, others are free to calculate
2551  unsigned int free_dofs = 0;
2552  for (std::size_t i=0; i != side_dofs.size(); ++i)
2553  if (!dof_is_fixed[side_dofs[i]])
2554  free_dof[free_dofs++] = i;
2555 
2556  // There may be nothing to project
2557  if (!free_dofs)
2558  continue;
2559 
2560  Ke.resize (free_dofs, free_dofs); Ke.zero();
2561  Fe.resize (free_dofs); Fe.zero();
2562  // The new side coefficients
2563  DenseVector<Number> Uside(free_dofs);
2564 
2565  // Initialize FE data on the side
2566  fe->attach_quadrature_rule (qsiderule.get());
2567  fe->reinit (elem, s);
2568  const unsigned int n_qp = qsiderule->n_points();
2569 
2570  // Loop over the quadrature points
2571  for (unsigned int qp=0; qp<n_qp; qp++)
2572  {
2573  // solution at the quadrature point
2574  Number fineval = f->component(var_component,
2575  xyz_values[qp],
2576  system.time);
2577  // solution grad at the quadrature point
2578  Gradient finegrad;
2579  if (cont == C_ONE)
2580  finegrad = g->component(var_component,
2581  xyz_values[qp],
2582  system.time);
2583 
2584  // Form side projection matrix
2585  for (std::size_t sidei=0, freei=0;
2586  sidei != side_dofs.size(); ++sidei)
2587  {
2588  unsigned int i = side_dofs[sidei];
2589  // fixed DoFs aren't test functions
2590  if (dof_is_fixed[i])
2591  continue;
2592  for (std::size_t sidej=0, freej=0;
2593  sidej != side_dofs.size(); ++sidej)
2594  {
2595  unsigned int j = side_dofs[sidej];
2596  if (dof_is_fixed[j])
2597  Fe(freei) -= phi[i][qp] * phi[j][qp] *
2598  JxW[qp] * Ue(j);
2599  else
2600  Ke(freei,freej) += phi[i][qp] *
2601  phi[j][qp] * JxW[qp];
2602  if (cont == C_ONE)
2603  {
2604  if (dof_is_fixed[j])
2605  Fe(freei) -= ((*dphi)[i][qp] *
2606  (*dphi)[j][qp]) *
2607  JxW[qp] * Ue(j);
2608  else
2609  Ke(freei,freej) += ((*dphi)[i][qp] *
2610  (*dphi)[j][qp])
2611  * JxW[qp];
2612  }
2613  if (!dof_is_fixed[j])
2614  freej++;
2615  }
2616  Fe(freei) += (fineval * phi[i][qp]) * JxW[qp];
2617  if (cont == C_ONE)
2618  Fe(freei) += (finegrad * (*dphi)[i][qp]) *
2619  JxW[qp];
2620  freei++;
2621  }
2622  }
2623 
2624  Ke.cholesky_solve(Fe, Uside);
2625 
2626  // Transfer new side solutions to element
2627  for (unsigned int i=0; i != free_dofs; ++i)
2628  {
2629  Number & ui = Ue(side_dofs[free_dof[i]]);
2631  std::abs(ui - Uside(i)) < TOLERANCE);
2632  ui = Uside(i);
2633  dof_is_fixed[side_dofs[free_dof[i]]] = true;
2634  }
2635  }
2636 
2637  const dof_id_type
2638  first = new_vector.first_local_index(),
2639  last = new_vector.last_local_index();
2640 
2641  // Lock the new_vector since it is shared among threads.
2642  {
2643  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
2644 
2645  for (unsigned int i = 0; i < n_dofs; i++)
2646  if (dof_is_fixed[i] &&
2647  (dof_indices[i] >= first) &&
2648  (dof_indices[i] < last))
2649  {
2650  new_vector.set(dof_indices[i], Ue(i));
2651  }
2652  }
2653  } // end elem loop
2654  } // end variables loop
2655 }
2656 
2657 
2658 } // namespace libMesh
This class implements the loops to other projection operations.
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...
Definition: fe_interface.C:510
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
virtual bool is_node_on_side(const unsigned int n, const unsigned int s) const =0
FEFamily family
The type of finite element.
Definition: fe_type.h:203
Real time
For time-dependent problems, this is the time t at the beginning of the current timestep.
Definition: system.h:1545
This class provides a wrapper with which to evaluate a (libMesh-style) function pointer in a Function...
dof_id_type end_old_dof(const processor_id_type proc) const
Definition: dof_map.h:600
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:117
double abs(double a)
const FEType & type() const
Definition: variable.h:119
OldSolutionValue(const OldSolutionValue &in)
A Node is like a Point, but with more information.
Definition: node.h:52
virtual void pre_fe_reinit(const System &, const Elem *e)
Reinitializes local data vectors/matrices on the current geometric element.
Definition: fem_context.C:1582
virtual numeric_index_type size() const =0
const NumericVector< Number > & old_solution
virtual bool is_vertex_on_parent(unsigned int c, unsigned int n) const
Definition: elem.C:3019
virtual numeric_index_type last_local_index() const =0
virtual ElemType type() const =0
BuildProjectionList(const System &system_in)
void get_edge_fe(unsigned int var, FEGenericBase< OutputShape > *&fe) const
Accessor for edge (3D only!) finite element object for variable var.
Definition: fem_context.h:1239
virtual void zero() libmesh_override
Set every element in the matrix to 0.
Definition: dense_matrix.h:792
This class provides the ability to map between arbitrary, user-defined strings and several data types...
Definition: parameters.h:63
void eval_old_dofs(const FEMContext &c, unsigned int var, std::vector< Output > &values)
virtual bool is_edge_on_side(const unsigned int e, const unsigned int s) const =0
unsigned int n_vars(const unsigned int s, const unsigned int vg) const
Definition: dof_object.h:745
virtual void init(const numeric_index_type, const numeric_index_type, const bool=false, const ParallelType=AUTOMATIC)=0
Change the dimension of the vector to N.
unsigned int p_level() const
Definition: elem.h:2422
virtual bool is_child_on_edge(const unsigned int c, const unsigned int e) const
Definition: elem.C:1665
void parallel_for(const Range &range, const Body &body)
Execute the provided function object in parallel on the specified range.
Definition: threads_none.h:73
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...
Definition: fe_interface.C:495
unsigned int dim
ImplicitSystem & sys
IntRange< unsigned short > side_index_range() const
Definition: elem.h:2083
virtual unsigned int n_edges() const =0
BoundaryProjectSolution(const std::set< boundary_id_type > &b_in, const std::vector< unsigned int > &variables_in, const System &system_in, FunctionBase< Number > *f_in, FunctionBase< Gradient > *g_in, const Parameters &parameters_in, NumericVector< Number > &new_v_in)
UniquePtr< FunctionBase< Number > > f
Dummy "splitting object" used to distinguish splitting constructors from copy constructors.
Definition: threads_none.h:63
ElemType
Defines an enum for geometric element types.
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:350
virtual void zero() libmesh_override
Set every element in the vector to 0.
Definition: dense_vector.h:374
const Elem * parent() const
Definition: elem.h:2346
std::vector< T > & get_values()
Definition: dense_vector.h:251
FEMFunctionWrapper(const FEMFunctionWrapper< Output > &fw)
void boundary_project_vector(const std::set< boundary_id_type > &b, const std::vector< unsigned int > &variables, NumericVector< Number > &new_vector, FunctionBase< Number > *f, FunctionBase< Gradient > *g=libmesh_nullptr, int is_adjoint=-1) const
Projects arbitrary boundary functions onto a vector of degree of freedom values for the current syste...
const std::vector< dof_id_type > & get_dof_indices() const
Accessor for element dof indices.
Definition: diff_context.h:366
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
This class implements projecting an arbitrary boundary function to the current mesh.
virtual Real hmin() const
Definition: elem.C:458
FEMFunctionWrapper(const FEMFunctionBase< Output > &f)
const class libmesh_nullptr_t libmesh_nullptr
dof_id_type dof_number(const unsigned int s, const unsigned int var, const unsigned int comp) const
Definition: dof_object.h:810
OrderWrapper order
The approximation order of the element.
Definition: fe_type.h:197
void get_element_fe(unsigned int var, FEGenericBase< OutputShape > *&fe) const
Accessor for interior finite element object for variable var for the largest dimension in the mesh...
Definition: fem_context.h:262
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:946
static const Real TOLERANCE
The StoredRange class defines a contiguous, divisible set of objects.
Definition: stored_range.h:52
void project_solution(FunctionBase< Number > *f, FunctionBase< Gradient > *g=libmesh_nullptr) const
Projects arbitrary functions onto the current solution.
void old_dof_indices(const Elem *const elem, std::vector< dof_id_type > &di, const unsigned int vn=libMesh::invalid_uint) const
After a mesh is refined and repartitioned it is possible that the _send_list will need to be augmente...
Definition: dof_map.C:2378
const Variable & variable(const unsigned int c) const
Definition: dof_map.h:1667
const unsigned int n_vars
Definition: tecplot_io.C:68
UniquePtr< FunctionBase< Gradient > > g
The libMesh namespace provides an interface to certain functionality in the library.
const std::vector< unsigned int > & variables
const_iterator begin() const
Beginning of the range.
Definition: stored_range.h:261
virtual unsigned int size() const libmesh_override
Definition: dense_vector.h:87
FEType get_fe_type() const
Definition: fe_abstract.h:421
NumericVector< Val > & target_vector
VectorSetAction(NumericVector< Val > &target_vec)
const Elem & get_elem() const
Accessor for current Elem object.
Definition: fem_context.h:871
virtual bool is_child_on_side(const unsigned int c, const unsigned int s) const =0
Output eval_at_node(const FEMContext &c, unsigned int i, unsigned int elem_dim, const Node &n, Real=0.)
IntRange< unsigned short > edge_index_range() const
Definition: elem.h:2074
SimpleRange< ChildRefIter > child_ref_range()
Returns a range with all children of a parent element, usable in range-based for loops.
Definition: elem.h:1699
Number fptr(const Point &p, const Parameters &, const std::string &libmesh_dbg_var(sys_name), const std::string &unknown_name)
Definition: projection.C:78
libmesh_assert(j)
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
virtual bool infinite() const =0
static UniquePtr< NumericVector< T > > build(const Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package())
Builds a NumericVector on the processors in communicator comm using the linear solver package specifi...
This class provides a wrapper with which to evaluate a (libMesh-style) function pointer in a Function...
virtual unsigned int n_nodes() const =0
std::vector< boundary_id_type > boundary_ids(const Node *node) const
BuildProjectionList(BuildProjectionList &other, Threads::split)
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:167
unsigned char edge
Current edge for edge_* to examine.
Definition: fem_context.h:982
Output eval_at_node(const FEMContext &c, unsigned int i, unsigned int, const Node &n, const Real time)
This class builds the send_list of old dof indices whose coefficients are needed to perform a project...
const std::set< unsigned char > & elem_dimensions() const
Definition: fem_context.h:913
const dof_id_type n_nodes
Definition: tecplot_io.C:67
const MeshBase & get_mesh() const
Definition: system.h:2014
static void get_shape_outputs(FEBase &fe)
spin_mutex spin_mtx
A convenient spin mutex object which can be used for obtaining locks.
Definition: threads.C:29
This class defines the notion of a variable in the system.
Definition: variable.h:49
const DofMap & get_dof_map() const
Definition: system.h:2030
const std::vector< Real > & get_JxW() const
Definition: fe_abstract.h:237
virtual void side_fe_reinit()
Reinitializes side FE objects on the current geometric element.
Definition: fem_context.C:1405
dof_id_type numeric_index_type
Definition: id_types.h:92
void eval_old_dofs(const FEMContext &, unsigned int, std::vector< Output >)
NumericVector< Number > & new_vector
This is the base class for classes which contain information related to any physical process that mig...
Definition: system.h:76
UniquePtr< FEMFunctionBase< Output > > _f
const std::vector< std::vector< OutputShape > > & get_phi() const
Definition: fe_base.h:208
static const Real out_of_elem_tol
The BoundaryInfo class contains information relevant to boundary conditions including storing faces...
Definition: boundary_info.h:56
bool active_on_subdomain(subdomain_id_type sid) const
Definition: variable.h:136
const Elem * child_ptr(unsigned int i) const
Definition: elem.h:2445
const Node & node_ref(const unsigned int i) const
Definition: elem.h:1896
void init_context(FEMContext &c)
This class provides all data required for a physics package (e.g.
Definition: fem_context.h:61
void init_context(FEMContext &c)
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)
Definition: fe_interface.C:569
RefinementState p_refinement_flag() const
Definition: elem.h:2521
subdomain_id_type subdomain_id() const
Definition: elem.h:1951
virtual unsigned int n_sides() const =0
void join(const BuildProjectionList &other)
const_iterator end() const
End of the range.
Definition: stored_range.h:266
const ProjectionAction & master_action
virtual void elem_fe_reinit(const std::vector< Point > *const pts=libmesh_nullptr)
Reinitializes interior FE objects on the current geometric element.
Definition: fem_context.C:1382
void insert(const FEMContext &c, unsigned int var_num, const DenseVector< Val > &Ue)
virtual unsigned int n_children() const =0
BoundaryProjectSolution(const BoundaryProjectSolution &in)
virtual void close()=0
Calls the NumericVector&#39;s internal assembly routines, ensuring that the values are consistent across ...
const std::vector< unsigned int > & variables
FEContinuity
defines an enum for finite element types to libmesh_assert a certain level (or type? Hcurl?) of continuity.
virtual numeric_index_type local_size() const =0
std::vector< object_type >::const_iterator const_iterator
Allows an StoredRange to behave like an STL container.
Definition: stored_range.h:58
virtual void edge_fe_reinit()
Reinitializes edge FE objects on the current geometric element.
Definition: fem_context.C:1425
SimpleRange< NodeRefIter > node_ref_range()
Returns a range with all nodes of an element, usable in range-based for loops.
Definition: elem.h:2047
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
unsigned int which_child_am_i(const Elem *e) const
Definition: elem.h:2487
virtual Real hmax() const
Definition: elem.C:475
void SCALAR_dof_indices(std::vector< dof_id_type > &di, const unsigned int vn, const bool old_dofs=false) const
Fills the vector di with the global degree of freedom indices corresponding to the SCALAR variable vn...
Definition: dof_map.C:2288
const std::vector< std::vector< OutputGradient > > & get_dphi() const
Definition: fe_base.h:216
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void operator()(const ConstElemRange &range)
DofObject * old_dof_object
This object on the last mesh.
Definition: dof_object.h:79
const Point & point(const unsigned int i) const
Definition: elem.h:1809
unsigned char side
Current side for side_* to examine.
Definition: fem_context.h:977
virtual bool is_vertex(const unsigned int i) const =0
void check_old_context(const FEMContext &c)
Output eval_at_point(const FEMContext &c, unsigned int i, const Point &p, Real=0.)
std::vector< dof_id_type > send_list
unsigned int number() const
Definition: system.h:2006
GenericProjector(const System &system_in, const FFunctor &f_in, const GFunctor *g_in, const ProjectionAction &act_in, const std::vector< unsigned int > &variables_in)
X_input swap(X_system)
Output eval_at_point(const FEMContext &c, unsigned int i, const Point &n, const Real time)
UniquePtr< QBase > default_quadrature_rule(const unsigned int dim, const int extraorder=0) const
Definition: fe_type.C:30
virtual void get(const std::vector< numeric_index_type > &index, T *values) const
Access multiple components at once.
unsigned int mesh_dimension() const
Definition: mesh_base.C:148
virtual Output component(const FEMContext &, unsigned int i, const Point &p, Real time=0.)
void resize(const unsigned int new_m, const unsigned int new_n)
Resize the matrix.
Definition: dense_matrix.h:776
ParallelType type() const
Gradient gptr(const Point &p, const Parameters &, const std::string &libmesh_dbg_var(sys_name), const std::string &unknown_name)
Definition: projection.C:93
virtual unsigned int dim() const =0
The DofObject defines an abstract base class for objects that have degrees of freedom associated with...
Definition: dof_object.h:51
Order
defines an enum for polynomial orders.
Definition: enum_order.h:32
RefinementState refinement_flag() const
Definition: elem.h:2505
void parallel_reduce(const Range &range, Body &body)
Execute the provided reduction operation in parallel on the specified range.
Definition: threads_none.h:101
void boundary_project_solution(const std::set< boundary_id_type > &b, const std::vector< unsigned int > &variables, FunctionBase< Number > *f, FunctionBase< Gradient > *g=libmesh_nullptr)
Projects arbitrary boundary functions onto a vector of degree of freedom values for the current syste...
virtual numeric_index_type first_local_index() const =0
virtual void clear()
Restores the NumericVector<T> to a pristine state.
unsigned int variable_scalar_number(const std::string &var, unsigned int component) const
Definition: system.h:2145
OldSolutionValue(const libMesh::System &sys_in, const NumericVector< Number > &old_sol)
bool check_old_context(const FEMContext &c, const Point &p)
virtual FEContinuity get_continuity() const =0
unsigned int n_vars() const
Definition: system.h:2086
virtual void set(const numeric_index_type i, const T value)=0
Sets v(i) = value.
Defines a dense vector for use in Finite Element-type computations.
void project_vector(NumericVector< Number > &new_vector, FunctionBase< Number > *f, FunctionBase< Gradient > *g=libmesh_nullptr, int is_adjoint=-1) const
Projects arbitrary functions onto a vector of degree of freedom values for the current system...
This action class can be used with a GenericProjector to set projection values (which must be of type...
unsigned int n_comp(const unsigned int s, const unsigned int var) const
Definition: dof_object.h:780
const std::set< boundary_id_type > & b
void operator()(const ConstElemRange &range) const
void get_side_fe(unsigned int var, FEGenericBase< OutputShape > *&fe) const
Accessor for edge/face (2D/3D) finite element object for variable var for the largest dimension in th...
Definition: fem_context.h:299
FEMFunctionBase is a base class from which users can derive in order to define "function-like" object...
dof_id_type first_old_dof(const processor_id_type proc) const
Definition: dof_map.h:545
processor_id_type processor_id()
Definition: libmesh_base.h:96
void cholesky_solve(const DenseVector< T2 > &b, DenseVector< T2 > &x)
For symmetric positive definite (SPD) matrices.
Definition: dense_matrix.C:911
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
virtual void localize(std::vector< T > &v_local) const =0
Creates a copy of the global vector in the local vector v_local.
sys get_dof_map()
void ErrorVector unsigned int
Definition: adjoints_ex3.C:360
const std::vector< Point > & get_xyz() const
Definition: fe_abstract.h:230
const Elem & get(const ElemType type_in)
This class forms the foundation from which generic finite elements may be derived.
GenericProjector(const GenericProjector &in)
uint8_t dof_id_type
Definition: id_types.h:64
virtual UniquePtr< NumericVector< T > > clone() const =0
processor_id_type n_processors()
Definition: libmesh_base.h:88
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.
Definition: dof_map.C:1917
void operator()(const ConstElemRange &range) const
static UniquePtr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.