tests/unit/unit_TermMatrix.cpp

/*
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_TermMatrix.cpp
	\author E. Lunéville
	\since 5 october 2012
	\date 11 september 2013

	Low level tests of TestMatrix 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 "../../src/term/SymbolicTermMatrix.hpp"

using namespace xlifepp;

namespace unit_TermMatrix {

// scalar 2D function
Real un(const Point& P, Parameters& pa = defaultParameters)
{
  return 1.;
}

Real xfy(const Point& P, Parameters& pa = defaultParameters)
{
  return P(1)*P(2);
}

Real G(const Point& P, const Point& Q,Parameters& pa = defaultParameters)
{
    return 1.;
}

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

  verboseLevel(30);

  // create a mesh and Domains
  std::vector<String> sidenames(4,"");
  sidenames[0]="Gamma_1";
  sidenames[3]="Gamma_2";
  verboseLevel(10);
  Mesh mesh2d(Rectangle(_xmin=0,_xmax=1,_ymin=0,_ymax=1,_nnodes=3,_side_names=sidenames),_triangle,1,_structured,"P1 mesh of [0,1]x[0,1]");
  Domain omega=mesh2d.domain("Omega");
  Domain gamma1=mesh2d.domain("Gamma_1");
//  Domain gamma2=mesh2d.domain("Gamma_2");

  // create Lagrange P1 space and unknown
  Interpolation *LagrangeP1=findInterpolation(Lagrange,standard,1,H1);
  Space V1(omega,*LagrangeP1,"V1",false);
  Unknown u(V1,"u");
  TestFunction v(u,"v");
  Function f1(un,"one");
  Function fxy(xfy,"xy");
  cpuTime("time for initialisation, mesh and space definition");

  // useful TermVector
  TermVector un(u,omega,1.,"un");
  TermVector unv(v,omega,1.,"unv");

  // test some volumic mass bilinear forms
  BilinearForm uv=intg(omega,u*v);
  TermMatrix M(uv,"UV");
  cpuTime("time for computation of M");
  M.viewStorage(out);
  out<<M;
  out<<"----> M*un|un="<<realPart(M*un|un)<<"\n\n";

  TermMatrix Md(uv,_csDualStorage,"UV");
  cpuTime("time for computation of Md");
  Md.viewStorage(out);
  out<<Md;
  out<<"----> Md*un|un="<<realPart(Md*un|un)<<"\n\n";
  clear(Md);

  BilinearForm uv1=intg(omega,f1*u*v);
  TermMatrix M1(uv1,"f1*UV");
  cpuTime("time for computation of M1");
  out<<M1;
  //clear(M1);

  BilinearForm xyuv=intg(omega,xfy*u*v);
  TermMatrix Mxy(xyuv,"xy*UV");
  cpuTime("time for computation of Mxy");
  out<<Mxy;
  out<<"----> Mxy*un|un="<<realPart(Mxy*un|un)<<"\n\n";
  clear(Mxy);

  BilinearForm uxyv=intg(omega,u*xfy*v);
  TermMatrix Mxy2(uxyv,"U*xy*V");
  cpuTime("time for computation of Mxy2");
  out<<Mxy2;
  clear(Mxy2);

  BilinearForm uvxy=intg(omega,u*(v*xfy));
  TermMatrix Mxy3(uvxy,"UV*xy");
  cpuTime("time for computation of Mxy3");
  out<<Mxy3;
  clear(Mxy3);

  // test some volumic rigidity bilinear forms
  BilinearForm gugv=intg(omega,grad(u)|grad(v));
  TermMatrix K(gugv,"K");
  cpuTime("time for computation of K");
  out<<K;
  out<<"----> K*un="<<K*un<<"\n\n";

  BilinearForm gugv1=intg(omega,f1*grad(u)|grad(v));
  TermMatrix K1(gugv1,"K1");
  cpuTime("time for computation of K1");
  out<<K1;
  clear(K1);

  BilinearForm gugvxy=intg(omega,xfy*grad(u)|grad(v));
  TermMatrix Kxy(gugvxy,"Kxy");
  cpuTime("time for computation of Kxy");
  out<<Kxy;
  clear(Kxy);

  BilinearForm dxuv=intg(omega,dx(u)*v);
  TermMatrix Dxuv(dxuv,"DxUV");
  cpuTime("time for computation of Dxuv");
  out<<Dxuv;
  out<<"----> DxUV*un="<<Dxuv*un<<"\n\n";
  clear(Dxuv);

//  ============  problem of deallocation (clear) with this guy ================
  TermMatrix Dxuvd(dxuv,"_skylineDualStorage, DxUVd");
  cpuTime("time for computation of Dxuvd");
  out<<Dxuvd;
  out<<"----> DxUVd*un="<<Dxuvd*un<<"\n\n";
  clear(Dxuvd);

  BilinearForm udxv=intg(omega,u*dx(v));
  TermMatrix Dxvu(udxv,"DxVU");
  cpuTime("time for computation of Dxvu");
  out<<Dxvu;
  out<<"----> DxVU*un="<<Dxvu*un<<"\n\n";

  // test some surfacic mass bilinear forms
  BilinearForm mg1=intg(gamma1,u*v);
  TermMatrix Mg1(mg1,"Mg1");
  cpuTime("time for computation of Mg1");
  out<<Mg1;
  TermVector ung(u,gamma1,1.,"un_gamma1");
  out<<"----> Mg1*ung|ung="<<realPart(Mg1*ung|ung)<<"\n\n";

  TermMatrix Mg1d(mg1,_skylineSymStorage,"Mg1d");
  cpuTime("time for computation of Mg1d");
  out<<Mg1d;
  out<<"----> Mg1d*ung|ung="<<realPart(Mg1d*ung|ung)<<"\n\n";

  // test some surfacic bilinear forms requiring domain extension
  TermMatrix Dxuv1(intg(gamma1,dx(u)*v),"Dxuv1");
  cpuTime("time for computation of Dxuv1");
  out<<Dxuv1;
  TermMatrix Dxvu1(intg(gamma1,u*dx(v)),"Dxvu1");
  cpuTime("time for computation of Dxvu1");
  out<<Dxvu1;
  TermMatrix Dyuv1(intg(gamma1,dy(u)*v),"Dyuv1");
  cpuTime("time for computation of Dyuv1");
  out<<Dyuv1;
  TermMatrix Dnuv1(intg(gamma1,ndotgrad(u)*v),"Dnuv1");
  cpuTime("time for computation of Dnuv1");
  out<<Dnuv1;
  TermMatrix Dtudtv(intg(gamma1,gradS(u)|gradS(v)),"Dtudtv");
  cpuTime("time for computation of Dtudtv");
  out<<Dtudtv;
  TermMatrix DtuGdtv(intg(gamma1,gamma1,gradS(u)|(G*gradS(v))),"DtuGdtv");
  cpuTime("time for computation of DtuGdtv");
  out<<DtuGdtv;

  //test some combination of bilinear forms
  Real k=1,alpha=1.;
  BilinearForm h=intg(omega,grad(u)|grad(v))-k*k*intg(omega,u*v)+alpha*intg(gamma1,u*v);
  TermMatrix H(h,"H");
  cpuTime("time for computation of H");
  out<<H;

  // test storage conversion
  H.setStorage(_skyline,_sym);
  out<<H;
  H.setStorage(_skyline,_dual);
  out<<H;
  H.setStorage(_dense,_row);
  out<<H;
  H.setStorage(_cs,_dual);
  out<<H;

  //test of linear combination of TermMatrix's
  TermMatrix Hlc=K-(k*k)*M+alpha*Mg1;
  //Hlc.name()="Hlc";
  cpuTime("time for computation of Hlc");
  out<<Hlc;

  clear(Hlc);

  //test some combination of bilinear forms
  BilinearForm h2=intg(omega,grad(u)|grad(v))+intg(omega,u*dx(v))+alpha*intg(gamma1,u*v);
  TermMatrix H2(h2,"H2");
  cpuTime("time for computation of H2");
  out<<H2;

  //test of linear combination of TermMatrix's
  TermMatrix H2lc=K+Dxvu+alpha*Mg1;
  //H2lc.name()="H2lc";
  cpuTime("time for computation of H2lc");
  out<<H2lc;
  //H2lc=H2lc-H2;
  H2lc-=H2;
  out<<H2lc;
  clear(H2lc);
  clear(H2);

  // test product and inverse
  TermMatrix M12=M1*M1;
  thePrintStream<<"M12=M1*M1="<<M12;
  TermMatrix invM1M1f=directSolve(M1,M1,_keep);
  out<<"invM1M1f=directSolve(M1,M1)="<<invM1M1f;
  TermMatrix Id(M1,_idMatrix,"Id");
  out<<Id;
  out<<"directSolve(M1,Id,_keep)="<<directSolve(M1,Id,_keep);
  TermMatrix invM1=inverse(M1);
  out<<"invM1=inverse(M1)="<<invM1;
  out<<"invM1*M1="<<invM1*M1;
  out<<"invM1*M12="<<invM1*M12;
  TermMatrix MT=invM1*M12-M1;
  out<<"norm(invM1*M12-M1)="<<norm2(MT)<<eol;
  out<<"norm(invM1*M12-M1)="<<norm2(invM1*M12-M1)<<eol;
  out<<"out<<invM1*M12-M1="<<invM1*M12-M1<<eol;

  thePrintStream<<"===== avant Mpq ====="<<eol;
  printMatrixStorages(thePrintStream);

  //test SymbolicTermMatrix
  TermMatrix Mf;
  factorize(M,Mf);
  SymbolicTermMatrix S=~M+2*inv(M)*M*(~M+~M);
  theCout<<"S="<<S<<eol;
  TermVector X(u,omega,1.);
  TermVector SX=S*X;
  TermVector S5=5*M*X;
  theCout<<"|S5-SX|="<<norm2(S5-SX)<<eol;

//  //other interpolations
//  Interpolation LagrangeP2(Lagrange,standard,2,H1);
//  Space V2(omega,LagrangeP2,"V2",false);
//  Unknown p(V2,"p");
//  TestFunction q(p,"q");
//  BilinearForm pq=intg(omega,p*q);
//  TermMatrix Mpq(pq,"pq");
//  cpuTime("time for computation of Mpq");
//  out<<Mpq;
//  TermVector pun(p,omega,1.,"pun");
//  TermVector qun(q,omega,1.,"qun");
//  out<<"\nMpq*un|un="<<realPart(Mpq*pun|qun)<<"\n\n";
//  clear(Mpq);
//
//  thePrintStream<<"===== avant Kpq ====="<<eol;
//  printMatrixStorages(thePrintStream);
//
//  BilinearForm gpgq=intg(omega,grad(p)|grad(q));
//  TermMatrix Kpq(gpgq,"Kpq");
//  cpuTime("time for computation of Kpq");
//  out<<Kpq;
//  out<<"\n"<<Kpq*pun<<"\n\n";
//  clear(Kpq);
//
//  thePrintStream<<"===== après Kpq ====="<<eol;
//  printMatrixStorages(thePrintStream);thePrintStream.flush();
//
//  Interpolation LagrangeP3(Lagrange,standard,3,H1);
//  Space V3(omega,LagrangeP3,"V3",false);
//  Unknown s(V3,"s");
//  TestFunction t(s,"t");
//  BilinearForm st=intg(omega,s*t);
//  TermMatrix Mst(st,"st");
//  cpuTime("time for computation of Mst");
//  out<<Mst;
//  TermVector sun(s,omega,1.,"pun");
//  TermVector tun(t,omega,1.,"qun");
//  out<<"\nMst*un|un="<<realPart(Mst*sun|tun)<<"\n\n";
//  clear(Mst);
//
//  BilinearForm gsgt=intg(omega,grad(s)|grad(t));
//  TermMatrix Kst(gsgt,"Kst");
//  cpuTime("time for computation of Kst");
//  out<<Kst;
//  out<<"\n"<<Kst*sun<<"\n";
//  clear(Kst);
//
//  Interpolation LagrangeP4(Lagrange,standard,4,H1);
//  Space V4(omega,LagrangeP4,"V4",false);
//  Unknown e(V4,"e");
//  TestFunction f(e,"f");
//  BilinearForm ef=intg(omega,e*f);
//  TermMatrix Mef(ef,"ef");
//  cpuTime("time for computation of Mef");
//  out<<Mef;
//  TermVector eun(e,omega,1.,"eun");
//  TermVector fun(f,omega,1.,"fun");
//  out<<"\nMef*un|un="<<realPart(Mef*eun|fun)<<"\n\n";
//  clear(Mef);
//
//  BilinearForm gegf=intg(omega,grad(e)|grad(f));
//  TermMatrix Kef(gegf,"Kef");
//  cpuTime("time for computation of Kef");
//  out<<Kef;
//  out<<"\n"<<Kef*eun<<"\n";
//  clear(Kef);

  //------------------------------------------------------------------------------------
  // 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); }

}

}