Index: include/physics/diff_qoi.h =================================================================== --- include/physics/diff_qoi.h (revision 6054) +++ include/physics/diff_qoi.h (working copy) @@ -25,6 +25,7 @@ #include "libmesh/diff_context.h" #include "libmesh/qoi_set.h" #include "libmesh/auto_ptr.h" +#include "libmesh/parallel.h" // C++ includes @@ -62,6 +63,12 @@ virtual ~DifferentiableQoI () {} /** + * Initialize system qoi. By default, does nothing to maintain backward compatibility for FEMSystem + * applications that control qoi. + */ + virtual void init_qoi( std::vector& /*sys_qoi*/, const QoISet& /*qoi_indices*/ ){}; + + /** * Clear all the data structures associated with * the QoI. */ @@ -156,6 +163,21 @@ * Copy of this object. User should override to copy any needed state. */ virtual AutoPtr clone() =0; + + /** + * Method to populate system qoi data structure with thread-local qoi. By default, simply + * sums thread qois into system qoi. Note that this is called by TBB, so we have to receive + * the QoIContributions functor class. + */ + virtual void thread_join( std::vector& sys_qoi, const std::vector& other_qoi, + const QoISet& qoi_indices ); + + /** + * Method to populate system qoi data structure with process-local qoi. By default, simply + * sums process qois into system qoi. + */ + virtual void parallel_op( std::vector& sys_qoi, std::vector& local_qoi, + const QoISet& qoi_indices ); }; } // namespace libMesh Index: include/systems/diff_system.h =================================================================== --- include/systems/diff_system.h (revision 6054) +++ include/systems/diff_system.h (working copy) @@ -151,8 +151,10 @@ /** * Attach external QoI object. */ - void attach_qoi( DifferentiableQoI* qoi ) - { this->diff_qoi = (qoi->clone()).release(); } + void attach_qoi( DifferentiableQoI* qoi_in, const QoISet& qoi_indices ) + { this->diff_qoi = (qoi_in->clone()).release(); + // User needs to resize qoi system qoi accordingly + this->diff_qoi->init_qoi( this->qoi, qoi_indices );} /** * A pointer to the solver object we're going to use. Index: src/physics/diff_qoi.C =================================================================== --- src/physics/diff_qoi.C (revision 6054) +++ src/physics/diff_qoi.C (working copy) @@ -21,14 +21,31 @@ namespace libMesh { - - DifferentiableQoI::DifferentiableQoI () : postprocess_sides(false), assemble_qoi_sides(false) { } +void DifferentiableQoI::thread_join( std::vector& sys_qoi, const std::vector& other_qoi, + const QoISet& ) +{ + for (unsigned int i=0; i != sys_qoi.size(); ++i) + sys_qoi[i] += other_qoi[i]; + return; +} +void DifferentiableQoI::parallel_op( std::vector& sys_qoi, std::vector& local_qoi, + const QoISet& ) +{ + // Sum everything into local_qoi + Parallel::sum(local_qoi); + + // Now put into system qoi + sys_qoi = local_qoi; + + return; +} + } // namespace libMesh Index: src/systems/fem_system.C =================================================================== --- src/systems/fem_system.C (revision 6054) +++ src/systems/fem_system.C (working copy) @@ -396,11 +396,9 @@ void join (const QoIContributions& other) { - const unsigned int my_size = this->qoi.size(); - libmesh_assert_equal_to (my_size, other.qoi.size()); - - for (unsigned int i=0; i != my_size; ++i) - this->qoi[i] += other.qoi[i]; + libmesh_assert_equal_to (this->qoi.size(), other.qoi.size()); + this->_diff_qoi.thread_join( this->qoi, other.qoi, _qoi_indices ); + return; } std::vector qoi; @@ -736,10 +734,8 @@ mesh.active_local_elements_end()), qoi_contributions); - Parallel::sum(qoi_contributions.qoi); + this->diff_qoi->parallel_op( this->qoi, qoi_contributions.qoi, qoi_indices ); - this->qoi = qoi_contributions.qoi; - STOP_LOG("assemble_qoi()", "FEMSystem"); } Index: examples/adjoints/adjoints_ex2/L-qoi.C =================================================================== --- examples/adjoints/adjoints_ex2/L-qoi.C (revision 6054) +++ examples/adjoints/adjoints_ex2/L-qoi.C (working copy) @@ -2,6 +2,13 @@ using namespace libMesh; +void LaplaceQoI::init_qoi( std::vector& sys_qoi, const QoISet& /*qoi_indices*/) +{ + //Only 1 qoi to worry about + sys_qoi.resize(1); + return; +} + // 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, Index: examples/adjoints/adjoints_ex2/L-shaped.h =================================================================== --- examples/adjoints/adjoints_ex2/L-shaped.h (revision 6054) +++ examples/adjoints/adjoints_ex2/L-shaped.h (working copy) @@ -18,7 +18,7 @@ const unsigned int number) : FEMSystem(es, name, number), _fe_family("LAGRANGE"), _fe_order(1), - _analytic_jacobians(true) { qoi.resize(1); } + _analytic_jacobians(true) { } std::string & fe_family() { return _fe_family; } unsigned int & fe_order() { return _fe_order; } Index: examples/adjoints/adjoints_ex2/L-qoi.h =================================================================== --- examples/adjoints/adjoints_ex2/L-qoi.h (revision 6054) +++ examples/adjoints/adjoints_ex2/L-qoi.h (working copy) @@ -19,6 +19,8 @@ LaplaceQoI(){} virtual ~LaplaceQoI(){} + virtual void init_qoi( std::vector& sys_qoi, const QoISet& qoi_indices); + virtual void postprocess( ){} virtual void element_qoi_derivative(DiffContext &context, const QoISet & qois); Index: examples/adjoints/adjoints_ex2/adjoints_ex2.C =================================================================== --- examples/adjoints/adjoints_ex2/adjoints_ex2.C (revision 6054) +++ examples/adjoints/adjoints_ex2/adjoints_ex2.C (working copy) @@ -304,10 +304,18 @@ // Build the FEMSystem LaplaceSystem &system = equation_systems.add_system ("LaplaceSystem"); + QoISet qois; + + std::vector qoi_indices; + qoi_indices.push_back(0); + qois.add_indices(qoi_indices); + + qois.set_weight(0, 0.5); + // Put some scope here to test that the cloning is working right { LaplaceQoI qoi; - system.attach_qoi( &qoi ); + system.attach_qoi( &qoi, qois ); } // Set its parameters @@ -343,21 +351,7 @@ // Get a pointer to the primal solution vector NumericVector &primal_solution = *system.solution; - - // Declare a QoISet object, we need this object to set weights for our QoI error contributions - QoISet qois; - // Declare a qoi_indices vector, each index will correspond to a QoI - std::vector qoi_indices; - qoi_indices.push_back(0); - qoi_indices.push_back(1); - qois.add_indices(qoi_indices); - - // Set weights for each index, these will weight the contribution of each QoI in the final error - // estimate to be used for flagging elements for refinement - qois.set_weight(0, 0.5); - qois.set_weight(1, 0.5); - // A SensitivityData object to hold the qois and parameters SensitivityData sensitivities(qois, system, system.get_parameter_vector()); @@ -471,17 +465,7 @@ write_output(equation_systems, a_step, "primal"); NumericVector &primal_solution = *system.solution; - - QoISet qois; - std::vector qoi_indices; - qoi_indices.push_back(0); - qoi_indices.push_back(1); - qois.add_indices(qoi_indices); - - qois.set_weight(0, 0.5); - qois.set_weight(1, 0.5); - SensitivityData sensitivities(qois, system, system.get_parameter_vector()); system.assemble_qoi_sides = true;