libMesh
L-shaped.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/fe_base.h"
6 #include "libmesh/quadrature.h"
7 #include "libmesh/string_to_enum.h"
8 #include "libmesh/parallel.h"
9 #include "libmesh/fem_context.h"
10 
11 // Local includes
12 #include "L-shaped.h"
13 
14 // Bring in everything from the libMesh namespace
15 using namespace libMesh;
16 
18 {
19  unsigned int T_var =
20  this->add_variable ("T", static_cast<Order>(_fe_order),
21  Utility::string_to_enum<FEFamily>(_fe_family));
22 
23  GetPot infile("l-shaped.in");
24  exact_QoI[0] = infile("QoI_0", 0.0);
25  exact_QoI[1] = infile("QoI_1", 0.0);
26 
27  // Do the parent's initialization after variables are defined
29 
30  // The temperature is evolving, with a first order time derivative
31  this->time_evolving(T_var, 1);
32 }
33 
35 {
36  FEMContext & c = cast_ref<FEMContext &>(context);
37 
38  // Now make sure we have requested all the data
39  // we need to build the linear system.
40  FEBase * elem_fe = libmesh_nullptr;
41  c.get_element_fe(0, elem_fe);
42  elem_fe->get_JxW();
43  elem_fe->get_phi();
44  elem_fe->get_dphi();
45  elem_fe->get_xyz();
46 
47  FEBase * side_fe = libmesh_nullptr;
48  c.get_side_fe(0, side_fe);
49 
50  side_fe->get_JxW();
51  side_fe->get_phi();
52  side_fe->get_dphi();
53  side_fe->get_xyz();
54 }
55 
56 #define optassert(X) {if (!(X)) libmesh_error_msg("Assertion " #X " failed.");}
57 
58 // Assemble the element contributions to the stiffness matrix
59 bool LaplaceSystem::element_time_derivative (bool request_jacobian,
60  DiffContext & context)
61 {
62  // Are the jacobians specified analytically ?
63  bool compute_jacobian = request_jacobian && _analytic_jacobians;
64 
65  FEMContext & c = cast_ref<FEMContext &>(context);
66 
67  // First we get some references to cell-specific data that
68  // will be used to assemble the linear system.
69  FEBase * elem_fe = libmesh_nullptr;
70  c.get_element_fe(0, elem_fe);
71 
72  // Element Jacobian * quadrature weights for interior integration
73  const std::vector<Real> & JxW = elem_fe->get_JxW();
74 
75  // Element basis functions
76  const std::vector<std::vector<RealGradient>> & dphi = elem_fe->get_dphi();
77 
78  // The number of local degrees of freedom in each variable
79  const unsigned int n_T_dofs = c.get_dof_indices(0).size();
80 
81  // The subvectors and submatrices we need to fill:
84 
85  // Now we will build the element Jacobian and residual.
86  // Constructing the residual requires the solution and its
87  // gradient from the previous timestep. This must be
88  // calculated at each quadrature point by summing the
89  // solution degree-of-freedom values by the appropriate
90  // weight functions.
91  unsigned int n_qpoints = c.get_element_qrule().n_points();
92 
93  for (unsigned int qp=0; qp != n_qpoints; qp++)
94  {
95  // Compute the solution gradient at the Newton iterate
96  Gradient grad_T = c.interior_gradient(0, qp);
97 
98  // The residual contribution from this element
99  for (unsigned int i=0; i != n_T_dofs; i++)
100  F(i) += JxW[qp] * (grad_T * dphi[i][qp]);
101  if (compute_jacobian)
102  for (unsigned int i=0; i != n_T_dofs; i++)
103  for (unsigned int j=0; j != n_T_dofs; ++j)
104  // The analytic jacobian
105  K(i,j) += JxW[qp] * (dphi[i][qp] * dphi[j][qp]);
106  } // end of the quadrature point qp-loop
107 
108  return compute_jacobian;
109 }
110 
111 // Set Dirichlet bcs, side contributions to global stiffness matrix
112 bool LaplaceSystem::side_constraint (bool request_jacobian,
113  DiffContext & context)
114 {
115  // Are the jacobians specified analytically ?
116  bool compute_jacobian = request_jacobian && _analytic_jacobians;
117 
118  FEMContext & c = cast_ref<FEMContext &>(context);
119 
120  // First we get some references to cell-specific data that
121  // will be used to assemble the linear system.
122  FEBase * side_fe = libmesh_nullptr;
123  c.get_side_fe(0, side_fe);
124 
125  // Element Jacobian * quadrature weights for interior integration
126  const std::vector<Real> & JxW = side_fe->get_JxW();
127 
128  // Side basis functions
129  const std::vector<std::vector<Real>> & phi = side_fe->get_phi();
130 
131  // Side Quadrature points
132  const std::vector<Point > & qside_point = side_fe->get_xyz();
133 
134  // The number of local degrees of freedom in each variable
135  const unsigned int n_T_dofs = c.get_dof_indices(0).size();
136 
137  // The subvectors and submatrices we need to fill:
140 
141  unsigned int n_qpoints = c.get_side_qrule().n_points();
142 
143  const Real penalty = 1./(TOLERANCE*TOLERANCE);
144 
145  for (unsigned int qp=0; qp != n_qpoints; qp++)
146  {
147  // Compute the solution at the old Newton iterate
148  Number T = c.side_value(0, qp);
149 
150  // We get the Dirichlet bcs from the exact solution
151  Number u_dirichlet = exact_solution (qside_point[qp]);
152 
153  // The residual from the boundary terms, penalize non-zero temperature
154  for (unsigned int i=0; i != n_T_dofs; i++)
155  F(i) += JxW[qp] * penalty * (T - u_dirichlet) * phi[i][qp];
156 
157  if (compute_jacobian)
158  for (unsigned int i=0; i != n_T_dofs; i++)
159  for (unsigned int j=0; j != n_T_dofs; ++j)
160  // The analytic jacobian
161  K(i,j) += JxW[qp] * penalty * phi[i][qp] * phi[j][qp];
162 
163  } // end of the quadrature point qp-loop
164 
165  return compute_jacobian;
166 }
167 
168 // Override the default DiffSystem postprocess function to compute the
169 // approximations to the QoIs
171 {
172  // Reset the array holding the computed QoIs
173  computed_QoI[0] = 0.0;
174  computed_QoI[1] = 0.0;
175 
177 
178  this->comm().sum(computed_QoI[0]);
179 
180  this->comm().sum(computed_QoI[1]);
181 
182 }
183 
184 // The exact solution to the singular problem,
185 // u_exact = r^(2/3)*sin(2*theta/3). We use this to set the Dirichlet boundary conditions
187 {
188  const Real x1 = p(0);
189  const Real x2 = p(1);
190 
191  Real theta = atan2(x2, x1);
192 
193  // Make sure 0 <= theta <= 2*pi
194  if (theta < 0)
195  theta += 2. * libMesh::pi;
196 
197  return pow(x1*x1 + x2*x2, 1./3.)*sin(2./3.*theta);
198 
199 }
virtual void init_context(DiffContext &context)
Definition: L-shaped.C:34
This class provides all data required for a physics package (e.g.
Definition: diff_context.h:54
Number(* exact_solution)(const Point &p, const Parameters &, const std::string &, const std::string &)
const QBase & get_side_qrule() const
Accessor for element side quadrature rule for the dimension of the current _elem. ...
Definition: fem_context.h:772
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
virtual bool element_time_derivative(bool request_jacobian, DiffContext &context)
Adds the time derivative contribution on elem to elem_residual.
Definition: L-shaped.C:58
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
static const Real TOLERANCE
The libMesh namespace provides an interface to certain functionality in the library.
const DenseVector< Number > & get_elem_residual() const
Const accessor for element residual.
Definition: diff_context.h:248
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::vector< Real > & get_JxW() const
Definition: fe_abstract.h:237
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
double pow(double a, int b)
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
Number side_value(unsigned int var, unsigned int qp) const
Definition: fem_context.C:575
const QBase & get_element_qrule() const
Accessor for element interior quadrature rule for the dimension of the current _elem.
Definition: fem_context.h:765
LaplaceExactSolution exact_solution
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
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
virtual bool side_constraint(bool request_jacobian, DiffContext &context)
Adds the constraint contribution on side of elem to elem_residual.
Definition: L-shaped.C:111
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
virtual void init_data()
Initializes the member data fields associated with the system, so that, e.g., assemble() may be used...
Definition: L-shaped.C:17
unsigned int n_points() const
Definition: quadrature.h:113
const std::vector< Point > & get_xyz() const
Definition: fe_abstract.h:230
virtual void postprocess()
Runs a postprocessing loop over all elements, and if postprocess_sides is true over all sides...
Definition: L-shaped.C:168
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
const Real pi
.
Definition: libmesh.h:172