fluid_problem.c

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

#define n_fields (DIMENSION+1)

#define newton_max_it 20
#define newton_atol 1e-12
#define newton_rtol 1e-5

#if DIMENSION == 2

#define N_SF 3
#define N_N 3
#define N_QP 3
static double QP[][2] = {{1./6,1./6},{1./6,2./3},{2./3,1./6}};
static double QW[] = {1./6,1./6,1./6};
static void shape_functions(double *uvw, double *sf)
{
  sf[0] = 1-uvw[0]-uvw[1];
  sf[1] = uvw[0];
  sf[2] = uvw[1];
}
static void grad_shape_functions(const double dxidx[2][2], double gsf[3][2])
{
  for (int i = 0;  i < 2; ++i) {
    gsf[0][i] = -dxidx[0][i]-dxidx[1][i];
    gsf[1][i] = dxidx[0][i];
    gsf[2][i] = dxidx[1][i];
  }
}
static double detDD(const double m[2][2])
{
  return m[0][0]*m[1][1] - m[0][1]*m[1][0];
}
static double invDD(const double m[2][2], double inv[2][2])
{
  double d = detDD(m);
  inv[0][0] = m[1][1]/d;
  inv[0][1] = -m[0][1]/d;
  inv[1][0] = -m[1][0]/d;
  inv[1][1] = m[0][0]/d;
  return d;
}
#else

#define N_SF 4
#define N_N 4
#define N_QP 5
static double QP[][3] = {
  {0.25, 0.25, 0.25},
  {0.5, 1./6, 1./6},
  {1./6, 0.5, 1./6},
  {1./6, 1./6, 0.5},
  {1./6, 1./6, 1./6}
};
static double QW[] = {-16./120, 9./120, 9./120, 9./120, 9./120};
static void shape_functions(double *uvw, double *sf)
{
  sf[0] = 1-uvw[0]-uvw[1]-uvw[2];
  sf[1] = uvw[0];
  sf[2] = uvw[1];
  sf[3] = uvw[2];
}
static void grad_shape_functions(const double dxidx[3][3], double gsf[4][3])
{
  for (int i = 0; i < 3; ++i) {
    gsf[0][i] = -dxidx[0][i]-dxidx[1][i]-dxidx[2][i];
    gsf[1][i] = dxidx[0][i];
    gsf[2][i] = dxidx[1][i];
    gsf[3][i] = dxidx[2][i];
  }
};
static double detDD(const double m[3][3])
{
  return m[0][0]*(m[1][1]*m[2][2]-m[2][1]*m[1][2])
        -m[1][0]*(m[0][1]*m[2][2]-m[2][1]*m[0][2])
        +m[2][0]*(m[0][1]*m[1][2]-m[1][1]*m[0][2]);
}
static double invDD(const double m[3][3], double inv[3][3]) {
  double d = detDD(m);
  inv[0][0] = (m[1][1]*m[2][2]-m[2][1]*m[1][2])/d;
  inv[0][1] = -(m[0][1]*m[2][2]-m[2][1]*m[0][2])/d;
  inv[0][2] = (m[0][1]*m[1][2]-m[1][1]*m[0][2])/d;
  inv[1][0] = -(m[1][0]*m[2][2]-m[2][0]*m[1][2])/d;
  inv[1][1] = (m[0][0]*m[2][2]-m[2][0]*m[0][2])/d;
  inv[1][2] = -(m[0][0]*m[1][2]-m[1][0]*m[0][2])/d;
  inv[2][0] = (m[1][0]*m[2][1]-m[2][0]*m[1][1])/d;
  inv[2][1] = -(m[0][0]*m[2][1]-m[2][0]*m[0][1])/d;
  inv[2][2] = (m[0][0]*m[1][1]-m[1][0]*m[0][1])/d;
  return d;
};

#endif

static int check_word_in_file(FILE *f, const char *w) {
  char word[256];
  if(fscanf(f,"%255s", word) != 1){
    printf("Invalid mesh file (\"%s\" expected).\n", w);
    return -1;
  }
  if(strcmp(word, w) != 0){
    printf("Invalid mesh file (\"%s\" expected).\n", w);
    return -1;
  }
  return 0;
}

int fluid_problem_export(const FluidProblem *problem, const char *output_dir, double t, int iter)
{
  char file_name[1024];
  sprintf(file_name,"%s/fluid_%05i.msh",output_dir, iter);
  FILE *f = fopen(file_name, "w");
  if (!f){
    printf("Cannot open file \"%s\" for writing.\n", file_name);
    return -1;
  }
  if(mesh_write_msh(problem->mesh, f))
    return -1;
  if(mesh_write_msh_scalar(problem->mesh, f, "porosity", t, iter, problem->porosity, 1, 0))
    return -1;
  int comp[]={0,1,2};
  if(mesh_write_msh_vector(problem->mesh, f, "uv", t, iter, problem->solution, n_fields, comp))
    return -1;
   if(mesh_write_msh_vector(problem->mesh, f, "force", t, iter, problem->node_force, DIMENSION, comp))
    return -1;
  if(mesh_write_msh_scalar(problem->mesh, f, "p", t, iter, problem->solution, n_fields, DIMENSION))
    return -1;
  fclose(f);
  return 0;
}

