fluid_problem.c 31.1 KB
Newer Older
1
#include "tools.h"
2
3
#include <stdlib.h>
#include <stdio.h>
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
4
#include <float.h>
5
6
#include <math.h>
#include "fluid_problem.h"
7
#include "mesh_find.h"
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
8
#include "mbxml.h"
9

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
10
#define n_fields (DIMENSION+1)
11

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
12
#define newton_max_it 10
13
#define newton_atol 2.5e-6
14
#define newton_rtol 1e-5
15

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
16
#if DIMENSION==2
17
18
19

#define N_SF 3
#define N_N 3
20
21
22
#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};
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
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;
}
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
50
#else
51
52
53

#define N_SF 4
#define N_N 4
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
54
55
56
57
58
59
60
61
62
#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};
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
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;
};

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
99
100
#endif

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
101
102
#if DIMENSION==2
static double particle_drag_coeff(double u[2], double mu, double rho, double vol, double c)
103
{
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
104
  double d = 2*sqrt(vol/M_PI);
105
106
  double normu = hypot(u[0],u[1]);
  //Reynolds/|u_p-u_f|
Matthieu Constant's avatar
Matthieu Constant committed
107
  double Re_O_u = sqrt(d)*c/mu*rho;
108
109
110
  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);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
111
  return f*Cd_u*rho/2*d;
112
}
113
#else
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
114
static double particle_drag_coeff(double u[3], double mu, double rho, double vol, double c)
115
{
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
116
  double d = 2*cbrt(3./4.*vol/M_PI);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
117
  double normu = sqrt(u[0]*u[0]+u[1]*u[1]+u[2]*u[2]);
118
119
120
121
122
123
  //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.;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
124
  return f*Cd_u*rho/2*(M_PI*r*r);
125
126
}
#endif
127

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
128
void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt, double *particle_force) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
129
130
131
132
133
134
135
136
137
138
139
140
141
  for (int i = 0; i < problem->n_particles*DIMENSION; ++i) {
    particle_force[i] = problem->particle_force[i];
  }
}

