/*
 * MigFlow - Copyright (C) <2010-2018>
 * <Universite catholique de Louvain (UCL), Belgium
 *  Universite de Montpellier, France>
 * 	
 * List of the contributors to the development of MigFlow: see AUTHORS file.
 * Description and complete License: see LICENSE file.
 * 	
 * This program (MigFlow) is free software: 
 * you can redistribute it and/or modify it under the terms of the GNU Lesser General 
 * Public License as published by the Free Software Foundation, either version
 * 3 of the License, or (at your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Lesser General Public License for more details.
 * 
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program (see COPYING and COPYING.LESSER files).  If not, 
 * see <http://www.gnu.org/licenses/>.
 */

#include "tools.h"
#include <stdlib.h>
#include <stdio.h>
#include <float.h>
#include <math.h>
#include "fluid_problem.h"
#include "mesh_find.h"
#include "fluid_problem_fem.h"

#define P  D
#define U 0
int global_count = 0;
int global_ind = 0;
#define D DIMENSION


#define LOCAL_MATRIX(i,j,a,b) local_matrix[((i)*n_fields*N_SF+N_SF*(a)+(j))*n_fields+(b)]

static void particle_force_f(FluidProblem *problem, double *f, double *dfdu, double *dfddp, double *sol, double *dsol, const double *solold, const double c, const double *dc, double a, double dt, int iel, int ip);
void fluid_problem_interpolate_rho_and_nu(const FluidProblem *problem, double a, double *rho, double *mu) {
  if (problem->n_fluids == 1){
    *rho = problem->rho[0];
    *mu = problem->mu[0];
  }
  else{
    //*rho = 1/(1/problem->rho[0]*a + 1/problem->rho[1]*(1-a));
    *rho = problem->rho[0]*a + problem->rho[1]*(1-a);
    //*mu = problem->mu[0]*a + problem->mu[1]*(1-a);
    *mu = (problem->mu[0]/problem->rho[0]*a + problem->mu[1]/problem->rho[1]*(1-a))*(*rho);
  }
}

inline static double mesh_dxidx(const Mesh *mesh, int iel, double dxidx[D][D]) {
  double dxdxi[D][D];
  const int *el = &mesh->elements[iel*N_N];
  for (int i = 0; i < D; ++i)
    for (int j = 0; j < D; ++j)
      dxdxi[i][j] = mesh->x[el[j+1]*3+i] - mesh->x[el[0]*3+i];
  return invDD(dxdxi, dxidx);
}

inline static double element_volume_from_detj(double detj) {
#if D == 2
    return detj/2;
#else        
    return detj/6;
#endif    
}

inline static double node_volume_from_detj(double detj) {
#if DIMENSION == 2
    return detj/6;
#else
    return detj/24;
#endif
}

int fluid_problem_n_fields(const FluidProblem *p) {return D+1;}


static void node_force_volume(FluidProblem *problem, const double *solution_old, double dt, double *all_local_vector, double *all_local_matrix)
{
  const Mesh *mesh = problem->mesh;
  const double *porosity = problem->porosity;
  const double *solution = problem->solution;
  int drag_in_stab = problem->drag_in_stab;
  int n_fields = fluid_problem_n_fields(problem);
  size_t local_size = N_SF*n_fields;
  double *s = malloc(sizeof(double)*n_fields);
  double *sold = malloc(sizeof(double)*n_fields);
  double *ds = malloc(sizeof(double)*n_fields*D);
  for (int ip=0; ip < problem->n_particles; ++ip) {
    int iel = problem->particle_element_id[ip];
    if(iel < 0){
      for (int d = 0; d < D; ++d)
        problem->particle_force[ip*D+d] = 0;
      problem->particle_force[ip*D+1] = problem->g*problem->particle_mass[ip];
      continue;
    }
    double *xi = problem->particle_uvw + D*ip;
    double phi[N_SF];
    shape_functions(xi,phi);
    const int *el = &mesh->elements[iel*N_N];
    double dxidx[D][D], dphi[N_SF][D];
    double detj = mesh_dxidx(mesh, iel, dxidx);
    grad_shape_functions(dxidx, dphi);
    double *local_vector = all_local_vector ? &all_local_vector[local_size*iel] : NULL;
    double *local_matrix = all_local_matrix ? &all_local_matrix[local_size*local_size*iel] : NULL;
    for (int i = 0; i < n_fields; ++i) {
      s[i] = 0;
      sold[i] = 0;
      for (int j = 0; j < D; ++j) {
        ds[i*D+j] = 0;
      }
    }
    double c = 0;
    double a = 0;
    double cold = 0;
    double dc[D] = {0};
    for (int i = 0; i < N_SF; ++i) {
      c += problem->porosity[el[i]]*phi[i];
      cold += problem->oldporosity[el[i]]*phi[i];
      if (problem->n_fluids == 2) {
        a += problem->concentration[iel*N_N+i]*phi[i];
      }
      for (int j = 0; j < n_fields; ++j) {
        double dof = solution[el[i]*n_fields+j];
        double dofold = solution_old[el[i]*n_fields+j];
        s[j] += phi[i]*dof;
        sold[j] += phi[i]*dofold;
        for (int k = 0; k < D; ++k)
          ds[j*D+k] += dphi[i][k]*dof;
      }
      for (int k=0; k<D; ++k){
        dc[k] += problem->porosity[el[i]]*dphi[i][k];
      }
    }
    double f[D],dfdu,dfddp;
    particle_force_f(problem,f,&dfdu,&dfddp,s,ds,sold,c,dc,a,dt,iel,ip);
    if (!local_vector)
      continue;
    double rho,nu;
    fluid_problem_interpolate_rho_and_nu(problem,a, &rho,&nu);
    double pspg = drag_in_stab?problem->taup[iel]/rho:0;
    for (int iD = 0; iD < D; ++iD) {
      for (int iphi = 0; iphi < N_SF; ++iphi){
        local_vector[iphi+N_SF*(U+iD)] += phi[iphi]*f[iD];
        local_vector[iphi+N_SF*P] += pspg*dphi[iphi][iD]*f[iD];
        for (int jD = 0; jD < D; ++jD) {
          double supg = drag_in_stab?sold[U+jD]*problem->taup[iel]:0;
          local_vector[iphi+N_SF*(U+iD)] += supg*dphi[iphi][jD]*f[iD];
        }
      }
    }
    if (!local_matrix)
      continue;
    for (int iD = 0; iD < D; ++iD){
      for (int iphi = 0; iphi < N_SF; ++iphi){
        for (int jphi = 0; jphi < N_SF; ++jphi){
          LOCAL_MATRIX(iphi,jphi,U+iD,U+iD) += phi[iphi]*phi[jphi]*dfdu;
          LOCAL_MATRIX(iphi,jphi,U+iD,P) += phi[iphi]*dphi[jphi][iD]*dfddp;
          LOCAL_MATRIX(iphi,jphi,P,U+iD) += pspg*dphi[iphi][iD]*phi[jphi]*dfdu;
          LOCAL_MATRIX(iphi,jphi,P,P) += pspg*dphi[iphi][iD]*dphi[jphi][iD]*dfddp;
          for (int jD = 0; jD < D; ++jD) {
            double supg = drag_in_stab?sold[U+jD]*problem->taup[iel]:0;
            LOCAL_MATRIX(iphi,jphi,U+iD,U+iD) += supg*dphi[iphi][jD]*phi[jphi]*dfdu;
            LOCAL_MATRIX(iphi,jphi,U+iD,P) += supg*dphi[iphi][jD]*dphi[jphi][iD]*dfddp;
          }
        }
      }
    }
  }
  free(s);
  free(ds);
  free(sold);
}

static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n, double *f0,const double *s,const double *ds, const double *sold, const double *dsold, const double c, const double *dc, const double rho, const double mu, const double dt, int eid, const double *data, double *f00, double *f01)
{
  const int vid = wbnd->vid;
  const int pid = wbnd->pid;
  const int cid = wbnd->cid;
  const int aid = wbnd->aid;
  double p = s[P];
  const int n_fields = fluid_problem_n_fields(problem);
  double u[D], c_du_o_c[D][D], dp[D];
  double s_c = 0;
  double unold = 0;
  double unext = 0;
  double dcn = 0;
  for (int iD = 0; iD < D; ++iD) {
    u[iD] = s[U+iD];
    dp[iD] = ds[P*D+iD];
    unold += sold[U+iD]*n[iD];
    if(vid>=0)
    unext += data[vid+iD]*n[iD];
    dcn += dc[iD]*n[iD];
    s_c += vid<0?0:(u[iD]-data[vid+iD])*n[iD];
    for (int jD = 0; jD < D; ++jD)
      c_du_o_c[iD][jD] = ds[(U+iD)*D+jD] -u[iD]/c*dc[jD];
  }
  double h = problem->element_size[eid];
  double sigma = (1+D)/(D*h)*mu*N_N;
  if (wbnd->type == BND_WALL && pid >= 0) {
      f0[P] = -(s[P]-data[pid])/h*problem->taup[eid];
      f00[P*n_fields+P] += -1/h*problem->taup[eid];
  }
  if(wbnd->type == BND_OPEN) {
    double unext = 0;
    for (int i = 0; i < D; ++i) {
      unext += (vid<0?s[U+i]:data[vid+i])*n[i];
      f00[P*n_fields+U+i] += vid<0?n[i]:0;
    }
    f0[P] = unext;
  }
  for (int id = 0; id < D; ++id) {
    f0[U+id] = (pid<0?0:data[pid]-p)*n[id] + (vid<0?0:sigma*(u[id]-data[vid+id]+s_c*n[id]));
    f00[(U+id)*n_fields+P] += (pid<0?0:-n[id]);
    f00[(U+id)*n_fields+U+id] += (vid<0?0:sigma);
    for (int jd = 0; jd <D; ++jd) {
      f0[U+id] += mu*(c_du_o_c[id][jd]+c_du_o_c[jd][id])*n[jd];
      f00[(U+id)*n_fields+U+id] += (vid<0?0:n[id]*n[jd]*sigma) - mu/c*dc[jd]*n[jd];
      f00[(U+id)*n_fields+U+jd] -= mu/c*dc[id]*n[jd];
      f01[(U+id)*n_fields*D+(U+id)*D+jd] += mu*n[jd];
      f01[(U+id)*n_fields*D+(U+jd)*D+id] += mu*n[jd];
    }
  }
  double unbnd = unold;
   if (wbnd->type == BND_WALL) unbnd = unold/2;
   else if(vid>0) unbnd = (unold+unext)/2;

   if (unbnd<0) {
     for (int id = 0; id < D; ++id) {
       f0[U+id] += unbnd*((vid<0 ? 0 : data[vid+id])-u[id])*rho/c;
       f00[(U+id)*n_fields+U+id] -= unbnd*rho/c;
     }
   }
}

static double pow2(double a) {return a*a;}
static void compute_stab_parameters(FluidProblem *problem, double dt) {
  const Mesh *mesh = problem->mesh;

  problem->element_size = realloc(problem->element_size,sizeof(double)*mesh->n_elements);
  problem->taup = realloc(problem->taup,sizeof(double)*mesh->n_elements);
  problem->tauc = realloc(problem->tauc,sizeof(double)*mesh->n_elements);

  double maxtaup = 0;
  double mintaup = DBL_MAX;
  const int n_fields = fluid_problem_n_fields(problem);
  double maxtaa = 0;
  double mintaa = DBL_MAX;
  for (int iel = 0; iel < mesh->n_elements; ++iel) {
    double normu = {0.};
    const int *el = &mesh->elements[iel*N_N];
    double a = 0;
    double amax = 0.5;
    double amin = 0.5;
    double c = 0;
    for (int i=0; i< N_N; ++i) {
      double normup = 0;
      c += problem->porosity[el[i]]/N_N;
      if(problem->n_fluids == 2){
        a += problem->concentration[iel*N_N+i];
        amax = fmax(amax,problem->concentration[iel*N_N+i]);
        amin = fmin(amin,problem->concentration[iel*N_N+i]);
      }
      for(int j = 0; j < DIMENSION; ++j) normup += pow2(problem->solution[el[i]*n_fields+j]);
      normu = fmax(normu,sqrt(normup));
    }
    double dxidx[DIMENSION][DIMENSION];
    double detj = mesh_dxidx(mesh,iel,dxidx);
#if DIMENSION == 2
    const double h = 2*sqrt(detj/(2*M_PI));
#else
    const double h = 2*cbrt(detj/(8*M_PI));
#endif
    problem->element_size[iel] = h;
    double rho, mu;
    a /= N_N;
    a = fmin(1.,fmax(0.,a));
    fluid_problem_interpolate_rho_and_nu(problem,a, &rho, &mu);
    double nu = mu/rho;
    problem->taup[iel] = 1/sqrt(pow2(2/dt)+pow2(2*normu/h)+pow2(4*nu/pow2(h)));
    problem->tauc[iel] = h*normu*fmin(h*normu/(6*nu),0.5);
  }
}

