libMesh
assembly.h
Go to the documentation of this file.
1 #ifndef ASSEMBLY_H
2 #define ASSEMBLY_H
3 
4 #include "libmesh/sparse_matrix.h"
5 #include "libmesh/numeric_vector.h"
6 #include "libmesh/dense_matrix.h"
7 #include "libmesh/dense_vector.h"
8 #include "libmesh/fe.h"
9 #include "libmesh/fe_interface.h"
10 #include "libmesh/fe_base.h"
11 #include "libmesh/elem_assembly.h"
12 #include "libmesh/quadrature_gauss.h"
13 
14 // rbOOmit includes
15 #include "libmesh/rb_theta.h"
16 #include "libmesh/rb_assembly_expansion.h"
17 
18 #define damping_epsilon 0.001
19 #define R_rad 12.0
20 
21 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
22 
23 // Bring in bits from the libMesh namespace.
24 // Just the bits we're using, since this is a header.
28 using libMesh::Number;
29 using libMesh::Point;
32 using libMesh::RBTheta;
34 using libMesh::Real;
36 using libMesh::MeshBase;
37 using libMesh::FEBase;
38 
39 // Functors for the parameter-dependent part of the affine decomposition of the PDE
40 struct ThetaA0 : RBTheta
41 {
42  virtual Number evaluate(const RBParameters & mu)
43  {
44  return Number(1., mu.get_value("frequency")*damping_epsilon);
45  }
46 };
47 
48 struct ThetaA1 : RBTheta
49 {
50  virtual Number evaluate(const RBParameters & mu)
51  {
52  return Number(-mu.get_value("frequency")*mu.get_value("frequency"), 0.);
53  }
54 };
55 
56 struct ThetaA2 : RBTheta
57 {
58  virtual Number evaluate(const RBParameters & mu)
59  {
60  return Number(0., mu.get_value("frequency"));
61  }
62 };
63 
64 struct ThetaA3 : RBTheta
65 {
66  virtual Number evaluate(const RBParameters & mu)
67  {
68  return Number(0.5/R_rad, mu.get_value("frequency"));
69  }
70 };
71 
72 struct ThetaF0 : RBTheta
73 {
74  virtual Number evaluate(const RBParameters & mu)
75  {
76  return Number(0., 2.*mu.get_value("frequency"));
77  }
78 };
79 
81 {
82  virtual Number evaluate(const RBParameters &)
83  {
84  return Number(1., 0.);
85  }
86 };
87 
89 {
90  virtual void interior_assembly(FEMContext & c)
91  {
92  const unsigned int p_var = 0;
93 
94  FEBase * elem_fe = libmesh_nullptr;
95  c.get_element_fe(p_var, elem_fe);
96 
97  const std::vector<Real> & JxW = elem_fe->get_JxW();
98 
99  const std::vector<std::vector<Real>> & phi = elem_fe->get_phi();
100 
101  const std::vector<std::vector<RealGradient>> & dphi = elem_fe->get_dphi();
102 
103  // The number of local degrees of freedom in each variable
104  const unsigned int n_p_dofs = c.get_dof_indices(p_var).size();
105 
106  // Now we will build the affine operator
107  unsigned int n_qpoints = c.get_element_qrule().n_points();
108 
109  // We don't need to conjugate phi or dphi, since basis functions are real-valued
110  for (unsigned int qp=0; qp != n_qpoints; qp++)
111  for (unsigned int i=0; i != n_p_dofs; i++)
112  for (unsigned int j=0; j != n_p_dofs; j++)
113  c.get_elem_jacobian()(i,j) += JxW[qp] * ((dphi[j][qp]*dphi[i][qp]) +
114  (phi[j][qp]*phi[i][qp]));
115  }
116 };
117 
118 struct A0 : ElemAssembly
119 {
120  virtual void interior_assembly(FEMContext & c)
121  {
122  const unsigned int p_var = 0;
123 
124  FEBase * elem_fe = libmesh_nullptr;
125  c.get_element_fe(p_var, elem_fe);
126 
127  const std::vector<Real> & JxW = elem_fe->get_JxW();
128 
129  // The velocity shape function gradients at interior
130  // quadrature points.
131  const std::vector<std::vector<RealGradient>> & dphi = elem_fe->get_dphi();
132 
133  // The number of local degrees of freedom in each variable
134  const unsigned int n_p_dofs = c.get_dof_indices(p_var).size();
135 
136  // Now we will build the affine operator
137  unsigned int n_qpoints = c.get_element_qrule().n_points();
138 
139  for (unsigned int qp=0; qp != n_qpoints; qp++)
140  for (unsigned int i=0; i != n_p_dofs; i++)
141  for (unsigned int j=0; j != n_p_dofs; j++)
142  c.get_elem_jacobian()(i,j) += JxW[qp] * (dphi[j][qp]*dphi[i][qp]);
143  }
144 };
145 
146 
147 struct A1 : ElemAssembly
148 {
149  virtual void interior_assembly(FEMContext & c)
150  {
151  const unsigned int p_var = 0;
152 
153  FEBase * elem_fe = libmesh_nullptr;
154  c.get_element_fe(p_var, elem_fe);
155 
156  const std::vector<Real> & JxW = elem_fe->get_JxW();
157 
158  const std::vector<std::vector<Real>> & phi = elem_fe->get_phi();
159 
160  // The number of local degrees of freedom in each variable
161  const unsigned int n_p_dofs = c.get_dof_indices(p_var).size();
162 
163  // Now we will build the affine operator
164  unsigned int n_qpoints = c.get_element_qrule().n_points();
165 
166  for (unsigned int qp=0; qp != n_qpoints; qp++)
167  for (unsigned int i=0; i != n_p_dofs; i++)
168  for (unsigned int j=0; j != n_p_dofs; j++)
169  c.get_elem_jacobian()(i,j) += JxW[qp] * (phi[j][qp]*phi[i][qp]);
170  }
171 };
172 
173 struct A2 : ElemAssembly
174 {
175  virtual void boundary_assembly(FEMContext & c)
176  {
177  if (c.has_side_boundary_id(1)) // Forcing on the horn "inlet"
178  {
179  const unsigned int p_var = 0;
180 
181  FEBase * side_fe = libmesh_nullptr;
182  c.get_side_fe(p_var, side_fe);
183 
184  const std::vector<Real> & JxW_face = side_fe->get_JxW();
185 
186  const std::vector<std::vector<Real>> & phi_face = side_fe->get_phi();
187 
188  // The number of local degrees of freedom in each variable
189  const unsigned int n_p_dofs = c.get_dof_indices(p_var).size();
190 
191  // Now we will build the affine operator
192  unsigned int n_sidepoints = c.get_side_qrule().n_points();
193 
194  for (unsigned int qp=0; qp != n_sidepoints; qp++)
195  for (unsigned int i=0; i != n_p_dofs; i++)
196  for (unsigned int j=0; j != n_p_dofs; j++)
197  c.get_elem_jacobian()(i,j) += JxW_face[qp] * phi_face[j][qp] * phi_face[i][qp];
198  }
199  }
200 };
201 
202 struct A3 : ElemAssembly
203 {
204  virtual void boundary_assembly(FEMContext & c)
205  {
206  if (c.has_side_boundary_id(2)) // Radiation condition on the "bubble"
207  {
208  const unsigned int p_var = 0;
209 
210  FEBase * side_fe = libmesh_nullptr;
211  c.get_side_fe(p_var, side_fe);
212 
213  const std::vector<Real> & JxW_face = side_fe->get_JxW();
214 
215  const std::vector<std::vector<Real>> & phi_face = side_fe->get_phi();
216 
217  // The number of local degrees of freedom in each variable
218  const unsigned int n_p_dofs = c.get_dof_indices(p_var).size();
219 
220  // Now we will build the affine operator
221  unsigned int n_sidepoints = c.get_side_qrule().n_points();
222 
223  for (unsigned int qp=0; qp != n_sidepoints; qp++)
224  for (unsigned int i=0; i != n_p_dofs; i++)
225  for (unsigned int j=0; j != n_p_dofs; j++)
226  c.get_elem_jacobian()(i,j) += JxW_face[qp] * phi_face[j][qp] * phi_face[i][qp];
227  }
228  }
229 };
230 
231 struct F0 : ElemAssembly
232 {
233  virtual void boundary_assembly(FEMContext & c)
234  {
235  if (c.has_side_boundary_id(1)) // Output is calculated on the horn "inlet"
236  {
237  const unsigned int p_var = 0;
238 
239  FEBase * side_fe = libmesh_nullptr;
240  c.get_side_fe(p_var, side_fe);
241 
242  const std::vector<Real> & JxW_face = side_fe->get_JxW();
243 
244  const std::vector<std::vector<Real>> & phi_face = side_fe->get_phi();
245 
246  // The number of local degrees of freedom in each variable
247  const unsigned int n_p_dofs = c.get_dof_indices(p_var).size();
248 
249  // Now we will build the affine operator
250  unsigned int n_sidepoints = c.get_side_qrule().n_points();
251 
252  for (unsigned int qp=0; qp != n_sidepoints; qp++)
253  for (unsigned int i=0; i != n_p_dofs; i++)
254  c.get_elem_residual()(i) += JxW_face[qp] * phi_face[i][qp];
255  }
256  }
257 };
258 
260 {
261  virtual void boundary_assembly(FEMContext & c)
262  {
263  if (c.has_side_boundary_id(1)) // Forcing on the horn "inlet"
264  {
265  const unsigned int p_var = 0;
266 
267  FEBase * side_fe = libmesh_nullptr;
268  c.get_side_fe(p_var, side_fe);
269 
270  const std::vector<Real> & JxW_face = side_fe->get_JxW();
271 
272  const std::vector<std::vector<Real>> & phi_face = side_fe->get_phi();
273 
274  // The number of local degrees of freedom in each variable
275  const unsigned int n_p_dofs = c.get_dof_indices(p_var).size();
276 
277  // Now we will build the affine operator
278  unsigned int n_sidepoints = c.get_side_qrule().n_points();
279 
280  for (unsigned int qp=0; qp != n_sidepoints; qp++)
281  for (unsigned int i=0; i != n_p_dofs; i++)
282  c.get_elem_residual()(i) += JxW_face[qp] * phi_face[i][qp];
283  }
284  }
285 };
286 
288 {
293  {
294  // set up the RBThetaExpansion object
295  attach_A_theta(&theta_a_0); // Attach the lhs theta
296  attach_A_theta(&theta_a_1);
297  attach_A_theta(&theta_a_2);
298  attach_A_theta(&theta_a_3);
299 
300  attach_F_theta(&theta_f_0); // Attach the rhs theta
301 
302  attach_output_theta(&theta_output_0);
303  }
304 
305  // The RBTheta member variables
312 };
313 
314 // Define an RBAssemblyExpansion class for this PDE
316 {
321  {
322  // And set up the RBAssemblyExpansion object
323  attach_A_assembly(&A0_assembly); // Attach the lhs assembly
324  attach_A_assembly(&A1_assembly);
325  attach_A_assembly(&A2_assembly);
326  attach_A_assembly(&A3_assembly);
327 
328  attach_F_assembly(&F0_assembly); // Attach the rhs assembly
329 
330  attach_output_assembly(&Output0_assembly);
331  }
332 
333  // The ElemAssembly objects
341 };
342 
343 #endif
344 
345 #endif
RealVectorValue RealGradient
Definition: assembly.h:202
ThetaOutput0 theta_output_0
Definition: assembly.h:311
Definition: assembly.h:98
virtual void boundary_assembly(FEMContext &c)
Perform the element boundary assembly.
Definition: assembly.h:261
Real get_value(const std::string &param_name) const
Get the value of the specific parameter.
Definition: rb_parameters.C:45
virtual Number evaluate(const RBParameters &mu)
Evaluate the functor object for the given parameter.
Definition: assembly.h:50
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 Number evaluate(const RBParameters &mu)
Evaluate the functor object for the given parameter.
Definition: assembly.h:42
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
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 Number evaluate(const RBParameters &mu)
Evaluate the functor object for the given parameter.
Definition: assembly.h:74
AcousticsRBThetaExpansion()
Constructor.
Definition: assembly.h:292
virtual Number evaluate(const RBParameters &mu)
Evaluate the functor object for the given parameter.
Definition: assembly.h:66
const DenseVector< Number > & get_elem_residual() const
Const accessor for element residual.
Definition: diff_context.h:248
This class stores the set of RBTheta functor objects that define the "parameter-dependent expansion" ...
This is the MeshBase class.
Definition: mesh_base.h:68
const DenseMatrix< Number > & get_elem_jacobian() const
Const accessor for element Jacobian.
Definition: diff_context.h:282
virtual Number evaluate(const RBParameters &)
Evaluate the functor object for the given parameter.
Definition: assembly.h:82
virtual void interior_assembly(FEMContext &c)
Perform the element interior assembly.
Definition: assembly.h:149
AcousticsRBAssemblyExpansion()
Constructor.
Definition: assembly.h:320
virtual Number evaluate(const RBParameters &mu)
Evaluate the functor object for the given parameter.
Definition: assembly.h:58
This class provides all data required for a physics package (e.g.
Definition: fem_context.h:61
This class is part of the rbOOmit framework.
Definition: rb_parameters.h:42
FEGenericBase< Real > FEBase
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.h:90
Definition: assembly.h:127
bool has_side_boundary_id(boundary_id_type id) const
Reports if the boundary id is found on the current side.
Definition: fem_context.C:209
Definition: assembly.h:39
virtual void boundary_assembly(FEMContext &c)
Perform the element boundary assembly.
Definition: assembly.h:204
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
This class stores the set of ElemAssembly functor objects that define the "parameter-independent expa...
ElemAssembly provides a per-element (interior and boundary) assembly functionality.
Definition: elem_assembly.h:38
const QBase & get_element_qrule() const
Accessor for element interior quadrature rule for the dimension of the current _elem.
Definition: fem_context.h:765
AcousticsInnerProduct acoustics_inner_product
Definition: assembly.h:340
virtual void boundary_assembly(FEMContext &c)
Perform the element boundary assembly.
Definition: assembly.h:175
virtual void interior_assembly(FEMContext &c)
Perform the element interior assembly.
Definition: assembly.h:120
This class is part of the rbOOmit framework.
Definition: rb_theta.h:46
virtual void boundary_assembly(FEMContext &c)
Perform the element boundary assembly.
Definition: assembly.h:233
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
Definition: assembly.h:69
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
unsigned int n_points() const
Definition: quadrature.h:113