Changeset 5357

Timestamp:

08/18/10 10:13:38 (15 years ago)

Author:

Eric.Larour

Message:

Extended TriaSearch to deal with several grids

Location:

issm/trunk

Files:

• src/c/modules/TriaSearchx/TriaSearchx.cpp (modified) (2 diffs)
• src/c/modules/TriaSearchx/TriaSearchx.h (modified) (1 diff)
• src/m/classes/public/ismodelselfconsistent.m (modified) (1 diff)
• src/m/classes/public/solveparallel.m (modified) (1 diff)
• src/m/utils/Exp/flowlines.m (modified) (1 diff)
• src/mex/TriaSearch/TriaSearch.cpp (modified) (4 diffs)
• test/Archives/Archive235.mat (added)
• test/Archives/Archive236.mat (added)

issm/trunk/src/c/modules/TriaSearchx/TriaSearchx.cpp

@@ -13 +13 @@
 using namespace std;
 
-void TriaSearchx(double* ptria,double* index,int nel, double* x, double* y, int nods,double x0, double y0){
+void TriaSearchx(double** ptria,double* index,int nel, double* x, double* y, int nods,double* x0, double* y0,int numberofgrids){
 
         /*Output*/
-        double tria;
+        double* tria=NULL;
+
+        /*allocate: */
+        tria=(double*)xmalloc(numberofgrids*sizeof(double));
 
         /*Intermediary*/
@@ -34 +37 @@
         Th.CreateSingleVertexToTriangleConnectivity();
 
-        //Get current point coordinates
-        r.x=x0; r.y=y0;
-        I=Th.R2ToI2(r);
+        for(i=0;i<numberofgrids;i++){
 
-        //Find triangle holding r/I
-        Triangle &tb=*Th.TriangleFindFromCoord(I,dete);
+                //Get current point coordinates
+                r.x=x0[i]; r.y=y0[i];
+                
+                I=Th.R2ToI2(r);
 
-        // internal point 
-        if (tb.det>0){ 
-                //Area coordinate
-                areacoord[0]= (double) dete[0]/ tb.det;
-                areacoord[1]= (double) dete[1] / tb.det;
-                areacoord[2]= (double) dete[2] / tb.det;
-                //3 vertices of the triangle
-                i0=Th.GetId(tb[0]);
-                i1=Th.GetId(tb[1]);
-                i2=Th.GetId(tb[2]);
-                //triangle number
-                tria=(double)Th.GetId(tb);
+                //Find triangle holding r/I
+                Triangle &tb=*Th.TriangleFindFromCoord(I,dete);
+
+                // internal point 
+                if (tb.det>0){ 
+                        //Area coordinate
+                        areacoord[0]= (double) dete[0]/ tb.det;
+                        areacoord[1]= (double) dete[1] / tb.det;
+                        areacoord[2]= (double) dete[2] / tb.det;
+                        //3 vertices of the triangle
+                        i0=Th.GetId(tb[0]);
+                        i1=Th.GetId(tb[1]);
+                        i2=Th.GetId(tb[2]);
+                        //triangle number
+                        tria[i]=(double)Th.GetId(tb);
+                }
+                //external point
+                else tria[i]=NAN;
         }
-        //external point
-        else tria=NAN;
-        
+
         /*Assign output pointers:*/
         *ptria=tria;

issm/trunk/src/c/modules/TriaSearchx/TriaSearchx.h

@@ -9 +9 @@
 
 /* local prototypes: */
-void TriaSearchx(double* ptria,double* index,int nel, double* x, double* y, int nods,double x0, double y0);
+void TriaSearchx(double** ptria,double* index,int nel, double* x, double* y, int nods,double* x0, double* y0,int numberofgrids);
 
 #endif

issm/trunk/src/m/classes/public/ismodelselfconsistent.m

@@ -351 +351 @@
 
         %LENGTH CONTROL FIELDS
-        fields={'maxiter','optscal','cm_responses','cm_jump'};
+        %fields={'maxiter','optscal','cm_responses','cm_jump'};
+        fields={'maxiter','optscal','cm_jump'};
         checksize(md,fields,[md.nsteps 1]);
 

issm/trunk/src/m/classes/public/solveparallel.m

@@ -42 +42 @@
         end
 end
+
+%post processes qmu results if necessary
+if md.qmu_analysis,
+        system(['rm -rf qmu' num2str(GetPId)]);
+end
+
+

issm/trunk/src/m/utils/Exp/flowlines.m

@@ -38 +38 @@
 
 %check seed points
-tria=tsearch(x,y,index,x0,y0);
-pos=find(isnan(tria));
-%seems to be a bug here: must do the test several times...
-while ~isempty(pos)
-        x0(pos)=[];
-        y0(pos)=[];
-        tria=tsearch(x,y,index,x0,y0);
-        pos=find(isnan(tria));
+tria=TriaSearch(index,x,y,x0,y0);
+if isnan(tria),
+        error('flowlines error message: seed point is not inside mesh!');
 end
 

issm/trunk/src/mex/TriaSearch/TriaSearch.cpp

@@ -14 +14 @@
 
         double* x=NULL;
+        double* y=NULL;
         int     nods;
 
-        double* y=NULL;
+        double* x0=NULL;
+        double* y0=NULL;
+        int     numberofgrids;
 
-        double  x0,y0;
-        
         /* output: */
-        double  tria;
+        double*  tria=NULL;
 
         /*Boot module: */
@@ -33 +34 @@
         FetchData(&x,&nods,XHANDLE);
         FetchData(&y,&nods,YHANDLE);
-        FetchData(&x0,X0HANDLE);
-        FetchData(&y0,Y0HANDLE);
+        FetchData(&x0,&numberofgrids,X0HANDLE);
+        FetchData(&y0,&numberofgrids,Y0HANDLE);
 
         /* Echo: {{{1*/
@@ -41 +42 @@
 
         /* Run core computations: */
-        TriaSearchx(&tria,index,nel,x,y,nods,x0,y0);
+        TriaSearchx(&tria,index,nel,x,y,nods,x0,y0,numberofgrids);
 
         /* c to matlab: */
-        tria++;
+        for(i=0;i<numberofgrids;i++)tria[i]++;
 
         /*Write data: */
-        WriteData(TRIA,tria);
+        WriteData(TRIA,tria,numberofgrids);
 
         /*end module: */
@@ -61 +62 @@
         _printf_("      index,x,y: mesh triangulatrion\n");
         _printf_("      x0,y0: coordinates of the point for which we are trying to find a triangle\n");
+        _printf_("      x0,y0 can be an array of points\n");
         _printf_("\n");
 }