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
init_data()17 void LaplaceSystem::init_data ()
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
28 FEMSystem::init_data();
29
30 // The temperature is evolving, with a first order time derivative
31 this->time_evolving(T_var, 1);
32 }
33
init_context(DiffContext & context)34 void LaplaceSystem::init_context(DiffContext & context)
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 = 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 = 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
element_time_derivative(bool request_jacobian,DiffContext & context)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 = 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.n_dof_indices(0);
80
81 // The subvectors and submatrices we need to fill:
82 DenseSubMatrix<Number> & K = c.get_elem_jacobian(0, 0);
83 DenseSubVector<Number> & F = c.get_elem_residual(0);
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
side_constraint(bool request_jacobian,DiffContext & context)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 = 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.n_dof_indices(0);
136
137 // The subvectors and submatrices we need to fill:
138 DenseSubMatrix<Number> & K = c.get_elem_jacobian(0, 0);
139 DenseSubVector<Number> & F = c.get_elem_residual(0);
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
postprocess()170 void LaplaceSystem::postprocess()
171 {
172 // Reset the array holding the computed QoIs
173 computed_QoI[0] = 0.0;
174 computed_QoI[1] = 0.0;
175
176 FEMSystem::postprocess();
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
exact_solution(const Point & p)186 Number LaplaceSystem::exact_solution(const Point & p)
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 }
200