tests/unit/unit_Mesh.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_Mesh.cpp
\authors E. Lunéville, Y. Lafranche
\since 10 may 2012
\date  23 oct 2012

Low level tests of Mesh 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
*/

//===============================================================================
//library dependences
#include "xlife++-libs.h"
#include "testUtils.hpp"

//===============================================================================
//stl dependences
#include <iostream>
#include <fstream>
#include <vector>

//===============================================================================

using namespace xlifepp;

namespace unit_Mesh {

void test_SubdvMesh(std::stringstream& out);
void test_MeshFromFile(std::stringstream& out);

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

  std::vector<String> sidenames(2, "");
  sidenames[0] = "Sigma_1";

  //test P1 2D mesh
  Mesh mesh1dP1(Segment(_xmin=0, _xmax=1, _nnodes=11, _side_names=sidenames), 1, _structured, "P1 mesh of [0,1], step=0.01");
  out << "------------------- P1 mesh of a segment -------------------------\n";
  out << mesh1dP1 << std::endl;
  mesh1dP1.saveToFile("mesh1dP1.vtk", _vtk, true);
  out << "\n------------------ P1 mesh of a rectangle ------------------------\n";
  sidenames.resize(4);
  sidenames[0] = "Gamma_1";
  sidenames[3] = "Gamma_2";
  Mesh mesh2dP1(Rectangle(_xmin=0, _xmax=1, _ymin=1, _ymax=3, _nnodes=Numbers(3, 5), _side_names=sidenames), _triangle, 1, _structured, "P1 mesh of [0,1]x[1,3]");
  out << mesh2dP1 << std::endl;
  out << "\n------------ P1 mesh of a rectangle with sides----------------------\n";
  mesh2dP1.buildSides();
  mesh2dP1.buildVertexElements();
  verboseLevel(9);
  out << mesh2dP1 << std::endl;
  mesh2dP1.saveToFile("mesh2dP1.vtk", _vtk, true);

  //test Q1 2D mesh
  sidenames[0] = "Sigma_1";
  sidenames[2] = "Sigma_2";
  sidenames[1] = "";
  Mesh mesh2dQ1(Rectangle(_xmin=1, _xmax=2, _ymin=1, _ymax=3, _nnodes=Numbers(3, 5), _side_names=sidenames), _quadrangle, 1, _structured, "Q1 mesh of [1,2]x[1,3]");
  mesh2dQ1.domain("Omega").rename("OmegaQ");
  out << mesh2dQ1 << std::endl;
  out << "\n------------ Q1 mesh of a rectangle with sides----------------------\n";
  mesh2dQ1.buildSides();
  mesh2dQ1.buildVertexElements();
  out << mesh2dQ1 << std::endl;
  mesh2dQ1.saveToFile("mesh2dQ1.vtk", _vtk, true);

  //merging meshes P1/Q1
  Mesh mesh2dP1Q1(mesh2dP1);
  out << "\n-------------------------- copy of P1 mesh --------------------------\n";
  out << mesh2dP1Q1 << std::endl;
  mesh2dP1Q1.merge(mesh2dQ1);
  verboseLevel(24);
  out << "\n----------------------- merge P1 and Q1 mesh ------------------------\n";
  out << mesh2dP1Q1 << std::endl;
  mesh2dP1Q1.saveToFile("mesh2dP1Q1.vtk", _vtk, true);
  verboseLevel(5);

  out << "\n------------------ Q1 mesh of a hexahedron ------------------------\n";
  sidenames.clear(); sidenames.resize(6);
  sidenames[0] = "z=1";
  sidenames[1] = "z=5";
  sidenames[2] = "y=1";
  sidenames[3] = "y=3";
  sidenames[4] = "x=0";
  sidenames[5] = "x=1";
  Mesh mesh3dQ1(Cuboid(_xmin=0, _xmax=1, _ymin=1, _ymax=3, _zmin=1, _zmax=5, _nnodes=Numbers(3, 5, 9), _side_names=sidenames), _hexahedron, 1, _structured, "Q1 mesh of [0,1]x[1,3]x[1,5]");
  out << mesh3dQ1 << std::endl;
  mesh3dQ1.saveToFile("mesh3dQ1.vtk", _vtk, true);

  out << "\n--------- convert hexahedron mesh to pyramid mesh -------------\n";
  Mesh mesh3DPyramid(mesh3dQ1,_pyramid);
  out << mesh3DPyramid << std::endl;
  mesh3DPyramid.saveToFile("mesh3DPyramid.msh", _msh, false);

