Newer
Older
Thomas Steinreiter
committed
#include "heatsystem.h"
#include "libmesh/dirichlet_boundaries.h"
#include "libmesh/dof_map.h"
#include "libmesh/elem.h"
#include "libmesh/fe_base.h"
#include "libmesh/fe_interface.h"
#include "libmesh/fem_context.h"
#include "libmesh/getpot.h"
#include "libmesh/mesh.h"
#include "libmesh/quadrature.h"
#include "libmesh/string_to_enum.h"
#include "libmesh/zero_function.h"
using namespace libMesh;
// user defined literal for libmesh's Real
constexpr auto operator"" _r(long double n) { return libMesh::Real(n); }
Thomas Steinreiter
committed
void HeatSystem::init_data() {
T_var = this->add_variable("T", static_cast<Order>(_fe_order),
Utility::string_to_enum<FEFamily>(_fe_family));
ZeroFunction<Number> zero;
this->get_dof_map().add_dirichlet_boundary(
Thomas Steinreiter
committed
// Do the parent's initialization after variables are defined
FEMSystem::init_data();
this->time_evolving(0);
}
void HeatSystem::init_context(DiffContext& context) {
Thomas Steinreiter
committed
Thomas Steinreiter
committed
for (const auto& dim : elem_dims) {
FEBase* fe = nullptr;
Thomas Steinreiter
committed
c.get_element_fe(T_var, fe, dim);
fe->get_JxW(); // For integration
fe->get_dphi(); // For bilinear form
fe->get_xyz(); // For forcing
fe->get_phi(); // For forcing
}
FEMSystem::init_context(context);
}
bool HeatSystem::element_time_derivative(bool request_jacobian,
DiffContext& context) {
Thomas Steinreiter
committed
const auto& mesh_dim = c.get_system().get_mesh().mesh_dimension();
Thomas Steinreiter
committed
// First we get some references to cell-specific data that
// will be used to assemble the linear system.
const auto& dim = c.get_elem().dim();
FEBase* fe = nullptr;
Thomas Steinreiter
committed
c.get_element_fe(T_var, fe, dim);
// Element Jacobian * quadrature weights for interior integration
Thomas Steinreiter
committed
Thomas Steinreiter
committed
Thomas Steinreiter
committed
Thomas Steinreiter
committed
// The number of local degrees of freedom in each variable
Thomas Steinreiter
committed
// The subvectors and submatrices we need to fill:
auto& K = c.get_elem_jacobian(T_var, T_var);
auto& F = c.get_elem_residual(T_var);
Thomas Steinreiter
committed
// Now we will build the element Jacobian and residual.
// Constructing the residual requires the solution and its
// gradient from the previous timestep. This must be
// calculated at each quadrature point by summing the
// solution degree-of-freedom values by the appropriate
// weight functions.
Thomas Steinreiter
committed
Thomas Steinreiter
committed
// Compute the solution gradient at the Newton iterate
Gradient grad_T = c.interior_gradient(T_var, qp);
Thomas Steinreiter
committed
Thomas Steinreiter
committed
// Only apply forcing to interior elements
Thomas Steinreiter
committed
(dim == mesh_dim) ? -k * u_exact * (dim * libMesh::pi * libMesh::pi)
: 0;
Thomas Steinreiter
committed
Thomas Steinreiter
committed
F(i) += JxWxNK * (grad_T * dphi[i][qp] + forcing * phi[i][qp]);
if (request_jacobian) {
Thomas Steinreiter
committed
JxWxNK * context.get_elem_solution_derivative();
for (std::size_t i{0}; i != n_T_dofs; i++)
for (std::size_t j{0}; j != n_T_dofs; ++j)
Thomas Steinreiter
committed
K(i, j) += JxWxNKxD * (dphi[i][qp] * dphi[j][qp]);
}
} // end of the quadrature point qp-loop
return request_jacobian;
}