libMesh
assembly.h
Go to the documentation of this file.
1 #ifndef ASSEMBLY_H
2 #define ASSEMBLY_H
3 
4 #if defined(LIBMESH_HAVE_SLEPC) && defined(LIBMESH_HAVE_GLPK)
5 
6 #include "libmesh/sparse_matrix.h"
7 #include "libmesh/numeric_vector.h"
8 #include "libmesh/dense_matrix.h"
9 #include "libmesh/dense_vector.h"
10 #include "libmesh/fe.h"
11 #include "libmesh/fe_interface.h"
12 #include "libmesh/fe_base.h"
13 #include "libmesh/elem_assembly.h"
14 #include "libmesh/quadrature_gauss.h"
15 
16 // rbOOmit includes
17 #include "libmesh/rb_theta.h"
18 #include "libmesh/rb_assembly_expansion.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 using libMesh::FEBase;
31 
32 // Functors for the parameter-dependent part of the affine decomposition of the PDE
33 // The RHS and outputs just require a constant value of 1, so use a default RBTheta object there
34 struct ThetaA0 : RBTheta { virtual Number evaluate(const RBParameters & mu) { return mu.get_value("mu_0"); } };
35 struct ThetaA1 : RBTheta { virtual Number evaluate(const RBParameters & mu) { return mu.get_value("mu_1"); } };
36 struct ThetaA2 : RBTheta { virtual Number evaluate(const RBParameters & mu) { return mu.get_value("mu_2"); } };
37 
38 struct B : ElemAssembly
39 {
40  // Assemble the H1 inner product. This will be used as the inner product
41  // for this problem.
42  virtual void interior_assembly(FEMContext & c)
43  {
44  const unsigned int u_var = 0;
45 
46  FEBase * elem_fe = libmesh_nullptr;
47  c.get_element_fe(u_var, elem_fe);
48 
49  const std::vector<Real> & JxW = elem_fe->get_JxW();
50 
51  const std::vector<std::vector<Real>> & phi = elem_fe->get_phi();
52 
53  // The velocity shape function gradients at interior
54  // quadrature points.
55  const std::vector<std::vector<RealGradient>> & dphi = elem_fe->get_dphi();
56 
57  // The number of local degrees of freedom in each variable
58  const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();
59 
60  // Now we will build the affine operator
61  unsigned int n_qpoints = c.get_element_qrule().n_points();
62 
63  for (unsigned int qp=0; qp != n_qpoints; qp++)
64  for (unsigned int i=0; i != n_u_dofs; i++)
65  for (unsigned int j=0; j != n_u_dofs; j++)
66  c.get_elem_jacobian()(i,j) += JxW[qp] * (phi[j][qp]*phi[i][qp] + dphi[j][qp]*dphi[i][qp]);
67  }
68 };
69 
70 
71 struct A0 : ElemAssembly
72 {
73  // Assemble the Laplacian operator
74  virtual void interior_assembly(FEMContext & c)
75  {
76  const unsigned int u_var = 0;
77 
78  FEBase * elem_fe = libmesh_nullptr;
79  c.get_element_fe(u_var, elem_fe);
80 
81  const std::vector<Real> & JxW = elem_fe->get_JxW();
82 
83  // The velocity shape function gradients at interior
84  // quadrature points.
85  const std::vector<std::vector<RealGradient>> & dphi = elem_fe->get_dphi();
86 
87  // The number of local degrees of freedom in each variable
88  const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();
89 
90  // Now we will build the affine operator
91  unsigned int n_qpoints = c.get_element_qrule().n_points();
92 
93  Real
94  min_x=0.,
95  max_x=0.5;
96 
97  Point centroid = c.get_elem().centroid();
98  if ((min_x <= centroid(0)) && (centroid(0) < max_x))
99  for (unsigned int qp=0; qp != n_qpoints; qp++)
100  for (unsigned int i=0; i != n_u_dofs; i++)
101  for (unsigned int j=0; j != n_u_dofs; j++)
102  c.get_elem_jacobian()(i,j) += JxW[qp] * dphi[j][qp]*dphi[i][qp];
103  }
104 };
105 
106 struct A1 : ElemAssembly
107 {
108  // Assemble the Laplacian operator
109  virtual void interior_assembly(FEMContext & c)
110  {
111  const unsigned int u_var = 0;
112 
113  FEBase * elem_fe = libmesh_nullptr;
114  c.get_element_fe(u_var, elem_fe);
115 
116  const std::vector<Real> & JxW = elem_fe->get_JxW();
117 
118  // The velocity shape function gradients at interior
119  // quadrature points.
120  const std::vector<std::vector<RealGradient>> & dphi = elem_fe->get_dphi();
121 
122  // The number of local degrees of freedom in each variable
123  const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();
124 
125  // Now we will build the affine operator
126  unsigned int n_qpoints = c.get_element_qrule().n_points();
127 
128  Real
129  min_x=0.5,
130  max_x=1.;
131 
132  Point centroid = c.get_elem().centroid();
133  if ((min_x <= centroid(0)) && (centroid(0) <= max_x))
134  for (unsigned int qp=0; qp != n_qpoints; qp++)
135  for (unsigned int i=0; i != n_u_dofs; i++)
136  for (unsigned int j=0; j != n_u_dofs; j++)
137  c.get_elem_jacobian()(i,j) += JxW[qp] * dphi[j][qp]*dphi[i][qp];
138  }
139 };
140 
141 struct A2 : ElemAssembly
142 {
143  // Convection in the x-direction
144  virtual void interior_assembly(FEMContext & c)
145  {
146  const unsigned int u_var = 0;
147 
148  FEBase * elem_fe = libmesh_nullptr;
149  c.get_element_fe(u_var, elem_fe);
150 
151  const std::vector<Real> & JxW = elem_fe->get_JxW();
152 
153  const std::vector<std::vector<Real>> & phi = elem_fe->get_phi();
154 
155  const std::vector<std::vector<RealGradient>> & dphi = elem_fe->get_dphi();
156 
157  // The number of local degrees of freedom in each variable
158  const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();
159 
160  // Now we will build the affine operator
161  unsigned int n_qpoints = c.get_element_qrule().n_points();
162 
163  for (unsigned int qp=0; qp != n_qpoints; qp++)
164  for (unsigned int i=0; i != n_u_dofs; i++)
165  for (unsigned int j=0; j != n_u_dofs; j++)
166  c.get_elem_jacobian()(i,j) -= JxW[qp] * dphi[i][qp](0)*phi[j][qp];
167  }
168 };
169 
170 
171 struct F0 : ElemAssembly
172 {
173  // Source term, 1 throughout the domain
174  virtual void interior_assembly(FEMContext & c)
175  {
176  const unsigned int u_var = 0;
177 
178  FEBase * elem_fe = libmesh_nullptr;
179  c.get_element_fe(u_var, elem_fe);
180 
181  const std::vector<Real> & JxW = elem_fe->get_JxW();
182 
183  const std::vector<std::vector<Real>> & phi = elem_fe->get_phi();
184 
185  // The number of local degrees of freedom in each variable
186  const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();
187 
188  // Now we will build the affine operator
189  unsigned int n_qpoints = c.get_element_qrule().n_points();
190 
191  for (unsigned int qp=0; qp != n_qpoints; qp++)
192  for (unsigned int i=0; i != n_u_dofs; i++)
193  c.get_elem_residual()(i) += JxW[qp] * (1.*phi[i][qp]);
194  }
195 };
196 
198 {
199  OutputAssembly(Real min_x_in, Real max_x_in,
200  Real min_y_in, Real max_y_in) :
201  min_x(min_x_in),
202  max_x(max_x_in),
203  min_y(min_y_in),
204  max_y(max_y_in)
205  {}
206 
207  // Output: Average value over the region [min_x,max_x]x[min_y,max_y]
208  virtual void interior_assembly(FEMContext & c)
209  {
210  const unsigned int u_var = 0;
211 
212  FEBase * elem_fe = libmesh_nullptr;
213  c.get_element_fe(u_var, elem_fe);
214 
215  const std::vector<Real> & JxW = elem_fe->get_JxW();
216 
217  const std::vector<std::vector<Real>> & phi = elem_fe->get_phi();
218 
219  // The number of local degrees of freedom in each variable
220  const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();
221 
222  // Now we will build the affine operator
223  unsigned int n_qpoints = c.get_element_qrule().n_points();
224 
225  Real output_area = (max_x-min_x) * (max_y-min_y);
226 
227  Point centroid = c.get_elem().centroid();
228  if ((min_x <= centroid(0)) && (centroid(0) <= max_x) &&
229  (min_y <= centroid(1)) && (centroid(1) <= max_y))
230  for (unsigned int qp=0; qp != n_qpoints; qp++)
231  for (unsigned int i=0; i != n_u_dofs; i++)
232  c.get_elem_residual()(i) += JxW[qp] * (1.*phi[i][qp]) / output_area;
233  }
234 
235  // Member variables that define the output region in 2D
236  Real min_x, max_x, min_y, max_y;
237 };
238 
239 // Define an RBThetaExpansion class for this PDE
240 struct Ex02RBThetaExpansion : RBThetaExpansion
241 {
242 
247  {
248  // set up the RBThetaExpansion object
249  attach_A_theta(&theta_a_0); // Attach the lhs theta
250  attach_A_theta(&theta_a_1);
251  attach_A_theta(&theta_a_2);
252 
253  attach_F_theta(&rb_theta); // Attach the rhs theta
254 
255  attach_output_theta(&rb_theta); // Attach output 0 theta
256  attach_output_theta(&rb_theta); // Attach output 1 theta
257  attach_output_theta(&rb_theta); // Attach output 2 theta
258  attach_output_theta(&rb_theta); // Attach output 3 theta
259  }
260 
261  // The RBTheta member variables
265  RBTheta rb_theta; // Default RBTheta object, just returns 1.
266 };
267 
268 // Define an RBAssemblyExpansion class for this PDE
269 struct Ex02RBAssemblyExpansion : RBAssemblyExpansion
270 {
271 
276  :
277  L0(0.72, 0.88, 0.72, 0.88), // We make sure these output regions conform to the mesh
278  L1(0.12, 0.28, 0.72, 0.88),
279  L2(0.12, 0.28, 0.12, 0.28),
280  L3(0.72, 0.88, 0.12, 0.28)
281  {
282  // And set up the RBAssemblyExpansion object
283  attach_A_assembly(&A0_assembly); // Attach the lhs assembly
284  attach_A_assembly(&A1_assembly);
285  attach_A_assembly(&A2_assembly);
286 
287  attach_F_assembly(&F0_assembly); // Attach the rhs assembly
288 
289  attach_output_assembly(&L0); // Attach output 0 assembly
290  attach_output_assembly(&L1); // Attach output 1 assembly
291  attach_output_assembly(&L2); // Attach output 2 assembly
292  attach_output_assembly(&L3); // Attach output 3 assembly
293  }
294 
295  // The ElemAssembly objects
305 };
306 
307 #endif // LIBMESH_HAVE_SLEPC && LIBMESH_HAVE_GLPK
308 
309 #endif
RealVectorValue RealGradient
virtual void interior_assembly(FEMContext &c)
Perform the element interior assembly.
Definition: assembly.h:42
Ex02RBAssemblyExpansion()
Constructor.
Definition: assembly.h:275
Definition: assembly.h:98
virtual Number evaluate(const RBParameters &mu)
Definition: assembly.h:35
virtual Number evaluate(const RBParameters &mu)
Definition: assembly.h:34
Ex02RBThetaExpansion()
Constructor.
Definition: assembly.h:246
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:199
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
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
Definition: assembly.h:38
const DenseMatrix< Number > & get_elem_jacobian() const
Const accessor for element Jacobian.
Definition: diff_context.h:282
OutputAssembly L3
Definition: assembly.h:304
virtual void interior_assembly(FEMContext &c)
Perform the element interior assembly.
Definition: assembly.h:109
virtual Number evaluate(const RBParameters &mu)
Definition: assembly.h:36
This class provides all data required for a physics package (e.g.
Definition: fem_context.h:61
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:208
OutputAssembly L0
Definition: assembly.h:301
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
OutputAssembly L1
Definition: assembly.h:302
const QBase & get_element_qrule() const
Accessor for element interior quadrature rule for the dimension of the current _elem.
Definition: fem_context.h:765
OutputAssembly L2
Definition: assembly.h:303
virtual void interior_assembly(FEMContext &c)
Perform the element interior assembly.
Definition: assembly.h:144
virtual Point centroid() const
Definition: elem.C:446
virtual void interior_assembly(FEMContext &c)
Perform the element interior assembly.
Definition: assembly.h:74
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:174
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