Commit 893ac506 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

smooth porosity

parent 4d279a1a
Pipeline #5457 passed with stage
in 54 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")
......
......@@ -1143,20 +1143,94 @@ static void fluid_problem_compute_node_volume(FluidProblem *problem) {
}
}
#include "vector.h"
static void fluid_problem_compute_porosity(FluidProblem *problem)
{
Mesh *mesh = problem->mesh;
int *touched_el_flag = malloc(sizeof(int)*mesh->n_elements);
double L2 = pow2(0.002);
int *nodes = NULL;
int *stack = NULL;
double *nodes_w = NULL;
int *neighbours = malloc(sizeof(int)*mesh->n_elements*(D+1));
for (int i = 0; i< mesh->n_elements; ++i){
touched_el_flag[i] = 0;
for (int d = 0; d < D+1; ++d)
neighbours[i*(D+1)+d] = -1;
}
for (int i = 0; i < mesh->n_interfaces; ++i) {
int *inter = mesh->interfaces + i*4;
if (inter[2] < 0) continue;
neighbours[inter[0]*(D+1)+inter[1]%(D+1)] = inter[2];
neighbours[inter[2]*(D+1)+inter[3]%(D+1)] = inter[0];
}
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;
double sf[N_SF];
shape_functions(&problem->particle_uvw[i*D],sf);
const int *el = &mesh->elements[problem->particle_element_id[i]*N_N];
for (int j=0; j<N_N;++j)
problem->porosity[el[j]] += sf[j]*volume[i];
double *xp = problem->particle_position+i*D;
vectorClear(nodes);
vectorClear(nodes_w);
vectorClear(stack);
int p = 0;
double wtot = 0;
*vectorPush(&stack) = problem->particle_element_id[i];
touched_el_flag[problem->particle_element_id[i]] = 1;
while (vectorSize(stack) > p) {
int elid = stack[p];
int close = 0;
for (int j = 0; j < N_N; ++j) {
int jnode = mesh->elements[elid*N_N+j];
double *x = mesh->x + jnode*3;
double d2 = 0;
for (int k = 0; k < D; ++k)
d2 += pow2(x[k]-xp[k]);
if (d2 < L2) {
close = 1;
int found = 0;
for (int l = 0; l < vectorSize(nodes); ++l){
if (nodes[l] == jnode) {
found = 1;
break;
}
}
if (found) continue;
*vectorPush(&nodes) = jnode;
double w = pow2(sqrt(L2)-sqrt(d2))*problem->node_volume[jnode];
*vectorPush(&nodes_w) = w;
wtot += w;
}
}
if (close) {
for (int j = 0; j < N_N; ++j) {
int nid = neighbours[elid*(D+1)+j];
if (nid == -1 || touched_el_flag[nid]!=0) continue;
touched_el_flag[nid] = 1;
*vectorPush(&stack) = nid;
}
}
p += 1;
}
if (vectorSize(nodes) >= 3) {
for (size_t j = 0; j < vectorSize(nodes); ++j) {
problem->porosity[nodes[j]] += volume[i]*nodes_w[j]/wtot;
}
}
else {
double sf[N_SF];
shape_functions(&problem->particle_uvw[i*D],sf);
const int *el = &mesh->elements[problem->particle_element_id[i]*N_N];
for (int j=0; j<N_N;++j){
problem->porosity[el[j]] += sf[j]*volume[i];
}
}
for (int j = 0; j < vectorSize(stack); ++j) {
touched_el_flag[stack[j]] = 0;
}
}
for (int i = 0; i < mesh->n_nodes; ++i){
if(problem->node_volume[i]==0.){
......@@ -1166,6 +1240,10 @@ static void fluid_problem_compute_porosity(FluidProblem *problem)
problem->porosity[i] = 1-problem->porosity[i]/problem->node_volume[i];
}
}
free(neighbours);
vectorFree(nodes);
vectorFree(nodes_w);
vectorFree(stack);
}
FluidProblem *fluid_problem_new(double g, int n_fluids, double *mu, double *rho, double sigma, double coeffStab, double volume_drag, const char *petsc_solver_type, int drag_in_stab)
......
......@@ -302,6 +302,9 @@ class FluidProblem :
self._lib.fluid_problem_set_reduced_gravity(self._fp,c_int(reduced_gravity))
self._lib.fluid_problem_set_stab_param(self._fp,c_double(stab_param))
nit = self._lib.fluid_problem_implicit_euler(self._fp, c_double(dt), c_double(newton_atol), c_double(newton_rtol), c_int(newton_max_it))
return nit
def advance_concentration(self,dt):
if self._n_fluids == 2 :
nv = np.linalg.norm(self.solution()[:,:self._dim],axis=1)
nvmax = np.max(nv[self.elements()],axis=1,keepdims=True)
......@@ -312,7 +315,6 @@ class FluidProblem :
print("sub-iterating advection for cfl : %i sub-iterations"%nsub)
for i in range(nsub) :
self._lib.fluid_problem_advance_concentration(self._fp,c_double(dt/nsub))
return nit
def set_particles(self, mass, volume, position, velocity, contact, init=False, reload = False):
"""Set location of the grains in the mesh and compute the porosity in each cell.
......
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