fluid_problem.c 22.4 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"
8
9
10
11
12

#define n_fields 3

#define newton_max_it 20
#define newton_atol 1e-12
13
#define newton_rtol 1e-5
14
15
16
17
18

#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};

19

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
20
int fluid_problem_export(const FluidProblem *problem, const char *output_dir, double t, int iter)
21
{
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
22
23
24
25
26
  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);
27
    return -1;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
28
29
30
31
  }
  if(mesh_write_msh(problem->mesh, f))
    return -1;
  if(mesh_write_msh_scalar(problem->mesh, f, "porosity", t, iter, problem->porosity, 1, 0))
32
    return -1;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
33
  if(mesh_write_msh_vector(problem->mesh, f, "uv", t, iter, problem->solution, n_fields, 0, 1))
34
    return -1;
Matthieu Constant's avatar
Matthieu Constant committed
35
36
   if(mesh_write_msh_vector(problem->mesh, f, "force", t, iter, problem->node_force, 2, 0, 1))
    return -1;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
37
38
39
  if(mesh_write_msh_scalar(problem->mesh, f, "p", t, iter, problem->solution, n_fields, 2))
    return -1;
  fclose(f);
40
41
42
  return 0;
}

43
static void particle_drag(double u[2], double mu, double rho, double d, double c, double drag[2], double dt, double mass, double rhop)
44
45
46
47
48
49
50
51
{
  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);
  double GoU = f*Cd_u*rho/2*d;
52
53
  drag[0] = GoU*u[0]/(1+(rhop+c/(1-c)*rho)*dt/(c/(1-c)*rho*mass)*GoU);
  drag[1] = GoU*u[1]/(1+(rhop+c/(1-c)*rho)*dt/(c/(1-c)*rho*mass)*GoU);
54
55
}

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
56
void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt, double *particle_force) {
57
58
  double *porosity = problem->porosity;
  Mesh *mesh = problem->mesh;
Matthieu Constant's avatar
Matthieu Constant committed
59
60
61
62
  const double *solution = problem->solution;
  for (int i = 0; i < mesh->n_nodes*2; ++i) {
    problem->node_force[i] = 0.0;
  }
63
64
65
66
  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];
67
68
69
70
71
72
73
74
    double vol = problem->particle_volume[i];
    double mass = problem->particle_mass[i];
    particle_force[i*2+0] = 0;
    particle_force[i*2+1] = -9.81*(mass-problem->rho*vol);
    if(iel < 0){
      continue;
    }

75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
    int *tri = problem->mesh->triangles+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
    };
110
    double rhop = mass/vol;
111
112
    double d = 2*sqrt(vol/M_PI);
    double drag[2];
113
    particle_drag(du,problem->mu,problem->rho,d,c,drag, dt, mass, rhop);
114
115
    particle_force[i*2+0] += -drag[0]-vol*gradp[0];
    particle_force[i*2+1] += -drag[1]-vol*gradp[1];
116
    for (int iphi = 0; iphi < 3; ++iphi) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
117
118
      problem->node_force[tri[iphi]*2+0] += drag[0]*phi[iphi];
      problem->node_force[tri[iphi]*2+1] += drag[1]*phi[iphi];
119
120
    }
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
121
122
123
124
  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];
  }
125
126
}

127
128
129
130
131
132
static void fluid_problem_assemble_system(FluidProblem *problem, const double *solution_old, double dt)
{
  LinearSystem *lsys = problem->linear_system;
  const Mesh *mesh = problem->mesh;
  const double *solution = problem->solution;
  const double *porosity = problem->porosity;
133
  const double *old_porosity = problem->old_porosity;
Matthieu Constant's avatar
Matthieu Constant committed
134
  const double *node_force = problem->node_force;
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
  double *local_vector = malloc(sizeof(double)*n_fields*3);
  double *local_matrix = malloc(sizeof(double)*n_fields*n_fields*9);
  const double mu = problem->mu;
  const double rho = problem->rho;
  const double epsilon = problem->epsilon;
  for (int iel=0; iel < mesh->n_triangles; ++iel) {
    const int *tri = &mesh->triangles[iel*3];
    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]}
    };
