/*
  Copyright (C) 2011, 2012 by the authors of the ASPECT code.

  This file is part of ASPECT.

  ASPECT is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; either version 2, or (at your option)
  any later version.

  ASPECT is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with ASPECT; see the file doc/COPYING.  If not see
  <http://www.gnu.org/licenses/>.
*/


#include <aspect/postprocess/energetics.h>
#include <aspect/simulator_access.h>

#include <deal.II/base/quadrature_lib.h>
#include <deal.II/fe/fe_values.h>


namespace aspect
{
  namespace Postprocess
  {
    template <int dim>
    std::pair<std::string,std::string>
    Energetics<dim>::execute (TableHandler &statistics)
    {
      // create a quadrature formula based on the temperature element alone.
      // be defensive about determining that what we think is the temperature
      // element is it in fact
      Assert (this->get_fe().n_base_elements() == 3+(this->n_compositional_fields()>0 ? 1 : 0),
              ExcNotImplemented());
      const QGauss<dim> quadrature_formula (this->get_fe().base_element(0).degree+1);
      const unsigned int n_q_points = quadrature_formula.size();

      FEValues<dim> fe_values (this->get_mapping(),
                               this->get_fe(),
                               quadrature_formula,
                               update_values   |  update_gradients | 
                               update_quadrature_points |
                               update_JxW_values);

      std::vector<double> temperature_values(n_q_points);
      std::vector<double> pressure_values(n_q_points);
      std::vector<Tensor<1,dim> > velocity_values(n_q_points);
      std::vector<SymmetricTensor<2,dim> > strain_rates(n_q_points);
      std::vector<Point<dim> > quadrature_points(n_q_points);

      std::vector<std::vector<double> > composition_values(this->n_compositional_fields());
      for ( unsigned int c=0; c<this->n_compositional_fields(); ++c) composition_values[c].resize(n_q_points, 0.0);

      typename DoFHandler<dim>::active_cell_iterator
      cell = this->get_dof_handler().begin_active(),
      endc = this->get_dof_handler().end();

      double local_adiabatic_heating = 0;
      double local_friction_heating = 0;

      // compute the integral quantities by quadrature
      for (; cell!=endc; ++cell)
        if (cell->is_locally_owned())
          {
            //get the relevant quantities from the solution vector
            fe_values.reinit (cell);
            fe_values[this->introspection().extractors.temperature].get_function_values (this->get_solution(), temperature_values);
            fe_values[this->introspection().extractors.pressure].get_function_values (this->get_solution(), pressure_values);
            fe_values[this->introspection().extractors.velocities].get_function_values (this->get_solution(), velocity_values);
            fe_values[this->introspection().extractors.velocities].get_function_symmetric_gradients (this->get_solution(), strain_rates);
            for (unsigned int c=0; c<this->n_compositional_fields(); ++c)
            {
                fe_values[this->introspection().extractors.compositional_fields[c]].get_function_values (this->get_solution(),
                    composition_values[c]);
            }
            
            //loop over quadrature pts
            quadrature_points = fe_values.get_quadrature_points();
            for (unsigned int q=0; q<n_q_points; ++q)
            {
              std::vector<double> composition(this->n_compositional_fields());
              for (unsigned int c=0; c<this->n_compositional_fields(); ++c)
                composition[c] = composition_values[c][q];

              Tensor<1,dim> grav = this->get_gravity_model().gravity_vector(quadrature_points[q]);
              double viscosity =   this->get_material_model().viscosity(temperature_values[q],
                                                                        pressure_values[q],
                                                                        composition,
                                                                        strain_rates[q],
                                                                        quadrature_points[q]);
              double alpha =      this->get_material_model().thermal_expansion_coefficient(temperature_values[q],
                                                                        pressure_values[q],
                                                                        composition,
                                                                        quadrature_points[q]);
              double density =    this->get_material_model().density(temperature_values[q],
                                                                        pressure_values[q],
                                                                        composition,
                                                                        quadrature_points[q]);

              local_friction_heating  += 2.0 * viscosity * (strain_rates[q] * strain_rates[q]) * fe_values.JxW(q);

              
              //Formula derived from momentum equation without subtracting reference state
              //local_adiabatic_heating += (velocity_values[q] * grav) * density * fe_values.JxW(q);

              //Formula derived from momentum equation with subtracting reference state
              //local_adiabatic_heating += -(velocity_values[q] * grav) * alpha * this->get_material_model().reference_density() * temperature_values[q] * fe_values.JxW(q);

              //This should be, I think, the appropriate equation.
              local_adiabatic_heating += (velocity_values[q] * grav) * (density - this->get_material_model().reference_density()) * fe_values.JxW(q);

              //Formula used in assembly.cc
              //local_adiabatic_heating += -(velocity_values[q] * grav) * alpha * density * temperature_values[q] * fe_values.JxW(q);
 
            }
          }

      //get the global values
      const double global_friction_heating
        = Utilities::MPI::sum (local_friction_heating, this->get_mpi_communicator());
      const double global_adiabatic_heating
        = Utilities::MPI::sum (local_adiabatic_heating, this->get_mpi_communicator());

      statistics.add_value ("Adiabatic heating",
                            global_adiabatic_heating);
      statistics.add_value ("Frictional heating",
                            global_friction_heating );

      {
        const char *columns[] = { "Adiabatic heating",
                                  "Frictional heating"
                                };
        for (unsigned int i=0; i<sizeof(columns)/sizeof(columns[0]); ++i)
          {
            statistics.set_precision (columns[i], 8);
            statistics.set_scientific (columns[i], true);
          }
      }

      std::ostringstream output;
      output.precision(4);
      output << global_friction_heating << "\t" 
             << global_adiabatic_heating <<"\t"<<1.0 - global_friction_heating/global_adiabatic_heating ;

      return std::pair<std::string, std::string> ("Energetics: friction/adiabatic/error:",
                                                  output.str());
    }
  }
}


// explicit instantiations
namespace aspect
{
  namespace Postprocess
  {
    ASPECT_REGISTER_POSTPROCESSOR(Energetics,
                                  "energetics",
                                  "A postprocessor that computes some statistics about "
                                  "the temperature field.")
  }
}