fluid_problem.c 34.8 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
// 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);
    }
Matthieu Constant's avatar
Matthieu Constant committed
244
    f[j] = rho*((u[j]-uold[j])/dt + u[j]/q*divu + utau/q);
245
246
247
  }
}

248
249
250
251
252
253
static void f_p_0(double *f, const double du[D][D], const double q, const double qold, const double dt)
{
  double divu = 0;
  for (int j = 0; j < D; ++j) {
    divu += du[j][j];
  }
Matthieu Constant's avatar
Matthieu Constant committed
254
  *f = divu+(q-qold)/dt;
255
256
}

Matthieu Constant's avatar
Matthieu Constant committed
257
258
259
260
261
static void f_q_0(double *f, const double c, const double q)
{
  *f = c-q;
}

262
263
264
265
266
267

//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) {
268
      tau[i][j] = du[i][j]-u[i]*dq[j]/q;
269
270
271
272
    }
  }
  for (int j = 0; j < D; ++j) {
    for (int k = 0; k < D; ++k) {
Matthieu Constant's avatar
Matthieu Constant committed
273
      f[j][k] = mu*(tau[k][j]+tau[j][k]) + (k==j ? -p : 0);
274
275
276
277
    }
  }
}

278
279
280
281
282
283
284
285
286
static void f_p_1(double f[D], const double dp[D], const double epsilon)
{
  for (int j = 0; j < D; ++j) {
    f[j] = epsilon*dp[j];
  }
}



Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
287
288
289
290
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;
291
  //a generaliser
Matthieu Constant's avatar
Matthieu Constant committed
292
  double deps = 1e-8;
293
  
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
294
295
296
297
298
  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
299
  const double epsilon = problem->epsilon;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
300
301
302
303
304
305
306
307
  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;
308
  //compute_node_force_implicit(problem, dt, all_local_vector, all_local_matrix);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
309
  for (int iel=0; iel < mesh->n_elements; ++iel) {
310
    const unsigned int *el = &mesh->elements[iel*N_N];
311
312
313
    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)
314
315
316
        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
317
318
    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];
Matthieu Constant's avatar
Matthieu Constant committed
319
320
321
322
    double dq[D]={0}, du[D][D]={{0}}, dp[D]={0};

    int P = D;
    int Q = D+1;
