libMesh
assembly.h
Go to the documentation of this file.
1 #ifndef ASSEMBLY_H
2 #define ASSEMBLY_H
3 
4 // libMesh includes
5 #include "libmesh/libmesh.h"
6 #include "libmesh/mesh.h"
7 #include "libmesh/equation_systems.h"
8 #include "libmesh/fe.h"
9 #include "libmesh/quadrature.h"
10 #include "libmesh/dof_map.h"
11 #include "libmesh/dense_matrix.h"
12 #include "libmesh/dense_vector.h"
13 #include "libmesh/fe_interface.h"
14 #include "libmesh/elem.h"
15 #include "libmesh/utility.h"
16 #include "libmesh/print_trace.h"
17 #include "libmesh/int_range.h"
18 #include "libmesh/bounding_box.h"
19 
20 // rbOOmit includes
21 #include "libmesh/rb_assembly_expansion.h"
22 #include "libmesh/rb_eim_theta.h"
23 #include "libmesh/rb_parametrized_function.h"
24 
25 // C++ includes
26 #include <cmath>
27 
28 // Bring in bits from the libMesh namespace.
29 // Just the bits we're using, since this is a header.
32 using libMesh::Number;
33 using libMesh::Point;
39 using libMesh::RBTheta;
41 using libMesh::Real;
43 using libMesh::Elem;
44 using libMesh::FEBase;
50 
52 {
53  unsigned int get_n_components() const override
54  {
55  return 1;
56  }
57 
58  virtual std::vector<Number>
59  evaluate(const RBParameters & mu,
60  const Point & p,
61  dof_id_type /*elem_id*/,
62  unsigned int /*qp*/,
63  subdomain_id_type /*subdomain_id*/,
64  const std::vector<Point> & /*p_perturb*/,
65  const std::vector<Real> & /*phi_i_qp*/) override
66  {
67  // // Old way, there is only 1 entry in the return vector
68  // Real center_x = mu.get_value("center_x");
69  // Real center_y = mu.get_value("center_y");
70  // return std::vector<Number> { std::exp(-2. * (pow<2>(center_x - p(0)) + pow<2>(center_y - p(1)))) };
71 
72  // New way, there are get_n_components() * mu.n_samples() entries in
73  // the return vector. In debug mode, we verify that the same
74  // number of samples are provided for all parameters when
75  // RBParameters::n_samples() is called.
76  std::vector<Number> ret(this->get_n_components() * mu.n_samples());
77  for (auto i : make_range(mu.n_samples()))
78  {
79  Real center_x = mu.get_sample_value("center_x", i);
80  Real center_y = mu.get_sample_value("center_y", i);
81  ret[i] = std::exp(-2. * (pow<2>(center_x - p(0)) + pow<2>(center_y - p(1))));
82  }
83  return ret;
84  }
85 };
86 
87 // A simple Theta(mu) function which just returns a constant value
88 // (hence does not depend on mu). The constant must be specified when
89 // calling the constructor.
91 {
95  ThetaConstant(Number val) : _val(val) {}
96 
104  virtual Number evaluate(const RBParameters & mu) override
105  {
106  libmesh_error_msg_if(mu.n_samples() > 1,
107  "You should only call the evaluate_vec() API when using multi-sample RBParameters objects.");
108 
109  return _val;
110  }
111 
119  virtual std::vector<Number> evaluate_vec(const std::vector<RBParameters> & mus) override
120  {
121  // Compute the number of values to be returned in the vector. For
122  // single-sample RBParameters objects, there would be mus.size()
123  // values returned in the vector. For multi-sample RBParameters
124  // objects, there are:
125  // sum_i mus[i].n_samples()
126  // total Thetas, i.e. one Theta per sample.
127  unsigned int count = 0;
128  for (const auto & mu : mus)
129  count += mu.n_samples();
130 
131  return std::vector<Number>(count, _val);
132  }
133 
134 private:
136 };
137 
138 struct A0 : ElemAssembly
139 {
140  // Assemble the Laplacian operator
141  virtual void interior_assembly(FEMContext & c)
142  {
143  const unsigned int u_var = 0;
144 
145  FEBase * elem_fe = nullptr;
146  c.get_element_fe(u_var, elem_fe);
147 
148  const std::vector<Real> & JxW = elem_fe->get_JxW();
149 
150  // The velocity shape function gradients at interior
151  // quadrature points.
152  const std::vector<std::vector<RealGradient>> & dphi = elem_fe->get_dphi();
153 
154  // The number of local degrees of freedom in each variable
155  const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();
156 
157  // Now we will build the affine operator
158  unsigned int n_qpoints = c.get_element_qrule().n_points();
159 
160  for (unsigned int qp=0; qp != n_qpoints; qp++)
161  for (unsigned int i=0; i != n_u_dofs; i++)
162  for (unsigned int j=0; j != n_u_dofs; j++)
163  c.get_elem_jacobian()(i,j) += JxW[qp] * dphi[j][qp]*dphi[i][qp];
164  }
165 };
166 
167 
169 {
170  // Use the L2 norm to find the best fit
171  virtual void interior_assembly(FEMContext & c)
172  {
173  const unsigned int u_var = 0;
174 
175  FEBase * elem_fe = nullptr;
176  c.get_element_fe(u_var, elem_fe);
177 
178  const std::vector<Real> & JxW = elem_fe->get_JxW();
179 
180  const std::vector<std::vector<Real>> & phi = elem_fe->get_phi();
181 
182  const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();
183 
184  unsigned int n_qpoints = c.get_element_qrule().n_points();
185 
186  for (unsigned int qp=0; qp != n_qpoints; qp++)
187  for (unsigned int i=0; i != n_u_dofs; i++)
188  for (unsigned int j=0; j != n_u_dofs; j++)
189  c.get_elem_jacobian()(i,j) += JxW[qp] * phi[j][qp]*phi[i][qp];
190  }
191 };
192 
194 {
195  EIM_F(RBEIMConstruction & rb_eim_con_in,
196  unsigned int basis_function_index_in) :
197  RBEIMAssembly(rb_eim_con_in,
198  basis_function_index_in)
199  {}
200 
201  virtual void interior_assembly(FEMContext & c)
202  {
203  // PDE variable number
204  const unsigned int u_var = 0;
205 
206  FEBase * elem_fe = nullptr;
207  c.get_element_fe(u_var, elem_fe);
208 
209  // EIM variable number
210  const unsigned int eim_var = 0;
211 
212  const std::vector<Real> & JxW = elem_fe->get_JxW();
213 
214  const std::vector<std::vector<Real>> & phi = elem_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  std::vector<Number> eim_values;
221  eim_var,
222  eim_values);
223 
224  for (unsigned int qp=0; qp != c.get_element_qrule().n_points(); qp++)
225  for (unsigned int i=0; i != n_u_dofs; i++)
226  c.get_elem_residual()(i) += JxW[qp] * (eim_values[qp]*phi[i][qp]);
227  }
228 };
229 
236 {
237  OutputAssembly(const BoundingBox & bbox_in) :
238  bbox(bbox_in)
239  {}
240 
241  // Output: Average value over the region [min_x,max_x]x[min_y,max_y]
242  virtual void interior_assembly(FEMContext & c)
243  {
244  const unsigned int u_var = 0;
245 
246  FEBase * fe = nullptr;
247  c.get_element_fe(u_var, fe);
248 
249  const std::vector<Real> & JxW = fe->get_JxW();
250  const std::vector<std::vector<Real>> & phi = fe->get_phi();
251 
252  // The number of local degrees of freedom in each variable
253  const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();
254 
255  // Now we will build the affine operator
256  unsigned int n_qpoints = c.get_element_qrule().n_points();
257 
258  // TODO: BoundingBox should be able to compute and return its area/volume
259  Real output_area = (bbox.max()(0) - bbox.min()(0)) * (bbox.max()(1) - bbox.min()(1));
260 
262  for (unsigned int qp=0; qp != n_qpoints; qp++)
263  for (unsigned int i=0; i != n_u_dofs; i++)
264  c.get_elem_residual()(i) += JxW[qp] * phi[i][qp] / output_area;
265  }
266 
267  // Member variables that define the output region in 2D
269 };
270 
271 // Define an RBThetaExpansion class for this PDE
273 {
278  theta_a_0(0.05),
279  output_theta(1.0)
280  {
282 
283  // Note: there are no Thetas associated with the RHS since we use
284  // an EIM approximation for the forcing term.
285 
286  // Attach an RBTheta object for each output. Here we just use the
287  // ThetaConstant class again, but this time with a value of 1.0.
290  }
291 
292  // The RBTheta member variables
295 };
296 
297 // Define an RBAssemblyExpansion class for this PDE
299 {
304  L0(BoundingBox(/*min=*/Point(-0.2, -0.2), /*max=*/Point(0.0, 0.0))),
305  L1(BoundingBox(/*min=*/Point(0.0, 0.0), /*max=*/Point(0.2, 0.2)))
306  {
308 
309  // Attach output assembly objects
312  }
313 
314  // A0 assembly object
316 
317  // Assembly objects associated with the output functionals
320 };
321 
322 #endif
RealVectorValue RealGradient
const DenseMatrix< Number > & get_elem_jacobian() const
Const accessor for element Jacobian.
Definition: diff_context.h:274
bool contains_point(const Point &) const
Definition: bounding_box.C:35
const Elem & get_elem() const
Accessor for current Elem object.
Definition: fem_context.h:908
A simple functor class that provides a RBParameter-dependent function.
TestClass subdomain_id_type
Based on the 4-byte comment warning above, this probably doesn&#39;t work with exodusII at all...
Definition: id_types.h:43
This class is part of the rbOOmit framework.
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
void evaluate_basis_function(dof_id_type elem_id, unsigned int var, std::vector< Number > &values)
Return the basis function values for all quadrature points for variable var on element elem_id...
ThetaConstant theta_a_0
Definition: assembly.h:293
EimTestRBAssemblyExpansion()
Constructor.
Definition: assembly.h:303
virtual Number evaluate(const RBParameters &mu) override
Evaluate theta for a single scalar-valued RBParameters object.
Definition: assembly.h:104
This class stores the set of RBTheta functor objects that define the "parameter-dependent expansion" ...
EimTestRBThetaExpansion()
Constructor.
Definition: assembly.h:277
BoundingBox bbox
Definition: assembly.h:268
This class provides functionality required to define an assembly object that arises from an "Empirica...
virtual void interior_assembly(FEMContext &c)
Perform the element interior assembly.
Definition: assembly.h:171
EIM_F(RBEIMConstruction &rb_eim_con_in, unsigned int basis_function_index_in)
Definition: assembly.h:195
ThetaConstant output_theta
Definition: assembly.h:294
T pow(const T &x)
Definition: utility.h:328
dof_id_type id() const
Definition: dof_object.h:823
Number _val
Definition: assembly.h:135
const Point & min() const
Definition: bounding_box.h:76
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.
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
OutputAssembly(const BoundingBox &bbox_in)
Definition: assembly.h:237
virtual void interior_assembly(FEMContext &c)
Perform the element interior assembly.
Definition: assembly.h:242
void attach_A_assembly(ElemAssembly *Aq_assembly)
Attach ElemAssembly object for the left-hand side (both interior and boundary assembly).
Defines a Cartesian bounding box by the two corner extremum.
Definition: bounding_box.h:40
unsigned int n_samples() const
Returns the number of samples stored for all parameters.
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
Real get_sample_value(const std::string &param_name, std::size_t sample_idx) const
Get the value of the specified parameter at the specified sample, throwing an error if it does not ex...
Definition: rb_parameters.C:96
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
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:134
virtual std::vector< Number > evaluate_vec(const std::vector< RBParameters > &mus) override
Evaluate theta for multiple mu values, each of which may have multiple "samples". ...
Definition: assembly.h:119
unsigned int get_n_components() const override
Specify the number of components in this parametrized function.
Definition: assembly.h:53
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:141
This class is part of the rbOOmit framework.
Definition: rb_theta.h:46
ThetaConstant(Number val)
Constructor.
Definition: assembly.h:95
const Point & max() const
Definition: bounding_box.h:85
virtual std::vector< Number > evaluate(const RBParameters &mu, const Point &p, dof_id_type, unsigned int, subdomain_id_type, const std::vector< Point > &, const std::vector< Real > &) override
Definition: assembly.h:59
Output assembly object which computes the average value of the solution variable inside a user-provid...
Definition: assembly.h:153
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
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
uint8_t dof_id_type
Definition: id_types.h:67
virtual void interior_assembly(FEMContext &c)
Perform the element interior assembly.
Definition: assembly.h:201
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.