CDMATH/tests/cdmath/MeshTests.cxx

   1 /*
   2  * meshtests.cxx
   3  *
   4  *  Created on: 24 janv. 2012
   5  *      Authors: CDMATH
   6  */
   7
   8 #include "MeshTests.hxx"
   9 #include "MEDLoader.hxx"
  10 #include "Face.hxx"
  11 #include "Cell.hxx"
  12 #include "Node.hxx"
  13
  14 #include <string>
  15 #include <cmath>
  16
  17 using namespace MEDCoupling;
  18 using namespace std;
  19
  20 void
  21 MeshTests::testNormals(Mesh mesh)
  22 {
  23         double eps=1.E-10;
  24     int dim     = mesh.getSpaceDimension();
  25     int nbCells = mesh.getNumberOfCells();
  26     int nbFaces;//local number of faces around a cell
  27     Face Fk;
  28     Cell Cj;
  29     int indexFace;
  30
  31     vector<double> test(dim,0);
  32     double norm;
  33
  34     for(int j=0; j<nbCells; j++)
  35     {
  36         Cj = mesh.getCell(j);
  37         nbFaces = Cj.getNumberOfFaces();
  38
  39         for(int i=0; i<dim; i++)
  40             test[i]=0;
  41
  42         for(int k=0; k<nbFaces; k++)
  43         {
  44             indexFace = Cj.getFacesId()[k];
  45             Fk = mesh.getFace(indexFace);
  46
  47             norm=0;
  48             for(int i=0; i<dim; i++)
  49             {
  50                 test[i] += Cj.getNormalVector(k, i)*Fk.getMeasure();
  51                 norm += Cj.getNormalVector(k, i) * Cj.getNormalVector(k, i);
  52             }
  53             CPPUNIT_ASSERT_DOUBLES_EQUAL( 1., norm, eps );
  54         }
  55
  56         for(int i=0; i<dim; i++)
  57             CPPUNIT_ASSERT_DOUBLES_EQUAL( 0., test[i], eps  );
  58     }
  59 }
  60 //----------------------------------------------------------------------
  61 void
  62 MeshTests::testClassMesh( void )
  63 //----------------------------------------------------------------------
  64 {
  65         double eps=1.E-10;
  66
  67         // Testing Mesh(xmin, xMax, nx)
  68         Mesh M1(0.0,4.0,4);
  69         CPPUNIT_ASSERT_EQUAL( 1, M1.getSpaceDimension() );
  70         CPPUNIT_ASSERT_EQUAL( 5, M1.getNumberOfNodes() );
  71         CPPUNIT_ASSERT_EQUAL( 4, M1.getNumberOfCells() );
  72         CPPUNIT_ASSERT_EQUAL( 5, M1.getNumberOfFaces() );
  73         CPPUNIT_ASSERT_EQUAL( 4, M1.getNumberOfEdges() );
  74         CPPUNIT_ASSERT_EQUAL( 0., M1.getFace(0).x() );
  75         CPPUNIT_ASSERT_EQUAL( 0., M1.getNode(0).x() );
  76         CPPUNIT_ASSERT_EQUAL( 1., M1.getFace(1).x() );
  77         CPPUNIT_ASSERT_EQUAL( 1., M1.getNode(1).x() );
  78         CPPUNIT_ASSERT_EQUAL( 2., M1.getFace(2).x() );
  79         CPPUNIT_ASSERT_EQUAL( 2., M1.getNode(2).x() );
  80         CPPUNIT_ASSERT_EQUAL( 3., M1.getFace(3).x() );
  81         CPPUNIT_ASSERT_EQUAL( 3., M1.getNode(3).x() );
  82         CPPUNIT_ASSERT_EQUAL( 4., M1.getFace(4).x() );
  83         CPPUNIT_ASSERT_EQUAL( 4., M1.getNode(4).x() );
  84         double x11=M1.getCells()[1].x();
  85         double y11=M1.getCells()[1].y();
  86         CPPUNIT_ASSERT_EQUAL( x11, 1.5 );
  87         CPPUNIT_ASSERT_EQUAL( y11, 0.0 );
  88         M1.setGroupAtFaceByCoords(0.,0.,0.,1.E-14,"LeftEdge") ;
  89         M1.setGroupAtFaceByCoords(4.,0.,0.,1.E-14,"RightEdge") ;
  90         CPPUNIT_ASSERT(M1.getFace(0).isBorder()==true);
  91         CPPUNIT_ASSERT(M1.getFace(1).isBorder()==false);
  92         CPPUNIT_ASSERT(M1.getFace(2).isBorder()==false);
  93         CPPUNIT_ASSERT(M1.getFace(3).isBorder()==false);
  94         CPPUNIT_ASSERT(M1.getFace(4).isBorder()==true);
  95     CPPUNIT_ASSERT(M1.getNameOfFaceGroups().size() == 3);//There is a default group named "Boundary" that is created by the mesh class
  96         CPPUNIT_ASSERT(M1.getNameOfFaceGroups()[0].compare("RightEdge")==0);
  97         CPPUNIT_ASSERT(M1.getNameOfFaceGroups()[1].compare("LeftEdge")==0);
  98         CPPUNIT_ASSERT(M1.getNameOfFaceGroups()[2].compare("Boundary")==0);
  99     double dx1=M1.minRatioVolSurf();
 100     CPPUNIT_ASSERT_DOUBLES_EQUAL( dx1, 1., eps  );
 101
 102     cout<<endl<<"Test mesh M1 normals"<<endl;
 103     testNormals(M1);
 104
 105         // Testing Mesh(xmin, xMax, nx, ymin, yMax, ny)
 106         double xmin=0.0;
 107         double xmax=4.0;
 108         double ymin=0.0;
 109         double ymax=4.0;
 110         Mesh M2(xmin,xmax,4,ymin,ymax,4);
 111         CPPUNIT_ASSERT_EQUAL( 4, M2.getNx() );
 112         CPPUNIT_ASSERT_EQUAL( 4, M2.getNy() );
 113         CPPUNIT_ASSERT_EQUAL( 2, M2.getSpaceDimension() );
 114         CPPUNIT_ASSERT_EQUAL( 25, M2.getNumberOfNodes() );
 115         CPPUNIT_ASSERT_EQUAL( 16, M2.getNumberOfCells() );
 116         CPPUNIT_ASSERT_EQUAL( 40, M2.getNumberOfFaces() );
 117         CPPUNIT_ASSERT_EQUAL( 40, M2.getNumberOfEdges() );
 118         CPPUNIT_ASSERT(M2.isQuadrangular());
 119         int nbCellsM2 = M2.getNumberOfCells();
 120         double areaM2=0;
 121         for(int i=0; i<nbCellsM2; i++)
 122                 areaM2+=M2.getCell(i).getMeasure();
 123         CPPUNIT_ASSERT_DOUBLES_EQUAL( 16., areaM2, eps );
 124         double x1=M2.getCells()[4].x();
 125         double y1=M2.getCells()[4].y();
 126         CPPUNIT_ASSERT_EQUAL( x1, 0.5 );
 127         CPPUNIT_ASSERT_EQUAL( y1, 1.5 );
 128
 129         double x2=M2.getNodes()[24].x();
 130         double y2=M2.getNodes()[24].y();
 131         CPPUNIT_ASSERT_EQUAL( x2, 4. );
 132         CPPUNIT_ASSERT_EQUAL( y2, 4. );
 133
 134         M2.setGroupAtPlan(xmax,0,eps,"RightEdge");
 135         M2.setGroupAtPlan(xmin,0,eps,"LeftEdge");
 136         M2.setGroupAtPlan(ymin,1,eps,"BottomEdge");
 137         M2.setGroupAtPlan(ymax,1,eps,"TopEdge");
 138         CPPUNIT_ASSERT_EQUAL( 5, int(M2.getNameOfFaceGroups().size()) );//There is a default group named "Boundary" that is created by the mesh class
 139     CPPUNIT_ASSERT(M2.getNameOfFaceGroups().size() == 5);
 140         CPPUNIT_ASSERT(M2.getNameOfFaceGroups()[1].compare("BottomEdge")==0);
 141         int nbFaces=M2.getNumberOfFaces();
 142     M2.setPeriodicFaces();
 143         std::map<int,int> indexFaces=M2.getIndexFacePeriodic();
 144         for (int i=0;i<nbFaces;i++)
 145         {
 146                 double x=M2.getFaces()[i].x();
 147                 double y=M2.getFaces()[i].y();
 148                 if (y==0. && x==0.5)
 149                 {
 150                         int indexFace=M2.getIndexFacePeriodic(i);
 151                         double xi=M2.getFace(indexFace).x();
 152                         double yi=M2.getFace(indexFace).y();
 153                         CPPUNIT_ASSERT_EQUAL( xi, x );
 154                         CPPUNIT_ASSERT_EQUAL( yi, ymax );
 155                         CPPUNIT_ASSERT_EQUAL( true, M2.getFace(indexFace).isBorder() );
 156                         CPPUNIT_ASSERT_EQUAL( indexFace, indexFaces[i] );
 157                 }
 158         }
 159     double dx2=M2.minRatioVolSurf();
 160     CPPUNIT_ASSERT_DOUBLES_EQUAL( dx2, 1./4, eps  );
 161
 162     cout<<"Test mesh M2 normals"<<endl;
 163     testNormals(M2);
 164
 165     // Testing 2D simplexization (regular triangle mesh)
 166     int splittingPolicy =0;
 167         Mesh M2Triangle(xmin,xmax,4,ymin,ymax,4,splittingPolicy);
 168         CPPUNIT_ASSERT_EQUAL( 4, M2Triangle.getNx() );
 169         CPPUNIT_ASSERT_EQUAL( 4, M2Triangle.getNy() );
 170         CPPUNIT_ASSERT_EQUAL( 2, M2Triangle.getSpaceDimension() );
 171         CPPUNIT_ASSERT_EQUAL( 25, M2Triangle.getNumberOfNodes() );
 172         CPPUNIT_ASSERT_EQUAL( 32, M2Triangle.getNumberOfCells() );
 173         CPPUNIT_ASSERT_EQUAL( 40+16, M2Triangle.getNumberOfFaces() );
 174         CPPUNIT_ASSERT_EQUAL( 40+16, M2Triangle.getNumberOfEdges() );
 175         CPPUNIT_ASSERT(M2Triangle.isTriangular());
 176         int nbCellsM2Triangle = M2Triangle.getNumberOfCells();
 177         double areaM2Triangle=0;
 178         for(int i=0; i<nbCellsM2Triangle; i++)
 179                 areaM2Triangle+=M2Triangle.getCell(i).getMeasure();
 180         CPPUNIT_ASSERT_DOUBLES_EQUAL( 16., areaM2Triangle, eps );
 181
 182         M2Triangle.setGroupAtPlan(xmax,0,eps,"RightEdge");
 183         M2Triangle.setGroupAtPlan(xmin,0,eps,"LeftEdge");
 184         M2Triangle.setGroupAtPlan(ymin,1,eps,"BottomEdge");
 185         M2Triangle.setGroupAtPlan(ymax,1,eps,"TopEdge");
 186         CPPUNIT_ASSERT_EQUAL( 5, int(M2Triangle.getNameOfFaceGroups().size()) );//There is a default group named "Boundary" that is created by the mesh class
 187         CPPUNIT_ASSERT(M2Triangle.getNameOfFaceGroups()[1].compare("BottomEdge")==0);
 188         std::map<int,int> indexFacesTriangle=M2Triangle.getIndexFacePeriodic();
 189
 190     double dx2Triangle=M2Triangle.minRatioVolSurf();
 191     CPPUNIT_ASSERT_DOUBLES_EQUAL( dx2Triangle, 0.5/3.414, 0.01  );
 192
 193     cout<<"Test mesh M2Triangle normals"<<endl;
 194     testNormals(M2Triangle);
 195
 196         // Testing Mesh(xmin, xmax, nx, ymin, ymax, ny, zmin, zmax, nz) (hexaèdres)
 197         xmin=0.0;
 198         xmax=1.0;
 199         ymin=0.0;
 200         ymax=1.0;
 201         double zmin=0.0;
 202         double zmax=1.0;
 203     Mesh M3(xmin,xmax,4,ymin,ymax,4,zmin,zmax,4);
 204     CPPUNIT_ASSERT_EQUAL( 3, M3.getSpaceDimension() );
 205         CPPUNIT_ASSERT_EQUAL( 5*5*5, M3.getNumberOfNodes() );
 206         CPPUNIT_ASSERT_EQUAL( 4*4*4, M3.getNumberOfCells() );
 207         CPPUNIT_ASSERT_EQUAL( 5*4*4*3, M3.getNumberOfFaces() );
 208         CPPUNIT_ASSERT_EQUAL( 5*5*4*3, M3.getNumberOfEdges() );
 209     CPPUNIT_ASSERT(M3.isHexahedral());
 210     int nbCellsM3 = M3.getNumberOfCells();
 211     double volM3=0;
 212     for(int i=0; i<nbCellsM3; i++)
 213         volM3+=M3.getCell(i).getMeasure();
 214     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1., volM3, eps );
 215
 216         M3.setGroupAtPlan(xmax,0,eps,"RightEdge");
 217         M3.setGroupAtPlan(xmin,0,eps,"LeftEdge");
 218         M3.setGroupAtPlan(ymin,1,eps,"BottomEdge");
 219         M3.setGroupAtPlan(ymax,1,eps,"TopEdge");
 220         M3.setGroupAtPlan(zmin,2,eps,"DownEdge");
 221         M3.setGroupAtPlan(zmax,2,eps,"UpEdge");
 222         CPPUNIT_ASSERT_EQUAL( 7, int(M3.getNameOfFaceGroups().size()) );//There is a default group named "Boundary" that is created by the mesh class
 223         CPPUNIT_ASSERT(M3.getNameOfFaceGroups()[1].compare("DownEdge")==0);
 224         nbFaces=M3.getNumberOfFaces();
 225     M3.setPeriodicFaces();
 226         indexFaces=M3.getIndexFacePeriodic();
 227         for (int i=0;i<nbFaces;i++)
 228         {
 229                 double x=M3.getFaces()[i].x();
 230                 double y=M3.getFaces()[i].y();
 231                 double z=M3.getFaces()[i].z();
 232                 if (z==0. && x==1./8 && y==1./8)
 233                 {
 234                         int indexFace=M3.getIndexFacePeriodic(i);
 235                         double xi=M3.getFace(indexFace).x();
 236                         double yi=M3.getFace(indexFace).y();
 237                         double zi=M3.getFace(indexFace).z();
 238                         CPPUNIT_ASSERT_EQUAL( xi, x );
 239                         CPPUNIT_ASSERT_EQUAL( yi, y );
 240                         CPPUNIT_ASSERT_EQUAL( zi, zmax );
 241                         CPPUNIT_ASSERT_EQUAL( true, M3.getFace(indexFace).isBorder() );
 242                         CPPUNIT_ASSERT_EQUAL( indexFace, indexFaces[i] );
 243                 }
 244
 245                 if (z==0.5 && y==0. && x==1.)
 246                         CPPUNIT_ASSERT_EQUAL( -1, M3.getIndexFacePeriodic(i) );
 247         }
 248
 249     double dx3=M3.minRatioVolSurf();
 250     CPPUNIT_ASSERT_DOUBLES_EQUAL( dx3, 0.25/6, eps  );
 251
 252     cout<<"Test mesh M3 normals"<<endl;
 253     testNormals(M3);
 254
 255     // Testing copies
 256     Mesh Mcopy1(M1);
 257     CPPUNIT_ASSERT_EQUAL( 1, Mcopy1.getSpaceDimension() );
 258     CPPUNIT_ASSERT_EQUAL( 5, Mcopy1.getNumberOfNodes() );
 259     CPPUNIT_ASSERT_EQUAL( 4, Mcopy1.getNumberOfCells() );
 260     CPPUNIT_ASSERT_EQUAL( 5, Mcopy1.getNumberOfFaces() );
 261     CPPUNIT_ASSERT_EQUAL( 4, Mcopy1.getNumberOfEdges() );
 262
 263     Mcopy1=M2;
 264     CPPUNIT_ASSERT_EQUAL( 2, Mcopy1.getSpaceDimension() );
 265     CPPUNIT_ASSERT_EQUAL( 25, Mcopy1.getNumberOfNodes() );
 266     CPPUNIT_ASSERT_EQUAL( 16, Mcopy1.getNumberOfCells() );
 267     CPPUNIT_ASSERT_EQUAL( 40, Mcopy1.getNumberOfFaces() );
 268     CPPUNIT_ASSERT_EQUAL( 40, Mcopy1.getNumberOfEdges() );
 269
 270     Mesh Mcopy2;
 271     Mcopy2=Mcopy1;
 272     CPPUNIT_ASSERT_EQUAL( 2, Mcopy2.getSpaceDimension() );
 273     CPPUNIT_ASSERT_EQUAL( 25, Mcopy2.getNumberOfNodes() );
 274     CPPUNIT_ASSERT_EQUAL( 16, Mcopy2.getNumberOfCells() );
 275     CPPUNIT_ASSERT_EQUAL( 40, Mcopy2.getNumberOfFaces() );
 276     CPPUNIT_ASSERT_EQUAL( 40, Mcopy2.getNumberOfEdges() );
 277
 278
 279     // Connection with MED
 280     string fileNameVTK="TestMesh";
 281     string fileNameMED="TestMesh";
 282
 283     M2.writeMED(fileNameMED);
 284     Mesh M22(fileNameMED + ".med");
 285     CPPUNIT_ASSERT_EQUAL( 2, M22.getSpaceDimension() );
 286     CPPUNIT_ASSERT_EQUAL( 25, M22.getNumberOfNodes() );
 287     CPPUNIT_ASSERT_EQUAL( 16, M22.getNumberOfCells() );
 288     CPPUNIT_ASSERT_EQUAL( 40, M22.getNumberOfFaces() );
 289     CPPUNIT_ASSERT_EQUAL( 40, M22.getNumberOfEdges() );
 290
 291     cout<<"Test mesh M22 normals "<<endl;
 292     testNormals(M22);
 293
 294     // Testing 3D simplexization (regular tetrahedra mesh)
 295     splittingPolicy = 0;
 296     Mesh M3Tetra(xmin,xmax,4,ymin,ymax,4,zmin,zmax,4,splittingPolicy);
 297     CPPUNIT_ASSERT_EQUAL( 3, M3Tetra.getSpaceDimension() );
 298     CPPUNIT_ASSERT(M3Tetra.isTetrahedral());
 299     int nbCellsM3Tetra = M3Tetra.getNumberOfCells();
 300     double volM3Tetra=0;
 301     for(int i=0; i<nbCellsM3Tetra; i++)
 302         volM3Tetra+=M3Tetra.getCell(i).getMeasure();
 303     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1., volM3Tetra, eps );
 304
 305         M3Tetra.setGroupAtPlan(xmax,0,eps,"RightEdge");
 306         M3Tetra.setGroupAtPlan(xmin,0,eps,"LeftEdge");
 307         M3Tetra.setGroupAtPlan(ymin,1,eps,"BottomEdge");
 308         M3Tetra.setGroupAtPlan(ymax,1,eps,"TopEdge");
 309         M3Tetra.setGroupAtPlan(zmin,2,eps,"DownEdge");
 310         M3Tetra.setGroupAtPlan(zmax,2,eps,"UpEdge");
 311         CPPUNIT_ASSERT_EQUAL( 7, int(M3Tetra.getNameOfFaceGroups().size()) );//There is a default group named "Boundary" that is created by the mesh class
 312         CPPUNIT_ASSERT(M3Tetra.getNameOfFaceGroups()[1].compare("DownEdge")==0);
 313         indexFaces=M3Tetra.getIndexFacePeriodic();
 314
 315     cout<<"Test mesh M3Tetra normals"<<endl;
 316     testNormals(M3Tetra);
 317
 318     //Testing a 2D unstructured mesh (triangles)
 319     Mesh M23("meshSquare.med");
 320     CPPUNIT_ASSERT(M23.getNameOfFaceGroups().size() == 5);//There is a default group named "Boundary" that is created by the mesh class;
 321     CPPUNIT_ASSERT(M23.getNameOfFaceGroups()[3].compare("Bottom")==0);
 322     CPPUNIT_ASSERT(M23.getNameOfFaceGroups()[2].compare("Left")==0);
 323     CPPUNIT_ASSERT(M23.getNameOfFaceGroups()[1].compare("Right")==0);
 324     CPPUNIT_ASSERT(M23.getNameOfFaceGroups()[0].compare("Top")==0);
 325     CPPUNIT_ASSERT(M23.isTriangular());
 326     int nbCellsM23 = M23.getNumberOfCells();
 327     double areaM23=0;
 328     for(int i=0; i<nbCellsM23; i++)
 329         areaM23+=M23.getCell(i).getMeasure();
 330     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1., areaM23 , eps);
 331
 332     cout<<"Test mesh M23 normals"<<endl;
 333     testNormals(M23);
 334
 335     Mcopy2.writeVTK(fileNameVTK);
 336     Mcopy2.writeMED(fileNameMED);
 337     Mesh M6(fileNameMED + ".med");
 338     CPPUNIT_ASSERT_EQUAL( 2, M6.getSpaceDimension() );
 339     CPPUNIT_ASSERT_EQUAL( 25, M6.getNumberOfNodes() );
 340     CPPUNIT_ASSERT_EQUAL( 16, M6.getNumberOfCells() );
 341     CPPUNIT_ASSERT_EQUAL( 40, M6.getNumberOfFaces() );
 342     CPPUNIT_ASSERT_EQUAL( 40, M6.getNumberOfEdges() );
 343
 344     /*
 345     const MEDCouplingMesh* M1MEDMesh = M2.getMEDCouplingMesh();
 346      */
 347
 348     //Test of a mesh with spaceDim=3 different from meshDim=2 (triangles)
 349     Mesh M4("meshSphere.med");
 350     CPPUNIT_ASSERT(M4.isTriangular());
 351     int nbCellsM4 = M4.getNumberOfCells();
 352     double areaM4=0;
 353     for(int i=0; i<nbCellsM4; i++)
 354         areaM4+=M4.getCell(i).getMeasure();
 355     CPPUNIT_ASSERT_DOUBLES_EQUAL( 4*3.14, areaM4, 1 );
 356
 357     cout<<"Test mesh M4 normals"<<endl;
 358     testNormals(M4);
 359
 360     //Testing a 3D unstructured mesh (tétraèdres)
 361     Mesh M5("meshCube.med");
 362     CPPUNIT_ASSERT(M5.isTetrahedral());
 363     int nbCellsM5 = M5.getNumberOfCells();
 364     double volM5=0;
 365     for(int i=0; i<nbCellsM5; i++)
 366         volM5+=M5.getCell(i).getMeasure();
 367     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1., volM5, eps );
 368
 369     cout<<"Test mesh M5 normals"<<endl;
 370     testNormals(M5);
 371
 372     //Testing Mesh( std::vector<double> points, std::string meshName )
 373     int nbCellsM7 = 2*11;
 374     int nbNodes = nbCellsM7+1;
 375     vector<double> points (nbNodes);
 376     xmin=0;
 377     xmax=1;
 378     double dx_min = (xmax-xmin)*2/nbCellsM7/3;
 379     double dx_max = 2*dx_min;
 380     points[0]=0;
 381     for(int i=0; i<nbNodes-1; i++)
 382     {
 383         if(i%2==0)
 384             points[i+1] = (dx_min+dx_max)*(i/2) + dx_min;
 385         else
 386             points[i+1] = (dx_min+dx_max)*(i/2) + dx_min + dx_max;
 387     }
 388     CPPUNIT_ASSERT_DOUBLES_EQUAL(points[0],        xmin,eps);
 389     CPPUNIT_ASSERT_DOUBLES_EQUAL(points[nbNodes-1],xmax,eps);
 390
 391     Mesh M7(points, "Checkerboard mesh");
 392
 393     double volM7=0;
 394     for(int i=0; i<nbCellsM7; i++)
 395         volM7+=M7.getCell(i).getMeasure();
 396     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1., volM7, eps );
 397 }