static void compute_node_force_implicit(FluidProblem *problem, double dt, double *all_local_vector, double *all_local_matrix) {
  double *porosity = problem->porosity;
  Mesh *mesh = problem->mesh;
  for (int i = 0; i < problem->n_particles; ++i) {
    int iel = problem->particle_element_id[i];
    if(iel < 0){
      continue;
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
142
143
144
145
146
    double *xi = problem->particle_uvw + DIMENSION*i;
    double phi[N_SF];
    shape_functions(xi,phi);
    uint32_t *el=problem->mesh->elements+iel*N_SF;
    double *se = problem->solution_explicit;
147
    double c=0, vf[DIMENSION]={0}, vfe[DIMENSION]={0}, gradp[DIMENSION]={0};
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
148
149
150
151
152
153
154
    double dxdxi[DIMENSION][DIMENSION],dxidx[DIMENSION][DIMENSION];
    for (int k=0; k<DIMENSION; ++k)
      for (int j=0; j<DIMENSION; ++j)
        dxdxi[k][j] = mesh->x[el[j+1]*3+k]-mesh->x[el[0]*3+k];
    invDD(dxdxi,dxidx);
    double dphi[N_SF][DIMENSION];
    grad_shape_functions(dxidx, dphi);
155
    double *si = problem->solution;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
156
157
158
    for (int iphi=0; iphi < N_SF; ++iphi) {
      c += porosity[el[iphi]]*phi[iphi];
      for (int j=0; j < DIMENSION; ++j) {
159
        vf[j]  += si[el[iphi]*n_fields+j]*phi[iphi];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
160
        vfe[j] += se[el[iphi]*n_fields+j]*phi[iphi];
161
        gradp[j] += dphi[iphi][j]*si[el[iphi]*n_fields+DIMENSION];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
162
163
      }
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
164
    double due[DIMENSION];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
165
166
167
    for (int j = 0; j < DIMENSION; ++j) {
      due[j] = problem->particle_velocity[i*DIMENSION+j]-vfe[j]/c;
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
168
169
    double vol = problem->particle_volume[i];
    double mass = problem->particle_mass[i];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
170
171
    double rhop = mass/vol;
    double drho = rhop -problem->rho;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
172
173
174
    double gamma = particle_drag_coeff(due,problem->mu,problem->rho,vol,c);
    double g = problem->g*drho/rhop;
    double fcoeff = mass/(mass+dt*gamma);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
175
176
177
178
179
180
181
182
183
184
185

    double drag[DIMENSION];
    for (int j = 0; j < DIMENSION; ++j) {
      problem->particle_new_velocity[i*DIMENSION+j] =
        fcoeff * (
          problem->particle_velocity[i*DIMENSION+j]
          +gamma*vf[j]*dt/c/mass-vol*gradp[j]*dt/mass
          +(j == 1 ? g*dt : 0.)
          );
      drag[j] = gamma*(problem->particle_new_velocity[i*DIMENSION+j]-vf[j]/c);
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
186
187
    
    for (int d = 0; d < DIMENSION; ++d)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
188
189
      problem->particle_force[i*DIMENSION+d] = -drag[d] -vol*gradp[d];
    problem->particle_force[i*DIMENSION+1] += g*mass;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
190
    
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
191
192
    double *local_vector = all_local_vector+iel*n_fields*N_SF;
    double *local_matrix = all_local_matrix+iel*n_fields*n_fields*N_SF*N_SF;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
193
194
    for (int iphi = 0; iphi < N_SF; ++iphi) {
      for (int d = 0; d < DIMENSION; ++d)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
195
        local_vector[iphi+N_SF*d] -= drag[d]*phi[iphi];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
196
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
197
198
    for (int iphi = 0; iphi < N_SF; ++iphi) {
      for (int jphi = 0; jphi < N_SF; ++jphi) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
199
200
201
        int N2 = n_fields*N_SF;
        int IDX = (iphi*N2+jphi)*n_fields;
#define LOCAL_MATRIX(a,b) local_matrix[IDX+N2*a+b]
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
202
203
        double f = fcoeff*phi[iphi];
        for (int d = 0; d < DIMENSION; ++d){
204
205
          LOCAL_MATRIX(d,d) -= -f/c*phi[jphi]*gamma;
          LOCAL_MATRIX(d,DIMENSION) -= -f*gamma*dt/mass*vol*dphi[jphi][d];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
206
        }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
207
208
      }
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
209
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
  //compute new porosity
  double *new_position = malloc(sizeof(double)*problem->n_particles*DIMENSION);
  double *new_uvw = malloc(sizeof(double)*problem->n_particles*DIMENSION);
  int *new_eid = malloc(sizeof(int)*problem->n_particles);
  for (int i = 0; i < problem->n_particles; ++i) {
    new_eid[i] = problem->particle_element_id[i];
    for (int d = 0; d < DIMENSION; ++d) {
      new_position[i*DIMENSION+d] = problem->particle_position[i*DIMENSION+d]+dt*problem->particle_new_velocity[i*DIMENSION+d];
    }
  }
  mesh_tree_particle_to_mesh(problem->mesh_tree,problem->n_particles, new_position, new_eid, new_uvw);
  {
    double *volume = problem->particle_volume;
    for (int i = 0; i < mesh->n_nodes; ++i){
      problem->new_porosity[i] = 0;
    }
    for (int i = 0; i < problem->n_particles; ++i) {
      if(new_eid[i] == -1) continue;
      double sf[N_SF];
      shape_functions(&new_uvw[i*DIMENSION],sf);
      const uint32_t *el = &mesh->elements[new_eid[i]*N_N];
      for (int j=0; j<N_N;++j)
        problem->new_porosity[el[j]] += sf[j]*volume[i];
    }
    for (int i = 0; i < mesh->n_nodes; ++i){
      problem->new_porosity[i] = (1-problem->new_porosity[i]/problem->node_volume[i]);
    }
  }
  free(new_eid);
  free(new_uvw);
  free(new_position);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
241
242
243
244
245
246
247
248
249
}

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 mu = problem->mu;
  const double rho = problem->rho;
Matthieu Constant's avatar
Matthieu Constant committed
250
  const double epsilon = problem->epsilon;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
251
252
253
254
255
256
257
258
259
  size_t local_size = N_SF*n_fields;
  size_t N_E = mesh->n_elements;
  double *all_local_vector = malloc(N_E*local_size*sizeof(double));
  double *all_local_matrix = 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;
  compute_node_force_implicit(problem, dt, all_local_vector, all_local_matrix);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
260
  for (int iel=0; iel < mesh->n_elements; ++iel) {
261
262
263
264
265
266
267
    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);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
268
269
    double *local_vector = &all_local_vector[N_SF*n_fields*iel];
    double *local_matrix = &all_local_matrix[N_SF*N_SF*n_fields*n_fields*iel];
270
271
272
    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) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
273
        dc[j] += dphi[i][j]*problem->porosity[el[i]];
274
275
276
277
278
        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];
      }
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
279
    for (int iqp = 0; iqp< N_QP; ++iqp) {
280
281
      double phi[N_SF];
      shape_functions(QP[iqp],phi);
282
      double u[DIMENSION]={0}, dudt[DIMENSION]={0}, p=0, c=0, dcdt=0;
283
284
285
286
287
288
      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];
        }
        p += solution[el[i]*n_fields+DIMENSION]*phi[i];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
289
290
        c += problem->porosity[el[i]]*phi[i];
        dcdt += (problem->new_porosity[el[i]]-problem->porosity[el[i]])/dt*phi[i];
291
      }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
292
      const double jw = QW[iqp]*detj;
293
      for (int iphi = 0; iphi < N_SF; ++iphi) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
294
        double phii = phi[iphi];
295
296
297
298
299
300
301
302
303
        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) {
304
305
            //utau += u[k]*dphii[k]*u[j]/c;
            utau += u[k]*tau[j][k];
306
307
308
309
            dphisigma += 2*mu*dphii[k]*0.5*(tau[k][j]+tau[j][k]);
          }
          local_vector[iphi+N_SF*j] += jw*(
              dphisigma
310
311
312
              //+rho*phii*dudt[j]-rho*utau
              +rho*phii*(dudt[j]-u[j]/c*dcdt+utau/c)
              -p*dphii[j]
313
314
315
316
317
318
319
              );
        }
        double dphidp = 0, divu = 0;
        for (int k = 0; k < DIMENSION; ++k){
          dphidp += dphii[k]*dp[k];
          divu += du[k][k];
        }
