libMesh
assembly.C
Go to the documentation of this file.
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 // Bring in bits from the libMesh namespace.
21 // Just the bits we're using, since this is a header.
25 using libMesh::Number;
26 using libMesh::Point;
27 using libMesh::RBTheta;
28 using libMesh::Real;
30 
31 // Kronecker delta function
32 inline Real kronecker_delta(unsigned int i,
33  unsigned int j)
34 {
35  return i == j ? 1. : 0.;
36 }
37 
38 Real elasticity_tensor(unsigned int i,
39  unsigned int j,
40  unsigned int k,
41  unsigned int l)
42 {
43  // Define the Poisson ratio and Young's modulus
44  const Real nu = 0.3;
45  const Real E = 1.;
46 
47  // Define the Lame constants (lambda_1 and lambda_2) based on nu and E
48  const Real lambda_1 = E * nu / ((1. + nu) * (1. - 2.*nu));
49  const Real lambda_2 = 0.5 * E / (1. + nu);
50 
51  return lambda_1 * kronecker_delta(i, j) * kronecker_delta(k, l)
52  + lambda_2 * (kronecker_delta(i, k) * kronecker_delta(j, l) + kronecker_delta(i, l) * kronecker_delta(j, k));
53 }
54 
56 {
57  const unsigned int n_components = rb_sys.n_vars();
58 
59  // make sure we have three components
60  libmesh_assert_equal_to (n_components, 3);
61 
62  const unsigned int u_var = rb_sys.u_var;
63  const unsigned int v_var = rb_sys.v_var;
64  const unsigned int w_var = rb_sys.w_var;
65 
66  FEBase * elem_fe = libmesh_nullptr;
67  c.get_element_fe(u_var, elem_fe);
68 
69  const std::vector<Real> & JxW = elem_fe->get_JxW();
70 
71  // The velocity shape function gradients at interior
72  // quadrature points.
73  const std::vector<std::vector<RealGradient>> & dphi = elem_fe->get_dphi();
74 
75  // Now we will build the affine operator
76  unsigned int n_qpoints = c.get_element_qrule().n_points();
77 
78  std::vector<unsigned int> n_var_dofs(n_components);
79  n_var_dofs[u_var] = c.get_dof_indices(u_var).size();
80  n_var_dofs[v_var] = c.get_dof_indices(v_var).size();
81  n_var_dofs[w_var] = c.get_dof_indices(w_var).size();
82 
83  for (unsigned int C_i = 0; C_i < n_components; C_i++)
84  {
85  unsigned int C_j = 0;
86  for (unsigned int C_k = 0; C_k < n_components; C_k++)
87  for (unsigned int C_l = 1; C_l < n_components; C_l++)
88  {
89  Real C_ijkl = elasticity_tensor(C_i, C_j, C_k, C_l);
90  for (unsigned int qp=0; qp<n_qpoints; qp++)
91  for (unsigned int i=0; i<n_var_dofs[C_i]; i++)
92  for (unsigned int j=0; j<n_var_dofs[C_k]; j++)
93  (c.get_elem_jacobian(C_i,C_k))(i,j) +=
94  JxW[qp]*(C_ijkl * dphi[i][qp](C_j)*dphi[j][qp](C_l));
95  }
96  }
97 
98  for (unsigned int C_i = 0; C_i < n_components; C_i++)
99  for (unsigned int C_j = 1; C_j < n_components; C_j++)
100  for (unsigned int C_k = 0; C_k < n_components; C_k++)
101  {
102  unsigned int C_l = 0;
103 
104  Real C_ijkl = elasticity_tensor(C_i, C_j, C_k, C_l);
105  for (unsigned int qp=0; qp<n_qpoints; qp++)
106  for (unsigned int i=0; i<n_var_dofs[C_i]; i++)
107  for (unsigned int j=0; j<n_var_dofs[C_k]; j++)
108  (c.get_elem_jacobian(C_i,C_k))(i,j) +=
109  JxW[qp]*(C_ijkl * dphi[i][qp](C_j)*dphi[j][qp](C_l));
110  }
111 
112 }
113 
115 {
116  const unsigned int n_components = rb_sys.n_vars();
117 
118  // make sure we have three components
119  libmesh_assert_equal_to (n_components, 3);
120 
121  const unsigned int u_var = rb_sys.u_var;
122  const unsigned int v_var = rb_sys.v_var;
123  const unsigned int w_var = rb_sys.w_var;
124 
125  FEBase * elem_fe = libmesh_nullptr;
126  c.get_element_fe(u_var, elem_fe);
127 
128  const std::vector<Real> & JxW = elem_fe->get_JxW();
129 
130  // The velocity shape function gradients at interior
131  // quadrature points.
132  const std::vector<std::vector<RealGradient>> & dphi = elem_fe->get_dphi();
133 
134  // Now we will build the affine operator
135  unsigned int n_qpoints = c.get_element_qrule().n_points();
136 
137  std::vector<unsigned int> n_var_dofs(n_components);
138  n_var_dofs[u_var] = c.get_dof_indices(u_var).size();
139  n_var_dofs[v_var] = c.get_dof_indices(v_var).size();
140  n_var_dofs[w_var] = c.get_dof_indices(w_var).size();
141 
142  for (unsigned int C_i = 0; C_i < n_components; C_i++)
143  for (unsigned int C_j = 1; C_j < n_components; C_j++)
144  for (unsigned int C_k = 0; C_k < n_components; C_k++)
145  for (unsigned int C_l = 1; C_l < n_components; C_l++)
146  {
147  Real C_ijkl = elasticity_tensor(C_i,C_j,C_k,C_l);
148  for (unsigned int qp=0; qp<n_qpoints; qp++)
149  for (unsigned int i=0; i<n_var_dofs[C_i]; i++)
150  for (unsigned int j=0; j<n_var_dofs[C_k]; j++)
151  (c.get_elem_jacobian(C_i,C_k))(i,j) +=
152  JxW[qp]*(C_ijkl * dphi[i][qp](C_j)*dphi[j][qp](C_l));
153  }
154 }
155 
157 {
158  const unsigned int n_components = rb_sys.n_vars();
159 
160  // make sure we have three components
161  libmesh_assert_equal_to (n_components, 3);
162 
163  const unsigned int u_var = rb_sys.u_var;
164  const unsigned int v_var = rb_sys.v_var;
165  const unsigned int w_var = rb_sys.w_var;
166 
167  FEBase * elem_fe = libmesh_nullptr;
168  c.get_element_fe(u_var, elem_fe);
169 
170  const std::vector<Real> & JxW = elem_fe->get_JxW();
171 
172  // The velocity shape function gradients at interior
173  // quadrature points.
174  const std::vector<std::vector<RealGradient>> & dphi = elem_fe->get_dphi();
175 
176  // Now we will build the affine operator
177  unsigned int n_qpoints = c.get_element_qrule().n_points();
178 
179  std::vector<unsigned int> n_var_dofs(n_components);
180  n_var_dofs[u_var] = c.get_dof_indices(u_var).size();
181  n_var_dofs[v_var] = c.get_dof_indices(v_var).size();
182  n_var_dofs[w_var] = c.get_dof_indices(w_var).size();
183 
184  for (unsigned int C_i = 0; C_i < n_components; C_i++)
185  {
186  unsigned int C_j = 0;
187 
188  for (unsigned int C_k = 0; C_k < n_components; C_k++)
189  {
190  unsigned int C_l = 0;
191 
192  Real C_ijkl = elasticity_tensor(C_i,C_j,C_k,C_l);
193  for (unsigned int qp=0; qp<n_qpoints; qp++)
194  for (unsigned int i=0; i<n_var_dofs[C_i]; i++)
195  for (unsigned int j=0; j<n_var_dofs[C_k]; j++)
196  (c.get_elem_jacobian(C_i,C_k))(i,j) +=
197  JxW[qp]*(C_ijkl * dphi[i][qp](C_j)*dphi[j][qp](C_l));
198  }
199  }
200 }
201 
203 {
205  (&c.get_elem(), c.side, BOUNDARY_ID_MAX_X))
206  {
207  const unsigned int u_var = 0;
208 
209  FEBase * side_fe = libmesh_nullptr;
210  c.get_side_fe(u_var, side_fe);
211 
212  const std::vector<Real> & JxW_side = side_fe->get_JxW();
213 
214  const std::vector<std::vector<Real>> & phi_side = side_fe->get_phi();
215 
216  // The number of local degrees of freedom in each variable
217  const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();
218 
219  // Now we will build the affine operator
220  unsigned int n_qpoints = c.get_side_qrule().n_points();
221  DenseSubVector<Number> & Fu = c.get_elem_residual(u_var);
222 
223  for (unsigned int qp=0; qp < n_qpoints; qp++)
224  for (unsigned int i=0; i < n_u_dofs; i++)
225  Fu(i) += JxW_side[qp] * (1. * phi_side[i][qp]);
226  }
227 }
228 
230 {
232  (&c.get_elem(), c.side, BOUNDARY_ID_MAX_X))
233  {
234  const unsigned int u_var = 0;
235  const unsigned int v_var = 1;
236 
237  FEBase * side_fe = libmesh_nullptr;
238  c.get_side_fe(u_var, side_fe);
239 
240  const std::vector<Real> & JxW_side = side_fe->get_JxW();
241 
242  const std::vector<std::vector<Real>> & phi_side = side_fe->get_phi();
243 
244  // The number of local degrees of freedom in each variable
245  const unsigned int n_v_dofs = c.get_dof_indices(u_var).size();
246 
247  // Now we will build the affine operator
248  unsigned int n_qpoints = c.get_side_qrule().n_points();
249  DenseSubVector<Number> & Fv = c.get_elem_residual(v_var);
250 
251  for (unsigned int qp=0; qp < n_qpoints; qp++)
252  for (unsigned int i=0; i < n_v_dofs; i++)
253  Fv(i) += JxW_side[qp] * (1. * phi_side[i][qp]);
254  }
255 }
256 
258 {
260  (&c.get_elem(), c.side, BOUNDARY_ID_MAX_X))
261  {
262  const unsigned int u_var = 0;
263  const unsigned int w_var = 2;
264 
265  FEBase * side_fe = libmesh_nullptr;
266  c.get_side_fe(u_var, side_fe);
267 
268  const std::vector<Real> & JxW_side = side_fe->get_JxW();
269 
270  const std::vector<std::vector<Real>> & phi_side = side_fe->get_phi();
271 
272  // The number of local degrees of freedom in each variable
273  const unsigned int n_w_dofs = c.get_dof_indices(w_var).size();
274 
275  // Now we will build the affine operator
276  unsigned int n_qpoints = c.get_side_qrule().n_points();
277  DenseSubVector<Number> & Fw = c.get_elem_residual(w_var);
278 
279  for (unsigned int qp=0; qp < n_qpoints; qp++)
280  for (unsigned int i=0; i < n_w_dofs; i++)
281  Fw(i) += JxW_side[qp] * (1. * phi_side[i][qp]);
282  }
283 }
284 
285 void AssemblyPointLoadX::get_nodal_rhs_values(std::map<numeric_index_type, Number> & values,
286  const System & sys,
287  const Node & node)
288 {
289  // First clear the values map
290  values.clear();
291 
292  if (sys.get_mesh().get_boundary_info().has_boundary_id
293  (&node, NODE_BOUNDARY_ID))
294  {
295  numeric_index_type dof_index =
296  node.dof_number(sys.number(), sys.variable_number("u"), 0);
297  values[dof_index] = 1.;
298  }
299 }
300 
301 void AssemblyPointLoadY::get_nodal_rhs_values(std::map<numeric_index_type, Number> & values,
302  const System & sys,
303  const Node & node)
304 {
305  // First clear the values map
306  values.clear();
307 
308  if (sys.get_mesh().get_boundary_info().has_boundary_id
309  (&node, NODE_BOUNDARY_ID))
310  {
311  numeric_index_type dof_index =
312  node.dof_number(sys.number(), sys.variable_number("v"), 0);
313  values[dof_index] = 1.;
314  }
315 }
316 
317 void AssemblyPointLoadZ::get_nodal_rhs_values(std::map<numeric_index_type, Number> & values,
318  const System & sys,
319  const Node & node)
320 {
321  // First clear the values map
322  values.clear();
323 
324  if (sys.get_mesh().get_boundary_info().has_boundary_id
325  (&node, NODE_BOUNDARY_ID))
326  {
327  numeric_index_type dof_index =
328  node.dof_number(sys.number(), sys.variable_number("w"), 0);
329  values[dof_index] = 1.;
330  }
331 }
332 
334 {
335  const unsigned int u_var = rb_sys.u_var;
336  const unsigned int v_var = rb_sys.v_var;
337  const unsigned int w_var = rb_sys.w_var;
338 
339  FEBase * elem_fe = libmesh_nullptr;
340  c.get_element_fe(u_var, elem_fe);
341 
342  const std::vector<Real> & JxW = elem_fe->get_JxW();
343 
344  // The velocity shape function gradients at interior
345  // quadrature points.
346  const std::vector<std::vector<RealGradient>>& dphi = elem_fe->get_dphi();
347 
348  // The number of local degrees of freedom in each variable
349  const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();
350  const unsigned int n_v_dofs = c.get_dof_indices(v_var).size();
351  const unsigned int n_w_dofs = c.get_dof_indices(w_var).size();
352 
353  // Now we will build the affine operator
354  unsigned int n_qpoints = c.get_element_qrule().n_points();
355 
356  DenseSubMatrix<Number> & Kuu = c.get_elem_jacobian(u_var, u_var);
357  DenseSubMatrix<Number> & Kvv = c.get_elem_jacobian(v_var, v_var);
358  DenseSubMatrix<Number> & Kww = c.get_elem_jacobian(w_var, w_var);
359 
360  for (unsigned int qp=0; qp<n_qpoints; qp++)
361  {
362  for (unsigned int i=0; i<n_u_dofs; i++)
363  for (unsigned int j=0; j<n_u_dofs; j++)
364  Kuu(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
365 
366  for (unsigned int i=0; i<n_v_dofs; i++)
367  for (unsigned int j=0; j<n_v_dofs; j++)
368  Kvv(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
369 
370  for (unsigned int i=0; i<n_w_dofs; i++)
371  for (unsigned int j=0; j<n_w_dofs; j++)
372  Kww(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
373  }
374 }
RealVectorValue RealGradient
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:117
virtual void get_nodal_rhs_values(std::map< numeric_index_type, Number > &values, const System &sys, const Node &node)
Definition: assembly.C:317
ImplicitSystem & sys
virtual void get_nodal_rhs_values(std::map< numeric_index_type, Number > &values, const System &sys, const Node &node)
Definition: assembly.C:301
const QBase & get_side_qrule() const
Accessor for element side quadrature rule for the dimension of the current _elem. ...
Definition: fem_context.h:772
virtual void boundary_assembly(FEMContext &c)
Perform the element boundary assembly.
Definition: assembly.C:229
unsigned int u_var
Variable numbers.
Definition: rb_classes.h:121
ElasticityRBConstruction & rb_sys
The ElasticityRBConstruction object that will use this assembly.
Definition: assembly.h:56
const std::vector< dof_id_type > & get_dof_indices() const
Accessor for element dof indices.
Definition: diff_context.h:366
Real elasticity_tensor(unsigned int i, unsigned int j, unsigned int k, unsigned int l)
Definition: assembly.C:38
const class libmesh_nullptr_t libmesh_nullptr
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
virtual void boundary_assembly(FEMContext &c)
Perform the element boundary assembly.
Definition: assembly.C:257
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
virtual void boundary_assembly(FEMContext &c)
Perform the element boundary assembly.
Definition: assembly.C:202
const DenseMatrix< Number > & get_elem_jacobian() const
Const accessor for element Jacobian.
Definition: diff_context.h:282
virtual void interior_assembly(FEMContext &c)
Perform the element interior assembly.
Definition: assembly.C:55
const MeshBase & get_mesh() const
Definition: system.h:2014
dof_id_type numeric_index_type
Definition: id_types.h:92
virtual void interior_assembly(FEMContext &c)
Perform the element interior assembly.
Definition: assembly.C:156
This class provides all data required for a physics package (e.g.
Definition: fem_context.h:61
Real kronecker_delta(unsigned int i, unsigned int j)
Definition: assembly.C:32
FEGenericBase< Real > FEBase
virtual void get_nodal_rhs_values(std::map< numeric_index_type, Number > &values, const System &sys, const Node &node)
Definition: assembly.C:285
This class provides an encapsulated access to all static public member functions of finite element cl...
Definition: fe_interface.h:63
virtual void interior_assembly(FEMContext &c)
Perform the element interior assembly.
Definition: assembly.C:333
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
unsigned char side
Current side for side_* to examine.
Definition: fem_context.h:977
ElemAssembly provides a per-element (interior and boundary) assembly functionality.
Definition: elem_assembly.h:38
bool has_boundary_id(const Node *const node, const boundary_id_type id) const
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 void interior_assembly(FEMContext &c)
Perform the element interior assembly.
Definition: assembly.C:114
unsigned int n_vars() const
Definition: system.h:2086
This class is part of the rbOOmit framework.
Definition: rb_theta.h:46
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