tests/unit/unit_TermVector.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_TermVector.cpp
	\author E. Lun�ville
	\since 5 oct 2012
	\date 21 aug 2013

	Low level tests of TestVector 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"

using namespace xlifepp;

namespace unit_TermVector {

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

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

Real fn(const Point& P, Parameters& pa = defaultParameters)
{
	Vector<Real>& n = getN();  //new method
	return dot(P,Point(n));
}

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

 // 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=11,_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
  Space V1(omega,_Lagrange,1,"V1",false);
  Unknown u(V1,"u"); TestFunction v(u,"v");

  //create TermVector from linear forms
  Function f1(one,"one");
  Function fxy(xy,"xy");
  TermVector un(v,omega,1.,"un");
  LinearForm lf=intg(omega,f1*v);
  TermVector B1(lf,"B");
  out<<B1;
  out<<"\n B1|un="<<realPart(B1|un)<<"\n\n";
  std::cout<<"fin compute B1\n";
  LinearForm lfxy=intg(omega,xy*v);
  TermVector Bxy(lfxy,"Bxy");
  out<<Bxy;
  out<<"\n Bxy|un="<<realPart(Bxy|un)<<"\n\n";
  std::cout<<"fin compute Bxy\n";
  LinearForm lfb1=intg(gamma2,v);
  TermVector C1(lfb1,"C1");
  out<<C1;
  //saveToFile("C1", C1, vtk);

  //test product
  out<<TermVector(2*intg(omega,f1*v));
  out<<TermVector(pi_*intg(omega,f1*v));
  out<<TermVector(i_*intg(omega,f1*v));
  out<<TermVector(Complex(3,0)*intg(omega,f1*v));
  out<<TermVector(Complex(0,3)*intg(omega,f1*v));

  //create TermVector from function
  TermVector F1(u,omega,fxy,"F1");
  out<<"TermVector from Function"<<F1;

  //create TermVector from symbolic function
  TermVector F2(u,omega,x_1*x_2,"F2");
  out<<"TermVector from SymbolicFunction"<<F2;

  //create TermVector from function involving normal vector
  //Parameters pa(0,"_n");
  Function f_fn(fn);f_fn.requireNx=true;
  TermVector Fn(u,gamma2,f_fn,"Fn");
  out<<"TermVector from Function involving normal vector"<<Fn;

  //examples of inner product
  TermVector ung2(v,gamma2,1.,"ung2");
  TermVector ung1(v,gamma1,1.,"ung1");
  out<<"\n C1|ung2="<<realPart(C1|ung2)<<"\n";
  out<<"\n C1|ung1="<<realPart(C1|ung1)<<"\n";
  out<<"\n C1|un="<<realPart(C1|un)<<"\n";
  std::cout<<"fin compute C1\n";

  //create Lagrange P2 space and unknown
  Interpolation* LagrangeP2=findInterpolation(Lagrange,standard,2,H1);
  Space V2(omega,*LagrangeP2,"V2",false);
  Unknown u2(V2,"u2");
  TestFunction v2(u2,"v2");

  TermVector un2(v2,omega,1.,"un2");
  LinearForm lf2=intg(omega,f1*v2);
  TermVector B2(lf2,"B2");
  out<<B2;
  out<<"\n B2|un2="<<realPart(B2|un2)<<"\n\n";
  std::cout<<"fin compute B2\n";

  LinearForm lfxy2=intg(omega,xy*v2);
  TermVector Bxy2(lfxy2,"Bxy2");
  out<<Bxy2;
  out<<"\n Bxy2|un2="<<realPart(Bxy2|un2)<<"\n\n";
  std::cout<<"fin compute Bxy2\n";

  LinearForm lfb2=intg(gamma2,v2);
  TermVector C2(lfb2,"C2");
  out<<C2;
  std::cout<<"fin compute C2\n";

  //examples of inner product
  TermVector ung22(v2,gamma2,1.,"ung22");
  TermVector ung12(v2,gamma1,1.,"ung12");
  out<<"\n C2|ung22="<<realPart(C2|ung22)<<"\n";
  out<<"\n C2|ung12="<<realPart(C2|ung12)<<"\n";
  out<<"\n C2|un2="<<realPart(C2|un2)<<"\n";