int fluid_problem_export_vtk(const FluidProblem *problem, const char *output_dir, double t, int iter)
{
  Mesh *mesh = problem->mesh;
  char file_name[1024];
  sprintf(file_name,"%s/fluid_%05i.vtu",output_dir, iter);
  FILE *f = fopen(file_name, "w");
  if (!f){
    printf("Cannot open file \"%s\" for writing.\n", file_name);
    return -1;
  }
  fprintf(f,"<VTKFile type=\"UnstructuredGrid\" format=\"0.1\">\n");
  fprintf(f,"<UnstructuredGrid>\n");
  fprintf(f,"<Piece NumberOfPoints=\"%i\" NumberOfCells=\"%i\">\n", mesh->n_nodes, mesh->n_elements);
  fprintf(f,"<Points>\n");
  fprintf(f,"<DataArray NumberOfComponents=\"3\" type=\"Float64\" format=\"ascii\">\n");
  for (int i = 0; i < 3*mesh->n_nodes; ++i) {
    fprintf(f,"%g ", mesh->x[i]);
  }
  fprintf(f,"\n</DataArray>\n");
  fprintf(f,"</Points>\n");
  fprintf(f,"<Cells>\n");
  fprintf(f,"<DataArray  type=\"Int32\" Name=\"connectivity\" format=\"ascii\">\n");
  for (int i = 0; i < mesh->n_elements*(DIMENSION+1); ++i) {
    fprintf(f,"%i ", mesh->elements[i]);
  }
  fprintf(f,"\n</DataArray>\n");
  fprintf(f,"<DataArray  type=\"Int32\" Name=\"offsets\" format=\"ascii\">\n");
  for(int i = 0; i < mesh->n_elements;++i) {
    fprintf(f,"%i ", (i+1)*(DIMENSION+1));
  }
#if DIMENSION == 2
  fprintf(f,"\n</DataArray>\n");
  fprintf(f,"<DataArray  type=\"Int32\" Name=\"types\" format=\"ascii\">\n");
  for (int i = 0; i < mesh->n_elements; ++i) {
    fprintf(f,"%i ", 5); // 10 for 3d
  }
#else
  fprintf(f,"\n</DataArray>\n");
  fprintf(f,"<DataArray  type=\"Int32\" Name=\"types\" format=\"ascii\">\n");
  for (int i = 0; i < mesh->n_elements; ++i) {
    fprintf(f,"%i ", 10); // 10 for 3d
  }
#endif
  fprintf(f,"\n</DataArray>\n");
  fprintf(f,"</Cells>\n");
  //fprintf(f,"<PointData Scalars=\"Pressure Porosity\" Vectors=\"Velocity\">\n");
  fprintf(f,"<PointData Scalars=\"Porosity Pressure\" Vectors=\"Velocity\">\n");
  fprintf(f,"<DataArray Name=\"Porosity\" NumberOfComponents = \"1\" type=\"Float64\" format=\"ascii\">\n");
  for (int i = 0; i < mesh->n_nodes; ++i) {
    fprintf(f, "%g ", problem->porosity[i]);
  }
  fprintf(f,"\n</DataArray>\n");
  fprintf(f,"<DataArray Name=\"Pressure\" NumberOfComponents = \"1\" type=\"Float64\" format=\"ascii\">\n");
  for (int i = 0; i < mesh->n_nodes; ++i) {
    fprintf(f, "%g ", problem->solution[i*(DIMENSION+1)+DIMENSION]);
  }
  fprintf(f,"\n</DataArray>\n");
  fprintf(f,"<DataArray Name=\"Velocity\" NumberOfComponents = \"3\" type=\"Float64\" format=\"ascii\">\n");
  for (int i = 0; i < mesh->n_nodes; ++i) {
    for (int j = 0; j < DIMENSION; ++j)
      fprintf(f, "%g ", problem->solution[i*(DIMENSION+1)+j]);
    if(DIMENSION==2)
      fprintf(f, "0 ");
  }
  fprintf(f,"\n</DataArray>\n");
  fprintf(f,"</PointData>\n");
  fprintf(f,"</Piece>\n");
  fprintf(f,"</UnstructuredGrid>\n");
  fprintf(f,"</VTKFile>\n");

          //t[:] = 5 if dim == 2 else 10
/*
output = open("%s/f_%05i.vtu" %(odir, oiter), "w")
          vertices = np.array(self._mesh.vertices)[:,:3]
          elements = np.array([[v[3] for v in e.vertices] for e in self._dof.getElements(0)], np.int32) -1
            v = np.zeros((vertices.shape[0], 3))
            for d in range (dim) :
                v[elements, d] = self._jac.x[:,:,d]

            output.write("<Cells>\n")
            output.write("<DataArray  type=\"Int32\" Name=\"connectivity\" format=\"ascii\">\n")
            np.savetxt(output, elements, fmt="%i")
            output.write("</DataArray>\n")
            output.write("<DataArray  type=\"Int32\" Name=\"offsets\" format=\"ascii\">\n")
            np.savetxt(output, (np.array(range(elements.shape[0]), np.int32)+1)*(dim+1), fmt="%i", newline=" ")
            output.write("</DataArray>\n")
            output.write("<DataArray  type=\"Int32\" Name=\"types\" format=\"ascii\">\n")
            t = np.ndarray((elements.shape[0]), np.int32)
            np.savetxt(output, t, fmt="%i", newline=" ")
            output.write("</DataArray>\n")
            output.write("</Cells>\n")
            output.write("<PointData Scalars=\"Pressure Porosity\" Vectors=\"Velocity\">\n")
            output.write("<DataArray Name=\"Porosity\" NumberOfComponents = \"1\" type=\"Float64\" format=\"ascii\">")
            p = np.ndarray((vertices.shape[0]))
            p[np.transpose(elements)] = self._jac.dofManager.getField(0, self._pf.porosity)
            np.savetxt(output, p)
            output.write("</DataArray>\n")
            output.write("<DataArray Name=\"Pressure\" NumberOfComponents = \"1\" type=\"Float64\" format=\"ascii\">")
            p[np.transpose(elements)] = self._dof.getField(dim, self._sol)
            np.savetxt(output, p)
            output.write("</DataArray>\n")
            output.write("<DataArray Name=\"Velocity\" NumberOfComponents = \"3\" type=\"Float64\" format=\"ascii\">")
            v = np.zeros((vertices.shape[0], 3))
            for d in range (self._jac.dim) :
                v[np.transpose(elements), d] = self._dof.getField(d, self._sol)[:self._jac.dim+1,:]
            np.savetxt(output, v)
            output.write("</DataArray>\n")
            output.write("</PointData>\n")
            output.write("</Piece>\n")
            output.write("</UnstructuredGrid>\n")
            output.write("</VTKFile>\n")

    */
  fclose(f);
  return 0;
}

