1 // local includes
2 #include "assembly.h"
3 #include "rb_classes.h"
4
5 // libMesh includes
6 #include "libmesh/sparse_matrix.h"
7 #include "libmesh/numeric_vector.h"
8 #include "libmesh/dense_matrix.h"
9 #include "libmesh/dense_submatrix.h"
10 #include "libmesh/dense_vector.h"
11 #include "libmesh/dense_subvector.h"
12 #include "libmesh/fe.h"
13 #include "libmesh/fe_interface.h"
14 #include "libmesh/fe_base.h"
15 #include "libmesh/elem_assembly.h"
16 #include "libmesh/quadrature_gauss.h"
17 #include "libmesh/boundary_info.h"
18 #include "libmesh/node.h"
19
20 // C++ includes
21 #include <functional>
22
23 #ifdef LIBMESH_ENABLE_DIRICHLET
24
25 // Bring in bits from the libMesh namespace.
26 // Just the bits we're using, since this is a header.
27 using libMesh::ElemAssembly;
28 using libMesh::FEInterface;
29 using libMesh::FEMContext;
30 using libMesh::Number;
31 using libMesh::Point;
32 using libMesh::RBTheta;
33 using libMesh::Real;
34 using libMesh::RealGradient;
35
36 // Kronecker delta function
kronecker_delta(unsigned int i,unsigned int j)37 inline Real kronecker_delta(unsigned int i,
38 unsigned int j)
39 {
40 return i == j ? 1. : 0.;
41 }
42
elasticity_tensor(unsigned int i,unsigned int j,unsigned int k,unsigned int l)43 Real elasticity_tensor(unsigned int i,
44 unsigned int j,
45 unsigned int k,
46 unsigned int l)
47 {
48 // Define the Poisson ratio and Young's modulus
49 const Real nu = 0.3;
50 const Real E = 1.;
51
52 // Define the Lame constants (lambda_1 and lambda_2) based on nu and E
53 const Real lambda_1 = E * nu / ((1. + nu) * (1. - 2.*nu));
54 const Real lambda_2 = 0.5 * E / (1. + nu);
55
56 return lambda_1 * kronecker_delta(i, j) * kronecker_delta(k, l)
57 + lambda_2 * (kronecker_delta(i, k) * kronecker_delta(j, l) + kronecker_delta(i, l) * kronecker_delta(j, k));
58 }
59
interior_assembly(FEMContext & c)60 void AssemblyA0::interior_assembly(FEMContext & c)
61 {
62 const unsigned int n_components = rb_sys.n_vars();
63
64 // make sure we have three components
65 libmesh_assert_equal_to (n_components, 3);
66
67 const unsigned int u_var = rb_sys.u_var;
68 const unsigned int v_var = rb_sys.v_var;
69 const unsigned int w_var = rb_sys.w_var;
70
71 FEBase * elem_fe = nullptr;
72 c.get_element_fe(u_var, elem_fe);
73
74 const std::vector<Real> & JxW = elem_fe->get_JxW();
75
76 // The velocity shape function gradients at interior
77 // quadrature points.
78 const std::vector<std::vector<RealGradient>> & dphi = elem_fe->get_dphi();
79
80 // Now we will build the affine operator
81 unsigned int n_qpoints = c.get_element_qrule().n_points();
82
83 std::vector<unsigned int> n_var_dofs(n_components);
84 n_var_dofs[u_var] = c.get_dof_indices(u_var).size();
85 n_var_dofs[v_var] = c.get_dof_indices(v_var).size();
86 n_var_dofs[w_var] = c.get_dof_indices(w_var).size();
87
88 for (unsigned int C_i = 0; C_i < n_components; C_i++)
89 {
90 unsigned int C_j = 0;
91 for (unsigned int C_k = 0; C_k < n_components; C_k++)
92 for (unsigned int C_l = 1; C_l < n_components; C_l++)
93 {
94 Real C_ijkl = elasticity_tensor(C_i, C_j, C_k, C_l);
95 for (unsigned int qp=0; qp<n_qpoints; qp++)
96 for (unsigned int i=0; i<n_var_dofs[C_i]; i++)
97 for (unsigned int j=0; j<n_var_dofs[C_k]; j++)
98 (c.get_elem_jacobian(C_i,C_k))(i,j) +=
99 JxW[qp]*(C_ijkl * dphi[i][qp](C_j)*dphi[j][qp](C_l));
100 }
101 }
102
103 for (unsigned int C_i = 0; C_i < n_components; C_i++)
104 for (unsigned int C_j = 1; C_j < n_components; C_j++)
105 for (unsigned int C_k = 0; C_k < n_components; C_k++)
106 {
107 unsigned int C_l = 0;
108
109 Real C_ijkl = elasticity_tensor(C_i, C_j, C_k, C_l);
110 for (unsigned int qp=0; qp<n_qpoints; qp++)
111 for (unsigned int i=0; i<n_var_dofs[C_i]; i++)
112 for (unsigned int j=0; j<n_var_dofs[C_k]; j++)
113 (c.get_elem_jacobian(C_i,C_k))(i,j) +=
114 JxW[qp]*(C_ijkl * dphi[i][qp](C_j)*dphi[j][qp](C_l));
115 }
116
117 }
118
interior_assembly(FEMContext & c)119 void AssemblyA1::interior_assembly(FEMContext & c)
120 {
121 const unsigned int n_components = rb_sys.n_vars();
122
123 // make sure we have three components
124 libmesh_assert_equal_to (n_components, 3);
125
126 const unsigned int u_var = rb_sys.u_var;
127 const unsigned int v_var = rb_sys.v_var;
128 const unsigned int w_var = rb_sys.w_var;
129
130 FEBase * elem_fe = nullptr;
131 c.get_element_fe(u_var, elem_fe);
132
133 const std::vector<Real> & JxW = elem_fe->get_JxW();
134
135 // The velocity shape function gradients at interior
136 // quadrature points.
137 const std::vector<std::vector<RealGradient>> & dphi = elem_fe->get_dphi();
138
139 // Now we will build the affine operator
140 unsigned int n_qpoints = c.get_element_qrule().n_points();
141
142 std::vector<unsigned int> n_var_dofs(n_components);
143 n_var_dofs[u_var] = c.get_dof_indices(u_var).size();
144 n_var_dofs[v_var] = c.get_dof_indices(v_var).size();
145 n_var_dofs[w_var] = c.get_dof_indices(w_var).size();
146
147 for (unsigned int C_i = 0; C_i < n_components; C_i++)
148 for (unsigned int C_j = 1; C_j < n_components; C_j++)
149 for (unsigned int C_k = 0; C_k < n_components; C_k++)
150 for (unsigned int C_l = 1; C_l < n_components; C_l++)
151 {
152 Real C_ijkl = elasticity_tensor(C_i,C_j,C_k,C_l);
153 for (unsigned int qp=0; qp<n_qpoints; qp++)
154 for (unsigned int i=0; i<n_var_dofs[C_i]; i++)
155 for (unsigned int j=0; j<n_var_dofs[C_k]; j++)
156 (c.get_elem_jacobian(C_i,C_k))(i,j) +=
157 JxW[qp]*(C_ijkl * dphi[i][qp](C_j)*dphi[j][qp](C_l));
158 }
159 }
160
interior_assembly(FEMContext & c)161 void AssemblyA2::interior_assembly(FEMContext & c)
162 {
163 const unsigned int n_components = rb_sys.n_vars();
164
165 // make sure we have three components
166 libmesh_assert_equal_to (n_components, 3);
167
168 const unsigned int u_var = rb_sys.u_var;
169 const unsigned int v_var = rb_sys.v_var;
170 const unsigned int w_var = rb_sys.w_var;
171
172 FEBase * elem_fe = nullptr;
173 c.get_element_fe(u_var, elem_fe);
174
175 const std::vector<Real> & JxW = elem_fe->get_JxW();
176
177 // The velocity shape function gradients at interior
178 // quadrature points.
179 const std::vector<std::vector<RealGradient>> & dphi = elem_fe->get_dphi();
180
181 // Now we will build the affine operator
182 unsigned int n_qpoints = c.get_element_qrule().n_points();
183
184 std::vector<unsigned int> n_var_dofs(n_components);
185 n_var_dofs[u_var] = c.get_dof_indices(u_var).size();
186 n_var_dofs[v_var] = c.get_dof_indices(v_var).size();
187 n_var_dofs[w_var] = c.get_dof_indices(w_var).size();
188
189 for (unsigned int C_i = 0; C_i < n_components; C_i++)
190 {
191 unsigned int C_j = 0;
192
193 for (unsigned int C_k = 0; C_k < n_components; C_k++)
194 {
195 unsigned int C_l = 0;
196
197 Real C_ijkl = elasticity_tensor(C_i,C_j,C_k,C_l);
198 for (unsigned int qp=0; qp<n_qpoints; qp++)
199 for (unsigned int i=0; i<n_var_dofs[C_i]; i++)
200 for (unsigned int j=0; j<n_var_dofs[C_k]; j++)
201 (c.get_elem_jacobian(C_i,C_k))(i,j) +=
202 JxW[qp]*(C_ijkl * dphi[i][qp](C_j)*dphi[j][qp](C_l));
203 }
204 }
205 }
206
boundary_assembly(FEMContext & c)207 void AssemblyF0::boundary_assembly(FEMContext & c)
208 {
209 if (rb_sys.get_mesh().get_boundary_info().has_boundary_id
210 (&c.get_elem(), c.side, BOUNDARY_ID_MAX_X))
211 {
212 const unsigned int u_var = 0;
213
214 FEBase * side_fe = nullptr;
215 c.get_side_fe(u_var, side_fe);
216
217 const std::vector<Real> & JxW_side = side_fe->get_JxW();
218
219 const std::vector<std::vector<Real>> & phi_side = side_fe->get_phi();
220
221 // The number of local degrees of freedom in each variable
222 const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();
223
224 // Now we will build the affine operator
225 unsigned int n_qpoints = c.get_side_qrule().n_points();
226 DenseSubVector<Number> & Fu = c.get_elem_residual(u_var);
227
228 for (unsigned int qp=0; qp < n_qpoints; qp++)
229 for (unsigned int i=0; i < n_u_dofs; i++)
230 Fu(i) += JxW_side[qp] * (1. * phi_side[i][qp]);
231 }
232 }
233
boundary_assembly(FEMContext & c)234 void AssemblyF1::boundary_assembly(FEMContext & c)
235 {
236 if (rb_sys.get_mesh().get_boundary_info().has_boundary_id
237 (&c.get_elem(), c.side, BOUNDARY_ID_MAX_X))
238 {
239 const unsigned int u_var = 0;
240 const unsigned int v_var = 1;
241
242 FEBase * side_fe = nullptr;
243 c.get_side_fe(u_var, side_fe);
244
245 const std::vector<Real> & JxW_side = side_fe->get_JxW();
246
247 const std::vector<std::vector<Real>> & phi_side = side_fe->get_phi();
248
249 // The number of local degrees of freedom in each variable
250 const unsigned int n_v_dofs = c.get_dof_indices(u_var).size();
251
252 // Now we will build the affine operator
253 unsigned int n_qpoints = c.get_side_qrule().n_points();
254 DenseSubVector<Number> & Fv = c.get_elem_residual(v_var);
255
256 for (unsigned int qp=0; qp < n_qpoints; qp++)
257 for (unsigned int i=0; i < n_v_dofs; i++)
258 Fv(i) += JxW_side[qp] * (1. * phi_side[i][qp]);
259 }
260 }
261
boundary_assembly(FEMContext & c)262 void AssemblyF2::boundary_assembly(FEMContext & c)
263 {
264 if (rb_sys.get_mesh().get_boundary_info().has_boundary_id
265 (&c.get_elem(), c.side, BOUNDARY_ID_MAX_X))
266 {
267 const unsigned int u_var = 0;
268 const unsigned int w_var = 2;
269
270 FEBase * side_fe = nullptr;
271 c.get_side_fe(u_var, side_fe);
272
273 const std::vector<Real> & JxW_side = side_fe->get_JxW();
274
275 const std::vector<std::vector<Real>> & phi_side = side_fe->get_phi();
276
277 // The number of local degrees of freedom in each variable
278 const unsigned int n_w_dofs = c.get_dof_indices(w_var).size();
279
280 // Now we will build the affine operator
281 unsigned int n_qpoints = c.get_side_qrule().n_points();
282 DenseSubVector<Number> & Fw = c.get_elem_residual(w_var);
283
284 for (unsigned int qp=0; qp < n_qpoints; qp++)
285 for (unsigned int i=0; i < n_w_dofs; i++)
286 Fw(i) += JxW_side[qp] * (1. * phi_side[i][qp]);
287 }
288 }
289
get_nodal_rhs_values(std::map<numeric_index_type,Number> & values,const System & sys,const Node & node)290 void AssemblyPointLoadX::get_nodal_rhs_values(std::map<numeric_index_type, Number> & values,
291 const System & sys,
292 const Node & node)
293 {
294 // First clear the values map
295 values.clear();
296
297 if (sys.get_mesh().get_boundary_info().has_boundary_id
298 (&node, NODE_BOUNDARY_ID))
299 {
300 numeric_index_type dof_index =
301 node.dof_number(sys.number(), sys.variable_number("u"), 0);
302 values[dof_index] = 1.;
303 }
304 }
305
get_nodal_rhs_values(std::map<numeric_index_type,Number> & values,const System & sys,const Node & node)306 void AssemblyPointLoadY::get_nodal_rhs_values(std::map<numeric_index_type, Number> & values,
307 const System & sys,
308 const Node & node)
309 {
310 // First clear the values map
311 values.clear();
312
313 if (sys.get_mesh().get_boundary_info().has_boundary_id
314 (&node, NODE_BOUNDARY_ID))
315 {
316 numeric_index_type dof_index =
317 node.dof_number(sys.number(), sys.variable_number("v"), 0);
318 values[dof_index] = 1.;
319 }
320 }
321
get_nodal_rhs_values(std::map<numeric_index_type,Number> & values,const System & sys,const Node & node)322 void AssemblyPointLoadZ::get_nodal_rhs_values(std::map<numeric_index_type, Number> & values,
323 const System & sys,
324 const Node & node)
325 {
326 // First clear the values map
327 values.clear();
328
329 if (sys.get_mesh().get_boundary_info().has_boundary_id
330 (&node, NODE_BOUNDARY_ID))
331 {
332 numeric_index_type dof_index =
333 node.dof_number(sys.number(), sys.variable_number("w"), 0);
334 values[dof_index] = 1.;
335 }
336 }
337
interior_assembly(FEMContext & c)338 void InnerProductAssembly::interior_assembly(FEMContext & c)
339 {
340 const unsigned int u_var = rb_sys.u_var;
341 const unsigned int v_var = rb_sys.v_var;
342 const unsigned int w_var = rb_sys.w_var;
343
344 FEBase * elem_fe = nullptr;
345 c.get_element_fe(u_var, elem_fe);
346
347 const std::vector<Real> & JxW = elem_fe->get_JxW();
348
349 // The velocity shape function gradients at interior
350 // quadrature points.
351 const std::vector<std::vector<RealGradient>>& dphi = elem_fe->get_dphi();
352
353 // (u_var, v_var, w_var) all have the same number of dofs
354 const unsigned int n_dofs = c.get_dof_indices(u_var).size();
355
356 // Now we will build the affine operator
357 unsigned int n_qpoints = c.get_element_qrule().n_points();
358
359 // Get references to stiffness matrix diagonal blocks
360 std::reference_wrapper<DenseSubMatrix<Number>> Kdiag[3] =
361 {
362 c.get_elem_jacobian(u_var, u_var),
363 c.get_elem_jacobian(v_var, v_var),
364 c.get_elem_jacobian(w_var, w_var)
365 };
366
367 for (unsigned int qp=0; qp<n_qpoints; qp++)
368 for (unsigned int d=0; d<3; ++d)
369 for (unsigned int i=0; i<n_dofs; i++)
370 for (unsigned int j=0; j<n_dofs; j++)
371 Kdiag[d](i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
372 }
373
374 #endif // LIBMESH_ENABLE_DIRICHLET
375