separation Ã2 fluides

src/fluid_problem.c

@@ -248,6 +248,7 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d
       const int iel = bnd[0];
       const int i0 = bnd[1];
       const int i1 = (bnd[1]+1) % 3;
+      const int i2 = (bnd[1]+2) % 3;
       const uint32_t *el = &mesh->elements[iel*N_N];
       const int nodes[2] = {el[i0],el[i1]};
       const double *x[2] = {&mesh->x[nodes[0]*3], &mesh->x[nodes[1]*3]};
@@ -287,8 +288,10 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d
 
           local_matrix[(i0*n_fields+Q)*local_size + (i0*n_fields+P)] += dphi[i0][k]*n[k]*problem->epsilon/2.*l;
           local_matrix[(i0*n_fields+Q)*local_size + (i1*n_fields+P)] += dphi[i1][k]*n[k]*problem->epsilon/2.*l;
+          local_matrix[(i0*n_fields+Q)*local_size + (i2*n_fields+P)] += dphi[i2][k]*n[k]*problem->epsilon/2.*l;
           local_matrix[(i1*n_fields+Q)*local_size + (i0*n_fields+P)] += dphi[i0][k]*n[k]*problem->epsilon/2.*l;
           local_matrix[(i1*n_fields+Q)*local_size + (i1*n_fields+P)] += dphi[i1][k]*n[k]*problem->epsilon/2.*l;
+          local_matrix[(i1*n_fields+Q)*local_size + (i2*n_fields+P)] += dphi[i2][k]*n[k]*problem->epsilon/2.*l;
           local_vector[i0+Q*N_SF] += dp[k]*n[k]*problem->epsilon/2.*l;
           local_vector[i1+Q*N_SF] += dp[k]*n[k]*problem->epsilon/2.*l;
           

testcases/separation/sep.py

+# Marblesbag - Copyright (C) <2010-2018>
+# <Universite catholique de Louvain (UCL), Belgium
+#  Universite de Montpellier, France>
+# 	
+# List of the contributors to the development of Marblesbag: see AUTHORS file.
+# Description and complete License: see LICENSE file.
+# 	
+# This program (Marblesbag) is free software: 
+# you can redistribute it and/or modify it under the terms of the GNU Lesser General 
+# 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
+# GNU Lesser General Public License for more details.
+# 
+# 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, 
+# see <http://www.gnu.org/licenses/>.
+
+#!/usr/bin/env python
+from marblesbag import fluid as fluid
+from marblesbag import scontact2
+
+import numpy as np
+import os
+import time
+import shutil
+import random
+#Physical parameters for the drops are the ones presented by Metzger et al. (2007) "Falling clouds of particles in viscous fluids"
+
+
+outputdir = "output"
+if not os.path.isdir(outputdir) :
+    os.makedirs(outputdir)
+
+t = 0
+ii = 0
+
+
+#physical parameters
+g =  -9.81                                      # gravity
+rho0 = 1000                                      # fluid density
+rho1 = 1000 
+nu0 = 1e-3                                   # kinematic viscosity
+nu1 = 1e-3                                   # kinematic viscosity
+tEnd = 100000                                     # final time
+
+#numerical parameters
+lcmin = .1                                  # mesh size
+dt = 1                                       # time step
+alpha = 1e-4                                    # stabilization coefficient
+epsilon = alpha*lcmin**2 /min(nu0,nu1)                    # stabilization parametre
+print('epsilon',epsilon)
+
+shutil.copy("mesh.msh", outputdir +"/mesh.msh")
+
+outf = 1                                       # number of iterations between output files
+
+#Object fluid creation + Boundary condition of the fluid (field 0 is horizontal velocity; field 1 is vertical velocity; field 2 is pressure)
+#Format: strong_boundaries = [(Boundary tag, Fluid field, Value)
+u0,v0,q0 = 0,1,2
+u1,v1,q1 = 3,4,5
+p = 6
+strong_boundaries = [
+        ("Bottom",u0,u0,0),
+        ("Bottom",v0,v0,0.),
+        ("Bottom",u1,u1,0),
+        ("Bottom",v1,v1,0.),
+        ("Top",u0,u0,0),
+        ("Top",v0,v0,0.),
+        ("Top",u1,u1,0),
+        ("Top",v1,v1,0.),
+        ("Top",q0,p,0.),
+        ("Left",u0,u0,0),
+        ("Left",v0,v0,0),
+        ("Left",u1,u1,0),
+        ("Left",v1,v1,0),
+        ("Right",u0,v0,0),
+        ("Right",u1,v1,0),
+        ("Right",v0,v0,0),
+        ("Right",v1,v1,0),
+        ]
+
+fluid = fluid.fluid_problem("mesh.msh",g,[nu0*rho0,nu1*rho1],[rho0,rho1],epsilon,strong_boundaries,2)
+ii = 0
+t = 0
+
+#set initial_condition
+x = fluid.coordinates()
+s = fluid.solution()
+a = (x[:,0]+5)/10
+s[:,2] = 0.5
+s[:,5] = 1-s[:,2]
+
+
+
+
+fluid.export_vtk(outputdir,0,0)
+
+ii = 0
+tic = time.clock()
+while t < tEnd : 
+    #Fluid solver
+    fluid.implicit_euler(dt)
+    t += dt
+    s = fluid.solution()
+    x = fluid.coordinates()
+    vel = (s[:,0]+s[:,3]-1/(20*nu0*rho0)*x[:,1]*(1-x[:, 1]))**2
+    print('Error', (vel.sum())**.5)
+    #Output files writting
+    if ii %outf == 0 :
+        ioutput = int(ii/outf) + 1
+        fluid.export_vtk(outputdir, t, ioutput)
+    ii += 1
+    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))
\ No newline at end of file