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