Changeset 22086

Timestamp:

09/15/17 18:50:36 (8 years ago)

Author:

Eric.Larour

Message:

CHG: new bamg code

Location:

issm/branches/trunk-larour-NatGeoScience2016/src

Files:

• c/modules/Bamgx/Bamgx.cpp (modified) (1 diff)
• c/modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp (modified) (5 diffs)
• m/modules/ContourToMesh.m (modified) (2 diffs)
• m/modules/ContourToNodes.m (modified) (2 diffs)
• m/modules/MeshProfileIntersection.m (modified) (2 diffs)
• wrappers/matlab/io/FetchMatlabData.cpp (modified) (5 diffs)

issm/branches/trunk-larour-NatGeoScience2016/src/c/modules/Bamgx/Bamgx.cpp

@@ -149 +149 @@
                 }
 
-                //Add geometry metric if provided
-                if(bamgopts->geometricalmetric) BTh.AddGeometryMetric(bamgopts);
-
                 //Smoothe metric
                 BTh.SmoothMetric(bamgopts->gradation);

issm/branches/trunk-larour-NatGeoScience2016/src/c/modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp

@@ -24 +24 @@
         R2     r;
         I2     I;
-        int    i,j;
+        int    i,j,k;
         int    it;
         int    i0,i1,i2;
         double areacoord[3];
         double aa,bb;
-        Icoor2 dete[3];
+        long long dete[3];
 
         /*Checks*/
@@ -66 +66 @@
         }
 
+        /*Create Single vertex to element connectivity*/
+        int* connectivity = xNew<int>(nods_data);
+        for(i=0;i<nels_data;i++){
+                for(j=0;j<3;j++){
+                        k = index_data[i*3+j]-1;
+                        _assert_(k>=0 & k<nods_data);
+                        connectivity[k]=i;
+                }
+        }
+
         /*Loop over output nodes*/
         for(i=0;i<N_interp;i++){
@@ -88 +98 @@
 
                         /*Area coordinates*/
-                        areacoord[0]= (double) dete[0]/tb.det;
-                        areacoord[1]= (double) dete[1]/tb.det;
-                        areacoord[2]= (double) dete[2]/tb.det;
+                        areacoord[0]= reCast<double>(dete[0])/reCast<double>(tb.det);
+                        areacoord[1]= reCast<double>(dete[1])/reCast<double>(tb.det);
+                        areacoord[2]= reCast<double>(dete[2])/reCast<double>(tb.det);
                         /*3 vertices of the triangle*/
                         i0=Th->GetId(tb[0]);
@@ -137 +147 @@
                         /*If we fall outside of the convex or outside of the mesh, return NaN*/
                         if(tb.det<0 || reft[it]<0){
-                                for (j=0;j<N_data;j++){
-                                        data_interp[i*N_data+j]=NAN;
-                                }
+                                _assert_(i0>=0 & i0<nods_data);
+                                it=connectivity[i0]; //or i1 or i2
+                                _assert_(it>=0 && it<nels_data);
+                                for(j=0;j<N_data;j++) data_interp[i*N_data+j]=data[N_data*it+j];
                         }
                         else{
+                                /*Inside the mesh!*/
                                 if(it<0 || it>=nels_data){
                                         _error_("Triangle number " << it << " not in [0 " << nels_data
                                                                 << "], report bug to developers (interpolation point: " <<x_interp[i]<<" "<<y_interp[i]<<")");
                                 }
-                                for (j=0;j<N_data;j++){
-
-                                        data_interp[i*N_data+j]=data[N_data*it+j];
-                                }
+                                for (j=0;j<N_data;j++) data_interp[i*N_data+j]=data[N_data*it+j];
                         }
                 }
@@ -158 +167 @@
         delete Th;
         xDelete<long>(reft);
+        xDelete<int>(connectivity);
         *pdata_interp=data_interp;
         return 1;

issm/branches/trunk-larour-NatGeoScience2016/src/m/modules/ContourToMesh.m

@@ -6 +6 @@
 %       
 %   index,x,y: mesh triangulation
