Changeset 3267

Timestamp:

03/12/10 07:27:15 (15 years ago)

Author:

Mathieu Morlighem

Message:

minor

Location:

issm/trunk/src/c/Bamgx

Files:

• Bamgx.cpp (modified) (1 diff)
• objects/BamgObjects.h (modified) (2 diffs)
• objects/QuadTree.cpp (modified) (2 diffs)
• objects/QuadTree.h (modified) (1 diff)
• objects/Triangles.cpp (modified) (15 diffs)
• objects/Triangles.h (modified) (2 diffs)
• objects/VertexOnGeom.h (modified) (2 diffs)

issm/trunk/src/c/Bamgx/Bamgx.cpp

@@ -113 +113 @@
                 if (verbosity>0) printf("Anisotropic mesh adaptation\n");
                 if (verbosity>1) printf("   Processing initial mesh and geometry...\n");
-                Triangles BTh(bamgmesh_in,bamgopts); 
+                Triangles BTh(bamggeom_in,bamgmesh_in,bamgopts); 
                 hmin = Max(hmin,BTh.MinimalHmin());
                 hmax = Min(hmax,BTh.MaximalHmax());

issm/trunk/src/c/Bamgx/objects/BamgObjects.h

@@ -2 +2 @@
 #define _BAMGOBEJECTS_H_
 
-/*Bamg objects (Must be compiled in this order!!*/
+/*Bamg objects (Must be compiled in this order!!)*/
 #include "./Metric.h"
 #include "./DoubleAndInt.h"
@@ -22 +22 @@
 #include "./Triangles.h"
 #include "./Geometry.h"
+#include "./QuadTree.h"
+#include "./SetOfE4.h"
 
 #endif

issm/trunk/src/c/Bamgx/objects/QuadTree.cpp

@@ -220 +220 @@
                 register long n0;
                 register QuadTreeBox* b;
-                IntQuad  h=MaxISize,h0;
-                IntQuad  hb=MaxISize;
+                long     h=MaxISize,h0;
+                long     hb=MaxISize;
                 Icoor1   i0=0,j0=0;
                 //iplus= i projected in [0,MaxISize-1] (example: if i<0, i=0)
@@ -356 +356 @@
                 int l; // level
                 QuadTreeBox * b;
-                IntQuad  h=MaxISize,h0;
-                IntQuad hb =  MaxISize;
+                long     h =MaxISize,h0;
+                long     hb=MaxISize;
                 Icoor1  i0=0,j0=0;
                 Icoor1  iplus( i<MaxISize?(i<0?0:i):MaxISize-1);

issm/trunk/src/c/Bamgx/objects/QuadTree.h

@@ -10 +10 @@
 namespace bamg {
 
-        typedef long IntQuad;
-
-        const int MaxDeep = 30;
-        const IntQuad MaxISize = ( 1L << MaxDeep); 
+        const int  MaxDeep = 30;
+        const long MaxISize = ( 1L << MaxDeep); 
 
         class Triangles;

issm/trunk/src/c/Bamgx/objects/Triangles.cpp

@@ -3 +3 @@
 #include <cmath>
 #include <ctime>
+#include <assert.h>
 
 #include "BamgObjects.h"
 #include "../shared/shared.h"
-#include "QuadTree.h"
-#include "SetOfE4.h"
 
 #undef __FUNCT__ 
@@ -18 +17 @@
 
         /*Constructors/Destructors*/
-        /*FUNCTION Triangles::Triangles(BamgMesh* bamgmesh, BamgOpts* bamgopts){{{1*/
-        Triangles::Triangles(BamgMesh* bamgmesh, BamgOpts* bamgopts):Gh(*(new Geometry())),BTh(*this){ 
+        /*FUNCTION Triangles::Triangles(BamgGeom* bamggeom,BamgMesh* bamgmesh, BamgOpts* bamgopts){{{1*/
+        Triangles::Triangles(BamgGeom* bamggeom,BamgMesh* bamgmesh, BamgOpts* bamgopts):Gh(*(new Geometry())),BTh(*this){ 
+
+                /*Initialize fields*/
                 PreInit(0);
+
+                /*Read background mesh*/
                 ReadMesh(bamgmesh,bamgopts);
+
+                /*Read Geometry*/
+                if(bamggeom->NumEdges==0) {
+                        /*Recreate geometry if needed*/
+                        printf("WARNING: mesh present but no geometry found. Reconstructing...\n");
+                        BuildGeometryFromMesh(bamgopts);
+                        Gh.AfterRead();
+                }
+                else{
+                        Gh.ReadGeometry(bamggeom,bamgopts);
+                        Gh.AfterRead();
+                }
+
+                /*Set integer coordinates*/
                 SetIntCoor();
-                GenerateMeshProperties();
+
+                /*Fill holes and generate mesh properties*/
+                ReconstructExistingMesh();
         }
         /*}}}1*/
@@ -32 +51 @@
                 ReadMesh(index,x,y,nods,nels);
                 SetIntCoor();
-                GenerateMeshProperties();
+                ReconstructExistingMesh();
         }
         /*}}}1*/
@@ -127 +146 @@
                   Gh.AfterRead(); 
                   SetIntCoor();
-                  GenerateMeshProperties();
+                  ReconstructExistingMesh();
 
                   if (!NbSubDomains){
@@ -523 +542 @@
                 }
 
-                /*Recreate geometry if needed*/
-                if(1!=0) {
-                        printf("WARNING: mesh present but no geometry found. Reconstructing...\n");
-                        /*Recreate geometry: */
-                        BuildGeometryFromMesh(bamgopts);
-                        Gh.AfterRead();
-                }
         }
         /*}}}1*/