#if D==2
static double particle_drag_coeff(double u[2], double mu, double rho, double vol, double c)
{
  double d = 2*sqrt(vol/M_PI);
  double normu = hypot(u[0],u[1]);
  //Reynolds/|u_p-u_f|
  double Re_O_u = d*c/mu*rho;
  double Cd_u = 0.63*sqrt(normu)+4.8/sqrt(Re_O_u);
  Cd_u = Cd_u*Cd_u;
  double f = pow(c,-(1.8-0.65*exp(-(1.5-log(Re_O_u*normu))*(1.5-log(Re_O_u*normu))/2.)));
  return f*Cd_u*rho/2*d;
}
#else
static double particle_drag_coeff(double u[3], double mu, double rho, double vol, double c)
{
  double d = 2*cbrt(3./4.*vol/M_PI);
  double normu = sqrt(u[0]*u[0]+u[1]*u[1]+u[2]*u[2]);
  //Reynolds/|u_p-u_f|
  double Re_O_u = d*c/mu*rho;
  double Cd_u = 0.63*sqrt(normu)+4.8/sqrt(Re_O_u);
  Cd_u = Cd_u*Cd_u;
  double f = pow(c,-1.8);
  double r = d/2.;
  return f*Cd_u*rho/2*(M_PI*r*r);
}
#endif

void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt, double *particle_force) {
  for (int i = 0; i < problem->n_particles*D; ++i) {
    particle_force[i] = problem->particle_force[i];
  }
}

static void particle_force_f(FluidProblem *problem, double *f, double *dfdu, double *dfddp, double *sol, double *dsol, const double *solold, const double c, const double *dc, double a, double dt, int iel, int ip) {
  const double *contact = &problem->particle_contact[ip*D];
  double u[D], uold[D], dp[D],G[D]={0};
  for (int iD = 0; iD < D; ++iD) {
    u[iD] = sol[U+iD];
    uold[iD] = solold[U+iD];
    dp[iD] = dsol[P*D+iD];
  }
  G[1] = problem->g;
  double mu, rho;
  fluid_problem_interpolate_rho_and_nu(problem,a, &rho, &mu);
  double Due[D];
  const double *up = &problem->particle_velocity[ip*D];
  for (int j = 0; j < D; ++j) {
    Due[j] = problem->particle_velocity[ip*D+j]-uold[j]/c;
  }
  double mass = problem->particle_mass[ip];
  double vol = problem->particle_volume[ip];
  double g = problem->g;
  if (problem->reduced_gravity){
    double rhop = mass/vol;
    g = g*(rhop-problem->rho[0])/rhop;
  }
  double gamma;
  
  if(problem->n_fluids == 2){
    double gamma0 = particle_drag_coeff(Due,problem->mu[0],problem->rho[0],vol,c);
    double gamma1 = particle_drag_coeff(Due,problem->mu[1],problem->rho[1],vol,c);
    gamma = gamma0*a + gamma1*(1-a);
  }
  else{
    gamma = particle_drag_coeff(Due,mu,rho,vol,c);
  }
  double gammastar = gamma/(1+dt/mass*gamma);
  for (int i = 0; i < D; ++i) {
    double upstar = up[i]+dt/mass*(contact[i]+mass*G[i]-vol*dp[i]);
    double drag = gammastar*(upstar-u[i]/c);
    problem->particle_force[ip*D+i] = -drag+G[i]*mass-vol*dp[i];
    f[U+i] = -drag-vol*dp[i];
  }
  *dfdu = gammastar/c;
  *dfddp = gammastar*dt/mass*vol-vol;
}

static void fluid_problem_f(const FluidProblem *problem, const double *sol, double *dsol, const double *solold, const double *dsolold, const double c, const double *dc, const double cold, const double rho, const double mu, double dt, int iel, double *f0, double *f1, double *f00, double *f10, double *f01, double *f11) {
  double p = sol[P];
  double taup = problem->taup[iel];
  double tauc = problem->tauc[iel];
  const int n_fields = fluid_problem_n_fields(problem);
  double u[D], uold[D], du[D][D], duold[D][D],dp[D];
  if(problem->stab_param != 0)
    taup = tauc = 0;
  double divu = 0;
  for (int iD = 0; iD < D; ++iD) {
    u[iD] = sol[U+iD];
    uold[iD] = solold[U+iD];
    dp[iD] = dsol[P*D+iD];
    for (int jD = 0; jD < D; ++jD) {
      du[iD][jD] = dsol[(U+iD)*D+jD];
      duold[iD][jD] = dsolold[(U+iD)*D+jD];
    }
    divu += du[iD][iD];
  }
  f0[P] = (c-cold)/dt;// + (sol[P]-solold[P])/dt*0.1;
  //f00[P*n_fields+P] = 1/dt*0.1;
  double G[D] = {0};
  G[1] = -problem->g;
  double rhoreduced = (problem->reduced_gravity ? (rho-problem->rho[0]) : rho);
  for (int i = 0; i < D; ++i) {
    f0[U+i] = 
      rho*(u[i]-uold[i])/dt 
      + dp[i] 
      + G[i]*rhoreduced*c
      + problem->volume_drag*u[i]*mu;
    f00[(U+i)*n_fields+U+i] = rho/dt + problem->volume_drag*mu;
    f01[(U+i)*n_fields*D+P*D+i] = 1;
    // advection :
    // dj(uj ui /c) = uj/c dj(ui) + ui/c dj(uj) - uj ui /(c*c) dj(c)
    //              = ui/c * (dj(uj)- uj/c*dj(c)) + uj/c*dj(ui)
    for (int j = 0; j < D; ++j) {
      f0[U+i] += rho/c*(uold[j]*du[i][j] + u[i]*(duold[j][j]-uold[j]/c*dc[j]));
      f00[(U+i)*n_fields+U+i] += rho/c*(duold[j][j]-uold[j]/c*dc[j]);
      f01[(U+i)*n_fields*D+(U+i)*D+j] += rho/c*uold[j];
    }
    for (int j = 0; j < D; ++j) {
      f1[(U+i)*D+j] = mu*(du[i][j]-u[i]*dc[j]/c+ du[j][i]-u[j]*dc[i]/c);
      f10[((U+i)*D+j)*n_fields+U+j] += -mu*dc[i]/c;
      f10[((U+i)*D+j)*n_fields+U+i] += -mu*dc[j]/c;
      f11[((U+i)*D+j)*n_fields*D+(U+j)*D+i] += mu; 
      f11[((U+i)*D+j)*n_fields*D+(U+i)*D+j] += mu;
      // SUPG
      double supg = uold[j]*taup;
      f1[(U+i)*D+j] += supg*f0[U+i];
      f10[((U+i)*D+j)*n_fields+U+i] += supg*f00[(U+i)*n_fields+U+i];
      f11[((U+i)*D+j)*(n_fields*D)+P*D+i] += supg*f01[(U+i)*n_fields*D+P*D+i]; 
      for (int k = 0; k < D; ++k)
        f11[((U+i)*D+j)*n_fields*D+(U+i)*D+k] += supg*f01[(U+i)*n_fields*D+(U+i)*D+k];
    }
    double lsic = tauc*rho;
    f1[(U+i)*D+i] += ((c-cold)/dt+divu)*lsic;
    for (int j = 0; j < D; ++j) {
      f11[((U+i)*D+i)*(n_fields*D)+(U+j)*D+j] += lsic;
    }
    f1[P*D+i] = -u[i];
    f10[(P*D+i)*n_fields+U+i] += -1;
    // PSPG
    double pspg = taup/rho;
    f1[P*D+i] += pspg*f0[U+i] + problem->stab_param*dp[i];
    f11[(P*D+i)*(n_fields*D)+P*D+i] += pspg*f01[(U+i)*n_fields*D+P*D+i]+problem->stab_param;
    f10[(P*D+i)*n_fields+U+i] += pspg*f00[(U+i)*n_fields+(U+i)];
    for (int j = 0; j < D; ++j) {
      f11[(P*D+i)*n_fields*D+(U+i)*D+j] = pspg*f01[(U+i)*n_fields*D+(U+i)*D+j];
    }
  }
}

static void compute_weak_boundary_conditions(FluidProblem *problem, const double *solution_old, double dt, double *all_local_vector, double *all_local_matrix)
{
  const Mesh *mesh = problem->mesh;
  const double *solution = problem->solution;
  const int n_fields = fluid_problem_n_fields(problem);
  const size_t local_size = N_SF*n_fields;
  double *s = malloc(sizeof(double)*n_fields);
  double *sold = malloc(sizeof(double)*n_fields);
  double *ds = malloc(sizeof(double)*n_fields*D);
  double *dsold = malloc(sizeof(double)*n_fields*D);
  double *f0 = malloc(sizeof(double)*n_fields);
  double *f00 = malloc(sizeof(double)*n_fields*n_fields);
  double *f01 = malloc(sizeof(double)*n_fields*n_fields*D);
  double i_bnd = 0;
  double iv_bnd_up = 0;
  double iv_bnd_down = 0;
  for (int ibnd = 0; ibnd < problem->n_weak_boundaries; ++ibnd){
    WeakBoundary *wbnd = problem->weak_boundaries + ibnd;
    int bndid= -1;
    for (int i = 0; i < mesh->n_boundaries; ++i) {
      if (strcmp(mesh->boundary_names[i],wbnd->tag) == 0){
        bndid = i;
      }
    }
    if (bndid == -1) {
      printf("Mesh has no boundary with name \"%s\".", wbnd->tag);
    }
    MeshBoundary *bnd = &problem->boundaries[bndid];
    int n_value = weak_boundary_n_values(wbnd);
    double *data = malloc(n_value*bnd->n_interfaces*N_LQP*sizeof(double));
    weak_boundary_values(mesh, bnd, wbnd, data);

    for (int iinterface = 0; iinterface < bnd->n_interfaces; ++iinterface) {
      const int *interface = &mesh->interfaces[bnd->interfaces[iinterface]*4];
      const int iel = interface[0];
      const int icl = interface[1];
      const int *cl = elbnd[icl];
      int nodes[N_LSF];
      double *x[N_LSF];
      const int *el = &mesh->elements[iel*N_N];
      for (int i = 0; i < D; ++i){
        nodes[i] = el[cl[i]];
        x[i] = &mesh->x[nodes[i]*3];
      }
      double dxidx[D][D], dphi[N_N][D];
      mesh_dxidx(mesh,iel,dxidx);
      double n[D],detbnd;
      get_normal_and_det(x,n,&detbnd);
      grad_shape_functions(dxidx, dphi);
      double dc[D] = {0};
      for (int j = 0; j < n_fields*D; ++j) {
        ds[j] = 0;
      }
      for (int i = 0; i < N_SF; ++i) {
        for (int j = 0; j < n_fields; ++j) {
          double dof = solution[el[i]*n_fields+j];
          double dofold = solution_old[el[i]*n_fields+j];
          for (int k = 0; k < D; ++k){
            ds[j*D+k] += dphi[i][k]*dof;
            dsold[j*D+k] += dphi[i][k]*dofold;
          }
        }
        for (int k=0; k<D; ++k){
          dc[k] += problem->porosity[el[i]]*dphi[i][k];
        }
      }
      for (int iqp = 0; iqp < N_LQP; ++iqp) {
        double *dataqp = &data[n_value*(N_LQP*iinterface+iqp)];
        double dp[D]={0};
        double *local_vector = &all_local_vector[local_size*iel];
        double *local_matrix = &all_local_matrix[local_size*local_size*iel];
        double phi[DIMENSION];
        l_shape_functions(LQP[iqp],phi);
        for (int i = 0; i < n_fields; ++i) {
          s[i] = 0;
          sold[i] = 0;
          f0[i] = 0;
        }
        for (int i = 0; i < n_fields*n_fields; ++i) {
          f00[i] = 0;
        }
        for (int i = 0; i < n_fields*n_fields*D; ++i) {
          f01[i] = 0;
        }
        double c = 0;
        double a = 0;
        for (int i = 0; i < DIMENSION; ++i) {
          c += problem->porosity[nodes[i]]*phi[i];
          if (problem->n_fluids == 2) {
            a += problem->concentration[iel*N_N+cl[i]]*phi[i];
          }
          for (int j = 0; j < n_fields; ++j) {
            double dof = solution[nodes[i]*n_fields+j];
            double dofold = solution_old[nodes[i]*n_fields+j];
            s[j] += phi[i]*dof;
            sold[j] += phi[i]*dofold;
          }
        }
        double rho, mu;
        fluid_problem_interpolate_rho_and_nu(problem,a, &rho, &mu);
        const double jw = LQW[iqp]*detbnd;
        if (wbnd->type != BND_WALL){
          for (int i = 0; i < D; ++i){
            if (wbnd->vid<0) {
              i_bnd -= a*s[U+i]*n[i]*jw*dt;
              iv_bnd_up -= s[U+i]*n[i]*jw;
            }
            else {
              iv_bnd_down -= dataqp[wbnd->vid+i]*n[i]*jw;
              i_bnd -= a*dataqp[wbnd->vid+i]*n[i]*jw*dt;
            }
          }
        }
        f_boundary(wbnd,problem,n,f0,s,ds,sold,dsold,c,dc,rho,mu,dt,iel,dataqp,f00,f01);
        for (int ifield = 0; ifield < n_fields; ++ifield) {
          for (int iphi = 0; iphi < D; ++iphi){
            local_vector[cl[iphi]+N_SF*ifield] += phi[iphi]*f0[ifield]*jw;
          }
        }
        for (int jfield = 0; jfield < n_fields; ++jfield) {
          for (int ifield = 0; ifield < n_fields; ++ifield){
            for (int iphi = 0; iphi < N_LSF; ++iphi){
              for (int jphi = 0; jphi < N_LSF; ++jphi){
                double d = phi[jphi]*f00[ifield*n_fields+jfield];
                LOCAL_MATRIX(cl[iphi],cl[jphi],ifield,jfield) += jw*phi[iphi]*d;
              }
              for (int jphi = 0; jphi < N_SF; ++jphi){
                double d = 0;
                for (int jD = 0; jD < D; ++jD) {
                  d += dphi[jphi][jD]*f01[ifield*n_fields*D+jfield*D+jD];
                }
                LOCAL_MATRIX(cl[iphi],jphi,ifield,jfield) += jw*phi[iphi]*d;
              }
            }
          }
        }
      }
    }
    free(data);
  }
  free(s);
  free(ds);
  free(sold);
  free(dsold);
  free(f0);
  free(f00);
}

