/**
# Contact angles on an embedded body

This file is used to impose contact angles on arbitrary embedded bodies
described using [EMBED] (embed.h).

*/

//#include "fractions.h"
#include "embed.h"
#include "geometrical_tools.h"

#define Pi 3.14159265358979323846

extern scalar tag, triple_cell, f1, Alpha;

/** We compute the normal corresponding to the angle we want to impose
on the contact cell knowing the initial volume fraction in that cell. */
coord normal_contact (coord ns, coord nf1, double angle){
  double det = ns.x*nf1.y - ns.y*nf1.x; // 2D cross product
  coord n;
  if ( det > 0){
    n.x = ns.x*cos(angle)-ns.y*sin(angle); // new fluid normal with the right contact angle
    n.y = ns.x*sin(angle)+ns.y*cos(angle);
  }
  else{
    n.x =  ns.x*cos(angle)+ns.y*sin(angle);
    n.y = -ns.x*sin(angle)+ns.y*cos(angle);
  }
  return n;
}              


void triple_cell_dectection(double angle, scalar f, scalar cs, vector ns, vector nf1, vector nf){
  scalar alpha1[],alpha2[];
  foreach(){
    tag[] = triple_cell[] = 0;
    int n=0;
    foreach_neighbor(){
      if(cs[]<1. && cs[]>0. && f[]<1.  && f[]>0. ){
        /** We locate the exact cell containing the triple point */
        coord b, ntemp, ntemp1, ntemp2, p_vof[2], p_embed[2];
        
        /** vof geometry */
        ntemp1 = interface_normal (point, f); // normal to the fluid interface
        foreach_dimension() nf1.x[] = ntemp1.x ; 
        //alpha1[] = line_alpha(f[], ntemp1);
				alpha1[] = line_alpha(f[],(coord){nf1.x[], nf1.y[]});
        /** embed solid geometry */
        
        embed_geometry (point, &b, &ntemp2);//ntemp = facet_normal(point, cs, fs);
        //embed_geometry (point, &b, (coord*){ns.x[], ns.y[]});
        //embed_geometry (point, &b, &ns[]);

       	foreach_dimension() ns.x[] = -ntemp2.x ; 
        alpha2[] = line_alpha(cs[], ntemp2);

        /** Does VOF interface intersects the embed solid in this cell ? */
        if ( facets (ntemp1, alpha1[], p_vof)==2 && facets (ntemp2, alpha2[], p_embed)==2  ){
          if (Intersect(p_vof[0],p_vof[1],p_embed[0],p_embed[1])) {
              /** Then the triple point lies in here */
              triple_cell[]=1; n++;

              double piangle = (angle*Pi)/180. ;

              //foreach_dimension() ntemp2.x = ns.x[] ; 
              ntemp = normal_contact((coord){ns.x[], ns.y[]}, (coord){nf1.x[], nf1.y[]}, piangle);

              foreach_dimension() nf.x[] = ntemp.x;

              Alpha[] = line_alpha(f[], ntemp); // fictive alpha is stored
          }
          else nf.x[] = nf.y[] = ns.x[] = ns.y[] = nodata;
        }
      }
    }
    if(n && cs[]==0.) tag[] =1.;
  }
}

void fraction_reconstruction(scalar f, vector nf){
	foreach(){
    if(tag[]){  
      int stencil=2; 
      for (int i1 = -stencil; i1 <= stencil; i1++){
        for (int j1 = -stencil; j1 <= stencil; j1++){
          if(triple_cell[i1,j1] ){
            coord a, b;

            a.x = -0.5 - i1; 
            a.y = -0.5 - j1;
            b.x =  0.5 - i1;
            b.y =  0.5 - j1;
            
          	f1[] = rectangle_fraction ((coord){nf.x[i1,j1], nf.y[i1,j1]}, Alpha[i1,j1], a, b);

          }
        }
      }  
  	}
 	}
 	boundary({f1});

  foreach(){
    if (tag[]){  
      f[]=f1[];
    }
  }
  boundary({f});
}