323
    for (int i = 0; i< N_SF; ++i) {
324
      for (int j = 0; j < D; ++j) {
Matthieu Constant's avatar
Matthieu Constant committed
325
326
        dq[j] += dphi[i][j]*solution[el[i]*n_fields+Q];
        dp[j] += dphi[i][j]*solution[el[i]*n_fields+P];
327
        for (int k = 0; k < D; ++k)
328
329
330
          du[k][j] += dphi[i][j]*solution[el[i]*n_fields+k];
      }
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
331
    for (int iqp = 0; iqp< N_QP; ++iqp) {
332
333
      double phi[N_SF];
      shape_functions(QP[iqp],phi);
334
      double u[D]={0}, uold[D]={0}, p=0, c=0, qold=0, q=0;
335
      
336
      for (int i = 0; i < N_SF; ++i) {
337
        for (int j = 0; j < D; ++j){
338
          u[j] += solution[el[i]*n_fields+j]*phi[i];
Matthieu Constant's avatar
Matthieu Constant committed
339
          uold[j] += solution_old[el[i]*n_fields+j]*phi[i];
340
        }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
341
342
        p += solution[el[i]*n_fields+P]*phi[i];
        q +=  solution[el[i]*n_fields+Q]*phi[i];
343
        qold += solution_old[el[i]*n_fields+Q]*phi[i];
344
345
        c += porosity[el[i]]*phi[i];
      }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
346
      const double jw = QW[iqp]*detj;
347

Matthieu Constant's avatar
Matthieu Constant committed
348
      double fu0[D],fu0e[D],fu1[D][D],fu1e[D][D],fp0,fp0e,fp1[D],fp1e[D],fq0,fq0e;
349
350
      f_u_0(fu0,u,du,q,dq,uold,dt,rho);
      f_u_1(fu1,u,du,q,dq,p,mu);
351
352
      f_p_0(&fp0,du,q,qold,dt);
      f_p_1(fp1,dp,epsilon);
Matthieu Constant's avatar
Matthieu Constant committed
353
      f_q_0(&fq0,c,q);
354

Matthieu Constant's avatar
Matthieu Constant committed
355
      double fu0_u[D][D],fu1_u[D][D][D],fu0_du[D][D][D],fu1_du[D][D][D][D],fu1_p[D][D],fu0_q[D],fu1_q[D][D],fu0_dq[D][D],fu1_dq[D][D][D];
Matthieu Constant's avatar
Matthieu Constant committed
356
      
Matthieu Constant's avatar
Matthieu Constant committed
357
      double fp0_du[D][D],fp1_dp[D][D],fp0_q,fq0_q;
Matthieu Constant's avatar
Matthieu Constant committed
358
      
Matthieu Constant's avatar
Matthieu Constant committed
359
360

      //discrete derivatives
Matthieu Constant's avatar
Matthieu Constant committed
361
362
363
364
365
      //k:equation; i:dphii; l:variable; j:dphij
      double buf;
      for (int l = 0; l < D; ++l){
        buf = u[l];
        u[l] += deps;
366
367
        f_u_0(fu0e,u,du,q,dq,uold,dt,rho);
        f_u_1(fu1e,u,du,q,dq,p,mu);
Matthieu Constant's avatar
Matthieu Constant committed
368
369
370
371
372
        u[l] = buf;
        for (int k = 0; k < D; ++k){
          fu0_u[k][l] = (fu0e[k]-fu0[k])/deps;
          for (int i = 0; i < D; ++i){
            fu1_u[k][i][l] = (fu1e[k][i]-fu1[k][i])/deps;
373
374
          }
        }
375
        for (int j = 0; j < D; ++j){
Matthieu Constant's avatar
Matthieu Constant committed
376
377
378
          buf = du[l][j];
          du[l][j] += deps;
          f_p_0(&fp0e,du,q,qold,dt);
379
380
          f_u_0(fu0e,u,du,q,dq,uold,dt,rho);
          f_u_1(fu1e,u,du,q,dq,p,mu);
Matthieu Constant's avatar
Matthieu Constant committed
381
382
          du[l][j] = buf;
          fp0_du[l][j] = (fp0e-fp0)/deps;
383
          for (int k = 0; k < D; ++k){
Matthieu Constant's avatar
Matthieu Constant committed
384
385
386
            fu0_du[k][l][j] = (fu0e[k]-fu0[k])/deps;
            for (int i = 0; i< D; ++i){
              fu1_du[k][i][l][j] = (fu1e[k][i]-fu1[k][i])/deps;
387
388
389
390
            }
          }
        }
      }
Matthieu Constant's avatar
Matthieu Constant committed
391
392
393
394
395
396
397
      buf = p;
      p += deps;
      f_u_1(fu1e,u,du,q,dq,p,mu);
      p = buf;
      for (int k = 0; k < D; ++k){
        for (int i = 0; i < D; ++i){
          fu1_p[k][i] = (fu1e[k][i]-fu1[k][i])/deps;
398
399
        }
      }
Matthieu Constant's avatar
Matthieu Constant committed
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
      for (int j = 0; j < D; ++j){
        buf = dp[j];
        dp[j] += deps;
        f_p_1(fp1e,dp,epsilon);
        dp[j] = buf;
        for (int i = 0; i < D; ++i){
          fp1_dp[i][j] = (fp1e[i]-fp1[i])/deps;
        }
      }
      buf = q;
      q += deps;
      f_q_0(&fq0e,c,q);
      f_u_0(fu0e,u,du,q,dq,uold,dt,rho);
      f_u_1(fu1e,u,du,q,dq,p,mu);
      f_p_0(&fp0e,du,q,qold,dt);
      q = buf;
      fq0_q = (fq0e-fq0)/deps;
      fp0_q = (fp0e-fp0)/deps;
      for (int k = 0; k < D;++k){
        fu0_q[k] = (fu0e[k]-fu0[k])/deps;
        for (int i = 0; i < D; ++i){
          fu1_q[k][i] = (fu1e[k][i]-fu1[k][i])/deps;
        }
      }
      for (int j = 0; j < D; ++j){
        buf = dq[j];
        dq[j] += deps;
        f_u_0(fu0e,u,du,q,dq,uold,dt,rho);
        f_u_1(fu1e,u,du,q,dq,p,mu);
        dq[j] = buf;
        for (int k = 0; k < D; ++k){
          fu0_dq[k][j] = (fu0e[k]-fu0[k])/deps;
          for (int i = 0; i < D; ++i){
            fu1_dq[k][i][j] = (fu1e[k][i]-fu1[k][j])/deps;
          }
435
436
437
        }
      }

Matthieu Constant's avatar
Matthieu Constant committed
438
      //filling
439
      for (int iphi = 0; iphi < N_SF; ++iphi){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
440
        double phii = phi[iphi];
441
        const double *dphii = dphi[iphi];
Matthieu Constant's avatar
Matthieu Constant committed
442
443
444
445
        for (int k = 0; k < D; ++k){
          double f = fu0[k]*phii;
          for (int i = 0; i < D; ++i){
            f += fu1[k][i]*dphii[i];
446
          }
Matthieu Constant's avatar
Matthieu Constant committed
447
          local_vector[iphi+N_SF*k] += jw*f;
448
        }
449
        double ff = fp0*phii;
Matthieu Constant's avatar
Matthieu Constant committed
450
451
        for (int i = 0; i < D; ++i){
          ff += fp1[i]*dphii[i];
452
        }
453
        local_vector[iphi+N_SF*P] += jw*ff;
Matthieu Constant's avatar
Matthieu Constant committed
454
        local_vector[iphi+N_SF*Q] += jw*fq0*phii;
Matthieu Constant's avatar
Matthieu Constant committed
455
        
456
        for (int jphi = 0; jphi < N_SF; ++jphi){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
457
          double phij = phi[jphi];
458
          const double *dphij = dphi[jphi];
459

460
461
          int N2 = n_fields*N_SF;
          int IDX = (iphi*N2+jphi)*n_fields;
462
#define LOCAL_MATRIX(a,b) local_matrix[IDX+N2*a+b]
Matthieu Constant's avatar
Matthieu Constant committed
463
464
          for (int k = 0; k < D; ++k){
            int U = k;
Matthieu Constant's avatar
Matthieu Constant committed
465
            double mu_p = 0;
Matthieu Constant's avatar
Matthieu Constant committed
466
467
468
469
            double mu_q = fu0_q[k]*phii*phij;
            for (int l = 0; l < D; ++l){
              int V = l;
              double mu_u = fu0_u[k][l]*phii*phij;
470
              for (int i = 0; i < D; ++i){
Matthieu Constant's avatar
Matthieu Constant committed
471
472
473
                mu_u += fu1_u[k][i][l]*dphii[i]*phij;
                for (int j = 0; j < D; ++j){
                  mu_u += fu1_du[k][i][l][j]*dphii[i]*dphij[j];
474
475
                }
              }
Matthieu Constant's avatar
Matthieu Constant committed
476
477
478
              for (int j = 0; j < D; ++j){
                mu_u += fu0_du[k][l][j]*phii*dphij[j];
              }
479
              LOCAL_MATRIX(U,V) += jw*mu_u;
480
            }
Matthieu Constant's avatar
Matthieu Constant committed
481
482
483
484
485
486
487
488
489
490
            for (int i = 0; i < D; ++i){
              for (int j = 0; j < D; ++j){
                mu_q += fu1_dq[k][i][j]*dphii[i]*dphij[j];
              }
              mu_p += fu1_p[k][i]*phij*dphii[i];
              mu_q += fu1_q[k][i]*phij*dphii[i];
            }
            for (int j = 0; j < D; ++j){
              mu_q += fu0_dq[k][j]*phii*dphij[j];
            }
491
492
            LOCAL_MATRIX(U,P) += jw*mu_p;
            LOCAL_MATRIX(U,Q) += jw*mu_q;
493
          }
Matthieu Constant's avatar
Matthieu Constant committed
494
495
496
497
498
499
500
501
502
503
504
505
506
507
          for (int l = 0; l < D; ++l){
            int V = l;
            double mp_u = 0;
            for (int j = 0; j < D; ++j){
              mp_u += fp0_du[l][j]*dphij[j]*phii;
            }
            LOCAL_MATRIX(P,V) += jw*mp_u;
          }
          double mp_p = 0;
          for (int i = 0; i < D; ++i){
            for (int j = 0; j < D; ++j){
              mp_p += fp1_dp[i][j]*dphii[i]*dphij[j];
            }
          }
508
          LOCAL_MATRIX(P,P) += jw*mp_p;
Matthieu Constant's avatar
Matthieu Constant committed
509
          double mp_q = fp0_q*phii*phij;
510
          LOCAL_MATRIX(P,Q) += jw*mp_q;
Matthieu Constant's avatar
Matthieu Constant committed
511
          double mq_q = fq0_q*phij*phii;
512
          LOCAL_MATRIX(Q,Q) += jw*mq_q;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
513
514
515
516
517
        }
      }
    }
    hxtLinearSystemAddToMatrix(lsys,iel,iel,local_matrix);
    hxtLinearSystemAddToRhs(lsys,rhs,iel,local_vector);
518
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
519
520
  free(all_local_matrix);
  free(all_local_vector);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
521

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
522
523
524
525
526
527
528
529
530
531
  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
532
533
      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
534
535
    }
  }
