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/transient_rb_theta_expansion.h"
16 #include "libmesh/transient_rb_assembly_expansion.h"
17 
18 // Bring in bits from the libMesh namespace.
19 // Just the bits we're using, since this is a header.
24 using libMesh::Number;
25 using libMesh::Point;
27 using libMesh::RBTheta;
28 using libMesh::Real;
32 using libMesh::FEBase;
33 
34 // Functors for the parameter-dependent part of the affine decomposition of the PDE
35 // The RHS and outputs just require a constant value of 1, so use a default RBTheta object there
36 struct ThetaA0 : RBTheta { virtual Number evaluate(const RBParameters &) { return 0.05; } };
37 struct ThetaA1 : RBTheta { virtual Number evaluate(const RBParameters & mu) { return mu.get_value("x_vel"); } };
38 struct ThetaA2 : RBTheta { virtual Number evaluate(const RBParameters & mu) { return mu.get_value("y_vel"); } };
39 
40 struct M0 : ElemAssembly
41 {
42  // L2 matrix
43  virtual void interior_assembly(FEMContext & c)
44  {
45  const unsigned int u_var = 0;
46 
47  FEBase * elem_fe = libmesh_nullptr;
48  c.get_element_fe(u_var, elem_fe);
49 
50  const std::vector<Real> & JxW = elem_fe->get_JxW();
51 
52  const std::vector<std::vector<Real>> & phi = elem_fe->get_phi();
53 
54  // The number of local degrees of freedom in each variable
55  const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();
56 
57  // Now we will build the affine operator
58  unsigned int n_qpoints = c.get_element_qrule().n_points();
59 
60  for (unsigned int qp=0; qp != n_qpoints; qp++)
61  for (unsigned int i=0; i != n_u_dofs; i++)
62  for (unsigned int j=0; j != n_u_dofs; j++)
63  c.get_elem_jacobian()(i,j) += JxW[qp] *phi[j][qp]*phi[i][qp];
64  }
65 };
66 
67 struct A0 : ElemAssembly
68 {
69  // Assemble the Laplacian operator
70  virtual void interior_assembly(FEMContext & c)
71  {
72  const unsigned int u_var = 0;
73 
74  FEBase * elem_fe = libmesh_nullptr;
75  c.get_element_fe(u_var, elem_fe);
76 
77  const std::vector<Real> & JxW = elem_fe->get_JxW();
78 
79  // The velocity shape function gradients at interior
80  // quadrature points.
81  const std::vector<std::vector<RealGradient>> & dphi = elem_fe->get_dphi();
82 
83  // The number of local degrees of freedom in each variable
84  const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();
85 
86  // Now we will build the affine operator
87  unsigned int n_qpoints = c.get_element_qrule().n_points();
88 
89  for (unsigned int qp=0; qp != n_qpoints; qp++)
90  for (unsigned int i=0; i != n_u_dofs; i++)
91  for (unsigned int j=0; j != n_u_dofs; j++)
92  c.get_elem_jacobian()(i,j) += JxW[qp] * dphi[j][qp]*dphi[i][qp];
93  }
94 };
95 
96 
97 struct A1 : ElemAssembly
98 {
99  // Convection in the x-direction
100  virtual void interior_assembly(FEMContext & c)
101  {
102  const unsigned int u_var = 0;
103 
104  FEBase * elem_fe = libmesh_nullptr;
105  c.get_element_fe(u_var, elem_fe);
106 
107  const std::vector<Real> & JxW = elem_fe->get_JxW();
108 
109  const std::vector<std::vector<Real>> & phi = elem_fe->get_phi();
110 
111  const std::vector<std::vector<RealGradient>> & dphi = elem_fe->get_dphi();
112 
113  // The number of local degrees of freedom in each variable
114  const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();
115 
116  // Now we will build the affine operator
117  unsigned int n_qpoints = c.get_element_qrule().n_points();
118 
119  for (unsigned int qp=0; qp != n_qpoints; qp++)
120  for (unsigned int i=0; i != n_u_dofs; i++)
121  for (unsigned int j=0; j != n_u_dofs; j++)
122  c.get_elem_jacobian()(i,j) += JxW[qp] * dphi[j][qp](0)*phi[i][qp];
123  }
124 };
125 
126 struct A2 : ElemAssembly
127 {
128  // Convection in the y-direction
129  virtual void interior_assembly(FEMContext & c)
130  {
131  const unsigned int u_var = 0;
132 
133  FEBase * elem_fe = libmesh_nullptr;
134  c.get_element_fe(u_var, elem_fe);
135 
136  const std::vector<Real> & JxW = elem_fe->get_JxW();
137 
138  const std::vector<std::vector<Real>> & phi = elem_fe->get_phi();
139 
140  const std::vector<std::vector<RealGradient>> & dphi = elem_fe->get_dphi();
141 
142  // The number of local degrees of freedom in each variable
143  const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();
144 
145  // Now we will build the affine operator
146  unsigned int n_qpoints = c.get_element_qrule().n_points();
147 
148  for (unsigned int qp=0; qp != n_qpoints; qp++)
149  for (unsigned int i=0; i != n_u_dofs; i++)
150  for (unsigned int j=0; j != n_u_dofs; j++)
151  c.get_elem_jacobian()(i,j) += JxW[qp] * dphi[j][qp](1)*phi[i][qp];
152  }
153 };
154 
155 struct F0 : ElemAssembly
156 {
157  // Source term, 1 throughout the domain
158  virtual void interior_assembly(FEMContext & c)
159  {
160  const unsigned int u_var = 0;
161 
162  FEBase * elem_fe = libmesh_nullptr;
163  c.get_element_fe(u_var, elem_fe);
164 
165  const std::vector<Real> & JxW = elem_fe->get_JxW();
166 
167  const std::vector<std::vector<Real>> & phi = elem_fe->get_phi();
168 
169  // The number of local degrees of freedom in each variable
170  const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();
171 
172  // Now we will build the affine operator
173  unsigned int n_qpoints = c.get_element_qrule().n_points();
174 
175  for (unsigned int qp=0; qp != n_qpoints; qp++)
176  for (unsigned int i=0; i != n_u_dofs; i++)
177  c.get_elem_residual()(i) += JxW[qp] * (1.*phi[i][qp]);
178  }
179 };
180 
182 {
183  OutputAssembly(Real min_x_in, Real max_x_in,
184  Real min_y_in, Real max_y_in)
185  :
186  min_x(min_x_in),
187  max_x(max_x_in),
188  min_y(min_y_in),
189  max_y(max_y_in)
190  {}
191 
192  // Output: Average value over the region [min_x,max_x]x[min_y,max_y]
193  virtual void interior_assembly(FEMContext & c)
194  {
195  const unsigned int u_var = 0;
196 
197  FEBase * elem_fe = libmesh_nullptr;
198  c.get_element_fe(u_var, elem_fe);
199 
200  const std::vector<Real> & JxW = elem_fe->get_JxW();
201 
202  const std::vector<std::vector<Real>> & phi = elem_fe->get_phi();
203 
204  // The number of local degrees of freedom in each variable
205  const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();
206 
207  // Now we will build the affine operator
208  unsigned int n_qpoints = c.get_element_qrule().n_points();
209 
210  Real output_area = (max_x-min_x) * (max_y-min_y);
211 
212  Point centroid = c.get_elem().centroid();
213  if ((min_x <= centroid(0)) && (centroid(0) <= max_x) &&
214  (min_y <= centroid(1)) && (centroid(1) <= max_y))
215  for (unsigned int qp=0; qp != n_qpoints; qp++)
216  for (unsigned int i=0; i != n_u_dofs; i++)
217  c.get_elem_residual()(i) += JxW[qp] * (1.*phi[i][qp]) / output_area;
218  }
219 
220  // Member variables that define the output region in 2D
221  Real min_x, max_x, min_y, max_y;
222 };
223 
224 // Define an RBThetaExpansion class for this PDE
226 {
227 
232  {
233  // set up the RBThetaExpansion object
234  attach_M_theta(&rb_theta); // Attach the time-derivative theta
235 
236  attach_A_theta(&theta_a_0); // Attach the lhs theta
237  attach_A_theta(&theta_a_1);
238  attach_A_theta(&theta_a_2);
239 
240  attach_F_theta(&rb_theta); // Attach the rhs theta
241 
242  attach_output_theta(&rb_theta); // Attach output 0 theta
243  attach_output_theta(&rb_theta); // Attach output 1 theta
244  attach_output_theta(&rb_theta); // Attach output 2 theta
245  attach_output_theta(&rb_theta); // Attach output 3 theta
246  }
247 
248  // The RBTheta member variables
249  ThetaA0 theta_a_0;
250  ThetaA1 theta_a_1;
251  ThetaA2 theta_a_2;
252  RBTheta rb_theta; // Default RBTheta object, just returns 1.
253 };
254 
255 // Define an RBAssemblyExpansion class for this PDE
257 {
258 
263  :
264  L0(0.72, 0.88, 0.72, 0.88), // We make sure these output regions conform to the mesh
265  L1(0.12, 0.28, 0.72, 0.88),
266  L2(0.12, 0.28, 0.12, 0.28),
267  L3(0.72, 0.88, 0.12, 0.28)
268  {
269  // And set up the RBAssemblyExpansion object
270  attach_M_assembly(&M0_assembly); // Attach the time-derivative assembly
271  attach_A_assembly(&A0_assembly); // Attach the lhs assembly
272  attach_A_assembly(&A1_assembly);
273  attach_A_assembly(&A2_assembly);
274 
275  attach_F_assembly(&F0_assembly); // Attach the rhs assembly
276 
277  attach_output_assembly(&L0); // Attach output 0 assembly
278  attach_output_assembly(&L1); // Attach output 1 assembly
279  attach_output_assembly(&L2); // Attach output 2 assembly
280  attach_output_assembly(&L3); // Attach output 3 assembly
281  }
282 
283  // The ElemAssembly objects
285  A0 A0_assembly;
286  A1 A1_assembly;
287  A2 A2_assembly;
288  F0 F0_assembly;
289  OutputAssembly L0;
292  OutputAssembly L3;
293 };
294 
295 #endif
RealVectorValue RealGradient
virtual void interior_assembly(FEMContext &c)
Perform the element interior assembly.
Definition: assembly.h:43
Definition: assembly.h:98
This extends RBAssemblyExpansion to provide an assembly expansion for the case of time-dependent PDEs...
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:37
CDRBAssemblyExpansion()
Constructor.
Definition: assembly.h:262
const std::vector< dof_id_type > & get_dof_indices() const
Accessor for element dof indices.
Definition: diff_context.h:366
OutputAssembly(Real min_x_in, Real max_x_in, Real min_y_in, Real max_y_in)
Definition: assembly.h:183
const class libmesh_nullptr_t libmesh_nullptr
This class allows one to associate Dirichlet boundary values with a given set of mesh boundary ids an...
This class stores the set of RBTheta functor objects that define the "parameter-dependent expansion" ...
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
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
CDRBThetaExpansion()
Constructor.
Definition: assembly.h:231
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.h:100
Definition: assembly.h:40
virtual Number evaluate(const RBParameters &mu)
Evaluate the functor object for the given parameter.
Definition: assembly.h:38
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:193
Definition: assembly.h:127
Definition: assembly.h:39
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
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
virtual void interior_assembly(FEMContext &c)
Perform the element interior assembly.
Definition: assembly.h:129
virtual Point centroid() const
Definition: elem.C:446
virtual void interior_assembly(FEMContext &c)
Perform the element interior assembly.
Definition: assembly.h:70
This class is part of the rbOOmit framework.
Definition: rb_theta.h:46
virtual void interior_assembly(FEMContext &c)
Perform the element interior assembly.
Definition: assembly.h:158
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
virtual Number evaluate(const RBParameters &)
Evaluate the functor object for the given parameter.
Definition: assembly.h:36