double fluid_problem_a_integ_volume(FluidProblem *problem)
{
  const Mesh *mesh = problem->mesh;
  double I_a = 0;
  for (int iel=0; iel < mesh->n_elements; ++iel) {
    const int *el = &mesh->elements[iel*N_N];
    double dxidx[D][D];
    const double detj = mesh_dxidx(mesh,iel,dxidx);
    for (int i = 0; i < N_N; ++i)
      for (int j = 0; j< N_N; ++j) 
        I_a += detj * problem->porosity[el[i]]*problem->concentration[iel*N_N+j]*mass_matrix[i][j];
  }
  return I_a;
}

double fluid_problem_volume_flux(FluidProblem *problem, const char *tagname)
{
  const Mesh *mesh = problem->mesh;
  const double *solution = problem->solution;
  double Q = 0;
  double wtotal = 0;
  int n_fields = D+1;
  int found = 0;
  for (int i = 0; i < N_LQP; ++i) wtotal += LQW[i];
  for (int ibnd = 0; ibnd < mesh->n_boundaries; ++ibnd) {
    if (strcmp(mesh->boundary_names[ibnd],tagname) != 0)
      continue;
    found = 1;
    MeshBoundary *bnd = &problem->boundaries[ibnd];
    for (int iinterface = 0; iinterface < bnd->n_interfaces; ++iinterface) {
      const int *interface = &mesh->interfaces[bnd->interfaces[iinterface]*4];
      const int iel = interface[0];
      const int icl = interface[1];
      const int *cl = elbnd[icl];
      int nodes[N_LSF];
      double *x[N_LSF];
      const int *el = &mesh->elements[iel*N_N];
      for (int i = 0; i < D; ++i){
        nodes[i] = el[cl[i]];
        x[i] = &mesh->x[nodes[i]*3];
      }
      double dxidx[D][D], dphi[N_N][D];
      mesh_dxidx(mesh,iel,dxidx);
      double n[D],detbnd;
      get_normal_and_det(x,n,&detbnd);
      double vnmean = 0;
      for (int i = 0; i < D; ++i) {
        for (int iD = 0; iD < D; ++iD) {
          vnmean += n[iD]*problem->solution[nodes[i]*n_fields+U+iD];
        }
      }
      Q += vnmean*detbnd*wtotal/D;
    }
  }
  if(found == 0) {
    printf("boundary '%s' not found\n", tagname);
    exit(1);
  }
  printf("Q = %g wtotal = %g\n", Q,wtotal);
  return Q;
}

double fluid_problem_a_integ_bnd(FluidProblem *problem, double dt)
{
  const Mesh *mesh = problem->mesh;
  const double *solution = problem->solution;
  const int n_fields = fluid_problem_n_fields(problem);
  const size_t local_size = N_SF*n_fields;
  double *s = malloc(sizeof(double)*n_fields);
  double i_bnd = 0;
  for (int ibnd = 0; ibnd < problem->n_weak_boundaries; ++ibnd){
    WeakBoundary *wbnd = problem->weak_boundaries + ibnd;
    int bndid= -1;
    for (int i = 0; i < mesh->n_boundaries; ++i) {
      if (strcmp(mesh->boundary_names[i],wbnd->tag) == 0){
        bndid = i;
      }
    }
    if (bndid == -1) {
      printf("Mesh has no boundary with name \"%s\".", wbnd->tag);
    }
    MeshBoundary *bnd = &problem->boundaries[bndid];
    int n_value = weak_boundary_n_values(wbnd);
    double *data = malloc(n_value*bnd->n_interfaces*N_LQP*sizeof(double));
    weak_boundary_values(mesh, bnd, wbnd, data);

    for (int iinterface = 0; iinterface < bnd->n_interfaces; ++iinterface) {
      const int *interface = &mesh->interfaces[bnd->interfaces[iinterface]*4];
      const int iel = interface[0];
      const int icl = interface[1];
      const int *cl = elbnd[icl];
      int nodes[N_LSF];
      double *x[N_LSF];
      const int *el = &mesh->elements[iel*N_N];
      for (int i = 0; i < D; ++i){
        nodes[i] = el[cl[i]];
        x[i] = &mesh->x[nodes[i]*3];
      }
      double dxidx[D][D], dphi[N_N][D];
      mesh_dxidx(mesh,iel,dxidx);
      double n[D],detbnd;
      get_normal_and_det(x,n,&detbnd);
      grad_shape_functions(dxidx, dphi);
      
      for (int iqp = 0; iqp < N_LQP; ++iqp) {
        double *dataqp = &data[n_value*(N_LQP*iinterface+iqp)];
        double phi[DIMENSION];
        l_shape_functions(LQP[iqp],phi);
        for (int i = 0; i < n_fields; ++i) {
          s[i] = 0;
        }
        double a = 0;
        for (int i = 0; i < DIMENSION; ++i) {
          if (problem->n_fluids == 2) {
            a += problem->concentration[iel*N_N+cl[i]]*phi[i];
          }
          for (int j = 0; j < n_fields; ++j) {
            double dof = solution[nodes[i]*n_fields+j];
            s[j] += phi[i]*dof;
          }
        }
        const double jw = LQW[iqp]*detbnd;
        if (wbnd->type != BND_WALL){
          for (int i = 0; i < D; ++i){
            if (wbnd->vid<0) {
              i_bnd -= a*s[U+i]*n[i]*jw*dt;
            }
            else {
              i_bnd -= a*dataqp[wbnd->vid+i]*n[i]*jw*dt;
            }
          }
        }
      }
    }
  }
  return i_bnd;
}

static void fluid_problem_volume(FluidProblem *problem, const double *solution_old, double dt, double *all_local_vector, double *all_local_matrix)
{
  const Mesh *mesh = problem->mesh;
  const double *porosity = problem->porosity;

  const double *solution = problem->solution;
  int n_fields = fluid_problem_n_fields(problem);
  size_t local_size = N_SF*n_fields;
  double *s = malloc(sizeof(double)*n_fields);
  double *sold = malloc(sizeof(double)*n_fields);
  double *ds = malloc(sizeof(double)*n_fields*D);
  double *dsold = malloc(sizeof(double)*n_fields*D);
  double *f0 = malloc(sizeof(double)*n_fields);
  double *f1 = malloc(sizeof(double)*n_fields*D);
  double *f00 = malloc(sizeof(double)*n_fields*n_fields);
  double *f10 = malloc(sizeof(double)*n_fields*n_fields*D);
  double *f01 = malloc(sizeof(double)*n_fields*n_fields*D);
  double *f11 = malloc(sizeof(double)*n_fields*n_fields*D*D);
  problem->kinetic_energy = 0;
  for (int iel=0; iel < mesh->n_elements; ++iel) {
    const int *el = &mesh->elements[iel*N_N];
    double dxidx[D][D], dphi[N_N][D];
    const double detj = mesh_dxidx(mesh,iel,dxidx);
    grad_shape_functions(dxidx, dphi);
    double *local_vector = &all_local_vector[local_size*iel];
    double *local_matrix = &all_local_matrix[local_size*local_size*iel];
    for (int i = 0; i < n_fields; ++i) {
      for (int j = 0; j < D; ++j) {
        ds[i*D+j] = 0;
        dsold[i*D+j] = 0;
      }
    }
    double dc[D] = {0}, da[D]={0};
    for (int i = 0; i < N_SF; ++i) {
      for (int j = 0; j < n_fields; ++j) {
        double dof = solution[el[i]*n_fields+j];
        double dofold = solution_old[el[i]*n_fields+j];
        for (int k = 0; k < D; ++k) {
          ds[j*D+k] += dphi[i][k]*dof;
          dsold[j*D+k] += dphi[i][k]*dofold;
        }
      }
      for (int k=0; k<D; ++k){
        dc[k] += problem->porosity[el[i]]*dphi[i][k];
      }
      if (problem->n_fluids==2) {
        for (int k=0; k<D; ++k){
          da[k] += problem->concentration[iel*N_N+i]*dphi[i][k];
        }
      }
    }
    for (int iqp = 0; iqp< N_QP; ++iqp) {
      double phi[N_SF];
      shape_functions(QP[iqp],phi);
      for (int i = 0; i < n_fields; ++i) {
        s[i] = 0;
        sold[i] = 0;
      }
      for (int i = 0; i < n_fields*n_fields; ++i) {
        f00[i] = 0;
      }
      for (int i = 0; i < D*n_fields*n_fields; ++i) {
        f10[i] = 0;
        f01[i] = 0;
      }
      for (int i = 0; i < D*D*n_fields*n_fields; ++i) {
        f11[i] = 0;
      }
      double c = 0, a = 0;
      double cold = 0, aold = 0;
      for (int i = 0; i < N_SF; ++i) {
        c += problem->porosity[el[i]]*phi[i];
        cold += problem->oldporosity[el[i]]*phi[i];
        for (int j = 0; j < n_fields; ++j) {
          double dof = solution[el[i]*n_fields+j];
          double dofold = solution_old[el[i]*n_fields+j];
          s[j] += phi[i]*dof;
          sold[j] += phi[i]*dofold;
        }
        if (problem->n_fluids==2) {
          a += problem->concentration[iel*N_N+i]*phi[i];
        }
      }
      double mu;
      double rho;
      fluid_problem_interpolate_rho_and_nu(problem,a, &rho, &mu);
      const double jw = QW[iqp]*detj;
      problem->kinetic_energy += rho*sqrt(s[U]*s[U]+s[U+1]*s[U+1])*sqrt(s[U]*s[U]+s[U+1]*s[U+1])/2*jw;

      fluid_problem_f(problem,s,ds,sold,dsold,c,dc,cold,rho,mu,dt,iel,f0,f1,f00,f10,f01,f11);
      for (int ifield = 0; ifield < n_fields; ++ifield) {
        for (int iphi = 0; iphi < N_SF; ++iphi){
          local_vector[iphi+N_SF*ifield] += phi[iphi]*f0[ifield]*jw;
          for (int id = 0; id < D; ++id) {
            local_vector[iphi+N_SF*ifield] += dphi[iphi][id]*f1[ifield*D+id]*jw;
          }
        }
      }
      for (int jfield = 0; jfield < n_fields; ++jfield) {
        for (int ifield = 0; ifield < n_fields; ++ifield){
          for (int iphi = 0; iphi < N_SF; ++iphi){
            for (int jphi = 0; jphi < N_SF; ++jphi){
              double m = jw*phi[iphi]*phi[jphi]*f00[ifield*n_fields+jfield];
              for (int id = 0; id < D; ++id) {
                m += jw*dphi[iphi][id]*phi[jphi]*f10[(ifield*D+id)*n_fields+jfield];
                for (int jd = 0; jd < D; ++jd) {
                  m += jw*dphi[iphi][id]*dphi[jphi][jd]*f11[(ifield*D+id)*n_fields*D+jfield*D+jd];
                }
              }
              for (int jd = 0; jd < D; ++jd) {
                m += jw*phi[iphi]*dphi[jphi][jd]*f01[ifield*(n_fields*D)+jfield*D+jd];
              }
              LOCAL_MATRIX(iphi,jphi,ifield,jfield) += m;
            }
          }
        }
      }
    }
  }
  free(s);
  free(ds);
  free(dsold);
  free(sold);
  free(f0);
  free(f1);
  free(f00);
  free(f01);
  free(f10);
  free(f11);
}