536
}
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
537
538


539

540
static void apply_boundary_conditions(const Mesh *mesh, double *solution, int nBnd, StrongBoundary *bnds)
541
{
542
543
  for (int ibnd = 0; ibnd < nBnd; ++ibnd) {
    StrongBoundary *bnd = bnds + ibnd;
544
545
546
547
548
549
550
551
    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);
    }
552
    double *x = malloc(mesh->phys_n_nodes[iphys]*sizeof(double)*D);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
553
554
555
    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];
556
557
      for (int k = 0; k < D; ++k)
        x[i*D+k] = mesh->x[j*3+k];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
558
559
    }
    bnd->apply(mesh->phys_n_nodes[iphys], x, v);
560
    for (int i = 0; i < mesh->phys_n_nodes[iphys]; ++i)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
561
562
563
      solution[mesh->phys_nodes[iphys][i]*n_fields+bnd->field] = v[i];
    free(x);
    free(v);
564
565
566
  }
}

567
int fluid_problem_implicit_euler(FluidProblem *problem, double dt)
568
{
569
570
  static double totaltime = 0;
  struct timespec tic, toc;
571
  const Mesh *mesh = problem->mesh;
572
  apply_boundary_conditions(mesh, problem->solution, problem->n_strong_boundaries, problem->strong_boundaries);
573
  double *solution_old = malloc(mesh->n_nodes*n_fields*sizeof(double));
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
574
  double *rhs = malloc(mesh->n_nodes*n_fields*sizeof(double));
575
576
  for (int i=0; i<mesh->n_nodes*n_fields; ++i)
    solution_old[i] = problem->solution[i];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
577
  problem->solution_explicit = solution_old;
578
579
  double *dx = malloc(sizeof(double)*mesh->n_nodes*n_fields);
  double norm0 = 0;
580
581
  int i;
  for (i=0; i<newton_max_it; ++i) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
582
583
    double norm = 0;
    hxtLinearSystemZeroMatrix(problem->linear_system);
584
585

    for (int i=0; i<mesh->n_nodes*n_fields; ++i){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
586
      rhs[i] = 0;
587
588
589
      //printf("solution=%g\n",problem->solution[i]);
    }
    //break;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
590
591
    fluid_problem_assemble_system(problem, rhs, solution_old, dt);
    hxtLinearSystemGetRhsNorm(problem->linear_system,rhs,&norm);
592
593
594
595
596
597
598
    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;
599
  clock_get(&tic);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
600
    hxtLinearSystemSolve(problem->linear_system, rhs, dx);
601
  clock_get(&toc);
602
  //break;
603
  totaltime += clock_diff(&tic, &toc);
604
605
606
607
    for (int j=0; j<mesh->n_nodes*n_fields; ++j) {
      problem->solution[j] -= dx[j];
    }
  }
