libMesh
poisson.C
Go to the documentation of this file.
1 /* This is where we define the assembly of the Laplace system */
2 
3 // General libMesh includes
4 #include "libmesh/getpot.h"
5 #include "libmesh/boundary_info.h"
6 #include "libmesh/dirichlet_boundaries.h"
7 #include "libmesh/dof_map.h"
8 #include "libmesh/fe_base.h"
9 #include "libmesh/fe_interface.h"
10 #include "libmesh/fem_context.h"
11 #include "libmesh/fem_system.h"
12 #include "libmesh/quadrature.h"
13 #include "libmesh/string_to_enum.h"
14 #include "libmesh/mesh.h"
15 #include "libmesh/parallel.h"
16 #include "libmesh/quadrature.h"
17 #include "libmesh/fem_context.h"
18 #include "libmesh/zero_function.h"
19 
20 // Local includes
21 #include "poisson.h"
22 
23 // Bring in everything from the libMesh namespace
24 using namespace libMesh;
25 
26 // Function to set the Dirichlet boundary function for the adjoint dirichlet boundary
27 class BdyFunction : public FunctionBase<Number>
28 {
29 public:
30  BdyFunction (unsigned int T_var)
31  : _T_var(T_var)
32  { this->_initialized = true; }
33 
34  virtual Number operator() (const Point &, const Real = 0)
35  { libmesh_not_implemented(); }
36 
37  virtual void operator() (const Point & p,
38  const Real,
39  DenseVector<Number> & output)
40  {
41  output.resize(2);
42  output.zero();
43  const Real x=p(0);
44  // Set the parabolic weighting
45  output(_T_var) = -x * (1 - x);
46  }
47 
49  { return UniquePtr<FunctionBase<Number>> (new BdyFunction(_T_var)); }
50 
51 private:
52  const unsigned int _T_var;
53 };
54 
56 {
57  T_var =
58  this->add_variable ("T", static_cast<Order>(_fe_order),
59  Utility::string_to_enum<FEFamily>(_fe_family));
60 
61  GetPot infile("poisson.in");
62  exact_QoI[0] = infile("QoI_0", 0.0);
63  alpha = infile("alpha", 100.0);
64 
65  // The temperature is evolving, with a first order time derivative
66  this->time_evolving(T_var, 1);
67 
68  // Now we will set the Dirichlet boundary conditions
69 
70  // Get boundary ids for all the boundaries
71  const boundary_id_type all_bdry_id[4] = {0, 1, 2, 3};
72  std::set<boundary_id_type> all_bdy(all_bdry_id, all_bdry_id+4);
73 
74  // For the adjoint problem, we will only set the bottom boundary to a non-zero value
75  const boundary_id_type bottom_bdry_id = 0;
76  std::set<boundary_id_type> bottom_bdry;
77  bottom_bdry.insert(bottom_bdry_id);
78 
79  // The T_only identifier for setting the boundary conditions for T
80  std::vector<unsigned int> T_only(1, T_var);
81 
82  // The zero function pointer for the primal all bdry bcs
84  // Boundary function for bottom bdry adjoint condition
85  BdyFunction bottom_adjoint(T_var);
86 
87  this->get_dof_map().add_dirichlet_boundary(DirichletBoundary (all_bdy, T_only, &zero));
88 
89  this->get_dof_map().add_adjoint_dirichlet_boundary(DirichletBoundary (bottom_bdry, T_only, &bottom_adjoint), 0);
90 
91  // Do the parent's initialization after variables are defined
93 }
94 
96 {
97  FEMContext & c = cast_ref<FEMContext &>(context);
98 
99  // Now make sure we have requested all the data
100  // we need to build the linear system.
101  FEBase* elem_fe = NULL;
102  c.get_element_fe( 0, elem_fe );
103  elem_fe->get_JxW();
104  elem_fe->get_phi();
105  elem_fe->get_dphi();
106  elem_fe->get_xyz();
107 
108  FEBase* side_fe = NULL;
109  c.get_side_fe( 0, side_fe );
110 
111  side_fe->get_JxW();
112  side_fe->get_phi();
113  side_fe->get_dphi();
114  side_fe->get_xyz();
115 
116 }
117 
118 #define optassert(X) {if (!(X)) libmesh_error_msg("Assertion " #X " failed.");}
119 
120 // Assemble the element contributions to the stiffness matrix
121 bool PoissonSystem::element_time_derivative (bool request_jacobian,
122  DiffContext & context)
123 {
124  // Are the jacobians specified analytically ?
125  bool compute_jacobian = request_jacobian && _analytic_jacobians;
126 
127  FEMContext & c = cast_ref<FEMContext &>(context);
128 
129  // First we get some references to cell-specific data that
130  // will be used to assemble the linear system.
131  FEBase* elem_fe = NULL;
132  c.get_element_fe( 0, elem_fe );
133 
134  // Element Jacobian * quadrature weights for interior integration
135  const std::vector<Real> & JxW = elem_fe->get_JxW();
136 
137  // Element basis functions
138  const std::vector<std::vector<Real>> & phi = elem_fe->get_phi();
139  const std::vector<std::vector<RealGradient>> & dphi = elem_fe->get_dphi();
140 
141  // Quadrature point locations
142  const std::vector<Point > & q_point = elem_fe->get_xyz();
143 
144  // The number of local degrees of freedom in each variable
145  const unsigned int n_T_dofs = c.get_dof_indices(0).size();
146 
147  // The subvectors and submatrices we need to fill:
150 
151  // Now we will build the element Jacobian and residual.
152  // Constructing the residual requires the solution and its
153  // gradient from the previous timestep. This must be
154  // calculated at each quadrature point by summing the
155  // solution degree-of-freedom values by the appropriate
156  // weight functions.
157  unsigned int n_qpoints = c.get_element_qrule().n_points();
158 
159  for (unsigned int qp=0; qp != n_qpoints; qp++)
160  {
161  // Compute the solution gradient at the Newton iterate
162  Gradient grad_T = c.interior_gradient(0, qp);
163 
164  // Location of the current qp
165  const Real x = q_point[qp](0);
166  const Real y = q_point[qp](1);
167 
168  // Forcing function
169  Real f = -alpha * ( ( (- 4 * alpha * alpha) * exp(-alpha*x) * y * (1 - y) ) + ( -8 + ( 8 * exp(-alpha*x) ) + ( 8 * ( 1 - exp(-alpha) )* x) ) );
170 
171  // The residual contribution from this element
172  for (unsigned int i=0; i != n_T_dofs; i++)
173  F(i) += JxW[qp] * ( (f * phi[i][qp]) - alpha*(grad_T * dphi[i][qp]) ) ;
174  if (compute_jacobian)
175  for (unsigned int i=0; i != n_T_dofs; i++)
176  for (unsigned int j=0; j != n_T_dofs; ++j)
177  // The analytic jacobian
178  K(i,j) += JxW[qp] * ( -alpha*(dphi[i][qp] * dphi[j][qp]) );
179  } // end of the quadrature point qp-loop
180 
181  return compute_jacobian;
182 }
183 
184 // Override the default DiffSystem postprocess function to compute the
185 // approximations to the QoIs
187 {
188  // Reset the array holding the computed QoIs
189  computed_QoI[0] = 0.0;
190 
192 
193  this->comm().sum(computed_QoI[0]);
194 }
virtual void init_context(DiffContext &context)
Definition: poisson.C:95
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
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:350
virtual void zero() libmesh_override
Set every element in the vector to 0.
Definition: dense_vector.h:374
const std::vector< dof_id_type > & get_dof_indices() const
Accessor for element dof indices.
Definition: diff_context.h:366
virtual void postprocess(void)
Runs a postprocessing loop over all elements, and if postprocess_sides is true over all sides...
Definition: poisson.C:186
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
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
Defines a dense subvector for use in finite element computations.
virtual bool element_time_derivative(bool request_jacobian, DiffContext &context)
Adds the time derivative contribution on elem to elem_residual.
Definition: poisson.C:121
const DenseMatrix< Number > & get_elem_jacobian() const
Const accessor for element Jacobian.
Definition: diff_context.h:282
PetscErrorCode Vec x
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_data()
Initializes the member data fields associated with the system, so that, e.g., assemble() may be used...
Definition: poisson.C:55
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
const unsigned int _T_var
Definition: poisson.C:52
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
virtual UniquePtr< FunctionBase< Number > > clone() const
Definition: poisson.C:48
const std::vector< std::vector< OutputGradient > > & get_dphi() const
Definition: fe_base.h:216
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const QBase & get_element_qrule() const
Accessor for element interior quadrature rule for the dimension of the current _elem.
Definition: fem_context.h:765
BdyFunction(unsigned int T_var)
Definition: poisson.C:30
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
This is the base class for functor-like classes.
void get_side_fe(unsigned int var, FEGenericBase< OutputShape > *&fe) const
Accessor for edge/face (2D/3D) finite element object for variable var for the largest dimension in th...
Definition: fem_context.h:299
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.
virtual void postprocess() libmesh_override
Runs a postprocessing loop over all elements, and if postprocess_sides is true over all sides...
Definition: fem_system.C:1107