libMesh::GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction > Class Template Reference

#include <generic_projector.h>

Public Member Functions

 GenericProjector (const System &system_in, const FFunctor &f_in, const GFunctor *g_in, const ProjectionAction &act_in, const std::vector< unsigned int > &variables_in)
 
 GenericProjector (const GenericProjector &in)
 
 ~GenericProjector ()
 
void operator() (const ConstElemRange &range) const
 

Private Attributes

const Systemsystem
 
const FFunctor & master_f
 
const GFunctor * master_g
 
bool g_was_copied
 
const ProjectionAction & master_action
 
const std::vector< unsigned int > & variables
 

Detailed Description

template<typename FFunctor, typename GFunctor, typename FValue, typename ProjectionAction>
class libMesh::GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction >

The GenericProjector class implements the core of other projection operations, using two input functors to read values to be projected and an output functor to set degrees of freedom in the result.

This may be executed in parallel on multiple threads.

Author
Roy H. Stogner
Date
2016

Definition at line 54 of file generic_projector.h.

Constructor & Destructor Documentation

template<typename FFunctor , typename GFunctor , typename FValue , typename ProjectionAction >
libMesh::GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction >::GenericProjector ( const System system_in,
const FFunctor &  f_in,
const GFunctor *  g_in,
const ProjectionAction &  act_in,
const std::vector< unsigned int > &  variables_in 
)
inline

Definition at line 68 of file generic_projector.h.

72  :
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  {}
const ProjectionAction & master_action
const std::vector< unsigned int > & variables
template<typename FFunctor , typename GFunctor , typename FValue , typename ProjectionAction >
libMesh::GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction >::GenericProjector ( const GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction > &  in)
inline

Definition at line 81 of file generic_projector.h.

81  :
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  {}
const class libmesh_nullptr_t libmesh_nullptr
const ProjectionAction & master_action
const std::vector< unsigned int > & variables
template<typename FFunctor , typename GFunctor , typename FValue , typename ProjectionAction >
libMesh::GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction >::~GenericProjector ( )
inline

Member Function Documentation

template<typename FFunctor , typename GFunctor , typename FValue , typename ProjectionAction >
void libMesh::GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction >::operator() ( const ConstElemRange range) const

Function definitions

Definition at line 567 of file generic_projector.h.

References std::abs(), libMesh::Variable::active_on_subdomain(), libMesh::FEGenericBase< OutputType >::build(), libMesh::C_ONE, libMesh::C_ZERO, libMesh::Elem::child_ptr(), libMesh::Elem::child_ref_range(), libMesh::DenseMatrix< T >::cholesky_solve(), libMesh::CONSTANT, libMesh::FEType::default_quadrature_rule(), libMesh::Elem::dim(), libMesh::DISCONTINUOUS, libMesh::FEInterface::dofs_on_edge(), libMesh::FEInterface::dofs_on_side(), libMesh::FEMContext::edge, libMesh::FEMContext::edge_fe_reinit(), libMesh::Elem::edge_index_range(), libMesh::FEMContext::elem_dimensions(), libMesh::FEMContext::elem_fe_reinit(), libMesh::FEType::family, libMesh::FEAbstract::get_continuity(), libMesh::DiffContext::get_dof_indices(), libMesh::System::get_dof_map(), libMesh::FEGenericBase< OutputType >::get_dphi(), libMesh::FEMContext::get_edge_fe(), libMesh::FEMContext::get_element_fe(), libMesh::FEAbstract::get_fe_type(), libMesh::FEAbstract::get_JxW(), libMesh::FEGenericBase< OutputType >::get_phi(), libMesh::FEMContext::get_side_fe(), libMesh::DenseVector< T >::get_values(), libMesh::FEAbstract::get_xyz(), libMesh::HERMITE, libMesh::Elem::hmin(), libMesh::Elem::infinite(), libMesh::FEInterface::inverse_map(), libMesh::Elem::is_child_on_edge(), libMesh::Elem::is_child_on_side(), libMesh::Elem::is_vertex(), libMesh::Elem::is_vertex_on_parent(), libMesh::Elem::JUST_COARSENED, libMesh::Elem::JUST_REFINED, libMesh::LAGRANGE, libmesh_nullptr, libMesh::GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction >::master_action, libMesh::GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction >::master_f, libMesh::GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction >::master_g, libMesh::MONOMIAL, libMesh::Elem::n_children(), libMesh::FEInterface::n_dofs_at_node(), n_nodes, libMesh::Elem::n_nodes(), libMesh::Elem::node_ref(), libMesh::DofObject::old_dof_object, libMesh::FEType::order, libMesh::Elem::p_level(), libMesh::Elem::p_refinement_flag(), libMesh::Elem::parent(), libMesh::Elem::point(), libMesh::FEMContext::pre_fe_reinit(), libMesh::Real, libMesh::Elem::refinement_flag(), libMesh::DenseVector< T >::resize(), libMesh::DenseMatrix< T >::resize(), libMesh::SCALAR, libMesh::FEMContext::side, libMesh::FEMContext::side_fe_reinit(), libMesh::Elem::side_index_range(), libMesh::Elem::subdomain_id(), libMesh::GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction >::system, libMesh::System::time, libMesh::TOLERANCE, libMesh::Variable::type(), libMesh::Elem::type(), libMesh::DofMap::variable(), libMesh::System::variable_scalar_number(), libMesh::GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction >::variables, libMesh::Elem::which_child_am_i(), libMesh::DenseMatrix< T >::zero(), and libMesh::DenseVector< T >::zero().

