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

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

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

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

Member Function Documentation

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

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::DenseMatrix< T >::cholesky_solve(), libMesh::CONSTANT, libMesh::FEType::default_quadrature_rule(), libMesh::Elem::dim(), libMesh::DISCONTINUOUS, libMesh::FEInterface::dofs_on_edge(), libMesh::FEInterface::dofs_on_side(), libMesh::FEMContext::edge, libMesh::FEMContext::edge_fe_reinit(), libMesh::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(), libMesh::Elem::n_edges(), n_nodes, libMesh::Elem::n_nodes(), libMesh::Elem::n_sides(), 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::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
1339  (i_am_child, n);
1340 
1341  // The hermite element vertex shape functions are weird
1342  if (fe_type.family == HERMITE)
1343  {
1344  Ue(current_dof) =
1345  f.eval_at_node(context,
1346  var_component,
1347  dim,
1348  elem->node_ref(n),
1349  system.time);
1350  dof_is_fixed[current_dof] = true;
1351  current_dof++;
1352  VectorValue<FValue> grad =
1353  is_parent_vertex ?
1354  g->eval_at_node(context,
1355  var_component,
1356  dim,
1357  elem->node_ref(n),
1358  system.time) :
1359  g->eval_at_point(context,
1360  var_component,
1361  elem->node_ref(n),
1362  system.time);
1363  // x derivative
1364  Ue(current_dof) = grad(0);
1365  dof_is_fixed[current_dof] = true;
1366  current_dof++;
1367  if (dim > 1)
1368  {
1369  // We'll finite difference mixed derivatives
1370  Real delta_x = TOLERANCE * elem->hmin();
1371 
1372  Point nxminus = elem->point(n),
1373  nxplus = elem->point(n);
1374  nxminus(0) -= delta_x;
1375  nxplus(0) += delta_x;
1376  VectorValue<FValue> gxminus =
1377  g->eval_at_point(context,
1378  var_component,
1379  nxminus,
1380  system.time);
1381  VectorValue<FValue> gxplus =
1382  g->eval_at_point(context,
1383  var_component,
1384  nxplus,
1385  system.time);
1386  // y derivative
1387  Ue(current_dof) = grad(1);
1388  dof_is_fixed[current_dof] = true;
1389  current_dof++;
1390  // xy derivative
1391  Ue(current_dof) = (gxplus(1) - gxminus(1))
1392  / 2. / delta_x;
1393  dof_is_fixed[current_dof] = true;
1394  current_dof++;
1395 
1396  if (dim > 2)
1397  {
1398  // z derivative
1399  Ue(current_dof) = grad(2);
1400  dof_is_fixed[current_dof] = true;
1401  current_dof++;
1402  // xz derivative
1403  Ue(current_dof) = (gxplus(2) - gxminus(2))
1404  / 2. / delta_x;
1405  dof_is_fixed[current_dof] = true;
1406  current_dof++;
1407  // We need new points for yz
1408  Point nyminus = elem->point(n),
1409  nyplus = elem->point(n);
1410  nyminus(1) -= delta_x;
1411  nyplus(1) += delta_x;
1412  VectorValue<FValue> gyminus =
1413  g->eval_at_point(context,
1414  var_component,
1415  nyminus,
1416  system.time);
1417  VectorValue<FValue> gyplus =
1418  g->eval_at_point(context,
1419  var_component,
1420  nyplus,
1421  system.time);
1422  // xz derivative
1423  Ue(current_dof) = (gyplus(2) - gyminus(2))
1424  / 2. / delta_x;
1425  dof_is_fixed[current_dof] = true;
1426  current_dof++;
1427  // Getting a 2nd order xyz is more tedious
1428  Point nxmym = elem->point(n),
1429  nxmyp = elem->point(n),
1430  nxpym = elem->point(n),
1431  nxpyp = elem->point(n);
1432  nxmym(0) -= delta_x;
1433  nxmym(1) -= delta_x;
1434  nxmyp(0) -= delta_x;
1435  nxmyp(1) += delta_x;
1436  nxpym(0) += delta_x;
1437  nxpym(1) -= delta_x;
1438  nxpyp(0) += delta_x;
1439  nxpyp(1) += delta_x;
1440  VectorValue<FValue> gxmym =
1441  g->eval_at_point(context,
1442  var_component,
1443  nxmym,
1444  system.time);
1445  VectorValue<FValue> gxmyp =
1446  g->eval_at_point(context,
1447  var_component,
1448  nxmyp,
1449  system.time);
1450  VectorValue<FValue> gxpym =
1451  g->eval_at_point(context,
1452  var_component,
1453  nxpym,
1454  system.time);
1455  VectorValue<FValue> gxpyp =
1456  g->eval_at_point(context,
1457  var_component,
1458  nxpyp,
1459  system.time);
1460  FValue gxzplus = (gxpyp(2) - gxmyp(2))
1461  / 2. / delta_x;
1462  FValue gxzminus = (gxpym(2) - gxmym(2))
1463  / 2. / delta_x;
1464  // xyz derivative
1465  Ue(current_dof) = (gxzplus - gxzminus)
1466  / 2. / delta_x;
1467  dof_is_fixed[current_dof] = true;
1468  current_dof++;
1469  }
1470  }
1471  }
1472  else
1473  {
1474  // Assume that other C_ONE elements have a single nodal
1475  // value shape function and nodal gradient component
1476  // shape functions
1477  libmesh_assert_equal_to (nc, 1 + dim);
1478  Ue(current_dof) = f.eval_at_node(context,
1479  var_component,
1480  dim,
1481  elem->node_ref(n),
1482  system.time);
1483  dof_is_fixed[current_dof] = true;
1484  current_dof++;
1485  VectorValue<FValue> grad =
1486  is_parent_vertex ?
1487  g->eval_at_node(context,
1488  var_component,
1489  dim,
1490  elem->node_ref(n),
1491  system.time) :
1492  g->eval_at_point(context,
1493  var_component,
1494  elem->node_ref(n),
1495  system.time);
1496  for (unsigned int i=0; i!= dim; ++i)
1497  {
1498  Ue(current_dof) = grad(i);
1499  dof_is_fixed[current_dof] = true;
1500  current_dof++;
1501  }
1502  }
1503  }
1504  else
1505  libmesh_error_msg("Unknown continuity " << cont);
1506  }
1507 
1508  STOP_LOG ("project_nodes","GenericProjector");
1509 
1510  START_LOG ("project_edges","GenericProjector");
1511 
1512  // In 3D with non-LAGRANGE, project any edge values next
1513  if (dim > 2 &&
1514  cont != DISCONTINUOUS &&
1515  (fe_type.family != LAGRANGE
1516 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1517  || elem->infinite()
1518 #endif
1519  ))
1520  {
1521  // If we're JUST_COARSENED we'll need a custom
1522  // evaluation, not just the standard edge FE
1523  const std::vector<Point> & xyz_values =
1524 #ifdef LIBMESH_ENABLE_AMR
1525  (elem->refinement_flag() == Elem::JUST_COARSENED) ?
1526  fe->get_xyz() :
1527 #endif
1528  edge_fe->get_xyz();
1529  const std::vector<Real> & JxW =
1530 #ifdef LIBMESH_ENABLE_AMR
1531  (elem->refinement_flag() == Elem::JUST_COARSENED) ?
1532  fe->get_JxW() :
1533 #endif
1534  edge_fe->get_JxW();
1535 
1536  const std::vector<std::vector<Real> > & phi =
1537 #ifdef LIBMESH_ENABLE_AMR
1538  (elem->refinement_flag() == Elem::JUST_COARSENED) ?
1539  fe->get_phi() :
1540 #endif
1541  edge_fe->get_phi();
1542  const std::vector<std::vector<RealGradient> > * dphi = libmesh_nullptr;
1543  if (cont == C_ONE)
1544  dphi =
1545 #ifdef LIBMESH_ENABLE_AMR
1546  (elem->refinement_flag() == Elem::JUST_COARSENED) ?
1547  &(fe->get_dphi()) :
1548 #endif
1549  &(edge_fe->get_dphi());
1550 
1551  for (unsigned char e=0; e != elem->n_edges(); ++e)
1552  {
1553  context.edge = e;
1554 
1555 #ifdef LIBMESH_ENABLE_AMR
1556  if (elem->refinement_flag() == Elem::JUST_COARSENED)
1557  {
1558  std::vector<Point> fine_points;
1559 
1560  UniquePtr<FEBase> fine_fe
1561  (FEBase::build (dim, base_fe_type));
1562 
1563  UniquePtr<QBase> qrule
1564  (base_fe_type.default_quadrature_rule(1));
1565  fine_fe->attach_quadrature_rule(qrule.get());
1566 
1567  const std::vector<Point> & child_xyz =
1568  fine_fe->get_xyz();
1569 
1570  for (unsigned int c = 0;
1571  c != elem->n_children(); ++c)
1572  {
1573  if (!elem->is_child_on_edge(c, e))
1574  continue;
1575 
1576  fine_fe->edge_reinit(elem->child_ptr(c), e);
1577  fine_points.insert(fine_points.end(),
1578  child_xyz.begin(),
1579  child_xyz.end());
1580  }
1581 
1582  std::vector<Point> fine_qp;
1583  FEInterface::inverse_map (dim, base_fe_type, elem,
1584  fine_points, fine_qp);
1585 
1586  context.elem_fe_reinit(&fine_qp);
1587  }
1588  else
1589 #endif // LIBMESH_ENABLE_AMR
1590  context.edge_fe_reinit();
1591 
1592  const unsigned int n_qp = xyz_values.size();
1593 
1594  FEInterface::dofs_on_edge(elem, dim, base_fe_type,
1595  e, side_dofs);
1596 
1597  // Some edge dofs are on nodes and already
1598  // fixed, others are free to calculate
1599  unsigned int free_dofs = 0;
1600  for (std::size_t i=0; i != side_dofs.size(); ++i)
1601  if (!dof_is_fixed[side_dofs[i]])
1602  free_dof[free_dofs++] = i;
1603 
1604  // There may be nothing to project
1605  if (!free_dofs)
1606  continue;
1607 
1608  Ke.resize (free_dofs, free_dofs); Ke.zero();
1609  Fe.resize (free_dofs); Fe.zero();
1610  // The new edge coefficients
1611  DenseVector<FValue> Uedge(free_dofs);
1612 
1613  // Loop over the quadrature points
1614  for (unsigned int qp=0; qp<n_qp; qp++)
1615  {
1616  // solution at the quadrature point
1617  FValue fineval = f.eval_at_point(context,
1618  var_component,
1619  xyz_values[qp],
1620  system.time);
1621  // solution grad at the quadrature point
1622  VectorValue<FValue> finegrad;
1623  if (cont == C_ONE)
1624  finegrad = g->eval_at_point(context,
1625  var_component,
1626  xyz_values[qp],
1627  system.time);
1628 
1629  // Form edge projection matrix
1630  for (std::size_t sidei=0, freei=0;
1631  sidei != side_dofs.size(); ++sidei)
1632  {
1633  unsigned int i = side_dofs[sidei];
1634  // fixed DoFs aren't test functions
1635  if (dof_is_fixed[i])
1636  continue;
1637  for (std::size_t sidej=0, freej=0;
1638  sidej != side_dofs.size(); ++sidej)
1639  {
1640  unsigned int j = side_dofs[sidej];
1641  if (dof_is_fixed[j])
1642  Fe(freei) -= phi[i][qp] * phi[j][qp] *
1643  JxW[qp] * Ue(j);
1644  else
1645  Ke(freei,freej) += phi[i][qp] *
1646  phi[j][qp] * JxW[qp];
1647  if (cont == C_ONE)
1648  {
1649  if (dof_is_fixed[j])
1650  Fe(freei) -= ( (*dphi)[i][qp] *
1651  (*dphi)[j][qp] ) *
1652  JxW[qp] * Ue(j);
1653  else
1654  Ke(freei,freej) += ( (*dphi)[i][qp] *
1655  (*dphi)[j][qp] )
1656  * JxW[qp];
1657  }
1658  if (!dof_is_fixed[j])
1659  freej++;
1660  }
1661  Fe(freei) += phi[i][qp] * fineval * JxW[qp];
1662  if (cont == C_ONE)
1663  Fe(freei) += (finegrad * (*dphi)[i][qp] ) *
1664  JxW[qp];
1665  freei++;
1666  }
1667  }
1668 
1669  Ke.cholesky_solve(Fe, Uedge);
1670 
1671  // Transfer new edge solutions to element
1672  for (unsigned int i=0; i != free_dofs; ++i)
1673  {
1674  FValue & ui = Ue(side_dofs[free_dof[i]]);
1676  std::abs(ui - Uedge(i)) < TOLERANCE);
1677  ui = Uedge(i);
1678  dof_is_fixed[side_dofs[free_dof[i]]] = true;
1679  }
1680  }
1681  } // end if (dim > 2, !DISCONTINUOUS, !LAGRANGE)
1682 
1683  STOP_LOG ("project_edges","GenericProjector");
1684 
1685  START_LOG ("project_sides","GenericProjector");
1686 
1687  // With non-LAGRANGE, project any side values (edges in 2D,
1688  // faces in 3D) next.
1689  if (dim > 1 &&
1690  cont != DISCONTINUOUS &&
1691  (fe_type.family != LAGRANGE
1692 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1693  || elem->infinite()
1694 #endif
1695  ))
1696  {
1697  // If we're JUST_COARSENED we'll need a custom
1698  // evaluation, not just the standard side FE
1699  const std::vector<Point> & xyz_values =
1700 #ifdef LIBMESH_ENABLE_AMR
1701  (elem->refinement_flag() == Elem::JUST_COARSENED) ?
1702  fe->get_xyz() :
1703 #endif // LIBMESH_ENABLE_AMR
1704  side_fe->get_xyz();
1705  const std::vector<Real> & JxW =
1706 #ifdef LIBMESH_ENABLE_AMR
1707  (elem->refinement_flag() == Elem::JUST_COARSENED) ?
1708  fe->get_JxW() :
1709 #endif // LIBMESH_ENABLE_AMR
1710  side_fe->get_JxW();
1711  const std::vector<std::vector<Real> > & phi =
1712 #ifdef LIBMESH_ENABLE_AMR
1713  (elem->refinement_flag() == Elem::JUST_COARSENED) ?
1714  fe->get_phi() :
1715 #endif // LIBMESH_ENABLE_AMR
1716  side_fe->get_phi();
1717  const std::vector<std::vector<RealGradient> > * dphi = libmesh_nullptr;
1718  if (cont == C_ONE)
1719  dphi =
1720 #ifdef LIBMESH_ENABLE_AMR
1721  (elem->refinement_flag() == Elem::JUST_COARSENED) ?
1722  &(fe->get_dphi()) :
1723 #endif // LIBMESH_ENABLE_AMR
1724  &(side_fe->get_dphi());
1725 
1726  for (unsigned char s=0; s != elem->n_sides(); ++s)
1727  {
1728  FEInterface::dofs_on_side(elem, dim, base_fe_type,
1729  s, side_dofs);
1730 
1731  // Some side dofs are on nodes/edges and already
1732  // fixed, others are free to calculate
1733  unsigned int free_dofs = 0;
1734  for (std::size_t i=0; i != side_dofs.size(); ++i)
1735  if (!dof_is_fixed[side_dofs[i]])
1736  free_dof[free_dofs++] = i;
1737 
1738  // There may be nothing to project
1739  if (!free_dofs)
1740  continue;
1741 
1742  Ke.resize (free_dofs, free_dofs); Ke.zero();
1743  Fe.resize (free_dofs); Fe.zero();
1744  // The new side coefficients
1745  DenseVector<FValue> Uside(free_dofs);
1746 
1747  context.side = s;
1748 
1749 #ifdef LIBMESH_ENABLE_AMR
1750  if (elem->refinement_flag() == Elem::JUST_COARSENED)
1751  {
1752  std::vector<Point> fine_points;
1753 
1754  UniquePtr<FEBase> fine_fe
1755  (FEBase::build (dim, base_fe_type));
1756 
1757  UniquePtr<QBase> qrule
1758  (base_fe_type.default_quadrature_rule(dim-1));
1759  fine_fe->attach_quadrature_rule(qrule.get());
1760 
1761  const std::vector<Point> & child_xyz =
1762  fine_fe->get_xyz();
1763 
1764  for (unsigned int c = 0;
1765  c != elem->n_children(); ++c)
1766  {
1767  if (!elem->is_child_on_side(c, s))
1768  continue;
1769 
1770  fine_fe->reinit(elem->child_ptr(c), s);
1771  fine_points.insert(fine_points.end(),
1772  child_xyz.begin(),
1773  child_xyz.end());
1774  }
1775 
1776  std::vector<Point> fine_qp;
1777  FEInterface::inverse_map (dim, base_fe_type, elem,
1778  fine_points, fine_qp);
1779 
1780  context.elem_fe_reinit(&fine_qp);
1781  }
1782  else
1783 #endif // LIBMESH_ENABLE_AMR
1784  context.side_fe_reinit();
1785 
1786  const unsigned int n_qp = xyz_values.size();
1787 
1788  // Loop over the quadrature points
1789  for (unsigned int qp=0; qp<n_qp; qp++)
1790  {
1791  // solution at the quadrature point
1792  FValue fineval = f.eval_at_point(context,
1793  var_component,
1794  xyz_values[qp],
1795  system.time);
1796  // solution grad at the quadrature point
1797  VectorValue<FValue> finegrad;
1798  if (cont == C_ONE)
1799  finegrad = g->eval_at_point(context,
1800  var_component,
1801  xyz_values[qp],
1802  system.time);
1803 
1804  // Form side projection matrix
1805  for (std::size_t sidei=0, freei=0;
1806  sidei != side_dofs.size(); ++sidei)
1807  {
1808  unsigned int i = side_dofs[sidei];
1809  // fixed DoFs aren't test functions
1810  if (dof_is_fixed[i])
1811  continue;
1812  for (std::size_t sidej=0, freej=0;
1813  sidej != side_dofs.size(); ++sidej)
1814  {
1815  unsigned int j = side_dofs[sidej];
1816  if (dof_is_fixed[j])
1817  Fe(freei) -= phi[i][qp] * phi[j][qp] *
1818  JxW[qp] * Ue(j);
1819  else
1820  Ke(freei,freej) += phi[i][qp] *
1821  phi[j][qp] * JxW[qp];
1822  if (cont == C_ONE)
1823  {
1824  if (dof_is_fixed[j])
1825  Fe(freei) -= ( (*dphi)[i][qp] *
1826  (*dphi)[j][qp] ) *
1827  JxW[qp] * Ue(j);
1828  else
1829  Ke(freei,freej) += ( (*dphi)[i][qp] *
1830  (*dphi)[j][qp] )
1831  * JxW[qp];
1832  }
1833  if (!dof_is_fixed[j])
1834  freej++;
1835  }
1836  Fe(freei) += (fineval * phi[i][qp]) * JxW[qp];
1837  if (cont == C_ONE)
1838  Fe(freei) += (finegrad * (*dphi)[i][qp] ) *
1839  JxW[qp];
1840  freei++;
1841  }
1842  }
1843 
1844  Ke.cholesky_solve(Fe, Uside);
1845 
1846  // Transfer new side solutions to element
1847  for (unsigned int i=0; i != free_dofs; ++i)
1848  {
1849  FValue & ui = Ue(side_dofs[free_dof[i]]);
1851  std::abs(ui - Uside(i)) < TOLERANCE);
1852  ui = Uside(i);
1853  dof_is_fixed[side_dofs[free_dof[i]]] = true;
1854  }
1855  }
1856  } // end if (dim > 1, !DISCONTINUOUS, !LAGRANGE)
1857 
1858  STOP_LOG ("project_sides","GenericProjector");
1859 
1860  START_LOG ("project_interior","GenericProjector");
1861 
1862  // Project the interior values, finally
1863 
1864  // Some interior dofs are on nodes/edges/sides and
1865  // already fixed, others are free to calculate
1866  unsigned int free_dofs = 0;
1867  for (unsigned int i=0; i != n_dofs; ++i)
1868  if (!dof_is_fixed[i])
1869  free_dof[free_dofs++] = i;
1870 
1871  // Project any remaining (interior) dofs in the non-LAGRANGE
1872  // case.
1873  if (free_dofs &&
1874  (fe_type.family != LAGRANGE
1875 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1876  || elem->infinite()
1877 #endif
1878  ))
1879  {
1880  const std::vector<Point> & xyz_values = fe->get_xyz();
1881  const std::vector<Real> & JxW = fe->get_JxW();
1882 
1883  const std::vector<std::vector<Real> > & phi = fe->get_phi();
1884  const std::vector<std::vector<RealGradient> > * dphi = libmesh_nullptr;
1885  if (cont == C_ONE)
1886  dphi = &(fe->get_dphi());
1887 
1888 #ifdef LIBMESH_ENABLE_AMR
1889  if (elem->refinement_flag() == Elem::JUST_COARSENED)
1890  {
1891  std::vector<Point> fine_points;
1892 
1893  UniquePtr<FEBase> fine_fe
1894  (FEBase::build (dim, base_fe_type));
1895 
1896  UniquePtr<QBase> qrule
1897  (base_fe_type.default_quadrature_rule(dim));
1898  fine_fe->attach_quadrature_rule(qrule.get());
1899 
1900  const std::vector<Point> & child_xyz =
1901  fine_fe->get_xyz();
1902 
1903  for (unsigned int c = 0;
1904  c != elem->n_children(); ++c)
1905  {
1906  fine_fe->reinit(elem->child_ptr(c));
1907  fine_points.insert(fine_points.end(),
1908  child_xyz.begin(),
1909  child_xyz.end());
1910  }
1911 
1912  std::vector<Point> fine_qp;
1913  FEInterface::inverse_map (dim, base_fe_type, elem,
1914  fine_points, fine_qp);
1915 
1916  context.elem_fe_reinit(&fine_qp);
1917  }
1918  else
1919 #endif // LIBMESH_ENABLE_AMR
1920  context.elem_fe_reinit();
1921 
1922  const unsigned int n_qp = xyz_values.size();
1923 
1924  Ke.resize (free_dofs, free_dofs); Ke.zero();
1925  Fe.resize (free_dofs); Fe.zero();
1926  // The new interior coefficients
1927  DenseVector<FValue> Uint(free_dofs);
1928 
1929  // Loop over the quadrature points
1930  for (unsigned int qp=0; qp<n_qp; qp++)
1931  {
1932  // solution at the quadrature point
1933  FValue fineval = f.eval_at_point(context,
1934  var_component,
1935  xyz_values[qp],
1936  system.time);
1937  // solution grad at the quadrature point
1938  VectorValue<FValue> finegrad;
1939  if (cont == C_ONE)
1940  finegrad = g->eval_at_point(context,
1941  var_component,
1942  xyz_values[qp],
1943  system.time);
1944 
1945  // Form interior projection matrix
1946  for (unsigned int i=0, freei=0; i != n_dofs; ++i)
1947  {
1948  // fixed DoFs aren't test functions
1949  if (dof_is_fixed[i])
1950  continue;
1951  for (unsigned int j=0, freej=0; j != n_dofs; ++j)
1952  {
1953  if (dof_is_fixed[j])
1954  Fe(freei) -= phi[i][qp] * phi[j][qp] * JxW[qp]
1955  * Ue(j);
1956  else
1957  Ke(freei,freej) += phi[i][qp] * phi[j][qp] *
1958  JxW[qp];
1959  if (cont == C_ONE)
1960  {
1961  if (dof_is_fixed[j])
1962  Fe(freei) -= ( (*dphi)[i][qp] *
1963  (*dphi)[j][qp] ) * JxW[qp] *
1964  Ue(j);
1965  else
1966  Ke(freei,freej) += ( (*dphi)[i][qp] *
1967  (*dphi)[j][qp] ) *
1968  JxW[qp];
1969  }
1970  if (!dof_is_fixed[j])
1971  freej++;
1972  }
1973  Fe(freei) += phi[i][qp] * fineval * JxW[qp];
1974  if (cont == C_ONE)
1975  Fe(freei) += (finegrad * (*dphi)[i][qp] ) * JxW[qp];
1976  freei++;
1977  }
1978  }
1979  Ke.cholesky_solve(Fe, Uint);
1980 
1981  // Transfer new interior solutions to element
1982  for (unsigned int i=0; i != free_dofs; ++i)
1983  {
1984  FValue & ui = Ue(free_dof[i]);
1986  std::abs(ui - Uint(i)) < TOLERANCE);
1987  ui = Uint(i);
1988  dof_is_fixed[free_dof[i]] = true;
1989  }
1990 
1991  } // if there are free interior dofs
1992 
1993  STOP_LOG ("project_interior","GenericProjector");
1994 
1995  // Make sure every DoF got reached!
1996  for (unsigned int i=0; i != n_dofs; ++i)
1997  libmesh_assert(dof_is_fixed[i]);
1998 
1999  action.insert(context, var, Ue);
2000  } // end variables loop
2001  } // end elements loop
2002 }
static void dofs_on_edge(const Elem *const elem, const unsigned int dim, const FEType &fe_t, unsigned int e, std::vector< unsigned int > &di)
Definition: fe_interface.C:497
double abs(double a)
static void dofs_on_side(const Elem *const elem, const unsigned int dim, const FEType &fe_t, unsigned int s, std::vector< unsigned int > &di)
Definition: fe_interface.C:482
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:2019
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:556
const ProjectionAction & master_action
FEGenericBase< Real > FEBase
const std::vector< unsigned int > & variables
std::vector< object_type >::const_iterator const_iterator
Definition: stored_range.h:58
static unsigned int n_dofs_at_node(const unsigned int dim, const FEType &fe_t, const ElemType t, const unsigned int n)
Definition: fe_interface.C:436
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
unsigned int variable_scalar_number(const std::string &var, unsigned int component) const
Definition: system.h:2134
static UniquePtr< FEGenericBase > build(const unsigned int dim, const FEType &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: