fluid_problem.c 33.7 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


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

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
38
#if DIMENSION==2
39
40
41

#define N_SF 3
#define N_N 3
42
43
44
#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};
45
46
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
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
72
#else
73
74
75

#define N_SF 4
#define N_N 4
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
76
77
78
79
80
81
82
83
84
#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};
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
110
111
112
113
114
115
116
117
118
119
120
static void shape_functions(double *uvw, double *sf)
{
  sf[0] = 1-uvw[0]-uvw[1]-uvw[2];
  sf[1] = uvw[0];
  sf[2] = uvw[1];
  sf[3] = uvw[2];
}
static void grad_shape_functions(const double dxidx[3][3], double gsf[4][3])
{
  for (int i = 0; i < 3; ++i) {
    gsf[0][i] = -dxidx[0][i]-dxidx[1][i]-dxidx[2][i];
    gsf[1][i] = dxidx[0][i];
    gsf[2][i] = dxidx[1][i];
    gsf[3][i] = dxidx[2][i];
  }
};
static double detDD(const double m[3][3])
{
  return m[0][0]*(m[1][1]*m[2][2]-m[2][1]*m[1][2])
        -m[1][0]*(m[0][1]*m[2][2]-m[2][1]*m[0][2])
        +m[2][0]*(m[0][1]*m[1][2]-m[1][1]*m[0][2]);
}
static double invDD(const double m[3][3], double inv[3][3]) {
  double d = detDD(m);
  inv[0][0] = (m[1][1]*m[2][2]-m[2][1]*m[1][2])/d;
  inv[0][1] = -(m[0][1]*m[2][2]-m[2][1]*m[0][2])/d;
  inv[0][2] = (m[0][1]*m[1][2]-m[1][1]*m[0][2])/d;
  inv[1][0] = -(m[1][0]*m[2][2]-m[2][0]*m[1][2])/d;
  inv[1][1] = (m[0][0]*m[2][2]-m[2][0]*m[0][2])/d;
  inv[1][2] = -(m[0][0]*m[1][2]-m[1][0]*m[0][2])/d;
  inv[2][0] = (m[1][0]*m[2][1]-m[2][0]*m[1][1])/d;
  inv[2][1] = -(m[0][0]*m[2][1]-m[2][0]*m[0][1])/d;
  inv[2][2] = (m[0][0]*m[1][1]-m[1][0]*m[0][1])/d;
  return d;
};

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
121
122
#endif

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

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
150
void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt, double *particle_force) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
151
152
153
154
155
156
  for (int i = 0; i < problem->n_particles*DIMENSION; ++i) {
    particle_force[i] = problem->particle_force[i];
  }
}

