The source file coupled_system.h with comments:

        #ifndef __coupled_system_h__
        #define __coupled_system_h__
        
DiffSystem framework files
        #include "libmesh/fem_function_base.h"
        #include "libmesh/fem_system.h"
        #include "libmesh/libmesh_common.h"
        #include "libmesh/parameter_vector.h"
        
        using namespace libMesh;
        
The Navier-Stokes system class. FEMSystem, TimeSolver and NewtonSolver will handle most tasks, but we must specify element residuals
        class CoupledSystem : public FEMSystem
        {
        public:
Constructor
          CoupledSystem(EquationSystems& es,
                       const std::string& name_in,
                       const unsigned int number_in)
            : FEMSystem(es, name_in, number_in), Peclet(1.) {qoi.resize(1);}
        
Function to get computed QoI values

          Number &get_QoI_value()
            {
              return computed_QoI;
            }
        
          Number &get_parameter_value(unsigned int parameter_index)
            {
              return parameters[parameter_index];
            }
        
          ParameterVector &get_parameter_vector()
            {
              parameter_vector.resize(parameters.size());
              for(unsigned int i = 0; i != parameters.size(); ++i)
        	{
        	  parameter_vector[i] = ¶meters[i];
        	}
        
              return parameter_vector;
            }
        
          Real &get_Pe()
            {
              return Peclet;
            }
        
         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 element_constraint (bool request_jacobian,
                                           DiffContext& context);
        
Postprocessed output
          virtual void postprocess ();
        
Parameters associated with the system
          std::vector<Number> parameters;
        
Indices for each variable;
          unsigned int p_var, u_var, v_var, C_var;
        
The ParameterVector object that will contain pointers to the system parameters
          ParameterVector parameter_vector;
        
The Peclet number for the species transport
          Real Peclet;
        
The functionals to be computed as QoIs
          Number computed_QoI;
        
        };
        
        
        class CoupledFEMFunctionsx : public FEMFunctionBase<Number>
        {
        public:
Constructor
          CoupledFEMFunctionsx(System& /* sys */, unsigned int var_number) {var = var_number;}
        
Destructor
          virtual ~CoupledFEMFunctionsx () {}
        
          virtual AutoPtr<FEMFunctionBase<Number> > clone () const
          {return AutoPtr<FEMFunctionBase<Number> >( new CoupledFEMFunctionsx(*this) ); }
        
          virtual void operator() (const FEMContext&, const Point&,
        			   const Real, DenseVector<Number>&)
          {libmesh_error();}
        
          virtual Number operator() (const FEMContext&, const Point& p,
        			     const Real time = 0.);
        
         private:
        
          unsigned int var;
        
        };
        
        
        class CoupledFEMFunctionsy : public FEMFunctionBase<Number>
        {
        public:
Constructor
          CoupledFEMFunctionsy(System& /* sys */, unsigned int var_number) {var = var_number;}
        
Destructor
          virtual ~CoupledFEMFunctionsy () {}
        
          virtual AutoPtr<FEMFunctionBase<Number> > clone () const
          {return AutoPtr<FEMFunctionBase<Number> >( new CoupledFEMFunctionsy(*this) ); }
        
          virtual void operator() (const FEMContext&, const Point&,
        			   const Real,
        			   DenseVector<Number>&)
          {libmesh_error();}
        
          virtual Number operator() (const FEMContext&, const Point& p,
        			     const Real time = 0.);
        
         private:
        
          unsigned int var;
        
        };
        
        #endif //__coupled_system_h__



The source file domain.h with comments:

        #include "femparameters.h"
        #include "libmesh/mesh.h"
        
Bring in everything from the libMesh namespace
        class FEMParameters;
        
        void build_domain (Mesh &mesh, FEMParameters ¶m);



The source file femparameters.h with comments:

        #ifndef __fem_parameters_h__
        #define __fem_parameters_h__
        
        #include <limits>
        
        #include "libmesh/libmesh_common.h"
        #include "libmesh/dof_map.h"
        #include "libmesh/enum_norm_type.h"
        #include "libmesh/getpot.h"
        
Bring in everything from the libMesh namespace
        using namespace libMesh;
        
        class FEMParameters
        {
        public:
            FEMParameters() :
              initial_timestep(0), n_timesteps(100),
              transient(true),
              deltat_reductions(0),
              timesolver_core("euler"),
              end_time(std::numeric_limits<Real>::max()),
              deltat(0.0001), timesolver_theta(0.5),
              timesolver_maxgrowth(0.), timesolver_tolerance(0.),
              timesolver_upper_tolerance(0.),
              steadystate_tolerance(0.),
              timesolver_norm(0,L2),
              dimension(2),
        	domaintype("square"), domainfile("mesh.xda"), elementtype("quad"),
        	fine_mesh_file_primal("fine_mesh.xda"), fine_mesh_soln_primal("fine_mesh_soln.xda"),
        	fine_mesh_file_adjoint("fine_mesh.xda"), fine_mesh_soln_adjoint("fine_mesh_soln.xda"),
              elementorder(2),
              domain_xmin(0.0), domain_ymin(0.0), domain_zmin(0.0),
              domain_edge_width(1.0), domain_edge_length(1.0), domain_edge_height(1.0),
              coarsegridx(1), coarsegridy(1), coarsegridz(1),
        	coarserefinements(0), coarsecoarsenings(0), extrarefinements(0),
              use_petsc_snes(false),
              time_solver_quiet(true), 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),
              initial_adaptivesteps(0), initial_sobolev_order(1),
        	initial_extra_quadrature(0),
        	indicator_type("kelly"), patch_reuse(true), sobolev_order(1), adjoint_residual_type("patchpatch"), alternate_with_uniform_steps("false"), alternate_step_number(2),
        	component_wise_error("false"), compare_to_fine_solution_primal(false), compare_to_fine_solution_adjoint(false), compute_sensitivities(false), do_forward_sensitivity(true), do_adjoint_sensitivity(false),
        	postprocess_adjoint(false),
              write_interval(10),
              write_gmv_error(false), write_tecplot_error(false),
              output_xda(false), output_xdr(false),
              output_bz2(true), output_gz(true),
              output_gmv(false), output_tecplot(false),
              run_simulation(true), run_postprocess(false),
              fe_family(1, "LAGRANGE"), fe_order(1, 1),
              extra_quadrature_order(0),
              analytic_jacobians(true), verify_analytic_jacobians(0.0),
              print_solution_norms(false), print_solutions(false),
              print_residual_norms(false), print_residuals(false),
              print_jacobian_norms(false), print_jacobians(false) {}
        
            void read(GetPot &input);
        
        
        
            unsigned int initial_timestep, n_timesteps;
            bool transient;
            unsigned int deltat_reductions;
            std::string timesolver_core;
            Real end_time, deltat, timesolver_theta,
        	 timesolver_maxgrowth, timesolver_tolerance,
        	 timesolver_upper_tolerance, steadystate_tolerance;
            std::vector<FEMNormType> timesolver_norm;
        
            unsigned int dimension;
            std::string domaintype, domainfile, elementtype, fine_mesh_file_primal, fine_mesh_soln_primal, fine_mesh_file_adjoint, fine_mesh_soln_adjoint;
            Real elementorder;
            Real domain_xmin, domain_ymin, domain_zmin;
            Real domain_edge_width, domain_edge_length, domain_edge_height;
            unsigned int coarsegridx, coarsegridy, coarsegridz;
            unsigned int coarserefinements, coarsecoarsenings, extrarefinements;
        
            bool use_petsc_snes;
            bool time_solver_quiet, 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;
            unsigned int initial_adaptivesteps;
            unsigned int initial_sobolev_order;
            unsigned int initial_extra_quadrature;
        
            std::string indicator_type;
            bool patch_reuse;
            unsigned int sobolev_order;
            std::string adjoint_residual_type;
            bool alternate_with_uniform_steps;
            unsigned int alternate_step_number;
            bool component_wise_error;
            bool compare_to_fine_solution_primal, compare_to_fine_solution_adjoint;
        
            bool compute_sensitivities;
            bool do_forward_sensitivity;
            bool do_adjoint_sensitivity;
        
            bool postprocess_adjoint;
        
            unsigned int write_interval;
            bool write_gmv_error, write_tecplot_error, output_xda, output_xdr,
        	 output_bz2, output_gz, output_gmv, output_tecplot;
            bool run_simulation, run_postprocess;
        
            std::vector<std::string> fe_family;
            std::vector<unsigned int> fe_order;
            int extra_quadrature_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 H-qoi.h with comments:

        #ifndef H_QOI_H
        #define H_QOI_H
        
        #include "libmesh/libmesh_common.h"
        #include "libmesh/elem.h"
        #include "libmesh/fe_base.h"
        #include "libmesh/fem_context.h"
        #include "libmesh/point.h"
        #include "libmesh/quadrature.h"
        #include "libmesh/diff_qoi.h"
        
Bring in everything from the libMesh namespace
        using namespace libMesh;
        
        class CoupledSystemQoI : public DifferentiableQoI
        {
        public:
          CoupledSystemQoI(){}
          virtual ~CoupledSystemQoI(){}
        
          virtual void init_qoi( std::vector<Number>& sys_qoi);
          virtual void postprocess( ){}
        
          virtual void side_qoi_derivative(DiffContext &context, const QoISet & qois);
        
          virtual void side_qoi(DiffContext &context, const QoISet & qois);
        
          virtual AutoPtr<DifferentiableQoI> clone( )
          { AutoPtr<DifferentiableQoI> my_clone( new CoupledSystemQoI );
            *my_clone = *this;
            return my_clone;
          }
        
        };
        #endif // H_QOI_H



The source file initial.h with comments:

        #include "libmesh/parameters.h"
        #include "libmesh/point.h"
        #include "libmesh/vector_value.h"
        
Bring in everything from the libMesh namespace
        using namespace libMesh;
        
        void read_initial_parameters();
        void finish_initialization();
        
        Number initial_value(const Point& p,
                             const Parameters&,
                             const std::string&,
                             const std::string&);
        
        Gradient initial_grad(const Point& p,
                              const Parameters&,
                              const std::string&,
                              const std::string&);



The source file mysystems.h with comments:

        #include "coupled_system.h"
        #include "femparameters.h"
        
Bring in everything from the libMesh namespace
        using namespace libMesh;
        
        FEMSystem &build_system(EquationSystems &es, GetPot &, FEMParameters& /* param */)
        {
          CoupledSystem &system = es.add_system<CoupledSystem> ("CoupledSystem");
        
          return system;
        }
        
        Number exact_value(const Point& /* p */,       // xyz location
                           const Parameters& /* param */,  // EquationSystem parameters
                           const std::string&, // sys_name
                           const std::string&) // unknown_name
        {
          std::cout<<"Warning ! No exact value function specified ! Returning 0." << std::endl;
        
          return 0;
        }
        
        Gradient exact_grad(const Point& /* p */,       // xyz location
                            const Parameters& /* param */,  // EquationSystems parameters
                            const std::string&, // sys_name
                            const std::string&) // unknown_name
        {
          std::cout<<"Warning ! No exact value function specified ! Returning 0." << std::endl;
        
          return 0;
        }



The source file output.h with comments:

        #include <string>
        
Bring in everything from the libMesh namespace
        using namespace libMesh;
        
        void start_output(unsigned int timesteps,
                          std::string filename,
                          std::string varname);



The source file adjoints_ex3.C with comments:

        #include <iostream>
        #include <sys/time.h>
        #include <iomanip>
        
Libmesh includes
        #include "libmesh/equation_systems.h"
        
        #include "libmesh/twostep_time_solver.h"
        #include "libmesh/euler_solver.h"
        #include "libmesh/euler2_solver.h"
        #include "libmesh/steady_solver.h"
        
        #include "libmesh/newton_solver.h"
        #include "libmesh/numeric_vector.h"
        #include "libmesh/petsc_diff_solver.h"
        
        #include "libmesh/mesh.h"
        #include "libmesh/mesh_tools.h"
        #include "libmesh/mesh_base.h"
        #include "libmesh/mesh_refinement.h"
        
        #include "libmesh/system.h"
        #include "libmesh/system_norm.h"
        
        #include "libmesh/adjoint_residual_error_estimator.h"
        #include "libmesh/const_fem_function.h"
        #include "libmesh/error_vector.h"
        #include "libmesh/fem_function_base.h"
        #include "libmesh/getpot.h"
        #include "libmesh/gmv_io.h"
        #include "libmesh/kelly_error_estimator.h"
        #include "libmesh/parameter_vector.h"
        #include "libmesh/patch_recovery_error_estimator.h"
        #include "libmesh/petsc_vector.h"
        #include "libmesh/sensitivity_data.h"
        #include "libmesh/tecplot_io.h"
        #include "libmesh/uniform_refinement_estimator.h"
        #include "libmesh/qoi_set.h"
        #include "libmesh/weighted_patch_recovery_error_estimator.h"
        
Local includes
        #include "coupled_system.h"
        #include "domain.h"
        #include "initial.h"
        #include "femparameters.h"
        #include "mysystems.h"
        #include "output.h"
        #include "H-qoi.h"
        
We solve a coupled Stokes + Convection Diffusion system in an H channel geometry with 2 inlets and 2 outlets. The QoI is the species flux from left outlet

Channel Geometry: Wall ---------------------------------------------------------------- 0 1 ---------------------------- ------------------------------- | | Wall | | Wall | | ----------------------------- ------------------------------- 2 2 ----------------------------------------------------------------- Wall

The equations governing this flow are: Stokes: -VectorLaplacian(velocity) + grad(pressure) = vector(0) Convection-Diffusion: - dot(velocity, grad(concentration) ) + Laplacian(concentration) = 0

The boundary conditions are: u_1(0) = -(y-2)*(y-3), u_2(0) = 0 ; u_1(1) = (y-2)*(y-3), u_2(1) = 0 ; u_1(walls) = 0, u_2(walls) = 0; C(0) = 1 ; C(1) = 0; grad(C) dot n (walls) = 0 ; grad(C) dot n (2) = 0 ; grad(C) dot n (3) = 0

The QoI is: Q((u,v), C) = integral_{left_outlet} - u * C ds

The complete equal order adjoint QoI error estimate is: (Notation for derivatives: grad(C) = C,1 + C,2) Q(u) - Q(u_h) \leq |e(u_1)|_{H1} |e(u_1^*)|_{H1} + |e(u_2)|_{H1} |e(u_2^*)|_{H1} + (1/Pe) * |e(C)|_{H1} |e(C^*)|_{H1} + ||e(u_1,1)||_{L2} ||e(p^*)||_{L2} + ||e(u_2,2)||_{L2} ||e(p^*)||_{L2} + ||e(u_1,1^*)||_{L2} ||e(p)||_{L2} + ||e(u_2,2^*)||_{L2} ||e(p)||_{L2} + ||e((u_1)_h C,1)||_{L2} ||e(C^*)||_{L2} + ||e(u1 C,1_h)||_{L2} ||e(C^*)||_{L2} ||e((u_2)_h C,2)||_{L2} ||e(C^*)||_{L2} + ||e(u2 C,2_h)||_{L2} ||e(C^*)||_{L2} = error_non_pressure + error_with_pressure + error_convection_diffusion_x + error_convection_diffusion_y

Bring in everything from the libMesh namespace
        using namespace libMesh;
        
Number output files

        std::string numbered_filename(unsigned int t_step, // The timestep count
                                      unsigned int a_step, // The adaptive step count
                                      std::string solution_type, // primal or adjoint solve
                                      std::string type,
                                      std::string extension,
                                      FEMParameters ¶m)
        {
          std::ostringstream file_name;
          file_name << solution_type
                    << ".out."
                    << type
                    << '.'
                    << std::setw(3)
                    << std::setfill('0')
                    << std::right
                    << t_step
                    << '.'
                    << std::setw(2)
                    << std::setfill('0')
                    << std::right
                    << a_step
                    << '.'
                    << extension;
        
          if (param.output_bz2)
            file_name << ".bz2";
          else if (param.output_gz)
            file_name << ".gz";
          return file_name.str();
        }
        
