Ice Sheet System Model  4.18
Code documentation
Functions
triangle.h File Reference
#include <stdio.h>
#include <math.h>

Go to the source code of this file.

Functions

int AssociateSegmentToElement (int **psegments, int nseg, int *index, int nel)
 
int OrderSegments (int **psegments, int nseg, int *index, int nel)
 
int GridInsideHole (double *px0, double *py0, int n, double *x, double *y)
 
int FindElement (int A, int B, int *index, int nel)
 
int SplitMeshForRifts (int *pnel, int **pindex, int *pnods, double **px, double **py, int *pnsegs, int **psegments, int **psegmentmarkerlist)
 
int IsGridOnRift (int *riftsegments, int nriftsegs, int node)
 
int GridElementsList (int **pGridElements, int *pNumGridElements, int node, double *index, int nel)
 
int IsNeighbor (int el1, int el2, int *index)
 
int IsOnRift (int el, int nriftsegs, int *riftsegments)
 
void RiftSegmentsFromSegments (int *pnriftsegs, int **priftsegments, int nel, int *index, int nsegs, int *segments)
 
int DetermineGridElementListOnOneSideOfRift (int *pNumGridElementListOnOneSideOfRift, int **pGridElementListOnOneSideOfRift, int segmentnumber, int nriftsegs, int *riftsegments, int node, int *index, int nel)
 
int UpdateSegments (int **psegments, int **psegmentmarkerlist, int *pnsegs, int *index, double *x, double *y, int *riftsegments, int nriftsegs, int nods, int nel)
 
int FindElement (double A, double B, int *index, int nel)
 
int IsRiftPresent (int *priftflag, int *pnumrifts, int *segmentmarkerlist, int nsegs)
 
int SplitRiftSegments (int **psegments, int **psegmentmarkerlist, int *pnumsegs, int *pnumrifts, int **priftsnumsegs, int ***priftssegments, int numrifts, int nods, int nels)
 
int OrderRifts (int **priftstips, int **riftssegments, int **riftspairs, int numrifts, int *riftsnumsegments, double *x, double *y, int nods, int nels)
 
int PenaltyPairs (double ***priftspenaltypairs, int **priftsnumpenaltypairs, int numrifts, int **riftssegments, int *riftsnumsegments, int **riftspairs, int *riftstips, double *x, double *y)
 
int RemoveCornersFromRifts (int **pindex, int *pnel, double **px, double **py, int *pnods, int *segments, int *segmentmarkers, int num_seg)
 
int PairRiftElements (int **priftsnumpairs, int ***priftspairs, int numrifts, int *riftsnumsegments, int **riftssegments, double *x, double *y)
 

Function Documentation

◆ AssociateSegmentToElement()

int AssociateSegmentToElement ( int **  psegments,
int  nseg,
int *  index,
int  nel 
)

Definition at line 7 of file AssociateSegmentToElement.cpp.

7  {
8 
9  /*node indices: */
10  int A,B;
11 
12  /*Recover segments: */
13  int* segments=*psegments;
14 
15  for(int i=0;i<nseg;i++){
16  A=segments[3*i+0];
17  B=segments[3*i+1];
18  segments[3*i+2]=FindElement(A,B,index,nel)+1; //matlab indexing.
19  }
20 
21  /*Assign output pointers: */
22  *psegments=segments;
23  return 1;
24 }

◆ OrderSegments()

int OrderSegments ( int **  psegments,
int  nseg,
int *  index,
int  nel 
)

Definition at line 7 of file OrderSegments.cpp.

7  {
8 
9  /*vertex indices: */
10  int A,B;
11 
12  /*element index*/
13  int el;
14 
15  /*Recover segments: */
16  int* segments=*psegments;
17 
18  for(int i=0;i<nseg;i++){
19  A=segments[3*i+0];
20  B=segments[3*i+1];
21  el=segments[3*i+2]-1; //after AssociateSegmentToElement, el was a matlab index, we need the c index now.
22 
23  if (index[3*el+0]==A){
24  if (index[3*el+2]==B){
25  segments[3*i+0]=B;
26  segments[3*i+1]=A;
27  }
28  }
29  else if (index[3*el+1]==A){
30  if (index[3*el+0]==B){
31  segments[3*i+0]=B;
32  segments[3*i+1]=A;
33  }
34  }
35  else{
36  if (index[3*el+1]==B){
37  segments[3*i+0]=B;
38  segments[3*i+1]=A;
39  }
40  }
41  }
42 
43  /*Assign output pointers: */
44  *psegments=segments;
45  return 1;
46 }

◆ GridInsideHole()

int GridInsideHole ( double *  px0,
double *  py0,
int  n,
double *  x,
double *  y 
)

Definition at line 14 of file GridInsideHole.cpp.

14  {
15 
16  double flag=0.0;
17  double xA,xB,xC,xD,xE;
18  double yA,yB,yC,yD,yE;
19 
20  /*Take first and last vertices: */
21  xA=x[0];
22  yA=y[0];
23  xB=x[n-1];
24  yB=y[n-1];
25 
26  /*Figure out middle of segment [A B]: */
27  xC=(xA+xB)/2;
28  yC=(yA+yB)/2;
29 
30  /*D and E are on each side of segment [A B], on the median line between segment [A B],
31  *at an angle of 10 degree (less than the minimum 30 enforced by the quality of the mesh: */
32  xD=xC+tan(10./180.*M_PI)*(yC-yA);
33  yD=yC+tan(10./180.*M_PI)*(xA-xC);
34  xE=xC-tan(10./180.*M_PI)*(yC-yA);
35  yE=yC-tan(10./180.*M_PI)*(xA-xC);
36 
37  /*Either E or D is inside profile (x,y): */
38  IsInPolySerial(&flag,&xD,&yD,1,x,y,n,2);
39  /*FIXME: used to be 'flag' and not '!flag', check*/
40  if(!flag){
41  /*D is inside the poly: */
42  *px0=xD;
43  *py0=yD;
44  }
45  else{
46  /*E is inside the poly: */
47  *px0=xE;
48  *py0=yE;
49  }
50  return 1;
51 }

◆ FindElement() [1/2]

int FindElement ( int  A,
int  B,
int *  index,
int  nel 
)

Definition at line 339 of file TriangleUtils.cpp.

339  {/*{{{*/
340 
341  int el=-1;
342  for (int n=0;n<nel;n++){
343  if(((index[3*n+0]==A) || (index[3*n+1]==A) || (index[3*n+2]==A)) && ((index[3*n+0]==B) || (index[3*n+1]==B) || (index[3*n+2]==B))){
344  el=n;
345  break;
346  }
347  }
348  return el;
349 }/*}}}*/

◆ SplitMeshForRifts()

int SplitMeshForRifts ( int *  pnel,
int **  pindex,
int *  pnods,
double **  px,
double **  py,
int *  pnsegs,
int **  psegments,
int **  psegmentmarkerlist 
)

Definition at line 7 of file SplitMeshForRifts.cpp.

7  {
8 
9  /*Some notes on dimensions:
10  index of size nelx3
11  x and y of size nodsx1
12  segments of size nsegsx3*/
13 
14  int i,j,k,l;
15  int node;
16  int el;
17  int nriftsegs;
18  int* riftsegments=NULL;
19  int* flags=NULL;
20  int NumGridElementListOnOneSideOfRift;
21  int* GridElementListOnOneSideOfRift=NULL;
22 
23  /*Recover input: */
24  int nel = *pnel;
25  int *index = *pindex;
26  int nods = *pnods;
27  double *x = *px;
28  double *y = *py;
29  int nsegs = *pnsegs;
30  int *segments = *psegments;
31  int *segmentmarkerlist = *psegmentmarkerlist;
32 
33  /*Establish list of segments that belong to a rift: */
34  /*riftsegments of size nriftsegsx4 (4 for first element on segment,second element,first node and second snode)*/
35  RiftSegmentsFromSegments(&nriftsegs,&riftsegments,nel,index,nsegs,segments);
36 
37  /*Go through all nodes of the rift segments, and start splitting the mesh: */
38  flags=xNewZeroInit<int>(nods); //to make sure we don't split the same nodes twice!
39  for (i=0;i<nriftsegs;i++){
40  for (j=0;j<2;j++){
41 
42  node=riftsegments[4*i+j+2];
43  if(flags[node-1]){
44  /*This node was already split, skip:*/
45  continue;
46  }
47  else{
48  flags[node-1]=1;
49  }
50 
51  if(IsGridOnRift(riftsegments,nriftsegs,node)){
52 
53  DetermineGridElementListOnOneSideOfRift(&NumGridElementListOnOneSideOfRift,&GridElementListOnOneSideOfRift,i,nriftsegs,riftsegments,node,index,nel);
54 
55  /*Summary: we have for node, a list of elements
56  * (GridElementListOnOneSideOfRift, of size
57  * NumGridElementListOnOneSideOfRift) that all contain node
58  *and that are on the same side of the rift. For all these
59  elements, we clone node into another node, and we swap all
60  instances of node in the triangulation *for those elements, to the
61  new node.*/
62 
63  //create new node
64  x=xReNew<double>(x,nods,nods+1);
65  y=xReNew<double>(y,nods,nods+1);
66  x[nods]=x[node-1]; //matlab indexing
67  y[nods]=y[node-1]; //matlab indexing
68 
69  //augment number of nodes
70  nods++;
71 
72  //change elements owning this node
73  for (k=0;k<NumGridElementListOnOneSideOfRift;k++){
74  el=GridElementListOnOneSideOfRift[k];
75  for (l=0;l<3;l++){
76  if (index[3*el+l]==node) index[3*el+l]=nods; //again, matlab indexing.
77  }
78  }
79  }// if(IsGridOnRift(riftsegments,nriftsegs,node))
80  } //for(j=0;j<2;j++)
81  } //for (i=0;i<nriftsegs;i++)
82 
83  /*update segments: they got modified completely by adding new nodes.*/
84  UpdateSegments(&segments,&segmentmarkerlist, &nsegs,index,x,y,riftsegments,nriftsegs,nods,nel);
85 
86  /*Assign output pointers: */
87  *pnel=nel;
88  *pindex=index;
89  *pnods=nods;
90  *px=x;
91  *py=y;
92  *pnsegs=nsegs;
93  *psegments=segments;
94  *psegmentmarkerlist=segmentmarkerlist;
95  return 1;
96 }

◆ IsGridOnRift()

int IsGridOnRift ( int *  riftsegments,
int  nriftsegs,
int  node 
)

Definition at line 12 of file TriangleUtils.cpp.

12  {/*{{{*/
13 
14  /*Does this node belong to 4 elements, or just 2? If it belongs to 4 elements, it is inside a rift,
15  *if it belongs to 2 elements, it is on the tip of a rift, or it has already been split across the rift (see below).*/
16 
17  int i;
18  int j;
19  int count;
20 
21  count=0;
22  for (i=0;i<nriftsegs;i++){
23  for (j=0;j<2;j++){
24  if ((*(riftsegments+4*i+2+j))==node) count++;
25  }
26  }
27  if (count==2){
28  return 1;
29  }
30  else{
31  return 0;
32  }
33 }/*}}}*/

◆ GridElementsList()

int GridElementsList ( int **  pGridElements,
int *  pNumGridElements,
int  node,
double *  index,
int  nel 
)

◆ IsNeighbor()

int IsNeighbor ( int  el1,
int  el2,
int *  index 
)

Definition at line 87 of file TriangleUtils.cpp.

87  {/*{{{*/
88  /*From a triangulation held in index, figure out if elements 1 and 2 have two nodes in common: */
89  int i,j;
90  int count=0;
91  for (i=0;i<3;i++){
92  for (j=0;j<3;j++){
93  if (index[3*el1+i]==index[3*el2+j])count++;
94  }
95  }
96  if (count==2){
97  return 1;
98  }
99  else{
100  return 0;
101  }
102 }/*}}}*/

◆ IsOnRift()

int IsOnRift ( int  el,
int  nriftsegs,
int *  riftsegments 
)

Definition at line 103 of file TriangleUtils.cpp.

103  {/*{{{*/
104  /*From a list of elements segments, figure out if el belongs to it: */
105  int i;
106  for (i=0;i<nriftsegs;i++){
107  if ((*(riftsegments+4*i+0)==el) || (*(riftsegments+4*i+1)==el)){
108  return 1;
109  }
110  }
111  return 0;
112 }/*}}}*/

◆ RiftSegmentsFromSegments()

void RiftSegmentsFromSegments ( int *  pnriftsegs,
int **  priftsegments,
int  nel,
int *  index,
int  nsegs,
int *  segments 
)

Definition at line 113 of file TriangleUtils.cpp.

113  {/*{{{*/
114 
115  int i,counter;
116  int el,el2;
117 
118  int nriftsegs;
119  int* riftsegments=NULL;
120  int* riftsegments_uncompressed=NULL;
121  int element_nodes[3];
122 
123  /*Allocate segmentflags: */
124  riftsegments_uncompressed=xNewZeroInit<int>(nsegs*5);
125 
126  /*Find the segments that belong to a rift: they are the ones that see two elements. The other ones belong to a boundary
127  *or a hole: */
128  nriftsegs=0;
129  for (i=0;i<nsegs;i++){
130  el=(int)*(segments+3*i+2)-1; //element found in AssociateSegmentToElements
131  /*Temporarily set nodes belonging to the segments to -1 in the triangulation index, and
132  *then proceed to find another element that owns the segment. If we don't find it, we know
133  *we are dealing with a boundary or hole, otherwise, we are dealing with a rift: */
134  element_nodes[0]=*(index+3*el+0);
135  element_nodes[1]=*(index+3*el+1);
136  element_nodes[2]=*(index+3*el+2);
137 
138  index[3*el+0]=-1;
139  index[3*el+1]=-1;
140  index[3*el+2]=-1;
141 
142  el2=FindElement(*(segments+3*i+0),*(segments+3*i+1),index,nel);
143 
144  /*Restore index: */
145  index[3*el+0]=element_nodes[0];
146  index[3*el+1]=element_nodes[1];
147  index[3*el+2]=element_nodes[2];
148 
149  if (el2!=-1){
150  /*el and el2 are on a segment rift, facing one another, plug them into riftsegments_uncompressed: */
151  riftsegments_uncompressed[5*i+0]=1;
152  riftsegments_uncompressed[5*i+1]=el;
153  riftsegments_uncompressed[5*i+2]=el2;
154  riftsegments_uncompressed[5*i+3]=segments[3*i+0];
155  riftsegments_uncompressed[5*i+4]=segments[3*i+1];
156  nriftsegs++;
157  }
158  }
159 
160  /*Compress riftsegments_uncompressed:*/
161  riftsegments=xNew<int>(nriftsegs*4);
162  counter=0;
163  for (i=0;i<nsegs;i++){
164  if (riftsegments_uncompressed[5*i+0]){
165  riftsegments[counter*4+0]=riftsegments_uncompressed[5*i+1];
166  riftsegments[counter*4+1]=riftsegments_uncompressed[5*i+2];
167  riftsegments[counter*4+2]=riftsegments_uncompressed[5*i+3];
168  riftsegments[counter*4+3]=riftsegments_uncompressed[5*i+4];
169  counter++;
170  }
171  }
172  xDelete<int>(riftsegments_uncompressed);
173 
174  /*Assign output pointers: */
175  *priftsegments=riftsegments;
176  *pnriftsegs=nriftsegs;
177 }/*}}}*/

◆ DetermineGridElementListOnOneSideOfRift()

int DetermineGridElementListOnOneSideOfRift ( int *  pNumGridElementListOnOneSideOfRift,
int **  pGridElementListOnOneSideOfRift,
int  segmentnumber,
int  nriftsegs,
int *  riftsegments,
int  node,
int *  index,
int  nel 
)

Definition at line 178 of file TriangleUtils.cpp.

178  {/*{{{*/
179 
180  int noerr=1;
181  int k,l,counter;
182  int newel;
183 
184  int* GridElements=NULL;
185  int NumGridElements;
186 
187  /*Output: */
188  int NumGridElementListOnOneSideOfRift;
189  int* GridElementListOnOneSideOfRift=NULL;
190 
191  /*Build a list of all the elements connected to this node: */
192  GridElementsList(&GridElements,&NumGridElements,node,index,nel);
193 
194  /*Figure out the list of elements that are on the same side of the rift. To do so, we start from one
195  * side of the rift and keep rotating in the same direction:*/
196  GridElementListOnOneSideOfRift=xNew<int>(NumGridElements);
197  //bootstrap the GridElementListOnOneSideOfRift by filling elements from riftsegments: */
198  GridElementListOnOneSideOfRift[0]=*(riftsegments+4*segmentnumber+0); /*this one does not belong to the same side, but is just there
199  for a rotation direction, we 'll take it out when we are
200  done rotating*/
201  GridElementListOnOneSideOfRift[1]=*(riftsegments+4*segmentnumber+1);
202  counter=1;
203  for (;;){
204  /*Find neighbour of element GridElementListOnOneSideOfRift[counter], not
205  * equal to GridElementListOnOneSideOfRift[counter-1]*/
206  for (k=0;k<NumGridElements;k++){
207  if(IsNeighbor(GridElements[k],GridElementListOnOneSideOfRift[counter],index)){
208  /*Verify this element is not already in our list of element on the same side of the rift: */
209  newel=1;
210  for (l=0;l<=counter;l++){
211  if (GridElements[k]==GridElementListOnOneSideOfRift[l]){
212  newel=0;
213  break;
214  }
215  }
216  if (newel){
217  counter++;
218  GridElementListOnOneSideOfRift[counter]=GridElements[k];
219  if (IsOnRift(GridElements[k],nriftsegs,riftsegments)){
220  break;
221  }
222  k=-1;
223  }
224  }
225  }
226  /*Reduce counter by 1 and get rift of first element in GridElementListOnOneSideOfRift:*/
227  NumGridElementListOnOneSideOfRift=counter;
228  for (l=0;l<NumGridElementListOnOneSideOfRift;l++){
229  GridElementListOnOneSideOfRift[l]=GridElementListOnOneSideOfRift[l+1];
230  }
231  break;
232  }// for (;;)
233 
234  /*Free ressources: */
235  xDelete<int>(GridElements);
236  /*Assign output pointers: */
237  *pNumGridElementListOnOneSideOfRift=NumGridElementListOnOneSideOfRift;
238  *pGridElementListOnOneSideOfRift=GridElementListOnOneSideOfRift;
239  return noerr;
240 }/*}}}*/

◆ UpdateSegments()

int UpdateSegments ( int **  psegments,
int **  psegmentmarkerlist,
int *  pnsegs,
int *  index,
double *  x,
double *  y,
int *  riftsegments,
int  nriftsegs,
int  nods,
int  nel 
)

Definition at line 241 of file TriangleUtils.cpp.

241  {/*{{{*/
242 
243  int noerr=1;
244  int i,j,k;
245  int el1,el2;
246 
247  int *segments = NULL;
248  int *segmentmarkerlist = NULL;
249  int nsegs;
250 
251  /*Recover input: */
252  segments = *psegments;
253  segmentmarkerlist = *psegmentmarkerlist;
254  nsegs = *pnsegs;
255 
256  /*Reallocate segments: */
257  segments =xReNew<int>(segments, nsegs*3,(nsegs+nriftsegs)*3);
258  segmentmarkerlist=xReNew<int>(segmentmarkerlist,nsegs,(nsegs+nriftsegs));
259 
260  /*First, update the existing segments to the new nodes :*/
261  for (i=0;i<nriftsegs;i++){
262  el1=riftsegments[4*i+0];
263  el2=riftsegments[4*i+1];
264  for (j=0;j<nsegs;j++){
265  if (segments[3*j+2]==(el1+1)){
266  /*segment j is the same as rift segment i.Let's update segments[j][:] using element el1 and the corresponding rift segment.
267  *Because riftsegments does not represent a list of rift segments anymore (it got heavily modified in SplitElementsForRifts,
268  *we can only rely on the position (x,y) of the rift nodes to create a segment:*/
269  for (k=0;k<3;k++){
270  if ((x[*(index+el1*3+k)-1]==x[*(segments+3*j+0)-1]) && (y[*(index+el1*3+k)-1]==y[*(segments+3*j+0)-1])){
271  *(segments+3*j+0)=*(index+el1*3+k); _assert_(segments[3*j+0]<nods+1);
272  break;
273  }
274  }
275  for (k=0;k<3;k++){
276  if ((x[*(index+el1*3+k)-1]==x[*(segments+3*j+1)-1]) && (y[*(index+el1*3+k)-1]==y[*(segments+3*j+1)-1])){
277  *(segments+3*j+1)=*(index+el1*3+k); _assert_(segments[3*j+1]<nods+1);
278  break;
279  }
280  }
281  /*Deal with el2: */
282  *(segments+3*(nsegs+i)+2)=el2+1;
283  *(segmentmarkerlist+(nsegs+i))=*(segmentmarkerlist+j);
284  for (k=0;k<3;k++){
285  if ((x[*(index+el2*3+k)-1]==x[*(segments+3*j+0)-1]) && (y[*(index+el2*3+k)-1]==y[*(segments+3*j+0)-1])){
286  *(segments+3*(nsegs+i)+0)=*(index+el2*3+k); _assert_(segments[3*(nsegs+i)+0]<nods+1);
287  break;
288  }
289  }
290  for (k=0;k<3;k++){
291  if ((x[*(index+el2*3+k)-1]==x[*(segments+3*j+1)-1]) && (y[*(index+el2*3+k)-1]==y[*(segments+3*j+1)-1])){
292  *(segments+3*(nsegs+i)+1)=*(index+el2*3+k); _assert_(segments[3*(nsegs+i)+1]<nods+1);
293  break;
294  }
295  }
296  }
297  if (*(segments+3*j+2)==(el2+1)){
298  /*segment j is the same as rift segment i.*/
299  /*Let's update segments[j][:] using element el2 and the corresponding rift segment: */
300  for (k=0;k<3;k++){
301  if ((x[*(index+el2*3+k)-1]==x[*(segments+3*j+0)-1]) && (y[*(index+el2*3+k)-1]==y[*(segments+3*j+0)-1])){
302  *(segments+3*j+0)=*(index+el2*3+k); _assert_(segments[3*j+0]<nods+1);
303  break;
304  }
305  }
306  for (k=0;k<3;k++){
307  if ((x[*(index+el2*3+k)-1]==x[*(segments+3*j+1)-1]) && (y[*(index+el2*3+k)-1]==y[*(segments+3*j+1)-1])){
308  *(segments+3*j+1)=*(index+el2*3+k);_assert_(segments[3*j+1]<nods+1);
309  break;
310  }
311  }
312  /*Deal with el1: */
313  *(segments+3*(nsegs+i)+2)=el1+1;
314  *(segmentmarkerlist+(nsegs+i))=*(segmentmarkerlist+j);
315  for (k=0;k<3;k++){
316  if ((x[*(index+el1*3+k)-1]==x[*(segments+3*j+0)-1]) && (y[*(index+el1*3+k)-1]==y[*(segments+3*j+0)-1])){
317  *(segments+3*(nsegs+i)+0)=*(index+el1*3+k);_assert_(segments[3*(nsegs+i)+0]<nods+1);
318  break;
319  }
320  }
321  for (k=0;k<3;k++){
322  if ((x[*(index+el1*3+k)-1]==x[*(segments+3*j+1)-1]) && (y[*(index+el1*3+k)-1]==y[*(segments+3*j+1)-1])){
323  *(segments+3*(nsegs+i)+1)=*(index+el1*3+k);_assert_(segments[3*(nsegs+i)+1]<nods+1);
324  break;
325  }
326  }
327  }
328  }
329  }
330  nsegs+=nriftsegs;
331 
332  /*Assign output pointers: */
333  *psegments=segments;
334  *psegmentmarkerlist=segmentmarkerlist;
335  *pnsegs=nsegs;
336 
337  return noerr;
338 }/*}}}*/

◆ FindElement() [2/2]

int FindElement ( double  A,
double  B,
int *  index,
int  nel 
)

◆ IsRiftPresent()

int IsRiftPresent ( int *  priftflag,
int *  pnumrifts,
int *  segmentmarkerlist,
int  nsegs 
)

Definition at line 482 of file TriangleUtils.cpp.

482  {/*{{{*/
483 
484  int i;
485  int noerr=1;
486 
487  /*output: */
488  int riftflag=0;
489  int numrifts=0;
490 
491  int maxmark=1; //default marker for regular segments
492 
493  /*Any marker >=2 indicates a certain rift: */
494  numrifts=0;
495  for (i=0;i<nsegs;i++){
496  if (segmentmarkerlist[i]>maxmark){
497  numrifts++;
498  maxmark=segmentmarkerlist[i];
499  }
500  }
501  if(numrifts)riftflag=1;
502 
503  /*Assign output pointers:*/
504  *priftflag=riftflag;
505  *pnumrifts=numrifts;
506  return noerr;
507 }/*}}}*/

