Index: /issm/trunk-jpl/src/c/modules/ContourToMeshx/ContourToMeshx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ContourToMeshx/ContourToMeshx.cpp (revision 12126) +++ /issm/trunk-jpl/src/c/modules/ContourToMeshx/ContourToMeshx.cpp (revision 12127) @@ -11,5 +11,5 @@ #include "./ContourToMeshx.h" -int ContourToMeshx( Vector** pin_nod,Vector** pin_elem, double* index, double* x, double* y,Contour** contours,int numcontours,char* interptype,int nel,int nods, int edgevalue) { +int ContourToMeshx( Vector** pin_nod,Vector** pin_elem, double* index, double* x, double* y,DataSet* contours,char* interptype,int nel,int nods, int edgevalue) { int noerr=1; @@ -37,5 +37,4 @@ /*initialize thread parameters: */ - gate.numcontours=numcontours; gate.contours=contours; gate.nods=nods; Index: /issm/trunk-jpl/src/c/modules/ContourToMeshx/ContourToMeshx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/ContourToMeshx/ContourToMeshx.h (revision 12126) +++ /issm/trunk-jpl/src/c/modules/ContourToMeshx/ContourToMeshx.h (revision 12127) @@ -13,6 +13,5 @@ typedef struct{ - int numcontours; - Contour** contours; + DataSet* contours; int nods; int edgevalue; @@ -25,5 +24,5 @@ /* local prototypes: */ -int ContourToMeshx( Vector** pin_nods,Vector** pin_elem, double* index, double* x, double* y,Contour** contours,int numcontours,char* interptype,int nel,int nods, int edgevalue); +int ContourToMeshx( Vector** pin_nods,Vector** pin_elem, double* index, double* x, double* y,DataSet* contours,char* interptype,int nel,int nods, int edgevalue); void* ContourToMeshxt(void* vContourToMeshxThreadStruct); Index: /issm/trunk-jpl/src/c/modules/ContourToMeshx/ContourToMeshxt.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ContourToMeshx/ContourToMeshxt.cpp (revision 12126) +++ /issm/trunk-jpl/src/c/modules/ContourToMeshx/ContourToMeshxt.cpp (revision 12127) @@ -26,13 +26,7 @@ /*Contour:*/ - Contour* contouri=NULL; - int numnodes; - double* xc=NULL; - double* yc=NULL; - + DataSet* contours=NULL; /*parameters: */ - int numcontours; - Contour** contours=NULL; int nods; int edgevalue; @@ -49,5 +43,4 @@ /*recover parameters :*/ - numcontours=gate->numcontours; contours=gate->contours; nods=gate->nods; @@ -61,10 +54,7 @@ /*Loop through all contours: */ - for (i=0;inods; - xc=contouri->x; - yc=contouri->y; - IsInPoly(in_nod,xc,yc,numnodes,x,y,i0,i1,edgevalue); + for (i=0;iSize();i++){ + Contour* contour=(Contour*)contours->GetObjectByOffset(i); + IsInPoly(in_nod,contour->x,contour->y,contour->nods,x,y,i0,i1,edgevalue); } Index: /issm/trunk-jpl/src/m/utils/BC/SetIceShelfBC.py =================================================================== --- /issm/trunk-jpl/src/m/utils/BC/SetIceShelfBC.py (revision 12127) +++ /issm/trunk-jpl/src/m/utils/BC/SetIceShelfBC.py (revision 12127) @@ -0,0 +1,93 @@ +def SetIceShelfBC(md,icefrontfile='') + #SETICESHELFBC - Create the boundary conditions for diagnostic and thermal models for a Ice Shelf with Ice Front + # + # Neumann BC are used on the ice front (an ANRGUS contour around the ice front + # must be given in input) + # Dirichlet BC are used elsewhere for diagnostic + # + # Usage: + # md=SetIceShelfBC(md,varargin) + # + # Example: + # md=SetIceShelfBC(md); + # md=SetIceShelfBC(md,'Front.exp'); + # + # See also: SETICESHEETBC, SETMARINEICESHEETBC + + #node on Dirichlet (boundary and ~icefront) + if not icefrontfile: + nodeonicefront=zeros(md.mesh.numberofvertices) + else: + if not os.path.isfile(icefrontfile): + print 'SetIceShelfBC error message: ice front file '+icefrontfile+ ' not found' + return [] + nodeinsideicefront=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,icefrontfile,'node',2); + nodeonicefront=double(md.mesh.vertexonboundary & nodeinsideicefront); + + pos=find(md.mesh.vertexonboundary & ~nodeonicefront); + md.diagnostic.spcvx=NaN*ones(md.mesh.numberofvertices,1); + md.diagnostic.spcvy=NaN*ones(md.mesh.numberofvertices,1); + md.diagnostic.spcvz=NaN*ones(md.mesh.numberofvertices,1); + md.diagnostic.spcvx(pos)=0; + md.diagnostic.spcvy(pos)=0; + md.diagnostic.spcvz(pos)=0; + md.diagnostic.referential=NaN*ones(md.mesh.numberofvertices,6); + + #Dirichlet Values + if (length(md.inversion.vx_obs)==md.mesh.numberofvertices & length(md.inversion.vy_obs)==md.mesh.numberofvertices) + disp(' boundary conditions for diagnostic model: spc set as observed velocities'); + md.diagnostic.spcvx(pos)=md.inversion.vx_obs(pos); + md.diagnostic.spcvy(pos)=md.inversion.vy_obs(pos); + else + disp(' boundary conditions for diagnostic model: spc set as zero'); + end + + #segment on Ice Front + #segment on Neumann (Ice Front) + pos=find(nodeonicefront(md.mesh.segments(:,1)) | nodeonicefront(md.mesh.segments(:,2))); + if (md.mesh.dimension==2) + pressureload=md.mesh.segments(pos,:); + elseif md.mesh.dimension==3 + pressureload_layer1=[md.mesh.segments(pos,1:2) md.mesh.segments(pos,2)+md.mesh.numberofvertices2d md.mesh.segments(pos,1)+md.mesh.numberofvertices2d md.mesh.segments(pos,3)]; + pressureload=[]; + for i=1:md.mesh.numberoflayers-1, + pressureload=[pressureload ;pressureload_layer1(:,1:4)+(i-1)*md.mesh.numberofvertices2d pressureload_layer1(:,5)+(i-1)*md.mesh.numberofelements2d ]; + end + end + + #Add water or air enum depending on the element + pressureload=[pressureload 1*md.mask.elementonfloatingice(pressureload(:,end))]; + + #plug onto model + md.diagnostic.icefront=pressureload; + + #Create zeros basalforcings and surfaceforcings + if isnan(md.surfaceforcings.precipitation), + md.surfaceforcings.precipitation=zeros(md.mesh.numberofvertices,1); + disp(' no surfaceforcings.precipitation specified: values set as zero'); + end + if isnan(md.surfaceforcings.mass_balance), + md.surfaceforcings.mass_balance=zeros(md.mesh.numberofvertices,1); + disp(' no surfaceforcings.mass_balance specified: values set as zero'); + end + if isnan(md.basalforcings.melting_rate), + md.basalforcings.melting_rate=zeros(md.mesh.numberofvertices,1); + disp(' no basalforcings.melting_rate specified: values set as zero'); + end + if isnan(md.balancethickness.thickening_rate), + md.balancethickness.thickening_rate=zeros(md.mesh.numberofvertices,1); + disp(' no balancethickness.thickening_rate specified: values set as zero'); + end + + md.prognostic.spcthickness=NaN*ones(md.mesh.numberofvertices,1); + md.balancethickness.spcthickness=NaN*ones(md.mesh.numberofvertices,1); + + if (length(md.initialization.temperature)==md.mesh.numberofvertices), + md.thermal.spctemperature=NaN*ones(md.mesh.numberofvertices,1); + pos=find(md.mesh.vertexonsurface); md.thermal.spctemperature(pos)=md.initialization.temperature(pos); #impose observed temperature on surface + if (length(md.basalforcings.geothermalflux)~=md.mesh.numberofvertices), + md.basalforcings.geothermalflux=zeros(md.mesh.numberofvertices,1); + end + else + disp(' no thermal boundary conditions created: no observed temperature found'); + end Index: /issm/trunk-jpl/src/modules/ContourToMesh/ContourToMesh.cpp =================================================================== --- /issm/trunk-jpl/src/modules/ContourToMesh/ContourToMesh.cpp (revision 12126) +++ /issm/trunk-jpl/src/modules/ContourToMesh/ContourToMesh.cpp (revision 12127) @@ -21,5 +21,5 @@ #include "./ContourToMesh.h" -void mexFunction( int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) { +WRAPPER(ContourToMesh){ int i,j; @@ -28,69 +28,46 @@ int edgevalue; double *index = NULL; + int nel; double *x = NULL; + int nods; double *y = NULL; char *interptype = NULL; + char* contourname = NULL; + DataSet* contours = NULL; /* output: */ Vector *in_nod = NULL; - int nods; Vector *in_elem = NULL; - int nel; - - //contours - mxArray *matlabstructure = NULL; - int numcontours; - Contour **contours = NULL; - Contour *contouri = NULL; /*Boot module: */ MODULEBOOT(); - /*checks on arguments on the matlab side: */ - //CheckNumMatlabArguments(nlhs,NLHS,nrhs,NRHS,__FUNCT__,&ContourToMeshUsage); Cant' use it here, as we have variable outputs. - if((nlhs!=1 && nlhs!=2) || (nrhs!=NRHS)){ + /*checks on output arguments on the matlab side: */ + #ifdef _HAVE_MATLAB_MODULES_ + if(nlhs!=1 && nlhs!=2){ + ContourToMeshUsage(); + _error_(" usage. See above"); + } + #endif + /*check on input arguments: */ + if(nrhs!=NRHS){ ContourToMeshUsage(); _error_(" usage. See above"); } - /*First, call expread on filename to build a contour array in the matlab workspace: */ - mexCallMATLAB( 1, &matlabstructure, 1, (mxArray**)&FILENAME, "expread"); /*Fetch inputs: */ - FetchData(&index,&nel,NULL,INDEXHANDLE); - FetchData(&x,&nods,NULL,XHANDLE); - FetchData(&y,NULL,NULL,YHANDLE); - FetchData(&edgevalue,EDGEVALUEHANDLE); - - //Fetch contours - numcontours=mxGetNumberOfElements(matlabstructure); - contours=(Contour**)xmalloc(numcontours*sizeof(Contour*)); - for(i=0;inods=(int)mxGetScalar(mxGetField(matlabstructure,i,"nods")); - //set pointers. - contouri->x=mxGetPr(mxGetField(matlabstructure,i,"x")); - contouri->y=mxGetPr(mxGetField(matlabstructure,i,"y")); - *(contours+i)=contouri; - } + FetchData(&index,&nel,NULL,INDEX); + FetchData(&x,&nods,NULL,X); + FetchData(&y,NULL,NULL,Y); + FetchData(&edgevalue,EDGEVALUE); + FetchData(&contourname,CONTOURNAME); + contours=DomainOutlineRead(contourname); /*Fetch interptype: */ - FetchData(&interptype,INTERPTYPEHANDLE); - - /* Debugging of contours :{{{1*/ - /*for(i=0;inods); - for (j=0;jnods;j++){ - printf(" %lf %lf\n",*(contouri->x+j),*(contouri->y+j)); - } - }*/ - /*}}}*/ + FetchData(&interptype,INTERPTYPE); /*Run interpolation routine: */ - ContourToMeshx( &in_nod,&in_elem,index,x,y,contours,numcontours,interptype,nel,nods,edgevalue); + ContourToMeshx( &in_nod,&in_elem,index,x,y,contours,interptype,nel,nods,edgevalue); /* output: */ @@ -112,5 +89,5 @@ } -void ContourToMeshUsage(void) +void ContourToMeshUsage(void)//{{{1 { printf("CONTOURTOMESH - Flag the elements or nodes inside a contour\n"); @@ -134,2 +111,3 @@ printf("\n"); } +//}}} Index: /issm/trunk-jpl/src/modules/ContourToMesh/ContourToMesh.h =================================================================== --- /issm/trunk-jpl/src/modules/ContourToMesh/ContourToMesh.h (revision 12126) +++ /issm/trunk-jpl/src/modules/ContourToMesh/ContourToMesh.h (revision 12127) @@ -8,32 +8,55 @@ #define _CONTOURTOMESH_H -/* local prototypes: */ -void ContourToMeshUsage(void); +#ifdef HAVE_CONFIG_H + #include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +/*Very important definition in case we are compiling a python module!: needs to come before header files inclusion*/ +#ifdef _HAVE_PYTHON_ +#define PY_ARRAY_UNIQUE_SYMBOL PythonIOSymbol +#endif + +/*Header files: */ #include "../../c/include/globals.h" +#include "../../c/toolkits/toolkits.h" +#include "../../c/include/include.h" #include "../../c/modules/modules.h" #include "../../c/Container/Container.h" #include "../../c/shared/shared.h" #include "../../c/issm-binding.h" +#include "../../c/io/io.h" +#include "../../c/EnumDefinitions/EnumDefinitions.h" + +#ifdef _HAVE_MATLAB_MODULES_ +/* serial input macros: */ +#define INDEX (mxArray *)prhs[0] +#define X (mxArray *)prhs[1] +#define Y (mxArray *)prhs[2] +#define CONTOURNAME (mxArray *)prhs[3] +#define INTERPTYPE (mxArray *)prhs[4] +#define EDGEVALUE (mxArray *)prhs[5] +/* serial output macros: */ +#define PLHS0 (mxArray**)&plhs[0] +#define PLHS1 (mxArray**)&plhs[1] +#endif + +#ifdef _HAVE_PYTHON_MODULES_ +/* serial input macros: */ +#define INDEX PyTuple_GetItem(args,0) +#define X PyTuple_GetItem(args,1) +#define Y PyTuple_GetItem(args,2) +#define CONTOURNAME PyTuple_GetItem(args,3) +#define INTERPTYPE PyTuple_GetItem(args,4) +#define EDGEVALUE PyTuple_GetItem(args,5) +/* serial output macros: */ +#define PLHS0 output,0 +#define PLHS1 output,1 +#endif #undef __FUNCT__ #define __FUNCT__ "ContourToMesh" - - -#ifndef ALL -#define ALL 0 -#endif - -/* input macros: */ -#define INDEXHANDLE prhs[0] -#define XHANDLE prhs[1] -#define YHANDLE prhs[2] -#define FILENAME prhs[3] -#define INTERPTYPEHANDLE prhs[4] -#define EDGEVALUEHANDLE prhs[5] - -/* serial output macros: */ -#define PLHS0 (mxArray**)&plhs[0] -#define PLHS1 (mxArray**)&plhs[1] /* serial arg counts: */ @@ -43,4 +66,6 @@ #define NRHS 6 +/* local prototypes: */ +void ContourToMeshUsage(void); #endif /* _CONTOURTOMESH_H */