static void fluid_problem_surface_tension_bnd(FluidProblem *problem, double *a_cg, double *all_local_vector){
  const Mesh *mesh = problem->mesh;
  double sigma = problem->sigma;
  const int n_fields = fluid_problem_n_fields(problem);
  const size_t local_size = N_SF*n_fields;
  for (int ibnd = 0; ibnd < problem->n_weak_boundaries; ++ibnd){
    WeakBoundary *wbnd = problem->weak_boundaries + ibnd;
    int bndid= -1;
    for (int i = 0; i < mesh->n_boundaries; ++i) {
      if (strcmp(mesh->boundary_names[i],wbnd->tag) == 0){
        bndid = i;
      }
    }
    if (bndid == -1) {
      printf("Mesh has no boundary with name \"%s\".", wbnd->tag);
    }
    MeshBoundary *bnd = &problem->boundaries[bndid];
    int n_value = weak_boundary_n_values(wbnd);
    for (int iinterface = 0; iinterface < bnd->n_interfaces; ++iinterface) {
      const int *interface = &mesh->interfaces[bnd->interfaces[iinterface]*4];
      const int iel = interface[0];
      const int icl = interface[1];
      const int *cl = elbnd[icl];
      int nodes[D];
      double *x[D];
      const int *el = &mesh->elements[iel*N_N];
      for (int i = 0; i < D; ++i){
        nodes[i] = el[elbnd[icl][i]];
        x[i] = &mesh->x[nodes[i]*3];
      }
      double dxidx[D][D], dphi[N_N][D];
      mesh_dxidx(mesh,iel,dxidx);
      double n[D],detbnd;
      get_normal_and_det(x,n,&detbnd);
      grad_shape_functions(dxidx, dphi);
      double da[D] = {0};
      double nda = 0;
      for (int k=0; k<D; ++k){
        for (int i = 0; i < N_SF; ++i) {
          da[k] += a_cg[el[i]]*dphi[i][k];
        }
        nda += da[k]*da[k];
      }
      nda = fmax(sqrt(nda),1e-8);
      for (int iqp = 0; iqp < N_LQP; ++iqp) {
        double *local_vector = &all_local_vector[local_size*iel];
        double phi[DIMENSION];
        l_shape_functions(LQP[iqp],phi);
        const double jw = LQW[iqp]*detbnd;
        double a = 0;
        double c = 0;
        for (int i = 0; i < D; ++i){
          a += a_cg[nodes[i]]*phi[i];
          c += problem->porosity[nodes[i]]*phi[i];
        }
        for (int iphi = 0; iphi < D; ++iphi) {
          for (int iD = 0; iD < D; ++iD) {
            for (int jD = 0; jD < D; ++jD){
              local_vector[cl[iphi]+N_SF*iD] += phi[iphi]*jw*c*da[iD]*da[jD]*n[jD]/nda*sigma;
            }
            local_vector[cl[iphi]+N_SF*iD] -= phi[iphi]*jw*c*nda*n[iD]*sigma;
          }
        }
      }
    }
  }
}

static void fluid_problem_surface_tension(FluidProblem *problem, double *a_cg, double *all_local_vector){
  if (problem->n_fluids == 1) return;
  const Mesh *mesh = problem->mesh;
  double sigma = problem->sigma;
  int n_fields = fluid_problem_n_fields(problem);
  for (int iel = 0; iel < mesh->n_elements; ++iel) {
    const int *el = &mesh->elements[iel*N_N];
    double dxidx[D][D], dphi[N_N][D];
    const double det = mesh_dxidx(mesh,iel,dxidx);
    const double vol = element_volume_from_detj(det);
    grad_shape_functions(dxidx, dphi);
    double da[D] = {0};
    double nda = 0;
    double c = 0;
    for (int i = 0; i < N_N; ++i){
      c += problem->porosity[el[i]]/N_N;
    }
    for (int k=0; k<D; ++k){
      for (int i = 0; i < N_SF; ++i) {
        //da[k] += a_cg[el[i]]*dphi[i][k];
        da[k] += problem->concentration[iel*N_N+i]*dphi[i][k];
      }
      nda += da[k]*da[k];
    }
    nda = fmax(sqrt(nda),1e-8);
    double *local_vector = &all_local_vector[iel*n_fields*N_N];
    for (int iphi = 0; iphi < N_N; ++iphi) {
      for (int id = 0; id < D; ++id) {
        for (int jd = 0; jd < D; ++jd) {
          local_vector[(U+id)*N_N+iphi] -= c*vol*da[id]*da[jd]*dphi[iphi][jd]/nda*sigma;
        }
        local_vector[(U+id)*N_N+iphi] += c*vol*nda*dphi[iphi][id]*sigma;
      }
    }
  }
}


static void fluid_problem_dg_to_cg(FluidProblem *problem, const double *dg, double *cg) {
  const Mesh *mesh = problem->mesh;
  for (int i = 0; i< mesh->n_nodes; ++i) {
    cg[i] = 0;
  }
  for (int iel = 0; iel < mesh->n_elements; ++iel) {
    const int *el = &mesh->elements[iel*N_N];
    double dxidx[D][D];
    const double det = mesh_dxidx(mesh,iel,dxidx);
    const double vol = node_volume_from_detj(det);
    for (int i = 0; i< N_N; ++i) {
      cg[el[i]] += vol * dg[iel*N_N+i];
    }
  }
  for (int i = 0; i< mesh->n_nodes; ++i) {
    cg[i] /= problem->node_volume[i];
  }
}

static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, const double *solution_old, double dt)
{
  HXTLinearSystem *lsys = problem->linear_system;
  const Mesh *mesh = problem->mesh;

  int n_fields = fluid_problem_n_fields(problem);
  
  size_t local_size = N_SF*n_fields;
  size_t N_E = mesh->n_elements;
  double *all_local_vector = (double*) malloc(N_E*local_size*sizeof(double));
  double *all_local_matrix = (double*) malloc(N_E*local_size*local_size*sizeof(double));
  for (size_t i = 0; i < N_E*local_size; ++i)
    all_local_vector[i] = 0;
  for (size_t i = 0; i < N_E*local_size*local_size; ++i)
    all_local_matrix[i] = 0;

  char *forced_row = malloc(sizeof(char)*n_fields*mesh->n_nodes);
  for (int i = 0; i < n_fields*mesh->n_nodes; ++i)
    forced_row[i] = -1;
  for (int ibnd = 0; ibnd < problem->n_strong_boundaries; ++ibnd) {
    const StrongBoundary *bnd = problem->strong_boundaries + ibnd;
    int iphys;
    for (iphys = 0; iphys < mesh->n_boundaries; ++iphys) {
      if (strcmp(bnd->tag, mesh->boundary_names[iphys]) == 0)
        break;
    }
    if (iphys == mesh->n_boundaries){
      printf("Boundary tag \"%s\" not found.\n", bnd->tag);
      exit(1);
    }
    for (int i = 0; i < problem->boundaries[iphys].n_nodes; ++i){
      forced_row[problem->boundaries[iphys].nodes[i]*n_fields+bnd->row] = bnd->field;
    }
  }
  node_force_volume(problem, solution_old, dt, all_local_vector, all_local_matrix);
  compute_weak_boundary_conditions(problem, solution_old, dt, all_local_vector, all_local_matrix);
  fluid_problem_volume(problem, solution_old, dt, all_local_vector, all_local_matrix);
  if (problem->n_fluids == 2){
    double *a_cg = malloc(sizeof(double)*mesh->n_nodes);
    fluid_problem_dg_to_cg(problem, problem->concentration, a_cg);
    fluid_problem_surface_tension(problem, a_cg, all_local_vector);
    fluid_problem_surface_tension_bnd(problem, a_cg, all_local_vector);
    free(a_cg);
  }
  for (int iel=0; iel < mesh->n_elements; ++iel) {
    double *local_vector = &all_local_vector[local_size*iel];
    double *local_matrix = &all_local_matrix[local_size*local_size*iel];
    for (int inode = 0; inode < N_SF; ++inode) {
      for(int ifield = 0; ifield < n_fields; ++ifield) {
        const int *el = &mesh->elements[iel*N_N];
        char forced_field = forced_row[el[inode]*n_fields+ifield];
        if (forced_field == -1) continue;
        int i = inode*n_fields + ifield;
        for (int jnode = 0; jnode< N_SF; ++jnode){
          for (int jfield= 0; jfield< n_fields; ++jfield){
            local_matrix[(inode*n_fields+ifield)*local_size+(jnode*n_fields+jfield)] =
              (inode==jnode && jfield == forced_field) ? 1. : 0.;
          }
        }
        local_vector[inode+ifield*N_SF] = 0;
      }
    }
    hxtLinearSystemAddToMatrix(lsys,iel,iel,local_matrix);
    hxtLinearSystemAddToRhs(lsys,rhs,iel,local_vector);
  }
  free(forced_row);
  free(all_local_matrix);
  free(all_local_vector);
}

static void apply_boundary_conditions(FluidProblem *problem)
{
  const Mesh *mesh = problem->mesh;
  int n_fields = fluid_problem_n_fields(problem);
  for (int ibnd = 0; ibnd < problem->n_strong_boundaries; ++ibnd) {
    StrongBoundary *bnd = problem->strong_boundaries + ibnd;
    int iphys;
    for (iphys = 0; iphys < mesh->n_boundaries; ++iphys) {
      if (strcmp(bnd->tag, mesh->boundary_names[iphys]) == 0)
        break;
    }
    if (iphys == mesh->n_boundaries){
      printf("Boundary tag \"%s\" not found.\n", bnd->tag);
      exit(1);
    }
    MeshBoundary *mbnd = &problem->boundaries[iphys];
    double *x = (double*)malloc(mbnd->n_nodes*sizeof(double)*D);
    double *v = (double*)malloc(mbnd->n_nodes*sizeof(double));
    for (int i = 0; i < mbnd->n_nodes; ++i){
      int j = mbnd->nodes[i];
      for (int k = 0; k < D; ++k)
        x[i*D+k] = mesh->x[j*3+k];
    }
    bnd->apply(mbnd->n_nodes, x, v);
    for (int i = 0; i < mbnd->n_nodes; ++i)
      problem->solution[mbnd->nodes[i]*n_fields+bnd->field] = v[i];
    free(x);
    free(v);
  }
}


int fluid_problem_implicit_euler(FluidProblem *problem, double dt, double newton_atol, double newton_rtol, int newton_max_it)
{
  static double totaltime = 0;
  struct timespec tic, toc;
  const Mesh *mesh = problem->mesh;
  int n_fields = fluid_problem_n_fields(problem);
  apply_boundary_conditions(problem);
  double *solution_old = (double*)malloc(mesh->n_nodes*n_fields*sizeof(double));
  double *rhs = (double*)malloc(mesh->n_nodes*n_fields*sizeof(double));
  for (int i=0; i<mesh->n_nodes*n_fields; ++i)
    solution_old[i] = problem->solution[i];
  problem->solution_explicit = solution_old;
  double *dx = (double*)malloc(sizeof(double)*mesh->n_nodes*n_fields);
  double norm0 = 0;
  int i;
  compute_stab_parameters(problem,dt);
  for (i=0; i<newton_max_it; ++i) {
    double norm = 0;
    hxtLinearSystemZeroMatrix(problem->linear_system);

    for (int i=0; i<mesh->n_nodes*n_fields; ++i){
      rhs[i] = 0;
    }
    fluid_problem_assemble_system(problem, rhs, solution_old, dt);
    hxtLinearSystemGetRhsNorm(problem->linear_system,rhs,&norm);
    printf("iter %i norm = %.16g\n", i, norm);
    if (norm < newton_atol)
      break;
    if (i == 0)
      norm0 = norm;
    else if(norm/norm0 < newton_rtol)
      break;
    clock_get(&tic);
    hxtLinearSystemSolve(problem->linear_system, rhs, dx);
    clock_get(&toc);
    totaltime += clock_diff(&tic, &toc);
    for (int j=0; j<mesh->n_nodes*n_fields; ++j) {
      problem->solution[j] -= dx[j];
    }
    node_force_volume(problem,problem->solution_explicit,dt,NULL,NULL);
    break; /// equation is linear
  }
  
  printf("total solve time : %g\n", totaltime);
  free(dx);
  free(rhs);
  free(solution_old);
  return i;
}

static void fluid_problem_compute_node_volume(FluidProblem *problem) {
  free(problem->node_volume);
  Mesh *mesh = problem->mesh;
  problem->node_volume = (double*)malloc(sizeof(double)*mesh->n_nodes);
  for (int i = 0; i < problem->mesh->n_nodes; ++i){
    problem->node_volume[i] = 0;
  }
  for (int iel = 0; iel < mesh->n_elements; ++iel) {
    double dxidx[D][D];
    const int *el = &mesh->elements[iel*N_N];
    const double detj = mesh_dxidx(mesh, iel, dxidx);
    const double vol = node_volume_from_detj(detj);
    for (int i = 0; i < N_N; ++i)
      problem->node_volume[el[i]] += vol;
  }
}