◆ SplitRiftSegments()

int SplitRiftSegments ( int **  psegments,
int **  psegmentmarkerlist,
int *  pnumsegs,
int *  pnumrifts,
int **  priftsnumsegs,
int ***  priftssegments,
int  numrifts,
int  nods,
int  nels 
)

Definition at line 350 of file TriangleUtils.cpp.

350  {/*{{{*/
351 
352  /*Using segment markers, wring out the rift segments from the segments. Rift markers are
353  *of the form 2+i where i=0 to number of rifts */
354 
355  int noerr=1;
356  int i,j,counter;
357 
358  /*input: */
359  int *segments = NULL;
360  int *segmentmarkerlist = NULL;
361  int numsegs;
362 
363  /*output: */
364  int new_numsegs;
365  int *riftsnumsegs = NULL;
366  int **riftssegments = NULL;
367  int *new_segments = NULL;
368  int *new_segmentmarkers = NULL;
369 
370  /*intermediary: */
371  int* riftsegment=NULL;
372 
373  /*Recover input arguments: */
374  segments = *psegments;
375  numsegs = *pnumsegs;
376  segmentmarkerlist = *psegmentmarkerlist;
377 
378  /*First, figure out how many segments will be left in 'segments': */
379  counter=0;
380  for (i=0;i<numsegs;i++){
381  if (segmentmarkerlist[i]==1)counter++; //1 is default marker for non-rifts;
382  }
383  /*Allocate new segments: */
384  new_numsegs=counter;
385  new_segments=xNew<int>(new_numsegs*3);
386  new_segmentmarkers=xNew<int>(new_numsegs);
387 
388  /*Copy new segments info : */
389  counter=0;
390  for (i=0;i<numsegs;i++){
391  if (segmentmarkerlist[i]==1){
392  new_segments[3*counter+0]=segments[3*i+0];
393  new_segments[3*counter+1]=segments[3*i+1];
394  new_segments[3*counter+2]=segments[3*i+2];
395  new_segmentmarkers[counter]=segmentmarkerlist[i];
396  counter++;
397  }
398  }
399 
400  /*Now deal with rift segments: */
401  riftsnumsegs=xNew<int>(numrifts);
402  riftssegments=xNew<int*>(numrifts);
403  for (i=0;i<numrifts;i++){
404  /*Figure out how many segments for rift i: */
405  counter=0;
406  for (j=0;j<numsegs;j++){
407  if (segmentmarkerlist[j]==2+i)counter++;
408  }
409  riftsnumsegs[i]=counter;
410  riftsegment=xNew<int>(counter*3);
411  /*Copy new segments info :*/
412  counter=0;
413  for (j=0;j<numsegs;j++){
414  if (segmentmarkerlist[j]==(2+i)){
415  riftsegment[3*counter+0]=segments[3*j+0];_assert_(riftsegment[3*counter+0]<nods+1);
416  riftsegment[3*counter+1]=segments[3*j+1];_assert_(riftsegment[3*counter+1]<nods+1);
417  riftsegment[3*counter+2]=segments[3*j+2];_assert_(riftsegment[3*counter+2]<nel+1);
418  counter++;
419  }
420  }
421  *(riftssegments+i)=riftsegment;
422  }
423 
424  /*Free ressources: */
425  xDelete<int>(segments);
426 
427  /*Assign output pointers: */
428  *psegments=new_segments;
429  *psegmentmarkerlist=new_segmentmarkers;
430  *pnumsegs=new_numsegs;
431  *pnumrifts=numrifts;
432  *priftssegments=riftssegments;
433  *priftsnumsegs=riftsnumsegs;
434  return noerr;
435 }/*}}}*/

◆ OrderRifts()

int OrderRifts ( int **  priftstips,
int **  riftssegments,
int **  riftspairs,
int  numrifts,
int *  riftsnumsegments,
double *  x,
double *  y,
int  nods,
int  nels 
)

Definition at line 508 of file TriangleUtils.cpp.