608
  printf("total solve time : %g\n", totaltime);
609
  free(dx);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
610
  free(rhs);
611
  free(solution_old);
612
  return i;
613
614
}

Matthieu Constant's avatar
Matthieu Constant committed
615
616
617
618
619
620
621
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
622
  for (int iel = 0; iel < mesh->n_elements; ++iel) {
623
    double dxdxi[D][D];
624
    const int *el = &mesh->elements[iel*N_N];
625
626
    for (int i=0; i<D; ++i)
      for (int j=0; j<D; ++j)
627
        dxdxi[i][j] = mesh->x[el[j+1]*3+i]-mesh->x[el[0]*3+i];
628
#if D == 2
629
    const double vol = detDD(dxdxi)/6;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
630
#else
631
    const double vol = detDD(dxdxi)/24;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
632
#endif
633
634
    for (int i = 0; i < N_N; ++i)
      problem->node_volume[el[i]] += vol;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
635
636
637
638
639
640
641
642
643
644
645
646
  }
}

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;
647
    double sf[N_SF];
648
    shape_functions(&problem->particle_uvw[i*D],sf);
649
    const uint32_t *el = &mesh->elements[problem->particle_element_id[i]*N_N];
650
651
    for (int j=0; j<N_N;++j)
      problem->porosity[el[j]] += sf[j]*volume[i];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