Matthieu Constant's avatar
Matthieu Constant committed
320
        local_vector[iphi+N_SF*DIMENSION] += jw*(epsilon*dphidp+phii*(divu+dcdt));
321
        for (int jphi = 0; jphi < N_SF; ++jphi) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
322
          double phij = phi[jphi];
323
          const double *dphij = dphi[jphi];
324
325
326
327
328
          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;
329
330
331
332
333
334
335
336
#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];
337
            udtau += u[j]*dtau[j];
338
339
340
            dphiidtau += dphii[j]*dtau[j];
          }
          for (int j = 0; j < DIMENSION; ++j) {
341
            int U = j;
342
343
344
            LOCAL_MATRIX(U,P) += -jw*dphii[j]*phij;
            LOCAL_MATRIX(P,U) +=  jw*phii*dphij[j];
            for (int k = 0; k < DIMENSION; ++k) {
345
              int V = k;
346
347
              //LOCAL_MATRIX(U,V) += -jw*rho*u[j]*dphii[k]*phij/c+jw*2*mu*dphii[k]*dtau[j]*0.5;
              LOCAL_MATRIX(U,V) += jw*rho*phii*phij*tau[j][k]/c+jw*2*mu*dphii[k]*dtau[j]*0.5;
348
            }
349
350
            //LOCAL_MATRIX(U,U) += jw*mu*2*0.5*dphiidtau + jw*rho*(phii*phij/dt-phij*dphii[j]*u[j]/c);
            LOCAL_MATRIX(U,U) += jw*mu*2*0.5*dphiidtau + jw*rho*phii*(phij/dt-phij/c*dcdt+udtau/c);
351
          }
Matthieu Constant's avatar
Matthieu Constant committed
352
          LOCAL_MATRIX(P,P) +=  jw*epsilon*dphiidphij;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