159
160
161
162
    for (int i=0; i < 3*n_fields; ++i)
      local_vector[i] = 0;
    for (int i=0; i < 9*n_fields*n_fields; ++i)
      local_matrix[i] = 0;
163
    double C[] = {porosity[tri[0]],porosity[tri[1]],porosity[tri[2]]};
164
165
    double OLDC[] = {old_porosity[tri[0]],old_porosity[tri[1]],old_porosity[tri[2]]};
    double DCDT[] = {(C[0]-OLDC[0])/dt,(C[1]-OLDC[1])/dt,(C[2]-OLDC[2])/dt};
166
167
    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]};
Matthieu Constant's avatar
Matthieu Constant committed
168
169
    double FU[] = {node_force[tri[0]*2+0], node_force[tri[1]*2+0], node_force[tri[2]*2+0]};
    double FV[] = {node_force[tri[0]*2+1], node_force[tri[1]*2+1], node_force[tri[2]*2+1]};
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
    double U_old[] = {solution_old[tri[0]*3+0],solution_old[tri[1]*3+0],solution_old[tri[2]*3+0]};
    double V_old[] = {solution_old[tri[0]*3+1],solution_old[tri[1]*3+1],solution_old[tri[2]*3+1]};
    double P[] = {solution[tri[0]*3+2],solution[tri[1]*3+2],solution[tri[2]*3+2]};
    double dc[2] = {
      dphi[0][0]*C[0]+dphi[1][0]*C[1]+dphi[2][0]*C[2],
      dphi[0][1]*C[0]+dphi[1][1]*C[1]+dphi[2][1]*C[2]
    };
    double du[2] = {
      dphi[0][0]*U[0]+dphi[1][0]*U[1]+dphi[2][0]*U[2],
      dphi[0][1]*U[0]+dphi[1][1]*U[1]+dphi[2][1]*U[2]
    };
    double dv[2] = {
      dphi[0][0]*V[0]+dphi[1][0]*V[1]+dphi[2][0]*V[2],
      dphi[0][1]*V[0]+dphi[1][1]*V[1]+dphi[2][1]*V[2]
    };
    double dp[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]
    };
    for (int iqp = 0; iqp< N_QP; ++iqp) {
      const double xi = QP[iqp][0], eta = QP[iqp][1];
      const double phi[] = {1-xi-eta, xi, eta};
      double u = phi[0]*U[0]+phi[1]*U[1]+phi[2]*U[2];
      double v = phi[0]*V[0]+phi[1]*V[1]+phi[2]*V[2];
Matthieu Constant's avatar
Matthieu Constant committed
194
195
      double fu = phi[0]*FU[0]+phi[1]*FU[1]+phi[2]*FU[2];
      double fv = phi[0]*FV[0]+phi[1]*FV[1]+phi[2]*FV[2];
196
197
198
199
      double u_old = phi[0]*U_old[0]+phi[1]*U_old[1]+phi[2]*U_old[2];
      double v_old = phi[0]*V_old[0]+phi[1]*V_old[1]+phi[2]*V_old[2];
      double p = phi[0]*P[0]+phi[1]*P[1]+phi[2]*P[2];
      double c = phi[0]*C[0]+phi[1]*C[1]+phi[2]*C[2];
200
      double dcdt = phi[0]*DCDT[0]+phi[1]*DCDT[1]+phi[2]*DCDT[2];
201
202
203
204
205
206
207
      const double jw = QW[iqp]*detj;
      for (int iphi = 0; iphi < 3; ++iphi) {
        double phii = phi[iphi];
        double dphii[2] = {dphi[iphi][0],dphi[iphi][1]};
        double tau[2][2] = {{du[0]-u*dc[0]/c,du[1]-u*dc[1]/c},
                            {dv[0]-v*dc[0]/c,dv[1]-v*dc[1]/c}};
        local_vector[iphi+0] += jw*(
Matthieu Constant's avatar
Matthieu Constant committed
208
           2*mu*(dphii[0]*tau[0][0]+dphii[1]*(tau[0][1]+tau[1][0])*0.5)
209
          +rho*phii*((u-u_old)/dt+(u*tau[0][0]+v*tau[0][1])/c)
Matthieu Constant's avatar
Matthieu Constant committed
210
          -p*dphii[0]-fu*phii
211
212
        );
        local_vector[iphi+3] += jw*(
Matthieu Constant's avatar
Matthieu Constant committed
213
           2*mu*(dphii[0]*(tau[0][1]+tau[1][0])*0.5+dphii[1]*tau[1][1])
214
          +rho*phii*((v-v_old)/dt+(u*tau[1][0]+v*tau[1][1])/c)
Matthieu Constant's avatar
Matthieu Constant committed
215
          -p*dphii[1]-fv*phii
216
        );
Matthieu Constant's avatar
Matthieu Constant committed
217

218
219
        local_vector[iphi+6] += jw*(
           epsilon*(dphii[0]*dp[0]+dphii[1]*dp[1])
220
           +phii*((du[0]+dv[1])+dcdt)
221
222
223
224
225
        );
        for (int jphi = 0; jphi < 3; ++jphi) {
          double phij = phi[jphi];
          double dphij[2] = {dphi[jphi][0],dphi[jphi][1]};
          double dtau[2] = {dphij[0]-phij*dc[0]/c,dphij[1]-phij*dc[1]/c};
Matthieu Constant's avatar
Matthieu Constant committed
226
          local_matrix[(iphi+0)*9+jphi+0] += jw*2*mu*(dphii[0]*dtau[0]+dphii[1]*dtau[1]*0.5)  // U U
227
                                            +jw*rho*phii*(phij/dt+(phij*tau[0][0]+u*dtau[0]+v*dtau[1])/c);
Matthieu Constant's avatar
Matthieu Constant committed
228
          local_matrix[(iphi+0)*9+jphi+3] += jw*2*mu*dphii[1]*dtau[0]*0.5                     // U V
229
230
                                            +jw*rho*phii*phij*tau[0][1]/c;
          local_matrix[(iphi+0)*9+jphi+6] +=-jw*dphii[0]*phij;                                // U P
Matthieu Constant's avatar
Matthieu Constant committed
231
          local_matrix[(iphi+3)*9+jphi+0] += jw*2*mu*dphii[0]*dtau[1]*0.5                     // V U
232
                                            +jw*rho*phii*phij*tau[1][0]/c;
Matthieu Constant's avatar
Matthieu Constant committed
233
          local_matrix[(iphi+3)*9+jphi+3] += jw*2*mu*(dphii[0]*dtau[0]*0.5+dphii[1]*dtau[1])  // V V
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
                                            +jw*rho*phii*(phij/dt+(phij*tau[1][1]+u*dtau[0]+v*dtau[1])/c);
          local_matrix[(iphi+3)*9+jphi+6] +=-jw*dphii[1]*phij;                                // V P
          local_matrix[(iphi+6)*9+jphi+0] += jw*phii*dphij[0];                                // P U
          local_matrix[(iphi+6)*9+jphi+3] += jw*phii*dphij[1];                                // P V
          local_matrix[(iphi+6)*9+jphi+6] += jw*epsilon*(dphii[0]*dphij[0]+dphii[1]*dphij[1]);// P P
        }
      }
    }
    linear_system_add_to_matrix(lsys,iel,iel,local_matrix);
    linear_system_add_to_rhs(lsys,iel,local_vector);
  }
  free(local_vector);
  free(local_matrix);
}

