fluid_problem.c 31.2 KB
Newer Older
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
1
2
/*
 * Marblesbag - Copyright (C) <2010-2018>
3
4
 * <Universite catholique de Louvain (UCL), Belgium
 *  Universite de Montpellier, France>
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
5
 * 	
6
 * List of the contributors to the development of Marblesbag: see AUTHORS file.
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
7
8
9
 * Description and complete License: see LICENSE file.
 * 	
 * This program (Marblesbag) is free software: 
10
 * you can redistribute it and/or modify it under the terms of the GNU Lesser General 
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
11
12
13
14
15
16
 * Public License as published by the Free Software Foundation, either version
 * 3 of the License, or (at your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
17
 * GNU Lesser General Public License for more details.
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
18
 * 
19
20
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program (see COPYING and COPYING.LESSER files).  If not, 
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
21
22
23
 * see <http://www.gnu.org/licenses/>.
 */

24
#include "tools.h"
25
26
#include <stdlib.h>
#include <stdio.h>
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
27
#include <float.h>
28
29
#include <math.h>
#include "fluid_problem.h"
30
#include "mesh_find.h"
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
31
#include "mbxml.h"
32

33
#define n_fields (D+2)
34

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
35
#define newton_max_it 10
Matthieu Constant's avatar
Matthieu Constant committed
36
#define newton_atol 5e-7
37
#define newton_rtol 1e-5
38
#define D DIMENSION
39

40
#if D==2
41
42
43

#define N_SF 3
#define N_N 3
44
45
46
#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};
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
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
74
#else
75
76
77

#define N_SF 4
#define N_N 4
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
78
79
80
81
82
83
84
85
86
#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};
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
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
123
124
#endif

125
#if D==2
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
126
static double particle_drag_coeff(double u[2], double mu, double rho, double vol, double c)
127
{
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
128
  double d = 2*sqrt(vol/M_PI);
129
130
  double normu = hypot(u[0],u[1]);
  //Reynolds/|u_p-u_f|
Matthieu Constant's avatar
Matthieu Constant committed
131
  double Re_O_u = sqrt(d)*c/mu*rho;
132
133
134
  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
135
  return f*Cd_u*rho/2*d;
136
}
137
#else
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
138
static double particle_drag_coeff(double u[3], double mu, double rho, double vol, double c)
139
{
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
140
  double d = 2*cbrt(3./4.*vol/M_PI);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
141
  double normu = sqrt(u[0]*u[0]+u[1]*u[1]+u[2]*u[2]);
142
143
144
145
146
147
  //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
148
  return f*Cd_u*rho/2*(M_PI*r*r);
149
150
}
#endif
151

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
152
void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt, double *particle_force) {
153
  for (int i = 0; i < problem->n_particles*D; ++i) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
154
155
156
157
158
159
160
161
162
    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];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
163
    double mass = problem->particle_mass[i];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
164
    if(iel < 0){
165
166
167
      for (int d = 0; d < D; ++d)
        problem->particle_force[i*D+d] = 0;
      problem->particle_force[i*D+1] = problem->g*mass;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
168
169
      continue;
    }
170
    double *xi = problem->particle_uvw + D*i;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
171
172
173
174
    double phi[N_SF];
    shape_functions(xi,phi);
    uint32_t *el=problem->mesh->elements+iel*N_SF;
    double *se = problem->solution_explicit;
175
176
177
178
    double c=0, vf[D]={0}, vfe[D]={0}, gradp[D]={0};
    double dxdxi[D][D],dxidx[D][D];
    for (int k=0; k<D; ++k)
      for (int j=0; j<D; ++j)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
179
180
        dxdxi[k][j] = mesh->x[el[j+1]*3+k]-mesh->x[el[0]*3+k];
    invDD(dxdxi,dxidx);
181
    double dphi[N_SF][D];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
182
    grad_shape_functions(dxidx, dphi);
183
    double *si = problem->solution;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