  out << "\n--------- mesh of cube using pyramids -------------\n";
  Mesh mesh3DPyramid2(Parallelepiped(_v1=Point(0,1,1), _v2=Point(1,1,1), _v4=Point(0,3,1),_v5=Point(0,1,5), _nnodes=Numbers(3, 5, 9), _side_names=sidenames),
                      _pyramid, 1, _structured, "Q1 mesh of [0,1]x[1,3]x[1,5]");
  out << mesh3DPyramid2 << std::endl;
  mesh3DPyramid2.saveToFile("mesh3DPyramid2.msh", _msh, false);

  out << "\n------------------ extrusion of 1D mesh ------------------------\n";
  Mesh mesh1dExtrude(mesh1dP1,Point(0.,0.),Point(0.,1.),10,1,1,0);
  out << mesh1dExtrude << std::endl;
  mesh1dExtrude.saveToFile("mesh1dExtrude.vtk", _vtk, true);

  out << "\n------------------ extrusion of 2D mesh ------------------------\n";
  Mesh mesh2dExtrudeP1(mesh2dP1,Point(0.,0.,0.),Point(0.,0.,1.),10,1,1,1);
  out << mesh2dExtrudeP1 << std::endl;
  mesh2dExtrudeP1.saveToFile("mesh2dExtrudeP1.vtk", _vtk, true);

  out << "\n------------------ extrusion of 2D mesh ------------------------\n";
  Mesh mesh2dExtrudeP1b(mesh2dP1,Point(0.,0.,0.),Point(0.,0.,1.),10,2,2,2);
  out << mesh2dExtrudeP1b << std::endl;
  mesh2dExtrudeP1b.saveToFile("mesh2dExtrudeP1b.vtk", _vtk, true);

  out << "\n------------------ extrusion of 2D mesh ------------------------\n";
  Mesh mesh2dExtrudeQ1(mesh2dQ1,Point(0.,0.,0.),Point(0.,0.,1.),10,1,1,1);
  out << mesh2dExtrudeQ1 << std::endl;
  mesh2dExtrudeQ1.saveToFile("mesh2dExtrudeQ1.vtk", _vtk, true);

  out << "\n-------------- parallelepiped mesh with prisms -------------------\n";
  Mesh mesh3dPrism(Parallelepiped(_v1=Point(0.,0.,0.), _v2=Point(0.,0.,1.), _v4=Point(1.,0.,0.), _v5=Point(0.,1.,0.),
                _nnodes=10, _side_names="Gamma"), _prism, 1, _structured, "Prism mesh of cube");
  out << mesh3dPrism << std::endl;
  mesh3dPrism.saveToFile("mesh3dPrism.vtk", _vtk, false);

  //    out<<"\n------------ P1 mesh of a rectangle with sides----------------------\n";
  //    mesh3dP1.buildSides();
  //    mesh3dP1.buildVertexElements();
  //    verboseLevel(9);
  //    out<<mesh3dP1;

  out << "\n------------------ creation of underlying P1 mesh ------------------------\n";
  //test creation of underlying P1 mesh
  Mesh* mesh2dQ1_P1= mesh2dQ1.createFirstOrderMesh();
  out << *mesh2dQ1_P1 << std::endl;

  //test creation of underlying first order mesh
  Mesh& mesh2dQ1ToP1= *mesh2dQ1.createFirstOrderMesh();
  out << mesh2dQ1ToP1 << std::endl;

  Mesh mesh1dFromBasic(Segment(_v1=Point(0.,0.,0.),_v2=Point(1.,2.,3.),_nnodes=12),1,_structured);
  Mesh mesh2dP1FromBasic(Rectangle(_xmin=0,_xmax=2,_ymin=0,_ymax=4,_nnodes=12),_triangle,1,_structured);
  Mesh mesh2dQ1FromBasic(Rectangle(_xmin=0.,_xmax=2.,_ymin=0.,_ymax=4.,_nnodes=12),_quadrangle,1,_structured);
  Mesh mesh3dQ1FromBasic(Cuboid(_xmin=0.,_xmax=2.,_ymin=0.,_ymax=4.,_zmin=0.,_zmax=1.,_nnodes=12),_hexahedron,1,_structured);

  //test of meshes built using a subdivision algorithm
  test_SubdvMesh(out);

