Index: include/physics/diff_qoi.h =================================================================== --- include/physics/diff_qoi.h (revision 6026) +++ include/physics/diff_qoi.h (working copy) @@ -143,13 +143,18 @@ const QoISet&) {} - /* + /** * Prepares the result of a build_context() call for use. * * Most FEMSystem-based problems will need to reimplement this in order to * call FE::get_*() as their particular QoI requires. */ virtual void init_context(DiffContext &) {} + + /** + * Copy of this object. User should override to copy any needed state. + */ + virtual DifferentiableQoI* clone() =0; }; } // namespace libMesh Index: include/systems/diff_system.h =================================================================== --- include/systems/diff_system.h (revision 6026) +++ include/systems/diff_system.h (working copy) @@ -132,6 +132,27 @@ * solvers, this will integrate dx/dt = F(x). */ virtual void solve (); + + /** + * Force the user to override clone for DifferentiableQoI + */ + virtual DifferentiableQoI* clone() + { libmesh_error(); + // dummy + return this; } + + /** + * Returns const reference to DifferentiableQoI object. Note that if no external + * QoI object is attached, the default is this. + */ + const DifferentiableQoI* get_qoi() + { return this->diff_qoi; } + + /** + * Attach external QoI object. + */ + void attach_qoi( DifferentiableQoI* qoi ) + { this->diff_qoi = qoi->clone(); } /** * A pointer to the solver object we're going to use. @@ -198,14 +219,16 @@ * users should create separate physics objects. */ DifferentiablePhysics *diff_physics; + +protected: /** * Pointer to object to use for quantity of interest assembly * evaluations. Defaults to \p this for backwards compatibility; in * the future users should create separate physics objects. */ - DifferentiableQoI *diff_qoi; -protected: + DifferentiableQoI* diff_qoi; + /** * Initializes the member data fields associated with * the system, so that, e.g., \p assemble() may be used. Index: src/systems/fem_system.C =================================================================== --- src/systems/fem_system.C (revision 6026) +++ src/systems/fem_system.C (working copy) @@ -288,7 +288,8 @@ * constructor to set context */ explicit - PostprocessContributions(FEMSystem &sys) : _sys(sys) {} + PostprocessContributions(FEMSystem &sys, + DifferentiableQoI &qoi) : _sys(sys), _qoi(qoi) {} /** * operator() for use with Threads::parallel_for(). @@ -309,14 +310,14 @@ if (_sys.fe_reinit_during_postprocess) _femcontext.elem_fe_reinit(); - _sys.element_postprocess(_femcontext); + _qoi.element_postprocess(_femcontext); for (_femcontext.side = 0; _femcontext.side != _femcontext.elem->n_sides(); ++_femcontext.side) { // Don't compute on non-boundary sides unless requested - if (!_sys.postprocess_sides || + if (!_qoi.postprocess_sides || (!_sys.compute_internal_sides && _femcontext.elem->neighbor(_femcontext.side) != NULL)) continue; @@ -325,7 +326,7 @@ if (_sys.fe_reinit_during_postprocess) _femcontext.side_fe_reinit(); - _sys.side_postprocess(_femcontext); + _qoi.side_postprocess(_femcontext); } } } @@ -333,6 +334,7 @@ private: FEMSystem& _sys; + DifferentiableQoI& _qoi; }; class QoIContributions @@ -342,15 +344,15 @@ * constructor to set context */ explicit - QoIContributions(FEMSystem &sys) : - qoi(sys.qoi.size(), 0.), _sys(sys) {} + QoIContributions(FEMSystem &sys, DifferentiableQoI &diff_qoi) : + qoi(sys.qoi.size(), 0.), _sys(sys), _diff_qoi(diff_qoi) {} /** * splitting constructor */ QoIContributions(const QoIContributions &other, Threads::split) : - qoi(other._sys.qoi.size(), 0.), _sys(other._sys) {} + qoi(other._sys.qoi.size(), 0.), _sys(other._sys), _diff_qoi(other._diff_qoi) {} /** * operator() for use with Threads::parallel_reduce(). @@ -369,21 +371,21 @@ _femcontext.pre_fe_reinit(_sys, el); _femcontext.elem_fe_reinit(); - _sys.element_qoi(_femcontext, _qoi_indices); + _diff_qoi.element_qoi(_femcontext, _qoi_indices); for (_femcontext.side = 0; _femcontext.side != _femcontext.elem->n_sides(); ++_femcontext.side) { // Don't compute on non-boundary sides unless requested - if (!_sys.assemble_qoi_sides || + if (!_diff_qoi.assemble_qoi_sides || (!_sys.compute_internal_sides && _femcontext.elem->neighbor(_femcontext.side) != NULL)) continue; _femcontext.side_fe_reinit(); - _sys.side_qoi(_femcontext, _qoi_indices); + _diff_qoi.side_qoi(_femcontext, _qoi_indices); } } @@ -404,6 +406,7 @@ private: FEMSystem& _sys; + DifferentiableQoI& _diff_qoi; const QoISet _qoi_indices; }; @@ -414,8 +417,9 @@ /** * constructor to set context */ - QoIDerivativeContributions(FEMSystem &sys, const QoISet& qoi_indices) : - _sys(sys), _qoi_indices(qoi_indices) {} + QoIDerivativeContributions(FEMSystem &sys, const QoISet& qoi_indices, + DifferentiableQoI &qoi ) : + _sys(sys), _qoi_indices(qoi_indices), _qoi(qoi) {} /** * operator() for use with Threads::parallel_for(). @@ -434,21 +438,21 @@ _femcontext.pre_fe_reinit(_sys, el); _femcontext.elem_fe_reinit(); - _sys.element_qoi_derivative(_femcontext, _qoi_indices); + _qoi.element_qoi_derivative(_femcontext, _qoi_indices); for (_femcontext.side = 0; _femcontext.side != _femcontext.elem->n_sides(); ++_femcontext.side) { // Don't compute on non-boundary sides unless requested - if (!_sys.assemble_qoi_sides || + if (!_qoi.assemble_qoi_sides || (!_sys.compute_internal_sides && _femcontext.elem->neighbor(_femcontext.side) != NULL)) continue; _femcontext.side_fe_reinit(); - _sys.side_qoi_derivative(_femcontext, _qoi_indices); + _qoi.side_qoi_derivative(_femcontext, _qoi_indices); } // We need some unmodified indices to use for constraining @@ -473,8 +477,8 @@ private: FEMSystem& _sys; - const QoISet& _qoi_indices; + DifferentiableQoI& _qoi; }; @@ -700,7 +704,7 @@ // Loop over every active mesh element on this processor Threads::parallel_for(elem_range.reset(mesh.active_local_elements_begin(), mesh.active_local_elements_end()), - PostprocessContributions(*this)); + PostprocessContributions(*this, *(this->diff_qoi))); STOP_LOG("postprocess()", "FEMSystem"); } @@ -723,7 +727,7 @@ // Create a non-temporary qoi_contributions object, so we can query // its results after the reduction - QoIContributions qoi_contributions(*this); + QoIContributions qoi_contributions(*this, *(this->diff_qoi)); // Loop over every active mesh element on this processor Threads::parallel_reduce(elem_range.reset(mesh.active_local_elements_begin(), @@ -756,7 +760,8 @@ // Loop over every active mesh element on this processor Threads::parallel_for(elem_range.reset(mesh.active_local_elements_begin(), mesh.active_local_elements_end()), - QoIDerivativeContributions(*this, qoi_indices)); + QoIDerivativeContributions(*this, qoi_indices, + *(this->diff_qoi))); STOP_LOG("assemble_qoi_derivative()", "FEMSystem"); } Index: examples/adjoints/adjoints_ex2/L-qoi.C =================================================================== --- examples/adjoints/adjoints_ex2/L-qoi.C (revision 0) +++ examples/adjoints/adjoints_ex2/L-qoi.C (working copy) @@ -0,0 +1,93 @@ +#include "L-qoi.h" + +using namespace libMesh; + +// We only have one QoI, so we don't bother checking the qois argument +// to see if it was requested from us +void LaplaceQoI::element_qoi (DiffContext &context, + const QoISet & /* qois */ ) + +{ + FEMContext &c = libmesh_cast_ref(context); + + // Element Jacobian * quadrature weights for interior integration + const std::vector &JxW = c.element_fe_var[0]->get_JxW(); + + const std::vector &xyz = c.element_fe_var[0]->get_xyz(); + + unsigned int n_qpoints = (c.get_element_qrule())->n_points(); + + Number dQoI_0 = 0.; + + // Loop over quadrature points + + for (unsigned int qp = 0; qp != n_qpoints; qp++) + { + // Get co-ordinate locations of the current quadrature point + const Real xf = xyz[qp](0); + const Real yf = xyz[qp](1); + + // If in the sub-domain omega, add the contribution to the integral R + if(fabs(xf - 0.875) <= 0.125 && fabs(yf - 0.125) <= 0.125) + { + // Get the solution value at the quadrature point + Number T = c.interior_value(0, qp); + + // Update the elemental increment dR for each qp + dQoI_0 += JxW[qp] * T; + } + + } + + // Update the computed value of the global functional R, by adding the contribution from this element + + c.elem_qoi[0] = c.elem_qoi[0] + dQoI_0; + +} + +// We only have one QoI, so we don't bother checking the qois argument +// to see if it was requested from us +void LaplaceQoI::element_qoi_derivative (DiffContext &context, + const QoISet & /* qois */) +{ + FEMContext &c = libmesh_cast_ref(context); + + // First we get some references to cell-specific data that + // will be used to assemble the linear system. + + // Element Jacobian * quadrature weights for interior integration + const std::vector &JxW = c.element_fe_var[0]->get_JxW(); + + // The basis functions for the element + const std::vector > &phi = c.element_fe_var[0]->get_phi(); + + // The element quadrature points + const std::vector &q_point = c.element_fe_var[0]->get_xyz(); + + // The number of local degrees of freedom in each variable + const unsigned int n_T_dofs = c.dof_indices_var[0].size(); + unsigned int n_qpoints = (c.get_element_qrule())->n_points(); + + // Fill the QoI RHS corresponding to this QoI. Since this is the 0th QoI + // we fill in the [0][i] subderivatives, i corresponding to the variable index. + // Our system has only one variable, so we only have to fill the [0][0] subderivative + DenseSubVector &Q = *c.elem_qoi_subderivatives[0][0]; + + // Loop over the qps + for (unsigned int qp=0; qp != n_qpoints; qp++) + { + const Real x = q_point[qp](0); + const Real y = q_point[qp](1); + + // If in the sub-domain over which QoI 0 is supported, add contributions + // to the adjoint rhs + if(fabs(x - 0.875) <= 0.125 && fabs(y - 0.125) <= 0.125) + { + for (unsigned int i=0; i != n_T_dofs; i++) + Q(i) += JxW[qp] *phi[i][qp] ; + } + + } // end of the quadrature point qp-loop +} + + Index: examples/adjoints/adjoints_ex2/L-qoi.h =================================================================== --- examples/adjoints/adjoints_ex2/L-qoi.h (revision 0) +++ examples/adjoints/adjoints_ex2/L-qoi.h (working copy) @@ -0,0 +1,32 @@ + +#ifndef L_QOI_H +#define L_QOI_H + +#include "libmesh/libmesh_common.h" +#include "libmesh/elem.h" +#include "libmesh/fe_base.h" +#include "libmesh/fem_context.h" +#include "libmesh/point.h" +#include "libmesh/quadrature.h" +#include "libmesh/diff_qoi.h" + +// Bring in everything from the libMesh namespace +using namespace libMesh; + +class LaplaceQoI : public DifferentiableQoI +{ +public: + LaplaceQoI(){} + virtual ~LaplaceQoI(){} + + virtual void postprocess( ){ } + + virtual void element_qoi_derivative(DiffContext &context, const QoISet & qois); + + virtual void element_qoi (DiffContext &context, const QoISet & qois); + + virtual DifferentiableQoI* clone( ) + { return new LaplaceQoI(); } + +}; +#endif // L_QOI_H Index: examples/adjoints/adjoints_ex2/L-shaped.C =================================================================== --- examples/adjoints/adjoints_ex2/L-shaped.C (revision 6026) +++ examples/adjoints/adjoints_ex2/L-shaped.C (working copy) @@ -23,9 +23,6 @@ parameters.push_back(infile("alpha_1", 1.0)); parameters.push_back(infile("alpha_2", 1.0)); - exact_QoI[0] = infile("QoI_0", 0.0); - exact_QoI[1] = infile("QoI_1", 0.0); - // Do the parent's initialization after variables are defined FEMSystem::init_data(); @@ -154,22 +151,6 @@ return compute_jacobian; } -// Override the default DiffSystem postprocess function to compute the -// approximations to the QoIs -void LaplaceSystem::postprocess() -{ - // Reset the array holding the computed QoIs - computed_QoI[0] = 0.0; - computed_QoI[1] = 0.0; - - FEMSystem::postprocess(); - - Parallel::sum(computed_QoI[0]); - - Parallel::sum(computed_QoI[1]); - -} - // The exact solution to the singular problem, // u_exact = r^(2/3)*sin(2*theta/3). We use this to set the Dirichlet boundary conditions Number LaplaceSystem::exact_solution(const Point& p)// xyz location Index: examples/adjoints/adjoints_ex2/L-shaped.h =================================================================== --- examples/adjoints/adjoints_ex2/L-shaped.h (revision 6026) +++ examples/adjoints/adjoints_ex2/L-shaped.h (working copy) @@ -23,23 +23,7 @@ std::string & fe_family() { return _fe_family; } unsigned int & fe_order() { return _fe_order; } bool & analytic_jacobians() { return _analytic_jacobians; } - - // Postprocessing function which we are going to override for this application - virtual void postprocess(void); - - Number &get_QoI_value(std::string type, unsigned int QoI_index) - { - if(type == "exact") - { - return exact_QoI[QoI_index]; - } - else - { - return computed_QoI[QoI_index]; - } - } - Number &get_parameter_value(unsigned int parameter_index) { return parameters[parameter_index]; @@ -71,16 +55,6 @@ // Constraint parts virtual bool side_constraint (bool request_jacobian, DiffContext &context); - - // Overloading the qoi function on elements - - virtual void element_qoi_derivative - (DiffContext &context, - const QoISet & qois); - - virtual void element_qoi - (DiffContext &context, - const QoISet & qois); Number exact_solution (const Point&); @@ -90,15 +64,7 @@ // The ParameterVector object that will contain pointers to // the system parameters ParameterVector parameter_vector; - - // Variables to hold the computed QoIs - Number computed_QoI[2]; - - // Variables to read in the exact QoIs from l-shaped.in - - Number exact_QoI[2]; - // The FE type to use std::string _fe_family; unsigned int _fe_order; Index: examples/adjoints/adjoints_ex2/adjoints_ex2.C =================================================================== --- examples/adjoints/adjoints_ex2/adjoints_ex2.C (revision 6026) +++ examples/adjoints/adjoints_ex2/adjoints_ex2.C (working copy) @@ -89,6 +89,7 @@ // Local includes #include "femparameters.h" #include "L-shaped.h" +#include "L-qoi.h" // Bring in everything from the libMesh namespace using namespace libMesh; @@ -303,6 +304,12 @@ // Build the FEMSystem LaplaceSystem &system = equation_systems.add_system ("LaplaceSystem"); + // Put some scope here to test that the cloning is working right + { + LaplaceQoI qoi; + system.attach_qoi( &qoi ); + } + // Set its parameters set_system_parameters(system, param); Index: examples/adjoints/adjoints_ex2/element_qoi.C =================================================================== --- examples/adjoints/adjoints_ex2/element_qoi.C (revision 6026) +++ examples/adjoints/adjoints_ex2/element_qoi.C (working copy) @@ -1,58 +0,0 @@ -// General libMesh includes -#include "libmesh/libmesh_common.h" -#include "libmesh/elem.h" -#include "libmesh/fe_base.h" -#include "libmesh/fem_context.h" -#include "libmesh/point.h" -#include "libmesh/quadrature.h" - -// Local includes -#include "L-shaped.h" - -// Bring in everything from the libMesh namespace -using namespace libMesh; - -// Function to assemble and compute the actual QoI, needed to compute parameter sensitivities - -// We only have one QoI, so we don't bother checking the qois argument -// to see if it was requested from us -void LaplaceSystem::element_qoi (DiffContext &context, - const QoISet & /* qois */ ) - -{ - FEMContext &c = libmesh_cast_ref(context); - - // Element Jacobian * quadrature weights for interior integration - const std::vector &JxW = c.element_fe_var[0]->get_JxW(); - - const std::vector &xyz = c.element_fe_var[0]->get_xyz(); - - unsigned int n_qpoints = (c.get_element_qrule())->n_points(); - - Number dQoI_0 = 0.; - - // Loop over quadrature points - - for (unsigned int qp = 0; qp != n_qpoints; qp++) - { - // Get co-ordinate locations of the current quadrature point - const Real xf = xyz[qp](0); - const Real yf = xyz[qp](1); - - // If in the sub-domain omega, add the contribution to the integral R - if(fabs(xf - 0.875) <= 0.125 && fabs(yf - 0.125) <= 0.125) - { - // Get the solution value at the quadrature point - Number T = c.interior_value(0, qp); - - // Update the elemental increment dR for each qp - dQoI_0 += JxW[qp] * T; - } - - } - - // Update the computed value of the global functional R, by adding the contribution from this element - - c.elem_qoi[0] = c.elem_qoi[0] + dQoI_0; - -} Index: examples/adjoints/adjoints_ex2/element_qoi_derivative.C =================================================================== --- examples/adjoints/adjoints_ex2/element_qoi_derivative.C (revision 6026) +++ examples/adjoints/adjoints_ex2/element_qoi_derivative.C (working copy) @@ -1,58 +0,0 @@ -// General libMesh includes -#include "libmesh/libmesh_common.h" -#include "libmesh/elem.h" -#include "libmesh/fe_base.h" -#include "libmesh/fem_context.h" -#include "libmesh/point.h" -#include "libmesh/quadrature.h" - -// Local includes -#include "L-shaped.h" - -// Bring in everything from the libMesh namespace -using namespace libMesh; - -// We only have one QoI, so we don't bother checking the qois argument -// to see if it was requested from us -void LaplaceSystem::element_qoi_derivative (DiffContext &context, - const QoISet & /* qois */) -{ - FEMContext &c = libmesh_cast_ref(context); - - // First we get some references to cell-specific data that - // will be used to assemble the linear system. - - // Element Jacobian * quadrature weights for interior integration - const std::vector &JxW = c.element_fe_var[0]->get_JxW(); - - // The basis functions for the element - const std::vector > &phi = c.element_fe_var[0]->get_phi(); - - // The element quadrature points - const std::vector &q_point = c.element_fe_var[0]->get_xyz(); - - // The number of local degrees of freedom in each variable - const unsigned int n_T_dofs = c.dof_indices_var[0].size(); - unsigned int n_qpoints = (c.get_element_qrule())->n_points(); - - // Fill the QoI RHS corresponding to this QoI. Since this is the 0th QoI - // we fill in the [0][i] subderivatives, i corresponding to the variable index. - // Our system has only one variable, so we only have to fill the [0][0] subderivative - DenseSubVector &Q = *c.elem_qoi_subderivatives[0][0]; - - // Loop over the qps - for (unsigned int qp=0; qp != n_qpoints; qp++) - { - const Real x = q_point[qp](0); - const Real y = q_point[qp](1); - - // If in the sub-domain over which QoI 0 is supported, add contributions - // to the adjoint rhs - if(fabs(x - 0.875) <= 0.125 && fabs(y - 0.125) <= 0.125) - { - for (unsigned int i=0; i != n_T_dofs; i++) - Q(i) += JxW[qp] *phi[i][qp] ; - } - - } // end of the quadrature point qp-loop -}