Write tecplot, gmv and xda/xdr output

        void write_output(EquationSystems &es,
                          unsigned int t_step, // The timestep count
                          unsigned int a_step, // The adaptive step count
                          std::string solution_type, // primal or adjoint solve
                          FEMParameters ¶m)
        {
          MeshBase &mesh = es.get_mesh();
        
        #ifdef LIBMESH_HAVE_GMV
          if (param.output_gmv)
            {
              std::ostringstream file_name_gmv;
              file_name_gmv << solution_type
                            << ".out.gmv."
                            << std::setw(3)
                            << std::setfill('0')
                            << std::right
                            << t_step
                            << '.'
                            << std::setw(2)
                            << std::setfill('0')
                            << std::right
                            << a_step;
        
              GMVIO(mesh).write_equation_systems
                (file_name_gmv.str(), es);
            }
        #endif
        
        #ifdef LIBMESH_HAVE_TECPLOT_API
          if (param.output_tecplot)
            {
              std::ostringstream file_name_tecplot;
              file_name_tecplot << solution_type
                                << ".out."
                                << std::setw(3)
                                << std::setfill('0')
                                << std::right
                                << t_step
                                << '.'
                                << std::setw(2)
                                << std::setfill('0')
                                << std::right
                                << a_step
                                << ".plt";
        
              TecplotIO(mesh).write_equation_systems
                (file_name_tecplot.str(), es);
            }
        #endif
        
          if (param.output_xda || param.output_xdr)
            mesh.renumber_nodes_and_elements();
          if (param.output_xda)
            {
              mesh.write(numbered_filename(t_step, a_step, solution_type, "mesh", "xda", param));
              es.write(numbered_filename(t_step, a_step, solution_type, "soln", "xda", param),
                       libMeshEnums::WRITE, EquationSystems::WRITE_DATA |
                       EquationSystems::WRITE_ADDITIONAL_DATA);
            }
          if (param.output_xdr)
            {
              mesh.write(numbered_filename(t_step, a_step, solution_type, "mesh", "xdr", param));
              es.write(numbered_filename(t_step, a_step, solution_type, "soln", "xdr", param),
                       libMeshEnums::WRITE, EquationSystems::WRITE_DATA |
                       EquationSystems::WRITE_ADDITIONAL_DATA);
            }
        }
        
        void write_output_headers(FEMParameters ¶m)
        {
Only one processor needs to take care of headers.
          if (libMesh::processor_id() != 0)
            return;
        
          start_output(param.initial_timestep, "out_clocktime.m", "vector_clocktime");
        
          if (param.run_simulation)
            {
              start_output(param.initial_timestep, "out_activemesh.m", "table_activemesh");
              start_output(param.initial_timestep, "out_solvesteps.m", "table_solvesteps");
        
              if (param.timesolver_tolerance)
                {
                  start_output(param.initial_timestep, "out_time.m", "vector_time");
                  start_output(param.initial_timestep, "out_timesteps.m", "vector_timesteps");
                }
              if (param.steadystate_tolerance)
                start_output(param.initial_timestep, "out_changerate.m", "vector_changerate");
            }
        }
        
        void write_output_solvedata(EquationSystems &es,
                                    unsigned int a_step,
                                    unsigned int newton_steps,
                                    unsigned int krylov_steps,
                                    unsigned int tv_sec,
                                    unsigned int tv_usec)
        {
          MeshBase &mesh = es.get_mesh();
          unsigned int n_active_elem = mesh.n_active_elem();
          unsigned int n_active_dofs = es.n_active_dofs();
        
          if (libMesh::processor_id() == 0)
            {
Write out the number of elements/dofs used
              std::ofstream activemesh ("out_activemesh.m",
                std::ios_base::app | std::ios_base::out);
              activemesh.precision(17);
              activemesh << (a_step + 1) << ' '
                         << n_active_elem << ' '
                         << n_active_dofs << std::endl;
        
Write out the number of solver steps used
              std::ofstream solvesteps ("out_solvesteps.m",
                std::ios_base::app | std::ios_base::out);
              solvesteps.precision(17);
              solvesteps << newton_steps << ' '
                         << krylov_steps << std::endl;
        
Write out the clock time
              std::ofstream clocktime ("out_clocktime.m",
                std::ios_base::app | std::ios_base::out);
              clocktime.precision(17);
              clocktime << tv_sec << '.' << tv_usec << std::endl;
            }
        }
        
        void write_output_footers(FEMParameters ¶m)
        {
Write footers on output .m files
          if (libMesh::processor_id() == 0)
            {
              std::ofstream clocktime ("out_clocktime.m",
                std::ios_base::app | std::ios_base::out);
              clocktime << "];" << std::endl;
        
              if (param.run_simulation)
                {
                  std::ofstream activemesh ("out_activemesh.m",
                    std::ios_base::app | std::ios_base::out);
                  activemesh << "];" << std::endl;
        
                  std::ofstream solvesteps ("out_solvesteps.m",
                    std::ios_base::app | std::ios_base::out);
                  solvesteps << "];" << std::endl;
        
                  if (param.timesolver_tolerance)
                    {
                      std::ofstream times ("out_time.m",
                        std::ios_base::app | std::ios_base::out);
                      times << "];" << std::endl;
                      std::ofstream timesteps ("out_timesteps.m",
                        std::ios_base::app | std::ios_base::out);
                      timesteps << "];" << std::endl;
                    }
                  if (param.steadystate_tolerance)
                    {
                      std::ofstream changerate ("out_changerate.m",
                        std::ios_base::app | std::ios_base::out);
                      changerate << "];" << std::endl;
                    }
                }
        
We're done, so write out a file (for e.g. Make to check)
              std::ofstream complete ("complete");
              complete << "complete" << std::endl;
            }
        }
        
        #if defined(LIBMESH_HAVE_GMV) || defined(LIBMESH_HAVE_TECPLOT_API)
        void write_error(EquationSystems &es,
                         ErrorVector &error,
                         unsigned int t_number,
                         unsigned int a_number,
                         FEMParameters ¶m,
                         std::string error_type)
        #else
        void write_error(EquationSystems &,
                         ErrorVector &,
                         unsigned int,
                         unsigned int,
                         FEMParameters &,
                         std::string)
        #endif
        {
        #ifdef LIBMESH_HAVE_GMV
          if (param.write_gmv_error)
            {
              std::ostringstream error_gmv;
              error_gmv << "error.gmv."
                        << std::setw(3)
                        << std::setfill('0')
                        << std::right
                        << a_number
                        << "."
                        << std::setw(2)
                        << std::setfill('0')
                        << std::right
                        << t_number
                        << error_type;
        
              error.plot_error(error_gmv.str(), es.get_mesh());
            }
        #endif
        
        #ifdef LIBMESH_HAVE_TECPLOT_API
          if (param.write_tecplot_error)
            {
              std::ostringstream error_tecplot;
              error_tecplot << "error.plt."
                            << std::setw(3)
                            << std::setfill('0')
                            << std::right
                            << a_number
                            << "."
                            << std::setw(2)
                            << std::setfill('0')
                            << std::right
                            << t_number
                            << error_type;
        
              error.plot_error(error_tecplot.str(), es.get_mesh());
            }
        #endif
        }
        
        void read_output(EquationSystems &es,
                         unsigned int t_step,
                         unsigned int a_step,
                         std::string solution_type,
                         FEMParameters ¶m)
        {
          MeshBase &mesh = es.get_mesh();
        
          std::string file_name_mesh, file_name_soln;
Look for ASCII files first
          if (param.output_xda)
            {
              file_name_mesh = numbered_filename(t_step, a_step, solution_type, "mesh", "xda", param);
              file_name_soln = numbered_filename(t_step, a_step, solution_type, "soln", "xda", param);
            }
          else if (param.output_xdr)
            {
              file_name_mesh = numbered_filename(t_step, a_step, solution_type, "mesh", "xdr", param);
              file_name_soln = numbered_filename(t_step, a_step, solution_type, "soln", "xdr", param);
            }
        
Read in the mesh
          mesh.read(file_name_mesh);
        
And the stored solution
          es.read(file_name_soln, libMeshEnums::READ,
                  EquationSystems::READ_HEADER |
                  EquationSystems::READ_DATA |
                  EquationSystems::READ_ADDITIONAL_DATA);
        
Put systems in a consistent state
          for (unsigned int i = 0; i != es.n_systems(); ++i)
            es.get_system<FEMSystem>(i).update();
        
Figure out the current time
          Real current_time = 0., current_timestep = 0.;
        
          if (param.timesolver_tolerance)
            {
              std::ifstream times ("out_time.m");
              std::ifstream timesteps ("out_timesteps.m");
              if (times.is_open() && timesteps.is_open())
                {
Read headers
                  const unsigned int headersize = 25;
                  char header[headersize];
                  timesteps.getline (header, headersize);
                  if (strcmp(header, "vector_timesteps = [") != 0)
                    {
                      std::cout << "Bad header in out_timesteps.m:" << std::endl
                                << header
                                << std::endl;
                      libmesh_error();
                    }
        
                  times.getline (header, headersize);
                  if (strcmp(header, "vector_time = [") != 0)
                    {
                      std::cout << "Bad header in out_time.m:" << std::endl
                                << header
                                << std::endl;
                      libmesh_error();
                    }
        
Read each timestep
                  for (unsigned int i = 0; i != t_step; ++i)
                    {
                      if (!times.good())
                        libmesh_error();
                      times >> current_time;
                      timesteps >> current_timestep;
                    }
Remember to increment the last timestep; out_times.m lists each *start* time
                  current_time += current_timestep;
                }
              else
                libmesh_error();
            }
          else
            current_time = t_step * param.deltat;
        
          for (unsigned int i = 0; i != es.n_systems(); ++i)
            es.get_system<FEMSystem>(i).time = current_time;
        }
        
        void set_system_parameters(FEMSystem &system, FEMParameters ¶m)
        {
Verify analytic jacobians against numerical ones?
          system.verify_analytic_jacobians = param.verify_analytic_jacobians;
        
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;
        
Solve this as a time-dependent or steady system
          if (param.transient)
            {
              UnsteadySolver *innersolver;
              if (param.timesolver_core == "euler2")
                {
                  Euler2Solver *euler2solver =
                    new Euler2Solver(system);
        
                  euler2solver->theta = param.timesolver_theta;
                  innersolver = euler2solver;
                }
              else if (param.timesolver_core == "euler")
                {
                  EulerSolver *eulersolver =
                    new EulerSolver(system);
        
                  eulersolver->theta = param.timesolver_theta;
                  innersolver = eulersolver;
                }
              else
                {
                  std::cerr << "Don't recognize core TimeSolver type: "
                            << param.timesolver_core << std::endl;
                  libmesh_error();
                }
        
              if (param.timesolver_tolerance)
                {
                  TwostepTimeSolver *timesolver =
                    new TwostepTimeSolver(system);
        
                  timesolver->max_growth       = param.timesolver_maxgrowth;
                  timesolver->target_tolerance = param.timesolver_tolerance;
                  timesolver->upper_tolerance  = param.timesolver_upper_tolerance;
                  timesolver->component_norm   = SystemNorm(param.timesolver_norm);
        
                  timesolver->core_time_solver =
                    AutoPtr<UnsteadySolver>(innersolver);
                  system.time_solver =
                    AutoPtr<UnsteadySolver>(timesolver);
                }
              else
                system.time_solver =
                  AutoPtr<TimeSolver>(innersolver);
            }
          else
            system.time_solver =
              AutoPtr<TimeSolver>(new SteadySolver(system));
        
          system.time_solver->reduce_deltat_on_diffsolver_failure =
                                                param.deltat_reductions;
          system.time_solver->quiet           = param.time_solver_quiet;
        
Set the time stepping options
          system.deltat = param.deltat;
        
And the integration options
          system.extra_quadrature_order = param.extra_quadrature_order;
        
And the nonlinear solver options
          if (param.use_petsc_snes)
            {
        #ifdef LIBMESH_HAVE_PETSC
              PetscDiffSolver *solver = new PetscDiffSolver(system);
              system.time_solver->diff_solver() = AutoPtr<DiffSolver>(solver);
        #else
              libmesh_error();
        #endif
            }
          else
            {
              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;
            }
        }
        
        #ifdef LIBMESH_ENABLE_AMR
        
        AutoPtr<MeshRefinement> build_mesh_refinement(MeshBase &mesh,
                                                      FEMParameters ¶m)
        {
          AutoPtr<MeshRefinement> mesh_refinement(new MeshRefinement(mesh));
          mesh_refinement->coarsen_by_parents() = true;
          mesh_refinement->absolute_global_tolerance() = param.global_tolerance;
          mesh_refinement->nelem_target()      = param.nelem_target;
          mesh_refinement->refine_fraction()   = param.refine_fraction;
          mesh_refinement->coarsen_fraction()  = param.coarsen_fraction;
          mesh_refinement->coarsen_threshold() = param.coarsen_threshold;
        
          return mesh_refinement;
        }
        
        #endif // LIBMESH_ENABLE_AMR
        
This function builds the Kelly error indicator. This indicator can be used for comparisons of adjoint and non-adjoint based error indicators
        AutoPtr<ErrorEstimator> build_error_estimator(FEMParameters& /* param */)
        {
          AutoPtr<ErrorEstimator> error_estimator;
        
          error_estimator.reset(new KellyErrorEstimator);
        
          return error_estimator;
        }
        
Functions to build the adjoint based error indicators The error_non_pressure and error_pressure constributions are estimated using the build_error_estimator_component_wise function below
        AutoPtr<ErrorEstimator>
        build_error_estimator_component_wise
          (FEMParameters ¶m,
           std::vector<std::vector<Real> > &term_weights,
           std::vector<libMeshEnums::FEMNormType> &primal_error_norm_type,
           std::vector<libMeshEnums::FEMNormType> &dual_error_norm_type)
        {
          AutoPtr<ErrorEstimator> error_estimator;
        
          AdjointResidualErrorEstimator *adjoint_residual_estimator = new AdjointResidualErrorEstimator;
        
          error_estimator.reset (adjoint_residual_estimator);
        
Both the primal and dual weights are going to be estimated using the patch recovery error estimator
          PatchRecoveryErrorEstimator *p1 =
            new PatchRecoveryErrorEstimator;
          adjoint_residual_estimator->primal_error_estimator().reset(p1);
        
          PatchRecoveryErrorEstimator *p2 =
            new PatchRecoveryErrorEstimator;
          adjoint_residual_estimator->dual_error_estimator().reset(p2);
        
Set the boolean for specifying whether we are reusing patches while building the patch recovery estimates
          p1->set_patch_reuse(param.patch_reuse);
          p2->set_patch_reuse(param.patch_reuse);
        
Using the user filled error norm type vector, we pass the type of norm to be used for the error in each variable, we can have different types of norms for the primal and dual variables
          unsigned int size = primal_error_norm_type.size();
        
          libmesh_assert_equal_to (size, dual_error_norm_type.size());
          for (unsigned int i = 0; i != size; ++i)
            {
              adjoint_residual_estimator->primal_error_estimator()->error_norm.set_type(i, primal_error_norm_type[i]);
              adjoint_residual_estimator->dual_error_estimator()->error_norm.set_type(i, dual_error_norm_type[i]);
            }
        
Now we set the right weights for each term in the error estimate, using the user provided term_weights matrix
          libmesh_assert_equal_to (size, term_weights.size());
          for (unsigned int i = 0; i != term_weights.size(); ++i)
            {
              libmesh_assert_equal_to (size, term_weights[i].size());
              adjoint_residual_estimator->error_norm.set_weight(i, term_weights[i][i]);
              for (unsigned int j = 0; j != size; ++j)
                if (i != j)
                  adjoint_residual_estimator->error_norm.set_off_diagonal_weight(i, j, term_weights[i][j]);
            }
        
          return error_estimator;
        }
        
The error_convection_diffusion_x and error_convection_diffusion_y are the nonlinear contributions which are computed using the build_weighted_error_estimator_component_wise below
        AutoPtr<ErrorEstimator>
        build_weighted_error_estimator_component_wise
          (FEMParameters ¶m,
           std::vector<std::vector<Real> > &term_weights,
           std::vector<libMeshEnums::FEMNormType> &primal_error_norm_type,
           std::vector<libMeshEnums::FEMNormType> &dual_error_norm_type,
           std::vector<FEMFunctionBase<Number>*> coupled_system_weight_functions)
        {
          AutoPtr<ErrorEstimator> error_estimator;
        
          AdjointResidualErrorEstimator *adjoint_residual_estimator = new AdjointResidualErrorEstimator;
        
          error_estimator.reset (adjoint_residual_estimator);
        
Using the user filled error norm type vector, we pass the type of norm to be used for the error in each variable, we can have different types of norms for the primal and dual variables

          WeightedPatchRecoveryErrorEstimator *p1 =
            new WeightedPatchRecoveryErrorEstimator;
          adjoint_residual_estimator->primal_error_estimator().reset(p1);
        
          PatchRecoveryErrorEstimator *p2 =
            new PatchRecoveryErrorEstimator;
          adjoint_residual_estimator->dual_error_estimator().reset(p2);
        
          p1->set_patch_reuse(param.patch_reuse);
          p2->set_patch_reuse(param.patch_reuse);
        
This is the critical difference with the build_error_estimate_component_wise
          p1->weight_functions.clear();
        
We pass the pointers to the user specified weight functions to the patch recovery error estimator objects declared above
          unsigned int size = primal_error_norm_type.size();
        
          libmesh_assert(coupled_system_weight_functions.size() == size);
          libmesh_assert(dual_error_norm_type.size() == size);
        
          for (unsigned int i = 0; i != size; ++i)
            {
              p1->weight_functions.push_back(coupled_system_weight_functions[i]);
              adjoint_residual_estimator->primal_error_estimator()->error_norm.set_type(i, primal_error_norm_type[i]);
              adjoint_residual_estimator->dual_error_estimator()->error_norm.set_type(i, dual_error_norm_type[i]);
            }
        
Now we set the right weights for each term in the error estimate, using the user provided term_weights matrix
          libmesh_assert_equal_to (size, term_weights.size());
          for (unsigned int i = 0; i != term_weights.size(); ++i)
            {
              libmesh_assert_equal_to (size, term_weights[i].size());
              adjoint_residual_estimator->error_norm.set_weight(i, term_weights[i][i]);
              for (unsigned int j = 0; j != size; ++j)
                if (i != j)
                  adjoint_residual_estimator->error_norm.set_off_diagonal_weight(i, j, term_weights[i][j]);
            }
        
          return error_estimator;
        }
        