static void fluid_problem_compute_porosity(FluidProblem *problem)
{
  Mesh *mesh = problem->mesh;
  double *volume = problem->particle_volume;
  for (int i = 0; i < mesh->n_nodes; ++i){
    problem->porosity[i] = 0;
  }
  for (int i = 0; i < problem->n_particles; ++i) {
    if(problem->particle_element_id[i] == -1) continue;
    double sf[N_SF];
    shape_functions(&problem->particle_uvw[i*D],sf);
    const int *el = &mesh->elements[problem->particle_element_id[i]*N_N];
    for (int j=0; j<N_N;++j)
      problem->porosity[el[j]] += sf[j]*volume[i];
  }
  for (int i = 0; i < mesh->n_nodes; ++i){
    if(problem->node_volume[i]==0.){
      problem->porosity[i] = 0.;
    }
    else{
      problem->porosity[i] = 1-problem->porosity[i]/problem->node_volume[i];
    }
  }
}

FluidProblem *fluid_problem_new(double g, int n_fluids, double *mu, double *rho, double sigma, double coeffStab, double volume_drag, const char *petsc_solver_type, int drag_in_stab)
{
  if (n_fluids != 1 && n_fluids != 2) {
    printf("error : only 1 or 2 fluids are supported\n");
  }
  static int initialized = 0;
  if (!initialized) {
    hxtInitializeLinearSystems(0, NULL);
    initialized = 1;
#ifdef HXT_HAVE_PETSC
    hxtPETScInsertOptions(petsc_solver_type, "fluid");
#endif
  }
  FluidProblem *problem = (FluidProblem*)malloc(sizeof(FluidProblem));
  problem->n_fluids = n_fluids;
  problem->sigma = sigma;
  problem->stab_param = 0;
  problem->drag_in_stab = drag_in_stab;
  problem->reduced_gravity = 0;
  problem->strong_boundaries = NULL;
  problem->n_strong_boundaries = 0;
  problem->weak_boundaries = NULL;
  problem->n_weak_boundaries = 0;
  problem->mu = (double*)malloc(sizeof(double)*problem->n_fluids);
  problem->rho = (double*)malloc(sizeof(double)*problem->n_fluids);
  for (int i = 0; i < problem->n_fluids; ++i){
    problem->rho[i] = rho[i];
    problem->mu[i] = mu[i];
  }
  problem->coeffStab = coeffStab;
  problem->g = g;
  problem->node_volume = NULL;
  problem->mesh_tree = NULL;
  problem->mesh = NULL;
  problem->taup = NULL;
  problem->tauc = NULL;
  problem->element_size = NULL;
  problem->porosity = NULL;
  problem->oldporosity = NULL;
  problem->solution = NULL;
  problem->concentration = NULL;
  problem->linear_system = NULL;
  problem->n_particles = 0;
  problem->particle_uvw = NULL;
  problem->particle_element_id = NULL;
  problem->particle_position = NULL;
  problem->particle_mass = NULL;
  problem->particle_force = NULL;
  problem->particle_volume = NULL;
  problem->particle_velocity = NULL;
  problem->particle_contact = NULL;
  problem->volume_drag = volume_drag;
  problem->boundaries = NULL;
  return problem;
}

void fluid_problem_set_stab_param(FluidProblem *p, double stab_param) {
  p->stab_param = stab_param;
}

void fluid_problem_set_reduced_gravity(FluidProblem *p, int reduced_gravity) {
  p->reduced_gravity = reduced_gravity;
}

void fluid_problem_free(FluidProblem *problem)
{
  free(problem->mu);
  free(problem->rho);
  free(problem->solution);
  free(problem->concentration);
  free(problem->porosity);
  free(problem->oldporosity);
  free(problem->node_volume);
  hxtLinearSystemDelete(&problem->linear_system);
  free(problem->particle_uvw);
  free(problem->particle_element_id);
  free(problem->particle_position);
  free(problem->particle_mass);
  free(problem->particle_force);
  free(problem->particle_volume);
  free(problem->particle_velocity);
  free(problem->particle_contact);
  free(problem->taup);
  free(problem->tauc);
  free(problem->element_size);
  for (int i = 0; i < problem->n_strong_boundaries; ++i)
    free(problem->strong_boundaries[i].tag);
  for (int i = 0; i < problem->n_weak_boundaries; ++i)
    free(problem->weak_boundaries[i].tag);
  free(problem->weak_boundaries);
  free(problem->strong_boundaries);
  mesh_tree_free(problem->mesh_tree);
  mesh_free(problem->mesh);
  mesh_boundaries_free(problem->boundaries,problem->mesh->n_boundaries);
  free(problem);
}

double fluid_problem_divergence_measure(FluidProblem *problem){
    double divu_max=0.0; double divu_mean=0.0;
    Mesh *mesh = problem->mesh;
    double *solution = problem->solution;
    int n_fields = fluid_problem_n_fields(problem);
    double *porosity = problem->porosity;
    
    for(int iel = 0; iel<mesh->n_elements; iel++){
        const unsigned int *el = &mesh->elements[iel*N_N];
        double C[N_N],Un[D][N_N];
        
        for (int i = 0; i < N_N; ++i){ // the for-loop is done over the N_N nodes of the elements iel
            C[i] = porosity[el[i]];
            for (int j=0; j<D; ++j) {
                Un[j][i] = solution[el[i]*n_fields+j]; // extract the velocity
            }
        }
        
        double dxdxi[D][D],dxidx[D][D];
        // Matrix for the linear mapping between the refence triangular element and any triangular elements of the mesh. This is also true in 3D.
        for (int i=0; i<D; ++i)
            for (int j=0; j<D; ++j)
                dxdxi[i][j] = mesh->x[el[j+1]*3+i]-mesh->x[el[0]*3+i];
    
        invDD(dxdxi, dxidx);
        double dphi[N_SF][D];
        grad_shape_functions(dxidx, dphi); // compute the gradient of the shape functions
        double divu = 0;
        for (int i=0; i<D; ++i){
            divu += dphi[i][i]*Un[i][i]/C[i]; // compute the divergence over this particular element
        }
        divu_max = fmax(divu_max, fabs(divu));
        divu_mean += divu;
    }
    divu_mean = divu_mean/mesh->n_elements;
    printf("The mean divergence is %.4e\n", divu_mean);
    printf("The max divergence is %.4e\n", divu_max);
    return divu_mean;
}
/*
 The function gradient_recovery_estimator aims to compute the average gardient jump over each edge in order to later on compute the linear interpolation.
 The result is stored in the two vectors gradient_x and gradient_y.
 
 @inputs: - problem
 - edges_vec  ==> an array of Edge structure (see mesh.h for more information about this structure)
 - n_edges    ==> an interger given the number of edges in the mesh following the Euler's formula
 - gradient_x ==> an array of length 3*mesh->n_elements
 - gradient_y ==> an array of length 3*mesh->n_elements
 
 @ouputs: the function returns nothing but it write in arrays gradient_x and gradient_y
 */
void gradient_recovery_estimator(FluidProblem *problem, double *gradient_x, double *gradient_y, int bool_pressure, int bool_velocity_U, int bool_velocity_V){
    if (bool_velocity_U) {
        printf("-------------------- START computing the jump in the gradient of U-velocity field --------------------\n");
    }else if(bool_pressure){
        printf("-------------------- START computing the jump in the gradient of pressure field --------------------\n");
    }else if(bool_velocity_V){
        printf("-------------------- START computing the jump in the gradient of V-velocity field --------------------\n");
    }else{
        printf("Error \n");
        exit(1);
    }
    
    Mesh *mesh = problem->mesh;
    double *porosity = problem->porosity;
    double *solution = problem->solution;
    int n_edges = mesh->n_interfaces;
    
    // Computing the jump in x and in y over each edge
    double Pn1[N_N], Pn2[N_N];
    double C1[N_N], C2[N_N];
    double Un1[D][N_N], Un2[D][N_N];
    int n_fields = fluid_problem_n_fields(problem);
    
    // Initialized a vector that count the edges for every tiangle
    // Each time we write into the output gradient_x and gradient_y vector
    // one of the entry is increment. In the end, kk should be 3 everywhere.
    int *kk = (int*) calloc(mesh->n_elements, sizeof(int));
    
    for (int ii=0; ii<n_edges; ii++){
        int *interface = &mesh->interfaces[ii*4];
        int el0 = interface[0];
        int cl0 = interface[1];
        int el1 = interface[2];
        int cl1 = interface[3];
        
        // The jump in the gradient is only computed for inner triangle
        if (el1>-1) {
            const int *el_T1 = &mesh->elements[el0*N_N];
            const int *el_T2 = &mesh->elements[el1*N_N];
            
            // Extracting pressure, porosity and velocity field from the solution
            for (int i=0;i<N_N; i++) {
                // Porosity field
                C1[i] = porosity[el_T1[i]];
                C2[i] = porosity[el_T2[i]];
                
                //Pressure field
                Pn1[i] = solution[el_T1[i]*n_fields+D]-problem->g*problem->rho[0]*mesh->x[el_T1[i]*3+1];
                Pn2[i] = solution[el_T2[i]*n_fields+D]-problem->g*problem->rho[0]*mesh->x[el_T2[i]*3+1];
                
                //Velocity field
                for (int jj = 0; jj < D; ++jj) {
                    Un1[jj][i] = solution[el_T1[i]*n_fields+jj];
                    Un2[jj][i] = solution[el_T2[i]*n_fields+jj];
                }

            }
            
            double dxdxi_T1[D][D], dxidx_T1[D][D];
            double dxdxi_T2[D][D], dxidx_T2[D][D];
            for(int i=0; i<D; i++){
                for(int j=0; j<D; j++){
                    dxdxi_T1[i][j] = mesh->x[el_T1[j+1]*3+i]-mesh->x[el_T1[0]*3+i];
                    dxdxi_T2[i][j] = mesh->x[el_T2[j+1]*3+i]-mesh->x[el_T2[0]*3+i];
                }
            }
            const double vol_T1 = detDD(dxdxi_T1)/2;
            const double vol_T2 = detDD(dxdxi_T2)/2;
            
            // Computing the pressure gradient, the porosity gradient and the velocity gradient
            invDD(dxdxi_T1, dxidx_T1); invDD(dxdxi_T2, dxidx_T2);
            double dphi_T1[N_SF][D], dphi_T2[N_SF][D];
            grad_shape_functions(dxidx_T1, dphi_T1); grad_shape_functions(dxidx_T2, dphi_T2);
            double dPn1dx = 0.0, dPn1dy = 0.0, dPn2dx = 0.0, dPn2dy = 0.0;
            //double dC1dx = 0.0, dC1dy = 0.0, dC2dx = 0.0, dC2dy = 0.0;
            double dU1dx = 0.0, dU1dy = 0.0, dU2dx = 0.0, dU2dy = 0.0;
            double dV1dx = 0.0, dV1dy = 0.0, dV2dx = 0.0, dV2dy = 0.0;
            
            for (int k=0; k<N_SF; ++k){
                // Porosity gradient
                //dC1dx += dphi_T1[k][0]*C1[k]*vol_T1;
                //dC1dy += dphi_T1[k][1]*C1[k]*vol_T1;
                //dC2dx += dphi_T2[k][0]*C2[k]*vol_T2;
                //dC2dy += dphi_T2[k][1]*C2[k]*vol_T2;
                
                // Pressure gradient
                dPn1dx += dphi_T1[k][0]*Pn1[k]*vol_T1;
                dPn1dy += dphi_T1[k][1]*Pn1[k]*vol_T1;
                dPn2dx += dphi_T2[k][0]*Pn2[k]*vol_T2;
                dPn2dy += dphi_T2[k][1]*Pn2[k]*vol_T2;
                
                // Velocity gradient
                dU1dx += dphi_T1[k][0]*Un1[0][k]/C1[k]*vol_T1; dU1dy += dphi_T1[k][1]*Un1[0][k]/C1[k]*vol_T1;
                dV1dx += dphi_T1[k][0]*Un1[1][k]/C1[k]*vol_T1; dV1dy += dphi_T1[k][1]*Un1[1][k]/C1[k]*vol_T1;
                dU2dx += dphi_T2[k][0]*Un2[0][k]/C2[k]*vol_T2; dU2dy += dphi_T2[k][1]*Un2[0][k]/C2[k]*vol_T2;
                dV2dx += dphi_T2[k][0]*Un2[1][k]/C2[k]*vol_T2; dV2dy += dphi_T2[k][1]*Un2[1][k]/C2[k]*vol_T2;
            }
            // Remember the jump in the edge structure
            //edges_vec[ii].jumpX = 0.5*(dU1dx+dU2dx)*bool_velocity_U + 0.5*(dPn1dx+dPn2dx)*bool_pressure + 0.5*(dV1dx+dV2dx)*bool_velocity_V;
            //edges_vec[ii].jumpY = 0.5*(dU1dy+dU2dy)*bool_velocity_U + 0.5*(dPn1dy+dPn2dy)*bool_pressure + 0.5*(dV1dy+dV2dy)*bool_velocity_V;
            
            // Return jump pressure gradient for each triangle
            gradient_x[el0*3+kk[el0]] = 0.5*(dU1dx+dU2dx)*bool_velocity_U + 0.5*(dPn1dx+dPn2dx)*bool_pressure + 0.5*(dV1dx+dV2dx)*bool_velocity_V;
            gradient_y[el0*3+kk[el0]] = 0.5*(dU1dy+dU2dy)*bool_velocity_U + 0.5*(dPn1dy+dPn2dy)*bool_pressure + 0.5*(dV1dy+dV2dy)*bool_velocity_V;
            
            gradient_x[el1*3+kk[el1]] = 0.5*(dU1dx+dU2dx)*bool_velocity_U + 0.5*(dPn1dx+dPn2dx)*bool_pressure + 0.5*(dV1dx+dV2dx)*bool_velocity_V;
            gradient_y[el1*3+kk[el1]] = 0.5*(dU1dy+dU2dy)*bool_velocity_U + 0.5*(dPn1dy+dPn2dy)*bool_pressure + 0.5*(dV1dy+dV2dy)*bool_velocity_V;
            kk[el0]++; kk[el1]++;
            
            
        }
    }
    
    free(kk);
}