@@ -2604 +2616 @@
         }
         /*}}}1*/
-        /*FUNCTION Triangles::GenerateMeshProperties{{{1*/
-        void Triangles::GenerateMeshProperties() {
+        /*FUNCTION Triangles::ReconstructExistingMesh{{{1*/
+        void Triangles::ReconstructExistingMesh() {
                 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/FillHoleInMesh)*/
+
+                /*This routine reconstruct an existing mesh to make it CONVEX:
+                 * -all the holes are filled
+                 * -concave boundaries are filled
+                 * A convex mesh is required for a lot of operations. This is why every mesh
+                 * goes through this process.
+                 * This routine also generates mesh properties such as adjencies,...
+                 */
 
                 int verbosity=0;
@@ -2614 +2634 @@
                 // find extrema coordinates of vertices pmin,pmax
                 long i;
-                if(verbosity>2) printf("      Filling holes in mesh of %i vertices\n",nbv); 
+                if(verbosity>2) printf("      Completing Triangles structure for a mesh of %i vertices\n",nbv); 
 
                 //initialize ordre
-                if (!ordre){
-                        throw ErrorException(__FUNCT__,exprintf("!ordre"));
-                }
+                assert(ordre);
                 for (i=0;i<nbv;i++) ordre[i]=0;
 
@@ -2651 +2669 @@
                                 }
                                 else if(st[k]>=0) {
-                                        if (triangles[i].TriangleAdj(j) || triangles[st[k] / 3].TriangleAdj((int) (st[k]%3))){
-                                                throw ErrorException(__FUNCT__,exprintf("(triangles[i].TriangleAdj(j) || triangles[st[k] / 3].TriangleAdj((int) (st[k]%3)))"));
-                                        }
+                                        assert(!triangles[i].TriangleAdj(j));
+                                        assert(!triangles[st[k]/3].TriangleAdj((int) (st[k]%3)));
 
                                         triangles[i].SetAdj2(j,triangles + st[k] / 3,(int) (st[k]%3));
@@ -2681 +2698 @@
                 for (i=0;i<edge4->nb();i++){
                         if (st[i] >=0){ // edge alone 
-                                if (i < nbe) {
+                                if (i<nbe){
                                         long i0=edge4->i(i);
                                         ordre[i0] = vertices+i0;
@@ -2689 +2706 @@
                                 else {
                                         k=k+1;
-                                        if (k <20) {
-                                                throw ErrorException(__FUNCT__,exprintf("Lost boundary edges %i : %i %i",i,edge4->i(i),edge4->j(i)));
+                                        if (k<10) {
+                                                //print only 10 edges
+                                                printf("Lost boundary edges %i : %i %i\n",i,edge4->i(i),edge4->j(i));
                                         }
-                                }
-                        }
-                }
-                if(k != 0) {
-                        throw ErrorException(__FUNCT__,exprintf("%i boundary edges  are not defined as edges",k));
+                                        else if (k==10){
+                                                printf("Other lost boundary edges not shown...\n");
+                                        }
+                                }
+                        }
+                }
+                if(k) {
+                        throw ErrorException(__FUNCT__,exprintf("%i boundary edges (from the geometry) are not defined as mesh edges",k));
                 }
 
@@ -2705 +2726 @@
                         vertices[i].vint=0;
                         if (ordre[i]){ 
-                                ordre[nbvb++] = ordre[i];
-                        }
-                }
-
-                Triangle* savetriangles= triangles;
+                                ordre[nbvb++]=ordre[i];
+                        }
+                }
+
+                Triangle* savetriangles=triangles;
                 long savenbt=nbt;
                 long savenbtx=nbtx;
-                SubDomain * savesubdomains = subdomains;
-                subdomains = 0;
+                SubDomain* savesubdomains=subdomains;
+                subdomains=0;
 
                 long  Nbtriafillhole = 2*nbvb;
