Changeset 3126

Timestamp:

02/24/10 15:57:47 (15 years ago)

Author:

Mathieu Morlighem

Message:

Increased InterpFromMeshToMesh2dx speed using Bamg

Location:

issm/trunk/src/c

Files:

• Bamgx/Bamgx.cpp (modified) (1 diff)
• Bamgx/Mesh2.h (modified) (3 diffs)
• Bamgx/objects/Triangles.cpp (modified) (4 diffs)
• Bamgx/objects/Vertex.cpp (modified) (1 diff)
• InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp (modified) (1 diff)
• InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.h (modified) (1 diff)
• InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dxt.cpp (deleted)
• Makefile.am (modified) (2 diffs)
• issm.h (modified) (1 diff)

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

@@ -166 +166 @@
                 Th.ReNumberingTheTriangleBySubDomain();
                 if(NbSmooth>0) Th.SmoothingVertex(NbSmooth,omega);
-                Th.BTh.ReMakeTriangleContainingTheVertex();
 
                 //display info

issm/trunk/src/c/Bamgx/Mesh2.h

@@ -648 +648 @@
                         Triangles(const char * ,Real8=-1) ;
                         Triangles(BamgMesh* bamgmesh,BamgOpts* bamgopts);
+                        Triangles(double* index,double* x,double* y,int nods,int nels);
 
                         Triangles(Int4 nbvx,Triangles & BT,int keepBackVertices=1) :Gh(BT.Gh),BTh(BT) {
@@ -724 +725 @@
                         Triangle * FindTriangleContaining(const I2 & ,Icoor2 [3],Triangle *tstart=0) const;
 
+                        void ReadMesh(double* index,double* x,double* y,int nods,int nels);
                         void ReadMesh(BamgMesh* bamgmesh, BamgOpts* bamgopts);
                         void WriteMesh(BamgMesh* bamgmesh,BamgOpts* bamgopts);
@@ -738 +740 @@
                         int UnCrack();
 
-                        void BuildGeometryFromMesh(BamgOpts* bamgopts);
+                        void BuildGeometryFromMesh(BamgOpts* bamgopts=NULL);
                         void FillHoleInMesh() ;
                         int  CrackMesh();

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

@@ -27 +27 @@
                 PreInit(0);
                 ReadMesh(bamgmesh,bamgopts);
+                SetIntCoor();
+                FillHoleInMesh();
+        }
+        /*}}}1*/
+        /*FUNCTION Triangles::Triangles(double* index,double* x,double* y,int nods,int nels){{{1*/
+        Triangles::Triangles(double* index,double* x,double* y,int nods,int nels):Gh(*(new Geometry())),BTh(*this){
+
+                PreInit(0);
+                ReadMesh(index,x,y,nods,nels);
                 SetIntCoor();
                 FillHoleInMesh();
@@ -226 +235 @@
 
         /*IO*/
-        /*FUNCTION Triangles::ReadMesh{{{1*/
+        /*FUNCTION Triangles::ReadMesh(double* index,double* x,double* y,int nods,int nels){{{1*/
+        void Triangles::ReadMesh(double* index,double* x,double* y,int nods,int nels){
+
+                Real8 Hmin = HUGE_VAL;// the infinie value 
+                Int4 i1,i2,i3,iref;
+                Int4 i,j;
+                Int4 hvertices =0;
+                Metric M1(1);
+                int Version=1,dim=2;
+                int verbose=0;
+
+                nbv=nods;
+                nbvx=nbv;
+                nbt=nels;
+                nbiv=nbvx;
+
+                //Vertices
+                if (verbose) printf("Reading vertices (%i)\n",nbv);
+                vertices=(Vertex*)xmalloc(nbv*sizeof(Vertex));
+                ordre=(Vertex**)xmalloc(nbv*sizeof(Vertex*));
+                for (i=0;i<nbv;i++){
+                        vertices[i].r.x=x[i];
+                        vertices[i].r.y=y[i];
+                        vertices[i].ReferenceNumber=1;
+                        vertices[i].DirOfSearch =NoDirOfSearch;
+                        vertices[i].m=M1;
+                }
+                nbtx=2*nbvx-2; // for filling The Holes and quadrilaterals 
+
+                //Triangles
+                if (verbose) printf("Reading triangles (%i)\n",nbt);
+                triangles =new Triangle[nbtx]; //we cannot allocate only nbt triangles since 
+                //other triangles will be added for each edge
+                for (i=0;i<nbt;i++){
+                        Triangle & t = triangles[i];
+                        i1=(Int4)index[i*3+0]-1; //for C indexing
+                        i2=(Int4)index[i*3+1]-1; //for C indexing
+                        i3=(Int4)index[i*3+2]-1; //for C indexing
+                        t=Triangle(this,i1,i2,i3);
+                        t.color=i;
+                }
+
+                /*Recreate geometry: */
+                if (verbose) printf("Building Geometry\n");
+                BuildGeometryFromMesh();
+                if (verbose) printf("Completing geometry\n");
+                Gh.AfterRead();
+        }
+        /*}}}1*/
+        /*FUNCTION Triangles::ReadMesh(BamgMesh* bamgmesh, BamgOpts* bamgopts){{{1*/
         void Triangles::ReadMesh(BamgMesh* bamgmesh, BamgOpts* bamgopts){
 
@@ -997 +1055 @@
                 /*Intermediary*/
                 int i,j,k,kk,it,jt;
+                int    verbosity=0;
+                double cutoffradian=10*Pi/180;
 
                 /*Recover options*/
-                int    verbosity=bamgopts->verbose;
-                double cutoffradian=bamgopts->MaximalAngleOfCorner*Pi/180;
+                if (bamgopts){
+                        verbosity=bamgopts->verbose;
+                        cutoffradian=bamgopts->MaximalAngleOfCorner*Pi/180;
+                }
 
                 //display info
@@ -2716 +2778 @@
                         t = a->t;
                         if (t<triangles || t>=triangles+nbt){
+                                a->Echo();
                                 throw ErrorException(__FUNCT__,exprintf("t<triangles || t>=triangles+nbt"));
                         }

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

@@ -207 +207 @@
                 printf("  Euclidean coordinates r.x: %g, r.y: %g\n",r.x,r.y);
                 printf("  ReferenceNumber = %i\n",ReferenceNumber);
+                m.Echo();
 
                 return;

issm/trunk/src/c/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp

@@ -1 @@
-/*!\file:  InterpFromMeshToMesh2dx.cpp
- * \brief  "c" core code for interpolating values from a structured grid.
- */ 
-
-
-#ifdef HAVE_CONFIG_H
-        #include "config.h"
-#else
-#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
-#endif
+/*!\file InterpFromMeshToMesh2dx
+ */
 
 #include "./InterpFromMeshToMesh2dx.h"
-#include "../shared/shared.h"
 
 #undef __FUNCT__ 
 #define __FUNCT__ "InterpFromMeshToMesh2dx"
-int InterpFromMeshToMesh2dx( Vec* pdata_prime,double* index_data, double* x_data, double* y_data, int nods_data,int nels_data, double* data, int data_length, double* x_prime, double* y_prime, int nods_prime,double default_value){
+
+#include "../shared/shared.h"
+#include "../include/macros.h"
+#include "../toolkits/toolkits.h"
+
+/*Bamg: */
+#include <unistd.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include "../Bamgx/Mesh2.h"
+#include "../Bamgx/QuadTree.h"
+
+using namespace bamg;
+using namespace std;
+
+int InterpFromMeshToMesh2dx(double** pdata_interp,double* index_data,double* x_data,double* y_data,int nods_data,int nels_data,
+                        double* data,int data_rows,int data_cols,double* x_interp,double* y_interp,int nods_interp,double default_value){
 
         /*Output*/
-        Vec data_prime=NULL;
+        double* data_interp=NULL;
+        int noerr;
 
         /*Intermediary*/
-        int i;
-        int interpolation_type;
-        bool debug;
-        double x_prime_min,x_prime_max;
-        double y_prime_min,y_prime_max;
+        R2     r;
+        I2     I;
+        int    i,j,k;
+        int    it;
+        int    i0,i1,i2;
+        double areacoord[3];
+        double aa,bb;
+        double data_value;
+        Icoor2 dete[3];
+        int verbose=0;
 
-        /*threading: */
-        InterpFromMeshToMesh2dxThreadStruct gate;
-        int num=1;
-        #ifdef _MULTITHREADING_
-        num=_NUMTHREADS_;
-        #endif
-
-        /*some checks*/
-        if (nels_data<1 || nods_data<3 || nods_prime==0){
-                throw ErrorException(__FUNCT__,"nothing to be done according to the mesh given in input");
+        /*Checks*/
+        if (data_cols<=0){
+                throw ErrorException(__FUNCT__,exprintf("data provided has a negative number of columns"));
         }
-
-        /*Set debug to 1 if there are lots of elements*/
-        debug=(bool)((double)nels_data*(double)nods_prime >= pow((double)10,(double)9));
-
-        /*figure out what kind of interpolation is needed*/
-        if (data_length==nods_data){
-                interpolation_type=1;
-        }
-        else if (data_length==nels_data){
-                interpolation_type=2;
-        }
-        else{
-                throw ErrorException(__FUNCT__,"length of vector data not supported yet. It should be of length (number of nodes) or (number of elements)!");
-        }
-
-        /*Get prime mesh extrema coordinates*/
-        x_prime_min=x_prime[0]; x_prime_max=x_prime[0];y_prime_min=y_prime[0]; y_prime_max=y_prime[0];
-        for (i=1;i<nods_prime;i++){
-                if (x_prime[i]<x_prime_min) x_prime_min=x_prime[i];
-                if (x_prime[i]>x_prime_max) x_prime_max=x_prime[i];
-                if (y_prime[i]<y_prime_min) y_prime_min=y_prime[i];
-                if (y_prime[i]>y_prime_max) y_prime_max=y_prime[i];
+        if (data_rows!=nods_data && data_rows!=nels_data){
+                throw ErrorException(__FUNCT__,exprintf("data provided should have either %i or %i lines (not %i)",nods_data,nels_data,data_rows));
         }
 
         /*Initialize output*/
-        data_prime=NewVec(nods_prime);
-        for (i=0;i<nods_prime;i++) VecSetValue(data_prime,i,default_value,INSERT_VALUES);
+        if (verbose) printf("Initializing output vector\n");
+        data_interp=(double*)xmalloc(nods_interp*data_cols*sizeof(double));
 
-        /*initialize thread parameters: */
-        gate.debug=debug;
-        gate.nods_prime=nods_prime;
-        gate.nels_data=nels_data;
-        gate.x_data=x_data;
-        gate.y_data=y_data;
-        gate.index_data=index_data;
-        gate.x_prime=x_prime;
-        gate.y_prime=y_prime;
-        gate.data=data;
-        gate.interpolation_type=interpolation_type;
-        gate.default_value=default_value;
-        gate.data_prime=data_prime;
-        gate.x_prime_min=x_prime_min;
-        gate.x_prime_max=x_prime_max;
-        gate.y_prime_min=y_prime_min;
-        gate.y_prime_max=y_prime_max;
+        // read background mesh 
+        if (verbose) printf("Reading mesh\n");
+        Triangles Th(index_data,x_data,y_data,nods_data,nels_data); 
+        Th.ReMakeTriangleContainingTheVertex();
 
-        /*launch the thread manager with InterpFromMeshToMesh2dxt as a core: */
-        LaunchThread(InterpFromMeshToMesh2dxt,(void*)&gate,num);
+        //Loop over output nodes
+        if (verbose) printf("Loop over the nodes\n");
+        for(i=0;i<nods_interp;i++){
+
+                //Get current point coordinates
+                r.x=x_interp[i]; r.y=y_interp[i];
+                I2 I=Th.toI2(r);
+
+                //Find triangle holding r/I
+                Triangle &tb=*Th.FindTriangleContaining(I,dete);
+
+                // internal point 
+                if (tb.det>0){ 
+                        //Area coordinate
+                        areacoord[0]= (Real8) dete[0]/ tb.det;
+                        areacoord[1]= (Real8) dete[1] / tb.det;
+                        areacoord[2]= (Real8) dete[2] / tb.det;
+                        //3 vertices of the triangle
+                        i0=Th.Number(tb[0]);
+                        i1=Th.Number(tb[1]);
+                        i2=Th.Number(tb[2]);
+                        //triangle number
+                        it=Th.Number(tb);
+                }
+
+                //external point
+                else {
+                        //Get closest adjacent triangle (inside the mesh)
+                        TriangleAdjacent ta=CloseBoundaryEdge(I,&tb,aa,bb).Adj();
+                        int k=ta;
+                        Triangle &tc=*(Triangle*)ta;
+                        //Area coordinate
+                        areacoord[VerticesOfTriangularEdge[k][1]] = aa;
+                        areacoord[VerticesOfTriangularEdge[k][0]] = bb;
+                        areacoord[OppositeVertex[k]] = 1 - aa -bb;
+                        //3 vertices of the triangle
+                        i0=Th.Number(tc[0]);
+                        i1=Th.Number(tc[1]);
+                        i2=Th.Number(tc[2]);
+                        //triangle number
+                        it=Th.Number(tc);
+                }
+
+                /*Last step, P1 interpolation*/
+                if (data_rows==nods_data){
+                        for (j=0;j<data_cols;j++){
+                                data_interp[i*data_cols+j]=areacoord[0]*data[data_cols*i0+j]+areacoord[1]*data[data_cols*i1+j]+areacoord[2]*data[data_cols*i2+j];
+                        }
+                }
+                else{
+                        for (j=0;j<data_cols;j++){
+                                if (it<0 || it>=nels_data){
+                                        throw ErrorException(__FUNCT__,exprintf("Triangle number %i not in [0 %i], because not correctly implemented yet... interpolate on grid first",it,nels_data));
+                                }
+                                data_interp[i*data_cols+j]=data[data_cols*it+j];
+                        }
+                }
+
+        }
 
         /*Assign output pointers:*/
-        *pdata_prime=data_prime;
+        if (verbose) printf("Assigning output\n");
+        *pdata_interp=data_interp;
+
+        /*No error return*/
+        return noerr;
+
 }
-
-
-

issm/trunk/src/c/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.h

@@ -1 @@
-/*!\file InterpFromMeshToMesh2dx.h
- * \brief: header file for Data interpolation routines.
- */
+/*!\file:  InterpFromMeshToMesh2dx.h
+ * \brief header file for Bamg module
+ */ 
 
 #ifndef _INTERPFROMMESHTOMESH2DX_H
 #define _INTERPFROMMESHTOMESH2DX_H
 
-#include "../toolkits/toolkits.h"
+#include "../objects/objects.h"
 
+/* local prototypes: */
+int InterpFromMeshToMesh2dx(double** pdata_interp,double* index_data,double* x_data,double* y_data,int nods_data,int nels_data,
+                        double* data,int data_rows,int data_cols,double* x_interp,double* y_interp,int nods_interp,double default_value);
 
-
-/*threading: */
-typedef struct{
-
-        int debug;
-        int nels_data;
-        double* x_data;
-        double* y_data;
-        double* index_data;
-        int nods_prime;
-        double* x_prime;
-        double* y_prime;
-        double* data;
-        int     interpolation_type;
-        double  default_value;
-        Vec     data_prime;
-        double  x_prime_min;
-        double  x_prime_max;
-        double  y_prime_min;
-        double  y_prime_max;
-
-} InterpFromMeshToMesh2dxThreadStruct;
-
-
-int InterpFromMeshToMesh2dx( Vec* pdata_prime,double* index_data, double* x_data, double* y_data, int nods_data,int nels_data, double* data, int data_length, double* x_prime, double* y_prime, int nods_prime,double default_value);
-void* InterpFromMeshToMesh2dxt(void* vInterpFromMeshToMesh2dxThreadStruct);
-
-
-#endif /* _INTERPFROMMESHTOMESH2DX_H */
-
+#endif

issm/trunk/src/c/Makefile.am

@@ -260 +260 @@
                                         ./InterpFromGridToMeshx/InterpFromGridToMeshxt.cpp\
                                         ./InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp\
-                                        ./InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dxt.cpp\
                                         ./InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.h\
                                         ./InterpFromMeshToMesh3dx/InterpFromMeshToMesh3dx.cpp\
@@ -587 +586 @@
                                         ./InterpFromGridToMeshx/InterpFromGridToMeshxt.cpp\
                                         ./InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp\
-                                        ./InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dxt.cpp\
                                         ./InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.h\
                                         ./InterpFromMeshToMesh3dx/InterpFromMeshToMesh3dx.cpp\

issm/trunk/src/c/issm.h

@@ -63 +63 @@
 #include "./MassFluxx/MassFluxx.h"
 #include "./Bamgx/Bamgx.h"
+#include "./Bamgx/InterpBamgx.h"