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),
        	reuse_preconditioner(false),
              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"), patch_reuse(true),
              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);
        
            std::string domainfile;
            unsigned int coarserefinements;
        
            bool solver_quiet, reuse_preconditioner, 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;
            bool patch_reuse;
        
            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-shaped.h with comments:

        #include "libmesh/enum_fe_family.h"
        #include "libmesh/fem_system.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) { qoi.resize(2); }
        
          std::string & fe_family() { return _fe_family;  }
          unsigned int & fe_order() { return _fe_order;  }
          bool & analytic_jacobians() { return _analytic_jacobians; }
        
Postprocessing function which we are going to override for this application

          virtual void postprocess(void);
        
          Number &get_QoI_value(std::string type, unsigned int QoI_index)
            {
              if(type == "exact")
        	{
        	  return exact_QoI[QoI_index];
        	}
              else
        	{
        	  return computed_QoI[QoI_index];
        	}
            }
        
          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);
        
Overloading the postprocess function

          virtual void element_postprocess(DiffContext &context);
        
          virtual void side_postprocess(DiffContext &context);
        
Overloading the qoi function on elements

          virtual void element_qoi_derivative
            (DiffContext &context,
             const QoISet & qois);
        
Overloading the qoi function on sides

          virtual void side_qoi_derivative
            (DiffContext &context,
             const QoISet & qois);
        
          Number exact_solution (const Point&);
        
Variables to hold the computed QoIs

          Number computed_QoI[2];
        
Variables to read in the exact QoIs from l-shaped.in

          Number exact_QoI[2];
        
The FE type to use
          std::string _fe_family;
          unsigned int _fe_order;
        
Calculate Jacobians analytically or not?
          bool _analytic_jacobians;
        };



The source file adjoints_ex4.C with comments:

Adjoints Example 1 - Laplace Equation in the L-Shaped Domain with Adjoint based mesh refinement



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). The kelly and adjoint residual error estimators are used to develop error indicators and guide mesh adaptation. Since we use the adjoint capabilities of libMesh in this example, we use the DiffSystem framework. This file (adjoints_ex1.C) contains the declaration of mesh and equation system objects, L-shaped.C contains the assembly of the system, element_qoi_derivative.C and side_qoi_derivative.C contain the RHS for the adjoint systems. Postprocessing to compute the QoIs is done in element_postprocess.C and side_postprocess.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/linear_solver.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"
        
Adjoint Related includes
        #include "libmesh/qoi_set.h"
        #include "libmesh/adjoint_refinement_estimator.h"
        
libMesh I/O includes
        #include "libmesh/getpot.h"
        #include "libmesh/gmv_io.h"
        
Local includes
        #include "femparameters.h"
        #include "L-shaped.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;
        }
        
This is where declare the adjoint refined error estimator. This estimator builds an error bound for Q(u) - Q(u_h), by solving the adjoint problem on a finer Finite Element space. For more details see the description of the Adjoint Refinement Error Estimator in adjoint_refinement_error_estimator.C
        AutoPtr<AdjointRefinementEstimator> build_adjoint_refinement_error_estimator(QoISet &qois)
        {
          AutoPtr<AdjointRefinementEstimator> error_estimator;
        
          std::cout<<"Computing the error estimate using the Adjoint Refinement Error Estimator"<<std::endl<<std::endl;
        
          AdjointRefinementEstimator *adjoint_refinement_estimator = new AdjointRefinementEstimator;
        
          error_estimator.reset (adjoint_refinement_estimator);
        
          adjoint_refinement_estimator->qoi_set() = qois;
        
We enrich the FE space for the dual problem by doing 2 uniform h refinements
          adjoint_refinement_estimator->number_h_refinements = 2;
        
          return error_estimator;
        }
        
        #endif // LIBMESH_ENABLE_AMR
        
        
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");
        
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();
          LinearSolver<Number> *linear_solver = system.get_linear_solver();
        
          {
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);
        
        	linear_solver->reuse_preconditioner(false);
        
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;
        
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<unsigned int> 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);
        
Make sure we get the contributions to the adjoint RHS from the sides
                system.assemble_qoi_sides = true;
        
We are about to solve the adjoint system, but before we do this we see the same preconditioner flag to reuse the preconditioner from the forward solver
                linear_solver->reuse_preconditioner(param.reuse_preconditioner);
        
Solve the adjoint system. This takes the transpose of the stiffness matrix and then solves the resulting system
                system.adjoint_solve();
        
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);
        
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);
        
Get a pointer to the solution vector of the adjoint problem for QoI 0
                NumericVector<Number> &dual_solution_1 = system.get_adjoint_solution(1);
        
Swap again
                primal_solution.swap(dual_solution_1);
        	write_output(equation_systems, a_step, "adjoint_1");
        
Swap back again
                primal_solution.swap(dual_solution_1);
        
        	std::cout << "Adaptive step " << a_step << ", we have " << mesh.n_active_elem()
                              << " active elements and "
                              << equation_systems.n_active_dofs()
        		  << " active dofs." << std::endl ;
        
