source: issm/trunk-jpl/src/modules/InternalFront/InternalFront.cpp@ 13236

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

CHG: preparing files for python inclusion
File size: 2.8 KB
Line 
1/*\file InternalFront.c
2 *\brief: build pressureload
3 */
4
5#include "./InternalFront.h"
6
7void InternalFrontUsage(void) {/*{{{*/
8 _pprintLine_("");
9 _pprintLine_(" usage: icefront = " << __FUNCT__ << "(md);");
10 _pprintLine_("");
11}/*}}}*/
12WRAPPER(InternalFront){
13
14 bool* elementonwater=NULL;
15 int* elements=NULL;
16 int* connectivity=NULL;
17 int* elementconnectivity=NULL;
18 int* front=NULL;
19 double* front2=NULL;
20 bool found;
21 int numberofelements,numberofsegments;
22 int N,M;
23 int i,j,ii,jj,id;
24
25 /*Boot module: */
26 MODULEBOOT();
27
28 /*checks on arguments on the matlab side: */
29 CheckNumMatlabArguments(nlhs,NLHS,nrhs,NRHS,__FUNCT__,&InternalFrontUsage);
30
31 /*Fetch required fields*/
32 FetchData(&numberofelements,mxGetAssignedField(MODEL,0,"numberofelements"));
33 if(numberofelements<=0) _error_("No elements found in the model");
34 FetchData(&elements,&M,&N,mxGetAssignedField(MODEL,0,"elements"));
35 if(M!=numberofelements || N!=3) _error_("Field 'elements' should be of size [md.numberofelements 3]");
36 FetchData(&elementonwater,&M,&N,mxGetAssignedField(MODEL,0,"elementonwater"));
37 if(M!=numberofelements || N!=1) _error_("Field 'elementonwater' should be of size [md.numberofelements 1]");
38 FetchData(&elementconnectivity,&M,&N,mxGetAssignedField(MODEL,0,"elementconnectivity"));
39 if(M!=numberofelements || N!=3) _error_("Field 'elementconnectivity' should be of size [md.numberofelements 3]");
40
41 /*Allocate and initialize all variables*/
42 numberofsegments=0;
43 front=xNew<int>(3*numberofelements*4);
44
45 /*Loop over all elements on water*/
46 for(i=0;i<numberofelements;i++){
47
48 /*Skip if on water*/
49 if(!elementonwater[i]) continue;
50
51 /*Loop over all three adjacent elements*/
52 for(j=0;j<3;j++){
53
54 /*Skip if adjacent element does not exist or is on water*/
55 id=elementconnectivity[i*3+j];
56 if(id==0) continue;
57 if(elementonwater[id-1]) continue;
58
59 /*We have an ice front to add here, let's go!*/
60 for(ii=0;ii<3;ii++){
61
62 found=false;
63 for(jj=0;jj<3;jj++){
64 if(elements[(id-1)*3+ii]==elements[i*3+jj]){
65 found=true;
66 break;
67 }
68 }
69
70 /*OK, we just found the node of id, which is not in i. We have the segment*/
71 if(!found){
72 front[numberofsegments*4+0]=elements[(id-1)*3+(ii+1)%3];
73 front[numberofsegments*4+1]=elements[(id-1)*3+(ii+2)%3];
74 front[numberofsegments*4+2]=id;
75 front[numberofsegments*4+3]=IceEnum;
76 numberofsegments+=1;
77 break;
78 }
79 }
80
81 /*In debugging mode, check that the segment has been found*/
82 _assert_(!found);
83 }
84 }
85
86 /*Now that we know how many segments there is we can allocate the final matrix*/
87 if(numberofsegments){
88 front2=xNew<double>(4*numberofsegments);
89 for(i=0;i<4*numberofsegments;i++) front2[i]=(double)front[i];
90 }
91 xfree((void**)&front);
92
93 /*write output datasets: */
94 WriteData(FRONT,front2,numberofsegments,4);
95
96 /*end module: */
97 MODULEEND();
98}
Note: See TracBrowser for help on using the repository browser.