  //test of meshes read from a file
  test_MeshFromFile(out);

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

//===============================================================================

/*!
   Test functions for subdivision algorithm.
*/

//! Tests for the sphere
void test_Subdv_Sph(std::stringstream& out, dimen_t nboctants = 8) {
  number_t order;
  Ball sphu; // unit sphere by default

  order=1;
  Mesh mesh3dP1S0VolSph(sphu, _tetrahedron, order, _subdiv, "mesh3dP1S0VolSph");
  out << "\n---------------- Volume P1 mesh of the unit sphere (subdiv 0) -------------------\n";
  out << mesh3dP1S0VolSph << std::endl;

  order=1;
  sphu = Ball(_center=Point(0.,0.,0.),_radius=1,_nboctants=nboctants,_nnodes=3);
  Mesh mesh3dP1S1VolSph(sphu, _tetrahedron, order, _subdiv, "mesh3dP1S1VolSph");
  out << "\n---------------- Volume P1 mesh of the unit sphere (subdiv 1) -------------------\n";
  out << mesh3dP1S1VolSph << std::endl;
  //mesh3dP1S1VolSph.saveToFile("mesh3dP1S1VolSph.vtk", _vtk, true);


  order=2;
  Mesh mesh3dP2S1VolSph(sphu, _tetrahedron, order, _subdiv, "mesh3dP2S1VolSph");
  out << "\n---------------- Volume P2 mesh of the unit sphere (subdiv 1) -------------------\n";
  out << mesh3dP2S1VolSph << std::endl;


  order=1;
  sphu = Ball(_center=Point(0.,0.,0.),_radius=1.,_nboctants=nboctants,_nnodes=2);
  Mesh mesh3dP1S0SurfSph(sphu, _triangle, order, _subdiv, "mesh3dP1S0SurfSph");
  out << "\n---------------- Surface P1 mesh of the unit sphere (subdiv 0) -------------------\n";
  out << mesh3dP1S0SurfSph << std::endl;

  order=1;
  sphu = Ball(_center=Point(0.,0.,0.),_radius=1.,_nboctants=nboctants,_nnodes=3);
  Mesh mesh3dP1S1SurfSph(sphu, _triangle, order, _subdiv, "mesh3dP1S1SurfSph");
  out << "\n---------------- Surface P1 mesh of the unit sphere (subdiv 1) -------------------\n";
  out << mesh3dP1S1SurfSph << std::endl;

  order=2;
  Mesh mesh3dP2S1SurfSph(sphu, _triangle, order, _subdiv, "mesh3dP2S1SurfSph");
  out << "\n---------------- Surface P2 mesh of the unit sphere (subdiv 1) -------------------\n";
  out << mesh3dP2S1SurfSph << std::endl;
}

//! Tests for the cone
void test_Subdv_Cone(std::stringstream& out) {
  std::string fname;
  number_t order=1;

  fname = "P1VolMeshTetCone";
  number_t n=5;
  Real radius=1.;
  RevCone cone(_center=Point(0.,0.,2.), _radius=radius, _apex=Point(-1.,-1.,0.), _nnodes=n);
  Mesh mesh3dP1S2VolConeT(cone, _tetrahedron, order, _subdiv, fname);
  out << "\n---------------- Volume P1 mesh of the cone (subdiv 2) with tetrahedrons -------------------\n";
  out << mesh3dP1S2VolConeT << std::endl;
  fname = "P1VolMeshTetTrunk";
  Real radius1=0.6, radius2=1.;
  RevTrunk trunk1(_center1=Point(-1.,-1.,0.), _radius1=radius1, _center2=Point(0.,0.,2.), _radius2=radius2,
                 _end1_shape=_gesCone, _end1_distance=1.5, _end2_shape=_gesEllipsoid, _end2_distance=0.7, _nnodes=n);
  Mesh mesh3dP1S2VolTrunkTE(trunk1, _tetrahedron, order, _subdiv, fname);
  out << "\n---------------- Volume P1 mesh of the trunk (subdiv 2) with tetrahedrons -------------------\n";
  out << mesh3dP1S2VolTrunkTE << std::endl;

  fname = "P1VolMeshHexTrunk";
  radius1=0.5;//radius1=1.e-15;
  RevTrunk trunk2(_center1=Point(-1.,-1.,0.), _radius1=radius1, _center2=Point(0.,0.,2.), _radius2=radius2, _nnodes=n);
  Mesh mesh3dP1S2VolTrunkH(trunk2, _hexahedron, order, _subdiv, fname);
  out << "\n---------------- Volume P1 mesh of the trunk (subdiv 2) with hexahedrons -------------------\n";
  out << mesh3dP1S2VolTrunkH << std::endl;

  fname = "P1SurfMeshQuaTrunk";
  RevTrunk trunk3(_center1=Point(-1.,-1.,0.), _radius1=radius1, _center2=Point(0.,0.,2.), _radius2=radius2, _nnodes=n);
  {std::stringstream ss; ss << fname << trunk3.nbSubdomains() << ".tex"; trunk3.teXFilename(ss.str());}
  Mesh mesh3dP1S2SurfConeQ(trunk3, _quadrangle, order, _subdiv, fname);
  out << "\n---------------- Surface P1 mesh of the cone (subdiv 2) with quadrangles -------------------\n";
  out << mesh3dP1S2SurfConeQ << std::endl;

  fname = "P1SurfMeshTriTrunk";
  {std::stringstream ss; ss << fname << trunk3.nbSubdomains() << ".tex"; trunk3.teXFilename(ss.str());}
  Mesh mesh3dP1S2SurfConeT(trunk3, _triangle, order, _subdiv, fname);
  out << "\n---------------- Surface P1 mesh of the cone (subdiv 2) with triangles -------------------\n";
  out << mesh3dP1S2SurfConeT << std::endl;
}

//! Tests for the cube
void test_Subdv_Cub(std::stringstream& out, dimen_t nboctants = 8) {
  number_t order=1;
  Cube cube(_center=Point(0.,0.,0.),_length=2.,_nboctants=nboctants, _nnodes=2); cube.rotate3d(Point(0.,0.,0.), 1.,0.,0., pi_);

  Mesh mesh3dP1S0VolCubT(cube, _tetrahedron,  order, _subdiv, "mesh3dP1S0VolCubT");
  out << "\n---------------- Volume P1 mesh of the cube (subdiv 0) with tetrahedrons -------------------\n";
  out << mesh3dP1S0VolCubT << std::endl;

  Cube cube2(_center=Point(0.,0.,0.),_length=1., _nboctants=2, _nnodes=2);
  std::string fname("mesh3dP1S0VolCubH");
  Mesh mesh3dP1S0VolCubH(cube2, _hexahedron,  3, _subdiv, fname);
  out << "\n---------------- Volume P1 mesh of the cube (subdiv 0) with hexahedrons -------------------\n";
  out << mesh3dP1S0VolCubH << std::endl;

  fname="mesh3dQ3S0SurfCubQ";
  Mesh mesh3dQ3S0SurfCubQ(cube2, _quadrangle,  3, _subdiv, fname);
  out << "\n---------------- Volume Q3 mesh of the cube (subdiv 0) with quadrangles -------------------\n";
  out << mesh3dQ3S0SurfCubQ << std::endl;
}

//! Tests for the cylinder
void test_Subdv_Cyl(std::stringstream& out) {
  number_t order=1;
  string_t fname;
  RevCylinder cyl;

  fname="mesh3dP1S0VolCylT";
  Mesh mesh3dP1S0VolCylT(cyl, _tetrahedron, order, _subdiv, fname);
  out << "\n---------------- Volume P1 mesh of the cylinder (subdiv 0) with tetrahedrons -------------------\n";
  out << mesh3dP1S0VolCylT << std::endl;

  fname="mesh3dP1S0VolCylH";
  Real radius=1.;
  number_t n=3;
  RevCylinder cyl2(_center1=Point(0.,0.,0.), _center2=Point(0.,0.,1.), _radius=radius, _nnodes=n);
  Mesh mesh3dP1S0VolCylH(cyl2, _hexahedron, order, _subdiv, fname);
  out << "\n---------------- Volume P1 mesh of the cylinder (subdiv 0) with hexahedrons -------------------\n";
  out << mesh3dP1S0VolCylH << std::endl;

  fname="mesh3dP1S0SurfCylH";
  Mesh mesh3dP1S0SurfCylH(cyl2, _quadrangle, order, _subdiv, fname);
  out << "\n---------------- Volume P1 mesh of the cylinder (subdiv 0) with quadrangles -------------------\n";
  out << mesh3dP1S0SurfCylH << std::endl;

  fname="mesh3dP1S0SurfCylT";
  Mesh mesh3dP1S0SurfCylT(cyl2, _triangle, order, _subdiv, fname);
  out << "\n---------------- Volume P1 mesh of the cylinder (subdiv 1) with triangles -------------------\n";
  out << mesh3dP1S0SurfCylT << std::endl;
}

//! Tests for the set of elements
void test_Subdv_TrSet(std::stringstream& out) {
  std::vector<Point> pts;
  pts.push_back(Point(0.,0.,0.));
  pts.push_back(Point(1.,0.,0.));
  pts.push_back(Point(0.,1.,0.3));
  pts.push_back(Point(0.,-1.,0.3));

  std::vector<std::vector<number_t> > elems;
  { number_t T[]={1,2,3}; elems.push_back(std::vector<number_t>(T,T+sizeof(T)/sizeof(number_t))); }
  { number_t T[]={1,4,2}; elems.push_back(std::vector<number_t>(T,T+sizeof(T)/sizeof(number_t))); }
  std::vector<std::vector<number_t> > bounds;
  std::vector<String> names;
  { number_t L[]={1,4}; bounds.push_back(std::vector<number_t>(L,L+sizeof(L)/sizeof(number_t))); names.push_back("S1"); }
  { number_t L[]={4,2}; bounds.push_back(std::vector<number_t>(L,L+sizeof(L)/sizeof(number_t))); names.push_back("S2"); }
  { number_t L[]={2,3}; bounds.push_back(std::vector<number_t>(L,L+sizeof(L)/sizeof(number_t))); names.push_back("S3"); }
  { number_t L[]={1,3}; bounds.push_back(std::vector<number_t>(L,L+sizeof(L)/sizeof(number_t))); names.push_back("S4"); }

  number_t nbsubdiv(1), order(3);

  SetOfElems sot(pts, elems, bounds, names, _triangle, nbsubdiv);
  sot.teXFilename("P1SurfMeshTriSet.tex");
  Mesh meshTriSet(sot, _triangle, order, _subdiv, "meshTriSet");
  out << "\n---------------- Surface P3 mesh starting from an initial mesh of triangles (subdiv 1) -------------------\n";
  out << meshTriSet << std::endl;

  pts.clear();
  pts.push_back(Point(0.,0.));
  pts.push_back(Point(0.,1.));
  pts.push_back(Point(1.,0.));
  pts.push_back(Point(1.,1.));
  pts.push_back(Point(1.5,0.5));
  pts.push_back(Point(1.5,1.5));
  pts.push_back(Point(2.,0.));
  pts.push_back(Point(2.,1.));
  pts.push_back(Point(3.,0.));
  pts.push_back(Point(3.,1.));
  pts.push_back(Point(0.,-1.));
  pts.push_back(Point(1.,-1.));
  pts.push_back(Point(1.5,-0.5));
  pts.push_back(Point(1.5,-1.5));
  pts.push_back(Point(2.,-1.));
  pts.push_back(Point(3.,-1.));

  elems.clear();
  { number_t T[]={1,3,4,2}; elems.push_back(std::vector<number_t>(T,T+sizeof(T)/sizeof(number_t))); }
  { number_t T[]={3,5,6,4}; elems.push_back(std::vector<number_t>(T,T+sizeof(T)/sizeof(number_t))); }
  { number_t T[]={5,7,8,6}; elems.push_back(std::vector<number_t>(T,T+sizeof(T)/sizeof(number_t))); }
  { number_t T[]={7,9,10,8}; elems.push_back(std::vector<number_t>(T,T+sizeof(T)/sizeof(number_t))); }
  { number_t T[]={11,12,3,1}; elems.push_back(std::vector<number_t>(T,T+sizeof(T)/sizeof(number_t))); }
  { number_t T[]={12,14,13,3}; elems.push_back(std::vector<number_t>(T,T+sizeof(T)/sizeof(number_t))); }
  { number_t T[]={14,15,7,13}; elems.push_back(std::vector<number_t>(T,T+sizeof(T)/sizeof(number_t))); }
  { number_t T[]={15,16,9,7}; elems.push_back(std::vector<number_t>(T,T+sizeof(T)/sizeof(number_t))); }
  bounds.clear();
  names.clear();
  { number_t L[]={1,2,11}; bounds.push_back(std::vector<number_t>(L,L+sizeof(L)/sizeof(number_t))); }
  { number_t L[]={11,12,14,15,16}; bounds.push_back(std::vector<number_t>(L,L+sizeof(L)/sizeof(number_t))); }
  { number_t L[]={16,9,10}; bounds.push_back(std::vector<number_t>(L,L+sizeof(L)/sizeof(number_t))); }
  { number_t L[]={2,4,6,8,10}; bounds.push_back(std::vector<number_t>(L,L+sizeof(L)/sizeof(number_t))); }
  { number_t L[]={3,5,7,13}; bounds.push_back(std::vector<number_t>(L,L+sizeof(L)/sizeof(number_t))); }

  SetOfElems soq(pts, elems, bounds, names, _quadrangle, nbsubdiv);
  soq.teXFilename("P1SurfMeshQuaSet.tex");
  Mesh meshQuaSet(soq, _quadrangle, order, _subdiv, "meshQuaSet");
  out << "\n---------------- Surface Q3 mesh starting from an initial mesh of quadrangles (subdiv 1) -------------------\n";
  out << meshQuaSet << std::endl;
}

//! Tests for the disk
void test_Subdv_Disk_Sub(std::stringstream& out) {
  real_t radius=1.;
  number_t order=1;
  Disk pdisk(_center=Point(0.,1.), _radius=radius, _angle1=10., _angle2=300., _nnodes=5);
  pdisk.teXFilename("P1SurfMeshTriDisk.tex");
  Mesh meshTriDisk(pdisk, _triangle, order, _subdiv, "meshTriDisk");
  out << "\n---------------- Surface mesh of a portion of disk (subdiv 1) -------------------\n";
  out << meshTriDisk << std::endl;

  pdisk.teXFilename("P1SurfMeshQuaDisk.tex");
  Mesh meshQuaDisk(pdisk, _quadrangle, order, _subdiv, "meshQuaDisk");
  out << "\n---------------- Surface mesh of a portion of disk (subdiv 1) -------------------\n";
  out << meshQuaDisk << std::endl;

  Mesh submeshQuaDisk(meshQuaDisk,1);
  out << "\n---------------- Subdivision of the previous one -------------------\n";
  out << submeshQuaDisk << std::endl;
}

//! Global test function
void test_SubdvMesh(std::stringstream& out) {
  test_Subdv_Sph(out,8); // nboctants = 8;
  test_Subdv_Cub(out,7); // nboctants = 7;
  test_Subdv_Cyl(out);
  test_Subdv_TrSet(out);
  test_Subdv_Disk_Sub(out);
  test_Subdv_Cone(out);
}
//===============================================================================

/*!
   Test functions for mesh read from a file.
*/
void test_MeshFromFile(std::stringstream& out) {
  using std::vector;
  // Mesh in Gmsh format
  String fn("concentric_sphere.msh");
  Mesh mesh3DGmsh=Mesh(inputsPathTo(fn),"mesh3DGmsh",msh);
  out << "\n---------------- Read file in Gmsh format (file " << fn << ") -------------------\n";
  out << mesh3DGmsh << std::endl;

  // Mesh in Melina format
  vector<String> vfn;
  vfn.push_back("tr4x4");
  vfn.push_back("qu4x4");
  for (vector<String>::const_iterator it_vfn=vfn.begin(); it_vfn < vfn.end(); it_vfn++) {
    Mesh mesh2DMel=Mesh(inputsPathTo(*it_vfn),"mesh2DMel",mel);
    out << "\n---------------- Read file in Melina format (file " << *it_vfn << ") -------------------\n";
    out << mesh2DMel << std::endl;
  }
  fn="unitCube_1e3.ply";
  out << "\n---------------- Read file in PLY format (file " << fn << ") -------------------\n";
  Mesh mesh2DPly(inputsPathTo(fn), "mesh PLY",ply);
  out << mesh2DPly << std::endl;
  fn="square.vtk";
  out << "\n---------------- Read file in VTK format (file " << fn << ") -------------------\n";
  Mesh meshSquareVtk(inputsPathTo(fn), "mesh square VTK",vtk);
  out << meshSquareVtk << std::endl;
  fn="screen_4.msh";
  out << "\n---------------- Read file in Msh format (file " << fn << ") -------------------\n";
  Mesh mesh2DMsh(inputsPathTo(fn), "mesh Msh",msh);
  out << mesh2DMsh << std::endl;
  fn="screen_4.mesh";
  out << "\n---------------- Read file in Medit format (file " << fn << ") -------------------\n";
  Mesh mesh2DMedit(inputsPathTo(fn), "mesh Medit",medit);
  out << mesh2DMedit << std::endl;
}

}