#if DIMENSION == 2
static void particle_drag(double u[2], double mu, double rho, double d, double c, double drag[2], double dt, double mass, double rhop, double gradp[2])
{
  double normu = hypot(u[0],u[1]);
  //Reynolds/|u_p-u_f|
  double Re_O_u = sqrt(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.;
  double GoU = f*Cd_u*rho/2*(2*r);

  double a = 0;
  double h = (1+(1-c)*rhop/(c*rho)+a)/(a+1);
  double F = (mass*GoU)/(mass+h*dt*GoU);

  drag[0] = (u[0]-gradp[0]*dt/rhop)*F;
  drag[1] = (u[1]-9.81*dt*(1-rho/rhop)-gradp[1]*dt/rhop)*F;
}
#else
static void particle_drag(double u[3], double mu, double rho, double d, double c, double drag[3], double dt, double mass, double rhop, double gradp[3])
{
  double normu = pow(u[0]*u[0]+u[1]*u[1]+u[2]*u[2],0.5);
  //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.;
  double GoU = f*Cd_u*rho/2*(M_PI*r*r);

  double a = 0;
  double h = (1+(1-c)*rhop/(c*rho)+a)/(a+1);
  double F = (mass*GoU)/(mass+h*dt*GoU);

  drag[0] = (u[0]-gradp[0]*dt/rhop)*F;
  drag[1] = (u[1]-9.81*dt*(1-rho/rhop)-gradp[1]*dt/rhop)*F;
  drag[2] = (u[2]-gradp[2]*dt/rhop)*F;

}
#endif

#if DIMENSION == 2
void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt, double *particle_force) {
  double *porosity = problem->porosity;
  Mesh *mesh = problem->mesh;
  const double *solution = problem->solution;
  for (int i = 0; i < mesh->n_nodes*DIMENSION; ++i) {
    problem->node_force[i] = 0.0;
  }
  for (int i = 0; i < problem->n_particles; ++i) {
    double *xi = problem->particle_uvw + 2*i;
    double phi[3] = {1-xi[0]-xi[1],xi[0],xi[1]};
    int iel = problem->particle_element_id[i];
    double vol = problem->particle_volume[i];
    double mass = problem->particle_mass[i];
    particle_force[i*2+0] = 0;
    particle_force[i*2+1] = problem->g*(mass-problem->rho*vol);
    if(iel < 0){
      continue;
    }

    int *tri = problem->mesh->elements+iel*3;
    double C[] = {porosity[tri[0]],porosity[tri[1]],porosity[tri[2]]};
    double U[] = {solution[tri[0]*3+0],solution[tri[1]*3+0],solution[tri[2]*3+0]};
    double V[] = {solution[tri[0]*3+1],solution[tri[1]*3+1],solution[tri[2]*3+1]};
    double P[] = {solution[tri[0]*3+2],solution[tri[1]*3+2],solution[tri[2]*3+2]};
    double vf[2] = {
      phi[0]*U[0]+phi[1]*U[1]+phi[2]*U[2],
      phi[0]*V[0]+phi[1]*V[1]+phi[2]*V[2]
    };
    double c = phi[0]*C[0]+phi[1]*C[1]+phi[2]*C[2];
    const double x[3][2] = {
      {mesh->x[tri[0]*3+0],mesh->x[tri[0]*3+1]},
      {mesh->x[tri[1]*3+0],mesh->x[tri[1]*3+1]},
      {mesh->x[tri[2]*3+0],mesh->x[tri[2]*3+1]}};
    const double dxdxi[2][2] = {
      {x[1][0]-x[0][0],x[2][0]-x[0][0]},
      {x[1][1]-x[0][1],x[2][1]-x[0][1]}};
    const double detj = dxdxi[0][0]*dxdxi[1][1]-dxdxi[0][1]*dxdxi[1][0];
    const double dxidx[2][2] = {
      {dxdxi[1][1]/detj,-dxdxi[0][1]/detj},
      {-dxdxi[1][0]/detj,dxdxi[0][0]/detj}
    };
    const double dphi[][2] = {
      {-dxidx[0][0]-dxidx[1][0],-dxidx[0][1]-dxidx[1][1]},
      {dxidx[0][0],dxidx[0][1]},
      {dxidx[1][0],dxidx[1][1]}
    };
    double gradp[2] = {
      dphi[0][0]*P[0]+dphi[1][0]*P[1]+dphi[2][0]*P[2],
      dphi[0][1]*P[0]+dphi[1][1]*P[1]+dphi[2][1]*P[2]
    };
    double du[2] = {
      problem->particle_velocity[i*2+0]-vf[0]/c,
      problem->particle_velocity[i*2+1]-vf[1]/c
    };
    double rhop = mass/vol;
    double d = 2*sqrt(vol/M_PI);
    double drag[2];
    particle_drag(du,problem->mu,problem->rho,d,c,drag, dt, mass, rhop, gradp);
    /*FILE *log = fopen("DragForce.txt", "at");
    if (!log) log = fopen("DragForce.txt", "wt");
    fprintf(log, "%d, %e, %e, %e, %e, %e, %e\n", i, du[0], du[1], drag[0], drag[1], mass, vol);
    fclose(log);*/

    particle_force[i*2+0] += -drag[0]-vol*gradp[0];
    particle_force[i*2+1] += -drag[1]-vol*gradp[1];
    //printf("%e,%e,%e,%e,%e,%e\n",particle_force[i*2+0],particle_force[i*2+1],-drag[0],-drag[1],0.,problem->g*(mass-problem->rho*vol));
    for (int iphi = 0; iphi < 3; ++iphi) {
      problem->node_force[tri[iphi]*2+0] += drag[0]*phi[iphi];
      problem->node_force[tri[iphi]*2+1] += drag[1]*phi[iphi];
    }
  }
  for (int i = 0; i < mesh->n_nodes; ++i) {
    problem->node_force[i*2+0] /= problem->node_volume[i];
    problem->node_force[i*2+1] /= problem->node_volume[i];
  }
}