Postprocess, compute the approximate QoIs and write them out to the console
                std::cout << "Postprocessing: " << std::endl;
        	system.postprocess_sides = true;
        	system.postprocess();
        	Number QoI_0_computed = system.get_QoI_value("computed", 0);
        	Number QoI_0_exact = system.get_QoI_value("exact", 0);
        	Number QoI_1_computed = system.get_QoI_value("computed", 1);
        	Number QoI_1_exact = system.get_QoI_value("exact", 1);
        
        	std::cout<< "The relative error in QoI 0 is " << std::setprecision(17)
                         << std::abs(QoI_0_computed - QoI_0_exact) /
                            std::abs(QoI_0_exact) << std::endl;
        
        	std::cout<< "The relative error in QoI 1 is " << std::setprecision(17)
                         << std::abs(QoI_1_computed - QoI_1_exact) /
                            std::abs(QoI_1_exact) << std::endl << std::endl;
        
We will declare an error vector for passing to the adjoint refinement error estimator
                ErrorVector QoI_elementwise_error;
        
Build an adjoint refinement error estimator object
                AutoPtr<AdjointRefinementEstimator> adjoint_refinement_error_estimator =
        	  build_adjoint_refinement_error_estimator(qois);
        
Estimate the error in each element using the Adjoint Refinement estimator
                adjoint_refinement_error_estimator->estimate_error(system, QoI_elementwise_error);
        
Print out the computed error estimate, note that we access the global error estimates using an accessor function, right now sum(QoI_elementwise_error) != global_QoI_error_estimate
                std::cout<< "The computed relative error in QoI 0 is " << std::setprecision(17)
                         << std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(0)) /
                            std::abs(QoI_0_exact) << std::endl;
        
        	std::cout<< "The computed relative error in QoI 1 is " << std::setprecision(17)
                         << std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(1)) /
        	  std::abs(QoI_1_exact) << std::endl << std::endl;
        
Also print out effecitivity indices (estimated error/true error)
                std::cout<< "The effectivity index for the computed error in QoI 0 is " << std::setprecision(17)
                         << std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(0)) /
                            std::abs(QoI_0_computed - QoI_0_exact) << std::endl;
        
        	std::cout<< "The effectivity index for the computed error in QoI 1 is " << std::setprecision(17)
                         << std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(1)) /
        	  std::abs(QoI_1_computed - QoI_1_exact) << std::endl << std::endl;
        
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)
        	  {
        	    mesh_refinement->uniformly_refine(1);
        	  }
Adaptively refine based on reaching an error tolerance
                else if(param.global_tolerance >= 0. && param.nelem_target == 0.)
        	  {
        	    mesh_refinement->flag_elements_by_error_tolerance (QoI_elementwise_error);
        
        	    mesh_refinement->refine_and_coarsen_elements();
        	  }
Adaptively refine based on reaching a target number of elements
                else
        	  {
        	    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 (QoI_elementwise_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)
              {
        	linear_solver->reuse_preconditioner(false);
        	system.solve();
        
        	write_output(equation_systems, a_step, "primal");
        
        	NumericVector<Number> &primal_solution = *system.solution;
        
        	QoISet qois;
        	std::vector<unsigned int> 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);
        
        	system.assemble_qoi_sides = true;
        	linear_solver->reuse_preconditioner(param.reuse_preconditioner);
        	system.adjoint_solve();
        
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);
        
        	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);
        
        	NumericVector<Number> &dual_solution_1 = system.get_adjoint_solution(1);
        
        	primal_solution.swap(dual_solution_1);
        	write_output(equation_systems, a_step, "adjoint_1");
        
        	primal_solution.swap(dual_solution_1);
        
        	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 << "Postprocessing: " << std::endl;
        	system.postprocess_sides = true;
        	system.postprocess();
        
        	Number QoI_0_computed = system.get_QoI_value("computed", 0);
        	Number QoI_0_exact = system.get_QoI_value("exact", 0);
        	Number QoI_1_computed = system.get_QoI_value("computed", 1);
        	Number QoI_1_exact = system.get_QoI_value("exact", 1);
        
        	std::cout<< "The relative error in QoI 0 is " << std::setprecision(17)
                         << std::abs(QoI_0_computed - QoI_0_exact) /
                            std::abs(QoI_0_exact) << std::endl;
        
        	std::cout<< "The relative error in QoI 1 is " << std::setprecision(17)
                         << std::abs(QoI_1_computed - QoI_1_exact) /
                            std::abs(QoI_1_exact) << std::endl << std::endl;
        
We will declare an error vector for passing to the adjoint refinement error estimator Right now, only the first entry of this vector will be filled (with the global QoI error estimate) Later, each entry of the vector will contain elementwise error that the user can sum to get the total error
                ErrorVector QoI_elementwise_error;
        
Build an adjoint refinement error estimator object
                AutoPtr<AdjointRefinementEstimator> adjoint_refinement_error_estimator =
        	  build_adjoint_refinement_error_estimator(qois);
        
Estimate the error in each element using the Adjoint Refinement estimator
                adjoint_refinement_error_estimator->estimate_error(system, QoI_elementwise_error);
        
Print out the computed error estimate, note that we access the global error estimates using an accessor function, right now sum(QoI_elementwise_error) != global_QoI_error_estimate
                std::cout<< "The computed relative error in QoI 0 is " << std::setprecision(17)
                         << std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(0)) /
                            std::abs(QoI_0_exact) << std::endl;
        
        	std::cout<< "The computed relative error in QoI 1 is " << std::setprecision(17)
                         << std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(1)) /
        	  std::abs(QoI_1_exact) << std::endl << std::endl;
        
Also print out effecitivity indices (estimated error/true error)
                std::cout<< "The effectivity index for the computed error in QoI 0 is " << std::setprecision(17)
                         << std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(0)) /
                            std::abs(QoI_0_computed - QoI_0_exact) << std::endl;
        
        	std::cout<< "The effectivity index for the computed error in QoI 1 is " << std::setprecision(17)
                         << std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(1)) /
        	  std::abs(QoI_1_computed - QoI_1_exact) << std::endl << std::endl;
        
              }
          }
        
          std::cerr << '[' << libMesh::processor_id()
                    << "] Completing output." << std::endl;
        
        #endif // #ifndef LIBMESH_ENABLE_AMR
        
All done.
          return 0;
        }



The source file element_postprocess.C with comments:

        #include "libmesh/libmesh_common.h"
        #include "libmesh/elem.h"
        #include "libmesh/fe_base.h"
        #include "libmesh/fem_context.h"
        #include "libmesh/point.h"
        #include "libmesh/quadrature.h"
        
Local includes
        #include "L-shaped.h"
        
Bring in everything from the libMesh namespace
        using namespace libMesh;
        
Last modified: Feb 13, 2009

Define the postprocess function to compute QoI 0, the integral of the the solution over a subdomain

        void LaplaceSystem::element_postprocess (DiffContext &context)
        
        {
          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();
        
The number of local degrees of freedom in each variable

          unsigned int n_qpoints = (c.get_element_qrule())->n_points();
        
The function R = int_{omega} T dR omega is a subset of Omega (the whole domain), omega = [0.75, 1.0] x [0.0, 0.25]

          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 x = xyz[qp](0);
              const Real y = xyz[qp](1);
        
If in the sub-domain omega, add the contribution to the integral R
              if(fabs(x - 0.875) <= 0.125 && fabs(y - 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

          computed_QoI[0] = computed_QoI[0] + dQoI_0;
        
        }



The source file element_qoi_derivative.C with comments:

        #include "libmesh/libmesh_common.h"
        #include "libmesh/elem.h"
        #include "libmesh/fe_base.h"
        #include "libmesh/fem_context.h"
        #include "libmesh/point.h"
        #include "libmesh/quadrature.h"
        
Local includes
        #include "L-shaped.h"
        
Bring in everything from the libMesh namespace
        using namespace libMesh;
        
We only have one QoI, so we don't bother checking the qois argument to see if it was requested from us
        void LaplaceSystem::element_qoi_derivative (DiffContext &context,
                                                    const QoISet & /* qois */)
        {
          FEMContext &c = libmesh_cast_ref<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 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(reuse_preconditioner);
            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_INPUT(patch_reuse);
        
            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-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");
          exact_QoI[0] = infile("QoI_0", 0.0);
          exact_QoI[1] = infile("QoI_1", 0.0);
        
Do the parent's initialization after variables are defined
          FEMSystem::init_data();
        
          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 functions
          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] * ( 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] * ( 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;
        }
        
Override the default DiffSystem postprocess function to compute the approximations to the QoIs
        void LaplaceSystem::postprocess()
        {
Reset the array holding the computed QoIs
          computed_QoI[0] = 0.0;
          computed_QoI[1] = 0.0;
        
          FEMSystem::postprocess();
        
          this->comm().sum(computed_QoI[0]);
        
          this->comm().sum(computed_QoI[1]);
        
        }
        
The exact solution to the singular problem, u_exact = r^(2/3)*sin(2*theta/3). We use this to set the Dirichlet boundary conditions
        Number LaplaceSystem::exact_solution(const Point& p)// xyz location
        {
          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 pow(x1*x1 + x2*x2, 1./3.)*sin(2./3.*theta);
        
        }



The source file side_postprocess.C with comments:

        #include "libmesh/libmesh_common.h"
        #include "libmesh/elem.h"
        #include "libmesh/fe_base.h"
        #include "libmesh/fem_context.h"
        #include "libmesh/point.h"
        #include "libmesh/quadrature.h"
        
Local includes
        #include "L-shaped.h"
        
Bring in everything from the libMesh namespace
        using namespace libMesh;
        
Define the postprocess function to compute QoI 1, the integral of the the normal derivative of the solution over part of the boundary

        void LaplaceSystem::side_postprocess(DiffContext &context)
        {
        
          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();
        
          const std::vector<Point > &q_point = c.side_fe_var[0]->get_xyz();
        
          const std::vector<Point> &face_normals = c.side_fe_var[0]->get_normals();
        
          unsigned int n_qpoints = (c.get_side_qrule())->n_points();
        
          Number dQoI_1 = 0. ;
        
Loop over qp's, compute the function at each qp and add to get the QoI
          for (unsigned int qp=0; qp != n_qpoints; qp++)
            {
              const Real x = q_point[qp](0);
              const Real y = q_point[qp](1);
        
              const Real TOL = 1.e-5;
        
If on the bottom horizontal bdry (y = -1)
              if(fabs(y - 1.0) <= TOL && x > 0.0)
        	{
Get the value of the gradient at this point
                  const Gradient grad_T = c.side_gradient(0,qp);
        
Add the contribution of this qp to the integral QoI
                  dQoI_1 += JxW[qp] * (grad_T * face_normals[qp] * x * (x - 1.)) ;
        	}
        
            } // end of the quadrature point qp-loop
        
          computed_QoI[1] = computed_QoI[1] + dQoI_1;
        }



The source file side_qoi_derivative.C with comments:

        #include "libmesh/libmesh_common.h"
        #include "libmesh/elem.h"
        #include "libmesh/fe_base.h"
        #include "libmesh/fem_context.h"
        #include "libmesh/point.h"
        #include "libmesh/quadrature.h"
        
Local includes
        #include "L-shaped.h"
        
Bring in everything from the libMesh namespace
        using namespace libMesh;
        
We only have one QoI, so we don't bother checking the qois argument to see if it was requested from us
        void LaplaceSystem::side_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.side_fe_var[0]->get_JxW();
        
Basis Functions
          const std::vector<std::vector<RealGradient> > &dphi = c.side_fe_var[0]->get_dphi();
        
The side quadrature points
          const std::vector<Point > &q_point = c.side_fe_var[0]->get_xyz();
        
Get the normal to the side at each qp
          const std::vector<Point> &face_normals = c.side_fe_var[0]->get_normals();
        
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_side_qrule())->n_points();
        
Fill the QoI RHS corresponding to this QoI. Since this is QoI 1 we fill in the [1][i] subderivatives, i corresponding to the variable index. Our system has only one variable, so we only have to fill the [1][0] subderivative
          DenseSubVector<Number> &Q1 = *c.elem_qoi_subderivatives[1][0];
        
          const Real TOL = 1.e-5;
        
          for (unsigned int qp=0; qp != n_qpoints; qp++)
            {
              const Real x = q_point[qp](0);
              const Real y = q_point[qp](1);
        
If on the sides where the boundary QoI is supported, add contributions to the adjoint rhs
              if(fabs(y - 1.0) <= TOL && x > 0.0)
          	{
          	  for (unsigned int i=0; i != n_T_dofs; i++)
        	    {
        	      Q1(i) += JxW[qp] * (dphi[i][qp] * face_normals[qp]) * x * (x - 1.);
        	    }
          	}
        
            } // end of the quadrature point qp-loop
        }



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),
  	reuse_preconditioner(false),
        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"), patch_reuse(true),
        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);
  
      std::string domainfile;
      unsigned int coarserefinements;
  
      bool solver_quiet, reuse_preconditioner, 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;
      bool patch_reuse;
  
      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-shaped.h without comments:

 
  #include "libmesh/enum_fe_family.h"
  #include "libmesh/fem_system.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) { qoi.resize(2); }
  
    std::string & fe_family() { return _fe_family;  }
    unsigned int & fe_order() { return _fe_order;  }
    bool & analytic_jacobians() { return _analytic_jacobians; }
  
  
    virtual void postprocess(void);
  
    Number &get_QoI_value(std::string type, unsigned int QoI_index)
      {
        if(type == "exact")
  	{
  	  return exact_QoI[QoI_index];
  	}
        else
  	{
  	  return computed_QoI[QoI_index];
  	}
      }
  
    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);
  
  
    virtual void element_postprocess(DiffContext &context);
  
    virtual void side_postprocess(DiffContext &context);
  
  
    virtual void element_qoi_derivative
      (DiffContext &context,
       const QoISet & qois);
  
  
    virtual void side_qoi_derivative
      (DiffContext &context,
       const QoISet & qois);
  
    Number exact_solution (const Point&);
  
  
    Number computed_QoI[2];
  
  
    Number exact_QoI[2];
  
    std::string _fe_family;
    unsigned int _fe_order;
  
    bool _analytic_jacobians;
  };