The main program.
        int main (int argc, char** argv)
        {
Initialize libMesh.
          LibMeshInit init (argc, argv);
        
Skip adaptive examples on a non-adaptive libMesh build
        #ifndef LIBMESH_ENABLE_AMR
          libmesh_example_assert(false, "--enable-amr");
        #else
        
This doesn't converge with Eigen BICGSTAB for some reason...
          libmesh_example_assert(libMesh::default_solver_package() != EIGEN_SOLVERS, "--enable-petsc");
        
          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 higher-dimensional examples on a lower-dimensional libMesh build
          libmesh_example_assert(2 <= LIBMESH_DIM, "2D support");
        
Create a mesh with the given dimension, distributed across the default MPI communicator.
          Mesh mesh(init.comm(), param.dimension);
        
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 << "Building mesh" << std::endl;
        
          if (!param.initial_timestep && param.run_simulation)
            build_domain(mesh, param);
        
          std::cout << "Building system" << std::endl;
        
          FEMSystem &system = build_system(equation_systems, infile, param);
        
          set_system_parameters(system, param);
        
          std::cout << "Initializing systems" << std::endl;
        
Create a QoI object, we will need this to compute the QoI
          {
            CoupledSystemQoI qoi;
Our QoI is computed on the side
            qoi.assemble_qoi_sides = true;
            system.attach_qoi(&qoi);
          }
        
          equation_systems.init ();
        
Print information about the mesh and system to the screen.
          mesh.print_info();
          equation_systems.print_info();
        
          std::cout<<"Starting adaptive loop"<<std::endl<<std::endl;
        
Count the number of steps used
          unsigned int a_step = 0;
        
Get the linear solver object to set the preconditioner reuse flag
          LinearSolver<Number> *linear_solver = system.get_linear_solver();
        
Adaptively solve the timestep
          for (; a_step <= param.max_adaptivesteps; ++a_step)
            {
We have to ensure that we are refining to either an error tolerance or to a target number of elements, and, not to a non-zero tolerance and a target number of elements simultaneously
              if(param.global_tolerance > 0 && param.nelem_target > 0)
                {
                  std::cout<<"We cant refine to both a non-zero tolerance and a target number of elements, EXITING adaptive loop. "<<std::endl<<std::endl;
                  break;
                }
        
We dont want to reuse the preconditioner for the forward solve
              linear_solver->reuse_preconditioner(false);
        
              std::cout<< "Adaptive step " << a_step << std::endl;
        
              std::cout<< "Solving the forward problem" <<std::endl;
        
              std::cout << "We have " << mesh.n_active_elem()
                            << " active elements and " << equation_systems.n_active_dofs()
                        << " active dofs." << std::endl << std::endl;
        
Solve the forward system
              system.solve();
        
Write the output
              write_output(equation_systems, 0, a_step, "primal", param);
        
Compute the QoI
              system.assemble_qoi();
        
We just call a postprocess here to set the variable computed_QoI to the value computed by assemble_qoi
              system.postprocess();
        
Get the value of the computed_QoI variable of the CoupledSystem class
              Number QoI_0_computed = (dynamic_cast<CoupledSystem&>(system)).get_QoI_value();
        
              std::cout<< "The boundary QoI is " << std::setprecision(17) << QoI_0_computed << std::endl << std::endl;
        
You dont need to compute error estimates and refine at the last adaptive step, only before that
              if(a_step!=param.max_adaptivesteps)
                {
                  NumericVector<Number> &primal_solution = (*system.solution);
        
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
                  std::cout<< "Solving the adjoint problem" <<std::endl;
                  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);
        
To plot the adjoint solution, we swap it with the primal solution and use the write_output function
                  NumericVector<Number> &dual_solution = system.get_adjoint_solution();
                  primal_solution.swap(dual_solution);
        
                  write_output(equation_systems, 0, a_step, "adjoint", param);
        
Swap back
                  primal_solution.swap(dual_solution);
        
We need the values of the parameters Pe from the system for the adjoint error estimate
                  Real Pe = (dynamic_cast<CoupledSystem&>(system)).get_Pe();
        
The total error is the sum: error = error_non_pressure + error_with_pressure + ... error_estimator_convection_diffusion_x + error_estimator_convection_diffusion_y
                  ErrorVector error;
        
We first construct the non-pressure contributions
                  ErrorVector error_non_pressure;
        
First we build the norm_type_vector_non_pressure vectors and weights_matrix_non_pressure matrix for the non-pressure term error contributions
                  std::vector<libMeshEnums::FEMNormType>
                    primal_norm_type_vector_non_pressure;
                  primal_norm_type_vector_non_pressure.push_back(H1_SEMINORM);
                  primal_norm_type_vector_non_pressure.push_back(H1_SEMINORM);
                  primal_norm_type_vector_non_pressure.push_back(L2);
                  primal_norm_type_vector_non_pressure.push_back(H1_SEMINORM);
        
                  std::vector<libMeshEnums::FEMNormType>
                    dual_norm_type_vector_non_pressure;
                  dual_norm_type_vector_non_pressure.push_back(H1_SEMINORM);
                  dual_norm_type_vector_non_pressure.push_back(H1_SEMINORM);
                  dual_norm_type_vector_non_pressure.push_back(L2);
                  dual_norm_type_vector_non_pressure.push_back(H1_SEMINORM);
        
                  std::vector<std::vector<Real> >
                    weights_matrix_non_pressure(system.n_vars(),
                      std::vector<Real>(system.n_vars(), 0.0));
                  weights_matrix_non_pressure[0][0] = 1.;
                  weights_matrix_non_pressure[1][1] = 1.;
                  weights_matrix_non_pressure[3][3] = 1./Pe;
        
We build the error estimator to estimate the contributions to the QoI error from the non pressure term
                  AutoPtr<ErrorEstimator> error_estimator_non_pressure =
                    build_error_estimator_component_wise
                      (param, weights_matrix_non_pressure,
                       primal_norm_type_vector_non_pressure,
                       dual_norm_type_vector_non_pressure);
        
Estimate the contributions to the QoI error from the non pressure terms
                  error_estimator_non_pressure->estimate_error(system, error_non_pressure);
        
Plot the estimated error from the non_pressure terms
                  write_error(equation_systems, error_non_pressure, 0, a_step, param, "_non_pressure");
        
Now for the pressure contributions
                  ErrorVector error_with_pressure;
        
Next we build the norm_type_vector_with_pressure vectors and weights_matrix_with_pressure matrix for the pressure term error contributions
                  std::vector<libMeshEnums::FEMNormType>
                    primal_norm_type_vector_with_pressure;
                  primal_norm_type_vector_with_pressure.push_back(H1_X_SEMINORM);
                  primal_norm_type_vector_with_pressure.push_back(H1_Y_SEMINORM);
                  primal_norm_type_vector_with_pressure.push_back(L2);
                  primal_norm_type_vector_with_pressure.push_back(L2);
        
                  std::vector<libMeshEnums::FEMNormType>
                    dual_norm_type_vector_with_pressure;
                  dual_norm_type_vector_with_pressure.push_back(H1_X_SEMINORM);
                  dual_norm_type_vector_with_pressure.push_back(H1_Y_SEMINORM);
                  dual_norm_type_vector_with_pressure.push_back(L2);
                  dual_norm_type_vector_with_pressure.push_back(L2);
        
                  std::vector<std::vector<Real> >
                    weights_matrix_with_pressure
                      (system.n_vars(),
                       std::vector<Real>(system.n_vars(), 0.0));
                  weights_matrix_with_pressure[0][2] = 1.;
        
                  weights_matrix_with_pressure[1][2] = 1.;
        
                  weights_matrix_with_pressure[2][0] = 1.;
                  weights_matrix_with_pressure[2][1] = 1.;
        
We build the error estimator to estimate the contributions to the QoI error from the pressure term
                  AutoPtr<ErrorEstimator> error_estimator_with_pressure =
                    build_error_estimator_component_wise
                      (param, weights_matrix_with_pressure,
                       primal_norm_type_vector_with_pressure,
                       dual_norm_type_vector_with_pressure);
        
Estimate the contributions to the QoI error from the pressure terms
                  error_estimator_with_pressure->estimate_error(system, error_with_pressure);
        
Plot the error due to the pressure terms
                  write_error(equation_systems, error_with_pressure, 0, a_step, param, "_with_pressure");
        
Now for the convection diffusion term errors (in the x and y directions)

                  ErrorVector error_convection_diffusion_x;
        
The norm type vectors and weights matrix for the convection_diffusion_x errors
                  std::vector<libMeshEnums::FEMNormType>
                    primal_norm_type_vector_convection_diffusion_x;
                  primal_norm_type_vector_convection_diffusion_x.push_back(L2);
                  primal_norm_type_vector_convection_diffusion_x.push_back(L2);
                  primal_norm_type_vector_convection_diffusion_x.push_back(L2);
                  primal_norm_type_vector_convection_diffusion_x.push_back(H1_X_SEMINORM);
        
                  std::vector<libMeshEnums::FEMNormType>
                    dual_norm_type_vector_convection_diffusion_x;
                  dual_norm_type_vector_convection_diffusion_x.push_back(L2);
                  dual_norm_type_vector_convection_diffusion_x.push_back(L2);
                  dual_norm_type_vector_convection_diffusion_x.push_back(L2);
Note that we need the error of the dual concentration in L2
                  dual_norm_type_vector_convection_diffusion_x.push_back(L2);
        
                  std::vector<std::vector<Real> >
                    weights_matrix_convection_diffusion_x
                      (system.n_vars(),
                       std::vector<Real>(system.n_vars(), 0.0));
        	  weights_matrix_convection_diffusion_x[0][3] = 1.;
                  weights_matrix_convection_diffusion_x[3][3] = 1.;
        
We will also have to build and pass the weight functions to the weighted patch recovery estimators

We pass the identity function as weights to error entries that the above matrix will scale to 0.
                  ConstFEMFunction<Number> identity(1);
        
Declare object of class CoupledFEMFunctionsx, the definition of the function contains the weight to be applied to the relevant terms For ||e(u1 C,1_h)||_{L2} ||e(C^*)||_{L2} term, returns C,1_h
                  CoupledFEMFunctionsx convdiffx0(system, 0);
For ||e((u_1)_h C,1)||_{L2} ||e(C^*)||_{L2} term, returns (u_1)_h
                  CoupledFEMFunctionsx convdiffx3(system, 3);
        
Make a vector of pointers to these objects
                  std::vector<FEMFunctionBase<Number> *> coupled_system_weight_functions_x;
                  coupled_system_weight_functions_x.push_back(&convdiffx0);
                  coupled_system_weight_functions_x.push_back(&identity);
                  coupled_system_weight_functions_x.push_back(&identity);
                  coupled_system_weight_functions_x.push_back(&convdiffx3);
        
Build the error estimator to estimate the contributions to the QoI error from the convection diffusion x term
                  AutoPtr<ErrorEstimator> error_estimator_convection_diffusion_x =
                  build_weighted_error_estimator_component_wise
                    (param, weights_matrix_convection_diffusion_x,
                     primal_norm_type_vector_convection_diffusion_x,
                     dual_norm_type_vector_convection_diffusion_x,
                     coupled_system_weight_functions_x);
        
Estimate the contributions to the QoI error from the convection diffusion x term
                  error_estimator_convection_diffusion_x->estimate_error
                    (system, error_convection_diffusion_x);
        
Plot this error
                  write_error(equation_systems, error_convection_diffusion_x,
                              0, a_step, param, "_convection_diffusion_x");
        
Now for the y direction terms
                  ErrorVector error_convection_diffusion_y;
        
The norm type vectors and weights matrix for the convection_diffusion_x errors
                  std::vector<libMeshEnums::FEMNormType>
                    primal_norm_type_vector_convection_diffusion_y;
                  primal_norm_type_vector_convection_diffusion_y.push_back(L2);
                  primal_norm_type_vector_convection_diffusion_y.push_back(L2);
                  primal_norm_type_vector_convection_diffusion_y.push_back(L2);
                  primal_norm_type_vector_convection_diffusion_y.push_back(H1_Y_SEMINORM);
        
                  std::vector<libMeshEnums::FEMNormType>
                    dual_norm_type_vector_convection_diffusion_y;
                  dual_norm_type_vector_convection_diffusion_y.push_back(L2);
                  dual_norm_type_vector_convection_diffusion_y.push_back(L2);
                  dual_norm_type_vector_convection_diffusion_y.push_back(L2);
Note that we need the error of the dual concentration in L2
                  dual_norm_type_vector_convection_diffusion_y.push_back(L2);
        
                  std::vector<std::vector<Real> >
                    weights_matrix_convection_diffusion_y
                      (system.n_vars(), std::vector<Real>(system.n_vars(), 0.0));
        	  weights_matrix_convection_diffusion_y[1][3] = 1.;
                  weights_matrix_convection_diffusion_y[3][3] = 1.;
        
For ||e(u2 C,2_h)||_{L2} ||e(C^*)||_{L2} term, returns C,2_h
                  CoupledFEMFunctionsy convdiffy1(system, 1);
For ||e((u_2)_h C,2)||_{L2} ||e(C^*)||_{L2} term, returns (u_2)_h
                  CoupledFEMFunctionsy convdiffy3(system, 3);
        
Make a vector of pointers to these objects
                  std::vector<FEMFunctionBase<Number> *> coupled_system_weight_functions_y;
                  coupled_system_weight_functions_y.push_back(&identity);
                  coupled_system_weight_functions_y.push_back(&convdiffy1);
                  coupled_system_weight_functions_y.push_back(&identity);
                  coupled_system_weight_functions_y.push_back(&convdiffy3);
        
Build the error estimator to estimate the contributions to the QoI error from the convection diffsion y term
                  AutoPtr<ErrorEstimator> error_estimator_convection_diffusion_y =
                    build_weighted_error_estimator_component_wise
                      (param, weights_matrix_convection_diffusion_y,
                       primal_norm_type_vector_convection_diffusion_y,
                       dual_norm_type_vector_convection_diffusion_y,
                       coupled_system_weight_functions_y);
        
Estimate the contributions to the QoI error from the convection diffusion y terms
                  error_estimator_convection_diffusion_y->estimate_error(system, error_convection_diffusion_y);
        
Plot this error
                  write_error(equation_systems, error_convection_diffusion_y, 0, a_step, param, "_convection_diffusion_y");
        
                  if(param.indicator_type == "adjoint_residual")
                    {
                      error.resize(error_non_pressure.size());
        
Now combine the contribs from the pressure and non pressure terms to get the complete estimate
                      for(unsigned int i = 0; i < error.size(); i++)
                        {
                          error[i] = error_non_pressure[i] + error_with_pressure[i] + error_convection_diffusion_x[i] + error_convection_diffusion_y[i];
                        }
                    }
                  else
                    {
                      std::cout<<"Using Kelly Estimator"<<std::endl;
Build the Kelly error estimator
                      AutoPtr<ErrorEstimator> error_estimator =  build_error_estimator(param);
        
Estimate the error
                      error_estimator->estimate_error(system, error);
                    }
        
                  write_error(equation_systems, error, 0, a_step, param, "_total");
        
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

                  if(param.refine_uniformly)
                    {
                      mesh_refinement->uniformly_refine(1);
                    }
                  else if(param.global_tolerance >= 0. && param.nelem_target == 0.) // Refine based on reaching an error tolerance
                    {
                      mesh_refinement->flag_elements_by_error_tolerance (error);
        
Carry out the adaptive mesh refinement/coarsening
                      mesh_refinement->refine_and_coarsen_elements();
                    }
                  else // Refine based on reaching a target number of elements
                    {
If we have reached the desired number of elements, we dont do any more adaptive steps
                      if (mesh.n_active_elem() >= param.nelem_target)
                        {
                          std::cout<<"We reached the target number of elements."<<std::endl <<std::endl;
                          break;
                        }
        
                      mesh_refinement->flag_elements_by_nelem_target (error);
        
Carry out the adaptive mesh refinement/coarsening
                      mesh_refinement->refine_and_coarsen_elements();
                    }
        
                  equation_systems.reinit();
        
                }
        
            }
        
        
          write_output_footers(param);
        
        #endif // #ifndef LIBMESH_ENABLE_AMR
        
All done.
          return 0;
        
        }



The source file coupled_system.C with comments:

        #include "libmesh/boundary_info.h"
        #include "libmesh/dirichlet_boundaries.h"
        #include "libmesh/dof_map.h"
        #include "libmesh/fe_base.h"
        #include "libmesh/fe_interface.h"
        #include "libmesh/fem_context.h"
        #include "libmesh/getpot.h"
        #include "libmesh/mesh.h"
        #include "libmesh/parallel.h"
        #include "libmesh/quadrature.h"
        #include "libmesh/string_to_enum.h"
        #include "libmesh/zero_function.h"
        
        #include "coupled_system.h"
        
Bring in everything from the libMesh namespace
        using namespace libMesh;
        
Function to set the Dirichlet boundary function for the inlet boundary velocities
        class BdyFunction : public FunctionBase<Number>
        {
        public:
          BdyFunction (unsigned int u_var, unsigned int v_var, int sign)
            : _u_var(u_var), _v_var(v_var), _sign(sign)
            { this->_initialized = true; }
        
          virtual Number operator() (const Point&, const Real = 0)
            { libmesh_not_implemented(); }
        
          virtual void operator() (const Point& p,
                                   const Real,
                                   DenseVector<Number>& output)
            {
              output.resize(2);
              output.zero();
              const Real y=p(1);
Set the parabolic inflow boundary conditions at stations 0 & 1
              output(_u_var) = (_sign)*((y-2) * (y-3));
              output(_v_var) = 0;
            }
        
          virtual AutoPtr<FunctionBase<Number> > clone() const
            { return AutoPtr<FunctionBase<Number> > (new BdyFunction(_u_var, _v_var, _sign)); }
        
        private:
          const unsigned int _u_var, _v_var;
          const Real _sign;
        };
        
        
        void CoupledSystem::init_data ()
        {
Check the input file for Reynolds number, application type, approximation type
          GetPot infile("coupled_system.in");
          Peclet = infile("Peclet", 1.);
          unsigned int pressure_p = infile("pressure_p", 1);
          std::string fe_family = infile("fe_family", std::string("LAGRANGE"));
        
LBB needs better-than-quadratic velocities for better-than-linear pressures, and libMesh needs non-Lagrange elements for better-than-quadratic velocities.
          libmesh_assert((pressure_p == 1) || (fe_family != "LAGRANGE"));
        
          FEFamily fefamily = Utility::string_to_enum<FEFamily>(fe_family);
        
Add the velocity components "u" & "v". They will be approximated using second-order approximation.
          u_var = this->add_variable ("u", static_cast<Order>(pressure_p+1),
        			      fefamily);
          v_var = this->add_variable ("v", static_cast<Order>(pressure_p+1),
        			      fefamily);
        
Add the pressure variable "p". This will be approximated with a first-order basis, providing an LBB-stable pressure-velocity pair.
          p_var = this->add_variable ("p", static_cast<Order>(pressure_p),
        			      fefamily);
        
Add the Concentration variable "C". They will be approximated using second-order approximation, the same as the velocity components
          C_var = this->add_variable ("C", static_cast<Order>(pressure_p+1),
        			      fefamily);
        
Tell the system to march velocity forward in time, but leave p as a constraint only
          this->time_evolving(u_var);
          this->time_evolving(v_var);
          this->time_evolving(C_var);
        
Useful debugging options
          this->verify_analytic_jacobians = infile("verify_analytic_jacobians", 0.);
          this->print_jacobians = infile("print_jacobians", false);
          this->print_element_jacobians = infile("print_element_jacobians", false);
        
Set Dirichlet boundary conditions
          const boundary_id_type left_inlet_id = 0;
          std::set<boundary_id_type> left_inlet_bdy;
          left_inlet_bdy.insert(left_inlet_id);
        
          const boundary_id_type right_inlet_id = 1;
          std::set<boundary_id_type> right_inlet_bdy;
          right_inlet_bdy.insert(right_inlet_id);
        
          const boundary_id_type outlets_id = 2;
          std::set<boundary_id_type> outlets_bdy;
          outlets_bdy.insert(outlets_id);
        
          const boundary_id_type wall_id = 3;
          std::set<boundary_id_type> wall_bdy;
          wall_bdy.insert(wall_id);
        
The uv identifier for the setting the inlet and wall velocity boundary conditions
          std::vector<unsigned int> uv(1, u_var);
          uv.push_back(v_var);
The C_only identifier for setting the concentrations at the inlets
          std::vector<unsigned int> C_only(1, C_var);
        
The zero and constant functions
          ZeroFunction<Number> zero;
          ConstFunction<Number> one(1);
        
We need two boundary functions for the inlets, because the signs on the velocities will be different
          int velocity_sign = 1;
          BdyFunction inflow_left(u_var, v_var, -velocity_sign);
          BdyFunction inflow_right(u_var, v_var, velocity_sign);
        
On the walls we will apply the no slip and no penetration boundary condition, u=0, v=0
          this->get_dof_map().add_dirichlet_boundary
                (DirichletBoundary (wall_bdy, uv, &zero));
        
On the inlet (left), we apply parabolic inflow boundary conditions for the velocity, u = - (y-2)*(y-3), v=0 and set C = 1
          this->get_dof_map().add_dirichlet_boundary
            (DirichletBoundary (left_inlet_bdy, uv, &inflow_left));
          this->get_dof_map().add_dirichlet_boundary
            (DirichletBoundary (left_inlet_bdy, C_only, &one));
        
On the inlet (right), we apply parabolic inflow boundary conditions for the velocity, u = (y-2)*(y-3), v=0 and set C = 0
          this->get_dof_map().add_dirichlet_boundary
            (DirichletBoundary (right_inlet_bdy, uv, &inflow_right));
          this->get_dof_map().add_dirichlet_boundary
            (DirichletBoundary (right_inlet_bdy, C_only, &zero));
        
Note that the remaining boundary conditions are the natural boundary conditions for the concentration on the wall (grad(c) dot n = 0) and natural boundary conditions for the velocity and the concentration on the outlets ((grad(velocity) dot n - p n) dot t = 0, grad(C) dot n = 0)

Do the parent's initialization after variables and boundary constraints are defined
          FEMSystem::init_data();
        }
        
        void CoupledSystem::init_context(DiffContext &context)
        {
          FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
        
We should prerequest all the data we will need to build the linear system. Note that the concentration and velocity components use the same basis.
          c.element_fe_var[u_var]->get_JxW();
          c.element_fe_var[u_var]->get_phi();
          c.element_fe_var[u_var]->get_dphi();
          c.element_fe_var[u_var]->get_xyz();
        
          c.element_fe_var[p_var]->get_phi();
        
          c.side_fe_var[u_var]->get_JxW();
          c.side_fe_var[u_var]->get_phi();
          c.side_fe_var[u_var]->get_xyz();
        }
        
        
        bool CoupledSystem::element_time_derivative (bool request_jacobian,
                                                    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.element_fe_var[u_var]->get_JxW();
        
The velocity shape functions at interior quadrature points.
          const std::vector<std::vector<Real> >& phi =
            c.element_fe_var[u_var]->get_phi();
        
The velocity shape function gradients at interior quadrature points.
          const std::vector<std::vector<RealGradient> >& dphi =
            c.element_fe_var[u_var]->get_dphi();
        
The pressure shape functions at interior quadrature points.
          const std::vector<std::vector<Real> >& psi =
            c.element_fe_var[p_var]->get_phi();
        
The number of local degrees of freedom in each variable
          const unsigned int n_p_dofs = c.dof_indices_var[p_var].size();
          const unsigned int n_u_dofs = c.dof_indices_var[u_var].size();
          libmesh_assert_equal_to (n_u_dofs, c.dof_indices_var[v_var].size());
        
The subvectors and submatrices we need to fill:
          DenseSubMatrix<Number> &Kuu = *c.elem_subjacobians[u_var][u_var];
          DenseSubMatrix<Number> &Kup = *c.elem_subjacobians[u_var][p_var];
          DenseSubVector<Number> &Fu = *c.elem_subresiduals[u_var];
        
          DenseSubMatrix<Number> &Kvv = *c.elem_subjacobians[v_var][v_var];
          DenseSubMatrix<Number> &Kvp = *c.elem_subjacobians[v_var][p_var];
          DenseSubVector<Number> &Fv = *c.elem_subresiduals[v_var];
        
          DenseSubMatrix<Number> &KCu = *c.elem_subjacobians[C_var][u_var];
          DenseSubMatrix<Number> &KCv = *c.elem_subjacobians[C_var][v_var];
          DenseSubMatrix<Number> &KCC = *c.elem_subjacobians[C_var][C_var];
          DenseSubVector<Number> &FC = *c.elem_subresiduals[C_var];
        
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.element_qrule->n_points();
        
          for (unsigned int qp=0; qp != n_qpoints; qp++)
            {
Compute the solution & its gradient at the old Newton iterate
              Number p = c.interior_value(p_var, qp),
        	u = c.interior_value(u_var, qp),
        	v = c.interior_value(v_var, qp);
              Gradient grad_u = c.interior_gradient(u_var, qp),
        	grad_v = c.interior_gradient(v_var, qp),
        	grad_C = c.interior_gradient(C_var, qp);
        
Definitions for convenience. It is sometimes simpler to do a dot product if you have the full vector at your disposal.
              NumberVectorValue U     (u,     v);
              const Number C_x = grad_C(0);
              const Number C_y = grad_C(1);
        
First, an i-loop over the velocity degrees of freedom. We know that n_u_dofs == n_v_dofs so we can compute contributions for both at the same time.
              for (unsigned int i=0; i != n_u_dofs; i++)
                {
Stokes equations residuals
                  Fu(i) += JxW[qp] *
                           (p*dphi[i][qp](0) -                // pressure term
        		    (grad_u*dphi[i][qp]));            // diffusion term
        
                  Fv(i) += JxW[qp] *
                           (p*dphi[i][qp](1) -                // pressure term
        		    (grad_v*dphi[i][qp]));            // diffusion term
        
Concentration Equation Residual
                  FC(i) += JxW[qp] *
        	           ( (U*grad_C)*phi[i][qp] +                // convection term
        	            (1./Peclet)*(grad_C*dphi[i][qp]) );     // diffusion term
        
Note that the Fp block is identically zero unless we are using some kind of artificial compressibility scheme...

                  if (request_jacobian && c.elem_solution_derivative)
                    {
                      libmesh_assert (c.elem_solution_derivative == 1.0);
        
Matrix contributions for the uu and vv couplings.
                      for (unsigned int j=0; j != n_u_dofs; j++)
                        {
                          Kuu(i,j) += JxW[qp] * (-(dphi[i][qp]*dphi[j][qp])); /* diffusion term  */
        
                          Kvv(i,j) += JxW[qp] * (-(dphi[i][qp]*dphi[j][qp])); /* diffusion term  */
        
        		  KCu(i,j) += JxW[qp]* ( (phi[j][qp]*C_x)*phi[i][qp] ); /* convection term */
        
        		  KCv(i,j) += JxW[qp]*( (phi[j][qp]*C_y)*phi[i][qp] );  /* convection term */
        
        		  KCC(i,j) += JxW[qp]*
        		              ( (U*dphi[j][qp])*phi[i][qp] +      /* nonlinear term (convection) */
        		              (1./Peclet)*(dphi[j][qp]*dphi[i][qp]) ); /* diffusion term */
        		}
        
Matrix contributions for the up and vp couplings.
                      for (unsigned int j=0; j != n_p_dofs; j++)
        		{
        		  Kup(i,j) += JxW[qp]*psi[j][qp]*dphi[i][qp](0);
        		  Kvp(i,j) += JxW[qp]*psi[j][qp]*dphi[i][qp](1);
        		}
        	    }
        	}
        
            } // end of the quadrature point qp-loop
        
              return request_jacobian;
        }
        
        
        bool CoupledSystem::element_constraint (bool request_jacobian,
                                               DiffContext &context)
        {
          FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
        
Here we define some references to cell-specific data that will be used to assemble the linear system.

Element Jacobian * quadrature weight for interior integration
          const std::vector<Real> &JxW = c.element_fe_var[u_var]->get_JxW();
        
The velocity shape function gradients at interior quadrature points.
          const std::vector<std::vector<RealGradient> >& dphi =
            c.element_fe_var[u_var]->get_dphi();
        
The pressure shape functions at interior quadrature points.
          const std::vector<std::vector<Real> >& psi =
            c.element_fe_var[p_var]->get_phi();
        
The number of local degrees of freedom in each variable
          const unsigned int n_u_dofs = c.dof_indices_var[u_var].size();
          const unsigned int n_p_dofs = c.dof_indices_var[p_var].size();
        
The subvectors and submatrices we need to fill:
          DenseSubMatrix<Number> &Kpu = *c.elem_subjacobians[p_var][u_var];
          DenseSubMatrix<Number> &Kpv = *c.elem_subjacobians[p_var][v_var];
          DenseSubVector<Number> &Fp = *c.elem_subresiduals[p_var];
        
Add the constraint given by the continuity equation
          unsigned int n_qpoints = c.element_qrule->n_points();
          for (unsigned int qp=0; qp != n_qpoints; qp++)
            {
Compute the velocity gradient at the old Newton iterate
              Gradient grad_u = c.interior_gradient(u_var, qp),
        	grad_v = c.interior_gradient(v_var, qp);
        
Now a loop over the pressure degrees of freedom. This computes the contributions of the continuity equation.
              for (unsigned int i=0; i != n_p_dofs; i++)
                {
                  Fp(i) += JxW[qp] * psi[i][qp] *
                           (grad_u(0) + grad_v(1));
        
                  if (request_jacobian && c.elem_solution_derivative)
                    {
                      libmesh_assert (c.elem_solution_derivative == 1.0);
        
                      for (unsigned int j=0; j != n_u_dofs; j++)
                        {
                          Kpu(i,j) += JxW[qp]*psi[i][qp]*dphi[j][qp](0);
                          Kpv(i,j) += JxW[qp]*psi[i][qp]*dphi[j][qp](1);
                        }
                    }
                }
            } // end of the quadrature point qp-loop
        
          return request_jacobian;
        }
        
        void CoupledSystem::postprocess()
        {
We need to overload the postprocess function to set the computed_QoI variable of the CoupledSystem class to the qoi value stored in System::qoi[0]

          computed_QoI = 0.0;
        
          computed_QoI = System::qoi[0];
        }
        