184
185
    for (int iphi=0; iphi < N_SF; ++iphi) {
      c += porosity[el[iphi]]*phi[iphi];
186
      for (int j=0; j < D; ++j) {
187
        vf[j]  += si[el[iphi]*n_fields+j]*phi[iphi];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
188
        vfe[j] += se[el[iphi]*n_fields+j]*phi[iphi];
189
        gradp[j] += dphi[iphi][j]*si[el[iphi]*n_fields+D];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
190
191
      }
    }
192
193
194
195
    double du[D], due[D];
    for (int j = 0; j < D; ++j) {
      du[j] = problem->particle_velocity[i*D+j]-vf[j]/c;
      due[j] = problem->particle_velocity[i*D+j]-vfe[j]/c;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
196
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
197
    double vol = problem->particle_volume[i];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
198
199
    double rhop = mass/vol;
    double drho = rhop -problem->rho;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
200
    double g = problem->g*drho/rhop;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
201
    double gamma = particle_drag_coeff(due,problem->mu,problem->rho,vol,c);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
202
203
    double fcoeff = mass/(mass+dt*gamma);
    
204
205
206
    for (int d = 0; d < D; ++d)
      problem->particle_force[i*D+d] = fcoeff*(-gamma*du[d]-vol*gradp[d]);
    problem->particle_force[i*D+1] += fcoeff*g*mass;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
207
    
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
208
209
    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
210
    for (int iphi = 0; iphi < N_SF; ++iphi) {
211
      for (int d = 0; d < D; ++d)
212
213
214
        local_vector[iphi+N_SF*d] -= fcoeff*gamma*(du[d]-dt/mass*vol*gradp[d])*phi[iphi];
      local_vector[iphi+N_SF*1] -= fcoeff*gamma*dt*g*phi[iphi];
      
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
215
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
216
217
    for (int iphi = 0; iphi < N_SF; ++iphi) {
      for (int jphi = 0; jphi < N_SF; ++jphi) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
218
219
220
        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
221
        double f = fcoeff*phi[iphi];
222
        for (int d = 0; d < D; ++d){
223
          LOCAL_MATRIX(d,d) -= -f/c*phi[jphi]*gamma;
224
          LOCAL_MATRIX(d,D) -= -f*gamma*dt/mass*vol*dphi[jphi][d];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
225
        }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
226
227
      }
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
228
229
230
  }
}

231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
// source term
static void f_u_0(
    double f[D], const double u[D], const double du[D][D], const double q, const double dq[D],
    const double *uold, const double dt, const double rho) {
  double divu = 0;
  for (int j = 0; j < D; ++j) {
    divu += du[j][j];
  }
  for (int j = 0; j < D; ++j) {
    double utau = 0;
    for (int k = 0; k < D; ++k) {
      utau += u[k]*(du[j][k]-u[j]*dq[k]/q);
    }
    f[j] = rho*((u[j]-uold[j])/dt + u[j]/q*divu + utau/q);
  }
}


//flux
static void f_u_1(double f[D][D], double u[D], const double du[D][D], const double q, const double dq[D], const double p, const double mu) {
  double tau[D][D];
  for (int i = 0; i < D; ++i) {
    for (int j = 0; j < D; ++j) {
      tau[i][j] = du[i][j]-u[i]*dq[j]/q;
    }
  }
  for (int j = 0; j < D; ++j) {
    for (int k = 0; k < D; ++k) {
      f[j][k] =  mu*(tau[k][j]+tau[j][k]) + (k==j ? -p : 0);
    }
  }
}

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
264
265
266
267
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;
268
269
  //a generaliser
  double deps = 1e-8;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
270
271
272
273
274
  const double *solution = problem->solution;
  const double *porosity = problem->porosity;
  const double *old_porosity = problem->old_porosity;
  const double mu = problem->mu;
  const double rho = problem->rho;
Matthieu Constant's avatar
Matthieu Constant committed
275
  const double epsilon = problem->epsilon;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
276
277
278
279
280
281
282
283
284
  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
285
  for (int iel=0; iel < mesh->n_elements; ++iel) {
286
    const unsigned int *el = &mesh->elements[iel*N_N];
287
288
289
    double dxdxi[D][D], dxidx[D][D], dphi[N_N][D];
    for (int i = 0; i < D; ++i)
      for (int j = 0; j < D; ++j)
290
291
292
        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
293
294
    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];