The source file adjoints_ex4.C without comments:

 
  
  
  #include <iostream>
  #include <iomanip>
  
  #include "libmesh/equation_systems.h"
  #include "libmesh/linear_solver.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/qoi_set.h"
  #include "libmesh/adjoint_refinement_estimator.h"
  
  #include "libmesh/getpot.h"
  #include "libmesh/gmv_io.h"
  
  #include "femparameters.h"
  #include "L-shaped.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;
  }
  
  AutoPtr<AdjointRefinementEstimator> build_adjoint_refinement_error_estimator(QoISet &qois)
  {
    AutoPtr<AdjointRefinementEstimator> error_estimator;
  
    std::cout<<"Computing the error estimate using the Adjoint Refinement Error Estimator"<<std::endl<<std::endl;
  
    AdjointRefinementEstimator *adjoint_refinement_estimator = new AdjointRefinementEstimator;
  
    error_estimator.reset (adjoint_refinement_estimator);
  
    adjoint_refinement_estimator->qoi_set() = qois;
  
    adjoint_refinement_estimator->number_h_refinements = 2;
  
    return error_estimator;
  }
  
  #endif // LIBMESH_ENABLE_AMR
  
  
  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");
  
    set_system_parameters(system, param);
  
    std::cout << "Initializing systems" << std::endl;
  
    equation_systems.init ();
  
    mesh.print_info();
    equation_systems.print_info();
    LinearSolver<Number> *linear_solver = system.get_linear_solver();
  
    {
      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);
  
  	linear_solver->reuse_preconditioner(false);
  
  	system.solve();
  
  	write_output(equation_systems, a_step, "primal");
  
  	NumericVector<Number> &primal_solution = *system.solution;
  
  	QoISet qois;
  
  	std::vector<unsigned int> 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);
  
  	system.assemble_qoi_sides = true;
  
  	linear_solver->reuse_preconditioner(param.reuse_preconditioner);
  
  	system.adjoint_solve();
  
  	system.set_adjoint_already_solved(true);
  
  	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);
  
  	NumericVector<Number> &dual_solution_1 = system.get_adjoint_solution(1);
  
  	primal_solution.swap(dual_solution_1);
  	write_output(equation_systems, a_step, "adjoint_1");
  
  	primal_solution.swap(dual_solution_1);
  
  	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 << "Postprocessing: " << std::endl;
  	system.postprocess_sides = true;
  	system.postprocess();
  	Number QoI_0_computed = system.get_QoI_value("computed", 0);
  	Number QoI_0_exact = system.get_QoI_value("exact", 0);
  	Number QoI_1_computed = system.get_QoI_value("computed", 1);
  	Number QoI_1_exact = system.get_QoI_value("exact", 1);
  
  	std::cout<< "The relative error in QoI 0 is " << std::setprecision(17)
                   << std::abs(QoI_0_computed - QoI_0_exact) /
                      std::abs(QoI_0_exact) << std::endl;
  
  	std::cout<< "The relative error in QoI 1 is " << std::setprecision(17)
                   << std::abs(QoI_1_computed - QoI_1_exact) /
                      std::abs(QoI_1_exact) << std::endl << std::endl;
  
  	ErrorVector QoI_elementwise_error;
  
  	AutoPtr<AdjointRefinementEstimator> adjoint_refinement_error_estimator =
  	  build_adjoint_refinement_error_estimator(qois);
  
  	adjoint_refinement_error_estimator->estimate_error(system, QoI_elementwise_error);
  
  	std::cout<< "The computed relative error in QoI 0 is " << std::setprecision(17)
                   << std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(0)) /
                      std::abs(QoI_0_exact) << std::endl;
  
  	std::cout<< "The computed relative error in QoI 1 is " << std::setprecision(17)
                   << std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(1)) /
  	  std::abs(QoI_1_exact) << std::endl << std::endl;
  
  	std::cout<< "The effectivity index for the computed error in QoI 0 is " << std::setprecision(17)
                   << std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(0)) /
                      std::abs(QoI_0_computed - QoI_0_exact) << std::endl;
  
  	std::cout<< "The effectivity index for the computed error in QoI 1 is " << std::setprecision(17)
                   << std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(1)) /
  	  std::abs(QoI_1_computed - QoI_1_exact) << std::endl << std::endl;
  
  
  	if(param.refine_uniformly)
  	  {
  	    mesh_refinement->uniformly_refine(1);
  	  }
  	else if(param.global_tolerance >= 0. && param.nelem_target == 0.)
  	  {
  	    mesh_refinement->flag_elements_by_error_tolerance (QoI_elementwise_error);
  
  	    mesh_refinement->refine_and_coarsen_elements();
  	  }
  	else
  	  {
  	    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 (QoI_elementwise_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)
        {
  	linear_solver->reuse_preconditioner(false);
  	system.solve();
  
  	write_output(equation_systems, a_step, "primal");
  
  	NumericVector<Number> &primal_solution = *system.solution;
  
  	QoISet qois;
  	std::vector<unsigned int> 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);
  
  	system.assemble_qoi_sides = true;
  	linear_solver->reuse_preconditioner(param.reuse_preconditioner);
  	system.adjoint_solve();
  
  	system.set_adjoint_already_solved(true);
  
  	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);
  
  	NumericVector<Number> &dual_solution_1 = system.get_adjoint_solution(1);
  
  	primal_solution.swap(dual_solution_1);
  	write_output(equation_systems, a_step, "adjoint_1");
  
  	primal_solution.swap(dual_solution_1);
  
  	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 << "Postprocessing: " << std::endl;
  	system.postprocess_sides = true;
  	system.postprocess();
  
  	Number QoI_0_computed = system.get_QoI_value("computed", 0);
  	Number QoI_0_exact = system.get_QoI_value("exact", 0);
  	Number QoI_1_computed = system.get_QoI_value("computed", 1);
  	Number QoI_1_exact = system.get_QoI_value("exact", 1);
  
  	std::cout<< "The relative error in QoI 0 is " << std::setprecision(17)
                   << std::abs(QoI_0_computed - QoI_0_exact) /
                      std::abs(QoI_0_exact) << std::endl;
  
  	std::cout<< "The relative error in QoI 1 is " << std::setprecision(17)
                   << std::abs(QoI_1_computed - QoI_1_exact) /
                      std::abs(QoI_1_exact) << std::endl << std::endl;
  
  	ErrorVector QoI_elementwise_error;
  
  	AutoPtr<AdjointRefinementEstimator> adjoint_refinement_error_estimator =
  	  build_adjoint_refinement_error_estimator(qois);
  
  	adjoint_refinement_error_estimator->estimate_error(system, QoI_elementwise_error);
  
  	std::cout<< "The computed relative error in QoI 0 is " << std::setprecision(17)
                   << std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(0)) /
                      std::abs(QoI_0_exact) << std::endl;
  
  	std::cout<< "The computed relative error in QoI 1 is " << std::setprecision(17)
                   << std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(1)) /
  	  std::abs(QoI_1_exact) << std::endl << std::endl;
  
  	std::cout<< "The effectivity index for the computed error in QoI 0 is " << std::setprecision(17)
                   << std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(0)) /
                      std::abs(QoI_0_computed - QoI_0_exact) << std::endl;
  
  	std::cout<< "The effectivity index for the computed error in QoI 1 is " << std::setprecision(17)
                   << std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(1)) /
  	  std::abs(QoI_1_computed - QoI_1_exact) << std::endl << std::endl;
  
        }
    }
  
    std::cerr << '[' << libMesh::processor_id()
              << "] Completing output." << std::endl;
  
  #endif // #ifndef LIBMESH_ENABLE_AMR
  
    return 0;
  }