#else
void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt, double *particle_force) {
  double *porosity = problem->porosity;
  Mesh *mesh = problem->mesh;
  const double *solution = problem->solution;
  for (int i = 0; i < mesh->n_nodes*3; ++i) {
    problem->node_force[i] = 0.0;
  }
  for (int i = 0; i < problem->n_particles; ++i) {
    double *xi = problem->particle_uvw + 3*i;
    double phi[4] = {1-xi[0]-xi[1]-xi[2],xi[0],xi[1],xi[2]};
    int iel = problem->particle_element_id[i];
    double vol = problem->particle_volume[i];
    double mass = problem->particle_mass[i];
    particle_force[i*3+0] = 0;
    particle_force[i*3+1] = problem->g*(mass-problem->rho*vol);
    particle_force[i*3+2] = 0;
    if(iel < 0){
      continue;
    }

    int *tet = problem->mesh->elements+iel*4;
    double C[] = {porosity[tet[0]],porosity[tet[1]],porosity[tet[2]],porosity[tet[3]]};
    double U[] = {solution[tet[0]*4+0],solution[tet[1]*4+0],solution[tet[2]*4+0],solution[tet[3]*4+0]};
    double V[] = {solution[tet[0]*4+1],solution[tet[1]*4+1],solution[tet[2]*4+1],solution[tet[3]*4+1]};
    double W[] = {solution[tet[0]*4+2],solution[tet[1]*4+2],solution[tet[2]*4+2],solution[tet[3]*4+2]};
    double P[] = {solution[tet[0]*4+3],solution[tet[1]*4+3],solution[tet[2]*4+3],solution[tet[3]*4+3]};
    double vf[3] = {
      phi[0]*U[0]+phi[1]*U[1]+phi[2]*U[2]+phi[3]*U[3],
      phi[0]*V[0]+phi[1]*V[1]+phi[2]*V[2]+phi[3]*V[3],
      phi[0]*W[0]+phi[1]*W[1]+phi[2]*W[2]+phi[3]*W[3]
    };
    double c = phi[0]*C[0]+phi[1]*C[1]+phi[2]*C[2]+phi[3]*C[3];
    /*printf("c=%e\n",c);
    printf("phi[0]=%e\n",phi[0]);
    printf("phi[1]=%e\n",phi[1]);
    printf("phi[2]=%e\n",phi[2]);
    printf("phi[3]=%e\n",phi[3]);*/
    const double x[4][3] = {
      {mesh->x[tet[0]*3+0],mesh->x[tet[0]*3+1],mesh->x[tet[0]*3+2]},
      {mesh->x[tet[1]*3+0],mesh->x[tet[1]*3+1],mesh->x[tet[1]*3+2]},
      {mesh->x[tet[2]*3+0],mesh->x[tet[2]*3+1],mesh->x[tet[2]*3+2]},
      {mesh->x[tet[3]*3+0],mesh->x[tet[3]*3+1],mesh->x[tet[3]*3+2]}};
    const double dxdxi[3][3] = {
      {x[1][0]-x[0][0],x[2][0]-x[0][0],x[3][0]-x[0][0]},
      {x[1][1]-x[0][1],x[2][1]-x[0][1],x[3][1]-x[0][1]},
      {x[1][2]-x[0][2],x[2][2]-x[0][2],x[3][2]-x[0][2]}};
    const double detj = dxdxi[0][0]*(dxdxi[1][1]*dxdxi[2][2]-dxdxi[2][1]*dxdxi[1][2])-dxdxi[1][0]*(dxdxi[0][1]*dxdxi[2][2]-dxdxi[2][1]*dxdxi[0][2])+dxdxi[2][0]*(dxdxi[0][1]*dxdxi[1][2]-dxdxi[1][1]*dxdxi[0][2]);
    const double dxidx[3][3] = {
      {(dxdxi[1][1]*dxdxi[2][2]-dxdxi[2][1]*dxdxi[1][2])/detj,-(dxdxi[0][1]*dxdxi[2][2]-dxdxi[2][1]*dxdxi[0][2])/detj, (dxdxi[0][1]*dxdxi[1][2]-dxdxi[1][1]*dxdxi[0][2])/detj},
      {-(dxdxi[1][0]*dxdxi[2][2]-dxdxi[2][0]*dxdxi[1][2])/detj,(dxdxi[0][0]*dxdxi[2][2]-dxdxi[2][0]*dxdxi[0][2])/detj,-(dxdxi[0][0]*dxdxi[1][2]-dxdxi[1][0]*dxdxi[0][2])/detj },
      {(dxdxi[1][0]*dxdxi[2][1]-dxdxi[2][0]*dxdxi[1][1])/detj,-(dxdxi[0][0]*dxdxi[2][1]-dxdxi[2][0]*dxdxi[0][1])/detj, (dxdxi[0][0]*dxdxi[1][1]-dxdxi[1][0]*dxdxi[0][1])/detj}
    };
    const double dphi[][3] = {
      {-dxidx[0][0]-dxidx[1][0]-dxidx[2][0],-dxidx[0][1]-dxidx[1][1]-dxidx[2][1],-dxidx[0][2]-dxidx[1][2]-dxidx[2][2]},
      {dxidx[0][0],dxidx[0][1],dxidx[0][2]},
      {dxidx[1][0],dxidx[1][1],dxidx[1][2]},
      {dxidx[2][0],dxidx[2][1],dxidx[2][2]}
    };
    double gradp[3] = {
      dphi[0][0]*P[0]+dphi[1][0]*P[1]+dphi[2][0]*P[2]+dphi[3][0]*P[3],
      dphi[0][1]*P[0]+dphi[1][1]*P[1]+dphi[2][1]*P[2]+dphi[3][1]*P[3],
      dphi[0][2]*P[0]+dphi[1][2]*P[1]+dphi[2][2]*P[2]+dphi[3][2]*P[3]
    };
    double du[3] = {
      problem->particle_velocity[i*3+0]-vf[0]/c,
      problem->particle_velocity[i*3+1]-vf[1]/c,
      problem->particle_velocity[i*3+2]-vf[2]/c
    };
    double rhop = mass/vol;
    double d = 2*pow(3./4.*vol/M_PI,1./3.);
    double drag[3];
    particle_drag(du,problem->mu,problem->rho,d,c,drag, dt, mass, rhop,gradp);
    particle_force[i*3+0] += -drag[0]-vol*gradp[0];
    particle_force[i*3+1] += -drag[1]-vol*gradp[1];
    particle_force[i*3+2] += -drag[2]-vol*gradp[2];
    //printf("drag[0]=%e\n",drag[0]);
    //printf("drag[1]=%e\n",drag[1]);
    //printf("drag[2]=%e\n",drag[2]);
    for (int iphi = 0; iphi < 4; ++iphi) {
      problem->node_force[tet[iphi]*3+0] += drag[0]*phi[iphi];
      problem->node_force[tet[iphi]*3+1] += drag[1]*phi[iphi];
      problem->node_force[tet[iphi]*3+2] += drag[2]*phi[iphi];
    }
  }
  for (int i = 0; i < mesh->n_nodes; ++i) {
    problem->node_force[i*3+0] /= problem->node_volume[i];
    problem->node_force[i*3+1] /= problem->node_volume[i];
    problem->node_force[i*3+2] /= problem->node_volume[i];
  }
}
#endif

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;
  const double *solution = problem->solution;
  const double *porosity = problem->porosity;
  const double *old_porosity = problem->old_porosity;
  const double *node_force = problem->node_force;
  const double mu = problem->mu;
  const double rho = problem->rho;
  const double epsilon = problem->epsilon;
  for (int iel=0; iel < mesh->n_elements; ++iel) {
    const unsigned int *el = &mesh->elements[iel*N_N];
    double dxdxi[DIMENSION][DIMENSION], dxidx[DIMENSION][DIMENSION], dphi[N_N][DIMENSION];
    for (int i = 0; i < DIMENSION; ++i)
      for (int j = 0; j < DIMENSION; ++j)
        dxdxi[i][j] = mesh->x[el[j+1]*3+i] - mesh->x[el[0]*3+i];
    const double detj = invDD(dxdxi, dxidx);
    grad_shape_functions(dxidx, dphi);
    double local_vector[n_fields*N_SF] = {0};
    double local_matrix[n_fields*n_fields*N_SF*N_SF] = {0};
    double dc[DIMENSION]={0}, du[DIMENSION][DIMENSION]={{0}}, dp[DIMENSION]={0};
    for (int i = 0; i< N_SF; ++i) {
      for (int j = 0; j < DIMENSION; ++j) {
        dc[j] += dphi[i][j]*porosity[el[i]];
        dp[j] += dphi[i][j]*solution[el[i]*n_fields+DIMENSION];
        for (int k = 0; k < DIMENSION; ++k)
          du[k][j] += dphi[i][j]*solution[el[i]*n_fields+k];
      }
    }
    for (int iqp = 0; iqp< N_QP; ++iqp) {
      double phi[N_SF];
      shape_functions(QP[iqp],phi);
      double u[DIMENSION]={0}, dudt[DIMENSION]={0}, fu[DIMENSION]={0}, p=0, c=0, dcdt=0;
      for (int i = 0; i < N_SF; ++i) {
        for (int j = 0; j < DIMENSION; ++j){
          u[j] += solution[el[i]*n_fields+j]*phi[i];
          dudt[j] += (solution[el[i]*n_fields+j]-solution_old[el[i]*n_fields+j])/dt*phi[i];
          fu[j] += node_force[el[i]*DIMENSION+j]*phi[i];
        }
        p += solution[el[i]*n_fields+DIMENSION]*phi[i];
        c += porosity[el[i]]*phi[i];
        dcdt += (porosity[el[i]]-old_porosity[el[i]])/dt*phi[i];
      }
      const double jw = QW[iqp]*detj;
      for (int iphi = 0; iphi < N_SF; ++iphi) {
        double phii = phi[iphi];
        const double *dphii = dphi[iphi];
        double tau[DIMENSION][DIMENSION];
        for (int i = 0; i < DIMENSION; ++i)
          for (int j = 0; j < DIMENSION; ++j)
            tau[i][j] = du[i][j]-u[i]*dc[j]/c;
        for (int j = 0; j < DIMENSION; ++j) {
          double utau = 0;
          double dphisigma=0;
          for (int k = 0; k < DIMENSION; ++k) {
            utau += u[k]*tau[j][k];
            dphisigma += 2*mu*dphii[k]*0.5*(tau[k][j]+tau[j][k]);
          }
          local_vector[iphi+N_SF*j] += jw*(
              dphisigma
              +rho*phii*(dudt[j]-u[j]/c*dcdt+utau/c)
              -p*dphii[j]-fu[j]*phii
              );
        }
        double dphidp = 0, divu = 0;
        for (int k = 0; k < DIMENSION; ++k){
          dphidp += dphii[k]*dp[k];
          divu += du[k][k];
        }
        local_vector[iphi+N_SF*DIMENSION] += jw*(epsilon*dphidp+phii*(divu+dcdt));
        for (int jphi = 0; jphi < N_SF; ++jphi) {
          double phij = phi[jphi];
          const double *dphij = dphi[jphi];
          int P = DIMENSION;
          //int N2 = N_SF*N_SF;
          //int IDX = iphi*N2+jphi;
          int N2 = n_fields*N_SF;
          int IDX = (iphi*N2+jphi)*n_fields;
#define LOCAL_MATRIX(a,b) local_matrix[IDX+N2*a+b]
          double dphiidphij = 0;
          double dtau[DIMENSION];
          double udtau = 0;
          double dphiidtau = 0;
          for (int j = 0; j < DIMENSION; ++j) {
            dtau[j] = dphij[j]-phij*dc[j]/c;
            dphiidphij += dphii[j]*dphij[j];
            udtau += u[j]*dtau[j];
            dphiidtau += dphii[j]*dtau[j];
          }
          for (int j = 0; j < DIMENSION; ++j) {
            int U = j;
            LOCAL_MATRIX(U,P) += -jw*dphii[j]*phij;
            LOCAL_MATRIX(P,U) +=  jw*phii*dphij[j];
            for (int k = 0; k < DIMENSION; ++k) {
              int V = k;
              LOCAL_MATRIX(U,V) += jw*rho*phii*phij*tau[j][k]/c+jw*2*mu*dphii[k]*dtau[j]*0.5;
            }
            LOCAL_MATRIX(U,U) += jw*mu*2*0.5*dphiidtau + jw*rho*phii*(phij/dt-phij/c*dcdt+udtau/c);
          }
          LOCAL_MATRIX(P,P) +=  jw*epsilon*dphiidphij;
        }
      }
    }
    hxtLinearSystemAddToMatrix(lsys,iel,iel,local_matrix);
    hxtLinearSystemAddToRhs(lsys,rhs,iel,local_vector);
  }

  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_phys; ++iphys) {
      if (strcmp(bnd->tag, mesh->phys_name[iphys]) == 0)
        break;
    }
    if (iphys == mesh->n_phys)
      printf("Boundary tag \"%s\" not found.", bnd->tag);
    for (int i = 0; i < mesh->phys_n_nodes[iphys]; ++i){
      hxtLinearSystemSetMatrixRowIdentity(lsys, mesh->phys_nodes[iphys][i], bnd->field);
      hxtLinearSystemSetRhsEntry(lsys, rhs,mesh->phys_nodes[iphys][i], bnd->field, 0.);
    }
  }
}



