The source file femparameters.h with comments:

        #ifndef __fem_parameters_h__
        #define __fem_parameters_h__
        
        #include <limits>
        #include <string>
        
        #include "libmesh/libmesh_common.h"
        #include "libmesh/getpot.h"
        
Bring in everything from the libMesh namespace
        using namespace libMesh;
        
        class FEMParameters
        {
        public:
            FEMParameters() :
              domainfile("lshaped.xda"),
              coarserefinements(0),
              solver_quiet(true),
              require_residual_reduction(true),
              min_step_length(1e-5),
              max_linear_iterations(200000), max_nonlinear_iterations(20),
              relative_step_tolerance(1.e-7), relative_residual_tolerance(1.e-10),
              initial_linear_tolerance(1.e-3), minimum_linear_tolerance(TOLERANCE*TOLERANCE),
              linear_tolerance_multiplier(1.e-3),
              nelem_target(30000), global_tolerance(0.0),
              refine_fraction(0.3), coarsen_fraction(0.3), coarsen_threshold(10),
              refine_uniformly(false),
              max_adaptivesteps(1),
              indicator_type("kelly"),
              fe_family(1, "LAGRANGE"), fe_order(1, 1),
              analytic_jacobians(true), verify_analytic_jacobians(0.0),
              print_solution_norms(false), print_solutions(false),
              print_residual_norms(false), print_residuals(false),
              print_jacobian_norms(false), print_jacobians(false) {}
        
            void read(GetPot &input);
        
            unsigned int dimension;
            Real elementorder;
            std::string domainfile;
            unsigned int coarserefinements;
        
            bool solver_quiet, require_residual_reduction;
            Real min_step_length;
            unsigned int max_linear_iterations, max_nonlinear_iterations;
            Real relative_step_tolerance, relative_residual_tolerance,
        	 initial_linear_tolerance, minimum_linear_tolerance,
        	 linear_tolerance_multiplier;
        
            unsigned int nelem_target;
            Real global_tolerance;
            Real refine_fraction, coarsen_fraction, coarsen_threshold;
            bool refine_uniformly;
            unsigned int max_adaptivesteps;
        
            std::string indicator_type;
        
            std::vector<std::string> fe_family;
            std::vector<unsigned int> fe_order;
        
            bool analytic_jacobians;
            Real verify_analytic_jacobians;
        
            bool print_solution_norms, print_solutions,
                 print_residual_norms, print_residuals,
                 print_jacobian_norms, print_jacobians;
        };
        
        #endif // __fem_parameters_h__



The source file L-qoi.h with comments:

        #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 init_qoi( std::vector<Number>& sys_qoi );
        
          virtual void postprocess( ){}
        
          virtual void element_qoi_derivative(DiffContext &context, const QoISet & qois);
        
          virtual void element_qoi (DiffContext &context, const QoISet & qois);
        
          virtual AutoPtr<DifferentiableQoI> clone( ) {
            return AutoPtr<DifferentiableQoI> ( new LaplaceQoI(*this) );
          }
        
        };
        #endif // L_QOI_H



The source file L-shaped.h with comments:

        #include "libmesh/enum_fe_family.h"
        #include "libmesh/fem_system.h"
        #include "libmesh/parameter_vector.h"
        #include "libmesh/qoi_set.h"
        #include "libmesh/system.h"
        
Bring in everything from the libMesh namespace
        using namespace libMesh;
        
FEMSystem, TimeSolver and NewtonSolver will handle most tasks, but we must specify element residuals
        class LaplaceSystem : public FEMSystem
        {
        public:
Constructor
          LaplaceSystem(EquationSystems& es,
                       const std::string& name_in,
                       const unsigned int number_in)
          : FEMSystem(es, name_in, number_in),
            _fe_family("LAGRANGE"), _fe_order(1),
            _analytic_jacobians(true) { }
        
          std::string & fe_family() { return _fe_family;  }
          unsigned int & fe_order() { return _fe_order;  }
          bool & analytic_jacobians() { return _analytic_jacobians; }
        
          Number &get_parameter_value(unsigned int parameter_index)
            {
              return parameters[parameter_index];
            }
        
          ParameterVector &get_parameter_vector()
            {
              parameter_vector.resize(parameters.size());
              for(unsigned int i = 0; i != parameters.size(); ++i)
        	{
        	  parameter_vector[i] = ¶meters[i];
        	}
        
              return parameter_vector;
            }
        
          protected:
System initialization
          virtual void init_data ();
        
Context initialization
          virtual void init_context (DiffContext &context);
        
Element residual and jacobian calculations Time dependent parts
          virtual bool element_time_derivative (bool request_jacobian,
        					DiffContext &context);
        
Constraint parts
          virtual bool side_constraint (bool request_jacobian,
        				DiffContext &context);
        
          Number exact_solution (const Point&);
        
Parameters associated with the system
          std::vector<Number> parameters;
        
The ParameterVector object that will contain pointers to the system parameters
          ParameterVector parameter_vector;
        
The FE type to use
          std::string _fe_family;
          unsigned int _fe_order;
        
Calculate Jacobians analytically or not?
          bool _analytic_jacobians;
        };



The source file adjoints_ex2.C with comments:



This example solves the Laplace equation on the classic "L-shaped" domain with adaptive mesh refinement. The exact solution is u(r,\theta) = r^{2/3} * \sin ( (2/3) * \theta). We scale this exact solution by a combination of parameters, (\alpha_{1} + 2 *\alpha_{2}) to get u = (\alpha_{1} + 2 *\alpha_{2}) * r^{2/3} * \sin ( (2/3) * \theta), which again satisfies the Laplace Equation. We define the Quantity of Interest in element_qoi.C, and compute the sensitivity of the QoI to \alpha_{1} and \alpha_{2} using the adjoint sensitivity method. Since we use the adjoint capabilities of libMesh in this example, we use the DiffSystem framework. This file (main.C) contains the declaration of mesh and equation system objects, L-shaped.C contains the assembly of the system, element_qoi_derivative.C contains the RHS for the adjoint system. Postprocessing to compute the the QoIs is done in element_qoi.C

The initial mesh contains three QUAD9 elements which represent the standard quadrants I, II, and III of the domain [-1,1]x[-1,1], i.e. Element 0: [-1,0]x[ 0,1] Element 1: [ 0,1]x[ 0,1] Element 2: [-1,0]x[-1,0] The mesh is provided in the standard libMesh ASCII format file named "lshaped.xda". In addition, an input file named "general.in" is provided which allows the user to set several parameters for the solution so that the problem can be re-run without a re-compile. The solution technique employed is to have a refinement loop with a linear (forward and adjoint) solve inside followed by a refinement of the grid and projection of the solution to the new grid In the final loop iteration, there is no additional refinement after the solve. In the input file "general.in", the variable "max_adaptivesteps" controls the number of refinement steps, and "refine_fraction" / "coarsen_fraction" determine the number of elements which will be refined / coarsened at each step.

C++ includes
        #include <iostream>
        #include <iomanip>
        
General libMesh includes
        #include "libmesh/equation_systems.h"
        #include "libmesh/error_vector.h"
        #include "libmesh/mesh.h"
        #include "libmesh/mesh_refinement.h"
        #include "libmesh/newton_solver.h"
        #include "libmesh/numeric_vector.h"
        #include "libmesh/steady_solver.h"
        #include "libmesh/system_norm.h"
        
Sensitivity Calculation related includes
        #include "libmesh/parameter_vector.h"
        #include "libmesh/sensitivity_data.h"
        
Error Estimator includes
        #include "libmesh/kelly_error_estimator.h"
        #include "libmesh/patch_recovery_error_estimator.h"
        
Adjoint Related includes
        #include "libmesh/adjoint_residual_error_estimator.h"
        #include "libmesh/qoi_set.h"
        
libMesh I/O includes
        #include "libmesh/getpot.h"
        #include "libmesh/gmv_io.h"
        
Local includes
        #include "femparameters.h"
        #include "L-shaped.h"
        #include "L-qoi.h"
        
Bring in everything from the libMesh namespace
        using namespace libMesh;
        
Local function declarations

Number output files, the files are give a prefix of primal or adjoint_i depending on whether the output is the primal solution or the dual solution for the ith QoI

Write gmv output

        void write_output(EquationSystems &es,
        		  unsigned int a_step, // The adaptive step count
        		  std::string solution_type) // primal or adjoint solve
        {
          MeshBase &mesh = es.get_mesh();
        
        #ifdef LIBMESH_HAVE_GMV
        
          std::ostringstream file_name_gmv;
          file_name_gmv << solution_type
                        << ".out.gmv."
                        << std::setw(2)
                        << std::setfill('0')
                        << std::right
                        << a_step;
        
          GMVIO(mesh).write_equation_systems
            (file_name_gmv.str(), es);
        
        #endif
        }
        