652
653
  }
  for (int i = 0; i < mesh->n_nodes; ++i){
Matthieu Constant's avatar
Matthieu Constant committed
654
655
656
657
    if(problem->node_volume[i]==0.){
      problem->porosity[i] = 1.;
    }
    else{
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
658
      problem->porosity[i] = (1-problem->porosity[i]/problem->node_volume[i]);
Matthieu Constant's avatar
Matthieu Constant committed
659
    }
Matthieu Constant's avatar
Matthieu Constant committed
660
661
662
  }
}

663
FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double mu, double rho, double epsilon, int n_strong_boundaries, StrongBoundary *strong_boundaries, int noGrains)
664
{
665
666
667
668
669
  static int initialized = 0;
  if (!initialized) {
    hxtInitializeLinearSystems(0, NULL);
    initialized = 1;
#ifdef HXT_HAVE_PETSC
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
670
    hxtPETScInsertOptions("-pc_type ilu -ksp_max_it 30", "fluid");
671
672
#endif
  }
673
  FluidProblem *problem = malloc(sizeof(FluidProblem));
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
674
675
676
  Mesh *mesh = mesh_load_msh(mesh_file_name);
  if (!mesh)
    return NULL;
677
  problem->mesh_tree = mesh_tree_create(mesh);
678
  problem->rho = rho;
Matthieu Constant's avatar
Matthieu Constant committed
679

680
  problem->mu = mu;
681
  problem->g = g;
682
  problem->mesh = mesh;
Matthieu Constant's avatar
Matthieu Constant committed
683
  problem->epsilon = epsilon;
684
685
686
687
  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
688
    problem->strong_boundaries[i].apply = strong_boundaries[i].apply;
689
690
    problem->strong_boundaries[i].field = strong_boundaries[i].field;
  }