295
    double dc[D]={0}, du[D][D]={{0}}, dp[D]={0};
296
    for (int i = 0; i< N_SF; ++i) {
297
      for (int j = 0; j < D; ++j) {
298
        dc[j] += dphi[i][j]*porosity[el[i]];
299
300
        dp[j] += dphi[i][j]*solution[el[i]*n_fields+D];
        for (int k = 0; k < D; ++k)
301
302
303
          du[k][j] += dphi[i][j]*solution[el[i]*n_fields+k];
      }
    }
304
305
    int P = D;
    int Q = D+1;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
306
    for (int iqp = 0; iqp< N_QP; ++iqp) {
307
308
      double phi[N_SF];
      shape_functions(QP[iqp],phi);
309
310
      double u[D]={0}, dudt[D]={0}, p=0, c=0, dcdt=0,q=0;
      
311
      for (int i = 0; i < N_SF; ++i) {
312
        for (int j = 0; j < D; ++j){
313
314
315
          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];
        }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
316
317
        p += solution[el[i]*n_fields+P]*phi[i];
        q +=  solution[el[i]*n_fields+Q]*phi[i];
318
        c += porosity[el[i]]*phi[i];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
319
320
        dcdt += (solution[el[i]*n_fields+Q]-solution_old[el[i]*n_fields+Q])/dt*phi[i];
        //dcdt += (problem->porosity[el[i]]-problem->old_porosity[el[i]])/dt*phi[i];
321
      }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
322
      const double jw = QW[iqp]*detj;
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356

      double fu0[D],fu1[D][D];
      f_u_0(fu0,u,du,q,dq,uold,dt,rho);
      f_u_1(fu1,u,du,q,dq,p,mu);

      double fu0_u[D][D],fu1_u[D][D][D],fu0_du[D][D][D],fu1_du[D][D][D][D];
      //fu0_u[f][u],fu1_u[f][u][x],fu0_du[f][u][sfx],fu1_du[f][u][x][sfx]
      for (int i = 0; i < D; ++i){
        double buf = u[i];
        u[i] += deps;
        f_u_0(fu0e,u,du,q,dq,uold,dt,rho);
        f_u_1(fu1e,u,du,q,dq,uold,dt,rho);
        u[i] = buf;
        for (int j = 0; j < D; ++j){
          fu0_u[j][i] = (fu0e[j]-fu0[j])/deps;
          for (int k = 0; k < D; ++k){
            fu1_u[j][i][k] = (fu1e[j][k]-fu1[j][k])/deps;
          }
        }
        for (int j = 0; j < D; ++j){
          buf = du[i][j];
          du[i][j] += deps;
          f_u_0(fu0e,u,du,q,dq,uold,dt,rho);
          f_u_1(fu1e,u,du,q,dq,p,mu);
          du[i][j] = buf;
          for (int k = 0; k < D; ++k){
            fu0_du[k][i][j] = (fu0e[k]-fu0[k])/deps;
            for (int l = 0; l < D; ++l){
              fu1_du[k][i][j][l] = (fu1e[k][l]-fu1[k][l])/deps;
            }
          }
        }
      }
      for (int iphi = 0; iphi < N_SF; ++iphi){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
357
        double phii = phi[iphi];
358
        const double *dphii = dphi[iphi];
359
360
361
362
        for (int j = 0; j < D; ++j){
          double f = fu0[j]*phii;
          for (int k = 0; k < D; ++k){
            f += fu1[j][k]*dphii[k];
363
          }
364
          local_vector[iphi+N_SF*j] += jw*f;
365
        }
366
367

        for (int jphi = 0; jphi < N_SF; ++jphi){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
368
          double phij = phi[jphi];
369
          const double *dphij = dphi[jphi];
370

371
372
          int N2 = n_fields*N_SF;
          int IDX = (iphi*N2+jphi)*n_fields;
373
#define LOCAL_MATRIX(a,b) local_matrix[IDX+N2*a+b]
374
          for (int j = 0; j < D; ++j){
375
            int U = j;
376
            for (int k = 0; k < D; ++k){
377
              int V = k;
378
379
380
381
382
383
384
385
              double m = fu0_u[j][k]*phii*phij;
              for (int i = 0; i < D; ++i){
                m += fu0_du[j][k][i]*dphij[i]*phii+fu1_u[j][k][i]*phij*dphii[i];
                for (int l = 0; l < D; ++l){
                  m += fu1_du[j][k][i][l]*dphij[i]*dphii[l];
                }
              }
              LOCAL_MATRIX(U,V) = jw*m;
386
387
            }
          }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
388
389
390
391
392
        }
      }
    }
    hxtLinearSystemAddToMatrix(lsys,iel,iel,local_matrix);
    hxtLinearSystemAddToRhs(lsys,rhs,iel,local_vector);