Set the parameters for the nonlinear and linear solvers to be used during the simulation

        void set_system_parameters(LaplaceSystem &system, FEMParameters ¶m)
        {
Use analytical jacobians?
          system.analytic_jacobians() = param.analytic_jacobians;
        
Verify analytic jacobians against numerical ones?
          system.verify_analytic_jacobians = param.verify_analytic_jacobians;
        
Use the prescribed FE type
          system.fe_family() = param.fe_family[0];
          system.fe_order() = param.fe_order[0];
        
More desperate debugging options
          system.print_solution_norms = param.print_solution_norms;
          system.print_solutions      = param.print_solutions;
          system.print_residual_norms = param.print_residual_norms;
          system.print_residuals      = param.print_residuals;
          system.print_jacobian_norms = param.print_jacobian_norms;
          system.print_jacobians      = param.print_jacobians;
        
No transient time solver
          system.time_solver =
              AutoPtr<TimeSolver>(new SteadySolver(system));
        
Nonlinear solver options
          {
            NewtonSolver *solver = new NewtonSolver(system);
            system.time_solver->diff_solver() = AutoPtr<DiffSolver>(solver);
        
            solver->quiet                       = param.solver_quiet;
            solver->max_nonlinear_iterations    = param.max_nonlinear_iterations;
            solver->minsteplength               = param.min_step_length;
            solver->relative_step_tolerance     = param.relative_step_tolerance;
            solver->relative_residual_tolerance = param.relative_residual_tolerance;
            solver->require_residual_reduction  = param.require_residual_reduction;
            solver->linear_tolerance_multiplier = param.linear_tolerance_multiplier;
            if (system.time_solver->reduce_deltat_on_diffsolver_failure)
              {
        	solver->continue_after_max_iterations = true;
        	solver->continue_after_backtrack_failure = true;
              }
        
And the linear solver options
            solver->max_linear_iterations       = param.max_linear_iterations;
            solver->initial_linear_tolerance    = param.initial_linear_tolerance;
            solver->minimum_linear_tolerance    = param.minimum_linear_tolerance;
          }
        }
        
Build the mesh refinement object and set parameters for refining/coarsening etc

        #ifdef LIBMESH_ENABLE_AMR
        
        AutoPtr<MeshRefinement> build_mesh_refinement(MeshBase &mesh,
                                                      FEMParameters ¶m)
        {
          AutoPtr<MeshRefinement> mesh_refinement(new MeshRefinement(mesh));
          mesh_refinement->coarsen_by_parents() = true;
          mesh_refinement->absolute_global_tolerance() = param.global_tolerance;
          mesh_refinement->nelem_target()      = param.nelem_target;
          mesh_refinement->refine_fraction()   = param.refine_fraction;
          mesh_refinement->coarsen_fraction()  = param.coarsen_fraction;
          mesh_refinement->coarsen_threshold() = param.coarsen_threshold;
        
          return mesh_refinement;
        }
        
        #endif // LIBMESH_ENABLE_AMR
        
This is where we declare the error estimators to be built and used for mesh refinement. The adjoint residual estimator needs two estimators. One for the forward component of the estimate and one for the adjoint weighting factor. Here we use the Patch Recovery indicator to estimate both the forward and adjoint weights. The H1 seminorm component of the error is used as dictated by the weak form the Laplace equation.

        AutoPtr<ErrorEstimator> build_error_estimator(FEMParameters ¶m)
        {
          AutoPtr<ErrorEstimator> error_estimator;
        
          if (param.indicator_type == "kelly")
            {
              std::cout<<"Using Kelly Error Estimator"<<std::endl;
        
              error_estimator.reset(new KellyErrorEstimator);
            }
          else if (param.indicator_type == "adjoint_residual")
            {
              std::cout<<"Using Adjoint Residual Error Estimator with Patch Recovery Weights"<<std::endl<<std::endl;
        
              AdjointResidualErrorEstimator *adjoint_residual_estimator = new AdjointResidualErrorEstimator;
        
              error_estimator.reset (adjoint_residual_estimator);
        
              adjoint_residual_estimator->error_plot_suffix = "error.gmv";
        
              PatchRecoveryErrorEstimator *p1 =
        	new PatchRecoveryErrorEstimator;
              adjoint_residual_estimator->primal_error_estimator().reset(p1);
        
              PatchRecoveryErrorEstimator *p2 =
        	new PatchRecoveryErrorEstimator;
              adjoint_residual_estimator->dual_error_estimator().reset(p2);
        
              adjoint_residual_estimator->primal_error_estimator()->error_norm.set_type(0, H1_SEMINORM);
        
              adjoint_residual_estimator->dual_error_estimator()->error_norm.set_type(0, H1_SEMINORM);
            }
          else
            {
              std::cerr << "Unknown indicator_type" << std::endl;
              libmesh_error();
            }
          return error_estimator;
        }
        
The main program.
        int main (int argc, char** argv)
        {
Initialize libMesh.
          LibMeshInit init (argc, argv);
        
Skip adaptive examples on a non-adaptive libMesh build
        #ifndef LIBMESH_ENABLE_AMR
          libmesh_example_assert(false, "--enable-amr");
        #else
        
          std::cout << "Started " << argv[0] << std::endl;
        
Make sure the general input file exists, and parse it
          {
            std::ifstream i("general.in");
            if (!i)
              {
                std::cerr << '[' << libMesh::processor_id()
                          << "] Can't find general.in; exiting early."
                          << std::endl;
                libmesh_error();
              }
          }
          GetPot infile("general.in");
        
Read in parameters from the input file
          FEMParameters param;
          param.read(infile);
        
Skip this default-2D example if libMesh was compiled as 1D-only.
          libmesh_example_assert(2 <= LIBMESH_DIM, "2D support");
        
Create a mesh, with dimension to be overridden later, distributed across the default MPI communicator.
          Mesh mesh(init.comm());
        
And an object to refine it
          AutoPtr<MeshRefinement> mesh_refinement =
            build_mesh_refinement(mesh, param);
        
And an EquationSystems to run on it
          EquationSystems equation_systems (mesh);
        
          std::cout << "Reading in and building the mesh" << std::endl;
        
Read in the mesh
          mesh.read(param.domainfile.c_str());
Make all the elements of the mesh second order so we can compute with a higher order basis
          mesh.all_second_order();
        
Create a mesh refinement object to do the initial uniform refinements on the coarse grid read in from lshaped.xda
          MeshRefinement initial_uniform_refinements(mesh);
          initial_uniform_refinements.uniformly_refine(param.coarserefinements);
        
          std::cout << "Building system" << std::endl;
        
Build the FEMSystem
          LaplaceSystem &system = equation_systems.add_system<LaplaceSystem> ("LaplaceSystem");
        
          QoISet qois;
        
          std::vector<unsigned int> 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 );
          }
        
Set its parameters
          set_system_parameters(system, param);
        
          std::cout << "Initializing systems" << std::endl;
        
          equation_systems.init ();
        