These functions supply the nonlinear weighting for the adjoint residual error estimate which arise due to the convection term in the convection-diffusion equation: ||e((u_1)_h C,1)||_{L2} ||e(C^*)||_{L2} + ||e(u1 C,1_h)||_{L2} ||e(C^*)||_{L2} ||e((u_2)_h C,2)||_{L2} ||e(C^*)||_{L2} + ||e(u2 C,2_h)||_{L2} ||e(C^*)||_{L2} These functions compute (u_1)_h or C,1_h , and (u_2)_h or C,2_h , and supply it to the weighted patch recovery error estimator In CoupledFEMFunctionsx, the object built with var = 0, returns the (u_1)_h weight, while the object built with var = 1, returns the C,1_h weight. The switch statment distinguishes the behavior of the two objects Same thing for CoupledFEMFunctionsy
        Number CoupledFEMFunctionsx::operator()(const FEMContext& c, const Point& p,
        					const Real /* time */)
        {
          Number weight = 0.0;
        
          switch(var)
            {
            case 0:
              {
        	Gradient grad_C = c.point_gradient(3, p);
        
        	weight = grad_C(0);
              }
              break;
        
            case 3:
              {
        	Number u = c.point_value(0, p);
        
        	weight = u;
              }
              break;
        
            default:
              {
        	std::cout<<"Wrong variable number"<<var<<" passed to CoupledFEMFunctionsx object ! Quitting !"<<std::endl;
        	libmesh_error();
              }
        
            }
        
          return weight;
        }
        
        Number CoupledFEMFunctionsy::operator()(const FEMContext& c, const Point& p,
        					const Real /* time */)
        {
          Number weight = 0.0;
        
          switch(var)
            {
            case 1:
              {
        	Gradient grad_C = c.point_gradient(3, p);
        
        	weight = grad_C(1);
              }
              break;
        
            case 3:
              {
        	Number v = c.point_value(1, p);
        
        	weight = v;
              }
              break;
        
            default:
              {
        	std::cout<<"Wrong variable number "<<var<<" passed to CoupledFEMFunctionsy object ! Quitting !"<<std::endl;
        	libmesh_error();
              }
            }
        
          return weight;
        }



The source file domain.C with comments:

        #include "libmesh/mesh.h"
        #include "libmesh/serial_mesh.h"
        #include "libmesh/mesh_generation.h"
        #include "libmesh/mesh_modification.h"
        #include "libmesh/mesh_refinement.h"
        
        #include "domain.h"
        
Bring in everything from the libMesh namespace
        using namespace libMesh;
        
        
        
Mesh construction
        void build_domain (Mesh &mesh, FEMParameters ¶m)
        {
          mesh.read(param.domainfile);
        
          std::cout<<"Making elements 2nd order"<<std::endl;
        
Right now we are setting approximation orders in the code, rather than reading them in That needs to be fixed and the second ordering should be done only if one of the approximation orders is greater than 1
          mesh.all_second_order();
        
        }



The source file femparameters.C with comments:

        #include "femparameters.h"
        
        using namespace libMesh;
        
        #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(initial_timestep);
            GETPOT_INT_INPUT(n_timesteps);
            GETPOT_INPUT(transient);
            GETPOT_INT_INPUT(deltat_reductions);
            GETPOT_INPUT(timesolver_core);
            GETPOT_INPUT(end_time);
            GETPOT_INPUT(deltat);
            GETPOT_INPUT(timesolver_theta);
            GETPOT_INPUT(timesolver_maxgrowth);
            GETPOT_INPUT(timesolver_upper_tolerance);
            GETPOT_INPUT(timesolver_tolerance);
            GETPOT_INPUT(steadystate_tolerance);
        
            GETPOT_REGISTER(timesolver_norm);
            const unsigned int n_timesolver_norm = input.vector_variable_size("timesolver_norm");
            timesolver_norm.resize(n_timesolver_norm, L2);
            for (unsigned int i=0; i != n_timesolver_norm; ++i)
              {
                int current_norm = 0; // L2
                if (timesolver_norm[i] == H1)
                  current_norm = 1;
                if (timesolver_norm[i] == H2)
                  current_norm = 2;
                current_norm = input("timesolver_norm", current_norm, i);
                if (current_norm == 0)
                  timesolver_norm[i] = L2;
                else if (current_norm == 1)
                  timesolver_norm[i] = H1;
                else if (current_norm == 2)
                  timesolver_norm[i] = H2;
                else
                  timesolver_norm[i] = DISCRETE_L2;
              }
        
            GETPOT_INT_INPUT(dimension);
            GETPOT_INPUT(domaintype);
            GETPOT_INPUT(domainfile);
            GETPOT_INPUT(elementtype);
            GETPOT_INPUT(fine_mesh_file_primal);
            GETPOT_INPUT(fine_mesh_soln_primal);
            GETPOT_INPUT(fine_mesh_file_adjoint);
            GETPOT_INPUT(fine_mesh_soln_adjoint);
            GETPOT_INPUT(elementorder);
            GETPOT_INPUT(domain_xmin);
            GETPOT_INPUT(domain_ymin);
            GETPOT_INPUT(domain_zmin);
            GETPOT_INPUT(domain_edge_width);
            GETPOT_INPUT(domain_edge_length);
            GETPOT_INPUT(domain_edge_height);
            GETPOT_INT_INPUT(coarsegridx);
            GETPOT_INT_INPUT(coarsegridy);
            GETPOT_INT_INPUT(coarsegridz);
            GETPOT_INT_INPUT(coarserefinements);
            GETPOT_INT_INPUT(coarsecoarsenings);
            GETPOT_INT_INPUT(extrarefinements);
            GETPOT_INPUT(use_petsc_snes);
            GETPOT_INPUT(time_solver_quiet);
            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_INT_INPUT(initial_adaptivesteps);
            GETPOT_INT_INPUT(initial_sobolev_order);
            GETPOT_INT_INPUT(initial_extra_quadrature);
            GETPOT_INPUT(refine_uniformly);
            GETPOT_INPUT(indicator_type);
            GETPOT_INPUT(adjoint_residual_type);
            GETPOT_INPUT(patch_reuse);
            GETPOT_INT_INPUT(sobolev_order);
            GETPOT_INPUT(alternate_with_uniform_steps);
            GETPOT_INPUT(alternate_step_number);
            GETPOT_INPUT(component_wise_error);
            GETPOT_INPUT(compute_sensitivities);
            GETPOT_INPUT(compare_to_fine_solution_primal);
            GETPOT_INPUT(compare_to_fine_solution_adjoint);
            GETPOT_INPUT(do_forward_sensitivity);
            GETPOT_INPUT(do_adjoint_sensitivity);
            GETPOT_INPUT(postprocess_adjoint);
            GETPOT_INT_INPUT(write_interval);
            GETPOT_INPUT(output_xda);
            GETPOT_INPUT(output_xdr);
        
            GETPOT_INPUT(output_gz);
        #ifndef LIBMESH_HAVE_GZSTREAM
            output_gz                   = false;
        #endif
        
            GETPOT_INPUT(output_bz2);
        #ifndef LIBMESH_HAVE_BZ2
            output_bz2                  = false;
        #endif
        
            GETPOT_INPUT(output_gmv);
            GETPOT_INPUT(write_gmv_error);
        
        #ifndef LIBMESH_HAVE_GMV
            output_gmv                  = false;
            write_gmv_error             = false;
        #endif
        
            GETPOT_INPUT(output_tecplot);
            GETPOT_INPUT(write_tecplot_error);
        #ifndef LIBMESH_HAVE_TECPLOT_API
            output_tecplot              = false;
            write_tecplot_error         = false;
        #endif
        
            GETPOT_INPUT(run_simulation);
            GETPOT_INPUT(run_postprocess);
        
            run_simulation              = input("run_simulation", run_simulation);
            run_postprocess             = input("run_postprocess", run_postprocess);
        
            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(extra_quadrature_order);
            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 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 H-qoi.C with comments:

        #include "H-qoi.h"
        
Here we define the functions to compute the QoI (side_qoi) and supply the right hand side for the associated adjoint problem (side_qoi_derivative)

        using namespace libMesh;
        
        void CoupledSystemQoI::init_qoi( std::vector<Number>& sys_qoi)
        {
Only 1 qoi to worry about
          sys_qoi.resize(1);
          return;
        }
        
This function supplies the right hand side for the adjoint problem We only have one QoI, so we don't bother checking the qois argument to see if it was requested from us
        void CoupledSystemQoI::side_qoi_derivative (DiffContext &context,
        					 const QoISet & /* qois */)
        {
          FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
        
Element Jacobian * quadrature weights for interior integration
          const std::vector<Real> &JxW = c.side_fe_var[0]->get_JxW();
        
Get velocity basis functions phi
          const std::vector<std::vector<Real> > &phi = c.side_fe_var[0]->get_phi();
        
          const std::vector<Point > &q_point = c.side_fe_var[0]->get_xyz();
        
The number of local degrees of freedom in each variable
          const unsigned int n_u_dofs = c.dof_indices_var[1].size();
        
          DenseSubVector<Number> &Qu = *c.elem_qoi_subderivatives[0][0];
          DenseSubVector<Number> &QC = *c.elem_qoi_subderivatives[0][3];
        
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.side_qrule->n_points();
        
          Number u = 0. ;
          Number C = 0. ;
        
// If side is on outlet
          if(c.has_side_boundary_id(2))
            {
Loop over all the qps on this side
              for (unsigned int qp=0; qp != n_qpoints; qp++)
          	{
        	  Real x = q_point[qp](0);
        
If side is on left outlet
                  if(x < 0.)
        	    {
Get u at the qp
                      u = c.side_value(0,qp);
        	      C = c.side_value(3,qp);
        
Add the contribution from each basis function
                      for (unsigned int i=0; i != n_u_dofs; i++)
        		{
        		  Qu(i) += JxW[qp] * -phi[i][qp] * C;
        		  QC(i) += JxW[qp] * phi[i][qp] * -u;
        		}
        	    } // end if
        
          	} // end quadrature loop
        
            } // end if on outlet
        }
        
This function computes the actual QoI
        void CoupledSystemQoI::side_qoi(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();
        
          const std::vector<Point > &q_point = c.side_fe_var[0]->get_xyz();
        
Loop over qp's, compute the function at each qp and add to get the QoI

          unsigned int n_qpoints = c.side_qrule->n_points();
        
          Number dQoI_0 = 0. ;
          Number u = 0. ;
          Number C = 0. ;
        
If side is on the left outlet
          if(c.has_side_boundary_id(2))
            {
Loop over all the qps on this side
              for (unsigned int qp=0; qp != n_qpoints; qp++)
        	{
        	  Real x = q_point[qp](0);
        
        	  if(x < 0.)
        	    {
Get u and C at the qp
                      u = c.side_value(0,qp);
        	      C = c.side_value(3,qp);
        
        	      dQoI_0 += JxW[qp] * -u * C;
        	    } // end if
        
          	} // end quadrature loop
        
            } // end if on bdry
        
          c.elem_qoi[0] += dQoI_0;
        
        }



The source file initial.C with comments:

        #include "initial.h"
        #include "femparameters.h"
        
Bring in everything from the libMesh namespace
        using namespace libMesh;
        
        void read_initial_parameters()
        {
        }
        
        void finish_initialization()
        {
        }
        
        
        
Initial conditions
        Number initial_value(const Point& /* p */,
                             const Parameters& /* param */,
                             const std::string&,
                             const std::string&)
        {
        
          return 1.;
        
        }
        
        
        
        Gradient initial_grad(const Point& /* p */,
                              const Parameters& /* param */,
                              const std::string&,
                              const std::string&)
        {
          return 0.;
        }



The source file output.C with comments:

        #include <fstream>
        #include <unistd.h>
        
        #include "libmesh/libmesh_common.h"
        
        #include "output.h"
        
Bring in everything from the libMesh namespace
        using namespace libMesh;
        
        void start_output(unsigned int timesteps,
                          std::string filename,
                          std::string varname)
        {
Truncate then append if we're doing a restart
          if (timesteps)
            {
              std::ifstream infile (filename.c_str());
              char buf[1025];
              if (!infile.is_open())
                libmesh_error();
        
Count off the lines we want to save, including the original header
              for (unsigned int i=0; i != timesteps+1; ++i)
                infile.getline(buf, 1024);
              if (!infile.good())
                libmesh_error();
              unsigned int length = infile.tellg();
              infile.close();
        
Then throw away the rest
              int err = truncate(filename.c_str(), length);
              if (err != 0)
                libmesh_error();
            }
Write new headers if we're starting a new simulation
          else
            {
              std::ofstream outfile (filename.c_str(), std::ios_base::out);
              outfile << varname << " = [" << std::endl;
            }
        }



