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

This class implements the loops to other projection operations. More...

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 >

This class implements the loops to other projection operations.

This may be executed in parallel on multiple threads.

Definition at line 54 of file system_projection.C.

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 
)

Definition at line 68 of file system_projection.C.

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)

Definition at line 81 of file system_projection.C.

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 ( )

Member Function Documentation

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

Definition at line 1077 of file system_projection.C.

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(), dim, 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::libmesh_assert(), 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::BoundaryProjectSolution::BoundaryProjectSolution(), libMesh::BuildProjectionList::BuildProjectionList(), and libMesh::GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction >::~GenericProjector().

1078 {
1079  LOG_SCOPE ("operator()", "GenericProjector");
1080 
1081  ProjectionAction action(master_action);
1082  FFunctor f(master_f);
1083  UniquePtr<GFunctor> g;
1084  if (master_g)
1085  g.reset(new GFunctor(*master_g));
1086 
1087  // The DofMap for this system
1088  const DofMap & dof_map = system.get_dof_map();
1089 
1090  // The element matrix and RHS for projections.
1091  // Note that Ke is always real-valued, whereas
1092  // Fe may be complex valued if complex number
1093  // support is enabled
1094  DenseMatrix<Real> Ke;
1095  DenseVector<FValue> Fe;
1096  // The new element degree of freedom coefficients
1097  DenseVector<FValue> Ue;
1098 
1099  // Context objects to contain all our required FE objects
1100  FEMContext context( system );
1101 
1102  // Loop over all the variables we've been requested to project, to
1103  // pre-request
1104  for (std::size_t v=0; v!=variables.size(); v++)
1105  {
1106  const unsigned int var = variables[v];
1107 
1108  // FIXME: Need to generalize this to vector-valued elements. [PB]
1109  FEBase * fe = libmesh_nullptr;
1110  FEBase * side_fe = libmesh_nullptr;
1111  FEBase * edge_fe = libmesh_nullptr;
1112 
1113  const std::set<unsigned char> & elem_dims =
1114  context.elem_dimensions();
1115 
1116  for (std::set<unsigned char>::const_iterator dim_it =
1117  elem_dims.begin(); dim_it != elem_dims.end(); ++dim_it)
1118  {
1119  const unsigned char dim = *dim_it;
1120 
1121  context.get_element_fe( var, fe, dim );
1122  if (fe->get_fe_type().family == SCALAR)
1123  continue;
1124  if (dim > 1)
1125  context.get_side_fe( var, side_fe, dim );
1126  if (dim > 2)
1127  context.get_edge_fe( var, edge_fe );
1128 
1129  fe->get_xyz();
1130 
1131  fe->get_phi();
1132  if (dim > 1)
1133  side_fe->get_phi();
1134  if (dim > 2)
1135  edge_fe->get_phi();
1136 
1137  const FEContinuity cont = fe->get_continuity();
1138  if (cont == C_ONE)
1139  {
1140  // Our C1 elements need gradient information
1141  libmesh_assert(g);
1142 
1143  fe->get_dphi();
1144  if (dim > 1)
1145  side_fe->get_dphi();
1146  if (dim > 2)
1147  edge_fe->get_dphi();
1148  }
1149  }
1150  }
1151 
1152  // Now initialize any user requested shape functions, xyz vals, etc
1153  f.init_context(context);
1154  if (g.get())
1155  g->init_context(context);
1156 
1157  // this->init_context(context);
1158 
1159  // Iterate over all the elements in the range
1160  for (ConstElemRange::const_iterator elem_it=range.begin(); elem_it != range.end();
1161  ++elem_it)
1162  {
1163  const Elem * elem = *elem_it;
1164 
1165  unsigned int dim = elem->dim();
1166 
1167  context.pre_fe_reinit(system, elem);
1168 
1169  // If we're doing AMR, this might be a grid projection with a cheap
1170  // early exit.
1171 #ifdef LIBMESH_ENABLE_AMR
1172  // If this element doesn't have an old_dof_object, but it
1173  // wasn't just refined or just coarsened into activity, then
1174  // it must be newly added, so the user is responsible for
1175  // setting the new dofs on it during a grid projection.
1176  if (!elem->old_dof_object &&
1177  elem->refinement_flag() != Elem::JUST_REFINED &&
1178  elem->refinement_flag() != Elem::JUST_COARSENED &&
1179  f.is_grid_projection())
1180  continue;
1181 #endif // LIBMESH_ENABLE_AMR
1182 
1183  // Loop over all the variables we've been requested to project, to
1184  // do the projection
1185  for (std::size_t v=0; v!=variables.size(); v++)
1186  {
1187  const unsigned int var = variables[v];
1188 
1189  const Variable & variable = dof_map.variable(var);
1190 
1191  const FEType & base_fe_type = variable.type();
1192 
1193  FEType fe_type = base_fe_type;
1194 
1195  // This may be a p refined element
1196  fe_type.order =
1197  libMesh::Order (fe_type.order + elem->p_level());
1198 
1199  if (fe_type.family == SCALAR)
1200  continue;
1201 
1202  // Per-subdomain variables don't need to be projected on
1203  // elements where they're not active
1204  if (!variable.active_on_subdomain(elem->subdomain_id()))
1205  continue;
1206 
1207  const std::vector<dof_id_type> & dof_indices =
1208  context.get_dof_indices(var);
1209 
1210  // The number of DOFs on the element
1211  const unsigned int n_dofs =
1212  cast_int<unsigned int>(dof_indices.size());
1213 
1214  const unsigned int var_component =
1216 
1217  // Zero the interpolated values
1218  Ue.resize (n_dofs); Ue.zero();
1219 
1220  // If we're projecting from an old grid
1221 #ifdef LIBMESH_ENABLE_AMR
1222  if (f.is_grid_projection() &&
1223  // And either this is an unchanged element
1224  ((elem->refinement_flag() != Elem::JUST_REFINED &&
1225  elem->refinement_flag() != Elem::JUST_COARSENED &&
1226  elem->p_refinement_flag() != Elem::JUST_REFINED &&
1227  elem->p_refinement_flag() != Elem::JUST_COARSENED) ||
1228  // Or this is a low order monomial element which has merely
1229  // been h refined
1230  (fe_type.family == MONOMIAL &&
1231  fe_type.order == CONSTANT &&
1232  elem->p_level() == 0 &&
1233  elem->refinement_flag() != Elem::JUST_COARSENED &&
1234  elem->p_refinement_flag() != Elem::JUST_COARSENED))
1235  )
1236  // then we can simply copy its old dof
1237  // values to new indices.
1238  {
1239  LOG_SCOPE ("copy_dofs", "GenericProjector");
1240 
1241  f.eval_old_dofs(context, var_component, Ue.get_values());
1242 
1243  action.insert(context, var, Ue);
1244 
1245  continue;
1246  }
1247 #endif // LIBMESH_ENABLE_AMR
1248 
1249  FEBase * fe = libmesh_nullptr;
1250  FEBase * side_fe = libmesh_nullptr;
1251  FEBase * edge_fe = libmesh_nullptr;
1252 
1253  context.get_element_fe( var, fe, dim );
1254  if (fe->get_fe_type().family == SCALAR)
1255  continue;
1256  if (dim > 1)
1257  context.get_side_fe( var, side_fe, dim );
1258  if (dim > 2)
1259  context.get_edge_fe( var, edge_fe );
1260 
1261  const FEContinuity cont = fe->get_continuity();
1262 
1263  std::vector<unsigned int> side_dofs;
1264 
1265  // Fixed vs. free DoFs on edge/face projections
1266  std::vector<char> dof_is_fixed(n_dofs, false); // bools
1267  std::vector<int> free_dof(n_dofs, 0);
1268 
1269  // The element type
1270  const ElemType elem_type = elem->type();
1271 
1272  // The number of nodes on the new element
1273  const unsigned int n_nodes = elem->n_nodes();
1274 
1275  START_LOG ("project_nodes","GenericProjector");
1276 
1277  // When interpolating C1 elements we need to know which
1278  // vertices were also parent vertices; we'll cache an
1279  // intermediate fact outside the nodes loop.
1280  int i_am_child = -1;
1281 #ifdef LIBMESH_ENABLE_AMR
1282  if (elem->parent())
1283  i_am_child = elem->parent()->which_child_am_i(elem);
1284 #endif
1285  // In general, we need a series of
1286  // projections to ensure a unique and continuous
1287  // solution. We start by interpolating nodes, then
1288  // hold those fixed and project edges, then
1289  // hold those fixed and project faces, then
1290  // hold those fixed and project interiors
1291  //
1292  // In the LAGRANGE case, we will save a lot of solution
1293  // evaluations (at a slight cost in accuracy) by simply
1294  // interpolating all nodes rather than projecting.
1295 
1296  // Interpolate vertex (or for LAGRANGE, all node) values first.
1297  unsigned int current_dof = 0;
1298  for (unsigned int n=0; n!= n_nodes; ++n)
1299  {
1300  // FIXME: this should go through the DofMap,
1301  // not duplicate dof_indices code badly!
1302  const unsigned int nc =
1303  FEInterface::n_dofs_at_node (dim, fe_type, elem_type, n);
1304 
1305  if (!elem->is_vertex(n) &&
1306  fe_type.family != LAGRANGE)
1307  {
1308  current_dof += nc;
1309  continue;
1310  }
1311 
1312  if (cont == DISCONTINUOUS)
1313  {
1314  libmesh_assert_equal_to (nc, 0);
1315  }
1316  else if (!nc)
1317  {
1318  // This should only occur for first-order LAGRANGE
1319  // FE on non-vertices of higher-order elements
1320  libmesh_assert (!elem->is_vertex(n));
1321  libmesh_assert_equal_to(fe_type.family, LAGRANGE);
1322  }
1323  // Assume that C_ZERO elements have a single nodal
1324  // value shape function at vertices
1325  else if (cont == C_ZERO)
1326  {
1327  Ue(current_dof) = f.eval_at_node(context,
1328  var_component,
1329  dim,
1330  elem->node_ref(n),
1331  system.time);
1332  dof_is_fixed[current_dof] = true;
1333  current_dof++;
1334  }
1335  else if (cont == C_ONE)
1336  {
1337  const bool is_parent_vertex = (i_am_child == -1) ||
1338  elem->parent()->is_vertex_on_parent(i_am_child, n);
1339 
1340  // The hermite element vertex shape functions are weird
1341  if (fe_type.family == HERMITE)
1342  {
1343  Ue(current_dof) =
1344  f.eval_at_node(context,
1345  var_component,
1346  dim,
1347  elem->node_ref(n),
1348  system.time);
1349  dof_is_fixed[current_dof] = true;
1350  current_dof++;
1351  VectorValue<FValue> grad =
1352  is_parent_vertex ?
1353  g->eval_at_node(context,
1354  var_component,
1355  dim,
1356  elem->node_ref(n),
1357  system.time) :
1358  g->eval_at_point(context,
1359  var_component,
1360  elem->node_ref(n),
1361  system.time);
1362  // x derivative
1363  Ue(current_dof) = grad(0);
1364  dof_is_fixed[current_dof] = true;
1365  current_dof++;
1366  if (dim > 1)
1367  {
1368  // We'll finite difference mixed derivatives
1369  Real delta_x = TOLERANCE * elem->hmin();
1370 
1371  Point nxminus = elem->point(n),
1372  nxplus = elem->point(n);
1373  nxminus(0) -= delta_x;
1374  nxplus(0) += delta_x;
1375  VectorValue<FValue> gxminus =
1376  g->eval_at_point(context,
1377  var_component,
1378  nxminus,
1379  system.time);
1380  VectorValue<FValue> gxplus =
1381  g->eval_at_point(context,
1382  var_component,
1383  nxplus,
1384  system.time);
1385  // y derivative
1386  Ue(current_dof) = grad(1);
1387  dof_is_fixed[current_dof] = true;
1388  current_dof++;
1389  // xy derivative
1390  Ue(current_dof) = (gxplus(1) - gxminus(1))
1391  / 2. / delta_x;
1392  dof_is_fixed[current_dof] = true;
1393  current_dof++;
1394 
1395  if (dim > 2)
1396  {
1397  // z derivative
1398  Ue(current_dof) = grad(2);
1399  dof_is_fixed[current_dof] = true;
1400  current_dof++;
1401  // xz derivative
1402  Ue(current_dof) = (gxplus(2) - gxminus(2))
1403  / 2. / delta_x;
1404  dof_is_fixed[current_dof] = true;
1405  current_dof++;
1406  // We need new points for yz
1407  Point nyminus = elem->point(n),
1408  nyplus = elem->point(n);
1409  nyminus(1) -= delta_x;
1410  nyplus(1) += delta_x;
1411  VectorValue<FValue> gyminus =
1412  g->eval_at_point(context,
1413  var_component,
1414  nyminus,
1415  system.time);
1416  VectorValue<FValue> gyplus =
1417  g->eval_at_point(context,
1418  var_component,
1419  nyplus,
1420  system.time);
1421  // xz derivative
1422  Ue(current_dof) = (gyplus(2) - gyminus(2))
1423  / 2. / delta_x;
1424  dof_is_fixed[current_dof] = true;
1425  current_dof++;
1426  // Getting a 2nd order xyz is more tedious
1427  Point nxmym = elem->point(n),
1428  nxmyp = elem->point(n),
1429  nxpym = elem->point(n),
1430  nxpyp = elem->point(n);
1431  nxmym(0) -= delta_x;
1432  nxmym(1) -= delta_x;
1433  nxmyp(0) -= delta_x;
1434  nxmyp(1) += delta_x;
1435  nxpym(0) += delta_x;
1436  nxpym(1) -= delta_x;
1437  nxpyp(0) += delta_x;
1438  nxpyp(1) += delta_x;
1439  VectorValue<FValue> gxmym =
1440  g->eval_at_point(context,
1441  var_component,
1442  nxmym,
1443  system.time);
1444  VectorValue<FValue> gxmyp =
1445  g->eval_at_point(context,
1446  var_component,
1447  nxmyp,
1448  system.time);
1449  VectorValue<FValue> gxpym =
1450  g->eval_at_point(context,
1451  var_component,
1452  nxpym,
1453  system.time);
1454  VectorValue<FValue> gxpyp =
1455  g->eval_at_point(context,
1456  var_component,
1457  nxpyp,
1458  system.time);
1459  FValue gxzplus = (gxpyp(2) - gxmyp(2))
1460  / 2. / delta_x;
1461  FValue gxzminus = (gxpym(2) - gxmym(2))
1462  / 2. / delta_x;
1463  // xyz derivative
1464  Ue(current_dof) = (gxzplus - gxzminus)
1465  / 2. / delta_x;
1466  dof_is_fixed[current_dof] = true;
1467  current_dof++;
1468  }
1469  }
1470  }
1471  else
1472  {
1473  // Assume that other C_ONE elements have a single nodal
1474  // value shape function and nodal gradient component
1475  // shape functions
1476  libmesh_assert_equal_to (nc, 1 + dim);
1477  Ue(current_dof) = f.eval_at_node(context,
1478  var_component,
1479  dim,
1480  elem->node_ref(n),
1481  system.time);
1482  dof_is_fixed[current_dof] = true;
1483  current_dof++;
1484  VectorValue<FValue> grad =
1485  is_parent_vertex ?
1486  g->eval_at_node(context,
1487  var_component,
1488  dim,
1489  elem->node_ref(n),
1490  system.time) :
1491  g->eval_at_point(context,
1492  var_component,
1493  elem->node_ref(n),
1494  system.time);
1495  for (unsigned int i=0; i!= dim; ++i)
1496  {
1497  Ue(current_dof) = grad(i);
1498  dof_is_fixed[current_dof] = true;
1499  current_dof++;
1500  }
1501  }
1502  }
1503  else
1504  libmesh_error_msg("Unknown continuity " << cont);
1505  }
1506 
1507  STOP_LOG ("project_nodes","GenericProjector");
1508 
1509  START_LOG ("project_edges","GenericProjector");
1510 
1511  // In 3D with non-LAGRANGE, project any edge values next
1512  if (dim > 2 &&
1513  cont != DISCONTINUOUS &&
1514  (fe_type.family != LAGRANGE
1515 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1516  || elem->infinite()
1517 #endif
1518  ))
1519  {
1520  // If we're JUST_COARSENED we'll need a custom
1521  // evaluation, not just the standard edge FE
1522  const std::vector<Point> & xyz_values =
1523 #ifdef LIBMESH_ENABLE_AMR
1524  (elem->refinement_flag() == Elem::JUST_COARSENED) ?
1525  fe->get_xyz() :
1526 #endif
1527  edge_fe->get_xyz();
1528  const std::vector<Real> & JxW =
1529 #ifdef LIBMESH_ENABLE_AMR
1530  (elem->refinement_flag() == Elem::JUST_COARSENED) ?
1531  fe->get_JxW() :
1532 #endif
1533  edge_fe->get_JxW();
1534 
1535  const std::vector<std::vector<Real>> & phi =
1536 #ifdef LIBMESH_ENABLE_AMR
1537  (elem->refinement_flag() == Elem::JUST_COARSENED) ?
1538  fe->get_phi() :
1539 #endif
1540  edge_fe->get_phi();
1541  const std::vector<std::vector<RealGradient>> * dphi = libmesh_nullptr;
1542  if (cont == C_ONE)
1543  dphi =
1544 #ifdef LIBMESH_ENABLE_AMR
1545  (elem->refinement_flag() == Elem::JUST_COARSENED) ?
1546  &(fe->get_dphi()) :
1547 #endif
1548  &(edge_fe->get_dphi());
1549 
1550  for (auto e : elem->edge_index_range())
1551  {
1552  context.edge = e;
1553 
1554 #ifdef LIBMESH_ENABLE_AMR
1555  if (elem->refinement_flag() == Elem::JUST_COARSENED)
1556  {
1557  std::vector<Point> fine_points;
1558 
1559  UniquePtr<FEBase> fine_fe
1560  (FEBase::build (dim, base_fe_type));
1561 
1562  UniquePtr<QBase> qrule
1563  (base_fe_type.default_quadrature_rule(1));
1564  fine_fe->attach_quadrature_rule(qrule.get());
1565 
1566  const std::vector<Point> & child_xyz =
1567  fine_fe->get_xyz();
1568 
1569  for (unsigned int c = 0, nc = elem->n_children();
1570  c != nc; ++c)
1571  {
1572  if (!elem->is_child_on_edge(c, e))
1573  continue;
1574 
1575  fine_fe->edge_reinit(elem->child_ptr(c), e);
1576  fine_points.insert(fine_points.end(),
1577  child_xyz.begin(),
1578  child_xyz.end());
1579  }
1580 
1581  std::vector<Point> fine_qp;
1582  FEInterface::inverse_map (dim, base_fe_type, elem,
1583  fine_points, fine_qp);
1584 
1585  context.elem_fe_reinit(&fine_qp);
1586  }
1587  else
1588 #endif // LIBMESH_ENABLE_AMR
1589  context.edge_fe_reinit();
1590 
1591  const unsigned int n_qp = xyz_values.size();
1592 
1593  FEInterface::dofs_on_edge(elem, dim, base_fe_type,
1594  e, side_dofs);
1595 
1596  // Some edge dofs are on nodes and already
1597  // fixed, others are free to calculate
1598  unsigned int free_dofs = 0;
1599  for (std::size_t i=0; i != side_dofs.size(); ++i)
1600  if (!dof_is_fixed[side_dofs[i]])
1601  free_dof[free_dofs++] = i;
1602 
1603  // There may be nothing to project
1604  if (!free_dofs)
1605  continue;
1606 
1607  Ke.resize (free_dofs, free_dofs); Ke.zero();
1608  Fe.resize (free_dofs); Fe.zero();
1609  // The new edge coefficients
1610  DenseVector<FValue> Uedge(free_dofs);
1611 
1612  // Loop over the quadrature points
1613  for (unsigned int qp=0; qp<n_qp; qp++)
1614  {
1615  // solution at the quadrature point
1616  FValue fineval = f.eval_at_point(context,
1617  var_component,
1618  xyz_values[qp],
1619  system.time);
1620  // solution grad at the quadrature point
1621  VectorValue<FValue> finegrad;
1622  if (cont == C_ONE)
1623  finegrad = g->eval_at_point(context,
1624  var_component,
1625  xyz_values[qp],
1626  system.time);
1627 
1628  // Form edge projection matrix
1629  for (std::size_t sidei=0, freei=0;
1630  sidei != side_dofs.size(); ++sidei)
1631  {
1632  unsigned int i = side_dofs[sidei];
1633  // fixed DoFs aren't test functions
1634  if (dof_is_fixed[i])
1635  continue;
1636  for (std::size_t sidej=0, freej=0;
1637  sidej != side_dofs.size(); ++sidej)
1638  {
1639  unsigned int j = side_dofs[sidej];
1640  if (dof_is_fixed[j])
1641  Fe(freei) -= phi[i][qp] * phi[j][qp] *
1642  JxW[qp] * Ue(j);
1643  else
1644  Ke(freei,freej) += phi[i][qp] *
1645  phi[j][qp] * JxW[qp];
1646  if (cont == C_ONE)
1647  {
1648  if (dof_is_fixed[j])
1649  Fe(freei) -= ( (*dphi)[i][qp] *
1650  (*dphi)[j][qp] ) *
1651  JxW[qp] * Ue(j);
1652  else
1653  Ke(freei,freej) += ( (*dphi)[i][qp] *
1654  (*dphi)[j][qp] )
1655  * JxW[qp];
1656  }
1657  if (!dof_is_fixed[j])
1658  freej++;
1659  }
1660  Fe(freei) += phi[i][qp] * fineval * JxW[qp];
1661  if (cont == C_ONE)
1662  Fe(freei) += (finegrad * (*dphi)[i][qp] ) *
1663  JxW[qp];
1664  freei++;
1665  }
1666  }
1667 
1668  Ke.cholesky_solve(Fe, Uedge);
1669 
1670  // Transfer new edge solutions to element
1671  for (unsigned int i=0; i != free_dofs; ++i)
1672  {
1673  FValue & ui = Ue(side_dofs[free_dof[i]]);
1675  std::abs(ui - Uedge(i)) < TOLERANCE);
1676  ui = Uedge(i);
1677  dof_is_fixed[side_dofs[free_dof[i]]] = true;
1678  }
1679  }
1680  } // end if (dim > 2, !DISCONTINUOUS, !LAGRANGE)
1681 
1682  STOP_LOG ("project_edges","GenericProjector");
1683 
1684  START_LOG ("project_sides","GenericProjector");
1685 
1686  // With non-LAGRANGE, project any side values (edges in 2D,
1687  // faces in 3D) next.
1688  if (dim > 1 &&
1689  cont != DISCONTINUOUS &&
1690  (fe_type.family != LAGRANGE
1691 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1692  || elem->infinite()
1693 #endif
1694  ))
1695  {
1696  // If we're JUST_COARSENED we'll need a custom
1697  // evaluation, not just the standard side FE
1698  const std::vector<Point> & xyz_values =
1699 #ifdef LIBMESH_ENABLE_AMR
1700  (elem->refinement_flag() == Elem::JUST_COARSENED) ?
1701  fe->get_xyz() :
1702 #endif // LIBMESH_ENABLE_AMR
1703  side_fe->get_xyz();
1704  const std::vector<Real> & JxW =
1705 #ifdef LIBMESH_ENABLE_AMR
1706  (elem->refinement_flag() == Elem::JUST_COARSENED) ?
1707  fe->get_JxW() :
1708 #endif // LIBMESH_ENABLE_AMR
1709  side_fe->get_JxW();
1710  const std::vector<std::vector<Real>> & phi =
1711 #ifdef LIBMESH_ENABLE_AMR
1712  (elem->refinement_flag() == Elem::JUST_COARSENED) ?
1713  fe->get_phi() :
1714 #endif // LIBMESH_ENABLE_AMR
1715  side_fe->get_phi();
1716  const std::vector<std::vector<RealGradient>> * dphi = libmesh_nullptr;
1717  if (cont == C_ONE)
1718  dphi =
1719 #ifdef LIBMESH_ENABLE_AMR
1720  (elem->refinement_flag() == Elem::JUST_COARSENED) ? &(fe->get_dphi()) :
1721 #endif // LIBMESH_ENABLE_AMR
1722  &(side_fe->get_dphi());
1723 
1724  for (auto s : elem->side_index_range())
1725  {
1726  FEInterface::dofs_on_side(elem, dim, base_fe_type,
1727  s, side_dofs);
1728 
1729  // Some side dofs are on nodes/edges and already
1730  // fixed, others are free to calculate
1731  unsigned int free_dofs = 0;
1732  for (std::size_t i=0; i != side_dofs.size(); ++i)
1733  if (!dof_is_fixed[side_dofs[i]])
1734  free_dof[free_dofs++] = i;
1735 
1736  // There may be nothing to project
1737  if (!free_dofs)
1738  continue;
1739 
1740  Ke.resize (free_dofs, free_dofs); Ke.zero();
1741  Fe.resize (free_dofs); Fe.zero();
1742  // The new side coefficients
1743  DenseVector<FValue> Uside(free_dofs);
1744 
1745  context.side = s;
1746 
1747 #ifdef LIBMESH_ENABLE_AMR
1748  if (elem->refinement_flag() == Elem::JUST_COARSENED)
1749  {
1750  std::vector<Point> fine_points;
1751 
1752  UniquePtr<FEBase> fine_fe
1753  (FEBase::build (dim, base_fe_type));
1754 
1755  UniquePtr<QBase> qrule
1756  (base_fe_type.default_quadrature_rule(dim-1));
1757  fine_fe->attach_quadrature_rule(qrule.get());
1758 
1759  const std::vector<Point> & child_xyz =
1760  fine_fe->get_xyz();
1761 
1762  for (unsigned int c = 0, nc = elem->n_children();
1763  c != nc; ++c)
1764  {
1765  if (!elem->is_child_on_side(c, s))
1766  continue;
1767 
1768  fine_fe->reinit(elem->child_ptr(c), s);
1769  fine_points.insert(fine_points.end(),
1770  child_xyz.begin(),
1771  child_xyz.end());
1772  }
1773 
1774  std::vector<Point> fine_qp;
1775  FEInterface::inverse_map (dim, base_fe_type, elem,
1776  fine_points, fine_qp);
1777 
1778  context.elem_fe_reinit(&fine_qp);
1779  }
1780  else
1781 #endif // LIBMESH_ENABLE_AMR
1782  context.side_fe_reinit();
1783 
1784  const unsigned int n_qp = xyz_values.size();
1785 
1786  // Loop over the quadrature points
1787  for (unsigned int qp=0; qp<n_qp; qp++)
1788  {
1789  // solution at the quadrature point
1790  FValue fineval = f.eval_at_point(context,
1791  var_component,
1792  xyz_values[qp],
1793  system.time);
1794  // solution grad at the quadrature point
1795  VectorValue<FValue> finegrad;
1796  if (cont == C_ONE)
1797  finegrad = g->eval_at_point(context,
1798  var_component,
1799  xyz_values[qp],
1800  system.time);
1801 
1802  // Form side projection matrix
1803  for (std::size_t sidei=0, freei=0;
1804  sidei != side_dofs.size(); ++sidei)
1805  {
1806  unsigned int i = side_dofs[sidei];
1807  // fixed DoFs aren't test functions
1808  if (dof_is_fixed[i])
1809  continue;
1810  for (std::size_t sidej=0, freej=0;
1811  sidej != side_dofs.size(); ++sidej)
1812  {
1813  unsigned int j = side_dofs[sidej];
1814  if (dof_is_fixed[j])
1815  Fe(freei) -= phi[i][qp] * phi[j][qp] *
1816  JxW[qp] * Ue(j);
1817  else
1818  Ke(freei,freej) += phi[i][qp] *
1819  phi[j][qp] * JxW[qp];
1820  if (cont == C_ONE)
1821  {
1822  if (dof_is_fixed[j])
1823  Fe(freei) -= ( (*dphi)[i][qp] *
1824  (*dphi)[j][qp] ) *
1825  JxW[qp] * Ue(j);
1826  else
1827  Ke(freei,freej) += ( (*dphi)[i][qp] *
1828  (*dphi)[j][qp] )
1829  * JxW[qp];
1830  }
1831  if (!dof_is_fixed[j])
1832  freej++;
1833  }
1834  Fe(freei) += (fineval * phi[i][qp]) * JxW[qp];
1835  if (cont == C_ONE)
1836  Fe(freei) += (finegrad * (*dphi)[i][qp] ) *
1837  JxW[qp];
1838  freei++;
1839  }
1840  }
1841 
1842  Ke.cholesky_solve(Fe, Uside);
1843 
1844  // Transfer new side solutions to element
1845  for (unsigned int i=0; i != free_dofs; ++i)
1846  {
1847  FValue & ui = Ue(side_dofs[free_dof[i]]);
1849  std::abs(ui - Uside(i)) < TOLERANCE);
1850  ui = Uside(i);
1851  dof_is_fixed[side_dofs[free_dof[i]]] = true;
1852  }
1853  }
1854  } // end if (dim > 1, !DISCONTINUOUS, !LAGRANGE)
1855 
1856  STOP_LOG ("project_sides","GenericProjector");
1857 
1858  START_LOG ("project_interior","GenericProjector");
1859 
1860  // Project the interior values, finally
1861 
1862  // Some interior dofs are on nodes/edges/sides and
1863  // already fixed, others are free to calculate
1864  unsigned int free_dofs = 0;
1865  for (unsigned int i=0; i != n_dofs; ++i)
1866  if (!dof_is_fixed[i])
1867  free_dof[free_dofs++] = i;
1868 
1869  // Project any remaining (interior) dofs in the non-LAGRANGE
1870  // case.
1871  if (free_dofs &&
1872  (fe_type.family != LAGRANGE
1873 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1874  || elem->infinite()
1875 #endif
1876  ))
1877  {
1878  const std::vector<Point> & xyz_values = fe->get_xyz();
1879  const std::vector<Real> & JxW = fe->get_JxW();
1880 
1881  const std::vector<std::vector<Real>> & phi = fe->get_phi();
1882  const std::vector<std::vector<RealGradient>> * dphi = libmesh_nullptr;
1883  if (cont == C_ONE)
1884  dphi = &(fe->get_dphi());
1885 
1886 #ifdef LIBMESH_ENABLE_AMR
1887  if (elem->refinement_flag() == Elem::JUST_COARSENED)
1888  {
1889  std::vector<Point> fine_points;
1890 
1891  UniquePtr<FEBase> fine_fe
1892  (FEBase::build (dim, base_fe_type));
1893 
1894  UniquePtr<QBase> qrule
1895  (base_fe_type.default_quadrature_rule(dim));
1896  fine_fe->attach_quadrature_rule(qrule.get());
1897 
1898  const std::vector<Point> & child_xyz =
1899  fine_fe->get_xyz();
1900 
1901  for (auto & child : elem->child_ref_range())
1902  {
1903  fine_fe->reinit(&child);
1904  fine_points.insert(fine_points.end(),
1905  child_xyz.begin(),
1906  child_xyz.end());
1907  }
1908 
1909  std::vector<Point> fine_qp;
1910  FEInterface::inverse_map (dim, base_fe_type, elem,
1911  fine_points, fine_qp);
1912 
1913  context.elem_fe_reinit(&fine_qp);
1914  }
1915  else
1916 #endif // LIBMESH_ENABLE_AMR
1917  context.elem_fe_reinit();
1918 
1919  const unsigned int n_qp = xyz_values.size();
1920 
1921  Ke.resize (free_dofs, free_dofs); Ke.zero();
1922  Fe.resize (free_dofs); Fe.zero();
1923  // The new interior coefficients
1924  DenseVector<FValue> Uint(free_dofs);
1925 
1926  // Loop over the quadrature points
1927  for (unsigned int qp=0; qp<n_qp; qp++)
1928  {
1929  // solution at the quadrature point
1930  FValue fineval = f.eval_at_point(context,
1931  var_component,
1932  xyz_values[qp],
1933  system.time);
1934  // solution grad at the quadrature point
1935  VectorValue<FValue> finegrad;
1936  if (cont == C_ONE)
1937  finegrad = g->eval_at_point(context,
1938  var_component,
1939  xyz_values[qp],
1940  system.time);
1941 
1942  // Form interior projection matrix
1943  for (unsigned int i=0, freei=0; i != n_dofs; ++i)
1944  {
1945  // fixed DoFs aren't test functions
1946  if (dof_is_fixed[i])
1947  continue;
1948  for (unsigned int j=0, freej=0; j != n_dofs; ++j)
1949  {
1950  if (dof_is_fixed[j])
1951  Fe(freei) -= phi[i][qp] * phi[j][qp] * JxW[qp]
1952  * Ue(j);
1953  else
1954  Ke(freei,freej) += phi[i][qp] * phi[j][qp] *
1955  JxW[qp];
1956  if (cont == C_ONE)
1957  {
1958  if (dof_is_fixed[j])
1959  Fe(freei) -= ( (*dphi)[i][qp] *
1960  (*dphi)[j][qp] ) * JxW[qp] *
1961  Ue(j);
1962  else
1963  Ke(freei,freej) += ( (*dphi)[i][qp] *
1964  (*dphi)[j][qp] ) *
1965  JxW[qp];
1966  }
1967  if (!dof_is_fixed[j])
1968  freej++;
1969  }
1970  Fe(freei) += phi[i][qp] * fineval * JxW[qp];
1971  if (cont == C_ONE)
1972  Fe(freei) += (finegrad * (*dphi)[i][qp] ) * JxW[qp];
1973  freei++;
1974  }
1975  }
1976  Ke.cholesky_solve(Fe, Uint);
1977 
1978  // Transfer new interior solutions to element
1979  for (unsigned int i=0; i != free_dofs; ++i)
1980  {
1981  FValue & ui = Ue(free_dof[i]);
1983  std::abs(ui - Uint(i)) < TOLERANCE);
1984  ui = Uint(i);
1985  dof_is_fixed[free_dof[i]] = true;
1986  }
1987 
1988  } // if there are free interior dofs
1989 
1990  STOP_LOG ("project_interior","GenericProjector");
1991 
1992  // Make sure every DoF got reached!
1993  for (unsigned int i=0; i != n_dofs; ++i)
1994  libmesh_assert(dof_is_fixed[i]);
1995 
1996  action.insert(context, var, Ue);
1997  } // end variables loop
1998  } // end elements loop
1999 }
static void dofs_on_edge(const Elem *const elem, const unsigned int dim, const FEType &fe_t, unsigned int e, std::vector< unsigned int > &di)
Fills the vector di with the local degree of freedom indices associated with edge e of element elem A...
Definition: fe_interface.C:510
Real time
For time-dependent problems, this is the time t at the beginning of the current timestep.
Definition: system.h:1545
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)
Fills the vector di with the local degree of freedom indices associated with side s of element elem A...
Definition: fe_interface.C:495
unsigned int dim
ElemType
Defines an enum for geometric element types.
const class libmesh_nullptr_t libmesh_nullptr
static const Real TOLERANCE
libmesh_assert(j)
const dof_id_type n_nodes
Definition: tecplot_io.C:67
const DofMap & get_dof_map() const
Definition: system.h:2030
static Point inverse_map(const unsigned int dim, const FEType &fe_t, const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true)
Definition: fe_interface.C:569
const ProjectionAction & master_action
FEGenericBase< Real > FEBase
const std::vector< unsigned int > & variables
FEContinuity
defines an enum for finite element types to libmesh_assert a certain level (or type? Hcurl?) of continuity.
std::vector< object_type >::const_iterator const_iterator
Allows an StoredRange to behave like an STL container.
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:436
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Order
defines an enum for polynomial orders.
Definition: enum_order.h:32
unsigned int variable_scalar_number(const std::string &var, unsigned int component) const
Definition: system.h:2145
static UniquePtr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.

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 system_projection.C.

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: