heatsystem.cpp 3.45 KB
Newer Older
Thomas Steinreiter's avatar
Thomas Steinreiter committed
1
#include "heatsystem.hpp"
2
3
4
5
6
7
8
9
10
11
12
13
14
15

#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;
Thomas Steinreiter's avatar
Thomas Steinreiter committed
16
17
// user defined literal for libmesh's Real
constexpr auto operator"" _r(long double n) { return libMesh::Real(n); }
18
19
20
21
22
23
24
25

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's avatar
Thomas Steinreiter committed
26
	    DirichletBoundary({0, 1, 2}, {T_var}, &zero));
27
28
29
30
31
32
33
34

	// Do the parent's initialization after variables are defined
	FEMSystem::init_data();

	this->time_evolving(0);
}

void HeatSystem::init_context(DiffContext& context) {
Thomas Steinreiter's avatar
Thomas Steinreiter committed
35
	auto& c = libmesh_cast_ref<FEMContext&>(context);
36

Thomas Steinreiter's avatar
Thomas Steinreiter committed
37
	const auto& elem_dims = c.elem_dimensions();
38

Thomas Steinreiter's avatar
Thomas Steinreiter committed
39
40
	for (const auto& dim : elem_dims) {
		FEBase* fe = nullptr;
41
42
43
44
45
46
47
48
49
50
51
52
53
		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's avatar
Thomas Steinreiter committed
54
	auto& c = libmesh_cast_ref<FEMContext&>(context);
55

Thomas Steinreiter's avatar
Thomas Steinreiter committed
56
	const auto& mesh_dim = c.get_system().get_mesh().mesh_dimension();
57
58
59

	// First we get some references to cell-specific data that
	// will be used to assemble the linear system.
Thomas Steinreiter's avatar
Thomas Steinreiter committed
60
61
	const auto& dim = c.get_elem().dim();
	FEBase* fe = nullptr;
62
63
64
	c.get_element_fe(T_var, fe, dim);

	// Element Jacobian * quadrature weights for interior integration
Thomas Steinreiter's avatar
Thomas Steinreiter committed
65
	const auto& JxW = fe->get_JxW();
66

Thomas Steinreiter's avatar
Thomas Steinreiter committed
67
	const auto& xyz = fe->get_xyz();
68

Thomas Steinreiter's avatar
Thomas Steinreiter committed
69
	const auto& phi = fe->get_phi();
70

Thomas Steinreiter's avatar
Thomas Steinreiter committed
71
	const auto& dphi = fe->get_dphi();
72
73

	// The number of local degrees of freedom in each variable
Thomas Steinreiter's avatar
Thomas Steinreiter committed
74
	const auto& n_T_dofs = c.get_dof_indices(T_var).size();
75
76

	// The subvectors and submatrices we need to fill:
Thomas Steinreiter's avatar
Thomas Steinreiter committed
77
78
	auto& K = c.get_elem_jacobian(T_var, T_var);
	auto& F = c.get_elem_residual(T_var);
79
80
81
82
83
84
85

	// 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's avatar
Thomas Steinreiter committed
86
	auto n_qpoints = c.get_element_qrule().n_points();
87

Thomas Steinreiter's avatar
Thomas Steinreiter committed
88
	for (std::size_t qp{0}; qp != n_qpoints; qp++) {
89
90
91
		// Compute the solution gradient at the Newton iterate
		Gradient grad_T = c.interior_gradient(T_var, qp);

Thomas Steinreiter's avatar
Thomas Steinreiter committed
92
		const auto& k = _k[dim];
93

Thomas Steinreiter's avatar
Thomas Steinreiter committed
94
		const auto& u_exact = +1.0_r;
95
96

		// Only apply forcing to interior elements
Thomas Steinreiter's avatar
Thomas Steinreiter committed
97
		const auto& forcing =
98
99
100
		    (dim == mesh_dim) ? -k * u_exact * (dim * libMesh::pi * libMesh::pi)
		                      : 0;

Thomas Steinreiter's avatar
Thomas Steinreiter committed
101
		const auto& JxWxNK = JxW[qp] * -k;
102

Thomas Steinreiter's avatar
Thomas Steinreiter committed
103
		for (std::size_t i{0}; i != n_T_dofs; i++)
104
105
			F(i) += JxWxNK * (grad_T * dphi[i][qp] + forcing * phi[i][qp]);
		if (request_jacobian) {
Thomas Steinreiter's avatar
Thomas Steinreiter committed
106
			const auto& JxWxNKxD =
107
108
			    JxWxNK * context.get_elem_solution_derivative();

Thomas Steinreiter's avatar
Thomas Steinreiter committed
109
110
			for (std::size_t i{0}; i != n_T_dofs; i++)
				for (std::size_t j{0}; j != n_T_dofs; ++j)
111
112
113
114
115
116
					K(i, j) += JxWxNKxD * (dphi[i][qp] * dphi[j][qp]);
		}
	} // end of the quadrature point qp-loop

	return request_jacobian;
}