Homec4science
Diffusion Gear (master) Browse

Gear/src/tessel.c
4b9ba9553668master

Last Change
Raw File

tessel.c
View Options

#ifdef TESSEL
#include <math.h>
#include <string.h>
#include <stdio.h>
#include <allvars.h>
struct Face /* a list of edges (only for 3d tesselation) */
{
};
struct Triangle
{
struct particle_data Pt1[3];
struct particle_data Pt2[3];
struct particle_data Pt3[3];
};
struct TriangleInList
{
int idx; /* index of current triangle (used for checks) */
struct particle_data* P[3]; /* pointers towards the 3 point */
struct TriangleInList* T[3]; /* pointers towards the 3 triangles */
int idxe[3]; /* index of point in the first triangle, opposite to the common edge */
struct Median* Med[3]; /* pointers towards 3 medians */
};
struct Median
{
double a; /* params for the equation of the segemnt */
double b; /* params for the equation of the segemnt */
double c; /* params for the equation of the segemnt */
struct vPoint *vPs; /* pointer towards starting vPoint of the segment */
struct vPoint *vPe; /* pointer towards stoping vPoint of the segment */
struct particle_data *Pa; /* pointer towards point A*/
struct particle_data *Pb; /* pointer towards point B*/
};
/* a voronoi point */
struct vPoint
{
double Pos[3];
int next;
};
/*
LOCAL VARIABLES
*/
#define MAXNUMTRIANGLES 1000000
#define MAXVPOINTS 5000000
#define MAXMEDIANS 5000000 /* obsolete */
int nT=0,numTinStack=0; /* number of triangles in the list */
struct TriangleInList Triangles[MAXNUMTRIANGLES]; /* list of triangles */
struct TriangleInList *TStack[MAXNUMTRIANGLES]; /* index of triangles to check */
struct Median MediansList[3*MAXNUMTRIANGLES][3];
int nvPoints=0; /* number of Voronoi Points */
int nMedians=0; /* number of Medians (obsolete) */
struct vPoint vPoints[MAXVPOINTS];
struct Median Medians[MAXMEDIANS];
double tesselDomainRadius,tesselDomainCenter[3];
struct particle_data Pe[3]; /* edges */
int td=0;
char tdname[50];
void setup_searching_radius()
{
int i;
for(i=0;i<NumPart;i++)
P[i].rSearch=0.1/10000;
}
void lines_intersections(double a0, double b0, double c0, double a1, double b1, double c1, double *x, double *y)
{
*x = (c1*b0 - c0*b1)/(a0*b1 - a1*b0);
*y = (c1*a0 - c0*a1)/(a1*b0 - a0*b1);
}
/*!
*/
struct Triangle TriangleInList2Triangle(struct TriangleInList Tl)
{
struct Triangle T;
T.Pt1->Pos[0] = Tl.P[0]->Pos[0];
T.Pt1->Pos[1] = Tl.P[0]->Pos[1];
T.Pt2->Pos[0] = Tl.P[1]->Pos[0];
T.Pt2->Pos[1] = Tl.P[1]->Pos[1];
T.Pt3->Pos[0] = Tl.P[2]->Pos[0];
T.Pt3->Pos[1] = Tl.P[2]->Pos[1];
return T;
}
/*! For a set of three points, construct a triangle
*/
struct Triangle MakeTriangleFromPoints(struct particle_data Pt1,struct particle_data Pt2,struct particle_data Pt3)
{
struct Triangle T;
T.Pt1->Pos[0] = Pt1.Pos[0];
T.Pt1->Pos[1] = Pt1.Pos[1];
T.Pt2->Pos[0] = Pt2.Pos[0];
T.Pt2->Pos[1] = Pt2.Pos[1];
T.Pt3->Pos[0] = Pt3.Pos[0];
T.Pt3->Pos[1] = Pt3.Pos[1];
return T;
}
/*! For a set of three points, this function computes the 3 medians.
*/
void TriangleMedians(struct particle_data Pt1,struct particle_data Pt2,struct particle_data Pt3,struct particle_data *Pmm1,struct particle_data *Pmm2,struct particle_data *Pmm3,struct particle_data *Pme1,struct particle_data *Pme2,struct particle_data *Pme3)
{
double ma1,mb1,mc1;
double ma2,mb2,mc2;
double ma3,mb3,mc3;
/* median 0-1 */
ma1 = 2*(Pt2.Pos[0] - Pt1.Pos[0]);
mb1 = 2*(Pt2.Pos[1] - Pt1.Pos[1]);
mc1 = (Pt1.Pos[0]*Pt1.Pos[0]) - (Pt2.Pos[0]*Pt2.Pos[0]) + (Pt1.Pos[1]*Pt1.Pos[1]) - (Pt2.Pos[1]*Pt2.Pos[1]);
/* median 1-2 */
ma2 = 2*(Pt3.Pos[0] - Pt2.Pos[0]);
mb2 = 2*(Pt3.Pos[1] - Pt2.Pos[1]);
mc2 = (Pt2.Pos[0]*Pt2.Pos[0]) - (Pt3.Pos[0]*Pt3.Pos[0]) + (Pt2.Pos[1]*Pt2.Pos[1]) - (Pt3.Pos[1]*Pt3.Pos[1]);
/* median 2-0 */
ma3 = 2*(Pt1.Pos[0] - Pt3.Pos[0]);
mb3 = 2*(Pt1.Pos[1] - Pt3.Pos[1]);
mc3 = (Pt3.Pos[0]*Pt3.Pos[0]) - (Pt1.Pos[0]*Pt1.Pos[0]) + (Pt3.Pos[1]*Pt3.Pos[1]) - (Pt1.Pos[1]*Pt1.Pos[1]);
/* intersection m0-1 -- m1-2 */
Pmm1->Pos[0] = (mc2*mb1 - mc1*mb2)/(ma1*mb2 - ma2*mb1);
Pmm1->Pos[1] = (mc2*ma1 - mc1*ma2)/(ma2*mb1 - ma1*mb2);
/* intersection m1-2 -- m2-0 */
Pmm2->Pos[0] = (mc2*mb1 - mc1*mb2)/(ma1*mb2 - ma2*mb1);
Pmm2->Pos[1] = (mc2*ma1 - mc1*ma2)/(ma2*mb1 - ma1*mb2);
/* intersection m2-0 -- m0-1 */
Pmm3->Pos[0] = (mc2*mb1 - mc1*mb2)/(ma1*mb2 - ma2*mb1);
Pmm3->Pos[1] = (mc2*ma1 - mc1*ma2)/(ma2*mb1 - ma1*mb2);
/* intersection m1-2 -- e1-2 */
Pme1->Pos[0] = 0.5*(Pt1.Pos[0] + Pt2.Pos[0]);
Pme1->Pos[1] = 0.5*(Pt1.Pos[1] + Pt2.Pos[1]);
/* intersection m2-3 -- e3-1 */
Pme2->Pos[0] = 0.5*(Pt2.Pos[0] + Pt3.Pos[0]);
Pme2->Pos[1] = 0.5*(Pt2.Pos[1] + Pt3.Pos[1]);
/* intersection m3-1 -- e1-2 */
Pme3->Pos[0] = 0.5*(Pt3.Pos[0] + Pt1.Pos[0]);
Pme3->Pos[1] = 0.5*(Pt3.Pos[1] + Pt1.Pos[1]);
}
/*! For a set of three points, this function computes their cirum-circle.
* Its radius is return, while the center is return using pointers.
*/
double CircumCircleProperties(struct particle_data Pt1,struct particle_data Pt2,struct particle_data Pt3, double *xc, double *yc)
{
double r;
double x21,x32,y21,y32;
double x12mx22,y12my22,x22mx32,y22my32;
double c1,c2;
x21 = Pt2.Pos[0]-Pt1.Pos[0];
x32 = Pt3.Pos[0]-Pt2.Pos[0];
y21 = Pt2.Pos[1]-Pt1.Pos[1];
y32 = Pt3.Pos[1]-Pt2.Pos[1];
x12mx22 = (Pt1.Pos[0]*Pt1.Pos[0])-(Pt2.Pos[0]*Pt2.Pos[0]);
y12my22 = (Pt1.Pos[1]*Pt1.Pos[1])-(Pt2.Pos[1]*Pt2.Pos[1]);
x22mx32 = (Pt2.Pos[0]*Pt2.Pos[0])-(Pt3.Pos[0]*Pt3.Pos[0]);
y22my32 = (Pt2.Pos[1]*Pt2.Pos[1])-(Pt3.Pos[1]*Pt3.Pos[1]);
c1 = x12mx22 + y12my22;
c2 = x22mx32 + y22my32;
*xc = (y32*c1 - y21*c2)/2.0/( x32*y21 - x21*y32 );
*yc = (x32*c1 - x21*c2)/2.0/( x21*y32 - x32*y21 );
r = sqrt( (Pt1.Pos[0]-*xc)*(Pt1.Pos[0]-*xc) + (Pt1.Pos[1]-*yc)*(Pt1.Pos[1]-*yc) ) ;
return r;
}
/*! For a given triangle T, the routine tells if the point P4
is in the circum circle of the triangle or not.
*/
int InCircumCircle(struct Triangle T,struct particle_data Pt4)
{
double a,b,c;
double d,e,f;
double g,h,i;
double det;
/*
a = T.Pt1->Pos[0] - Pt4.Pos[0];
b = T.Pt1->Pos[1] - Pt4.Pos[1];
c = (T.Pt1->Pos[0]*T.Pt1->Pos[0] - Pt4.Pos[0]*Pt4.Pos[0]) + (T.Pt1->Pos[1]*T.Pt1->Pos[1] - Pt4.Pos[1]*Pt4.Pos[1]);
d = T.Pt2->Pos[0] - Pt4.Pos[0];
e = T.Pt2->Pos[1] - Pt4.Pos[1];
f = (T.Pt2->Pos[0]*T.Pt2->Pos[0] - Pt4.Pos[0]*Pt4.Pos[0]) + (T.Pt2->Pos[1]*T.Pt2->Pos[1] - Pt4.Pos[1]*Pt4.Pos[1]);
g = T.Pt3->Pos[0] - Pt4.Pos[0];
h = T.Pt3->Pos[1] - Pt4.Pos[1];
i = (T.Pt3->Pos[0]*T.Pt3->Pos[0] - Pt4.Pos[0]*Pt4.Pos[0]) + (T.Pt3->Pos[1]*T.Pt3->Pos[1] - Pt4.Pos[1]*Pt4.Pos[1]);
*/
/*
Volker Formula
*/
a = T.Pt2->Pos[0] - T.Pt1->Pos[0];
b = T.Pt2->Pos[1] - T.Pt1->Pos[1];
c = a*a + b*b;
d = T.Pt3->Pos[0] - T.Pt1->Pos[0];
e = T.Pt3->Pos[1] - T.Pt1->Pos[1];
f = d*d + e*e;
g = Pt4.Pos[0] - T.Pt1->Pos[0];
h = Pt4.Pos[1] - T.Pt1->Pos[1];
i = g*g + h*h;
det = a*e*i - a*f*h - b*d*i + b*f*g + c*d*h - c*e*g;
if (det<0)
return 1; /* inside */
else
return 0; /* outside */
}
/*! For a given triangle T, the routine tells if the point P4
lie inside the triangle or not.
*/
int InTriangle(struct Triangle T,struct particle_data Pt4)
{
double c1,c2,c3;
/* here, we use the cross product */
c1 = (T.Pt2->Pos[0]-T.Pt1->Pos[0])*(Pt4.Pos[1]-T.Pt1->Pos[1]) - (T.Pt2->Pos[1]-T.Pt1->Pos[1])*(Pt4.Pos[0]-T.Pt1->Pos[0]);
c2 = (T.Pt3->Pos[0]-T.Pt2->Pos[0])*(Pt4.Pos[1]-T.Pt2->Pos[1]) - (T.Pt3->Pos[1]-T.Pt2->Pos[1])*(Pt4.Pos[0]-T.Pt2->Pos[0]);
c3 = (T.Pt1->Pos[0]-T.Pt3->Pos[0])*(Pt4.Pos[1]-T.Pt3->Pos[1]) - (T.Pt1->Pos[1]-T.Pt3->Pos[1])*(Pt4.Pos[0]-T.Pt3->Pos[0]);
if ( (c1>0) && (c2>0) && (c3>0) ) /* inside */
return 1;
else
return 0;
}
int InTriangleOrOutside(struct Triangle T,struct particle_data Pt4)
{
double c1,c2,c3;
c1 = (T.Pt2->Pos[0]-T.Pt1->Pos[0])*(Pt4.Pos[1]-T.Pt1->Pos[1]) - (T.Pt2->Pos[1]-T.Pt1->Pos[1])*(Pt4.Pos[0]-T.Pt1->Pos[0]);
if (c1<0)
return 2; /* to triangle T[2] */
c2 = (T.Pt3->Pos[0]-T.Pt2->Pos[0])*(Pt4.Pos[1]-T.Pt2->Pos[1]) - (T.Pt3->Pos[1]-T.Pt2->Pos[1])*(Pt4.Pos[0]-T.Pt2->Pos[0]);
if (c2<0)
return 0; /* to triangle T[1] */
c3 = (T.Pt1->Pos[0]-T.Pt3->Pos[0])*(Pt4.Pos[1]-T.Pt3->Pos[1]) - (T.Pt1->Pos[1]-T.Pt3->Pos[1])*(Pt4.Pos[0]-T.Pt3->Pos[0]);
if (c3<0)
return 1; /* to triangle T[0] */
return -1; /* the point is inside */
}
/*! For a given triangle, orient it positively.
*/
struct Triangle OrientTriangle(struct Triangle T)
{
double a,b,c,d;
double det;
struct particle_data Ptsto;
a = T.Pt2->Pos[0] - T.Pt1->Pos[0];
b = T.Pt2->Pos[1] - T.Pt1->Pos[1];
c = T.Pt3->Pos[0] - T.Pt1->Pos[0];
d = T.Pt3->Pos[1] - T.Pt1->Pos[1];
det = (a*d) - (b*c);
if (det<0)
{
Ptsto.Pos[0] = T.Pt1->Pos[0];
Ptsto.Pos[1] = T.Pt1->Pos[1];
T.Pt1->Pos[0] = T.Pt3->Pos[0];
T.Pt1->Pos[1] = T.Pt3->Pos[1];
T.Pt3->Pos[0] = Ptsto.Pos[0];
T.Pt3->Pos[1] = Ptsto.Pos[1];
T = OrientTriangle(T);
}
return T;
}
/*! For a given triangle, orient it positively.
*/
struct TriangleInList OrientTriangleInList(struct TriangleInList T)
{
double a,b,c,d;
double det;
struct particle_data Ptsto;
a = T.P[1]->Pos[0] - T.P[0]->Pos[0];
b = T.P[1]->Pos[1] - T.P[0]->Pos[1];
c = T.P[2]->Pos[0] - T.P[0]->Pos[0];
d = T.P[2]->Pos[1] - T.P[0]->Pos[1];
det = (a*d) - (b*c);
if (det<0)
{
Ptsto.Pos[0] = T.P[0]->Pos[0];
Ptsto.Pos[1] = T.P[0]->Pos[1];
T.P[0]->Pos[0] = T.P[2]->Pos[0];
T.P[0]->Pos[1] = T.P[2]->Pos[1];
T.P[2]->Pos[0] = Ptsto.Pos[0];
T.P[2]->Pos[1] = Ptsto.Pos[1];
T = OrientTriangleInList(T);
}
return T;
}
/*! This function computes the extension of the domain.
* It computes:
* len : max-min
* tesselDomainCenter : min + 0.5*len
* tesselDomainRadius = len*1.5;
*/
void FindTesselExtent()
{
int i,j;
double xmin[3], xmax[3],len;
/* determine local extension */
for(j = 0; j < 3; j++)
{
xmin[j] = MAX_REAL_NUMBER;
xmax[j] = -MAX_REAL_NUMBER;
}
for(i = 0; i < NumPart; i++)
{
for(j = 0; j < 3; j++)
{
if(xmin[j] > P[i].Pos[j])
xmin[j] = P[i].Pos[j];
if(xmax[j] < P[i].Pos[j])
xmax[j] = P[i].Pos[j];
}
}
len = 0;
for(j = 0; j < 3; j++)
{
if(xmax[j] - xmin[j] > len)
len = xmax[j] - xmin[j];
}
for(j = 0; j < 3; j++)
tesselDomainCenter[j] = xmin[j] + len/2.;
tesselDomainRadius = len*2;
printf("tesselDomainRadius = %g\n",tesselDomainRadius);
printf("tesselDomainCenter = (%g %g)\n",tesselDomainCenter[0],tesselDomainCenter[1]);
}
int FindSegmentInTriangle(struct TriangleInList *T,double v,struct particle_data P[3])
{
double v0,v1,v2;
double x0,x1,x2;
double y0,y1,y2;
double f;
double x,y;
int iP;
/* if the triangle as an edge point, do nothing */
if ( (T->P[0]==&Pe[0]) || (T->P[1]==&Pe[0]) || (T->P[2]==&Pe[0]) )
return 0;
/* if the triangle as an edge point, do nothing */
if ( (T->P[0]==&Pe[1]) || (T->P[1]==&Pe[1]) || (T->P[2]==&Pe[1]) )
return 0;
/* if the triangle as an edge point, do nothing */
if ( (T->P[0]==&Pe[2]) || (T->P[1]==&Pe[2]) || (T->P[2]==&Pe[2]) )
return 0;
iP = 0;
v0 = T->P[0]->Mass;
v1 = T->P[1]->Mass;
v2 = T->P[2]->Mass;
//printf("Triangle %d : %g %g %g\n",T->idx,v0,v1,v2);
/* we could also use the sign v-v0 * v-v1 ??? */
if (( ((v>v0)&&(v<v1)) || ((v>v1)&&(v<v0)) )&& (v0 != v1)) /* in 0-1 */
{
x0 = T->P[0]->Pos[0];
y0 = T->P[0]->Pos[1];
x1 = T->P[1]->Pos[0];
y1 = T->P[1]->Pos[1];
f = (v-v0)/(v1-v0);
P[iP].Pos[0] = f*(x1-x0) + x0;
P[iP].Pos[1] = f*(y1-y0) + y0;
iP++;
}
if (( ((v>v1)&&(v<v2)) || ((v>v2)&&(v<v1)) )&& (v1 != v2)) /* in 1-2 */
{
x0 = T->P[1]->Pos[0];
y0 = T->P[1]->Pos[1];
x1 = T->P[2]->Pos[0];
y1 = T->P[2]->Pos[1];
f = (v-v1)/(v2-v1);
P[iP].Pos[0] = f*(x1-x0) + x0;
P[iP].Pos[1] = f*(y1-y0) + y0;
iP++;
}
if (( ((v>v2)&&(v<v0)) || ((v>v0)&&(v<v2)) )&& (v2 != v0)) /* in 2-0 */
{
x0 = T->P[2]->Pos[0];
y0 = T->P[2]->Pos[1];
x1 = T->P[0]->Pos[0];
y1 = T->P[0]->Pos[1];
f = (v-v2)/(v0-v2);
P[iP].Pos[0] = f*(x1-x0) + x0;
P[iP].Pos[1] = f*(y1-y0) + y0;
iP++;
}
return iP;
}
void CheckTriangles(void)
{
int iT;
struct TriangleInList *T,*Te;
for (iT=0;iT<nT;iT++)
{
T = &Triangles[iT];
Te = T->T[0];
if (Te!=NULL)
{
if ((Te->T[0]!=NULL)&&(Te->T[0] == T))
{
}
else
if ((Te->T[1]!=NULL)&&(Te->T[1] == T))
{
}
else
if ((Te->T[2]!=NULL)&&(Te->T[2] == T))
{
}
else
{
printf("Triangle %d does not point towards %d, while T->T2=%d\n",Te->idx,T->idx,T->T[0]->idx);
exit(-1);
}
}
Te = T->T[1];
if (Te!=NULL)
{
if ((Te->T[0]!=NULL)&&(Te->T[0] == T))
{
}
else
if ((Te->T[1]!=NULL)&&(Te->T[1] == T))
{
}
else
if ((Te->T[2]!=NULL)&&(Te->T[2] == T))
{
}
else
{
printf("Triangle %d does not point towards %d, while T->T2=%d\n",Te->idx,T->idx,T->T[1]->idx);
exit(-1);
}
}
Te = T->T[2];
if (Te!=NULL)
{
if ((Te->T[0]!=NULL)&&(Te->T[0] == T))
{
}
else
if ((Te->T[1]!=NULL)&&(Te->T[1] == T))
{
}
else
if ((Te->T[2]!=NULL)&&(Te->T[2] == T))
{
}
else
{
printf("Triangle %d does not point towards %d, while T->T2=%d\n",Te->idx,T->idx,T->T[2]->idx);
exit(-1);
}
}
}
}
void CheckPoints(int n)
{
int i,p,iT;
for (i=0;i<n;i++)
{
/* find p */
iT = P[i].iT;
p=-1;
if (Triangles[iT].P[0]==&P[i])
p=0;
else
if (Triangles[iT].P[1]==&P[i])
p=1;
else
if (Triangles[iT].P[2]==&P[i])
p=2;
// printf(" ---> i=%d (id=%d) iT=%d \n",i,P[i].ID,iT);
// printf(" P.ID=%d%09d\n", (int) (Triangles[iT].P[0]->ID / 1000000000), (int) (Triangles[iT].P[0]->ID % 1000000000));
// printf(" P.ID=%d%09d\n", (int) (Triangles[iT].P[1]->ID / 1000000000), (int) (Triangles[iT].P[1]->ID % 1000000000));
// printf(" P.ID=%d%09d\n", (int) (Triangles[iT].P[2]->ID / 1000000000), (int) (Triangles[iT].P[2]->ID % 1000000000));
if (p==-1)
endrun(-123555);
//printf("%d ok\n",i);
}
}
/*! Flip two triangles.
Te = T.T[i]
*/
void FlipTriangle(int i,struct TriangleInList *T,struct TriangleInList *Te,struct TriangleInList *T1,struct TriangleInList *T2)
{
struct TriangleInList Ts1,Ts2;
int i0,i1,i2;
int j0,j1,j2;
int j;
Ts1 = *T; /* save the content of the pointed triangle */
Ts2 = *Te; /* save the content of the pointed triangle */
j = T->idxe[i]; /* index of point opposite to i */
i0= i;
i1= (i+1) % 3;
i2= (i+2) % 3;
j0= j;
j1= (j+1) % 3;
j2= (j+2) % 3;
/* triangle 1 */
T1->P[0] = Ts1.P[i0];
T1->P[1] = Ts1.P[i1];
T1->P[2] = Ts2.P[j0];
T1->T[0] = Ts2.T[j1];
T1->T[1] = T2;
T1->T[2] = Ts1.T[i2];
T1->idxe[0] = Ts2.idxe[j1];
T1->idxe[1] = 1;
T1->idxe[2] = Ts1.idxe[i2];
/* triangle 2 */
T2->P[0] = Ts2.P[j0];
T2->P[1] = Ts2.P[j1];
T2->P[2] = Ts1.P[i0];
/* restore links with points */
// if (Ts1.P[i0]->iT==Ts1.idx)
// Ts1.P[i0]->iT=T1->idx;
// if (Ts1.P[i1]->iT==Ts1.idx)
// Ts1.P[i1]->iT=T1->idx;
if (Ts1.P[i2]->iT==Ts1.idx)
Ts1.P[i2]->iT=T2->idx;
/* if (Ts2.P[j0]->iT==Ts2.idx)
Ts2.P[j0]->iT=T2->idx;
if (Ts2.P[j1]->iT==Ts2.idx)
Ts2.P[j1]->iT=T2->idx; */
if (Ts2.P[j2]->iT==Ts2.idx)
Ts2.P[j2]->iT=T1->idx;
T2->T[0] = Ts1.T[i1];
T2->T[1] = T1;
T2->T[2] = Ts2.T[j2];
T2->idxe[0] = Ts1.idxe[i1];
T2->idxe[1] = 1;
T2->idxe[2] = Ts2.idxe[j2];
/* restore links with adjacents triangles */
if (Ts1.T[i1]!=NULL)
{
Ts1.T[i1]->T[ Ts1.idxe[i1] ] = T2;
Ts1.T[i1]->idxe[ Ts1.idxe[i1] ] = 0;
}
if (Ts1.T[i2] !=NULL)
{
Ts1.T[i2]->T[ Ts1.idxe[i2] ] = T1;
Ts1.T[i2]->idxe[ Ts1.idxe[i2] ] = 2;
}
if (Ts2.T[j1] !=NULL)
{
Ts2.T[j1]->T[ Ts2.idxe[j1] ] = T1;
Ts2.T[j1]->idxe[ Ts2.idxe[j1] ] = 0;
}
if (Ts2.T[j2] !=NULL)
{
Ts2.T[j2]->T[ Ts2.idxe[j2] ] = T2;
Ts2.T[j2]->idxe[ Ts2.idxe[j2] ] = 2;
}
}
void DoTrianglesInStack(void)
{
struct TriangleInList *T,*Te,*T1,*T2,*Tee;
struct TriangleInList Ts1,Ts2;
struct particle_data Pt;
int istack;
int idx1,idx2;
int i;
istack=0;
while(numTinStack>0)
{
int insphere=0;
T = TStack[istack];
//printf(" DoInStack T=%d (istack=%d, numTinStack=%d)\n",T->idx,istack,numTinStack);
/* find the opposite point of the 3 adjacent triangles */
/*******************/
/* triangle 1 */
/*******************/
i = 0;
Te = T->T[i];
if (Te!=NULL)
{
/* index of opposite point */
Pt = *Te->P[T->idxe[i]];
insphere = InCircumCircle(TriangleInList2Triangle(*T),Pt);
if (insphere)
{
//printf("insphere (1)... %g %g %g in T=%d\n",Pt.Pos[0],Pt.Pos[1],Pt.Pos[2],T->idx);
/* index of the new triangles */
idx1 = T->idx;
idx2 = Te->idx;
T1 = &Triangles[idx1];
T2 = &Triangles[idx2];
FlipTriangle(i,T,Te,T1,T2);
/* add triangles in stack */
if (numTinStack+1>MAXNUMTRIANGLES)
{
printf("\nNo more memory ! i=0\n");
printf("numTinStack+1=%d > MAXNUMTRIANGLES=%d\n",numTinStack+1,MAXNUMTRIANGLES);
printf("You should increase MAXNUMTRIANGLES\n\n");
exit(-1);
}
TStack[istack ] = T1;
TStack[istack+numTinStack] = T2;
numTinStack++;
continue;
}
}
/*******************/
/* triangle 2 */
/*******************/
i = 1;
Te = T->T[i];
if (Te!=NULL)
{
/* index of opposite point */
Pt = *Te->P[T->idxe[i]];
insphere = InCircumCircle(TriangleInList2Triangle(*T),Pt);
if (insphere)
{
//printf("insphere (2)... %g %g %g in T=%d\n",Pt.Pos[0],Pt.Pos[1],Pt.Pos[2],T->idx);
/* index of the new triangles */
idx1 = T->idx;
idx2 = Te->idx;
T1 = &Triangles[idx1];
T2 = &Triangles[idx2];
FlipTriangle(i,T,Te,T1,T2);
/* add triangles in stack */
if (numTinStack+1>MAXNUMTRIANGLES)
{
printf("\nNo more memory ! i=1\n");
printf("numTinStack+1=%d > MAXNUMTRIANGLES=%d\n",numTinStack+1,MAXNUMTRIANGLES);
printf("You should increase MAXNUMTRIANGLES\n\n");
exit(-1);
}
TStack[istack ] = T1;
TStack[istack+numTinStack] = T2;
numTinStack++;
continue;
}
}
/*******************/
/* triangle 3 */
/*******************/
i = 2;
Te = T->T[i];
if (Te!=NULL)
{
/* index of opposite point */
Pt = *Te->P[T->idxe[i]];
insphere = InCircumCircle(TriangleInList2Triangle(*T),Pt);
if (insphere)
{
//printf("insphere (3)... %g %g %g in T=%d\n",Pt.Pos[0],Pt.Pos[1],Pt.Pos[2],T->idx);
/* index of the new triangles */
idx1 = T->idx;
idx2 = Te->idx;
T1 = &Triangles[idx1];
T2 = &Triangles[idx2];
FlipTriangle(i,T,Te,T1,T2);
/* add triangles in stack */
if (numTinStack+1>MAXNUMTRIANGLES)
{
printf("\nNo more memory ! i=2\n");
printf("numTinStack+1=%d > MAXNUMTRIANGLES=%d\n",numTinStack+1,MAXNUMTRIANGLES);
printf("You should increase MAXNUMTRIANGLES\n\n");
exit(-1);
}
TStack[istack ] = T1;
TStack[istack+numTinStack] = T2;
numTinStack++;
continue;
}
}
numTinStack--;
istack++;
//printf("one triangle less...(istack=%d numTinStack=%d)\n",istack,numTinStack);
}
}
void Check(void)
{
int iT;
printf("===========================\n");
for(iT=0;iT<nT;iT++)
{
printf("* T %d\n",Triangles[iT].idx);
printf("pt1 %g %g %g\n",Triangles[iT].P[0]->Pos[0],Triangles[iT].P[0]->Pos[1],Triangles[iT].P[0]->Pos[2]);
printf("pt2 %g %g %g\n",Triangles[iT].P[1]->Pos[0],Triangles[iT].P[1]->Pos[1],Triangles[iT].P[1]->Pos[2]);
printf("pt3 %g %g %g\n",Triangles[iT].P[2]->Pos[0],Triangles[iT].P[2]->Pos[1],Triangles[iT].P[2]->Pos[2]);
if (Triangles[iT].T[0]!=NULL)
printf("T1 %d\n",Triangles[iT].T[0]->idx);
else
printf("T1 x\n");
if (Triangles[iT].T[1]!=NULL)
printf("T2 %d\n",Triangles[iT].T[1]->idx);
else
printf("T2 x\n");
if (Triangles[iT].T[2]!=NULL)
printf("T3 %d\n",Triangles[iT].T[2]->idx);
else
printf("T3 x\n");
}
printf("===========================\n");
}
/*! Split a triangle in 3, using the point P inside it.
Update the global list.
*/
void SplitTriangle(struct TriangleInList *pT,struct particle_data *Pt)
{
struct TriangleInList T,*T0,*T1,*T2,*Te;
int idx,idx0,idx1,idx2;
T = *pT; /* save the content of the pointed triangle */
idx = T.idx;
/* index of the new triangles */
idx0 = idx;
idx1 = nT;
idx2 = nT+1;
/* increment counter */
nT=nT+2;
/* check memory */
if (nT>MAXNUMTRIANGLES)
{
printf("\nNo more memory !\n");
printf("nT=%d > MAXNUMTRIANGLES=%d\n",nT,MAXNUMTRIANGLES);
printf("You should increase MAXNUMTRIANGLES\n\n");
exit(-1);
}
/* create pointers towards the triangles */
T0 = &Triangles[idx0];
T1 = &Triangles[idx1];
T2 = &Triangles[idx2];
/* first */
T0->idx = idx0;
T0->P[0] = T.P[0];
T0->P[1] = T.P[1];
T0->P[2] = Pt;
/* second */
T1->idx = idx1;
T1->P[0] = T.P[1];
T1->P[1] = T.P[2];
T1->P[2] = Pt;
/* third */
T2->idx = idx2;
T2->P[0] = T.P[2];
T2->P[1] = T.P[0];
T2->P[2] = Pt;
/* link points with triangles */
/* the new point Pt belongs to the
3 new triangle, link it to the first */
Pt->iT = idx0;
/* if P[0] points towards the old Triangle, link it with T0 */
if (T.P[0]->iT==T.idx)
{
T.P[0]->iT=T0->idx;
}
/* if P[1] points towards the old Triangle, link it with T1 */
if (T.P[1]->iT==T.idx)
{
T.P[1]->iT=T1->idx;
}
/* if P[2] points towards the old Triangle, link it with T2 */
if (T.P[2]->iT==T.idx)
{
T.P[2]->iT=T2->idx;
}
/* add adjacents */
T0->T[0] = T1;
T0->T[1] = T2;
T0->T[2] = T.T[2];
T1->T[0] = T2;
T1->T[1] = T0;
T1->T[2] = T.T[0];
T2->T[0] = T0;
T2->T[1] = T1;
T2->T[2] = T.T[1];
/* add ext point */
T0->idxe[0] = 1;
T0->idxe[1] = 0;
T0->idxe[2] = T.idxe[2];
T1->idxe[0] = 1;
T1->idxe[1] = 0;
T1->idxe[2] = T.idxe[0];
T2->idxe[0] = 1;
T2->idxe[1] = 0;
T2->idxe[2] = T.idxe[1];
/* restore links with adgacents triangles */
Te = T0->T[2];
if (Te!=NULL)
{
Te->T[ T0->idxe[2]] = T0;
Te->idxe[T0->idxe[2]] = 2;
}
Te = T1->T[2];
if (Te!=NULL)
{
Te->T[ T1->idxe[2]] = T1;
Te->idxe[T1->idxe[2]] = 2;
}
Te = T2->T[2];
if (Te!=NULL)
{
Te->T[ T2->idxe[2]] = T2;
Te->idxe[T2->idxe[2]] = 2;
}
/* add the new triangles in the stack */
TStack[numTinStack] = T0;
numTinStack++;
TStack[numTinStack] = T1;
numTinStack++;
TStack[numTinStack] = T2;
numTinStack++;
//printf("--> add in stack %d %d %d\n",T0->idx,T1->idx,T2->idx);
}
int FindTriangle(struct particle_data *Pt)
{
int iT;
/* find triangle containing the point */
for(iT=0;iT<nT;iT++) /* loop over all triangles */
{
if (InTriangle(TriangleInList2Triangle( Triangles[iT] ),*Pt))
break;
}
return iT;
}
int NewFindTriangle(struct particle_data *Pt)
{
int iT;
struct TriangleInList *T;
int e;
iT = 0; /* star with first triangle */
T = &Triangles[iT];
while (1)
{
/* test position of the point relative to the triangle */
e = InTriangleOrOutside(TriangleInList2Triangle( *T ),*Pt);
//printf("T=%d e=%d Te=%d\n",T->idx,e,T->T[e]->idx);
if (e==-1) /* the point is inside */
break;
T = T->T[e];
if (T==NULL)
{
printf("point lie outside the limits.\n");
printf("Pt x=%g y=%g z=%g\n",Pt->Pos[0],Pt->Pos[1],Pt->Pos[2]);
endrun(-1);
}
}
//printf("done with find triangle (T=%d)\n",T->idx);
return T->idx;
}
/*! Add a new point in the tesselation
*/
void AddPoint(struct particle_data *Pt)
{
int iT;
/* find the triangle that contains the point P */
//iT= FindTriangle(Pt);
//printf("Find triangle...\n");
//sprintf(tdname,"%s%s_%04d.dat","./","tria",td++);
//dump_triangles(tdname);
iT= NewFindTriangle(Pt);
/* create the new triangles */
//printf("Split triangle...\n");
SplitTriangle(&Triangles[iT],Pt);
//sprintf(tdname,"%s%s_%04d.dat","./","trib",td++);
//dump_triangles(tdname);
/* test the new triangles and divide and modify if necessary */
//printf("Do triangles in stack... nT=%d\n",nT);
DoTrianglesInStack();
/* check */
//CheckTriangles();
}
void AddGhostPoints(int istart,int nadd)
{
int i;
if (ThisTask==0)
printf("%d new ghost points\n",nadd);
/* loop over all ghostpoints */
for (i=istart;i<istart+nadd;i++)
{
AddPoint(&P[i]);
//printf("check after %d\n",i);
#ifdef CHECK_IN_TESSEL
CheckTriangles();
#endif
}
//if (ThisTask==0)
// printf("%d triangles\n",nT);
#ifdef CHECK_IN_TESSEL
CheckTriangles();
#endif
}
/*! Add gost points to the tesselation
*/
void ExtendWithGhostPoints()
{
size_t bytes;
/* all particles try to find ngbs particles that lie in another proc */
if (ThisTask==0)
printf("ExtendWithGhostPoints...\n");
All.MaxgPart=All.MaxPart;
NumgPart=0;
/* if(!(gP = malloc(bytes = All.MaxgPart * sizeof(struct ghost_particle_data))))
{
printf("failed to allocate memory for `gP' (%g MB).\n", bytes / (1024.0 * 1024.0));
endrun(1);
} */
ghost();
/* free(gP);*/
if (ThisTask==0)
printf("ExtendWithGhostPoints done.\n");
}
/*! Compute all medians properties (a,b,c)
* For each triangle, for each edge, the function computes the
* median properties which is stored in MediansList
*/
void ComputeMediansProperties()
{
int iT;
/* loop over all triangles */
for(iT=0;iT<nT;iT++)
{
struct particle_data Pt0,Pt1,Pt2;
Pt0.Pos[0] = Triangles[iT].P[0]->Pos[0];
Pt0.Pos[1] = Triangles[iT].P[0]->Pos[1];
Pt1.Pos[0] = Triangles[iT].P[1]->Pos[0];
Pt1.Pos[1] = Triangles[iT].P[1]->Pos[1];
Pt2.Pos[0] = Triangles[iT].P[2]->Pos[0];
Pt2.Pos[1] = Triangles[iT].P[2]->Pos[1];
/* median 0-1 */
MediansList[iT][2].a = 2*(Pt1.Pos[0] - Pt0.Pos[0]);
MediansList[iT][2].b = 2*(Pt1.Pos[1] - Pt0.Pos[1]);
MediansList[iT][2].c = (Pt0.Pos[0]*Pt0.Pos[0]) - (Pt1.Pos[0]*Pt1.Pos[0]) + (Pt0.Pos[1]*Pt0.Pos[1]) - (Pt1.Pos[1]*Pt1.Pos[1]);
/* median 1-2 */
MediansList[iT][0].a = 2*(Pt2.Pos[0] - Pt1.Pos[0]);
MediansList[iT][0].b = 2*(Pt2.Pos[1] - Pt1.Pos[1]);
MediansList[iT][0].c = (Pt1.Pos[0]*Pt1.Pos[0]) - (Pt2.Pos[0]*Pt2.Pos[0]) + (Pt1.Pos[1]*Pt1.Pos[1]) - (Pt2.Pos[1]*Pt2.Pos[1]);
/* median 2-0 */
MediansList[iT][1].a = 2*(Pt0.Pos[0] - Pt2.Pos[0]);
MediansList[iT][1].b = 2*(Pt0.Pos[1] - Pt2.Pos[1]);
MediansList[iT][1].c = (Pt2.Pos[0]*Pt2.Pos[0]) - (Pt0.Pos[0]*Pt0.Pos[0]) + (Pt2.Pos[1]*Pt2.Pos[1]) - (Pt0.Pos[1]*Pt0.Pos[1]);
/* link The triangle with the MediansList */
Triangles[iT].Med[0] = &MediansList[iT][0]; /* median 1-2 */
Triangles[iT].Med[1] = &MediansList[iT][1]; /* median 2-0 */
Triangles[iT].Med[2] = &MediansList[iT][2]; /* median 0-1 */
}
}
/*! Compute the intersetions of medians around a point of index p (index of the point in the triangle T)
*
*/
void ComputeMediansAroundPoint(struct TriangleInList *Tstart,int iPstart)
{
/*
Tstart : pointer to first triangle
iPstart : index of master point relative to triangle Tstart
if p = 0:
T1 = T0->T[iTn]; pn=1
if p = 1:
T1 = T0->T[iTn]; pn=2
if p = 0:
T1 = T0->T[iTn]; pn=3
iTn = (p+1) % 3;
*/
double x,y;
struct TriangleInList *T0,*T1;
int iP0,iP1;
int iT1;
struct particle_data *initialPoint;
int iM0,iM1;
int next_Median=-1; /* index towards next median */
int number_of_vPoints=0;
int number_of_Medians=0;
int next;
int initialvPoint;
T0 = Tstart;
iP0 = iPstart;
initialPoint = T0->P[iP0];
initialPoint->ivPoint=nvPoints;
initialvPoint=nvPoints;
//printf("\n--> rotating around T=%d p=%d\n",T0->idx,iP0);
/* rotate around the point */
while (1)
{
/* next triangle */
iT1= (iP0+1) % 3;
T1 = T0->T[iT1];
if (T1==NULL)
{
//printf("reach an edge\n");
T0->P[iP0]->IsDone=2;
//printf("%g %g\n",T0->P[iP0]->Pos[0],T0->P[iP0]->Pos[1]);
return;
}
//printf(" next triangle = %d\n",T1->idx);
/* index of point in the triangle */
iP1 = T0->idxe[iT1]; /* index of point opposite to iTn */
iP1 = (iP1+1) % 3; /* next index of point opposite to iTn */
//printf(" initial point=%g %g current point =%g %g iP1=%d\n",initialPoint->Pos[0],initialPoint->Pos[1],T1->P[iP1]->Pos[0],T1->P[iP1]->Pos[1],iP1);
/* check */
if (initialPoint!=T1->P[iP1])
{
printf(" problem : initial point=%g %g current point =%g %g iP1=%d\n",initialPoint->Pos[0],initialPoint->Pos[1],T1->P[iP1]->Pos[0],T1->P[iP1]->Pos[1],iP1);
exit(-1);
}
/* compute the intersection of the two medians */
iM0 = (iP0+1) % 3;
iM1 = (iP1+1) % 3;
lines_intersections(T0->Med[iM0]->a,T0->Med[iM0]->b,T0->Med[iM0]->c,T1->Med[iM1]->a,T1->Med[iM1]->b,T1->Med[iM1]->c,&x,&y);
/* create a new vPoint and put it to the vPoints list */
vPoints[nvPoints].Pos[0] = x;
vPoints[nvPoints].Pos[1] = y;
vPoints[nvPoints].Pos[2] = 0.5;
vPoints[nvPoints].next = nvPoints+1;
/* end point for T0 */
T0->Med[iM0]->vPe = &vPoints[nvPoints];
/* here, we could add the point to T0->Med[(iM0+2) % 3] as Ps */
/* start point for T0 */
T1->Med[iM1]->vPs = &vPoints[nvPoints];
/* here, we could add the point to T1->Med[(iM1+2) % 3] as Ps */
/* increment vPoints */
nvPoints++;
number_of_vPoints++;
if (nvPoints==MAXVPOINTS)
{
printf("you need to increase ");
endrun(-55655);
}
if (T1==Tstart) /* end of loop */
{
//printf(" end of loop\n");
initialPoint->nvPoints = number_of_vPoints;
if (nvPoints-2<initialvPoint)
{
printf("this is bad... please check.\n");
endrun(-72354);
}
initialPoint->ivPoint = nvPoints-2; /* the current point points towards the previous-1 one */
vPoints[nvPoints-1].next = initialvPoint; /* close the loop : the last points towards the first */
/* create the median list */
/* first vPoint */
// next = initialPoint->ivPoint;
//
// for (j = 0; j < initialPoint->nvPoints; j++)
// {
// Medians[nMedians].vPe = vPoints[prev];
// Medians[nMedians].vPs = vPoints[next];
// next = vPoints[next].next;
// }
break;
}
T0 = T1;
iP0 = iP1;
}
}
/*! Compute all medians intersections and define Ps and Pe
* For each triangle, compute the medians
*/
// void oldComputeMediansIntersections()
// {
// int i,p,iT;
//
// for (i=0;i<NumPart+NumgPart;i++)
// P[i].IsDone = 0;
//
// /* loop over all triangles */
// for(iT=0;iT<nT;iT++)
// {
// /* loop over points in triangle */
// for(p=0;p<3;p++)
// {
// if (!(Triangles[iT].P[p]->IsDone))
// {
// //printf("in Triangle T %d do point %d\n",iT,p);
// Triangles[iT].P[p]->IsDone = 1;
// ComputeMediansAroundPoint(&Triangles[iT],p);
// }
// }
// }
// }
void ComputeMediansIntersections()
{
int i,p,iT;
for (i=0;i<NumPart+NumgPart;i++)
{
/* find p */
iT = P[i].iT;
p=-1;
if (Triangles[iT].P[0]==&P[i])
p=0;
else
if (Triangles[iT].P[1]==&P[i])
p=1;
else
if (Triangles[iT].P[2]==&P[i])
p=2;
if (p==-1)
endrun(-12354);
ComputeMediansAroundPoint(&Triangles[P[i].iT],p);
}
// /* do the same for edges */
// for (i=0;i<3;i++)
// {
// /* find p */
// iT = Pe[i].iT;
// p=-1;
// if (Triangles[iT].P[0]->ID==Pe[i].ID)
// p=0;
// else
// if (Triangles[iT].P[1]->ID==Pe[i].ID)
// p=1;
// else
// if (Triangles[iT].P[2]->ID==Pe[i].ID)
// p=2;
//
//
// if (p==-1)
// endrun(-12354);
//
// ComputeMediansAroundPoint(&Triangles[Pe[i].iT],p);
// }
}
/*! Compute the density for all particles
*/
int ComputeDensity()
{
int i,j;
int next; /* next voronoi point */
int np;
double x0,y0,x1,y1;
double area;
/* do the normal points first */
for(i=0;i<NumPart;i++)
{
next = P[i].ivPoint;
np = P[i].nvPoints;
x0 = 0; /* this ensure the first loop to give 0 */
y0 = 0; /* this ensure the first loop to give 0 */
area = 0;
for (j = 0; j < np; j++)
{
x1 = vPoints[next].Pos[0];
y1 = vPoints[next].Pos[1];
area = area + (x0*y1 - x1*y0);
x0 = x1;
y0 = y1;
next = vPoints[next].next;
}
/* connect the last with the first */
next = P[i].ivPoint;
x1 = vPoints[next].Pos[0];
y1 = vPoints[next].Pos[1];
area = area + (x0*y1 - x1*y0);
/* */
area = 0.5*fabs(area);
P[i].Volume=area;
P[i].Density=P[i].Mass/area;
P[i].Pressure = (P[i].Entropy) * pow(P[i].Density, GAMMA); /* SphP[i].Entropy this is bad, as we use SphP */
}
/* do the ghost points now */
for(i=NumPart;i<NumPart+NumgPart;i++)
{
P[i].Volume = P[P[i].iPref].Volume;
P[i].Density= P[P[i].iPref].Density;
P[i].Pressure = P[P[i].iPref].Pressure;
}
return 0;
}
/*! Compute the density for all particles
*/
int CheckCompletenessForThisPoint(int i)
{
/* loop over the triangles which contains the point */
/* first triangle */
int iT,iTstart;
int p,iP0,iP1,iPnext0;
struct TriangleInList *T,*Tnext,*Tstart;
double rc,xc,yc,rcmax;
iTstart = P[i].iT;
Tstart=&Triangles[iTstart];
/* find index in the triangle */
p=-1;
if (Tstart->P[0]==&P[i])
p=0;
else
if (Tstart->P[1]==&P[i])
p=1;
else
if (Tstart->P[2]==&P[i])
p=2;
if (p==-1)
endrun(-312354);
iT = iTstart; /* index of triangle */
iP0 = p; /* index of the ref point "0" in the triangle */
T = Tstart;
/* visit all triangles around the point */
rcmax=0;
while (1)
{
/* compute circumcircle properties of the current tirangle T */
rc = CircumCircleProperties(*T->P[0],*T->P[1],*T->P[2], &xc, &yc);
if (rc>rcmax)
rcmax=rc;
iP1 = (iP0+1) % 3 ; /* point 1 in the current triangle */
Tnext = T->T[iP1]; /* triangle opposit to iP1 */
if (Tnext==NULL) /* we reach an edge */
{
printf("mmmh... why are we here ?");
endrun(-31235445);
return -1; /* need to increase the size for this point */
}
/* now, find the index of the ref. point in the next triangle */
iPnext0 = (T->idxe[iP1] + 1) % 3 ;
/* check */
if (&P[i]!=Tnext->P[iPnext0])
//if (P[i].ID!=Tnext->P[iPnext0]->ID)
{
printf(" problem : initial point %d (id=%d x=%g y=%g z=%g) next-> (id=%d x=%g y=%g z=%g)\n",i,P[i].ID,P[i].Pos[0],P[i].Pos[1],P[i].Pos[2],Tnext->P[iPnext0]->ID,Tnext->P[i]->Pos[0],Tnext->P[i]->Pos[1],Tnext->P[i]->Pos[2]);
endrun(-3123544);
}
/* end of loop */
if (Tnext==Tstart)
break;
iP0 = iPnext0;
T = Tnext;
}
if (P[i].rSearch>2*rcmax)
return 0; /* ok, no need to */
else
{
/* need to increase hi */
//printf("need to increase hi for i=%d %g %g \n",i,P[i].rSearch,2*rcmax);
P[i].rSearch *= 1.5;
return 1;
}
}
/*! Construct the Delaunay tesselation
*/
int ConstructDelaunay()
{
int i,j;
nT=0;
if (ThisTask==0)
printf("Start ConstructDelaunay...\n");
/* find domain extent */
FindTesselExtent();
/* set edges Pe, the 3 points are in an equilateral triangle around all particles */
for (j=0;j<3;j++)
{
Pe[j].Pos[0] = tesselDomainCenter[0] + tesselDomainRadius * cos(2./3.*PI*j)*2;
Pe[j].Pos[1] = tesselDomainCenter[1] + tesselDomainRadius * sin(2./3.*PI*j)*2;
Pe[j].Pos[2] = 0;
Pe[j].Mass = 0;
}
/* Triangle list */
Triangles[0].idx = 0;
Triangles[0].P[0] = &Pe[0];
Triangles[0].P[1] = &Pe[1];
Triangles[0].P[2] = &Pe[2];
Triangles[0].T[0] = NULL;
Triangles[0].T[1] = NULL;
Triangles[0].T[2] = NULL;
Triangles[0].idxe[0] = -1;
Triangles[0].idxe[1] = 1;
Triangles[0].idxe[2] = -1;
nT++;
OrientTriangleInList(Triangles[0]);
if (ThisTask==0)
printf("Start local AddPoint in Delaunay.\n");
/* loop over all points */
for (i=0;i<NumPart;i++)
{
AddPoint(&P[i]);
#ifdef CHECK_IN_TESSEL
CheckPoints(i-1);
#endif
}
//if (ThisTask==0)
// printf("%d triangles\n",nT);
/* add ghost points */
#ifdef CHECK_IN_TESSEL
CheckTriangles();
#endif
//printf("check after AddPoint local is ok\n");
ExtendWithGhostPoints();
/* check */
#ifdef CHECK_IN_TESSEL
CheckTriangles();
CheckPoints(NumPart+NumgPart);
#endif
if (ThisTask==0)
printf("ConstructDelaunay done.\n");
}
/*! Compute the Voronoi tesselation from the Delonay one
*/
int ComputeVoronoi()
{
if (ThisTask==0)
printf("Start ComputeVoronoi...\n");
/* reset nvPoints */
nvPoints=0;
ComputeMediansProperties();
ComputeMediansIntersections();
ComputeDensity();
if (ThisTask==0)
printf("ComputeVoronoi done.\n");
}
/*! Compute accelerations
*/
void tessel_convert_energy_to_entropy()
{
#ifndef ISOTHERM_EQS
int i;
double a3=1;
if(All.ComovingIntegrationOn)
a3 = (All.Time * All.Time * All.Time);
else
a3 = 1;
if(header.flag_entropy_instead_u == 0)
for(i = 0; i < N_gas; i++)
{
P[i].Entropy = GAMMA_MINUS1 * SphP[i].Entropy / pow(P[i].Density / a3, GAMMA_MINUS1);
P[i].Pressure = (P[i].Entropy) * pow(P[i].Density, GAMMA);
}
#endif
}
/*! Compute accelerations
*/
void tessel_compute_accelerations()
{
int i;
int np;
/* first triangle */
int iT,iTstart;
int p,iP0,iP1,iPnext0;
struct TriangleInList *T,*Tnext,*Tstart;
double acc[3];
double dx, dy, dz;
double adx,ady,adz;
double Rij,Aij;
double eij[3],mpij[3],cmij[3],cij[3];
struct particle_data *Pi,*Pj;
struct vPoint *vP1,*vP2;
int nextvp1,nextvp2;
double PjpPi,PjmPi;
for(i=0;i<N_gas;i++)
{
/* some init */
acc[0] = acc[1] = acc[2] = 0;
Pi=&P[i];
np=0;
iTstart = P[i].iT;
Tstart=&Triangles[iTstart];
/* find index in the triangle */
p=-1;
if (Tstart->P[0]==&P[i])
p=0;
else
if (Tstart->P[1]==&P[i])
p=1;
else
if (Tstart->P[2]==&P[i])
p=2;
if (p==-1)
endrun(-712354);
iT = iTstart; /* index of triangle */
iP0 = p; /* index of the ref point "0" in the triangle */
T = Tstart;
np = 0;
nextvp1 = Pi->ivPoint;
/* visit all triangles around the point */
while (1)
{
iP1 = (iP0+1) % 3 ; /* point 1 in the current triangle */
Tnext = T->T[iP1]; /* triangle opposit to iP1 */
if (Tnext==NULL) /* we reach an edge */
{
printf("mmmh... why are we here ?");
endrun(-71235445);
return -1; /* need to increase the size for this point */
}
/* we have a ngb point */
np++;
Pj = T->P[iP1];
nextvp2 = vPoints[nextvp1].next;
vP1 = &vPoints[nextvp1];
vP2 = &vPoints[nextvp2];
dx = Pi->Pos[0] - Pj->Pos[0];
dy = Pi->Pos[1] - Pj->Pos[1];
dz = Pi->Pos[2] - Pj->Pos[2];
Rij = sqrt(dx * dx + dy * dy + dz * dz);
eij[0] = (Pj->Pos[0]-Pi->Pos[0])/Rij ;
eij[1] = (Pj->Pos[1]-Pi->Pos[1])/Rij ;
eij[2] = (Pj->Pos[2]-Pi->Pos[2])/Rij ;
/* mid point between i and j */
mpij[0] = 0.5*(Pj->Pos[0]+Pi->Pos[0]) ;
mpij[1] = 0.5*(Pj->Pos[1]+Pi->Pos[1]) ;
mpij[2] = 0.5*(Pj->Pos[2]+Pi->Pos[2]) ;
#ifdef TWODIMS
/* center of mass of the face */
cmij[0] = 0.5*(vP1->Pos[0]+vP2->Pos[0]);
cmij[1] = 0.5*(vP1->Pos[1]+vP2->Pos[1]);
cmij[2] = 0.5*(vP1->Pos[2]+vP2->Pos[2]);
/* area of the face */
adx = vP2->Pos[0]-vP1->Pos[0];
ady = vP2->Pos[1]-vP1->Pos[1];
adz = vP2->Pos[2]-vP1->Pos[2];
Aij = sqrt(adx * adx + ady * ady + adz * adz);
#else
need to compute mass center here
#endif
/* cij (mid point between ij, towards mass center of the face) cmij - mpij */
cij[0] = cmij[0]-mpij[0];
cij[1] = cmij[1]-mpij[1];
cij[2] = cmij[2]-mpij[2];
/* compute pressure */
PjpPi = Pj->Pressure + Pi->Pressure;
PjmPi = Pj->Pressure - Pi->Pressure;
/* now, the accelection */
acc[0] += -Aij*( 0.5*PjpPi*eij[0] + PjmPi*cij[0]/Rij);
acc[1] += -Aij*( 0.5*PjpPi*eij[1] + PjmPi*cij[1]/Rij);
acc[2] += -Aij*( 0.5*PjpPi*eij[2] + PjmPi*cij[2]/Rij);
/* now, find the index of the ref. point in the next triangle */
iPnext0 = (T->idxe[iP1] + 1) % 3 ;
/* check */
if (&P[i]!=Tnext->P[iPnext0])
//if (P[i].ID!=Tnext->P[iPnext0]->ID)
{
printf(" problem : initial point %d (id=%d x=%g y=%g z=%g) next-> (id=%d x=%g y=%g z=%g)\n",i,P[i].ID,P[i].Pos[0],P[i].Pos[1],P[i].Pos[2],Tnext->P[iPnext0]->ID,Tnext->P[i]->Pos[0],Tnext->P[i]->Pos[1],Tnext->P[i]->Pos[2]);
endrun(-7123544);
}
/* end of loop */
if (Tnext==Tstart)
break;
iP0 = iPnext0;
T = Tnext;
nextvp1 = nextvp2;
}
P[i].tesselAccel[0] = acc[0]/P[i].Mass;
P[i].tesselAccel[1] = acc[1]/P[i].Mass;
P[i].tesselAccel[2] = acc[2]/P[i].Mass;
}
}
void tessel_kick(float dt_kick)
{
int i,j;
for(i = 0; i < NumPart; i++)
{
for(j = 0; j < 3; j++)
P[i].Vel[j] += P[i].tesselAccel[j] * dt_kick;
}
}
void tessel_drift(float dt_drift)
{
int i,j;
for(i = 0; i < NumPart; i++)
{
for(j = 0; j < 3; j++)
P[i].Pos[j] += P[i].Vel[j] * dt_drift;
}
}
double tessel_get_timestep()
{
double ax, ay, az, ac, csnd, hsml;
int i;
double dt,dtmin;
dtmin = 1e32;
for(i = 0; i < NumPart; i++)
{
//ax = P[i].tesselAccel[0];
//ay = P[i].tesselAccel[1];
//az = P[i].tesselAccel[2];
//ac = sqrt(ax * ax + ay * ay + az * az);
csnd = sqrt(GAMMA * P[i].Pressure / P[i].Density);
#ifdef TWODIMS
hsml = sqrt(P[i].Volume/PI);
#else
endrun("-787899");
#endif
dt = 2 * All.CourantFac * hsml / csnd;
if (dt<dtmin)
dtmin = dt;
}
return dtmin;
}
/************************************************************/
/* */
/* IO ROUTINES */
/* */
/************************************************************/
#define SKIP {my_fwrite(&blksize,sizeof(int),1,fd);}
void dump_triangles(char *fname)
{
int iT;
FILE *fd = 0;
int blksize;
double f;
int idx;
//sprintf(fname,"%s%s_%04d","./","tria",1);
//printf("------> dump_triangle %d (%s)\n",nT,fname);
if(!(fd = fopen(fname, "w")))
{
printf("can't open file `%s' for writing triangles.\n", fname);
}
else
{
/* number of triangles */
blksize=sizeof(int);
SKIP;
my_fwrite(&nT, sizeof(nT), 1, fd);
SKIP;
/* triangles */
blksize=sizeof(f)*nT*3*3;
SKIP;
for (iT=0;iT<nT;iT++)
{
f = Triangles[iT].P[0]->Pos[0];
my_fwrite(&f, sizeof(f), 1, fd);
f = Triangles[iT].P[0]->Pos[1];
my_fwrite(&f, sizeof(f), 1, fd);
f = Triangles[iT].P[0]->Pos[2];
my_fwrite(&f, sizeof(f), 1, fd);
f = Triangles[iT].P[1]->Pos[0];
my_fwrite(&f, sizeof(f), 1, fd);
f = Triangles[iT].P[1]->Pos[1];
my_fwrite(&f, sizeof(f), 1, fd);
f = Triangles[iT].P[1]->Pos[2];
my_fwrite(&f, sizeof(f), 1, fd);
f = Triangles[iT].P[2]->Pos[0];
my_fwrite(&f, sizeof(f), 1, fd);
f = Triangles[iT].P[2]->Pos[1];
my_fwrite(&f, sizeof(f), 1, fd);
f = Triangles[iT].P[2]->Pos[2];
my_fwrite(&f, sizeof(f), 1, fd);
}
SKIP;
/* triangles */
blksize=sizeof(idx)*nT;
SKIP;
for (iT=0;iT<nT;iT++)
{
idx = Triangles[iT].idx;
my_fwrite(&idx, sizeof(idx), 1, fd);
}
SKIP;
close(fd);
}
}
void dump_voronoi(char *fname)
{
int i,j;
FILE *fd = 0;
int blksize;
double f;
int np;
int n;
int next;
//sprintf(fname,"%s%s_%04d","./","tria",1);
if(!(fd = fopen(fname, "w")))
{
printf("can't open file `%s' for writing voronoi.\n", fname);
}
else
{
n = NumPart+NumgPart;
/* number of points */
blksize=sizeof(int);
SKIP;
my_fwrite(&n, sizeof(n), 1, fd);
SKIP;
for (i=0;i<n;i++)
{
next = P[i].ivPoint;
np = P[i].nvPoints;
/* number of points in the cell */
blksize=sizeof(int);
SKIP;
my_fwrite(&np, sizeof(np), 1, fd);
SKIP;
blksize=sizeof(f)*np*3;
SKIP;
for (j = 0; j < np; j++)
{
f = vPoints[next].Pos[0];
my_fwrite(&f, sizeof(f), 1, fd);
f = vPoints[next].Pos[1];
my_fwrite(&f, sizeof(f), 1, fd);
f = 0.5;
my_fwrite(&f, sizeof(f), 1, fd);
next = vPoints[next].next;
}
SKIP;
}
close(fd);
}
}
/************************************************************/
/* */
/* PYTHON INTERFACE */
/* */
/************************************************************/
#ifdef PY_INTERFACE
#include <Python.h>
#include <numpy/arrayobject.h>
PyObject *gadget_GetAllDelaunayTriangles(self, args)
PyObject *self;
PyObject *args;
{
import_array(); /* needed but strange no ? */
PyObject *OutputList;
PyObject *OutputDict;
PyArrayObject *tri = NULL;
npy_intp dim[2];
int iT;
/* loop over all triangles */
OutputList = PyList_New(0);
for (iT=0;iT<nT;iT++)
{
/* 3x3 vector */
dim[0]=3;
dim[1]=3;
tri = (PyArrayObject *) PyArray_SimpleNew(2,dim,PyArray_DOUBLE);
*(double *) (tri->data + 0*(tri->strides[0]) + 0*tri->strides[1]) = Triangles[iT].P[0]->Pos[0];
*(double *) (tri->data + 0*(tri->strides[0]) + 1*tri->strides[1]) = Triangles[iT].P[0]->Pos[1];
*(double *) (tri->data + 0*(tri->strides[0]) + 2*tri->strides[1]) = 0;
*(double *) (tri->data + 1*(tri->strides[0]) + 0*tri->strides[1]) = Triangles[iT].P[1]->Pos[0];
*(double *) (tri->data + 1*(tri->strides[0]) + 1*tri->strides[1]) = Triangles[iT].P[1]->Pos[1];
*(double *) (tri->data + 1*(tri->strides[0]) + 2*tri->strides[1]) = 0;
*(double *) (tri->data + 2*(tri->strides[0]) + 0*tri->strides[1]) = Triangles[iT].P[2]->Pos[0];
*(double *) (tri->data + 2*(tri->strides[0]) + 1*tri->strides[1]) = Triangles[iT].P[2]->Pos[1];
*(double *) (tri->data + 2*(tri->strides[0]) + 2*tri->strides[1]) = 0;
OutputDict = PyDict_New();
PyDict_SetItem(OutputDict,PyString_FromString("id"),PyInt_FromLong(Triangles[iT].idx) );
PyDict_SetItem(OutputDict,PyString_FromString("coord"),(PyObject*)tri);
//(PyObject*)tri
PyList_Append(OutputList, OutputDict );
}
return Py_BuildValue("O",OutputList);
}
PyObject *gadget_GetAllvPoints(self, args)
PyObject *self;
PyObject *args;
{
import_array(); /* needed but strange no ? */
PyArrayObject *pos;
npy_intp ld[2];
int i;
ld[0] = nvPoints;
ld[1] = 3;
pos = (PyArrayObject *) PyArray_SimpleNew(2,ld,PyArray_FLOAT);
for (i = 0; i < pos->dimensions[0]; i++)
{
*(float *) (pos->data + i*(pos->strides[0]) + 0*pos->strides[1]) = vPoints[i].Pos[0];
*(float *) (pos->data + i*(pos->strides[0]) + 1*pos->strides[1]) = vPoints[i].Pos[1];
*(float *) (pos->data + i*(pos->strides[0]) + 2*pos->strides[1]) = 0.5;
}
return PyArray_Return(pos);
}
PyObject *gadget_GetAllvDensities(PyObject* self)
{
import_array(); /* needed but strange no ? */
PyArrayObject *rho;
npy_intp ld[1];
int i;
ld[0] = NumPart;
rho = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE);
for (i = 0; i < rho->dimensions[0]; i++)
{
*(double *) (rho->data + i*(rho->strides[0])) = P[i].Density;
}
return PyArray_Return(rho);
}
PyObject *gadget_GetAllvVolumes(PyObject* self)
{
import_array(); /* needed but strange no ? */
PyArrayObject *volume;
npy_intp ld[1];
int i;
ld[0] = NumPart;
volume = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE);
for (i = 0; i < volume->dimensions[0]; i++)
{
*(double *) (volume->data + i*(volume->strides[0])) = P[i].Volume;
}
return PyArray_Return(volume);
}
PyObject *gadget_GetAllvPressures(PyObject* self)
{
import_array(); /* needed but strange no ? */
PyArrayObject *pressure;
npy_intp ld[1];
int i;
ld[0] = NumPart;
pressure = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE);
for (i = 0; i < pressure->dimensions[0]; i++)
{
*(double *) (pressure->data + i*(pressure->strides[0])) = P[i].Pressure;
}
return PyArray_Return(pressure);
}
PyObject *gadget_GetAllvEnergySpec(PyObject* self)
{
import_array(); /* needed but strange no ? */
PyArrayObject *energy;
npy_intp ld[1];
int i;
double a3inv;
double EnergySpec;
ld[0] = NumPart;
energy = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE);
if(All.ComovingIntegrationOn)
{
a3inv = 1 / (All.Time * All.Time * All.Time);
}
else
a3inv = 1;
for (i = 0; i < energy->dimensions[0]; i++)
{
#ifdef ISOTHERM_EQS
*(double *) (energy->data + i*(energy->strides[0])) = P[i].Entropy;
#else
EnergySpec=P[i].Entropy / GAMMA_MINUS1 * pow(P[i].Density * a3inv, GAMMA_MINUS1);
if (EnergySpec<All.MinEgySpec)
*(double *) (energy->data + i*(energy->strides[0])) = All.MinEgySpec;
else
*(double *) (energy->data + i*(energy->strides[0])) = P[i].Entropy / GAMMA_MINUS1 * pow(P[i].Density * a3inv, GAMMA_MINUS1);
#endif
}
return PyArray_Return(energy);
}
PyObject *gadget_GetAllvAccelerations(PyObject* self)
{
import_array(); /* needed but strange no ? */
PyArrayObject *acc;
npy_intp ld[2];
int i;
ld[0] = NumPart;
ld[1] = 3;
acc = (PyArrayObject *) PyArray_SimpleNew(2,ld,PyArray_FLOAT);
for (i = 0; i < acc->dimensions[0]; i++)
{
*(float *) (acc->data + i*(acc->strides[0]) + 0*acc->strides[1]) = P[i].tesselAccel[0];
*(float *) (acc->data + i*(acc->strides[0]) + 1*acc->strides[1]) = P[i].tesselAccel[1];
*(float *) (acc->data + i*(acc->strides[0]) + 2*acc->strides[1]) = P[i].tesselAccel[2];
}
return PyArray_Return(acc);
}
PyObject *gadget_GetvPointsForOnePoint(self, args)
PyObject *self;
PyObject *args;
{
import_array(); /* needed but strange no ? */
PyArrayObject *pos;
npy_intp ld[2];
int i;
int np=0;
if (!PyArg_ParseTuple(args,"i",&i))
return NULL;
int next;
next = P[i].ivPoint;
np = P[i].nvPoints;
ld[0] = np;
ld[1] = 3;
pos = (PyArrayObject *) PyArray_SimpleNew(2,ld,PyArray_FLOAT);
for (i = 0; i < pos->dimensions[0]; i++)
{
*(float *) (pos->data + i*(pos->strides[0]) + 0*pos->strides[1]) = vPoints[next].Pos[0];
*(float *) (pos->data + i*(pos->strides[0]) + 1*pos->strides[1]) = vPoints[next].Pos[1];
*(float *) (pos->data + i*(pos->strides[0]) + 2*pos->strides[1]) = 0.5;
next = vPoints[next].next;
}
next = vPoints[next].next;
// if (next!=0)
// {
// printf("error in tessel_get_vPointsForOnePoint\n");
// return NULL;
// }
return PyArray_Return(pos);
}
PyObject *gadget_GetNgbPointsForOnePoint(self, args)
PyObject *self;
PyObject *args;
{
import_array(); /* needed but strange no ? */
PyArrayObject *pos;
npy_intp ld[2];
int i,k;
int np=0;
int d=1; /* 1=counterclockwise 2=clockwise */
/* first triangle */
int iT,iTstart;
int p,iP0,iP1,iPnext0;
struct TriangleInList *T,*Tnext,*Tstart;
if (!PyArg_ParseTuple(args,"i",&i))
return NULL;
/*
first loop to count np (number of ngbs)
*/
iTstart = P[i].iT;
Tstart=&Triangles[iTstart];
/* find index in the triangle */
p=-1;
if (Tstart->P[0]==&P[i])
p=0;
else
if (Tstart->P[1]==&P[i])
p=1;
else
if (Tstart->P[2]==&P[i])
p=2;
if (p==-1)
endrun(-712354);
iT = iTstart; /* index of triangle */
iP0 = p; /* index of the ref point "0" in the triangle */
T = Tstart;
np = 0;
/* visit all triangles around the point */
while (1)
{
iP1 = (iP0+d) % 3 ; /* point 1 in the current triangle */
Tnext = T->T[iP1]; /* triangle opposit to iP1 */
if (Tnext==NULL) /* we reach an edge */
{
printf("mmmh... why are we here ?");
endrun(-71235445);
return -1; /* need to increase the size for this point */
}
/* we have a ngb point */
np++;
/* now, find the index of the ref. point in the next triangle */
iPnext0 = (T->idxe[iP1] + d) % 3 ;
/* check */
if (&P[i]!=Tnext->P[iPnext0])
//if (P[i].ID!=Tnext->P[iPnext0]->ID)
{
printf(" problem : initial point %d (id=%d x=%g y=%g z=%g) next-> (id=%d x=%g y=%g z=%g)\n",i,P[i].ID,P[i].Pos[0],P[i].Pos[1],P[i].Pos[2],Tnext->P[iPnext0]->ID,Tnext->P[i]->Pos[0],Tnext->P[i]->Pos[1],Tnext->P[i]->Pos[2]);
endrun(-7123544);
}
/* end of loop */
if (Tnext==Tstart)
break;
iP0 = iPnext0;
T = Tnext;
}
ld[0] = np;
ld[1] = 3;
pos = (PyArrayObject *) PyArray_SimpleNew(2,ld,PyArray_FLOAT);
/*
second loop to record positions
*/
iTstart = P[i].iT;
Tstart=&Triangles[iTstart];
k=0;
/* find index in the triangle */
p=-1;
if (Tstart->P[0]==&P[i])
p=0;
else
if (Tstart->P[1]==&P[i])
p=1;
else
if (Tstart->P[2]==&P[i])
p=2;
if (p==-1)
endrun(-712354);
iT = iTstart; /* index of triangle */
iP0 = p; /* index of the ref point "0" in the triangle */
T = Tstart;
np = 0;
/* visit all triangles around the point */
while (1)
{
iP1 = (iP0+d) % 3 ; /* point 1 in the current triangle */
Tnext = T->T[iP1]; /* triangle opposit to iP1 */
if (Tnext==NULL) /* we reach an edge */
{
printf("mmmh... why are we here ?");
endrun(-71235445);
return -1; /* need to increase the size for this point */
}
/* we have a ngb point */
if (T->P[iP1]->iPref==-1)
{
*(float *) (pos->data + k*(pos->strides[0]) + 0*pos->strides[1]) = T->P[iP1]->Pos[0];
*(float *) (pos->data + k*(pos->strides[0]) + 1*pos->strides[1]) = T->P[iP1]->Pos[1];
*(float *) (pos->data + k*(pos->strides[0]) + 2*pos->strides[1]) = 0;
}
else /* a ghost point */
{
*(float *) (pos->data + k*(pos->strides[0]) + 0*pos->strides[1]) = P[T->P[iP1]->iPref].Pos[0];
*(float *) (pos->data + k*(pos->strides[0]) + 1*pos->strides[1]) = P[T->P[iP1]->iPref].Pos[1];
*(float *) (pos->data + k*(pos->strides[0]) + 2*pos->strides[1]) = 0;
}
k++;
/* now, find the index of the ref. point in the next triangle */
iPnext0 = (T->idxe[iP1] + d) % 3 ;
/* check */
if (&P[i]!=Tnext->P[iPnext0])
//if (P[i].ID!=Tnext->P[iPnext0]->ID)
{
printf(" problem : initial point %d (id=%d x=%g y=%g z=%g) next-> (id=%d x=%g y=%g z=%g)\n",i,P[i].ID,P[i].Pos[0],P[i].Pos[1],P[i].Pos[2],Tnext->P[iPnext0]->ID,Tnext->P[i]->Pos[0],Tnext->P[i]->Pos[1],Tnext->P[i]->Pos[2]);
endrun(-7123544);
}
/* end of loop */
if (Tnext==Tstart)
break;
iP0 = iPnext0;
T = Tnext;
}
return PyArray_Return(pos);
}
PyObject *gadget_GetNgbPointsAndFacesForOnePoint(self, args)
PyObject *self;
PyObject *args;
{
import_array(); /* needed but strange no ? */
PyArrayObject *list;
PyArrayObject *elt;
PyArrayObject *pos1,*pos2,*pos3;
npy_intp ld[2];
int i;
int np=0;
int d=1; /* 1=counterclockwise 2=clockwise */
/* first triangle */
int iT,iTstart;
int p,iP0,iP1,iPnext0;
struct TriangleInList *T,*Tnext,*Tstart;
int nextv,nextv2;
if (!PyArg_ParseTuple(args,"i",&i))
return NULL;
list = PyList_New(0);
iTstart = P[i].iT;
Tstart=&Triangles[iTstart];
/* find index in the triangle */
p=-1;
if (Tstart->P[0]==&P[i])
p=0;
else
if (Tstart->P[1]==&P[i])
p=1;
else
if (Tstart->P[2]==&P[i])
p=2;
if (p==-1)
endrun(-712354);
iT = iTstart; /* index of triangle */
iP0 = p; /* index of the ref point "0" in the triangle */
T = Tstart;
np = 0;
nextv = P[i].ivPoint;
/* visit all triangles around the point */
while (1)
{
iP1 = (iP0+d) % 3 ; /* point 1 in the current triangle */
Tnext = T->T[iP1]; /* triangle opposit to iP1 */
if (Tnext==NULL) /* we reach an edge */
{
printf("mmmh... why are we here ?");
endrun(-71235445);
return -1; /* need to increase the size for this point */
}
/* we have a ngb point */
elt = PyTuple_New(3);
pos1 = PyTuple_New(3);
pos2 = PyTuple_New(3);
pos3 = PyTuple_New(3);
PyTuple_SetItem(pos1, (Py_ssize_t)0, PyFloat_FromDouble((double)T->P[iP1]->Pos[0]));
PyTuple_SetItem(pos1, (Py_ssize_t)1, PyFloat_FromDouble((double)T->P[iP1]->Pos[1]));
PyTuple_SetItem(pos1, (Py_ssize_t)2, PyFloat_FromDouble((double)0.5));
PyTuple_SetItem(pos2, (Py_ssize_t)0, PyFloat_FromDouble((double)vPoints[nextv].Pos[0]));
PyTuple_SetItem(pos2, (Py_ssize_t)1, PyFloat_FromDouble((double)vPoints[nextv].Pos[1]));
PyTuple_SetItem(pos2, (Py_ssize_t)2, PyFloat_FromDouble((double)0.5));
nextv2 = vPoints[nextv].next;
PyTuple_SetItem(pos3, (Py_ssize_t)0, PyFloat_FromDouble((double)vPoints[nextv2].Pos[0]));
PyTuple_SetItem(pos3, (Py_ssize_t)1, PyFloat_FromDouble((double)vPoints[nextv2].Pos[1]));
PyTuple_SetItem(pos3, (Py_ssize_t)2, PyFloat_FromDouble((double)0.5));
PyTuple_SetItem(elt,(Py_ssize_t)0,pos1);
PyTuple_SetItem(elt,(Py_ssize_t)1,pos2);
PyTuple_SetItem(elt,(Py_ssize_t)2,pos3);
PyList_Append(list,elt);
nextv = vPoints[nextv].next;
/* now, find the index of the ref. point in the next triangle */
iPnext0 = (T->idxe[iP1] + d) % 3 ;
/* check */
if (&P[i]!=Tnext->P[iPnext0])
//if (P[i].ID!=Tnext->P[iPnext0]->ID)
{
printf(" problem : initial point %d (id=%d x=%g y=%g z=%g) next-> (id=%d x=%g y=%g z=%g)\n",i,P[i].ID,P[i].Pos[0],P[i].Pos[1],P[i].Pos[2],Tnext->P[iPnext0]->ID,Tnext->P[i]->Pos[0],Tnext->P[i]->Pos[1],Tnext->P[i]->Pos[2]);
endrun(-7123544);
}
/* end of loop */
if (Tnext==Tstart)
break;
iP0 = iPnext0;
T = Tnext;
}
return Py_BuildValue("O",list);
}
PyObject *gadget_GetAllGhostPositions(PyObject* self)
{
PyArrayObject *pos;
npy_intp ld[2];
int i;
import_array(); /* needed but strange no ? */
ld[0] = NumgPart;
ld[1] = 3;
pos = (PyArrayObject *) PyArray_SimpleNew(2,ld,PyArray_FLOAT);
for (i = 0; i < pos->dimensions[0]; i++)
{
*(float *) (pos->data + i*(pos->strides[0]) + 0*pos->strides[1]) = P[NumPart+i].Pos[0];
*(float *) (pos->data + i*(pos->strides[0]) + 1*pos->strides[1]) = P[NumPart+i].Pos[1];
*(float *) (pos->data + i*(pos->strides[0]) + 2*pos->strides[1]) = P[NumPart+i].Pos[2];
}
return PyArray_Return(pos);
}
PyObject *gadget_GetAllGhostvDensities(PyObject* self)
{
import_array(); /* needed but strange no ? */
PyArrayObject *rho;
npy_intp ld[1];
int i;
ld[0] = NumgPart;
rho = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE);
for (i = 0; i < rho->dimensions[0]; i++)
{
*(double *) (rho->data + i*(rho->strides[0])) = P[NumPart+i].Density;
}
return PyArray_Return(rho);
}
PyObject *gadget_GetAllGhostvVolumes(PyObject* self)
{
import_array(); /* needed but strange no ? */
PyArrayObject *volume;
npy_intp ld[1];
int i;
ld[0] = NumgPart;
volume = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE);
for (i = 0; i < volume->dimensions[0]; i++)
{
*(double *) (volume->data + i*(volume->strides[0])) = P[NumPart+i].Volume;
}
return PyArray_Return(volume);
}
#endif /* PY_INTERFACE */
#endif
c4science · Help