353
354
355
356
357
        }
      }
    }
    hxtLinearSystemAddToMatrix(lsys,iel,iel,local_matrix);
    hxtLinearSystemAddToRhs(lsys,rhs,iel,local_vector);
358
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
359
360
  free(all_local_matrix);
  free(all_local_vector);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
361

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
362
363
364
365
366
367
368
369
370
371
  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){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
372
373
      hxtLinearSystemSetMatrixRowIdentity(lsys, mesh->phys_nodes[iphys][i], bnd->field);
      hxtLinearSystemSetRhsEntry(lsys, rhs,mesh->phys_nodes[iphys][i], bnd->field, 0.);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
374
375
    }
  }
376
}
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
377
378


379

380
static void apply_boundary_conditions(const Mesh *mesh, double *solution, int nBnd, StrongBoundary *bnds)
381
{
382
383
  for (int ibnd = 0; ibnd < nBnd; ++ibnd) {
    StrongBoundary *bnd = bnds + ibnd;
384
385
386
387
388
389
390
391
    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);
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
392
    double *x = malloc(mesh->phys_n_nodes[iphys]*sizeof(double)*DIMENSION);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
393
394
395
    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];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
396
397
      for (int k = 0; k < DIMENSION; ++k)
        x[i*DIMENSION+k] = mesh->x[j*3+k];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
398
399
    }
    bnd->apply(mesh->phys_n_nodes[iphys], x, v);
400
    for (int i = 0; i < mesh->phys_n_nodes[iphys]; ++i)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
401
402
403
      solution[mesh->phys_nodes[iphys][i]*n_fields+bnd->field] = v[i];
    free(x);
    free(v);
404
405
406
  }
}

407
int fluid_problem_implicit_euler(FluidProblem *problem, double dt)
408
{
409
410
  static double totaltime = 0;
  struct timespec tic, toc;
411
  const Mesh *mesh = problem->mesh;
412
  apply_boundary_conditions(mesh, problem->solution, problem->n_strong_boundaries, problem->strong_boundaries);
413
  double *solution_old = malloc(mesh->n_nodes*n_fields*sizeof(double));
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
414
  double *rhs = malloc(mesh->n_nodes*n_fields*sizeof(double));
415
416
  for (int i=0; i<mesh->n_nodes*n_fields; ++i)
    solution_old[i] = problem->solution[i];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
417
  problem->solution_explicit = solution_old;
418
419
  double *dx = malloc(sizeof(double)*mesh->n_nodes*n_fields);
  double norm0 = 0;
420
421
  int i;
  for (i=0; i<newton_max_it; ++i) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
422
423
424
425
426
427
    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);
428
429
430
431
432
433
434
    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;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
435
    clock_get(&tic);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
436
    hxtLinearSystemSolve(problem->linear_system, rhs, dx);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
437
438
    clock_get(&toc);
    totaltime += clock_diff(&tic, &toc);
439
440
441
442
    for (int j=0; j<mesh->n_nodes*n_fields; ++j) {
      problem->solution[j] -= dx[j];
    }
  }
443
  printf("total solve time : %g\n", totaltime);
444
  free(dx);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
445
  free(rhs);
446
  free(solution_old);
447
  return i;
448
449
}

Matthieu Constant's avatar
Matthieu Constant committed
450
451
452
453
454
455
456
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;
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
457
  for (int iel = 0; iel < mesh->n_elements; ++iel) {
458
459
460
461
462
    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];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
463
#if DIMENSION == 2
464
    const double vol = detDD(dxdxi)/6;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
465
#else
466
    const double vol = detDD(dxdxi)/24;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
467
#endif
468
469
    for (int i = 0; i < N_N; ++i)
      problem->node_volume[el[i]] += vol;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
470
471
472
473
474
475
476
477
478
479
480
481
  }
}

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;
482
483
    double sf[N_SF];
    shape_functions(&problem->particle_uvw[i*DIMENSION],sf);
484
    const uint32_t *el = &mesh->elements[problem->particle_element_id[i]*N_N];
485
486
    for (int j=0; j<N_N;++j)
      problem->porosity[el[j]] += sf[j]*volume[i];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
