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 #include "libmesh/boundary_info.h"
14 
15 // rbOOmit includes
16 #include "libmesh/rb_theta.h"
17 #include "libmesh/rb_assembly_expansion.h"
18 #include "libmesh/rb_parametrized_function.h"
19 #include "libmesh/rb_eim_construction.h"
20 #include "libmesh/rb_eim_theta.h"
21 
22 // Bring in bits from the libMesh namespace.
23 // Just the bits we're using, since this is a header.
29 using libMesh::Number;
30 using libMesh::Point;
35 using libMesh::RBTheta;
41 using libMesh::Real;
43 using libMesh::Elem;
44 using libMesh::FEBase;
45 
47 {
49 };
50 
51 // The "x component" of the function we're approximating with EIM
52 struct Gx : public RBParametrizedFunction
53 {
54  virtual Number evaluate(const RBParameters & mu,
55  const Point & p,
56  const Elem &)
57  {
58  Real curvature = mu.get_value("curvature");
59  return 1. + curvature*p(0);
60  }
61 };
62 
63 // The "y component" of the function we're approximating with EIM
64 struct Gy : public RBParametrizedFunction
65 {
66  virtual Number evaluate(const RBParameters & mu,
67  const Point & p,
68  const Elem &)
69  {
70  Real curvature = mu.get_value("curvature");
71  return 1. + curvature*p(0);
72  }
73 };
74 
75 // The "z component" of the function we're approximating with EIM
76 struct Gz : public RBParametrizedFunction
77 {
78  virtual Number evaluate(const RBParameters & mu,
79  const Point & p,
80  const Elem &)
81  {
82  Real curvature = mu.get_value("curvature");
83  return 1./(1. + curvature*p(0));
84  }
85 };
86 
87 struct ThetaA0 : RBTheta
88 {
89  virtual Number evaluate(const RBParameters & mu)
90  {
91  return mu.get_value("kappa") * mu.get_value("Bi");
92  }
93 };
94 
96 {
97  virtual void boundary_assembly(FEMContext & c)
98  {
99  std::vector<boundary_id_type> bc_ids;
101  for (std::vector<boundary_id_type>::const_iterator b =
102  bc_ids.begin(); b != bc_ids.end(); ++b)
103  if (*b == 1 || *b == 2 || *b == 3 || *b == 4)
104  {
105  const unsigned int u_var = 0;
106 
107  FEBase * side_fe = nullptr;
108  c.get_side_fe(u_var, side_fe);
109 
110  const std::vector<Real> & JxW_side = side_fe->get_JxW();
111 
112  const std::vector<std::vector<Real>> & phi_side = side_fe->get_phi();
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_sidepoints = c.get_side_qrule().n_points();
119 
120  for (unsigned int qp=0; qp != n_sidepoints; 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_side[qp] * phi_side[j][qp]*phi_side[i][qp];
124 
125  break;
126  }
127  }
128 };
129 
130 struct ThetaA1 : RBTheta
131 {
132  virtual Number evaluate(const RBParameters & mu)
133  {
134  return mu.get_value("kappa") * mu.get_value("Bi") * mu.get_value("curvature");
135  }
136 };
137 
139 {
140  virtual void boundary_assembly(FEMContext & c)
141  {
142  std::vector<boundary_id_type> bc_ids;
144  for (std::vector<boundary_id_type>::const_iterator b =
145  bc_ids.begin(); b != bc_ids.end(); ++b)
146  if (*b == 1 || *b == 3) // y == -0.2, y == 0.2
147  {
148  const unsigned int u_var = 0;
149 
150  FEBase * side_fe = nullptr;
151  c.get_side_fe(u_var, side_fe);
152 
153  const std::vector<Real> & JxW_side = side_fe->get_JxW();
154 
155  const std::vector<std::vector<Real>> & phi_side = side_fe->get_phi();
156 
157  const std::vector<Point> & xyz = side_fe->get_xyz();
158 
159  // The number of local degrees of freedom in each variable
160  const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();
161 
162  // Now we will build the affine operator
163  unsigned int n_sidepoints = c.get_side_qrule().n_points();
164 
165  for (unsigned int qp=0; qp != n_sidepoints; qp++)
166  {
167  Real x_hat = xyz[qp](0);
168 
169  for (unsigned int i=0; i != n_u_dofs; i++)
170  for (unsigned int j=0; j != n_u_dofs; j++)
171  c.get_elem_jacobian()(i,j) += JxW_side[qp] * x_hat * phi_side[j][qp]*phi_side[i][qp];
172  }
173 
174  break;
175  }
176  }
177 };
178 
179 struct ThetaA2 : RBTheta {
180  virtual Number evaluate(const RBParameters & mu)
181  {
182  return 0.2*mu.get_value("kappa") * mu.get_value("Bi") * mu.get_value("curvature");
183  }
184 };
186 {
187  virtual void boundary_assembly(FEMContext & c)
188  {
189  std::vector<boundary_id_type> bc_ids;
191  for (std::vector<boundary_id_type>::const_iterator b =
192  bc_ids.begin(); b != bc_ids.end(); ++b)
193  if (*b == 2 || *b == 4) // x == 0.2, x == -0.2
194  {
195  const unsigned int u_var = 0;
196 
197  FEBase * side_fe = nullptr;
198  c.get_side_fe(u_var, side_fe);
199 
200  const std::vector<Real> & JxW_side = side_fe->get_JxW();
201 
202  const std::vector<std::vector<Real>> & phi_side = side_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_sidepoints = c.get_side_qrule().n_points();
209 
210  if (*b==2)
211  {
212  for (unsigned int qp=0; qp != n_sidepoints; qp++)
213  {
214  for (unsigned int i=0; i != n_u_dofs; i++)
215  for (unsigned int j=0; j != n_u_dofs; j++)
216  c.get_elem_jacobian()(i,j) += JxW_side[qp] * phi_side[j][qp]*phi_side[i][qp];
217  }
218  }
219 
220  if (*b==4)
221  {
222  for (unsigned int qp=0; qp != n_sidepoints; qp++)
223  {
224  for (unsigned int i=0; i != n_u_dofs; i++)
225  for (unsigned int j=0; j != n_u_dofs; j++)
226  c.get_elem_jacobian()(i,j) -= JxW_side[qp] * phi_side[j][qp]*phi_side[i][qp];
227  }
228  }
229  }
230  }
231 };
232 
234 {
235  ThetaEIM(RBEIMEvaluation & rb_eim_eval_in, unsigned int index_in) :
236  RBEIMTheta(rb_eim_eval_in, index_in)
237  {}
238 
239  virtual Number evaluate(const RBParameters & mu)
240  {
241  return mu.get_value("kappa") * RBEIMTheta::evaluate(mu);
242  }
243 };
244 
246 {
248  unsigned int basis_function_index_in) :
249  RBEIMAssembly(rb_eim_con_in,
250  basis_function_index_in)
251  {}
252 
253  virtual void interior_assembly(FEMContext & c)
254  {
255  // PDE variable numbers
256  const unsigned int u_var = 0;
257 
258  // EIM variable numbers
259  const unsigned int Gx_var = 0;
260  const unsigned int Gy_var = 1;
261  const unsigned int Gz_var = 2;
262 
263  FEBase * elem_fe = nullptr;
264  c.get_element_fe(u_var, elem_fe);
265 
266  const std::vector<Real> & JxW = elem_fe->get_JxW();
267 
268  const std::vector<std::vector<RealGradient>> & dphi = elem_fe->get_dphi();
269 
270  // The number of local degrees of freedom in each variable
271  const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();
272 
273  std::vector<Number> eim_values_Gx;
275  c.get_elem(),
277  eim_values_Gx);
278 
279  std::vector<Number> eim_values_Gy;
281  c.get_elem(),
283  eim_values_Gy);
284 
285  std::vector<Number> eim_values_Gz;
287  c.get_elem(),
289  eim_values_Gz);
290 
291  for (unsigned int qp=0; qp != c.get_element_qrule().n_points(); qp++)
292  for (unsigned int i=0; i != n_u_dofs; i++)
293  for (unsigned int j=0; j != n_u_dofs; j++)
294  c.get_elem_jacobian()(i,j) += JxW[qp] * (eim_values_Gx[qp]*dphi[i][qp](0)*dphi[j][qp](0) +
295  eim_values_Gy[qp]*dphi[i][qp](1)*dphi[j][qp](1) +
296  eim_values_Gz[qp]*dphi[i][qp](2)*dphi[j][qp](2));
297  }
298 };
299 
300 
301 struct ThetaF0 : RBTheta
302 {
303  virtual Number evaluate(const RBParameters &) { return 1.; }
304 };
305 
306 struct AssemblyF0 : ElemAssembly
307 {
308 
309  virtual void interior_assembly(FEMContext & c)
310  {
311  const unsigned int u_var = 0;
312 
313  FEBase * elem_fe = nullptr;
314  c.get_element_fe(u_var, elem_fe);
315 
316  const std::vector<Real> & JxW = elem_fe->get_JxW();
317 
318  const std::vector<std::vector<Real>> & phi = elem_fe->get_phi();
319 
320  // The number of local degrees of freedom in each variable
321  const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();
322 
323  // Now we will build the affine operator
324  unsigned int n_qpoints = c.get_element_qrule().n_points();
325 
326  for (unsigned int qp=0; qp != n_qpoints; qp++)
327  for (unsigned int i=0; i != n_u_dofs; i++)
328  c.get_elem_residual()(i) += JxW[qp] * (1.*phi[i][qp]);
329  }
330 };
331 
332 struct ThetaF1 : RBTheta
333 {
334  virtual Number evaluate(const RBParameters & mu)
335  {
336  return mu.get_value("curvature");
337  }
338 };
339 
340 struct AssemblyF1 : ElemAssembly
341 {
342 
343  virtual void interior_assembly(FEMContext & c)
344  {
345  const unsigned int u_var = 0;
346 
347  FEBase * elem_fe = nullptr;
348  c.get_element_fe(u_var, elem_fe);
349 
350  const std::vector<Real> & JxW = elem_fe->get_JxW();
351 
352  const std::vector<std::vector<Real>> & phi = elem_fe->get_phi();
353 
354  const std::vector<Point> & xyz = elem_fe->get_xyz();
355 
356  // The number of local degrees of freedom in each variable
357  const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();
358 
359  // Now we will build the affine operator
360  unsigned int n_qpoints = c.get_element_qrule().n_points();
361 
362  for (unsigned int qp=0; qp != n_qpoints; qp++)
363  {
364  Real x_hat = xyz[qp](0);
365 
366  for (unsigned int i=0; i != n_u_dofs; i++)
367  c.get_elem_residual()(i) += JxW[qp] * (1.*x_hat*phi[i][qp]);
368  }
369  }
370 };
371 
373 {
374  // Assemble the Laplacian operator
375  virtual void interior_assembly(FEMContext & c)
376  {
377  const unsigned int u_var = 0;
378 
379  FEBase * elem_fe = nullptr;
380  c.get_element_fe(u_var, elem_fe);
381 
382  const std::vector<Real> & JxW = elem_fe->get_JxW();
383 
384  const std::vector<std::vector<RealGradient>> & dphi = elem_fe->get_dphi();
385 
386  // The number of local degrees of freedom in each variable
387  const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();
388 
389  // Now we will build the affine operator
390  unsigned int n_qpoints = c.get_element_qrule().n_points();
391 
392  for (unsigned int qp=0; qp != n_qpoints; qp++)
393  for (unsigned int i=0; i != n_u_dofs; i++)
394  for (unsigned int j=0; j != n_u_dofs; j++)
395  c.get_elem_jacobian()(i,j) += JxW[qp] * dphi[j][qp]*dphi[i][qp];
396  }
397 };
398 
400 {
401  // Use the L2 inner product to find the best fit
402  virtual void interior_assembly(FEMContext & c)
403  {
404  FEBase * elem_fe = nullptr;
405  c.get_element_fe(0, elem_fe);
406 
407  const std::vector<Real> & JxW = elem_fe->get_JxW();
408 
409  const std::vector<std::vector<Real>> & phi = elem_fe->get_phi();
410 
411  const unsigned int n_dofs = c.get_dof_indices().size();
412 
413  unsigned int n_qpoints = c.get_element_qrule().n_points();
414 
415  for (unsigned int qp=0; qp != n_qpoints; qp++)
416  for (unsigned int i=0; i != n_dofs; i++)
417  for (unsigned int j=0; j != n_dofs; j++)
418  {
419  c.get_elem_jacobian()(i,j) += JxW[qp] * phi[j][qp]*phi[i][qp];
420  }
421  }
422 };
423 
424 // Define an RBThetaExpansion class for this PDE
425 // The A terms depend on EIM, so we deal with them later
427 {
428 
433  {
437  attach_F_theta(&theta_f0); // Attach the rhs theta
439  }
440 
441  // The RBTheta member variables
447 };
448 
449 // Define an RBAssemblyExpansion class for this PDE
450 // The A terms depend on EIM, so we deal with them later
452 {
453 
458  {
459  // Point to the RBConstruction object
460  assembly_a0.rb_con = &rb_con;
461  assembly_a1.rb_con = &rb_con;
462  assembly_a2.rb_con = &rb_con;
463 
467  attach_F_assembly(&assembly_f0); // Attach the rhs assembly
469  }
470 
471  // The ElemAssembly objects
477 };
478 
479 #endif
Real get_value(const std::string &param_name) const
Get the value of the specific parameter.
RealVectorValue RealGradient
const DenseMatrix< Number > & get_elem_jacobian() const
Const accessor for element Jacobian.
Definition: diff_context.h:283
ThetaF1 theta_f1
Definition: assembly.h:446
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
virtual void boundary_assembly(FEMContext &c)
Perform the element boundary assembly.
Definition: assembly.h:140
const Elem & get_elem() const
Accessor for current Elem object.
Definition: fem_context.h:883
Definition: assembly.h:64
A simple functor class that provides a RBParameter-dependent function.
This class provides functionality required to define an RBTheta object that arises from an "Empirical...
Definition: rb_eim_theta.h:42
void attach_F_assembly(ElemAssembly *Fq_assembly)
Attach ElemAssembly object for the right-hand side (both interior and boundary assembly).
virtual Number evaluate(const RBParameters &mu)
Evaluate the functor object for the given parameter.
Definition: assembly.h:132
virtual Number evaluate(const RBParameters &mu)
Evaluate the functor object for the given parameter.
Definition: assembly.h:89
This class is part of the rbOOmit framework.
virtual Number evaluate(const RBParameters &mu, const Point &p, const Elem &)
Evaluate this parametrized function for the parameter value mu at the point p.
Definition: assembly.h:54
This is the base class from which all geometric element types are derived.
Definition: elem.h:101
ThetaF0 theta_f0
Definition: assembly.h:445
virtual void interior_assembly(FEMContext &c)
Perform the element interior assembly.
Definition: assembly.h:402
virtual void evaluate_basis_function(unsigned int var, const Elem &element, const std::vector< Point > &points, std::vector< Number > &values)
Evaluate variable var_number of this object&#39;s EIM basis function at the points points, where the points are in reference coordinates.
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:131
Ex6ThetaExpansion()
Constructor.
Definition: assembly.h:432
const MeshBase & get_mesh() const
Definition: system.h:2067
virtual Number evaluate(const RBParameters &mu, const Point &p, const Elem &)
Evaluate this parametrized function for the parameter value mu at the point p.
Definition: assembly.h:66
This class stores the set of RBTheta functor objects that define the "parameter-dependent expansion" ...
AssemblyF1 assembly_f1
Definition: assembly.h:476
virtual void boundary_assembly(FEMContext &c)
Perform the element boundary assembly.
Definition: assembly.h:187
Definition: assembly.h:76
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.
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:309
void attach_A_assembly(ElemAssembly *Aq_assembly)
Attach ElemAssembly object for the left-hand side (both interior and boundary assembly).
int8_t boundary_id_type
Definition: id_types.h:51
ThetaEIM(RBEIMEvaluation &rb_eim_eval_in, unsigned int index_in)
Definition: assembly.h:235
virtual void interior_assembly(FEMContext &c)
Perform the element interior assembly.
Definition: assembly.h:253
Defines a dense submatrix for use in Finite Element-type computations.
virtual Number evaluate(const RBParameters &mu)
Evaluate the functor object for the given parameter.
Definition: assembly.h:180
This class provides all data required for a physics package (e.g.
Definition: fem_context.h:61
const std::vector< dof_id_type > & get_dof_indices() const
Accessor for element dof indices.
Definition: diff_context.h:367
unsigned int n_points() const
Definition: quadrature.h:127
AssemblyF0 assembly_f0
Definition: assembly.h:475
AssemblyA0 assembly_a0
Definition: assembly.h:472
This class is part of the rbOOmit framework.
Definition: rb_parameters.h:42
virtual Number evaluate(const RBParameters &)
Evaluate the functor object for the given parameter.
Definition: assembly.h:303
FEGenericBase< Real > FEBase
This class provides an encapsulated access to all static public member functions of finite element cl...
Definition: fe_interface.h:65
Ex6AssemblyExpansion(RBConstruction &rb_con)
Constructor.
Definition: assembly.h:457
virtual Number evaluate(const RBParameters &mu)
Evaluate this RBEIMTheta object at the parameter mu.
Definition: assembly.h:239
AssemblyA2 assembly_a2
Definition: assembly.h:474
virtual void boundary_assembly(FEMContext &c)
Perform the element boundary assembly.
Definition: assembly.h:97
const DenseVector< Number > & get_elem_residual() const
Const accessor for element residual.
Definition: diff_context.h:249
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
ThetaA0 theta_a0
Definition: assembly.h:442
unsigned char side
Current side for side_* to examine.
Definition: fem_context.h:991
ThetaA1 theta_a1
Definition: assembly.h:443
const std::vector< Point > & get_points() const
Definition: quadrature.h:142
virtual Number evaluate(const RBParameters &mu)
Evaluate the functor object for the given parameter.
Definition: assembly.h:334
RBConstruction * rb_con
Definition: assembly.h:48
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
AssemblyEIM(RBEIMConstruction &rb_eim_con_in, unsigned int basis_function_index_in)
Definition: assembly.h:247
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
This class is part of the rbOOmit framework.
virtual void interior_assembly(FEMContext &c)
Perform the element interior assembly.
Definition: assembly.h:343
AssemblyA1 assembly_a1
Definition: assembly.h:473
ThetaA2 theta_a2
Definition: assembly.h:444
This class is part of the rbOOmit framework.
Definition: rb_theta.h:46
This class is part of the rbOOmit framework.
std::vector< boundary_id_type > boundary_ids(const Node *node) const
virtual void interior_assembly(FEMContext &c)
Perform the element interior assembly.
Definition: assembly.h:375
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
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.
virtual Number evaluate(const RBParameters &mu, const Point &p, const Elem &)
Evaluate this parametrized function for the parameter value mu at the point p.
Definition: assembly.h:78
Definition: assembly.h:52
const QBase & get_element_qrule() const
Accessor for element interior quadrature rule for the dimension of the current _elem.
Definition: fem_context.h:777
const QBase & get_side_qrule() const
Accessor for element side quadrature rule for the dimension of the current _elem. ...
Definition: fem_context.h:784