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
Line 
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;
29 const char* fnames[10];
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
49 /* state: */
50 double* state=NULL;
51
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 _error_("bad usage");
69 }
70
71 /*Fetch index_in: */
72 if(mxIsDouble(prhs[0])){
73 nel=mxGetM(prhs[0]);
74 tindex_in=mxGetPr(prhs[0]);
75 index_in=(double*)xmalloc(nel*3*sizeof(double));
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{
83 _error_("first argument should be the element list");
84 }
85
86 /*Fetch x_in: */
87 if(mxIsDouble(prhs[1])){
88 nods=mxGetM(prhs[1]);
89 x_inm=mxGetPr(prhs[1]);
90 x_in=(double*)xmalloc(nods*sizeof(double));
91 for (i=0;i<nods;i++){
92 x_in[i]=x_inm[i];
93 }
94 }
95 else{
96 _error_("second argument should be the x corrdinate list");
97 }
98
99 /*Fetch y_in: */
100 if(mxIsDouble(prhs[2])){
101 y_inm=mxGetPr(prhs[2]);
102 y_in=(double*)xmalloc(nods*sizeof(double));
103 for (i=0;i<nods;i++){
104 y_in[i]=y_inm[i];
105 }
106 }
107 else{
108 _error_("third argument should be the y corrdinate list");
109 }
110
111 /*Fetch segments_in: */
112 if(mxIsDouble(prhs[3])){
113 num_seg=mxGetM(prhs[3]);
114 tsegments_in=mxGetPr(prhs[3]);
115 segments_in=(double*)xmalloc(num_seg*3*sizeof(double));
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{
123 _error_("fourth argument should be the segments list");
124 }
125
126 /*Fetch segment markers: */
127 if(mxIsDouble(prhs[4])){
128 tsegmentmarkers_in=mxGetPr(prhs[4]);
129 segmentmarkers_in=(double*)xmalloc(num_seg*sizeof(double));
130 for (i=0;i<num_seg;i++){
131 segmentmarkers_in[i]=tsegmentmarkers_in[i];
132 }
133 }
134 else{
135 _error_("fourth argument should be the segmentmarkers list");
136 }
137
138 /*First, do some fixing on the existing mesh: we do not want any element belonging entirely to the segment list (ie:
139 *all the nodes of this element belong to the segments (tends to happen when there are corners: */
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: */
150 OrderSegments(&segments_in,num_seg, index_in,nel);
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:*/
156 SplitRiftSegments(&segments_in,&segmentmarkers_in,&num_seg,&out_numrifts,&out_riftsnumsegments,&out_riftssegments,numrifts,nods,nel);
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: */
162 OrderRifts(&out_riftstips, out_riftssegments,out_riftspairs,numrifts,out_riftsnumsegments,x_in,y_in,nods,nel);
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 : */
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");
176
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");
183
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");
190
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");
197
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");
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";
215 fnames[7] = "fraction";
216 fnames[8] = "fractionincrement";
217 fnames[9] = "state";
218
219 dimensions[0]=out_numrifts;
220
221 pmxa_array=mxCreateStructArray( ndim,dimensions,10,fnames);
222
223 for (i=0;i<out_numrifts;i++){
224
225 /*Segments: */
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");
232 mxSetField(pmxa_array,i,"segments",pmxa_array3);
233 mxSetField(pmxa_array,i,"numsegs",mxCreateDoubleScalar((double)out_riftsnumsegments[i]));
234
235 /*Element pairs: */
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");
242 mxSetField(pmxa_array,i,"pairs",pmxa_array3);
243
244 /*Tips: */
245 pmxa_array2= mxCreateDoubleMatrix(0,0,mxREAL);
246 mxSetM(pmxa_array2,1);
247 pair=(double*)mxMalloc(2*sizeof(double));
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: */
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");
261 mxSetField(pmxa_array,i,"penaltypairs",pmxa_array3);
262
263 /*Friction fraction, fractionincrement and fill: */
264 mxSetField(pmxa_array,i,"friction",mxCreateDoubleScalar(0));
265 mxSetField(pmxa_array,i,"fill",mxCreateDoubleScalar(IceEnum)); //default is ice
266 mxSetField(pmxa_array,i,"fraction",mxCreateDoubleScalar(0)); //default is ice
267 mxSetField(pmxa_array,i,"fractionincrement",mxCreateDoubleScalar(0.1));
268
269 /*State: */
270 state=(double*)mxMalloc(out_riftsnumpenaltypairs[i]*sizeof(double));
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);
279 }
280 }
281 else{
282 /*output NaN :*/
283 pNaN=(double*)mxMalloc(sizeof(double));
284 *pNaN=NAN;
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{
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.");
303}
Note: See TracBrowser for help on using the repository browser.