487
488
489
  }
  for (int i = 0; i < mesh->n_nodes; ++i){
    problem->porosity[i] = (1-problem->porosity[i]/problem->node_volume[i]);
Matthieu Constant's avatar
Matthieu Constant committed
490
491
492
  }
}

Matthieu Constant's avatar
Matthieu Constant committed
493
FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double mu, double rho, double epsilon, int n_strong_boundaries, StrongBoundary *strong_boundaries)
494
{
495
496
497
498
499
  static int initialized = 0;
  if (!initialized) {
    hxtInitializeLinearSystems(0, NULL);
    initialized = 1;
#ifdef HXT_HAVE_PETSC
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
500
    hxtPETScInsertOptions("-pc_type ilu -ksp_max_it 30", "fluid");
501
502
#endif
  }
503
  FluidProblem *problem = malloc(sizeof(FluidProblem));
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
504
505
506
  Mesh *mesh = mesh_load_msh(mesh_file_name);
  if (!mesh)
    return NULL;
507
  problem->mesh_tree = mesh_tree_create(mesh);
508
  problem->rho = rho;
Matthieu Constant's avatar
Matthieu Constant committed
509

510
  problem->mu = mu;
511
  problem->g = g;
512
  problem->mesh = mesh;
Matthieu Constant's avatar
Matthieu Constant committed
513
  problem->epsilon = epsilon;
514
515
516
517
  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);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
518
    problem->strong_boundaries[i].apply = strong_boundaries[i].apply;
519
520
    problem->strong_boundaries[i].field = strong_boundaries[i].field;
  }
521
  problem->porosity = malloc(mesh->n_nodes*sizeof(double));
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
522
  problem->new_porosity = malloc(mesh->n_nodes*sizeof(double));
523
  problem->solution = malloc(mesh->n_nodes*n_fields*sizeof(double));
Matthieu Constant's avatar
Matthieu Constant committed
524
525
526
527
  for (int i = 0; i < problem->mesh->n_nodes; ++i)
  {
    problem->porosity[i] = 1.;
  }
528
  // begin to remove
529
  for (int i = 0; i < problem->mesh->n_nodes*N_N; ++i){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
530
    problem->solution[i] = 0.;
531
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
532
  problem->node_volume = NULL;
Matthieu Constant's avatar
Matthieu Constant committed
533
  fluid_problem_compute_node_volume(problem);
534
  // end to remove
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
535
  problem->linear_system = NULL;
536
#ifdef HXT_HAVE_PETSC
Matthieu Constant's avatar
Matthieu Constant committed
537
    hxtLinearSystemCreatePETSc(&problem->linear_system,mesh->n_elements,N_N,n_fields,mesh->elements,"fluid");
538
#else
539
    hxtLinearSystemCreateLU(&problem->linear_system,mesh->n_elements,N_N,n_fields,mesh->elements);
540
#endif
541
  problem->n_particles = 0;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
542
543
544
545
546
547
548
549
  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_new_velocity = NULL;
550
551
552
553
554
555
556
  return problem;
}

void fluid_problem_free(FluidProblem *problem)
{
  free(problem->solution);
  free(problem->porosity);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
557
  free(problem->new_porosity);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
558
  hxtLinearSystemDelete(&problem->linear_system);
559
560
561
562
  free(problem->particle_uvw);
  free(problem->particle_element_id);
  free(problem->particle_position);
  free(problem->particle_mass);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
563
  free(problem->particle_force);
564
565
566
567
568
  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);
569
  mesh_tree_free(problem->mesh_tree);
570
571
572
  mesh_free(problem->mesh);
}

Matthieu Constant's avatar
Matthieu Constant committed
573
void fluid_problem_adapt_gen_mesh(FluidProblem *problem, double lcmax, double lcmin, double n_el)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
574
{
Matthieu Constant's avatar
Matthieu Constant committed
575
  double ngradM = 0.;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
576
577
578
  Mesh *mesh = problem->mesh;
  double *solution = problem->solution;
  double *new_mesh_size = malloc(sizeof(double)*mesh->n_nodes);
Matthieu Constant's avatar
Matthieu Constant committed
579
580
581
582
583
  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.;
  }