691
  problem->porosity = malloc(mesh->n_nodes*sizeof(double));
692
  problem->old_porosity = malloc(mesh->n_nodes*sizeof(double));
693
694
695
696
697
  if(noGrains){
    for(int i = 0; i < mesh->n_nodes; ++i){
      problem->porosity[i] = 1.;
      problem->old_porosity[i] = 1.;
    }
Matthieu Constant's avatar
Matthieu Constant committed
698
  }
699
  problem->solution = malloc(mesh->n_nodes*n_fields*sizeof(double));
700
  // begin to remove
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
701
  for (int i = 0; i < problem->mesh->n_nodes*n_fields; ++i){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
702
    problem->solution[i] = 0.;
703
  }
704
705
706
  for (int i = 0; i < mesh->n_nodes; ++i){
    problem->solution[i*n_fields+(D+1)] = 1.;
  }
Matthieu Constant's avatar
Matthieu Constant committed
707
708
  problem->node_volume = malloc(0);
  fluid_problem_compute_node_volume(problem);
709
  // end to remove
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
710
  problem->linear_system = NULL;
711
#ifdef HXT_HAVE_PETSC
Matthieu Constant's avatar
Matthieu Constant committed
712
    hxtLinearSystemCreatePETSc(&problem->linear_system,mesh->n_elements,N_N,n_fields,mesh->elements,"fluid");
713
#else
714
    hxtLinearSystemCreateLU(&problem->linear_system,mesh->n_elements,N_N,n_fields,mesh->elements);
715
#endif
716
717
718
719
720
  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
721
  problem->particle_force = malloc(0);
722
723
  problem->particle_volume = malloc(0);
  problem->particle_velocity = malloc(0);
724
725
726
727
728
729
730
  return problem;
}

void fluid_problem_free(FluidProblem *problem)
{
  free(problem->solution);
  free(problem->porosity);
731
  free(problem->old_porosity);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
732
  hxtLinearSystemDelete(&problem->linear_system);
733
734
735
736
  free(problem->particle_uvw);
  free(problem->particle_element_id);
  free(problem->particle_position);
  free(problem->particle_mass);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
737
  free(problem->particle_force);
738
739
740
741
742
  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);
743
  mesh_tree_free(problem->mesh_tree);
744
745
746
  mesh_free(problem->mesh);
}

