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

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

merged trunk-jpl and trunk for revision 12326M
File size: 9.9 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();
68 printf(" ");
69 mexErrMsgTxt(" ");
70 }
71
72 /*Fetch index_in: */
73 if(mxIsDouble(prhs[0])){
74 nel=mxGetM(prhs[0]);
75 tindex_in=mxGetPr(prhs[0]);
[12060]76 index_in=(double*)xmalloc(nel*3*sizeof(double));
[1]77 for (i=0;i<nel;i++){
78 for (j=0;j<3;j++){
79 *(index_in+3*i+j)=*(tindex_in+nel*j+i);
80 }
81 }
82 }
83 else{
84 printf("%s%s\n",__FUNCT__," error message: first argument should be the element list!");
85 mexErrMsgTxt(" ");
86 }
87
88 /*Fetch x_in: */
89 if(mxIsDouble(prhs[1])){
90 nods=mxGetM(prhs[1]);
91 x_inm=mxGetPr(prhs[1]);
[12060]92 x_in=(double*)xmalloc(nods*sizeof(double));
[1]93 for (i=0;i<nods;i++){
94 x_in[i]=x_inm[i];
95 }
96 }
97 else{
98 printf("%s%s\n",__FUNCT__," error message: second argument should be the x corrdinate list!");
99 mexErrMsgTxt(" ");
100 }
101
102 /*Fetch y_in: */
103 if(mxIsDouble(prhs[2])){
104 y_inm=mxGetPr(prhs[2]);
[12060]105 y_in=(double*)xmalloc(nods*sizeof(double));
[1]106 for (i=0;i<nods;i++){
107 y_in[i]=y_inm[i];
108 }
109 }
110 else{
111 printf("%s%s\n",__FUNCT__," error message: third argument should be the y corrdinate list!");
112 mexErrMsgTxt(" ");
113 }
114
115 /*Fetch segments_in: */
116 if(mxIsDouble(prhs[3])){
117 num_seg=mxGetM(prhs[3]);
118 tsegments_in=mxGetPr(prhs[3]);
[12060]119 segments_in=(double*)xmalloc(num_seg*3*sizeof(double));
[1]120 for (i=0;i<num_seg;i++){
121 for (j=0;j<3;j++){
122 *(segments_in+3*i+j)=*(tsegments_in+num_seg*j+i);
123 }
124 }
125 }
126 else{
127 printf("%s%s\n",__FUNCT__," error message: fourth argument should be the segments list!");
128 mexErrMsgTxt(" ");
129 }
130
131 /*Fetch segment markers: */
132 if(mxIsDouble(prhs[4])){
133 tsegmentmarkers_in=mxGetPr(prhs[4]);
[12060]134 segmentmarkers_in=(double*)xmalloc(num_seg*sizeof(double));
[1]135 for (i=0;i<num_seg;i++){
136 segmentmarkers_in[i]=tsegmentmarkers_in[i];
137 }
138 }
139 else{
140 printf("%s%s\n",__FUNCT__," error message: fourth argument should be the segmentmarkers list!");
141 mexErrMsgTxt(" ");
142 }
143
[3755]144 /*
[1]145 printf("Index: \n");
146 for (i=0;i<nel;i++){
147 for(j=0;j<3;j++){
148 printf("%lf ",*(index_in+3*i+j));
149 }
150 printf("\n");
151 }
152 printf("x,y: \n");
153 for (i=0;i<nods;i++){
154 printf("%16.16lf %16.16lf\n",x_in[i],y_in[i]);
155 }
156 printf("segments:\n");
157 for (i=0;i<num_seg;i++){
158 for(j=0;j<3;j++){
159 printf("%lf ",*(segments_in+3*i+j));
160 }
161 printf("%lf ",segmentmarkers_in[i]);
162 printf("\n");
163 }
[3755]164 */
[1]165
166 /*First, do some fixing on the existing mesh: we do not want any element belonging entirely to the segment list (ie:
[8306]167 *all the nodes of this element belong to the segments (tends to happen when there are corners: */
[1]168 RemoveCornersFromRifts(&index_in,&nel,&x_in,&y_in,&nods,segments_in,segmentmarkers_in,num_seg);
169
170 /*Figure out if we have rifts, and how many: */
171 IsRiftPresent(&riftflag,&numrifts,segmentmarkers_in,num_seg);
172
173 if(riftflag){
174 SplitMeshForRifts(&nel,&index_in,&nods,&x_in,&y_in,&num_seg,&segments_in,&segmentmarkers_in);
175 }
176
177 /*Order segments so that their normals point outside the domain: */
[4453]178 OrderSegments(&segments_in,num_seg, index_in,nel);
[1]179
180 if(riftflag){
181
182 /*We do not want to output segments mixed with rift segments: wring out the rifts from the segments, using the
183 *segmentmarkerlist:*/
[12081]184 SplitRiftSegments(&segments_in,&segmentmarkers_in,&num_seg,&out_numrifts,&out_riftsnumsegments,&out_riftssegments,numrifts,nods,nel);
[1]185
186 /*Using rift segments, associate rift faces in pairs, each pair face representing opposite flanks of the rifts facing one another directly: */
187 PairRiftElements(&out_riftsnumpairs,&out_riftspairs,out_numrifts,out_riftsnumsegments,out_riftssegments,x_in,y_in);
188
189 /*Order rifts so that they start from one tip, go to the other tip, and back: */
[12071]190 OrderRifts(&out_riftstips, out_riftssegments,out_riftspairs,numrifts,out_riftsnumsegments,x_in,y_in,nods,nel);
[1]191
192 /*Create penalty pairs, used by Imp: */
193 PenaltyPairs(&out_riftspenaltypairs,&out_riftsnumpenaltypairs,numrifts,out_riftssegments,out_riftsnumsegments,out_riftspairs,out_riftstips,x_in,y_in);
194 }
195
196
197 /*Output : */
[12085]198 WriteData(&plhs[0],index_in,nel,3);
199 //pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
200 //mxSetM(pmxa_array,3);
201 //mxSetN(pmxa_array,nel);
202 //mxSetPr(pmxa_array,index_in);
203 //mexCallMATLAB( 1, &plhs[0], 1, &pmxa_array, "transpose");
[1]204
[12085]205 WriteData(&plhs[1],x_in,nods,1);
206 //pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
207 //mxSetM(pmxa_array,1);
208 //mxSetN(pmxa_array,nods);
209 //mxSetPr(pmxa_array,x_in);
210 //mexCallMATLAB( 1, &plhs[1], 1, &pmxa_array, "transpose");
[1]211
[12085]212 WriteData(&plhs[2],y_in,nods,1);
213 //pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
214 //mxSetM(pmxa_array,1);
215 //mxSetN(pmxa_array,nods);
216 //mxSetPr(pmxa_array,y_in);
217 //mexCallMATLAB( 1, &plhs[2], 1, &pmxa_array, "transpose");
[1]218
[12085]219 WriteData(&plhs[3],segments_in,num_seg,3);
220//pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
221//mxSetM(pmxa_array,3);
222//mxSetN(pmxa_array,num_seg);
223//mxSetPr(pmxa_array,segments_in);
224//mexCallMATLAB( 1, &plhs[3], 1, &pmxa_array, "transpose");
[1]225
[12085]226 WriteData(&plhs[4],segmentmarkers_in,num_seg,1);
227 //pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
228 //mxSetM(pmxa_array,1);
229 //mxSetN(pmxa_array,num_seg);
230 //mxSetPr(pmxa_array,segmentmarkers_in);
231 //mexCallMATLAB( 1, &plhs[4], 1, &pmxa_array, "transpose");
[1]232
233 if(riftflag){
234 /*Create a structure rifts where if i is a rift number, we have the following fields rifts(i).segments and rifts(i).numsegs: */
235
236 fnames[0] = "numsegs";
237 fnames[1] = "segments";
238 fnames[2] = "pairs";
239 fnames[3] = "tips";
240 fnames[4] = "penaltypairs";
241 fnames[5] = "fill";
242 fnames[6] = "friction";
[1784]243 fnames[7] = "fraction";
[1805]244 fnames[8] = "fractionincrement";
[6748]245 fnames[9] = "state";
[1]246
247 dimensions[0]=out_numrifts;
248
[6748]249 pmxa_array=mxCreateStructArray( ndim,dimensions,10,fnames);
[1]250
251 for (i=0;i<out_numrifts;i++){
[12085]252
[1]253 /*Segments: */
[12085]254 WriteData(&pmxa_array3,out_riftssegments[i],out_riftsnumsegments[i],3);
255 //pmxa_array2= mxCreateDoubleMatrix(0,0,mxREAL);
256 //mxSetM(pmxa_array2,3);
257 //mxSetN(pmxa_array2,out_riftsnumsegments[i]);
258 //mxSetPr(pmxa_array2,out_riftssegments[i]);
259 //mexCallMATLAB( 1, &pmxa_array3, 1, &pmxa_array2, "transpose");
[1]260 mxSetField(pmxa_array,i,"segments",pmxa_array3);
[8963]261 mxSetField(pmxa_array,i,"numsegs",mxCreateDoubleScalar((double)out_riftsnumsegments[i]));
[1]262
263 /*Element pairs: */
[12085]264 WriteData(&pmxa_array3,out_riftspairs[i],out_riftsnumpairs[i],2);
265 //pmxa_array2= mxCreateDoubleMatrix(0,0,mxREAL);
266 //mxSetM(pmxa_array2,2);
267 //mxSetN(pmxa_array2,out_riftsnumpairs[i]);
268 //mxSetPr(pmxa_array2,out_riftspairs[i]);
269 //mexCallMATLAB( 1, &pmxa_array3, 1, &pmxa_array2, "transpose");
[1]270 mxSetField(pmxa_array,i,"pairs",pmxa_array3);
271
272 /*Tips: */
273 pmxa_array2= mxCreateDoubleMatrix(0,0,mxREAL);
274 mxSetM(pmxa_array2,1);
[12085]275 pair=(double*)mxMalloc(2*sizeof(double));
[1]276 pair[0]=*(out_riftstips+2*i+0);
277 pair[1]=*(out_riftstips+2*i+1);
278 mxSetN(pmxa_array2,2);
279 mxSetPr(pmxa_array2,pair);
280 mxSetField(pmxa_array,i,"tips",pmxa_array2);
281
282 /*Penalty pairs: */
[12085]283 WriteData(&pmxa_array3,out_riftspenaltypairs[i],out_riftsnumpenaltypairs[i],7);
284 //pmxa_array2= mxCreateDoubleMatrix(0,0,mxREAL);
285 //mxSetM(pmxa_array2,7);
286 //mxSetN(pmxa_array2,out_riftsnumpenaltypairs[i]);
287 //mxSetPr(pmxa_array2,out_riftspenaltypairs[i]);
288 //mexCallMATLAB( 1, &pmxa_array3, 1, &pmxa_array2, "transpose");
[1]289 mxSetField(pmxa_array,i,"penaltypairs",pmxa_array3);
290
[1805]291 /*Friction fraction, fractionincrement and fill: */
[1]292 mxSetField(pmxa_array,i,"friction",mxCreateDoubleScalar(0));
[3567]293 mxSetField(pmxa_array,i,"fill",mxCreateDoubleScalar(IceEnum)); //default is ice
[1784]294 mxSetField(pmxa_array,i,"fraction",mxCreateDoubleScalar(0)); //default is ice
[1805]295 mxSetField(pmxa_array,i,"fractionincrement",mxCreateDoubleScalar(0.1));
[6748]296
297 /*State: */
[12085]298 state=(double*)mxMalloc(out_riftsnumpenaltypairs[i]*sizeof(double));
[6748]299 for(j=0;j<out_riftsnumpenaltypairs[i];j++)state[j]=FreeEnum;
300 pmxa_array2= mxCreateDoubleMatrix(0,0,mxREAL);
301 mxSetM(pmxa_array2,1);
302 mxSetN(pmxa_array2,out_riftsnumpenaltypairs[i]);
303 mxSetPr(pmxa_array2,state);
304 mexCallMATLAB( 1, &pmxa_array3, 1, &pmxa_array2, "transpose");
305
306 mxSetField(pmxa_array,i,"state",pmxa_array3);
[1]307 }
308 }
309 else{
310 /*output NaN :*/
[12085]311 pNaN=(double*)mxMalloc(sizeof(double));
[3567]312 *pNaN=NAN;
[1]313 pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
314 mxSetM(pmxa_array,1);
315 mxSetN(pmxa_array,1);
316 mxSetPr(pmxa_array,pNaN);
317
318 }
319 plhs[5]=pmxa_array;
320
321 return;
322}
323
324
325void TriMeshProcessRiftsUsage(void)
326{
327 printf("\n");
328 printf(" usage: [index2,x2,y2,segments2,segmentmarkers2,rifts2]=TriMeshProcessrifts(index1,x1,y1,segments1,segmentmarkers1) \n");
329 printf(" where: (index1,x1,y1,segments1,segmentmarkers1) is an initial triangulation.\n");
330 printf(" index2,x2,y2,segments2,segmentmarkers2,rifts2 is the resulting triangulation where rifts have been processed.\n");
331}
Note: See TracBrowser for help on using the repository browser.