/*
 The function fluid_problem_adapt_gen_mesh_in_advance aims to memorize the smallest element size at each iteration to produce in the end an "in adavance" adapted mesh. This function
 applies the gradient recovery estimator over a given file (porosity, velocity or pressure). It computes for each element the estimator and the new mesh size. Finally, it write the result
 in an attribute min_h_target which is a list of minimum h value of the structure Mesh. We do not need any file to transport form one iteration to the other.
 
 The file mesh.c contains the function that initializes the min_h_target list. Notice that this operation should be done at the beginning of the "in advance" strategy and each time the 
 mesh is physical modified by an lc.pos.
 
 @inputs: - problem  ==> a structure containing the solution and the mesh
          - lcmax    ==> a double being the maximum value that an element can have locally
          - lcmin    ==> a double being the minimum value that an element can have locally (very important for the scale separation hypothesis)
          - cmax     ==> a double ...
          - cmin     ==> a double ...
          - e_target ==> the target error over one element
          - U_0      ==> the characteristic velocity of the problem
          - P_0      ==> the characteristic pressure of the problem
 */
void fluid_problem_adapt_gen_mesh_in_advance(FluidProblem *problem, Mesh *mesh, double lcmax, double lcmin, double n_el, double cmax, double cmin, double e_target, double U_0, double P_0){
    // These definition should be more general
    double h_max = lcmax;
    double h_min = lcmin;
    double alpha = 2;
    double h_k = 0.0, h_target = 0.0;
    double tol = 1e-12;
    double gradC_x = 0.0, gradC_y = 0.0;
    double u_bar_x = 0.0, u_bar_y = 0.0;
    double v_bar_x = 0.0, v_bar_y = 0.0;
    double p_bar_x = 0.0, p_bar_y = 0.0;
    int bool1, bool2, bool3;
    
    // Definition of points and weights for the quadrature
    double PQ[4][2] ={{1./3, 1./3},{3./5, 1./5},{1./5, 3./5},{1./5, 1./5}};
    double w[4]={-0.281250000000000, +0.260416666666667, +0.260416666666667, +0.260416666666667};
    
    // Problem and solution definition
    double ngradM = 0.;
    //Mesh *mesh = problem->mesh;
    double *solution = problem->solution;
    
    double J, eta;
    
    // Let's call the function to handle the edges
    //int n_edges = (mesh->n_elements+1) + (mesh->n_nodes) - 2; // Euler's formula for planar graph
    //Edge *edges_vec = (Edge*) calloc(n_edges, sizeof(Edge));
    //initialized_edges_array(mesh, edges_vec, n_edges);
    // Horizontal velocity jump gradient
    double *grad_u_x = (double*) calloc(3*mesh->n_elements, sizeof(double));
    double *grad_u_y = (double*) calloc(3*mesh->n_elements, sizeof(double));
    gradient_recovery_estimator(problem, grad_u_x, grad_u_y, 0, 1, 0);
    
    // Vertical velocity jump gradient
    double *grad_v_x = (double*) calloc(3*mesh->n_elements, sizeof(double));
    double *grad_v_y = (double*) calloc(3*mesh->n_elements, sizeof(double));
    gradient_recovery_estimator(problem, grad_v_x, grad_v_y, 0, 0, 1);
    
    // Pressure jump gradient
    double *grad_p_x = (double*) calloc(3*mesh->n_elements, sizeof(double));
    double *grad_p_y = (double*) calloc(3*mesh->n_elements, sizeof(double));
    gradient_recovery_estimator(problem, grad_p_x, grad_p_y, 1, 0, 0);
    
    
    int n_fields = fluid_problem_n_fields(problem);
    double *porosity = problem->porosity;
    // Loop over each elements in the domain
    for(int iel = 0; iel<mesh->n_elements; iel++){
        const unsigned int *el = &mesh->elements[iel*N_N];
        double C[N_N],Pn[N_N],Un[D][N_N];
        // C -> porosity (VECTOR)
        // Pn -> nodal pressure (VECTOR)
        // Un -> nodal velocity in 2D or 3D (MATRIX)
        for (int i = 0; i < N_N; ++i){ // the for-loop is done over the N_N nodes of the elements iel
            C[i] = porosity[el[i]];
            Pn[i] = solution[el[i]*n_fields+D]-problem->g*problem->rho[0]*mesh->x[el[i]*3+1];
            //printf("Pn[%d]=%.6f\n", i, Pn[i]);
            //printf("p = %.6f\n", solution[el[i]*n_fields+D]);
            //printf("P0 = %.6f\n", P_0);
            for (int jj = 0; jj < D; ++jj) {
                Un[jj][i] = solution[el[i]*n_fields+jj];
            }
        }
        
        double dxdxi[D][D], dxidx[D][D];
        for(int i=0; i<D; i++){
            for(int j=0; j<D; j++){
                dxdxi[i][j] = mesh->x[el[j+1]*3+i]-mesh->x[el[0]*3+i];
            }
        }
#if D==2
        const double vol = detDD(dxdxi)/2;
        // We need to find the diameter of the circle circumscribing the triangle which id given by diam = 2R = abc/(2S)
        double X1 = mesh->x[el[0]*3+0], Y1 = mesh->x[el[0]*3+1];
        double X2 = mesh->x[el[1]*3+0], Y2 = mesh->x[el[1]*3+1];
        double X3 = mesh->x[el[2]*3+0], Y3 = mesh->x[el[2]*3+1];
        J = (X2-X1)*(Y3-Y1) - (X3-X1)*(Y2-Y1);
        double a = sqrt((X1-X2)*(X1-X2)+(Y1-Y2)*(Y1-Y2));
        double b = sqrt((X2-X3)*(X2-X3)+(Y2-Y3)*(Y2-Y3));
        double c = sqrt((X3-X1)*(X3-X1)+(Y3-Y1)*(Y3-Y1));
        h_k = 0.5*a*b*c/vol;
        
#else
        const double vol = detDD(dxdxi)/6;
        h_k = 0.0;
#endif
        // computation of pressure, porosity and velocity gradient
        invDD(dxdxi, dxidx);
        double dphi[N_SF][D];
        grad_shape_functions(dxidx, dphi);
        double dPndx = 0.0, dPndy = 0.0;
        double dCdx  = 0.0, dCdy  = 0.0;
        double dUdx  = 0.0, dUdy  = 0.0;
        double dVdx  = 0.0, dVdy  = 0.0;
        for (int k=0; k<N_SF; ++k){
            // Pressure gradient
            dPndx += dphi[k][0]*Pn[k]*vol;
            dPndy += dphi[k][1]*Pn[k]*vol;
            // Porosity gradient
            dCdx  += dphi[k][0]*C[k]*vol;
            dCdy  += dphi[k][1]*C[k]*vol;
            // Velocity gardient in u(x,y) and v(x,y) only for 2D problem
            dUdx += dphi[k][0]*Un[0][k]/C[k]*vol; dUdy += dphi[k][1]*Un[0][k]/C[k]*vol;
            dVdx += dphi[k][0]*Un[1][k]/C[k]*vol; dVdy += dphi[k][1]*Un[1][k]/C[k]*vol;
        }
        
        // Computation of the integrale over the reference triangle
        double I_h = 0.0; double eta_pressure=0.0;
        for (int i=0; i<4; i++) {
            // ---------------------------------------------------------------------------------------------
            // Linear interpolation of the middle jump value over the triangle, the shape function are,
            // f1 = 2x + 2y - 1, f2 = 1 - 2x and f3 = 1 - 2y
            // ---------------------------------------------------------------------------------------------
            
            // HORIZONTAL VELOCITY
            u_bar_x = (2*PQ[i][0]+2*PQ[i][1]-1)*grad_u_x[iel*3+0]+ (1-2*PQ[i][0])*grad_u_x[iel*3+1] + (1-2*PQ[i][1])*grad_u_x[iel*3+2];
            u_bar_y = (2*PQ[i][0]+2*PQ[i][1]-1)*grad_u_y[iel*3+0]+ (1-2*PQ[i][0])*grad_u_y[iel*3+1] + (1-2*PQ[i][1])*grad_u_y[iel*3+2];
            
            // VERTICAL VELOCITY
            v_bar_x = (2*PQ[i][0]+2*PQ[i][1]-1)*grad_v_x[iel*3+0]+ (1-2*PQ[i][0])*grad_v_x[iel*3+1] + (1-2*PQ[i][1])*grad_v_x[iel*3+2];
            v_bar_y = (2*PQ[i][0]+2*PQ[i][1]-1)*grad_v_y[iel*3+0]+ (1-2*PQ[i][0])*grad_v_y[iel*3+1] + (1-2*PQ[i][1])*grad_v_y[iel*3+2];
            
            // PRESSURE
            p_bar_x = (2*PQ[i][0]+2*PQ[i][1]-1)*grad_p_x[iel*3+0]+ (1-2*PQ[i][0])*grad_p_x[iel*3+1] + (1-2*PQ[i][1])*grad_p_x[iel*3+2];
            p_bar_y = (2*PQ[i][0]+2*PQ[i][1]-1)*grad_p_y[iel*3+0]+ (1-2*PQ[i][0])*grad_p_y[iel*3+1] + (1-2*PQ[i][1])*grad_p_y[iel*3+2];
            
            // INTEGRATION USING GAUSS QUADRATURE
            I_h += w[i]*((u_bar_x-dUdx)*(u_bar_x-dUdx)+(u_bar_y-dUdy)*(u_bar_y-dUdy))*J/(U_0*U_0);
            
            I_h += w[i]*((v_bar_x-dVdx)*(v_bar_x-dVdx)+(v_bar_y-dVdy)*(v_bar_y-dVdy))*J/(U_0*U_0);
            
            I_h += w[i]*((p_bar_x-dPndx)*(p_bar_x-dPndx)+(p_bar_y-dPndy)*(p_bar_y-dPndy))*J/(P_0*P_0);
        }
        eta = sqrt(I_h);
        h_target = h_k*sqrt(e_target/eta);
        bool1 = h_target<=h_k/alpha;  // then h_k = h_k/alpha;
        bool2 = h_target>=h_k*alpha;  // then h_k = h_k*alpha;
        bool3 = (h_target<h_k*alpha) && (h_target>h_k/alpha);
        h_target = bool1*fmax(h_k/alpha, h_min) + bool2*fmin(h_k*alpha, h_max) + bool3*h_k;
        mesh->min_h_target[iel] = fmin(h_target, 1.2*mesh->min_h_target[iel]);
    }
    //free(edges_vec);
    free(grad_u_x); free(grad_u_y);
    free(grad_v_x); free(grad_v_y);
    free(grad_p_x); free(grad_p_y);
    global_count++;
}

/*
 The function fluid_problem_adapt_gen_mesh aims to adapt the present mesh contained in the structure FluidProblem under a bunch of constraints
 and it write the adaptation inside the file filename under a adapt.pos file. 
 
 @inputs: - problem  ==> a structure containing the solution and the mesh
          - filename ==> a file to write the .pos in it, creating the metric map for GMSH
 
 @ouputs: return nothing but write the .pos file for GMSH purpose.
 */