static void compute_node_force_implicit(FluidProblem *problem, double dt, double *all_local_vector, double *all_local_matrix) {
Matthieu Constant's avatar
Matthieu Constant committed
157
  double *compacity = problem->compacity;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
158
  Mesh *mesh = problem->mesh;
159
  int n_fields = problem->nFields;
Matthieu Constant's avatar
debug    
Matthieu Constant committed
160
  int N_E = problem->mesh->n_elements;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
161
162
  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){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
165
166
167
      for (int d = 0; d < DIMENSION; ++d)
        problem->particle_force[i*DIMENSION+d] = 0;
      problem->particle_force[i*DIMENSION+1] = problem->g*mass;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
168
169
      continue;
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
170
171
172
173
174
    double *xi = problem->particle_uvw + DIMENSION*i;
    double phi[N_SF];
    shape_functions(xi,phi);
    uint32_t *el=problem->mesh->elements+iel*N_SF;
    double *se = problem->solution_explicit;
Matthieu Constant's avatar
Matthieu Constant committed
175
    
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
176
177
178
179
180
181
182
    double dxdxi[DIMENSION][DIMENSION],dxidx[DIMENSION][DIMENSION];
    for (int k=0; k<DIMENSION; ++k)
      for (int j=0; j<DIMENSION; ++j)
        dxdxi[k][j] = mesh->x[el[j+1]*3+k]-mesh->x[el[0]*3+k];
    invDD(dxdxi,dxidx);
    double dphi[N_SF][DIMENSION];
    grad_shape_functions(dxidx, dphi);
183
    double *si = problem->solution;
Matthieu Constant's avatar
debug    
Matthieu Constant committed
184
    for (int ff=0;ff<problem->nFluids;++ff){
Matthieu Constant's avatar
Matthieu Constant committed
185
186
    double c=0,vf[DIMENSION]={0}, vfe[DIMENSION]={0}, gradp[DIMENSION]={0};
    double beta = 0;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
187
    for (int iphi=0; iphi < N_SF; ++iphi) {
Matthieu Constant's avatar
Matthieu Constant committed
188
      c += compacity[el[iphi]]*phi[iphi];
Matthieu Constant's avatar
Matthieu Constant committed
189
      beta += si[el[iphi]*n_fields+DIMENSION+1+ff*(DIMENSION+1)];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
190
      for (int j=0; j < DIMENSION; ++j) {
Matthieu Constant's avatar
Matthieu Constant committed
191
192
        vf[j]  += si[el[iphi]*n_fields+j+ff*(DIMENSION+1)]*phi[iphi];
        vfe[j] += se[el[iphi]*n_fields+j+ff*(DIMENSION+1)]*phi[iphi];
Matthieu Constant's avatar
Matthieu Constant committed
193
        gradp[j] += dphi[iphi][j]*si[el[iphi]*n_fields+n_fields-2];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
194
195
196
197
      }
    }
    double du[DIMENSION], due[DIMENSION];
    for (int j = 0; j < DIMENSION; ++j) {
Matthieu Constant's avatar
Matthieu Constant committed
198
199
      du[j] = problem->particle_velocity[i*DIMENSION+j]-vf[j]/(1-c);
      due[j] = problem->particle_velocity[i*DIMENSION+j]-vfe[j]/(1-c);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
200
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
201
    double vol = problem->particle_volume[i];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
202
    double rhop = mass/vol;
Matthieu Constant's avatar
debug    
Matthieu Constant committed
203
    double drho = rhop -problem->rho[ff];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
204
    double g = problem->g*drho/rhop;
Matthieu Constant's avatar
Matthieu Constant committed
205
    double gamma = particle_drag_coeff(due,problem->mu[ff],problem->rho[ff],vol,1-c);
Matthieu Constant's avatar
debug    
Matthieu Constant committed
206
    double fcoeff = beta*mass/(mass+dt*gamma);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
207
208
209
210
211
    
    for (int d = 0; d < DIMENSION; ++d)
      problem->particle_force[i*DIMENSION+d] = fcoeff*(-gamma*du[d]-vol*gradp[d]);
    problem->particle_force[i*DIMENSION+1] += fcoeff*g*mass;
    
Matthieu Constant's avatar
Matthieu Constant committed
212

Matthieu Constant's avatar
Matthieu Constant committed
213
214
    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];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
215
216
    for (int iphi = 0; iphi < N_SF; ++iphi) {
      for (int d = 0; d < DIMENSION; ++d)
Matthieu Constant's avatar
Matthieu Constant committed
217
218
        local_vector[iphi+N_SF*d+ff*(DIMENSION+1)] -= fcoeff*gamma*(du[d]-dt/mass*vol*gradp[d])*phi[iphi];
      local_vector[iphi+N_SF*1+ff*(DIMENSION+1)] -= fcoeff*gamma*dt*g*phi[iphi];
219
      
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
220
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
221
222
    for (int iphi = 0; iphi < N_SF; ++iphi) {
      for (int jphi = 0; jphi < N_SF; ++jphi) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
223
224
225
        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
226
227
        double f = fcoeff*phi[iphi];
        for (int d = 0; d < DIMENSION; ++d){
Matthieu Constant's avatar
Matthieu Constant committed
228
229
          LOCAL_MATRIX(d+ff*(DIMENSION+1),d+ff*(DIMENSION+1)) -= -beta*f/(1-c)*phi[jphi]*gamma;
          LOCAL_MATRIX(d+ff*(DIMENSION+1),DIMENSION) -= -beta*f*gamma*dt/mass*vol*dphi[jphi][d];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
230
        }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
231
232
      }
    }
Matthieu Constant's avatar
Matthieu Constant committed
233
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
234
235
236
237
238
239
240
241
  }
}

