sandbox/lchirco/signature.h
We want to detect thin sheets and ligaments in 2 phase flows. First, we compute the quadratic form f (x, x) = x_i x_j T_{ij}, where T_{ij} are the quadratic moments, by integrating over a shell of arbitrary thickness and radius. Finally, we compute the signature s of the quadratic form.
int find_moments_level(scalar f, double length, double target_delta){
int level = depth();
bool found = false;
while (level >= 0 && !found){
double num = 1<<level;
if (length/num > target_delta) found = true;
level--;
}
return level;
}
To compute the eigenvalues used in the signature method in 3D we use the GNU scientific library.
#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>
void eigen ( double data[], double tau[], double dim){
gsl_matrix_view m
= gsl_matrix_view_array (data, dim, dim);
gsl_vector_complex *eval = gsl_vector_complex_alloc (dim);
gsl_matrix_complex *evec = gsl_matrix_complex_alloc (dim, dim);
gsl_eigen_nonsymmv_workspace * w =
gsl_eigen_nonsymmv_alloc (dim);
gsl_eigen_nonsymmv (&m.matrix, eval, evec, w);
gsl_eigen_nonsymmv_free (w);
gsl_eigen_nonsymmv_sort (eval, evec,
GSL_EIGEN_SORT_ABS_DESC);
for (int i = 0; i < dim; i++)
{
gsl_complex eval_i
= gsl_vector_complex_get (eval, i);
gsl_vector_complex_view evec_i
= gsl_matrix_complex_column (evec, i);
// printf ("eigenvalue = %g + %gi\n", GSL_REAL(eval_i), GSL_IMAG(eval_i));
tau[i]=GSL_REAL(eval_i);
}
gsl_vector_complex_free(eval);
gsl_matrix_complex_free(evec);
}
This function is used to compute the signature at a given level of refinement lev
. It can be found using the function find_moments_level
.
void compute_signature_neigh_level(scalar f, scalar phii, scalar s, int lev){
scalar int_xx[], int_xy[], int_yy[];
#if dimension == 3
scalar int_xz[], int_yz[], int_zz[];
#endif
vector Tij[], sign[];
foreach_level(lev){
int_xx[] = 0;
int_xy[] = 0;
int_yy[] = 0;
#if dimension == 3
int_xz[] = 0;
int_yz[] = 0;
int_zz[] = 0;
#endif
}
boundary ({phii});
We integrate over the cells of the 5x5 stencil and subtract the values of those in the 3x3 stencil. We can say the thickness of the shell is equal to the grid-size and the radius is twice the grid-size
double mom_xx, mom_yy, mom_xy;
#if dimension == 3
double mom_xz, mom_yz, mom_zz;
#endif
foreach_level(lev){
mom_xx = mom_xy = mom_yy = 0.;
#if dimension == 3
mom_xz = mom_yz = mom_zz = 0.;
#endif
// if (phii[] > -0.99){
double xp = x;
double yp = y;
#if dimension == 3
double zp = z;
#endif
foreach_neighbor() {
mom_xx += sq(x - xp)*phii[];
mom_yy += sq(y - yp)*phii[];
mom_xy += (x - xp)*(y - yp)*phii[];
#if dimension == 3
mom_xz += (x - xp)*(z - zp)*phii[];
mom_yz += (y - yp)*(z - zp)*phii[];
mom_zz += sq(z - zp)*phii[];
#endif
}
foreach_neighbor(1) {
mom_xx -= sq(x - xp)*phii[];
mom_yy -= sq(y - yp)*phii[];
mom_xy -= (x - xp)*(y - yp)*phii[];
#if dimension == 3
mom_xz -= (x - xp)*(z - zp)*phii[];
mom_yz -= (y - yp)*(z - zp)*phii[];
mom_zz -= sq(z - zp)*phii[];
#endif
}
int_xx[] = mom_xx;
int_yy[] = mom_yy;
int_xy[] = mom_xy;
#if dimension == 3
int_xz[] = mom_xz;
int_yz[] = mom_yz;
int_zz[] = mom_zz;
#endif
}
foreach_level(lev)
foreach_dimension()
Tij.x[] = 0;
We now compute the eigenvalues of the quadratic form. For a 2x2 symmetric matrix the eigenvalues can be found analytically. In the 3D case we call the function eigen
that uses the GSL library.
#if dimension == 2
double Txx, Tyy, Txy;
double tau[2]={0.};
foreach_level(lev){
Txx = int_xx[]/sq(Delta);
Tyy = int_yy[]/sq(Delta);
Txy = int_xy[]/sq(Delta);
double data[] = { Txx, Txy,
Txy, Tyy };
eigen(data, tau, 2);
Tij.x[] = tau[0];
Tij.y[] = tau[1];
}
#else //dimension == 3
double Txx, Tyy, Tzz, Txy, Txz, Tyz=0.;
double tau[3]={0.};
foreach_level(lev){
Txx = int_xx[]/sq(Delta);
Tyy = int_yy[]/sq(Delta);
Txy = int_xy[]/sq(Delta);
Tzz = int_zz[]/sq(Delta);
Txz = int_xz[]/sq(Delta);
Tyz = int_yz[]/sq(Delta);
double data[] = { Txx, Txy, Txz,
Txy, Tyy, Tyz,
Txz, Tyz, Tzz };
eigen(data, tau, 3);
Tij.x[] = tau[0];
Tij.y[] = tau[1];
Tij.z[] = tau[2];
}
#endif
The signature of the quadratic form are the signs of the eigenvalues.
foreach_level(lev){
foreach_dimension(){
sign.x[] = fabs(Tij.x[])/(Tij.x[] + 1.e-10);
if (fabs(Tij.x[]) < 10.) sign.x[] = 0;
}
}
We identify thin ligaments by looking at the signature. In 2D we can have the signatures:
- (+,+) bulk of phase
- (+,-) thin ligament
- (+,0) or (-,0) interface
- (-,-) outside of phase
#if dimension ==2
foreach_level(lev){
if (sign.x[]*sign.y[] > 0.999 && sign.x[] > 0.999) s[] = 2; //bulk of phase
if (sign.x[]*sign.y[] > 0.999 && sign.x[] < -0.999) s[] = -1; //outside of phase
if (fabs(sign.x[]*sign.y[]) < 0.01) s[] = 0; //interface
if (sign.x[]*sign.y[] < -0.999) s[] = 1; //ligament
}
#else //dimension == 3
foreach_level(lev){
s[] = 1; //thin
if (sign.x[] == 0. || sign.y[] == 0. || sign.z[] == 0.) s[] = 0; //interface
if (sign.x[] > 0.999 && sign.y[] > 0.999 && sign.z[] > 0.999) s[] = 2; //bulk of phase
if (sign.x[] < -0.999 && sign.y[] < -0.999 && sign.z[] < -0.999) s[] = -1; //outside of phase
}
#endif
}
This function is used to perforate thin sheets. The scalar s
must be filled using the compute_signature_neigh_level
function.
void change_topology (scalar f, scalar s, scalar M, int lev, const int max_change, bool large){
double f_avg[max_change]; // average f in the neighbor
int num = 0.; int i_change = 0;
srand(time(NULL)); // Initialization of random seed.
int r;
Cache to_perf = {0};
#ifndef SQUARES
int i;
scalar hole[];
double xc[max_change], yc[max_change], zc[max_change];
#endif
foreach_level(lev){
M[] = 0.;
r = rand() % 50;
if (s[] > 0.99 && s[] < 1.01 && r==0 && i_change < max_change){
f_avg[i_change] = 0.;
num = 0;
#ifndef SQUARES
xc[i_change] = x; yc[i_change] = y; zc[i_change] = z;
#endif
#ifdef SQUARES
if (large == false){
foreach_neighbor(1){ //compute average f
f_avg[i_change] += f[];
cache_append (&to_perf, point, i_change);
num += 1;
}
}
else{
foreach_neighbor(){ //compute average f
f_avg[i_change] += f[];
cache_append (&to_perf, point, i_change);
num += 1;
}
}
f_avg[i_change]/=num;
cache_shrink (&to_perf);
#endif
i_change++;
}
}
#ifndef SQUARES
printf("Holes are spheres \n");
vertex scalar phi[];
double Ddelta = 0.001/64.;
double R = 1.2*2.5*Ddelta;
foreach_vertex() {
phi[] = HUGE;
for (i = 0; i < max_change; i++){
phi[] = intersection (phi[], (sq(x - xc[i])/sq(1)+sq(y -yc[i])/sq(1)+sq(z -zc[i])/sq(1) - sq(R)));
// phi[] = intersection (phi[], (/*sq(x - xc[i]) +*/ (sq(y - yc[i]) + sq(z - zc[i]) - sq(R))));
}
}
boundary ({phi});
fractions (phi, hole);
foreach(){
if (f[] > 1e-4 && hole[] < 1 - 1.e-4){
f[] = hole[];
}
}
boundary({f});
#endif
#ifdef SQUARES
printf("Holes are squares \n");
foreach_cache(to_perf){
s[] = -2.;
M[] = f[] - (1 - f_avg[_flags]);
if (f_avg[_flags] < 0.3) {
s[] = -3.;
f[] = 1.;
}
else
f[] = 0.;
}
#endif
for (int ilev = lev; ilev < depth(); ilev++)
foreach_level(ilev){
s.prolongation = M.prolongation = f.prolongation = refine_injection;
if(is_refined(cell) && s[] == -2){
s.prolongation (point, s);
M.prolongation (point, M);
f.prolongation (point, f);
}
}
free(to_perf.p);
boundary(all);
}