Print information about the mesh and system to the screen.
          mesh.print_info();
          equation_systems.print_info();
        
          {
Adaptively solve the timestep
            unsigned int a_step = 0;
            for (; a_step != param.max_adaptivesteps; ++a_step)
              {
We can't adapt to both a tolerance and a target mesh size
                if (param.global_tolerance != 0.)
        	  libmesh_assert_equal_to (param.nelem_target, 0);
If we aren't adapting to a tolerance we need a target mesh size
                    else
                      libmesh_assert_greater (param.nelem_target, 0);
        
Solve the forward problem
                system.solve();
        
Write out the computed primal solution
                write_output(equation_systems, a_step, "primal");
        
Get a pointer to the primal solution vector
                NumericVector<Number> &primal_solution = *system.solution;
        
A SensitivityData object to hold the qois and parameters
                SensitivityData sensitivities(qois, system, system.get_parameter_vector());
        
Make sure we get the contributions to the adjoint RHS from the sides
                system.assemble_qoi_sides = true;
        
Here we solve the adjoint problem inside the adjoint_qoi_parameter_sensitivity function, so we have to set the adjoint_already_solved boolean to false
                system.set_adjoint_already_solved(false);
        
Compute the sensitivities
                system.adjoint_qoi_parameter_sensitivity(qois, system.get_parameter_vector(), sensitivities);
        
Now that we have solved the adjoint, set the adjoint_already_solved boolean to true, so we dont solve unneccesarily in the error estimator
                system.set_adjoint_already_solved(true);
        
        	GetPot infile_l_shaped("l-shaped.in");
        
        	Number sensitivity_QoI_0_0_computed = sensitivities[0][0];
        	Number sensitivity_QoI_0_0_exact = infile_l_shaped("sensitivity_0_0", 0.0);
        	Number sensitivity_QoI_0_1_computed = sensitivities[0][1];
        	Number sensitivity_QoI_0_1_exact = infile_l_shaped("sensitivity_0_1", 0.0);
        
        	std::cout << "Adaptive step " << a_step << ", we have " << mesh.n_active_elem()
                              << " active elements and "
                              << equation_systems.n_active_dofs()
        		  << " active dofs." << std::endl ;
        
        	std::cout<<"Sensitivity of QoI one to Parameter one is "<<sensitivity_QoI_0_0_computed<<std::endl;
        	std::cout<<"Sensitivity of QoI one to Parameter two is "<<sensitivity_QoI_0_1_computed<<std::endl;
        
        	std::cout<< "The relative error in sensitivity QoI_0_0 is "
        		 << std::setprecision(17)
                         << std::abs(sensitivity_QoI_0_0_computed - sensitivity_QoI_0_0_exact) /
                            std::abs(sensitivity_QoI_0_0_exact) << std::endl;
        	std::cout<< "The relative error in sensitivity QoI_0_1 is "
                         << std::setprecision(17)
                         << std::abs(sensitivity_QoI_0_1_computed - sensitivity_QoI_0_1_exact) /
                            std::abs(sensitivity_QoI_0_1_exact) << std::endl << std::endl;
        
Get a pointer to the solution vector of the adjoint problem for QoI 0
                NumericVector<Number> &dual_solution_0 = system.get_adjoint_solution(0);
        
Swap the primal and dual solutions so we can write out the adjoint solution
                primal_solution.swap(dual_solution_0);
        	write_output(equation_systems, a_step, "adjoint_0");
        
Swap back
                primal_solution.swap(dual_solution_0);
        
We have to refine either based on reaching an error tolerance or a number of elements target, which should be verified above Otherwise we flag elements by error tolerance or nelem target

Uniform refinement
                if(param.refine_uniformly)
        	  {
        	    std::cout<<"Refining Uniformly"<<std::endl<<std::endl;
        
        	    mesh_refinement->uniformly_refine(1);
        	  }
Adaptively refine based on reaching an error tolerance
                else if(param.global_tolerance >= 0. && param.nelem_target == 0.)
        	  {
Now we construct the data structures for the mesh refinement process
                    ErrorVector error;
        
Build an error estimator object
                    AutoPtr<ErrorEstimator> error_estimator =
        	      build_error_estimator(param);
        
Estimate the error in each element using the Adjoint Residual or Kelly error estimator
                    error_estimator->estimate_error(system, error);
        
        	    mesh_refinement->flag_elements_by_error_tolerance (error);
        
        	    mesh_refinement->refine_and_coarsen_elements();
        	  }
Adaptively refine based on reaching a target number of elements
                else
        	  {
Now we construct the data structures for the mesh refinement process
                    ErrorVector error;
        
Build an error estimator object
                    AutoPtr<ErrorEstimator> error_estimator =
        	      build_error_estimator(param);
        
Estimate the error in each element using the Adjoint Residual or Kelly error estimator
                    error_estimator->estimate_error(system, error);
        
        	    if (mesh.n_active_elem() >= param.nelem_target)
        	      {
        		std::cout<<"We reached the target number of elements."<<std::endl <<std::endl;
        		break;
        	      }
        
        	    mesh_refinement->flag_elements_by_nelem_target (error);
        
        	    mesh_refinement->refine_and_coarsen_elements();
        	  }
        
Dont forget to reinit the system after each adaptive refinement !
                equation_systems.reinit();
        
                std::cout << "Refined mesh to "
                          << mesh.n_active_elem()
                          << " active elements and "
                          << equation_systems.n_active_dofs()
                          << " active dofs." << std::endl;
              }
        
Do one last solve if necessary
            if (a_step == param.max_adaptivesteps)
              {
        	system.solve();
        
        	write_output(equation_systems, a_step, "primal");
        
        	NumericVector<Number> &primal_solution = *system.solution;
        
        	SensitivityData sensitivities(qois, system, system.get_parameter_vector());
        
        	system.assemble_qoi_sides = true;
        
Here we solve the adjoint problem inside the adjoint_qoi_parameter_sensitivity function, so we have to set the adjoint_already_solved boolean to false
                system.set_adjoint_already_solved(false);
        
        	system.adjoint_qoi_parameter_sensitivity(qois, system.get_parameter_vector(), sensitivities);
        
Now that we have solved the adjoint, set the adjoint_already_solved boolean to true, so we dont solve unneccesarily in the error estimator
                system.set_adjoint_already_solved(true);
        
        	GetPot infile_l_shaped("l-shaped.in");
        
        	Number sensitivity_QoI_0_0_computed = sensitivities[0][0];
        	Number sensitivity_QoI_0_0_exact = infile_l_shaped("sensitivity_0_0", 0.0);
        	Number sensitivity_QoI_0_1_computed = sensitivities[0][1];
        	Number sensitivity_QoI_0_1_exact = infile_l_shaped("sensitivity_0_1", 0.0);
        
        	std::cout << "Adaptive step " << a_step << ", we have " << mesh.n_active_elem()
        		  << " active elements and "
        		  << equation_systems.n_active_dofs()
        		  << " active dofs." << std::endl ;
        
        	std::cout<<"Sensitivity of QoI one to Parameter one is "<<sensitivity_QoI_0_0_computed<<std::endl;
        	std::cout<<"Sensitivity of QoI one to Parameter two is "<<sensitivity_QoI_0_1_computed<<std::endl;
        
        	std::cout<< "The error in sensitivity QoI_0_0 is "
                         << std::setprecision(17)
                         << std::abs(sensitivity_QoI_0_0_computed - sensitivity_QoI_0_0_exact)/sensitivity_QoI_0_0_exact << std::endl;
        	std::cout<< "The error in sensitivity QoI_0_1 is "
                         << std::setprecision(17)
                         << std::abs(sensitivity_QoI_0_1_computed - sensitivity_QoI_0_1_exact)/sensitivity_QoI_0_1_exact << std::endl << std::endl;
        
        	NumericVector<Number> &dual_solution_0 = system.get_adjoint_solution(0);
        
        	primal_solution.swap(dual_solution_0);
        	write_output(equation_systems, a_step, "adjoint_0");
        
        	primal_solution.swap(dual_solution_0);
              }
          }
        
          std::cerr << '[' << libMesh::processor_id()
                    << "] Completing output." << std::endl;
        
        #endif // #ifndef LIBMESH_ENABLE_AMR
        
All done.
          return 0;
        }



