20 #ifndef GENERIC_PROJECTOR_H 21 #define GENERIC_PROJECTOR_H 27 #include "libmesh/dense_matrix.h" 28 #include "libmesh/dense_vector.h" 29 #include "libmesh/dof_map.h" 30 #include "libmesh/elem.h" 31 #include "libmesh/fe_base.h" 32 #include "libmesh/fe_interface.h" 33 #include "libmesh/int_range.h" 34 #include "libmesh/libmesh_logging.h" 35 #include "libmesh/mesh_tools.h" 36 #include "libmesh/numeric_vector.h" 37 #include "libmesh/quadrature.h" 38 #include "libmesh/system.h" 39 #include "libmesh/threads.h" 40 #include "libmesh/tensor_tools.h" 41 #include "libmesh/libmesh_common.h" 67 template <
typename SendT,
typename T>
69 { converted = received; }
81 template <
typename FFunctor,
typename GFunctor,
82 typename FValue,
typename ProjectionAction>
90 typedef std::unordered_map<dof_id_type, std::vector<dof_id_type>>
NodesToElemMap;
125 ProjectionAction & act_in,
126 const std::vector<unsigned int> & variables_in,
135 if (!nodes_to_elem_in)
161 std::unordered_map<dof_id_type, typename FFunctor::ValuePushType>
ids_to_save;
182 std::pair<
typename FFunctor::ValuePushType,
192 unsigned short node_num,
197 template <
typename InsertInput,
198 typename std::enable_if<
203 template <
typename InsertInput,
204 typename std::enable_if<
209 template <
typename InsertInput,
210 typename std::enable_if<
213 void insert_ids(
const std::vector<dof_id_type> & ids,
214 const std::vector<InsertInput> & vals,
217 template <
typename InsertInput,
218 typename std::enable_if<
221 void insert_ids(
const std::vector<dof_id_type> & ids,
222 const std::vector<InsertInput> & vals,
246 typedef std::pair<const Node *, std::tuple<const Elem *, unsigned short, var_set>>
267 std::unique_ptr<GFunctor>
g;
270 (
const std::vector<dof_id_type> & dof_indices_var,
271 const std::vector<unsigned int> & involved_dofs,
272 unsigned int var_component,
300 typedef std::unordered_multimap<const Node *, std::tuple<const Elem *, unsigned short, var_set>>
ves_multimap;
388 template <
typename Value>
390 (std::unordered_map<
dof_id_type, std::pair<Value, processor_id_type>> & ids_to_push,
391 ProjectionAction & action)
const;
404 template <
typename Val>
427 void insert(
const std::vector<dof_id_type> & dof_indices,
434 unsigned int size = Ue.
size();
436 libmesh_assert_equal_to(size, dof_indices.size());
442 for (
unsigned int i = 0; i != size; ++i)
443 if ((dof_indices[i] >= first) && (dof_indices[i] < last))
459 template <
typename Output>
473 _f(fw.
_f->clone()) {}
483 {
return _f->component(c, i, n, time); }
490 {
return _f->component(c, i, n, time); }
496 std::vector<Output> & )
504 std::vector<dof_id_type> &,
505 std::vector<Output> &)
512 std::vector<dof_id_type> &,
513 std::vector<Output> &)
517 std::unique_ptr<FEMFunctionBase<Output>>
_f;
531 #ifdef LIBMESH_ENABLE_AMR 532 template <
typename Output,
545 const std::vector<unsigned int> * vars) :
548 old_context(sys_in, vars, false)
552 auto make_lookup = [
this](
unsigned int v)
554 const unsigned int vcomp = sys.variable_scalar_number(v,0);
555 if (vcomp >= component_to_var.size())
556 component_to_var.resize(vcomp+1, static_cast<unsigned int>(-1));
557 component_to_var[vcomp] = v;
571 old_context(sys, in.old_context.active_vars(), false),
572 component_to_var(in.component_to_var)
576 static void get_shape_outputs(
FEAbstract & fe);
584 const std::set<unsigned char> & elem_dims =
585 old_context.elem_dimensions();
588 for (
const auto &
dim : elem_dims)
591 if (old_context.active_vars())
592 for (
auto var : *old_context.active_vars())
594 old_context.get_element_fe(var, fe,
dim);
595 get_shape_outputs(*fe);
600 old_context.get_element_fe(var, fe,
dim);
601 get_shape_outputs(*fe);
611 LOG_SCOPE (
"check_old_context(c)",
"OldSolutionBase");
613 if (last_elem != &elem)
617 old_context.pre_fe_reinit(sys, elem.
parent());
630 old_context.pre_fe_reinit(sys, &elem);
644 LOG_SCOPE (
"check_old_context(c,p)",
"OldSolutionBase");
646 if (last_elem != &elem)
650 old_context.pre_fe_reinit(sys, elem.
parent());
661 const Real master_tol = out_of_elem_tol / elem.
hmax() * 2;
664 if (child.close_to_point(p, master_tol))
666 old_context.pre_fe_reinit(sys, &child);
671 (old_context.get_elem().close_to_point(p, master_tol));
678 old_context.pre_fe_reinit(sys, &elem);
687 const Real master_tol = out_of_elem_tol / elem.
hmax() * 2;
689 if (!old_context.get_elem().close_to_point(p, master_tol))
691 libmesh_assert_equal_to
695 if (child.close_to_point(p, master_tol))
697 old_context.pre_fe_reinit(sys, &child);
702 (old_context.get_elem().close_to_point(p, master_tol));
727 template <
typename Output,
742 const std::vector<unsigned int> * vars) :
744 old_solution(old_sol)
747 this->old_context.set_custom_solution(&old_solution);
751 OldSolutionBase<Output, point_output>(in.sys, in.old_context.active_vars()),
752 old_solution(in.old_solution)
755 this->old_context.set_custom_solution(&old_solution);
761 unsigned int elem_dim,
770 bool skip_context_check)
772 LOG_SCOPE (
"eval_at_point()",
"OldSolutionValue");
774 if (!skip_context_check)
775 if (!this->check_old_context(c, p))
779 libmesh_assert_less(i, this->component_to_var.size());
780 unsigned int var = this->component_to_var[i];
783 (this->old_context.*point_output)(var, p, n, this->out_of_elem_tol);
787 template <
typename T = Output,
793 std::vector<Output> & derivs)
795 LOG_SCOPE (
"eval_mixed_derivatives",
"OldSolutionValue");
798 libmesh_assert_less(c.get_elem().get_node_index(&n),
799 c.get_elem().n_vertices());
802 libmesh_assert_less(i, this->component_to_var.size());
803 unsigned int var = this->component_to_var[i];
806 const unsigned int n_mixed = (
dim-1) * (
dim-1);
807 derivs.resize(n_mixed);
812 if (old_dof_object &&
813 old_dof_object->
n_vars(this->sys.number()) &&
814 old_dof_object->
n_comp(this->sys.number(), var))
818 std::vector<dof_id_type> old_ids(n_mixed);
819 std::iota(old_ids.begin(), old_ids.end(), first_old_id);
820 old_solution.get(old_ids, derivs);
824 std::fill(derivs.begin(), derivs.end(), 0);
828 template <
typename T = Output,
831 const FEMContext &,
unsigned int,
unsigned int,
const Node &, std::vector<Output> &)
833 libmesh_error_msg(
"eval_mixed_derivatives should only be applicable for Hermite finite element " 834 "types. I don't know how you got here");
838 unsigned int node_num,
839 unsigned int var_num,
840 std::vector<dof_id_type> & indices,
841 std::vector<DofValueType> & values)
843 LOG_SCOPE (
"eval_old_dofs(node)",
"OldSolutionValue");
849 this->sys.get_dof_map().dof_indices(elem, node_num, indices, var_num);
851 std::vector<dof_id_type> old_indices;
853 this->sys.get_dof_map().old_dof_indices(elem, node_num, old_indices, var_num);
855 libmesh_assert_equal_to (old_indices.size(), indices.size());
859 bool invalid_old_index =
false;
860 for (
const auto & di : old_indices)
862 invalid_old_index =
true;
864 values.resize(old_indices.size());
865 if (invalid_old_index)
873 values[i] = old_solution(di);
877 old_solution.get(old_indices, values);
883 unsigned int sys_num,
884 unsigned int var_num,
885 std::vector<dof_id_type> & indices,
886 std::vector<DofValueType> & values)
888 LOG_SCOPE (
"eval_old_dofs(elem)",
"OldSolutionValue");
892 const Elem & old_elem =
909 libmesh_assert_greater(elem.
n_systems(), sys_num);
911 const std::pair<unsigned int, unsigned int>
913 const unsigned int vg = vg_and_offset.first;
914 const unsigned int vig = vg_and_offset.second;
916 unsigned int n_comp = old_dof_object.
n_comp_group(sys_num,vg);
917 n_comp = std::min(n_comp, nc);
919 std::vector<dof_id_type> old_dof_indices(n_comp);
921 for (
unsigned int i=0; i != n_comp; ++i)
924 old_dof_object.
dof_number(sys_num, vg, vig, i, n_comp);
930 old_dof_indices[i] = d_old;
934 values.resize(n_comp);
935 old_solution.get(old_dof_indices, values);
937 for (
unsigned int i=n_comp; i != nc; ++i)
945 values.resize(nc, 0);
986 #ifdef LIBMESH_USE_COMPLEX_NUMBERS 1016 #endif // LIBMESH_USE_COMPLEX_NUMBERS 1027 bool extra_hanging_dofs,
1030 LOG_SCOPE (
"Number eval_at_node()",
"OldSolutionValue");
1037 libmesh_assert_less(i, this->component_to_var.size());
1038 unsigned int var = this->component_to_var[i];
1053 if (old_dof_object &&
1054 (!extra_hanging_dofs ||
1057 old_dof_object->
n_vars(sys.number()) &&
1058 old_dof_object->
n_comp(sys.number(), var))
1061 old_dof_object->
dof_number(sys.number(), var, 0);
1062 return old_solution(old_id);
1065 return this->eval_at_point(c, i, n, 0,
false);
1076 bool extra_hanging_dofs,
1079 LOG_SCOPE (
"Number eval_at_node()",
"OldSolutionValue");
1086 libmesh_assert_less(i, this->component_to_var.size());
1087 unsigned int var = this->component_to_var[i];
1104 if (old_dof_object &&
1105 (!extra_hanging_dofs ||
1108 old_dof_object->
n_vars(sys.number()) &&
1109 old_dof_object->
n_comp(sys.number(), var))
1117 return_val(
dim) = old_solution(old_id);
1123 return this->eval_at_point(c, i, n, 0,
false);
1134 unsigned int elem_dim,
1136 bool extra_hanging_dofs,
1139 LOG_SCOPE (
"Gradient eval_at_node()",
"OldSolutionValue");
1146 libmesh_assert_less(i, this->component_to_var.size());
1147 unsigned int var = this->component_to_var[i];
1162 if (old_dof_object &&
1163 (!extra_hanging_dofs ||
1166 old_dof_object->
n_vars(sys.number()) &&
1167 old_dof_object->
n_comp(sys.number(), var))
1170 for (
unsigned int d = 0; d != elem_dim; ++d)
1173 old_dof_object->
dof_number(sys.number(), var, d+1);
1174 g(d) = old_solution(old_id);
1179 return this->eval_at_point(c, i, n, 0,
false);
1193 libmesh_error_msg(
"You shouldn't need to call eval_at_node for the gradient " 1194 "functor for a vector-valued finite element type");
1198 template <
typename Output,
1205 #endif // LIBMESH_ENABLE_AMR 1211 template <
typename FFunctor,
typename GFunctor,
1212 typename FValue,
typename ProjectionAction>
1216 LOG_SCOPE (
"project",
"GenericProjector");
1220 done_saving_ids =
false;
1224 ProjectionAction action(master_action);
1227 std::unordered_map<dof_id_type, std::pair<typename FFunctor::ValuePushType, processor_id_type>>
1235 std::vector<node_projection> vertices(sort_work.
vertices.begin(),
1238 done_saving_ids = sort_work.
edges.empty() &&
1248 done_saving_ids = sort_work.
sides.empty() && sort_work.
interiors.empty();
1250 this->send_and_insert_dof_values(ids_to_push, action);
1253 std::vector<node_projection> edges(sort_work.
edges.begin(), sort_work.
edges.end());
1260 done_saving_ids = sort_work.
interiors.empty();
1262 this->send_and_insert_dof_values(ids_to_push, action);
1265 std::vector<node_projection> sides(sort_work.
sides.begin(), sort_work.
sides.end());
1272 done_saving_ids =
true;
1273 this->send_and_insert_dof_values(ids_to_push, action);
1283 template <
typename FFunctor,
typename GFunctor,
1284 typename FValue,
typename ProjectionAction>
1296 for (
const auto & var : this->projector.variables)
1303 const std::set<unsigned char> & elem_dims =
1304 context.elem_dimensions();
1306 for (
const auto &
dim : elem_dims)
1308 context.get_element_fe( var, fe,
dim );
1312 context.get_side_fe( var, side_fe,
dim );
1314 context.get_edge_fe( var, edge_fe );
1335 this->conts[var] = cont;
1350 f.init_context(context);
1354 template <
typename FFunctor,
typename GFunctor,
1355 typename FValue,
typename ProjectionAction>
1362 for (
const auto & var : this->projector.variables)
1363 if (this->conts[var] ==
C_ONE)
1366 g = std::make_unique<GFunctor>(*p.
master_g);
1367 g->init_context(context);
1372 template <
typename FFunctor,
typename GFunctor,
typename FValue,
typename ProjectionAction>
1373 template <
typename InsertInput,
1374 typename std::enable_if<
1381 libmesh_error_msg(
"ID insertion should only occur when the provided input matches that " 1382 "expected by the ProjectionAction");
1385 template <
typename FFunctor,
typename GFunctor,
typename FValue,
typename ProjectionAction>
1386 template <
typename InsertInput,
1387 typename std::enable_if<
1399 auto iter = new_ids_to_push.find(
id);
1400 if (iter == new_ids_to_push.end())
1401 action.insert(
id, val);
1405 iter->second = std::make_pair(val, pid);
1407 if (!this->projector.done_saving_ids)
1410 new_ids_to_save[id] = val;
1414 template <
typename FFunctor,
typename GFunctor,
typename FValue,
typename ProjectionAction>
1415 template <
typename InsertInput,
1416 typename std::enable_if<
1421 const std::vector<dof_id_type> &,
1422 const std::vector<InsertInput> &,
1425 libmesh_error_msg(
"ID insertion should only occur when the provided input matches that " 1426 "expected by the ProjectionAction");
1429 template <
typename FFunctor,
typename GFunctor,
typename FValue,
typename ProjectionAction>
1430 template <
typename InsertInput,
1431 typename std::enable_if<
1436 const std::vector<dof_id_type> & ids,
1437 const std::vector<InsertInput> & vals,
1440 libmesh_assert_equal_to(ids.size(), vals.size());
1450 const InsertInput & val = vals[i];
1452 auto iter = new_ids_to_push.find(
id);
1453 if (iter == new_ids_to_push.end())
1454 action.insert(
id, val);
1458 iter->second = std::make_pair(val, pid);
1460 if (!this->projector.done_saving_ids)
1463 new_ids_to_save[id] = val;
1468 template <
typename FFunctor,
typename GFunctor,
1469 typename FValue,
typename ProjectionAction>
1478 template <
typename FFunctor,
typename GFunctor,
1479 typename FValue,
typename ProjectionAction>
1483 LOG_SCOPE (
"SortAndCopy::operator()",
"GenericProjector");
1509 std::unordered_map<const Node *, std::pair<var_set, var_set>> copied_nodes;
1516 std::vector<unsigned short> extra_hanging_dofs;
1517 bool all_extra_hanging_dofs =
true;
1518 for (
auto v_num : this->projector.variables)
1520 if (extra_hanging_dofs.size() <= v_num)
1521 extra_hanging_dofs.resize(v_num+1,
false);
1522 extra_hanging_dofs[v_num] =
1525 if (!extra_hanging_dofs[v_num])
1526 all_extra_hanging_dofs =
false;
1529 for (
const auto & elem : range)
1533 bool copy_this_elem =
false;
1535 #ifdef LIBMESH_ENABLE_AMR 1537 if (f.is_grid_projection())
1546 if (!old_dof_object &&
1556 copy_this_elem =
true;
1559 bool reinitted =
false;
1561 const unsigned int p_level = elem->p_level();
1565 const bool copy_possible =
1570 std::vector<typename FFunctor::ValuePushType> Ue(1);
1571 std::vector<dof_id_type> elem_dof_ids(1);
1573 for (
auto v_num : this->projector.variables)
1587 context.pre_fe_reinit(
system, elem);
1590 f.eval_old_dofs(*elem, fe_type, sys_num, v_num,
1593 action.insert(elem_dof_ids[0], Ue[0]);
1598 #endif // LIBMESH_ENABLE_AMR 1600 const int dim = elem->dim();
1602 const unsigned int n_vertices = elem->n_vertices();
1603 const unsigned int n_edges = elem->n_edges();
1604 const unsigned int n_nodes = elem->n_nodes();
1607 const unsigned int n_sides = (
dim > 1) * elem->n_sides();
1610 var_set vertex_vars, edge_vars, side_vars;
1614 const bool has_edge_nodes = (
n_nodes > n_vertices &&
dim > 2);
1617 const bool has_side_nodes =
1618 (
n_nodes > n_vertices + ((
dim > 2) * n_edges));
1622 const bool has_interior_nodes =
1623 (
n_nodes > n_vertices + ((
dim > 2) * n_edges) + n_sides);
1625 for (
auto v_num : this->projector.variables)
1631 const auto & dof_map = this->projector.system.get_dof_map();
1632 const auto vg = dof_map.var_group_from_var_number(v_num);
1633 const bool add_p_level = dof_map.should_p_refine(vg);
1645 vertex_vars.insert(vertex_vars.end(), v_num);
1651 edge_vars.insert(edge_vars.end(), v_num);
1658 side_vars.insert(side_vars.end(), v_num);
1663 for (
unsigned int n = 0; n !=
n_nodes; ++n)
1664 if (elem->is_face(n))
1667 side_vars.insert(side_vars.end(), v_num);
1673 (has_interior_nodes &&
1676 #ifdef LIBMESH_ENABLE_AMR 1679 if ((f.is_grid_projection() &&
1682 elem->p_level() == 0 &&
1688 #endif // LIBMESH_ENABLE_AMR 1691 if (interiors.empty() || interiors.back() != elem)
1692 interiors.push_back(elem);
1699 auto erase_covered_vars = []
1702 auto covered_range = covered.equal_range(node);
1703 for (
const auto & v_ent :
as_range(covered_range))
1704 for (
const unsigned int var_covered :
1705 std::get<2>(v_ent.second))
1706 remaining.erase(var_covered);
1709 auto erase_nonhanging_vars = [&extra_hanging_dofs]
1712 auto covered_range = covered.equal_range(node);
1713 for (
const auto & v_ent :
as_range(covered_range))
1714 for (
const unsigned int var_covered :
1715 std::get<2>(v_ent.second))
1716 if (!extra_hanging_dofs[var_covered])
1717 remaining.erase(var_covered);
1720 auto erase_copied_vars = [&copied_nodes, &extra_hanging_dofs]
1721 (
const Node * node,
bool is_vertex,
var_set & remaining)
1723 auto copying_range = copied_nodes.equal_range(node);
1724 for (
const auto & v_ent :
as_range(copying_range))
1726 for (
const unsigned int var_covered :
1728 if (is_vertex || !extra_hanging_dofs[var_covered])
1729 remaining.erase(var_covered);
1731 for (
const unsigned int var_covered :
1732 v_ent.second.second)
1733 if (!is_vertex || !extra_hanging_dofs[var_covered])
1734 remaining.erase(var_covered);
1738 for (
unsigned int v=0; v != n_vertices; ++v)
1740 const Node * node = elem->node_ptr(v);
1742 auto remaining_vars = vertex_vars;
1744 erase_covered_vars(node, remaining_vars, vertices);
1746 if (remaining_vars.empty())
1749 if (!all_extra_hanging_dofs)
1751 erase_nonhanging_vars(node, remaining_vars, edges);
1752 if (remaining_vars.empty())
1755 erase_nonhanging_vars(node, remaining_vars, sides);
1756 if (remaining_vars.empty())
1760 erase_copied_vars(node,
true, remaining_vars);
1761 if (remaining_vars.empty())
1766 std::vector<dof_id_type> node_dof_ids;
1767 std::vector<typename FFunctor::ValuePushType> values;
1769 for (
auto var : remaining_vars)
1771 f.eval_old_dofs(*elem, v, var, node_dof_ids, values);
1772 insert_ids(node_dof_ids, values, node->
processor_id());
1774 copied_nodes[node].first.insert(remaining_vars.begin(),
1775 remaining_vars.end());
1776 this->find_dofs_to_send(*node, *elem, v, remaining_vars);
1780 (node, std::make_tuple(elem, v, std::move(remaining_vars)));
1785 for (
unsigned int e=0; e != n_edges; ++e)
1787 const Node * node = elem->node_ptr(n_vertices+e);
1789 auto remaining_vars = edge_vars;
1791 erase_covered_vars(node, remaining_vars, edges);
1792 if (remaining_vars.empty())
1795 erase_covered_vars(node, remaining_vars, sides);
1796 if (remaining_vars.empty())
1799 if (!all_extra_hanging_dofs)
1801 erase_nonhanging_vars(node, remaining_vars, vertices);
1802 if (remaining_vars.empty())
1806 erase_copied_vars(node,
false, remaining_vars);
1807 if (remaining_vars.empty())
1812 std::vector<dof_id_type> edge_dof_ids;
1813 std::vector<typename FFunctor::ValuePushType> values;
1815 for (
auto var : remaining_vars)
1817 f.eval_old_dofs(*elem, n_vertices+e, var, edge_dof_ids, values);
1818 insert_ids(edge_dof_ids, values, node->
processor_id());
1821 copied_nodes[node].second.insert(remaining_vars.begin(),
1822 remaining_vars.end());
1823 this->find_dofs_to_send(*node, *elem, n_vertices+e, remaining_vars);
1827 (node, std::make_tuple(elem, e, std::move(remaining_vars)));
1833 for (
unsigned int side=0; side != n_sides; ++side)
1835 const Node * node =
nullptr;
1836 unsigned short node_num = n_vertices+(
dim>2)*n_edges+side;
1838 node = elem->node_ptr(node_num);
1842 for (
unsigned int n = 0; n !=
n_nodes; ++n)
1844 if (!elem->is_face(n))
1847 if (elem->is_node_on_side(n, side))
1850 node = elem->node_ptr(node_num);
1859 auto remaining_vars = side_vars;
1861 erase_covered_vars(node, remaining_vars, edges);
1862 if (remaining_vars.empty())
1865 erase_covered_vars(node, remaining_vars, sides);
1866 if (remaining_vars.empty())
1869 if (!all_extra_hanging_dofs)
1871 erase_nonhanging_vars(node, remaining_vars, vertices);
1872 if (remaining_vars.empty())
1876 erase_copied_vars(node,
false, remaining_vars);
1877 if (remaining_vars.empty())
1882 std::vector<dof_id_type> side_dof_ids;
1883 std::vector<typename FFunctor::ValuePushType> values;
1885 for (
auto var : remaining_vars)
1887 f.eval_old_dofs(*elem, node_num, var, side_dof_ids, values);
1888 insert_ids(side_dof_ids, values, node->
processor_id());
1891 copied_nodes[node].second.insert(remaining_vars.begin(),
1892 remaining_vars.end());
1893 this->find_dofs_to_send(*node, *elem, node_num, remaining_vars);
1897 (node, std::make_tuple(elem, side, std::move(remaining_vars)));
1904 std::vector<typename FFunctor::ValuePushType> U;
1905 std::vector<dof_id_type> dof_ids;
1907 for (
auto v_num : this->projector.variables)
1914 f.eval_old_dofs(*elem, fe_type, sys_num, v_num,
1916 action.insert(dof_ids, U);
1918 if (has_interior_nodes)
1920 f.eval_old_dofs(*elem,
n_nodes-1, v_num, dof_ids, U);
1921 action.insert(dof_ids, U);
1929 template <
typename FFunctor,
typename GFunctor,
1930 typename FValue,
typename ProjectionAction>
1936 for (
const auto & pair : other_mm)
1938 const Node * node = pair.first;
1939 auto remaining_vars = std::get<2>(pair.second);
1941 auto my_range =
self.equal_range(node);
1942 for (
const auto & v_ent :
as_range(my_range))
1943 for (
const unsigned int var_covered :
1944 std::get<2>(v_ent.second))
1945 remaining_vars.erase(var_covered);
1947 if (!remaining_vars.empty())
1949 (node, std::make_tuple(std::get<0>(pair.second),
1950 std::get<1>(pair.second),
1951 std::move(remaining_vars)));
1955 merge_multimaps(vertices, other.
vertices);
1956 merge_multimaps(edges, other.
edges);
1957 merge_multimaps(sides, other.
sides);
1959 std::sort(interiors.begin(), interiors.end());
1960 std::vector<const Elem *> other_interiors = other.
interiors;
1961 std::sort(other_interiors.begin(), other_interiors.end());
1962 std::vector<const Elem *> merged_interiors;
1963 std::set_union(interiors.begin(), interiors.end(),
1964 other_interiors.begin(), other_interiors.end(),
1965 std::back_inserter(merged_interiors));
1966 interiors.swap(merged_interiors);
1973 template <
typename Output,
typename Input>
1975 raw_value(
const Input & input,
unsigned int )
1980 template <
typename Output,
template <
typename>
class WrapperClass,
typename T>
1984 raw_value(
const WrapperClass<T> & input,
unsigned int index)
1986 return input.slice(index);
1989 template <
typename T>
1990 typename T::value_type
1991 grad_component(
const T &,
unsigned int);
1993 template <
typename T>
1995 grad_component(
const VectorValue<T> & grad,
unsigned int component)
1997 return grad(component);
2000 template <
typename T>
2002 grad_component(
const TensorValue<T> & grad,
unsigned int component)
2004 libmesh_error_msg(
"This function should only be called for Hermites. " 2005 "I don't know how you got here");
2006 return grad(component, component);
2012 template <
typename FFunctor,
typename GFunctor,
2013 typename FValue,
typename ProjectionAction>
2014 void GenericProjector<FFunctor, GFunctor, FValue, ProjectionAction>::ProjectVertices::operator()
2017 LOG_SCOPE (
"project_vertices",
"GenericProjector");
2023 std::vector<unsigned short> extra_hanging_dofs;
2024 for (
auto v_num : this->projector.variables)
2026 if (extra_hanging_dofs.size() <= v_num)
2027 extra_hanging_dofs.resize(v_num+1,
false);
2028 extra_hanging_dofs[v_num] =
2032 for (
const auto & v_pair : range)
2034 const Node & vertex = *v_pair.first;
2035 const Elem & elem = *std::get<0>(v_pair.second);
2036 const unsigned int n = std::get<1>(v_pair.second);
2037 const var_set & vertex_vars = std::get<2>(v_pair.second);
2039 context.pre_fe_reinit(
system, &elem);
2041 this->find_dofs_to_send(vertex, elem, n, vertex_vars);
2045 for (
const auto & var : vertex_vars)
2048 const FEType & base_fe_type = variable.
type();
2049 const unsigned int var_component =
2056 const FEFieldType field_type = this->field_types[var];
2060 libmesh_assert_equal_to(vertex.
n_comp(sys_num, var), 0);
2062 else if (cont ==
C_ZERO ||
2068 libmesh_assert_equal_to(vertex.
n_comp(sys_num, var), 0);
2072 const FValue val = f.eval_at_node
2073 (context, var_component, 0,
2074 vertex, extra_hanging_dofs[var],
system.
time);
2078 libmesh_assert_equal_to(vertex.
n_comp(sys_num, var), elem.
dim());
2087 const auto insert_val =
2088 raw_value<typename ProjectionAction::InsertInput>(val, i);
2105 else if (cont ==
C_ONE)
2115 const int dim = elem.
dim();
2122 for (
const auto e_id : (*this->projector.nodes_to_elem)[vertex.
id()])
2125 libmesh_assert_equal_to(
dim, e.
dim());
2128 #ifdef LIBMESH_ENABLE_AMR 2129 bool is_old_vertex =
true;
2132 const int i_am_child =
2138 const bool is_old_vertex =
false;
2145 f.eval_at_node(context,
2149 extra_hanging_dofs[var],
2153 typename GFunctor::FunctorValue grad =
2155 g->eval_at_node(context,
2159 extra_hanging_dofs[var],
2161 g->eval_at_point(context,
2167 insert_id(first_id+1, grad.slice(0),
2170 if (
dim > 1 && is_old_vertex && f.is_grid_projection())
2172 for (
int i = 1; i <
dim; ++i)
2173 insert_id(first_id+i+1, grad.slice(i),
2177 std::vector<FValue> derivs;
2178 f.eval_mixed_derivatives
2179 (context, var_component,
dim, vertex, derivs);
2181 insert_id(first_id+
dim+1+i, derivs[i],
2193 nxplus = elem.
point(n);
2194 nxminus(0) -= delta_x;
2195 nxplus(0) += delta_x;
2196 typename GFunctor::FunctorValue gxminus =
2197 g->eval_at_point(context,
2202 typename GFunctor::FunctorValue gxplus =
2203 g->eval_at_point(context,
2209 insert_id(first_id+2, grad.slice(1),
2212 insert_id(first_id+3,
2213 (grad_component(gxplus, 1) - grad_component(gxminus, 1)) / 2. / delta_x,
2220 insert_id(first_id+4, grad.slice(2),
2223 insert_id(first_id+5,
2224 (grad_component(gxplus, 2) - grad_component(gxminus, 2)) / 2. / delta_x,
2229 nyplus = elem.
point(n);
2230 nyminus(1) -= delta_x;
2231 nyplus(1) += delta_x;
2232 typename GFunctor::FunctorValue gyminus =
2233 g->eval_at_point(context,
2238 typename GFunctor::FunctorValue gyplus =
2239 g->eval_at_point(context,
2245 insert_id(first_id+6,
2246 (grad_component(gyplus, 2) - grad_component(gyminus, 2)) / 2. / delta_x,
2250 nxmyp = elem.
point(n),
2251 nxpym = elem.
point(n),
2252 nxpyp = elem.
point(n);
2253 nxmym(0) -= delta_x;
2254 nxmym(1) -= delta_x;
2255 nxmyp(0) -= delta_x;
2256 nxmyp(1) += delta_x;
2257 nxpym(0) += delta_x;
2258 nxpym(1) -= delta_x;
2259 nxpyp(0) += delta_x;
2260 nxpyp(1) += delta_x;
2261 typename GFunctor::FunctorValue gxmym =
2262 g->eval_at_point(context,
2267 typename GFunctor::FunctorValue gxmyp =
2268 g->eval_at_point(context,
2273 typename GFunctor::FunctorValue gxpym =
2274 g->eval_at_point(context,
2279 typename GFunctor::FunctorValue gxpyp =
2280 g->eval_at_point(context,
2285 FValue gxzplus = (grad_component(gxpyp, 2) - grad_component(gxmyp, 2))
2287 FValue gxzminus = (grad_component(gxpym, 2) - grad_component(gxmym, 2))
2290 insert_id(first_id+7,
2291 (gxzplus - gxzminus) / 2. / delta_x,
2294 #endif // LIBMESH_DIM > 2 2296 #endif // LIBMESH_DIM > 1 2303 libmesh_assert_equal_to(
2308 this->projector.system.get_dof_map().should_p_refine(
2309 this->projector.system.get_dof_map().var_group_from_var_number(var))),
2310 (
unsigned int)(1 +
dim));
2313 f.eval_at_node(context, var_component,
dim,
2314 vertex, extra_hanging_dofs[var],
2317 typename GFunctor::FunctorValue grad =
2319 g->eval_at_node(context, var_component,
dim,
2320 vertex, extra_hanging_dofs[var],
2322 g->eval_at_point(context, var_component, vertex,
2324 for (
int i=0; i!=
dim; ++i)
2325 insert_id(first_id + i + 1, grad.slice(i),
2330 libmesh_error_msg(
"Unknown continuity " << cont);
2336 template <
typename FFunctor,
typename GFunctor,
2337 typename FValue,
typename ProjectionAction>
2341 LOG_SCOPE (
"project_edges",
"GenericProjector");
2345 for (
const auto & e_pair : range)
2347 const Elem & elem = *std::get<0>(e_pair.second);
2351 #ifdef LIBMESH_ENABLE_AMR 2352 if (f.is_grid_projection() &&
2358 #endif // LIBMESH_ENABLE_AMR 2360 const Node & edge_node = *e_pair.first;
2361 const int dim = elem.
dim();
2362 const var_set & edge_vars = std::get<2>(e_pair.second);
2364 const unsigned short edge_num = std::get<1>(e_pair.second);
2365 const unsigned short node_num = elem.
n_vertices() + edge_num;
2366 context.edge = edge_num;
2367 context.pre_fe_reinit(
system, &elem);
2369 this->find_dofs_to_send(edge_node, elem, node_num, edge_vars);
2373 for (
const auto & var : edge_vars)
2376 const FEType & base_fe_type = variable.
type();
2377 const unsigned int var_component =
2383 FEType fe_type = base_fe_type;
2393 if (fe_type.
order > 1)
2397 FValue fval = f.eval_at_point
2398 (context, var_component, edge_node,
system.
time,
2410 const auto insert_val =
2411 raw_value<typename ProjectionAction::InsertInput>(fval, i);
2413 insert_id(dof_id, insert_val, pid);
2420 insert_id(dof_id, fval, pid);
2426 #ifdef LIBMESH_ENABLE_AMR 2434 #endif // LIBMESH_ENABLE_AMR 2438 context.get_element_fe( var, fe,
dim );
2440 context.get_edge_fe( var, edge_fe );
2445 #ifdef LIBMESH_ENABLE_AMR 2451 #ifdef LIBMESH_ENABLE_AMR 2454 std::vector<Point> fine_points;
2456 std::unique_ptr<FEGenericBase<typename FFunctor::RealType>> fine_fe
2459 std::unique_ptr<QBase> qrule
2461 fine_fe->attach_quadrature_rule(qrule.get());
2463 const std::vector<Point> & child_xyz =
2466 for (
unsigned int c = 0, nc = elem.
n_children();
2472 fine_fe->edge_reinit(elem.
child_ptr(c), context.edge);
2473 fine_points.insert(fine_points.end(),
2478 std::vector<Point> fine_qp;
2481 context.elem_fe_reinit(&fine_qp);
2484 #endif // LIBMESH_ENABLE_AMR 2485 context.edge_fe_reinit();
2487 const std::vector<dof_id_type> & dof_indices =
2488 context.get_dof_indices(var);
2490 std::vector<unsigned int> edge_dofs;
2492 context.edge, edge_dofs);
2494 this->construct_projection
2495 (dof_indices, edge_dofs, var_component,
2496 &edge_node, proj_fe);
2502 template <
typename FFunctor,
typename GFunctor,
2503 typename FValue,
typename ProjectionAction>
2507 LOG_SCOPE (
"project_sides",
"GenericProjector");
2511 for (
const auto & s_pair : range)
2513 const Elem & elem = *std::get<0>(s_pair.second);
2517 #ifdef LIBMESH_ENABLE_AMR 2518 if (f.is_grid_projection() &&
2524 #endif // LIBMESH_ENABLE_AMR 2526 const Node & side_node = *s_pair.first;
2527 const int dim = elem.
dim();
2528 const var_set & side_vars = std::get<2>(s_pair.second);
2530 const unsigned int side_num = std::get<1>(s_pair.second);
2531 unsigned short node_num = elem.
n_vertices()+side_num;
2546 context.side = side_num;
2547 context.pre_fe_reinit(
system, &elem);
2549 this->find_dofs_to_send(side_node, elem, node_num, side_vars);
2553 for (
const auto & var : side_vars)
2556 const FEType & base_fe_type = variable.
type();
2557 const unsigned int var_component =
2563 FEType fe_type = base_fe_type;
2580 if (fe_type.
order > 1 &&
2581 side_node.
n_comp(sys_num, var))
2585 FValue fval = f.eval_at_point
2586 (context, var_component, side_node,
system.
time,
2598 const auto insert_val =
2599 raw_value<typename ProjectionAction::InsertInput>(fval, i);
2601 insert_id(dof_id, insert_val, pid);
2608 insert_id(dof_id, fval, pid);
2614 #ifdef LIBMESH_ENABLE_AMR 2622 #endif // LIBMESH_ENABLE_AMR 2626 context.get_element_fe( var, fe,
dim );
2628 context.get_side_fe( var, side_fe );
2633 #ifdef LIBMESH_ENABLE_AMR 2639 #ifdef LIBMESH_ENABLE_AMR 2642 std::vector<Point> fine_points;
2644 std::unique_ptr<FEGenericBase<typename FFunctor::RealType>> fine_fe
2647 std::unique_ptr<QBase> qrule
2649 fine_fe->attach_quadrature_rule(qrule.get());
2651 const std::vector<Point> & child_xyz =
2654 for (
unsigned int c = 0, nc = elem.
n_children();
2660 fine_fe->reinit(elem.
child_ptr(c), context.side);
2661 fine_points.insert(fine_points.end(),
2666 std::vector<Point> fine_qp;
2669 context.elem_fe_reinit(&fine_qp);
2672 #endif // LIBMESH_ENABLE_AMR 2673 context.side_fe_reinit();
2675 const std::vector<dof_id_type> & dof_indices =
2676 context.get_dof_indices(var);
2678 std::vector<unsigned int> side_dofs;
2680 context.side, side_dofs, add_p_level);
2682 this->construct_projection
2683 (dof_indices, side_dofs, var_component,
2684 &side_node, proj_fe);
2690 template <
typename FFunctor,
typename GFunctor,
2691 typename FValue,
typename ProjectionAction>
2695 LOG_SCOPE (
"project_interiors",
"GenericProjector");
2700 for (
const auto & elem : range)
2702 unsigned char dim = cast_int<unsigned char>(elem->dim());
2704 context.pre_fe_reinit(
system, elem);
2708 for (
const auto & var : this->projector.variables)
2715 const FEType & base_fe_type = variable.
type();
2721 context.get_element_fe( var, fe,
dim );
2723 FEType fe_type = base_fe_type;
2726 const bool add_p_level = dof_map.should_p_refine(vg);
2729 fe_type.
order = fe_type.
order + add_p_level * elem->p_level();
2731 const unsigned int var_component =
2739 if (fe_type.
order > 1)
2741 const unsigned int first_interior_node =
2742 (elem->n_vertices() +
2743 ((elem->dim() > 2) * elem->n_edges()) +
2744 ((elem->dim() > 1) * elem->n_sides()));
2745 const unsigned int n_nodes = elem->n_nodes();
2748 for (
unsigned int n = first_interior_node; n <
n_nodes; ++n)
2750 const Node & interior_node = elem->node_ref(n);
2752 FValue fval = f.eval_at_point
2753 (context, var_component, interior_node,
2761 for (
unsigned int i = 0; i < elem->dim(); ++i)
2767 const auto insert_val =
2768 raw_value<typename ProjectionAction::InsertInput>(fval, i);
2770 insert_id(dof_id, insert_val, pid);
2777 insert_id(dof_id, fval, pid);
2784 #ifdef LIBMESH_ENABLE_AMR 2787 std::vector<Point> fine_points;
2789 std::unique_ptr<FEGenericBase<typename FFunctor::RealType>> fine_fe
2792 std::unique_ptr<QBase> qrule
2794 fine_fe->attach_quadrature_rule(qrule.get());
2796 const std::vector<Point> & child_xyz =
2799 for (
auto & child : elem->child_ref_range())
2801 fine_fe->reinit(&child);
2802 fine_points.insert(fine_points.end(),
2807 std::vector<Point> fine_qp;
2810 context.elem_fe_reinit(&fine_qp);
2813 #endif // LIBMESH_ENABLE_AMR 2814 context.elem_fe_reinit();
2816 const std::vector<dof_id_type> & dof_indices =
2817 context.get_dof_indices(var);
2819 const unsigned int n_dofs =
2820 cast_int<unsigned int>(dof_indices.size());
2822 std::vector<unsigned int> all_dofs(n_dofs);
2823 std::iota(all_dofs.begin(), all_dofs.end(), 0);
2825 this->construct_projection
2826 (dof_indices, all_dofs, var_component,
2834 template <
typename FFunctor,
typename GFunctor,
2835 typename FValue,
typename ProjectionAction>
2838 ProjectionAction>::SubFunctor::find_dofs_to_send
2853 std::vector<dof_id_type> node_dof_ids, patch_dof_ids;
2854 for (
const auto & var : vars)
2861 dof_map.
dof_indices(elem, node_num, node_dof_ids, var);
2864 node_dof_ids.end()));
2865 const std::vector<dof_id_type> & patch =
2866 (*this->projector.nodes_to_elem)[node.
id()];
2867 for (
const auto & elem_id : patch)
2869 const Elem & patch_elem =
mesh.elem_ref(elem_id);
2873 std::sort(patch_dof_ids.begin(), patch_dof_ids.end());
2877 std::vector<dof_id_type> diff_ids(node_dof_ids.size());
2878 auto it = std::set_difference(node_dof_ids.begin(), node_dof_ids.end(),
2879 patch_dof_ids.begin(), patch_dof_ids.end(), diff_ids.begin());
2880 diff_ids.resize(it-diff_ids.begin());
2881 node_dof_ids.swap(diff_ids);
2882 if (node_dof_ids.empty())
2887 for (
auto id : node_dof_ids)
2894 template <
typename FFunctor,
typename GFunctor,
2895 typename FValue,
typename ProjectionAction>
2896 template <
typename Value>
2900 (std::unordered_map<
dof_id_type, std::pair<Value, processor_id_type>> & ids_to_push,
2901 ProjectionAction & action)
const 2903 LOG_SCOPE (
"send_and_insert_dof_values",
"GenericProjector");
2907 std::unordered_map<processor_id_type, std::vector<dof_id_type>>
2908 pushed_dof_ids, received_dof_ids;
2909 std::unordered_map<processor_id_type, std::vector<typename TypeToSend<Value>::type>>
2910 pushed_dof_values, received_dof_values;
2911 for (
auto & id_val_pid : ids_to_push)
2916 pushed_dof_ids[pid].push_back(id_val_pid.first);
2917 pushed_dof_values[pid].push_back(
convert_to_send(id_val_pid.second.first));
2923 auto ids_action_functor =
2924 [&action, &received_dof_ids, &received_dof_values]
2926 const std::vector<dof_id_type> & data)
2928 auto iter = received_dof_values.find(pid);
2929 if (iter == received_dof_values.end())
2932 received_dof_ids.end());
2933 received_dof_ids[pid] = data;
2937 auto & vals = iter->second;
2938 std::size_t vals_size = vals.size();
2939 libmesh_assert_equal_to(vals_size, data.size());
2940 for (std::size_t i=0; i != vals_size; ++i)
2944 action.insert(data[i], val);
2946 received_dof_values.erase(iter);
2952 auto values_action_functor =
2953 [&action, &received_dof_ids, &received_dof_values]
2955 const std::vector<typename TypeToSend<Value>::type> & data)
2957 auto iter = received_dof_ids.find(pid);
2958 if (iter == received_dof_ids.end())
2962 received_dof_values.end());
2963 received_dof_values[pid] = data;
2967 auto & ids = iter->second;
2968 std::size_t ids_size = ids.size();
2969 libmesh_assert_equal_to(ids_size, data.size());
2970 for (std::size_t i=0; i != ids_size; ++i)
2974 action.insert(ids[i], val);
2976 received_dof_ids.erase(iter);
2980 Parallel::push_parallel_vector_data
2981 (
system.
comm(), pushed_dof_ids, ids_action_functor);
2983 Parallel::push_parallel_vector_data
2984 (
system.
comm(), pushed_dof_values, values_action_functor);
2994 template <
typename FFunctor,
typename GFunctor,
2995 typename FValue,
typename ProjectionAction>
2998 (
const std::vector<dof_id_type> & dof_indices_var,
2999 const std::vector<unsigned int> & involved_dofs,
3000 unsigned int var_component,
3004 const auto & JxW = fe.
get_JxW();
3005 const auto & phi = fe.
get_phi();
3006 const std::vector<std::vector<typename FEGenericBase<typename FFunctor::RealType>::OutputGradient>> * dphi =
nullptr;
3007 const std::vector<Point> & xyz_values = fe.
get_xyz();
3009 const std::unordered_map<dof_id_type, typename FFunctor::ValuePushType> &
ids_to_save =
3010 this->projector.ids_to_save;
3015 const unsigned int n_involved_dofs =
3016 cast_int<unsigned int>(involved_dofs.size());
3018 std::vector<dof_id_type> free_dof_ids;
3020 std::vector<char> dof_is_fixed(n_involved_dofs,
false);
3024 const dof_id_type id = dof_indices_var[involved_dofs[i]];
3027 free_dof_ids.push_back(
id);
3030 dof_is_fixed[i] =
true;
3031 Uinvolved(i) = iter->second;
3035 const unsigned int free_dofs = free_dof_ids.
size();
3050 const unsigned int n_qp =
3051 cast_int<unsigned int>(xyz_values.size());
3054 for (
unsigned int qp=0; qp<n_qp; qp++)
3057 FValue fineval = f.eval_at_point(context,
3063 typename GFunctor::FunctorValue finegrad;
3065 finegrad = g->eval_at_point(context,
3072 for (
unsigned int sidei=0, freei=0;
3073 sidei != n_involved_dofs; ++sidei)
3075 unsigned int i = involved_dofs[sidei];
3077 if (dof_is_fixed[sidei])
3079 for (
unsigned int sidej=0, freej=0;
3080 sidej != n_involved_dofs; ++sidej)
3082 unsigned int j = involved_dofs[sidej];
3083 if (dof_is_fixed[sidej])
3084 Fe(freei) -= phi[i][qp] * phi[j][qp] *
3085 JxW[qp] * Uinvolved(sidej);
3087 Ke(freei,freej) += phi[i][qp] *
3088 phi[j][qp] * JxW[qp];
3091 if (dof_is_fixed[sidej])
3094 JxW[qp] * Uinvolved(sidej);
3100 if (!dof_is_fixed[sidej])
3103 Fe(freei) += phi[i][qp] * fineval * JxW[qp];
3117 insert_ids(free_dof_ids, Ufree.
get_values(), pid);
3123 #endif // GENERIC_PROJECTOR_H The GenericProjector class implements the core of other projection operations, using two input functo...
static unsigned int n_dofs_per_elem(const unsigned int dim, const FEType &fe_t, const ElemType t)
std::unordered_map< dof_id_type, std::pair< typename FFunctor::ValuePushType, processor_id_type > > new_ids_to_push
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
void send_and_insert_dof_values(std::unordered_map< dof_id_type, std::pair< Value, processor_id_type >> &ids_to_push, ProjectionAction &action) const
virtual bool is_vertex_on_parent(unsigned int c, unsigned int n) const
FEFamily family
The type of finite element.
void check_old_context(const FEMContext &c)
Real time
For time-dependent problems, this is the time t at the beginning of the current timestep.
OldSolutionBase(const libMesh::System &sys_in, const std::vector< unsigned int > *vars)
RefinementState refinement_flag() const
dof_id_type dof_number(const unsigned int s, const unsigned int var, const unsigned int comp) const
const Elem * parent() const
unsigned int variable_scalar_number(std::string_view var, unsigned int component) const
void eval_old_dofs(const Elem &elem, const FEType &fe_type, unsigned int sys_num, unsigned int var_num, std::vector< dof_id_type > &indices, std::vector< DofValueType > &values)
SubFunctor(GenericProjector &p)
OldSolutionValue(const OldSolutionValue &in)
A Node is like a Point, but with more information.
const Variable & variable(unsigned int var) const
Return a constant reference to Variable var.
bool check_old_context(const FEMContext &c, const Point &p)
const NumericVector< Number > & old_solution
void eval_old_dofs(const Elem &elem, unsigned int node_num, unsigned int var_num, std::vector< dof_id_type > &indices, std::vector< DofValueType > &values)
unsigned int get_node_index(const Node *node_ptr) const
unsigned int n_comp(const unsigned int s, const unsigned int var) const
void operator()(const node_range &range)
void construct_projection(const std::vector< dof_id_type > &dof_indices_var, const std::vector< unsigned int > &involved_dofs, unsigned int var_component, const Node *node, const FEGenericBase< typename FFunctor::RealType > &fe)
std::unique_ptr< GFunctor > master_g_deepcopy
Needed for C1 type elements only.
virtual bool is_face(const unsigned int i) const =0
NodesToElemMap * nodes_to_elem
void init_context(FEMContext &c)
~GenericProjector()=default
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
static constexpr Real TOLERANCE
const Elem & get_elem() const
Accessor for current Elem object.
static Point inverse_map(const unsigned int dim, const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true, const bool extra_checks=true)
ProjectInteriors(GenericProjector &p)
ProjectSides(ProjectSides &p_s, Threads::split)
void libmesh_merge_move(T &target, T &source)
The OldSolutionBase input functor abstract base class is the root of the OldSolutionValue and OldSolu...
virtual bool is_child_on_side(const unsigned int c, const unsigned int s) const =0
void eval_mixed_derivatives(const FEMContext &, unsigned int, unsigned int, const Node &, std::vector< Output > &)
void eval_mixed_derivatives(const FEMContext &, unsigned int, unsigned int, const Node &, std::vector< Output > &)
Dummy "splitting object" used to distinguish splitting constructors from copy constructors.
RefinementState p_refinement_flag() const
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...
ProjectVertices(GenericProjector &p)
std::vector< T > & get_values()
FEMFunctionWrapper(const FEMFunctionWrapper< Output > &fw)
ProjectSides(GenericProjector &p)
std::unordered_multimap< const Node *, std::tuple< const Elem *, unsigned short, var_set > > ves_multimap
ProjectionAction & master_action
This is the base class from which all geometric element types are derived.
FEMFunctionWrapper(const FEMFunctionBase< Output > &f)
static FEFieldType field_type(const FEType &fe_type)
RefinementState
Enumeration of possible element refinement states.
virtual bool is_node_on_side(const unsigned int n, const unsigned int s) const =0
const Parallel::Communicator & comm() const
SortAndCopy(GenericProjector &p)
virtual unsigned int n_children() const =0
unsigned int p_level() const
OrderWrapper order
The approximation order of the element.
void set_algebraic_type(const AlgebraicType atype)
Setting which determines whether to initialize algebraic structures (elem_*) on each element and set ...
ProjectVertices(ProjectVertices &p_v, Threads::split)
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.
The OldSolutionValue input functor class can be used with GenericProjector to read values from a solu...
void eval_old_dofs(const Elem &, const FEType &, unsigned int, unsigned int, std::vector< dof_id_type > &, std::vector< Output > &)
virtual Real hmax() const
NumericVector< Val > & target_vector
VectorSetAction(NumericVector< Val > &target_vec)
const MeshBase & get_mesh() const
void insert_id(dof_id_type id, const InsertInput &val, processor_id_type pid)
void operator()(const interior_range &range)
void join(const SubFunctor &other)
OldSolutionBase(const OldSolutionBase &in)
void iota(ForwardIter first, ForwardIter last, T value)
Utility::iota was created back when std::iota was just an SGI STL extension.
virtual void request_phi() const =0
request phi calculations
uint8_t processor_id_type
This is the MeshBase class.
SimpleRange< ChildRefIter > child_ref_range()
Returns a range with all children of a parent element, usable in range-based for loops.
virtual FEContinuity get_continuity() const =0
bool should_p_refine_var(unsigned int var) const
Whether the given variable should be p-refined.
unsigned int var_group_from_var_number(unsigned int var_num) const
FEType get_fe_type() const
std::vector< const Elem * > interiors
static void get_shape_outputs(FEAbstract &fe)
This class handles the numbering of degrees of freedom on a mesh.
std::unique_ptr< GFunctor > g
std::unique_ptr< QBase > default_quadrature_rule(const unsigned int dim, const int extraorder=0) const
const std::vector< std::vector< OutputGradient > > & get_dphi() const
bool is_grid_projection()
void eval_old_dofs(const Elem &, unsigned int, unsigned int, std::vector< dof_id_type > &, std::vector< Output > &)
const dof_id_type n_nodes
unsigned int number() const
std::unordered_map< dof_id_type, typename FFunctor::ValuePushType > ids_to_save
This class defines the notion of a variable in the system.
dof_id_type numeric_index_type
OldSolutionValue(const libMesh::System &sys_in, const NumericVector< Number > &old_sol, const std::vector< unsigned int > *vars)
virtual Real hmin() const
void project(const ConstElemRange &range)
Function definitions.
static bool extra_hanging_dofs(const FEType &fe_t)
TensorTools::MakeBaseNumber< Output >::type DofValueType
virtual unsigned int n_nodes() const =0
void insert_ids(const std::vector< dof_id_type > &ids, const std::vector< InsertInput > &vals, processor_id_type pid)
bool is_sorted(const std::vector< KeyType > &v)
void operator()(const ConstElemRange &range)
static const processor_id_type invalid_processor_id
An invalid processor_id to distinguish DoFs that have not been assigned to a processor.
unsigned int n_vars(const unsigned int s, const unsigned int vg) const
StoredRange< std::vector< node_projection >::const_iterator, node_projection > node_range
Manages consistently variables, degrees of freedom, and coefficient vectors.
unsigned int n_systems() const
std::unordered_map< dof_id_type, std::vector< dof_id_type > > NodesToElemMap
Convenience typedef for the Node-to-attached-Elem mapping that may be passed in to the constructor...
SimpleRange< IndexType > as_range(const std::pair< IndexType, IndexType > &p)
Helper function that allows us to treat a homogenous pair as a range.
void insert(dof_id_type id, Val val)
ProjectEdges(ProjectEdges &p_e, Threads::split)
std::vector< FEFieldType > field_types
SubProjector(GenericProjector &p)
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
This class provides all data required for a physics package (e.g.
void init_context(FEMContext &c)
std::unique_ptr< FEMFunctionBase< Output > > _f
virtual_for_inffe const std::vector< Point > & get_xyz() const
DofObject * get_old_dof_object()
Pointer accessor for previously public old_dof_object.
bool active_on_subdomain(subdomain_id_type sid) const
void insert(const std::vector< dof_id_type > &dof_indices, const DenseVector< Val > &Ue)
const std::vector< unsigned int > & variables
void eval_mixed_derivatives(const FEMContext &libmesh_dbg_var(c), unsigned int i, unsigned int dim, const Node &n, std::vector< Output > &derivs)
unsigned int which_child_am_i(const Elem *e) const
DofObject & get_old_dof_object_ref()
As above, but do not use in situations where the old_dof_object may be nullptr, since this function a...
Output eval_at_point(const FEMContext &c, unsigned int i, const Point &n, const Real time, bool)
virtual void request_dphi() const =0
request dphi calculations
DofValueType ValuePushType
static unsigned int n_dofs_at_node(const unsigned int dim, const FEType &fe_t, const ElemType t, const unsigned int n)
void operator()(const node_range &range)
TensorTools::MakeReal< Output >::type RealType
virtual numeric_index_type first_local_index() const =0
GenericProjector(const System &system_in, FFunctor &f_in, GFunctor *g_in, ProjectionAction &act_in, const std::vector< unsigned int > &variables_in, NodesToElemMap *nodes_to_elem_in=nullptr)
void find_dofs_to_send(const Node &node, const Elem &elem, unsigned short node_num, const var_set &vars)
GenericProjector & projector
virtual unsigned int n_vertices() const =0
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static const Real out_of_elem_tol
void convert_from_receive(SendT &received, T &converted)
subdomain_id_type subdomain_id() const
DofValueType ValuePushType
virtual unsigned short dim() const =0
ProjectInteriors(ProjectInteriors &p_i, Threads::split)
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
TensorTools::MakeReal< Output >::type RealType
For ease of communication, we allow users to translate their own value types to a more easily computa...
std::unordered_map< dof_id_type, typename FFunctor::ValuePushType > new_ids_to_save
virtual const Elem & elem_ref(const dof_id_type i) const
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...
TensorTools::MakeBaseNumber< Output >::type DofValueType
The DofObject defines an abstract base class for objects that have degrees of freedom associated with...
std::pair< const Node *, std::tuple< const Elem *, unsigned short, var_set > > node_projection
std::set< unsigned int > var_set
ProjectEdges(GenericProjector &p)
void parallel_reduce(const Range &range, Body &body)
Execute the provided reduction operation in parallel on the specified range.
virtual unsigned int size() const override final
std::pair< unsigned int, unsigned int > var_to_vg_and_offset(const unsigned int s, const unsigned int var) const
SortAndCopy(SortAndCopy &other, Threads::split)
virtual void set(const numeric_index_type i, const T value)=0
Sets v(i) = value.
std::vector< unsigned int > component_to_var
Defines a dense vector for use in Finite Element-type computations.
void cholesky_solve(const DenseVector< T2 > &b, DenseVector< T2 > &x)
For symmetric positive definite (SPD) matrices.
The VectorSetAction output functor class can be used with a GenericProjector to set projection values...
StoredRange< std::vector< const Elem * >::const_iterator, const Elem * > interior_range
unsigned int n_vars() const
This class forms the foundation from which generic finite elements may be derived.
const TypeToSend< T >::type convert_to_send(const T &in)
unsigned int n_comp_group(const unsigned int s, const unsigned int vg) const
processor_id_type processor_id() const
void operator()(const node_range &range)
Output eval_at_node(const FEMContext &c, unsigned int i, unsigned int elem_dim, const Node &n, bool, Real=0.)
FEMFunctionBase is a base class from which users can derive in order to define "function-like" object...
const DofMap & get_dof_map() const
processor_id_type processor_id() const
Output eval_at_node(const FEMContext &c, unsigned int i, unsigned int, const Node &n, bool, const Real time)
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...
The FEMFunctionWrapper input functor class can be used with a GenericProjector to read values from an...
A Point defines a location in LIBMESH_DIM dimensional Real space.
const Point & point(const unsigned int i) const
void join(const SortAndCopy &other)
virtual numeric_index_type last_local_index() const =0
void ErrorVector unsigned int
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
virtual bool is_child_on_edge(const unsigned int c, const unsigned int e) const
const Elem & get(const ElemType type_in)
This class forms the foundation from which generic finite elements may be derived.
std::vector< FEContinuity > conts
NodesToElemMap nodes_to_elem_ourcopy
nodes_to_elem is either a shallow copy of a map passed in to the constructor, or points to nodes_to_e...
GenericProjector(const GenericProjector &in)
bool is_grid_projection()
const Elem * child_ptr(unsigned int i) const
This class defines a tensor in LIBMESH_DIM dimensional Real or Complex space.
Output eval_at_point(const FEMContext &c, unsigned int i, const Point &p, Real, bool skip_context_check)
const std::vector< std::vector< OutputShape > > & get_phi() const
const FEType & type() const
spin_mutex spin_mtx
A convenient spin mutex object which can be used for obtaining locks.
FEFieldType
defines an enum for finite element field types - i.e.