The source file coupled_system.h without comments:

 
  #ifndef __coupled_system_h__
  #define __coupled_system_h__
  
  #include "libmesh/fem_function_base.h"
  #include "libmesh/fem_system.h"
  #include "libmesh/libmesh_common.h"
  #include "libmesh/parameter_vector.h"
  
  using namespace libMesh;
  
  class CoupledSystem : public FEMSystem
  {
  public:
    CoupledSystem(EquationSystems& es,
                 const std::string& name_in,
                 const unsigned int number_in)
      : FEMSystem(es, name_in, number_in), Peclet(1.) {qoi.resize(1);}
  
  
    Number &get_QoI_value()
      {
        return computed_QoI;
      }
  
    Number &get_parameter_value(unsigned int parameter_index)
      {
        return parameters[parameter_index];
      }
  
    ParameterVector &get_parameter_vector()
      {
        parameter_vector.resize(parameters.size());
        for(unsigned int i = 0; i != parameters.size(); ++i)
  	{
  	  parameter_vector[i] = &parameters[i];
  	}
  
        return parameter_vector;
      }
  
    Real &get_Pe()
      {
        return Peclet;
      }
  
   protected:
  
    virtual void init_data ();
  
    virtual void init_context(DiffContext &context);
  
    virtual bool element_time_derivative (bool request_jacobian,
                                          DiffContext& context);
  
    virtual bool element_constraint (bool request_jacobian,
                                     DiffContext& context);
  
    virtual void postprocess ();
  
    std::vector<Number> parameters;
  
    unsigned int p_var, u_var, v_var, C_var;
  
    ParameterVector parameter_vector;
  
    Real Peclet;
  
    Number computed_QoI;
  
  };
  
  
  class CoupledFEMFunctionsx : public FEMFunctionBase<Number>
  {
  public:
    CoupledFEMFunctionsx(System& /* sys */, unsigned int var_number) {var = var_number;}
  
    virtual ~CoupledFEMFunctionsx () {}
  
    virtual AutoPtr<FEMFunctionBase<Number> > clone () const
    {return AutoPtr<FEMFunctionBase<Number> >( new CoupledFEMFunctionsx(*this) ); }
  
    virtual void operator() (const FEMContext&, const Point&,
  			   const Real, DenseVector<Number>&)
    {libmesh_error();}
  
    virtual Number operator() (const FEMContext&, const Point& p,
  			     const Real time = 0.);
  
   private:
  
    unsigned int var;
  
  };
  
  
  class CoupledFEMFunctionsy : public FEMFunctionBase<Number>
  {
  public:
    CoupledFEMFunctionsy(System& /* sys */, unsigned int var_number) {var = var_number;}
  
    virtual ~CoupledFEMFunctionsy () {}
  
    virtual AutoPtr<FEMFunctionBase<Number> > clone () const
    {return AutoPtr<FEMFunctionBase<Number> >( new CoupledFEMFunctionsy(*this) ); }
  
    virtual void operator() (const FEMContext&, const Point&,
  			   const Real,
  			   DenseVector<Number>&)
    {libmesh_error();}
  
    virtual Number operator() (const FEMContext&, const Point& p,
  			     const Real time = 0.);
  
   private:
  
    unsigned int var;
  
  };
  
  #endif //__coupled_system_h__



The source file domain.h without comments:

 
  #include "femparameters.h"
  #include "libmesh/mesh.h"
  
  class FEMParameters;
  
  void build_domain (Mesh &mesh, FEMParameters &param);



The source file femparameters.h without comments:

 
  #ifndef __fem_parameters_h__
  #define __fem_parameters_h__
  
  #include <limits>
  
  #include "libmesh/libmesh_common.h"
  #include "libmesh/dof_map.h"
  #include "libmesh/enum_norm_type.h"
  #include "libmesh/getpot.h"
  
  using namespace libMesh;
  
  class FEMParameters
  {
  public:
      FEMParameters() :
        initial_timestep(0), n_timesteps(100),
        transient(true),
        deltat_reductions(0),
        timesolver_core("euler"),
        end_time(std::numeric_limits<Real>::max()),
        deltat(0.0001), timesolver_theta(0.5),
        timesolver_maxgrowth(0.), timesolver_tolerance(0.),
        timesolver_upper_tolerance(0.),
        steadystate_tolerance(0.),
        timesolver_norm(0,L2),
        dimension(2),
  	domaintype("square"), domainfile("mesh.xda"), elementtype("quad"),
  	fine_mesh_file_primal("fine_mesh.xda"), fine_mesh_soln_primal("fine_mesh_soln.xda"),
  	fine_mesh_file_adjoint("fine_mesh.xda"), fine_mesh_soln_adjoint("fine_mesh_soln.xda"),
        elementorder(2),
        domain_xmin(0.0), domain_ymin(0.0), domain_zmin(0.0),
        domain_edge_width(1.0), domain_edge_length(1.0), domain_edge_height(1.0),
        coarsegridx(1), coarsegridy(1), coarsegridz(1),
  	coarserefinements(0), coarsecoarsenings(0), extrarefinements(0),
        use_petsc_snes(false),
        time_solver_quiet(true), 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),
        initial_adaptivesteps(0), initial_sobolev_order(1),
  	initial_extra_quadrature(0),
  	indicator_type("kelly"), patch_reuse(true), sobolev_order(1), adjoint_residual_type("patchpatch"), alternate_with_uniform_steps("false"), alternate_step_number(2),
  	component_wise_error("false"), compare_to_fine_solution_primal(false), compare_to_fine_solution_adjoint(false), compute_sensitivities(false), do_forward_sensitivity(true), do_adjoint_sensitivity(false),
  	postprocess_adjoint(false),
        write_interval(10),
        write_gmv_error(false), write_tecplot_error(false),
        output_xda(false), output_xdr(false),
        output_bz2(true), output_gz(true),
        output_gmv(false), output_tecplot(false),
        run_simulation(true), run_postprocess(false),
        fe_family(1, "LAGRANGE"), fe_order(1, 1),
        extra_quadrature_order(0),
        analytic_jacobians(true), verify_analytic_jacobians(0.0),
        print_solution_norms(false), print_solutions(false),
        print_residual_norms(false), print_residuals(false),
        print_jacobian_norms(false), print_jacobians(false) {}
  
      void read(GetPot &input);
  
  
  
      unsigned int initial_timestep, n_timesteps;
      bool transient;
      unsigned int deltat_reductions;
      std::string timesolver_core;
      Real end_time, deltat, timesolver_theta,
  	 timesolver_maxgrowth, timesolver_tolerance,
  	 timesolver_upper_tolerance, steadystate_tolerance;
      std::vector<FEMNormType> timesolver_norm;
  
      unsigned int dimension;
      std::string domaintype, domainfile, elementtype, fine_mesh_file_primal, fine_mesh_soln_primal, fine_mesh_file_adjoint, fine_mesh_soln_adjoint;
      Real elementorder;
      Real domain_xmin, domain_ymin, domain_zmin;
      Real domain_edge_width, domain_edge_length, domain_edge_height;
      unsigned int coarsegridx, coarsegridy, coarsegridz;
      unsigned int coarserefinements, coarsecoarsenings, extrarefinements;
  
      bool use_petsc_snes;
      bool time_solver_quiet, 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;
      unsigned int initial_adaptivesteps;
      unsigned int initial_sobolev_order;
      unsigned int initial_extra_quadrature;
  
      std::string indicator_type;
      bool patch_reuse;
      unsigned int sobolev_order;
      std::string adjoint_residual_type;
      bool alternate_with_uniform_steps;
      unsigned int alternate_step_number;
      bool component_wise_error;
      bool compare_to_fine_solution_primal, compare_to_fine_solution_adjoint;
  
      bool compute_sensitivities;
      bool do_forward_sensitivity;
      bool do_adjoint_sensitivity;
  
      bool postprocess_adjoint;
  
      unsigned int write_interval;
      bool write_gmv_error, write_tecplot_error, output_xda, output_xdr,
  	 output_bz2, output_gz, output_gmv, output_tecplot;
      bool run_simulation, run_postprocess;
  
      std::vector<std::string> fe_family;
      std::vector<unsigned int> fe_order;
      int extra_quadrature_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 H-qoi.h without comments:

 
  #ifndef H_QOI_H
  #define H_QOI_H
  
  #include "libmesh/libmesh_common.h"
  #include "libmesh/elem.h"
  #include "libmesh/fe_base.h"
  #include "libmesh/fem_context.h"
  #include "libmesh/point.h"
  #include "libmesh/quadrature.h"
  #include "libmesh/diff_qoi.h"
  
  using namespace libMesh;
  
  class CoupledSystemQoI : public DifferentiableQoI
  {
  public:
    CoupledSystemQoI(){}
    virtual ~CoupledSystemQoI(){}
  
    virtual void init_qoi( std::vector<Number>& sys_qoi);
    virtual void postprocess( ){}
  
    virtual void side_qoi_derivative(DiffContext &context, const QoISet & qois);
  
    virtual void side_qoi(DiffContext &context, const QoISet & qois);
  
    virtual AutoPtr<DifferentiableQoI> clone( )
    { AutoPtr<DifferentiableQoI> my_clone( new CoupledSystemQoI );
      *my_clone = *this;
      return my_clone;
    }
  
  };
  #endif // H_QOI_H



The source file initial.h without comments:

 
  #include "libmesh/parameters.h"
  #include "libmesh/point.h"
  #include "libmesh/vector_value.h"
  
  using namespace libMesh;
  
  void read_initial_parameters();
  void finish_initialization();
  
  Number initial_value(const Point& p,
                       const Parameters&,
                       const std::string&,
                       const std::string&);
  
  Gradient initial_grad(const Point& p,
                        const Parameters&,
                        const std::string&,
                        const std::string&);



The source file mysystems.h without comments:

 
  #include "coupled_system.h"
  #include "femparameters.h"
  
  using namespace libMesh;
  
  FEMSystem &build_system(EquationSystems &es, GetPot &, FEMParameters& /* param */)
  {
    CoupledSystem &system = es.add_system<CoupledSystem> ("CoupledSystem");
  
    return system;
  }
  
  Number exact_value(const Point& /* p */,       // xyz location
                     const Parameters& /* param */,  // EquationSystem parameters
                     const std::string&, // sys_name
                     const std::string&) // unknown_name
  {
    std::cout<<"Warning ! No exact value function specified ! Returning 0." << std::endl;
  
    return 0;
  }
  
  Gradient exact_grad(const Point& /* p */,       // xyz location
                      const Parameters& /* param */,  // EquationSystems parameters
                      const std::string&, // sys_name
                      const std::string&) // unknown_name
  {
    std::cout<<"Warning ! No exact value function specified ! Returning 0." << std::endl;
  
    return 0;
  }



The source file output.h without comments:

 
  #include <string>
  
  using namespace libMesh;
  
  void start_output(unsigned int timesteps,
                    std::string filename,
                    std::string varname);



The source file adjoints_ex3.C without comments:

 
  #include <iostream>
  #include <sys/time.h>
  #include <iomanip>
  
  #include "libmesh/equation_systems.h"
  
  #include "libmesh/twostep_time_solver.h"
  #include "libmesh/euler_solver.h"
  #include "libmesh/euler2_solver.h"
  #include "libmesh/steady_solver.h"
  
  #include "libmesh/newton_solver.h"
  #include "libmesh/numeric_vector.h"
  #include "libmesh/petsc_diff_solver.h"
  
  #include "libmesh/mesh.h"
  #include "libmesh/mesh_tools.h"
  #include "libmesh/mesh_base.h"
  #include "libmesh/mesh_refinement.h"
  
  #include "libmesh/system.h"
  #include "libmesh/system_norm.h"
  
  #include "libmesh/adjoint_residual_error_estimator.h"
  #include "libmesh/const_fem_function.h"
  #include "libmesh/error_vector.h"
  #include "libmesh/fem_function_base.h"
  #include "libmesh/getpot.h"
  #include "libmesh/gmv_io.h"
  #include "libmesh/kelly_error_estimator.h"
  #include "libmesh/parameter_vector.h"
  #include "libmesh/patch_recovery_error_estimator.h"
  #include "libmesh/petsc_vector.h"
  #include "libmesh/sensitivity_data.h"
  #include "libmesh/tecplot_io.h"
  #include "libmesh/uniform_refinement_estimator.h"
  #include "libmesh/qoi_set.h"
  #include "libmesh/weighted_patch_recovery_error_estimator.h"
  
  #include "coupled_system.h"
  #include "domain.h"
  #include "initial.h"
  #include "femparameters.h"
  #include "mysystems.h"
  #include "output.h"
  #include "H-qoi.h"
  
  
  
  
  
  
  
  using namespace libMesh;
  
  
  std::string numbered_filename(unsigned int t_step, // The timestep count
                                unsigned int a_step, // The adaptive step count
                                std::string solution_type, // primal or adjoint solve
                                std::string type,
                                std::string extension,
                                FEMParameters &param)
  {
    std::ostringstream file_name;
    file_name << solution_type
              << ".out."
              << type
              << '.'
              << std::setw(3)
              << std::setfill('0')
              << std::right
              << t_step
              << '.'
              << std::setw(2)
              << std::setfill('0')
              << std::right
              << a_step
              << '.'
              << extension;
  
    if (param.output_bz2)
      file_name << ".bz2";
    else if (param.output_gz)
      file_name << ".gz";
    return file_name.str();
  }
  
  
  void write_output(EquationSystems &es,
                    unsigned int t_step, // The timestep count
                    unsigned int a_step, // The adaptive step count
                    std::string solution_type, // primal or adjoint solve
                    FEMParameters &param)
  {
    MeshBase &mesh = es.get_mesh();
  
  #ifdef LIBMESH_HAVE_GMV
    if (param.output_gmv)
      {
        std::ostringstream file_name_gmv;
        file_name_gmv << solution_type
                      << ".out.gmv."
                      << std::setw(3)
                      << std::setfill('0')
                      << std::right
                      << t_step
                      << '.'
                      << std::setw(2)
                      << std::setfill('0')
                      << std::right
                      << a_step;
  
        GMVIO(mesh).write_equation_systems
          (file_name_gmv.str(), es);
      }
  #endif
  
  #ifdef LIBMESH_HAVE_TECPLOT_API
    if (param.output_tecplot)
      {
        std::ostringstream file_name_tecplot;
        file_name_tecplot << solution_type
                          << ".out."
                          << std::setw(3)
                          << std::setfill('0')
                          << std::right
                          << t_step
                          << '.'
                          << std::setw(2)
                          << std::setfill('0')
                          << std::right
                          << a_step
                          << ".plt";
  
        TecplotIO(mesh).write_equation_systems
          (file_name_tecplot.str(), es);
      }
  #endif
  
    if (param.output_xda || param.output_xdr)
      mesh.renumber_nodes_and_elements();
    if (param.output_xda)
      {
        mesh.write(numbered_filename(t_step, a_step, solution_type, "mesh", "xda", param));
        es.write(numbered_filename(t_step, a_step, solution_type, "soln", "xda", param),
                 libMeshEnums::WRITE, EquationSystems::WRITE_DATA |
                 EquationSystems::WRITE_ADDITIONAL_DATA);
      }
    if (param.output_xdr)
      {
        mesh.write(numbered_filename(t_step, a_step, solution_type, "mesh", "xdr", param));
        es.write(numbered_filename(t_step, a_step, solution_type, "soln", "xdr", param),
                 libMeshEnums::WRITE, EquationSystems::WRITE_DATA |
                 EquationSystems::WRITE_ADDITIONAL_DATA);
      }
  }
  
  void write_output_headers(FEMParameters &param)
  {
    if (libMesh::processor_id() != 0)
      return;
  
    start_output(param.initial_timestep, "out_clocktime.m", "vector_clocktime");
  
    if (param.run_simulation)
      {
        start_output(param.initial_timestep, "out_activemesh.m", "table_activemesh");
        start_output(param.initial_timestep, "out_solvesteps.m", "table_solvesteps");
  
        if (param.timesolver_tolerance)
          {
            start_output(param.initial_timestep, "out_time.m", "vector_time");
            start_output(param.initial_timestep, "out_timesteps.m", "vector_timesteps");
          }
        if (param.steadystate_tolerance)
          start_output(param.initial_timestep, "out_changerate.m", "vector_changerate");
      }
  }
  
  void write_output_solvedata(EquationSystems &es,
                              unsigned int a_step,
                              unsigned int newton_steps,
                              unsigned int krylov_steps,
                              unsigned int tv_sec,
                              unsigned int tv_usec)
  {
    MeshBase &mesh = es.get_mesh();
    unsigned int n_active_elem = mesh.n_active_elem();
    unsigned int n_active_dofs = es.n_active_dofs();
  
    if (libMesh::processor_id() == 0)
      {
        std::ofstream activemesh ("out_activemesh.m",
          std::ios_base::app | std::ios_base::out);
        activemesh.precision(17);
        activemesh << (a_step + 1) << ' '
                   << n_active_elem << ' '
                   << n_active_dofs << std::endl;
  
        std::ofstream solvesteps ("out_solvesteps.m",
          std::ios_base::app | std::ios_base::out);
        solvesteps.precision(17);
        solvesteps << newton_steps << ' '
                   << krylov_steps << std::endl;
  
        std::ofstream clocktime ("out_clocktime.m",
          std::ios_base::app | std::ios_base::out);
        clocktime.precision(17);
        clocktime << tv_sec << '.' << tv_usec << std::endl;
      }
  }
  
  void write_output_footers(FEMParameters &param)
  {
    if (libMesh::processor_id() == 0)
      {
        std::ofstream clocktime ("out_clocktime.m",
          std::ios_base::app | std::ios_base::out);
        clocktime << "];" << std::endl;
  
        if (param.run_simulation)
          {
            std::ofstream activemesh ("out_activemesh.m",
              std::ios_base::app | std::ios_base::out);
            activemesh << "];" << std::endl;
  
            std::ofstream solvesteps ("out_solvesteps.m",
              std::ios_base::app | std::ios_base::out);
            solvesteps << "];" << std::endl;
  
            if (param.timesolver_tolerance)
              {
                std::ofstream times ("out_time.m",
                  std::ios_base::app | std::ios_base::out);
                times << "];" << std::endl;
                std::ofstream timesteps ("out_timesteps.m",
                  std::ios_base::app | std::ios_base::out);
                timesteps << "];" << std::endl;
              }
            if (param.steadystate_tolerance)
              {
                std::ofstream changerate ("out_changerate.m",
                  std::ios_base::app | std::ios_base::out);
                changerate << "];" << std::endl;
              }
          }
  
        std::ofstream complete ("complete");
        complete << "complete" << std::endl;
      }
  }
  
  #if defined(LIBMESH_HAVE_GMV) || defined(LIBMESH_HAVE_TECPLOT_API)
  void write_error(EquationSystems &es,
                   ErrorVector &error,
                   unsigned int t_number,
                   unsigned int a_number,
                   FEMParameters &param,
                   std::string error_type)
  #else
  void write_error(EquationSystems &,
                   ErrorVector &,
                   unsigned int,
                   unsigned int,
                   FEMParameters &,
                   std::string)
  #endif
  {
  #ifdef LIBMESH_HAVE_GMV
    if (param.write_gmv_error)
      {
        std::ostringstream error_gmv;
        error_gmv << "error.gmv."
                  << std::setw(3)
                  << std::setfill('0')
                  << std::right
                  << a_number
                  << "."
                  << std::setw(2)
                  << std::setfill('0')
                  << std::right
                  << t_number
                  << error_type;
  
        error.plot_error(error_gmv.str(), es.get_mesh());
      }
  #endif
  
  #ifdef LIBMESH_HAVE_TECPLOT_API
    if (param.write_tecplot_error)
      {
        std::ostringstream error_tecplot;
        error_tecplot << "error.plt."
                      << std::setw(3)
                      << std::setfill('0')
                      << std::right
                      << a_number
                      << "."
                      << std::setw(2)
                      << std::setfill('0')
                      << std::right
                      << t_number
                      << error_type;
  
        error.plot_error(error_tecplot.str(), es.get_mesh());
      }
  #endif
  }
  
  void read_output(EquationSystems &es,
                   unsigned int t_step,
                   unsigned int a_step,
                   std::string solution_type,
                   FEMParameters &param)
  {
    MeshBase &mesh = es.get_mesh();
  
    std::string file_name_mesh, file_name_soln;
    if (param.output_xda)
      {
        file_name_mesh = numbered_filename(t_step, a_step, solution_type, "mesh", "xda", param);
        file_name_soln = numbered_filename(t_step, a_step, solution_type, "soln", "xda", param);
      }
    else if (param.output_xdr)
      {
        file_name_mesh = numbered_filename(t_step, a_step, solution_type, "mesh", "xdr", param);
        file_name_soln = numbered_filename(t_step, a_step, solution_type, "soln", "xdr", param);
      }
  
    mesh.read(file_name_mesh);
  
    es.read(file_name_soln, libMeshEnums::READ,
            EquationSystems::READ_HEADER |
            EquationSystems::READ_DATA |
            EquationSystems::READ_ADDITIONAL_DATA);
  
    for (unsigned int i = 0; i != es.n_systems(); ++i)
      es.get_system<FEMSystem>(i).update();
  
    Real current_time = 0., current_timestep = 0.;
  
    if (param.timesolver_tolerance)
      {
        std::ifstream times ("out_time.m");
        std::ifstream timesteps ("out_timesteps.m");
        if (times.is_open() && timesteps.is_open())
          {
            const unsigned int headersize = 25;
            char header[headersize];
            timesteps.getline (header, headersize);
            if (strcmp(header, "vector_timesteps = [") != 0)
              {
                std::cout << "Bad header in out_timesteps.m:" << std::endl
                          << header
                          << std::endl;
                libmesh_error();
              }
  
            times.getline (header, headersize);
            if (strcmp(header, "vector_time = [") != 0)
              {
                std::cout << "Bad header in out_time.m:" << std::endl
                          << header
                          << std::endl;
                libmesh_error();
              }
  
            for (unsigned int i = 0; i != t_step; ++i)
              {
                if (!times.good())
                  libmesh_error();
                times >> current_time;
                timesteps >> current_timestep;
              }
            current_time += current_timestep;
          }
        else
          libmesh_error();
      }
    else
      current_time = t_step * param.deltat;
  
    for (unsigned int i = 0; i != es.n_systems(); ++i)
      es.get_system<FEMSystem>(i).time = current_time;
  }
  
  void set_system_parameters(FEMSystem &system, FEMParameters &param)
  {
    system.verify_analytic_jacobians = param.verify_analytic_jacobians;
  
    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;
  
    if (param.transient)
      {
        UnsteadySolver *innersolver;
        if (param.timesolver_core == "euler2")
          {
            Euler2Solver *euler2solver =
              new Euler2Solver(system);
  
            euler2solver->theta = param.timesolver_theta;
            innersolver = euler2solver;
          }
        else if (param.timesolver_core == "euler")
          {
            EulerSolver *eulersolver =
              new EulerSolver(system);
  
            eulersolver->theta = param.timesolver_theta;
            innersolver = eulersolver;
          }
        else
          {
            std::cerr << "Don't recognize core TimeSolver type: "
                      << param.timesolver_core << std::endl;
            libmesh_error();
          }
  
        if (param.timesolver_tolerance)
          {
            TwostepTimeSolver *timesolver =
              new TwostepTimeSolver(system);
  
            timesolver->max_growth       = param.timesolver_maxgrowth;
            timesolver->target_tolerance = param.timesolver_tolerance;
            timesolver->upper_tolerance  = param.timesolver_upper_tolerance;
            timesolver->component_norm   = SystemNorm(param.timesolver_norm);
  
            timesolver->core_time_solver =
              AutoPtr<UnsteadySolver>(innersolver);
            system.time_solver =
              AutoPtr<UnsteadySolver>(timesolver);
          }
        else
          system.time_solver =
            AutoPtr<TimeSolver>(innersolver);
      }
    else
      system.time_solver =
        AutoPtr<TimeSolver>(new SteadySolver(system));
  
    system.time_solver->reduce_deltat_on_diffsolver_failure =
                                          param.deltat_reductions;
    system.time_solver->quiet           = param.time_solver_quiet;
  
    system.deltat = param.deltat;
  
    system.extra_quadrature_order = param.extra_quadrature_order;
  
    if (param.use_petsc_snes)
      {
  #ifdef LIBMESH_HAVE_PETSC
        PetscDiffSolver *solver = new PetscDiffSolver(system);
        system.time_solver->diff_solver() = AutoPtr<DiffSolver>(solver);
  #else
        libmesh_error();
  #endif
      }
    else
      {
        NewtonSolver *solver = new NewtonSolver(system);
        system.time_solver->diff_solver() = AutoPtr<DiffSolver>(solver);
  
        solver->quiet                       = param.solver_quiet;
        solver->max_nonlinear_iterations    = param.max_nonlinear_iterations;
        solver->minsteplength               = param.min_step_length;
        solver->relative_step_tolerance     = param.relative_step_tolerance;
        solver->relative_residual_tolerance = param.relative_residual_tolerance;
        solver->require_residual_reduction  = param.require_residual_reduction;
        solver->linear_tolerance_multiplier = param.linear_tolerance_multiplier;
        if (system.time_solver->reduce_deltat_on_diffsolver_failure)
          {
            solver->continue_after_max_iterations = true;
            solver->continue_after_backtrack_failure = true;
          }
  
        solver->max_linear_iterations       = param.max_linear_iterations;
        solver->initial_linear_tolerance    = param.initial_linear_tolerance;
        solver->minimum_linear_tolerance    = param.minimum_linear_tolerance;
      }
  }
  
  #ifdef LIBMESH_ENABLE_AMR
  
  AutoPtr<MeshRefinement> build_mesh_refinement(MeshBase &mesh,
                                                FEMParameters &param)
  {
    AutoPtr<MeshRefinement> mesh_refinement(new MeshRefinement(mesh));
    mesh_refinement->coarsen_by_parents() = true;
    mesh_refinement->absolute_global_tolerance() = param.global_tolerance;
    mesh_refinement->nelem_target()      = param.nelem_target;
    mesh_refinement->refine_fraction()   = param.refine_fraction;
    mesh_refinement->coarsen_fraction()  = param.coarsen_fraction;
    mesh_refinement->coarsen_threshold() = param.coarsen_threshold;
  
    return mesh_refinement;
  }
  
  #endif // LIBMESH_ENABLE_AMR
  
  AutoPtr<ErrorEstimator> build_error_estimator(FEMParameters& /* param */)
  {
    AutoPtr<ErrorEstimator> error_estimator;
  
    error_estimator.reset(new KellyErrorEstimator);
  
    return error_estimator;
  }
  
  AutoPtr<ErrorEstimator>
  build_error_estimator_component_wise
    (FEMParameters &param,
     std::vector<std::vector<Real> > &term_weights,
     std::vector<libMeshEnums::FEMNormType> &primal_error_norm_type,
     std::vector<libMeshEnums::FEMNormType> &dual_error_norm_type)
  {
    AutoPtr<ErrorEstimator> error_estimator;
  
    AdjointResidualErrorEstimator *adjoint_residual_estimator = new AdjointResidualErrorEstimator;
  
    error_estimator.reset (adjoint_residual_estimator);
  
    PatchRecoveryErrorEstimator *p1 =
      new PatchRecoveryErrorEstimator;
    adjoint_residual_estimator->primal_error_estimator().reset(p1);
  
    PatchRecoveryErrorEstimator *p2 =
      new PatchRecoveryErrorEstimator;
    adjoint_residual_estimator->dual_error_estimator().reset(p2);
  
    p1->set_patch_reuse(param.patch_reuse);
    p2->set_patch_reuse(param.patch_reuse);
  
    unsigned int size = primal_error_norm_type.size();
  
    libmesh_assert_equal_to (size, dual_error_norm_type.size());
    for (unsigned int i = 0; i != size; ++i)
      {
        adjoint_residual_estimator->primal_error_estimator()->error_norm.set_type(i, primal_error_norm_type[i]);
        adjoint_residual_estimator->dual_error_estimator()->error_norm.set_type(i, dual_error_norm_type[i]);
      }
  
    libmesh_assert_equal_to (size, term_weights.size());
    for (unsigned int i = 0; i != term_weights.size(); ++i)
      {
        libmesh_assert_equal_to (size, term_weights[i].size());
        adjoint_residual_estimator->error_norm.set_weight(i, term_weights[i][i]);
        for (unsigned int j = 0; j != size; ++j)
          if (i != j)
            adjoint_residual_estimator->error_norm.set_off_diagonal_weight(i, j, term_weights[i][j]);
      }
  
    return error_estimator;
  }
  
  AutoPtr<ErrorEstimator>
  build_weighted_error_estimator_component_wise
    (FEMParameters &param,
     std::vector<std::vector<Real> > &term_weights,
     std::vector<libMeshEnums::FEMNormType> &primal_error_norm_type,
     std::vector<libMeshEnums::FEMNormType> &dual_error_norm_type,
     std::vector<FEMFunctionBase<Number>*> coupled_system_weight_functions)
  {
    AutoPtr<ErrorEstimator> error_estimator;
  
    AdjointResidualErrorEstimator *adjoint_residual_estimator = new AdjointResidualErrorEstimator;
  
    error_estimator.reset (adjoint_residual_estimator);
  
  
    WeightedPatchRecoveryErrorEstimator *p1 =
      new WeightedPatchRecoveryErrorEstimator;
    adjoint_residual_estimator->primal_error_estimator().reset(p1);
  
    PatchRecoveryErrorEstimator *p2 =
      new PatchRecoveryErrorEstimator;
    adjoint_residual_estimator->dual_error_estimator().reset(p2);
  
    p1->set_patch_reuse(param.patch_reuse);
    p2->set_patch_reuse(param.patch_reuse);
  
    p1->weight_functions.clear();
  
    unsigned int size = primal_error_norm_type.size();
  
    libmesh_assert(coupled_system_weight_functions.size() == size);
    libmesh_assert(dual_error_norm_type.size() == size);
  
    for (unsigned int i = 0; i != size; ++i)
      {
        p1->weight_functions.push_back(coupled_system_weight_functions[i]);
        adjoint_residual_estimator->primal_error_estimator()->error_norm.set_type(i, primal_error_norm_type[i]);
        adjoint_residual_estimator->dual_error_estimator()->error_norm.set_type(i, dual_error_norm_type[i]);
      }
  
    libmesh_assert_equal_to (size, term_weights.size());
    for (unsigned int i = 0; i != term_weights.size(); ++i)
      {
        libmesh_assert_equal_to (size, term_weights[i].size());
        adjoint_residual_estimator->error_norm.set_weight(i, term_weights[i][i]);
        for (unsigned int j = 0; j != size; ++j)
          if (i != j)
            adjoint_residual_estimator->error_norm.set_off_diagonal_weight(i, j, term_weights[i][j]);
      }
  
    return error_estimator;
  }
  
  int main (int argc, char** argv)
  {
    LibMeshInit init (argc, argv);
  
  #ifndef LIBMESH_ENABLE_AMR
    libmesh_example_assert(false, "--enable-amr");
  #else
  
    libmesh_example_assert(libMesh::default_solver_package() != EIGEN_SOLVERS, "--enable-petsc");
  
    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(), param.dimension);
  
    AutoPtr<MeshRefinement> mesh_refinement =
      build_mesh_refinement(mesh, param);
  
    EquationSystems equation_systems (mesh);
  
    std::cout << "Building mesh" << std::endl;
  
    if (!param.initial_timestep && param.run_simulation)
      build_domain(mesh, param);
  
    std::cout << "Building system" << std::endl;
  
    FEMSystem &system = build_system(equation_systems, infile, param);
  
    set_system_parameters(system, param);
  
    std::cout << "Initializing systems" << std::endl;
  
    {
      CoupledSystemQoI qoi;
      qoi.assemble_qoi_sides = true;
      system.attach_qoi(&qoi);
    }
  
    equation_systems.init ();
  
    mesh.print_info();
    equation_systems.print_info();
  
    std::cout<<"Starting adaptive loop"<<std::endl<<std::endl;
  
    unsigned int a_step = 0;
  
    LinearSolver<Number> *linear_solver = system.get_linear_solver();
  
    for (; a_step <= param.max_adaptivesteps; ++a_step)
      {
        if(param.global_tolerance > 0 && param.nelem_target > 0)
          {
            std::cout<<"We cant refine to both a non-zero tolerance and a target number of elements, EXITING adaptive loop. "<<std::endl<<std::endl;
            break;
          }
  
        linear_solver->reuse_preconditioner(false);
  
        std::cout<< "Adaptive step " << a_step << std::endl;
  
        std::cout<< "Solving the forward problem" <<std::endl;
  
        std::cout << "We have " << mesh.n_active_elem()
                      << " active elements and " << equation_systems.n_active_dofs()
                  << " active dofs." << std::endl << std::endl;
  
        system.solve();
  
        write_output(equation_systems, 0, a_step, "primal", param);
  
        system.assemble_qoi();
  
        system.postprocess();
  
        Number QoI_0_computed = (dynamic_cast<CoupledSystem&>(system)).get_QoI_value();
  
        std::cout<< "The boundary QoI is " << std::setprecision(17) << QoI_0_computed << std::endl << std::endl;
  
        if(a_step!=param.max_adaptivesteps)
          {
            NumericVector<Number> &primal_solution = (*system.solution);
  
            system.assemble_qoi_sides = true;
  
            linear_solver->reuse_preconditioner(param.reuse_preconditioner);
  
            std::cout<< "Solving the adjoint problem" <<std::endl;
            system.adjoint_solve();
  
  	system.set_adjoint_already_solved(true);
  
            NumericVector<Number> &dual_solution = system.get_adjoint_solution();
            primal_solution.swap(dual_solution);
  
            write_output(equation_systems, 0, a_step, "adjoint", param);
  
            primal_solution.swap(dual_solution);
  
            Real Pe = (dynamic_cast<CoupledSystem&>(system)).get_Pe();
  
            ErrorVector error;
  
            ErrorVector error_non_pressure;
  
            std::vector<libMeshEnums::FEMNormType>
              primal_norm_type_vector_non_pressure;
            primal_norm_type_vector_non_pressure.push_back(H1_SEMINORM);
            primal_norm_type_vector_non_pressure.push_back(H1_SEMINORM);
            primal_norm_type_vector_non_pressure.push_back(L2);
            primal_norm_type_vector_non_pressure.push_back(H1_SEMINORM);
  
            std::vector<libMeshEnums::FEMNormType>
              dual_norm_type_vector_non_pressure;
            dual_norm_type_vector_non_pressure.push_back(H1_SEMINORM);
            dual_norm_type_vector_non_pressure.push_back(H1_SEMINORM);
            dual_norm_type_vector_non_pressure.push_back(L2);
            dual_norm_type_vector_non_pressure.push_back(H1_SEMINORM);
  
            std::vector<std::vector<Real> >
              weights_matrix_non_pressure(system.n_vars(),
                std::vector<Real>(system.n_vars(), 0.0));
            weights_matrix_non_pressure[0][0] = 1.;
            weights_matrix_non_pressure[1][1] = 1.;
            weights_matrix_non_pressure[3][3] = 1./Pe;
  
            AutoPtr<ErrorEstimator> error_estimator_non_pressure =
              build_error_estimator_component_wise
                (param, weights_matrix_non_pressure,
                 primal_norm_type_vector_non_pressure,
                 dual_norm_type_vector_non_pressure);
  
            error_estimator_non_pressure->estimate_error(system, error_non_pressure);
  
            write_error(equation_systems, error_non_pressure, 0, a_step, param, "_non_pressure");
  
            ErrorVector error_with_pressure;
  
            std::vector<libMeshEnums::FEMNormType>
              primal_norm_type_vector_with_pressure;
            primal_norm_type_vector_with_pressure.push_back(H1_X_SEMINORM);
            primal_norm_type_vector_with_pressure.push_back(H1_Y_SEMINORM);
            primal_norm_type_vector_with_pressure.push_back(L2);
            primal_norm_type_vector_with_pressure.push_back(L2);
  
            std::vector<libMeshEnums::FEMNormType>
              dual_norm_type_vector_with_pressure;
            dual_norm_type_vector_with_pressure.push_back(H1_X_SEMINORM);
            dual_norm_type_vector_with_pressure.push_back(H1_Y_SEMINORM);
            dual_norm_type_vector_with_pressure.push_back(L2);
            dual_norm_type_vector_with_pressure.push_back(L2);
  
            std::vector<std::vector<Real> >
              weights_matrix_with_pressure
                (system.n_vars(),
                 std::vector<Real>(system.n_vars(), 0.0));
            weights_matrix_with_pressure[0][2] = 1.;
  
            weights_matrix_with_pressure[1][2] = 1.;
  
            weights_matrix_with_pressure[2][0] = 1.;
            weights_matrix_with_pressure[2][1] = 1.;
  
            AutoPtr<ErrorEstimator> error_estimator_with_pressure =
              build_error_estimator_component_wise
                (param, weights_matrix_with_pressure,
                 primal_norm_type_vector_with_pressure,
                 dual_norm_type_vector_with_pressure);
  
            error_estimator_with_pressure->estimate_error(system, error_with_pressure);
  
            write_error(equation_systems, error_with_pressure, 0, a_step, param, "_with_pressure");
  
  
            ErrorVector error_convection_diffusion_x;
  
            std::vector<libMeshEnums::FEMNormType>
              primal_norm_type_vector_convection_diffusion_x;
            primal_norm_type_vector_convection_diffusion_x.push_back(L2);
            primal_norm_type_vector_convection_diffusion_x.push_back(L2);
            primal_norm_type_vector_convection_diffusion_x.push_back(L2);
            primal_norm_type_vector_convection_diffusion_x.push_back(H1_X_SEMINORM);
  
            std::vector<libMeshEnums::FEMNormType>
              dual_norm_type_vector_convection_diffusion_x;
            dual_norm_type_vector_convection_diffusion_x.push_back(L2);
            dual_norm_type_vector_convection_diffusion_x.push_back(L2);
            dual_norm_type_vector_convection_diffusion_x.push_back(L2);
            dual_norm_type_vector_convection_diffusion_x.push_back(L2);
  
            std::vector<std::vector<Real> >
              weights_matrix_convection_diffusion_x
                (system.n_vars(),
                 std::vector<Real>(system.n_vars(), 0.0));
  	  weights_matrix_convection_diffusion_x[0][3] = 1.;
            weights_matrix_convection_diffusion_x[3][3] = 1.;
  
  
            ConstFEMFunction<Number> identity(1);
  
            CoupledFEMFunctionsx convdiffx0(system, 0);
  	  CoupledFEMFunctionsx convdiffx3(system, 3);
  
            std::vector<FEMFunctionBase<Number> *> coupled_system_weight_functions_x;
            coupled_system_weight_functions_x.push_back(&convdiffx0);
            coupled_system_weight_functions_x.push_back(&identity);
            coupled_system_weight_functions_x.push_back(&identity);
            coupled_system_weight_functions_x.push_back(&convdiffx3);
  
            AutoPtr<ErrorEstimator> error_estimator_convection_diffusion_x =
            build_weighted_error_estimator_component_wise
              (param, weights_matrix_convection_diffusion_x,
               primal_norm_type_vector_convection_diffusion_x,
               dual_norm_type_vector_convection_diffusion_x,
               coupled_system_weight_functions_x);
  
            error_estimator_convection_diffusion_x->estimate_error
              (system, error_convection_diffusion_x);
  
            write_error(equation_systems, error_convection_diffusion_x,
                        0, a_step, param, "_convection_diffusion_x");
  
            ErrorVector error_convection_diffusion_y;
  
            std::vector<libMeshEnums::FEMNormType>
              primal_norm_type_vector_convection_diffusion_y;
            primal_norm_type_vector_convection_diffusion_y.push_back(L2);
            primal_norm_type_vector_convection_diffusion_y.push_back(L2);
            primal_norm_type_vector_convection_diffusion_y.push_back(L2);
            primal_norm_type_vector_convection_diffusion_y.push_back(H1_Y_SEMINORM);
  
            std::vector<libMeshEnums::FEMNormType>
              dual_norm_type_vector_convection_diffusion_y;
            dual_norm_type_vector_convection_diffusion_y.push_back(L2);
            dual_norm_type_vector_convection_diffusion_y.push_back(L2);
            dual_norm_type_vector_convection_diffusion_y.push_back(L2);
            dual_norm_type_vector_convection_diffusion_y.push_back(L2);
  
            std::vector<std::vector<Real> >
              weights_matrix_convection_diffusion_y
                (system.n_vars(), std::vector<Real>(system.n_vars(), 0.0));
  	  weights_matrix_convection_diffusion_y[1][3] = 1.;
            weights_matrix_convection_diffusion_y[3][3] = 1.;
  
            CoupledFEMFunctionsy convdiffy1(system, 1);
  	  CoupledFEMFunctionsy convdiffy3(system, 3);
  
            std::vector<FEMFunctionBase<Number> *> coupled_system_weight_functions_y;
            coupled_system_weight_functions_y.push_back(&identity);
            coupled_system_weight_functions_y.push_back(&convdiffy1);
            coupled_system_weight_functions_y.push_back(&identity);
            coupled_system_weight_functions_y.push_back(&convdiffy3);
  
            AutoPtr<ErrorEstimator> error_estimator_convection_diffusion_y =
              build_weighted_error_estimator_component_wise
                (param, weights_matrix_convection_diffusion_y,
                 primal_norm_type_vector_convection_diffusion_y,
                 dual_norm_type_vector_convection_diffusion_y,
                 coupled_system_weight_functions_y);
  
            error_estimator_convection_diffusion_y->estimate_error(system, error_convection_diffusion_y);
  
            write_error(equation_systems, error_convection_diffusion_y, 0, a_step, param, "_convection_diffusion_y");
  
            if(param.indicator_type == "adjoint_residual")
              {
                error.resize(error_non_pressure.size());
  
                for(unsigned int i = 0; i < error.size(); i++)
                  {
                    error[i] = error_non_pressure[i] + error_with_pressure[i] + error_convection_diffusion_x[i] + error_convection_diffusion_y[i];
                  }
              }
            else
              {
                std::cout<<"Using Kelly Estimator"<<std::endl;
                AutoPtr<ErrorEstimator> error_estimator =  build_error_estimator(param);
  
                error_estimator->estimate_error(system, error);
              }
  
            write_error(equation_systems, error, 0, a_step, param, "_total");
  
  
  
            if(param.refine_uniformly)
              {
                mesh_refinement->uniformly_refine(1);
              }
            else if(param.global_tolerance >= 0. && param.nelem_target == 0.) // Refine based on reaching an error tolerance
              {
                mesh_refinement->flag_elements_by_error_tolerance (error);
  
                mesh_refinement->refine_and_coarsen_elements();
              }
            else // Refine based on reaching a target number of elements
              {
                if (mesh.n_active_elem() >= param.nelem_target)
                  {
                    std::cout<<"We reached the target number of elements."<<std::endl <<std::endl;
                    break;
                  }
  
                mesh_refinement->flag_elements_by_nelem_target (error);
  
                mesh_refinement->refine_and_coarsen_elements();
              }
  
            equation_systems.reinit();
  
          }
  
      }
  
  
    write_output_footers(param);
  
  #endif // #ifndef LIBMESH_ENABLE_AMR
  
    return 0;
  
  }