  //linear algebra on TermVector
  TermVector t1(u, omega, xy), t2(u,gamma1,one);
  TermVector t=t1+t2;t.name()="";
  out<<"\n t1 = "<<t1<<"\n";
  out<<"\n t2 = "<<t2<<"\n";
  out<<"\n t1 + t2 = "<<t<<"\n";

  t=2*t1+Complex(0.,1.)*t2;
  out<<"\n 2*t1+Complex(0.,1.)*t2 = "<<t<<"\n";

  t=2*t1-Complex(0.,1.)*t2;
  out<<"\n 2*t1-Complex(0.,1.)*t2 = "<<t<<"\n";

  out<<"\n abs(t) = "<<abs(t)<<"\n";
  out<<"\n imag(t) = "<<imag(t)<<"\n";
  out<<"\n real(t) = "<<real(t)<<"\n";
  out<<"\n conj(t) = "<<conj(t)<<"\n";

  Real dt=2.;
  t=(t1-t2)/dt;
  out<<"\n (t1-t2)/dt = "<<t<<"\n";
  out<<"\n on the fly abs((t1-t2)/dt) = "<<abs((t1-t2)/dt)<<"\n";
  out<<"\n on the fly out<<(t1-t2)/dt = "<<(t1-t2)/dt<<"\n";

  //interpolation of TermVector
  Point p(0.5,0.5), q(0.11,0.33);
  Real r=0.;
  out<<"t1.subVector(u)(p,r) = "<<t1.subVector(u)(p,r)<<"\n";
  out<<"t1.subVector(u)(q,r) = "<<t1.subVector(u)(q,r)<<"\n";
  out<<"t1(p,r) = "<<t1(p,r)<<"\n";
  out<<"t1(q,r) = "<<t1(q,r)<<"\n";
  out<<"t1(u,p,r) = "<<t1(u,p,r)<<"\n";
  out<<"t1(u,q,r) = "<<t1(u,q,r)<<"\n";

  //assign constant value to a part of a TermVector
  TermVector ts(u2,omega,fxy);
  ts.setValue(gamma1,1.);
  ts.setValue(gamma2,2.);
  Domain gamma=merge(gamma1,gamma2,"Gamma");
  ts.setValue(gamma,0.);
  ts.setValue(gamma1,x); //assign function values to a part of a TermVector
  ts.setValue(gamma,0.);
  TermVector tg(u2,gamma1,x);
  ts.setValue(tg);      //assign TermVector values to a part of a TermVector
  //saveToFile("ts",ts,vtu);

  //interpolation on a domain (TermVector::mapTo)
  Mesh mdisk(Disk(_center = Point(0.5,0.5),_radius=0.5, _nnodes=11,_domain_name="Disk"),_triangle,1,_gmsh);
  Domain disk=mdisk.domain("Disk");
  Space Vd(disk,P1,"Vd"); Unknown vd(Vd,"vd");
  TermVector td=t1.mapTo(disk,vd, false);
  //saveToFile("t1",t1,_vtk);
  //saveToFile("td",td,_vtk);

  //non linear operation on TermVector
  out<<"by fun t1*t1 = "<<TermVector(t1,t1,fun_productR)<<eol;      //using explicit function
  out<<"by symbolic function t1^2 = "<<TermVector(t1,x_1^2)<<eol;   //using symbolic function
  out<<"t1*t1 = "<<t1*t1<<eol;                                      //using direct expression
  TermVector st1 = TermVector(t1,sin(x_1));
  out<<"st1 = "<<st1<<eol;
  out<<"abs(st1+t1) = "<<TermVector(st1,t1,abs(x_1+x_2))<<eol;
  out<<"abs(sin(t1)+t1)= "<<abs(sin(t1)+t1)<<eol;

  //linear form from TermVector
  out<<"TermVector(intg(omega,F1*v)) : "<<TermVector(intg(omega,F1*v))<<eol;

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

}