/* transform.c
 * by Francois Labelle (http://www.cs.berkeley.edu/~flab/)
 * December 3, 1999
 *
 * Linear transformation functions used for anisotropic meshing. */

#include <math.h>
#include "mesh.h"

/*= TOOLS ===================================================================*/

void Eigen(struct EllipticityTensor *T, real c, real s, real l1, real l2)

/* Specify a real symmetric matrix by its eigenvectors instead
 *
 * (c, s): direction of the first eigenvector (no need to normalize)
 * l1:     first eigenvalue, i.e. compression along (c,s)
 * l2:     second eigenvalue, i.e. compression orthogonal to (c,s)  */

{
    real d2, r;
    
    /* normalize c and s */
    d2 = c*c + s*s;
    if (d2 == 0) {
        c = 1;      // in fact, should do this as soon as IEEE
        s = 0;      // arithmetic starts its "graceful" underflow
    }
    else {
        r = 1.0/sqrt(d2);
        c *= r;
        s *= r;
    }
    
    /* compute ellipticity tensor */
    T->ix = (c*c) * l1 + (s*s) * l2;
    T->jy = (s*s) * l1 + (c*c) * l2;
    T->jx = T->iy = (c*s) * (l1 - l2);
}

/*= SOME EXAMPLES ===========================================================*/

/* You can use either the absolute coordinates (x, y), or the normalized
 * coordinates (nx, ny) for your calculations.
 *
 * Normalized coordinates lie in the box [-1,1] x [-1,1] */

/*---------------------------------------------------------------------------*/

void T_hstretch(struct EllipticityTensor *T, real x, real y, real nx, real ny)
{
    T->ix = 0.5; T->jx = 0.0;
    T->iy = 0.0; T->jy = 2.0;
}

/*---------------------------------------------------------------------------*/

void T_vstretch(struct EllipticityTensor *T, real x, real y, real nx, real ny)
{
    T->ix = 2.0; T->jx = 0.0;
    T->iy = 0.0; T->jy = 0.5;
}

/*---------------------------------------------------------------------------*/

void T_dstretch(struct EllipticityTensor *T, real x, real y, real nx, real ny)
{
    T->ix =  1.25; T->jx = -0.75;
    T->iy = -0.75; T->jy =  1.25;
}

/*---------------------------------------------------------------------------*/

void T_sink(struct EllipticityTensor *T, real x, real y, real nx, real ny)
{
    real s = sqrt(sqrt(nx*nx+ny*ny)) + 1.0;
    
    Eigen(T, nx, ny, 1.0/s, s);
}

/*---------------------------------------------------------------------------*/

void T_swirl(struct EllipticityTensor *T, real x, real y, real nx, real ny)
{
    real s = sqrt(sqrt(nx*nx+ny*ny)) + 1.0;
    
    Eigen(T, nx, ny, s, 1.0/s);
}

/*---------------------------------------------------------------------------*/

void T_hline(struct EllipticityTensor *T, real x, real y, real nx, real ny)
{
    real s = 1.0/(fabs(ny) + 0.4);
    
    T->ix = s; T->jx = 0;
    T->iy = 0; T->jy = s;
}

/*---------------------------------------------------------------------------*/

void T_vline(struct EllipticityTensor *T, real x, real y, real nx, real ny)
{
    real s = 1.0/(fabs(nx) + 0.4);
    
    T->ix = s; T->jx = 0;
    T->iy = 0; T->jy = s;
}

/*---------------------------------------------------------------------------*/

void T_dline(struct EllipticityTensor *T, real x, real y, real nx, real ny)
{
    real s = 1.0/(0.707*fabs(nx-ny) + 0.4);
    
    T->ix = s; T->jx = 0;
    T->iy = 0; T->jy = s;
}

/*---------------------------------------------------------------------------*/

void T_center(struct EllipticityTensor *T, real x, real y, real nx, real ny)
{
    real s = 1.0/(sqrt(nx*nx+ny*ny) + 0.4);
    
    T->ix = s; T->jx = 0;
    T->iy = 0; T->jy = s;
}

/*---------------------------------------------------------------------------*/

void T_perimeter(struct EllipticityTensor *T, real x, real y, real nx, real ny)
{
    real s = sqrt(nx*nx+ny*ny) + 0.4;
    
    T->ix = s; T->jx = 0;
    T->iy = 0; T->jy = s;
}

/*---------------------------------------------------------------------------*/

void T_left(struct EllipticityTensor *T, real x, real y, real nx, real ny)
{
    real s = exp(-nx) + 0.6;
    
    T->ix = s; T->jx = 0;
    T->iy = 0; T->jy = s;
}

/*---------------------------------------------------------------------------*/

void T_right(struct EllipticityTensor *T, real x, real y, real nx, real ny)
{
    real s = exp(nx) + 0.6;
    
    T->ix = s; T->jx = 0;
    T->iy = 0; T->jy = s;
}

/*---------------------------------------------------------------------------*/

void T_top(struct EllipticityTensor *T, real x, real y, real nx, real ny)
{
    real s = exp(ny) + 0.6;
    
    T->ix = s; T->jx = 0;
    T->iy = 0; T->jy = s;
}

/*---------------------------------------------------------------------------*/

void T_bottom(struct EllipticityTensor *T, real x, real y, real nx, real ny)
{
    real s = exp(-ny) + 0.6;
    
    T->ix = s; T->jx = 0;
    T->iy = 0; T->jy = s;
}

/*= CUSTOM TRANSFORMATION  ==================================================*/

/* For each point (x, y) in the physical space, assign a linear
 * transformation "T".
 *
 * The determinant of T should be positive, otherwise it
 * means that you favor inverted triangles!
 *
 * The transformation should be continuous over the space, else the mesher
 * may keep generating triangles around discontinuities. */

void T_custom(struct EllipticityTensor *T, real x, real y, real nx, real ny)
{
    real angle = atan2(ny, nx);
    real dist = sqrt(nx*nx + ny*ny);
    real s = 1.0 + 7*(0.6 + 0.25*sin(5*angle) - dist);
    
    if (s < 0.4) s = 0.5;
    if (s > 2.5) s = 2.5;
    
    T->ix = s; T->jx = 0;
    T->iy = 0; T->jy = s;
}

/*---------------------------------------------------------------------------*/