-%   contourname: name of .exp file containing the contours.
+%   contourname: name of .exp or .shp file containing the contours.
 %   interptype: string defining type of interpolation ('element', or 'node').
 %   edgevalue: integer (0, 1, or 2) defining the value associated to the nodes on the edges of the polygons.
@@ -25 +25 @@
 end
 
+%Some conversion of files: 
+[path,name,ext]=fileparts(contourname); 
+if strcmpi(ext,'.shp'),
+        %read contour from shapefile:
+        contour=shpread(contourname); 
+
+        %write it to a temporary filename: 
+        contourname=tempname; 
+        expwrite(contour,contourname);
+end
+
 %Call mex module
 [in_nod,in_elem] = ContourToMesh_matlab(index,x,y,contourname,interptype,edgevalue);

issm/branches/trunk-larour-NatGeoScience2016/src/m/modules/ContourToNodes.m

@@ -6 +6 @@
 %
 %   x,y: list of nodes
-%   contourname: name of .exp file containing the contours, or resulting structure from call to expread
+%   contourname: name of Argus or Shape file containing the contours, or resulting structure from call to expread
 %   edgevalue: integer (0, 1 or 2) defining the value associated to the nodes on the edges of the polygons
 %   flags: vector of flags (0 or 1), of size nodes
@@ -16 +16 @@
 end
 
-%Some conversion of files: 
-[path,name,ext]=fileparts(contourname); 
-if strcmpi(ext,'.shp'),
-        %read contour from shapefile:
-        contour=shpread(contourname); 
+%Some conversion of files:  
+if ischar(contourname),
+        [path,name,ext]=fileparts(contourname); 
+        if strcmpi(ext,'.shp'),
+                %read contour from shapefile:
+                contour=shpread(contourname); 
 
-        %write it to a temporary filename: 
-        contourname=tempname; 
-        expwrite(contour,contourname);
+                %write it to a temporary filename: 
+                contourname=tempname; 
+                expwrite(contour,contourname);
+        end
 end
 

issm/branches/trunk-larour-NatGeoScience2016/src/m/modules/MeshProfileIntersection.m

@@ -7 +7 @@
 %          input:
 %               index,x,y is a triangulation
-%               filename: name of Argus style .exp file containing the segments (can be groups of disconnected segments)
+%               filename: name of Argus or Shape file containing the segments (can be groups of disconnected segments)
 %          output:
 %               segments: array made of x1,y1,x2,y2,element_id lines (x1,y1) and (x2,y2) are segment extremities for a segment 
@@ -19 +19 @@
 end
 
+[path,name,ext]=fileparts(filename); 
+if strcmpi(ext,'.shp'),
+        
+        %convert to expfile and store in a temporary directory: 
+        oldfilename=filename;
+        filename=[tempname '.exp'];
+        shp2exp(oldfilename,filename);
+end
+
 % Call mex module
 [segments] = MeshProfileIntersection_matlab(index,x,y,filename);

issm/branches/trunk-larour-NatGeoScience2016/src/wrappers/matlab/io/FetchMatlabData.cpp

@@ -537 +537 @@
         FetchData(&bamgopts->gradation,mxGetField(dataref,0,"gradation"));
         FetchData(&bamgopts->Hessiantype,mxGetField(dataref,0,"Hessiantype"));
-        FetchData(&bamgopts->MaxCornerAngle,mxGetField(dataref,0,"MaxCornerAngle"));
         FetchData(&bamgopts->maxnbv,mxGetField(dataref,0,"maxnbv"));
         FetchData(&bamgopts->maxsubdiv,mxGetField(dataref,0,"maxsubdiv"));
@@ -545 +544 @@
         FetchData(&bamgopts->omega,mxGetField(dataref,0,"omega"));
         FetchData(&bamgopts->power,mxGetField(dataref,0,"power"));
-        FetchData(&bamgopts->random,mxGetField(dataref,0,"random"));
         FetchData(&bamgopts->verbose,mxGetField(dataref,0,"verbose"));
 
         FetchData(&bamgopts->Crack,mxGetField(dataref,0,"Crack"));