508  {/*{{{*/
509 
510  int noerr=1;
511  int i,j,k,counter;
512 
513  /*intermediary: */
514  int *riftsegments = NULL;
515  int *riftpairs = NULL;
516  int numsegs;
517 
518  /*ordering and copy: */
519  int *order = NULL;
520  int *riftsegments_copy = NULL;
521  int *riftpairs_copy = NULL;
522 
523  /*node and element manipulation: */
524  int node1,node2,node3,node4,temp_node,tip1,tip2,node;
525  int el2;
526  int already_ordered=0;
527 
528  /*output: */
529  int* riftstips=NULL;
530 
531  /*Allocate byproduct of this routine, riftstips: */
532  riftstips=xNew<int>(numrifts*2);
533 
534  /*Go through all rifts: */
535  for (i=0;i<numrifts;i++){
536  riftsegments = riftssegments[i];
537  riftpairs = riftspairs[i];
538  numsegs = riftsnumsegments[i];
539 
540  /*Allocate copy of riftsegments and riftpairs,
541  *as well as ordering vector: */
542  riftsegments_copy=xNew<int>(numsegs*3);
543  riftpairs_copy=xNew<int>(numsegs*2);
544  order=xNew<int>(numsegs);
545 
546  /*First find the tips, using the pairs. If a pair of elements has one node in common, this node is a rift tip: */
547  tip1=-1;
548  tip2=-1;
549 
550  for (j=0;j<numsegs;j++){
551  el2=*(riftpairs+2*j+1);
552  node1=*(riftsegments+3*j+0);
553  node2=*(riftsegments+3*j+1);
554  /*Summary, el1 and el2 are facing one another across the rift. node1 and node2 belong to el1 and
555  *are located on the rift. Find node3 and node4, nodes belonging to el2 and located on the rift: */
556  for (k=0;k<numsegs;k++){
557  if (*(riftsegments+3*k+2)==el2){
558  node3=*(riftsegments+3*k+0);
559  node4=*(riftsegments+3*k+1);
560  break;
561  }
562  }
563  /* Make sure node3 faces node1 and node4 faces node2: */
564  _assert_(node1<nods+1 && node4<nods+1);
565  _assert_(node1>0 && node4>0);
566  if ((x[node1-1]==x[node4-1]) && (y[node1-1]==y[node4-1])){
567  /*Swap node3 and node4:*/
568  temp_node=node3;
569  node3=node4;
570  node4=temp_node;
571  }
572 
573  /*Figure out if a tip is on this element: */
574  if (node3==node1){
575  /*node1 is a tip*/
576  if (tip1==-1) {
577  tip1=node1;
578  continue;
579  }
580  if ((tip2==-1) && (node1!=tip1)){
581  tip2=node1;
582  break;
583  }
584  }
585 
586  if (node4==node2){
587  /*node2 is a tip*/
588  if (tip1==-1){
589  tip1=node2;
590  continue;
591  }
592  if ((tip2==-1) && (node2!=tip1)){
593  tip2=node2;
594  break;
595  }
596  }
597  }
598 
599  /*Record tips in riftstips: */
600  *(riftstips+2*i+0)=tip1;
601  *(riftstips+2*i+1)=tip2;
602 
603  /*We have the two tips for this rift. Go from tip1 to tip2, and figure out the order in which segments are sequential.
604  *Because two elements are connected to tip1, we chose one first, which defines the direction we are rotating along the rift. */
605  node=tip1;
606  for (counter=0;counter<numsegs;counter++){
607  for (j=0;j<numsegs;j++){
608  node1=*(riftsegments+3*j+0);
609  node2=*(riftsegments+3*j+1);
610 
611  if ((node1==node) || (node2==node)){
612  /*Ok, this segment is connected to node, plug its index into order, unless we already plugged it before: */
613  already_ordered=0;
614  for (k=0;k<counter;k++){
615  if(order[k]==j){
616  already_ordered=1;
617  break;
618  }
619  }
620  if (!already_ordered){
621  order[counter]=j;
622  if(node1==node){
623  node=node2;
624  }
625  else if(node2==node){
626  node=node1;
627  }
628  break;
629  }
630  }
631  }
632  }
633 
634  /*Using the order vector, and the riftsegments_copy and riftspairs_copy, reorder the segments and the pairs: */
635  for (j=0;j<numsegs;j++){
636  _assert_(order[j]<numsegs);
637  *(riftsegments_copy+3*j+0)=*(riftsegments+3*order[j]+0);
638  *(riftsegments_copy+3*j+1)=*(riftsegments+3*order[j]+1);
639  *(riftsegments_copy+3*j+2)=*(riftsegments+3*order[j]+2);
640  *(riftpairs_copy+2*j+0)=*(riftpairs+2*order[j]+0);
641  *(riftpairs_copy+2*j+1)=*(riftpairs+2*order[j]+1);
642  }
643 
644  for (j=0;j<numsegs;j++){
645  *(riftsegments+3*j+0)=*(riftsegments_copy+3*j+0);
646  *(riftsegments+3*j+1)=*(riftsegments_copy+3*j+1);
647  *(riftsegments+3*j+2)=*(riftsegments_copy+3*j+2);
648  *(riftpairs+2*j+0)=*(riftpairs_copy+2*j+0);
649  *(riftpairs+2*j+1)=*(riftpairs_copy+2*j+1);
650  }
651 
652  xDelete<int>(order);
653  xDelete<int>(riftsegments_copy);
654  xDelete<int>(riftpairs_copy);
655 
656  }
657 
658  /*Assign output pointer:*/
659  *priftstips=riftstips;
660  return noerr;
661 }/*}}}*/

◆ PenaltyPairs()

int PenaltyPairs ( double ***  priftspenaltypairs,
int **  priftsnumpenaltypairs,
int  numrifts,
int **  riftssegments,
int *  riftsnumsegments,
int **  riftspairs,
int *  riftstips,
double *  x,
double *  y 
)

Definition at line 662 of file TriangleUtils.cpp.