static void apply_boundary_conditions(const Mesh *mesh, double *solution, int nBnd, StrongBoundary *bnds)
{
  for (int ibnd = 0; ibnd < nBnd; ++ibnd) {
    StrongBoundary *bnd = bnds + ibnd;
    int iphys;
    for (iphys = 0; iphys < mesh->n_phys; ++iphys) {
      if (strcmp(bnd->tag, mesh->phys_name[iphys]) == 0)
        break;
    }
    if (iphys == mesh->n_phys){
      printf("Boundary tag \"%s\" not found.", bnd->tag);
    }
    double *x = malloc(mesh->phys_n_nodes[iphys]*sizeof(double)*DIMENSION);
    double *v = malloc(mesh->phys_n_nodes[iphys]*sizeof(double));
    for (int i = 0; i < mesh->phys_n_nodes[iphys]; ++i){
      int j = mesh->phys_nodes[iphys][i];
      for (int k = 0; k < DIMENSION; ++k)
        x[i*DIMENSION+k] = mesh->x[j*3+k];
    }
    bnd->apply(mesh->phys_n_nodes[iphys], x, v);
    for (int i = 0; i < mesh->phys_n_nodes[iphys]; ++i)
      solution[mesh->phys_nodes[iphys][i]*n_fields+bnd->field] = v[i];
    free(x);
    free(v);
  }
}

int fluid_problem_implicit_euler(FluidProblem *problem, double dt)
{
  static double totaltime = 0;
  struct timespec tic, toc;
  const Mesh *mesh = problem->mesh;
  apply_boundary_conditions(mesh, problem->solution, problem->n_strong_boundaries, problem->strong_boundaries);
  double *solution_old = malloc(mesh->n_nodes*n_fields*sizeof(double));
  double *rhs = 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];
  double *dx = malloc(sizeof(double)*mesh->n_nodes*n_fields);
  double norm0 = 0;
  int i;
  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 = %g\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];
    }
  }
  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 = 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 dxdxi[DIMENSION][DIMENSION];
    const int *el = &mesh->elements[iel*N_N];
    for (int i=0; i<DIMENSION; ++i)
      for (int j=0; j<DIMENSION; ++j)
        dxdxi[i][j] = mesh->x[el[j+1]*3+i]-mesh->x[el[0]*3+i];
#if DIMENSION == 2
    const double vol = detDD(dxdxi)/6;
#else
    const double vol = detDD(dxdxi)/24;
#endif
    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*DIMENSION],sf);
    const uint32_t *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){
    problem->porosity[i] = (1-problem->porosity[i]/problem->node_volume[i]);
  }
}

FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double mu, double rho, double epsilon, int n_strong_boundaries, StrongBoundary *strong_boundaries)
{
  static int initialized = 0;
  if (!initialized) {
    hxtInitializeLinearSystems(0, NULL);
    initialized = 1;
#ifdef HXT_HAVE_PETSC
    hxtPETScInsertOptions("-pc_type ilu", "fluid");
#endif
  }
  FluidProblem *problem = malloc(sizeof(FluidProblem));
  Mesh *mesh = mesh_load_msh(mesh_file_name);
  if (!mesh)
    return NULL;
  problem->rho = rho;
  problem->epsilon = epsilon;
  problem->mu = mu;
  problem->g = g;
  problem->mesh = mesh;
  problem->strong_boundaries = malloc(sizeof(StrongBoundary)*n_strong_boundaries);
  problem->n_strong_boundaries = n_strong_boundaries;
  for (int i = 0; i < n_strong_boundaries; ++i) {
    problem->strong_boundaries[i].tag = strdup(strong_boundaries[i].tag);
    problem->strong_boundaries[i].apply = strong_boundaries[i].apply;
    problem->strong_boundaries[i].field = strong_boundaries[i].field;
  }
  problem->porosity = malloc(mesh->n_nodes*sizeof(double));
  problem->node_force = malloc(mesh->n_nodes*DIMENSION*sizeof(double));
  for (int i = 0; i < mesh->n_nodes*DIMENSION; ++i) {
    problem->node_force[i] = 0.0;
  }
  problem->old_porosity = malloc(mesh->n_nodes*sizeof(double));
  problem->solution = malloc(mesh->n_nodes*n_fields*sizeof(double));
  // begin to remove
  for (int i = 0; i < problem->mesh->n_nodes*N_N; ++i){
    problem->solution[i] = 0.;
  }
  problem->node_volume = malloc(0);
  fluid_problem_compute_node_volume(problem);
  // end to remove
  problem->linear_system = NULL;
#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
  problem->n_particles = 0;
  problem->particle_uvw = malloc(0);
  problem->particle_element_id = malloc(0);
  problem->particle_position = malloc(0);
  problem->particle_mass = malloc(0);
  problem->particle_volume = malloc(0);
  problem->particle_velocity = malloc(0);
  return problem;
}

