/*
XLiFE++ is an extended library of finite elements written in C++
    Copyright (C) 2014  Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Yvon; Nguyen, Manh-Ha; Chambeyron, Colin

    This program 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 3 of the License, or
    (at your option) any later version.
    This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.
 */

/*!
	\file unit_operator.cpp
	\author E. Lunéville
	\since 23 feb 2012
	\date 11 may 2012

	Low level tests of Operator class and related classes.
	Almost functionalities are checked.
	This function may either creates a reference file storing the results (check=false)
	or compares results to those stored in the reference file (check=true)
	It returns reporting informations in a string
*/

#include "xlife++-libs.h"
#include "testUtils.hpp"

#include <iostream>
#include <fstream>

using namespace xlifepp;

namespace unit_operator {

// scalar spectral function
Real sinBasis(const Point& P, Parameters& pa = defaultParameters)
{
  Real x = P.x();
  Real h = pa("h");                 // get the parameter h (user definition)
  Int n = pa("basis index");        // get the index of function to compute
  return std::sqrt(2. / h) * std::sin(n * pi_ * x / h); //computation
}

//vector spectral function
Vector<Real> VecBasis(const Point& P, Parameters& pa = defaultParameters)
{
  Real x = P.x();
  Real h = pa("h");                 // get the parameter h (user definition)
  Int n = pa("basis index");        // get the index of function to compute
  Vector<Real> res(2);
  res(1) = std::sqrt(2. / h) * std::sin(n * pi_ * x / h); // computation
  res(2) = std::sqrt(2. / h) * std::cos(n * pi_ * x / h);
  return res;
}

//vector spectral function
Vector<Complex> VecBasis3(const Point& P, Parameters& pa = defaultParameters)
{
  Real x = P.x();
  Real h = pa("h");                 // get the parameter h (user definition)
  Int n = pa("basis index");        // get the index of function to compute
  Vector<Complex> res(3);
  res(1) = std::sqrt(2. / h) * std::sin(n * pi_ * x / h);
  res(2) = std::sqrt(2. / h) * std::cos(n * pi_ * x / h);
  res(3) = Complex(0,1.);
  return res;
}

Real scalF(const Point& P, Parameters& pa = defaultParameters)
{
  Real x = P.x();
  return x; //computation
}

Vector<Real> VecF(const Point& P, Parameters& pa = defaultParameters)
{
    return Vector<Real>(1,0.);
}


// matrix function
Matrix<Real> MatF(const Point& P, Parameters& pa = defaultParameters)
{
  Matrix<Real> res(2, 2);
  Real x = P.x();
  res(1, 1) = x;   res(1, 2) = 0.;
  res(2, 1) = 0.5; res(2, 2) = 2 * x;
  return res;
}

Matrix<Complex> MatFc(const Point& P, Parameters& pa = defaultParameters)
{
  Matrix<Complex> res(3, 3);
  Real x = P.x();
  res(1, 1) = Complex(0,1); res(2, 2) = x; res(3, 3) = 1.;
  return res;
}

// scalar kernel function
Real G(const Point& M, const Point& P, Parameters& pa = defaultParameters)
{
//  Real r = M.distance(P);
  return 1.;    // computation
}

String unit_operator(bool check)
{
  String rootname = "unit_operator";
  trace_p->push(rootname);
  std::stringstream out;                  // string stream receiving results
  out.precision(testPrec);

  // test of DifferentialOperator
  DifferentialOperator* do1_p = findDifferentialOperator(_id);
  out << "difop id : name=" << do1_p->name() << " order=" << do1_p->order()
      << " normal=" << words(do1_p->normalRequired())
      << " extension=" << words(do1_p->extensionRequired()) << "\n";
  findDifferentialOperator(_dt);
  findDifferentialOperator(_dx);
  findDifferentialOperator(_dy);
  findDifferentialOperator(_dz);
  findDifferentialOperator(_grad);
  findDifferentialOperator(_div);
  findDifferentialOperator(_curl);
  findDifferentialOperator(_gradS);
  findDifferentialOperator(_gradG);
  findDifferentialOperator(_divS);
  findDifferentialOperator(_divG);
  findDifferentialOperator(_curlS);
  findDifferentialOperator(_curlG);
  findDifferentialOperator(_epsilon);
  findDifferentialOperator(_epsilonG);
  findDifferentialOperator(_ntimes);
  findDifferentialOperator(_timesn);
  findDifferentialOperator(_ncrossntimes);
  findDifferentialOperator(_timesncrossn);
  findDifferentialOperator(_timesn);
  findDifferentialOperator(_ndot);
  findDifferentialOperator(_ndotgrad);
  findDifferentialOperator(_ndiv);
  findDifferentialOperator(_ncross);
  findDifferentialOperator(_ncrosscurl);
  findDifferentialOperator(_ncrossncross);
  verboseLevel(2);
  printListDiffOp(out);

  // test of Value
  Value::printValueTypeRTINames(out);
  Real s = 1.;
  Vector<Real> vr(3); Vector<Complex> vc(3, i_);
  Matrix<Real> mr(3, 3, 1.); Matrix<Complex> mc(3, 3, i_);
  Value Vs(s); out << "real Value Vs(s) :" << Vs;
  Value Vi(i_); out << "complex Value Vi(i) :" << Vi;
  Value Vvr(vr); out << "real vector Value Vvr(vr) :" << Vvr;
  Value Vvc(vc); out << "complex vector Value Vvc(vc) :" << Vvc;
  Value Vmr(mr); out << "real matrix Value Vmr(mr) :" << Vmr;
  Value Vmc(mc); out << "complex matrix Value Vmc(mc) :" << Vmc;
  out << "get real Value s= " << Vs.value<Real>() << "\n";
  out << "get complex Value i= " << Vi.value<Complex>() << "\n";
  out << "get real vector Value vr= " << Vvr.value<Vector<Real> >() << "\n";
  out << "get complex vector Value vc= " << Vvc.value<Vector<Complex> >() << "\n";
  out << "get real matrix Value mr= " << Vmr.value<Matrix<Real> >() << "\n";
  out << "get complex matrix Value mc= " << Vmc.value<Matrix<Complex> >() << "\n";

  // test of OperatorOnUnknown
  Mesh msh;  //fake mesh
  GeomDomain Omega(msh, "Omega", 3);
  Parameters ps;
  ps<<Parameter(1.,"h")<<Parameter(1,"basis index");
  Function F(sinBasis, "sinBasis", ps);
  Function VF(VecBasis, "(sin,cos)", ps);
  Function MF(MatF, "mat function", ps);
  Space V(Omega, Function(sinBasis, "sinBasis", ps), 10, "sinus basis");
  Unknown u(V,"u");
  out << "definition of unknown u : " << u;
  Unknown w(V, "w", 3);
  out << "definition of unknown w : " << w;

  // test constructor and assignment
  OperatorOnUnknown opu(u), opgu(grad(u));
  out << "test OperatorOnUnknown opu(u) : " << opu;
  out << "test OperatorOnUnknown opgu(grad(u)) : " << opgu;
  opu = d0(u);
  out << "test opu=d0(u) : " << opu;
  out << "test opgu=opu : " << opgu;
  // test all primary operator
  out << "test id(u) : " << id(u);
  out << "test d0(u) : " << d0(u);
  out << "test dt(u) : " << dt(u);
  out << "test d1(u) : " << d1(u);
  out << "test dx(u) : " << dx(u);
  out << "test d2(u) : " << d2(u);
  out << "test dy(u) : " << dy(u);
  out << "test d3(u) : " << d3(u);
  out << "test dz(u) : " << dz(u);
  out << "test grad(u) : " << grad(u);
  out << "test nabla(u) : " << nabla(u);
  out << "test div(w) : " << div(w);
  out << "test curl(w) : " << curl(w);
  out << "test rot(w) : " << rot(w);
//  out << "test gradS(u) : " << gradS(u);
//  out << "test nablaS(u) : " << nablaS(u);
//  out << "test divS(w) : " << divS(w);
//  out << "test curlS(w) : " << curlS(w);
//  out << "test rotS(w) : " << rotS(w);
  out << "test nx(u) : " << nx(u);
  out << "test ndot(w) : " << ndot(w);
  out << "test ncross(w) : " << ncross(w);
  out << "test ncrossncross(w) : " << ncrossncross(w);
  out << "test ndotgrad(u) : " << ndotgrad(u);
  out << "test ndiv(w) : " << ndiv(w);
  out << "test ncrosscurl(w) : " << ncrosscurl(w);
  out << "test epsilon(w) : " << epsilon(w);
  out << "test epsilonG(w) : " << epsilonG(w,1,0,0);
  out << "test epsilonG(w) : " << epsilonG(w,1,0,1,0);
//  out << "test gradG(u,1.) : " << gradG(u, 1.);
//  out << "test nablaG(u,1,2) : " << nablaG(u, 1, 2);
//  out << "test divG(w,1,2,3) : " << divG(w, 1, 2, 3);
//  out << "test curlG(w,1,2,3) : " << curlG(w, 1, 2, 3);
//  out << "test rotG(w,0,1) : " << rotG(w, 0, 1);

  // test aliases
  out << "test n*u : " << _n* u;
  out << "test u*n : " << u* _n;
  out << "test n|w : " << (_n | w);
  out << "test w|n : " << (w | _n);
  out << "test n^w : " << (_n ^ w);
  out << "test n^(n^w) : " << (_n ^ (_n ^ w));
  out << "test n|grad(u) : " << (_n | grad(u));
  out << "test grad(u)|n : " << (grad(u) | _n);
  out << "test n*div(w) : " << (_n * div(w));
  out << "test div(w)*n : " << (div(w)*_n);
  out << "test n^curl(w) : " << (_n ^ curl(w));
  out << "test n^rot(w) : " << (_n ^ rot(w));

  // test left and right operations with function
  out << "test F*opu : " << (F * opu);
  out << "test opu*F : " << (opu * F);
  out << "test F*opgu*F : " << (F * opgu * F);
  out << "test F*u : " << (F * u);
  out << "test u*F : " << (u * F);
  out << "test F*grad(u)*F : " << (F * grad(u)*F);
  out << "test VF|grad(u)*F : " << (VF | grad(u)*F);
  out << "test VF^grad(u)*VF : " << ((VF ^ grad(u)) * VF);  //cross product in 2D returns a scalar !
  out << "test MF*grad(u)*F : " << (MF * grad(u)*F);
  out << "test (MF*grad(u))|VF : " << ((MF * grad(u)) | VF);

  // test using explicit C++ functions
  out << "test u*scalF : " << (u * scalF);
  out << "test scalF*u : " << (scalF * u);
  out << "test grad(u)*scalF : " << (grad(u)*scalF);
  out << "test (adj(MatF)*grad(u))|VecF : " << ((adj(MatF)*grad(u)) | VecF);

  // test left and right operations with constant
  out << "test u*vr: " << (u * vr);
  out << "test s*grad(u) : " << (s * grad(u));
  out << "test i*grad(u) : " << (i_ * grad(u));
  out << "test vr|grad(u) : " << (vr | grad(u));
  out << "test vc|grad(u) : " << (vc | grad(u));
  out << "test mr*grad(u) : " << (mr * grad(u));
  out << "test mc*grad(u) : " << (mc * grad(u));
  out << "test s*grad(u)*i : " << (s * grad(u)*i_);
  out << "test i*grad(u)|vr : " << (i_ * grad(u) | vr);
  out << "test (vr|grad(u))*mr : " << ((vr | grad(u))*mr);

  // test conjugate/transpose and adjoint operator
  out << "test grad(conj(u)) : " << grad(conj(u));
  out << "test trans(mr)*grad(conj(u)) : " << (trans(mr)*grad(conj(u)));
  out << "test grad(conj(u))|conj(vc) : " << (grad(conj(u)) | conj(vc));
  out << "test (adj(mc)*grad(conj(u)))|conj(vc) : " << (adj(mc)*grad(conj(u)) | conj(vc));
  out << "test (adj(mc)*grad(conj(u)))|adj(vc) : " << (adj(mc)*grad(conj(u)) | adj(vc));

  // test expressions with component
  out << "test grad(w[1]) : " << grad(w[1]);
  out << "test grad(w[2])*s : " << (grad(w[2])*s);
  out << "test (adj(mc)*grad(conj(w[2])))|adj(vc) : " << (adj(mc)*grad(conj(w[2])) | adj(vc));
  out << "test adj(mc)*(grad(conj(w[2])))*adj(mc)) : " << (adj(mc) * (grad(conj(w[2]))*adj(mc)));

  //test kernel form
  Kernel ker=Helmholtz3dKernel(1.);
  //Kernel ker=Laplace3dKernel();
  out << "\ntest ker*u : " << ker*u;
  out << "\ntest u*ker : " << u*ker;
  out << "\ntest u*ker*u : " << u*ker*u;

  //out << "\ntest G*u : " << G*u;                      //not working on visual studio, ambiguity ?
  out << "\ntest u*G : " << u*G;
  out << "\ntest u*G*u : " << u*G*u;
  out << "\ntest u*(grad_x(ker)|_nx)*u : " << u*(grad_x(ker)|_nx)*u;
  out << "\ntest u*grad_y(grad_x(ker)|_nx)*u : " << u*grad_y(grad_x(ker)|_nx)*u;
  out << "\ntest u*grad_y(grad_x(ker)|_nx)|ny*u : " << u*(grad_y(grad_x(ker)|_nx)|_ny)*u;
  out << "\ntest u*(_nx^(_ny^ker))*u : " << u*(_nx^(_ny^ker))*u;
  out << "\ntest u*(grad_y(grad_x(ker))*u : " << u*(grad_y(grad_x(ker)))*u;
  //out << "\ntest u*(_nx^(_ny^G))*u : " << u*(_nx^(_ny^G))*u;  //ambiguity between ny^Function et ny^Kernel, G may be cast to both

  //test operator on function
  out << "\ntest id(F) : " << id(F);
  out << "\ntest ntimes(F) : " << ntimes(F);
  out << "\ntest ndot(VF) : " << ndot(VF);
  out << "\ntest ncross(VF) : " << ncross(VF);
  out << "\ntest ncrossncross(VF) : " << ncrossncross(VF);
  out << "\ntest n*F : " << _n*F;
  out << "\ntest F*n : " << F *_n;
  out << "\ntest (n^n)*F : " << ((_n^_n)*F);
  out << "\ntest n|VF : " << (_n|VF);
  out << "\ntest VF|*n : " << (VF |_n);
  out << "\ntest n^VF : " << (_n^VF);
  out << "\ntest n^(n^VF) : " << (_n^(_n^VF));
  out << "\ntest (n^n)^VF) : " << ((_n^_n)^VF);

  //evaluate operator on function
  Point M(1.);
  Real val;
  Vector<Real> v_val, n(2); n[0]=1.; n[1]=1.;
  out << "\ntest id(F)(M) : " << id(F).eval(M,val);
  out << "\ntest id(VF)(M) : " << id(VF).eval(M,v_val);
  out << "\ntest (VF|_n)(M) : " << (VF|_n).eval(M,val,&n);
  out << "\ntest (_n^VF)(M) : " << (_n^VF).eval(M,val,&n);  //2D ->scalar result
  out << "\ntest (MatF*_n)(M) : " << (MatF*_n).eval(M,v_val,&n);
  out << "\ntest (_n*MatF)(M) : " << (_n*MatF).eval(M,v_val,&n);

  Matrix<Real> A(2,2);A(1,1)=10;
  out << "\ntest (A*_n)(M) : " << (A*_n).eval(M,v_val,&n);

  Function H(VecBasis3, "(sin,cos,1)", ps);
  n.resize(3);n[2]=0.;
  Complex cval;
  Vector<Complex> cv_val;
  out << "\ntest id(H)(M) : " << id(H).eval(M,cv_val);
  out << "\ntest (H|_n)(M) : " << (H|_n).eval(M,cval,&n);
  out << "\ntest (_n^H)(M) : " << (_n^H).eval(M,cv_val,&n);
  out << "\ntest (MatFc*_n)(M) : " << (MatFc*_n).eval(M,cv_val,&n);

  //test operator involving TermVector
  verboseLevel(25);
  Mesh msq(Square(_origin=Point(0.,0.),_length=1.,_nnodes=5,_domain_name="Omega"),_triangle,1,_structured,"msq");
  Domain omega=msq.domain("Omega");
  Space Vh(omega,P1,"Vh");  Unknown uh(Vh,"uh");
  TermVector x1(uh,omega,_x1,"x1");
  out << "test x1*uh : " << x1*uh;
  out << "test x1*uh*x1 : " << x1*uh*x1;
  out << "test (x1^3)*uh : " << (x1^3)*uh;
  out << "test intg(omega,(x1^2)*uh) : " << intg(omega,(x1^2)*uh);


  //------------------------------------------------------------------------------------
  // save results in a file or compare results with some references value in a file
  //------------------------------------------------------------------------------------
  trace_p->pop();
  if (check) { return diffResFile(out, rootname); }
  else { return saveResToFile(out, rootname); }
}

}