19 #include "libmesh/dof_map.h" 22 #include "libmesh/boundary_info.h" 23 #include "libmesh/dense_matrix.h" 24 #include "libmesh/dense_vector.h" 25 #include "libmesh/dirichlet_boundaries.h" 26 #include "libmesh/elem.h" 27 #include "libmesh/elem_range.h" 28 #include "libmesh/fe_base.h" 29 #include "libmesh/fe_interface.h" 30 #include "libmesh/fe_type.h" 31 #include "libmesh/function_base.h" 32 #include "libmesh/int_range.h" 33 #include "libmesh/libmesh_logging.h" 34 #include "libmesh/linear_solver.h" 35 #include "libmesh/mesh_base.h" 36 #include "libmesh/null_output_iterator.h" 37 #include "libmesh/mesh_tools.h" 38 #include "libmesh/nonlinear_implicit_system.h" 39 #include "libmesh/numeric_vector.h" 40 #include "libmesh/parallel_algebra.h" 41 #include "libmesh/parallel_elem.h" 42 #include "libmesh/parallel_node.h" 43 #include "libmesh/periodic_boundaries.h" 44 #include "libmesh/periodic_boundary.h" 45 #include "libmesh/periodic_boundary_base.h" 46 #include "libmesh/point_locator_base.h" 47 #include "libmesh/quadrature.h" 48 #include "libmesh/raw_accessor.h" 49 #include "libmesh/sparse_matrix.h" 50 #include "libmesh/system.h" 51 #include "libmesh/tensor_tools.h" 52 #include "libmesh/threads.h" 53 #include "libmesh/enum_to_string.h" 54 #include "libmesh/coupling_matrix.h" 68 #include <unordered_set> 75 class ComputeConstraints
80 #ifdef LIBMESH_ENABLE_PERIODIC
84 const unsigned int variable_number) :
85 _constraints(constraints),
87 #ifdef LIBMESH_ENABLE_PERIODIC
88 _periodic_boundaries(periodic_boundaries),
91 _variable_number(variable_number)
96 const Variable & var_description = _dof_map.variable(_variable_number);
98 #ifdef LIBMESH_ENABLE_PERIODIC 99 std::unique_ptr<PointLocatorBase> point_locator;
100 const bool have_periodic_boundaries =
101 !_periodic_boundaries.empty();
102 if (have_periodic_boundaries && !range.
empty())
103 point_locator = _mesh.sub_point_locator();
106 for (
const auto & elem : range)
109 #ifdef LIBMESH_ENABLE_AMR 115 #ifdef LIBMESH_ENABLE_PERIODIC 119 if (have_periodic_boundaries)
122 _periodic_boundaries,
134 #ifdef LIBMESH_ENABLE_PERIODIC 138 const unsigned int _variable_number;
143 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS 144 class ComputeNodeConstraints
148 #ifdef LIBMESH_ENABLE_PERIODIC
152 _node_constraints(node_constraints),
153 #ifdef LIBMESH_ENABLE_PERIODIC
154 _periodic_boundaries(periodic_boundaries),
161 #ifdef LIBMESH_ENABLE_PERIODIC 162 std::unique_ptr<PointLocatorBase> point_locator;
163 bool have_periodic_boundaries = !_periodic_boundaries.empty();
164 if (have_periodic_boundaries && !range.
empty())
165 point_locator = _mesh.sub_point_locator();
168 for (
const auto & elem : range)
170 #ifdef LIBMESH_ENABLE_AMR 173 #ifdef LIBMESH_ENABLE_PERIODIC 177 if (have_periodic_boundaries)
179 _periodic_boundaries,
189 #ifdef LIBMESH_ENABLE_PERIODIC 194 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS 197 #ifdef LIBMESH_ENABLE_DIRICHLET 208 AddConstraint(
DofMap & dof_map_in) : dof_map(dof_map_in) {}
212 const Number constraint_rhs)
const = 0;
215 class AddPrimalConstraint :
public AddConstraint
218 AddPrimalConstraint(
DofMap & dof_map_in) : AddConstraint(dof_map_in) {}
222 const Number constraint_rhs)
const 224 if (!dof_map.is_constrained_dof(dof_number))
225 dof_map.add_constraint_row (dof_number, constraint_row,
226 constraint_rhs,
true);
230 class AddAdjointConstraint :
public AddConstraint
233 const unsigned int qoi_index;
236 AddAdjointConstraint(
DofMap & dof_map_in,
unsigned int qoi_index_in)
237 : AddConstraint(dof_map_in), qoi_index(qoi_index_in) {}
241 const Number constraint_rhs)
const 243 dof_map.add_adjoint_constraint_row
244 (qoi_index, dof_number, constraint_row, constraint_rhs,
256 class ConstrainDirichlet
264 const AddConstraint & add_fn;
278 return std::numeric_limits<Real>::quiet_NaN();
295 return std::numeric_limits<Number>::quiet_NaN();
308 struct SingleElemBoundaryInfo
311 const std::map<
boundary_id_type, std::set<std::pair<unsigned int, DirichletBoundary *>>> & ordered_map_in) :
313 boundary_id_to_ordered_dirichlet_boundaries(ordered_map_in),
321 const std::map<boundary_id_type, std::set<std::pair<unsigned int, DirichletBoundary *>>> & boundary_id_to_ordered_dirichlet_boundaries;
324 unsigned short n_sides;
325 unsigned short n_edges;
331 std::map<const DirichletBoundary *, std::vector<bool>> is_boundary_node_map;
332 std::map<const DirichletBoundary *, std::vector<bool>> is_boundary_side_map;
333 std::map<const DirichletBoundary *, std::vector<bool>> is_boundary_edge_map;
334 std::map<const DirichletBoundary *, std::vector<bool>> is_boundary_shellface_map;
336 std::map<const DirichletBoundary *, std::vector<bool>> is_boundary_nodeset_map;
340 std::set<std::pair<unsigned int, DirichletBoundary *>> ordered_dbs;
349 bool reinit(
const Elem * elem_in)
358 is_boundary_node_map.clear();
359 is_boundary_side_map.clear();
360 is_boundary_edge_map.clear();
361 is_boundary_shellface_map.clear();
362 is_boundary_nodeset_map.clear();
369 bool has_dirichlet_constraint =
false;
373 std::vector<boundary_id_type> ids_vec;
375 for (
unsigned char s = 0; s != n_sides; ++s)
380 bool do_this_side =
false;
381 for (
const auto & bc_id : ids_vec)
383 if (
auto it = boundary_id_to_ordered_dirichlet_boundaries.find(bc_id);
384 it != boundary_id_to_ordered_dirichlet_boundaries.end())
389 ordered_dbs.insert(it->second.begin(), it->second.end());
392 for (
const auto & db_pair : it->second)
399 auto pr = is_boundary_side_map.emplace(db_pair.second, std::vector<bool>(n_sides,
false));
400 pr.first->second[s] =
true;
408 has_dirichlet_constraint =
true;
411 for (
unsigned int n = 0; n !=
n_nodes; ++n)
419 for (
const auto & db_pair : ordered_dbs)
423 if (
auto side_it = is_boundary_side_map.find(db_pair.second);
424 side_it != is_boundary_side_map.end() && side_it->second[s])
426 auto pr = is_boundary_node_map.emplace(db_pair.second, std::vector<bool>(
n_nodes,
false));
427 pr.first->second[n] =
true;
433 for (
unsigned int e = 0; e != n_edges; ++e)
441 for (
const auto & db_pair : ordered_dbs)
445 if (
auto side_it = is_boundary_side_map.find(db_pair.second);
446 side_it != is_boundary_side_map.end() && side_it->second[s])
448 auto pr = is_boundary_edge_map.emplace(db_pair.second, std::vector<bool>(n_edges,
false));
449 pr.first->second[e] =
true;
457 for (
unsigned int n=0; n !=
n_nodes; ++n)
461 for (
const auto & bc_id : ids_vec)
463 if (
auto it = boundary_id_to_ordered_dirichlet_boundaries.find(bc_id);
464 it != boundary_id_to_ordered_dirichlet_boundaries.end())
467 ordered_dbs.insert(it->second.begin(), it->second.end());
470 for (
const auto & db_pair : it->second)
472 auto pr = is_boundary_node_map.emplace(db_pair.second, std::vector<bool>(
n_nodes,
false));
473 pr.first->second[n] =
true;
475 auto pr2 = is_boundary_nodeset_map.emplace(db_pair.second, std::vector<bool>(
n_nodes,
false));
476 pr2.first->second[n] =
true;
479 has_dirichlet_constraint =
true;
486 for (
unsigned short e=0; e != n_edges; ++e)
490 bool do_this_side =
false;
491 for (
const auto & bc_id : ids_vec)
493 if (
auto it = boundary_id_to_ordered_dirichlet_boundaries.find(bc_id);
494 it != boundary_id_to_ordered_dirichlet_boundaries.end())
499 ordered_dbs.insert(it->second.begin(), it->second.end());
502 for (
const auto & db_pair : it->second)
504 auto pr = is_boundary_edge_map.emplace(db_pair.second, std::vector<bool>(n_edges,
false));
505 pr.first->second[e] =
true;
513 has_dirichlet_constraint =
true;
516 for (
unsigned int n = 0; n !=
n_nodes; ++n)
524 for (
const auto & db_pair : ordered_dbs)
528 if (
auto edge_it = is_boundary_edge_map.find(db_pair.second);
529 edge_it != is_boundary_edge_map.end() && edge_it->second[e])
531 auto pr = is_boundary_node_map.emplace(db_pair.second, std::vector<bool>(
n_nodes,
false));
532 pr.first->second[n] =
true;
540 for (
unsigned short shellface=0; shellface != 2; ++shellface)
543 bool do_this_shellface =
false;
545 for (
const auto & bc_id : ids_vec)
547 if (
auto it = boundary_id_to_ordered_dirichlet_boundaries.find(bc_id);
548 it != boundary_id_to_ordered_dirichlet_boundaries.end())
550 has_dirichlet_constraint =
true;
551 do_this_shellface =
true;
554 ordered_dbs.insert(it->second.begin(), it->second.end());
557 for (
const auto & db_pair : it->second)
559 auto pr = is_boundary_shellface_map.emplace(db_pair.second, std::vector<bool>(2,
false));
560 pr.first->second[shellface] =
true;
565 if (do_this_shellface)
568 for (
unsigned int n = 0; n !=
n_nodes; ++n)
569 for (
const auto & db_pair : ordered_dbs)
573 if (
auto side_it = is_boundary_shellface_map.find(db_pair.second);
574 side_it != is_boundary_shellface_map.end() && side_it->second[shellface])
576 auto pr = is_boundary_node_map.emplace(db_pair.second, std::vector<bool>(
n_nodes,
false));
577 pr.first->second[n] =
true;
583 return has_dirichlet_constraint;
590 template<
typename OutputType>
591 void apply_lagrange_dirichlet_impl(
const SingleElemBoundaryInfo & sebi,
597 const Elem * elem = sebi.elem;
627 const unsigned int var_component =
631 const unsigned int var = variable.
number();
637 std::unique_ptr<FEMContext> context;
643 context = std::make_unique<FEMContext>
650 if (f_system && context.get())
651 context->pre_fe_reinit(*f_system, elem);
657 const std::vector<dof_id_type> & dof_indices =
661 const unsigned int n_dofs =
662 cast_int<unsigned int>(dof_indices.size());
665 std::vector<char> dof_is_fixed(n_dofs,
false);
677 unsigned int current_dof = 0;
678 for (
unsigned int n=0; n!= sebi.n_nodes; ++n)
683 const unsigned int nc =
693 auto is_boundary_node_it = sebi.is_boundary_node_map.find(&dirichlet);
694 const bool is_boundary_node =
695 (is_boundary_node_it != sebi.is_boundary_node_map.end() &&
696 is_boundary_node_it->second[n]);
699 auto is_boundary_nodeset_it = sebi.is_boundary_nodeset_map.find(&dirichlet);
700 const bool is_boundary_nodeset =
701 (is_boundary_nodeset_it != sebi.is_boundary_nodeset_map.end() &&
702 is_boundary_nodeset_it->second[n]);
705 if ( !(is_boundary_node || is_boundary_nodeset) )
712 libmesh_assert_equal_to (nc, n_vec_dim);
713 for (
unsigned int c = 0; c < n_vec_dim; c++)
716 f_component(f, f_fem, context.get(), var_component+c,
717 elem->
point(n), time);
718 dof_is_fixed[current_dof+c] =
true;
720 current_dof += n_vec_dim;
727 for (
unsigned int i = 0; i < n_dofs; i++)
731 add_fn (dof_indices[i], empty_row, Ue(i));
739 template<
typename OutputType>
740 void apply_dirichlet_impl(
const SingleElemBoundaryInfo & sebi,
746 const Elem * elem = sebi.elem;
753 typedef OutputType OutputShape;
794 const unsigned int var_component =
798 const unsigned int var = variable.
number();
820 std::unique_ptr<FEMContext> context;
826 context = std::make_unique<FEMContext>
847 if (f_system && context.get())
848 context->pre_fe_reinit(*f_system, elem);
854 const std::vector<dof_id_type> & dof_indices =
858 const unsigned int n_dofs =
859 cast_int<unsigned int>(dof_indices.size());
862 std::vector<char> dof_is_fixed(n_dofs,
false);
863 std::vector<int> free_dof(n_dofs, 0);
881 unsigned int current_dof = 0;
882 for (
unsigned int n=0; n!= sebi.n_nodes; ++n)
889 const unsigned int nc =
896 auto is_boundary_node_it = sebi.is_boundary_node_map.find(&dirichlet);
902 const bool not_boundary_node =
903 (is_boundary_node_it == sebi.is_boundary_node_map.end() ||
904 !is_boundary_node_it->second[n]);
906 if (!elem->
is_vertex(n) || not_boundary_node)
913 libmesh_assert_equal_to (nc, 0);
919 libmesh_assert_equal_to (nc, n_vec_dim);
920 for (
unsigned int c = 0; c < n_vec_dim; c++)
923 f_component(f, f_fem, context.get(), var_component+c,
924 elem->
point(n), time);
925 dof_is_fixed[current_dof+c] =
true;
927 current_dof += n_vec_dim;
933 f_component(f, f_fem, context.get(), var_component,
934 elem->
point(n), time);
935 dof_is_fixed[current_dof] =
true;
938 g_component(g, g_fem, context.get(), var_component,
939 elem->
point(n), time);
941 Ue(current_dof) = grad(0);
942 dof_is_fixed[current_dof] =
true;
948 nxplus = elem->
point(n);
952 g_component(g, g_fem, context.get(), var_component,
955 g_component(g, g_fem, context.get(), var_component,
958 Ue(current_dof) = grad(1);
959 dof_is_fixed[current_dof] =
true;
962 Ue(current_dof) = (gxplus(1) - gxminus(1))
964 dof_is_fixed[current_dof] =
true;
970 Ue(current_dof) = grad(2);
971 dof_is_fixed[current_dof] =
true;
974 Ue(current_dof) = (gxplus(2) - gxminus(2))
976 dof_is_fixed[current_dof] =
true;
980 nyplus = elem->
point(n);
984 g_component(g, g_fem, context.get(), var_component,
987 g_component(g, g_fem, context.get(), var_component,
990 Ue(current_dof) = (gyplus(2) - gyminus(2))
992 dof_is_fixed[current_dof] =
true;
996 nxmyp = elem->
point(n),
997 nxpym = elem->
point(n),
998 nxpyp = elem->
point(n);
1008 g_component(g, g_fem, context.get(), var_component,
1011 g_component(g, g_fem, context.get(), var_component,
1014 g_component(g, g_fem, context.get(), var_component,
1017 g_component(g, g_fem, context.get(), var_component,
1019 Number gxzplus = (gxpyp(2) - gxmyp(2))
1021 Number gxzminus = (gxpym(2) - gxmym(2))
1024 Ue(current_dof) = (gxzplus - gxzminus)
1026 dof_is_fixed[current_dof] =
true;
1034 else if (cont ==
C_ONE)
1036 libmesh_assert_equal_to (nc, 1 +
dim);
1038 f_component(f, f_fem, context.get(), var_component,
1039 elem->
point(n), time);
1040 dof_is_fixed[current_dof] =
true;
1043 g_component(g, g_fem, context.get(), var_component,
1044 elem->
point(n), time);
1045 for (
unsigned int i=0; i!=
dim; ++i)
1047 Ue(current_dof) = grad(i);
1048 dof_is_fixed[current_dof] =
true;
1053 libmesh_error_msg(
"Unknown continuity cont = " << cont);
1073 const std::vector<std::vector<OutputShape>> & phi = edge_fe->
get_phi();
1074 const std::vector<Point> & xyz_values = edge_fe->
get_xyz();
1075 const std::vector<Real> & JxW = edge_fe->
get_JxW();
1078 const std::vector<std::vector<OutputGradient>> * dphi =
nullptr;
1081 const std::vector<std::vector<OutputGradient>> & ref_dphi = edge_fe->
get_dphi();
1086 std::vector<unsigned int> edge_dofs;
1092 if (
const auto is_boundary_edge_it = sebi.is_boundary_edge_map.find(&dirichlet);
1093 is_boundary_edge_it != sebi.is_boundary_edge_map.end())
1095 for (
unsigned int e=0; e != sebi.n_edges; ++e)
1097 if (!is_boundary_edge_it->second[e])
1103 const unsigned int n_edge_dofs =
1104 cast_int<unsigned int>(edge_dofs.size());
1108 unsigned int free_dofs = 0;
1109 for (
unsigned int i=0; i != n_edge_dofs; ++i)
1110 if (!dof_is_fixed[edge_dofs[i]])
1111 free_dof[free_dofs++] = i;
1127 for (
unsigned int qp=0; qp<n_qp; qp++)
1130 OutputNumber fineval(0);
1133 for (
unsigned int c = 0; c < n_vec_dim; c++)
1135 f_component(f, f_fem, context.get(), var_component+c,
1136 xyz_values[qp], time);
1139 OutputNumberGradient finegrad;
1142 unsigned int g_rank;
1156 libmesh_error_msg(
"Unknown field type!");
1160 for (
unsigned int c = 0; c < n_vec_dim; c++)
1161 for (
unsigned int d = 0; d < g_rank; d++)
1162 g_accessor(c + d*
dim ) =
1163 g_component(g, g_fem, context.get(), var_component,
1164 xyz_values[qp], time)(c);
1167 for (
unsigned int sidei=0, freei=0; sidei != n_edge_dofs; ++sidei)
1169 unsigned int i = edge_dofs[sidei];
1171 if (dof_is_fixed[i])
1173 for (
unsigned int sidej=0, freej=0; sidej != n_edge_dofs; ++sidej)
1175 unsigned int j = edge_dofs[sidej];
1176 if (dof_is_fixed[j])
1177 Fe(freei) -= phi[i][qp] * phi[j][qp] *
1180 Ke(freei,freej) += phi[i][qp] *
1181 phi[j][qp] * JxW[qp];
1184 if (dof_is_fixed[j])
1185 Fe(freei) -= ((*dphi)[i][qp].contract((*dphi)[j][qp]) ) *
1188 Ke(freei,freej) += ((*dphi)[i][qp].contract((*dphi)[j][qp]))
1191 if (!dof_is_fixed[j])
1194 Fe(freei) += phi[i][qp] * fineval * JxW[qp];
1196 Fe(freei) += (finegrad.contract( (*dphi)[i][qp]) ) *
1205 for (
unsigned int i=0; i != free_dofs; ++i)
1207 Number & ui = Ue(edge_dofs[free_dof[i]]);
1211 dof_is_fixed[edge_dofs[free_dof[i]]] =
true;
1231 const std::vector<std::vector<OutputShape>> & phi = side_fe->
get_phi();
1232 const std::vector<Point> & xyz_values = side_fe->
get_xyz();
1233 const std::vector<Real> & JxW = side_fe->
get_JxW();
1236 const std::vector<std::vector<OutputGradient>> * dphi =
nullptr;
1239 const std::vector<std::vector<OutputGradient>> & ref_dphi = side_fe->
get_dphi();
1244 std::vector<unsigned int> side_dofs;
1250 if (
const auto is_boundary_side_it = sebi.is_boundary_side_map.find(&dirichlet);
1251 is_boundary_side_it != sebi.is_boundary_side_map.end())
1253 for (
unsigned int s=0; s != sebi.n_sides; ++s)
1255 if (!is_boundary_side_it->second[s])
1261 const unsigned int n_side_dofs =
1262 cast_int<unsigned int>(side_dofs.size());
1266 unsigned int free_dofs = 0;
1267 for (
unsigned int i=0; i != n_side_dofs; ++i)
1268 if (!dof_is_fixed[side_dofs[i]])
1269 free_dof[free_dofs++] = i;
1281 side_fe->
reinit(elem, s);
1285 for (
unsigned int qp=0; qp<n_qp; qp++)
1288 OutputNumber fineval(0);
1291 for (
unsigned int c = 0; c < n_vec_dim; c++)
1293 f_component(f, f_fem, context.get(), var_component+c,
1294 xyz_values[qp], time);
1297 OutputNumberGradient finegrad;
1300 unsigned int g_rank;
1314 libmesh_error_msg(
"Unknown field type!");
1318 for (
unsigned int c = 0; c < n_vec_dim; c++)
1319 for (
unsigned int d = 0; d < g_rank; d++)
1320 g_accessor(c + d*
dim ) =
1321 g_component(g, g_fem, context.get(), var_component,
1322 xyz_values[qp], time)(c);
1325 for (
unsigned int sidei=0, freei=0; sidei != n_side_dofs; ++sidei)
1327 unsigned int i = side_dofs[sidei];
1329 if (dof_is_fixed[i])
1331 for (
unsigned int sidej=0, freej=0; sidej != n_side_dofs; ++sidej)
1333 unsigned int j = side_dofs[sidej];
1334 if (dof_is_fixed[j])
1335 Fe(freei) -= phi[i][qp] * phi[j][qp] *
1338 Ke(freei,freej) += phi[i][qp] *
1339 phi[j][qp] * JxW[qp];
1342 if (dof_is_fixed[j])
1343 Fe(freei) -= ((*dphi)[i][qp].contract((*dphi)[j][qp])) *
1346 Ke(freei,freej) += ((*dphi)[i][qp].contract((*dphi)[j][qp]))
1349 if (!dof_is_fixed[j])
1352 Fe(freei) += (fineval * phi[i][qp]) * JxW[qp];
1354 Fe(freei) += (finegrad.contract((*dphi)[i][qp])) *
1363 for (
unsigned int i=0; i != free_dofs; ++i)
1365 Number & ui = Ue(side_dofs[free_dof[i]]);
1371 dof_is_fixed[side_dofs[free_dof[i]]] =
true;
1391 const std::vector<std::vector<OutputShape>> & phi = fe->
get_phi();
1392 const std::vector<Point> & xyz_values = fe->
get_xyz();
1393 const std::vector<Real> & JxW = fe->
get_JxW();
1396 const std::vector<std::vector<OutputGradient>> * dphi =
nullptr;
1399 const std::vector<std::vector<OutputGradient>> & ref_dphi = fe->
get_dphi();
1407 if (
const auto is_boundary_shellface_it = sebi.is_boundary_shellface_map.find(&dirichlet);
1408 is_boundary_shellface_it != sebi.is_boundary_shellface_map.end())
1410 for (
unsigned int shellface=0; shellface != 2; ++shellface)
1412 if (!is_boundary_shellface_it->second[shellface])
1416 std::vector<unsigned int> shellface_dofs(n_dofs);
1417 std::iota(shellface_dofs.begin(), shellface_dofs.end(), 0);
1421 unsigned int free_dofs = 0;
1422 for (
unsigned int i=0; i != n_dofs; ++i)
1423 if (!dof_is_fixed[shellface_dofs[i]])
1424 free_dof[free_dofs++] = i;
1440 for (
unsigned int qp=0; qp<n_qp; qp++)
1443 OutputNumber fineval(0);
1446 for (
unsigned int c = 0; c < n_vec_dim; c++)
1448 f_component(f, f_fem, context.get(), var_component+c,
1449 xyz_values[qp], time);
1452 OutputNumberGradient finegrad;
1455 unsigned int g_rank;
1469 libmesh_error_msg(
"Unknown field type!");
1473 for (
unsigned int c = 0; c < n_vec_dim; c++)
1474 for (
unsigned int d = 0; d < g_rank; d++)
1475 g_accessor(c + d*
dim ) =
1476 g_component(g, g_fem, context.get(), var_component,
1477 xyz_values[qp], time)(c);
1480 for (
unsigned int shellfacei=0, freei=0;
1481 shellfacei != n_dofs; ++shellfacei)
1483 unsigned int i = shellface_dofs[shellfacei];
1485 if (dof_is_fixed[i])
1487 for (
unsigned int shellfacej=0, freej=0;
1488 shellfacej != n_dofs; ++shellfacej)
1490 unsigned int j = shellface_dofs[shellfacej];
1491 if (dof_is_fixed[j])
1492 Fe(freei) -= phi[i][qp] * phi[j][qp] *
1495 Ke(freei,freej) += phi[i][qp] *
1496 phi[j][qp] * JxW[qp];
1499 if (dof_is_fixed[j])
1500 Fe(freei) -= ((*dphi)[i][qp].contract((*dphi)[j][qp])) *
1503 Ke(freei,freej) += ((*dphi)[i][qp].contract((*dphi)[j][qp]))
1506 if (!dof_is_fixed[j])
1509 Fe(freei) += (fineval * phi[i][qp]) * JxW[qp];
1511 Fe(freei) += (finegrad.contract((*dphi)[i][qp])) *
1520 for (
unsigned int i=0; i != free_dofs; ++i)
1522 Number & ui = Ue(shellface_dofs[free_dof[i]]);
1524 std::abs(ui - Ushellface(i)) <
TOLERANCE);
1526 dof_is_fixed[shellface_dofs[free_dof[i]]] =
true;
1536 for (
unsigned int i = 0; i < n_dofs; i++)
1540 add_fn (dof_indices[i], empty_row, Ue(i));
1546 ConstrainDirichlet (
DofMap & dof_map_in,
1550 const AddConstraint & add_in) :
1551 dof_map(dof_map_in),
1554 dirichlets(dirichlets_in),
1558 ConstrainDirichlet (ConstrainDirichlet &&) =
default;
1559 ConstrainDirichlet (
const ConstrainDirichlet &) =
default;
1563 ConstrainDirichlet & operator= (
const ConstrainDirichlet &) =
delete;
1564 ConstrainDirichlet & operator= (ConstrainDirichlet &&) =
delete;
1578 System * system =
nullptr;
1583 std::map<boundary_id_type, std::set<std::pair<unsigned int, DirichletBoundary *>>>
1584 boundary_id_to_ordered_dirichlet_boundaries;
1589 const auto & dirichlet = dirichlets[dirichlet_id];
1592 for (
const auto & b_id : dirichlet->
b)
1593 boundary_id_to_ordered_dirichlet_boundaries[b_id].emplace(dirichlet_id, dirichlet.get());
1595 for (
const auto & var : dirichlet->
variables)
1598 auto current_system = variable.
system();
1601 system = current_system;
1603 libmesh_error_msg_if(current_system != system,
1604 "All variables should be defined on the same System");
1613 libmesh_error_msg_if(!system,
"Valid System not found for any Variables.");
1620 auto fem_context = std::make_unique<FEMContext>
1621 (*system,
nullptr,
false);
1632 SingleElemBoundaryInfo sebi(boundary_info, boundary_id_to_ordered_dirichlet_boundaries);
1635 for (
const auto & elem : range)
1643 bool has_dirichlet_constraint = sebi.reinit(elem);
1646 if (!has_dirichlet_constraint)
1649 for (
const auto & db_pair : sebi.ordered_dbs)
1652 const auto & dirichlet = db_pair.second;
1655 for (
const auto & var : dirichlet->
variables)
1661 libmesh_assert_equal_to(variable.
number(), var);
1665 if (fe_type.family ==
SCALAR)
1676 this->apply_lagrange_dirichlet_impl<Real>( sebi, variable, *dirichlet, *fem_context );
1678 this->apply_dirichlet_impl<Real>( sebi, variable, *dirichlet, *fem_context );
1683 this->apply_dirichlet_impl<RealGradient>( sebi, variable, *dirichlet, *fem_context );
1687 libmesh_error_msg(
"Unknown field type!");
1697 #endif // LIBMESH_ENABLE_DIRICHLET 1710 #ifdef LIBMESH_ENABLE_CONSTRAINTS 1715 parallel_object_only();
1725 const DofConstraints::const_iterator lower =
1737 parallel_object_only();
1739 LOG_SCOPE(
"create_dof_constraints()",
"DofMap");
1754 this->
comm().
min(constraint_rows_empty);
1759 const bool possible_local_constraints =
false 1761 || !constraint_rows_empty
1762 #ifdef LIBMESH_ENABLE_AMR 1765 #ifdef LIBMESH_ENABLE_PERIODIC 1768 #ifdef LIBMESH_ENABLE_DIRICHLET 1774 bool possible_global_constraints = possible_local_constraints;
1775 #if defined(LIBMESH_ENABLE_PERIODIC) || defined(LIBMESH_ENABLE_DIRICHLET) || defined(LIBMESH_ENABLE_AMR) 1778 this->
comm().
max(possible_global_constraints);
1786 #ifdef LIBMESH_ENABLE_CONSTRAINTS 1791 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS 1795 if (!possible_global_constraints)
1805 mesh.local_elements_end());
1816 #ifdef LIBMESH_ENABLE_PERIODIC 1819 this->
comm().
max(need_point_locator);
1821 if (need_point_locator)
1825 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS 1828 #ifdef LIBMESH_ENABLE_PERIODIC
1832 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS 1838 for (
unsigned int variable_number=0; variable_number<
n_vars;
1839 ++variable_number, range.reset())
1843 #ifdef LIBMESH_ENABLE_PERIODIC
1849 #ifdef LIBMESH_ENABLE_DIRICHLET 1868 AddPrimalConstraint(*
this)));
1882 ConstrainDirichlet(*
this,
mesh, time, adb_q,
1883 AddAdjointConstraint(*
this, qoi_index)));
1887 #endif // LIBMESH_ENABLE_DIRICHLET 1891 if (!constraint_rows_empty)
1912 bool constraint_rows_empty = constraint_rows.empty();
1913 this->
comm().
min(constraint_rows_empty);
1919 #ifdef LIBMESH_ENABLE_PERIODIC 1921 "Periodic boundary conditions are not yet implemented for spline meshes");
1926 std::unique_ptr<SparsityPattern::Build> sp;
1927 std::unique_ptr<SparseMatrix<Number>> mat;
1929 const unsigned int n_adjoint_rhs =
1933 std::vector<std::unique_ptr<NumericVector<Number>>>
1934 solve_rhs(n_adjoint_rhs+1);
1939 std::set<dof_id_type> my_dirichlet_spline_dofs;
1942 std::unordered_set<dof_id_type> was_previously_constrained;
1944 const unsigned int sys_num = this->
sys_number();
1945 for (
auto & node_row : constraint_rows)
1947 const Node * node = node_row.first;
1969 for (
const auto & [pr, val] : node_row.second)
1971 const Elem * spline_elem = pr.first;
1974 const Node & spline_node =
1979 dc_row[spline_dof_id] = val;
1985 was_previously_constrained.insert(constrained_id);
1989 for (
auto & row_entry : dc_row)
1990 my_dirichlet_spline_dofs.insert(row_entry.first);
2009 if (this->
comm().size() > 1)
2013 their_dirichlet_spline_dofs;
2019 my_dirichlet_spline_dofs.end()));
2021 for (
auto d : my_dirichlet_spline_dofs)
2023 libmesh_assert_less(d, this->
end_dof(this->
comm().size()-1));
2024 while (d >= this->
end_dof(destination_pid))
2028 their_dirichlet_spline_dofs[destination_pid].push_back(d);
2031 auto receive_dof_functor =
2032 [& my_dirichlet_spline_dofs]
2034 const std::vector<dof_id_type> & dofs)
2036 my_dirichlet_spline_dofs.insert(dofs.begin(), dofs.end());
2039 Parallel::push_parallel_vector_data
2040 (this->
comm(), their_dirichlet_spline_dofs, receive_dof_functor);
2050 bool important_prior_constraints =
2051 !was_previously_constrained.empty();
2052 this->
comm().
max(important_prior_constraints);
2054 if (important_prior_constraints)
2070 mat->attach_dof_map(*
this);
2072 mat->attach_sparsity_pattern(*sp);
2075 for (
auto & node_row : constraint_rows)
2077 const Node * node = node_row.first;
2092 if (was_previously_constrained.count(constrained_id))
2098 std::vector<dof_id_type>
dof_indices(1, constrained_id);
2106 DofConstraintValueMap::const_iterator rhsit =
2107 vals.find(constrained_id);
2108 F(0) = (rhsit == vals.end()) ? 0 : rhsit->second;
2111 if (rhsit != vals.end())
2129 if (!was_previously_constrained.count(d) &&
2130 !my_dirichlet_spline_dofs.count(d))
2135 std::unique_ptr<LinearSolver<Number>> linear_solver =
2138 std::unique_ptr<NumericVector<Number>> projected_vals =
2145 for (
auto sd : my_dirichlet_spline_dofs)
2154 const unsigned int max_its = 5000;
2156 linear_solver->solve(*mat, *projected_vals,
2157 *(solve_rhs[q]), tol, max_its);
2163 for (
auto sd : my_dirichlet_spline_dofs)
2166 Number constraint_rhs = (*projected_vals)(sd);
2168 std::pair<DofConstraintValueMap::iterator, bool> rhs_it =
2169 vals.emplace(sd, constraint_rhs);
2171 rhs_it.first->second = constraint_rhs;
2181 const Number constraint_rhs,
2182 const bool forbid_constraint_overwrite)
2185 libmesh_error_msg_if(forbid_constraint_overwrite && this->
is_constrained_dof(dof_number),
2186 "ERROR: DOF " << dof_number <<
" was already constrained!");
2188 libmesh_assert_less(dof_number, this->
n_dofs());
2192 libmesh_assert_msg(!constraint_row.count(dof_number),
2193 "Error: constraint_row for dof " << dof_number <<
2194 " should not contain an entry for dof " << dof_number);
2197 for (
const auto & pr : constraint_row)
2198 libmesh_assert_less(pr.first, this->n_dofs());
2204 std::pair<DofConstraintValueMap::iterator, bool> rhs_it =
2207 rhs_it.first->second = constraint_rhs;
2214 const Number constraint_rhs,
2215 const bool forbid_constraint_overwrite)
2218 if (forbid_constraint_overwrite)
2221 "ERROR: DOF " << dof_number <<
" has no corresponding primal constraint!");
2253 bool print_nonlocal)
const 2255 parallel_object_only();
2257 std::string local_constraints =
2262 this->
comm().
send(0, local_constraints);
2266 os <<
"Processor 0:\n";
2267 os << local_constraints;
2272 os <<
"Processor " << p <<
":\n";
2273 os << local_constraints;
2280 std::ostringstream os;
2281 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS 2287 os <<
"Node Constraints:" 2293 if (!print_nonlocal &&
2298 const Point & offset = pr.second;
2300 os <<
"Constraints for Node id " << node->id()
2303 for (
const auto & [cnode, val] : row)
2304 os <<
" (" << cnode->id() <<
"," << val <<
")\t";
2306 os <<
"rhs: " << offset;
2310 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS 2317 os <<
"DoF Constraints:" 2326 DofConstraintValueMap::const_iterator rhsit =
2331 os <<
"Constraints for DoF " << i
2334 for (
const auto & item : row)
2335 os <<
" (" << item.first <<
"," << item.second <<
")\t";
2337 os <<
"rhs: " << rhs;
2341 for (
unsigned int qoi_index = 0,
2343 qoi_index != n_qois; ++qoi_index)
2345 os <<
"Adjoint " << qoi_index <<
" DoF rhs values:" 2350 for (
const auto & [i, rhs] : adjoint_map_it->second)
2356 os <<
"RHS for DoF " << i
2369 std::vector<dof_id_type> & elem_dofs,
2370 bool asymmetric_constraint_rows)
const 2372 libmesh_assert_equal_to (elem_dofs.size(), matrix.
m());
2373 libmesh_assert_equal_to (elem_dofs.size(), matrix.
n());
2385 LOG_SCOPE(
"constrain_elem_matrix()",
"DofMap");
2388 if ((C.
m() == matrix.
m()) &&
2389 (C.
n() == elem_dofs.size()))
2396 libmesh_assert_equal_to (matrix.
m(), matrix.
n());
2397 libmesh_assert_equal_to (matrix.
m(), elem_dofs.size());
2398 libmesh_assert_equal_to (matrix.
n(), elem_dofs.size());
2401 for (
unsigned int i=0,
2402 n_elem_dofs = cast_int<unsigned int>(elem_dofs.size());
2403 i != n_elem_dofs; i++)
2412 if (asymmetric_constraint_rows)
2414 DofConstraints::const_iterator
2427 for (
const auto & item : constraint_row)
2428 for (
unsigned int j=0; j != n_elem_dofs; j++)
2429 if (elem_dofs[j] == item.first)
2430 matrix(i,j) = -item.second;
2440 std::vector<dof_id_type> & elem_dofs,
2441 bool asymmetric_constraint_rows)
const 2443 libmesh_assert_equal_to (elem_dofs.size(), matrix.
m());
2444 libmesh_assert_equal_to (elem_dofs.size(), matrix.
n());
2445 libmesh_assert_equal_to (elem_dofs.size(), rhs.
size());
2457 LOG_SCOPE(
"cnstrn_elem_mat_vec()",
"DofMap");
2460 if ((C.
m() == matrix.
m()) &&
2461 (C.
n() == elem_dofs.size()))
2468 libmesh_assert_equal_to (matrix.
m(), matrix.
n());
2469 libmesh_assert_equal_to (matrix.
m(), elem_dofs.size());
2470 libmesh_assert_equal_to (matrix.
n(), elem_dofs.size());
2473 for (
unsigned int i=0,
2474 n_elem_dofs = cast_int<unsigned int>(elem_dofs.size());
2475 i != n_elem_dofs; i++)
2487 if (asymmetric_constraint_rows)
2489 DofConstraints::const_iterator
2499 for (
const auto & item : constraint_row)
2500 for (
unsigned int j=0; j != n_elem_dofs; j++)
2501 if (elem_dofs[j] == item.first)
2502 matrix(i,j) = -item.second;
2519 std::vector<dof_id_type> & elem_dofs,
2520 bool asymmetric_constraint_rows,
2521 int qoi_index)
const 2523 libmesh_assert_equal_to (elem_dofs.size(), matrix.
m());
2524 libmesh_assert_equal_to (elem_dofs.size(), matrix.
n());
2525 libmesh_assert_equal_to (elem_dofs.size(), rhs.
size());
2538 LOG_SCOPE(
"hetero_cnstrn_elem_mat_vec()",
"DofMap");
2541 if ((C.
m() == matrix.
m()) &&
2542 (C.
n() == elem_dofs.size()))
2550 rhs_values = &it->second;
2565 libmesh_assert_equal_to (matrix.
m(), matrix.
n());
2566 libmesh_assert_equal_to (matrix.
m(), elem_dofs.size());
2567 libmesh_assert_equal_to (matrix.
n(), elem_dofs.size());
2569 for (
unsigned int i=0,
2570 n_elem_dofs = cast_int<unsigned int>(elem_dofs.size());
2571 i != n_elem_dofs; i++)
2586 if (asymmetric_constraint_rows)
2588 DofConstraints::const_iterator
2595 for (
const auto & item : constraint_row)
2596 for (
unsigned int j=0; j != n_elem_dofs; j++)
2597 if (elem_dofs[j] == item.first)
2598 matrix(i,j) = -item.second;
2602 const DofConstraintValueMap::const_iterator valpos =
2603 rhs_values->find(dof_id);
2605 rhs(i) = (valpos == rhs_values->end()) ?
2621 std::vector<dof_id_type> & elem_dofs,
2624 libmesh_assert_equal_to (elem_dofs.size(), matrix.
m());
2625 libmesh_assert_equal_to (elem_dofs.size(), matrix.
n());
2626 libmesh_assert_equal_to (elem_dofs.size(), rhs.
size());
2632 if (this->_dof_constraints.empty())
2641 this->build_constraint_matrix_and_vector (C, H, elem_dofs);
2643 LOG_SCOPE(
"hetero_cnstrn_elem_jac_res()",
"DofMap");
2646 if ((C.
m() != matrix.
m()) ||
2647 (C.
n() != elem_dofs.size()))
2658 libmesh_assert_equal_to (matrix.
m(), matrix.
n());
2659 libmesh_assert_equal_to (matrix.
m(), elem_dofs.size());
2660 libmesh_assert_equal_to (matrix.
n(), elem_dofs.size());
2662 for (
unsigned int i=0,
2663 n_elem_dofs = cast_int<unsigned int>(elem_dofs.size());
2664 i != n_elem_dofs; i++)
2668 if (
auto pos = _dof_constraints.find(dof_id);
2669 pos != _dof_constraints.end())
2682 for (
const auto & item : constraint_row)
2683 for (
unsigned int j=0; j != n_elem_dofs; j++)
2684 if (elem_dofs[j] == item.first)
2685 matrix(i,j) = -item.second;
2687 const DofConstraintValueMap::const_iterator valpos =
2688 _primal_constraint_values.find(dof_id);
2690 Number & rhs_val = rhs(i);
2691 rhs_val = (valpos == _primal_constraint_values.end()) ?
2692 0 : -valpos->second;
2693 for (
const auto & [constraining_dof, coef] : constraint_row)
2694 rhs_val -= coef * solution_local(constraining_dof);
2695 rhs_val += solution_local(dof_id);
2703 std::vector<dof_id_type> & elem_dofs,
2706 libmesh_assert_equal_to (elem_dofs.size(), rhs.
size());
2712 if (this->_dof_constraints.empty())
2720 this->build_constraint_matrix_and_vector (C, H, elem_dofs);
2722 LOG_SCOPE(
"hetero_cnstrn_elem_res()",
"DofMap");
2725 if ((C.
m() != rhs.
size()) ||
2726 (C.
n() != elem_dofs.size()))
2733 for (
unsigned int i=0,
2734 n_elem_dofs = cast_int<unsigned int>(elem_dofs.size());
2735 i != n_elem_dofs; i++)
2739 if (
auto pos = _dof_constraints.find(dof_id);
2740 pos != _dof_constraints.end())
2747 const DofConstraintValueMap::const_iterator valpos =
2748 _primal_constraint_values.find(dof_id);
2750 Number & rhs_val = rhs(i);
2751 rhs_val = (valpos == _primal_constraint_values.end()) ?
2752 0 : -valpos->second;
2753 for (
const auto & [constraining_dof, coef] : constraint_row)
2754 rhs_val -= coef * solution_local(constraining_dof);
2755 rhs_val += solution_local(dof_id);
2763 std::vector<dof_id_type> & elem_dofs,
2766 libmesh_assert_equal_to (elem_dofs.size(), rhs.
size());
2772 if (this->_dof_constraints.empty())
2778 this->build_constraint_matrix (C, elem_dofs);
2780 LOG_SCOPE(
"cnstrn_elem_residual()",
"DofMap");
2783 if (C.
n() != elem_dofs.size())
2790 for (
unsigned int i=0,
2791 n_elem_dofs = cast_int<unsigned int>(elem_dofs.size());
2792 i != n_elem_dofs; i++)
2796 if (
auto pos = _dof_constraints.find(dof_id);
2797 pos != _dof_constraints.end())
2804 Number & rhs_val = rhs(i);
2806 for (
const auto & [constraining_dof, coef] : constraint_row)
2807 rhs_val -= coef * solution_local(constraining_dof);
2808 rhs_val += solution_local(dof_id);
2816 std::vector<dof_id_type> & elem_dofs,
2817 bool asymmetric_constraint_rows,
2818 int qoi_index)
const 2820 libmesh_assert_equal_to (elem_dofs.size(), matrix.
m());
2821 libmesh_assert_equal_to (elem_dofs.size(), matrix.
n());
2822 libmesh_assert_equal_to (elem_dofs.size(), rhs.
size());
2835 LOG_SCOPE(
"hetero_cnstrn_elem_vec()",
"DofMap");
2838 if ((C.
m() == matrix.
m()) &&
2839 (C.
n() == elem_dofs.size()))
2847 rhs_values = &it->second;
2858 for (
unsigned int i=0,
2859 n_elem_dofs = cast_int<unsigned int>(elem_dofs.size());
2860 i != n_elem_dofs; i++)
2869 if (asymmetric_constraint_rows && rhs_values)
2871 const DofConstraintValueMap::const_iterator valpos =
2872 rhs_values->find(dof_id);
2874 rhs(i) = (valpos == rhs_values->end()) ?
2889 std::vector<dof_id_type> & row_dofs,
2890 std::vector<dof_id_type> & col_dofs,
2891 bool asymmetric_constraint_rows)
const 2893 libmesh_assert_equal_to (row_dofs.size(), matrix.
m());
2894 libmesh_assert_equal_to (col_dofs.size(), matrix.
n());
2907 std::vector<dof_id_type> orig_row_dofs(row_dofs);
2908 std::vector<dof_id_type> orig_col_dofs(col_dofs);
2913 LOG_SCOPE(
"constrain_elem_matrix()",
"DofMap");
2915 row_dofs = orig_row_dofs;
2916 col_dofs = orig_col_dofs;
2918 bool constraint_found =
false;
2922 if ((R.
m() == matrix.
m()) &&
2923 (R.
n() == row_dofs.size()))
2926 constraint_found =
true;
2929 if ((C.
m() == matrix.
n()) &&
2930 (C.
n() == col_dofs.size()))
2933 constraint_found =
true;
2937 if (constraint_found)
2939 libmesh_assert_equal_to (matrix.
m(), row_dofs.size());
2940 libmesh_assert_equal_to (matrix.
n(), col_dofs.size());
2943 for (
unsigned int i=0,
2944 n_row_dofs = cast_int<unsigned int>(row_dofs.size());
2945 i != n_row_dofs; i++)
2950 if (row_dofs[i] != col_dofs[j])
2956 if (asymmetric_constraint_rows)
2958 DofConstraints::const_iterator
2967 for (
const auto & item : constraint_row)
2968 for (
unsigned int j=0,
2969 n_col_dofs = cast_int<unsigned int>(col_dofs.size());
2970 j != n_col_dofs; j++)
2971 if (col_dofs[j] == item.first)
2972 matrix(i,j) = -item.second;
2981 std::vector<dof_id_type> & row_dofs,
2984 libmesh_assert_equal_to (rhs.
size(), row_dofs.size());
2995 LOG_SCOPE(
"constrain_elem_vector()",
"DofMap");
2998 if ((R.
m() == rhs.
size()) &&
2999 (R.
n() == row_dofs.size()))
3005 libmesh_assert_equal_to (row_dofs.size(), rhs.
size());
3007 for (
unsigned int i=0,
3008 n_row_dofs = cast_int<unsigned int>(row_dofs.size());
3009 i != n_row_dofs; i++)
3024 std::vector<dof_id_type> & row_dofs,
3027 libmesh_assert_equal_to (v.
size(), row_dofs.size());
3028 libmesh_assert_equal_to (w.
size(), row_dofs.size());
3039 LOG_SCOPE(
"cnstrn_elem_dyad_mat()",
"DofMap");
3042 if ((R.
m() == v.
size()) &&
3043 (R.
n() == row_dofs.size()))
3053 libmesh_assert_equal_to (row_dofs.size(), v.
size());
3054 libmesh_assert_equal_to (row_dofs.size(), w.
size());
3058 for (
unsigned int i=0,
3059 n_row_dofs = cast_int<unsigned int>(row_dofs.size());
3060 i != n_row_dofs; i++)
3089 bool homogeneous)
const 3091 parallel_object_only();
3096 LOG_SCOPE(
"enforce_constraints_exactly()",
"DofMap");
3106 std::unique_ptr<NumericVector<Number>> v_built;
3114 i<v_built->last_local_index(); i++)
3115 v_built->set(i, (*v)(i));
3117 v_global = v_built.
get();
3130 v_local = v_built.
get();
3147 libmesh_assert_equal_to (
this, &(system.
get_dof_map()));
3161 for (
const auto & [dof, val] : constraint_row)
3182 bool homogeneous)
const 3184 parallel_object_only();
3195 std::unique_ptr<NumericVector<Number>> solution_built;
3197 solution_local = solution;
3204 solution_built->close();
3205 solution_local = solution_built.
get();
3213 libmesh_assert_equal_to (
this, &(system.
get_dof_map()));
3221 for (
const auto & [dof, val] : constraint_row)
3223 exact_value += (*solution_local)(constrained_dof);
3238 parallel_object_only();
3246 libmesh_assert_equal_to (
this, &(system.
get_dof_map()));
3253 for (
const auto & j : constraint_row)
3254 jac->
set(constrained_dof, j.first, -j.second);
3255 jac->
set(constrained_dof, constrained_dof, 1);
3261 unsigned int q)
const 3263 parallel_object_only();
3268 LOG_SCOPE(
"enforce_adjoint_constraints_exactly()",
"DofMap");
3272 std::unique_ptr<NumericVector<Number>> v_built;
3280 i<v_built->last_local_index(); i++)
3281 v_built->set(i, v(i));
3283 v_global = v_built.
get();
3295 v_local = v_built.
get();
3305 libmesh_error_msg(
"ERROR: Unknown v.type() == " << v.
type());
3314 const AdjointDofConstraintValues::const_iterator
3318 nullptr : &adjoint_constraint_map_it->second;
3328 if (
const auto adjoint_constraint_it = constraint_map->find(constrained_dof);
3329 adjoint_constraint_it != constraint_map->end())
3333 for (
const auto & j : constraint_row)
3353 std::pair<Real, Real>
3364 Real max_absolute_error = 0., max_relative_error = 0.;
3368 libmesh_assert_equal_to (
this, &(system.
get_dof_map()));
3371 std::vector<dof_id_type> local_dof_indices;
3373 for (
const auto & elem :
mesh.active_local_element_ptr_range())
3376 std::vector<dof_id_type> raw_dof_indices = local_dof_indices;
3387 libmesh_assert_equal_to (C.
m(), raw_dof_indices.size());
3388 libmesh_assert_equal_to (C.
n(), local_dof_indices.size());
3399 DofConstraints::const_iterator
3406 DofConstraintValueMap::const_iterator rhsit =
3413 if (local_dof_indices[j] != global_dof)
3415 vec(local_dof_indices[j]);
3418 max_absolute_error = std::max(max_absolute_error,
3420 max_relative_error = std::max(max_relative_error,
3427 return std::pair<Real, Real>(max_absolute_error, max_relative_error);
3433 std::vector<dof_id_type> & elem_dofs,
3434 const bool called_recursively)
const 3436 LOG_SCOPE_IF(
"build_constraint_matrix()",
"DofMap", !called_recursively);
3439 typedef std::set<dof_id_type> RCSet;
3442 bool we_have_constraints =
false;
3449 for (
const auto & dof : elem_dofs)
3452 we_have_constraints =
true;
3455 DofConstraints::const_iterator
3465 for (
const auto & item : constraint_row)
3466 dof_set.insert (item.first);
3471 if (!we_have_constraints)
3474 for (
const auto & dof : elem_dofs)
3475 dof_set.erase (dof);
3483 if (!dof_set.empty() ||
3484 !called_recursively)
3486 const unsigned int old_size =
3487 cast_int<unsigned int>(elem_dofs.size());
3490 elem_dofs.insert(elem_dofs.end(),
3491 dof_set.begin(), dof_set.end());
3496 cast_int<unsigned int>(elem_dofs.size()));
3499 for (
unsigned int i=0; i != old_size; i++)
3503 DofConstraints::const_iterator
3513 for (
const auto & item : constraint_row)
3514 for (
unsigned int j=0,
3515 n_elem_dofs = cast_int<unsigned int>(elem_dofs.size());
3516 j != n_elem_dofs; j++)
3517 if (elem_dofs[j] == item.first)
3518 C(i,j) = item.second;
3532 if ((C.
n() == Cnew.
m()) &&
3533 (Cnew.
n() == elem_dofs.size()))
3536 libmesh_assert_equal_to (C.
n(), elem_dofs.size());
3544 std::vector<dof_id_type> & elem_dofs,
3546 const bool called_recursively)
const 3548 LOG_SCOPE_IF(
"build_constraint_matrix_and_vector()",
"DofMap", !called_recursively);
3551 typedef std::set<dof_id_type> RCSet;
3554 bool we_have_constraints =
false;
3561 for (
const auto & dof : elem_dofs)
3564 we_have_constraints =
true;
3567 DofConstraints::const_iterator
3577 for (
const auto & item : constraint_row)
3578 dof_set.insert (item.first);
3583 if (!we_have_constraints)
3586 for (
const auto & dof : elem_dofs)
3587 dof_set.erase (dof);
3595 if (!dof_set.empty() ||
3596 !called_recursively)
3603 rhs_values = &it->second;
3605 const unsigned int old_size =
3606 cast_int<unsigned int>(elem_dofs.size());
3609 elem_dofs.insert(elem_dofs.end(),
3610 dof_set.begin(), dof_set.end());
3615 cast_int<unsigned int>(elem_dofs.size()));
3619 for (
unsigned int i=0; i != old_size; i++)
3623 DofConstraints::const_iterator
3633 for (
const auto & item : constraint_row)
3634 for (
unsigned int j=0,
3635 n_elem_dofs = cast_int<unsigned int>(elem_dofs.size());
3636 j != n_elem_dofs; j++)
3637 if (elem_dofs[j] == item.first)
3638 C(i,j) = item.second;
3642 if (
const auto rhsit = rhs_values->find(elem_dofs[i]);
3643 rhsit != rhs_values->end())
3644 H(i) = rhsit->second;
3661 if ((C.
n() == Cnew.
m()) &&
3662 (Cnew.
n() == elem_dofs.size()))
3671 libmesh_assert_equal_to (C.
n(), elem_dofs.size());
3679 parallel_object_only();
3688 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS 3690 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS 3692 this->
comm().
max(has_constraints);
3693 if (!has_constraints)
3698 const unsigned int max_qoi_num =
3702 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS 3704 std::vector<Parallel::Request> packed_range_sends;
3718 std::map<processor_id_type, std::set<dof_id_type>> pushed_ids;
3720 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS 3721 std::map<processor_id_type, std::set<dof_id_type>> pushed_node_ids;
3724 const unsigned int sys_num = this->
sys_number();
3727 for (
auto & elem :
as_range(
mesh.active_not_local_elements_begin(),
3728 mesh.active_not_local_elements_end()))
3744 for (
unsigned int v=0; v !=
n_vars; ++v)
3746 const unsigned int n_comp = elem->
n_comp(sys_num,v);
3747 for (
unsigned int c=0; c != n_comp; ++c)
3749 const unsigned int id =
3757 for (
unsigned short n = 0; n !=
n_nodes; ++n)
3761 for (
unsigned int v=0; v !=
n_vars; ++v)
3763 const unsigned int n_comp = node.
n_comp(sys_num,v);
3764 for (
unsigned int c=0; c != n_comp; ++c)
3766 const unsigned int id =
3774 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS 3775 for (
unsigned short n = 0; n !=
n_nodes; ++n)
3783 std::map<processor_id_type, std::vector<dof_id_type>>
3784 pushed_id_vecs, received_id_vecs;
3785 for (
auto & p : pushed_ids)
3786 pushed_id_vecs[p.first].assign(p.second.begin(), p.second.end());
3788 std::map<processor_id_type, std::vector<std::vector<std::pair<dof_id_type,Real>>>>
3789 pushed_keys_vals, received_keys_vals;
3790 std::map<processor_id_type, std::vector<std::vector<Number>>> pushed_rhss, received_rhss;
3791 for (
auto & p : pushed_id_vecs)
3793 auto & keys_vals = pushed_keys_vals[p.first];
3794 keys_vals.reserve(p.second.size());
3796 auto & rhss = pushed_rhss[p.first];
3797 rhss.reserve(p.second.size());
3798 for (
auto & pushed_id : p.second)
3801 keys_vals.emplace_back(row.begin(), row.end());
3803 rhss.push_back(std::vector<Number>(max_qoi_num+1));
3804 std::vector<Number> & rhs = rhss.back();
3805 DofConstraintValueMap::const_iterator rhsit =
3810 for (
unsigned int q = 0; q != max_qoi_num; ++q)
3812 AdjointDofConstraintValues::const_iterator adjoint_map_it =
3819 adjoint_map_it->second;
3821 DofConstraintValueMap::const_iterator adj_rhsit =
3822 constraint_map.find(pushed_id);
3825 (adj_rhsit == constraint_map.end()) ?
3826 0 : adj_rhsit->second;
3831 auto ids_action_functor =
3832 [& received_id_vecs]
3834 const std::vector<dof_id_type> & data)
3836 received_id_vecs[pid] = data;
3839 Parallel::push_parallel_vector_data
3840 (this->
comm(), pushed_id_vecs, ids_action_functor);
3842 auto keys_vals_action_functor =
3843 [& received_keys_vals]
3845 const std::vector<std::vector<std::pair<dof_id_type,Real>>> & data)
3847 received_keys_vals[pid] = data;
3850 Parallel::push_parallel_vector_data
3851 (this->
comm(), pushed_keys_vals, keys_vals_action_functor);
3853 auto rhss_action_functor =
3856 const std::vector<std::vector<Number>> & data)
3858 received_rhss[pid] = data;
3861 Parallel::push_parallel_vector_data
3862 (this->
comm(), pushed_rhss, rhss_action_functor);
3867 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS 3868 std::map<processor_id_type, std::vector<dof_id_type>>
3869 pushed_node_id_vecs, received_node_id_vecs;
3870 for (
auto & p : pushed_node_ids)
3871 pushed_node_id_vecs[p.first].assign(p.second.begin(), p.second.end());
3873 std::map<processor_id_type, std::vector<std::vector<std::pair<dof_id_type,Real>>>>
3874 pushed_node_keys_vals, received_node_keys_vals;
3875 std::map<processor_id_type, std::vector<Point>> pushed_offsets, received_offsets;
3877 for (
auto & p : pushed_node_id_vecs)
3883 std::set<const Node *> nodes_requested;
3885 auto & node_keys_vals = pushed_node_keys_vals[pid];
3886 node_keys_vals.reserve(p.second.size());
3888 auto & offsets = pushed_offsets[pid];
3889 offsets.reserve(p.second.size());
3891 for (
auto & pushed_node_id : p.second)
3895 const std::size_t row_size = row.size();
3896 node_keys_vals.push_back
3897 (std::vector<std::pair<dof_id_type,Real>>());
3898 std::vector<std::pair<dof_id_type,Real>> & this_node_kv =
3899 node_keys_vals.back();
3900 this_node_kv.reserve(row_size);
3901 for (
const auto & j : row)
3903 this_node_kv.emplace_back(j.first->id(), j.second);
3908 if (j.first->processor_id() != pid && dist_mesh)
3909 nodes_requested.insert(j.first);
3923 (pid, &
mesh, nodes_requested.begin(), nodes_requested.end(),
3924 packed_range_sends.back(), range_tag);
3928 auto node_ids_action_functor =
3929 [& received_node_id_vecs]
3931 const std::vector<dof_id_type> & data)
3933 received_node_id_vecs[pid] = data;
3936 Parallel::push_parallel_vector_data
3937 (this->
comm(), pushed_node_id_vecs, node_ids_action_functor);
3939 auto node_keys_vals_action_functor =
3940 [& received_node_keys_vals]
3942 const std::vector<std::vector<std::pair<dof_id_type,Real>>> & data)
3944 received_node_keys_vals[pid] = data;
3947 Parallel::push_parallel_vector_data
3948 (this->
comm(), pushed_node_keys_vals,
3949 node_keys_vals_action_functor);
3951 auto node_offsets_action_functor =
3952 [& received_offsets]
3954 const std::vector<Point> & data)
3956 received_offsets[pid] = data;
3959 Parallel::push_parallel_vector_data
3960 (this->
comm(), pushed_offsets, node_offsets_action_functor);
3965 for (
auto & [pid, pushed_ids_to_me] : received_id_vecs)
3969 const auto & pushed_keys_vals_to_me = received_keys_vals.at(pid);
3970 const auto & pushed_rhss_to_me = received_rhss.at(pid);
3972 libmesh_assert_equal_to (pushed_ids_to_me.size(),
3973 pushed_keys_vals_to_me.size());
3974 libmesh_assert_equal_to (pushed_ids_to_me.size(),
3975 pushed_rhss_to_me.size());
3986 for (
auto & kv : pushed_keys_vals_to_me[i])
3988 libmesh_assert_less(kv.first, this->n_dofs());
3989 row[kv.first] = kv.second;
3992 const Number primal_rhs = pushed_rhss_to_me[i][max_qoi_num];
3997 if (primal_rhs !=
Number(0))
4002 for (
unsigned int q = 0; q != max_qoi_num; ++q)
4004 AdjointDofConstraintValues::iterator adjoint_map_it =
4007 const Number adj_rhs = pushed_rhss_to_me[i][q];
4018 adjoint_map_it->second;
4020 if (adj_rhs !=
Number(0))
4021 constraint_map[constrained] = adj_rhs;
4023 constraint_map.erase(constrained);
4029 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS 4031 for (
auto & [pid, pushed_node_ids_to_me] : received_node_id_vecs)
4038 (
Node**)
nullptr, range_tag);
4042 const auto & pushed_node_keys_vals_to_me = received_node_keys_vals.at(pid);
4043 const auto & pushed_offsets_to_me = received_offsets.at(pid);
4045 libmesh_assert_equal_to (pushed_node_ids_to_me.size(),
4046 pushed_node_keys_vals_to_me.size());
4047 libmesh_assert_equal_to (pushed_node_ids_to_me.size(),
4048 pushed_offsets_to_me.size());
4052 dof_id_type constrained_id = pushed_node_ids_to_me[i];
4060 for (
auto & kv : pushed_node_keys_vals_to_me[i])
4064 row[key_node] = kv.second;
4070 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS 4077 typedef std::set<dof_id_type> DoF_RCSet;
4078 DoF_RCSet unexpanded_dofs;
4081 unexpanded_dofs.insert(i.first);
4087 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS 4088 typedef std::set<const Node *> Node_RCSet;
4089 Node_RCSet unexpanded_nodes;
4092 unexpanded_nodes.insert(i.first);
4096 bool unexpanded_set_nonempty = !unexpanded_nodes.empty();
4097 this->
comm().
max(unexpanded_set_nonempty);
4099 while (unexpanded_set_nonempty)
4102 parallel_object_only();
4105 Node_RCSet node_request_set;
4108 std::map<processor_id_type, std::vector<dof_id_type>>
4112 std::map<processor_id_type, dof_id_type> node_ids_on_proc;
4115 for (
const auto & i : unexpanded_nodes)
4118 for (
const auto & j : row)
4120 const Node *
const node = j.first;
4126 !_node_constraints.count(node))
4127 node_request_set.insert(node);
4133 unexpanded_nodes.clear();
4136 for (
const auto & node : node_request_set)
4139 libmesh_assert_less (node->processor_id(), this->
n_processors());
4140 node_ids_on_proc[node->processor_id()]++;
4143 for (
auto pair : node_ids_on_proc)
4144 requested_node_ids[pair.first].reserve(pair.second);
4147 for (
const auto & node : node_request_set)
4148 requested_node_ids[node->processor_id()].push_back(node->id());
4150 typedef std::vector<std::pair<dof_id_type, Real>> row_datum;
4152 auto node_row_gather_functor =
4156 & packed_range_sends,
4159 const std::vector<dof_id_type> & ids,
4160 std::vector<row_datum> & data)
4164 std::set<const Node *> nodes_requested;
4167 const std::size_t query_size = ids.size();
4169 data.resize(query_size);
4170 for (std::size_t i=0; i != query_size; ++i)
4177 std::size_t row_size = row.size();
4178 data[i].reserve(row_size);
4179 for (
const auto & j : row)
4181 const Node * node = j.first;
4182 data[i].emplace_back(node->
id(), j.second);
4188 nodes_requested.insert(node);
4214 (pid, &
mesh, nodes_requested.begin(), nodes_requested.end(),
4215 packed_range_sends.back(), range_tag);
4219 typedef Point node_rhs_datum;
4221 auto node_rhs_gather_functor =
4225 const std::vector<dof_id_type> & ids,
4226 std::vector<node_rhs_datum> & data)
4229 const std::size_t query_size = ids.size();
4231 data.resize(query_size);
4232 for (std::size_t i=0; i != query_size; ++i)
4239 data[i](0) = std::numeric_limits<Real>::quiet_NaN();
4243 auto node_row_action_functor =
4250 const std::vector<dof_id_type> & ids,
4251 const std::vector<row_datum> & data)
4258 (
Node**)
nullptr, range_tag);
4261 const std::size_t query_size = ids.size();
4263 for (std::size_t i=0; i != query_size; ++i)
4269 if (data[i].empty())
4280 for (
auto & pair : data[i])
4282 const Node * key_node =
4285 row[key_node] = pair.second;
4289 unexpanded_nodes.insert(constrained_node);
4294 auto node_rhs_action_functor =
4298 const std::vector<dof_id_type> & ids,
4299 const std::vector<node_rhs_datum> & data)
4302 const std::size_t query_size = ids.size();
4304 for (std::size_t i=0; i != query_size; ++i)
4317 row_datum * node_row_ex =
nullptr;
4318 Parallel::pull_parallel_vector_data
4319 (this->
comm(), requested_node_ids, node_row_gather_functor,
4320 node_row_action_functor, node_row_ex);
4323 node_rhs_datum * node_rhs_ex =
nullptr;
4324 Parallel::pull_parallel_vector_data
4325 (this->
comm(), requested_node_ids, node_rhs_gather_functor,
4326 node_rhs_action_functor, node_rhs_ex);
4331 unexpanded_set_nonempty = !unexpanded_nodes.empty();
4332 this->
comm().
max(unexpanded_set_nonempty);
4334 Parallel::wait(packed_range_sends);
4335 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS 4359 const unsigned int max_qoi_num =
4364 typedef std::set<dof_id_type> RCSet;
4365 RCSet unexpanded_set;
4368 unexpanded_set.insert(i.first);
4370 while (!unexpanded_set.empty())
4371 for (RCSet::iterator i = unexpanded_set.begin();
4372 i != unexpanded_set.end(); )
4375 DofConstraints::iterator
4382 DofConstraintValueMap::iterator rhsit =
4388 std::vector<DofConstraintValueMap::iterator> adjoint_rhs_iterators;
4389 adjoint_rhs_iterators.resize(max_qoi_num);
4392 std::vector<Number> adjoint_constraint_rhs(max_qoi_num, 0.0);
4397 const std::size_t q = adjoint_map.first;
4398 adjoint_rhs_iterators[q] = adjoint_map.second.find(*i);
4400 adjoint_constraint_rhs[q] =
4401 (adjoint_rhs_iterators[q] == adjoint_map.second.end()) ?
4402 0 : adjoint_rhs_iterators[q]->second;
4405 std::vector<dof_id_type> constraints_to_expand;
4407 for (
const auto & item : constraint_row)
4408 if (item.first != *i && this->is_constrained_dof(item.first))
4410 unexpanded_set.insert(item.first);
4411 constraints_to_expand.push_back(item.first);
4414 for (
const auto & expandable : constraints_to_expand)
4416 const Real this_coef = constraint_row[expandable];
4418 DofConstraints::const_iterator
4425 for (
const auto & item : subconstraint_row)
4429 constraint_row[item.first] += item.second * this_coef;
4434 constraint_rhs += subrhsit->second * this_coef;
4439 if (
auto adjoint_subrhsit = adjoint_map.second.find(expandable);
4440 adjoint_subrhsit != adjoint_map.second.end())
4441 adjoint_constraint_rhs[adjoint_map.first] += adjoint_subrhsit->second * this_coef;
4444 constraint_row.erase(expandable);
4449 if (constraint_rhs !=
Number(0))
4456 if (constraint_rhs !=
Number(0))
4457 rhsit->second = constraint_rhs;
4465 const std::size_t q = adjoint_map.first;
4467 if(adjoint_rhs_iterators[q] == adjoint_map.second.end())
4469 if (adjoint_constraint_rhs[q] !=
Number(0))
4470 (adjoint_map.second)[*i] = adjoint_constraint_rhs[q];
4472 adjoint_map.second.erase(*i);
4476 if (adjoint_constraint_rhs[q] !=
Number(0))
4477 adjoint_rhs_iterators[q]->second = adjoint_constraint_rhs[q];
4479 adjoint_map.second.erase(adjoint_rhs_iterators[q]);
4483 if (constraints_to_expand.empty())
4484 i = unexpanded_set.erase(i);
4501 #ifdef LIBMESH_ENABLE_CONSTRAINTS 4511 typedef std::set<dof_id_type> RCSet;
4512 RCSet unexpanded_set;
4518 for (
const auto & i : dof_constraints_copy)
4519 unexpanded_set.insert(i.first);
4521 while (!unexpanded_set.empty())
4522 for (RCSet::iterator i = unexpanded_set.begin();
4523 i != unexpanded_set.end(); )
4526 DofConstraints::iterator
4527 pos = dof_constraints_copy.find(*i);
4539 std::vector<dof_id_type> constraints_to_expand;
4541 for (
const auto & item : constraint_row)
4542 if (item.first != *i && this->is_constrained_dof(item.first))
4544 unexpanded_set.insert(item.first);
4545 constraints_to_expand.push_back(item.first);
4548 for (
const auto & expandable : constraints_to_expand)
4550 const Real this_coef = constraint_row[expandable];
4552 DofConstraints::const_iterator
4553 subpos = dof_constraints_copy.find(expandable);
4559 for (
const auto & item : subconstraint_row)
4561 libmesh_error_msg_if(item.first == expandable,
"Constraint loop detected");
4563 constraint_row[item.first] += item.second * this_coef;
4572 constraint_row.erase(expandable);
4591 if (constraints_to_expand.empty())
4592 i = unexpanded_set.erase(i);
4613 parallel_object_only();
4622 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS 4624 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS 4626 this->
comm().
max(has_constraints);
4627 if (!has_constraints)
4634 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS 4635 std::map<processor_id_type, std::set<dof_id_type>> pushed_node_ids;
4636 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS 4638 std::map<processor_id_type, std::set<dof_id_type>> pushed_ids;
4644 while (constrained >=
_end_df[constrained_proc_id])
4645 constrained_proc_id++;
4650 for (
auto & j : row)
4655 while (constraining >=
_end_df[constraining_proc_id])
4656 constraining_proc_id++;
4659 constraining_proc_id != constrained_proc_id)
4660 pushed_ids[constraining_proc_id].insert(constrained);
4667 std::vector<std::vector<std::pair<dof_id_type, Real>>>>
4668 pushed_keys_vals, pushed_keys_vals_to_me;
4670 std::map<processor_id_type, std::vector<std::pair<dof_id_type, Number>>>
4671 pushed_ids_rhss, pushed_ids_rhss_to_me;
4680 for (
const auto & [pid, pid_ids] : pushed_ids)
4682 const std::size_t ids_size = pid_ids.size();
4683 std::vector<std::vector<std::pair<dof_id_type, Real>>> &
4684 keys_vals = pushed_keys_vals[pid];
4685 std::vector<std::pair<dof_id_type,Number>> &
4686 ids_rhss = pushed_ids_rhss[pid];
4687 keys_vals.resize(ids_size);
4688 ids_rhss.resize(ids_size);
4691 std::set<dof_id_type>::const_iterator it;
4692 for (push_i = 0, it = pid_ids.begin();
4693 it != pid_ids.end(); ++push_i, ++it)
4697 keys_vals[push_i].assign(row.begin(), row.end());
4699 DofConstraintValueMap::const_iterator rhsit =
4701 ids_rhss[push_i].first = constrained;
4702 ids_rhss[push_i].second =
4711 auto ids_rhss_action_functor =
4712 [& pushed_ids_rhss_to_me]
4714 const std::vector<std::pair<dof_id_type, Number>> & data)
4716 pushed_ids_rhss_to_me[pid] = data;
4719 auto keys_vals_action_functor =
4720 [& pushed_keys_vals_to_me]
4722 const std::vector<std::vector<std::pair<dof_id_type, Real>>> & data)
4724 pushed_keys_vals_to_me[pid] = data;
4727 Parallel::push_parallel_vector_data
4728 (this->
comm(), pushed_ids_rhss, ids_rhss_action_functor);
4729 Parallel::push_parallel_vector_data
4730 (this->
comm(), pushed_keys_vals, keys_vals_action_functor);
4733 auto receive_dof_constraints =
4735 & pushed_ids_rhss_to_me,
4736 & pushed_keys_vals_to_me]
4739 for (
const auto & [pid, ids_rhss] : pushed_ids_rhss_to_me)
4741 const auto & keys_vals = pushed_keys_vals_to_me[pid];
4743 libmesh_assert_equal_to
4744 (ids_rhss.size(), keys_vals.size());
4756 for (
auto & key_val : keys_vals[i])
4758 libmesh_assert_less(key_val.first, this->n_dofs());
4759 row[key_val.first] = key_val.second;
4761 if (ids_rhss[i].second !=
Number(0))
4771 receive_dof_constraints();
4773 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS 4777 const Node * constrained = i.first;
4783 for (
auto & j : row)
4785 const Node * constraining = j.first;
4789 pushed_node_ids[constraining->
processor_id()].insert(constrained->
id());
4795 std::vector<std::vector<std::pair<dof_id_type,Real>>>>
4796 pushed_node_keys_vals, pushed_node_keys_vals_to_me;
4797 std::map<processor_id_type, std::vector<std::pair<dof_id_type, Point>>>
4798 pushed_node_ids_offsets, pushed_node_ids_offsets_to_me;
4799 std::map<processor_id_type, std::vector<const Node *>> pushed_node_vecs;
4801 for (
const auto & [pid, pid_ids]: pushed_node_ids)
4803 const std::size_t ids_size = pid_ids.size();
4804 std::vector<std::vector<std::pair<dof_id_type,Real>>> &
4805 keys_vals = pushed_node_keys_vals[pid];
4806 std::vector<std::pair<dof_id_type, Point>> &
4807 ids_offsets = pushed_node_ids_offsets[pid];
4808 keys_vals.resize(ids_size);
4809 ids_offsets.resize(ids_size);
4810 std::set<Node *> nodes;
4813 std::set<dof_id_type>::const_iterator it;
4814 for (push_i = 0, it = pid_ids.begin();
4815 it != pid_ids.end(); ++push_i, ++it)
4820 nodes.insert(constrained);
4823 std::size_t row_size = row.size();
4824 keys_vals[push_i].reserve(row_size);
4825 for (
const auto & j : row)
4827 Node * constraining =
const_cast<Node *
>(j.first);
4829 keys_vals[push_i].emplace_back(constraining->
id(), j.second);
4832 nodes.insert(constraining);
4835 ids_offsets[push_i].first = *it;
4841 auto & pid_nodes = pushed_node_vecs[pid];
4842 pid_nodes.assign(nodes.begin(), nodes.end());
4846 auto node_ids_offsets_action_functor =
4847 [& pushed_node_ids_offsets_to_me]
4849 const std::vector<std::pair<dof_id_type, Point>> & data)
4851 pushed_node_ids_offsets_to_me[pid] = data;
4854 auto node_keys_vals_action_functor =
4855 [& pushed_node_keys_vals_to_me]
4857 const std::vector<std::vector<std::pair<dof_id_type, Real>>> & data)
4859 pushed_node_keys_vals_to_me[pid] = data;
4863 Parallel::push_parallel_vector_data
4864 (this->
comm(), pushed_node_ids_offsets, node_ids_offsets_action_functor);
4865 Parallel::push_parallel_vector_data
4866 (this->
comm(), pushed_node_keys_vals, node_keys_vals_action_functor);
4872 auto null_node_functor = [](
processor_id_type,
const std::vector<const Node *> &){};
4875 Parallel::push_parallel_packed_range
4876 (this->
comm(), pushed_node_vecs, &
mesh, null_node_functor);
4878 for (
const auto & [pid, ids_offsets] : pushed_node_ids_offsets_to_me)
4880 const auto & keys_vals = pushed_node_keys_vals_to_me[pid];
4882 libmesh_assert_equal_to
4883 (ids_offsets.size(), keys_vals.size());
4888 dof_id_type constrained_id = ids_offsets[i].first;
4896 for (
auto & key_val : keys_vals[i])
4899 row[key_node] = key_val.second;
4902 ids_offsets[i].second;
4906 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS 4919 typedef std::map<dof_id_type, std::set<dof_id_type>> DofConstrainsMap;
4920 DofConstrainsMap dof_id_constrains;
4924 for (
const auto & j : row)
4929 while (constraining >=
_end_df[constraining_proc_id])
4930 constraining_proc_id++;
4933 dof_id_constrains[constraining].insert(constrained);
4941 for (
const auto & elem :
as_range(
mesh.active_not_local_elements_begin(),
4942 mesh.active_not_local_elements_end()))
4944 std::vector<dof_id_type> my_dof_indices;
4947 for (
const auto & dof : my_dof_indices)
4949 if (
auto dcmi = dof_id_constrains.find(dof);
4950 dcmi != dof_id_constrains.end())
4952 for (
const auto & constrained : dcmi->second)
4955 while (constrained >=
_end_df[the_constrained_proc_id])
4956 the_constrained_proc_id++;
4959 if (elemproc != the_constrained_proc_id)
4960 pushed_ids[elemproc].insert(constrained);
4966 pushed_ids_rhss.clear();
4967 pushed_ids_rhss_to_me.clear();
4968 pushed_keys_vals.clear();
4969 pushed_keys_vals_to_me.clear();
4974 Parallel::push_parallel_vector_data
4975 (this->
comm(), pushed_ids_rhss, ids_rhss_action_functor);
4976 Parallel::push_parallel_vector_data
4977 (this->
comm(), pushed_keys_vals, keys_vals_action_functor);
4979 receive_dof_constraints();
4991 (elements_to_couple,
4992 temporary_coupling_matrices,
4995 mesh.active_local_elements_begin(),
4996 mesh.active_local_elements_end(),
5000 std::set<dof_id_type> requested_dofs;
5002 for (
const auto & pr : elements_to_couple)
5004 const Elem * elem = pr.first;
5007 std::vector<dof_id_type> element_dofs;
5010 for (
auto dof : element_dofs)
5011 requested_dofs.insert(dof);
5019 std::set<dof_id_type> & unexpanded_dofs,
5022 typedef std::set<dof_id_type> DoF_RCSet;
5026 const unsigned int max_qoi_num =
5032 bool unexpanded_set_nonempty = !unexpanded_dofs.empty();
5033 this->
comm().
max(unexpanded_set_nonempty);
5035 while (unexpanded_set_nonempty)
5038 parallel_object_only();
5041 DoF_RCSet dof_request_set;
5044 std::map<processor_id_type, std::vector<dof_id_type>>
5048 std::map<processor_id_type, dof_id_type>
5052 for (
const auto & unexpanded_dof : unexpanded_dofs)
5061 dof_request_set.insert(unexpanded_dof);
5069 for (
const auto & j : row)
5077 dof_request_set.insert(constraining_dof);
5083 unexpanded_dofs.clear();
5087 for (
const auto & i : dof_request_set)
5091 dof_ids_on_proc[proc_id]++;
5094 for (
auto & pair : dof_ids_on_proc)
5096 requested_dof_ids[pair.first].reserve(pair.second);
5101 for (
const auto & i : dof_request_set)
5105 requested_dof_ids[proc_id].push_back(i);
5108 typedef std::vector<std::pair<dof_id_type, Real>> row_datum;
5110 typedef std::vector<Number> rhss_datum;
5112 auto row_gather_functor =
5115 const std::vector<dof_id_type> & ids,
5116 std::vector<row_datum> & data)
5119 const std::size_t query_size = ids.size();
5121 data.resize(query_size);
5122 for (std::size_t i=0; i != query_size; ++i)
5128 std::size_t row_size = row.size();
5129 data[i].reserve(row_size);
5130 for (
const auto & j : row)
5132 data[i].push_back(j);
5160 auto rhss_gather_functor =
5164 const std::vector<dof_id_type> & ids,
5165 std::vector<rhss_datum> & data)
5168 const std::size_t query_size = ids.size();
5170 data.resize(query_size);
5171 for (std::size_t i=0; i != query_size; ++i)
5177 DofConstraintValueMap::const_iterator rhsit =
5183 for (
unsigned int q = 0; q != max_qoi_num; ++q)
5185 AdjointDofConstraintValues::const_iterator adjoint_map_it =
5190 data[i].push_back(0);
5195 adjoint_map_it->second;
5197 DofConstraintValueMap::const_iterator adj_rhsit =
5198 constraint_map.find(constrained);
5200 ((adj_rhsit == constraint_map.end()) ?
5201 0 : adj_rhsit->second);
5207 auto row_action_functor =
5211 const std::vector<dof_id_type> & ids,
5212 const std::vector<row_datum> & data)
5215 const std::size_t query_size = ids.size();
5217 for (std::size_t i=0; i != query_size; ++i)
5223 if (data[i].empty())
5232 for (
auto & pair : data[i])
5234 libmesh_assert_less(pair.first, this->n_dofs());
5235 row[pair.first] = pair.second;
5239 unexpanded_dofs.insert(constrained);
5244 auto rhss_action_functor =
5248 const std::vector<dof_id_type> & ids,
5249 const std::vector<rhss_datum> & data)
5252 const std::size_t query_size = ids.size();
5254 for (std::size_t i=0; i != query_size; ++i)
5256 if (!data[i].empty())
5259 if (data[i][0] !=
Number(0))
5264 for (
unsigned int q = 0; q != max_qoi_num; ++q)
5266 AdjointDofConstraintValues::iterator adjoint_map_it =
5270 data[i][q+1] ==
Number(0))
5278 adjoint_map_it->second;
5280 if (data[i][q+1] !=
Number(0))
5281 constraint_map[constrained] =
5284 constraint_map.erase(constrained);
5292 row_datum * row_ex =
nullptr;
5293 Parallel::pull_parallel_vector_data
5294 (this->
comm(), requested_dof_ids, row_gather_functor,
5295 row_action_functor, row_ex);
5298 rhss_datum * rhs_ex =
nullptr;
5299 Parallel::pull_parallel_vector_data
5300 (this->
comm(), requested_dof_ids, rhss_gather_functor,
5301 rhss_action_functor, rhs_ex);
5305 unexpanded_set_nonempty = !unexpanded_dofs.empty();
5306 this->
comm().
max(unexpanded_set_nonempty);
5313 parallel_object_only();
5322 this->
comm().
max(has_constraints);
5323 if (!has_constraints)
5332 for (
const auto & j : constraint_row)
5347 #endif // LIBMESH_ENABLE_CONSTRAINTS 5350 #ifdef LIBMESH_ENABLE_AMR 5360 libmesh_assert_greater (elem->
p_level(), p);
5361 libmesh_assert_less (s, elem->
n_sides());
5363 const unsigned int sys_num = this->
sys_number();
5367 for (
unsigned int n = 0; n !=
n_nodes; ++n)
5371 const unsigned int low_nc =
5373 const unsigned int high_nc =
5385 for (
unsigned int i = low_nc; i != high_nc; ++i)
5393 const unsigned int total_dofs = node.
n_comp(sys_num, var);
5394 libmesh_assert_greater_equal (total_dofs, high_nc);
5397 for (
unsigned int j = low_nc; j != high_nc; ++j)
5399 const unsigned int i = total_dofs - j - 1;
5407 #endif // LIBMESH_ENABLE_AMR 5410 #ifdef LIBMESH_ENABLE_DIRICHLET 5418 unsigned int qoi_index)
5420 unsigned int old_size = cast_int<unsigned int>
5422 for (
unsigned int i = old_size; i <= qoi_index; ++i)
5427 (std::make_unique<DirichletBoundary>(dirichlet_boundary));
5451 unsigned int old_size = cast_int<unsigned int>
5453 for (
unsigned int i = old_size; i <= q; ++i)
5463 auto lam = [&boundary_to_remove](
const auto & bdy)
5464 {
return bdy->b == boundary_to_remove.
b && bdy->variables == boundary_to_remove.
variables;};
5475 unsigned int qoi_index)
5480 auto lam = [&boundary_to_remove](
const auto & bdy)
5481 {
return bdy->b == boundary_to_remove.
b && bdy->variables == boundary_to_remove.
variables;};
5496 const std::set<boundary_id_type>& mesh_side_bcids =
5498 const std::set<boundary_id_type>& mesh_edge_bcids =
5500 const std::set<boundary_id_type>& mesh_node_bcids =
5502 const std::set<boundary_id_type>& dbc_bcids = boundary.
b;
5507 for (
const auto & bc_id : dbc_bcids)
5512 bool found_bcid = (mesh_side_bcids.find(bc_id) != mesh_side_bcids.end() ||
5513 mesh_edge_bcids.find(bc_id) != mesh_edge_bcids.end() ||
5514 mesh_node_bcids.find(bc_id) != mesh_node_bcids.end());
5521 libmesh_error_msg_if(!found_bcid,
5522 "Could not find Dirichlet boundary id " << bc_id <<
" in mesh!");
5526 #endif // LIBMESH_ENABLE_DIRICHLET 5529 #ifdef LIBMESH_ENABLE_PERIODIC 5536 if (!existing_boundary)
5546 existing_boundary->
merge(periodic_boundary);
5551 inverse_boundary->
merge(periodic_boundary);
void remove_adjoint_dirichlet_boundary(const DirichletBoundary &dirichlet_boundary, unsigned int q)
Removes from the system the specified Dirichlet boundary for the adjoint equation defined by Quantity...
std::unique_ptr< FEMFunctionBase< Gradient > > g_fem
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
std::string get_local_constraints(bool print_nonlocal=false) const
Gets a string reporting all DoF and Node constraints local to this processor.
FEFamily family
The type of finite element.
void constrain_element_dyad_matrix(DenseVector< Number > &v, DenseVector< Number > &w, std::vector< dof_id_type > &row_dofs, bool asymmetric_constraint_rows=true) const
Constrains a dyadic element matrix B = v w'.
bool is_constrained_node(const Node *node) const
void add_adjoint_dirichlet_boundary(const DirichletBoundary &dirichlet_boundary, unsigned int q)
Adds a copy of the specified Dirichlet boundary to the system, corresponding to the adjoint problem d...
void check_dirichlet_bcid_consistency(const MeshBase &mesh, const DirichletBoundary &boundary) const
Check that all the ids in dirichlet_bcids are actually present in the mesh.
dof_id_type dof_number(const unsigned int s, const unsigned int var, const unsigned int comp) const
virtual Output component(unsigned int i, const Point &p, Real time=0.)
constraint_rows_type & get_constraint_rows()
Constraint rows accessors.
A Node is like a Point, but with more information.
std::set< GhostingFunctor * >::const_iterator coupling_functors_begin() const
Beginning of range of coupling functors.
void build_constraint_matrix_and_vector(DenseMatrix< Number > &C, DenseVector< Number > &H, std::vector< dof_id_type > &elem_dofs, int qoi_index=-1, const bool called_recursively=false) const
Build the constraint matrix C and the forcing vector H associated with the element degree of freedom ...
void get_side_fe(unsigned int var, FEGenericBase< OutputShape > *&fe) const
Accessor for edge/face (2D/3D) finite element object for variable var for the largest dimension in th...
virtual void zero() override final
Set every element in the vector to 0.
virtual void pre_fe_reinit(const System &, const Elem *e)
Reinitializes local data vectors/matrices on the current geometric element.
unsigned int n_comp(const unsigned int s, const unsigned int var) const
std::unique_ptr< PointLocatorBase > sub_point_locator() const
std::unique_ptr< FunctionBase< Number > > f
void scatter_constraints(MeshBase &)
Sends constraint equations to constraining processors.
void add_periodic_boundary(const PeriodicBoundaryBase &periodic_boundary)
Adds a copy of the specified periodic boundary to the system.
static void compute_node_constraints(NodeConstraints &constraints, const Elem *elem)
Computes the nodal constraint contributions (for non-conforming adapted meshes), using Lagrange geome...
bool _error_on_constraint_loop
This flag indicates whether or not we do an opt-mode check for the presence of constraint loops...
void add_adjoint_constraint_row(const unsigned int qoi_index, const dof_id_type dof_number, const DofConstraintRow &constraint_row, const Number constraint_rhs, const bool forbid_constraint_overwrite)
Adds a copy of the user-defined row to the constraint matrix, using an inhomogeneous right-hand-side ...
void get_edge_fe(unsigned int var, FEGenericBase< OutputShape > *&fe) const
Accessor for edge (3D only!) finite element object for variable var.
void process_mesh_constraint_rows(const MeshBase &mesh)
Adds any spline constraints from the Mesh to our DoF constraints.
The IntRange templated class is intended to make it easy to loop over integers which are indices of a...
void constrain_element_matrix_and_vector(DenseMatrix< Number > &matrix, DenseVector< Number > &rhs, std::vector< dof_id_type > &elem_dofs, bool asymmetric_constraint_rows=true) const
Constrains the element matrix and vector.
virtual void get(const std::vector< numeric_index_type > &index, T *values) const
Access multiple components at once.
static std::unique_ptr< LinearSolver< T > > build(const libMesh::Parallel::Communicator &comm_in, const SolverPackage solver_package=libMesh::default_solver_package())
Builds a LinearSolver using the linear solver package specified by solver_package.
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
virtual void zero() override final
Sets all elements of the matrix to 0 and resets any decomposition flag which may have been previously...
void parallel_for(const Range &range, const Body &body)
Execute the provided function object in parallel on the specified range.
dof_id_type n_dofs() const
static constexpr Real TOLERANCE
Real jacobian_tolerance
Defaults to zero, but can be set to a custom small negative value to try and avoid spurious zero (or ...
We're using a class instead of a typedef to allow forward declarations and future flexibility...
virtual bool is_edge_on_side(const unsigned int e, const unsigned int s) const =0
MessageTag get_unique_tag(int tagvalue=MessageTag::invalid_tag) const
dof_id_type n_local_constrained_dofs() const
virtual numeric_index_type size() const =0
void resize(const unsigned int n)
Resize the vector.
void remove_dirichlet_boundary(const DirichletBoundary &dirichlet_boundary)
Removes the specified Dirichlet boundary from the system.
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, const bool add_p_level=true)
Fills the vector di with the local degree of freedom indices associated with edge e of element elem A...
std::vector< dof_id_type > _send_list
A list containing all the global DOF indices that affect the solution on my processor.
const FEType & variable_type(const unsigned int c) const
std::unique_ptr< SparsityPattern::Build > build_sparsity(const MeshBase &mesh, bool calculate_constrained=false, bool use_condensed_system=false) const
Builds a sparsity pattern for matrices using the current degree-of-freedom numbering and coupling...
This is the base class from which all geometric element types are derived.
boundary_id_type pairedboundary
void gather_constraints(MeshBase &mesh, std::set< dof_id_type > &unexpanded_dofs, bool look_for_constrainees)
Helper function for querying about constraint equations on other processors.
static FEFieldType field_type(const FEType &fe_type)
virtual bool is_node_on_side(const unsigned int n, const unsigned int s) const =0
std::unique_ptr< DirichletBoundaries > _dirichlet_boundaries
Data structure containing Dirichlet functions.
NumericVector< Number > * rhs
The system matrix.
This class allows one to associate Dirichlet boundary values with a given set of mesh boundary ids an...
const Parallel::Communicator & comm() const
The Node constraint storage format.
unsigned int p_level() const
std::vector< unsigned int > variables
void set_algebraic_type(const AlgebraicType atype)
Setting which determines whether to initialize algebraic structures (elem_*) on each element and set ...
void shellface_boundary_ids(const Elem *const elem, const unsigned short int shellface, std::vector< boundary_id_type > &vec_to_fill) const
void boundary_ids(const Node *node, std::vector< boundary_id_type > &vec_to_fill) const
Fills a user-provided std::vector with the boundary ids associated with Node node.
std::map< const Elem *, const CouplingMatrix *, CompareDofObjectsByPIDAndThenID > map_type
What elements do we care about and what variables do we care about on each element?
This class defines a vector in LIBMESH_DIM dimensional Real or Complex space.
The StoredRange class defines a contiguous, divisible set of objects.
The libMesh namespace provides an interface to certain functionality in the library.
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
const MeshBase & get_mesh() const
void send_packed_range(const unsigned int dest_processor_id, const Context *context, Iter range_begin, const Iter range_end, const MessageTag &tag=no_tag, std::size_t approx_buffer_size=1000000) const
Real distance(const Point &p)
Number exact_value(const Point &p, const Parameters ¶meters, const std::string &, const std::string &)
static std::unique_ptr< SparseMatrix< T > > build(const Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package(), const MatrixBuildType matrix_build_type=MatrixBuildType::AUTOMATIC)
Builds a SparseMatrix<T> using the linear solver package specified by solver_package.
void iota(ForwardIter first, ForwardIter last, T value)
Utility::iota was created back when std::iota was just an SGI STL extension.
unsigned int sys_number() const
uint8_t processor_id_type
This is the MeshBase class.
const Variable & variable(const unsigned int c) const override
virtual void set(const numeric_index_type i, const numeric_index_type j, const T value)=0
Set the element (i,j) to value.
void enforce_adjoint_constraints_exactly(NumericVector< Number > &v, unsigned int q) const
Heterogeneously constrains the numeric vector v, which represents an adjoint solution defined on the ...
AdjointDofConstraintValues _adjoint_constraint_values
void vector_mult_transpose(DenseVector< T > &dest, const DenseVector< T > &arg) const
Performs the matrix-vector multiplication, dest := (*this)^T * arg.
dof_id_type end_dof() const
void check_for_constraint_loops()
static void compute_constraints(DofConstraints &constraints, DofMap &dof_map, const unsigned int variable_number, const Elem *elem)
Computes the constraint matrix contributions (for non-conforming adapted meshes) corresponding to var...
void merge(const PeriodicBoundaryBase &pb)
const std::set< boundary_id_type > & get_node_boundary_ids() const
void enforce_constraints_on_residual(const NonlinearImplicitSystem &system, NumericVector< Number > *rhs, NumericVector< Number > const *solution, bool homogeneous=true) const
This class handles the numbering of degrees of freedom on a mesh.
std::unique_ptr< FEMFunctionBase< Number > > f_fem
virtual bool is_node_on_edge(const unsigned int n, const unsigned int e) const =0
const std::vector< std::vector< OutputGradient > > & get_dphi() const
unsigned int first_scalar_number() const
processor_id_type n_processors() const
virtual bool is_serial() const
void set_jacobian_tolerance(Real tol)
Set the Jacobian tolerance used for determining when the mapping fails.
boundary_id_type myboundary
The boundary ID of this boundary and its counterpart.
const dof_id_type n_nodes
dof_id_type first_dof() const
Status receive(const unsigned int dest_processor_id, T &buf, const MessageTag &tag=any_tag) const
This class defines the notion of a variable in the system.
static void merge_ghost_functor_outputs(GhostingFunctor::map_type &elements_to_ghost, CouplingMatricesSet &temporary_coupling_matrices, const std::set< GhostingFunctor *>::iterator &gf_begin, const std::set< GhostingFunctor *>::iterator &gf_end, const MeshBase::const_element_iterator &elems_begin, const MeshBase::const_element_iterator &elems_end, processor_id_type p)
void vector_mult_add(DenseVector< T > &dest, const T factor, const DenseVector< T > &arg) const
Performs the scaled matrix-vector multiplication, dest += factor * (*this) * arg. ...
bool has_adjoint_dirichlet_boundaries(unsigned int q) const
const Node & node_ref(const unsigned int i) const
void min(const T &r, T &o, Request &req) const
void add_constraint_row(const dof_id_type dof_number, const DofConstraintRow &constraint_row, const Number constraint_rhs, const bool forbid_constraint_overwrite)
Adds a copy of the user-defined row to the constraint matrix, using an inhomogeneous right-hand-side ...
const System & get_system() const
Accessor for associated system.
virtual unsigned int n_nodes() const =0
std::set< std::unique_ptr< CouplingMatrix >, Utility::CompareUnderlying > CouplingMatricesSet
bool is_sorted(const std::vector< KeyType > &v)
void heterogeneously_constrain_element_vector(const DenseMatrix< Number > &matrix, DenseVector< Number > &rhs, std::vector< dof_id_type > &elem_dofs, bool asymmetric_constraint_rows=true, int qoi_index=-1) const
Constrains the element vector.
void receive_packed_range(const unsigned int dest_processor_id, Context *context, OutputIter out, const T *output_type, const MessageTag &tag=any_tag) const
unsigned int n_vars(const unsigned int s, const unsigned int vg) const
Manages consistently variables, degrees of freedom, and coefficient vectors.
unsigned int n_variables() const override
SimpleRange< IndexType > as_range(const std::pair< IndexType, IndexType > &p)
Helper function that allows us to treat a homogenous pair as a range.
virtual void init_context(const FEMContext &)
Prepares a context object for use.
std::unique_ptr< NumericVector< Number > > solution
Data structure to hold solution values.
We're using a class instead of a typedef to allow forward declarations and future flexibility...
void print_dof_constraints(std::ostream &os=libMesh::out, bool print_nonlocal=false) const
Prints (from processor 0) all DoF and Node constraints.
The BoundaryInfo class contains information relevant to boundary conditions including storing faces...
bool is_constrained_dof(const dof_id_type dof) const
Manages consistently variables, degrees of freedom, coefficient vectors, matrices and non-linear solv...
static const dof_id_type invalid_id
An invalid id to distinguish an uninitialized DofObject.
virtual_for_inffe const std::vector< Real > & get_JxW() const
void allgather_recursive_constraints(MeshBase &)
Gathers constraint equation dependencies from other processors.
This class provides all data required for a physics package (e.g.
virtual void reinit(const Elem *elem, const std::vector< Point > *const pts=nullptr, const std::vector< Real > *const weights=nullptr)=0
This is at the core of this class.
const std::vector< dof_id_type > & get_dof_indices() const
Accessor for element dof indices.
DofConstraints _dof_constraints
Data structure containing DOF constraints.
unsigned int n_points() const
virtual_for_inffe const std::vector< Point > & get_xyz() const
bool active_on_subdomain(subdomain_id_type sid) const
virtual unsigned int n_edges() const =0
virtual std::unique_ptr< PeriodicBoundaryBase > clone(TransformationType t=FORWARD) const =0
If we want the DofMap to be able to make copies of references and store them in the underlying map...
static unsigned int n_vec_dim(const MeshBase &mesh, const FEType &fe_type)
virtual void close()=0
Calls the NumericVector's internal assembly routines, ensuring that the values are consistent across ...
std::pair< Real, Real > max_constraint_error(const System &system, NumericVector< Number > *v=nullptr) const
Tests the constrained degrees of freedom on the numeric vector v, which represents a solution defined...
A do-nothing class for templated methods that expect output iterator arguments.
Storage for DofConstraint right hand sides for a particular problem.
std::string enum_to_string(const T e)
virtual bool closed() const
static FEContinuity get_continuity(const FEType &fe_type)
Returns the input FEType's FEContinuity based on the underlying FEFamily and potentially the Order...
const std::set< boundary_id_type > & get_boundary_ids() const
std::vector< std::unique_ptr< DirichletBoundaries > > _adjoint_dirichlet_boundaries
Data structure containing Dirichlet functions.
static unsigned int n_dofs_at_node(const unsigned int dim, const FEType &fe_t, const ElemType t, const unsigned int n)
virtual const Elem * elem_ptr(const dof_id_type i) const =0
virtual numeric_index_type first_local_index() const =0
virtual unsigned int n_sides() const =0
virtual void edge_reinit(const Elem *elem, const unsigned int edge, const Real tolerance=TOLERANCE, const std::vector< Point > *pts=nullptr, const std::vector< Real > *weights=nullptr)=0
Reinitializes all the physical element-dependent data based on the edge of the element elem...
ParallelType type() const
std::unique_ptr< PeriodicBoundaries > _periodic_boundaries
Data structure containing periodic boundaries.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void left_multiply_transpose(const DenseMatrix< T > &A)
Left multiplies by the transpose of the matrix A.
void process_constraints(MeshBase &)
Postprocesses any constrained degrees of freedom to be constrained only in terms of unconstrained dof...
subdomain_id_type subdomain_id() const
std::unique_ptr< FunctionBase< Gradient > > g
timpi_pure bool verify(const T &r) const
dof_id_type n_constrained_dofs() const
void constrain_element_vector(DenseVector< Number > &rhs, std::vector< dof_id_type > &dofs, bool asymmetric_constraint_rows=true) const
Constrains the element vector.
void max(const T &r, T &o, Request &req) const
FEContinuity
defines an enum for finite element types to libmesh_assert a certain level (or type? Hcurl?) of continuity.
const Node * node_ptr(const unsigned int i) const
const DirichletBoundaries * get_adjoint_dirichlet_boundaries(unsigned int q) const
virtual bool is_vertex(const unsigned int i) const =0
bool _verify_dirichlet_bc_consistency
Flag which determines whether we should do some additional checking of the consistency of the Dirichl...
void send(const unsigned int dest_processor_id, const T &buf, const MessageTag &tag=no_tag) const
void add_constraints_to_send_list()
Adds entries to the _send_list vector corresponding to DoFs which are dependencies for constraint equ...
virtual void right_multiply(const DenseMatrixBase< T > &M2) override final
Performs the operation: (*this) <- (*this) * M3.
SparseMatrix< Number > * matrix
The system matrix.
DofConstraintValueMap _primal_constraint_values
virtual numeric_index_type local_size() const =0
void heterogeneously_constrain_element_matrix_and_vector(DenseMatrix< Number > &matrix, DenseVector< Number > &rhs, std::vector< dof_id_type > &elem_dofs, bool asymmetric_constraint_rows=true, int qoi_index=-1) const
Constrains the element matrix and vector.
const std::set< boundary_id_type > & get_edge_boundary_ids() const
virtual Output component(const FEMContext &, unsigned int i, const Point &p, Real time=0.)
void resize(const unsigned int new_m, const unsigned int new_n)
Resizes the matrix to the specified size and calls zero().
std::map< dof_id_type, Real, std::less< dof_id_type >, Threads::scalable_allocator< std::pair< const dof_id_type, Real > > > DofConstraintRow
A row of the Dof constraint matrix.
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
static std::unique_ptr< NumericVector< T > > build(const Parallel::Communicator &comm, SolverPackage solver_package=libMesh::default_solver_package(), ParallelType parallel_type=AUTOMATIC)
Builds a NumericVector on the processors in communicator comm using the linear solver package specifi...
unsigned int mesh_dimension() const
dof_id_type n_local_dofs() const
virtual unsigned int size() const override final
std::unique_ptr< NumericVector< Number > > current_local_solution
All the values I need to compute my contribution to the simulation at hand.
void add_dirichlet_boundary(const DirichletBoundary &dirichlet_boundary)
Adds a copy of the specified Dirichlet boundary to the system.
void build_constraint_matrix(DenseMatrix< Number > &C, std::vector< dof_id_type > &elem_dofs, const bool called_recursively=false) const
Build the constraint matrix C associated with the element degree of freedom indices elem_dofs...
void get_element_fe(unsigned int var, FEGenericBase< OutputShape > *&fe) const
Accessor for interior finite element object for variable var for the largest dimension in the mesh...
void heterogeneously_constrain_element_residual(DenseVector< Number > &rhs, std::vector< dof_id_type > &elem_dofs, NumericVector< Number > &solution_local) const
Constrains the element residual.
The base class for defining periodic boundaries.
void create_dof_constraints(const MeshBase &, Real time=0)
Rebuilds the raw degree of freedom and DofObject constraints, based on attached DirichletBoundary obj...
std::set< GhostingFunctor * >::const_iterator coupling_functors_end() const
End of range of coupling functors.
std::set< boundary_id_type > b
std::map< const Node *, Real, std::less< const Node * >, Threads::scalable_allocator< std::pair< const Node *const, Real > > > NodeConstraintRow
A row of the Node constraint mapping.
virtual void set(const numeric_index_type i, const T value)=0
Sets v(i) = value.
void edge_boundary_ids(const Elem *const elem, const unsigned short int edge, std::vector< boundary_id_type > &vec_to_fill) const
void cholesky_solve(const DenseVector< T2 > &b, DenseVector< T2 > &x)
For symmetric positive definite (SPD) matrices.
const QBase & get_edge_qrule() const
Accessor for element edge quadrature rule.
void heterogeneously_constrain_element_jacobian_and_residual(DenseMatrix< Number > &matrix, DenseVector< Number > &rhs, std::vector< dof_id_type > &elem_dofs, NumericVector< Number > &solution_local) const
Constrains the element Jacobian and residual.
unsigned int number() const
const FEMap & get_fe_map() const
FEFamily
defines an enum for finite element families.
virtual dof_id_type n_elem() const =0
virtual const Node * node_ptr(const dof_id_type i) const =0
processor_id_type processor_id() const
void constrain_element_residual(DenseVector< Number > &rhs, std::vector< dof_id_type > &elem_dofs, NumericVector< Number > &solution_local) const
Constrains the element residual.
static void compute_periodic_constraints(DofConstraints &constraints, DofMap &dof_map, const PeriodicBoundaries &boundaries, const MeshBase &mesh, const PointLocatorBase *point_locator, const unsigned int variable_number, const Elem *elem)
Computes the constraint matrix contributions (for periodic boundary conditions) corresponding to vari...
const DofMap & get_dof_map() const
processor_id_type processor_id() const
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, const bool add_p_level=true)
Fills the vector di with the local degree of freedom indices associated with side s of element elem A...
A Point defines a location in LIBMESH_DIM dimensional Real space.
dof_id_type node_id(const unsigned int i) const
const Point & point(const unsigned int i) const
virtual numeric_index_type last_local_index() const =0
The constraint matrix storage format.
const std::vector< dof_id_type > & get_send_list() const
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
static void compute_periodic_node_constraints(NodeConstraints &constraints, const PeriodicBoundaries &boundaries, const MeshBase &mesh, const PointLocatorBase *point_locator, const Elem *elem)
Computes the node position constraint equation contributions (for meshes with periodic boundary condi...
This class forms the foundation from which generic finite elements may be derived.
void constrain_nothing(std::vector< dof_id_type > &dofs) const
Does not actually constrain anything, but modifies dofs in the same way as any of the constrain funct...
const QBase & get_element_qrule() const
Accessor for element interior quadrature rule for the dimension of the current _elem.
This class provides single index access to FieldType (i.e.
void vector_mult(DenseVector< T > &dest, const DenseVector< T > &arg) const
Performs the matrix-vector multiplication, dest := (*this) * arg.
void constrain_p_dofs(unsigned int var, const Elem *elem, unsigned int s, unsigned int p)
Constrains degrees of freedom on side s of element elem which correspond to variable number var and t...
void check_for_cyclic_constraints()
Throw an error if we detect any constraint loops, i.e.
const QBase & get_side_qrule() const
Accessor for element side quadrature rule for the dimension of the current _elem. ...
std::vector< dof_id_type > _end_df
Last DOF index (plus 1) on processor p.
bool local_index(dof_id_type dof_index) const
const std::vector< std::vector< OutputShape > > & get_phi() const
const FEType & type() const
NodeConstraints _node_constraints
Data structure containing DofObject constraints.
void enforce_constraints_exactly(const System &system, NumericVector< Number > *v=nullptr, bool homogeneous=false) const
Constrains the numeric vector v, which represents a solution defined on the mesh. ...
spin_mutex spin_mtx
A convenient spin mutex object which can be used for obtaining locks.
void constrain_element_matrix(DenseMatrix< Number > &matrix, std::vector< dof_id_type > &elem_dofs, bool asymmetric_constraint_rows=true) const
Constrains the element matrix.
virtual void localize(std::vector< T > &v_local) const =0
Creates a copy of the global vector in the local vector v_local.
void enforce_constraints_on_jacobian(const NonlinearImplicitSystem &system, SparseMatrix< Number > *jac) const