663  {
664 
665  int noerr=1;
666  int i,j,k,k0;
667 
668  double el1,el2,node1,node2,node3,node4;
669  double temp_node;
670 
671  /*output: */
672  double **riftspenaltypairs = NULL;
673  double *riftpenaltypairs = NULL;
674  int *riftsnumpenaltypairs = NULL;
675 
676  /*intermediary: */
677  int numsegs;
678  int* riftsegments=NULL;
679  int* riftpairs=NULL;
680  int counter;
681  double normal[2];
682  double length;
683  int k1,k2;
684 
685  /*Allocate: */
686  riftspenaltypairs=xNew<double*>(numrifts);
687  riftsnumpenaltypairs=xNew<int>(numrifts);
688 
689  for(i=0;i<numrifts;i++){
690  numsegs=riftsnumsegs[i];
691  riftsegments=riftssegments[i];
692  riftpairs=riftspairs[i];
693 
694  /*allocate riftpenaltypairs, and riftnumpenaltypairs: */
695  if((numsegs/2-1)!=0)riftpenaltypairs=xNewZeroInit<double>((numsegs/2-1)*RIFTPENALTYPAIRSWIDTH);
696 
697  /*Go through only one flank of the rifts, not counting the tips: */
698  counter=0;
699  for(j=0;j<(numsegs/2);j++){
700  el1=*(riftpairs+2*j+0);
701  el2=*(riftpairs+2*j+1);
702  node1=*(riftsegments+3*j+0);
703  node2=*(riftsegments+3*j+1);
704  /*Find segment index to recover node3 and node4, facing node1 and node2: */
705  k0=-1;
706  for(k=0;k<numsegs;k++){
707  if(*(riftsegments+3*k+2)==el2){
708  k0=k;
709  break;
710  }
711  }
712  node3=*(riftsegments+3*k0+0);
713  node4=*(riftsegments+3*k0+1);
714 
715  /* Make sure node3 faces node1 and node4 faces node2: */
716  if ((x[(int)node1-1]==x[(int)node4-1]) && (y[(int)node1-1]==y[(int)node4-1])){
717  /*Swap node3 and node4:*/
718  temp_node=node3;
719  node3=node4;
720  node4=temp_node;
721  }
722  /*Ok, we have node1 facing node3, and node2 facing node4. Compute the normal to
723  *this segment, and its length: */
724  normal[0]=cos(atan2(x[(int)node1-1]-x[(int)node2-1],y[(int)node2-1]-y[(int)node1-1]));
725  normal[1]=sin(atan2(x[(int)node1-1]-x[(int)node2-1],y[(int)node2-1]-y[(int)node1-1]));
726  length=sqrt(pow(x[(int)node2-1]-x[(int)node1-1],(double)2)+pow(y[(int)node2-1]-y[(int)node1-1],(double)2));
727 
728  /*Be careful here, we want penalty loads on each node, not on each segment. This means we cannot plug node1,
729  * node2, node3 and node4 directly into riftpenaltypairs. We need to include node1, node2, node3 and node4,
730  * only once. We'll add the normals and the lengths : */
731 
732  if(node1!=node3){ //exclude tips from loads
733  k1=-1;
734  for(k=0;k<counter;k++){
735  if( (*(riftpenaltypairs+k*7+0))==node1){
736  k1=k;
737  break;
738  }
739  }
740  if(k1==-1){
741  *(riftpenaltypairs+counter*7+0)=node1;
742  *(riftpenaltypairs+counter*7+1)=node3;
743  *(riftpenaltypairs+counter*7+2)=el1;
744  *(riftpenaltypairs+counter*7+3)=el2;
745  *(riftpenaltypairs+counter*7+4)=normal[0];
746  *(riftpenaltypairs+counter*7+5)=normal[1];
747  *(riftpenaltypairs+counter*7+6)=length/2;
748  counter++;
749  }
750  else{
751  *(riftpenaltypairs+k1*7+4)+=normal[0];
752  *(riftpenaltypairs+k1*7+5)+=normal[1];
753  *(riftpenaltypairs+k1*7+6)+=length/2;
754  }
755  }
756  if(node2!=node4){
757  k2=-1;
758  for(k=0;k<counter;k++){
759  if( (*(riftpenaltypairs+k*7+0))==node2){
760  k2=k;
761  break;
762  }
763  }
764  if(k2==-1){
765  *(riftpenaltypairs+counter*7+0)=node2;
766  *(riftpenaltypairs+counter*7+1)=node4;
767  *(riftpenaltypairs+counter*7+2)=el1;
768  *(riftpenaltypairs+counter*7+3)=el2;
769  *(riftpenaltypairs+counter*7+4)=normal[0];
770  *(riftpenaltypairs+counter*7+5)=normal[1];
771  *(riftpenaltypairs+counter*7+6)=length/2;
772  counter++;
773  }
774  else{
775  *(riftpenaltypairs+k2*7+4)+=normal[0];
776  *(riftpenaltypairs+k2*7+5)+=normal[1];
777  *(riftpenaltypairs+k2*7+6)+=length/2;
778  }
779  }
780  }
781  /*Renormalize normals: */
782  for(j=0;j<counter;j++){
783  double magnitude=sqrt(pow( double(riftpenaltypairs[j*7+4]),2) + pow( double(riftpenaltypairs[j*7+5]),2) );
784  *(riftpenaltypairs+j*7+4)=*(riftpenaltypairs+j*7+4)/magnitude;
785  *(riftpenaltypairs+j*7+5)=*(riftpenaltypairs+j*7+5)/magnitude;
786  }
787 
788  riftspenaltypairs[i]=riftpenaltypairs;
789  riftsnumpenaltypairs[i]=(numsegs/2-1);
790  }
791 
792  /*Assign output pointers: */
793  *priftspenaltypairs=riftspenaltypairs;
794  *priftsnumpenaltypairs=riftsnumpenaltypairs;
795  return noerr;
796 }/*}}}*/

◆ RemoveCornersFromRifts()

int RemoveCornersFromRifts ( int **  pindex,
int *  pnel,
double **  px,
double **  py,
int *  pnods,
int *  segments,
int *  segmentmarkers,
int  num_seg 
)

Definition at line 797 of file TriangleUtils.cpp.

797  {/*{{{*/
798 
799  int noerr=1;
800  int i,j,k;
801  int node1,node2,node3;
802  int el;
803  double pair[2];
804  int pair_count=0;
805  int triple=0;
806 
807  /*Recover input: */
808  int *index = *pindex;
809  int nel = *pnel;
810  double *x = *px;
811  double *y = *py;
812  int nods = *pnods;
813 
814  for (i=0;i<num_seg;i++){
815  node1=*(segments+3*i+0);
816  node2=*(segments+3*i+1);
817  /*Find all elements connected to [node1 node2]: */
818  pair_count=0;
819  for (j=0;j<nel;j++){
820  if (*(index+3*j+0)==node1){
821  if ((*(index+3*j+1)==node2) || (*(index+3*j+2)==node2)){
822  pair[pair_count]=j;
823  pair_count++;
824  }
825  }
826  if (*(index+3*j+1)==node1){
827  if ((*(index+3*j+0)==node2) || (*(index+3*j+2)==node2)){
828  pair[pair_count]=j;
829  pair_count++;
830  }
831  }
832  if (*(index+3*j+2)==node1){
833  if ((*(index+3*j+0)==node2) || (*(index+3*j+1)==node2)){
834  pair[pair_count]=j;
835  pair_count++;
836  }
837  }
838  }
839  /*Ok, we have pair_count elements connected to this segment. For each of these elements,
840  *figure out if the third node also belongs to a segment: */
841  if ((pair_count==0) || (pair_count==1)){ //we only select the rift segments, which belong to 2 elements
842  continue;
843  }
844  else{
845  for (j=0;j<pair_count;j++){
846  el=(int)pair[j];
847  triple=0;
848  /*First find node3: */
849  if (*(index+3*el+0)==node1){
850  if (*(index+3*el+1)==node2)node3=*(index+3*el+2);
851  else node3=*(index+3*el+1);
852  }
853  if (*(index+3*el+1)==node1){
854  if (*(index+3*el+0)==node2)node3=*(index+3*el+2);
855  else node3=*(index+3*el+0);
856  }
857  if (*(index+3*el+2)==node1){
858  if (*(index+3*el+0)==node2)node3=*(index+3*el+1);
859  else node3=*(index+3*el+0);
860  }
861  /*Ok, we have node3. Does node3 belong to a segment? : */
862  for (k=0;k<num_seg;k++){
863  if ((node3==*(segments+3*k+0)) || (node3==*(segments+3*k+1))){
864  triple=1;
865  break;
866  }
867  }
868  if(triple==1){
869  /*el is a corner element: we need to split it in 3 triangles: */
870  x=xReNew<double>(x,nods,nods+1);
871  y=xReNew<double>(y,nods,nods+1);
872  x[nods]=(x[(int)node1-1]+x[(int)node2-1]+x[(int)node3-1])/3;
873  y[nods]=(y[(int)node1-1]+y[(int)node2-1]+y[(int)node3-1])/3;
874  index=xReNew<int>(index,nel*3,(nel+2*3));
875  /*First, reassign element el: */
876  *(index+3*el+0)=node1;
877  *(index+3*el+1)=node2;
878  *(index+3*el+2)=nods+1;
879  /*Other two elements: */
880  *(index+3*nel+0)=node2;
881  *(index+3*nel+1)=node3;
882  *(index+3*nel+2)=nods+1;
883 
884  *(index+3*(nel+1)+0)=node3;
885  *(index+3*(nel+1)+1)=node1;
886  *(index+3*(nel+1)+2)=nods+1;
887  /*we need to change the segment elements corresponding to el: */
888  for (k=0;k<num_seg;k++){
889  if (*(segments+3*k+2)==(el+1)){
890  if ( ((*(segments+3*k+0)==node1) && (*(segments+3*k+1)==node2)) || ((*(segments+3*k+0)==node2) && (*(segments+3*k+1)==node1))) *(segments+3*k+2)=el+1;
891  if ( ((*(segments+3*k+0)==node2) && (*(segments+3*k+1)==node3)) || ((*(segments+3*k+0)==node3) && (*(segments+3*k+1)==node2))) *(segments+3*k+2)=nel+1;
892  if ( ((*(segments+3*k+0)==node3) && (*(segments+3*k+1)==node1)) || ((*(segments+3*k+0)==node1) && (*(segments+3*k+1)==node3))) *(segments+3*k+2)=nel+2;
893  }
894  }
895 
896  nods=nods+1;
897  nel=nel+2;
898  i=0;
899  break;
900  }
901  } //for (j=0;j<pair_count;j++)
902  }
903  }// for (i=0;i<num_seg;i++)
904 
905  /*Assign output pointers: */
906  *pindex=index;
907  *pnel=nel;
908  *px=x;
909  *py=y;
910  *pnods=nods;
911  return noerr;
912 }/*}}}*/

