libMesh
heatsystem.C
Go to the documentation of this file.
1 #include "heatsystem.h"
2 
3 #include "libmesh/dirichlet_boundaries.h"
4 #include "libmesh/dof_map.h"
5 #include "libmesh/fe_base.h"
6 #include "libmesh/fe_interface.h"
7 #include "libmesh/fem_context.h"
8 #include "libmesh/getpot.h"
9 #include "libmesh/mesh.h"
10 #include "libmesh/quadrature.h"
11 #include "libmesh/string_to_enum.h"
12 #include "libmesh/zero_function.h"
13 #include "libmesh/elem.h"
14 
15 using namespace libMesh;
16 
18 {
19  T_var = this->add_variable("T", static_cast<Order>(_fe_order),
20  Utility::string_to_enum<FEFamily>(_fe_family));
21 
22  const unsigned int dim = this->get_mesh().mesh_dimension();
23 
24  // Add dirichlet boundaries on all but the boundary element side
25  const boundary_id_type all_ids[6] = {0, 1, 2, 3, 4, 5};
26  std::set<boundary_id_type> nonyplus_bdys(all_ids, all_ids+(dim*2));
27  const boundary_id_type yplus_id = (dim == 3) ? 3 : 2;
28  nonyplus_bdys.erase(yplus_id);
29 
30  std::vector<unsigned int> T_only(1, T_var);
32 
33  // Most DirichletBoundary users will want to supply a "locally
34  // indexed" functor
35  this->get_dof_map().add_dirichlet_boundary
36  (DirichletBoundary (nonyplus_bdys, T_only, zero, LOCAL_VARIABLE_ORDER));
37 
38  // Do the parent's initialization after variables are defined
40 
41  // The temperature is evolving, with a first-order time derivative
42  this->time_evolving(T_var, 1);
43 }
44 
45 
46 
48 {
49  FEMContext & c = cast_ref<FEMContext &>(context);
50 
51  const std::set<unsigned char> & elem_dims =
52  c.elem_dimensions();
53 
54  for (std::set<unsigned char>::const_iterator dim_it =
55  elem_dims.begin(); dim_it != elem_dims.end(); ++dim_it)
56  {
57  const unsigned char dim = *dim_it;
58 
59  FEBase * fe = libmesh_nullptr;
60 
61  c.get_element_fe(T_var, fe, dim);
62 
63  fe->get_JxW(); // For integration
64  fe->get_dphi(); // For bilinear form
65  fe->get_xyz(); // For forcing
66  fe->get_phi(); // For forcing
67  }
68 
69  FEMSystem::init_context(context);
70 }
71 
72 
73 bool HeatSystem::element_time_derivative (bool request_jacobian,
74  DiffContext & context)
75 {
76  FEMContext & c = cast_ref<FEMContext &>(context);
77 
78  const unsigned int mesh_dim =
80 
81  // First we get some references to cell-specific data that
82  // will be used to assemble the linear system.
83  const unsigned int dim = c.get_elem().dim();
84  FEBase * fe = libmesh_nullptr;
85  c.get_element_fe(T_var, fe, dim);
86 
87  // Element Jacobian * quadrature weights for interior integration
88  const std::vector<Real> & JxW = fe->get_JxW();
89 
90  const std::vector<Point> & xyz = fe->get_xyz();
91 
92  const std::vector<std::vector<Real>> & phi = fe->get_phi();
93 
94  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
95 
96  // The number of local degrees of freedom in each variable
97  const unsigned int n_T_dofs = c.get_dof_indices(T_var).size();
98 
99  // The subvectors and submatrices we need to fill:
100  DenseSubMatrix<Number> & K = c.get_elem_jacobian(T_var, T_var);
102 
103  // Now we will build the element Jacobian and residual.
104  // Constructing the residual requires the solution and its
105  // gradient from the previous timestep. This must be
106  // calculated at each quadrature point by summing the
107  // solution degree-of-freedom values by the appropriate
108  // weight functions.
109  unsigned int n_qpoints = c.get_element_qrule().n_points();
110 
111  for (unsigned int qp=0; qp != n_qpoints; qp++)
112  {
113  // Compute the solution gradient at the Newton iterate
114  Gradient grad_T = c.interior_gradient(T_var, qp);
115 
116  const Number k = _k[dim];
117 
118  const Point & p = xyz[qp];
119 
120  // solution + laplacian depend on problem dimension
121  const Number u_exact = (mesh_dim == 2) ?
122  std::sin(libMesh::pi*p(0)) * std::sin(libMesh::pi*p(1)) :
123  std::sin(libMesh::pi*p(0)) * std::sin(libMesh::pi*p(1)) *
124  std::sin(libMesh::pi*p(2));
125 
126  // Only apply forcing to interior elements
127  const Number forcing = (dim == mesh_dim) ?
128  -k * u_exact * (dim * libMesh::pi * libMesh::pi) : 0;
129 
130  const Number JxWxNK = JxW[qp] * -k;
131 
132  for (unsigned int i=0; i != n_T_dofs; i++)
133  F(i) += JxWxNK * (grad_T * dphi[i][qp] + forcing * phi[i][qp]);
134  if (request_jacobian)
135  {
136  const Number JxWxNKxD = JxWxNK *
138 
139  for (unsigned int i=0; i != n_T_dofs; i++)
140  for (unsigned int j=0; j != n_T_dofs; ++j)
141  K(i,j) += JxWxNKxD * (dphi[i][qp] * dphi[j][qp]);
142  }
143  } // end of the quadrature point qp-loop
144 
145  return request_jacobian;
146 }
This class provides all data required for a physics package (e.g.
Definition: diff_context.h:54
ConstFunction that simply returns 0.
Definition: zero_function.h:35
unsigned int dim
virtual bool element_time_derivative(bool request_jacobian, DiffContext &context)
Adds the time derivative contribution on elem to elem_residual.
Definition: heatsystem.C:110
const std::vector< dof_id_type > & get_dof_indices() const
Accessor for element dof indices.
Definition: diff_context.h:366
const class libmesh_nullptr_t libmesh_nullptr
This class allows one to associate Dirichlet boundary values with a given set of mesh boundary ids an...
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...
Definition: fem_context.h:262
The libMesh namespace provides an interface to certain functionality in the library.
const Number zero
.
Definition: libmesh.h:178
const DenseVector< Number > & get_elem_residual() const
Const accessor for element residual.
Definition: diff_context.h:248
const Elem & get_elem() const
Accessor for current Elem object.
Definition: fem_context.h:871
Defines a dense subvector for use in finite element computations.
const DenseMatrix< Number > & get_elem_jacobian() const
Const accessor for element Jacobian.
Definition: diff_context.h:282
const std::set< unsigned char > & elem_dimensions() const
Definition: fem_context.h:913
const MeshBase & get_mesh() const
Definition: system.h:2014
int8_t boundary_id_type
Definition: id_types.h:51
const std::vector< Real > & get_JxW() const
Definition: fe_abstract.h:237
virtual void init_context(DiffContext &context)
Definition: heatsystem.C:77
Defines a dense submatrix for use in Finite Element-type computations.
const std::vector< std::vector< OutputShape > > & get_phi() const
Definition: fe_base.h:208
This class provides all data required for a physics package (e.g.
Definition: fem_context.h:61
Gradient interior_gradient(unsigned int var, unsigned int qp) const
Definition: fem_context.C:381
const System & get_system() const
Accessor for associated system.
Definition: diff_context.h:104
const std::vector< std::vector< OutputGradient > > & get_dphi() const
Definition: fe_base.h:216
Real get_elem_solution_derivative() const
The derivative of the current elem_solution w.r.t.
Definition: diff_context.h:419
unsigned int mesh_dimension() const
Definition: mesh_base.C:148
const QBase & get_element_qrule() const
Accessor for element interior quadrature rule for the dimension of the current _elem.
Definition: fem_context.h:765
virtual unsigned int dim() const =0
virtual void init_data()
Initializes the member data fields associated with the system, so that, e.g., assemble() may be used...
Definition: heatsystem.C:36
virtual void init_data() libmesh_override
Initializes the member data fields associated with the system, so that, e.g., assemble() may be used...
Definition: fem_system.C:845
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
unsigned int n_points() const
Definition: quadrature.h:113
sys get_dof_map()
const std::vector< Point > & get_xyz() const
Definition: fe_abstract.h:230
This class forms the foundation from which generic finite elements may be derived.
const Real pi
.
Definition: libmesh.h:172
virtual void init_context(DiffContext &) libmesh_override
Definition: fem_system.C:1345