static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, const double *solution_old, double dt)
{
  HXTLinearSystem *lsys = problem->linear_system;
  const Mesh *mesh = problem->mesh;
  const double *solution = problem->solution;
242
243
244
  int n_fields = problem->nFields;
  const double *compacity = problem->compacity;
  const double *old_compacity = problem->old_compacity;
Matthieu Constant's avatar
Matthieu Constant committed
245
246
  const double *mu = problem->mu;
  const double *rho = problem->rho;
Matthieu Constant's avatar
Matthieu Constant committed
247
  const double epsilon = problem->epsilon;
248
  int n_fluids = problem->nFluids;
Matthieu Constant's avatar
Matthieu Constant committed
249
  size_t local_size = N_SF*n_fields;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
250
251
252
253
254
255
256
257
  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);
Matthieu Constant's avatar
Matthieu Constant committed
258
 
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
259
  for (int iel=0; iel < mesh->n_elements; ++iel) {
260
261
    const unsigned int *el = &mesh->elements[iel*N_N];
    double dxdxi[DIMENSION][DIMENSION], dxidx[DIMENSION][DIMENSION], dphi[N_N][DIMENSION];
Matthieu Constant's avatar
Matthieu Constant committed
262
263
    int P = DIMENSION;    
    double *local_vector = &all_local_vector[N_SF*n_fields*iel];
Matthieu Constant's avatar
Matthieu Constant committed
264
265
    double c=0,p=0, dp[DIMENSION]={0};
    for (int i = 0; i < DIMENSION; ++i)
266
267
        for (int j = 0; j < DIMENSION; ++j)
          dxdxi[i][j] = mesh->x[el[j+1]*3+i] - mesh->x[el[0]*3+i];
Matthieu Constant's avatar
debug    
Matthieu Constant committed
268
269
    const double detj = invDD(dxdxi, dxidx);
    grad_shape_functions(dxidx, dphi);
Matthieu Constant's avatar
Matthieu Constant committed
270
    for (int iqp = 0; iqp< N_QP; ++iqp) {
Matthieu Constant's avatar
debug    
Matthieu Constant committed
271
272
      const double jwl = QW[iqp]*detj;
      double phi[N_SF];
Matthieu Constant's avatar
Matthieu Constant committed
273
274
275
276
277
278
279
280
281
      shape_functions(QP[iqp],phi);
      for (int i = 0; i < N_SF; ++i) {
          double phii = phi[i];
          p += solution[el[i]*n_fields+P]*phi[i];
          c += compacity[el[i]]*phi[i];
          local_vector[i+N_SF*P] += jwl*phii*(c-1);
          for (int j = 0; j < DIMENSION; ++j)
            dp[j]   += dphi[i][j]*solution[el[i]*n_fields+P];
      }
Matthieu Constant's avatar
Matthieu Constant committed
282
    }//printf("c=%g,p=%g\n",c,p);
Matthieu Constant's avatar
Matthieu Constant committed
283
    
Matthieu Constant's avatar
Matthieu Constant committed
284
    double *local_matrix = &all_local_matrix[N_SF*N_SF*n_fields*n_fields*iel];
Matthieu Constant's avatar
Matthieu Constant committed
285
286
    for (int ff = 0; ff < n_fluids; ++ff)
    {
Matthieu Constant's avatar
Matthieu Constant committed
287
      int Q = n_fields-1;//DIMENSION+ff*(DIMENSION+1);
Matthieu Constant's avatar
Matthieu Constant committed
288
      double dq[DIMENSION]={0}, du[DIMENSION][DIMENSION]={{0}}, dqdq[DIMENSION]={0};
289
290
      for (int i = 0; i< N_SF; ++i) {
        for (int j = 0; j < DIMENSION; ++j) {
Matthieu Constant's avatar
Matthieu Constant committed
291
          dq[j]   += dphi[i][j]*solution[el[i]*n_fields+Q];
Matthieu Constant's avatar
Matthieu Constant committed
292
          dqdq[j] += dphi[i][j];
293
294
          for (int k = 0; k < DIMENSION; ++k)
            du[k][j] += dphi[i][j]*solution[el[i]*n_fields+k+ff*(DIMENSION+1)];
295
296
        }
      }
297
298
299
      for (int iqp = 0; iqp< N_QP; ++iqp) {
        double phi[N_SF];
        shape_functions(QP[iqp],phi);
Matthieu Constant's avatar
Matthieu Constant committed
300
        double u[DIMENSION]={0}, dudt[DIMENSION]={0}, dqdt=0, q=0;
301
302
303
        for (int i = 0; i < N_SF; ++i) {
          for (int j = 0; j < DIMENSION; ++j){
            u[j] += solution[el[i]*n_fields+j+ff*(DIMENSION+1)]*phi[i];
Matthieu Constant's avatar
debug    
Matthieu Constant committed
304
            dudt[j] += (solution[el[i]*n_fields+j+ff*(DIMENSION+1)]-solution_old[el[i]*n_fields+j+ff*(DIMENSION+1)])/dt*phi[i];
305
          }
306
          q +=  solution[el[i]*n_fields+Q]*phi[i];
Matthieu Constant's avatar
Matthieu Constant committed
307
          dqdt += (solution[el[i]*n_fields+Q]-solution_old[el[i]*n_fields+Q])*phi[i]/dt;
308
          //dcdt += (problem->porosity[el[i]]-problem->old_porosity[el[i]])/dt*phi[i];
309
        }
Matthieu Constant's avatar
Matthieu Constant committed
310
          printf("q=%g,dqdt=%g\n",q,dqdt);
311
312
313
314
315
316
317
        const double jw = QW[iqp]*detj;
        for (int iphi = 0; iphi < N_SF; ++iphi) {
          double phii = phi[iphi];
          const double *dphii = dphi[iphi];
          double tau[DIMENSION][DIMENSION];
          for (int i = 0; i < DIMENSION; ++i)
            for (int j = 0; j < DIMENSION; ++j)
Matthieu Constant's avatar
debug    
Matthieu Constant committed
318
              tau[i][j] = du[i][j]-u[i]*dq[j]/q;
319
320
321
322
          double dphidp = 0, divu = 0;
          for (int k = 0; k < DIMENSION; ++k){
            dphidp += dphii[k]*dp[k];
            divu += du[k][k];
323
324
          }
          for (int j = 0; j < DIMENSION; ++j) {
325
326
            double utau = 0;
            double dphisigma=0;
327
            for (int k = 0; k < DIMENSION; ++k) {
328
              //utau += u[k]*dphii[k]*u[j]/c;
Matthieu Constant's avatar
debut q    
Matthieu Constant committed
329
              utau += u[k]*tau[j][k];
Matthieu Constant's avatar
Matthieu Constant committed
330
              dphisigma += 2*mu[ff]*dphii[k]*0.5*(tau[k][j]+tau[j][k]);
331
332
333
334
            }
            local_vector[iphi+N_SF*j] += jw*(
                dphisigma
                //+rho*phii*dudt[j]-rho*utau
Matthieu Constant's avatar
Matthieu Constant committed
335
                +rho[ff]*phii*(dudt[j]+u[j]/q*divu+utau/q)
336
337
338
                -p*dphii[j]
                );
          }
Matthieu Constant's avatar
Matthieu Constant committed
339
          local_vector[iphi+N_SF*P] += jw*phii*(q);
Matthieu Constant's avatar
Matthieu Constant committed
340
          local_vector[iphi+N_SF*Q] += jw*phii*(dqdt+divu)+jw*epsilon*dphidp;
341
342
343
344
345
346
347
348
349
350
351
352
353
          for (int jphi = 0; jphi < N_SF; ++jphi) {
            double phij = phi[jphi];
            const double *dphij = dphi[jphi];
            //int N2 = N_SF*N_SF;
            //int IDX = iphi*N2+jphi;
            int N2 = n_fields*N_SF;
            int IDX = (iphi*N2+jphi)*n_fields;
  #define LOCAL_MATRIX(a,b) local_matrix[IDX+N2*a+b]
            double dphiidphij = 0;
            double dtau[DIMENSION];
            double udtau = 0;
            double dphiidtau = 0;
            for (int j = 0; j < DIMENSION; ++j) {
Matthieu Constant's avatar
Matthieu Constant committed
354
              dtau[j] = dphij[j]-phij*dq[j]/q;
355
356
357
358
359
              dphiidphij += dphii[j]*dphij[j];
              udtau += u[j]*dtau[j];
              dphiidtau += dphii[j]*dtau[j];
            }
            for (int j = 0; j < DIMENSION; ++j) {
Matthieu Constant's avatar
Matthieu Constant committed
360
              int U = j+ff*(DIMENSION+1);
361
              LOCAL_MATRIX(U,P) += -jw*dphii[j]*phij;
Matthieu Constant's avatar
Matthieu Constant committed
362
              LOCAL_MATRIX(Q,U) +=  jw*phii*dphij[j];
363
364
365
366
              double utau =0;
              double dphisigmadq = 0;
              double dutaudq = 0;
              for (int k = 0; k < DIMENSION; ++k) {
Matthieu Constant's avatar
Matthieu Constant committed
367
                int V = k+ff*(DIMENSION+1);
368
                utau += u[k]*tau[j][k];
Matthieu Constant's avatar
Matthieu Constant committed
369
370
                dutaudq += u[k]*u[j]*(dq[k]*phij/(q*q)-dqdq[k]/q);
                dphisigmadq += 2*mu[ff]*dphii[k]*0.5*(u[k]*dq[j]*phij/(q*q)+u[j]*dq[k]*phij/(q*q)-u[k]*dqdq[j]/q-u[j]*dqdq[k]/q);
371
372
                //LOCAL_MATRIX(U,V) += -jw*rho*u[j]*dphii[k]*phij/c+jw*2*mu*dphii[k]*dtau[j]*0.5;
                LOCAL_MATRIX(U,V) +=
Matthieu Constant's avatar
Matthieu Constant committed
373
374
375
                   jw*rho[ff]*phii*phij*tau[j][k]/q
                  +jw*2*mu[ff]*dphii[k]*dtau[j]*0.5
                  +jw*rho[ff]*phii*u[j]*dphij[k]/q;
376
377
              }
              //LOCAL_MATRIX(U,U) += jw*mu*2*0.5*dphiidtau + jw*rho*(phii*phij/dt-phij*dphii[j]*u[j]/c);
Matthieu Constant's avatar
Matthieu Constant committed
378
379
              LOCAL_MATRIX(U,U) += jw*mu[ff]*2*0.5*dphiidtau + jw*rho[ff]*phii*(phij/dt+phij/q*divu+udtau/q);
              LOCAL_MATRIX(U,Q) += jw*rho[ff]*phii*(-u[j]*divu*phij/(q*q)-utau*phij/(q*q)+dutaudq/q)+jw*dphisigmadq;
380
            }
Matthieu Constant's avatar
Matthieu Constant committed
381
382
383
            LOCAL_MATRIX(Q,P) += jw*epsilon*dphiidphij;
            LOCAL_MATRIX(P,Q) += -jw*phii*phij;
            LOCAL_MATRIX(Q,Q) += jw*phii*phij/dt;
384
          }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
385
386
        }
      }
387
388
      hxtLinearSystemAddToMatrix(lsys,iel,iel,local_matrix);
      hxtLinearSystemAddToRhs(lsys,rhs,iel,local_vector);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
389
    }
390
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
391
392
  free(all_local_matrix);
  free(all_local_vector);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
393

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


411

412
static void apply_boundary_conditions(const Mesh *mesh, int n_fields, double *solution, int nBnd, StrongBoundary *bnds)
413
{
414
415
  for (int ibnd = 0; ibnd < nBnd; ++ibnd) {
    StrongBoundary *bnd = bnds + ibnd;
416
417
418
419
420
421
422
423
    int iphys;
    for (iphys = 0; iphys < mesh->n_phys; ++iphys) {
      if (strcmp(bnd->tag, mesh->phys_name[iphys]) == 0)
        break;
    }
    if (iphys == mesh->n_phys){
      printf("Boundary tag \"%s\" not found.", bnd->tag);
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
424
    double *x = malloc(mesh->phys_n_nodes[iphys]*sizeof(double)*DIMENSION);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
425
426
427
    double *v = malloc(mesh->phys_n_nodes[iphys]*sizeof(double));
    for (int i = 0; i < mesh->phys_n_nodes[iphys]; ++i){
      int j = mesh->phys_nodes[iphys][i];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
428
429
      for (int k = 0; k < DIMENSION; ++k)
        x[i*DIMENSION+k] = mesh->x[j*3+k];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
430
431
    }
    bnd->apply(mesh->phys_n_nodes[iphys], x, v);
432
    for (int i = 0; i < mesh->phys_n_nodes[iphys]; ++i)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
433
434
435
      solution[mesh->phys_nodes[iphys][i]*n_fields+bnd->field] = v[i];
    free(x);
    free(v);
436
437
438
  }
}

439
int fluid_problem_implicit_euler(FluidProblem *problem, double dt)
440
{
441
442
  static double totaltime = 0;
  struct timespec tic, toc;
443
  const Mesh *mesh = problem->mesh;
444
445
  apply_boundary_conditions(mesh, problem->nFields, problem->solution, problem->n_strong_boundaries, problem->strong_boundaries);
  int n_fields = problem->nFields;
446
  double *solution_old = malloc(mesh->n_nodes*n_fields*sizeof(double));
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
447
  double *rhs = malloc(mesh->n_nodes*n_fields*sizeof(double));
Matthieu Constant's avatar
Matthieu Constant committed
448
  for (int i=0; i<mesh->n_nodes*n_fields; ++i){
449
    solution_old[i] = problem->solution[i];
Matthieu Constant's avatar
Matthieu Constant committed
450
451
    printf("old_solution=%g,solution=%g\n",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
    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);
Matthieu Constant's avatar
Matthieu Constant committed
462
    //exit(0);
Matthieu Constant's avatar
Matthieu Constant committed
463
   // if (i==1) break;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
464
    hxtLinearSystemGetRhsNorm(problem->linear_system,rhs,&norm);
465
466
467
468
469
470
471
    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;
472
  clock_get(&tic);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
473
    hxtLinearSystemSolve(problem->linear_system, rhs, dx);
474
475
  clock_get(&toc);
  totaltime += clock_diff(&tic, &toc);
476
477
478
    for (int j=0; j<mesh->n_nodes*n_fields; ++j) {
      problem->solution[j] -= dx[j];
    }
Matthieu Constant's avatar
Matthieu Constant committed
479
   /* for (int j=0; j<mesh->n_nodes*n_fields; ++j){
Matthieu Constant's avatar
Matthieu Constant committed
480
481
      printf("sol=%g\n",dx[j]);
    }
Matthieu Constant's avatar
Matthieu Constant committed
482
    exit(0);*/
483
  }
484
  printf("total solve time : %g\n", totaltime);
485
  free(dx);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
486
  free(rhs);
487
  free(solution_old);
488
  return i;
489
490
}

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

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

534
FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double *mu, double *rho, double epsilon, int n_fluids, int n_strong_boundaries, StrongBoundary *strong_boundaries)
535
{
536
537
538
539
540
  static int initialized = 0;
  if (!initialized) {
    hxtInitializeLinearSystems(0, NULL);
    initialized = 1;
#ifdef HXT_HAVE_PETSC
Matthieu Constant's avatar
Matthieu Constant committed
541
    hxtPETScInsertOptions("-pc_type lu -ksp_max_it 30", "fluid");
542
543
#endif
  }
544
  FluidProblem *problem = malloc(sizeof(FluidProblem));
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
545
546
547
  Mesh *mesh = mesh_load_msh(mesh_file_name);
  if (!mesh)
    return NULL;
548
  problem->mesh_tree = mesh_tree_create(mesh);
549
550
551
552
  problem->rho = malloc(n_fluids*sizeof(double));
  problem->mu = malloc(n_fluids*sizeof(double));
  for (int i = 0; i < n_fluids; ++i)
  {
Matthieu Constant's avatar
Matthieu Constant committed
553
    printf("mu=%g,rho=%g\n",mu[i],rho[i]);
554
555
556
557
558
559
    problem->rho[i] = rho[i];
    problem->mu[i] = mu[i];
  }
  problem->nFluids = n_fluids;
  problem->nFields = (DIMENSION+1)*n_fluids+1;
  int n_fields = problem->nFields;
560
  problem->g = g;
561
  problem->mesh = mesh;
Matthieu Constant's avatar
Matthieu Constant committed
562
  problem->epsilon = epsilon;
563
564
565
566
  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
567
    problem->strong_boundaries[i].apply = strong_boundaries[i].apply;
568
569
    problem->strong_boundaries[i].field = strong_boundaries[i].field;
  }
570
571
  problem->compacity = malloc(mesh->n_nodes*sizeof(double));
  problem->old_compacity = malloc(mesh->n_nodes*sizeof(double));
572
  problem->solution = malloc(mesh->n_nodes*n_fields*sizeof(double));
Matthieu Constant's avatar
Matthieu Constant committed
573
574
  for (int i = 0; i < problem->mesh->n_nodes; ++i)
  {
Matthieu Constant's avatar
debug    
Matthieu Constant committed
575
    problem->compacity[i] = 0.;
Matthieu Constant's avatar
Matthieu Constant committed
576
  }
577
  // begin to remove
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
578
  for (int i = 0; i < problem->mesh->n_nodes*n_fields; ++i){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
579
    problem->solution[i] = 0.;
580
  }
Matthieu Constant's avatar
Matthieu Constant committed
581
582
  problem->node_volume = malloc(0);
  fluid_problem_compute_node_volume(problem);
583
  // end to remove
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
584
  problem->linear_system = NULL;
585
#ifdef HXT_HAVE_PETSC
Matthieu Constant's avatar
Matthieu Constant committed
586
    hxtLinearSystemCreatePETSc(&problem->linear_system,mesh->n_elements,N_N,n_fields,mesh->elements,"fluid");
587
#else
588
    hxtLinearSystemCreateLU(&problem->linear_system,mesh->n_elements,N_N,n_fields,mesh->elements);
589
#endif
590
591
592
593
594
  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
595
  problem->particle_force = malloc(0);
596
597
  problem->particle_volume = malloc(0);
  problem->particle_velocity = malloc(0);
598
599
600
601
602
603
  return problem;
}

