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 1, the integral of the the normal 16 // derivative of the solution over part of the boundary side_postprocess(DiffContext & context)17void LaplaceSystem::side_postprocess(DiffContext & context) 18 { 19 FEMContext & c = cast_ref<FEMContext &>(context); 20 21 // First we get some references to cell-specific data that 22 // will be used to assemble the linear system. 23 FEBase * side_fe = nullptr; 24 c.get_side_fe(0, side_fe); 25 26 // Element Jacobian * quadrature weights for interior integration 27 const std::vector<Real> & JxW = side_fe->get_JxW(); 28 29 const std::vector<Point > & q_point = side_fe->get_xyz(); 30 31 const std::vector<Point> & face_normals = side_fe->get_normals(); 32 33 unsigned int n_qpoints = c.get_side_qrule().n_points(); 34 35 Number dQoI_1 = 0.; 36 37 // Loop over qp's, compute the function at each qp and add 38 // to get the QoI 39 for (unsigned int qp=0; qp != n_qpoints; qp++) 40 { 41 const Real x = q_point[qp](0); 42 const Real y = q_point[qp](1); 43 44 const Real TOL = 1.e-5; 45 46 // If on the bottom horizontal bdry (y = -1) 47 if (std::abs(y - 1.0) <= TOL && x > 0.0) 48 { 49 // Get the value of the gradient at this point 50 const Gradient grad_T = c.side_gradient(0, qp); 51 52 // Add the contribution of this qp to the integral QoI 53 dQoI_1 += JxW[qp] * (grad_T * face_normals[qp] * x * (x - 1.)); 54 } 55 56 } // end of the quadrature point qp-loop 57 58 computed_QoI[1] = computed_QoI[1] + dQoI_1; 59 } 60