void fluid_problem_free(FluidProblem *problem)
{
  free(problem->solution);
  free(problem->node_force);
  free(problem->porosity);
  free(problem->old_porosity);
  hxtLinearSystemDelete(&problem->linear_system);
  free(problem->particle_uvw);
  free(problem->particle_element_id);
  free(problem->particle_position);
  free(problem->particle_mass);
  free(problem->particle_volume);
  free(problem->particle_velocity);
  for (int i = 0; i < problem->n_strong_boundaries; ++i)
    free(problem->strong_boundaries[i].tag);
  free(problem->strong_boundaries);
  mesh_free(problem->mesh);
}

void fluid_problem_adapt_gen_mesh(FluidProblem *problem, double gradmin, double gradmax, double gradPmin, double gradPmax, double lcmin, double n_el)
{
  double ngradM = 0.;
  Mesh *mesh = problem->mesh;
  double *solution = problem->solution;
  double *new_mesh_size = malloc(sizeof(double)*mesh->n_nodes);
  double *num_lc = malloc(sizeof(double)*mesh->n_nodes);
  for (int i = 0; i < mesh->n_nodes; ++i){
    new_mesh_size[i] = 0.;
    num_lc[i] = 0.;
  }
  gradmax = sqrt(gradmax);
  gradmin = sqrt(gradmin);
  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],U[DIMENSION][N_N], P[N_N];
    for (int i = 0; i < N_N; ++i){
      C[i] = porosity[el[i]];
      P[i] = solution[el[i]*n_fields+DIMENSION];
      for (int j=0; j<DIMENSION; ++j) {
        U[j][i] = solution[el[i]*n_fields+j];
      }
    }
    double dxdxi[DIMENSION][DIMENSION],dxidx[DIMENSION][DIMENSION];
    for (int i=0; i<DIMENSION; ++i)
      for (int j=0; j<DIMENSION; ++j)
        dxdxi[i][j] = mesh->x[el[j+1]*3+i]-mesh->x[el[0]*3+i];
    invDD(dxdxi, dxidx);
    double dphi[N_SF][DIMENSION];
    grad_shape_functions(dxidx, dphi);
    double ngrad2 = 0;
    for (int i=0; i<DIMENSION; ++i){
      for (int j=0; j<DIMENSION; ++j){
        double gradU = 0;
        for (int k=0; k<N_SF; ++k){
          gradU += dphi[k][j]*U[i][k]/C[k];
        }
        ngrad2 += gradU*gradU;
      }
    }
/*    double ngrad = pow(ngrad2, 0.25); */
/*    for (int j=0; j<DIMENSION; ++j){*/
/*      double gradP = 0;*/
/*      for (int k=0; k<N_SF; ++k){*/
/*        gradP += dphi[k][j]*P[k];*/
/*      }*/
/*      ngrad2 += gradP*gradP;*/
/*    }*/
    double ngrad = sqrt(ngrad2);
    double lc = lcmin /fmin(1,fmax(ngrad/(gradmax),gradmin/(gradmax)));
    
    ngrad = pow(ngrad2, 0.25); 
    for (int j=0; j<DIMENSION; ++j){
      double gradP = 0;
      for (int k=0; k<N_SF; ++k){
        gradP += dphi[k][j]*P[k];
      }
      ngrad2 += gradP*gradP;
    }
    ngrad = sqrt(ngrad2);
    lc = fmin(lcmin /fmin(1,fmax(ngrad/(gradPmax),gradPmin/(gradPmax))),lc);
    
    
    for (int j = 0; j < N_N; ++j){
      new_mesh_size[el[j]] += lc;
      num_lc[el[j]] += 1.;
    }
  }
  double nOfEl = 0.0;