249
static void apply_boundary_conditions(const Mesh *mesh, double *solution, int nBnd, StrongBoundary *bnds)
250
{
251
252
  for (int ibnd = 0; ibnd < nBnd; ++ibnd) {
    StrongBoundary *bnd = bnds + ibnd;
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
    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)
      solution[mesh->phys_nodes[iphys][i]*n_fields+bnd->field] = bnd->v;
  }
}

void fluid_problem_implicit_euler(FluidProblem *problem, double dt)
{
268
269
  static double totaltime = 0;
  struct timespec tic, toc;
270
  const Mesh *mesh = problem->mesh;
271
  apply_boundary_conditions(mesh, problem->solution, problem->n_strong_boundaries, problem->strong_boundaries);
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
  double *solution_old = 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;
  for (int i=0; i<newton_max_it; ++i) {
    linear_system_zero_matrix_and_rhs(problem->linear_system);
    fluid_problem_assemble_system(problem, solution_old, dt);
    double norm = linear_system_get_rhs_norm(problem->linear_system);
    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;
288
  clock_get(&tic);
289
    linear_system_solve(problem->linear_system, dx);
290
291
  clock_get(&toc);
  totaltime += clock_diff(&tic, &toc);
292
293
294
295
    for (int j=0; j<mesh->n_nodes*n_fields; ++j) {
      problem->solution[j] -= dx[j];
    }
  }
296
  printf("total solve time : %g\n", totaltime);
297
298
299
300
  free(dx);
  free(solution_old);
}

