system_projection.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 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 #include <numeric> // std::iota
23 
24 // Local includes
25 #include "libmesh/libmesh_config.h"
26 
27 #ifdef LIBMESH_HAVE_METAPHYSICL
28 
29 // FIXME - having to do this with MetaPhysicL brings me shame - RHS
31 // Template specialization declarations in here need to *precede* code
32 // using them.
33 #include "metaphysicl/dynamicsparsenumberarray_decl.h"
35 
36 #include "libmesh/compare_types.h"
37 #include "libmesh/int_range.h"
38 
39 using MetaPhysicL::DynamicSparseNumberArray;
40 
41 namespace libMesh
42 {
43 // From the perspective of libMesh gradient vectors, a DSNA is a
44 // scalar component
45 template <typename T, typename I>
46 struct ScalarTraits<MetaPhysicL::DynamicSparseNumberArray<T,I> >
47 {
48  static const bool value = true;
49 };
50 
51 // And although MetaPhysicL knows how to combine DSNA with something
52 // else, we need to teach libMesh too.
53 template <typename T, typename I, typename T2>
54 struct CompareTypes<MetaPhysicL::DynamicSparseNumberArray<T,I>, T2>
55 {
56  typedef typename
57  MetaPhysicL::DynamicSparseNumberArray
59 };
60 }
61 
62 
63 #endif
64 
65 #include "libmesh/boundary_info.h"
66 #include "libmesh/dense_matrix.h"
67 #include "libmesh/dense_vector.h"
68 #include "libmesh/dof_map.h"
69 #include "libmesh/elem.h"
70 #include "libmesh/fe_base.h"
71 #include "libmesh/fe_interface.h"
74 #include "libmesh/mesh_base.h"
75 #include "libmesh/numeric_vector.h"
76 #include "libmesh/quadrature.h"
77 #include "libmesh/sparse_matrix.h"
78 #include "libmesh/system.h"
79 #include "libmesh/threads.h"
82 
83 
84 
85 #ifdef LIBMESH_HAVE_METAPHYSICL
86 // FIXME - wrapping MetaPhysicL is shameful - RHS
88 // Include MetaPhysicL definitions finally
89 #include "metaphysicl/dynamicsparsenumberarray.h"
91 
92 // And make sure we instantiate the methods we'll need to use on them.
94 
95 namespace libMesh {
96 typedef DynamicSparseNumberArray<Real, dof_id_type> DSNAN;
97 
98 template void
101 template void
103  DenseVector<DSNAN> &) const;
104 }
105 #endif
106 
107 
108 
109 namespace libMesh
110 {
111 
112 // ------------------------------------------------------------
113 // Helper class definitions
114 
115 #ifdef LIBMESH_ENABLE_AMR
116 
128 {
129 private:
130  const System & system;
131 
132 public:
133  BuildProjectionList (const System & system_in) :
134  system(system_in),
135  send_list()
136  {}
137 
139  system(other.system),
140  send_list()
141  {}
142 
143  void unique();
144  void operator()(const ConstElemRange & range);
145  void join (const BuildProjectionList & other);
146  std::vector<dof_id_type> send_list;
147 };
148 
149 #endif // LIBMESH_ENABLE_AMR
150 
151 
158 {
159 private:
160  const std::set<boundary_id_type> & b;
161  const std::vector<unsigned int> & variables;
162  const System & system;
163  std::unique_ptr<FunctionBase<Number>> f;
164  std::unique_ptr<FunctionBase<Gradient>> g;
167 
168 public:
169  BoundaryProjectSolution (const std::set<boundary_id_type> & b_in,
170  const std::vector<unsigned int> & variables_in,
171  const System & system_in,
172  FunctionBase<Number> * f_in,
173  FunctionBase<Gradient> * g_in,
174  const Parameters & parameters_in,
175  NumericVector<Number> & new_v_in) :
176  b(b_in),
177  variables(variables_in),
178  system(system_in),
179  f(f_in ? f_in->clone() : std::unique_ptr<FunctionBase<Number>>()),
180  g(g_in ? g_in->clone() : std::unique_ptr<FunctionBase<Gradient>>()),
181  parameters(parameters_in),
182  new_vector(new_v_in)
183  {
184  libmesh_assert(f.get());
185  f->init();
186  if (g.get())
187  g->init();
188  }
189 
191  b(in.b),
192  variables(in.variables),
193  system(in.system),
194  f(in.f.get() ? in.f->clone() : std::unique_ptr<FunctionBase<Number>>()),
195  g(in.g.get() ? in.g->clone() : std::unique_ptr<FunctionBase<Gradient>>()),
196  parameters(in.parameters),
197  new_vector(in.new_vector)
198  {
199  libmesh_assert(f.get());
200  f->init();
201  if (g.get())
202  g->init();
203  }
204 
205  void operator()(const ConstElemRange & range) const;
206 };
207 
208 
209 
210 // ------------------------------------------------------------
211 // System implementation
212 void System::project_vector (NumericVector<Number> & vector,
213  int is_adjoint) const
214 {
215  // Create a copy of the vector, which currently
216  // contains the old data.
217  std::unique_ptr<NumericVector<Number>>
218  old_vector (vector.clone());
219 
220  // Project the old vector to the new vector
221  this->project_vector (*old_vector, vector, is_adjoint);
222 }
223 
224 
230 void System::project_vector (const NumericVector<Number> & old_v,
231  NumericVector<Number> & new_v,
232  int
233 #ifdef LIBMESH_ENABLE_AMR
234  is_adjoint
235 #endif
236  ) const
237 {
238  LOG_SCOPE ("project_vector(old,new)", "System");
239 
246  new_v.clear();
247 
248 #ifdef LIBMESH_ENABLE_AMR
249 
250  // Resize the new vector and get a serial version.
251  NumericVector<Number> * new_vector_ptr = nullptr;
252  std::unique_ptr<NumericVector<Number>> new_vector_built;
253  NumericVector<Number> * local_old_vector;
254  std::unique_ptr<NumericVector<Number>> local_old_vector_built;
255  const NumericVector<Number> * old_vector_ptr = nullptr;
256 
257  ConstElemRange active_local_elem_range
258  (this->get_mesh().active_local_elements_begin(),
259  this->get_mesh().active_local_elements_end());
260 
261  // If the old vector was uniprocessor, make the new
262  // vector uniprocessor
263  if (old_v.type() == SERIAL)
264  {
265  new_v.init (this->n_dofs(), false, SERIAL);
266  new_vector_ptr = &new_v;
267  old_vector_ptr = &old_v;
268  }
269 
270  // Otherwise it is a parallel, distributed vector, which
271  // we need to localize.
272  else if (old_v.type() == PARALLEL)
273  {
274  // Build a send list for efficient localization
275  BuildProjectionList projection_list(*this);
276  Threads::parallel_reduce (active_local_elem_range,
277  projection_list);
278 
279  // Create a sorted, unique send_list
280  projection_list.unique();
281 
282  new_v.init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
283  new_vector_built = NumericVector<Number>::build(this->comm());
284  local_old_vector_built = NumericVector<Number>::build(this->comm());
285  new_vector_ptr = new_vector_built.get();
286  local_old_vector = local_old_vector_built.get();
287  new_vector_ptr->init(this->n_dofs(), false, SERIAL);
288  local_old_vector->init(old_v.size(), false, SERIAL);
289  old_v.localize(*local_old_vector, projection_list.send_list);
290  local_old_vector->close();
291  old_vector_ptr = local_old_vector;
292  }
293  else if (old_v.type() == GHOSTED)
294  {
295  // Build a send list for efficient localization
296  BuildProjectionList projection_list(*this);
297  Threads::parallel_reduce (active_local_elem_range,
298  projection_list);
299 
300  // Create a sorted, unique send_list
301  projection_list.unique();
302 
303  new_v.init (this->n_dofs(), this->n_local_dofs(),
304  this->get_dof_map().get_send_list(), false, GHOSTED);
305 
306  local_old_vector_built = NumericVector<Number>::build(this->comm());
307  new_vector_ptr = &new_v;
308  local_old_vector = local_old_vector_built.get();
309  local_old_vector->init(old_v.size(), old_v.local_size(),
310  projection_list.send_list, false, GHOSTED);
311  old_v.localize(*local_old_vector, projection_list.send_list);
312  local_old_vector->close();
313  old_vector_ptr = local_old_vector;
314  }
315  else // unknown old_v.type()
316  libmesh_error_msg("ERROR: Unknown old_v.type() == " << old_v.type());
317 
318  // Note that the above will have zeroed the new_vector.
319  // Just to be sure, assert that new_vector_ptr and old_vector_ptr
320  // were successfully set before trying to deref them.
321  libmesh_assert(new_vector_ptr);
322  libmesh_assert(old_vector_ptr);
323 
324  NumericVector<Number> & new_vector = *new_vector_ptr;
325  const NumericVector<Number> & old_vector = *old_vector_ptr;
326 
327  const unsigned int n_variables = this->n_vars();
328 
329  if (n_variables)
330  {
331  std::vector<unsigned int> vars(n_variables);
332  std::iota(vars.begin(), vars.end(), 0);
333 
334  // Use a typedef to make the calling sequence for parallel_for() a bit more readable
335  typedef
336  GenericProjector<OldSolutionValue<Number, &FEMContext::point_value>,
337  OldSolutionValue<Gradient, &FEMContext::point_gradient>,
338  Number, VectorSetAction<Number>> FEMProjector;
339 
340  OldSolutionValue<Number, &FEMContext::point_value> f(*this, old_vector);
341  OldSolutionValue<Gradient, &FEMContext::point_gradient> g(*this, old_vector);
342  VectorSetAction<Number> setter(new_vector);
343 
344  Threads::parallel_for (active_local_elem_range,
345  FEMProjector(*this, f, &g, setter, vars));
346 
347  // Copy the SCALAR dofs from old_vector to new_vector
348  // Note: We assume that all SCALAR dofs are on the
349  // processor with highest ID
350  if (this->processor_id() == (this->n_processors()-1))
351  {
352  const DofMap & dof_map = this->get_dof_map();
353  for (unsigned int var=0; var<this->n_vars(); var++)
354  if (this->variable(var).type().family == SCALAR)
355  {
356  // We can just map SCALAR dofs directly across
357  std::vector<dof_id_type> new_SCALAR_indices, old_SCALAR_indices;
358  dof_map.SCALAR_dof_indices (new_SCALAR_indices, var, false);
359  dof_map.SCALAR_dof_indices (old_SCALAR_indices, var, true);
360  for (auto i : index_range(new_SCALAR_indices))
361  new_vector.set(new_SCALAR_indices[i], old_vector(old_SCALAR_indices[i]));
362  }
363  }
364  }
365 
366  new_vector.close();
367 
368  // If the old vector was serial, we probably need to send our values
369  // to other processors
370  //
371  // FIXME: I'm not sure how to make a NumericVector do that without
372  // creating a temporary parallel vector to use localize! - RHS
373  if (old_v.type() == SERIAL)
374  {
375  std::unique_ptr<NumericVector<Number>> dist_v = NumericVector<Number>::build(this->comm());
376  dist_v->init(this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
377  dist_v->close();
378 
379  for (dof_id_type i=0; i!=dist_v->size(); i++)
380  if (new_vector(i) != 0.0)
381  dist_v->set(i, new_vector(i));
382 
383  dist_v->close();
384 
385  dist_v->localize (new_v, this->get_dof_map().get_send_list());
386  new_v.close();
387  }
388  // If the old vector was parallel, we need to update it
389  // and free the localized copies
390  else if (old_v.type() == PARALLEL)
391  {
392  // We may have to set dof values that this processor doesn't
393  // own in certain special cases, like LAGRANGE FIRST or
394  // HERMITE THIRD elements on second-order meshes
395  for (dof_id_type i=0; i!=new_v.size(); i++)
396  if (new_vector(i) != 0.0)
397  new_v.set(i, new_vector(i));
398  new_v.close();
399  }
400 
401  if (is_adjoint == -1)
402  this->get_dof_map().enforce_constraints_exactly(*this, &new_v);
403  else if (is_adjoint >= 0)
404  this->get_dof_map().enforce_adjoint_constraints_exactly(new_v,
405  is_adjoint);
406 
407 #else
408 
409  // AMR is disabled: simply copy the vector
410  new_v = old_v;
411 
412 #endif // #ifdef LIBMESH_ENABLE_AMR
413 }
414 
415 
416 #ifdef LIBMESH_ENABLE_AMR
417 #ifdef LIBMESH_HAVE_METAPHYSICL
418 
419 template <typename Output>
421 {
422 public:
423  typedef DynamicSparseNumberArray<Output, dof_id_type> type;
424 };
425 
426 template <typename InnerOutput>
427 class DSNAOutput<VectorValue<InnerOutput> >
428 {
429 public:
431 };
432 
442 template <typename Output,
443  void (FEMContext::*point_output) (unsigned int,
444  const Point &,
445  Output &,
446  const Real) const>
448 {
449 public:
450  typedef typename DSNAOutput<Output>::type DSNA;
451 
453  OldSolutionBase<Output, point_output>(sys_in)
454  {
455  this->old_context.set_algebraic_type(FEMContext::OLD_DOFS_ONLY);
456  }
457 
459  OldSolutionBase<Output, point_output>(in.sys)
460  {
461  this->old_context.set_algebraic_type(FEMContext::OLD_DOFS_ONLY);
462  }
463 
464  DSNA eval_at_node (const FEMContext & c,
465  unsigned int i,
466  unsigned int elem_dim,
467  const Node & n,
468  Real /* time */ = 0.);
469 
470  DSNA eval_at_point(const FEMContext & c,
471  unsigned int i,
472  const Point & p,
473  Real /* time */ = 0.);
474 
475  void eval_old_dofs (const FEMContext & c,
476  unsigned int var,
477  std::vector<DSNA> & values)
478  {
479  LOG_SCOPE ("eval_old_dofs()", "OldSolutionValue");
480 
481  this->check_old_context(c);
482 
483  const std::vector<dof_id_type> & old_dof_indices =
484  this->old_context.get_dof_indices(var);
485 
486  libmesh_assert_equal_to (old_dof_indices.size(), values.size());
487 
488  for (auto i : index_range(values))
489  {
490  values[i].resize(1);
491  values[i].raw_at(0) = 1;
492  values[i].raw_index(0) = old_dof_indices[i];
493  }
494  }
495 };
496 
497 
498 
499 template<>
500 inline
501 DynamicSparseNumberArray<Real, dof_id_type>
502 OldSolutionCoefs<Real, &FEMContext::point_value>::
503 eval_at_point(const FEMContext & c,
504  unsigned int i,
505  const Point & p,
506  Real /* time */)
507 {
508  LOG_SCOPE ("eval_at_point()", "OldSolutionCoefs");
509 
510  if (!this->check_old_context(c, p))
511  return 0;
512 
513  // Get finite element object
514  FEGenericBase<Real> * fe = nullptr;
515  this->old_context.get_element_fe<Real>
516  (i, fe, this->old_context.get_elem_dim());
517 
518  // Build a FE for calculating phi(p)
519  FEGenericBase<Real> * fe_new =
520  this->old_context.build_new_fe(fe, p);
521 
522  // Get the values and global indices of the shape functions
523  const std::vector<std::vector<Real> > & phi = fe_new->get_phi();
524  const std::vector<dof_id_type> & dof_indices =
525  this->old_context.get_dof_indices(i);
526 
527  const std::size_t n_dofs = phi.size();
528  libmesh_assert_equal_to(n_dofs, dof_indices.size());
529 
530  DynamicSparseNumberArray<Real, dof_id_type> returnval;
531  returnval.resize(n_dofs);
532 
533  for (auto j : index_range(phi))
534  {
535  returnval.raw_at(j) = phi[j][0];
536  returnval.raw_index(j) = dof_indices[j];
537  }
538 
539  return returnval;
540 }
541 
542 
543 
544 template<>
545 inline
549  unsigned int i,
550  const Point & p,
551  Real /* time */)
552 {
553  LOG_SCOPE ("eval_at_point()", "OldSolutionCoefs");
554 
555  if (!this->check_old_context(c, p))
556  return 0;
557 
558  // Get finite element object
559  FEGenericBase<Real> * fe = nullptr;
560  this->old_context.get_element_fe<Real>
561  (i, fe, this->old_context.get_elem_dim());
562 
563  // Build a FE for calculating phi(p)
564  FEGenericBase<Real> * fe_new =
565  this->old_context.build_new_fe(fe, p);
566 
567  // Get the values and global indices of the shape functions
568  const std::vector<std::vector<RealGradient> > & dphi = fe_new->get_dphi();
569  const std::vector<dof_id_type> & dof_indices =
570  this->old_context.get_dof_indices(i);
571 
572  const std::size_t n_dofs = dphi.size();
573  libmesh_assert_equal_to(n_dofs, dof_indices.size());
574 
576 
577  for (unsigned int d = 0; d != LIBMESH_DIM; ++d)
578  returnval(d).resize(n_dofs);
579 
580  for (auto j : index_range(dphi))
581  for (int d = 0; d != LIBMESH_DIM; ++d)
582  {
583  returnval(d).raw_at(j) = dphi[j][0](d);
584  returnval(d).raw_index(j) = dof_indices[j];
585  }
586 
587  return returnval;
588 }
589 
590 
591 template<>
592 inline
593 DynamicSparseNumberArray<Real, dof_id_type>
596  unsigned int i,
597  unsigned int /* elem_dim */,
598  const Node & n,
599  Real /* time */)
600 {
601  LOG_SCOPE ("Real eval_at_node()", "OldSolutionCoefs");
602 
603  // Optimize for the common case, where this node was part of the
604  // old solution.
605  //
606  // Be sure to handle cases where the variable wasn't defined on
607  // this node (due to changing subdomain support) or where the
608  // variable has no components on this node (due to Elem order
609  // exceeding FE order)
610  if (n.old_dof_object &&
611  n.old_dof_object->n_vars(sys.number()) &&
612  n.old_dof_object->n_comp(sys.number(), i))
613  {
614  DynamicSparseNumberArray<Real, dof_id_type> returnval;
615  const dof_id_type old_id =
616  n.old_dof_object->dof_number(sys.number(), i, 0);
617  returnval.resize(1);
618  returnval.raw_at(0) = 1;
619  returnval.raw_index(0) = old_id;
620  return returnval;
621  }
622 
623  return this->eval_at_point(c, i, n, 0);
624 }
625 
626 
627 
628 template<>
629 inline
633  unsigned int i,
634  unsigned int elem_dim,
635  const Node & n,
636  Real /* time */)
637 {
638  LOG_SCOPE ("RealGradient eval_at_node()", "OldSolutionCoefs");
639 
640  // Optimize for the common case, where this node was part of the
641  // old solution.
642  //
643  // Be sure to handle cases where the variable wasn't defined on
644  // this node (due to changing subdomain support) or where the
645  // variable has no components on this node (due to Elem order
646  // exceeding FE order)
647  if (n.old_dof_object &&
648  n.old_dof_object->n_vars(sys.number()) &&
649  n.old_dof_object->n_comp(sys.number(), i))
650  {
652  for (unsigned int d = 0; d != elem_dim; ++d)
653  {
654  const dof_id_type old_id =
655  n.old_dof_object->dof_number(sys.number(), i, d+1);
656  g(d).resize(1);
657  g(d).raw_at(0) = 1;
658  g(d).raw_index(0) = old_id;
659  }
660  return g;
661  }
662 
663  return this->eval_at_point(c, i, n, 0);
664 }
665 
666 
667 
676 template <typename ValIn, typename ValOut>
678 {
679 private:
681 
682 public:
684  target_matrix(target_mat) {}
685 
686  void insert(const FEMContext & c,
687  unsigned int var_num,
688  const DenseVector<DynamicSparseNumberArray<ValIn, dof_id_type> > & Ue)
689  {
690  const numeric_index_type
691  begin_dof = target_matrix.row_start(),
692  end_dof = target_matrix.row_stop();
693 
694  const std::vector<dof_id_type> & dof_indices =
695  c.get_dof_indices(var_num);
696 
697  unsigned int size = Ue.size();
698 
699  libmesh_assert_equal_to(size, dof_indices.size());
700 
701  // Lock the target matrix since it is shared among threads.
702  {
703  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
704 
705  for (unsigned int i = 0; i != size; ++i)
706  {
707  const dof_id_type dof_i = dof_indices[i];
708  if ((dof_i >= begin_dof) && (dof_i < end_dof))
709  {
710  const DynamicSparseNumberArray<ValIn,dof_id_type> & dnsa = Ue(i);
711  const std::size_t dnsa_size = dnsa.size();
712  for (unsigned int j = 0; j != dnsa_size; ++j)
713  {
714  const dof_id_type dof_j = dnsa.raw_index(j);
715  const ValIn dof_val = dnsa.raw_at(j);
716  target_matrix.set(dof_i, dof_j, dof_val);
717  }
718  }
719  }
720  }
721  }
722 };
723 
724 
725 
730 void System::projection_matrix (SparseMatrix<Number> & proj_mat) const
731 {
732  LOG_SCOPE ("projection_matrix()", "System");
733 
734  const unsigned int n_variables = this->n_vars();
735 
736  if (n_variables)
737  {
738  ConstElemRange active_local_elem_range
739  (this->get_mesh().active_local_elements_begin(),
740  this->get_mesh().active_local_elements_end());
741 
742  std::vector<unsigned int> vars(n_variables);
743  std::iota(vars.begin(), vars.end(), 0);
744 
745  // Use a typedef to make the calling sequence for parallel_for() a bit more readable
746  typedef OldSolutionCoefs<Real, &FEMContext::point_value> OldSolutionValueCoefs;
747  typedef OldSolutionCoefs<RealGradient, &FEMContext::point_gradient> OldSolutionGradientCoefs;
748 
749  typedef
750  GenericProjector<OldSolutionValueCoefs,
751  OldSolutionGradientCoefs,
752  DynamicSparseNumberArray<Real,dof_id_type>,
753  MatrixFillAction<Real, Number> > ProjMatFiller;
754 
755  OldSolutionValueCoefs f(*this);
756  OldSolutionGradientCoefs g(*this);
757  MatrixFillAction<Real, Number> setter(proj_mat);
758 
759  Threads::parallel_for (active_local_elem_range,
760  ProjMatFiller(*this, f, &g, setter, vars));
761 
762  // Set the SCALAR dof transfer entries too.
763  // Note: We assume that all SCALAR dofs are on the
764  // processor with highest ID
765  if (this->processor_id() == (this->n_processors()-1))
766  {
767  const DofMap & dof_map = this->get_dof_map();
768  for (unsigned int var=0; var<this->n_vars(); var++)
769  if (this->variable(var).type().family == SCALAR)
770  {
771  // We can just map SCALAR dofs directly across
772  std::vector<dof_id_type> new_SCALAR_indices, old_SCALAR_indices;
773  dof_map.SCALAR_dof_indices (new_SCALAR_indices, var, false);
774  dof_map.SCALAR_dof_indices (old_SCALAR_indices, var, true);
775  const unsigned int new_n_dofs =
776  cast_int<unsigned int>(new_SCALAR_indices.size());
777 
778  for (unsigned int i=0; i<new_n_dofs; i++)
779  {
780  proj_mat.set( new_SCALAR_indices[i],
781  old_SCALAR_indices[i], 1);
782  }
783  }
784  }
785  }
786 }
787 #endif // LIBMESH_HAVE_METAPHYSICL
788 #endif // LIBMESH_ENABLE_AMR
789 
790 
791 
796 void System::project_solution (ValueFunctionPointer fptr,
797  GradientFunctionPointer gptr,
798  const Parameters & parameters) const
799 {
800  WrappedFunction<Number> f(*this, fptr, &parameters);
801  WrappedFunction<Gradient> g(*this, gptr, &parameters);
802  this->project_solution(&f, &g);
803 }
804 
805 
810 void System::project_solution (FunctionBase<Number> * f,
811  FunctionBase<Gradient> * g) const
812 {
813  this->project_vector(*solution, f, g);
814 
815  solution->localize(*current_local_solution, _dof_map->get_send_list());
816 }
817 
818 
823 void System::project_solution (FEMFunctionBase<Number> * f,
824  FEMFunctionBase<Gradient> * g) const
825 {
826  this->project_vector(*solution, f, g);
827 
828  solution->localize(*current_local_solution, _dof_map->get_send_list());
829 }
830 
831 
836 void System::project_vector (ValueFunctionPointer fptr,
837  GradientFunctionPointer gptr,
838  const Parameters & parameters,
839  NumericVector<Number> & new_vector,
840  int is_adjoint) const
841 {
842  WrappedFunction<Number> f(*this, fptr, &parameters);
843  WrappedFunction<Gradient> g(*this, gptr, &parameters);
844  this->project_vector(new_vector, &f, &g, is_adjoint);
845 }
846 
851 void System::project_vector (NumericVector<Number> & new_vector,
854  int is_adjoint) const
855 {
856  LOG_SCOPE ("project_vector(FunctionBase)", "System");
857 
858  libmesh_assert(f);
859 
860  WrappedFunctor<Number> f_fem(*f);
861 
862  if (g)
863  {
864  WrappedFunctor<Gradient> g_fem(*g);
865 
866  this->project_vector(new_vector, &f_fem, &g_fem, is_adjoint);
867  }
868  else
869  this->project_vector(new_vector, &f_fem, nullptr, is_adjoint);
870 }
871 
872 
877 void System::project_vector (NumericVector<Number> & new_vector,
880  int is_adjoint) const
881 {
882  LOG_SCOPE ("project_fem_vector()", "System");
883 
884  libmesh_assert (f);
885 
886  ConstElemRange active_local_range
887  (this->get_mesh().active_local_elements_begin(),
888  this->get_mesh().active_local_elements_end() );
889 
890  VectorSetAction<Number> setter(new_vector);
891 
892  const unsigned int n_variables = this->n_vars();
893 
894  std::vector<unsigned int> vars(n_variables);
895  std::iota(vars.begin(), vars.end(), 0);
896 
897  // Use a typedef to make the calling sequence for parallel_for() a bit more readable
898  typedef
900  Number, VectorSetAction<Number>> FEMProjector;
901 
903 
904  if (g)
905  {
907 
909  (active_local_range,
910  FEMProjector(*this, fw, &gw, setter, vars));
911  }
912  else
914  (active_local_range,
915  FEMProjector(*this, fw, nullptr, setter, vars));
916 
917  // Also, load values into the SCALAR dofs
918  // Note: We assume that all SCALAR dofs are on the
919  // processor with highest ID
920  if (this->processor_id() == (this->n_processors()-1))
921  {
922  // FIXME: Do we want to first check for SCALAR vars before building this? [PB]
923  FEMContext context( *this );
924 
925  const DofMap & dof_map = this->get_dof_map();
926  for (unsigned int var=0; var<this->n_vars(); var++)
927  if (this->variable(var).type().family == SCALAR)
928  {
929  // FIXME: We reinit with an arbitrary element in case the user
930  // doesn't override FEMFunctionBase::component. Is there
931  // any use case we're missing? [PB]
932  Elem * el = const_cast<Elem *>(*(this->get_mesh().active_local_elements_begin()));
933  context.pre_fe_reinit(*this, el);
934 
935  std::vector<dof_id_type> SCALAR_indices;
936  dof_map.SCALAR_dof_indices (SCALAR_indices, var);
937  const unsigned int n_SCALAR_dofs =
938  cast_int<unsigned int>(SCALAR_indices.size());
939 
940  for (unsigned int i=0; i<n_SCALAR_dofs; i++)
941  {
942  const dof_id_type global_index = SCALAR_indices[i];
943  const unsigned int component_index =
944  this->variable_scalar_number(var,i);
945 
946  new_vector.set(global_index, f->component(context, component_index, Point(), this->time));
947  }
948  }
949  }
950 
951  new_vector.close();
952 
953 #ifdef LIBMESH_ENABLE_CONSTRAINTS
954  if (is_adjoint == -1)
955  this->get_dof_map().enforce_constraints_exactly(*this, &new_vector);
956  else if (is_adjoint >= 0)
957  this->get_dof_map().enforce_adjoint_constraints_exactly(new_vector,
958  is_adjoint);
959 #endif
960 }
961 
962 
968 void System::boundary_project_solution (const std::set<boundary_id_type> & b,
969  const std::vector<unsigned int> & variables,
970  ValueFunctionPointer fptr,
971  GradientFunctionPointer gptr,
972  const Parameters & parameters)
973 {
974  WrappedFunction<Number> f(*this, fptr, &parameters);
975  WrappedFunction<Gradient> g(*this, gptr, &parameters);
976  this->boundary_project_solution(b, variables, &f, &g);
977 }
978 
979 
985 void System::boundary_project_solution (const std::set<boundary_id_type> & b,
986  const std::vector<unsigned int> & variables,
989 {
990  this->boundary_project_vector(b, variables, *solution, f, g);
991 
992  solution->localize(*current_local_solution);
993 }
994 
995 
996 
997 
998 
1003 void System::boundary_project_vector (const std::set<boundary_id_type> & b,
1004  const std::vector<unsigned int> & variables,
1005  ValueFunctionPointer fptr,
1006  GradientFunctionPointer gptr,
1007  const Parameters & parameters,
1008  NumericVector<Number> & new_vector,
1009  int is_adjoint) const
1010 {
1011  WrappedFunction<Number> f(*this, fptr, &parameters);
1012  WrappedFunction<Gradient> g(*this, gptr, &parameters);
1013  this->boundary_project_vector(b, variables, new_vector, &f, &g,
1014  is_adjoint);
1015 }
1016 
1021 void System::boundary_project_vector (const std::set<boundary_id_type> & b,
1022  const std::vector<unsigned int> & variables,
1023  NumericVector<Number> & new_vector,
1026  int is_adjoint) const
1027 {
1028  LOG_SCOPE ("boundary_project_vector()", "System");
1029 
1031  (ConstElemRange (this->get_mesh().active_local_elements_begin(),
1032  this->get_mesh().active_local_elements_end() ),
1033  BoundaryProjectSolution(b, variables, *this, f, g,
1034  this->get_equation_systems().parameters,
1035  new_vector)
1036  );
1037 
1038  // We don't do SCALAR dofs when just projecting the boundary, so
1039  // we're done here.
1040 
1041  new_vector.close();
1042 
1043 #ifdef LIBMESH_ENABLE_CONSTRAINTS
1044  if (is_adjoint == -1)
1045  this->get_dof_map().enforce_constraints_exactly(*this, &new_vector);
1046  else if (is_adjoint >= 0)
1047  this->get_dof_map().enforce_adjoint_constraints_exactly(new_vector,
1048  is_adjoint);
1049 #endif
1050 }
1051 
1052 
1053 
1054 #ifdef LIBMESH_ENABLE_AMR
1055 void BuildProjectionList::unique()
1056 {
1057  // Sort the send list. After this duplicated
1058  // elements will be adjacent in the vector
1059  std::sort(this->send_list.begin(),
1060  this->send_list.end());
1061 
1062  // Now use std::unique to remove duplicate entries
1063  std::vector<dof_id_type>::iterator new_end =
1064  std::unique (this->send_list.begin(),
1065  this->send_list.end());
1066 
1067  // Remove the end of the send_list. Use the "swap trick"
1068  // from Effective STL
1069  std::vector<dof_id_type>
1070  (this->send_list.begin(), new_end).swap (this->send_list);
1071 }
1072 
1073 
1074 
1075 void BuildProjectionList::operator()(const ConstElemRange & range)
1076 {
1077  // The DofMap for this system
1078  const DofMap & dof_map = system.get_dof_map();
1079 
1080  const dof_id_type first_old_dof = dof_map.first_old_dof();
1081  const dof_id_type end_old_dof = dof_map.end_old_dof();
1082 
1083  // We can handle all the variables at once.
1084  // The old global DOF indices
1085  std::vector<dof_id_type> di;
1086 
1087  // Iterate over the elements in the range
1088  for (const auto & elem : range)
1089  {
1090  // If this element doesn't have an old_dof_object with dofs for the
1091  // current system, then it must be newly added, so the user
1092  // is responsible for setting the new dofs.
1093 
1094  // ... but we need a better way to test for that; the code
1095  // below breaks on any FE type for which the elem stores no
1096  // dofs.
1097  // if (!elem->old_dof_object || !elem->old_dof_object->has_dofs(system.number()))
1098  // continue;
1099 
1100  // Examining refinement flags instead should distinguish
1101  // between refinement-added and user-added elements lacking
1102  // old_dof_object
1103  if (!elem->old_dof_object &&
1104  elem->refinement_flag() != Elem::JUST_REFINED &&
1105  elem->refinement_flag() != Elem::JUST_COARSENED)
1106  continue;
1107 
1108  const Elem * parent = elem->parent();
1109 
1110  if (elem->refinement_flag() == Elem::JUST_REFINED)
1111  {
1112  libmesh_assert(parent);
1113 
1114  // We used to hack_p_level here, but that wasn't thread-safe
1115  // so now we take p refinement flags into account in
1116  // old_dof_indices
1117 
1118  dof_map.old_dof_indices (parent, di);
1119 
1120  for (auto & node : elem->node_ref_range())
1121  {
1122  const DofObject * old_dofs = node.old_dof_object;
1123 
1124  if (old_dofs)
1125  {
1126  const unsigned int sysnum = system.number();
1127  const unsigned int nvg = old_dofs->n_var_groups(sysnum);
1128 
1129  for (unsigned int vg=0; vg != nvg; ++vg)
1130  {
1131  const unsigned int nvig =
1132  old_dofs->n_vars(sysnum, vg);
1133  for (unsigned int vig=0; vig != nvig; ++vig)
1134  {
1135  const unsigned int n_comp =
1136  old_dofs->n_comp_group(sysnum, vg);
1137  for (unsigned int c=0; c != n_comp; ++c)
1138  di.push_back
1139  (old_dofs->dof_number(sysnum, vg, vig, c, n_comp));
1140  }
1141  }
1142  }
1143  }
1144 
1145  std::sort(di.begin(), di.end());
1146  std::vector<dof_id_type>::iterator new_end =
1147  std::unique(di.begin(), di.end());
1148  std::vector<dof_id_type>(di.begin(), new_end).swap(di);
1149  }
1150  else if (elem->refinement_flag() == Elem::JUST_COARSENED)
1151  {
1152  std::vector<dof_id_type> di_child;
1153  di.clear();
1154  for (auto & child : elem->child_ref_range())
1155  {
1156  dof_map.old_dof_indices (&child, di_child);
1157  di.insert(di.end(), di_child.begin(), di_child.end());
1158  }
1159  }
1160  else
1161  dof_map.old_dof_indices (elem, di);
1162 
1163  for (std::size_t i=0; i != di.size(); ++i)
1164  if (di[i] < first_old_dof || di[i] >= end_old_dof)
1165  this->send_list.push_back(di[i]);
1166  } // end elem loop
1167 }
1168 
1169 
1170 
1171 void BuildProjectionList::join(const BuildProjectionList & other)
1172 {
1173  // Joining simply requires I add the dof indices from the other object
1174  this->send_list.insert(this->send_list.end(),
1175  other.send_list.begin(),
1176  other.send_list.end());
1177 }
1178 #endif // LIBMESH_ENABLE_AMR
1179 
1180 
1181 
1182 void BoundaryProjectSolution::operator()(const ConstElemRange & range) const
1183 {
1184  // We need data to project
1185  libmesh_assert(f.get());
1186 
1194  // The dimensionality of the current mesh
1195  const unsigned int dim = system.get_mesh().mesh_dimension();
1196 
1197  // The DofMap for this system
1198  const DofMap & dof_map = system.get_dof_map();
1199 
1200  // Boundary info for the current mesh
1201  const BoundaryInfo & boundary_info =
1202  system.get_mesh().get_boundary_info();
1203 
1204  // The element matrix and RHS for projections.
1205  // Note that Ke is always real-valued, whereas
1206  // Fe may be complex valued if complex number
1207  // support is enabled
1208  DenseMatrix<Real> Ke;
1210  // The new element coefficients
1212 
1213 
1214  // Loop over all the variables we've been requested to project
1215  for (std::size_t v=0; v!=variables.size(); v++)
1216  {
1217  const unsigned int var = variables[v];
1218 
1219  const Variable & variable = dof_map.variable(var);
1220 
1221  const FEType & fe_type = variable.type();
1222 
1223  if (fe_type.family == SCALAR)
1224  continue;
1225 
1226  const unsigned int var_component =
1227  system.variable_scalar_number(var, 0);
1228 
1229  // Get FE objects of the appropriate type
1230  std::unique_ptr<FEBase> fe (FEBase::build(dim, fe_type));
1231 
1232  // Prepare variables for projection
1233  std::unique_ptr<QBase> qedgerule (fe_type.default_quadrature_rule(1));
1234  std::unique_ptr<QBase> qsiderule (fe_type.default_quadrature_rule(dim-1));
1235 
1236  // The values of the shape functions at the quadrature
1237  // points
1238  const std::vector<std::vector<Real>> & phi = fe->get_phi();
1239 
1240  // The gradients of the shape functions at the quadrature
1241  // points on the child element.
1242  const std::vector<std::vector<RealGradient>> * dphi = nullptr;
1243 
1244  const FEContinuity cont = fe->get_continuity();
1245 
1246  if (cont == C_ONE)
1247  {
1248  // We'll need gradient data for a C1 projection
1249  libmesh_assert(g.get());
1250 
1251  const std::vector<std::vector<RealGradient>> &
1252  ref_dphi = fe->get_dphi();
1253  dphi = &ref_dphi;
1254  }
1255 
1256  // The Jacobian * quadrature weight at the quadrature points
1257  const std::vector<Real> & JxW =
1258  fe->get_JxW();
1259 
1260  // The XYZ locations of the quadrature points
1261  const std::vector<Point> & xyz_values =
1262  fe->get_xyz();
1263 
1264  // The global DOF indices
1265  std::vector<dof_id_type> dof_indices;
1266  // Side/edge DOF indices
1267  std::vector<unsigned int> side_dofs;
1268 
1269  // Container to catch IDs passed back from BoundaryInfo.
1270  std::vector<boundary_id_type> bc_ids;
1271 
1272  // Iterate over all the elements in the range
1273  for (const auto & elem : range)
1274  {
1275  // Per-subdomain variables don't need to be projected on
1276  // elements where they're not active
1277  if (!variable.active_on_subdomain(elem->subdomain_id()))
1278  continue;
1279 
1280  const unsigned short n_nodes = elem->n_nodes();
1281  const unsigned short n_edges = elem->n_edges();
1282  const unsigned short n_sides = elem->n_sides();
1283 
1284  // Find out which nodes, edges and sides are on a requested
1285  // boundary:
1286  std::vector<bool> is_boundary_node(n_nodes, false),
1287  is_boundary_edge(n_edges, false),
1288  is_boundary_side(n_sides, false);
1289 
1290  // We also maintain a separate list of nodeset-based boundary nodes
1291  std::vector<bool> is_boundary_nodeset(n_nodes, false);
1292 
1293  for (unsigned char s=0; s != n_sides; ++s)
1294  {
1295  // First see if this side has been requested
1296  boundary_info.boundary_ids (elem, s, bc_ids);
1297  bool do_this_side = false;
1298  for (const auto & bc_id : bc_ids)
1299  if (b.count(bc_id))
1300  {
1301  do_this_side = true;
1302  break;
1303  }
1304  if (!do_this_side)
1305  continue;
1306 
1307  is_boundary_side[s] = true;
1308 
1309  // Then see what nodes and what edges are on it
1310  for (unsigned int n=0; n != n_nodes; ++n)
1311  if (elem->is_node_on_side(n,s))
1312  is_boundary_node[n] = true;
1313  for (unsigned int e=0; e != n_edges; ++e)
1314  if (elem->is_edge_on_side(e,s))
1315  is_boundary_edge[e] = true;
1316  }
1317 
1318  // We can also project on nodes, so we should also independently
1319  // check whether the nodes have been requested
1320  for (unsigned int n=0; n != n_nodes; ++n)
1321  {
1322  boundary_info.boundary_ids (elem->node_ptr(n), bc_ids);
1323 
1324  for (const auto & bc_id : bc_ids)
1325  if (b.count(bc_id))
1326  {
1327  is_boundary_node[n] = true;
1328  is_boundary_nodeset[n] = true;
1329  }
1330  }
1331 
1332  // We can also project on edges, so we should also independently
1333  // check whether the edges have been requested
1334  for (unsigned short e=0; e != n_edges; ++e)
1335  {
1336  boundary_info.edge_boundary_ids (elem, e, bc_ids);
1337 
1338  for (const auto & bc_id : bc_ids)
1339  if (b.count(bc_id))
1340  is_boundary_edge[e] = true;
1341  }
1342 
1343  // Update the DOF indices for this element based on
1344  // the current mesh
1345  dof_map.dof_indices (elem, dof_indices, var);
1346 
1347  // The number of DOFs on the element
1348  const unsigned int n_dofs =
1349  cast_int<unsigned int>(dof_indices.size());
1350 
1351  // Fixed vs. free DoFs on edge/face projections
1352  std::vector<char> dof_is_fixed(n_dofs, false); // bools
1353  std::vector<int> free_dof(n_dofs, 0);
1354 
1355  // The element type
1356  const ElemType elem_type = elem->type();
1357 
1358  // Zero the interpolated values
1359  Ue.resize (n_dofs); Ue.zero();
1360 
1361  // In general, we need a series of
1362  // projections to ensure a unique and continuous
1363  // solution. We start by interpolating boundary nodes, then
1364  // hold those fixed and project boundary edges, then hold
1365  // those fixed and project boundary faces,
1366 
1367  // Interpolate node values first
1368  unsigned int current_dof = 0;
1369  for (unsigned short n = 0; n != n_nodes; ++n)
1370  {
1371  // FIXME: this should go through the DofMap,
1372  // not duplicate dof_indices code badly!
1373  const unsigned int nc =
1374  FEInterface::n_dofs_at_node (dim, fe_type, elem_type,
1375  n);
1376  if ((!elem->is_vertex(n) || !is_boundary_node[n]) &&
1377  !is_boundary_nodeset[n])
1378  {
1379  current_dof += nc;
1380  continue;
1381  }
1382  if (cont == DISCONTINUOUS)
1383  {
1384  libmesh_assert_equal_to (nc, 0);
1385  }
1386  // Assume that C_ZERO elements have a single nodal
1387  // value shape function
1388  else if (cont == C_ZERO)
1389  {
1390  libmesh_assert_equal_to (nc, 1);
1391  Ue(current_dof) = f->component(var_component,
1392  elem->point(n),
1393  system.time);
1394  dof_is_fixed[current_dof] = true;
1395  current_dof++;
1396  }
1397  // The hermite element vertex shape functions are weird
1398  else if (fe_type.family == HERMITE)
1399  {
1400  Ue(current_dof) = f->component(var_component,
1401  elem->point(n),
1402  system.time);
1403  dof_is_fixed[current_dof] = true;
1404  current_dof++;
1405  Gradient grad = g->component(var_component,
1406  elem->point(n),
1407  system.time);
1408  // x derivative
1409  Ue(current_dof) = grad(0);
1410  dof_is_fixed[current_dof] = true;
1411  current_dof++;
1412  if (dim > 1)
1413  {
1414  // We'll finite difference mixed derivatives
1415  Point nxminus = elem->point(n),
1416  nxplus = elem->point(n);
1417  nxminus(0) -= TOLERANCE;
1418  nxplus(0) += TOLERANCE;
1419  Gradient gxminus = g->component(var_component,
1420  nxminus,
1421  system.time);
1422  Gradient gxplus = g->component(var_component,
1423  nxplus,
1424  system.time);
1425  // y derivative
1426  Ue(current_dof) = grad(1);
1427  dof_is_fixed[current_dof] = true;
1428  current_dof++;
1429  // xy derivative
1430  Ue(current_dof) = (gxplus(1) - gxminus(1))
1431  / 2. / TOLERANCE;
1432  dof_is_fixed[current_dof] = true;
1433  current_dof++;
1434 
1435  if (dim > 2)
1436  {
1437  // z derivative
1438  Ue(current_dof) = grad(2);
1439  dof_is_fixed[current_dof] = true;
1440  current_dof++;
1441  // xz derivative
1442  Ue(current_dof) = (gxplus(2) - gxminus(2))
1443  / 2. / TOLERANCE;
1444  dof_is_fixed[current_dof] = true;
1445  current_dof++;
1446  // We need new points for yz
1447  Point nyminus = elem->point(n),
1448  nyplus = elem->point(n);
1449  nyminus(1) -= TOLERANCE;
1450  nyplus(1) += TOLERANCE;
1451  Gradient gyminus = g->component(var_component,
1452  nyminus,
1453  system.time);
1454  Gradient gyplus = g->component(var_component,
1455  nyplus,
1456  system.time);
1457  // xz derivative
1458  Ue(current_dof) = (gyplus(2) - gyminus(2))
1459  / 2. / TOLERANCE;
1460  dof_is_fixed[current_dof] = true;
1461  current_dof++;
1462  // Getting a 2nd order xyz is more tedious
1463  Point nxmym = elem->point(n),
1464  nxmyp = elem->point(n),
1465  nxpym = elem->point(n),
1466  nxpyp = elem->point(n);
1467  nxmym(0) -= TOLERANCE;
1468  nxmym(1) -= TOLERANCE;
1469  nxmyp(0) -= TOLERANCE;
1470  nxmyp(1) += TOLERANCE;
1471  nxpym(0) += TOLERANCE;
1472  nxpym(1) -= TOLERANCE;
1473  nxpyp(0) += TOLERANCE;
1474  nxpyp(1) += TOLERANCE;
1475  Gradient gxmym = g->component(var_component,
1476  nxmym,
1477  system.time);
1478  Gradient gxmyp = g->component(var_component,
1479  nxmyp,
1480  system.time);
1481  Gradient gxpym = g->component(var_component,
1482  nxpym,
1483  system.time);
1484  Gradient gxpyp = g->component(var_component,
1485  nxpyp,
1486  system.time);
1487  Number gxzplus = (gxpyp(2) - gxmyp(2))
1488  / 2. / TOLERANCE;
1489  Number gxzminus = (gxpym(2) - gxmym(2))
1490  / 2. / TOLERANCE;
1491  // xyz derivative
1492  Ue(current_dof) = (gxzplus - gxzminus)
1493  / 2. / TOLERANCE;
1494  dof_is_fixed[current_dof] = true;
1495  current_dof++;
1496  }
1497  }
1498  }
1499  // Assume that other C_ONE elements have a single nodal
1500  // value shape function and nodal gradient component
1501  // shape functions
1502  else if (cont == C_ONE)
1503  {
1504  libmesh_assert_equal_to (nc, 1 + dim);
1505  Ue(current_dof) = f->component(var_component,
1506  elem->point(n),
1507  system.time);
1508  dof_is_fixed[current_dof] = true;
1509  current_dof++;
1510  Gradient grad = g->component(var_component,
1511  elem->point(n),
1512  system.time);
1513  for (unsigned int i=0; i!= dim; ++i)
1514  {
1515  Ue(current_dof) = grad(i);
1516  dof_is_fixed[current_dof] = true;
1517  current_dof++;
1518  }
1519  }
1520  else
1521  libmesh_error_msg("Unknown continuity " << cont);
1522  }
1523 
1524  // In 3D, project any edge values next
1525  if (dim > 2 && cont != DISCONTINUOUS)
1526  for (unsigned short e = 0; e != n_edges; ++e)
1527  {
1528  if (!is_boundary_edge[e])
1529  continue;
1530 
1531  FEInterface::dofs_on_edge(elem, dim, fe_type, e,
1532  side_dofs);
1533 
1534  const unsigned int n_side_dofs =
1535  cast_int<unsigned int>(side_dofs.size());
1536 
1537  // Some edge dofs are on nodes and already
1538  // fixed, others are free to calculate
1539  unsigned int free_dofs = 0;
1540  for (auto i : IntRange<unsigned int>(0, n_side_dofs))
1541  if (!dof_is_fixed[side_dofs[i]])
1542  free_dof[free_dofs++] = i;
1543 
1544  // There may be nothing to project
1545  if (!free_dofs)
1546  continue;
1547 
1548  Ke.resize (free_dofs, free_dofs); Ke.zero();
1549  Fe.resize (free_dofs); Fe.zero();
1550  // The new edge coefficients
1551  DenseVector<Number> Uedge(free_dofs);
1552 
1553  // Initialize FE data on the edge
1554  fe->attach_quadrature_rule (qedgerule.get());
1555  fe->edge_reinit (elem, e);
1556  const unsigned int n_qp = qedgerule->n_points();
1557 
1558  // Loop over the quadrature points
1559  for (unsigned int qp=0; qp<n_qp; qp++)
1560  {
1561  // solution at the quadrature point
1562  Number fineval = f->component(var_component,
1563  xyz_values[qp],
1564  system.time);
1565  // solution grad at the quadrature point
1566  Gradient finegrad;
1567  if (cont == C_ONE)
1568  finegrad = g->component(var_component,
1569  xyz_values[qp],
1570  system.time);
1571 
1572  // Form edge projection matrix
1573  for (unsigned int sidei=0, freei=0;
1574  sidei != n_side_dofs; ++sidei)
1575  {
1576  unsigned int i = side_dofs[sidei];
1577  // fixed DoFs aren't test functions
1578  if (dof_is_fixed[i])
1579  continue;
1580  for (unsigned int sidej=0, freej=0;
1581  sidej != n_side_dofs; ++sidej)
1582  {
1583  unsigned int j = side_dofs[sidej];
1584  if (dof_is_fixed[j])
1585  Fe(freei) -= phi[i][qp] * phi[j][qp] *
1586  JxW[qp] * Ue(j);
1587  else
1588  Ke(freei,freej) += phi[i][qp] *
1589  phi[j][qp] * JxW[qp];
1590  if (cont == C_ONE)
1591  {
1592  if (dof_is_fixed[j])
1593  Fe(freei) -= ((*dphi)[i][qp] *
1594  (*dphi)[j][qp]) *
1595  JxW[qp] * Ue(j);
1596  else
1597  Ke(freei,freej) += ((*dphi)[i][qp] *
1598  (*dphi)[j][qp])
1599  * JxW[qp];
1600  }
1601  if (!dof_is_fixed[j])
1602  freej++;
1603  }
1604  Fe(freei) += phi[i][qp] * fineval * JxW[qp];
1605  if (cont == C_ONE)
1606  Fe(freei) += (finegrad * (*dphi)[i][qp]) *
1607  JxW[qp];
1608  freei++;
1609  }
1610  }
1611 
1612  Ke.cholesky_solve(Fe, Uedge);
1613 
1614  // Transfer new edge solutions to element
1615  for (unsigned int i=0; i != free_dofs; ++i)
1616  {
1617  Number & ui = Ue(side_dofs[free_dof[i]]);
1618  libmesh_assert(std::abs(ui) < TOLERANCE ||
1619  std::abs(ui - Uedge(i)) < TOLERANCE);
1620  ui = Uedge(i);
1621  dof_is_fixed[side_dofs[free_dof[i]]] = true;
1622  }
1623  }
1624 
1625  // Project any side values (edges in 2D, faces in 3D)
1626  if (dim > 1 && cont != DISCONTINUOUS)
1627  for (unsigned short s = 0; s != n_sides; ++s)
1628  {
1629  if (!is_boundary_side[s])
1630  continue;
1631 
1632  FEInterface::dofs_on_side(elem, dim, fe_type, s,
1633  side_dofs);
1634 
1635  // Some side dofs are on nodes/edges and already
1636  // fixed, others are free to calculate
1637  unsigned int free_dofs = 0;
1638  for (auto i : index_range(side_dofs))
1639  if (!dof_is_fixed[side_dofs[i]])
1640  free_dof[free_dofs++] = i;
1641 
1642  // There may be nothing to project
1643  if (!free_dofs)
1644  continue;
1645 
1646  Ke.resize (free_dofs, free_dofs); Ke.zero();
1647  Fe.resize (free_dofs); Fe.zero();
1648  // The new side coefficients
1649  DenseVector<Number> Uside(free_dofs);
1650 
1651  // Initialize FE data on the side
1652  fe->attach_quadrature_rule (qsiderule.get());
1653  fe->reinit (elem, s);
1654  const unsigned int n_qp = qsiderule->n_points();
1655 
1656  const unsigned int n_side_dofs =
1657  cast_int<unsigned int>(side_dofs.size());
1658 
1659  // Loop over the quadrature points
1660  for (unsigned int qp=0; qp<n_qp; qp++)
1661  {
1662  // solution at the quadrature point
1663  Number fineval = f->component(var_component,
1664  xyz_values[qp],
1665  system.time);
1666  // solution grad at the quadrature point
1667  Gradient finegrad;
1668  if (cont == C_ONE)
1669  finegrad = g->component(var_component,
1670  xyz_values[qp],
1671  system.time);
1672 
1673  // Form side projection matrix
1674  for (unsigned int sidei=0, freei=0;
1675  sidei != n_side_dofs; ++sidei)
1676  {
1677  unsigned int i = side_dofs[sidei];
1678  // fixed DoFs aren't test functions
1679  if (dof_is_fixed[i])
1680  continue;
1681  for (unsigned int sidej=0, freej=0;
1682  sidej != n_side_dofs; ++sidej)
1683  {
1684  unsigned int j = side_dofs[sidej];
1685  if (dof_is_fixed[j])
1686  Fe(freei) -= phi[i][qp] * phi[j][qp] *
1687  JxW[qp] * Ue(j);
1688  else
1689  Ke(freei,freej) += phi[i][qp] *
1690  phi[j][qp] * JxW[qp];
1691  if (cont == C_ONE)
1692  {
1693  if (dof_is_fixed[j])
1694  Fe(freei) -= ((*dphi)[i][qp] *
1695  (*dphi)[j][qp]) *
1696  JxW[qp] * Ue(j);
1697  else
1698  Ke(freei,freej) += ((*dphi)[i][qp] *
1699  (*dphi)[j][qp])
1700  * JxW[qp];
1701  }
1702  if (!dof_is_fixed[j])
1703  freej++;
1704  }
1705  Fe(freei) += (fineval * phi[i][qp]) * JxW[qp];
1706  if (cont == C_ONE)
1707  Fe(freei) += (finegrad * (*dphi)[i][qp]) *
1708  JxW[qp];
1709  freei++;
1710  }
1711  }
1712 
1713  Ke.cholesky_solve(Fe, Uside);
1714 
1715  // Transfer new side solutions to element
1716  for (unsigned int i=0; i != free_dofs; ++i)
1717  {
1718  Number & ui = Ue(side_dofs[free_dof[i]]);
1719  libmesh_assert(std::abs(ui) < TOLERANCE ||
1720  std::abs(ui - Uside(i)) < TOLERANCE);
1721  ui = Uside(i);
1722  dof_is_fixed[side_dofs[free_dof[i]]] = true;
1723  }
1724  }
1725 
1726  const dof_id_type
1727  first = new_vector.first_local_index(),
1728  last = new_vector.last_local_index();
1729 
1730  // Lock the new_vector since it is shared among threads.
1731  {
1732  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1733 
1734  for (unsigned int i = 0; i < n_dofs; i++)
1735  if (dof_is_fixed[i] &&
1736  (dof_indices[i] >= first) &&
1737  (dof_indices[i] < last))
1738  {
1739  new_vector.set(dof_indices[i], Ue(i));
1740  }
1741  }
1742  } // end elem loop
1743  } // end variables loop
1744 }
1745 
1746 
1747 } // namespace libMesh
Manages the family, order, etc. parameters for a given FE.
Definition: fe_type.h:179
FEFamily family
Definition: fe_type.h:204
void eval_old_dofs(const FEMContext &c, unsigned int var, std::vector< DSNA > &values)
Wrap a libMesh-style function pointer into a FunctionBase object.
double abs(double a)
dof_id_type dof_number(const unsigned int s, const unsigned int var, const unsigned int comp) const
Definition: dof_object.h:833
const Elem * parent() const
Definition: elem.h:2479
A geometric point in (x,y,z) space associated with a DOF.
Definition: node.h:52
virtual void pre_fe_reinit(const System &, const Elem *e)
Definition: fem_context.C:1552
unsigned int n_comp(const unsigned int s, const unsigned int var) const
Definition: dof_object.h:803
BuildProjectionList(const System &system_in)
SparseMatrix< ValOut > & target_matrix
unsigned int n_var_groups(const unsigned int s) const
Definition: dof_object.h:758
virtual void get(const std::vector< numeric_index_type > &index, T *values) const
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Definition: dof_map.C:1930
void parallel_for(const Range &range, const Body &body)
Definition: threads_none.h:73
std::unique_ptr< FunctionBase< Number > > f
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)
virtual numeric_index_type size() const =0
void resize(const unsigned int n)
Definition: dense_vector.h:355
DynamicSparseNumberArray< Output, dof_id_type > type
The base class for all geometric element types.
Definition: elem.h:100
IntRange< std::size_t > index_range(const std::vector< T > &vec)
Definition: int_range.h:104
virtual std::unique_ptr< NumericVector< T > > clone() const =0
static const Real TOLERANCE
Utility class for defining generic ranges for threading.
Definition: stored_range.h:52
const unsigned int n_vars
Definition: tecplot_io.C:69
const std::vector< unsigned int > & variables
virtual void init(const numeric_index_type n, const numeric_index_type n_local, const bool fast=false, const ParallelType ptype=AUTOMATIC)=0
boundary_id_type bc_id
Definition: xdr_io.C:51
void old_dof_indices(const Elem *const elem, std::vector< dof_id_type > &di, const unsigned int vn=libMesh::invalid_uint) const
Definition: dof_map.C:2434
void iota(ForwardIter first, ForwardIter last, T value)
Definition: utility.h:57
void SCALAR_dof_indices(std::vector< dof_id_type > &di, const unsigned int vn, const bool old_dofs=false) const
Definition: dof_map.C:2348
virtual numeric_index_type row_stop() const =0
virtual void set(const numeric_index_type i, const numeric_index_type j, const T value)=0
std::vector< boundary_id_type > boundary_ids(const Node *node) const
BuildProjectionList(BuildProjectionList &other, Threads::split)
Manages the degrees of freedom (DOFs) in a simulation.
Definition: dof_map.h:176
std::unique_ptr< QBase > default_quadrature_rule(const unsigned int dim, const int extraorder=0) const
Definition: fe_type.C:31
const std::vector< std::vector< OutputGradient > > & get_dphi() const
Definition: fe_base.h:215
const dof_id_type n_nodes
Definition: tecplot_io.C:68
spin_mutex spin_mtx
Definition: threads.C:29
A variable which is solved for in a System of equations.
Definition: variable.h:49
const Variable & variable(const unsigned int c) const
Definition: dof_map.h:1762
dof_id_type numeric_index_type
Definition: id_types.h:92
OldSolutionCoefs(const libMesh::System &sys_in)
NumericVector< Number > & new_vector
unsigned int n_vars(const unsigned int s, const unsigned int vg) const
Definition: dof_object.h:768
VectorValue< DynamicSparseNumberArray< InnerOutput, dof_id_type > > type
Manages consistently variables, degrees of freedom, and coefficient vectors.
Definition: system.h:92
MetaPhysicL::DynamicSparseNumberArray< typename CompareTypes< T, T2 >::supertype, I > supertype
void _cholesky_back_substitute(const DenseVector< T2 > &b, DenseVector< T2 > &x) const
Used by the Mesh to keep track of boundary nodes and elements.
Definition: boundary_info.h:57
dof_id_type end_old_dof(const processor_id_type proc) const
Definition: dof_map.h:664
const std::vector< dof_id_type > & get_dof_indices() const
Definition: diff_context.h:367
bool active_on_subdomain(subdomain_id_type sid) const
Definition: variable.h:136
BoundaryProjectSolution(const BoundaryProjectSolution &in)
virtual void close()=0
static std::unique_ptr< NumericVector< Number > > build(const Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package())
OldSolutionCoefs(const OldSolutionCoefs &in)
ParallelType type() const
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
DofObject * old_dof_object
Definition: dof_object.h:79
DynamicSparseNumberArray< Real, dof_id_type > DSNAN
std::vector< dof_id_type > send_list
virtual numeric_index_type row_start() const =0
virtual void zero() override
Definition: dense_matrix.h:808
void swap(Iterator &lhs, Iterator &rhs)
virtual numeric_index_type local_size() const =0
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)
Definition: dense_matrix.h:792
void parallel_reduce(const Range &range, Body &body)
Definition: threads_none.h:101
MatrixFillAction(SparseMatrix< ValOut > &target_mat)
virtual void set(const numeric_index_type i, const T value)=0
void cholesky_solve(const DenseVector< T2 > &b, DenseVector< T2 > &x)
void insert(const FEMContext &c, unsigned int var_num, const DenseVector< DynamicSparseNumberArray< ValIn, dof_id_type > > &Ue)
const std::set< boundary_id_type > & b
std::unique_ptr< FunctionBase< Gradient > > g
unsigned int n_comp_group(const unsigned int s, const unsigned int vg) const
Definition: dof_object.h:816
static const bool value
Definition: compare_types.h:58
A matrix object used for finite element assembly and numerics.
Definition: dense_matrix.h:54
A geometric point in (x,y,z) space.
Definition: point.h:38
DSNAOutput< Output >::type DSNA
const Elem & get(const ElemType type_in)
dof_id_type first_old_dof(const processor_id_type proc) const
Definition: dof_map.h:609
virtual void zero() override
Definition: dense_vector.h:379
uint8_t dof_id_type
Definition: id_types.h:64
const std::vector< std::vector< OutputShape > > & get_phi() const
Definition: fe_base.h:207
const FEType & type() const
Definition: variable.h:119
std::vector< boundary_id_type > edge_boundary_ids(const Elem *const elem, const unsigned short int edge) const
virtual void localize(std::vector< T > &v_local) const =0