393
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
394
395
  free(all_local_matrix);
  free(all_local_vector);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
396

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
397
398
399
400
401
402
403
404
405
406
  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
407
408
      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
409
410
    }
  }
411
}
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
412
413


414

415
static void apply_boundary_conditions(const Mesh *mesh, double *solution, int nBnd, StrongBoundary *bnds)
416
{
417
418
  for (int ibnd = 0; ibnd < nBnd; ++ibnd) {
    StrongBoundary *bnd = bnds + ibnd;
419
420
421
422
423
424
425
426
    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);
    }
427
    double *x = malloc(mesh->phys_n_nodes[iphys]*sizeof(double)*D);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
428
429
430
    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];
431
432
      for (int k = 0; k < D; ++k)
        x[i*D+k] = mesh->x[j*3+k];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
433
434
    }
    bnd->apply(mesh->phys_n_nodes[iphys], x, v);
435
    for (int i = 0; i < mesh->phys_n_nodes[iphys]; ++i)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
436
437
438
      solution[mesh->phys_nodes[iphys][i]*n_fields+bnd->field] = v[i];
    free(x);
    free(v);
439
440
441
  }
}

442
int fluid_problem_implicit_euler(FluidProblem *problem, double dt)
443
{
444
445
  static double totaltime = 0;
  struct timespec tic, toc;
446
  const Mesh *mesh = problem->mesh;
447
  apply_boundary_conditions(mesh, problem->solution, problem->n_strong_boundaries, problem->strong_boundaries);
448
  double *solution_old = malloc(mesh->n_nodes*n_fields*sizeof(double));
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
449
  double *rhs = malloc(mesh->n_nodes*n_fields*sizeof(double));
450
451
  for (int i=0; i<mesh->n_nodes*n_fields; ++i)
    solution_old[i] = problem->solution[i];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
452
  problem->solution_explicit = solution_old;
453
454
  double *dx = malloc(sizeof(double)*mesh->n_nodes*n_fields);
  double norm0 = 0;
455
456
  int i;
  for (i=0; i<newton_max_it; ++i) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
457
458
459
460
461
462
    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);
463
464
465
466
467
468
469
    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;
470
  clock_get(&tic);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
471
    hxtLinearSystemSolve(problem->linear_system, rhs, dx);
472
473
  clock_get(&toc);
  totaltime += clock_diff(&tic, &toc);
474
475
476
477
    for (int j=0; j<mesh->n_nodes*n_fields; ++j) {
      problem->solution[j] -= dx[j];
    }
  }
478
  printf("total solve time : %g\n", totaltime);
479
  free(dx);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
480
  free(rhs);
481
  free(solution_old);
482
  return i;
483
484
}

Matthieu Constant's avatar
Matthieu Constant committed
485
486
487
488
489
490
491
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
492
  for (int iel = 0; iel < mesh->n_elements; ++iel) {
493
    double dxdxi[D][D];
494
    const int *el = &mesh->elements[iel*N_N];
495
496
    for (int i=0; i<D; ++i)
      for (int j=0; j<D; ++j)
497
        dxdxi[i][j] = mesh->x[el[j+1]*3+i]-mesh->x[el[0]*3+i];
498
#if D == 2
499
    const double vol = detDD(dxdxi)/6;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
500
#else
501
    const double vol = detDD(dxdxi)/24;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
502
#endif
503
504
    for (int i = 0; i < N_N; ++i)
      problem->node_volume[el[i]] += vol;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
505
506
507
508
509
510
511
512
513
514
515
516
  }
}

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;
517
    double sf[N_SF];
