Point start_point(scalar f,double xs, double ys){
double dist = sq(L0*pow(2.,0.5));
foreach(reduction(min:dist)){
if (f[]>0. && f[]<1.&&fabs(x-X0)>Delta&&fabs(x-(X0+L0))>Delta&&fabs(y-(Y0+L0))>Delta&&fabs(y-Y0)>Delta){
double dis = (sq(x-xs)+sq(y-ys));
if (dis<dist)
dist = dis;
}
}
double xp=(2*L0)+X0;// These coordinates cannot be located.
double yp=(2*L0)+Y0;
foreach(){
if (sq(x-xs)+sq(y-ys)==dist && f[]>0 && f[]<1.){ // Only true on a single PID.
xp=x;
yp=y;
}
}
return locate(xp,yp);
}
void find_facets(Point point,scalar f,double xy[4]){
coord n;
n = mycs (point, f);
double alpha = plane_alpha (f[], n);
coord segment[2];
if (facets (n, alpha, segment) == 2){
xy[0] = x + segment[0].x*Delta;
xy[1] = y + segment[0].y*Delta;
xy[2] = x + segment[1].x*Delta;
xy[3] = y + segment[1].y*Delta;
}else{
fprintf(stdout,"Warning:\nCould not find facets; expect unexpected behaviour.\n");
}
}
void next_point(Point point,scalar f,double xp,double yp,int dir[2],int * iter,double * xu,double * yu){
int nb=0;
foreach_dimension()
nb+=((f[-1]>0.)*(f[-1]<1.))+((f[1]>0.)*(f[1]<1.));
if (nb==2&&(fabs(dir[0])+fabs(dir[1])==1)){ // Only two neighbouring interfacial cells could be easy...
for(int ii = -1;ii<=1;ii++){
for(int jj = -1;jj<=1;jj++){
if ((abs(ii)+(abs(jj))==1) && (dir[0]!=-ii || dir[1]!=-jj) && f[ii,jj]>0. && f[ii,jj]<1.){// Simple neighbour
dir[0]=ii;
dir[1]=jj;
if (!is_refined(neighbor(ii,jj))){//Special care for refined next cells
*xu=(5*xp+x)/6.+((double)ii*Delta/3.);
*yu=(5*yp+y)/6.+((double)jj*Delta/3.);
return;
}else if(fine(f,(xp>x)+ii,(yp>y)+jj)>0 && fine(f,(xp>x)+ii,(yp>y)+jj)<1){//still OK
*xu=(5*xp+x)/6.+((double)ii*Delta/3.);
*yu=(5*yp+y)/6.+((double)jj*Delta/3.);
return;
}else if(ii==0){//y-dir
if(fine(f,(xp<x),((yp>y))+jj)>0 && fine(f,(xp<x),(yp>y)+jj)<1){ //poorly reconstructed, do not follow facet
*xu=x+(0.5*Delta*((xp<x)-0.5));
*yu=y+(0.75*Delta*jj) ;
return;
}
}else if(fine(f,(xp>x)+ii,(yp<y))>0 && fine(f,(xp>x)+ii,(yp<y))<1){ //x-dir, poorly reconstructed, do not follow facet
*yu=y+(0.5*Delta*((yp<y)-0.5));
*xu=x+(0.75*Delta*ii);
return;
}
}
}
}
}
if (fabs(xp-X0)<Delta||fabs(xp-(X0+L0))<=Delta||fabs(yp-Y0)<Delta||fabs(yp-(Y0+L0))<Delta){//domain boundary cell: stop tracing
*iter += (1<<25);// add large number
*xu=X0+L0/2.;
*yu=Y0+L0/2.;
return;
}
// If previous efforts have not worked, tryo to follow the interface
if (xp==x+Delta/2&&((f[1]<1.&&f[1]>0.)||(!is_leaf(neighbor(1,0))&&coarse(f,1,0)<1&&coarse(f,1,0)>0))){
dir[0]=1.;
dir[1]=0;
*xu=xp+Delta/10.;
*yu=yp;
return;
}
else if (xp==x-Delta/2.&&((f[-1]<1.&&f[-1]>0.)||(!is_leaf(neighbor(-1,0))&&coarse(f,-1,0)<1&&coarse(f,-1,0)>0))){
dir[0]=-1;
dir[1]=0;
*xu=xp-Delta/10.;
*yu=yp;
return;
}
else if (yp==y+Delta/2.&&((f[0,1]<1.&&f[0,1]>0.)||(!is_leaf(neighbor(0,1))&&coarse(f,0,1)<1&&coarse(f,0,1)>0))){
dir[0]=0;
dir[1]=1;
*xu=xp;
*yu=yp+Delta/10.;
return;
}
else if (yp==y-Delta/2.&&((f[0,-1]<1.&&f[0,-1]>0.)||(!is_leaf(neighbor(0,-1))&&coarse(f,0,-1)<1&&coarse(f,0,-1)>0))) {
dir[0]=0;
dir[1]=-1;
*xu=xp;
*yu=yp-Delta/10;
return;
}
//As a last resort, try a diagnonal neighbour, blindly. If this does not go well, we are in in trouble
double fac[4];
find_facets(point,f,fac);
double meanx = (fac[0]+fac[2])/2.;
double meany = (fac[1]+fac[3])/2.;
double dx,dy;
if (yp>meany){
dy= 0.3;
}else
dy = -0.3;
if (xp>meanx){
dx= 0.3;
}else
dx = -0.3;
dir[0]=0; dir[1]=0;
*xu=xp+Delta*dx;
*yu=yp+Delta*dy;
return;
}
void update_new_facet(Point point, double xyn[2],double xyf[4],int dir[2]){
if ((fabs(dir[0])+fabs(dir[1])==1)){ // Try to make use of the direction.
for (int g=-1;g<=1;g+=2){
if (dir[0]==g){ //x-dir
double dg = (double)g;
if(dg*xyf[0]==min(dg*xyf[0],dg*xyf[2])){
xyn[0]=xyf[2];xyn[1]=xyf[3];
return;
}else if(dg*xyf[2]==min(dg*xyf[0],dg*xyf[2])){
xyn[0]=xyf[0];xyn[1]=xyf[1];
return;
}
}else if(dir[1]==g){ //y-dir
double dg = (double)g;
if(dg*xyf[1]==min(dg*xyf[1],dg*xyf[3])){
xyn[0]=xyf[2];xyn[1]=xyf[3];
return;
}else if(dg*xyf[3]==min(dg*xyf[1],dg*xyf[3])){
xyn[0]=xyf[0];xyn[1]=xyf[1];
return;
}
}
}
}
// Else choose the one farest away from the previous. This is not very robust and is the main cause of 'bounce backs'.
if ((sq(xyf[0]-xyn[0])+sq(xyf[1]-xyn[1])) <= (sq(xyf[2]-xyn[0])+sq(xyf[3]-xyn[1]))){
xyn[0]=xyf[2]; xyn[1]=xyf[3];
}else{
xyn[0]=xyf[0]; xyn[1]=xyf[1];
}
}
void loop_interfacial_cells(FILE * fp, scalar f,scalar tag,double xynn[2]){
if (dimension!=2){
fprintf(stdout,"2D only!\n");
if (pid()==0)
fprintf(fp,"# 2D only!\n");
return;
}
boundary({f});//neighbors on trees etc.
double xyn[2];
xyn[0]=xynn[0]; xyn[1]=xynn[1];
int ifc= 0;
foreach(reduction(+:ifc)){
if ( f[]>0. && f[]<1.)
ifc++;
}
double xu,yu;
Point point;
point = start_point(f,xyn[0],xyn[1]);
// int i = 1;
//static FILE * fpp = popen ("gfsview2D tag.gfv -s","w"); //For trouble shooting
int igg=0;
double l =0,cur;
double xyf[4];
int dir[2]={0,0};
while (igg <= ifc){
xu=X0-L0;
yu=Y0-L0;
if (point.level>0){ // Local pid()
find_facets(point,f,xyf);
l+=pow(sq(xyf[0]-xyf[2])+sq(xyf[1]-xyf[3]),0.5);
cur = tag[];
//tag[]+=100.; // For trouble shooting
update_new_facet(point,xyn,xyf,dir); // For tracing the interface in the correct direction
next_point(point,f,xyn[0],xyn[1],dir,&igg,&xu,&yu); // Aims to Update xu and yu to be in next leaf cell along the interface
}else{ // Set small values for the tracing variables on the other threads
cur=-HUGE;
xyf[0]=xu; xyf[1]=yu; xyf[2]=xu; xyf[3]=yu;
xyn[0]=xu;xyn[1]=yu;
dir[0]=-2;dir[1]=-2;
}
@ if _MPI
// Set send buffers, so reduced variables can be recieved in the original ones.
double oldx=xu;double oldy=yu;double oldl = l; double xyfold[4];
xyfold[0]=xyf[0];xyfold[1]=xyf[1];xyfold[2]=xyf[2];xyfold[3]=xyf[3];
double oldcur=cur; double oldxyn[2];
oldxyn[0]=xyn[0];oldxyn[1]=xyn[1];
int olddir[2];
olddir[1]=dir[1];olddir[0]=dir[0];
int oldi=igg;
// Update all workers, just in case one has to take over, this should be improved for better performance.
MPI_Allreduce(&oldx, &xu, 1, MPI_DOUBLE, MPI_MAX,MPI_COMM_WORLD);
MPI_Allreduce(&oldy, &yu, 1, MPI_DOUBLE, MPI_MAX,MPI_COMM_WORLD);
MPI_Allreduce(&xyfold, &xyf, 4, MPI_DOUBLE, MPI_MAX,MPI_COMM_WORLD);
MPI_Allreduce(&oldxyn, &xyn, 2, MPI_DOUBLE, MPI_MAX,MPI_COMM_WORLD);
MPI_Allreduce(&oldl, &l, 1, MPI_DOUBLE, MPI_MAX,MPI_COMM_WORLD);
MPI_Allreduce(&oldcur, &cur, 1, MPI_DOUBLE, MPI_MAX,MPI_COMM_WORLD);
MPI_Allreduce(&olddir, &dir, 2, MPI_INT, MPI_MAX,MPI_COMM_WORLD);
MPI_Allreduce(&oldi, &igg, 1, MPI_INT, MPI_MAX,MPI_COMM_WORLD);
@ endif
if (pid()==0&&igg<ifc){ // Our favorite worker writes in the file
fprintf(fp,"%d\t%g\t%g\t%g\t%g\n",igg,(xyf[0]+xyf[2])/2,(xyf[1]+xyf[3])/2,l,cur);
fflush(fp);
}
// All threads try to find the new location.
point = locate(xu,yu);
igg++;
// i = igg;
//output_gfs(fpp); //For trouble shooting purposes:
}
}