The source file element_postprocess.C without comments:

 
  #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 "L-shaped.h"
  
  using namespace libMesh;
  
  
  
  void LaplaceSystem::element_postprocess (DiffContext &context)
  
  {
    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 x = xyz[qp](0);
        const Real y = xyz[qp](1);
  
        if(fabs(x - 0.875) <= 0.125 && fabs(y - 0.125) <= 0.125)
    	{
    	  Number T = c.interior_value(0, qp);
  
    	  dQoI_0 += JxW[qp] * T;
    	}
      }
  
  
    computed_QoI[0] = computed_QoI[0] + dQoI_0;
  
  }



The source file element_qoi_derivative.C without comments:

 
  #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 "L-shaped.h"
  
  using namespace libMesh;
  
  void LaplaceSystem::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 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(reuse_preconditioner);
      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_INPUT(patch_reuse);
  
      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-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");
    exact_QoI[0] = infile("QoI_0", 0.0);
    exact_QoI[1] = infile("QoI_1", 0.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] * ( 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] * ( 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;
  }
  
  void LaplaceSystem::postprocess()
  {
    computed_QoI[0] = 0.0;
    computed_QoI[1] = 0.0;
  
    FEMSystem::postprocess();
  
    this->comm().sum(computed_QoI[0]);
  
    this->comm().sum(computed_QoI[1]);
  
  }
  
  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 pow(x1*x1 + x2*x2, 1./3.)*sin(2./3.*theta);
  
  }



