libMesh
H-qoi.C
Go to the documentation of this file.
1 #include "H-qoi.h"
2 
3 // Here we define the functions to compute the QoI (side_qoi)
4 // and supply the right hand side for the associated adjoint problem (side_qoi_derivative)
5 
6 using namespace libMesh;
7 
8 void CoupledSystemQoI::init_qoi(std::vector<Number> & sys_qoi)
9 {
10  //Only 1 qoi to worry about
11  sys_qoi.resize(1);
12 }
13 
14 // This function supplies the right hand side for the adjoint problem
15 // We only have one QoI, so we don't bother checking the qois argument
16 // to see if it was requested from us
18  const QoISet & /* qois */)
19 {
20  FEMContext & c = cast_ref<FEMContext &>(context);
21 
22  FEBase * side_fe = libmesh_nullptr;
23  c.get_side_fe(0, side_fe);
24 
25  // Element Jacobian * quadrature weights for interior integration
26  const std::vector<Real> & JxW = side_fe->get_JxW();
27 
28  // Side basis functions
29  const std::vector<std::vector<Real>> & phi = side_fe->get_phi();
30  const std::vector<Point > & q_point = side_fe->get_xyz();
31 
32  // The number of local degrees of freedom in each variable
33  const unsigned int n_u_dofs = c.get_dof_indices(1).size();
34 
37 
38  // Now we will build the element Jacobian and residual.
39  // Constructing the residual requires the solution and its
40  // gradient from the previous timestep. This must be
41  // calculated at each quadrature point by summing the
42  // solution degree-of-freedom values by the appropriate
43  // weight functions.
44  unsigned int n_qpoints = c.get_side_qrule().n_points();
45 
46  Number u = 0.;
47  Number C = 0.;
48 
49  // // If side is on outlet
50  if (c.has_side_boundary_id(2))
51  {
52  // Loop over all the qps on this side
53  for (unsigned int qp=0; qp != n_qpoints; qp++)
54  {
55  Real x = q_point[qp](0);
56 
57  // If side is on left outlet
58  if (x < 0.)
59  {
60  // Get u at the qp
61  u = c.side_value(0, qp);
62  C = c.side_value(3, qp);
63 
64  // Add the contribution from each basis function
65  for (unsigned int i=0; i != n_u_dofs; i++)
66  {
67  Qu(i) += JxW[qp] * -phi[i][qp] * C;
68  QC(i) += JxW[qp] * phi[i][qp] * -u;
69  }
70  } // end if
71 
72  } // end quadrature loop
73 
74  } // end if on outlet
75 }
76 
77 // This function computes the actual QoI
79  const QoISet & /* qois */)
80 {
81 
82  FEMContext & c = cast_ref<FEMContext &>(context);
83 
84  // First we get some references to cell-specific data that
85  // will be used to assemble the linear system.
86  FEBase * side_fe = libmesh_nullptr;
87  c.get_side_fe(0, side_fe);
88 
89  // Element Jacobian * quadrature weights for interior integration
90  const std::vector<Real> & JxW = side_fe->get_JxW();
91  const std::vector<Point> & q_point = side_fe->get_xyz();
92 
93  // Loop over qp's, compute the function at each qp and add
94  // to get the QoI
95  unsigned int n_qpoints = c.get_side_qrule().n_points();
96 
97  Number dQoI_0 = 0.;
98  Number u = 0.;
99  Number C = 0.;
100 
101  // If side is on the left outlet
102  if (c.has_side_boundary_id(2))
103  {
104  //Loop over all the qps on this side
105  for (unsigned int qp=0; qp != n_qpoints; qp++)
106  {
107  Real x = q_point[qp](0);
108 
109  if (x < 0.)
110  {
111  // Get u and C at the qp
112  u = c.side_value(0, qp);
113  C = c.side_value(3, qp);
114 
115  dQoI_0 += JxW[qp] * -u * C;
116  } // end if
117  } // end quadrature loop
118  } // end if on bdry
119 
120  c.get_qois()[0] += dQoI_0;
121 }
virtual void side_qoi(DiffContext &context, const QoISet &qois)
Does any work that needs to be done on side of elem in a quantity of interest assembly loop...
Definition: H-qoi.C:78
This class provides all data required for a physics package (e.g.
Definition: diff_context.h:54
Data structure for specifying which Quantities of Interest should be calculated in an adjoint or a pa...
Definition: qoi_set.h:45
const QBase & get_side_qrule() const
Accessor for element side quadrature rule for the dimension of the current _elem. ...
Definition: fem_context.h:772
const std::vector< dof_id_type > & get_dof_indices() const
Accessor for element dof indices.
Definition: diff_context.h:366
const class libmesh_nullptr_t libmesh_nullptr
The libMesh namespace provides an interface to certain functionality in the library.
virtual void init_qoi(std::vector< Number > &sys_qoi)
Definition: H-qoi.C:8
Defines a dense subvector for use in finite element computations.
PetscErrorCode Vec x
const std::vector< DenseVector< Number > > & get_qoi_derivatives() const
Const accessor for QoI derivatives.
Definition: diff_context.h:330
const std::vector< Number > & get_qois() const
Const accessor for QoI vector.
Definition: diff_context.h:318
This class provides all data required for a physics package (e.g.
Definition: fem_context.h:61
bool has_side_boundary_id(boundary_id_type id) const
Reports if the boundary id is found on the current side.
Definition: fem_context.C:209
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual void side_qoi_derivative(DiffContext &context, const QoISet &qois)
Does any work that needs to be done on side of elem in a quantity of interest derivative assembly loo...
Definition: H-qoi.C:17
Number side_value(unsigned int var, unsigned int qp) const
Definition: fem_context.C:575
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
unsigned int n_points() const
Definition: quadrature.h:113
This class forms the foundation from which generic finite elements may be derived.