#if DIMENSION==2
  for (int i = 0; i < mesh->n_nodes; ++i){
    new_mesh_size[i] /= num_lc[i];
    double nOfElByNode = 4*problem->node_volume[i]/(new_mesh_size[i]*new_mesh_size[i]*sqrt(3));
    nOfEl += nOfElByNode;
  }
  double ratioEl = n_el/nOfEl;
  for (int i = 0; i < mesh->n_nodes; ++i){
    new_mesh_size[i] = fmax(new_mesh_size[i]/sqrt(ratioEl),lcmin);
  }
#else
  for (int i = 0; i < mesh->n_nodes; ++i){
    new_mesh_size[i] /= num_lc[i];
    double nOfElByNode = 12*problem->node_volume[i]/(new_mesh_size[i]*new_mesh_size[i]*new_mesh_size[i]*sqrt(2));
    nOfEl += nOfElByNode;
  }
  double ratioEl = n_el/nOfEl;
  for (int i = 0; i < mesh->n_nodes; ++i){
    new_mesh_size[i] = fmax(new_mesh_size[i]/cbrt(ratioEl),lcmin);
  }
#endif
  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, DIMENSION == 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(DIMENSION)" adapt/mesh.geo");
}


void fluid_problem_adapt_mesh(FluidProblem *problem, double gradmin, double gradmax, double gradPmin, double gradPmax, double lcmin, double n_el)
{
  fluid_problem_adapt_gen_mesh(problem, gradmin, gradmax, gradPmin, gradPmax, lcmin, n_el);
  Mesh *new_mesh = mesh_load_msh("adapt/mesh.msh");
  double *new_solution = malloc(sizeof(double)*new_mesh->n_nodes*N_N);
  double *new_xi = malloc(sizeof(double)*new_mesh->n_nodes*DIMENSION);
  int *new_eid = malloc(sizeof(int)*new_mesh->n_nodes);
  for (int i = 0; i < new_mesh->n_nodes; ++i)
    new_eid[i] = -1;
  double *newx = malloc(sizeof(double)*DIMENSION*new_mesh->n_nodes);
  for (int i = 0;i < new_mesh->n_nodes;++i) 
    for (int k = 0; k < DIMENSION; ++k)
      newx[i*DIMENSION+k] = new_mesh->x[i*3+k];
  mesh_particles_to_mesh(problem->mesh, 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*DIMENSION, 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);
  free(problem->solution);
  problem->solution = new_solution;
  free(problem->porosity);
  problem->porosity = malloc(sizeof(double)*new_mesh->n_nodes);
  free(problem->old_porosity);
  problem->old_porosity = malloc(sizeof(double)*new_mesh->n_nodes);
  free(problem->node_force);
  problem->node_force = malloc(sizeof(double)*DIMENSION*new_mesh->n_nodes);
  mesh_free(problem->mesh);
  problem->mesh = new_mesh;
  for (int i = 0; i < problem->n_particles; ++i)
    problem->particle_element_id[i] = -1;
  mesh_particles_to_mesh(problem->mesh, problem->n_particles, problem->particle_position, problem->particle_element_id, problem->particle_uvw);
  fluid_problem_compute_node_volume(problem);
  fluid_problem_compute_porosity(problem);
  hxtLinearSystemDelete(&problem->linear_system);
  #ifdef HXT_HAVE_PETSC
    hxtLinearSystemCreatePETSc(&problem->linear_system,new_mesh->n_elements,N_N,n_fields,new_mesh->elements,"fluid");
  #else
    hxtLinearSystemCreateLU(&problem->linear_system,new_mesh->n_elements,N_N,n_fields,new_mesh->elements);
  #endif 
}


void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity, long int *elid)
{
  
  struct timespec tic, toc;
  
  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_volume);
    free(problem->particle_velocity);
    problem->particle_uvw = malloc(sizeof(double)*DIMENSION*n);
    problem->particle_element_id = malloc(sizeof(int)*n);
    for (int i = 0; i < n; ++i) {
      problem->particle_element_id[i]=-1;
    }
    problem->particle_position = malloc(sizeof(double)*DIMENSION*n);
    problem->particle_mass = malloc(sizeof(double)*n);
    problem->particle_volume = malloc(sizeof(double)*n);
    problem->particle_velocity = malloc(sizeof(double)*n*DIMENSION);
  }
  if (elid) {
    for (int i = 0; i < n; ++i) {
      problem->particle_element_id[i]=elid[i];
    }
  }
  clock_get(&tic);
  mesh_particles_to_mesh(problem->mesh, n, position, problem->particle_element_id, problem->particle_uvw);
  clock_get(&toc);
  printf("Time mesh_particles_to_mesh : %.6e \n",clock_diff(&tic,&toc));
  for (int i = 0; i < n; ++i) {
    problem->particle_mass[i] = mass[i];
    problem->particle_volume[i] = volume[i];
    for (int k = 0; k < DIMENSION; ++k) {
      problem->particle_position[i*DIMENSION+k] = position[i*DIMENSION+k];
      problem->particle_velocity[i*DIMENSION+k] = velocity[i*DIMENSION+k];
    }
  }
  for (int i = 0; i < problem->mesh->n_nodes; ++i){
    problem->old_porosity[i] = problem->porosity[i];
  }
  fluid_problem_compute_porosity(problem);
}

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

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