The source file coupled_system.C without comments:

 
  #include "libmesh/boundary_info.h"
  #include "libmesh/dirichlet_boundaries.h"
  #include "libmesh/dof_map.h"
  #include "libmesh/fe_base.h"
  #include "libmesh/fe_interface.h"
  #include "libmesh/fem_context.h"
  #include "libmesh/getpot.h"
  #include "libmesh/mesh.h"
  #include "libmesh/parallel.h"
  #include "libmesh/quadrature.h"
  #include "libmesh/string_to_enum.h"
  #include "libmesh/zero_function.h"
  
  #include "coupled_system.h"
  
  using namespace libMesh;
  
  class BdyFunction : public FunctionBase<Number>
  {
  public:
    BdyFunction (unsigned int u_var, unsigned int v_var, int sign)
      : _u_var(u_var), _v_var(v_var), _sign(sign)
      { this->_initialized = true; }
  
    virtual Number operator() (const Point&, const Real = 0)
      { libmesh_not_implemented(); }
  
    virtual void operator() (const Point& p,
                             const Real,
                             DenseVector<Number>& output)
      {
        output.resize(2);
        output.zero();
        const Real y=p(1);
        output(_u_var) = (_sign)*((y-2) * (y-3));
        output(_v_var) = 0;
      }
  
    virtual AutoPtr<FunctionBase<Number> > clone() const
      { return AutoPtr<FunctionBase<Number> > (new BdyFunction(_u_var, _v_var, _sign)); }
  
  private:
    const unsigned int _u_var, _v_var;
    const Real _sign;
  };
  
  
  void CoupledSystem::init_data ()
  {
    GetPot infile("coupled_system.in");
    Peclet = infile("Peclet", 1.);
    unsigned int pressure_p = infile("pressure_p", 1);
    std::string fe_family = infile("fe_family", std::string("LAGRANGE"));
  
    libmesh_assert((pressure_p == 1) || (fe_family != "LAGRANGE"));
  
    FEFamily fefamily = Utility::string_to_enum<FEFamily>(fe_family);
  
    u_var = this->add_variable ("u", static_cast<Order>(pressure_p+1),
  			      fefamily);
    v_var = this->add_variable ("v", static_cast<Order>(pressure_p+1),
  			      fefamily);
  
    p_var = this->add_variable ("p", static_cast<Order>(pressure_p),
  			      fefamily);
  
    C_var = this->add_variable ("C", static_cast<Order>(pressure_p+1),
  			      fefamily);
  
    this->time_evolving(u_var);
    this->time_evolving(v_var);
    this->time_evolving(C_var);
  
    this->verify_analytic_jacobians = infile("verify_analytic_jacobians", 0.);
    this->print_jacobians = infile("print_jacobians", false);
    this->print_element_jacobians = infile("print_element_jacobians", false);
  
    const boundary_id_type left_inlet_id = 0;
    std::set<boundary_id_type> left_inlet_bdy;
    left_inlet_bdy.insert(left_inlet_id);
  
    const boundary_id_type right_inlet_id = 1;
    std::set<boundary_id_type> right_inlet_bdy;
    right_inlet_bdy.insert(right_inlet_id);
  
    const boundary_id_type outlets_id = 2;
    std::set<boundary_id_type> outlets_bdy;
    outlets_bdy.insert(outlets_id);
  
    const boundary_id_type wall_id = 3;
    std::set<boundary_id_type> wall_bdy;
    wall_bdy.insert(wall_id);
  
    std::vector<unsigned int> uv(1, u_var);
    uv.push_back(v_var);
    std::vector<unsigned int> C_only(1, C_var);
  
    ZeroFunction<Number> zero;
    ConstFunction<Number> one(1);
  
    int velocity_sign = 1;
    BdyFunction inflow_left(u_var, v_var, -velocity_sign);
    BdyFunction inflow_right(u_var, v_var, velocity_sign);
  
    this->get_dof_map().add_dirichlet_boundary
          (DirichletBoundary (wall_bdy, uv, &zero));
  
    this->get_dof_map().add_dirichlet_boundary
      (DirichletBoundary (left_inlet_bdy, uv, &inflow_left));
    this->get_dof_map().add_dirichlet_boundary
      (DirichletBoundary (left_inlet_bdy, C_only, &one));
  
    this->get_dof_map().add_dirichlet_boundary
      (DirichletBoundary (right_inlet_bdy, uv, &inflow_right));
    this->get_dof_map().add_dirichlet_boundary
      (DirichletBoundary (right_inlet_bdy, C_only, &zero));
  
  
    FEMSystem::init_data();
  }
  
  void CoupledSystem::init_context(DiffContext &context)
  {
    FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
  
    c.element_fe_var[u_var]->get_JxW();
    c.element_fe_var[u_var]->get_phi();
    c.element_fe_var[u_var]->get_dphi();
    c.element_fe_var[u_var]->get_xyz();
  
    c.element_fe_var[p_var]->get_phi();
  
    c.side_fe_var[u_var]->get_JxW();
    c.side_fe_var[u_var]->get_phi();
    c.side_fe_var[u_var]->get_xyz();
  }
  
  
  bool CoupledSystem::element_time_derivative (bool request_jacobian,
                                              DiffContext &context)
  {
    FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
  
  
    const std::vector<Real> &JxW =
      c.element_fe_var[u_var]->get_JxW();
  
    const std::vector<std::vector<Real> >& phi =
      c.element_fe_var[u_var]->get_phi();
  
    const std::vector<std::vector<RealGradient> >& dphi =
      c.element_fe_var[u_var]->get_dphi();
  
    const std::vector<std::vector<Real> >& psi =
      c.element_fe_var[p_var]->get_phi();
  
    const unsigned int n_p_dofs = c.dof_indices_var[p_var].size();
    const unsigned int n_u_dofs = c.dof_indices_var[u_var].size();
    libmesh_assert_equal_to (n_u_dofs, c.dof_indices_var[v_var].size());
  
    DenseSubMatrix<Number> &Kuu = *c.elem_subjacobians[u_var][u_var];
    DenseSubMatrix<Number> &Kup = *c.elem_subjacobians[u_var][p_var];
    DenseSubVector<Number> &Fu = *c.elem_subresiduals[u_var];
  
    DenseSubMatrix<Number> &Kvv = *c.elem_subjacobians[v_var][v_var];
    DenseSubMatrix<Number> &Kvp = *c.elem_subjacobians[v_var][p_var];
    DenseSubVector<Number> &Fv = *c.elem_subresiduals[v_var];
  
    DenseSubMatrix<Number> &KCu = *c.elem_subjacobians[C_var][u_var];
    DenseSubMatrix<Number> &KCv = *c.elem_subjacobians[C_var][v_var];
    DenseSubMatrix<Number> &KCC = *c.elem_subjacobians[C_var][C_var];
    DenseSubVector<Number> &FC = *c.elem_subresiduals[C_var];
  
    unsigned int n_qpoints = c.element_qrule->n_points();
  
    for (unsigned int qp=0; qp != n_qpoints; qp++)
      {
        Number p = c.interior_value(p_var, qp),
  	u = c.interior_value(u_var, qp),
  	v = c.interior_value(v_var, qp);
        Gradient grad_u = c.interior_gradient(u_var, qp),
  	grad_v = c.interior_gradient(v_var, qp),
  	grad_C = c.interior_gradient(C_var, qp);
  
        NumberVectorValue U     (u,     v);
        const Number C_x = grad_C(0);
        const Number C_y = grad_C(1);
  
        for (unsigned int i=0; i != n_u_dofs; i++)
          {
            Fu(i) += JxW[qp] *
                     (p*dphi[i][qp](0) -                // pressure term
  		    (grad_u*dphi[i][qp]));            // diffusion term
  
            Fv(i) += JxW[qp] *
                     (p*dphi[i][qp](1) -                // pressure term
  		    (grad_v*dphi[i][qp]));            // diffusion term
  
  	  FC(i) += JxW[qp] *
  	           ( (U*grad_C)*phi[i][qp] +                // convection term
  	            (1./Peclet)*(grad_C*dphi[i][qp]) );     // diffusion term
  
  
            if (request_jacobian && c.elem_solution_derivative)
              {
                libmesh_assert (c.elem_solution_derivative == 1.0);
  
                for (unsigned int j=0; j != n_u_dofs; j++)
                  {
                    Kuu(i,j) += JxW[qp] * (-(dphi[i][qp]*dphi[j][qp])); /* diffusion term  */
  
                    Kvv(i,j) += JxW[qp] * (-(dphi[i][qp]*dphi[j][qp])); /* diffusion term  */
  
  		  KCu(i,j) += JxW[qp]* ( (phi[j][qp]*C_x)*phi[i][qp] ); /* convection term */
  
  		  KCv(i,j) += JxW[qp]*( (phi[j][qp]*C_y)*phi[i][qp] );  /* convection term */
  
  		  KCC(i,j) += JxW[qp]*
  		              ( (U*dphi[j][qp])*phi[i][qp] +      /* nonlinear term (convection) */
  		              (1./Peclet)*(dphi[j][qp]*dphi[i][qp]) ); /* diffusion term */
  		}
  
  	      for (unsigned int j=0; j != n_p_dofs; j++)
  		{
  		  Kup(i,j) += JxW[qp]*psi[j][qp]*dphi[i][qp](0);
  		  Kvp(i,j) += JxW[qp]*psi[j][qp]*dphi[i][qp](1);
  		}
  	    }
  	}
  
      } // end of the quadrature point qp-loop
  
        return request_jacobian;
  }
  
  
  bool CoupledSystem::element_constraint (bool request_jacobian,
                                         DiffContext &context)
  {
    FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
  
  
    const std::vector<Real> &JxW = c.element_fe_var[u_var]->get_JxW();
  
    const std::vector<std::vector<RealGradient> >& dphi =
      c.element_fe_var[u_var]->get_dphi();
  
    const std::vector<std::vector<Real> >& psi =
      c.element_fe_var[p_var]->get_phi();
  
    const unsigned int n_u_dofs = c.dof_indices_var[u_var].size();
    const unsigned int n_p_dofs = c.dof_indices_var[p_var].size();
  
    DenseSubMatrix<Number> &Kpu = *c.elem_subjacobians[p_var][u_var];
    DenseSubMatrix<Number> &Kpv = *c.elem_subjacobians[p_var][v_var];
    DenseSubVector<Number> &Fp = *c.elem_subresiduals[p_var];
  
    unsigned int n_qpoints = c.element_qrule->n_points();
    for (unsigned int qp=0; qp != n_qpoints; qp++)
      {
        Gradient grad_u = c.interior_gradient(u_var, qp),
  	grad_v = c.interior_gradient(v_var, qp);
  
        for (unsigned int i=0; i != n_p_dofs; i++)
          {
            Fp(i) += JxW[qp] * psi[i][qp] *
                     (grad_u(0) + grad_v(1));
  
            if (request_jacobian && c.elem_solution_derivative)
              {
                libmesh_assert (c.elem_solution_derivative == 1.0);
  
                for (unsigned int j=0; j != n_u_dofs; j++)
                  {
                    Kpu(i,j) += JxW[qp]*psi[i][qp]*dphi[j][qp](0);
                    Kpv(i,j) += JxW[qp]*psi[i][qp]*dphi[j][qp](1);
                  }
              }
          }
      } // end of the quadrature point qp-loop
  
    return request_jacobian;
  }
  
  void CoupledSystem::postprocess()
  {
  
    computed_QoI = 0.0;
  
    computed_QoI = System::qoi[0];
  }
  
  Number CoupledFEMFunctionsx::operator()(const FEMContext& c, const Point& p,
  					const Real /* time */)
  {
    Number weight = 0.0;
  
    switch(var)
      {
      case 0:
        {
  	Gradient grad_C = c.point_gradient(3, p);
  
  	weight = grad_C(0);
        }
        break;
  
      case 3:
        {
  	Number u = c.point_value(0, p);
  
  	weight = u;
        }
        break;
  
      default:
        {
  	std::cout<<"Wrong variable number"<<var<<" passed to CoupledFEMFunctionsx object ! Quitting !"<<std::endl;
  	libmesh_error();
        }
  
      }
  
    return weight;
  }
  
  Number CoupledFEMFunctionsy::operator()(const FEMContext& c, const Point& p,
  					const Real /* time */)
  {
    Number weight = 0.0;
  
    switch(var)
      {
      case 1:
        {
  	Gradient grad_C = c.point_gradient(3, p);
  
  	weight = grad_C(1);
        }
        break;
  
      case 3:
        {
  	Number v = c.point_value(1, p);
  
  	weight = v;
        }
        break;
  
      default:
        {
  	std::cout<<"Wrong variable number "<<var<<" passed to CoupledFEMFunctionsy object ! Quitting !"<<std::endl;
  	libmesh_error();
        }
      }
  
    return weight;
  }



The source file domain.C without comments:

 
  #include "libmesh/mesh.h"
  #include "libmesh/serial_mesh.h"
  #include "libmesh/mesh_generation.h"
  #include "libmesh/mesh_modification.h"
  #include "libmesh/mesh_refinement.h"
  
  #include "domain.h"
  
  using namespace libMesh;
  
  
  
  void build_domain (Mesh &mesh, FEMParameters &param)
  {
    mesh.read(param.domainfile);
  
    std::cout<<"Making elements 2nd order"<<std::endl;
  
    mesh.all_second_order();
  
  }



The source file femparameters.C without comments:

 
  #include "femparameters.h"
  
  using namespace libMesh;
  
  #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(initial_timestep);
      GETPOT_INT_INPUT(n_timesteps);
      GETPOT_INPUT(transient);
      GETPOT_INT_INPUT(deltat_reductions);
      GETPOT_INPUT(timesolver_core);
      GETPOT_INPUT(end_time);
      GETPOT_INPUT(deltat);
      GETPOT_INPUT(timesolver_theta);
      GETPOT_INPUT(timesolver_maxgrowth);
      GETPOT_INPUT(timesolver_upper_tolerance);
      GETPOT_INPUT(timesolver_tolerance);
      GETPOT_INPUT(steadystate_tolerance);
  
      GETPOT_REGISTER(timesolver_norm);
      const unsigned int n_timesolver_norm = input.vector_variable_size("timesolver_norm");
      timesolver_norm.resize(n_timesolver_norm, L2);
      for (unsigned int i=0; i != n_timesolver_norm; ++i)
        {
          int current_norm = 0; // L2
          if (timesolver_norm[i] == H1)
            current_norm = 1;
          if (timesolver_norm[i] == H2)
            current_norm = 2;
          current_norm = input("timesolver_norm", current_norm, i);
          if (current_norm == 0)
            timesolver_norm[i] = L2;
          else if (current_norm == 1)
            timesolver_norm[i] = H1;
          else if (current_norm == 2)
            timesolver_norm[i] = H2;
          else
            timesolver_norm[i] = DISCRETE_L2;
        }
  
      GETPOT_INT_INPUT(dimension);
      GETPOT_INPUT(domaintype);
      GETPOT_INPUT(domainfile);
      GETPOT_INPUT(elementtype);
      GETPOT_INPUT(fine_mesh_file_primal);
      GETPOT_INPUT(fine_mesh_soln_primal);
      GETPOT_INPUT(fine_mesh_file_adjoint);
      GETPOT_INPUT(fine_mesh_soln_adjoint);
      GETPOT_INPUT(elementorder);
      GETPOT_INPUT(domain_xmin);
      GETPOT_INPUT(domain_ymin);
      GETPOT_INPUT(domain_zmin);
      GETPOT_INPUT(domain_edge_width);
      GETPOT_INPUT(domain_edge_length);
      GETPOT_INPUT(domain_edge_height);
      GETPOT_INT_INPUT(coarsegridx);
      GETPOT_INT_INPUT(coarsegridy);
      GETPOT_INT_INPUT(coarsegridz);
      GETPOT_INT_INPUT(coarserefinements);
      GETPOT_INT_INPUT(coarsecoarsenings);
      GETPOT_INT_INPUT(extrarefinements);
      GETPOT_INPUT(use_petsc_snes);
      GETPOT_INPUT(time_solver_quiet);
      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_INT_INPUT(initial_adaptivesteps);
      GETPOT_INT_INPUT(initial_sobolev_order);
      GETPOT_INT_INPUT(initial_extra_quadrature);
      GETPOT_INPUT(refine_uniformly);
      GETPOT_INPUT(indicator_type);
      GETPOT_INPUT(adjoint_residual_type);
      GETPOT_INPUT(patch_reuse);
      GETPOT_INT_INPUT(sobolev_order);
      GETPOT_INPUT(alternate_with_uniform_steps);
      GETPOT_INPUT(alternate_step_number);
      GETPOT_INPUT(component_wise_error);
      GETPOT_INPUT(compute_sensitivities);
      GETPOT_INPUT(compare_to_fine_solution_primal);
      GETPOT_INPUT(compare_to_fine_solution_adjoint);
      GETPOT_INPUT(do_forward_sensitivity);
      GETPOT_INPUT(do_adjoint_sensitivity);
      GETPOT_INPUT(postprocess_adjoint);
      GETPOT_INT_INPUT(write_interval);
      GETPOT_INPUT(output_xda);
      GETPOT_INPUT(output_xdr);
  
      GETPOT_INPUT(output_gz);
  #ifndef LIBMESH_HAVE_GZSTREAM
      output_gz                   = false;
  #endif
  
      GETPOT_INPUT(output_bz2);
  #ifndef LIBMESH_HAVE_BZ2
      output_bz2                  = false;
  #endif
  
      GETPOT_INPUT(output_gmv);
      GETPOT_INPUT(write_gmv_error);
  
  #ifndef LIBMESH_HAVE_GMV
      output_gmv                  = false;
      write_gmv_error             = false;
  #endif
  
      GETPOT_INPUT(output_tecplot);
      GETPOT_INPUT(write_tecplot_error);
  #ifndef LIBMESH_HAVE_TECPLOT_API
      output_tecplot              = false;
      write_tecplot_error         = false;
  #endif
  
      GETPOT_INPUT(run_simulation);
      GETPOT_INPUT(run_postprocess);
  
      run_simulation              = input("run_simulation", run_simulation);
      run_postprocess             = input("run_postprocess", run_postprocess);
  
      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(extra_quadrature_order);
      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);
  
  
  
  }



The source file H-qoi.C without comments:

 
  #include "H-qoi.h"
  
  
  using namespace libMesh;
  
  void CoupledSystemQoI::init_qoi( std::vector<Number>& sys_qoi)
  {
    sys_qoi.resize(1);
    return;
  }
  
  void CoupledSystemQoI::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<Real> > &phi = c.side_fe_var[0]->get_phi();
  
    const std::vector<Point > &q_point = c.side_fe_var[0]->get_xyz();
  
    const unsigned int n_u_dofs = c.dof_indices_var[1].size();
  
    DenseSubVector<Number> &Qu = *c.elem_qoi_subderivatives[0][0];
    DenseSubVector<Number> &QC = *c.elem_qoi_subderivatives[0][3];
  
    unsigned int n_qpoints = c.side_qrule->n_points();
  
    Number u = 0. ;
    Number C = 0. ;
  
    if(c.has_side_boundary_id(2))
      {
        for (unsigned int qp=0; qp != n_qpoints; qp++)
    	{
  	  Real x = q_point[qp](0);
  
  	  if(x < 0.)
  	    {
  	      u = c.side_value(0,qp);
  	      C = c.side_value(3,qp);
  
  	      for (unsigned int i=0; i != n_u_dofs; i++)
  		{
  		  Qu(i) += JxW[qp] * -phi[i][qp] * C;
  		  QC(i) += JxW[qp] * phi[i][qp] * -u;
  		}
  	    } // end if
  
    	} // end quadrature loop
  
      } // end if on outlet
  }
  
  void CoupledSystemQoI::side_qoi(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<Point > &q_point = c.side_fe_var[0]->get_xyz();
  
  
    unsigned int n_qpoints = c.side_qrule->n_points();
  
    Number dQoI_0 = 0. ;
    Number u = 0. ;
    Number C = 0. ;
  
    if(c.has_side_boundary_id(2))
      {
        for (unsigned int qp=0; qp != n_qpoints; qp++)
  	{
  	  Real x = q_point[qp](0);
  
  	  if(x < 0.)
  	    {
  	      u = c.side_value(0,qp);
  	      C = c.side_value(3,qp);
  
  	      dQoI_0 += JxW[qp] * -u * C;
  	    } // end if
  
    	} // end quadrature loop
  
      } // end if on bdry
  
    c.elem_qoi[0] += dQoI_0;
  
  }



The source file initial.C without comments:

 
  #include "initial.h"
  #include "femparameters.h"
  
  using namespace libMesh;
  
  void read_initial_parameters()
  {
  }
  
  void finish_initialization()
  {
  }
  
  
  
  Number initial_value(const Point& /* p */,
                       const Parameters& /* param */,
                       const std::string&,
                       const std::string&)
  {
  
    return 1.;
  
  }
  
  
  
  Gradient initial_grad(const Point& /* p */,
                        const Parameters& /* param */,
                        const std::string&,
                        const std::string&)
  {
    return 0.;
  }



The source file output.C without comments:

 
  #include <fstream>
  #include <unistd.h>
  
  #include "libmesh/libmesh_common.h"
  
  #include "output.h"
  
  using namespace libMesh;
  
  void start_output(unsigned int timesteps,
                    std::string filename,
                    std::string varname)
  {
    if (timesteps)
      {
        std::ifstream infile (filename.c_str());
        char buf[1025];
        if (!infile.is_open())
          libmesh_error();
  
        for (unsigned int i=0; i != timesteps+1; ++i)
          infile.getline(buf, 1024);
        if (!infile.good())
          libmesh_error();
        unsigned int length = infile.tellg();
        infile.close();
  
        int err = truncate(filename.c_str(), length);
        if (err != 0)
          libmesh_error();
      }
    else
      {
        std::ofstream outfile (filename.c_str(), std::ios_base::out);
        outfile << varname << " = [" << std::endl;
      }
  }



The console output of the program:

make[4]: Entering directory `/net/spark/workspace/roystgnr/libmesh/git/devel/examples/adjoints/adjoints_ex3'
***************************************************************
* Running Example adjoints_ex3:
*  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_ex3/.libs/lt-example-devel
Building mesh
Making elements 2nd order
Building system
Initializing systems
 Mesh Information:
  mesh_dimension()=2
  spatial_dimension()=3
  n_nodes()=99
    n_local_nodes()=27
  n_elem()=16
    n_local_elem()=4
    n_active_elem()=16
  n_subdomains()=1
  n_partitions()=4
  n_processors()=4
  n_threads()=1
  processor_id()=0

 EquationSystems
  n_systems()=1
   System #0, "CoupledSystem"
    Type "Implicit"
    Variables={ "u" "v" } "p" "C" 
    Finite Element Types="LAGRANGE", "JACOBI_20_00" "LAGRANGE", "JACOBI_20_00" "LAGRANGE", "JACOBI_20_00" 
    Infinite Element Mapping="CARTESIAN" "CARTESIAN" "CARTESIAN" 
    Approximation Orders="SECOND", "THIRD" "FIRST", "THIRD" "SECOND", "THIRD" 
    n_dofs()=331
    n_local_dofs()=91
    n_constrained_dofs()=138
    n_local_constrained_dofs()=34
    n_vectors()=1
    n_matrices()=1
    DofMap Sparsity
      Average  On-Processor Bandwidth <= 36.7855
      Average Off-Processor Bandwidth <= 5.91541
      Maximum  On-Processor Bandwidth <= 62
      Maximum Off-Processor Bandwidth <= 40
    DofMap Constraints
      Number of DoF Constraints = 138
      Number of Heterogenous Constraints= 5
      Average DoF Constraint Length= 0
      Number of Node Constraints = 0

Starting adaptive loop

Adaptive step 0
Solving the forward problem
We have 16 active elements and 193 active dofs.

  Nonlinear solver converged, step 1, residual reduction 9.03212e-14 < 1e-05
The boundary QoI is 0.083342124064158307

Solving the adjoint problem
Adaptive step 1
Solving the forward problem
We have 22 active elements and 259 active dofs.

  Nonlinear solver converged, step 0, residual reduction 5.6868816087085118e-07 < 1.0000000000000001e-05
The boundary QoI is 0.082286987672833115

Solving the adjoint problem
Adaptive step 2
Solving the forward problem
We have 31 active elements and 363 active dofs.

  Nonlinear solver converged, step 1, residual reduction 3.8583229881642115e-14 < 1.0000000000000001e-05
  Nonlinear solver converged, step 1, relative step size 7.5013647980960183e-07 < 1.0000000000000001e-05
The boundary QoI is 0.083169669189009962

Solving the adjoint problem
Adaptive step 3
Solving the forward problem
We have 43 active elements and 522 active dofs.

  Nonlinear solver converged, step 1, residual reduction 4.3999517315436118e-14 < 1.0000000000000001e-05
  Nonlinear solver converged, step 1, relative step size 6.565604269298289e-06 < 1.0000000000000001e-05
The boundary QoI is 0.083339483627459909

Solving the adjoint problem
Adaptive step 4
Solving the forward problem
We have 58 active elements and 671 active dofs.

  Nonlinear solver converged, step 0, residual reduction 3.9242900634249945e-06 < 1.0000000000000001e-05
The boundary QoI is 0.083044585890271971

Solving the adjoint problem
Adaptive step 5
Solving the forward problem
We have 76 active elements and 898 active dofs.

  Nonlinear solver converged, step 1, residual reduction 5.9607546528653657e-14 < 1.0000000000000001e-05
  Nonlinear solver converged, step 1, relative step size 2.0501238293058485e-07 < 1.0000000000000001e-05
The boundary QoI is 0.08333916169049338


 -------------------------------------------------------------------------------------------------------------------
| Processor id:   0                                                                                                 |
| Num Processors: 4                                                                                                 |
| Time:           Fri Apr 19 11:47:50 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.09678, Active time=1.97643                                                        |
 ---------------------------------------------------------------------------------------------------------------------
| 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   |
|---------------------------------------------------------------------------------------------------------------------|
|                                                                                                                     |
|                                                                                                                     |
| AdjointResidualErrorEstimator                                                                                       |
|   estimate_error()                      20        0.0014      0.000069    0.9382      0.046909    0.07     47.47    |
|                                                                                                                     |
| DofMap                                                                                                              |
|   add_neighbors_to_send_list()          31        0.0024      0.000079    0.0046      0.000149    0.12     0.23     |
|   build_constraint_matrix()             298       0.0030      0.000010    0.0030      0.000010    0.15     0.15     |
|   build_sparsity()                      6         0.0024      0.000406    0.0075      0.001243    0.12     0.38     |
|   cnstrn_elem_mat_vec()                 105       0.0072      0.000068    0.0072      0.000068    0.36     0.36     |
|   constrain_elem_matrix()               44        0.0013      0.000030    0.0013      0.000030    0.07     0.07     |
|   constrain_elem_vector()               149       0.0007      0.000005    0.0007      0.000005    0.04     0.04     |
|   create_dof_constraints()              31        0.0156      0.000503    0.0607      0.001959    0.79     3.07     |
|   distribute_dofs()                     31        0.0056      0.000180    0.0288      0.000928    0.28     1.46     |
|   dof_indices()                         25520     0.2032      0.000008    0.2032      0.000008    10.28    10.28    |
|   enforce_constraints_exactly()         37        0.0068      0.000184    0.0068      0.000184    0.34     0.34     |
|   old_dof_indices()                     590       0.0067      0.000011    0.0067      0.000011    0.34     0.34     |
|   prepare_send_list()                   31        0.0001      0.000003    0.0001      0.000003    0.01     0.01     |
|   reinit()                              31        0.0060      0.000193    0.0060      0.000193    0.30     0.30     |
|                                                                                                                     |
| EquationSystems                                                                                                     |
|   build_discontinuous_solution_vector() 25        0.0019      0.000078    0.0033      0.000134    0.10     0.17     |
|   build_solution_vector()               11        0.0012      0.000113    0.0099      0.000899    0.06     0.50     |
|                                                                                                                     |
| FE                                                                                                                  |
|   compute_shape_functions()             21372     0.1065      0.000005    0.1065      0.000005    5.39     5.39     |
|   init_shape_functions()                5032      0.0401      0.000008    0.0401      0.000008    2.03     2.03     |
|   inverse_map()                         10827     0.0364      0.000003    0.0364      0.000003    1.84     1.84     |
|                                                                                                                     |
| FEMSystem                                                                                                           |
|   assemble_qoi()                        6         0.0031      0.000518    0.0486      0.008096    0.16     2.46     |
|   assemble_qoi_derivative()             5         0.0266      0.005329    0.0438      0.008754    1.35     2.21     |
|   assembly()                            10        0.0316      0.003164    0.0721      0.007211    1.60     3.65     |
|   assembly(get_jacobian)                5         0.0132      0.002642    0.0295      0.005898    0.67     1.49     |
|   assembly(get_residual)                10        0.0094      0.000940    0.0430      0.004302    0.48     2.18     |
|                                                                                                                     |
| FEMap                                                                                                               |
|   compute_affine_map()                  21372     0.0648      0.000003    0.0648      0.000003    3.28     3.28     |
|   compute_face_map()                    1896      0.0142      0.000007    0.0399      0.000021    0.72     2.02     |
|   init_face_shape_functions()           812       0.0013      0.000002    0.0013      0.000002    0.06     0.06     |
|   init_reference_to_physical_map()      5032      0.0538      0.000011    0.0538      0.000011    2.72     2.72     |
|                                                                                                                     |
| GMVIO                                                                                                               |
|   write_nodal_data()                    11        0.0138      0.001252    0.0145      0.001321    0.70     0.74     |
|                                                                                                                     |
| ImplicitSystem                                                                                                      |
|   adjoint_solve()                       5         0.0027      0.000547    0.1677      0.033542    0.14     8.49     |
|                                                                                                                     |
| LocationMap                                                                                                         |
|   find()                                400       0.0005      0.000001    0.0005      0.000001    0.02     0.02     |
|   init()                                10        0.0006      0.000062    0.0006      0.000062    0.03     0.03     |
|                                                                                                                     |
| Mesh                                                                                                                |
|   all_first_order()                     25        0.0016      0.000064    0.0016      0.000064    0.08     0.08     |
|   all_second_order()                    1         0.0001      0.000142    0.0001      0.000142    0.01     0.01     |
|   contract()                            5         0.0000      0.000009    0.0001      0.000024    0.00     0.01     |
|   find_neighbors()                      57        0.0050      0.000088    0.0121      0.000212    0.25     0.61     |
|   renumber_nodes_and_elem()             69        0.0005      0.000007    0.0005      0.000007    0.03     0.03     |
|                                                                                                                     |
| MeshCommunication                                                                                                   |
|   compute_hilbert_indices()             32        0.0030      0.000093    0.0030      0.000093    0.15     0.15     |
|   find_global_indices()                 32        0.0011      0.000034    0.0125      0.000391    0.06     0.63     |
|   parallel_sort()                       32        0.0021      0.000065    0.0054      0.000168    0.10     0.27     |
|                                                                                                                     |
| MeshOutput                                                                                                          |
|   write_equation_systems()              11        0.0003      0.000023    0.0255      0.002315    0.01     1.29     |
|                                                                                                                     |
| MeshRefinement                                                                                                      |
|   _coarsen_elements()                   10        0.0001      0.000005    0.0002      0.000015    0.00     0.01     |
|   _refine_elements()                    10        0.0007      0.000070    0.0033      0.000330    0.04     0.17     |
|   add_point()                           400       0.0006      0.000002    0.0012      0.000003    0.03     0.06     |
|   make_coarsening_compatible()          20        0.0005      0.000027    0.0005      0.000027    0.03     0.03     |
|   make_flags_parallel_consistent()      15        0.0008      0.000055    0.0055      0.000367    0.04     0.28     |
|   make_refinement_compatible()          20        0.0001      0.000004    0.0002      0.000011    0.00     0.01     |
|                                                                                                                     |
| MetisPartitioner                                                                                                    |
|   partition()                           32        0.0107      0.000333    0.0273      0.000854    0.54     1.38     |
|                                                                                                                     |
| NewtonSolver                                                                                                        |
|   solve()                               6         0.0940      0.015661    0.5123      0.085382    4.75     25.92    |
|                                                                                                                     |
| Parallel                                                                                                            |
|   allgather()                           169       0.0067      0.000039    0.0071      0.000042    0.34     0.36     |
|   broadcast()                           22        0.0001      0.000003    0.0001      0.000003    0.00     0.00     |
|   max(bool)                             52        0.0018      0.000036    0.0018      0.000036    0.09     0.09     |
|   max(scalar)                           3484      0.0150      0.000004    0.0150      0.000004    0.76     0.76     |
|   max(vector)                           787       0.0047      0.000006    0.0143      0.000018    0.24     0.72     |
|   max(vector)                     10        0.0001      0.000007    0.0003      0.000029    0.00     0.01     |
|   min(bool)                             4177      0.0168      0.000004    0.0168      0.000004    0.85     0.85     |
|   min(scalar)                           3367      0.4617      0.000137    0.4617      0.000137    23.36    23.36    |
|   min(vector)                           792       0.0052      0.000007    0.0154      0.000019    0.26     0.78     |
|   probe()                               1128      0.0072      0.000006    0.0072      0.000006    0.37     0.37     |
|   receive()                             1128      0.0023      0.000002    0.0096      0.000009    0.11     0.49     |
|   send()                                1128      0.0011      0.000001    0.0011      0.000001    0.06     0.06     |
|   send_receive()                        1192      0.0033      0.000003    0.0151      0.000013    0.17     0.76     |
|   sum()                                 354       0.0050      0.000014    0.3842      0.001085    0.25     19.44    |
|                                                                                                                     |
| Parallel::Request                                                                                                   |
|   wait()                                1128      0.0008      0.000001    0.0008      0.000001    0.04     0.04     |
|                                                                                                                     |
| Partitioner                                                                                                         |
|   set_node_processor_ids()              58        0.0035      0.000060    0.0195      0.000337    0.18     0.99     |
|   set_parent_processor_ids()            32        0.0003      0.000010    0.0003      0.000010    0.02     0.02     |
|                                                                                                                     |
| PatchRecoveryErrorEstimator                                                                                         |
|   estimate_error()                      120       0.1552      0.001293    0.5901      0.004917    7.85     29.86    |
|                                                                                                                     |
| PetscLinearSolver                                                                                                   |
|   solve()                               15        0.3861      0.025737    0.3861      0.025737    19.53    19.53    |
|                                                                                                                     |
| ProjectVector                                                                                                       |
|   operator()                            10        0.0028      0.000280    0.0135      0.001348    0.14     0.68     |
|                                                                                                                     |
| System                                                                                                              |
|   project_vector()                      10        0.0122      0.001222    0.0321      0.003208    0.62     1.62     |
|                                                                                                                     |
| WeightedPatchRecoveryErrorEstimator                                                                                 |
|   estimate_error()                      40        0.0787      0.001968    0.3467      0.008667    3.98     17.54    |
|                                                                                                                     |
| XdrIO                                                                                                               |
|   read()                                1         0.0003      0.000273    0.0003      0.000305    0.01     0.02     |
 ---------------------------------------------------------------------------------------------------------------------
| Totals:                                 113746    1.9764                                          100.00            |
 ---------------------------------------------------------------------------------------------------------------------

 
***************************************************************
* Done Running Example adjoints_ex3:
*  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_ex3'

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

Hosted By:
SourceForge.net Logo