518
    shape_functions(&problem->particle_uvw[i*D],sf);
519
    const uint32_t *el = &mesh->elements[problem->particle_element_id[i]*N_N];
520
521
    for (int j=0; j<N_N;++j)
      problem->porosity[el[j]] += sf[j]*volume[i];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
522
523
  }
  for (int i = 0; i < mesh->n_nodes; ++i){
Matthieu Constant's avatar
Matthieu Constant committed
524
525
526
527
    if(problem->node_volume[i]==0.){
      problem->porosity[i] = 1.;
    }
    else{
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
528
      problem->porosity[i] = (1-problem->porosity[i]/problem->node_volume[i]);
Matthieu Constant's avatar
Matthieu Constant committed
529
    }
Matthieu Constant's avatar
Matthieu Constant committed
530
531
532
  }
}

Matthieu Constant's avatar
Matthieu Constant committed
533
FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double mu, double rho, double epsilon, int n_strong_boundaries, StrongBoundary *strong_boundaries)
534
{
535
536
537
538
539
  static int initialized = 0;
  if (!initialized) {
    hxtInitializeLinearSystems(0, NULL);
    initialized = 1;
#ifdef HXT_HAVE_PETSC
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
540
    hxtPETScInsertOptions("-pc_type ilu -ksp_max_it 30", "fluid");
541
542
#endif
  }
543
  FluidProblem *problem = malloc(sizeof(FluidProblem));
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
544
545
546
  Mesh *mesh = mesh_load_msh(mesh_file_name);
  if (!mesh)
    return NULL;
547
  problem->mesh_tree = mesh_tree_create(mesh);
548
  problem->rho = rho;
Matthieu Constant's avatar
Matthieu Constant committed
549

550
  problem->mu = mu;
551
  problem->g = g;
552
  problem->mesh = mesh;
Matthieu Constant's avatar
Matthieu Constant committed
553
  problem->epsilon = epsilon;
554
555
556
557
  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
558
    problem->strong_boundaries[i].apply = strong_boundaries[i].apply;
559
560
    problem->strong_boundaries[i].field = strong_boundaries[i].field;
  }
561
  problem->porosity = malloc(mesh->n_nodes*sizeof(double));
562
  problem->old_porosity = malloc(mesh->n_nodes*sizeof(double));
563
  problem->solution = malloc(mesh->n_nodes*n_fields*sizeof(double));
Matthieu Constant's avatar
Matthieu Constant committed
564
565
566
567
568
  for (int i = 0; i < problem->mesh->n_nodes; ++i)
  {
    problem->porosity[i] = 1.;
    problem->old_porosity[i] = 1.;
  }
569
  // begin to remove
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
570
  for (int i = 0; i < problem->mesh->n_nodes*n_fields; ++i){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
571
    problem->solution[i] = 0.;
572
  }
Matthieu Constant's avatar
Matthieu Constant committed
573
574
  problem->node_volume = malloc(0);
  fluid_problem_compute_node_volume(problem);
575
  // end to remove
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
576
  problem->linear_system = NULL;
577
#ifdef HXT_HAVE_PETSC
Matthieu Constant's avatar
Matthieu Constant committed
578
    hxtLinearSystemCreatePETSc(&problem->linear_system,mesh->n_elements,N_N,n_fields,mesh->elements,"fluid");
579
#else
580
    hxtLinearSystemCreateLU(&problem->linear_system,mesh->n_elements,N_N,n_fields,mesh->elements);
581
#endif
582
583
584
585
586
  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);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
587
  problem->particle_force = malloc(0);
588
589
  problem->particle_volume = malloc(0);
  problem->particle_velocity = malloc(0);
