1 // General libMesh includes
2 #include "libmesh/libmesh_common.h"
3 #include "libmesh/elem.h"
4 #include "libmesh/fe_base.h"
5 #include "libmesh/fem_context.h"
6 #include "libmesh/point.h"
7 #include "libmesh/quadrature.h"
8 
9 // Local includes
10 #include "L-shaped.h"
11 
12 // Bring in everything from the libMesh namespace
13 using namespace libMesh;
14 
15 // Define the postprocess function to compute QoI 0, the integral of the the solution
16 // over a subdomain
element_postprocess(DiffContext & context)17 void LaplaceSystem::element_postprocess (DiffContext & context)
18 {
19   FEMContext & c = cast_ref<FEMContext &>(context);
20 
21   FEBase * elem_fe = nullptr;
22   c.get_element_fe(0, elem_fe);
23 
24   // Element Jacobian * quadrature weights for interior integration
25   const std::vector<Real> & JxW = elem_fe->get_JxW();
26 
27   const std::vector<Point> & xyz = elem_fe->get_xyz();
28 
29   // The number of local degrees of freedom in each variable
30   unsigned int n_qpoints = c.get_element_qrule().n_points();
31 
32   // The function R = int_{omega} T dR
33   // omega is a subset of Omega (the whole domain), omega = [0.75, 1.0] x [0.0, 0.25]
34   Number dQoI_0 = 0.;
35 
36   // Loop over quadrature points
37   for (unsigned int qp = 0; qp != n_qpoints; qp++)
38     {
39       // Get co-ordinate locations of the current quadrature point
40       const Real x = xyz[qp](0);
41       const Real y = xyz[qp](1);
42 
43       // If in the sub-domain omega, add the contribution to the integral R
44       if (std::abs(x - 0.875) <= 0.125 && std::abs(y - 0.125) <= 0.125)
45         {
46           // Get the solution value at the quadrature point
47           Number T = c.interior_value(0, qp);
48 
49           // Update the elemental increment dR for each qp
50           dQoI_0 += JxW[qp] * T;
51         }
52     }
53 
54   // Update the computed value of the global functional R, by adding the contribution from this element
55   computed_QoI[0] = computed_QoI[0] + dQoI_0;
56 }
57