#ifndef _analytic_soln_hpp_
#define _analytic_soln_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 <cmath>

#ifndef MINIFE_SCALAR
#define MINIFE_SCALAR double;
#endif

namespace miniFE {

typedef MINIFE_SCALAR Scalar;

// The 'soln' function below computes the analytic solution for
// steady state temperature in a brick-shaped domain (formally called
// a rectangular parallelepiped). The inputs to the function are
// the x,y,z coordinates of the point at which temperature is to be
// computed, and the number of terms p,q in the series expansion.
//
// The equations used for the temperature solution are equations 9 and 10
// in section 6.2 of Carslaw & Jaeger, "Conduction of Heat in Solids".
//
// The paralellepiped being used is defined by this domain:
// 0 <= x <= 1.0
// 0 <= y <= 1.0
// 0 <= z <= 1.0
//
// With boundary conditions prescribing the temperature to be 1.0 on
// the x==1.0 face, and 0.0 on all other faces.
//
// Thus, in the equations from Carslaw & Jaeger, the following constants
// are used:
//
// a == b == c == 1.0  (the extents of the domain)
// v1 == 0.0           (temperature at x == 0.0)
// v2 == 1.0           (temperature at x == 1.0)
//

const Scalar PI = 3.141592653589793238462;
const Scalar PI_SQR = PI*PI;
const Scalar term0 = 16.0/(PI_SQR);

inline Scalar fcn_l(int p, int q)
{
  return std::sqrt((2*p+1)*(2*p+1)*PI_SQR + (2*q+1)*(2*q+1)*PI_SQR);
}

inline Scalar fcn(int n, Scalar u)
{
  return (2*n+1)*PI*u;
}

inline Scalar soln(Scalar x, Scalar y, Scalar z, int max_p, int max_q)
{
  Scalar sum = 0;
  for(int p=0; p<=max_p; ++p) {
    const Scalar p21y = fcn(p, y);
    const Scalar sin_py = std::sin(p21y)/(2*p+1);
    for(int q=0; q<=max_q; ++q) {
      const Scalar q21z = fcn(q, z);
      const Scalar sin_qz = std::sin(q21z)/(2*q+1);

      const Scalar l = fcn_l(p, q);

      const Scalar sinh1 = std::sinh(l*x);
      const Scalar sinh2 = std::sinh(l);

      const Scalar tmp = (sinh1*sin_py)*(sin_qz/sinh2);

      //if the scalar l gets too big, sinh(l) becomes inf.
      //if that happens, tmp is a NaN.
      //crude check for NaN:
      //if tmp != tmp, tmp is NaN
      if (tmp == tmp) {
        sum += tmp;
      }
      else {
        //if we got a NaN, break out of this inner loop and go to
        //the next iteration of the outer loop.
        break;
      }
    }
  }
  return term0*sum;
}

}//namespace miniFE

#endif /* _analytic_soln_hpp_ */