@@ -2723 +2744 @@
 
                 for (i=2 ; det( ordre[0]->i, ordre[1]->i, ordre[i]->i ) == 0;) 
-                 if  ( ++i >= nbvb) {
-                         throw ErrorException(__FUNCT__,exprintf("GenerateMeshProperties: All the vertices are aligned"));
+                 if  (++i>=nbvb) {
+                         throw ErrorException(__FUNCT__,exprintf("ReconstructExistingMesh: All the vertices are aligned"));
                  }
-                Exchange( ordre[2], ordre[i]);
+                Exchange(ordre[2], ordre[i]);
 
                 Vertex *  v0=ordre[0], *v1=ordre[1];
 
-
-                triangles[0](0) = 0; // sommet pour infini 
+                triangles[0](0) = NULL; // Infinite vertex
                 triangles[0](1) = v0;
                 triangles[0](2) = v1;
 
-                triangles[1](0) = 0;// sommet pour infini 
+                triangles[1](0) = NULL;// Infinite vertex
                 triangles[1](2) = v0;
                 triangles[1](1) = v1;
@@ -2745 +2765 @@
                 triangles[0].SetAdj2(e2, &triangles[1] ,e1);
 
-                triangles[0].det = -1;  // faux triangles
-                triangles[1].det = -1;  // faux triangles
+                triangles[0].det = -1;  // boundary triangles
+                triangles[1].det = -1;  // boundary triangles
 
                 triangles[0].SetTriangleContainingTheVertex();
@@ -2881 +2901 @@
                 SetVertexFieldOn();
 
-                for (i=0;i<nbe;i++)
-                 if(edges[i].onGeometry) 
-                  for(int j=0;j<2;j++)
-                        if (!edges[i].adj[j])
-                         if(!edges[i][j].onGeometry->IsRequiredVertex()) {
-                                 throw ErrorException(__FUNCT__,exprintf("adj and vertex required esge(?)"));
-                         }
+                /*Check requirements consistency*/
+                for (i=0;i<nbe;i++){
+                        /*If the current mesh edge is on Geometry*/
+                        if(edges[i].onGeometry){
+                                for(int j=0;j<2;j++){
+                                        if (!edges[i].adj[j]){
+                                                if(!edges[i][j].onGeometry->IsRequiredVertex()) {
+                                                        printf("Error: adj and vertex requires edges [%i %i]=%i on = %i\n",i,j,Number(edges[i][j]),Gh.Number(edges[i].onGeometry));
+                                                        if (edges[i][j].onGeometry->OnGeomVertex())
+                                                         printf(" Vertex %i\n",Gh.Number(edges[i][j].onGeometry->gv));
+                                                        else if (edges[i][j].onGeometry->OnGeomEdge())
+                                                         printf(" Edges %i \n",Gh.Number(edges[i][j].onGeometry->ge));
+                                                        else
+                                                         printf(" = %p\n",edges[i][j].onGeometry);
+                                                        throw ErrorException(__FUNCT__,exprintf("See above (might be cryptic...)"));
+                                                }
+                                        }
+                                }
+                        }
+                }
         }
         /*}}}1*/
@@ -5271 +5304 @@
                 //  ReMakeTriangleContainingTheVertex();
 
-                GenerateMeshProperties();
+                ReconstructExistingMesh();
 
                 if (verbosity>2){

issm/trunk/src/c/Bamgx/objects/Triangles.h

@@ -68 +68 @@
 
                         //Constructors/Destructors
-                        Triangles(BamgMesh* bamgmesh,BamgOpts* bamgopts);
+                        Triangles(BamgGeom* bamggeom,BamgMesh* bamgmesh,BamgOpts* bamgopts);
                         Triangles(double* index,double* x,double* y,int nods,int nels);
                         Triangles(Triangles &,Geometry * pGh=0,Triangles* pBTh=0,long nbvxx=0 ); //copy operator
@@ -145 +145 @@
                         int  UnCrack();
                         void BuildGeometryFromMesh(BamgOpts* bamgopts=NULL);
-                        void GenerateMeshProperties() ;
+                        void ReconstructExistingMesh();
                         int  CrackMesh();
 

issm/trunk/src/c/Bamgx/objects/VertexOnGeom.h

@@ -24 +24 @@
                         double abscisse;  
                         union{ 
-                                GeometricalVertex * gv; // if abscisse <0; 
-                                GeometricalEdge * ge;  // if abscisse in [0..1]
+                                GeometricalVertex* gv; // if abscisse <0; 
+                                GeometricalEdge*   ge; // if abscisse in [0..1]
                         };
 
@@ -46 +46 @@
 
                         //Inline methods
-                        //inline void Set(const Triangles &,long,Triangles &);
                         void Set(const VertexOnGeom&,const Triangles &,Triangles &);