Commit e5136752 authored by Matthieu Constant's avatar Matthieu Constant
Browse files

merge

parents ec66243b 560ba5d2
Pipeline #5531 failed with stage
in 1 minute and 38 seconds
from distutils.sysconfig import get_python_inc
import platform
import os
import ycm_core
# These are the compilation flags that will be used in case there's no
# compilation database set (by default, one is not set).
# CHANGE THIS LIST OF FLAGS. YES, THIS IS THE DROID YOU HAVE BEEN LOOKING FOR.
flags = [
'-Wno-unused-function',
'-std=gnu99',
'-Iscontact',
'-Ihxt',
'-Isrc',
'-DDIMENSION=2'
]
# Clang automatically sets the '-std=' flag to 'c++14' for MSVC 2015 or later,
# which is required for compiling the standard library, and to 'c++11' for older
# versions.
#if platform.system() != 'Windows':
# flags.append( '-std=c++11' )
# Set this to the absolute path to the folder (NOT the file!) containing the
# compile_commands.json file to use that instead of 'flags'. See here for
# more details: http://clang.llvm.org/docs/JSONCompilationDatabase.html
#
# You can get CMake to generate this file for you by adding:
# set( CMAKE_EXPORT_COMPILE_COMMANDS 1 )
# to your CMakeLists.txt file.
#
# Most projects will NOT need to set this to anything; you can just change the
# 'flags' list of compilation flags. Notice that YCM itself uses that approach.
compilation_database_folder = ''
if os.path.exists( compilation_database_folder ):
database = ycm_core.CompilationDatabase( compilation_database_folder )
else:
database = None
SOURCE_EXTENSIONS = [ '.cpp', '.cxx', '.cc', '.c', '.m', '.mm' ]
def DirectoryOfThisScript():
return os.path.dirname( os.path.abspath( __file__ ) )
def IsHeaderFile( filename ):
extension = os.path.splitext( filename )[ 1 ]
return extension in [ '.h', '.hxx', '.hpp', '.hh' ]
def FindCorrespondingSourceFile( filename ):
if IsHeaderFile( filename ):
basename = os.path.splitext( filename )[ 0 ]
for extension in SOURCE_EXTENSIONS:
replacement_file = basename + extension
if os.path.exists( replacement_file ):
return replacement_file
return filename
def FlagsForFile( filename, **kwargs ):
# If the file is a header, try to find the corresponding source file and
# retrieve its flags from the compilation database if using one. This is
# necessary since compilation databases don't have entries for header files.
# In addition, use this source file as the translation unit. This makes it
# possible to jump from a declaration in the header file to its definition in
# the corresponding source file.
filename = FindCorrespondingSourceFile( filename )
if not database:
return {
'flags': flags,
'include_paths_relative_to_dir': DirectoryOfThisScript(),
'override_filename': filename
}
compilation_info = database.GetCompilationInfoForFile( filename )
if not compilation_info.compiler_flags_:
return None
# Bear in mind that compilation_info.compiler_flags_ does NOT return a
# python list, but a "list-like" StringVec object.
final_flags = list( compilation_info.compiler_flags_ )
# NOTE: This is just for YouCompleteMe; it's highly likely that your project
# does NOT need to remove the stdlib flag. DO NOT USE THIS IN YOUR
# ycm_extra_conf IF YOU'RE NOT 100% SURE YOU NEED IT.
try:
final_flags.remove( '-stdlib=libc++' )
except ValueError:
pass
return {
'flags': final_flags,
'include_paths_relative_to_dir': compilation_info.compiler_working_dir_,
'override_filename': filename
}
......@@ -35,7 +35,7 @@ if(PETSC_FOUND)
list(APPEND SRC hxt_linear_system_petsc.c)
endif(PETSC_FOUND)
include_directories(. ../tools ${EXTERNAL_INCLUDES})
include_directories(. ${CMAKE_SOURCE_DIR}/tools ${EXTERNAL_INCLUDES})
add_library(mbfluid2 SHARED ${SRC})
target_link_libraries(mbfluid2 ${EXTERNAL_LIBRARIES} mbtools2)
target_compile_definitions(mbfluid2 PUBLIC "-DDIMENSION=2")
......
This diff is collapsed.
......@@ -87,6 +87,13 @@ struct FluidProblem {
double *particle_force;
double *grad_a_cg;
double *all_kappa;
double *particle_node_forces;
double *particle_node_forces_du;
double *particle_nodes_weight;
double *force_cst, *force_dp, *force_v;
double *gamma;
int *particle_nodes;
int *particle_nodes_n;
int n_fluids;
//stabilisation coefficients
double *taup;
......
......@@ -812,6 +812,8 @@ static void _particleProblemGenContacts2(ParticleProblem *p, const double alert)
vectorFree(cc);
}
static void _particleProblemGenContacts(ParticleProblem *p, const double alert)
{
//_particleProblemGenContacts2(p,alert);
......@@ -999,7 +1001,7 @@ static void _particleProblemSolveContacts(ParticleProblem *p, double dt, double
}
}
static void _particleProblemSolveContactsQueue(ParticleProblem *p, double dt, double tol) {
static void _genContactByParticle(ParticleProblem *p, size_t **contactByParticlePPtr, size_t **contactByParticlePtr) {
size_t *nContactByParticle = malloc(sizeof(size_t)*(vectorSize(p->particles)));
for (size_t i = 0; i < vectorSize(p->particles); ++i) {
nContactByParticle[i] = 0;
......@@ -1024,8 +1026,67 @@ static void _particleProblemSolveContactsQueue(ParticleProblem *p, double dt, do
contactByParticle[contactByParticleP[c->o1]+(nContactByParticle[c->o1]++)] = i;
}
free(nContactByParticle);
*contactByParticlePtr = contactByParticle;
*contactByParticlePPtr = contactByParticleP;
}
static void _particleProblemComputeStressTensor(ParticleProblem *p, const double *nodes, const double radius, const int n_nodes)
{
size_t *contactByParticle, *contactByParticleP;
_genContactByParticle(p, &contactByParticleP, &contactByParticle);
double bbmin[DIMENSION], bbmax[DIMENSION];
_particleProblemBoundingBox(p, bbmin, bbmax);
Cell *tree = cellNew(bbmin, bbmax);
double amin[DIMENSION], amax[DIMENSION];
int *found = NULL;
vectorPushN(&found, 100);
vectorClear(found);
// Particles
Contact *oldContacts = vectorDup(p->contacts), *cold;
vectorClear(p->contacts);
int ntot = 0;
for (size_t i = 0; i < vectorSize(p->particles); ++i) {
particleBoundingBox(&p->particles[i], &p->position[i * DIMENSION], amin, amax);
vectorClear(found);
cellAdd(tree, amin, amax, i, &found);
}
for (int i = 0; i < n_nodes; ++i){
for (int j = 0; j < DIMENSION; ++j){
amin[j] = nodes[i*DIMENSION+j]-radius;
amax[j] = nodes[i*DIMENSION+j]+radius;
}
vectorClear(found);
cellSearch(tree,amin,amax,&found);
for (size_t j = 0; j < vectorSize(found); ++j){
for (int k = contactByParticleP[found[j]]; k < contactByParticleP[found[j]+1]; ++k){
Contact *c = &p->contacts[contactByParticle[k]];
double *xg = &p->position[c->o1*DIMENSION];
double r = p->particles[c->o1].r;
double xc[DIMENSION];
double d = 0;
for (int l = 0; l < DIMENSION; ++l){
xc[l] = xg[l]-r*c->n[l];
d += (xc[l]-nodes[i*DIMENSION+l])*(xc[l]-nodes[i*DIMENSION+l]);
}
d = sqrt(d);
if (d<radius){
double F = 0;
}
else{
continue;
}
}
}
}
vectorFree(found);
cellFree(tree);
}
static void _particleProblemSolveContactsQueue(ParticleProblem *p, double dt, double tol) {
Fifo *queue = fifoNew(vectorSize(p->contacts));
size_t *contactByParticle, *contactByParticleP;
_genContactByParticle(p, &contactByParticleP, &contactByParticle);
int *activeContact = malloc(sizeof(int)*vectorSize(p->contacts));
for (size_t i = 0; i < vectorSize(p->contacts); ++i) {
fifoPush(queue, i);
......
......@@ -75,6 +75,8 @@ class Darcy(unittest.TestCase) :
n_test = 10
c_v = np.linspace(0.01,0.85,n_test)
c_v = [.5]
n_test=1
#c_v = np.array([0.1])
pM = np.zeros_like(c_v)
pm = np.zeros_like(c_v)
......@@ -84,13 +86,13 @@ class Darcy(unittest.TestCase) :
g = 0 # gravity
rho = 1000 # fluid density
rhop = 1500 # grains density
L = 2
L = 1
nu = 1e-3 # kinematic viscosity
mu = 1
r = 5e-3 # grains radius
#numerical parameters
dt = 0.001 # time step
tEnd = dt*100 # final time
tEnd = dt*200 # final time
outf = 10 # number of iterations between output files
for compacity in c_v:
a_g = np.pi*r**2
......@@ -108,8 +110,10 @@ class Darcy(unittest.TestCase) :
fluid.solution()[:,0] = V
fluid.set_open_boundary("Left",velocity=[V,0])
fluid.set_open_boundary("Right",pressure=0)
fluid.set_wall_boundary("Bottom")
fluid.set_wall_boundary("Top")
fluid.set_wall_boundary("Bottom",velocity=[0,0])
fluid.set_wall_boundary("Top",velocity=[0,0])
fluid.set_strong_boundary("Top",0,0)
fluid.set_strong_boundary("Bottom",0,0)
fp = np.zeros_like(fluid.porosity())
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces(),init=True)
fluid.export_vtk(outputdir,0,0)
......@@ -121,8 +125,8 @@ class Darcy(unittest.TestCase) :
#Computation loop
while t < tEnd :
#Fluid solver
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),-fluid.compute_node_force(dt),init=True)
it = fluid.implicit_euler(dt)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),-fluid.compute_node_force(dt))
t += dt
#Output files writting
if ii %outf == 0 :
......@@ -138,10 +142,12 @@ class Darcy(unittest.TestCase) :
Dp_drag[j] = gamma*N*V/(1-compacity)**2
pressure = fluid.solution()[:,2]
pM[j] = np.max(pressure)
pm[j] = np.min(pressure)
pM[j] = np.mean(fluid.solution()[fluid.coordinates()[:,0]<-.99,2])
pm[j] = np.mean(fluid.solution()[fluid.coordinates()[:,0]>.99,2])
j = j+1
dp = (pM-pm)/L
dp = (pM-pm)/(2*L)
print(dp)
print(Dp_drag)
error = np.sum(abs(dp-Dp_drag)/Dp_drag)/n_test
self.assertLess(error,1./10., "error is too large in darcy")
\ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment