Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 13633) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 13634) @@ -360,6 +360,7 @@ ./classes/objects/Options/GenericOption.h\ ./classes/objects/Options/OptionUtilities.cpp\ - ./classes/objects/Options/OptionUtilities.h - + ./classes/objects/Options/OptionUtilities.h\ + ./classes/RiftStruct.cpp\ + ./classes/RiftStruct.h #}}} #DAKOTA sources {{{ @@ -848,4 +849,6 @@ ./modules/TriMeshx/TriMeshx.h\ ./modules/TriMeshx/TriMeshx.cpp\ + ./modules/TriMeshProcessRiftsx/TriMeshProcessRiftsx.h\ + ./modules/TriMeshProcessRiftsx/TriMeshProcessRiftsx.cpp\ ./modules/PointCloudFindNeighborsx/PointCloudFindNeighborsx.h\ ./modules/PointCloudFindNeighborsx/PointCloudFindNeighborsx.cpp\ Index: /issm/trunk-jpl/src/c/classes/DofIndexing.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/DofIndexing.cpp (revision 13633) +++ /issm/trunk-jpl/src/c/classes/DofIndexing.cpp (revision 13634) @@ -10,6 +10,6 @@ #include +#include #include "./classes.h" -#include #include "../EnumDefinitions/EnumDefinitions.h" #include "../shared/shared.h" Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 13633) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 13634) @@ -20,4 +20,5 @@ class Materials; class Profiler; +class Elements; /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/RiftStruct.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/RiftStruct.cpp (revision 13634) +++ /issm/trunk-jpl/src/c/classes/RiftStruct.cpp (revision 13634) @@ -0,0 +1,95 @@ +/*!\file RiftStruct.c + * \brief: implementation of the RiftStruct object + */ + +#include "./classes.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../shared/shared.h" +#include "../include/include.h" + +/*RiftStruct constructors and destructor*/ +RiftStruct::RiftStruct(void){/*{{{*/ + + this->numrifts = 0; + this->riftsnumsegments = NULL; + this->riftssegments = NULL; + this->riftsnumpairs = NULL; + this->riftspairs = NULL; + this->riftsnumpenaltypairs = NULL; + this->riftspenaltypairs = NULL; + this->riftstips = NULL; + this->state = NULL; + +}/*}}}*/ +RiftStruct::RiftStruct(int numrifts_in,int *riftsnumsegments_in,double **riftssegments_in,int *riftsnumpairs_in,double **riftspairs_in,int *riftsnumpenaltypairs_in,double **riftspenaltypairs_in,double* riftstips_in){/*{{{*/ + + int i; + + /*numrifts*/ + this->numrifts = numrifts_in; + if(!numrifts_in) return; + + /*riftsnumsegments*/ + _assert_(riftsnumsegments_in); + this->riftsnumsegments=xNew(numrifts_in); + xMemCpy(this->riftsnumsegments,riftsnumsegments_in,numrifts_in); + + /*riftssegments*/ + _assert_(riftssegments_in); + this->riftssegments=xNew(numrifts_in); + for(i=0;iriftssegments[i]=xNew(riftsnumsegments_in[i]*3); + xMemCpy(this->riftssegments[i],riftssegments_in[i],riftsnumsegments_in[i]*3); + } + + /*riftsnumpairs*/ + _assert_(riftsnumpairs_in); + this->riftsnumpairs=xNew(numrifts_in); + xMemCpy(this->riftsnumpairs,riftsnumpairs_in,numrifts_in); + + /*riftspairs*/ + _assert_(riftspairs_in); + this->riftspairs=xNew(numrifts_in); + for(i=0;iriftspairs[i]=xNew(riftsnumpairs_in[i]*2); + xMemCpy(this->riftspairs[i],riftspairs_in[i],riftsnumpairs_in[i]*2); + } + + /*riftsnumpenaltypairs*/ + _assert_(riftsnumpenaltypairs_in); + this->riftsnumpenaltypairs=xNew(numrifts_in); + xMemCpy(this->riftsnumpenaltypairs,riftsnumpenaltypairs_in,numrifts_in); + + /*riftspenaltypairs*/ + _assert_(riftspenaltypairs_in); + this->riftspenaltypairs=xNew(numrifts_in); + for(i=0;iriftspenaltypairs[i]=xNew(riftsnumpenaltypairs_in[i]*7); + xMemCpy(this->riftspenaltypairs[i],riftspenaltypairs_in[i],riftsnumpenaltypairs_in[i]*7); + } + + /*riftstips*/ + _assert_(riftstips_in); + this->riftstips=xNew(2*numrifts_in); + xMemCpy(this->riftstips,riftstips_in,2*numrifts_in); + + /*state*/ + this->state=xNew(numrifts_in); + for(i=0;istate[i]=xNew(riftsnumpenaltypairs_in[i]); + for(int j=0;jstate[i])[j]=FreeEnum; + } + +}/*}}}*/ +RiftStruct::~RiftStruct(void){/*{{{*/ + + xDelete(this->riftsnumsegments); + xDelete(this->riftssegments); + xDelete(this->riftsnumpairs); + xDelete(this->riftspairs); + xDelete(this->riftsnumpenaltypairs); + xDelete(this->riftspenaltypairs); + xDelete(this->riftstips); + xDelete(this->state); + +}/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/RiftStruct.h =================================================================== --- /issm/trunk-jpl/src/c/classes/RiftStruct.h (revision 13634) +++ /issm/trunk-jpl/src/c/classes/RiftStruct.h (revision 13634) @@ -0,0 +1,26 @@ +/*!\file: RiftStruct.h + * \brief place holder for Rift Structure in mex module + */ + +#ifndef _RIFTSTRUCT_H_ +#define _RIFTSTRUCT_H_ + +class RiftStruct{ + + public: + int numrifts; + double **riftssegments; + int *riftsnumsegments; + double **riftspairs; + int *riftsnumpairs; + double **riftspenaltypairs; + int *riftsnumpenaltypairs; + double *riftstips; + double **state; + + RiftStruct(); + RiftStruct(int numrifts_in,int *riftsnumsegments_in,double **riftssegments_in,int *riftsnumpairs_in,double **riftspairs_in,int *riftsnumpenaltypairs_in,double **riftspenaltypairs_in,double* riftstips_in); + ~RiftStruct(); +}; + +#endif Index: /issm/trunk-jpl/src/c/classes/Update.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Update.h (revision 13633) +++ /issm/trunk-jpl/src/c/classes/Update.h (revision 13634) @@ -7,7 +7,6 @@ /*Headers:*/ -/*{{{*/ #include "../shared/shared.h" -/*}}}*/ +class IoModel; class Update{ Index: /issm/trunk-jpl/src/c/classes/classes.h =================================================================== --- /issm/trunk-jpl/src/c/classes/classes.h (revision 13633) +++ /issm/trunk-jpl/src/c/classes/classes.h (revision 13634) @@ -32,4 +32,5 @@ #include "./AdolcEdf.h" #include "./IssmComm.h" +#include "./RiftStruct.h" #endif Index: /issm/trunk-jpl/src/c/classes/objects/Vertex.h =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Vertex.h (revision 13633) +++ /issm/trunk-jpl/src/c/classes/objects/Vertex.h (revision 13634) @@ -9,11 +9,10 @@ /*{{{*/ #include "../classes.h" -class IoModel; -template class Vector; -class Parameters; #include "../../shared/Exceptions/exceptions.h" #include "../../toolkits/toolkits.h" #include "../../include/include.h" - +template class Vector; +class Parameters; +class IoModel; /*}}}*/ Index: /issm/trunk-jpl/src/c/modules/TriMeshProcessRiftsx/TriMeshProcessRiftsx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/TriMeshProcessRiftsx/TriMeshProcessRiftsx.cpp (revision 13634) +++ /issm/trunk-jpl/src/c/modules/TriMeshProcessRiftsx/TriMeshProcessRiftsx.cpp (revision 13634) @@ -0,0 +1,79 @@ +/*!\file: TriMeshProcessRifts.cpp + * \brief split a mesh where a rift (or fault) is present + */ + +#include "./TriMeshProcessRiftsx.h" +#include "../../classes/RiftStruct.h" +#include "../../shared/shared.h" +#include "../../include/include.h" +#include "../../toolkits/toolkits.h" + +void TriMeshProcessRiftsx(double** pindex, int* pnel,double** px,double** py,int* pnods,double** psegments,double** psegmentmarkers,int *pnum_seg,RiftStruct **priftstruct){ + + /*Output*/ + int numrifts,numrifts0; + int *riftsnumsegments = NULL; + double **riftssegments = NULL; + int *riftsnumpairs = NULL; + double **riftspairs = NULL; + double *riftstips = NULL; + double **riftspenaltypairs = NULL; + int *riftsnumpenaltypairs = NULL; + + /*Recover initial mesh*/ + int nel = *pnel; + double *index = *pindex; + double *x = *px; + double *y = *py; + int nods = *pnods; + double *segments = *psegments; + double *segmentmarkers = *psegmentmarkers; + int num_seg = *pnum_seg; + + /*Intermediary*/ + int riftflag; + + /*First, do some fixing on the existing mesh: we do not want any element belonging entirely to the segment list (ie: + *all the nodes of this element belong to the segments (tends to happen when there are corners: */ + RemoveCornersFromRifts(&index,&nel,&x,&y,&nods,segments,segmentmarkers,num_seg); + + /*Figure out if we have rifts, and how many: */ + IsRiftPresent(&riftflag,&numrifts0,segmentmarkers,num_seg); + + if(!riftflag) _error_("No rift present in mesh"); + + /*Split mesh*/ + SplitMeshForRifts(&nel,&index,&nods,&x,&y,&num_seg,&segments,&segmentmarkers); + + /*Order segments so that their normals point outside the domain: */ + OrderSegments(&segments,num_seg, index,nel); + + /*We do not want to output segments mixed with rift segments: wring out the rifts from the segments, using the + *segmentmarkerlist:*/ + SplitRiftSegments(&segments,&segmentmarkers,&num_seg,&numrifts,&riftsnumsegments,&riftssegments,numrifts0,nods,nel); + + /*Using rift segments, associate rift faces in pairs, each pair face representing opposite flanks of the rifts facing one another directly: */ + PairRiftElements(&riftsnumpairs,&riftspairs,numrifts,riftsnumsegments,riftssegments,x,y); + + /*Order rifts so that they start from one tip, go to the other tip, and back: */ + OrderRifts(&riftstips,riftssegments,riftspairs,numrifts,riftsnumsegments,x,y,nods,nel); + + /*Create penalty pairs, used by Imp: */ + PenaltyPairs(&riftspenaltypairs,&riftsnumpenaltypairs,numrifts,riftssegments,riftsnumsegments,riftspairs,riftstips,x,y); + + /*Create Riftstruct*/ + RiftStruct* riftstruct = new RiftStruct(numrifts,riftsnumsegments,riftssegments,riftsnumpairs,riftspairs,riftsnumpenaltypairs,riftspenaltypairs,riftstips); + + /*Assign output pointers for mesh*/ + *pnel = nel; + *pindex = index; + *px = x; + *py = y; + *pnods = nods; + *psegments = segments; + *psegmentmarkers = segmentmarkers; + *pnum_seg = num_seg; + + /*Assign output pointers for rifts*/ + *priftstruct = riftstruct; +} Index: /issm/trunk-jpl/src/c/modules/TriMeshProcessRiftsx/TriMeshProcessRiftsx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/TriMeshProcessRiftsx/TriMeshProcessRiftsx.h (revision 13634) +++ /issm/trunk-jpl/src/c/modules/TriMeshProcessRiftsx/TriMeshProcessRiftsx.h (revision 13634) @@ -0,0 +1,12 @@ +/*!\file: TriMeshProcessRiftsx.h + * \brief header file for TriMeshProcessRifts module + */ + +#ifndef _TRIMESHPROCESSRIFTX_H +#define _TRIMESHPROCESSRIFTX_H + +class RiftStruct; + +void TriMeshProcessRiftsx(double** pindex, int* pnel,double** px,double** py,int* pnods,double** psegments,double** psegmentmarkers,int *pnum_seg,RiftStruct **priftstruct); + +#endif /* _TRIMESHPROCESSRIFTX_H*/ Index: /issm/trunk-jpl/src/c/modules/modules.h =================================================================== --- /issm/trunk-jpl/src/c/modules/modules.h (revision 13633) +++ /issm/trunk-jpl/src/c/modules/modules.h (revision 13634) @@ -117,4 +117,5 @@ #include "./TriaSearchx/TriaSearchx.h" #include "./TriMeshx/TriMeshx.h" +#include "./TriMeshProcessRiftsx/TriMeshProcessRiftsx.h" #include "./ThicknessAbsMisfitx/ThicknessAbsMisfitx.h" #include "./ThicknessAbsGradientx/ThicknessAbsGradientx.h" Index: /issm/trunk-jpl/src/c/shared/TriMesh/TriMeshUtils.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/TriMesh/TriMeshUtils.cpp (revision 13633) +++ /issm/trunk-jpl/src/c/shared/TriMesh/TriMeshUtils.cpp (revision 13634) @@ -266,5 +266,5 @@ /*Reallocate segments: */ segments =xReNew(segments, nsegs*3,(nsegs+nriftsegs)*3); - segmentmarkerlist=xReNew(segmentmarkerlist,nsegs*3,(nsegs+nriftsegs)*3); + segmentmarkerlist=xReNew(segmentmarkerlist,nsegs,(nsegs+nriftsegs)); /*First, update the existing segments to the new nodes :*/ Index: /issm/trunk-jpl/src/c/shared/shared.h =================================================================== --- /issm/trunk-jpl/src/c/shared/shared.h (revision 13633) +++ /issm/trunk-jpl/src/c/shared/shared.h (revision 13634) @@ -6,18 +6,18 @@ #define _SHARED_H_ -#include "Alloc/alloc.h" -#include "Alloc/xNewDelete.h" -#include "Bamg/shared.h" -#include "Elements/elements.h" -#include "Exceptions/exceptions.h" -#include "Exp/exp.h" -#include "Matrix/matrix.h" -#include "MemOps/xMemCpy.h" -#include "Numerics/numerics.h" -#include "Sorting/sorting.h" -#include "String/sharedstring.h" -#include "Threads/issm_threads.h" -#include "TriMesh/trimesh.h" -#include "Wrapper/wrappershared.h" +#include "./Alloc/alloc.h" +#include "./Alloc/xNewDelete.h" +#include "./Bamg/shared.h" +#include "./Elements/elements.h" +#include "./Exceptions/exceptions.h" +#include "./Exp/exp.h" +#include "./Matrix/matrix.h" +#include "./MemOps/xMemCpy.h" +#include "./Numerics/numerics.h" +#include "./Sorting/sorting.h" +#include "./String/sharedstring.h" +#include "./Threads/issm_threads.h" +#include "./TriMesh/trimesh.h" +#include "./Wrapper/wrappershared.h" #endif Index: /issm/trunk-jpl/src/modules/TriMeshProcessRifts/TriMeshProcessRifts.cpp =================================================================== --- /issm/trunk-jpl/src/modules/TriMeshProcessRifts/TriMeshProcessRifts.cpp (revision 13633) +++ /issm/trunk-jpl/src/modules/TriMeshProcessRifts/TriMeshProcessRifts.cpp (revision 13634) @@ -11,19 +11,8 @@ _printLine_(" index2,x2,y2,segments2,segmentmarkers2,rifts2 is the resulting triangulation where rifts have been processed."); }/*}}}*/ -void TriMeshProcessRiftsx(double** pindex, int* pnel,double** px,double** py,int* pnods,double** psegments,double** psegmentmarkers,int *pnum_seg, - int *pnumrifts, int **priftsnumsegments, double ***priftssegments, - int **priftsnumpairs, double ***priftspairs, double **priftstips, - double ***priftspenaltypairs, int **priftsnumpenaltypairs); WRAPPER(TriMeshProcessRifts){ /* returned quantities: */ - int numrifts; - int *riftsnumsegments = NULL; - double **riftssegments = NULL; - int *riftsnumpairs = NULL; - double **riftspairs = NULL; - double *riftstips = NULL; - double **riftspenaltypairs = NULL; - int *riftsnumpenaltypairs = NULL; + RiftStruct *riftstruct = NULL; /* input: */ @@ -51,5 +40,5 @@ /*call x layer*/ TriMeshProcessRiftsx(&index,&nel,&x,&y,&nods,&segments,&segmentmarkers,&num_seg, - &numrifts,&riftsnumsegments,&riftssegments,&riftsnumpairs,&riftspairs,&riftstips,&riftspenaltypairs,&riftsnumpenaltypairs); + &riftstruct); /*Output : */ @@ -60,12 +49,11 @@ WriteData(SEGMENTMARKERSOUT,segmentmarkers,num_seg,1); - mxArray *pmxa_array = NULL; - mxArray *pmxa_array2 = NULL; - mxArray *pmxa_array3 = NULL; + mxArray *pmxa_array = NULL; + mxArray *pmxa_array2 = NULL; + mxArray *pmxa_array3 = NULL; const char *fnames[10]; mwSize ndim = 2; mwSize dimensions[2] = {1,1}; double *pair = NULL; - /*Create a structure rifts where if i is a rift number, we have the following fields rifts(i).segments and rifts(i).numsegs: */ @@ -81,16 +69,16 @@ fnames[9] = "state"; - dimensions[0]=numrifts; - pmxa_array=mxCreateStructArray( ndim,dimensions,10,fnames); + dimensions[0]=riftstruct->numrifts; + pmxa_array=mxCreateStructArray(ndim,dimensions,10,fnames); - for(int i=0;inumrifts;i++){ /*Segments: */ - WriteData(&pmxa_array3,riftssegments[i],riftsnumsegments[i],3); + WriteData(&pmxa_array3,riftstruct->riftssegments[i],riftstruct->riftsnumsegments[i],3); mxSetField(pmxa_array,i,"segments",pmxa_array3); - mxSetField(pmxa_array,i,"numsegs",mxCreateDoubleScalar((double)riftsnumsegments[i])); + mxSetField(pmxa_array,i,"numsegs",mxCreateDoubleScalar((double)riftstruct->riftsnumsegments[i])); /*Element pairs: */ - WriteData(&pmxa_array3,riftspairs[i],riftsnumpairs[i],2); + WriteData(&pmxa_array3,riftstruct->riftspairs[i],riftstruct->riftsnumpairs[i],2); mxSetField(pmxa_array,i,"pairs",pmxa_array3); @@ -99,6 +87,6 @@ mxSetM(pmxa_array2,1); pair=(double*)mxMalloc(2*sizeof(double)); - pair[0]=*(riftstips+2*i+0); - pair[1]=*(riftstips+2*i+1); + pair[0]=*(riftstruct->riftstips+2*i+0); + pair[1]=*(riftstruct->riftstips+2*i+1); mxSetN(pmxa_array2,2); mxSetPr(pmxa_array2,pair); @@ -106,5 +94,5 @@ /*Penalty pairs: */ - WriteData(&pmxa_array3,riftspenaltypairs[i],riftsnumpenaltypairs[i],7); + WriteData(&pmxa_array3,riftstruct->riftspenaltypairs[i],riftstruct->riftsnumpenaltypairs[i],7); mxSetField(pmxa_array,i,"penaltypairs",pmxa_array3); @@ -116,10 +104,8 @@ /*State: */ - double *state=(double*)mxMalloc(riftsnumpenaltypairs[i]*sizeof(double)); - for(int j=0;jriftsnumpenaltypairs[i]); + mxSetPr(pmxa_array2,riftstruct->state[i]); mexCallMATLAB( 1, &pmxa_array3, 1, &pmxa_array2, "transpose"); @@ -132,78 +118,2 @@ MODULEEND(); } -void TriMeshProcessRiftsx(double** pindex, int* pnel,double** px,double** py,int* pnods,double** psegments,double** psegmentmarkers,int *pnum_seg, - int *pnumrifts, int **priftsnumsegments, double ***priftssegments, - int **priftsnumpairs, double ***priftspairs, double **priftstips, - double ***priftspenaltypairs, int **priftsnumpenaltypairs){ - - /*Output*/ - int numrifts,numrifts0; - int *riftsnumsegments = NULL; - double **riftssegments = NULL; - int *riftsnumpairs = NULL; - double **riftspairs = NULL; - double *riftstips = NULL; - double **riftspenaltypairs = NULL; - int *riftsnumpenaltypairs = NULL; - - /*Recover initial mesh*/ - int nel = *pnel; - double *index = *pindex; - double *x = *px; - double *y = *py; - int nods = *pnods; - double *segments = *psegments; - double *segmentmarkers = *psegmentmarkers; - int num_seg = *pnum_seg; - - /*Intermediary*/ - int riftflag; - - /*First, do some fixing on the existing mesh: we do not want any element belonging entirely to the segment list (ie: - *all the nodes of this element belong to the segments (tends to happen when there are corners: */ - RemoveCornersFromRifts(&index,&nel,&x,&y,&nods,segments,segmentmarkers,num_seg); - - /*Figure out if we have rifts, and how many: */ - IsRiftPresent(&riftflag,&numrifts0,segmentmarkers,num_seg); - - if(!riftflag) _error_("No rift present in mesh"); - - /*Split mesh*/ - SplitMeshForRifts(&nel,&index,&nods,&x,&y,&num_seg,&segments,&segmentmarkers); - - /*Order segments so that their normals point outside the domain: */ - OrderSegments(&segments,num_seg, index,nel); - - /*We do not want to output segments mixed with rift segments: wring out the rifts from the segments, using the - *segmentmarkerlist:*/ - SplitRiftSegments(&segments,&segmentmarkers,&num_seg,&numrifts,&riftsnumsegments,&riftssegments,numrifts0,nods,nel); - - /*Using rift segments, associate rift faces in pairs, each pair face representing opposite flanks of the rifts facing one another directly: */ - PairRiftElements(&riftsnumpairs,&riftspairs,numrifts,riftsnumsegments,riftssegments,x,y); - - /*Order rifts so that they start from one tip, go to the other tip, and back: */ - OrderRifts(&riftstips,riftssegments,riftspairs,numrifts,riftsnumsegments,x,y,nods,nel); - - /*Create penalty pairs, used by Imp: */ - PenaltyPairs(&riftspenaltypairs,&riftsnumpenaltypairs,numrifts,riftssegments,riftsnumsegments,riftspairs,riftstips,x,y); - - /*Assign output pointers for mesh*/ - *pnel = nel; - *pindex = index; - *px = x; - *py = y; - *pnods = nods; - *psegments = segments; - *psegmentmarkers = segmentmarkers; - *pnum_seg = num_seg; - - /*Assign output pointers for rifts*/ - *pnumrifts = numrifts; - *priftsnumsegments = riftsnumsegments; - *priftssegments = riftssegments; - *priftsnumpairs = riftsnumpairs; - *priftspairs = riftspairs; - *priftstips = riftstips; - *priftspenaltypairs = riftspenaltypairs; - *priftsnumpenaltypairs = riftsnumpenaltypairs; -}