590
591
592
593
594
595
596
  return problem;
}

void fluid_problem_free(FluidProblem *problem)
{
  free(problem->solution);
  free(problem->porosity);
597
  free(problem->old_porosity);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
598
  hxtLinearSystemDelete(&problem->linear_system);
599
600
601
602
  free(problem->particle_uvw);
  free(problem->particle_element_id);
  free(problem->particle_position);
  free(problem->particle_mass);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
603
  free(problem->particle_force);
604
605
606
607
608
  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);
609
  mesh_tree_free(problem->mesh_tree);
610
611
612
  mesh_free(problem->mesh);
}

Matthieu Constant's avatar
Matthieu Constant committed
613
void fluid_problem_adapt_gen_mesh(FluidProblem *problem, double lcmax, double lcmin, double n_el)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
614
{
Matthieu Constant's avatar
Matthieu Constant committed
615
  double ngradM = 0.;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
616
617
618
  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
619
620
621
622
623
  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
624
625
626
627
628
629
630
  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
631
  double *porosity = problem->porosity;
Matthieu Constant's avatar
Matthieu Constant committed
632
633
  for (int iel = 0; iel < mesh->n_elements; ++iel) {
    const unsigned int *el = &mesh->elements[iel*N_N];
634
    double C[N_N],U[D][N_N], P[N_N];
Matthieu Constant's avatar
Matthieu Constant committed
635
636
    for (int i = 0; i < N_N; ++i){
      C[i] = porosity[el[i]];
637
638
      P[i] = solution[el[i]*n_fields+D];
      for (int j=0; j<D; ++j) {
Matthieu Constant's avatar
Matthieu Constant committed
639
640
641
        U[j][i] = solution[el[i]*n_fields+j];
      }
    }
642
643
644
    double dxdxi[D][D],dxidx[D][D];
    for (int i=0; i<D; ++i)
      for (int j=0; j<D; ++j)
Matthieu Constant's avatar
Matthieu Constant committed
645
        dxdxi[i][j] = mesh->x[el[j+1]*3+i]-mesh->x[el[0]*3+i];
646
#if D == 2
Matthieu Constant's avatar
Matthieu Constant committed
647
648
649
650
651
    const double vol = detDD(dxdxi)/2;
#else        
    const double vol = detDD(dxdxi)/6;
#endif    
    invDD(dxdxi, dxidx);
652
    double dphi[N_SF][D];
Matthieu Constant's avatar
Matthieu Constant committed
653
654
    grad_shape_functions(dxidx, dphi);
    double ngrad2 = 0;
655
656
    for (int i=0; i<D; ++i){
      for (int j=0; j<D; ++j){
Matthieu Constant's avatar
Matthieu Constant committed
657
658
659
660
661
662
663
664
665
666
667
668
        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;
669
    for (int j=0; j<D; ++j){
Matthieu Constant's avatar
Matthieu Constant committed
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
      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
688
  for (int iel = 0; iel < mesh->n_elements; ++iel) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
689
    const unsigned int *el = &mesh->elements[iel*N_N];
690
    double C[N_N],U[D][N_N], P[N_N];
691
692
    for (int i = 0; i < N_N; ++i){
      C[i] = porosity[el[i]];
693
694
      P[i] = solution[el[i]*n_fields+D];
      for (int j=0; j<D; ++j) {
695
696
697
        U[j][i] = solution[el[i]*n_fields+j];
      }
    }
698
699
700
    double dxdxi[D][D],dxidx[D][D];
    for (int i=0; i<D; ++i)
      for (int j=0; j<D; ++j)
701
702
        dxdxi[i][j] = mesh->x[el[j+1]*3+i]-mesh->x[el[0]*3+i];
    invDD(dxdxi, dxidx);
703
    double dphi[N_SF][D];
704
705
    grad_shape_functions(dxidx, dphi);
    double ngrad2 = 0;
706
707
    for (int i=0; i<D; ++i){
      for (int j=0; j<D; ++j){
708
709
710
711
712
713
714
        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
715
/*    double ngrad = pow(ngrad2, 0.25); */
716
/*    for (int j=0; j<D; ++j){*/
Matthieu Constant's avatar
Matthieu Constant committed
717
718
719
720
721
722
/*      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
723
    double ngrad = sqrt(ngrad2)/gradM;
Matthieu Constant's avatar
Matthieu Constant committed
724
    double lc = lcmin /fmin(1,fmax(ngrad/(gradmax),gradmin/(gradmax)));
Matthieu Constant's avatar
Matthieu Constant committed
725
    
Matthieu Constant's avatar
Matthieu Constant committed
726
    ngrad2 = 0.; 
727
    for (int j=0; j<D; ++j){
Matthieu Constant's avatar
Matthieu Constant committed
728
729
730
731
732
733
      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
734
735
    ngrad = sqrt(ngrad2)/gradPM;
    lc = fmin(fmin(lcmin /fmin(1,fmax(ngrad/(gradPmax),gradPmin/(gradPmax))),lc),lcmax);
Matthieu Constant's avatar
Matthieu Constant committed
736
737
    
    
Matthieu Constant's avatar
Matthieu Constant committed
738
739
740
741
742
    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
743
  double nOfEl = 0.0;
744
#if D==2
Matthieu Constant's avatar
Matthieu Constant committed
745
746
  for (int i = 0; i < mesh->n_nodes; ++i){
    new_mesh_size[i] /= num_lc[i];
Matthieu Constant's avatar
Matthieu Constant committed
747
748
749
750
751
752
753
    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
754
#else
Matthieu Constant's avatar
Matthieu Constant committed
755
  for (int i = 0; i < mesh->n_nodes; ++i){
Matthieu Constant's avatar
Matthieu Constant committed
756
757
    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
758
    nOfEl += nOfElByNode;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
759
  }
Matthieu Constant's avatar
Matthieu Constant committed
760
  double ratioEl = n_el/nOfEl;
Matthieu Constant's avatar
Matthieu Constant committed
761
  printf("RatioEl = %e\n",ratioEl);
Matthieu Constant's avatar
Matthieu Constant committed
762
763
764
765
  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
766
  FILE *f = fopen("adapt/lc.pos", "w");
767
768
  if(!f)
    printf("cannot open file adapt/lc.pos for writing\n");
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
769
  fprintf(f, "View \"lc\" {\n");
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
770
  for (int iel = 0; iel < mesh->n_elements; ++iel) {
771
    unsigned int *el = problem->mesh->elements+iel*N_N;
772
    fprintf(f, D == 2 ? "ST(" : "SS(");
773
774
775
776
    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
777
  }
778
  fprintf(f,"};\n");
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
779
780
  fclose(f);
  free(new_mesh_size);
Matthieu Constant's avatar
Matthieu Constant committed
781
  free(num_lc);
782
  system("gmsh -"xstr(D)" adapt/mesh.geo");
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
783
784
785
}


Matthieu Constant's avatar
Matthieu Constant committed
786
void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin, double n_el)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
787
{
Matthieu Constant's avatar
Matthieu Constant committed
788
789
790
  struct timespec tic, toc;
  fluid_problem_adapt_gen_mesh(problem, lcmax, lcmin, n_el);
  clock_get(&tic);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
791
  Mesh *new_mesh = mesh_load_msh("adapt/mesh.msh");
Matthieu Constant's avatar
Matthieu Constant committed
792
793
  clock_get(&toc);
  printf("Time mesh_load : %.6e \n",clock_diff(&tic,&toc));
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
794
  double *new_solution = malloc(sizeof(double)*new_mesh->n_nodes*n_fields);
795
  double *new_xi = malloc(sizeof(double)*new_mesh->n_nodes*D);
Matthieu Constant's avatar
Matthieu Constant committed
796
  clock_get(&tic);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
797
  int *new_eid = malloc(sizeof(int)*new_mesh->n_nodes);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
798
799
  for (int i = 0; i < new_mesh->n_nodes; ++i)
    new_eid[i] = -1;
Matthieu Constant's avatar
Matthieu Constant committed
800
801
802
  clock_get(&toc);
  printf("Time new_eid : %.6e \n",clock_diff(&tic,&toc));
  clock_get(&tic);
803
  double *newx = malloc(sizeof(double)*D*new_mesh->n_nodes);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
804
  for (int i = 0;i < new_mesh->n_nodes;++i) 
805
806
    for (int k = 0; k < D; ++k)
      newx[i*D+k] = new_mesh->x[i*3+k];
Matthieu Constant's avatar
Matthieu Constant committed
807
808
809
  clock_get(&toc);
  printf("Time newx : %.6e \n",clock_diff(&tic,&toc));
  clock_get(&tic);
810
  mesh_tree_particle_to_mesh(problem->mesh_tree,new_mesh->n_nodes, newx, new_eid, new_xi);
Matthieu Constant's avatar
Matthieu Constant committed
811
812
  clock_get(&toc);
  printf("Time find new point in old mesh : %.6e \n",clock_diff(&tic,&toc));
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
813
  for (int i = 0; i < new_mesh->n_nodes; ++i) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
814
815
816
817
818
    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);
    }
819
    double phi[N_SF];
820
    shape_functions(new_xi+i*D, phi);
821
822
823
    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
824
        new_solution[i*n_fields+j] += problem->solution[el[k]*n_fields+j]*phi[k];
825
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
826
827
828
829
830
831
832
833
834
835
  }
  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
836
837
838
839
840
  for (int i = 0; i < new_mesh->n_nodes; ++i)
  {
    problem->porosity[i] = 1.;
    problem->old_porosity[i] = 1.;
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
841
842
  mesh_free(problem->mesh);
  problem->mesh = new_mesh;
843
844
  mesh_tree_free(problem->mesh_tree);
  problem->mesh_tree = mesh_tree_create(problem->mesh);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
845
846
  for (int i = 0; i < problem->n_particles; ++i)
    problem->particle_element_id[i] = -1;
847
  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
848
  fluid_problem_compute_node_volume(problem);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
849
  fluid_problem_compute_porosity(problem);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
850
851
  hxtLinearSystemDelete(&problem->linear_system);
  #ifdef HXT_HAVE_PETSC
Matthieu Constant's avatar
Matthieu Constant committed
852
    hxtLinearSystemCreatePETSc(&problem->linear_system,new_mesh->n_elements,N_N,n_fields,new_mesh->elements,"fluid");
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
853
  #else
854
    hxtLinearSystemCreateLU(&problem->linear_system,new_mesh->n_elements,N_N,n_fields,new_mesh->elements);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
855
  #endif 
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
856
857
}

Matthieu Constant's avatar
Matthieu Constant committed
858
859
860
861
void fluid_problem_after_import(FluidProblem *problem)
{
  mesh_tree_free(problem->mesh_tree);
  problem->mesh_tree = mesh_tree_create(problem->mesh);
862
863
  free(problem->old_porosity);
  problem->old_porosity = malloc(sizeof(double)*problem->mesh->n_nodes);
Matthieu Constant's avatar
Matthieu Constant committed
864
865
866
867
868
869
870
871
872
873
874
875
876
  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
877

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
878
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity, long int *elid, int init)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
879
{
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
880
881
882
  
  struct timespec tic, toc;
  
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
883
884
885
886
887
888
  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
889
    free(problem->particle_force);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
890
891
    free(problem->particle_volume);
    free(problem->particle_velocity);
892
    problem->particle_uvw = malloc(sizeof(double)*D*n);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
893
    problem->particle_element_id = malloc(sizeof(int)*n);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
894
895
896
    for (int i = 0; i < n; ++i) {
      problem->particle_element_id[i]=-1;
    }
897
    problem->particle_position = malloc(sizeof(double)*D*n);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
898
    problem->particle_mass = malloc(sizeof(double)*n);
899
    problem->particle_force = malloc(sizeof(double)*n*D);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
900
    problem->particle_volume = malloc(sizeof(double)*n);
901
    problem->particle_velocity = malloc(sizeof(double)*n*D);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
902
903
904
905