-        FetchData(&bamgopts->geometricalmetric,mxGetField(dataref,0,"geometricalmetric"));
         FetchData(&bamgopts->KeepVertices,mxGetField(dataref,0,"KeepVertices"));
         FetchData(&bamgopts->splitcorners,mxGetField(dataref,0,"splitcorners"));
@@ -999 +996 @@
 mxArray* mxGetAssignedField(const mxArray* pmxa_array,int number,const char* field){
 
-        //output
-        mxArray* mxfield=NULL;
-
-        //input
-        mxArray    *inputs[2];
-        mxArray    *pindex      = NULL;
-        const char *fnames[2];
-        mwSize      ndim        = 2;
-        mwSize      onebyone[2] = {1,1};
-
-        //We want to call the subsasgn method, and get the returned array.This ensures that if we are running 
-        //large sized problems, the data is truly loaded from disk by the model subsasgn class method.
-        inputs[0]=(mxArray*)pmxa_array; //this is the model
-
-        //create index structure used in the assignment (index.type='.' and index.subs='x' for field x for ex)
-        fnames[0] = "type";
-        fnames[1] = "subs";
-        pindex=mxCreateStructArray( ndim,onebyone,2,fnames);
-        mxSetField( pindex, 0, "type",mxCreateString("."));
-        mxSetField( pindex, 0, "subs",mxCreateString(field));
-        inputs[1]=pindex;
-
-        mexCallMATLAB( 1, &mxfield, 2, (mxArray**)inputs, "subsref");
+        /*Output*/
+        mxArray *mxfield = NULL;
+
+        if(mxIsStruct(pmxa_array)){
+                mxfield = mxGetField(pmxa_array,number,field);
+        }
+        else{
+                /*This is an object, mxGetField returns NULL in old version of matlab (we do not have access to them)*/
+
+                /*Intermediaries*/
+                mxArray    *inputs[2];
+                mwSize      ndim        = 2;
+                mwSize      onebyone[2] = {1,1};
+
+                /*create index structure used in the assignment (index.type='.' and index.subs='x' for field x*/
+                const char *fnames[2];
+                fnames[0] = "type"; fnames[1] = "subs";
+                mxArray* pindex=mxCreateStructArray( ndim,onebyone,2,fnames);
+                mxSetField( pindex, 0, "type",mxCreateString("."));
+                mxSetField( pindex, 0, "subs",mxCreateString(field));
+                inputs[0]=(mxArray*)pmxa_array; //this is the model
+                inputs[1]=pindex;
+
+                mexCallMATLAB( 1, &mxfield, 2, (mxArray**)inputs, "subsref");
+        }
+
+        if(mxfield == NULL) _error_("Could not find field "<< field <<" in structure");
 
         return mxfield;
@@ -1091 +1092 @@
 
         return(ochar);
-}/*}}}*/
-GenericOption<Options**>* OptionStructParse( char* name, const mxArray* prhs[]){ /*{{{*/
-
-        int            i;
-        char           namei[161];
-        Option*                   option      = NULL;
-        GenericOption<Options**>  *ostruct    = NULL;
-        const mwSize  *ipt        = NULL;
-        const mxArray *structi;
-        mwIndex        sindex;
-
-        /*check and parse the name  */
-        ostruct=new GenericOption<Options**>();
-        ostruct->name =xNew<char>(strlen(name)+1);
-        memcpy(ostruct->name,name,(strlen(name)+1)*sizeof(char));
-
-        /*check and parse the value  */
-        if (!mxIsClass(prhs[0],"struct")){
-                _error_("Value of option \"" << ostruct->name  << "\" must be class \"struct\", not class \"" << mxGetClassName(prhs[0]) <<"\".");
-        }
-        ostruct->numel=mxGetNumberOfElements(prhs[0]);
-        ostruct->ndims=mxGetNumberOfDimensions(prhs[0]);
-        ipt           =mxGetDimensions(prhs[0]);
-        ostruct->size =xNew<int>(ostruct->ndims);
-        for (i=0; i<ostruct->ndims; i++) ostruct->size[i]=(int)ipt[i];
-        if (ostruct->numel) ostruct->value=xNew<Options*>(ostruct->numel);
-
-        /*loop through and process each element of the struct array  */
-        for (sindex=0; sindex<ostruct->numel; sindex++) {
-                ostruct->value[sindex]=new Options;
-
-                /*loop through and process each field for the element  */
-                for (i=0; i<mxGetNumberOfFields(prhs[0]); i++) {
-                        sprintf(namei,"%s.%s",name,mxGetFieldNameByNumber(prhs[0],i));
-                        structi=mxGetFieldByNumber(prhs[0],sindex,i);
-
-                        option=(Option*)OptionParse(namei,&structi);
-                        ostruct->value[sindex]->AddObject((Object*)option);
-                        option=NULL;
-                }
-        }
-
-        return(ostruct);
-}/*}}}*/
-GenericOption<Options*>* OptionCellParse( char* name, const mxArray* prhs[]){ /*{{{*/
-
-        int            i;
-        int           *dims;
-        char           namei[161];
-        char           cstr[81];
-        GenericOption<Options*> *ocell      = NULL;
-        Option        *option     = NULL;
-        const mwSize  *ipt        = NULL;
-        const mxArray *celli;
-        mwIndex        cindex;
-
-        /*check and parse the name  */
-        ocell=new GenericOption<Options*>();
-        ocell->name =xNew<char>(strlen(name)+1);
-        memcpy(ocell->name,name,(strlen(name)+1)*sizeof(char));
-
-        /*check and parse the value  */
-        if (!mxIsClass(prhs[0],"cell")){
-                _error_("Value of option \"" << ocell->name  << "\" must be class \"cell\", not class \"" << mxGetClassName(prhs[0]) <<"\".");
-        }
-
-        ocell->numel=mxGetNumberOfElements(prhs[0]);
-        ocell->ndims=mxGetNumberOfDimensions(prhs[0]);
-        ipt         =mxGetDimensions(prhs[0]);
-        ocell->size =xNew<int>(ocell->ndims);
-        for (i=0; i<ocell->ndims; i++) ocell->size[i]=(int)ipt[i];
-        ocell->value=new Options;
-
-        /*loop through and process each element of the cell array  */
-        dims=xNew<int>(ocell->ndims);
-        for (cindex=0; cindex<ocell->numel; cindex++) {
-                ColumnWiseDimsFromIndex(dims,(int)cindex,ocell->size,ocell->ndims);
-                StringFromDims(cstr,dims,ocell->ndims);
-                #ifdef _INTEL_WIN_
-                        _snprintf(namei,161,"%s%s",name,cstr);
-                #else
-                        snprintf(namei,161,"%s%s",name,cstr);
-                #endif
-                celli=mxGetCell(prhs[0],cindex);
-
-                option=(Option*)OptionParse(namei,&celli);
-                ocell->value->AddObject((Object*)option);
-                option=NULL;
-        }
-        xDelete<int>(dims);
-
-        return(ocell);
 }/*}}}*/
 Option* OptionParse(char* name, const mxArray* prhs[]){ /*{{{*/
@@ -1198 +1107 @@
         else if(mxIsClass(prhs[0],"char"))
          option=(Option*)OptionCharParse(name,prhs);
-        else if(mxIsClass(prhs[0],"struct"))
-         option=(Option*)OptionStructParse(name,prhs);
-        else if(mxIsClass(prhs[0],"cell"))
-         option=(Option*)OptionCellParse(name,prhs);
         else {
-                _printf0_("  Converting value of option \"" << name << "\" from unrecognized class \"" << mxGetClassName(prhs[0]) << "\" to class \"" << "struct" << "\".\n");
-                if (!mexCallMATLAB(1,lhs,1,(mxArray**)prhs,"struct")) {
-                        option=(Option*)OptionStructParse(name,(const mxArray**)lhs);
-                        mxDestroyArray(lhs[0]);
-                }
-                else _error_("Second argument value of option \""<< name <<"\" is of unrecognized class \""<< mxGetClassName(prhs[0]) <<"\".");
+                _error_("Second argument value of option \""<< name <<"\" is of unrecognized class \""<< mxGetClassName(prhs[0]) <<"\".");
         }