source: issm/trunk/src/modules/TriMeshProcessRifts/TriMeshProcessRifts.cpp@ 12706

Last change on this file since 12706 was 12706, checked in by Mathieu Morlighem, 13 years ago

merged trunk-jpl and trunk for revision 12703
File size: 9.2 KB
RevLine 
[1]1/*!\file: TriMeshProcessRifts.cpp
2 * \brief split a mesh where a rift (or fault) is present
3 */
4
5#include "./TriMeshProcessRifts.h"
6
7void mexFunction( int nlhs, mxArray* plhs[],
8 int nrhs, const mxArray* prhs[] ) {
9
10
11 /*Matlab arrays: */
12 mxArray* pmxa_array=NULL;
13 mxArray* pmxa_array2=NULL;
14 mxArray* pmxa_array3=NULL;
15 int i,j,k,counter;
16
17 /* returned quantities: */
18 int out_numrifts;
19 int* out_riftsnumsegments=NULL;
20 double** out_riftssegments=NULL;
21 int* out_riftsnumpairs=NULL;
22 double** out_riftspairs=NULL;
23 double* out_riftstips=NULL;
24 double** out_riftspenaltypairs=NULL;
25 int* out_riftsnumpenaltypairs=NULL;
26
27 /*empty rifts structure: */
28 double* pNaN=NULL;
[6748]29 const char* fnames[10];
[1]30 mwSize ndim=2;
31 mwSize dimensions[2] = {1,1};
32 double* pair=NULL;
33
34
35 /* input: */
36 double* tindex_in=NULL;
37 double* index_in=NULL;
38 int nel;
39 double* x_inm=NULL; //matlab vector
40 double* x_in=NULL; //copy of matlab vector
41 int nods;
42 double* y_inm=NULL;//matlab vector
43 double* y_in=NULL;//copy of matlab vector
44 double* tsegments_in=NULL;
45 double* segments_in=NULL;
46 double* tsegmentmarkers_in=NULL;
47 double* segmentmarkers_in=NULL;
48
[6748]49 /* state: */
50 double* state=NULL;
51
[1]52 int num_seg;
53
54 /*rifts: */
55 int riftflag;
56 int numrifts;
57
58 /* verify correct usage: */
59 if (nlhs==0 && nrhs==0) {
60 /* special case: */
61 TriMeshProcessRiftsUsage();
62 return;
63 }
64
65 if (!( (nlhs==6) || (nrhs==5))){
66 mexPrintf(" %s format error.\n", __FUNCT__);
67 TriMeshProcessRiftsUsage();
[12706]68 _error_("bad usage");
[1]69 }
70
71 /*Fetch index_in: */
72 if(mxIsDouble(prhs[0])){
73 nel=mxGetM(prhs[0]);
74 tindex_in=mxGetPr(prhs[0]);
[12060]75 index_in=(double*)xmalloc(nel*3*sizeof(double));
[1]76 for (i=0;i<nel;i++){
77 for (j=0;j<3;j++){
78 *(index_in+3*i+j)=*(tindex_in+nel*j+i);
79 }
80 }
81 }
82 else{
[12706]83 _error_("first argument should be the element list");
[1]84 }
85
86 /*Fetch x_in: */
87 if(mxIsDouble(prhs[1])){
88 nods=mxGetM(prhs[1]);
89 x_inm=mxGetPr(prhs[1]);
[12060]90 x_in=(double*)xmalloc(nods*sizeof(double));
[1]91 for (i=0;i<nods;i++){
92 x_in[i]=x_inm[i];
93 }
94 }
95 else{
[12706]96 _error_("second argument should be the x corrdinate list");
[1]97 }
98
99 /*Fetch y_in: */
100 if(mxIsDouble(prhs[2])){
101 y_inm=mxGetPr(prhs[2]);
[12060]102 y_in=(double*)xmalloc(nods*sizeof(double));
[1]103 for (i=0;i<nods;i++){
104 y_in[i]=y_inm[i];
105 }
106 }
107 else{
[12706]108 _error_("third argument should be the y corrdinate list");
[1]109 }
110
111 /*Fetch segments_in: */
112 if(mxIsDouble(prhs[3])){
113 num_seg=mxGetM(prhs[3]);
114 tsegments_in=mxGetPr(prhs[3]);
[12060]115 segments_in=(double*)xmalloc(num_seg*3*sizeof(double));
[1]116 for (i=0;i<num_seg;i++){
117 for (j=0;j<3;j++){
118 *(segments_in+3*i+j)=*(tsegments_in+num_seg*j+i);
119 }
120 }
121 }
122 else{
[12706]123 _error_("fourth argument should be the segments list");
[1]124 }
125
126 /*Fetch segment markers: */
127 if(mxIsDouble(prhs[4])){
128 tsegmentmarkers_in=mxGetPr(prhs[4]);
[12060]129 segmentmarkers_in=(double*)xmalloc(num_seg*sizeof(double));
[1]130 for (i=0;i<num_seg;i++){
131 segmentmarkers_in[i]=tsegmentmarkers_in[i];
132 }
133 }
134 else{
[12706]135 _error_("fourth argument should be the segmentmarkers list");
[1]136 }
137
138 /*First, do some fixing on the existing mesh: we do not want any element belonging entirely to the segment list (ie:
[8306]139 *all the nodes of this element belong to the segments (tends to happen when there are corners: */
[1]140 RemoveCornersFromRifts(&index_in,&nel,&x_in,&y_in,&nods,segments_in,segmentmarkers_in,num_seg);
141
142 /*Figure out if we have rifts, and how many: */
143 IsRiftPresent(&riftflag,&numrifts,segmentmarkers_in,num_seg);
144
145 if(riftflag){
146 SplitMeshForRifts(&nel,&index_in,&nods,&x_in,&y_in,&num_seg,&segments_in,&segmentmarkers_in);
147 }
148
149 /*Order segments so that their normals point outside the domain: */
[4453]150 OrderSegments(&segments_in,num_seg, index_in,nel);
[1]151
152 if(riftflag){
153
154 /*We do not want to output segments mixed with rift segments: wring out the rifts from the segments, using the
155 *segmentmarkerlist:*/
[12081]156 SplitRiftSegments(&segments_in,&segmentmarkers_in,&num_seg,&out_numrifts,&out_riftsnumsegments,&out_riftssegments,numrifts,nods,nel);
[1]157
158 /*Using rift segments, associate rift faces in pairs, each pair face representing opposite flanks of the rifts facing one another directly: */
159 PairRiftElements(&out_riftsnumpairs,&out_riftspairs,out_numrifts,out_riftsnumsegments,out_riftssegments,x_in,y_in);
160
161 /*Order rifts so that they start from one tip, go to the other tip, and back: */
[12071]162 OrderRifts(&out_riftstips, out_riftssegments,out_riftspairs,numrifts,out_riftsnumsegments,x_in,y_in,nods,nel);
[1]163
164 /*Create penalty pairs, used by Imp: */
165 PenaltyPairs(&out_riftspenaltypairs,&out_riftsnumpenaltypairs,numrifts,out_riftssegments,out_riftsnumsegments,out_riftspairs,out_riftstips,x_in,y_in);
166 }
167
168
169 /*Output : */
[12085]170 WriteData(&plhs[0],index_in,nel,3);
171 //pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
172 //mxSetM(pmxa_array,3);
173 //mxSetN(pmxa_array,nel);
174 //mxSetPr(pmxa_array,index_in);
175 //mexCallMATLAB( 1, &plhs[0], 1, &pmxa_array, "transpose");
[1]176
[12085]177 WriteData(&plhs[1],x_in,nods,1);
178 //pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
179 //mxSetM(pmxa_array,1);
180 //mxSetN(pmxa_array,nods);
181 //mxSetPr(pmxa_array,x_in);
182 //mexCallMATLAB( 1, &plhs[1], 1, &pmxa_array, "transpose");
[1]183
[12085]184 WriteData(&plhs[2],y_in,nods,1);
185 //pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
186 //mxSetM(pmxa_array,1);
187 //mxSetN(pmxa_array,nods);
188 //mxSetPr(pmxa_array,y_in);
189 //mexCallMATLAB( 1, &plhs[2], 1, &pmxa_array, "transpose");
[1]190
[12085]191 WriteData(&plhs[3],segments_in,num_seg,3);
192//pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
193//mxSetM(pmxa_array,3);
194//mxSetN(pmxa_array,num_seg);
195//mxSetPr(pmxa_array,segments_in);
196//mexCallMATLAB( 1, &plhs[3], 1, &pmxa_array, "transpose");
[1]197
[12085]198 WriteData(&plhs[4],segmentmarkers_in,num_seg,1);
199 //pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
200 //mxSetM(pmxa_array,1);
201 //mxSetN(pmxa_array,num_seg);
202 //mxSetPr(pmxa_array,segmentmarkers_in);
203 //mexCallMATLAB( 1, &plhs[4], 1, &pmxa_array, "transpose");
[1]204
205 if(riftflag){
206 /*Create a structure rifts where if i is a rift number, we have the following fields rifts(i).segments and rifts(i).numsegs: */
207
208 fnames[0] = "numsegs";
209 fnames[1] = "segments";
210 fnames[2] = "pairs";
211 fnames[3] = "tips";
212 fnames[4] = "penaltypairs";
213 fnames[5] = "fill";
214 fnames[6] = "friction";
[1784]215 fnames[7] = "fraction";
[1805]216 fnames[8] = "fractionincrement";
[6748]217 fnames[9] = "state";
[1]218
219 dimensions[0]=out_numrifts;
220
[6748]221 pmxa_array=mxCreateStructArray( ndim,dimensions,10,fnames);
[1]222
223 for (i=0;i<out_numrifts;i++){
[12085]224
[1]225 /*Segments: */
[12085]226 WriteData(&pmxa_array3,out_riftssegments[i],out_riftsnumsegments[i],3);
227 //pmxa_array2= mxCreateDoubleMatrix(0,0,mxREAL);
228 //mxSetM(pmxa_array2,3);
229 //mxSetN(pmxa_array2,out_riftsnumsegments[i]);
230 //mxSetPr(pmxa_array2,out_riftssegments[i]);
231 //mexCallMATLAB( 1, &pmxa_array3, 1, &pmxa_array2, "transpose");
[1]232 mxSetField(pmxa_array,i,"segments",pmxa_array3);
[8963]233 mxSetField(pmxa_array,i,"numsegs",mxCreateDoubleScalar((double)out_riftsnumsegments[i]));
[1]234
235 /*Element pairs: */
[12085]236 WriteData(&pmxa_array3,out_riftspairs[i],out_riftsnumpairs[i],2);
237 //pmxa_array2= mxCreateDoubleMatrix(0,0,mxREAL);
238 //mxSetM(pmxa_array2,2);
239 //mxSetN(pmxa_array2,out_riftsnumpairs[i]);
240 //mxSetPr(pmxa_array2,out_riftspairs[i]);
241 //mexCallMATLAB( 1, &pmxa_array3, 1, &pmxa_array2, "transpose");
[1]242 mxSetField(pmxa_array,i,"pairs",pmxa_array3);
243
244 /*Tips: */
245 pmxa_array2= mxCreateDoubleMatrix(0,0,mxREAL);
246 mxSetM(pmxa_array2,1);
[12085]247 pair=(double*)mxMalloc(2*sizeof(double));
[1]248 pair[0]=*(out_riftstips+2*i+0);
249 pair[1]=*(out_riftstips+2*i+1);
250 mxSetN(pmxa_array2,2);
251 mxSetPr(pmxa_array2,pair);
252 mxSetField(pmxa_array,i,"tips",pmxa_array2);
253
254 /*Penalty pairs: */
[12085]255 WriteData(&pmxa_array3,out_riftspenaltypairs[i],out_riftsnumpenaltypairs[i],7);
256 //pmxa_array2= mxCreateDoubleMatrix(0,0,mxREAL);
257 //mxSetM(pmxa_array2,7);
258 //mxSetN(pmxa_array2,out_riftsnumpenaltypairs[i]);
259 //mxSetPr(pmxa_array2,out_riftspenaltypairs[i]);
260 //mexCallMATLAB( 1, &pmxa_array3, 1, &pmxa_array2, "transpose");
[1]261 mxSetField(pmxa_array,i,"penaltypairs",pmxa_array3);
262
[1805]263 /*Friction fraction, fractionincrement and fill: */
[1]264 mxSetField(pmxa_array,i,"friction",mxCreateDoubleScalar(0));
[3567]265 mxSetField(pmxa_array,i,"fill",mxCreateDoubleScalar(IceEnum)); //default is ice
[1784]266 mxSetField(pmxa_array,i,"fraction",mxCreateDoubleScalar(0)); //default is ice
[1805]267 mxSetField(pmxa_array,i,"fractionincrement",mxCreateDoubleScalar(0.1));
[6748]268
269 /*State: */
[12085]270 state=(double*)mxMalloc(out_riftsnumpenaltypairs[i]*sizeof(double));
[6748]271 for(j=0;j<out_riftsnumpenaltypairs[i];j++)state[j]=FreeEnum;
272 pmxa_array2= mxCreateDoubleMatrix(0,0,mxREAL);
273 mxSetM(pmxa_array2,1);
274 mxSetN(pmxa_array2,out_riftsnumpenaltypairs[i]);
275 mxSetPr(pmxa_array2,state);
276 mexCallMATLAB( 1, &pmxa_array3, 1, &pmxa_array2, "transpose");
277
278 mxSetField(pmxa_array,i,"state",pmxa_array3);
[1]279 }
280 }
281 else{
282 /*output NaN :*/
[12085]283 pNaN=(double*)mxMalloc(sizeof(double));
[3567]284 *pNaN=NAN;
[1]285 pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
286 mxSetM(pmxa_array,1);
287 mxSetN(pmxa_array,1);
288 mxSetPr(pmxa_array,pNaN);
289
290 }
291 plhs[5]=pmxa_array;
292
293 return;
294}
295
296
297void TriMeshProcessRiftsUsage(void)
298{
[12706]299 _printLine_("");
300 _printLine_(" usage: [index2,x2,y2,segments2,segmentmarkers2,rifts2]=TriMeshProcessrifts(index1,x1,y1,segments1,segmentmarkers1) ");
301 _printLine_(" where: (index1,x1,y1,segments1,segmentmarkers1) is an initial triangulation.");
302 _printLine_(" index2,x2,y2,segments2,segmentmarkers2,rifts2 is the resulting triangulation where rifts have been processed.");
[1]303}
Note: See TracBrowser for help on using the repository browser.