This method projects an arbitrary boundary solution to the current mesh. The input function f
gives the arbitrary solution, while the new_vector
(which should already be correctly sized) gives the solution (to be computed) on the current mesh.
1492 const BoundaryInfo & boundary_info =
1500 DenseVector<Number> Fe;
1502 DenseVector<Number> Ue;
1510 const Variable & variable = dof_map.variable(var);
1512 const FEType & fe_type = variable.type();
1514 if (fe_type.family ==
SCALAR)
1517 const unsigned int var_component =
1524 std::unique_ptr<QBase> qedgerule (fe_type.default_quadrature_rule(1));
1525 std::unique_ptr<QBase> qsiderule (fe_type.default_quadrature_rule(
dim-1));
1529 const std::vector<std::vector<Real>> & phi = fe->get_phi();
1533 const std::vector<std::vector<RealGradient>> * dphi =
nullptr;
1542 const std::vector<std::vector<RealGradient>> &
1543 ref_dphi = fe->get_dphi();
1548 const std::vector<Real> & JxW =
1552 const std::vector<Point> & xyz_values =
1556 std::vector<dof_id_type> dof_indices;
1558 std::vector<unsigned int> side_dofs;
1561 std::vector<boundary_id_type> bc_ids;
1564 for (
const auto & elem : range)
1568 if (!variable.active_on_subdomain(elem->subdomain_id()))
1571 const unsigned short n_nodes = elem->n_nodes();
1572 const unsigned short n_edges = elem->n_edges();
1573 const unsigned short n_sides = elem->n_sides();
1577 std::vector<bool> is_boundary_node(
n_nodes,
false),
1578 is_boundary_edge(n_edges,
false),
1579 is_boundary_side(n_sides,
false);
1582 std::vector<bool> is_boundary_nodeset(
n_nodes,
false);
1584 for (
unsigned char s=0; s != n_sides; ++s)
1587 boundary_info.boundary_ids (elem, s, bc_ids);
1588 bool do_this_side =
false;
1589 for (
const auto & bc_id : bc_ids)
1592 do_this_side =
true;
1598 is_boundary_side[s] =
true;
1601 for (
unsigned int n=0; n !=
n_nodes; ++n)
1602 if (elem->is_node_on_side(n,s))
1603 is_boundary_node[n] =
true;
1604 for (
unsigned int e=0; e != n_edges; ++e)
1605 if (elem->is_edge_on_side(e,s))
1606 is_boundary_edge[e] =
true;
1611 for (
unsigned int n=0; n !=
n_nodes; ++n)
1613 boundary_info.boundary_ids (elem->node_ptr(n), bc_ids);
1615 for (
const auto & bc_id : bc_ids)
1618 is_boundary_node[n] =
true;
1619 is_boundary_nodeset[n] =
true;
1625 for (
unsigned short e=0; e != n_edges; ++e)
1627 boundary_info.edge_boundary_ids (elem, e, bc_ids);
1629 for (
const auto & bc_id : bc_ids)
1631 is_boundary_edge[e] =
true;
1636 dof_map.dof_indices (elem, dof_indices, var);
1639 const unsigned int n_dofs =
1640 cast_int<unsigned int>(dof_indices.size());
1643 std::vector<char> dof_is_fixed(n_dofs,
false);
1644 std::vector<int> free_dof(n_dofs, 0);
1647 Ue.resize (n_dofs); Ue.zero();
1656 unsigned int current_dof = 0;
1657 for (
unsigned short n = 0; n !=
n_nodes; ++n)
1663 const unsigned int nc =
1666 if ((!elem->is_vertex(n) || !is_boundary_node[n]) &&
1667 !is_boundary_nodeset[n])
1674 libmesh_assert_equal_to (nc, 0);
1680 libmesh_assert_equal_to (nc, 1);
1681 Ue(current_dof) =
f->component(var_component,
1684 dof_is_fixed[current_dof] =
true;
1688 else if (fe_type.family ==
HERMITE)
1690 Ue(current_dof) =
f->component(var_component,
1693 dof_is_fixed[current_dof] =
true;
1695 Gradient grad =
g->component(var_component,
1699 Ue(current_dof) = grad(0);
1700 dof_is_fixed[current_dof] =
true;
1706 Point nxminus = elem->point(n),
1707 nxplus = elem->point(n);
1710 Gradient gxminus =
g->component(var_component,
1713 Gradient gxplus =
g->component(var_component,
1717 Ue(current_dof) = grad(1);
1718 dof_is_fixed[current_dof] =
true;
1721 Ue(current_dof) = (gxplus(1) - gxminus(1))
1723 dof_is_fixed[current_dof] =
true;
1730 Ue(current_dof) = grad(2);
1731 dof_is_fixed[current_dof] =
true;
1734 Ue(current_dof) = (gxplus(2) - gxminus(2))
1736 dof_is_fixed[current_dof] =
true;
1739 Point nyminus = elem->point(n),
1740 nyplus = elem->point(n);
1743 Gradient gyminus =
g->component(var_component,
1746 Gradient gyplus =
g->component(var_component,
1750 Ue(current_dof) = (gyplus(2) - gyminus(2))
1752 dof_is_fixed[current_dof] =
true;
1755 Point nxmym = elem->point(n),
1756 nxmyp = elem->point(n),
1757 nxpym = elem->point(n),
1758 nxpyp = elem->point(n);
1767 Gradient gxmym =
g->component(var_component,
1770 Gradient gxmyp =
g->component(var_component,
1773 Gradient gxpym =
g->component(var_component,
1776 Gradient gxpyp =
g->component(var_component,
1779 Number gxzplus = (gxpyp(2) - gxmyp(2))
1781 Number gxzminus = (gxpym(2) - gxmym(2))
1784 Ue(current_dof) = (gxzplus - gxzminus)
1786 dof_is_fixed[current_dof] =
true;
1789 #endif // LIBMESH_DIM > 2 1791 #endif // LIBMESH_DIM > 1 1796 else if (cont ==
C_ONE)
1798 libmesh_assert_equal_to (nc, 1 +
dim);
1799 Ue(current_dof) =
f->component(var_component,
1802 dof_is_fixed[current_dof] =
true;
1804 Gradient grad =
g->component(var_component,
1807 for (
unsigned int i=0; i!=
dim; ++i)
1809 Ue(current_dof) = grad(i);
1810 dof_is_fixed[current_dof] =
true;
1815 libmesh_error_msg(
"Unknown continuity " << cont);
1820 for (
unsigned short e = 0; e != n_edges; ++e)
1822 if (!is_boundary_edge[e])
1828 const unsigned int n_side_dofs =
1829 cast_int<unsigned int>(side_dofs.size());
1833 unsigned int free_dofs = 0;
1835 if (!dof_is_fixed[side_dofs[i]])
1836 free_dof[free_dofs++] = i;
1842 Ke.resize (free_dofs, free_dofs); Ke.zero();
1843 Fe.resize (free_dofs); Fe.zero();
1845 DenseVector<Number> Uedge(free_dofs);
1848 fe->attach_quadrature_rule (qedgerule.get());
1849 fe->edge_reinit (elem, e);
1850 const unsigned int n_qp = qedgerule->n_points();
1853 for (
unsigned int qp=0; qp<n_qp; qp++)
1856 Number fineval =
f->component(var_component,
1862 finegrad =
g->component(var_component,
1867 for (
unsigned int sidei=0, freei=0;
1868 sidei != n_side_dofs; ++sidei)
1870 unsigned int i = side_dofs[sidei];
1872 if (dof_is_fixed[i])
1874 for (
unsigned int sidej=0, freej=0;
1875 sidej != n_side_dofs; ++sidej)
1877 unsigned int j = side_dofs[sidej];
1878 if (dof_is_fixed[j])
1879 Fe(freei) -= phi[i][qp] * phi[j][qp] *
1882 Ke(freei,freej) += phi[i][qp] *
1883 phi[j][qp] * JxW[qp];
1886 if (dof_is_fixed[j])
1887 Fe(freei) -= ((*dphi)[i][qp] *
1891 Ke(freei,freej) += ((*dphi)[i][qp] *
1895 if (!dof_is_fixed[j])
1898 Fe(freei) += phi[i][qp] * fineval * JxW[qp];
1900 Fe(freei) += (finegrad * (*dphi)[i][qp]) *
1906 Ke.cholesky_solve(Fe, Uedge);
1909 for (
unsigned int i=0; i != free_dofs; ++i)
1911 Number & ui = Ue(side_dofs[free_dof[i]]);
1915 dof_is_fixed[side_dofs[free_dof[i]]] =
true;
1921 for (
unsigned short s = 0; s != n_sides; ++s)
1923 if (!is_boundary_side[s])
1931 unsigned int free_dofs = 0;
1933 if (!dof_is_fixed[side_dofs[i]])
1934 free_dof[free_dofs++] = i;
1940 Ke.resize (free_dofs, free_dofs); Ke.zero();
1941 Fe.resize (free_dofs); Fe.zero();
1943 DenseVector<Number> Uside(free_dofs);
1946 fe->attach_quadrature_rule (qsiderule.get());
1947 fe->reinit (elem, s);
1948 const unsigned int n_qp = qsiderule->n_points();
1950 const unsigned int n_side_dofs =
1951 cast_int<unsigned int>(side_dofs.size());
1954 for (
unsigned int qp=0; qp<n_qp; qp++)
1957 Number fineval =
f->component(var_component,
1963 finegrad =
g->component(var_component,
1968 for (
unsigned int sidei=0, freei=0;
1969 sidei != n_side_dofs; ++sidei)
1971 unsigned int i = side_dofs[sidei];
1973 if (dof_is_fixed[i])
1975 for (
unsigned int sidej=0, freej=0;
1976 sidej != n_side_dofs; ++sidej)
1978 unsigned int j = side_dofs[sidej];
1979 if (dof_is_fixed[j])
1980 Fe(freei) -= phi[i][qp] * phi[j][qp] *
1983 Ke(freei,freej) += phi[i][qp] *
1984 phi[j][qp] * JxW[qp];
1987 if (dof_is_fixed[j])
1988 Fe(freei) -= ((*dphi)[i][qp] *
1992 Ke(freei,freej) += ((*dphi)[i][qp] *
1996 if (!dof_is_fixed[j])
1999 Fe(freei) += (fineval * phi[i][qp]) * JxW[qp];
2001 Fe(freei) += (finegrad * (*dphi)[i][qp]) *
2007 Ke.cholesky_solve(Fe, Uside);
2010 for (
unsigned int i=0; i != free_dofs; ++i)
2012 Number & ui = Ue(side_dofs[free_dof[i]]);
2016 dof_is_fixed[side_dofs[free_dof[i]]] =
true;
2028 for (
unsigned int i = 0; i < n_dofs; i++)
2029 if (dof_is_fixed[i] &&
2030 (dof_indices[i] >= first) &&
2031 (dof_indices[i] < last))
Real time
For time-dependent problems, this is the time t at the beginning of the current timestep.
unsigned int variable_scalar_number(std::string_view var, unsigned int component) const
static constexpr Real TOLERANCE
std::unique_ptr< FunctionBase< Number > > f
const std::vector< unsigned int > & variables
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
const MeshBase & get_mesh() 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, bool add_p_level=true)
Fills the vector di with the local degree of freedom indices associated with side s of element elem A...
ADRealEigenVector< T, D, asd > abs(const ADRealEigenVector< T, D, asd > &)
const dof_id_type n_nodes
NumericVector< Number > & new_vector
NumberVectorValue Gradient
static std::unique_ptr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.
template class LIBMESH_EXPORT DenseMatrix< Real >
static unsigned int n_dofs_at_node(const unsigned int dim, const FEType &fe_t, const ElemType t, const unsigned int n)
virtual numeric_index_type first_local_index() const =0
FEContinuity
defines an enum for finite element types to libmesh_assert a certain level (or type? Hcurl?) of continuity.
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...
unsigned int mesh_dimension() 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, bool add_p_level=true)
Fills the vector di with the local degree of freedom indices associated with edge e of element elem A...
virtual void set(const numeric_index_type i, const T value)=0
Sets v(i) = value.
const std::set< boundary_id_type > & b
std::unique_ptr< FunctionBase< Gradient > > g
const DofMap & get_dof_map() const
virtual numeric_index_type last_local_index() const =0
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
spin_mutex spin_mtx
A convenient spin mutex object which can be used for obtaining locks.