source: issm/trunk/src/c/shared/Triangle/TriangleUtils.cpp@ 27232

Last change on this file since 27232 was 27232, checked in by Mathieu Morlighem, 3 years ago

merged trunk-jpl and trunk for revision 27230
File size: 27.1 KB
RevLine 
[1]1/*
[22498]2 * TriangleUtils: mesh manipulation routines:
[1]3 */
4
[3252]5#include <stdio.h>
6
[22498]7#include "./triangle.h"
[3335]8#include "../Exceptions/exceptions.h"
[14904]9#include "../MemOps/MemOps.h"
[1]10
[6748]11#define RIFTPENALTYPAIRSWIDTH 8
[18063]12int IsGridOnRift(int* riftsegments, int nriftsegs, int node){/*{{{*/
[1]13
[8301]14 /*Does this node belong to 4 elements, or just 2? If it belongs to 4 elements, it is inside a rift,
[1]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).*/
[13622]16
[1]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++){
[8301]24 if ((*(riftsegments+4*i+2+j))==node) count++;
[1]25 }
26 }
27 if (count==2){
28 return 1;
29 }
30 else{
31 return 0;
32 }
[12080]33}/*}}}*/
[18063]34int GridElementsList(int** pGridElements, int* pNumGridElements,int node,int* index,int nel){/*{{{*/
[1]35
[8301]36 /*From a node, recover all the elements that are connected to it: */
[1]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
[8301]47 * the node. We start by allocating GridElements with that size, and realloc
[1]48 * more if needed.*/
49
50 current_size=max_number_elements;
51 NumGridElements=0;
[12436]52 GridElements=xNew<int>(max_number_elements);
[1]53
54 for (i=0;i<nel;i++){
55 for (j=0;j<3;j++){
[14222]56 if (index[3*i+j]==node){
[1]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: */
[12447]64 GridElementsRealloc=xReNew<int>(GridElements,current_size,(current_size+max_number_elements));
[1]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){
[12436]80 xDelete<int>(GridElements);
[1]81 }
82 /*Allocate return pointers: */
83 *pGridElements=GridElements;
84 *pNumGridElements=NumGridElements;
85 return noerr;
[12080]86}/*}}}*/
[18063]87int IsNeighbor(int el1,int el2,int* index){/*{{{*/
[8301]88 /*From a triangulation held in index, figure out if elements 1 and 2 have two nodes in common: */
[1]89 int i,j;
90 int count=0;
91 for (i=0;i<3;i++){
92 for (j=0;j<3;j++){
[14222]93 if (index[3*el1+i]==index[3*el2+j])count++;
[1]94 }
95 }
96 if (count==2){
97 return 1;
98 }
99 else{
100 return 0;
101 }
[12080]102}/*}}}*/
[18063]103int IsOnRift(int el,int nriftsegs,int* riftsegments){/*{{{*/
[1]104 /*From a list of elements segments, figure out if el belongs to it: */
105 int i;
106 for (i=0;i<nriftsegs;i++){
[16218]107 if ((*(riftsegments+4*i+0)==el) || (*(riftsegments+4*i+1)==el)){
[1]108 return 1;
109 }
110 }
111 return 0;
[12080]112}/*}}}*/
[18063]113void RiftSegmentsFromSegments(int* pnriftsegs, int** priftsegments, int nel,int* index,int nsegs,int* segments){/*{{{*/
[13622]114
[1]115 int i,counter;
116 int el,el2;
[13622]117
[1]118 int nriftsegs;
119 int* riftsegments=NULL;
120 int* riftsegments_uncompressed=NULL;
[14222]121 int element_nodes[3];
[1]122
123 /*Allocate segmentflags: */
[12446]124 riftsegments_uncompressed=xNewZeroInit<int>(nsegs*5);
[1]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
[8301]131 /*Temporarily set nodes belonging to the segments to -1 in the triangulation index, and
[1]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: */
[8301]134 element_nodes[0]=*(index+3*el+0);
135 element_nodes[1]=*(index+3*el+1);
136 element_nodes[2]=*(index+3*el+2);
[1]137
[14222]138 index[3*el+0]=-1;
139 index[3*el+1]=-1;
140 index[3*el+2]=-1;
[1]141
142 el2=FindElement(*(segments+3*i+0),*(segments+3*i+1),index,nel);
143
144 /*Restore index: */
[14222]145 index[3*el+0]=element_nodes[0];
146 index[3*el+1]=element_nodes[1];
147 index[3*el+2]=element_nodes[2];
[1]148
149 if (el2!=-1){
150 /*el and el2 are on a segment rift, facing one another, plug them into riftsegments_uncompressed: */
[14222]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++;
[1]157 }
158 }
159
160 /*Compress riftsegments_uncompressed:*/
[12436]161 riftsegments=xNew<int>(nriftsegs*4);
[1]162 counter=0;
163 for (i=0;i<nsegs;i++){
[14222]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];
[1]169 counter++;
170 }
171 }
[12436]172 xDelete<int>(riftsegments_uncompressed);
[13622]173
[1]174 /*Assign output pointers: */
175 *priftsegments=riftsegments;
176 *pnriftsegs=nriftsegs;
[12080]177}/*}}}*/
[18063]178int DetermineGridElementListOnOneSideOfRift(int* pNumGridElementListOnOneSideOfRift, int** pGridElementListOnOneSideOfRift, int segmentnumber, int nriftsegs, int* riftsegments, int node,int* index,int nel){/*{{{*/
[1]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
[8301]191 /*Build a list of all the elements connected to this node: */
192 GridElementsList(&GridElements,&NumGridElements,node,index,nel);
[3595]193
[1]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:*/
[12436]196 GridElementListOnOneSideOfRift=xNew<int>(NumGridElements);
[1]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
[27232]234 /*Free resources: */
[12436]235 xDelete<int>(GridElements);
[1]236 /*Assign output pointers: */
237 *pNumGridElementListOnOneSideOfRift=NumGridElementListOnOneSideOfRift;
238 *pGridElementListOnOneSideOfRift=GridElementListOnOneSideOfRift;
239 return noerr;
[12080]240}/*}}}*/
[18063]241int UpdateSegments(int** psegments,int** psegmentmarkerlist, int* pnsegs,int* index, double* x,double* y,int* riftsegments,int nriftsegs,int nods,int nel){/*{{{*/
[1]242
243 int noerr=1;
244 int i,j,k;
245 int el1,el2;
246
[14222]247 int *segments = NULL;
248 int *segmentmarkerlist = NULL;
249 int nsegs;
[1]250
251 /*Recover input: */
[13363]252 segments = *psegments;
253 segmentmarkerlist = *psegmentmarkerlist;
254 nsegs = *pnsegs;
[1]255
256 /*Reallocate segments: */
[14222]257 segments =xReNew<int>(segments, nsegs*3,(nsegs+nriftsegs)*3);
258 segmentmarkerlist=xReNew<int>(segmentmarkerlist,nsegs,(nsegs+nriftsegs));
[1]259
[8301]260 /*First, update the existing segments to the new nodes :*/
[1]261 for (i=0;i<nriftsegs;i++){
[12080]262 el1=riftsegments[4*i+0];
263 el2=riftsegments[4*i+1];
[1]264 for (j=0;j<nsegs;j++){
[12080]265 if (segments[3*j+2]==(el1+1)){
[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,
[8301]268 *we can only rely on the position (x,y) of the rift nodes to create a segment:*/
[1]269 for (k=0;k<3;k++){
[14222]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])){
[12086]271 *(segments+3*j+0)=*(index+el1*3+k); _assert_(segments[3*j+0]<nods+1);
[1]272 break;
273 }
274 }
275 for (k=0;k<3;k++){
[14222]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])){
[12086]277 *(segments+3*j+1)=*(index+el1*3+k); _assert_(segments[3*j+1]<nods+1);
[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++){
[14222]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])){
[12086]286 *(segments+3*(nsegs+i)+0)=*(index+el2*3+k); _assert_(segments[3*(nsegs+i)+0]<nods+1);
[1]287 break;
288 }
289 }
290 for (k=0;k<3;k++){
[14222]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])){
[12086]292 *(segments+3*(nsegs+i)+1)=*(index+el2*3+k); _assert_(segments[3*(nsegs+i)+1]<nods+1);
[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++){
[14222]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])){
[12086]302 *(segments+3*j+0)=*(index+el2*3+k); _assert_(segments[3*j+0]<nods+1);
[1]303 break;
304 }
305 }
306 for (k=0;k<3;k++){
[14222]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])){
[12086]308 *(segments+3*j+1)=*(index+el2*3+k);_assert_(segments[3*j+1]<nods+1);
[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++){
[14222]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])){
[12086]317 *(segments+3*(nsegs+i)+0)=*(index+el1*3+k);_assert_(segments[3*(nsegs+i)+0]<nods+1);
[1]318 break;
319 }
320 }
321 for (k=0;k<3;k++){
[14222]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])){
[12086]323 *(segments+3*(nsegs+i)+1)=*(index+el1*3+k);_assert_(segments[3*(nsegs+i)+1]<nods+1);
[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;
[13622]336
[1]337 return noerr;
[12080]338}/*}}}*/
[18063]339int FindElement(int A,int B,int* index,int nel){/*{{{*/
[1]340
341 int el=-1;
[14222]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))){
[1]344 el=n;
345 break;
346 }
347 }
348 return el;
[12080]349}/*}}}*/
[18063]350int SplitRiftSegments(int** psegments,int** psegmentmarkerlist, int* pnumsegs, int* pnumrifts,int** priftsnumsegs,int*** priftssegments,int numrifts,int nods,int nel){/*{{{*/
[1]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: */
[14222]359 int *segments = NULL;
360 int *segmentmarkerlist = NULL;
[1]361 int numsegs;
[13622]362
[1]363 /*output: */
[14222]364 int new_numsegs;
365 int *riftsnumsegs = NULL;
366 int **riftssegments = NULL;
367 int *new_segments = NULL;
368 int *new_segmentmarkers = NULL;
[1]369
370 /*intermediary: */
[14222]371 int* riftsegment=NULL;
[1]372
373 /*Recover input arguments: */
[14222]374 segments = *psegments;
375 numsegs = *pnumsegs;
376 segmentmarkerlist = *psegmentmarkerlist;
[1]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;
[14222]385 new_segments=xNew<int>(new_numsegs*3);
386 new_segmentmarkers=xNew<int>(new_numsegs);
[1]387
388 /*Copy new segments info : */
389 counter=0;
390 for (i=0;i<numsegs;i++){
391 if (segmentmarkerlist[i]==1){
[12080]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];
[1]395 new_segmentmarkers[counter]=segmentmarkerlist[i];
396 counter++;
397 }
398 }
399
400 /*Now deal with rift segments: */
[12436]401 riftsnumsegs=xNew<int>(numrifts);
[14222]402 riftssegments=xNew<int*>(numrifts);
[1]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;
[14222]410 riftsegment=xNew<int>(counter*3);
[1]411 /*Copy new segments info :*/
412 counter=0;
413 for (j=0;j<numsegs;j++){
414 if (segmentmarkerlist[j]==(2+i)){
[12086]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);
[1]418 counter++;
419 }
420 }
421 *(riftssegments+i)=riftsegment;
422 }
423
[27232]424 /*Free resources: */
[14222]425 xDelete<int>(segments);
[1]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;
[12080]435}/*}}}*/
[18063]436int PairRiftElements(int** priftsnumpairs,int*** priftspairs,int numrifts,int* riftsnumsegments,int** riftssegments,double* x,double* y){/*{{{*/
[1]437
438 int noerr=1;
439 int i,j,k;
440
441 /*output: */
[14222]442 int *riftsnumpairs = NULL;
443 int **riftspairs = NULL;
[1]444
445 /*intermediary :*/
[14222]446 int numsegs;
447 int* segments=NULL;
448 int* pairs=NULL;
449 int node1,node2,node3,node4;
[1]450
[12436]451 riftsnumpairs=xNew<int>(numrifts);
[14222]452 riftspairs=xNew<int*>(numrifts);
[1]453 for (i=0;i<numrifts;i++){
454 segments=riftssegments[i];
[14222]455 numsegs =riftsnumsegments[i];
[1]456 riftsnumpairs[i]=numsegs;
[14222]457 pairs=xNew<int>(2*numsegs);
[1]458 for (j=0;j<numsegs;j++){
[14222]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;
[1]461 /*Find element facing on other side of rift: */
462 for (k=0;k<numsegs;k++){
463 if (k==j)continue;
[14222]464 node3=segments[3*k+0]-1; node4=segments[3*k+1]-1;
[8301]465 /*We are trying to find 2 elements, where position of node3 == position of node1, and position of node4 == position of node2*/
[16521]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])) ){
[1]468 /*We found the corresponding element: */
[14222]469 pairs[2*j+1]=segments[3*k+2];
[1]470 break;
471 }
472 }
473 }
474 riftspairs[i]=pairs;
475 }
476
477 /*Assign output pointers: */
478 *priftsnumpairs=riftsnumpairs;
479 *priftspairs=riftspairs;
480 return noerr;
[12080]481}/*}}}*/
[18063]482int IsRiftPresent(int* priftflag,int* pnumrifts,int* segmentmarkerlist,int nsegs){/*{{{*/
[1]483
484 int i;
485 int noerr=1;
[13622]486
[1]487 /*output: */
488 int riftflag=0;
489 int numrifts=0;
490
[14222]491 int maxmark=1; //default marker for regular segments
[1]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 }
[14222]501 if(numrifts)riftflag=1;
[1]502
503 /*Assign output pointers:*/
504 *priftflag=riftflag;
505 *pnumrifts=numrifts;
506 return noerr;
[12080]507}/*}}}*/
[18063]508int OrderRifts(int** priftstips,int** riftssegments,int** riftspairs,int numrifts,int* riftsnumsegments,double* x,double* y,int nods,int nels){/*{{{*/
[13622]509
[1]510 int noerr=1;
511 int i,j,k,counter;
512
513 /*intermediary: */
[14222]514 int *riftsegments = NULL;
515 int *riftpairs = NULL;
[1]516 int numsegs;
517
518 /*ordering and copy: */
[14222]519 int *order = NULL;
520 int *riftsegments_copy = NULL;
521 int *riftpairs_copy = NULL;
[1]522
[8301]523 /*node and element manipulation: */
[14222]524 int node1,node2,node3,node4,temp_node,tip1,tip2,node;
525 int el2;
526 int already_ordered=0;
[1]527
528 /*output: */
[14222]529 int* riftstips=NULL;
[1]530
531 /*Allocate byproduct of this routine, riftstips: */
[14222]532 riftstips=xNew<int>(numrifts*2);
[1]533
534 /*Go through all rifts: */
535 for (i=0;i<numrifts;i++){
[14222]536 riftsegments = riftssegments[i];
537 riftpairs = riftspairs[i];
538 numsegs = riftsnumsegments[i];
[13622]539
[1]540 /*Allocate copy of riftsegments and riftpairs,
541 *as well as ordering vector: */
[14222]542 riftsegments_copy=xNew<int>(numsegs*3);
543 riftpairs_copy=xNew<int>(numsegs*2);
[12436]544 order=xNew<int>(numsegs);
[1]545
[8301]546 /*First find the tips, using the pairs. If a pair of elements has one node in common, this node is a rift tip: */
[1]547 tip1=-1;
548 tip2=-1;
549
550 for (j=0;j<numsegs;j++){
551 el2=*(riftpairs+2*j+1);
[14222]552 node1=*(riftsegments+3*j+0);
553 node2=*(riftsegments+3*j+1);
[8301]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: */
[1]556 for (k=0;k<numsegs;k++){
557 if (*(riftsegments+3*k+2)==el2){
[14222]558 node3=*(riftsegments+3*k+0);
559 node4=*(riftsegments+3*k+1);
[1]560 break;
561 }
562 }
[8301]563 /* Make sure node3 faces node1 and node4 faces node2: */
[12086]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])){
[8301]567 /*Swap node3 and node4:*/
568 temp_node=node3;
569 node3=node4;
570 node4=temp_node;
[1]571 }
572
573 /*Figure out if a tip is on this element: */
[8301]574 if (node3==node1){
575 /*node1 is a tip*/
[1]576 if (tip1==-1) {
[8301]577 tip1=node1;
[1]578 continue;
579 }
[8301]580 if ((tip2==-1) && (node1!=tip1)){
581 tip2=node1;
[1]582 break;
583 }
584 }
[13622]585
[8301]586 if (node4==node2){
587 /*node2 is a tip*/
[1]588 if (tip1==-1){
[8301]589 tip1=node2;
[1]590 continue;
591 }
[8301]592 if ((tip2==-1) && (node2!=tip1)){
593 tip2=node2;
[1]594 break;
595 }
596 }
597 }
598
599 /*Record tips in riftstips: */
[14222]600 *(riftstips+2*i+0)=tip1;
601 *(riftstips+2*i+1)=tip2;
[1]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. */
[8301]605 node=tip1;
[1]606 for (counter=0;counter<numsegs;counter++){
607 for (j=0;j<numsegs;j++){
[14222]608 node1=*(riftsegments+3*j+0);
609 node2=*(riftsegments+3*j+1);
[13622]610
[8301]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: */
[1]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;
[8301]622 if(node1==node){
623 node=node2;
[1]624 }
[8301]625 else if(node2==node){
626 node=node1;
[1]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++){
[6412]636 _assert_(order[j]<numsegs);
[1]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 }
[13622]643
[1]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
[12436]652 xDelete<int>(order);
[14222]653 xDelete<int>(riftsegments_copy);
654 xDelete<int>(riftpairs_copy);
[1]655
656 }
657
658 /*Assign output pointer:*/
659 *priftstips=riftstips;
660 return noerr;
[12080]661}/*}}}*/
[18063]662int PenaltyPairs(double*** priftspenaltypairs,int** priftsnumpenaltypairs,int numrifts,int** riftssegments,/*{{{*/
[14222]663 int* riftsnumsegs,int** riftspairs,int* riftstips,double* x,double* y){
[1]664
665 int noerr=1;
666 int i,j,k,k0;
667
[8301]668 double el1,el2,node1,node2,node3,node4;
[9320]669 double temp_node;
[1]670
671 /*output: */
[12436]672 double **riftspenaltypairs = NULL;
673 double *riftpenaltypairs = NULL;
674 int *riftsnumpenaltypairs = NULL;
[1]675
676 /*intermediary: */
677 int numsegs;
[14222]678 int* riftsegments=NULL;
679 int* riftpairs=NULL;
[1]680 int counter;
681 double normal[2];
682 double length;
683 int k1,k2;
684
685 /*Allocate: */
[12436]686 riftspenaltypairs=xNew<double*>(numrifts);
687 riftsnumpenaltypairs=xNew<int>(numrifts);
[1]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: */
[12446]695 if((numsegs/2-1)!=0)riftpenaltypairs=xNewZeroInit<double>((numsegs/2-1)*RIFTPENALTYPAIRSWIDTH);
[13622]696
[1]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);
[8301]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: */
[1]705 k0=-1;
706 for(k=0;k<numsegs;k++){
707 if(*(riftsegments+3*k+2)==el2){
708 k0=k;
709 break;
710 }
711 }
[8301]712 node3=*(riftsegments+3*k0+0);
713 node4=*(riftsegments+3*k0+1);
[1]714
[8301]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;
[1]721 }
[8301]722 /*Ok, we have node1 facing node3, and node2 facing node4. Compute the normal to
[1]723 *this segment, and its length: */
[8301]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));
[1]727
[8301]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,
[1]730 * only once. We'll add the normals and the lengths : */
731
[8301]732 if(node1!=node3){ //exclude tips from loads
[1]733 k1=-1;
734 for(k=0;k<counter;k++){
[8301]735 if( (*(riftpenaltypairs+k*7+0))==node1){
[1]736 k1=k;
737 break;
738 }
739 }
740 if(k1==-1){
[8301]741 *(riftpenaltypairs+counter*7+0)=node1;
742 *(riftpenaltypairs+counter*7+1)=node3;
[1]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 }
[8301]756 if(node2!=node4){
[1]757 k2=-1;
758 for(k=0;k<counter;k++){
[8301]759 if( (*(riftpenaltypairs+k*7+0))==node2){
[1]760 k2=k;
761 break;
762 }
763 }
764 if(k2==-1){
[8301]765 *(riftpenaltypairs+counter*7+0)=node2;
766 *(riftpenaltypairs+counter*7+1)=node4;
[1]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++){
[14222]783 double magnitude=sqrt(pow( double(riftpenaltypairs[j*7+4]),2) + pow( double(riftpenaltypairs[j*7+5]),2) );
[1]784 *(riftpenaltypairs+j*7+4)=*(riftpenaltypairs+j*7+4)/magnitude;
785 *(riftpenaltypairs+j*7+5)=*(riftpenaltypairs+j*7+5)/magnitude;
786 }
[13622]787
[1]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;
[14222]796}/*}}}*/
797int RemoveCornersFromRifts(int** pindex,int* pnel,double** px,double** py,int* pnods,int* segments,int* segmentmarkers,int num_seg){/*{{{*/
[1]798
799 int noerr=1;
800 int i,j,k;
[14222]801 int node1,node2,node3;
[1]802 int el;
803 double pair[2];
804 int pair_count=0;
805 int triple=0;
806
807 /*Recover input: */
[14222]808 int *index = *pindex;
809 int nel = *pnel;
810 double *x = *px;
811 double *y = *py;
812 int nods = *pnods;
[1]813
814 for (i=0;i<num_seg;i++){
[8301]815 node1=*(segments+3*i+0);
816 node2=*(segments+3*i+1);
817 /*Find all elements connected to [node1 node2]: */
[1]818 pair_count=0;
819 for (j=0;j<nel;j++){
[8301]820 if (*(index+3*j+0)==node1){
821 if ((*(index+3*j+1)==node2) || (*(index+3*j+2)==node2)){
[1]822 pair[pair_count]=j;
823 pair_count++;
824 }
825 }
[8301]826 if (*(index+3*j+1)==node1){
827 if ((*(index+3*j+0)==node2) || (*(index+3*j+2)==node2)){
[1]828 pair[pair_count]=j;
829 pair_count++;
830 }
831 }
[8301]832 if (*(index+3*j+2)==node1){
833 if ((*(index+3*j+0)==node2) || (*(index+3*j+1)==node2)){
[1]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,
[8301]840 *figure out if the third node also belongs to a segment: */
[1]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;
[8301]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);
[1]852 }
[8301]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);
[1]856 }
[8301]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);
[1]860 }
[8301]861 /*Ok, we have node3. Does node3 belong to a segment? : */
[1]862 for (k=0;k<num_seg;k++){
[8301]863 if ((node3==*(segments+3*k+0)) || (node3==*(segments+3*k+1))){
[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: */
[13248]870 x=xReNew<double>(x,nods,nods+1);
871 y=xReNew<double>(y,nods,nods+1);
[8301]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;
[14222]874 index=xReNew<int>(index,nel*3,(nel+2*3));
[1]875 /*First, reassign element el: */
[8301]876 *(index+3*el+0)=node1;
877 *(index+3*el+1)=node2;
[1]878 *(index+3*el+2)=nods+1;
879 /*Other two elements: */
[8301]880 *(index+3*nel+0)=node2;
881 *(index+3*nel+1)=node3;
[1]882 *(index+3*nel+2)=nods+1;
883
[8301]884 *(index+3*(nel+1)+0)=node3;
885 *(index+3*(nel+1)+1)=node1;
[1]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++){
[14222]889 if (*(segments+3*k+2)==(el+1)){
[14235]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;
[1]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;
[14222]912}/*}}}*/
Note: See TracBrowser for help on using the repository browser.