Ice Sheet System Model  4.18
Code documentation
Macros | Functions
TriangleUtils.cpp File Reference
#include <stdio.h>
#include "./triangle.h"
#include "../Exceptions/exceptions.h"
#include "../MemOps/MemOps.h"

Go to the source code of this file.

Macros

#define RIFTPENALTYPAIRSWIDTH   8
 

Functions

int IsGridOnRift (int *riftsegments, int nriftsegs, int node)
 
int GridElementsList (int **pGridElements, int *pNumGridElements, int node, int *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 (int A, int B, int *index, int nel)
 
int SplitRiftSegments (int **psegments, int **psegmentmarkerlist, int *pnumsegs, int *pnumrifts, int **priftsnumsegs, int ***priftssegments, int numrifts, int nods, int nel)
 
int PairRiftElements (int **priftsnumpairs, int ***priftspairs, int numrifts, int *riftsnumsegments, int **riftssegments, double *x, double *y)
 
int IsRiftPresent (int *priftflag, int *pnumrifts, int *segmentmarkerlist, int nsegs)
 
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 *riftsnumsegs, 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)
 

Macro Definition Documentation

◆ RIFTPENALTYPAIRSWIDTH

#define RIFTPENALTYPAIRSWIDTH   8

Definition at line 11 of file TriangleUtils.cpp.

Function Documentation

◆ 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,
int *  index,
int  nel 
)

Definition at line 34 of file TriangleUtils.cpp.

34  {/*{{{*/
35 
36  /*From a node, recover all the elements that are connected to it: */
37  int i,j;
38  int noerr=1;
39 
40  int max_number_elements=12;
41  int current_size;
42  int NumGridElements;
43  int* GridElements=NULL;
44  int* GridElementsRealloc=NULL;
45 
46  /*From a mesh with 30 degrees minimum angle, we get 12 possible elements that own
47  * the node. We start by allocating GridElements with that size, and realloc
48  * more if needed.*/
49 
50  current_size=max_number_elements;
51  NumGridElements=0;
52  GridElements=xNew<int>(max_number_elements);
53 
54  for (i=0;i<nel;i++){
55  for (j=0;j<3;j++){
56  if (index[3*i+j]==node){
57  if (NumGridElements<=(current_size-1)){
58  GridElements[NumGridElements]=i;
59  NumGridElements++;
60  break;
61  }
62  else{
63  /*Reallocate another max_number_elements slots in the GridElements: */
64  GridElementsRealloc=xReNew<int>(GridElements,current_size,(current_size+max_number_elements));
65  if (!GridElementsRealloc){
66  noerr=0;
67  goto cleanup_and_return;
68  }
69  current_size+=max_number_elements;
70  GridElements=GridElementsRealloc;
71  GridElements[NumGridElements]=i;
72  NumGridElements++;
73  break;
74  }
75  }
76  }
77  }
78  cleanup_and_return:
79  if(!noerr){
80  xDelete<int>(GridElements);
81  }
82  /*Allocate return pointers: */
83  *pGridElements=GridElements;
84  *pNumGridElements=NumGridElements;
85  return noerr;
86 }/*}}}*/

◆ 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()

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 }/*}}}*/

◆ SplitRiftSegments()

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

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 }/*}}}*/

◆ 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 }/*}}}*/

◆ 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 }/*}}}*/

◆ 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 *  riftsnumsegs,
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 }/*}}}*/
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
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
IsOnRift
int IsOnRift(int el, int nriftsegs, int *riftsegments)
Definition: TriangleUtils.cpp:103
RIFTPENALTYPAIRSWIDTH
#define RIFTPENALTYPAIRSWIDTH
Definition: TriangleUtils.cpp:11
Generated on Thu Jul 2 2020 08:09:25 for Ice Sheet System Model by   doxygen 1.8.19