The source file femparameters.C with comments:

        #include "femparameters.h"
        
        #define GETPOT_INPUT(A) { A = input(#A, A);\
          std::string stringval = input(#A, std::string());\
          variable_names.push_back(std::string(#A "=") + stringval); };
        #define GETPOT_INT_INPUT(A) { A = input(#A, (int)A);\
          std::string stringval = input(#A, std::string());\
          variable_names.push_back(std::string(#A "=") + stringval); };
        #define GETPOT_REGISTER(A) { \
          std::string stringval = input(#A, std::string());\
          variable_names.push_back(std::string(#A "=") + stringval); };
        
        
        void FEMParameters::read(GetPot &input)
        {
            std::vector<std::string> variable_names;
        
        
            GETPOT_INT_INPUT(coarserefinements);
            GETPOT_INPUT(domainfile);
        
            GETPOT_INPUT(solver_quiet);
            GETPOT_INPUT(require_residual_reduction);
            GETPOT_INPUT(min_step_length);
            GETPOT_INT_INPUT(max_linear_iterations);
            GETPOT_INT_INPUT(max_nonlinear_iterations);
            GETPOT_INPUT(relative_step_tolerance);
            GETPOT_INPUT(relative_residual_tolerance);
            GETPOT_INPUT(initial_linear_tolerance);
            GETPOT_INPUT(minimum_linear_tolerance);
            GETPOT_INPUT(linear_tolerance_multiplier);
            GETPOT_INT_INPUT(nelem_target);
            GETPOT_INPUT(global_tolerance);
            GETPOT_INPUT(refine_fraction);
            GETPOT_INPUT(coarsen_fraction);
            GETPOT_INPUT(coarsen_threshold);
            GETPOT_INT_INPUT(max_adaptivesteps);
            GETPOT_INPUT(refine_uniformly);
            GETPOT_INPUT(indicator_type);
        
            GETPOT_REGISTER(fe_family);
            const unsigned int n_fe_family =
              std::max(1u, input.vector_variable_size("fe_family"));
            fe_family.resize(n_fe_family, "LAGRANGE");
            for (unsigned int i=0; i != n_fe_family; ++i)
              fe_family[i]              = input("fe_family", fe_family[i].c_str(), i);
            GETPOT_REGISTER(fe_order);
            const unsigned int n_fe_order =
              input.vector_variable_size("fe_order");
            fe_order.resize(n_fe_order, 1);
            for (unsigned int i=0; i != n_fe_order; ++i)
              fe_order[i]               = input("fe_order", (int)fe_order[i], i);
        
            GETPOT_INPUT(analytic_jacobians);
            GETPOT_INPUT(verify_analytic_jacobians);
            GETPOT_INPUT(print_solution_norms);
            GETPOT_INPUT(print_solutions);
            GETPOT_INPUT(print_residual_norms);
            GETPOT_INPUT(print_residuals);
            GETPOT_INPUT(print_jacobian_norms);
            GETPOT_INPUT(print_jacobians);
        
          std::vector<std::string> bad_variables =
            input.unidentified_arguments(variable_names);
        
          if (libMesh::processor_id() == 0 && !bad_variables.empty())
            {
              std::cerr << "ERROR: Unrecognized variables:" << std::endl;
              for (unsigned int i = 0; i != bad_variables.size(); ++i)
                std::cerr << bad_variables[i] << std::endl;
              std::cerr << "not found among recognized variables." << std::endl;
              for (unsigned int i = 0; i != variable_names.size(); ++i)
                std::cerr << variable_names[i] << std::endl;
              libmesh_error();
            }
         }



The source file L-qoi.C with comments:

        #include "L-qoi.h"
        
        using namespace libMesh;
        
        void LaplaceQoI::init_qoi( std::vector<Number>& sys_qoi )
        {
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,
        			      const QoISet & /* qois */ )
        
        {
          FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
        
Element Jacobian * quadrature weights for interior integration
          const std::vector<Real> &JxW = c.element_fe_var[0]->get_JxW();
        
          const std::vector<Point> &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<FEMContext&>(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<Real> &JxW = c.element_fe_var[0]->get_JxW();
        
The basis functions for the element
          const std::vector<std::vector<Real> >          &phi = c.element_fe_var[0]->get_phi();
        
The element quadrature points
          const std::vector<Point > &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<Number> &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
        }



The source file L-shaped.C with comments:

        #include "libmesh/getpot.h"
        #include "libmesh/fe_base.h"
        #include "libmesh/quadrature.h"
        #include "libmesh/string_to_enum.h"
        #include "libmesh/parallel.h"
        #include "libmesh/fem_context.h"
        
Local includes
        #include "L-shaped.h"
        
Bring in everything from the libMesh namespace
        using namespace libMesh;
        
        void LaplaceSystem::init_data ()
        {
          this->add_variable ("T", static_cast<Order>(_fe_order),
                              Utility::string_to_enum<FEFamily>(_fe_family));
        
          GetPot infile("l-shaped.in");
          parameters.push_back(infile("alpha_1", 1.0));
          parameters.push_back(infile("alpha_2", 1.0));
        
Do the parent's initialization after variables are defined
          FEMSystem::init_data();
        
          this->time_evolving(0);
        }
        
        void LaplaceSystem::init_context(DiffContext &context)
        {
          FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
        
Now make sure we have requested all the data we need to build the linear system.
          c.element_fe_var[0]->get_JxW();
          c.element_fe_var[0]->get_phi();
          c.element_fe_var[0]->get_dphi();
        
          c.side_fe_var[0]->get_JxW();
          c.side_fe_var[0]->get_phi();
          c.side_fe_var[0]->get_dphi();
        }
        
        #define optassert(X) {if (!(X)) libmesh_error();}
        
Assemble the element contributions to the stiffness matrix
        bool LaplaceSystem::element_time_derivative (bool request_jacobian,
        					  DiffContext &context)
        {
Are the jacobians specified analytically ?
          bool compute_jacobian = request_jacobian && _analytic_jacobians;
        
          FEMContext &c = libmesh_cast_ref<FEMContext&>(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<Real> &JxW = c.element_fe_var[0]->get_JxW();
        
Element basis function gradients
          const std::vector<std::vector<RealGradient> > &dphi = c.element_fe_var[0]->get_dphi();
        
The number of local degrees of freedom in each variable
          const unsigned int n_T_dofs = c.dof_indices_var[0].size();
        
The subvectors and submatrices we need to fill:
          DenseSubMatrix<Number> &K = *c.elem_subjacobians[0][0];
          DenseSubVector<Number> &F = *c.elem_subresiduals[0];
        
Now we will build the element Jacobian and residual. Constructing the residual requires the solution and its gradient from the previous timestep. This must be calculated at each quadrature point by summing the solution degree-of-freedom values by the appropriate weight functions.
          unsigned int n_qpoints = (c.get_element_qrule())->n_points();
        
          for (unsigned int qp=0; qp != n_qpoints; qp++)
            {
Compute the solution gradient at the Newton iterate
              Gradient grad_T = c.interior_gradient(0, qp);
        
The residual contribution from this element
              for (unsigned int i=0; i != n_T_dofs; i++)
                F(i) += JxW[qp] * (parameters[0] + (2.*parameters[1])) * ( grad_T * dphi[i][qp] ) ;
              if (compute_jacobian)
                for (unsigned int i=0; i != n_T_dofs; i++)
                  for (unsigned int j=0; j != n_T_dofs; ++j)
The analytic jacobian
                    K(i,j) += JxW[qp] * (parameters[0] + (2.*parameters[1])) * ( dphi[i][qp] * dphi[j][qp] );
            } // end of the quadrature point qp-loop
        
          return compute_jacobian;
        }
        
Set Dirichlet bcs, side contributions to global stiffness matrix
        bool LaplaceSystem::side_constraint (bool request_jacobian,
        				  DiffContext &context)
        {
Are the jacobians specified analytically ?
          bool compute_jacobian = request_jacobian && _analytic_jacobians;
        
          FEMContext &c = libmesh_cast_ref<FEMContext&>(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<Real> &JxW = c.side_fe_var[0]->get_JxW();
        
Side basis functions
          const std::vector<std::vector<Real> > &phi = c.side_fe_var[0]->get_phi();
        
Side Quadrature points
          const std::vector<Point > &qside_point = c.side_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();
        
The subvectors and submatrices we need to fill:
          DenseSubMatrix<Number> &K = *c.elem_subjacobians[0][0];
          DenseSubVector<Number> &F = *c.elem_subresiduals[0];
        
          unsigned int n_qpoints = (c.get_side_qrule())->n_points();
        
          const Real penalty = 1./(TOLERANCE*TOLERANCE);
        
          for (unsigned int qp=0; qp != n_qpoints; qp++)
            {
Compute the solution at the old Newton iterate
              Number T = c.side_value(0, qp);
        
We get the Dirichlet bcs from the exact solution
              Number u_dirichlet = exact_solution (qside_point[qp]);
        
The residual from the boundary terms, penalize non-zero temperature
              for (unsigned int i=0; i != n_T_dofs; i++)
        	F(i) += JxW[qp] * penalty * ( T - u_dirichlet) * phi[i][qp];
              if (compute_jacobian)
        	for (unsigned int i=0; i != n_T_dofs; i++)
        	  for (unsigned int j=0; j != n_T_dofs; ++j)
The analytic jacobian
                    K(i,j) += JxW[qp] * penalty * phi[i][qp] * phi[j][qp];
        
            } // end of the quadrature point qp-loop
        
          return compute_jacobian;
        }
        
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
        {
          const Real x1 = p(0);
          const Real x2 = p(1);
        
          Real theta = atan2(x2,x1);
        
Make sure 0 <= theta <= 2*pi
          if (theta < 0)
            theta += 2. * libMesh::pi;
        
          return (parameters[0] + (2.*parameters[1])) * pow(x1*x1 + x2*x2, 1./3.)*sin(2./3.*theta);
        }



The source file femparameters.h without comments:

 
  #ifndef __fem_parameters_h__
  #define __fem_parameters_h__
  
  #include <limits>
  #include <string>
  
  #include "libmesh/libmesh_common.h"
  #include "libmesh/getpot.h"
  
  using namespace libMesh;
  
  class FEMParameters
  {
  public:
      FEMParameters() :
        domainfile("lshaped.xda"),
        coarserefinements(0),
        solver_quiet(true),
        require_residual_reduction(true),
        min_step_length(1e-5),
        max_linear_iterations(200000), max_nonlinear_iterations(20),
        relative_step_tolerance(1.e-7), relative_residual_tolerance(1.e-10),
        initial_linear_tolerance(1.e-3), minimum_linear_tolerance(TOLERANCE*TOLERANCE),
        linear_tolerance_multiplier(1.e-3),
        nelem_target(30000), global_tolerance(0.0),
        refine_fraction(0.3), coarsen_fraction(0.3), coarsen_threshold(10),
        refine_uniformly(false),
        max_adaptivesteps(1),
        indicator_type("kelly"),
        fe_family(1, "LAGRANGE"), fe_order(1, 1),
        analytic_jacobians(true), verify_analytic_jacobians(0.0),
        print_solution_norms(false), print_solutions(false),
        print_residual_norms(false), print_residuals(false),
        print_jacobian_norms(false), print_jacobians(false) {}
  
      void read(GetPot &input);
  
      unsigned int dimension;
      Real elementorder;
      std::string domainfile;
      unsigned int coarserefinements;
  
      bool solver_quiet, require_residual_reduction;
      Real min_step_length;
      unsigned int max_linear_iterations, max_nonlinear_iterations;
      Real relative_step_tolerance, relative_residual_tolerance,
  	 initial_linear_tolerance, minimum_linear_tolerance,
  	 linear_tolerance_multiplier;
  
      unsigned int nelem_target;
      Real global_tolerance;
      Real refine_fraction, coarsen_fraction, coarsen_threshold;
      bool refine_uniformly;
      unsigned int max_adaptivesteps;
  
      std::string indicator_type;
  
      std::vector<std::string> fe_family;
      std::vector<unsigned int> fe_order;
  
      bool analytic_jacobians;
      Real verify_analytic_jacobians;
  
      bool print_solution_norms, print_solutions,
           print_residual_norms, print_residuals,
           print_jacobian_norms, print_jacobians;
  };
  
  #endif // __fem_parameters_h__



The source file L-qoi.h without comments:

 
  #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"
  
  using namespace libMesh;
  
  class LaplaceQoI : public DifferentiableQoI
  {
  public:
    LaplaceQoI(){}
    virtual ~LaplaceQoI(){}
  
    virtual void init_qoi( std::vector<Number>& sys_qoi );
  
    virtual void postprocess( ){}
  
    virtual void element_qoi_derivative(DiffContext &context, const QoISet & qois);
  
    virtual void element_qoi (DiffContext &context, const QoISet & qois);
  
    virtual AutoPtr<DifferentiableQoI> clone( ) {
      return AutoPtr<DifferentiableQoI> ( new LaplaceQoI(*this) );
    }
  
  };
  #endif // L_QOI_H



The source file L-shaped.h without comments:

 
  #include "libmesh/enum_fe_family.h"
  #include "libmesh/fem_system.h"
  #include "libmesh/parameter_vector.h"
  #include "libmesh/qoi_set.h"
  #include "libmesh/system.h"
  
  using namespace libMesh;
  
  class LaplaceSystem : public FEMSystem
  {
  public:
    LaplaceSystem(EquationSystems& es,
                 const std::string& name_in,
                 const unsigned int number_in)
    : FEMSystem(es, name_in, number_in),
      _fe_family("LAGRANGE"), _fe_order(1),
      _analytic_jacobians(true) { }
  
    std::string & fe_family() { return _fe_family;  }
    unsigned int & fe_order() { return _fe_order;  }
    bool & analytic_jacobians() { return _analytic_jacobians; }
  
    Number &get_parameter_value(unsigned int parameter_index)
      {
        return parameters[parameter_index];
      }
  
    ParameterVector &get_parameter_vector()
      {
        parameter_vector.resize(parameters.size());
        for(unsigned int i = 0; i != parameters.size(); ++i)
  	{
  	  parameter_vector[i] = &parameters[i];
  	}
  
        return parameter_vector;
      }
  
    protected:
    virtual void init_data ();
  
    virtual void init_context (DiffContext &context);
  
    virtual bool element_time_derivative (bool request_jacobian,
  					DiffContext &context);
  
    virtual bool side_constraint (bool request_jacobian,
  				DiffContext &context);
  
    Number exact_solution (const Point&);
  
    std::vector<Number> parameters;
  
    ParameterVector parameter_vector;
  
    std::string _fe_family;
    unsigned int _fe_order;
  
    bool _analytic_jacobians;
  };



The source file adjoints_ex2.C without comments:

 
  
  
  #include <iostream>
  #include <iomanip>
  
  #include "libmesh/equation_systems.h"
  #include "libmesh/error_vector.h"
  #include "libmesh/mesh.h"
  #include "libmesh/mesh_refinement.h"
  #include "libmesh/newton_solver.h"
  #include "libmesh/numeric_vector.h"
  #include "libmesh/steady_solver.h"
  #include "libmesh/system_norm.h"
  
  #include "libmesh/parameter_vector.h"
  #include "libmesh/sensitivity_data.h"
  
  #include "libmesh/kelly_error_estimator.h"
  #include "libmesh/patch_recovery_error_estimator.h"
  
  #include "libmesh/adjoint_residual_error_estimator.h"
  #include "libmesh/qoi_set.h"
  
  #include "libmesh/getpot.h"
  #include "libmesh/gmv_io.h"
  
  #include "femparameters.h"
  #include "L-shaped.h"
  #include "L-qoi.h"
  
  using namespace libMesh;
  
  
  
  
  void write_output(EquationSystems &es,
  		  unsigned int a_step, // The adaptive step count
  		  std::string solution_type) // primal or adjoint solve
  {
    MeshBase &mesh = es.get_mesh();
  
  #ifdef LIBMESH_HAVE_GMV
  
    std::ostringstream file_name_gmv;
    file_name_gmv << solution_type
                  << ".out.gmv."
                  << std::setw(2)
                  << std::setfill('0')
                  << std::right
                  << a_step;
  
    GMVIO(mesh).write_equation_systems
      (file_name_gmv.str(), es);
  
  #endif
  }
  
  
  void set_system_parameters(LaplaceSystem &system, FEMParameters &param)
  {
    system.analytic_jacobians() = param.analytic_jacobians;
  
    system.verify_analytic_jacobians = param.verify_analytic_jacobians;
  
    system.fe_family() = param.fe_family[0];
    system.fe_order() = param.fe_order[0];
  
    system.print_solution_norms = param.print_solution_norms;
    system.print_solutions      = param.print_solutions;
    system.print_residual_norms = param.print_residual_norms;
    system.print_residuals      = param.print_residuals;
    system.print_jacobian_norms = param.print_jacobian_norms;
    system.print_jacobians      = param.print_jacobians;
  
    system.time_solver =
        AutoPtr<TimeSolver>(new SteadySolver(system));
  
    {
      NewtonSolver *solver = new NewtonSolver(system);
      system.time_solver->diff_solver() = AutoPtr<DiffSolver>(solver);
  
      solver->quiet                       = param.solver_quiet;
      solver->max_nonlinear_iterations    = param.max_nonlinear_iterations;
      solver->minsteplength               = param.min_step_length;
      solver->relative_step_tolerance     = param.relative_step_tolerance;
      solver->relative_residual_tolerance = param.relative_residual_tolerance;
      solver->require_residual_reduction  = param.require_residual_reduction;
      solver->linear_tolerance_multiplier = param.linear_tolerance_multiplier;
      if (system.time_solver->reduce_deltat_on_diffsolver_failure)
        {
  	solver->continue_after_max_iterations = true;
  	solver->continue_after_backtrack_failure = true;
        }
  
      solver->max_linear_iterations       = param.max_linear_iterations;
      solver->initial_linear_tolerance    = param.initial_linear_tolerance;
      solver->minimum_linear_tolerance    = param.minimum_linear_tolerance;
    }
  }
  
  
  #ifdef LIBMESH_ENABLE_AMR
  
  AutoPtr<MeshRefinement> build_mesh_refinement(MeshBase &mesh,
                                                FEMParameters &param)
  {
    AutoPtr<MeshRefinement> mesh_refinement(new MeshRefinement(mesh));
    mesh_refinement->coarsen_by_parents() = true;
    mesh_refinement->absolute_global_tolerance() = param.global_tolerance;
    mesh_refinement->nelem_target()      = param.nelem_target;
    mesh_refinement->refine_fraction()   = param.refine_fraction;
    mesh_refinement->coarsen_fraction()  = param.coarsen_fraction;
    mesh_refinement->coarsen_threshold() = param.coarsen_threshold;
  
    return mesh_refinement;
  }
  
  #endif // LIBMESH_ENABLE_AMR
  
  
  AutoPtr<ErrorEstimator> build_error_estimator(FEMParameters &param)
  {
    AutoPtr<ErrorEstimator> error_estimator;
  
    if (param.indicator_type == "kelly")
      {
        std::cout<<"Using Kelly Error Estimator"<<std::endl;
  
        error_estimator.reset(new KellyErrorEstimator);
      }
    else if (param.indicator_type == "adjoint_residual")
      {
        std::cout<<"Using Adjoint Residual Error Estimator with Patch Recovery Weights"<<std::endl<<std::endl;
  
        AdjointResidualErrorEstimator *adjoint_residual_estimator = new AdjointResidualErrorEstimator;
  
        error_estimator.reset (adjoint_residual_estimator);
  
        adjoint_residual_estimator->error_plot_suffix = "error.gmv";
  
        PatchRecoveryErrorEstimator *p1 =
  	new PatchRecoveryErrorEstimator;
        adjoint_residual_estimator->primal_error_estimator().reset(p1);
  
        PatchRecoveryErrorEstimator *p2 =
  	new PatchRecoveryErrorEstimator;
        adjoint_residual_estimator->dual_error_estimator().reset(p2);
  
        adjoint_residual_estimator->primal_error_estimator()->error_norm.set_type(0, H1_SEMINORM);
  
        adjoint_residual_estimator->dual_error_estimator()->error_norm.set_type(0, H1_SEMINORM);
      }
    else
      {
        std::cerr << "Unknown indicator_type" << std::endl;
        libmesh_error();
      }
    return error_estimator;
  }
  
  int main (int argc, char** argv)
  {
    LibMeshInit init (argc, argv);
  
  #ifndef LIBMESH_ENABLE_AMR
    libmesh_example_assert(false, "--enable-amr");
  #else
  
    std::cout << "Started " << argv[0] << std::endl;
  
    {
      std::ifstream i("general.in");
      if (!i)
        {
          std::cerr << '[' << libMesh::processor_id()
                    << "] Can't find general.in; exiting early."
                    << std::endl;
          libmesh_error();
        }
    }
    GetPot infile("general.in");
  
    FEMParameters param;
    param.read(infile);
  
    libmesh_example_assert(2 <= LIBMESH_DIM, "2D support");
  
    Mesh mesh(init.comm());
  
    AutoPtr<MeshRefinement> mesh_refinement =
      build_mesh_refinement(mesh, param);
  
    EquationSystems equation_systems (mesh);
  
    std::cout << "Reading in and building the mesh" << std::endl;
  
    mesh.read(param.domainfile.c_str());
    mesh.all_second_order();
  
    MeshRefinement initial_uniform_refinements(mesh);
    initial_uniform_refinements.uniformly_refine(param.coarserefinements);
  
    std::cout << "Building system" << std::endl;
  
    LaplaceSystem &system = equation_systems.add_system<LaplaceSystem> ("LaplaceSystem");
  
    QoISet qois;
  
    std::vector<unsigned int> qoi_indices;
    qoi_indices.push_back(0);
    qois.add_indices(qoi_indices);
  
    qois.set_weight(0, 0.5);
  
    {
      LaplaceQoI qoi;
      system.attach_qoi( &qoi );
    }
  
    set_system_parameters(system, param);
  
    std::cout << "Initializing systems" << std::endl;
  
    equation_systems.init ();
  
    mesh.print_info();
    equation_systems.print_info();
  
    {
      unsigned int a_step = 0;
      for (; a_step != param.max_adaptivesteps; ++a_step)
        {
  	if (param.global_tolerance != 0.)
  	  libmesh_assert_equal_to (param.nelem_target, 0);
              else
                libmesh_assert_greater (param.nelem_target, 0);
  
  	system.solve();
  
  	write_output(equation_systems, a_step, "primal");
  
  	NumericVector<Number> &primal_solution = *system.solution;
  
  	SensitivityData sensitivities(qois, system, system.get_parameter_vector());
  
  	system.assemble_qoi_sides = true;
  
  	system.set_adjoint_already_solved(false);
  
  	system.adjoint_qoi_parameter_sensitivity(qois, system.get_parameter_vector(), sensitivities);
  
  	system.set_adjoint_already_solved(true);
  
  	GetPot infile_l_shaped("l-shaped.in");
  
  	Number sensitivity_QoI_0_0_computed = sensitivities[0][0];
  	Number sensitivity_QoI_0_0_exact = infile_l_shaped("sensitivity_0_0", 0.0);
  	Number sensitivity_QoI_0_1_computed = sensitivities[0][1];
  	Number sensitivity_QoI_0_1_exact = infile_l_shaped("sensitivity_0_1", 0.0);
  
  	std::cout << "Adaptive step " << a_step << ", we have " << mesh.n_active_elem()
                        << " active elements and "
                        << equation_systems.n_active_dofs()
  		  << " active dofs." << std::endl ;
  
  	std::cout<<"Sensitivity of QoI one to Parameter one is "<<sensitivity_QoI_0_0_computed<<std::endl;
  	std::cout<<"Sensitivity of QoI one to Parameter two is "<<sensitivity_QoI_0_1_computed<<std::endl;
  
  	std::cout<< "The relative error in sensitivity QoI_0_0 is "
  		 << std::setprecision(17)
                   << std::abs(sensitivity_QoI_0_0_computed - sensitivity_QoI_0_0_exact) /
                      std::abs(sensitivity_QoI_0_0_exact) << std::endl;
  	std::cout<< "The relative error in sensitivity QoI_0_1 is "
                   << std::setprecision(17)
                   << std::abs(sensitivity_QoI_0_1_computed - sensitivity_QoI_0_1_exact) /
                      std::abs(sensitivity_QoI_0_1_exact) << std::endl << std::endl;
  
  	NumericVector<Number> &dual_solution_0 = system.get_adjoint_solution(0);
  
  	primal_solution.swap(dual_solution_0);
  	write_output(equation_systems, a_step, "adjoint_0");
  
  	primal_solution.swap(dual_solution_0);
  
  
  	if(param.refine_uniformly)
  	  {
  	    std::cout<<"Refining Uniformly"<<std::endl<<std::endl;
  
  	    mesh_refinement->uniformly_refine(1);
  	  }
  	else if(param.global_tolerance >= 0. && param.nelem_target == 0.)
  	  {
  	    ErrorVector error;
  
  	    AutoPtr<ErrorEstimator> error_estimator =
  	      build_error_estimator(param);
  
  	    error_estimator->estimate_error(system, error);
  
  	    mesh_refinement->flag_elements_by_error_tolerance (error);
  
  	    mesh_refinement->refine_and_coarsen_elements();
  	  }
  	else
  	  {
  	    ErrorVector error;
  
  	    AutoPtr<ErrorEstimator> error_estimator =
  	      build_error_estimator(param);
  
  	    error_estimator->estimate_error(system, error);
  
  	    if (mesh.n_active_elem() >= param.nelem_target)
  	      {
  		std::cout<<"We reached the target number of elements."<<std::endl <<std::endl;
  		break;
  	      }
  
  	    mesh_refinement->flag_elements_by_nelem_target (error);
  
  	    mesh_refinement->refine_and_coarsen_elements();
  	  }
  
  	equation_systems.reinit();
  
          std::cout << "Refined mesh to "
                    << mesh.n_active_elem()
                    << " active elements and "
                    << equation_systems.n_active_dofs()
                    << " active dofs." << std::endl;
        }
  
      if (a_step == param.max_adaptivesteps)
        {
  	system.solve();
  
  	write_output(equation_systems, a_step, "primal");
  
  	NumericVector<Number> &primal_solution = *system.solution;
  
  	SensitivityData sensitivities(qois, system, system.get_parameter_vector());
  
  	system.assemble_qoi_sides = true;
  
  	system.set_adjoint_already_solved(false);
  
  	system.adjoint_qoi_parameter_sensitivity(qois, system.get_parameter_vector(), sensitivities);
  
  	system.set_adjoint_already_solved(true);
  
  	GetPot infile_l_shaped("l-shaped.in");
  
  	Number sensitivity_QoI_0_0_computed = sensitivities[0][0];
  	Number sensitivity_QoI_0_0_exact = infile_l_shaped("sensitivity_0_0", 0.0);
  	Number sensitivity_QoI_0_1_computed = sensitivities[0][1];
  	Number sensitivity_QoI_0_1_exact = infile_l_shaped("sensitivity_0_1", 0.0);
  
  	std::cout << "Adaptive step " << a_step << ", we have " << mesh.n_active_elem()
  		  << " active elements and "
  		  << equation_systems.n_active_dofs()
  		  << " active dofs." << std::endl ;
  
  	std::cout<<"Sensitivity of QoI one to Parameter one is "<<sensitivity_QoI_0_0_computed<<std::endl;
  	std::cout<<"Sensitivity of QoI one to Parameter two is "<<sensitivity_QoI_0_1_computed<<std::endl;
  
  	std::cout<< "The error in sensitivity QoI_0_0 is "
                   << std::setprecision(17)
                   << std::abs(sensitivity_QoI_0_0_computed - sensitivity_QoI_0_0_exact)/sensitivity_QoI_0_0_exact << std::endl;
  	std::cout<< "The error in sensitivity QoI_0_1 is "
                   << std::setprecision(17)
                   << std::abs(sensitivity_QoI_0_1_computed - sensitivity_QoI_0_1_exact)/sensitivity_QoI_0_1_exact << std::endl << std::endl;
  
  	NumericVector<Number> &dual_solution_0 = system.get_adjoint_solution(0);
  
  	primal_solution.swap(dual_solution_0);
  	write_output(equation_systems, a_step, "adjoint_0");
  
  	primal_solution.swap(dual_solution_0);
        }
    }
  
    std::cerr << '[' << libMesh::processor_id()
              << "] Completing output." << std::endl;
  
  #endif // #ifndef LIBMESH_ENABLE_AMR
  
    return 0;
  }



The source file femparameters.C without comments:

 
  #include "femparameters.h"
  
  #define GETPOT_INPUT(A) { A = input(#A, A);\
    std::string stringval = input(#A, std::string());\
    variable_names.push_back(std::string(#A "=") + stringval); };
  #define GETPOT_INT_INPUT(A) { A = input(#A, (int)A);\
    std::string stringval = input(#A, std::string());\
    variable_names.push_back(std::string(#A "=") + stringval); };
  #define GETPOT_REGISTER(A) { \
    std::string stringval = input(#A, std::string());\
    variable_names.push_back(std::string(#A "=") + stringval); };
  
  
  void FEMParameters::read(GetPot &input)
  {
      std::vector<std::string> variable_names;
  
  
      GETPOT_INT_INPUT(coarserefinements);
      GETPOT_INPUT(domainfile);
  
      GETPOT_INPUT(solver_quiet);
      GETPOT_INPUT(require_residual_reduction);
      GETPOT_INPUT(min_step_length);
      GETPOT_INT_INPUT(max_linear_iterations);
      GETPOT_INT_INPUT(max_nonlinear_iterations);
      GETPOT_INPUT(relative_step_tolerance);
      GETPOT_INPUT(relative_residual_tolerance);
      GETPOT_INPUT(initial_linear_tolerance);
      GETPOT_INPUT(minimum_linear_tolerance);
      GETPOT_INPUT(linear_tolerance_multiplier);
      GETPOT_INT_INPUT(nelem_target);
      GETPOT_INPUT(global_tolerance);
      GETPOT_INPUT(refine_fraction);
      GETPOT_INPUT(coarsen_fraction);
      GETPOT_INPUT(coarsen_threshold);
      GETPOT_INT_INPUT(max_adaptivesteps);
      GETPOT_INPUT(refine_uniformly);
      GETPOT_INPUT(indicator_type);
  
      GETPOT_REGISTER(fe_family);
      const unsigned int n_fe_family =
        std::max(1u, input.vector_variable_size("fe_family"));
      fe_family.resize(n_fe_family, "LAGRANGE");
      for (unsigned int i=0; i != n_fe_family; ++i)
        fe_family[i]              = input("fe_family", fe_family[i].c_str(), i);
      GETPOT_REGISTER(fe_order);
      const unsigned int n_fe_order =
        input.vector_variable_size("fe_order");
      fe_order.resize(n_fe_order, 1);
      for (unsigned int i=0; i != n_fe_order; ++i)
        fe_order[i]               = input("fe_order", (int)fe_order[i], i);
  
      GETPOT_INPUT(analytic_jacobians);
      GETPOT_INPUT(verify_analytic_jacobians);
      GETPOT_INPUT(print_solution_norms);
      GETPOT_INPUT(print_solutions);
      GETPOT_INPUT(print_residual_norms);
      GETPOT_INPUT(print_residuals);
      GETPOT_INPUT(print_jacobian_norms);
      GETPOT_INPUT(print_jacobians);
  
    std::vector<std::string> bad_variables =
      input.unidentified_arguments(variable_names);
  
    if (libMesh::processor_id() == 0 && !bad_variables.empty())
      {
        std::cerr << "ERROR: Unrecognized variables:" << std::endl;
        for (unsigned int i = 0; i != bad_variables.size(); ++i)
          std::cerr << bad_variables[i] << std::endl;
        std::cerr << "not found among recognized variables." << std::endl;
        for (unsigned int i = 0; i != variable_names.size(); ++i)
          std::cerr << variable_names[i] << std::endl;
        libmesh_error();
      }
   }



The source file L-qoi.C without comments:

 
  #include "L-qoi.h"
  
  using namespace libMesh;
  
  void LaplaceQoI::init_qoi( std::vector<Number>& sys_qoi )
  {
    sys_qoi.resize(1);
    return;
  }
  
  void LaplaceQoI::element_qoi (DiffContext &context,
  			      const QoISet & /* qois */ )
  
  {
    FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
  
    const std::vector<Real> &JxW = c.element_fe_var[0]->get_JxW();
  
    const std::vector<Point> &xyz = c.element_fe_var[0]->get_xyz();
  
    unsigned int n_qpoints = (c.get_element_qrule())->n_points();
  
    Number dQoI_0 = 0.;
  
  
    for (unsigned int qp = 0; qp != n_qpoints; qp++)
      {
        const Real xf = xyz[qp](0);
        const Real yf = xyz[qp](1);
  
        if(fabs(xf - 0.875) <= 0.125 && fabs(yf - 0.125) <= 0.125)
        	{
        	  Number T = c.interior_value(0, qp);
  
        	  dQoI_0 += JxW[qp] * T;
        	}
  
      }
  
  
    c.elem_qoi[0] = c.elem_qoi[0] + dQoI_0;
  
  }
  
  void LaplaceQoI::element_qoi_derivative (DiffContext &context,
  					 const QoISet & /* qois */)
  {
    FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
  
  
    const std::vector<Real> &JxW = c.element_fe_var[0]->get_JxW();
  
    const std::vector<std::vector<Real> >          &phi = c.element_fe_var[0]->get_phi();
  
    const std::vector<Point > &q_point = c.element_fe_var[0]->get_xyz();
  
    const unsigned int n_T_dofs = c.dof_indices_var[0].size();
    unsigned int n_qpoints = (c.get_element_qrule())->n_points();
  
    DenseSubVector<Number> &Q = *c.elem_qoi_subderivatives[0][0];
  
    for (unsigned int qp=0; qp != n_qpoints; qp++)
      {
        const Real x = q_point[qp](0);
        const Real y = q_point[qp](1);
  
        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
  }



The source file L-shaped.C without comments:

 
  #include "libmesh/getpot.h"
  #include "libmesh/fe_base.h"
  #include "libmesh/quadrature.h"
  #include "libmesh/string_to_enum.h"
  #include "libmesh/parallel.h"
  #include "libmesh/fem_context.h"
  
  #include "L-shaped.h"
  
  using namespace libMesh;
  
  void LaplaceSystem::init_data ()
  {
    this->add_variable ("T", static_cast<Order>(_fe_order),
                        Utility::string_to_enum<FEFamily>(_fe_family));
  
    GetPot infile("l-shaped.in");
    parameters.push_back(infile("alpha_1", 1.0));
    parameters.push_back(infile("alpha_2", 1.0));
  
    FEMSystem::init_data();
  
    this->time_evolving(0);
  }
  
  void LaplaceSystem::init_context(DiffContext &context)
  {
    FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
  
    c.element_fe_var[0]->get_JxW();
    c.element_fe_var[0]->get_phi();
    c.element_fe_var[0]->get_dphi();
  
    c.side_fe_var[0]->get_JxW();
    c.side_fe_var[0]->get_phi();
    c.side_fe_var[0]->get_dphi();
  }
  
  #define optassert(X) {if (!(X)) libmesh_error();}
  
  bool LaplaceSystem::element_time_derivative (bool request_jacobian,
  					  DiffContext &context)
  {
    bool compute_jacobian = request_jacobian && _analytic_jacobians;
  
    FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
  
  
    const std::vector<Real> &JxW = c.element_fe_var[0]->get_JxW();
  
    const std::vector<std::vector<RealGradient> > &dphi = c.element_fe_var[0]->get_dphi();
  
    const unsigned int n_T_dofs = c.dof_indices_var[0].size();
  
    DenseSubMatrix<Number> &K = *c.elem_subjacobians[0][0];
    DenseSubVector<Number> &F = *c.elem_subresiduals[0];
  
    unsigned int n_qpoints = (c.get_element_qrule())->n_points();
  
    for (unsigned int qp=0; qp != n_qpoints; qp++)
      {
        Gradient grad_T = c.interior_gradient(0, qp);
  
        for (unsigned int i=0; i != n_T_dofs; i++)
          F(i) += JxW[qp] * (parameters[0] + (2.*parameters[1])) * ( grad_T * dphi[i][qp] ) ;
        if (compute_jacobian)
          for (unsigned int i=0; i != n_T_dofs; i++)
            for (unsigned int j=0; j != n_T_dofs; ++j)
              K(i,j) += JxW[qp] * (parameters[0] + (2.*parameters[1])) * ( dphi[i][qp] * dphi[j][qp] );
      } // end of the quadrature point qp-loop
  
    return compute_jacobian;
  }
  
  bool LaplaceSystem::side_constraint (bool request_jacobian,
  				  DiffContext &context)
  {
    bool compute_jacobian = request_jacobian && _analytic_jacobians;
  
    FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
  
  
    const std::vector<Real> &JxW = c.side_fe_var[0]->get_JxW();
  
    const std::vector<std::vector<Real> > &phi = c.side_fe_var[0]->get_phi();
  
    const std::vector<Point > &qside_point = c.side_fe_var[0]->get_xyz();
  
    const unsigned int n_T_dofs = c.dof_indices_var[0].size();
  
    DenseSubMatrix<Number> &K = *c.elem_subjacobians[0][0];
    DenseSubVector<Number> &F = *c.elem_subresiduals[0];
  
    unsigned int n_qpoints = (c.get_side_qrule())->n_points();
  
    const Real penalty = 1./(TOLERANCE*TOLERANCE);
  
    for (unsigned int qp=0; qp != n_qpoints; qp++)
      {
        Number T = c.side_value(0, qp);
  
        Number u_dirichlet = exact_solution (qside_point[qp]);
  
        for (unsigned int i=0; i != n_T_dofs; i++)
  	F(i) += JxW[qp] * penalty * ( T - u_dirichlet) * phi[i][qp];
        if (compute_jacobian)
  	for (unsigned int i=0; i != n_T_dofs; i++)
  	  for (unsigned int j=0; j != n_T_dofs; ++j)
  	    K(i,j) += JxW[qp] * penalty * phi[i][qp] * phi[j][qp];
  
      } // end of the quadrature point qp-loop
  
    return compute_jacobian;
  }
  
  Number LaplaceSystem::exact_solution(const Point& p)// xyz location
  {
    const Real x1 = p(0);
    const Real x2 = p(1);
  
    Real theta = atan2(x2,x1);
  
    if (theta < 0)
      theta += 2. * libMesh::pi;
  
    return (parameters[0] + (2.*parameters[1])) * pow(x1*x1 + x2*x2, 1./3.)*sin(2./3.*theta);
  }



The console output of the program:

make[4]: Entering directory `/net/spark/workspace/roystgnr/libmesh/git/devel/examples/adjoints/adjoints_ex2'
***************************************************************
* Running Example adjoints_ex2:
*  mpirun -np 4 example-devel  -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc
***************************************************************
 
Started /net/spark/workspace/roystgnr/libmesh/git/devel/examples/adjoints/adjoints_ex2/.libs/lt-example-devel
Reading in and building the mesh
Building system
Initializing systems
 Mesh Information:
  mesh_dimension()=2
  spatial_dimension()=3
  n_nodes()=225
    n_local_nodes()=63
  n_elem()=63
    n_local_elem()=18
    n_active_elem()=48
  n_subdomains()=1
  n_partitions()=4
  n_processors()=4
  n_threads()=1
  processor_id()=0

 EquationSystems
  n_systems()=1
   System #0, "LaplaceSystem"
    Type "Implicit"
    Variables="T" 
    Finite Element Types="LAGRANGE", "JACOBI_20_00" 
    Infinite Element Mapping="CARTESIAN" 
    Approximation Orders="SECOND", "THIRD" 
    n_dofs()=225
    n_local_dofs()=63
    n_constrained_dofs()=0
    n_local_constrained_dofs()=0
    n_vectors()=1
    n_matrices()=1
    DofMap Sparsity
      Average  On-Processor Bandwidth <= 12.64
      Average Off-Processor Bandwidth <= 2.09778
      Maximum  On-Processor Bandwidth <= 25
      Maximum Off-Processor Bandwidth <= 16
    DofMap Constraints
      Number of DoF Constraints = 0
      Number of Node Constraints = 0

  Nonlinear solver converged, step 0, residual reduction 2.23191e-14 < 1e-09
Adaptive step 0, we have 48 active elements and 225 active dofs.
Sensitivity of QoI one to Parameter one is 0.00543563
Sensitivity of QoI one to Parameter two is 0.0108713
The relative error in sensitivity QoI_0_0 is 0.00027993001869348769
The relative error in sensitivity QoI_0_1 is 0.0002799300307674672

Refining Uniformly

Refined mesh to 192 active elements and 833 active dofs.
  Nonlinear solver converged, step 0, residual reduction 3.9022539966880865e-12 < 1.0000000000000001e-09
Adaptive step 1, we have 192 active elements and 833 active dofs.
Sensitivity of QoI one to Parameter one is 0.0054364699404312834
Sensitivity of QoI one to Parameter two is 0.010872939878913211
The relative error in sensitivity QoI_0_0 is 0.00012462933890172227
The relative error in sensitivity QoI_0_1 is 0.00012462951816411252

Refining Uniformly

Refined mesh to 768 active elements and 3201 active dofs.
  Nonlinear solver converged, step 0, residual reduction 2.298909708194788e-11 < 1.0000000000000001e-09
Adaptive step 2, we have 768 active elements and 3201 active dofs.
Sensitivity of QoI one to Parameter one is 0.005436875146067748
Sensitivity of QoI one to Parameter two is 0.010873750291287549
The relative error in sensitivity QoI_0_0 is 5.0103931717751492e-05
The relative error in sensitivity QoI_0_1 is 5.010400969459807e-05

Refining Uniformly

Refined mesh to 3072 active elements and 12545 active dofs.
  Nonlinear solver converged, step 0, residual reduction 2.1301677837111685e-10 < 1.0000000000000001e-09
Adaptive step 3, we have 3072 active elements and 12545 active dofs.
Sensitivity of QoI one to Parameter one is 0.0054370391087449682
Sensitivity of QoI one to Parameter two is 0.010874078216664211
The error in sensitivity QoI_0_0 is 1.9947921575799245e-05
The error in sensitivity QoI_0_1 is 1.9947997509128744e-05


 -------------------------------------------------------------------------------------------------------------------
| Processor id:   0                                                                                                 |
| Num Processors: 4                                                                                                 |
| Time:           Fri Apr 19 11:47:13 2013                                                                          |
| OS:             Linux                                                                                             |
| HostName:       spark.ices.utexas.edu                                                                             |
| OS Release:     2.6.32-279.22.1.el6.x86_64                                                                        |
| OS Version:     #1 SMP Tue Feb 5 14:33:39 CST 2013                                                                |
| Machine:        x86_64                                                                                            |
| Username:       roystgnr                                                                                          |
| Configuration:  ../configure  '--enable-everything'                                                               |
|  'METHODS=devel'                                                                                                  |
|  '--prefix=/h2/roystgnr/libmesh-test'                                                                             |
|  'CXX=distcc /usr/bin/g++'                                                                                        |
|  'CC=distcc /usr/bin/gcc'                                                                                         |
|  'FC=distcc /usr/bin/gfortran'                                                                                    |
|  'F77=distcc /usr/bin/gfortran'                                                                                   |
|  'PETSC_DIR=/opt/apps/ossw/libraries/petsc/petsc-3.3-p2'                                                          |
|  'PETSC_ARCH=gcc-system-mkl-gf-10.3.12.361-mpich2-1.4.1p1-cxx-opt'                                                |
|  'SLEPC_DIR=/opt/apps/ossw/libraries/slepc/slepc-3.3-p2-petsc-3.3-p2-cxx-opt'                                     |
|  'TRILINOS_DIR=/opt/apps/ossw/libraries/trilinos/trilinos-10.12.2/sl6/gcc-system/mpich2-1.4.1p1/mkl-gf-10.3.12.361'|
|  'VTK_DIR=/opt/apps/ossw/libraries/vtk/vtk-5.10.0/sl6/gcc-system'                                                 |
|  'HDF5_DIR=/opt/apps/ossw/libraries/hdf5/hdf5-1.8.9/sl6/gcc-system'                                               |
 -------------------------------------------------------------------------------------------------------------------
 ----------------------------------------------------------------------------------------------------------------
| libMesh Performance: Alive time=3.09365, Active time=2.57486                                                   |
 ----------------------------------------------------------------------------------------------------------------
| Event                              nCalls    Total Time  Avg Time    Total Time  Avg Time    % of Active Time  |
|                                              w/o Sub     w/o Sub     With Sub    With Sub    w/o S    With S   |
|----------------------------------------------------------------------------------------------------------------|
|                                                                                                                |
|                                                                                                                |
| DofMap                                                                                                         |
|   add_neighbors_to_send_list()     4         0.0058      0.001459    0.0071      0.001776    0.23     0.28     |
|   build_sparsity()                 4         0.0043      0.001067    0.0189      0.004728    0.17     0.73     |
|   create_dof_constraints()         4         0.0042      0.001043    0.0042      0.001043    0.16     0.16     |
|   distribute_dofs()                4         0.0172      0.004304    0.0757      0.018931    0.67     2.94     |
|   dof_indices()                    33866     0.1913      0.000006    0.1913      0.000006    7.43     7.43     |
|   old_dof_indices()                4032      0.0220      0.000005    0.0220      0.000005    0.86     0.86     |
|   prepare_send_list()              4         0.0001      0.000014    0.0001      0.000014    0.00     0.00     |
|   reinit()                         4         0.0228      0.005689    0.0228      0.005689    0.88     0.88     |
|                                                                                                                |
| EquationSystems                                                                                                |
|   build_solution_vector()          8         0.0052      0.000654    0.0334      0.004174    0.20     1.30     |
|                                                                                                                |
| FE                                                                                                             |
|   compute_shape_functions()        30576     0.3762      0.000012    0.3762      0.000012    14.61    14.61    |
|   init_shape_functions()           2128      0.0154      0.000007    0.0154      0.000007    0.60     0.60     |
|   inverse_map()                    14570     0.0510      0.000003    0.0510      0.000003    1.98     1.98     |
|                                                                                                                |
| FEMSystem                                                                                                      |
|   assemble_qoi()                   20        0.0324      0.001620    0.5118      0.025590    1.26     19.88    |
|   assemble_qoi_derivative()        4         0.0747      0.018670    0.1248      0.031211    2.90     4.85     |
|   assembly()                       8         0.0574      0.007175    0.1593      0.019915    2.23     6.19     |
|   assembly(get_jacobian)           4         0.0280      0.006994    0.0787      0.019666    1.09     3.06     |
|   assembly(get_residual)           20        0.0671      0.003355    0.3191      0.015953    2.61     12.39    |
|                                                                                                                |
| FEMap                                                                                                          |
|   compute_affine_map()             30576     0.0868      0.000003    0.0868      0.000003    3.37     3.37     |
|   compute_face_map()               2016      0.0133      0.000007    0.0357      0.000018    0.52     1.39     |
|   init_face_shape_functions()      64        0.0002      0.000003    0.0002      0.000003    0.01     0.01     |
|   init_reference_to_physical_map() 2128      0.0173      0.000008    0.0173      0.000008    0.67     0.67     |
|                                                                                                                |
| GMVIO                                                                                                          |
|   write_nodal_data()               8         0.1193      0.014910    0.1199      0.014983    4.63     4.66     |
|                                                                                                                |
| ImplicitSystem                                                                                                 |
|   adjoint_solve()                  4         0.0009      0.000217    0.4419      0.110475    0.03     17.16    |
|                                                                                                                |
| LocationMap                                                                                                    |
|   find()                           20460     0.0235      0.000001    0.0235      0.000001    0.91     0.91     |
|   init()                           8         0.0063      0.000786    0.0063      0.000786    0.24     0.24     |
|                                                                                                                |
| Mesh                                                                                                           |
|   all_second_order()               1         0.0001      0.000090    0.0001      0.000090    0.00     0.00     |
|   contract()                       3         0.0018      0.000586    0.0057      0.001886    0.07     0.22     |
|   find_neighbors()                 6         0.0153      0.002550    0.0237      0.003951    0.59     0.92     |
|   renumber_nodes_and_elem()        3         0.0039      0.001299    0.0039      0.001299    0.15     0.15     |
|                                                                                                                |
| MeshCommunication                                                                                              |
|   broadcast()                      1         0.0002      0.000236    0.0005      0.000525    0.01     0.02     |
|   compute_hilbert_indices()        7         0.0096      0.001378    0.0096      0.001378    0.37     0.37     |
|   find_global_indices()            7         0.0020      0.000281    0.0221      0.003158    0.08     0.86     |
|   parallel_sort()                  7         0.0008      0.000114    0.0093      0.001333    0.03     0.36     |
|                                                                                                                |
| MeshOutput                                                                                                     |
|   write_equation_systems()         8         0.0002      0.000024    0.1543      0.019285    0.01     5.99     |
|                                                                                                                |
| MeshRefinement                                                                                                 |
|   _coarsen_elements()              3         0.0004      0.000146    0.0008      0.000261    0.02     0.03     |
|   _refine_elements()               8         0.0269      0.003363    0.1265      0.015806    1.04     4.91     |
|   add_point()                      20460     0.0272      0.000001    0.0525      0.000003    1.06     2.04     |
|   make_coarsening_compatible()     6         0.0096      0.001602    0.0096      0.001602    0.37     0.37     |
|   make_flags_parallel_consistent() 6         0.0045      0.000754    0.0096      0.001593    0.18     0.37     |
|   make_refinement_compatible()     6         0.0000      0.000002    0.0000      0.000006    0.00     0.00     |
|                                                                                                                |
| MetisPartitioner                                                                                               |
|   partition()                      6         0.0320      0.005337    0.0564      0.009404    1.24     2.19     |
|                                                                                                                |
| NewtonSolver                                                                                                   |
|   solve()                          4         0.1303      0.032566    0.5814      0.145342    5.06     22.58    |
|                                                                                                                |
| Parallel                                                                                                       |
|   allgather()                      28        0.0239      0.000852    0.0240      0.000857    0.93     0.93     |
|   broadcast()                      27        0.0002      0.000006    0.0001      0.000005    0.01     0.00     |
|   max(bool)                        19        0.0468      0.002461    0.0468      0.002461    1.82     1.82     |
|   max(scalar)                      937       0.0048      0.000005    0.0048      0.000005    0.19     0.19     |
|   max(vector)                      219       0.0014      0.000006    0.0044      0.000020    0.05     0.17     |
|   min(bool)                        1152      0.0120      0.000010    0.0120      0.000010    0.47     0.47     |
|   min(scalar)                      922       0.3139      0.000340    0.3139      0.000340    12.19    12.19    |
|   min(vector)                      219       0.0016      0.000007    0.0061      0.000028    0.06     0.24     |
|   probe()                          222       0.0122      0.000055    0.0122      0.000055    0.47     0.47     |
|   receive()                        222       0.0006      0.000003    0.0129      0.000058    0.02     0.50     |
|   send()                           222       0.0004      0.000002    0.0004      0.000002    0.02     0.02     |
|   send_receive()                   236       0.0008      0.000003    0.0144      0.000061    0.03     0.56     |
|   sum()                            105       0.0108      0.000103    0.2567      0.002445    0.42     9.97     |
|                                                                                                                |
| Parallel::Request                                                                                              |
|   wait()                           222       0.0002      0.000001    0.0002      0.000001    0.01     0.01     |
|                                                                                                                |
| Partitioner                                                                                                    |
|   set_node_processor_ids()         6         0.0074      0.001229    0.0349      0.005810    0.29     1.35     |
|   set_parent_processor_ids()       6         0.0015      0.000250    0.0015      0.000250    0.06     0.06     |
|                                                                                                                |
| PetscLinearSolver                                                                                              |
|   solve()                          8         0.5427      0.067839    0.5427      0.067839    21.08    21.08    |
|                                                                                                                |
| ProjectVector                                                                                                  |
|   operator()                       6         0.0168      0.002795    0.0693      0.011546    0.65     2.69     |
|                                                                                                                |
| System                                                                                                         |
|   project_vector()                 6         0.0695      0.011591    0.1508      0.025136    2.70     5.86     |
 ----------------------------------------------------------------------------------------------------------------
| Totals:                            165850    2.5749                                          100.00            |
 ----------------------------------------------------------------------------------------------------------------

 
***************************************************************
* Done Running Example adjoints_ex2:
*  mpirun -np 4 example-devel  -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc
***************************************************************
make[4]: Leaving directory `/net/spark/workspace/roystgnr/libmesh/git/devel/examples/adjoints/adjoints_ex2'

Site Created By: libMesh Developers
Last modified: April 23 2013 04:19:30 UTC

Hosted By:
SourceForge.net Logo