void fluid_problem_adapt_gen_mesh(FluidProblem *problem, double lcmax, double lcmin, double n_el, double cmax, double cmin, double e_target, const char *filename, double U_0, double P_0)
{
    // Extracting mesh information
    Mesh *mesh = problem->mesh;
    double *new_mesh_size = (double*)malloc(sizeof(double)*mesh->n_nodes);
    double *num_lc = (double*)malloc(sizeof(double)*mesh->n_nodes);
    for (int i = 0; i < mesh->n_nodes; ++i){
        new_mesh_size[i] = 0.;
        num_lc[i] = 0.;
    }
    
    fluid_problem_adapt_gen_mesh_in_advance(problem, mesh, lcmax, lcmin, n_el, cmax, cmin, e_target, U_0, P_0);
    
    for (int iel = 0; iel<mesh->n_elements; iel++) {
        const unsigned int *el = &mesh->elements[iel*N_N];
        for (int j = 0; j < N_N; ++j){
            new_mesh_size[el[j]] += mesh->min_h_target[iel];
            num_lc[el[j]]+=1.;
        }
    }
    
    // Last loop to achieve the mean of h_target for the nodal position
    for(int i=0; i<mesh->n_nodes; i++){
        new_mesh_size[i] = new_mesh_size[i]/num_lc[i];
        //printf("%d \n", new_mesh_size[i]>lcmin);
    }
    
    // writing in the adapt.pos file
    FILE *f = fopen("adapt/lc.pos", "w");
    if(!f)
        printf("cannot open file adapt/lc.pos for writing\n");
    fprintf(f, "View \"lc\" {\n");
    for (int iel = 0; iel < mesh->n_elements; ++iel) {
        unsigned int *el = problem->mesh->elements+iel*N_N;
        fprintf(f, D == 2 ? "ST(" : "SS(");
        for (int i = 0; i < N_N; ++i)
            fprintf(f, "%.8g, %.8g, %.8g%s", mesh->x[el[i]*3+0], mesh->x[el[i]*3+1],mesh->x[el[i]*3+2], i != N_N-1 ? "," : "){");
        for (int i = 0; i < N_N; ++i)
            fprintf(f, "%.8g%s", new_mesh_size[el[i]],i != N_N-1?",":"};\n");
    }
    fprintf(f,"};\n");
    fclose(f);
    
    free(new_mesh_size);
    free(num_lc);
    
    system("gmsh -" xstr(D) " adapt/mesh.geo");
    
}

// allocate all mesh-dependant structures
static void fluid_problem_set_mesh(FluidProblem *problem, Mesh *mesh) {
    printf("Here in set_mesh\n");
  if(problem->mesh)
    mesh_free(problem->mesh);
  problem->mesh = mesh;
  if(problem->mesh_tree)
    mesh_tree_free(problem->mesh_tree);
    
  printf("CHECK1\n");
  problem->mesh_tree = mesh_tree_create(mesh);
  printf("CHECK2\n");
  free(problem->porosity);
  problem->porosity = (double*)malloc(mesh->n_nodes*sizeof(double));
  free(problem->oldporosity);
  problem->oldporosity = (double*)malloc(mesh->n_nodes*sizeof(double));
  for(int i = 0; i < mesh->n_nodes; ++i){
    problem->porosity[i] = 1.;
    problem->oldporosity[i] = 1.;
  }
  free(problem->solution);
  int n_fields = fluid_problem_n_fields(problem);
  problem->solution = (double*)malloc(mesh->n_nodes*n_fields*sizeof(double));
  if (problem->n_fluids == 2) {
    free(problem->concentration);
    problem->concentration = (double*)malloc(mesh->n_elements*N_N*sizeof(double));
  }
    printf("Transfert solution \n");
  problem->taup = malloc(sizeof(double)*mesh->n_elements);
  for (int i = 0; i < problem->mesh->n_elements; ++i) problem->taup[i] = 0;
  for (int i = 0; i < problem->mesh->n_nodes*n_fields; ++i){
    problem->solution[i] = 0.;
  }
  if (problem->linear_system)
    hxtLinearSystemDelete(&problem->linear_system);
#ifdef HXT_HAVE_PETSC
    hxtLinearSystemCreatePETSc(&problem->linear_system,mesh->n_elements,N_N,n_fields,mesh->elements,"fluid");
#else
    hxtLinearSystemCreateLU(&problem->linear_system,mesh->n_elements,N_N,n_fields,mesh->elements);
#endif
  fluid_problem_compute_node_volume(problem);
  if (problem->boundaries) mesh_boundaries_free(problem->boundaries,problem->mesh->n_boundaries);
  problem->boundaries = mesh_boundaries_create(problem->mesh);
}

void fluid_problem_set_mesh_and_transfer_solution(FluidProblem *problem, Mesh *new_mesh)
{
  int n_fields = fluid_problem_n_fields(problem);
  double *new_solution = (double*)malloc(sizeof(double)*new_mesh->n_nodes*n_fields);
  double *new_xi = (double*)malloc(sizeof(double)*new_mesh->n_nodes*D);
  int *new_eid = (int*)malloc(sizeof(int)*new_mesh->n_nodes);
  for (int i = 0; i < new_mesh->n_nodes; ++i)
    new_eid[i] = -1;
  double *newx = (double*)malloc(sizeof(double)*D*new_mesh->n_nodes);
  for (int i = 0;i < new_mesh->n_nodes;++i)
    for (int k = 0; k < D; ++k)
      newx[i*D+k] = new_mesh->x[i*3+k];
  mesh_tree_particle_to_mesh(problem->mesh_tree,new_mesh->n_nodes, newx, new_eid, new_xi);
  for (int i = 0; i < new_mesh->n_nodes; ++i) {
    unsigned int *el = problem->mesh->elements+new_eid[i]*N_N;
    if(new_eid[i] < 0) {
      printf("new mesh point not found in old mesh\n");
      exit(1);
    }
    double phi[N_SF];
    shape_functions(new_xi+i*D, phi);
    for (int j=0; j<n_fields; ++j) {
      new_solution[i*n_fields+j] = 0;
      for (int k = 0; k < N_SF; ++k)
        new_solution[i*n_fields+j] += problem->solution[el[k]*n_fields+j]*phi[k];
    }
  }
  free(new_eid);
  free(new_xi);
  free(newx);
  fluid_problem_set_mesh(problem,new_mesh);
  free(problem->solution);
  problem->solution = new_solution;
  for (int i = 0; i < problem->n_particles; ++i)
    problem->particle_element_id[i] = -1;
  mesh_tree_particle_to_mesh(problem->mesh_tree, problem->n_particles, problem->particle_position, problem->particle_element_id, problem->particle_uvw);
  fluid_problem_compute_porosity(problem);
}

void fluid_problem_adapt_load_mesh(FluidProblem *problem, const char *filename)
{
  Mesh *new_mesh = mesh_load_msh("adapt/mesh.msh");
  fluid_problem_set_mesh_and_transfer_solution(problem, new_mesh);
}



void fluid_problem_adapt_mesh(FluidProblem *problem, double e_target, double lcmax, double lcmin, int n_el, double cmax, double cmin, double U_0, double P_0)
{
    Mesh *mesh = problem->mesh;
    int bool_porosity = 0, bool_pressure = 0, bool_velocity = 1;
    
    fluid_problem_adapt_gen_mesh(problem, lcmax, lcmin, n_el, cmax, cmin, e_target, "adapt/lc.pos", U_0, P_0);
    //fluid_problem_adapt_load_mesh(problem,"adapt/mesh.msh");
}

/*
 This boolean is linked to the function mesh_initialized_min_h_target in the mesh.c file.
 if boolean_init == true,  then we initialize the min_h_target attribut of Mesh structure
 if boolean_init == false, then we do not
 */
void fluid_problem_estimator_evaluation(FluidProblem *problem, double e_target, double lcmax, double lcmin, double n_el, double cmax, double cmin, int boolean_init, double U_0, double P_0)
{
    Mesh *mesh = problem->mesh;
    int bool_porosity = 0, bool_pressure = 0, bool_velocity = 1;
    
    if (boolean_init==1) {
        printf("==> Initialize the min_h_target list \n");
        mesh_initialized_min_h_target(mesh);
        problem->mesh = mesh;
    }else{
        printf("==> Call to the fluid_problem_adapt_gen_mesh_in_advance function\n");
        fluid_problem_adapt_gen_mesh_in_advance(problem, mesh, lcmax, lcmin, n_el, cmax, cmin, e_target, U_0, P_0);
        problem->mesh = mesh;
    }
}


/*
 This function aims to create the new mesh according to the estimator. If bool_init_tree is set to 1, it means that it is the first time the mesh needs to be refined thus the 
 Tree_cell *Tree structure has to be initialized and to be maintained all along the simulation. 
 
 Otherwise, the Tree_cell *Tree structure already exists and a mapping between children tree and mesh elements was established. Using this Tree and the mapping, we are able to
 refine the current mesh.
 */

void fluid_problem_local_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin, int bool_init_tree){
    Mesh *mesh = problem->mesh;
    
    int n_elements_new = 0;
    int n_nodes_new = 0;
    int n_boundaries_new = 0, delta_boundaries_new = 0;
    
    //printing(mesh);
    //exit(1);
    
    if (bool_init_tree) {
        // INITIALISE THE TREE
        printf("==> INITIALIZED THE TREE STRUCTURE \n");
        problem->Tree = (Tree_cell*) malloc(mesh->n_elements*sizeof(Tree_cell));
        problem->size_initial_mesh = mesh->n_elements;
        problem->size_initial_node = mesh->n_nodes;
        deep_tree_initialized(mesh, problem->Tree, mesh->n_elements);
        
        // INITIALISE THE MAPPING
        problem->map = malloc(problem->size_initial_mesh*sizeof(Tree_cell*));
        deep_tree_get_mapping(problem->Tree, problem->map, problem->size_initial_mesh);
        
    }else{
        printf("==> ADAPT ON THE MESH STRUCTURE \n");
        int n_nodes_new = 0;
        int n_elements_new = 0;
        int n_boundaries_new = 0;
        int n_nodes_removed = 0;
        int n_elements_removed = 0;
        // Let's create the correspondence table for the removed node
        int len_node = mesh->n_nodes - problem->size_initial_node;
        int *tab_node = (int*) malloc(len_node*sizeof(int));
        for (int ii=0; ii<len_node; ii++) {
            tab_node[ii]=-1;
        }
        // Let's create the correpondence table for the removed element
        int len_elt = mesh->n_elements - problem->size_initial_mesh;
        int *tab_elt = (int*) malloc(len_elt*sizeof(int));
        for (int ii=0; ii<len_elt; ii++) {
            tab_elt[ii]=-1;
        }
        deep_tree_refinement(mesh, problem->Tree, problem->map, problem->size_initial_mesh, lcmin, &n_nodes_new, &n_elements_new, &n_boundaries_new, len_node, tab_node, problem->size_initial_node, len_elt, tab_elt, problem->size_initial_mesh, &n_nodes_removed, &n_elements_removed);
        
        printf("There are %d nodes, %d elements in the initial mesh\n", problem->size_initial_node, problem->size_initial_mesh);
        printf("There are %d nodes, %d elements in the old mesh\n", mesh->n_nodes, mesh->n_elements);
        printf("There are %d nodes, %d elements and %d boundaries in the mesh\n", n_nodes_new, n_elements_new, n_boundaries_new);
        int *elements_new = (int*) malloc(3*n_elements_new*sizeof(int));
        double *x_new = (double*) malloc(3*n_nodes_new*sizeof(double));
        int *boundaries_new = (int*) calloc(2*n_boundaries_new, sizeof(int));
        int *boundary_tags_new = (int*) malloc(n_boundaries_new*sizeof(int));

        Tree_cell **map_new = malloc(n_elements_new*sizeof(Tree_cell*));
        if (len_node==0 && len_elt==0) {
            for (int n=0; n<mesh->n_nodes; n++) {
                for (int ii=0; ii<3; ii++) {
                    x_new[n*3+ii] = mesh->x[n*3+ii];
                }
            }
        }else{
            for (int n=0; n<mesh->n_nodes; n++) {
                int m = n;
                if (n>=problem->size_initial_node && n<len_node+problem->size_initial_node) {
                    m = tab_node[n-problem->size_initial_node];
                }else if(n>=len_node+problem->size_initial_node){
                    m = n-n_nodes_removed;
                }
                if (m>-1) {
                    for (int ii=0; ii<3; ii++) {
                        x_new[m*3+ii] = mesh->x[n*3+ii];
                    }
                }
            }
        }
        
        deep_tree_travel(mesh, problem->Tree, problem->map, map_new, x_new, elements_new, boundaries_new, boundary_tags_new, n_nodes_new, n_elements_new, n_boundaries_new, problem->size_initial_mesh, len_node, tab_node, problem->size_initial_node, len_elt, tab_elt, problem->size_initial_mesh, n_nodes_removed, n_elements_removed);
        
        /*printf("Print in post_processing files \n");
        printing_post_processing_files(mesh, elements_new, x_new, boundaries_new, n_elements_new, n_boundaries_new);
        
        if (mesh->n_nodes>855) {
            int tag1 = 855;
            printf("(%f, %f)\n", x_new[tag1*3+0], x_new[tag1*3+1]);
        }
        
        if (mesh->n_nodes>916) {
            int tag1 = 916;
            printf("(%f, %f)\n", x_new[tag1*3+0], x_new[tag1*3+1]);
        }
        
        if (mesh->n_nodes>1158) {
            int tag1 = 1158;
            printf("(%f, %f)\n", x_new[tag1*3+0], x_new[tag1*3+1]);
        }
        
        if (mesh->n_nodes>1223) {
            int tag1 = 1223;
            printf("(%f, %f)\n", x_new[tag1*3+0], x_new[tag1*3+1]);
        }
        
        FILE *file2 = fopen("/Users/margauxboxho/migflow/testcases/my_local_drop/element_new.dat", "w");
        if(file2 == NULL){
            printf("Error at opening file \n");
            exit(1);
        }
        for (int iel=0; iel<n_elements_new; iel++) {
            fprintf(file2, "%4d, %4d, %4d, ", elements_new[iel*3+0], elements_new[iel*3+1], elements_new[iel*3+2]);
        }
        fclose(file2);*/

        
        free(problem->map);
        problem->map = NULL;
        problem->map = map_new;
        
        // CHECK ORIENTATION
        for (int i= 0; i < n_elements_new; ++i) {
            double *x[3];
            for (int j = 0; j < 3; ++j) {
                x[j] = &x_new[elements_new[i*3+j]*3];
            }
            double dx0[2] = {x[1][0]-x[0][0],x[1][1]-x[0][1]};
            double dx1[2] = {x[2][0]-x[0][0],x[2][1]-x[0][1]};
            double det = dx0[0]*dx1[1]-dx0[1]*dx1[0];
            if (det < 0) {
                int sw = elements_new[i*3];
                elements_new[i*3] = elements_new[i*3+1];
                elements_new[i*3+1] = sw;
            }
        }
        
        fluid_problem_set_elements(problem, n_nodes_new, x_new, n_elements_new, elements_new, n_boundaries_new, boundaries_new, boundary_tags_new, mesh->n_boundaries, mesh->boundary_names,1);
        
        free(elements_new);
        free(x_new);
        free(boundaries_new);
        free(boundary_tags_new);
        free(tab_node);
        free(tab_elt);
    }
}