The source file side_postprocess.C without comments:

 
  #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 "L-shaped.h"
  
  using namespace libMesh;
  
  
  void LaplaceSystem::side_postprocess(DiffContext &context)
  {
  
    FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
  
  
    const std::vector<Real> &JxW = c.side_fe_var[0]->get_JxW();
  
    const std::vector<Point > &q_point = c.side_fe_var[0]->get_xyz();
  
    const std::vector<Point> &face_normals = c.side_fe_var[0]->get_normals();
  
    unsigned int n_qpoints = (c.get_side_qrule())->n_points();
  
    Number dQoI_1 = 0. ;
  
    for (unsigned int qp=0; qp != n_qpoints; qp++)
      {
        const Real x = q_point[qp](0);
        const Real y = q_point[qp](1);
  
        const Real TOL = 1.e-5;
  
        if(fabs(y - 1.0) <= TOL && x > 0.0)
  	{
  	  const Gradient grad_T = c.side_gradient(0,qp);
  
  	  dQoI_1 += JxW[qp] * (grad_T * face_normals[qp] * x * (x - 1.)) ;
  	}
  
      } // end of the quadrature point qp-loop
  
    computed_QoI[1] = computed_QoI[1] + dQoI_1;
  }



The source file side_qoi_derivative.C without comments:

 
  #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 "L-shaped.h"
  
  using namespace libMesh;
  
  void LaplaceSystem::side_qoi_derivative (DiffContext &context,
                                           const QoISet & /* qois */)
  {
    FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
  
  
    const std::vector<Real> &JxW = c.side_fe_var[0]->get_JxW();
  
    const std::vector<std::vector<RealGradient> > &dphi = c.side_fe_var[0]->get_dphi();
  
    const std::vector<Point > &q_point = c.side_fe_var[0]->get_xyz();
  
    const std::vector<Point> &face_normals = c.side_fe_var[0]->get_normals();
  
    const unsigned int n_T_dofs = c.dof_indices_var[0].size();
    unsigned int n_qpoints = (c.get_side_qrule())->n_points();
  
    DenseSubVector<Number> &Q1 = *c.elem_qoi_subderivatives[1][0];
  
    const Real TOL = 1.e-5;
  
    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(y - 1.0) <= TOL && x > 0.0)
    	{
    	  for (unsigned int i=0; i != n_T_dofs; i++)
  	    {
  	      Q1(i) += JxW[qp] * (dphi[i][qp] * face_normals[qp]) * x * (x - 1.);
  	    }
    	}
  
      } // end of the quadrature point qp-loop
  }



The console output of the program:

