libMesh
Public Member Functions | Private Attributes | List of all members
libMesh::GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction > Class Template Reference

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. More...

#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
 Function definitions. More...
 

Private Attributes

const System & system
 
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

◆ GenericProjector() [1/2]

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 
)

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  {}
libMesh::GenericProjector::master_action
const ProjectionAction & master_action
Definition: generic_projector.h:64
libMesh::GenericProjector::variables
const std::vector< unsigned int > & variables
Definition: generic_projector.h:65

◆ GenericProjector() [2/2]

template<typename FFunctor , typename GFunctor , typename FValue , typename ProjectionAction >
libMesh::GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction >::GenericProjector ( const GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction > &  in)

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) : nullptr),
85  g_was_copied(in.master_g),
86  master_action(in.master_action),
87  variables(in.variables)
88  {}
libMesh::GenericProjector::master_action
const ProjectionAction & master_action
Definition: generic_projector.h:64
libMesh::GenericProjector::variables
const std::vector< unsigned int > & variables
Definition: generic_projector.h:65

◆ ~GenericProjector()

template<typename FFunctor , typename GFunctor , typename FValue , typename ProjectionAction >
libMesh::GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction >::~GenericProjector ( )

Definition at line 90 of file generic_projector.h.

References libMesh::GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction >::g_was_copied, and libMesh::GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction >::master_g.

91  {
92  if (g_was_copied)
93  delete master_g;
94  }

Member Function Documentation

◆ operator()()

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 588 of file generic_projector.h.

References std::abs(), libMesh::Variable::active_on_subdomain(), libMesh::FEGenericBase< OutputType >::build(), libMesh::C_ONE, libMesh::C_ZERO, libMesh::DenseMatrix< T >::cholesky_solve(), libMesh::CONSTANT, libMesh::FEType::default_quadrature_rule(), libMesh::DISCONTINUOUS, libMesh::FEInterface::dofs_on_edge(), libMesh::FEInterface::dofs_on_side(), libMesh::FEMContext::edge, libMesh::FEMContext::edge_fe_reinit(), libMesh::FEMContext::elem_dimensions(), libMesh::FEMContext::elem_fe_reinit(), libMesh::FEType::family, libMesh::FEAbstract::get_continuity(), libMesh::DiffContext::get_dof_indices(), 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::FEInterface::inverse_map(), libMesh::Elem::JUST_COARSENED, libMesh::Elem::JUST_REFINED, libMesh::LAGRANGE, libMesh::MONOMIAL, libMesh::FEInterface::n_dofs_at_node(), libMesh::MeshTools::n_nodes(), libMesh::FEType::order, libMesh::FEMContext::pre_fe_reinit(), libMesh::Real, libMesh::DenseVector< T >::resize(), libMesh::DenseMatrix< T >::resize(), libMesh::SCALAR, libMesh::FEMContext::side, libMesh::FEMContext::side_fe_reinit(), libMesh::DiffContext::time, libMesh::TOLERANCE, libMesh::Variable::type(), libMesh::DofMap::variable(), libMesh::DenseMatrix< T >::zero(), and libMesh::DenseVector< T >::zero().

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

Member Data Documentation

◆ g_was_copied

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.

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

◆ master_action

template<typename FFunctor , typename GFunctor , typename FValue , typename ProjectionAction >
const ProjectionAction& libMesh::GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction >::master_action
private

Definition at line 64 of file generic_projector.h.

◆ master_f

template<typename FFunctor , typename GFunctor , typename FValue , typename ProjectionAction >
const FFunctor& libMesh::GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction >::master_f
private

Definition at line 61 of file generic_projector.h.

◆ master_g

template<typename FFunctor , typename GFunctor , typename FValue , typename ProjectionAction >
const GFunctor* libMesh::GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction >::master_g
private

Definition at line 62 of file generic_projector.h.

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

◆ system

template<typename FFunctor , typename GFunctor , typename FValue , typename ProjectionAction >
const System& libMesh::GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction >::system
private

Definition at line 57 of file generic_projector.h.

◆ variables

template<typename FFunctor , typename GFunctor , typename FValue , typename ProjectionAction >
const std::vector<unsigned int>& libMesh::GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction >::variables
private

Definition at line 65 of file generic_projector.h.

The documentation for this class was generated from the following file:
Generated on Thu Dec 13 2018 19:33:41 for libMesh by   doxygen 1.8.14