Matthieu Constant's avatar
Matthieu Constant committed
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
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_triangles; ++iel) {
    const int *tri = &mesh->triangles[iel*3];
    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];
    problem->node_volume[tri[0]] += detj/6;
    problem->node_volume[tri[1]] += detj/6;
    problem->node_volume[tri[2]] += detj/6;
  }
}

324
FluidProblem *fluid_problem_new(const char *mesh_file_name, double mu, double rho, double epsilon, int n_strong_boundaries, StrongBoundary *strong_boundaries)
325
326
327
328
329
330
331
332
333
{
  FluidProblem *problem = malloc(sizeof(FluidProblem));
  problem->epsilon = epsilon;
  problem->rho = rho;
  problem->mu = mu;
  Mesh *mesh = mesh_load_msh(mesh_file_name);
  if (!mesh)
    return NULL;
  problem->mesh = mesh;
334
335
336
337
338
339
340
  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].v = strong_boundaries[i].v;
    problem->strong_boundaries[i].field = strong_boundaries[i].field;
  }
341
  problem->porosity = malloc(mesh->n_nodes*sizeof(double));
Matthieu Constant's avatar
Matthieu Constant committed
342
343
344
345
  problem->node_force = malloc(mesh->n_nodes*2*sizeof(double));
  for (int i = 0; i < mesh->n_nodes*2; ++i) {
    problem->node_force[i] = 0.0;
  }
346
  problem->old_porosity = malloc(mesh->n_nodes*sizeof(double));
347
348
  problem->solution = malloc(mesh->n_nodes*n_fields*sizeof(double));
  // begin to remove
349
350
  //static StrongBoundary strong_boundaries[] = {{"Bottom",1,VEXT},{"Bottom",0,0.},{"Top",0,0.},{"PtFix",2,0.},{"Box",0,0.},{NULL, 0, 0.}};
  //static StrongBoundary strong_boundaries[] = {{"Top",2,0},{"Top",1,0.},{"Top",0,0.},{"Box",1,0.},{"Box",0,0.},{NULL, 0, 0.}};
351
352
353
354
355
  for (int i = 0; i < problem->mesh->n_nodes; ++i){
    problem->solution[i*3+0] = 0.;
    problem->solution[i*3+1] = 0;
    problem->solution[i*3+2] = 0.;
  }
Matthieu Constant's avatar
Matthieu Constant committed
356
357
  problem->node_volume = malloc(0);
  fluid_problem_compute_node_volume(problem);
358
  // end to remove
359
360
361
362
363
364
365
366
  problem->linear_system = linear_system_new(mesh, n_fields, problem->n_strong_boundaries, problem->strong_boundaries); 
  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);
367
368
369
370
371
372
  return problem;
}

