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 = libmesh_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 = libmesh_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 = libmesh_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 = libmesh_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  // Now we will build the affine operator
274  unsigned int n_qpoints = c.get_element_qrule().n_points();
275 
276  std::vector<Number> eim_values_Gx;
277  evaluate_basis_function(Gx_var,
278  c.get_elem(),
279  c.get_element_qrule(),
280  eim_values_Gx);
281 
282  std::vector<Number> eim_values_Gy;
283  evaluate_basis_function(Gy_var,
284  c.get_elem(),
285  c.get_element_qrule(),
286  eim_values_Gy);
287 
288  std::vector<Number> eim_values_Gz;
289  evaluate_basis_function(Gz_var,
290  c.get_elem(),
291  c.get_element_qrule(),
292  eim_values_Gz);
293 
294  for (unsigned int qp=0; qp != n_qpoints; qp++)
295  for (unsigned int i=0; i != n_u_dofs; i++)
296  for (unsigned int j=0; j != n_u_dofs; j++)
297  c.get_elem_jacobian()(i,j) += JxW[qp] * (eim_values_Gx[qp]*dphi[i][qp](0)*dphi[j][qp](0) +
298  eim_values_Gy[qp]*dphi[i][qp](1)*dphi[j][qp](1) +
299  eim_values_Gz[qp]*dphi[i][qp](2)*dphi[j][qp](2));
300  }
301 };
302 
303 
304 struct ThetaF0 : RBTheta
305 {
306  virtual Number evaluate(const RBParameters &) { return 1.; }
307 };
308 
309 struct AssemblyF0 : ElemAssembly
310 {
311 
312  virtual void interior_assembly(FEMContext & c)
313  {
314  const unsigned int u_var = 0;
315 
316  FEBase * elem_fe = libmesh_nullptr;
317  c.get_element_fe(u_var, elem_fe);
318 
319  const std::vector<Real> & JxW = elem_fe->get_JxW();
320 
321  const std::vector<std::vector<Real>> & phi = elem_fe->get_phi();
322 
323  // The number of local degrees of freedom in each variable
324  const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();
325 
326  // Now we will build the affine operator
327  unsigned int n_qpoints = c.get_element_qrule().n_points();
328 
329  for (unsigned int qp=0; qp != n_qpoints; qp++)
330  for (unsigned int i=0; i != n_u_dofs; i++)
331  c.get_elem_residual()(i) += JxW[qp] * (1.*phi[i][qp]);
332  }
333 };
334 
335 struct ThetaF1 : RBTheta
336 {
337  virtual Number evaluate(const RBParameters & mu)
338  {
339  return mu.get_value("curvature");
340  }
341 };
342 
343 struct AssemblyF1 : ElemAssembly
344 {
345 
346  virtual void interior_assembly(FEMContext & c)
347  {
348  const unsigned int u_var = 0;
349 
350  FEBase * elem_fe = libmesh_nullptr;
351  c.get_element_fe(u_var, elem_fe);
352 
353  const std::vector<Real> & JxW = elem_fe->get_JxW();
354 
355  const std::vector<std::vector<Real>> & phi = elem_fe->get_phi();
356 
357  const std::vector<Point> & xyz = elem_fe->get_xyz();
358 
359  // The number of local degrees of freedom in each variable
360  const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();
361 
362  // Now we will build the affine operator
363  unsigned int n_qpoints = c.get_element_qrule().n_points();
364 
365  for (unsigned int qp=0; qp != n_qpoints; qp++)
366  {
367  Real x_hat = xyz[qp](0);
368 
369  for (unsigned int i=0; i != n_u_dofs; i++)
370  c.get_elem_residual()(i) += JxW[qp] * (1.*x_hat*phi[i][qp]);
371  }
372  }
373 };
374 
376 {
377  // Assemble the Laplacian operator
378  virtual void interior_assembly(FEMContext & c)
379  {
380  const unsigned int u_var = 0;
381 
382  FEBase * elem_fe = libmesh_nullptr;
383  c.get_element_fe(u_var, elem_fe);
384 
385  const std::vector<Real> & JxW = elem_fe->get_JxW();
386 
387  const std::vector<std::vector<RealGradient>> & dphi = elem_fe->get_dphi();
388 
389  // The number of local degrees of freedom in each variable
390  const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();
391 
392  // Now we will build the affine operator
393  unsigned int n_qpoints = c.get_element_qrule().n_points();
394 
395  for (unsigned int qp=0; qp != n_qpoints; qp++)
396  for (unsigned int i=0; i != n_u_dofs; i++)
397  for (unsigned int j=0; j != n_u_dofs; j++)
398  c.get_elem_jacobian()(i,j) += JxW[qp] * dphi[j][qp]*dphi[i][qp];
399  }
400 };
401 
403 {
404  // Use the L2 inner product to find the best fit
405  virtual void interior_assembly(FEMContext & c)
406  {
407  FEBase * elem_fe = libmesh_nullptr;
408  c.get_element_fe(0, elem_fe);
409 
410  const std::vector<Real> & JxW = elem_fe->get_JxW();
411 
412  const std::vector<std::vector<Real>> & phi = elem_fe->get_phi();
413 
414  const unsigned int n_dofs = c.get_dof_indices().size();
415 
416  unsigned int n_qpoints = c.get_element_qrule().n_points();
417 
418  for (unsigned int qp=0; qp != n_qpoints; qp++)
419  for (unsigned int i=0; i != n_dofs; i++)
420  for (unsigned int j=0; j != n_dofs; j++)
421  {
422  c.get_elem_jacobian()(i,j) += JxW[qp] * phi[j][qp]*phi[i][qp];
423  }
424  }
425 };
426 
427 // Define an RBThetaExpansion class for this PDE
428 // The A terms depend on EIM, so we deal with them later
430 {
431 
436  {
437  attach_A_theta(&theta_a0);
438  attach_A_theta(&theta_a1);
439  attach_A_theta(&theta_a2);
440  attach_F_theta(&theta_f0); // Attach the rhs theta
441  attach_F_theta(&theta_f1);
442  }
443 
444  // The RBTheta member variables
450 };
451 
452 // Define an RBAssemblyExpansion class for this PDE
453 // The A terms depend on EIM, so we deal with them later
455 {
456 
461  {
462  // Point to the RBConstruction object
463  assembly_a0.rb_con = &rb_con;
464  assembly_a1.rb_con = &rb_con;
465  assembly_a2.rb_con = &rb_con;
466 
467  attach_A_assembly(&assembly_a0);
468  attach_A_assembly(&assembly_a1);
469  attach_A_assembly(&assembly_a2);
470  attach_F_assembly(&assembly_f0); // Attach the rhs assembly
471  attach_F_assembly(&assembly_f1);
472  }
473 
474  // The ElemAssembly objects
480 };
481 
482 #endif
RealVectorValue RealGradient
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:117
ThetaF1 theta_f1
Definition: assembly.h:449
virtual void boundary_assembly(FEMContext &c)
Perform the element boundary assembly.
Definition: assembly.h:140
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
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:132
const QBase & get_side_qrule() const
Accessor for element side quadrature rule for the dimension of the current _elem. ...
Definition: fem_context.h:772
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
const std::vector< dof_id_type > & get_dof_indices() const
Accessor for element dof indices.
Definition: diff_context.h:366
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
ThetaF0 theta_f0
Definition: assembly.h:448
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
virtual void interior_assembly(FEMContext &c)
Perform the element interior assembly.
Definition: assembly.h:405
const DenseVector< Number > & get_elem_residual() const
Const accessor for element residual.
Definition: diff_context.h:248
Ex6ThetaExpansion()
Constructor.
Definition: assembly.h:435
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
const Elem & get_elem() const
Accessor for current Elem object.
Definition: fem_context.h:871
This class stores the set of RBTheta functor objects that define the "parameter-dependent expansion" ...
AssemblyF1 assembly_f1
Definition: assembly.h:479
virtual void boundary_assembly(FEMContext &c)
Perform the element boundary assembly.
Definition: assembly.h:187
Definition: assembly.h:76
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:312
std::vector< boundary_id_type > boundary_ids(const Node *node) const
const DenseMatrix< Number > & get_elem_jacobian() const
Const accessor for element Jacobian.
Definition: diff_context.h:282
const MeshBase & get_mesh() const
Definition: system.h:2014
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
AssemblyF0 assembly_f0
Definition: assembly.h:478
AssemblyA0 assembly_a0
Definition: assembly.h:475
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:306
FEGenericBase< Real > FEBase
This class provides an encapsulated access to all static public member functions of finite element cl...
Definition: fe_interface.h:63
Ex6AssemblyExpansion(RBConstruction &rb_con)
Constructor.
Definition: assembly.h:460
virtual Number evaluate(const RBParameters &mu)
Evaluate this RBEIMTheta object at the parameter mu.
Definition: assembly.h:239
AssemblyA2 assembly_a2
Definition: assembly.h:477
virtual void boundary_assembly(FEMContext &c)
Perform the element boundary assembly.
Definition: assembly.h:97
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
ThetaA0 theta_a0
Definition: assembly.h:445
unsigned char side
Current side for side_* to examine.
Definition: fem_context.h:977
ThetaA1 theta_a1
Definition: assembly.h:446
virtual Number evaluate(const RBParameters &mu)
Evaluate the functor object for the given parameter.
Definition: assembly.h:337
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
const QBase & get_element_qrule() const
Accessor for element interior quadrature rule for the dimension of the current _elem.
Definition: fem_context.h:765
AssemblyEIM(RBEIMConstruction &rb_eim_con_in, unsigned int basis_function_index_in)
Definition: assembly.h:247
This class is part of the rbOOmit framework.
virtual void interior_assembly(FEMContext &c)
Perform the element interior assembly.
Definition: assembly.h:346
AssemblyA1 assembly_a1
Definition: assembly.h:476
ThetaA2 theta_a2
Definition: assembly.h:447
This class is part of the rbOOmit framework.
Definition: rb_theta.h:46
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
This class is part of the rbOOmit framework.
virtual void interior_assembly(FEMContext &c)
Perform the element interior assembly.
Definition: assembly.h:378
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 &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