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)17void 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