Matthieu Constant's avatar
Matthieu Constant committed
584
585
586
587
588
589
590
  double gradmax = 0.;
  double gradPmax = 0.;
  double gradmin = 1e8;
  double gradPmin = 1e8;
  double gradM = 0.;
  double gradPM = 0.;
  double volT = 0.;
Matthieu Constant's avatar
Matthieu Constant committed
591
  double *porosity = problem->porosity;
Matthieu Constant's avatar
Matthieu Constant committed
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
  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];
#if DIMENSION == 2
    const double vol = detDD(dxdxi)/2;
#else        
    const double vol = detDD(dxdxi)/6;
#endif    
    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;
      }
    }
    gradM += sqrt(ngrad2)*vol;
    volT += vol;
    gradmax = fmax(gradmax,ngrad2);
    gradmin = fmin(gradmin,ngrad2);
    ngrad2 = 0;
    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;
    }
    gradPM += sqrt(ngrad2)*vol;
    gradPmax = fmax(gradPmax,ngrad2);
    gradPmin = fmin(gradPmin,ngrad2);
  }
  gradPM /= volT;
  gradM /= volT;
  gradPmin = 20*sqrt(gradPmin)/gradPM;
  gradPmax = .1*sqrt(gradPmax)/gradPM;
  gradmin = 20*sqrt(gradmin)/gradM;
  gradmax = .1*sqrt(gradmax)/gradM;
  printf("gradPmax=%e\n",gradPmax);
  printf("gradPmin=%e\n",gradPmin);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
648
  for (int iel = 0; iel < mesh->n_elements; ++iel) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
649
    const unsigned int *el = &mesh->elements[iel*N_N];
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
    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;
Matthieu Constant's avatar
Matthieu Constant committed
666
    for (int i=0; i<DIMENSION; ++i){
667
668
669
670
671
672
673
674
      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;
      }
    }
Matthieu Constant's avatar
Matthieu Constant committed
675
676
677
678
679
680
681
682
/*    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;*/
/*    }*/
Matthieu Constant's avatar
Matthieu Constant committed
683
    double ngrad = sqrt(ngrad2)/gradM;
Matthieu Constant's avatar
Matthieu Constant committed
684
    double lc = lcmin /fmin(1,fmax(ngrad/(gradmax),gradmin/(gradmax)));
Matthieu Constant's avatar
Matthieu Constant committed
685
    
Matthieu Constant's avatar
Matthieu Constant committed
686
    ngrad2 = 0.; 
Matthieu Constant's avatar
Matthieu Constant committed
687
688
689
690
691
692
693
    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;
    }
Matthieu Constant's avatar
Matthieu Constant committed
694
695
    ngrad = sqrt(ngrad2)/gradPM;
    lc = fmin(fmin(lcmin /fmin(1,fmax(ngrad/(gradPmax),gradPmin/(gradPmax))),lc),lcmax);
Matthieu Constant's avatar
Matthieu Constant committed
696
697
    
    
Matthieu Constant's avatar
Matthieu Constant committed
698
699
700
701
702
    for (int j = 0; j < N_N; ++j){
      new_mesh_size[el[j]] += lc;
      num_lc[el[j]] += 1.;
    }
  }
Matthieu Constant's avatar
Matthieu Constant committed
703
  double nOfEl = 0.0;
Matthieu Constant's avatar
Matthieu Constant committed
704
#if DIMENSION==2
Matthieu Constant's avatar
Matthieu Constant committed
705
706
  for (int i = 0; i < mesh->n_nodes; ++i){
    new_mesh_size[i] /= num_lc[i];
Matthieu Constant's avatar
Matthieu Constant committed
707
708
709
710
711
712
713
    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);
  }
Matthieu Constant's avatar
Matthieu Constant committed
714
#else
Matthieu Constant's avatar
Matthieu Constant committed
715
  for (int i = 0; i < mesh->n_nodes; ++i){
Matthieu Constant's avatar
Matthieu Constant committed
716
717
    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));
Matthieu Constant's avatar
Matthieu Constant committed
718
    nOfEl += nOfElByNode;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
719
  }
Matthieu Constant's avatar
Matthieu Constant committed
720
  double ratioEl = n_el/nOfEl;
Matthieu Constant's avatar
Matthieu Constant committed
721
  printf("RatioEl = %e\n",ratioEl);