make[4]: Entering directory `/net/spark/workspace/roystgnr/libmesh/git/devel/examples/adjoints/adjoints_ex4'
***************************************************************
* Running Example adjoints_ex4:
*  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_ex4/.libs/lt-example-devel
Reading in and building the mesh
Building system
Initializing systems
 Mesh Information:
  mesh_dimension()=2
  spatial_dimension()=3
  n_nodes()=65
    n_local_nodes()=21
  n_elem()=15
    n_local_elem()=5
    n_active_elem()=12
  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()=65
    n_local_dofs()=21
    n_constrained_dofs()=0
    n_local_constrained_dofs()=0
    n_vectors()=1
    n_matrices()=1
    DofMap Sparsity
      Average  On-Processor Bandwidth <= 10.2923
      Average Off-Processor Bandwidth <= 3.47692
      Maximum  On-Processor Bandwidth <= 21
      Maximum Off-Processor Bandwidth <= 16
    DofMap Constraints
      Number of DoF Constraints = 0
      Number of Node Constraints = 0

  Nonlinear solver converged, step 0, residual reduction 8.99985e-11 < 1e-09
Adaptive step 0, we have 12 active elements and 65 active dofs.
Postprocessing: 
The relative error in QoI 0 is inf
The relative error in QoI 1 is inf

Computing the error estimate using the Adjoint Refinement Error Estimator

The computed relative error in QoI 0 is inf
The computed relative error in QoI 1 is inf

The effectivity index for the computed error in QoI 0 is 6.2029227321009616e-05
The effectivity index for the computed error in QoI 1 is 0.0045983013122305322

Refined mesh to 48 active elements and 225 active dofs.
  Nonlinear solver converged, step 0, residual reduction 6.3455558581624914e-11 < 1.0000000000000001e-09
Adaptive step 1, we have 48 active elements and 225 active dofs.
Postprocessing: 
The relative error in QoI 0 is inf
The relative error in QoI 1 is inf

Computing the error estimate using the Adjoint Refinement Error Estimator

The computed relative error in QoI 0 is inf
The computed relative error in QoI 1 is inf

The effectivity index for the computed error in QoI 0 is 0.00022987680218270143
The effectivity index for the computed error in QoI 1 is 0.0005945107583468955

Refined mesh to 192 active elements and 833 active dofs.
  Nonlinear solver converged, step 0, residual reduction 5.5183873316805855e-11 < 1.0000000000000001e-09
Adaptive step 2, we have 192 active elements and 833 active dofs.
Postprocessing: 
The relative error in QoI 0 is inf
The relative error in QoI 1 is inf

Computing the error estimate using the Adjoint Refinement Error Estimator

The computed relative error in QoI 0 is inf
The computed relative error in QoI 1 is inf

The effectivity index for the computed error in QoI 0 is 0.00010471304480264564
The effectivity index for the computed error in QoI 1 is 5.8173302203917334e-05


 -------------------------------------------------------------------------------------------------------------------
| Processor id:   0                                                                                                 |
| Num Processors: 4                                                                                                 |
| Time:           Fri Apr 19 11:48:21 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=2.00285, Active time=1.90319                                                   |
 ----------------------------------------------------------------------------------------------------------------
| 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()     15        0.0181      0.001207    0.0246      0.001641    0.95     1.29     |
|   build_sparsity()                 15        0.0123      0.000822    0.0320      0.002135    0.65     1.68     |
|   create_dof_constraints()         15        0.0089      0.000595    0.0089      0.000595    0.47     0.47     |
|   distribute_dofs()                15        0.0647      0.004315    0.1662      0.011082    3.40     8.73     |
|   dof_indices()                    15806     0.1603      0.000010    0.1603      0.000010    8.42     8.42     |
|   old_dof_indices()                11070     0.1055      0.000010    0.1055      0.000010    5.54     5.54     |
|   prepare_send_list()              15        0.0002      0.000015    0.0002      0.000015    0.01     0.01     |
|   reinit()                         15        0.0896      0.005974    0.0896      0.005974    4.71     4.71     |
|                                                                                                                |
| EquationSystems                                                                                                |
|   build_solution_vector()          9         0.0012      0.000132    0.0048      0.000538    0.06     0.25     |
|                                                                                                                |
| FE                                                                                                             |
|   compute_shape_functions()        7426      0.1281      0.000017    0.1281      0.000017    6.73     6.73     |
|   init_shape_functions()           796       0.0101      0.000013    0.0101      0.000013    0.53     0.53     |
|   inverse_map()                    19260     0.1136      0.000006    0.1136      0.000006    5.97     5.97     |
|                                                                                                                |
| FEMSystem                                                                                                      |
|   assemble_qoi_derivative()        6         0.0363      0.006049    0.1273      0.021222    1.91     6.69     |
|   assembly()                       3         0.0025      0.000835    0.0139      0.004646    0.13     0.73     |
|   assembly(get_jacobian)           6         0.0359      0.005977    0.2499      0.041653    1.88     13.13    |
|   assembly(get_residual)           6         0.0277      0.004614    0.1178      0.019635    1.45     6.19     |
|   numerical_elem_jacobian()        1134      0.1210      0.000107    0.1210      0.000107    6.36     6.36     |
|   numerical_side_jacobian()        132       0.0061      0.000046    0.0061      0.000046    0.32     0.32     |
|   postprocess()                    3         0.0017      0.000557    0.0095      0.003166    0.09     0.50     |
|                                                                                                                |
| FEMap                                                                                                          |
|   compute_affine_map()             7426      0.0427      0.000006    0.0427      0.000006    2.24     2.24     |
|   compute_face_map()               748       0.0086      0.000012    0.0227      0.000030    0.45     1.19     |
|   init_face_shape_functions()      48        0.0002      0.000004    0.0002      0.000004    0.01     0.01     |
|   init_reference_to_physical_map() 796       0.0109      0.000014    0.0109      0.000014    0.57     0.57     |
|                                                                                                                |
| GMVIO                                                                                                          |
|   write_nodal_data()               9         0.0253      0.002811    0.0258      0.002869    1.33     1.36     |
|                                                                                                                |
| ImplicitSystem                                                                                                 |
|   adjoint_solve()                  6         0.0008      0.000133    0.7429      0.123813    0.04     39.03    |
|                                                                                                                |
| LocationMap                                                                                                    |
|   find()                           26460     0.0422      0.000002    0.0422      0.000002    2.22     2.22     |
|   init()                           23        0.0165      0.000718    0.0165      0.000718    0.87     0.87     |
|                                                                                                                |
| Mesh                                                                                                           |
|   all_second_order()               1         0.0001      0.000095    0.0001      0.000095    0.00     0.00     |
|   contract()                       14        0.0079      0.000568    0.0241      0.001720    0.42     1.27     |
|   find_neighbors()                 17        0.0580      0.003411    0.0614      0.003615    3.05     3.23     |
|   renumber_nodes_and_elem()        14        0.0161      0.001152    0.0161      0.001152    0.85     0.85     |
|                                                                                                                |
| MeshCommunication                                                                                              |
|   broadcast()                      1         0.0002      0.000232    0.0005      0.000529    0.01     0.03     |
|   compute_hilbert_indices()        6         0.0008      0.000133    0.0008      0.000133    0.04     0.04     |
|   find_global_indices()            6         0.0003      0.000051    0.0029      0.000481    0.02     0.15     |
|   parallel_sort()                  6         0.0006      0.000101    0.0011      0.000181    0.03     0.06     |
|                                                                                                                |
| MeshOutput                                                                                                     |
|   write_equation_systems()         9         0.0002      0.000026    0.0315      0.003496    0.01     1.65     |
|                                                                                                                |
| MeshRefinement                                                                                                 |
|   _coarsen_elements()              20        0.0113      0.000565    0.0121      0.000603    0.59     0.63     |
|   _refine_elements()               23        0.0531      0.002307    0.1560      0.006785    2.79     8.20     |
|   add_point()                      26460     0.0521      0.000002    0.0973      0.000004    2.74     5.11     |
|   make_coarsening_compatible()     28        0.0309      0.001102    0.0309      0.001102    1.62     1.62     |
|   make_flags_parallel_consistent() 28        0.0159      0.000569    0.0248      0.000884    0.84     1.30     |
|   make_refinement_compatible()     28        0.0001      0.000002    0.0002      0.000006    0.00     0.01     |
|                                                                                                                |
| MetisPartitioner                                                                                               |
|   partition()                      5         0.0034      0.000674    0.0062      0.001242    0.18     0.33     |
|                                                                                                                |
| NewtonSolver                                                                                                   |
|   solve()                          3         0.0149      0.004959    0.0514      0.017149    0.78     2.70     |
|                                                                                                                |
| Parallel                                                                                                       |
|   allgather()                      59        0.0051      0.000086    0.0051      0.000087    0.27     0.27     |
|   broadcast()                      27        0.0001      0.000005    0.0001      0.000003    0.01     0.00     |
|   max(bool)                        73        0.0039      0.000054    0.0039      0.000054    0.21     0.21     |
|   max(scalar)                      2454      0.0096      0.000004    0.0096      0.000004    0.50     0.50     |
|   max(vector)                      594       0.0041      0.000007    0.0105      0.000018    0.22     0.55     |
|   min(bool)                        3077      0.0117      0.000004    0.0117      0.000004    0.61     0.61     |
|   min(scalar)                      2405      0.0192      0.000008    0.0192      0.000008    1.01     1.01     |
|   min(vector)                      594       0.0045      0.000008    0.0111      0.000019    0.24     0.58     |
|   probe()                          672       0.0035      0.000005    0.0035      0.000005    0.19     0.19     |
|   receive()                        672       0.0023      0.000003    0.0059      0.000009    0.12     0.31     |
|   send()                           672       0.0014      0.000002    0.0014      0.000002    0.07     0.07     |
|   send_receive()                   684       0.0029      0.000004    0.0111      0.000016    0.15     0.58     |
|   sum()                            124       0.0018      0.000015    0.0022      0.000018    0.09     0.12     |
|                                                                                                                |
| Parallel::Request                                                                                              |
|   wait()                           672       0.0007      0.000001    0.0007      0.000001    0.04     0.04     |
|                                                                                                                |
| Partitioner                                                                                                    |
|   set_node_processor_ids()         5         0.0006      0.000122    0.0020      0.000406    0.03     0.11     |
|   set_parent_processor_ids()       5         0.0001      0.000023    0.0001      0.000023    0.01     0.01     |
|                                                                                                                |
| PetscLinearSolver                                                                                              |
|   solve()                          15        0.3746      0.024975    0.3746      0.024975    19.68    19.68    |
|                                                                                                                |
| ProjectVector                                                                                                  |
|   operator()                       36        0.0652      0.001811    0.2591      0.007197    3.42     13.61    |
|                                                                                                                |
| System                                                                                                         |
|   project_vector()                 36        0.0348      0.000966    0.3606      0.010016    1.83     18.94    |
 ----------------------------------------------------------------------------------------------------------------
| Totals:                            130808    1.9032                                          100.00            |
 ----------------------------------------------------------------------------------------------------------------

 
***************************************************************
* Done Running Example adjoints_ex4:
*  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_ex4'

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

Hosted By:
SourceForge.net Logo