Changeset 8301
- Timestamp:
- 05/16/11 15:34:39 (14 years ago)
- Location:
- issm/trunk/src/c/shared
- Files:
-
- 10 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/shared/Exp/DomainOutlineRead.cpp
r6412 r8301 1 1 /*!\file: DomainOutlineRead.cpp 2 * \brief DomainOutlineRead.c: read the gridcoordinates defined in a domain2 * \brief DomainOutlineRead.c: read the vertex coordinates defined in a domain 3 3 * outline from Argus (.exp file). The first contour in the file is for 4 4 * the outside domain outline. The following contours represent holes in … … 11 11 #include "../Exceptions/exceptions.h" 12 12 13 int DomainOutlineRead(int* pnprof,int** pprofn grids,double*** ppprofx,double*** ppprofy,char* domainname){13 int DomainOutlineRead(int* pnprof,int** pprofnvertices,double*** ppprofx,double*** ppprofy,char* domainname){ 14 14 15 15 … … 26 26 /*output: */ 27 27 int nprof; //number of profiles in the domainname file 28 int* profn grids=NULL; //array holding the number of grids for the nprof profiles28 int* profnvertices=NULL; //array holding the number of vertices for the nprof profiles 29 29 double** pprofx=NULL; //array of profiles x coordinates 30 30 double** pprofy=NULL; //array of profiles y coordinates … … 61 61 62 62 /*Allocate and initialize all the profiles: */ 63 profn grids=(int*)xmalloc(nprof*sizeof(int));63 profnvertices=(int*)xmalloc(nprof*sizeof(int)); 64 64 pprofx=(double**)xmalloc(nprof*sizeof(double*)); 65 65 pprofy=(double**)xmalloc(nprof*sizeof(double*)); … … 80 80 fscanf(fid,"%s %s %s %s\n",chardummy,chardummy,chardummy,chardummy); 81 81 82 /*Get number of profile grids: */82 /*Get number of profile vertices: */ 83 83 fscanf(fid,"%u %s\n",&n,chardummy); 84 84 … … 86 86 fscanf(fid,"%s %s %s %s %s\n",chardummy,chardummy,chardummy,chardummy,chardummy); 87 87 88 /*Allocate grids: */88 /*Allocate vertices: */ 89 89 x=(double*)xmalloc(n*sizeof(double)); 90 90 y=(double*)xmalloc(n*sizeof(double)); 91 91 92 92 93 /*Read grids: */93 /*Read vertices: */ 94 94 for (i=0;i<n;i++){ 95 95 fscanf(fid,"%lf %lf\n",x+i,y+i); … … 102 102 103 103 /*Assign pointers: */ 104 profn grids[counter]=n;104 profnvertices[counter]=n; 105 105 pprofx[counter]=x; 106 106 pprofy[counter]=y; … … 120 120 /*Assign output pointers: */ 121 121 *pnprof=nprof; 122 *pprofn grids=profngrids;122 *pprofnvertices=profnvertices; 123 123 *ppprofx=pprofx; 124 124 *ppprofy=pprofy; -
issm/trunk/src/c/shared/Exp/IsInPoly.cpp
r8270 r8301 15 15 16 16 /*IsInPoly {{{1*/ 17 int IsInPoly(Vec in,double* xc,double* yc,int num grids,double* x,double* y,int i0,int i1, int edgevalue){17 int IsInPoly(Vec in,double* xc,double* yc,int numvertices,double* x,double* y,int i0,int i1, int edgevalue){ 18 18 19 19 int i; … … 26 26 27 27 /*Get extrema*/ 28 for (i=1;i<num grids;i++){28 for (i=1;i<numvertices;i++){ 29 29 if(xc[i]<xmin) xmin=xc[i]; 30 30 if(xc[i]>xmax) xmax=xc[i]; … … 33 33 } 34 34 35 /*Go through all grids of the mesh:*/35 /*Go through all vertices of the mesh:*/ 36 36 for (i=i0;i<i1;i++){ 37 37 … … 39 39 VecGetValues(in,1,&i,&value); 40 40 if (value){ 41 /*this gridalready is inside one of the contours, continue*/41 /*this vertex already is inside one of the contours, continue*/ 42 42 continue; 43 43 } 44 44 45 /*pick up grid(x[i],y[i]) and figure out if located inside contour (xc,yc)*/45 /*pick up vertex (x[i],y[i]) and figure out if located inside contour (xc,yc)*/ 46 46 x0=x[i]; y0=y[i]; 47 47 if(x0<xmin || x0>xmax || y0<ymin || y0>ymax){ … … 49 49 } 50 50 else{ 51 value=pnpoly(num grids,xc,yc,x0,y0,edgevalue);51 value=pnpoly(numvertices,xc,yc,x0,y0,edgevalue); 52 52 } 53 53 VecSetValues(in,1,&i,&value,INSERT_VALUES); … … 56 56 }/*}}}*/ 57 57 /*IsOutsidePoly {{{1*/ 58 int IsOutsidePoly(Vec in,double* xc,double* yc,int num grids,double* x,double* y,int i0,int i1, int edgevalue){58 int IsOutsidePoly(Vec in,double* xc,double* yc,int numvertices,double* x,double* y,int i0,int i1, int edgevalue){ 59 59 60 60 int i,j; … … 67 67 68 68 /*Get extrema*/ 69 for (i=1;i<num grids;i++){69 for (i=1;i<numvertices;i++){ 70 70 if(xc[i]<xmin) xmin=xc[i]; 71 71 if(xc[i]>xmax) xmax=xc[i]; … … 74 74 } 75 75 76 /*Go through all grids of the mesh:*/76 /*Go through all vertices of the mesh:*/ 77 77 for (i=i0;i<i1;i++){ 78 78 … … 80 80 VecGetValues(in,1,&i,&value); 81 81 if (value){ 82 /*this gridalready is inside one of the contours, continue*/82 /*this vertex already is inside one of the contours, continue*/ 83 83 continue; 84 84 } 85 85 86 /*pick up grid(x[i],y[i]) and figure out if located inside contour (xc,yc)*/86 /*pick up vertex (x[i],y[i]) and figure out if located inside contour (xc,yc)*/ 87 87 x0=x[i]; y0=y[i]; 88 88 if(x0<xmin || x0>xmax || y0<ymin || y0>ymax){ … … 90 90 } 91 91 else{ 92 value=1-pnpoly(num grids,xc,yc,x0,y0,edgevalue);92 value=1-pnpoly(numvertices,xc,yc,x0,y0,edgevalue); 93 93 } 94 94 VecSetValues(in,1,&i,&value,INSERT_VALUES); -
issm/trunk/src/c/shared/Exp/IsInPolySerial.cpp
r5016 r8301 9 9 10 10 11 int IsInPolySerial(double* in,double* xc,double* yc,int num grids,double* x,double* y,int nods, int edgevalue){11 int IsInPolySerial(double* in,double* xc,double* yc,int numvertices,double* x,double* y,int nods, int edgevalue){ 12 12 13 13 int i,j; 14 14 double x0,y0; 15 15 16 /*Go through all grids of the mesh:*/16 /*Go through all vertices of the mesh:*/ 17 17 for (i=0;i<nods;i++){ 18 18 if (in[i]){ 19 /*this gridalready is inside one of the contours, continue*/19 /*this vertex already is inside one of the contours, continue*/ 20 20 continue; 21 21 } 22 /*pick up grid: */22 /*pick up vertex: */ 23 23 x0=x[i]; 24 24 y0=y[i]; 25 if (pnpoly(num grids,xc,yc,x0,y0,edgevalue)){25 if (pnpoly(numvertices,xc,yc,x0,y0,edgevalue)){ 26 26 in[i]=1; 27 27 } -
issm/trunk/src/c/shared/Exp/exp.h
r8270 r8301 9 9 #include "../../toolkits/toolkits.h" 10 10 11 int IsInPoly(Vec in,double* xc,double* yc,int num grids,double* x,double* y,int i0,int i1, int edgevalue);12 int IsOutsidePoly(Vec in,double* xc,double* yc,int num grids,double* x,double* y,int i0,int i1, int edgevalue);13 int IsInPolySerial(double* in,double* xc,double* yc,int num grids,double* x,double* y,int nods, int edgevalue);14 int DomainOutlineRead(int* pnprof,int** pprofn grids,double*** ppprofx,double*** ppprofy,char* domainname);11 int IsInPoly(Vec in,double* xc,double* yc,int numvertices,double* x,double* y,int i0,int i1, int edgevalue); 12 int IsOutsidePoly(Vec in,double* xc,double* yc,int numvertices,double* x,double* y,int i0,int i1, int edgevalue); 13 int IsInPolySerial(double* in,double* xc,double* yc,int numvertices,double* x,double* y,int nods, int edgevalue); 14 int DomainOutlineRead(int* pnprof,int** pprofnvertices,double*** ppprofx,double*** ppprofy,char* domainname); 15 15 int pnpoly(int npol, double *xp, double *yp, double x, double y, int edgevalue); 16 16 -
issm/trunk/src/c/shared/TriMesh/AssociateSegmentToElement.cpp
r3332 r8301 14 14 double* segments=NULL; 15 15 16 /* gridindices: */16 /*node indices: */ 17 17 double A,B; 18 18 -
issm/trunk/src/c/shared/TriMesh/GridInsideHole.cpp
r3332 r8301 18 18 double yA,yB,yC,yD,yE; 19 19 20 /*Take first and last grids: */20 /*Take first and last vertices: */ 21 21 xA=x[0]; 22 22 yA=y[0]; -
issm/trunk/src/c/shared/TriMesh/OrderSegments.cpp
r3332 r8301 14 14 double* segments=NULL; 15 15 16 /* gridindices: */16 /*vertex indices: */ 17 17 double A,B; 18 18 /*element indices: */ -
issm/trunk/src/c/shared/TriMesh/SplitMeshForRifts.cpp
r5016 r8301 17 17 18 18 int i,j,k,l; 19 int grid;19 int node; 20 20 int el; 21 21 … … 50 50 /*Establish list of segments that belong to a rift: */ 51 51 RiftSegmentsFromSegments(&nriftsegs,&riftsegments,nel,index,nsegs,segments); /*riftsegments of size nriftsegsx4 (4 for first element on segment,second element, 52 first grid and second sgrid)*/52 first node and second snode)*/ 53 53 54 /*Go through all grids of the rift segments, and start splitting the mesh: */55 flags=(int*)xcalloc(nods,sizeof(int)); //to make sure we don't split the same grids twice!54 /*Go through all nodes of the rift segments, and start splitting the mesh: */ 55 flags=(int*)xcalloc(nods,sizeof(int)); //to make sure we don't split the same nodes twice! 56 56 for (i=0;i<nriftsegs;i++){ 57 57 for (j=0;j<2;j++){ 58 58 59 grid=*(riftsegments+4*i+j+2);60 if(flags[ grid-1]){61 /*This gridwas already split, skip:*/59 node=*(riftsegments+4*i+j+2); 60 if(flags[node-1]){ 61 /*This node was already split, skip:*/ 62 62 continue; 63 63 } 64 64 else{ 65 flags[ grid-1]=1;65 flags[node-1]=1; 66 66 } 67 67 68 if(IsGridOnRift(riftsegments,nriftsegs, grid)){68 if(IsGridOnRift(riftsegments,nriftsegs,node)){ 69 69 70 DetermineGridElementListOnOneSideOfRift(&NumGridElementListOnOneSideOfRift,&GridElementListOnOneSideOfRift,i,nriftsegs,riftsegments, grid,index,nel);70 DetermineGridElementListOnOneSideOfRift(&NumGridElementListOnOneSideOfRift,&GridElementListOnOneSideOfRift,i,nriftsegs,riftsegments,node,index,nel); 71 71 72 /*Summary: we have for grid, a list of elements (GridElementListOnOneSideOfRift, of size NumGridElementListOnOneSideOfRift) that all contain grid73 *and that are on the same side of the rift. For all these elements, we clone grid into another grid, and we swap all instances of gridin the triangulation74 *for those elements, to the new grid.*/72 /*Summary: we have for node, a list of elements (GridElementListOnOneSideOfRift, of size NumGridElementListOnOneSideOfRift) that all contain node 73 *and that are on the same side of the rift. For all these elements, we clone node into another node, and we swap all instances of node in the triangulation 74 *for those elements, to the new node.*/ 75 75 76 //augment number of grids76 //augment number of nodes 77 77 nods=nods+1; 78 //create new grid78 //create new node 79 79 x=(double*)xrealloc(x,nods*sizeof(double)); 80 80 y=(double*)xrealloc(y,nods*sizeof(double)); 81 x[nods-1]=x[ grid-1]; //matlab indexing82 y[nods-1]=y[ grid-1]; //matlab indexing81 x[nods-1]=x[node-1]; //matlab indexing 82 y[nods-1]=y[node-1]; //matlab indexing 83 83 84 //change elements owning this grid84 //change elements owning this node 85 85 for (k=0;k<NumGridElementListOnOneSideOfRift;k++){ 86 86 el=GridElementListOnOneSideOfRift[k]; 87 87 for (l=0;l<3;l++){ 88 if (*(index+3*el+l)== grid)*(index+3*el+l)=nods; //again, matlab indexing.88 if (*(index+3*el+l)==node)*(index+3*el+l)=nods; //again, matlab indexing. 89 89 } 90 90 } 91 }// if(IsGridOnRift(riftsegments,nriftsegs, grid))91 }// if(IsGridOnRift(riftsegments,nriftsegs,node)) 92 92 } //for(j=0;j<2;j++) 93 93 } //for (i=0;i<nriftsegs;i++) 94 94 95 /*update segments: they got modified completely by adding new grids.*/95 /*update segments: they got modified completely by adding new nodes.*/ 96 96 UpdateSegments(&segments,&segmentmarkerlist, &nsegs,index,x,y,riftsegments,nriftsegs); 97 97 -
issm/trunk/src/c/shared/TriMesh/TriMeshUtils.cpp
r6748 r8301 11 11 12 12 #define RIFTPENALTYPAIRSWIDTH 8 13 int IsGridOnRift(int* riftsegments, int nriftsegs, int grid){14 15 /*Does this gridbelong to 4 elements, or just 2? If it belongs to 4 elements, it is inside a rift,13 int IsGridOnRift(int* riftsegments, int nriftsegs, int node){ 14 15 /*Does this node belong to 4 elements, or just 2? If it belongs to 4 elements, it is inside a rift, 16 16 *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).*/ 17 17 … … 23 23 for (i=0;i<nriftsegs;i++){ 24 24 for (j=0;j<2;j++){ 25 if ((*(riftsegments+4*i+2+j))== grid) count++;25 if ((*(riftsegments+4*i+2+j))==node) count++; 26 26 } 27 27 } … … 35 35 36 36 37 int GridElementsList(int** pGridElements, int* pNumGridElements,int grid,double * index,int nel){38 39 /*From a grid, recover all the elements that are connected to it: */37 int GridElementsList(int** pGridElements, int* pNumGridElements,int node,double * index,int nel){ 38 39 /*From a node, recover all the elements that are connected to it: */ 40 40 int i,j; 41 41 int noerr=1; … … 48 48 49 49 /*From a mesh with 30 degrees minimum angle, we get 12 possible elements that own 50 * the grid. We start by allocating GridElements with that size, and realloc50 * the node. We start by allocating GridElements with that size, and realloc 51 51 * more if needed.*/ 52 52 … … 57 57 for (i=0;i<nel;i++){ 58 58 for (j=0;j<3;j++){ 59 if ((int)*(index+3*i+j)== grid){59 if ((int)*(index+3*i+j)==node){ 60 60 if (NumGridElements<=(current_size-1)){ 61 61 GridElements[NumGridElements]=i; … … 91 91 92 92 int IsNeighbor(int el1,int el2,double* index){ 93 /*From a triangulation held in index, figure out if elements 1 and 2 have two grids in common: */93 /*From a triangulation held in index, figure out if elements 1 and 2 have two nodes in common: */ 94 94 int i,j; 95 95 int count=0; … … 132 132 int* riftsegments=NULL; 133 133 int* riftsegments_uncompressed=NULL; 134 double element_ grids[3];134 double element_nodes[3]; 135 135 136 136 /*Allocate segmentflags: */ … … 142 142 for (i=0;i<nsegs;i++){ 143 143 el=(int)*(segments+3*i+2)-1; //element found in AssociateSegmentToElements 144 /*Temporarily set grids belonging to the segments to -1 in the triangulation index, and144 /*Temporarily set nodes belonging to the segments to -1 in the triangulation index, and 145 145 *then proceed to find another element that owns the segment. If we don't find it, we know 146 146 *we are dealing with a boundary or hole, otherwise, we are dealing with a rift: */ 147 element_ grids[0]=*(index+3*el+0);148 element_ grids[1]=*(index+3*el+1);149 element_ grids[2]=*(index+3*el+2);147 element_nodes[0]=*(index+3*el+0); 148 element_nodes[1]=*(index+3*el+1); 149 element_nodes[2]=*(index+3*el+2); 150 150 151 151 *(index+3*el+0)=-1; … … 156 156 157 157 /*Restore index: */ 158 *(index+3*el+0)=element_ grids[0];159 *(index+3*el+1)=element_ grids[1];160 *(index+3*el+2)=element_ grids[2];158 *(index+3*el+0)=element_nodes[0]; 159 *(index+3*el+1)=element_nodes[1]; 160 *(index+3*el+2)=element_nodes[2]; 161 161 162 162 if (el2!=-1){ … … 195 195 ******************************************************************************************************************************/ 196 196 197 int DetermineGridElementListOnOneSideOfRift(int* pNumGridElementListOnOneSideOfRift, int** pGridElementListOnOneSideOfRift, int segmentnumber, int nriftsegs, int* riftsegments, int grid,double* index,int nel){197 int DetermineGridElementListOnOneSideOfRift(int* pNumGridElementListOnOneSideOfRift, int** pGridElementListOnOneSideOfRift, int segmentnumber, int nriftsegs, int* riftsegments, int node,double* index,int nel){ 198 198 199 199 int noerr=1; … … 208 208 int* GridElementListOnOneSideOfRift=NULL; 209 209 210 /*Build a list of all the elements connected to this grid: */211 GridElementsList(&GridElements,&NumGridElements, grid,index,nel);210 /*Build a list of all the elements connected to this node: */ 211 GridElementsList(&GridElements,&NumGridElements,node,index,nel); 212 212 213 213 /*Figure out the list of elements that are on the same side of the rift. To do so, we start from one … … 282 282 segmentmarkerlist=(double*)xrealloc(segmentmarkerlist,(nsegs+nriftsegs)*sizeof(double)); 283 283 284 /*First, update the existing segments to the new grids :*/284 /*First, update the existing segments to the new nodes :*/ 285 285 for (i=0;i<nriftsegs;i++){ 286 286 el1=*(riftsegments+4*i+0); … … 290 290 /*segment j is the same as rift segment i.Let's update segments[j][:] using element el1 and the corresponding rift segment. 291 291 *Because riftsegments does not represent a list of rift segments anymore (it got heavily modified in SplitElementsForRifts, 292 *we can only rely on the position (x,y) of the rift grids to create a segment:*/292 *we can only rely on the position (x,y) of the rift nodes to create a segment:*/ 293 293 for (k=0;k<3;k++){ 294 294 if ((x[(int)*(index+el1*3+k)-1]==x[(int)*(segments+3*j+0)-1]) && (y[(int)*(index+el1*3+k)-1]==y[(int)*(segments+3*j+0)-1])){ … … 379 379 IsInPoly 380 380 ******************************************************************************************************************************/ 381 int IsInPoly(double* in,double* xc,double* yc,int num grids,double* x,double* y,int nods){381 int IsInPoly(double* in,double* xc,double* yc,int numnodes,double* x,double* y,int nods){ 382 382 383 383 int i,j; 384 384 double x0,y0; 385 385 386 /*Go through all grids of the mesh:*/386 /*Go through all nodes of the mesh:*/ 387 387 for (i=0;i<nods;i++){ 388 388 if (in[i]){ 389 /*this gridalready is inside one of the contours, continue*/389 /*this node already is inside one of the contours, continue*/ 390 390 continue; 391 391 } 392 /*pick up grid: */392 /*pick up node: */ 393 393 x0=x[i]; 394 394 y0=y[i]; 395 if (pnpoly(num grids,xc,yc,x0,y0)){395 if (pnpoly(numnodes,xc,yc,x0,y0)){ 396 396 in[i]=1; 397 397 } … … 524 524 double* segments=NULL; 525 525 double* pairs=NULL; 526 int grid1,grid2,grid3,grid4;526 int node1,node2,node3,node4; 527 527 528 528 riftsnumpairs=(int*)xmalloc(numrifts*sizeof(int)); … … 535 535 for (j=0;j<numsegs;j++){ 536 536 *(pairs+2*j+0)=*(segments+3*j+2); //retrieve element to which this segment belongs. 537 grid1=(int)*(segments+3*j+0)-1; grid2=(int)*(segments+3*j+1)-1;537 node1=(int)*(segments+3*j+0)-1; node2=(int)*(segments+3*j+1)-1; 538 538 /*Find element facing on other side of rift: */ 539 539 for (k=0;k<numsegs;k++){ 540 540 if (k==j)continue; 541 grid3=(int)*(segments+3*k+0)-1; grid4=(int)*(segments+3*k+1)-1;542 /*We are trying to find 2 elements, where position of grid3 == position of grid1, and position of grid4 == position of grid2*/543 if ( (x[ grid3]==x[grid1]) && (y[grid3]==y[grid1]) && (x[grid4]==x[grid2]) && (y[grid4]==y[grid2])544 || (x[ grid3]==x[grid2]) && (y[grid3]==y[grid2]) && (x[grid4]==x[grid1]) && (y[grid4]==y[grid1]) ){541 node3=(int)*(segments+3*k+0)-1; node4=(int)*(segments+3*k+1)-1; 542 /*We are trying to find 2 elements, where position of node3 == position of node1, and position of node4 == position of node2*/ 543 if ( (x[node3]==x[node1]) && (y[node3]==y[node1]) && (x[node4]==x[node2]) && (y[node4]==y[node2]) 544 || (x[node3]==x[node2]) && (y[node3]==y[node2]) && (x[node4]==x[node1]) && (y[node4]==y[node1]) ){ 545 545 /*We found the corresponding element: */ 546 546 *(pairs+2*j+1)=*(segments+3*k+2); … … 585 585 double* riftsegments=NULL; 586 586 double* riftpairs=NULL; 587 int grid1,grid2,grid3,grid4,temp_grid;587 int node1,node2,node3,node4,temp_node; 588 588 double el1,el2; 589 int newnods; //temporary # gridcounter.589 int newnods; //temporary # node counter. 590 590 double xmin,ymin; 591 591 double* xreal=NULL; 592 592 double* yreal=NULL; 593 int* grids=NULL;594 int* merging grids=NULL;593 int* nodes=NULL; 594 int* mergingnodes=NULL; 595 595 int max_size; 596 596 int redundant; … … 608 608 newnods=nods; 609 609 610 /*Figure out a unique value to flag x and y for gridremoval: */610 /*Figure out a unique value to flag x and y for node removal: */ 611 611 xmin=x[0]; 612 612 ymin=y[0]; … … 618 618 ymin=ymin-dabs(ymin); 619 619 620 /*Initialize two arrays, one for grids that are going to be merged, the other with corresponding grids being merge into: */620 /*Initialize two arrays, one for nodes that are going to be merged, the other with corresponding nodes being merge into: */ 621 621 max_size=0; 622 622 for (i=0;i<numrifts1;i++){ 623 623 max_size+=rifts1numsegs[i]; 624 624 } 625 grids=(int*)xmalloc(max_size*sizeof(int));626 merging grids=(int*)xmalloc(max_size*sizeof(int));627 628 /*Go through the rifts segments, and identify which gridwe are going to merge with its counterpart on the other side629 *of the rift. The way we identify this grid is by looking at the element pairs, and the corresponding grids: */625 nodes=(int*)xmalloc(max_size*sizeof(int)); 626 mergingnodes=(int*)xmalloc(max_size*sizeof(int)); 627 628 /*Go through the rifts segments, and identify which node we are going to merge with its counterpart on the other side 629 *of the rift. The way we identify this node is by looking at the element pairs, and the corresponding nodes: */ 630 630 counter=0; 631 631 for (i=0;i<numrifts1;i++){ … … 635 635 el1=*(riftpairs+2*j+0); 636 636 el2=*(riftpairs+2*j+1); 637 grid1=(int)*(riftsegments+3*j+0);638 grid2=(int)*(riftsegments+3*j+1);639 /*Summary, el1 and el2 are facing one another across the rift. grid1 and grid2 belong to el1 and640 *are located on the rift. Find grid3 and grid4, grids belonging to el2 and located on the rift: */637 node1=(int)*(riftsegments+3*j+0); 638 node2=(int)*(riftsegments+3*j+1); 639 /*Summary, el1 and el2 are facing one another across the rift. node1 and node2 belong to el1 and 640 *are located on the rift. Find node3 and node4, nodes belonging to el2 and located on the rift: */ 641 641 for (k=0;k<rifts1numsegs[i];k++){ 642 642 if (*(riftsegments+3*k+2)==el2){ 643 grid3=(int)*(riftsegments+3*k+0);644 grid4=(int)*(riftsegments+3*k+1);643 node3=(int)*(riftsegments+3*k+0); 644 node4=(int)*(riftsegments+3*k+1); 645 645 break; 646 646 } 647 647 } 648 /* Make sure grid3 faces grid1 and grid4 faces grid2: */649 if ((x[ grid1]==x[grid4]) && (y[grid1]==y[grid4])){650 /*Swap grid3 and grid4:*/651 temp_ grid=grid3;652 grid3=grid4;653 grid4=temp_grid;654 } 655 /* Is any of these two grid pairs on the tip of a rift, in which case, we don't include it in grids or merginggrids: */656 if (( grid1==grid3) || (grid2==grid4)){657 if( grid1!=grid3){658 /*Add grid1 and grid3 to grids and merginggrids if they have not already been added: */648 /* Make sure node3 faces node1 and node4 faces node2: */ 649 if ((x[node1]==x[node4]) && (y[node1]==y[node4])){ 650 /*Swap node3 and node4:*/ 651 temp_node=node3; 652 node3=node4; 653 node4=temp_node; 654 } 655 /* Is any of these two node pairs on the tip of a rift, in which case, we don't include it in nodes or mergingnodes: */ 656 if ((node1==node3) || (node2==node4)){ 657 if(node1!=node3){ 658 /*Add node1 and node3 to nodes and mergingnodes if they have not already been added: */ 659 659 redundant=0; 660 660 for (k=0;k<counter;k++){ 661 if ((merging grids[k]==grid1) || (grids[k]==grid1))redundant=1;661 if ((mergingnodes[k]==node1) || (nodes[k]==node1))redundant=1; 662 662 } 663 663 if(!redundant){ 664 /*Ok, add grid1 to grids, and grid3 to merginggrids: */665 grids[counter]=grid1;666 merging grids[counter]=grid3;664 /*Ok, add node1 to nodes, and node3 to mergingnodes: */ 665 nodes[counter]=node1; 666 mergingnodes[counter]=node3; 667 667 counter++; 668 668 } 669 669 } 670 if( grid2!=grid4){671 /*Add grid2 and grid4 to grids and merginggrids if they have not already been added: */670 if(node2!=node4){ 671 /*Add node2 and node4 to nodes and mergingnodes if they have not already been added: */ 672 672 redundant=0; 673 673 for (k=0;k<counter;k++){ 674 if ((merging grids[k]==grid2) || (grids[k]==grid2))redundant=1;674 if ((mergingnodes[k]==node2) || (nodes[k]==node2))redundant=1; 675 675 } 676 676 if(!redundant){ 677 /*Ok, add grid2 to grids, and grid4 to merginggrids: */678 grids[counter]=grid2;679 merging grids[counter]=grid4;677 /*Ok, add node2 to nodes, and node4 to mergingnodes: */ 678 nodes[counter]=node2; 679 mergingnodes[counter]=node4; 680 680 counter++; 681 681 } … … 683 683 } 684 684 else{ 685 /*Check that grid1 is not already present in the merginggrids: */685 /*Check that node1 is not already present in the mergingnodes: */ 686 686 redundant=0; 687 687 for (k=0;k<counter;k++){ 688 if ((merging grids[k]==grid1) || (grids[k]==grid1))redundant=1;688 if ((mergingnodes[k]==node1) || (nodes[k]==node1))redundant=1; 689 689 } 690 690 if(!redundant){ 691 /*Ok, add grid1 to grids, and grid3 to merginggrids: */692 grids[counter]=grid1;693 merging grids[counter]=grid3;691 /*Ok, add node1 to nodes, and node3 to mergingnodes: */ 692 nodes[counter]=node1; 693 mergingnodes[counter]=node3; 694 694 counter++; 695 695 } 696 /*Check that grid2 is not already present in the merginggrids: */696 /*Check that node2 is not already present in the mergingnodes: */ 697 697 redundant=0; 698 698 for (k=0;k<counter;k++){ 699 if ((merging grids[k]==grid1) || (grids[k]==grid1))redundant=1;699 if ((mergingnodes[k]==node1) || (nodes[k]==node1))redundant=1; 700 700 } 701 701 if(!redundant){ 702 /*Ok, add grid2 to grids, and grid4 to merginggrids: */703 grids[counter]=grid2;704 merging grids[counter]=grid4;702 /*Ok, add node2 to nodes, and node4 to mergingnodes: */ 703 nodes[counter]=node2; 704 mergingnodes[counter]=node4; 705 705 counter++; 706 706 } … … 709 709 } 710 710 711 /*Ok, we have counter pairs of grids (grids and merginggrids): start merging grids in the triangulation: */711 /*Ok, we have counter pairs of nodes (nodes and mergingnodes): start merging nodes in the triangulation: */ 712 712 newnods=nods; 713 713 for (i=0;i<counter;i++){ 714 grid1=grids[i];715 grid3=merginggrids[i];716 /*Merge grid3 into grid1: */717 x[ grid3]=xmin; //flag for later removal from x718 y[ grid3]=ymin; //flag for later removal from y714 node1=nodes[i]; 715 node3=mergingnodes[i]; 716 /*Merge node3 into node1: */ 717 x[node3]=xmin; //flag for later removal from x 718 y[node3]=ymin; //flag for later removal from y 719 719 newnods--; 720 720 for (k=0;k<nel;k++){ 721 if (*(index+3*k+0)== grid3)*(index+3*k+0)=grid1;722 if (*(index+3*k+1)== grid3)*(index+3*k+1)=grid1;723 if (*(index+3*k+2)== grid3)*(index+3*k+2)=grid1;721 if (*(index+3*k+0)==node3)*(index+3*k+0)=node1; 722 if (*(index+3*k+1)==node3)*(index+3*k+1)=node1; 723 if (*(index+3*k+2)==node3)*(index+3*k+2)=node1; 724 724 } 725 725 } … … 805 805 double* riftpairs_copy=NULL; 806 806 807 /* gridand element manipulation: */808 int grid1,grid2,grid3,grid4,temp_grid,tip1,tip2,grid;807 /*node and element manipulation: */ 808 int node1,node2,node3,node4,temp_node,tip1,tip2,node; 809 809 double el1,el2; 810 810 int already_ordered=0; … … 830 830 order=(int*)xmalloc(numsegs*sizeof(int)); 831 831 832 /*First find the tips, using the pairs. If a pair of elements has one grid in common, this gridis a rift tip: */832 /*First find the tips, using the pairs. If a pair of elements has one node in common, this node is a rift tip: */ 833 833 tip1=-1; 834 834 tip2=-1; … … 837 837 el1=*(riftpairs+2*j+0); 838 838 el2=*(riftpairs+2*j+1); 839 grid1=(int)*(riftsegments+3*j+0);840 grid2=(int)*(riftsegments+3*j+1);841 /*Summary, el1 and el2 are facing one another across the rift. grid1 and grid2 belong to el1 and842 *are located on the rift. Find grid3 and grid4, grids belonging to el2 and located on the rift: */839 node1=(int)*(riftsegments+3*j+0); 840 node2=(int)*(riftsegments+3*j+1); 841 /*Summary, el1 and el2 are facing one another across the rift. node1 and node2 belong to el1 and 842 *are located on the rift. Find node3 and node4, nodes belonging to el2 and located on the rift: */ 843 843 for (k=0;k<numsegs;k++){ 844 844 if (*(riftsegments+3*k+2)==el2){ 845 grid3=(int)*(riftsegments+3*k+0);846 grid4=(int)*(riftsegments+3*k+1);845 node3=(int)*(riftsegments+3*k+0); 846 node4=(int)*(riftsegments+3*k+1); 847 847 break; 848 848 } 849 849 } 850 /* Make sure grid3 faces grid1 and grid4 faces grid2: */851 if ((x[ grid1]==x[grid4]) && (y[grid1]==y[grid4])){852 /*Swap grid3 and grid4:*/853 temp_ grid=grid3;854 grid3=grid4;855 grid4=temp_grid;850 /* Make sure node3 faces node1 and node4 faces node2: */ 851 if ((x[node1]==x[node4]) && (y[node1]==y[node4])){ 852 /*Swap node3 and node4:*/ 853 temp_node=node3; 854 node3=node4; 855 node4=temp_node; 856 856 } 857 857 858 858 /*Figure out if a tip is on this element: */ 859 if ( grid3==grid1){860 /* grid1 is a tip*/859 if (node3==node1){ 860 /*node1 is a tip*/ 861 861 if (tip1==-1) { 862 tip1= grid1;862 tip1=node1; 863 863 continue; 864 864 } 865 if ((tip2==-1) && ( grid1!=tip1)){866 tip2= grid1;865 if ((tip2==-1) && (node1!=tip1)){ 866 tip2=node1; 867 867 break; 868 868 } 869 869 } 870 870 871 if ( grid4==grid2){872 /* grid2 is a tip*/871 if (node4==node2){ 872 /*node2 is a tip*/ 873 873 if (tip1==-1){ 874 tip1= grid2;874 tip1=node2; 875 875 continue; 876 876 } 877 if ((tip2==-1) && ( grid2!=tip1)){878 tip2= grid2;877 if ((tip2==-1) && (node2!=tip1)){ 878 tip2=node2; 879 879 break; 880 880 } … … 889 889 /*We have the two tips for this rift. Go from tip1 to tip2, and figure out the order in which segments are sequential. 890 890 *Because two elements are connected to tip1, we chose one first, which defines the direction we are rotating along the rift. */ 891 grid=tip1;891 node=tip1; 892 892 for (counter=0;counter<numsegs;counter++){ 893 893 for (j=0;j<numsegs;j++){ 894 grid1=(int)*(riftsegments+3*j+0);895 grid2=(int)*(riftsegments+3*j+1);894 node1=(int)*(riftsegments+3*j+0); 895 node2=(int)*(riftsegments+3*j+1); 896 896 897 if (( grid1==grid) || (grid2==grid)){898 /*Ok, this segment is connected to grid, plug its index into order, unless we already plugged it before: */897 if ((node1==node) || (node2==node)){ 898 /*Ok, this segment is connected to node, plug its index into order, unless we already plugged it before: */ 899 899 already_ordered=0; 900 900 for (k=0;k<counter;k++){ … … 906 906 if (!already_ordered){ 907 907 order[counter]=j; 908 if( grid1==grid){909 grid=grid2;908 if(node1==node){ 909 node=node2; 910 910 } 911 else if( grid2==grid){912 grid=grid1;911 else if(node2==node){ 912 node=node1; 913 913 } 914 914 break; … … 958 958 int i,j,k,k0; 959 959 960 double el1,el2, grid1,grid2,grid3,grid4;961 double tip1,tip2,temp_ grid;960 double el1,el2,node1,node2,node3,node4; 961 double tip1,tip2,temp_node; 962 962 963 963 /*output: */ … … 994 994 el1=*(riftpairs+2*j+0); 995 995 el2=*(riftpairs+2*j+1); 996 grid1=*(riftsegments+3*j+0);997 grid2=*(riftsegments+3*j+1);998 /*Find segment index to recover grid3 and grid4, facing grid1 and grid2: */996 node1=*(riftsegments+3*j+0); 997 node2=*(riftsegments+3*j+1); 998 /*Find segment index to recover node3 and node4, facing node1 and node2: */ 999 999 k0=-1; 1000 1000 for(k=0;k<numsegs;k++){ … … 1004 1004 } 1005 1005 } 1006 grid3=*(riftsegments+3*k0+0);1007 grid4=*(riftsegments+3*k0+1);1008 1009 /* Make sure grid3 faces grid1 and grid4 faces grid2: */1010 if ((x[(int) grid1-1]==x[(int)grid4-1]) && (y[(int)grid1-1]==y[(int)grid4-1])){1011 /*Swap grid3 and grid4:*/1012 temp_ grid=grid3;1013 grid3=grid4;1014 grid4=temp_grid;1006 node3=*(riftsegments+3*k0+0); 1007 node4=*(riftsegments+3*k0+1); 1008 1009 /* Make sure node3 faces node1 and node4 faces node2: */ 1010 if ((x[(int)node1-1]==x[(int)node4-1]) && (y[(int)node1-1]==y[(int)node4-1])){ 1011 /*Swap node3 and node4:*/ 1012 temp_node=node3; 1013 node3=node4; 1014 node4=temp_node; 1015 1015 } 1016 /*Ok, we have grid1 facing grid3, and grid2 facing grid4. Compute the normal to1016 /*Ok, we have node1 facing node3, and node2 facing node4. Compute the normal to 1017 1017 *this segment, and its length: */ 1018 normal[0]=cos(atan2(x[(int) grid1-1]-x[(int)grid2-1],y[(int)grid2-1]-y[(int)grid1-1]));1019 normal[1]=sin(atan2(x[(int) grid1-1]-x[(int)grid2-1],y[(int)grid2-1]-y[(int)grid1-1]));1020 length=sqrt(pow(x[(int) grid2-1]-x[(int)grid1-1],(double)2)+pow(y[(int)grid2-1]-y[(int)grid1-1],(double)2));1021 1022 /*Be careful here, we want penalty loads on each grid, not on each segment. This means we cannot plug grid1,1023 * grid2, grid3 and grid4 directly into riftpenaltypairs. We need to include grid1, grid2, grid3 and grid4,1018 normal[0]=cos(atan2(x[(int)node1-1]-x[(int)node2-1],y[(int)node2-1]-y[(int)node1-1])); 1019 normal[1]=sin(atan2(x[(int)node1-1]-x[(int)node2-1],y[(int)node2-1]-y[(int)node1-1])); 1020 length=sqrt(pow(x[(int)node2-1]-x[(int)node1-1],(double)2)+pow(y[(int)node2-1]-y[(int)node1-1],(double)2)); 1021 1022 /*Be careful here, we want penalty loads on each node, not on each segment. This means we cannot plug node1, 1023 * node2, node3 and node4 directly into riftpenaltypairs. We need to include node1, node2, node3 and node4, 1024 1024 * only once. We'll add the normals and the lengths : */ 1025 1025 1026 if( grid1!=grid3){ //exclude tips from loads1026 if(node1!=node3){ //exclude tips from loads 1027 1027 k1=-1; 1028 1028 for(k=0;k<counter;k++){ 1029 if( (*(riftpenaltypairs+k*7+0))== grid1){1029 if( (*(riftpenaltypairs+k*7+0))==node1){ 1030 1030 k1=k; 1031 1031 break; … … 1033 1033 } 1034 1034 if(k1==-1){ 1035 *(riftpenaltypairs+counter*7+0)= grid1;1036 *(riftpenaltypairs+counter*7+1)= grid3;1035 *(riftpenaltypairs+counter*7+0)=node1; 1036 *(riftpenaltypairs+counter*7+1)=node3; 1037 1037 *(riftpenaltypairs+counter*7+2)=el1; 1038 1038 *(riftpenaltypairs+counter*7+3)=el2; … … 1048 1048 } 1049 1049 } 1050 if( grid2!=grid4){1050 if(node2!=node4){ 1051 1051 k2=-1; 1052 1052 for(k=0;k<counter;k++){ 1053 if( (*(riftpenaltypairs+k*7+0))== grid2){1053 if( (*(riftpenaltypairs+k*7+0))==node2){ 1054 1054 k2=k; 1055 1055 break; … … 1057 1057 } 1058 1058 if(k2==-1){ 1059 *(riftpenaltypairs+counter*7+0)= grid2;1060 *(riftpenaltypairs+counter*7+1)= grid4;1059 *(riftpenaltypairs+counter*7+0)=node2; 1060 *(riftpenaltypairs+counter*7+1)=node4; 1061 1061 *(riftpenaltypairs+counter*7+2)=el1; 1062 1062 *(riftpenaltypairs+counter*7+3)=el2; … … 1101 1101 int noerr=1; 1102 1102 int i,j,k; 1103 double grid1,grid2,grid3;1103 double node1,node2,node3; 1104 1104 int el; 1105 1105 … … 1123 1123 1124 1124 for (i=0;i<num_seg;i++){ 1125 grid1=*(segments+3*i+0);1126 grid2=*(segments+3*i+1);1127 /*Find all elements connected to [ grid1 grid2]: */1125 node1=*(segments+3*i+0); 1126 node2=*(segments+3*i+1); 1127 /*Find all elements connected to [node1 node2]: */ 1128 1128 pair_count=0; 1129 1129 for (j=0;j<nel;j++){ 1130 if (*(index+3*j+0)== grid1){1131 if ((*(index+3*j+1)== grid2) || (*(index+3*j+2)==grid2)){1130 if (*(index+3*j+0)==node1){ 1131 if ((*(index+3*j+1)==node2) || (*(index+3*j+2)==node2)){ 1132 1132 pair[pair_count]=j; 1133 1133 pair_count++; 1134 1134 } 1135 1135 } 1136 if (*(index+3*j+1)== grid1){1137 if ((*(index+3*j+0)== grid2) || (*(index+3*j+2)==grid2)){1136 if (*(index+3*j+1)==node1){ 1137 if ((*(index+3*j+0)==node2) || (*(index+3*j+2)==node2)){ 1138 1138 pair[pair_count]=j; 1139 1139 pair_count++; 1140 1140 } 1141 1141 } 1142 if (*(index+3*j+2)== grid1){1143 if ((*(index+3*j+0)== grid2) || (*(index+3*j+1)==grid2)){1142 if (*(index+3*j+2)==node1){ 1143 if ((*(index+3*j+0)==node2) || (*(index+3*j+1)==node2)){ 1144 1144 pair[pair_count]=j; 1145 1145 pair_count++; … … 1148 1148 } 1149 1149 /*Ok, we have pair_count elements connected to this segment. For each of these elements, 1150 *figure out if the third gridalso belongs to a segment: */1150 *figure out if the third node also belongs to a segment: */ 1151 1151 if ((pair_count==0) || (pair_count==1)){ //we only select the rift segments, which belong to 2 elements 1152 1152 continue; … … 1156 1156 el=(int)pair[j]; 1157 1157 triple=0; 1158 /*First find grid3: */1159 if (*(index+3*el+0)== grid1){1160 if (*(index+3*el+1)== grid2)grid3=*(index+3*el+2);1161 else grid3=*(index+3*el+1);1162 } 1163 if (*(index+3*el+1)== grid1){1164 if (*(index+3*el+0)== grid2)grid3=*(index+3*el+2);1165 else grid3=*(index+3*el+0);1166 } 1167 if (*(index+3*el+2)== grid1){1168 if (*(index+3*el+0)== grid2)grid3=*(index+3*el+1);1169 else grid3=*(index+3*el+0);1170 } 1171 /*Ok, we have grid3. Does grid3 belong to a segment? : */1158 /*First find node3: */ 1159 if (*(index+3*el+0)==node1){ 1160 if (*(index+3*el+1)==node2)node3=*(index+3*el+2); 1161 else node3=*(index+3*el+1); 1162 } 1163 if (*(index+3*el+1)==node1){ 1164 if (*(index+3*el+0)==node2)node3=*(index+3*el+2); 1165 else node3=*(index+3*el+0); 1166 } 1167 if (*(index+3*el+2)==node1){ 1168 if (*(index+3*el+0)==node2)node3=*(index+3*el+1); 1169 else node3=*(index+3*el+0); 1170 } 1171 /*Ok, we have node3. Does node3 belong to a segment? : */ 1172 1172 for (k=0;k<num_seg;k++){ 1173 if (( grid3==*(segments+3*k+0)) || (grid3==*(segments+3*k+1))){1173 if ((node3==*(segments+3*k+0)) || (node3==*(segments+3*k+1))){ 1174 1174 triple=1; 1175 1175 break; … … 1180 1180 x=(double*)xrealloc(x,(nods+1)*sizeof(double)); 1181 1181 y=(double*)xrealloc(y,(nods+1)*sizeof(double)); 1182 x[nods]=(x[(int) grid1-1]+x[(int)grid2-1]+x[(int)grid3-1])/3;1183 y[nods]=(y[(int) grid1-1]+y[(int)grid2-1]+y[(int)grid3-1])/3;1182 x[nods]=(x[(int)node1-1]+x[(int)node2-1]+x[(int)node3-1])/3; 1183 y[nods]=(y[(int)node1-1]+y[(int)node2-1]+y[(int)node3-1])/3; 1184 1184 1185 1185 index=(double*)xrealloc(index,(nel+2)*3*sizeof(double)); 1186 1186 /*First, reassign element el: */ 1187 *(index+3*el+0)= grid1;1188 *(index+3*el+1)= grid2;1187 *(index+3*el+0)=node1; 1188 *(index+3*el+1)=node2; 1189 1189 *(index+3*el+2)=nods+1; 1190 1190 /*Other two elements: */ 1191 *(index+3*nel+0)= grid2;1192 *(index+3*nel+1)= grid3;1191 *(index+3*nel+0)=node2; 1192 *(index+3*nel+1)=node3; 1193 1193 *(index+3*nel+2)=nods+1; 1194 1194 1195 *(index+3*(nel+1)+0)= grid3;1196 *(index+3*(nel+1)+1)= grid1;1195 *(index+3*(nel+1)+0)=node3; 1196 *(index+3*(nel+1)+1)=node1; 1197 1197 *(index+3*(nel+1)+2)=nods+1; 1198 1198 /*we need to change the segment elements corresponding to el: */ 1199 1199 for (k=0;k<num_seg;k++){ 1200 1200 if (*(segments+3*k+2)==(double)(el+1)){ 1201 if ( ((*(segments+3*k+0)== grid1) && (*(segments+3*k+1)==grid2)) || ((*(segments+3*k+0)==grid2) && (*(segments+3*k+1)==grid1))) *(segments+3*k+2)=(double)(el+1);1202 if ( ((*(segments+3*k+0)== grid2) && (*(segments+3*k+1)==grid3)) || ((*(segments+3*k+0)==grid3) && (*(segments+3*k+1)==grid2))) *(segments+3*k+2)=(double)(nel+1);1203 if ( ((*(segments+3*k+0)== grid3) && (*(segments+3*k+1)==grid1)) || ((*(segments+3*k+0)==grid1) && (*(segments+3*k+1)==grid3))) *(segments+3*k+2)=(double)(nel+2);1201 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)=(double)(el+1); 1202 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)=(double)(nel+1); 1203 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)=(double)(nel+2); 1204 1204 } 1205 1205 } -
issm/trunk/src/c/shared/TriMesh/trimesh.h
r1439 r8301 22 22 int SplitMeshForRifts(int* pnel,double** pindex,int* pnods,double** px,double** py,int* pnsegs,double** psegments,double** psegmentmarkerlist); 23 23 24 int IsGridOnRift(int* riftsegments, int nriftsegs, int grid);25 int GridElementsList(int** pGridElements, int* pNumGridElements,int grid,double * index,int nel);24 int IsGridOnRift(int* riftsegments, int nriftsegs, int node); 25 int GridElementsList(int** pGridElements, int* pNumGridElements,int node,double * index,int nel); 26 26 int IsNeighbor(int el1,int el2,double* index); 27 27 int IsOnRift(int el,int nriftsegs,int* riftsegments); 28 28 int RiftSegmentsFromSegments(int* pnriftsegs, int** priftsegments, int nel, double* index, int nsegs,double* segments); 29 int DetermineGridElementListOnOneSideOfRift(int* pNumGridElementListOnOneSideOfRift, int** pGridElementListOnOneSideOfRift, int segmentnumber, int nriftsegs, int* riftsegments, int grid,double* index,int nel);29 int DetermineGridElementListOnOneSideOfRift(int* pNumGridElementListOnOneSideOfRift, int** pGridElementListOnOneSideOfRift, int segmentnumber, int nriftsegs, int* riftsegments, int node,double* index,int nel); 30 30 int UpdateSegments(double** psegments,double** psegmentmarkerlist, int* pnsegs, double* index, double* x,double* y,int* riftsegments,int nriftsegs); 31 31 int pnpoly(int npol, double *xp, double *yp, double x, double y);
Note:
See TracChangeset
for help on using the changeset viewer.