void fluid_problem_free(FluidProblem *problem)
{
  free(problem->solution);
Matthieu Constant's avatar
Matthieu Constant committed
373
  free(problem->node_force);
374
  free(problem->porosity);
375
  free(problem->old_porosity);
376
  linear_system_free(problem->linear_system);
377
378
379
380
381
382
383
384
385
  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);
386
387
388
  mesh_free(problem->mesh);
}

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
389
static void fluid_problem_compute_porosity(FluidProblem *problem)
390
{
391
  Mesh *mesh = problem->mesh;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
392
  double *volume = problem->particle_volume;
393
394
395
  for (int i = 0; i < mesh->n_nodes; ++i){
    problem->porosity[i] = 0;
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
396
  for (int i = 0; i < problem->n_particles; ++i) {
397
    if(problem->particle_element_id[i] == -1) continue;
398
399
400
401
402
403
404
    const int *tri = &mesh->triangles[problem->particle_element_id[i]*3];
    double u = problem->particle_uvw[i*2+0], v = problem->particle_uvw[i*2+1];
    problem->porosity[tri[0]] += (1-u-v)*volume[i];
    problem->porosity[tri[1]] += u*volume[i];
    problem->porosity[tri[2]] += v*volume[i];
  }
  for (int i = 0; i < mesh->n_nodes; ++i){
Matthieu Constant's avatar
Matthieu Constant committed
405
    problem->porosity[i] = (1-problem->porosity[i]/problem->node_volume[i]);
406
  }
407
}
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
408
409
410
411
412
413
414
415
416
417

void fluid_problem_adapt_mesh(FluidProblem *problem, double gradmin, double gradmax, double lcmin)
{
  Mesh *mesh = problem->mesh;
  double *solution = problem->solution;
  double *new_mesh_size = malloc(sizeof(double)*mesh->n_nodes);
  for (int i = 0; i < mesh->n_nodes; ++i)
    new_mesh_size[i] = DBL_MAX;
  gradmax = sqrt(gradmax);
  gradmin = sqrt(gradmin);
Matthieu Constant's avatar
Matthieu Constant committed
418
  double *porosity = problem->porosity;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
419
420
  for (int iel = 0; iel < mesh->n_triangles; ++iel) {
    int *tri = problem->mesh->triangles+iel*3;
Matthieu Constant's avatar
Matthieu Constant committed
421
    double C[] = {porosity[tri[0]], porosity[tri[1]], porosity[tri[2]]};
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
    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]};
    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 gradU[2] = {
Matthieu Constant's avatar
Matthieu Constant committed
442
443
      dphi[0][0]*U[0]/C[0]+dphi[1][0]*U[1]/C[1]+dphi[2][0]*U[2]/C[2],
      dphi[0][1]*U[0]/C[0]+dphi[1][1]*U[1]/C[1]+dphi[2][1]*U[2]/C[2]
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
444
445
    };
    double gradV[2] = {
Matthieu Constant's avatar
Matthieu Constant committed
446
447
      dphi[0][0]*V[0]/C[0]+dphi[1][0]*V[1]/C[1]+dphi[2][0]*V[2]/C[2],
      dphi[0][1]*V[0]/C[0]+dphi[1][1]*V[1]/C[1]+dphi[2][1]*V[2]/C[2]
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
448
449
450
451
452
453
454
    };
    double ngrad = pow(gradU[0]*gradU[0]+gradU[1]*gradU[1] + gradV[0]*gradV[0] + gradV[1]*gradV[1], 0.25);
    double lc = lcmin /fmin(1,fmax(ngrad/gradmax,gradmin/gradmax));
    for (int j = 0; j < 3; ++j)
      new_mesh_size[tri[j]] = fmin(new_mesh_size[tri[j]], lc);
  }
  FILE *f = fopen("adapt/lc.pos", "w");
455
456
457
458
  if(!f)
    printf("cannot open file adapt/lc.pos for writing\n");
  else
    printf("open ok\n");
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
459
  fprintf(f, "View \"lc\" {\n");