void fluid_problem_load_msh(FluidProblem *problem, const char *mesh_file_name) {
  Mesh *mesh = mesh_load_msh(mesh_file_name);
  if (!mesh){
    printf("cannot import mesh file '%s'\n", mesh_file_name);
  }
  fluid_problem_set_mesh(problem,mesh);
}


void fluid_problem_after_import(FluidProblem *problem)
{
  if(problem->mesh_tree)
    mesh_tree_free(problem->mesh_tree);
  problem->mesh_tree = mesh_tree_create(problem->mesh);
  int n_fields = fluid_problem_n_fields(problem);
  for (int i = 0; i < problem->n_particles; ++i)
    problem->particle_element_id[i] = -1;
  mesh_tree_particle_to_mesh(problem->mesh_tree, problem->n_particles, problem->particle_position, problem->particle_element_id, problem->particle_uvw);
  fluid_problem_compute_node_volume(problem);
  fluid_problem_compute_porosity(problem);
  if(problem->linear_system)
    hxtLinearSystemDelete(&problem->linear_system);
  #ifdef HXT_HAVE_PETSC
    hxtLinearSystemCreatePETSc(&problem->linear_system,problem->mesh->n_elements,N_N,n_fields,problem->mesh->elements,"fluid");
  #else
    hxtLinearSystemCreateLU(&problem->linear_system,problem->mesh->n_elements,N_N,n_fields,problem->mesh->elements);
  #endif 
}


void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity, double *contact, long int *elid, int init, int reload)
{
  
  struct timespec tic, toc;
  int n_fields = fluid_problem_n_fields(problem);
  
  if(problem->n_particles != n) {
    problem->n_particles = n;
    free(problem->particle_uvw);
    free(problem->particle_element_id);
    free(problem->particle_position);
    free(problem->particle_mass);
    free(problem->particle_force);
    free(problem->particle_volume);
    free(problem->particle_velocity);
    free(problem->particle_contact);
    problem->particle_uvw = (double*)malloc(sizeof(double)*D*n);
    problem->particle_element_id = (int*)malloc(sizeof(int)*n);
    for (int i = 0; i < n; ++i) {
      problem->particle_element_id[i]=-1;
    }
    problem->particle_position = (double*)malloc(sizeof(double)*D*n);
    problem->particle_mass = (double*)malloc(sizeof(double)*n);
    problem->particle_force = (double*)malloc(sizeof(double)*n*D);
    problem->particle_volume = (double*)malloc(sizeof(double)*n);
    problem->particle_velocity = (double*)malloc(sizeof(double)*n*D);
    problem->particle_contact = (double*)malloc(sizeof(double)*n*D);
  }
  if (elid) {
    for (int i = 0; i < n; ++i) {
      problem->particle_element_id[i]=elid[i];
    }
  }
  mesh_tree_particle_to_mesh(problem->mesh_tree, n, position, problem->particle_element_id, problem->particle_uvw);
  for (int i = 0; i < n; ++i) {
    problem->particle_mass[i] = mass[i];
    problem->particle_volume[i] = volume[i];
    for (int k = 0; k < D; ++k) {
      problem->particle_position[i*D+k] = position[i*D+k];
      problem->particle_velocity[i*D+k] = velocity[i*D+k];
      problem->particle_contact[i*D+k] = contact[i*D+k];
    }
  }
  if (!reload){
    for (int i=0; i<problem->mesh->n_nodes; ++i){
      problem->oldporosity[i] = problem->porosity[i];
    }
  }
  fluid_problem_compute_porosity(problem);
  if (init){
    for (int i=0; i<problem->mesh->n_nodes; ++i){
      problem->oldporosity[i] = problem->porosity[i];
    }
  }
}

int *fluid_problem_particle_element_id(FluidProblem *problem)
{
  return problem->particle_element_id;
}

int fluid_problem_n_particles(FluidProblem *problem)
{
  return problem->n_particles;
}

double *fluid_problem_solution(const FluidProblem *p)
{
  return p->solution;
}

double *fluid_problem_coordinates(const FluidProblem *p)
{
  return p->mesh->x;
}

int fluid_problem_n_nodes(const FluidProblem *p)
{
  return p->mesh->n_nodes;
}

void fluid_problem_set_strong_boundary(FluidProblem *problem, const char *tag, int field, int row, BoundaryCallback *apply)
{
  problem->n_strong_boundaries++;
  problem->strong_boundaries = realloc(problem->strong_boundaries,sizeof(StrongBoundary)*problem->n_strong_boundaries);
  int i = problem->n_strong_boundaries-1;
  problem->strong_boundaries[i].tag = strdup(tag);
  problem->strong_boundaries[i].apply = apply;
  problem->strong_boundaries[i].field = field;
  problem->strong_boundaries[i].row = row;
}

void fluid_problem_set_weak_boundary(FluidProblem *p, const char *tag, BoundaryType type, BoundaryCallback *cb, int vid, int pid, int cid, int aid)
{
  for (int i = 0; i < p->n_weak_boundaries; ++i) {
    if (strcmp(p->weak_boundaries[i].tag, tag) == 0){
      printf("Weak boundary is already defined for tag \"%s\".\n", tag);
      exit(1);
    }
  }
  p->n_weak_boundaries++;
  p->weak_boundaries = realloc(p->weak_boundaries,sizeof(WeakBoundary)*p->n_weak_boundaries);
  WeakBoundary *bnd = &p->weak_boundaries[p->n_weak_boundaries-1];
  bnd->tag = strdup(tag);
  bnd->vid = vid;
  bnd->pid = pid;
  bnd->cid = cid;
  bnd->aid = aid;
  bnd->type = type;
  bnd->field_cb = cb;
}

int fluid_problem_n_elements(const FluidProblem *p)
{
  return p->mesh->n_elements;
}

int *fluid_problem_elements(const FluidProblem *p)
{
  return p->mesh->elements;
}

double *fluid_problem_porosity(const FluidProblem *p)
{
  return p->porosity;
}

double *fluid_problem_stab_param(FluidProblem *problem)
{
  return problem->taup;
}

double *fluid_problem_old_porosity(const FluidProblem *p)
{
  return p->oldporosity;
}

int fluid_problem_n_boundaries(const FluidProblem *p){
  return p->mesh->n_boundaries;
}

void fluid_problem_boundary_name(const FluidProblem *p, int i, char **phys_name)
{
  *phys_name = p->mesh->boundary_names[i];
}

double *fluid_problem_element_size(const FluidProblem *p)
{
  return p->element_size;
}

void fluid_problem_set_elements(FluidProblem *p, int n_nodes, double *x, int n_elements, int *elements, int n_boundaries, int *boundaries, int *boundary_tags,int n_physicals, char **physicals, int transfer_solution)
{
    printf("Set mesh and transfert solution %i elements\n", n_elements);
  Mesh *m = malloc(sizeof(Mesh));
  m->n_elements = n_elements;
  m->elements = malloc(n_elements*sizeof(int)*N_N);
  memcpy(m->elements,elements,N_N*sizeof(int)*n_elements);
  m->n_nodes = n_nodes;
  m->n_elements = n_elements;
  m->x = malloc(n_nodes*3*sizeof(double));
  memcpy(m->x,x,sizeof(double)*3*n_nodes);
  m->elements = malloc(sizeof(int)*(DIMENSION+1)*n_elements);
  memcpy(m->elements,elements,(DIMENSION+1)*sizeof(int)*n_elements);
  m->n_boundaries = n_physicals;
  m->boundary_names = malloc(sizeof(char*)*n_physicals);
  for (int i = 0; i < n_physicals; ++i) {
    m->boundary_names[i] = strdup(physicals[i]);
  }
  printf("n_boundaries : %i\n",n_boundaries);
  /*for (int i = 0; i < n_boundaries; ++i) {
    printf("%i %i with tag %i \n",boundaries[i*2+0],boundaries[i*2+1],boundary_tags[i]);
  }*/
  mesh_gen_edges(m,n_boundaries, boundaries, boundary_tags);
    
  
  FILE *ftmp = fopen("test_mesh.msh","w");
  mesh_write_msh(m,ftmp);
  fclose(ftmp);
    
  if(transfer_solution)
    fluid_problem_set_mesh_and_transfer_solution(p,m);
  else
    fluid_problem_set_mesh(p,m);
}


int fluid_problem_mesh_n_boundaries(const FluidProblem *p) {
  return p->mesh->n_boundaries;
}

void fluid_problem_mesh_boundary_info(const FluidProblem *p, int bid, char **bname, int *bsize)
{
  *bname = p->mesh->boundary_names[bid];
  *bsize = p->boundaries[bid].n_interfaces;
}

void fluid_problem_mesh_boundary_interface_nodes(const FluidProblem *p, int bid, int *binterfaces)
{
  const MeshBoundary *bnd = &p->boundaries[bid];
  for (int i = 0; i < bnd->n_interfaces; ++i) {
    int *inter = &p->mesh->interfaces[bnd->interfaces[i]*4];
    int *el = &p->mesh->elements[inter[0]*N_N];
    int *cl = elbnd[inter[1]];
    for (int j= 0; j< N_LSF; ++j) {
      binterfaces[i*N_LSF+j] = el[cl[j]];
    }
  }
}

int weak_boundary_n_values(const WeakBoundary *wbnd) {
  return (wbnd->vid>=0)*DIMENSION+(wbnd->pid>=0)+(wbnd->cid>=0)+(wbnd->aid>=0);
}

void weak_boundary_values(const Mesh *mesh, const MeshBoundary *bnd, const WeakBoundary *wbnd, double *data) {
  int n_value = weak_boundary_n_values(wbnd);
  if (n_value != 0) {
    double *x = malloc(bnd->n_interfaces*sizeof(double)*N_LQP*D);
    for (int i = 0; i < bnd->n_interfaces; ++i){
      int *interface = &mesh->interfaces[bnd->interfaces[i]*4];
      const int *cl = elbnd[interface[1]];
      const int *el = &mesh->elements[interface[0]*N_N];
      for (int iqp = 0; iqp < N_LQP; ++iqp) {
        double phil[N_LQP];
        l_shape_functions(LQP[iqp],phil);
        for (int id = 0; id < D; ++id){
          double xx = 0;
          for (int iphi = 0; iphi < D; ++iphi) {
            xx += mesh->x[el[cl[iphi]]*3+id]*phil[iphi];
          }
          x[(i*N_LQP+iqp)*D+id] = xx; 
        }
      }
    }
    wbnd->field_cb(N_LQP*bnd->n_interfaces, x, data);
    free(x);
  }
}