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 // Bring in bits from the libMesh namespace.
19 // Just the bits we're using, since this is a header.
23 using libMesh::Number;
24 using libMesh::Point;
27 using libMesh::RBTheta;
29 using libMesh::Real;
31 using libMesh::FEBase;
32 
33 // Functors for the parameter-dependent part of the affine decomposition of the PDE
34 // The RHS and outputs just require a constant value of 1, so use a default RBTheta object there
35 struct ThetaA0 : RBTheta { virtual Number evaluate(const RBParameters &) { return 0.05; } };
36 struct ThetaA1 : RBTheta { virtual Number evaluate(const RBParameters & mu) { return mu.get_value("x_vel"); } };
37 struct ThetaA2 : RBTheta { virtual Number evaluate(const RBParameters & mu) { return mu.get_value("y_vel"); } };
38 
39 struct A0 : ElemAssembly
40 {
41  // Assemble the Laplacian operator
42  virtual void interior_assembly(FEMContext & c)
43  {
44  const unsigned int u_var = 0;
45 
46  FEBase * elem_fe = nullptr;
47  c.get_element_fe(u_var, elem_fe);
48 
49  const std::vector<Real> & JxW = elem_fe->get_JxW();
50 
51  // The velocity shape function gradients at interior
52  // quadrature points.
53  const std::vector<std::vector<RealGradient>> & dphi = elem_fe->get_dphi();
54 
55  // The number of local degrees of freedom in each variable
56  const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();
57 
58  // Now we will build the affine operator
59  unsigned int n_qpoints = c.get_element_qrule().n_points();
60 
61  for (unsigned int qp=0; qp != n_qpoints; qp++)
62  for (unsigned int i=0; i != n_u_dofs; i++)
63  for (unsigned int j=0; j != n_u_dofs; j++)
64  c.get_elem_jacobian()(i,j) += JxW[qp] * dphi[j][qp]*dphi[i][qp];
65  }
66 };
67 
68 
69 struct A1 : ElemAssembly
70 {
71  // Convection in the x-direction
72  virtual void interior_assembly(FEMContext & c)
73  {
74  const unsigned int u_var = 0;
75 
76  FEBase * elem_fe = nullptr;
77  c.get_element_fe(u_var, elem_fe);
78 
79  const std::vector<Real> & JxW = elem_fe->get_JxW();
80 
81  const std::vector<std::vector<Real>> & phi = elem_fe->get_phi();
82 
83  const std::vector<std::vector<RealGradient>> & dphi = elem_fe->get_dphi();
84 
85  // The number of local degrees of freedom in each variable
86  const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();
87 
88  // Now we will build the affine operator
89  unsigned int n_qpoints = c.get_element_qrule().n_points();
90 
91  for (unsigned int qp=0; qp != n_qpoints; qp++)
92  for (unsigned int i=0; i != n_u_dofs; i++)
93  for (unsigned int j=0; j != n_u_dofs; j++)
94  c.get_elem_jacobian()(i,j) += JxW[qp] * dphi[j][qp](0)*phi[i][qp];
95  }
96 };
97 
98 struct A2 : ElemAssembly
99 {
100  // Convection in the y-direction
101  virtual void interior_assembly(FEMContext & c)
102  {
103  const unsigned int u_var = 0;
104 
105  FEBase * elem_fe = nullptr;
106  c.get_element_fe(u_var, elem_fe);
107 
108  const std::vector<Real> & JxW = elem_fe->get_JxW();
109 
110  const std::vector<std::vector<Real>> & phi = elem_fe->get_phi();
111 
112  const std::vector<std::vector<RealGradient>> & dphi = elem_fe->get_dphi();
113 
114  // The number of local degrees of freedom in each variable
115  const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();
116 
117  // Now we will build the affine operator
118  unsigned int n_qpoints = c.get_element_qrule().n_points();
119 
120  for (unsigned int qp=0; qp != n_qpoints; qp++)
121  for (unsigned int i=0; i != n_u_dofs; i++)
122  for (unsigned int j=0; j != n_u_dofs; j++)
123  c.get_elem_jacobian()(i,j) += JxW[qp] * dphi[j][qp](1)*phi[i][qp];
124  }
125 };
126 
127 struct F0 : ElemAssembly
128 {
129  // Source term, 1 throughout the domain
130  virtual void interior_assembly(FEMContext & c)
131  {
132  const unsigned int u_var = 0;
133 
134  FEBase * elem_fe = nullptr;
135  c.get_element_fe(u_var, elem_fe);
136 
137  const std::vector<Real> & JxW = elem_fe->get_JxW();
138 
139  const std::vector<std::vector<Real>> & phi = elem_fe->get_phi();
140 
141  // The number of local degrees of freedom in each variable
142  const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();
143 
144  // Now we will build the affine operator
145  unsigned int n_qpoints = c.get_element_qrule().n_points();
146 
147  for (unsigned int qp=0; qp != n_qpoints; qp++)
148  for (unsigned int i=0; i != n_u_dofs; i++)
149  c.get_elem_residual()(i) += JxW[qp] * (1.*phi[i][qp]);
150  }
151 };
152 
154 {
155  OutputAssembly(Real min_x_in, Real max_x_in,
156  Real min_y_in, Real max_y_in)
157  :
158  min_x(min_x_in),
159  max_x(max_x_in),
160  min_y(min_y_in),
161  max_y(max_y_in)
162  {}
163 
164  // Output: Average value over the region [min_x,max_x]x[min_y,max_y]
165  virtual void interior_assembly(FEMContext & c)
166  {
167  const unsigned int u_var = 0;
168 
169  FEBase * elem_fe = nullptr;
170  c.get_element_fe(u_var, elem_fe);
171 
172  const std::vector<Real> & JxW = elem_fe->get_JxW();
173 
174  const std::vector<std::vector<Real>> & phi = elem_fe->get_phi();
175 
176  // The number of local degrees of freedom in each variable
177  const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();
178 
179  // Now we will build the affine operator
180  unsigned int n_qpoints = c.get_element_qrule().n_points();
181 
182  Real output_area = (max_x-min_x) * (max_y-min_y);
183 
184  Point avg = c.get_elem().vertex_average();
185  if ((min_x <= avg(0)) && (avg(0) <= max_x) &&
186  (min_y <= avg(1)) && (avg(1) <= max_y))
187  for (unsigned int qp=0; qp != n_qpoints; qp++)
188  for (unsigned int i=0; i != n_u_dofs; i++)
189  c.get_elem_residual()(i) += JxW[qp] * (1.*phi[i][qp]) / output_area;
190  }
191 
192  // Member variables that define the output region in 2D
194 };
195 
196 // Define an RBThetaExpansion class for this PDE
198 {
199 
204  {
205  // set up the RBThetaExpansion object
206  attach_A_theta(&theta_a_0); // Attach the lhs theta
209 
210  attach_F_theta(&rb_theta); // Attach the rhs theta
211 
212  attach_output_theta(&rb_theta); // Attach output 0 theta
213  attach_output_theta(&rb_theta); // Attach output 1 theta
214  attach_output_theta(&rb_theta); // Attach output 2 theta
215  attach_output_theta(&rb_theta); // Attach output 3 theta
216  }
217 
218  // The RBTheta member variables
222  RBTheta rb_theta; // Default RBTheta object, just returns 1.
223 };
224 
225 // Define an RBAssemblyExpansion class for this PDE
227 {
228 
233  L0(0.7, 0.8, 0.7, 0.8),
234  L1(0.2, 0.3, 0.7, 0.8),
235  L2(0.2, 0.3, 0.2, 0.3),
236  L3(0.7, 0.8, 0.2, 0.3)
237  {
238  // And set up the RBAssemblyExpansion object
239  attach_A_assembly(&A0_assembly); // Attach the lhs assembly
242 
243  attach_F_assembly(&F0_assembly); // Attach the rhs assembly
244 
245  attach_output_assembly(&L0); // Attach output 0 assembly
246  attach_output_assembly(&L1); // Attach output 1 assembly
247  attach_output_assembly(&L2); // Attach output 2 assembly
248  attach_output_assembly(&L3); // Attach output 3 assembly
249  }
250 
251  // The ElemAssembly objects
260 };
261 
262 #endif
Real get_value(const std::string &param_name) const
Get the value of the specified parameter, throw an error if it does not exist.
Definition: rb_parameters.C:64
RealVectorValue RealGradient
const DenseMatrix< Number > & get_elem_jacobian() const
Const accessor for element Jacobian.
Definition: diff_context.h:274
ThetaA2 theta_a_2
Definition: assembly.h:221
Definition: assembly.h:98
ThetaA1 theta_a_1
Definition: assembly.h:220
const Elem & get_elem() const
Accessor for current Elem object.
Definition: fem_context.h:908
virtual Number evaluate(const RBParameters &mu)
Evaluate the functor object for the given parameter.
Definition: assembly.h:36
CDRBAssemblyExpansion()
Constructor.
Definition: assembly.h:232
OutputAssembly(Real min_x_in, Real max_x_in, Real min_y_in, Real max_y_in)
Definition: assembly.h:155
void attach_F_assembly(ElemAssembly *Fq_assembly)
Attach ElemAssembly object for the right-hand side (both interior and boundary assembly).
This class stores the set of RBTheta functor objects that define the "parameter-dependent expansion" ...
CDRBThetaExpansion()
Constructor.
Definition: assembly.h:203
virtual void attach_F_theta(RBTheta *theta_q_f)
Attach a pointer to a functor object that defines one of the theta_q_a terms.
OutputAssembly L0
Definition: assembly.h:256
RBTheta rb_theta
Definition: assembly.h:222
OutputAssembly L1
Definition: assembly.h:257
virtual void interior_assembly(FEMContext &c)
Perform the element interior assembly.
Definition: assembly.h:72
virtual void attach_output_theta(std::vector< std::unique_ptr< RBTheta >> &theta_q_l)
Attach a vector of pointers to functor objects that define one of the outputs.
OutputAssembly L2
Definition: assembly.h:258
virtual Number evaluate(const RBParameters &mu)
Evaluate the functor object for the given parameter.
Definition: assembly.h:37
This class provides all data required for a physics package (e.g.
Definition: fem_context.h:62
const std::vector< dof_id_type > & get_dof_indices() const
Accessor for element dof indices.
Definition: diff_context.h:363
unsigned int n_points() const
Definition: quadrature.h:123
This class is part of the rbOOmit framework.
Definition: rb_parameters.h:52
FEGenericBase< Real > FEBase
This class provides an encapsulated access to all static public member functions of finite element cl...
Definition: fe_interface.h:65
virtual void interior_assembly(FEMContext &c)
Perform the element interior assembly.
Definition: assembly.h:165
void attach_A_assembly(ElemAssembly *Aq_assembly)
Attach ElemAssembly object for the left-hand side (both interior and boundary assembly).
Definition: assembly.h:127
const DenseVector< Number > & get_elem_residual() const
Const accessor for element residual.
Definition: diff_context.h:242
Definition: assembly.h:39
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
virtual void interior_assembly(FEMContext &c)
Perform the element interior assembly.
Definition: assembly.h:101
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:277
virtual void attach_output_assembly(std::vector< std::unique_ptr< ElemAssembly >> &output_assembly)
Attach ElemAssembly object for an output (both interior and boundary assembly).
virtual void interior_assembly(FEMContext &c)
Perform the element interior assembly.
Definition: assembly.h:42
This class is part of the rbOOmit framework.
Definition: rb_theta.h:46
ThetaA0 theta_a_0
Definition: assembly.h:219
OutputAssembly L3
Definition: assembly.h:259
virtual void interior_assembly(FEMContext &c)
Perform the element interior assembly.
Definition: assembly.h:130
Output assembly object which computes the average value of the solution variable inside a user-provid...
Definition: assembly.h:153
Definition: assembly.h:69
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
virtual Number evaluate(const RBParameters &)
Evaluate the functor object for the given parameter.
Definition: assembly.h:35
const QBase & get_element_qrule() const
Accessor for element interior quadrature rule for the dimension of the current _elem.
Definition: fem_context.h:802
Point vertex_average() const
Definition: elem.C:498
virtual void attach_A_theta(RBTheta *theta_q_a)
Attach a pointer to a functor object that defines one of the theta_q_a terms.