◆ PairRiftElements()

int PairRiftElements ( int **  priftsnumpairs,
int ***  priftspairs,
int  numrifts,
int *  riftsnumsegments,
int **  riftssegments,
double *  x,
double *  y 
)

Definition at line 436 of file TriangleUtils.cpp.

436  {/*{{{*/
437 
438  int noerr=1;
439  int i,j,k;
440 
441  /*output: */
442  int *riftsnumpairs = NULL;
443  int **riftspairs = NULL;
444 
445  /*intermediary :*/
446  int numsegs;
447  int* segments=NULL;
448  int* pairs=NULL;
449  int node1,node2,node3,node4;
450 
451  riftsnumpairs=xNew<int>(numrifts);
452  riftspairs=xNew<int*>(numrifts);
453  for (i=0;i<numrifts;i++){
454  segments=riftssegments[i];
455  numsegs =riftsnumsegments[i];
456  riftsnumpairs[i]=numsegs;
457  pairs=xNew<int>(2*numsegs);
458  for (j=0;j<numsegs;j++){
459  pairs[2*j+0]=segments[3*j+2]; //retrieve element to which this segment belongs.
460  node1=segments[3*j+0]-1; node2=segments[3*j+1]-1;
461  /*Find element facing on other side of rift: */
462  for (k=0;k<numsegs;k++){
463  if (k==j)continue;
464  node3=segments[3*k+0]-1; node4=segments[3*k+1]-1;
465  /*We are trying to find 2 elements, where position of node3 == position of node1, and position of node4 == position of node2*/
466  if ( ((x[node3]==x[node1]) && (y[node3]==y[node1]) && (x[node4]==x[node2]) && (y[node4]==y[node2]))
467  || ((x[node3]==x[node2]) && (y[node3]==y[node2]) && (x[node4]==x[node1]) && (y[node4]==y[node1])) ){
468  /*We found the corresponding element: */
469  pairs[2*j+1]=segments[3*k+2];
470  break;
471  }
472  }
473  }
474  riftspairs[i]=pairs;
475  }
476 
477  /*Assign output pointers: */
478  *priftsnumpairs=riftsnumpairs;
479  *priftspairs=riftspairs;
480  return noerr;
481 }/*}}}*/
IsGridOnRift
int IsGridOnRift(int *riftsegments, int nriftsegs, int node)
Definition: TriangleUtils.cpp:12
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
DetermineGridElementListOnOneSideOfRift
int DetermineGridElementListOnOneSideOfRift(int *pNumGridElementListOnOneSideOfRift, int **pGridElementListOnOneSideOfRift, int segmentnumber, int nriftsegs, int *riftsegments, int node, int *index, int nel)
Definition: TriangleUtils.cpp:178
IsNeighbor
int IsNeighbor(int el1, int el2, int *index)
Definition: TriangleUtils.cpp:87
GridElementsList
int GridElementsList(int **pGridElements, int *pNumGridElements, int node, int *index, int nel)
Definition: TriangleUtils.cpp:34
FindElement
int FindElement(int A, int B, int *index, int nel)
Definition: TriangleUtils.cpp:339
UpdateSegments
int UpdateSegments(int **psegments, int **psegmentmarkerlist, int *pnsegs, int *index, double *x, double *y, int *riftsegments, int nriftsegs, int nods, int nel)
Definition: TriangleUtils.cpp:241
RiftSegmentsFromSegments
void RiftSegmentsFromSegments(int *pnriftsegs, int **priftsegments, int nel, int *index, int nsegs, int *segments)
Definition: TriangleUtils.cpp:113
IsOnRift
int IsOnRift(int el, int nriftsegs, int *riftsegments)
Definition: TriangleUtils.cpp:103
RIFTPENALTYPAIRSWIDTH
#define RIFTPENALTYPAIRSWIDTH
Definition: TriangleUtils.cpp:11
IsInPolySerial
int IsInPolySerial(double *in, double *xc, double *yc, int numvertices, double *x, double *y, int nods, int edgevalue)
Definition: exp.cpp:45
FindElement
int FindElement(int A, int B, int *index, int nel)
Definition: TriangleUtils.cpp:339
M_PI
#define M_PI
Definition: GridInsideHole.cpp:12
Generated on Thu Jul 2 2020 08:09:25 for Ice Sheet System Model by   doxygen 1.8.19