Matthieu Constant's avatar
Matthieu Constant committed
747
void fluid_problem_adapt_gen_mesh(FluidProblem *problem, double lcmax, double lcmin, double n_el)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
748
{
Matthieu Constant's avatar
Matthieu Constant committed
749
  double ngradM = 0.;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
750
751
752
  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
753
754
755
756
757
  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
758
759
760
761
762
763
764
  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
765
  double *porosity = problem->porosity;
Matthieu Constant's avatar
Matthieu Constant committed
766
767
  for (int iel = 0; iel < mesh->n_elements; ++iel) {
    const unsigned int *el = &mesh->elements[iel*N_N];
768
    double C[N_N],U[D][N_N], P[N_N];
Matthieu Constant's avatar
Matthieu Constant committed
769
770
    for (int i = 0; i < N_N; ++i){
      C[i] = porosity[el[i]];
771
772
      P[i] = solution[el[i]*n_fields+D];
      for (int j=0; j<D; ++j) {
Matthieu Constant's avatar
Matthieu Constant committed
773
774
775
        U[j][i] = solution[el[i]*n_fields+j];
      }
    }
776
777
778
    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
779
        dxdxi[i][j] = mesh->x[el[j+1]*3+i]-mesh->x[el[0]*3+i];
780
#if D == 2
Matthieu Constant's avatar
Matthieu Constant committed
781
782
783
784
785
    const double vol = detDD(dxdxi)/2;
#else        
    const double vol = detDD(dxdxi)/6;
#endif    
    invDD(dxdxi, dxidx);
786
    double dphi[N_SF][D];
Matthieu Constant's avatar
Matthieu Constant committed
787
788
    grad_shape_functions(dxidx, dphi);
    double ngrad2 = 0;
789
790
    for (int i=0; i<D; ++i){
      for (int j=0; j<D; ++j){
Matthieu Constant's avatar
Matthieu Constant committed
791
792
793
794
795
796
797
798
799
800
801
802
        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;
803
    for (int j=0; j<D; ++j){
Matthieu Constant's avatar
Matthieu Constant committed
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
      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
822
  for (int iel = 0; iel < mesh->n_elements; ++iel) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
823
    const unsigned int *el = &mesh->elements[iel*N_N];
824
    double C[N_N],U[D][N_N], P[N_N];
825
826
    for (int i = 0; i < N_N; ++i){
      C[i] = porosity[el[i]];
827
828
      P[i] = solution[el[i]*n_fields+D];
      for (int j=0; j<D; ++j) {
829
830
831
        U[j][i] = solution[el[i]*n_fields+j];
      }
    }
832
833
834
    double dxdxi[D][D],dxidx[D][D];
    for (int i=0; i<D; ++i)
      for (int j=0; j<D; ++j)
835
836
        dxdxi[i][j] = mesh->x[el[j+1]*3+i]-mesh->x[el[0]*3+i];
    invDD(dxdxi, dxidx);
837
    double dphi[N_SF][D];
838
839
    grad_shape_functions(dxidx, dphi);
    double ngrad2 = 0;
840
841
    for (int i=0; i<D; ++i){
      for (int j=0; j<D; ++j){
842
843
844
845
846
847
848
        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
849
/*    double ngrad = pow(ngrad2, 0.25); */
850
/*    for (int j=0; j<D; ++j){*/
Matthieu Constant's avatar
Matthieu Constant committed
851
852
853
854
855
856
/*      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
857
    double ngrad = sqrt(ngrad2)/gradM;
Matthieu Constant's avatar
Matthieu Constant committed
858
    double lc = lcmin /fmin(1,fmax(ngrad/(gradmax),gradmin/(gradmax)));
Matthieu Constant's avatar
Matthieu Constant committed
859
    
Matthieu Constant's avatar
Matthieu Constant committed
860
    ngrad2 = 0.; 
861
    for (int j=0; j<D; ++j){
Matthieu Constant's avatar
Matthieu Constant committed
862
863
864
865
866
867
      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
868
869
    ngrad = sqrt(ngrad2)/gradPM;
    lc = fmin(fmin(lcmin /fmin(1,fmax(ngrad/(gradPmax),gradPmin/(gradPmax))),lc),lcmax);
Matthieu Constant's avatar
Matthieu Constant committed
870
871
    
    
Matthieu Constant's avatar
Matthieu Constant committed
872
873
874
875
876
    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
877
  double nOfEl = 0.0;
878
#if D==2
Matthieu Constant's avatar
Matthieu Constant committed
879
880
  for (int i = 0; i < mesh->n_nodes; ++i){
    new_mesh_size[i] /= num_lc[i];
Matthieu Constant's avatar
Matthieu Constant committed
881
882
883
884
885
886
887
    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
888
#else
Matthieu Constant's avatar
Matthieu Constant committed
889
  for (int i = 0; i < mesh->n_nodes; ++i){
Matthieu Constant's avatar
Matthieu Constant committed
890
891
    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
892
    nOfEl += nOfElByNode;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
893
  }
Matthieu Constant's avatar
Matthieu Constant committed
894
  double ratioEl = n_el/nOfEl;
Matthieu Constant's avatar
Matthieu Constant committed
895
  printf("RatioEl = %e\n",ratioEl);
Matthieu Constant's avatar
Matthieu Constant committed
896
897
898
899
  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
900
  FILE *f = fopen("adapt/lc.pos", "w");
901
902
  if(!f)
    printf("cannot open file adapt/lc.pos for writing\n");
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
903
  fprintf(f, "View \"lc\" {\n");