Referenced by libMesh::GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction >::~GenericProjector().

568 {
569  LOG_SCOPE ("operator()", "GenericProjector");
570 
571  ProjectionAction action(master_action);
572  FFunctor f(master_f);
573  std::unique_ptr<GFunctor> g;
574  if (master_g)
575  g.reset(new GFunctor(*master_g));
576 
577  // The DofMap for this system
578  const DofMap & dof_map = system.get_dof_map();
579 
580  // The element matrix and RHS for projections.
581  // Note that Ke is always real-valued, whereas
582  // Fe may be complex valued if complex number
583  // support is enabled
584  DenseMatrix<Real> Ke;
585  DenseVector<FValue> Fe;
586  // The new element degree of freedom coefficients
587  DenseVector<FValue> Ue;
588 
589  // Context objects to contain all our required FE objects
590  FEMContext context( system );
591 
592  // Loop over all the variables we've been requested to project, to
593  // pre-request
594  for (std::size_t v=0; v!=variables.size(); v++)
595  {
596  const unsigned int var = variables[v];
597 
598  // FIXME: Need to generalize this to vector-valued elements. [PB]
599  FEBase * fe = libmesh_nullptr;
600  FEBase * side_fe = libmesh_nullptr;
601  FEBase * edge_fe = libmesh_nullptr;
602 
603  const std::set<unsigned char> & elem_dims =
604  context.elem_dimensions();
605 
606  for (std::set<unsigned char>::const_iterator dim_it =
607  elem_dims.begin(); dim_it != elem_dims.end(); ++dim_it)
608  {
609  const unsigned char dim = *dim_it;
610 
611  context.get_element_fe( var, fe, dim );
612  if (fe->get_fe_type().family == SCALAR)
613  continue;
614  if (dim > 1)
615  context.get_side_fe( var, side_fe, dim );
616  if (dim > 2)
617  context.get_edge_fe( var, edge_fe );
618 
619  fe->get_xyz();
620 
621  fe->get_phi();
622  if (dim > 1)
623  side_fe->get_phi();
624  if (dim > 2)
625  edge_fe->get_phi();
626 
627  const FEContinuity cont = fe->get_continuity();
628  if (cont == C_ONE)
629  {
630  // Our C1 elements need gradient information
631  libmesh_assert(g);
632 
633  fe->get_dphi();
634  if (dim > 1)
635  side_fe->get_dphi();
636  if (dim > 2)
637  edge_fe->get_dphi();
638  }
639  }
640  }
641 
642  // Now initialize any user requested shape functions, xyz vals, etc
643  f.init_context(context);
644  if (g.get())
645  g->init_context(context);
646 
647  // this->init_context(context);
648 
649  // Iterate over all the elements in the range
650  for (ConstElemRange::const_iterator elem_it=range.begin(); elem_it != range.end();
651  ++elem_it)
652  {
653  const Elem * elem = *elem_it;
654 
655  unsigned int dim = elem->dim();
656 
657  context.pre_fe_reinit(system, elem);
658 
659  // If we're doing AMR, this might be a grid projection with a cheap
660  // early exit.
661 #ifdef LIBMESH_ENABLE_AMR
662  // If this element doesn't have an old_dof_object, but it
663  // wasn't just refined or just coarsened into activity, then
664  // it must be newly added, so the user is responsible for
665  // setting the new dofs on it during a grid projection.
666  if (!elem->old_dof_object &&
667  elem->refinement_flag() != Elem::JUST_REFINED &&
668  elem->refinement_flag() != Elem::JUST_COARSENED &&
669  f.is_grid_projection())
670  continue;
671 #endif // LIBMESH_ENABLE_AMR
672 
673  // Loop over all the variables we've been requested to project, to
674  // do the projection
675  for (std::size_t v=0; v!=variables.size(); v++)
676  {
677  const unsigned int var = variables[v];
678 
679  const Variable & variable = dof_map.variable(var);
680 
681  const FEType & base_fe_type = variable.type();
682 
683  FEType fe_type = base_fe_type;
684 
685  // This may be a p refined element
686  fe_type.order =
687  libMesh::Order (fe_type.order + elem->p_level());
688 
689  if (fe_type.family == SCALAR)
690  continue;
691 
692  // Per-subdomain variables don't need to be projected on
693  // elements where they're not active
694  if (!variable.active_on_subdomain(elem->subdomain_id()))
695  continue;
696 
697  const std::vector<dof_id_type> & dof_indices =
698  context.get_dof_indices(var);
699 
700  // The number of DOFs on the element
701  const unsigned int n_dofs =
702  cast_int<unsigned int>(dof_indices.size());
703 
704  const unsigned int var_component =
706 
707  // Zero the interpolated values
708  Ue.resize (n_dofs); Ue.zero();
709 
710  // If we're projecting from an old grid
711 #ifdef LIBMESH_ENABLE_AMR
712  if (f.is_grid_projection() &&
713  // And either this is an unchanged element
714  ((elem->refinement_flag() != Elem::JUST_REFINED &&
715  elem->refinement_flag() != Elem::JUST_COARSENED &&
716  elem->p_refinement_flag() != Elem::JUST_REFINED &&
717  elem->p_refinement_flag() != Elem::JUST_COARSENED) ||
718  // Or this is a low order monomial element which has merely
719  // been h refined
720  (fe_type.family == MONOMIAL &&
721  fe_type.order == CONSTANT &&
722  elem->p_level() == 0 &&
723  elem->refinement_flag() != Elem::JUST_COARSENED &&
724  elem->p_refinement_flag() != Elem::JUST_COARSENED))
725  )
726  // then we can simply copy its old dof
727  // values to new indices.
728  {
729  LOG_SCOPE ("copy_dofs", "GenericProjector");
730 
731  f.eval_old_dofs(context, var_component, Ue.get_values());
732 
733  action.insert(context, var, Ue);
734 
735  continue;
736  }
737 #endif // LIBMESH_ENABLE_AMR
738 
739  FEBase * fe = libmesh_nullptr;
740  FEBase * side_fe = libmesh_nullptr;
741  FEBase * edge_fe = libmesh_nullptr;
742 
743  context.get_element_fe( var, fe, dim );
744  if (fe->get_fe_type().family == SCALAR)
745  continue;
746  if (dim > 1)
747  context.get_side_fe( var, side_fe, dim );
748  if (dim > 2)
749  context.get_edge_fe( var, edge_fe );
750 
751  const FEContinuity cont = fe->get_continuity();
752 
753  std::vector<unsigned int> side_dofs;
754 
755  // Fixed vs. free DoFs on edge/face projections
756  std::vector<char> dof_is_fixed(n_dofs, false); // bools
757  std::vector<int> free_dof(n_dofs, 0);
758 
759  // The element type
760  const ElemType elem_type = elem->type();
761 
762  // The number of nodes on the new element
763  const unsigned int n_nodes = elem->n_nodes();
764 
765  START_LOG ("project_nodes","GenericProjector");
766 
767  // When interpolating C1 elements we need to know which
768  // vertices were also parent vertices; we'll cache an
769  // intermediate fact outside the nodes loop.
770  int i_am_child = -1;
771 #ifdef LIBMESH_ENABLE_AMR
772  if (elem->parent())
773  i_am_child = elem->parent()->which_child_am_i(elem);
774 #endif
775  // In general, we need a series of
776  // projections to ensure a unique and continuous
777  // solution. We start by interpolating nodes, then
778  // hold those fixed and project edges, then
779  // hold those fixed and project faces, then
780  // hold those fixed and project interiors
781  //
782  // In the LAGRANGE case, we will save a lot of solution
783  // evaluations (at a slight cost in accuracy) by simply
784  // interpolating all nodes rather than projecting.
785 
786  // Interpolate vertex (or for LAGRANGE, all node) values first.
787  unsigned int current_dof = 0;
788  for (unsigned int n=0; n!= n_nodes; ++n)
789  {
790  // FIXME: this should go through the DofMap,
791  // not duplicate dof_indices code badly!
792  const unsigned int nc =
793  FEInterface::n_dofs_at_node (dim, fe_type, elem_type, n);
794 
795  if (!elem->is_vertex(n) &&
796  fe_type.family != LAGRANGE)
797  {
798  current_dof += nc;
799  continue;
800  }
801 
802  if (cont == DISCONTINUOUS)
803  {
804  libmesh_assert_equal_to (nc, 0);
805  }
806  else if (!nc)
807  {
808  // This should only occur for first-order LAGRANGE
809  // FE on non-vertices of higher-order elements
810  libmesh_assert (!elem->is_vertex(n));
811  libmesh_assert_equal_to(fe_type.family, LAGRANGE);
812  }
813  // Assume that C_ZERO elements have a single nodal
814  // value shape function at vertices
815  else if (cont == C_ZERO)
816  {
817  Ue(current_dof) = f.eval_at_node(context,
818  var_component,
819  dim,
820  elem->node_ref(n),
821  system.time);
822  dof_is_fixed[current_dof] = true;
823  current_dof++;
824  }
825  else if (cont == C_ONE)
826  {
827  const bool is_parent_vertex = (i_am_child == -1) ||
828  elem->parent()->is_vertex_on_parent(i_am_child, n);
829 
830  // The hermite element vertex shape functions are weird
831  if (fe_type.family == HERMITE)
832  {
833  Ue(current_dof) =
834  f.eval_at_node(context,
835  var_component,
836  dim,
837  elem->node_ref(n),
838  system.time);
839  dof_is_fixed[current_dof] = true;
840  current_dof++;
841  VectorValue<FValue> grad =
842  is_parent_vertex ?
843  g->eval_at_node(context,
844  var_component,
845  dim,
846  elem->node_ref(n),
847  system.time) :
848  g->eval_at_point(context,
849  var_component,
850  elem->node_ref(n),
851  system.time);
852  // x derivative
853  Ue(current_dof) = grad(0);
854  dof_is_fixed[current_dof] = true;
855  current_dof++;
856  if (dim > 1)
857  {
858  // We'll finite difference mixed derivatives
859  Real delta_x = TOLERANCE * elem->hmin();
860 
861  Point nxminus = elem->point(n),
862  nxplus = elem->point(n);
863  nxminus(0) -= delta_x;
864  nxplus(0) += delta_x;
865  VectorValue<FValue> gxminus =
866  g->eval_at_point(context,
867  var_component,
868  nxminus,
869  system.time);
870  VectorValue<FValue> gxplus =
871  g->eval_at_point(context,
872  var_component,
873  nxplus,
874  system.time);
875  // y derivative
876  Ue(current_dof) = grad(1);
877  dof_is_fixed[current_dof] = true;
878  current_dof++;
879  // xy derivative
880  Ue(current_dof) = (gxplus(1) - gxminus(1))
881  / 2. / delta_x;
882  dof_is_fixed[current_dof] = true;
883  current_dof++;
884 
885  if (dim > 2)
886  {
887  // z derivative
888  Ue(current_dof) = grad(2);
889  dof_is_fixed[current_dof] = true;
890  current_dof++;
891  // xz derivative
892  Ue(current_dof) = (gxplus(2) - gxminus(2))
893  / 2. / delta_x;
894  dof_is_fixed[current_dof] = true;
895  current_dof++;
896  // We need new points for yz
897  Point nyminus = elem->point(n),
898  nyplus = elem->point(n);
899  nyminus(1) -= delta_x;
900  nyplus(1) += delta_x;
901  VectorValue<FValue> gyminus =
902  g->eval_at_point(context,
903  var_component,
904  nyminus,
905  system.time);
906  VectorValue<FValue> gyplus =
907  g->eval_at_point(context,
908  var_component,
909  nyplus,
910  system.time);
911  // xz derivative
912  Ue(current_dof) = (gyplus(2) - gyminus(2))
913  / 2. / delta_x;
914  dof_is_fixed[current_dof] = true;
915  current_dof++;
916  // Getting a 2nd order xyz is more tedious
917  Point nxmym = elem->point(n),
918  nxmyp = elem->point(n),
919  nxpym = elem->point(n),
920  nxpyp = elem->point(n);
921  nxmym(0) -= delta_x;
922  nxmym(1) -= delta_x;
923  nxmyp(0) -= delta_x;
924  nxmyp(1) += delta_x;
925  nxpym(0) += delta_x;
926  nxpym(1) -= delta_x;
927  nxpyp(0) += delta_x;
928  nxpyp(1) += delta_x;
929  VectorValue<FValue> gxmym =
930  g->eval_at_point(context,
931  var_component,
932  nxmym,
933  system.time);
934  VectorValue<FValue> gxmyp =
935  g->eval_at_point(context,
936  var_component,
937  nxmyp,
938  system.time);
939  VectorValue<FValue> gxpym =
940  g->eval_at_point(context,
941  var_component,
942  nxpym,
943  system.time);
944  VectorValue<FValue> gxpyp =
945  g->eval_at_point(context,
946  var_component,
947  nxpyp,
948  system.time);
949  FValue gxzplus = (gxpyp(2) - gxmyp(2))
950  / 2. / delta_x;
951  FValue gxzminus = (gxpym(2) - gxmym(2))
952  / 2. / delta_x;
953  // xyz derivative
954  Ue(current_dof) = (gxzplus - gxzminus)
955  / 2. / delta_x;
956  dof_is_fixed[current_dof] = true;
957  current_dof++;
958  }
959  }
960  }
961  else
962  {
963  // Assume that other C_ONE elements have a single nodal
964  // value shape function and nodal gradient component
965  // shape functions
966  libmesh_assert_equal_to (nc, 1 + dim);
967  Ue(current_dof) = f.eval_at_node(context,
968  var_component,
969  dim,
970  elem->node_ref(n),
971  system.time);
972  dof_is_fixed[current_dof] = true;
973  current_dof++;
974  VectorValue<FValue> grad =
975  is_parent_vertex ?
976  g->eval_at_node(context,
977  var_component,
978  dim,
979  elem->node_ref(n),
980  system.time) :
981  g->eval_at_point(context,
982  var_component,
983  elem->node_ref(n),
984  system.time);
985  for (unsigned int i=0; i!= dim; ++i)
986  {
987  Ue(current_dof) = grad(i);
988  dof_is_fixed[current_dof] = true;
989  current_dof++;
990  }
991  }
992  }
993  else
994  libmesh_error_msg("Unknown continuity " << cont);
995  }
996 
997  STOP_LOG ("project_nodes","GenericProjector");
998 
999  START_LOG ("project_edges","GenericProjector");
1000 
1001  // In 3D with non-LAGRANGE, project any edge values next
1002  if (dim > 2 &&
1003  cont != DISCONTINUOUS &&
1004  (fe_type.family != LAGRANGE
1005 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1006  || elem->infinite()
1007 #endif
1008  ))
1009  {
1010  // If we're JUST_COARSENED we'll need a custom
1011  // evaluation, not just the standard edge FE
1012  const std::vector<Point> & xyz_values =
1013 #ifdef LIBMESH_ENABLE_AMR
1014  (elem->refinement_flag() == Elem::JUST_COARSENED) ?
1015  fe->get_xyz() :
1016 #endif
1017  edge_fe->get_xyz();
1018  const std::vector<Real> & JxW =
1019 #ifdef LIBMESH_ENABLE_AMR
1020  (elem->refinement_flag() == Elem::JUST_COARSENED) ?
1021  fe->get_JxW() :
1022 #endif
1023  edge_fe->get_JxW();
1024 
1025  const std::vector<std::vector<Real>> & phi =
1026 #ifdef LIBMESH_ENABLE_AMR
1027  (elem->refinement_flag() == Elem::JUST_COARSENED) ?
1028  fe->get_phi() :
1029 #endif
1030  edge_fe->get_phi();
1031  const std::vector<std::vector<RealGradient>> * dphi = libmesh_nullptr;
1032  if (cont == C_ONE)
1033  dphi =
1034 #ifdef LIBMESH_ENABLE_AMR
1035  (elem->refinement_flag() == Elem::JUST_COARSENED) ?
1036  &(fe->get_dphi()) :
1037 #endif
1038  &(edge_fe->get_dphi());
1039 
1040  for (auto e : elem->edge_index_range())
1041  {
1042  context.edge = e;
1043 
1044 #ifdef LIBMESH_ENABLE_AMR
1045  if (elem->refinement_flag() == Elem::JUST_COARSENED)
1046  {
1047  std::vector<Point> fine_points;
1048 
1049  std::unique_ptr<FEBase> fine_fe
1050  (FEBase::build (dim, base_fe_type));
1051 
1052  std::unique_ptr<QBase> qrule
1053  (base_fe_type.default_quadrature_rule(1));
1054  fine_fe->attach_quadrature_rule(qrule.get());
1055 
1056  const std::vector<Point> & child_xyz =
1057  fine_fe->get_xyz();
1058 
1059  for (unsigned int c = 0, nc = elem->n_children();
1060  c != nc; ++c)
1061  {
1062  if (!elem->is_child_on_edge(c, e))
1063  continue;
1064 
1065  fine_fe->edge_reinit(elem->child_ptr(c), e);
1066  fine_points.insert(fine_points.end(),
1067  child_xyz.begin(),
1068  child_xyz.end());
1069  }
1070 
1071  std::vector<Point> fine_qp;
1072  FEInterface::inverse_map (dim, base_fe_type, elem,
1073  fine_points, fine_qp);
1074 
1075  context.elem_fe_reinit(&fine_qp);
1076  }
1077  else
1078 #endif // LIBMESH_ENABLE_AMR
1079  context.edge_fe_reinit();
1080 
1081  const unsigned int n_qp = xyz_values.size();
1082 
1083  FEInterface::dofs_on_edge(elem, dim, base_fe_type,
1084  e, side_dofs);
1085 
1086  // Some edge dofs are on nodes and already
1087  // fixed, others are free to calculate
1088  unsigned int free_dofs = 0;
1089  for (std::size_t i=0; i != side_dofs.size(); ++i)
1090  if (!dof_is_fixed[side_dofs[i]])
1091  free_dof[free_dofs++] = i;
1092 
1093  // There may be nothing to project
1094  if (!free_dofs)
1095  continue;
1096 
1097  Ke.resize (free_dofs, free_dofs); Ke.zero();
1098  Fe.resize (free_dofs); Fe.zero();
1099  // The new edge coefficients
1100  DenseVector<FValue> Uedge(free_dofs);
1101 
1102  // Loop over the quadrature points
1103  for (unsigned int qp=0; qp<n_qp; qp++)
1104  {
1105  // solution at the quadrature point
1106  FValue fineval = f.eval_at_point(context,
1107  var_component,
1108  xyz_values[qp],
1109  system.time);
1110  // solution grad at the quadrature point
1111  VectorValue<FValue> finegrad;
1112  if (cont == C_ONE)
1113  finegrad = g->eval_at_point(context,
1114  var_component,
1115  xyz_values[qp],
1116  system.time);
1117 
1118  // Form edge projection matrix
1119  for (std::size_t sidei=0, freei=0;
1120  sidei != side_dofs.size(); ++sidei)
1121  {
1122  unsigned int i = side_dofs[sidei];
1123  // fixed DoFs aren't test functions
1124  if (dof_is_fixed[i])
1125  continue;
1126  for (std::size_t sidej=0, freej=0;
1127  sidej != side_dofs.size(); ++sidej)
1128  {
1129  unsigned int j = side_dofs[sidej];
1130  if (dof_is_fixed[j])
1131  Fe(freei) -= phi[i][qp] * phi[j][qp] *
1132  JxW[qp] * Ue(j);
1133  else
1134  Ke(freei,freej) += phi[i][qp] *
1135  phi[j][qp] * JxW[qp];
1136  if (cont == C_ONE)
1137  {
1138  if (dof_is_fixed[j])
1139  Fe(freei) -= ( (*dphi)[i][qp] *
1140  (*dphi)[j][qp] ) *
1141  JxW[qp] * Ue(j);
1142  else
1143  Ke(freei,freej) += ( (*dphi)[i][qp] *
1144  (*dphi)[j][qp] )
1145  * JxW[qp];
1146  }
1147  if (!dof_is_fixed[j])
1148  freej++;
1149  }
1150  Fe(freei) += phi[i][qp] * fineval * JxW[qp];
1151  if (cont == C_ONE)
1152  Fe(freei) += (finegrad * (*dphi)[i][qp] ) *
1153  JxW[qp];
1154  freei++;
1155  }
1156  }
1157 
1158  Ke.cholesky_solve(Fe, Uedge);
1159 
1160  // Transfer new edge solutions to element
1161  for (unsigned int i=0; i != free_dofs; ++i)
1162  {
1163  FValue & ui = Ue(side_dofs[free_dof[i]]);
1164  libmesh_assert(bool(std::abs(ui) < TOLERANCE) ||
1165  bool(std::abs(ui - Uedge(i)) < TOLERANCE));
1166  ui = Uedge(i);
1167  dof_is_fixed[side_dofs[free_dof[i]]] = true;
1168  }
1169  }
1170  } // end if (dim > 2, !DISCONTINUOUS, !LAGRANGE)
1171 
1172  STOP_LOG ("project_edges","GenericProjector");
1173 
1174  START_LOG ("project_sides","GenericProjector");
1175 
1176  // With non-LAGRANGE, project any side values (edges in 2D,
1177  // faces in 3D) next.
1178  if (dim > 1 &&
1179  cont != DISCONTINUOUS &&
1180  (fe_type.family != LAGRANGE
1181 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1182  || elem->infinite()
1183 #endif
1184  ))
1185  {
1186  // If we're JUST_COARSENED we'll need a custom
1187  // evaluation, not just the standard side FE
1188  const std::vector<Point> & xyz_values =
1189 #ifdef LIBMESH_ENABLE_AMR
1190  (elem->refinement_flag() == Elem::JUST_COARSENED) ?
1191  fe->get_xyz() :
1192 #endif // LIBMESH_ENABLE_AMR
1193  side_fe->get_xyz();
1194  const std::vector<Real> & JxW =
1195 #ifdef LIBMESH_ENABLE_AMR
1196  (elem->refinement_flag() == Elem::JUST_COARSENED) ?
1197  fe->get_JxW() :
1198 #endif // LIBMESH_ENABLE_AMR
1199  side_fe->get_JxW();
1200  const std::vector<std::vector<Real>> & phi =
1201 #ifdef LIBMESH_ENABLE_AMR
1202  (elem->refinement_flag() == Elem::JUST_COARSENED) ?
1203  fe->get_phi() :
1204 #endif // LIBMESH_ENABLE_AMR
1205  side_fe->get_phi();
1206  const std::vector<std::vector<RealGradient>> * dphi = libmesh_nullptr;
1207  if (cont == C_ONE)
1208  dphi =
1209 #ifdef LIBMESH_ENABLE_AMR
1210  (elem->refinement_flag() == Elem::JUST_COARSENED) ? &(fe->get_dphi()) :
1211 #endif // LIBMESH_ENABLE_AMR
1212  &(side_fe->get_dphi());
1213 
1214  for (auto s : elem->side_index_range())
1215  {
1216  FEInterface::dofs_on_side(elem, dim, base_fe_type,
1217  s, side_dofs);
1218 
1219  // Some side dofs are on nodes/edges and already
1220  // fixed, others are free to calculate
1221  unsigned int free_dofs = 0;
1222  for (std::size_t i=0; i != side_dofs.size(); ++i)
1223  if (!dof_is_fixed[side_dofs[i]])
1224  free_dof[free_dofs++] = i;
1225 
1226  // There may be nothing to project
1227  if (!free_dofs)
1228  continue;
1229 
1230  Ke.resize (free_dofs, free_dofs); Ke.zero();
1231  Fe.resize (free_dofs); Fe.zero();
1232  // The new side coefficients
1233  DenseVector<FValue> Uside(free_dofs);
1234 
1235  context.side = s;
1236 
1237 #ifdef LIBMESH_ENABLE_AMR
1238  if (elem->refinement_flag() == Elem::JUST_COARSENED)
1239  {
1240  std::vector<Point> fine_points;
1241 
1242  std::unique_ptr<FEBase> fine_fe
1243  (FEBase::build (dim, base_fe_type));
1244 
1245  std::unique_ptr<QBase> qrule
1246  (base_fe_type.default_quadrature_rule(dim-1));
1247  fine_fe->attach_quadrature_rule(qrule.get());
1248 
1249  const std::vector<Point> & child_xyz =
1250  fine_fe->get_xyz();
1251 
1252  for (unsigned int c = 0, nc = elem->n_children();
1253  c != nc; ++c)
1254  {
1255  if (!elem->is_child_on_side(c, s))
1256  continue;
1257 
1258  fine_fe->reinit(elem->child_ptr(c), s);
1259  fine_points.insert(fine_points.end(),
1260  child_xyz.begin(),
1261  child_xyz.end());
1262  }
1263 
1264  std::vector<Point> fine_qp;
1265  FEInterface::inverse_map (dim, base_fe_type, elem,
1266  fine_points, fine_qp);
1267 
1268  context.elem_fe_reinit(&fine_qp);
1269  }
1270  else
1271 #endif // LIBMESH_ENABLE_AMR
1272  context.side_fe_reinit();
1273 
1274  const unsigned int n_qp = xyz_values.size();
1275 
1276  // Loop over the quadrature points
1277  for (unsigned int qp=0; qp<n_qp; qp++)
1278  {
1279  // solution at the quadrature point
1280  FValue fineval = f.eval_at_point(context,
1281  var_component,
1282  xyz_values[qp],
1283  system.time);
1284  // solution grad at the quadrature point
1285  VectorValue<FValue> finegrad;
1286  if (cont == C_ONE)
1287  finegrad = g->eval_at_point(context,
1288  var_component,
1289  xyz_values[qp],
1290  system.time);
1291 
1292  // Form side projection matrix
1293  for (std::size_t sidei=0, freei=0;
1294  sidei != side_dofs.size(); ++sidei)
1295  {
1296  unsigned int i = side_dofs[sidei];
1297  // fixed DoFs aren't test functions
1298  if (dof_is_fixed[i])
1299  continue;
1300  for (std::size_t sidej=0, freej=0;
1301  sidej != side_dofs.size(); ++sidej)
1302  {
1303  unsigned int j = side_dofs[sidej];
1304  if (dof_is_fixed[j])
1305  Fe(freei) -= phi[i][qp] * phi[j][qp] *
1306  JxW[qp] * Ue(j);
1307  else
1308  Ke(freei,freej) += phi[i][qp] *
1309  phi[j][qp] * JxW[qp];
1310  if (cont == C_ONE)
1311  {
1312  if (dof_is_fixed[j])
1313  Fe(freei) -= ( (*dphi)[i][qp] *
1314  (*dphi)[j][qp] ) *
1315  JxW[qp] * Ue(j);
1316  else
1317  Ke(freei,freej) += ( (*dphi)[i][qp] *
1318  (*dphi)[j][qp] )
1319  * JxW[qp];
1320  }
1321  if (!dof_is_fixed[j])
1322  freej++;
1323  }
1324  Fe(freei) += (fineval * phi[i][qp]) * JxW[qp];
1325  if (cont == C_ONE)
1326  Fe(freei) += (finegrad * (*dphi)[i][qp] ) *
1327  JxW[qp];
1328  freei++;
1329  }
1330  }
1331 
1332  Ke.cholesky_solve(Fe, Uside);
1333 
1334  // Transfer new side solutions to element
1335  for (unsigned int i=0; i != free_dofs; ++i)
1336  {
1337  FValue & ui = Ue(side_dofs[free_dof[i]]);
1338  libmesh_assert(bool(std::abs(ui) < TOLERANCE) ||
1339  bool(std::abs(ui - Uside(i)) < TOLERANCE));
1340  ui = Uside(i);
1341  dof_is_fixed[side_dofs[free_dof[i]]] = true;
1342  }
1343  }
1344  } // end if (dim > 1, !DISCONTINUOUS, !LAGRANGE)
1345 
1346  STOP_LOG ("project_sides","GenericProjector");
1347 
1348  START_LOG ("project_interior","GenericProjector");
1349 
1350  // Project the interior values, finally
1351 
1352  // Some interior dofs are on nodes/edges/sides and
1353  // already fixed, others are free to calculate
1354  unsigned int free_dofs = 0;
1355  for (unsigned int i=0; i != n_dofs; ++i)
1356  if (!dof_is_fixed[i])
1357  free_dof[free_dofs++] = i;
1358 
1359  // Project any remaining (interior) dofs in the non-LAGRANGE
1360  // case.
1361  if (free_dofs &&
1362  (fe_type.family != LAGRANGE
1363 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1364  || elem->infinite()
1365 #endif
1366  ))
1367  {
1368  const std::vector<Point> & xyz_values = fe->get_xyz();
1369  const std::vector<Real> & JxW = fe->get_JxW();
1370 
1371  const std::vector<std::vector<Real>> & phi = fe->get_phi();
1372  const std::vector<std::vector<RealGradient>> * dphi = libmesh_nullptr;
1373  if (cont == C_ONE)
1374  dphi = &(fe->get_dphi());
1375 
1376 #ifdef LIBMESH_ENABLE_AMR
1377  if (elem->refinement_flag() == Elem::JUST_COARSENED)
1378  {
1379  std::vector<Point> fine_points;
1380 
1381  std::unique_ptr<FEBase> fine_fe
1382  (FEBase::build (dim, base_fe_type));
1383 
1384  std::unique_ptr<QBase> qrule
1385  (base_fe_type.default_quadrature_rule(dim));
1386  fine_fe->attach_quadrature_rule(qrule.get());
1387 
1388  const std::vector<Point> & child_xyz =
1389  fine_fe->get_xyz();
1390 
1391  for (auto & child : elem->child_ref_range())
1392  {
1393  fine_fe->reinit(&child);
1394  fine_points.insert(fine_points.end(),
1395  child_xyz.begin(),
1396  child_xyz.end());
1397  }
1398 
1399  std::vector<Point> fine_qp;
1400  FEInterface::inverse_map (dim, base_fe_type, elem,
1401  fine_points, fine_qp);
1402 
1403  context.elem_fe_reinit(&fine_qp);
1404  }
1405  else
1406 #endif // LIBMESH_ENABLE_AMR
1407  context.elem_fe_reinit();
1408 
1409  const unsigned int n_qp = xyz_values.size();
1410 
1411  Ke.resize (free_dofs, free_dofs); Ke.zero();
1412  Fe.resize (free_dofs); Fe.zero();
1413  // The new interior coefficients
1414  DenseVector<FValue> Uint(free_dofs);
1415 
1416  // Loop over the quadrature points
1417  for (unsigned int qp=0; qp<n_qp; qp++)
1418  {
1419  // solution at the quadrature point
1420  FValue fineval = f.eval_at_point(context,
1421  var_component,
1422  xyz_values[qp],
1423  system.time);
1424  // solution grad at the quadrature point
1425  VectorValue<FValue> finegrad;
1426  if (cont == C_ONE)
1427  finegrad = g->eval_at_point(context,
1428  var_component,
1429  xyz_values[qp],
1430  system.time);
1431 
1432  // Form interior projection matrix
1433  for (unsigned int i=0, freei=0; i != n_dofs; ++i)
1434  {
1435  // fixed DoFs aren't test functions
1436  if (dof_is_fixed[i])
1437  continue;
1438  for (unsigned int j=0, freej=0; j != n_dofs; ++j)
1439  {
1440  if (dof_is_fixed[j])
1441  Fe(freei) -= phi[i][qp] * phi[j][qp] * JxW[qp]
1442  * Ue(j);
1443  else
1444  Ke(freei,freej) += phi[i][qp] * phi[j][qp] *
1445  JxW[qp];
1446  if (cont == C_ONE)
1447  {
1448  if (dof_is_fixed[j])
1449  Fe(freei) -= ( (*dphi)[i][qp] *
1450  (*dphi)[j][qp] ) * JxW[qp] *
1451  Ue(j);
1452  else
1453  Ke(freei,freej) += ( (*dphi)[i][qp] *
1454  (*dphi)[j][qp] ) *
1455  JxW[qp];
1456  }
1457  if (!dof_is_fixed[j])
1458  freej++;
1459  }
1460  Fe(freei) += phi[i][qp] * fineval * JxW[qp];
1461  if (cont == C_ONE)
1462  Fe(freei) += (finegrad * (*dphi)[i][qp] ) * JxW[qp];
1463  freei++;
1464  }
1465  }
1466  Ke.cholesky_solve(Fe, Uint);
1467 
1468  // Transfer new interior solutions to element
1469  for (unsigned int i=0; i != free_dofs; ++i)
1470  {
1471  FValue & ui = Ue(free_dof[i]);
1472  libmesh_assert(bool(std::abs(ui) < TOLERANCE) ||
1473  bool(std::abs(ui - Uint(i)) < TOLERANCE));
1474  ui = Uint(i);
1475  dof_is_fixed[free_dof[i]] = true;
1476  }
1477 
1478  } // if there are free interior dofs
1479 
1480  STOP_LOG ("project_interior","GenericProjector");
1481 
1482  // Make sure every DoF got reached!
1483  for (unsigned int i=0; i != n_dofs; ++i)
1484  libmesh_assert(dof_is_fixed[i]);
1485 
1486  action.insert(context, var, Ue);
1487  } // end variables loop
1488  } // end elements loop
1489 }
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)
Definition: fe_interface.C:493
double abs(double a)
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)
Definition: fe_interface.C:480
const class libmesh_nullptr_t libmesh_nullptr
static const Real TOLERANCE
const dof_id_type n_nodes
Definition: tecplot_io.C:67
const DofMap & get_dof_map() const
Definition: system.h:2030
static std::unique_ptr< FEGenericBase > build(const unsigned int dim, const FEType &type)
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:547
const ProjectionAction & master_action
FEGenericBase< Real > FEBase
const std::vector< unsigned int > & variables
std::vector< object_type >::const_iterator const_iterator
Definition: stored_range.h:58
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:430
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
unsigned int variable_scalar_number(const std::string &var, unsigned int component) const
Definition: system.h:2145

Member Data Documentation

template<typename FFunctor , typename GFunctor , typename FValue , typename ProjectionAction >
bool libMesh::GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction >::g_was_copied
private

Definition at line 63 of file generic_projector.h.

template<typename FFunctor , typename GFunctor , typename FValue , typename ProjectionAction >
const ProjectionAction& libMesh::GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction >::master_action
private
template<typename FFunctor , typename GFunctor , typename FValue , typename ProjectionAction >
const FFunctor& libMesh::GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction >::master_f
private
template<typename FFunctor , typename GFunctor , typename FValue , typename ProjectionAction >
const GFunctor* libMesh::GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction >::master_g
private
template<typename FFunctor , typename GFunctor , typename FValue , typename ProjectionAction >
const System& libMesh::GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction >::system
private
template<typename FFunctor , typename GFunctor , typename FValue , typename ProjectionAction >
const std::vector<unsigned int>& libMesh::GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction >::variables
private

The documentation for this class was generated from the following file: