#ifndef _verify_solution_hpp_
#define _verify_solution_hpp_

//@HEADER
// ************************************************************************
//
// MiniFE: Simple Finite Element Assembly and Solve
// Copyright (2006-2013) Sandia Corporation
//
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
//
// This library is free software; you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; either version 2.1 of the
// License, or (at your option) any later version.
//
// This library 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
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
// USA
//
// ************************************************************************
//@HEADER

#include <sstream>
#include <stdexcept>
#include <map>
#include <algorithm>

#include <simple_mesh_description.hpp>
#include <analytic_soln.hpp>
#include <box_utils.hpp>
#include <utils.hpp>

#ifdef HAVE_MPI
#include <mpi.h>
#endif

namespace miniFE {

template<typename Scalar>
struct err_info {
  Scalar err;
  Scalar computed;
  Scalar analytic;
  Scalar coords[3];
};

template<typename VectorType>
int
verify_solution(const simple_mesh_description<typename VectorType::GlobalOrdinalType>& mesh,
                const VectorType& x, double tolerance, bool verify_whole_domain = false)
{
  typedef typename VectorType::GlobalOrdinalType GlobalOrdinal;
  typedef typename VectorType::ScalarType Scalar;

  int global_nodes_x = mesh.global_box[0][1]+1;
  int global_nodes_y = mesh.global_box[1][1]+1;
  int global_nodes_z = mesh.global_box[2][1]+1;
  Box box;
  copy_box(mesh.local_box, box);

  //num-owned-nodes in each dimension is num-elems+1
  //only if num-elems > 0 in that dimension *and*
  //we are at the high end of the global range in that dimension:
  if (box[0][1] > box[0][0] && box[0][1] == mesh.global_box[0][1]) ++box[0][1];
  if (box[1][1] > box[1][0] && box[1][1] == mesh.global_box[1][1]) ++box[1][1];
  if (box[2][1] > box[2][0] && box[2][1] == mesh.global_box[2][1]) ++box[2][1];

  std::vector<GlobalOrdinal> rows;
  std::vector<Scalar> row_coords;

  int roffset = 0;
  for(int iz=box[2][0]; iz<box[2][1]; ++iz) {
   for(int iy=box[1][0]; iy<box[1][1]; ++iy) {
    for(int ix=box[0][0]; ix<box[0][1]; ++ix) {
      GlobalOrdinal row_id =
          get_id<GlobalOrdinal>(global_nodes_x, global_nodes_y, global_nodes_z,
                                ix, iy, iz);
      Scalar x, y, z;
      get_coords(row_id, global_nodes_x, global_nodes_y, global_nodes_z, x, y, z);

      bool verify_this_point = false;
      if (verify_whole_domain) verify_this_point = true;
      else if (std::abs(x - 0.5) < 0.05 && std::abs(y - 0.5) < 0.05 && std::abs(z - 0.5) < 0.05) {
        verify_this_point = true;
      }

      if (verify_this_point) {
        rows.push_back(roffset);
        row_coords.push_back(x);
        row_coords.push_back(y);
        row_coords.push_back(z);
      }

      ++roffset;
    }
   }
  }

  int return_code = 0;

  const int num_terms = 300;

  err_info<Scalar> max_error;
  max_error.err = 0.0;

  for(size_t i=0; i<rows.size(); ++i) {
    Scalar computed_soln = x.coefs[rows[i]];
    Scalar x = row_coords[i*3];
    Scalar y = row_coords[i*3+1];
    Scalar z = row_coords[i*3+2];
    Scalar analytic_soln = 0.0;
    //set exact boundary-conditions:
    if (x == 1.0) {
      //x==1 is first, we want soln to be 1 even around the edges
      //of the x==1 plane where y and/or z may be 0 or 1...
      analytic_soln = 1;
    }
    else if (x == 0.0 || y == 0.0 || z == 0.0) {
      analytic_soln = 0;
    }
    else if (y == 1.0 || z == 1.0) {
      analytic_soln = 0;
    }
    else {
      analytic_soln = soln(x, y, z, num_terms, num_terms);
    }

#ifdef MINIFE_DEBUG_VERBOSE
std::cout<<"("<<x<<","<<y<<","<<z<<") row "<<rows[i]<<": computed: "<<computed_soln<<",  analytic: "<<analytic_soln<<std::endl;
#endif
    Scalar err = std::abs(analytic_soln - computed_soln);
    if (err > max_error.err) {
      max_error.err = err;
      max_error.computed = computed_soln;
      max_error.analytic = analytic_soln;
      max_error.coords[0] = x;
      max_error.coords[1] = y;
      max_error.coords[2] = z;
    }
  }

  Scalar local_max_err = max_error.err;
  Scalar global_max_err = 0;
#ifdef HAVE_MPI
  MPI_Allreduce(&local_max_err, &global_max_err, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
#else
  global_max_err = local_max_err;
#endif

  if (local_max_err == global_max_err) {
    if (max_error.err > tolerance) {
      std::cout << "max absolute error is "<<max_error.err<<":"<<std::endl;
      std::cout << "   at position ("<<max_error.coords[0]<<","<<max_error.coords[1]<<","<<max_error.coords[2]<<"), "<<std::endl;
      std::cout << "   computed solution: "<<max_error.computed<<",  analytic solution: "<<max_error.analytic<<std::endl;
    }
    else {
      std::cout << "solution matches analytic solution to within "<<tolerance<<" or better."<<std::endl;
    }
  }

  if (global_max_err > tolerance) {
    return_code = 1;
  }

  return return_code;
}

}//namespace miniFE

#endif