Matthieu Constant's avatar
Matthieu Constant committed
722
723
724
725
  for (int i = 0; i < mesh->n_nodes; ++i){
    new_mesh_size[i] = fmax(new_mesh_size[i]/cbrt(ratioEl),lcmin);
  }
#endif
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
726
  FILE *f = fopen("adapt/lc.pos", "w");
727
728
  if(!f)
    printf("cannot open file adapt/lc.pos for writing\n");
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
729
  fprintf(f, "View \"lc\" {\n");
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
730
  for (int iel = 0; iel < mesh->n_elements; ++iel) {
731
732
733
734
735
736
    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");
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
737
  }
738
  fprintf(f,"};\n");
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
739
740
  fclose(f);
  free(new_mesh_size);
Matthieu Constant's avatar
Matthieu Constant committed
741
  free(num_lc);
742
  system("gmsh -"xstr(DIMENSION)" adapt/mesh.geo");
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
743
744
745
}


Matthieu Constant's avatar
Matthieu Constant committed
746
void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin, double n_el)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
747
{
Matthieu Constant's avatar
Matthieu Constant committed
748
749
750
  struct timespec tic, toc;
  fluid_problem_adapt_gen_mesh(problem, lcmax, lcmin, n_el);
  clock_get(&tic);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
751
  Mesh *new_mesh = mesh_load_msh("adapt/mesh.msh");
Matthieu Constant's avatar
Matthieu Constant committed
752
  clock_get(&toc);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
753
  //printf("Time mesh_load : %.6e \n",clock_diff(&tic,&toc));
754
  double *new_solution = malloc(sizeof(double)*new_mesh->n_nodes*N_N);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
755
  double *new_xi = malloc(sizeof(double)*new_mesh->n_nodes*DIMENSION);
Matthieu Constant's avatar
Matthieu Constant committed
756
  clock_get(&tic);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
757
  int *new_eid = malloc(sizeof(int)*new_mesh->n_nodes);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
758
759
  for (int i = 0; i < new_mesh->n_nodes; ++i)
    new_eid[i] = -1;
Matthieu Constant's avatar
Matthieu Constant committed
760
  clock_get(&toc);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
761
  //printf("Time new_eid : %.6e \n",clock_diff(&tic,&toc));
Matthieu Constant's avatar
Matthieu Constant committed
762
  clock_get(&tic);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
763
764
765
766
  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];
Matthieu Constant's avatar
Matthieu Constant committed
767
  clock_get(&toc);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
768
  //printf("Time newx : %.6e \n",clock_diff(&tic,&toc));
Matthieu Constant's avatar
Matthieu Constant committed
769
  clock_get(&tic);
770
  mesh_tree_particle_to_mesh(problem->mesh_tree,new_mesh->n_nodes, newx, new_eid, new_xi);
Matthieu Constant's avatar
Matthieu Constant committed
771
  clock_get(&toc);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
772
  //printf("Time find new point in old mesh : %.6e \n",clock_diff(&tic,&toc));
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
773
  for (int i = 0; i < new_mesh->n_nodes; ++i) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
774
775
776
777
778
    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);
    }
779
780
781
782
783
    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)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
784
        new_solution[i*n_fields+j] += problem->solution[el[k]*n_fields+j]*phi[k];
785
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
786
787
788
789
790
791
  }
  free(new_eid);
  free(new_xi);
  free(newx);
  free(problem->solution);
  problem->solution = new_solution;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
792
793
  free(problem->new_porosity);
  problem->new_porosity = malloc(sizeof(double)*new_mesh->n_nodes);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
794
795
  free(problem->porosity);
  problem->porosity = malloc(sizeof(double)*new_mesh->n_nodes);