void fluid_problem_free(FluidProblem *problem)
{
  free(problem->solution);
604
605
  free(problem->compacity);
  free(problem->old_compacity);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
606
  hxtLinearSystemDelete(&problem->linear_system);
607
608
609
610
  free(problem->particle_uvw);
  free(problem->particle_element_id);
  free(problem->particle_position);
  free(problem->particle_mass);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
611
  free(problem->particle_force);
612
613
614
615
616
  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);
617
  mesh_tree_free(problem->mesh_tree);
618
619
620
  mesh_free(problem->mesh);
}

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


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

Matthieu Constant's avatar
Matthieu Constant committed
868
869
870
871
void fluid_problem_after_import(FluidProblem *problem)
{
  mesh_tree_free(problem->mesh_tree);
  problem->mesh_tree = mesh_tree_create(problem->mesh);
Matthieu Constant's avatar
debug    
Matthieu Constant committed
872
873
  free(problem->old_compacity);
  problem->old_compacity = malloc(sizeof(double)*problem->mesh->n_nodes);
Matthieu Constant's avatar
Matthieu Constant committed
874
875
876
877
878
879
880
  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
Matthieu Constant's avatar
debug    
Matthieu Constant committed
881
    hxtLinearSystemCreatePETSc(&problem->linear_system,problem->mesh->n_elements,N_N,problem->nFields,problem->mesh->elements,"fluid");
Matthieu Constant's avatar
Matthieu Constant committed
882
883
884
885
886
  #else
    hxtLinearSystemCreateLU(&problem->linear_system,new_mesh->n_elements,N_N,n_fields,new_mesh->elements);
  #endif 
}

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
887

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
888
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
889
{
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
890
  
891
  int n_fields = problem->nFields;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
892
  
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
893
894
895
896
897
898
  if(problem->n_particles != n)