460
    printf("write ok\n");
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
  for (int iel = 0; iel < mesh->n_triangles; ++iel) {
    int *tri = problem->mesh->triangles+iel*3;
    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]}};
    double LC[3] = {new_mesh_size[tri[0]], new_mesh_size[tri[1]], new_mesh_size[tri[2]]};
    fprintf(f, "ST(%.8g, %.8g, 0, %.8g, %.8g, 0, %.8g, %.8g, 0){%.8g, %.8g, %.8g};\n",
        x[0][0],x[0][1],x[1][0],x[1][1],x[2][0],x[2][1],
        LC[0],LC[1],LC[2]);
  }
  fprintf(f, "};\n");
  fclose(f);
  free(new_mesh_size);
  system("gmsh -2 adapt/mesh.geo");

  Mesh *new_mesh = mesh_load_msh("adapt/mesh.msh");
  double *new_solution = malloc(sizeof(double)*new_mesh->n_nodes*3);
  double *new_xi = malloc(sizeof(double)*new_mesh->n_nodes*2);
  int *new_eid = malloc(sizeof(int)*new_mesh->n_nodes);
  double *newx = malloc(sizeof(double)*2*new_mesh->n_nodes);
  for (int i = 0;i < new_mesh->n_nodes;++i) {
    newx[i*2+0] = new_mesh->x[i*3+0];
    newx[i*2+1] = new_mesh->x[i*3+1];
  }
  mesh_particles_to_mesh(mesh, new_mesh->n_nodes, newx, new_eid, new_xi);
  for (int i = 0; i < new_mesh->n_nodes; ++i) {
    int *tri = problem->mesh->triangles+new_eid[i]*3;
489
490
    if(new_eid[i] < 0)
      continue;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
    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 *xi = new_xi+i*2;
    double phi[] = {1-xi[0]-xi[1],xi[0],xi[1]};
    new_solution[i*3+0] = U[0]*phi[0]+U[1]*phi[1]+U[2]*phi[2];
    new_solution[i*3+1] = V[0]*phi[0]+V[1]*phi[1]+V[2]*phi[2];
    new_solution[i*3+2] = 0.;
  }
  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);
Matthieu Constant's avatar
Matthieu Constant committed
508
509
  free(problem->node_force);
  problem->node_force = malloc(sizeof(double)*2*new_mesh->n_nodes);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
510
511
512
  mesh_free(problem->mesh);
  problem->mesh = new_mesh;
  mesh_particles_to_mesh(problem->mesh, problem->n_particles, problem->particle_position, problem->particle_element_id, problem->particle_uvw);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
513
  fluid_problem_compute_node_volume(problem);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
514
515
  fluid_problem_compute_porosity(problem);
  linear_system_free(problem->linear_system);
Matthieu Constant's avatar
Matthieu Constant committed
516
  problem->linear_system = linear_system_new(problem->mesh, n_fields, problem->n_strong_boundaries, problem->strong_boundaries);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
}

void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity)
{
  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)*2*n);
    problem->particle_element_id = malloc(sizeof(int)*n);
    problem->particle_position = malloc(sizeof(double)*2*n);
    problem->particle_mass = malloc(sizeof(double)*n);
    problem->particle_volume = malloc(sizeof(double)*n);
    problem->particle_velocity = malloc(sizeof(double)*n*2);
  }
  mesh_particles_to_mesh(problem->mesh, n, position, problem->particle_element_id, problem->particle_uvw);
  for (int i = 0; i < n; ++i) {
    problem->particle_position[i*2+0] = position[i*2+0];
    problem->particle_position[i*2+1] = position[i*2+1];
    problem->particle_mass[i] = mass[i];
    problem->particle_volume[i] = volume[i];
    problem->particle_velocity[i*2+0] = velocity[i*2+0];
    problem->particle_velocity[i*2+1] = velocity[i*2+1];
  }
  for (int i = 0; i < problem->mesh->n_nodes; ++i){
    problem->old_porosity[i] = problem->porosity[i];
  }
  fluid_problem_compute_porosity(problem);
}