Matthieu Constant's avatar
Matthieu Constant committed
796
797
798
799
  for (int i = 0; i < new_mesh->n_nodes; ++i)
  {
    problem->porosity[i] = 1.;
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
800
801
  mesh_free(problem->mesh);
  problem->mesh = new_mesh;
802
803
  mesh_tree_free(problem->mesh_tree);
  problem->mesh_tree = mesh_tree_create(problem->mesh);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
804
805
  for (int i = 0; i < problem->n_particles; ++i)
    problem->particle_element_id[i] = -1;
806
  mesh_tree_particle_to_mesh(problem->mesh_tree, problem->n_particles, problem->particle_position, problem->particle_element_id, problem->particle_uvw);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
807
  fluid_problem_compute_node_volume(problem);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
808
  fluid_problem_compute_porosity(problem);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
809
810
  hxtLinearSystemDelete(&problem->linear_system);
  #ifdef HXT_HAVE_PETSC
Matthieu Constant's avatar
Matthieu Constant committed
811
    hxtLinearSystemCreatePETSc(&problem->linear_system,new_mesh->n_elements,N_N,n_fields,new_mesh->elements,"fluid");
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
812
  #else
813
    hxtLinearSystemCreateLU(&problem->linear_system,new_mesh->n_elements,N_N,n_fields,new_mesh->elements);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
814
  #endif 
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
815
816
}

Matthieu Constant's avatar
Matthieu Constant committed
817
818
819
820
void fluid_problem_after_import(FluidProblem *problem)
{
  mesh_tree_free(problem->mesh_tree);
  problem->mesh_tree = mesh_tree_create(problem->mesh);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
821
822
  free(problem->new_porosity);
  problem->new_porosity = malloc(sizeof(double)*problem->mesh->n_nodes);
Matthieu Constant's avatar
Matthieu Constant committed
823
824
825
826
827
828
829
830
831
832
833
834
835
  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);
  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,new_mesh->n_elements,N_N,n_fields,new_mesh->elements);
  #endif 
}

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
836
837

void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity, long int *elid)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
838
{
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
839
840
841
  
  struct timespec tic, toc;
  
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
842
843
844
845
846
847
  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);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
848
    free(problem->particle_force);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
849
850
    free(problem->particle_volume);
    free(problem->particle_velocity);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
851
    free(problem->particle_new_velocity);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
852
    problem->particle_uvw = malloc(sizeof(double)*DIMENSION*n);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
853
    problem->particle_element_id = malloc(sizeof(int)*n);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
854
855
856
857
    for (int i = 0; i < n; ++i) {
      problem->particle_element_id[i]=-1;
    }
    problem->particle_position = malloc(sizeof(double)*DIMENSION*n);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
858
    problem->particle_mass = malloc(sizeof(double)*n);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
859
    problem->particle_force = malloc(sizeof(double)*n*DIMENSION);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
860
    problem->particle_volume = malloc(sizeof(double)*n);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
861
    problem->particle_velocity = malloc(sizeof(double)*n*DIMENSION);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
862
    problem->particle_new_velocity = malloc(sizeof(double)*n*DIMENSION);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
863
864
865
866
867
  }
  if (elid) {
    for (int i = 0; i < n; ++i) {
      problem->particle_element_id[i]=elid[i];
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
868
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
869
  clock_get(&tic);
870
  mesh_tree_particle_to_mesh(problem->mesh_tree, n, position, problem->particle_element_id, problem->particle_uvw);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
871
  clock_get(&toc);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
872
  //printf("Time mesh_particles_to_mesh : %.6e \n",clock_diff(&tic,&toc));
Matthieu Constant's avatar
Matthieu Constant committed
873
  clock_get(&tic);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
874
875
876
  for (int i = 0; i < n; ++i) {
    problem->particle_mass[i] = mass[i];
    problem->particle_volume[i] = volume[i];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
877
878
879
880
    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];
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
881
  }
Matthieu Constant's avatar
Matthieu Constant committed
882
  clock_get(&toc);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
883
  //printf("Time initialize particle_fluid : %.6e \n",clock_diff(&tic,&toc));
Matthieu Constant's avatar
Matthieu Constant committed
884
  clock_get(&tic);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
885
  fluid_problem_compute_porosity(problem);
Matthieu Constant's avatar
Matthieu Constant committed
886
  clock_get(&toc);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
887
  //printf("Time compute_porosity : %.6e \n",clock_diff(&tic,&toc));
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
888
889
890
891
892
893
894
895
896
897
}

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

int fluid_problem_n_particles(FluidProblem *problem)
{
  return problem->n_particles;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
898
}