Commit 0991333a authored by Michel Henry's avatar Michel Henry
Browse files

mapping to fix

parent c6bbcbb9
Pipeline #8313 failed with stages
in 1 minute and 5 seconds
...@@ -1142,7 +1142,6 @@ void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, const dou ...@@ -1142,7 +1142,6 @@ void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, const dou
printf("Boundary tag \"%s\" not found.\n", bnd->tag); printf("Boundary tag \"%s\" not found.\n", bnd->tag);
exit(1); exit(1);
} }
// To change for periodic conditions
for (int i = 0; i < problem->boundaries[iphys].n_nodes; ++i){ for (int i = 0; i < problem->boundaries[iphys].n_nodes; ++i){
forced_row[problem->boundaries[iphys].nodes[i]*n_fields+bnd->row] = bnd->field; forced_row[problem->boundaries[iphys].nodes[i]*n_fields+bnd->row] = bnd->field;
} }
...@@ -1164,18 +1163,14 @@ void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, const dou ...@@ -1164,18 +1163,14 @@ void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, const dou
const int *el = &mesh->elements[iel*N_N]; const int *el = &mesh->elements[iel*N_N];
char forced_field = forced_row[el[inode]*n_fields+ifield]; char forced_field = forced_row[el[inode]*n_fields+ifield];
if (forced_field == -1) continue; if (forced_field == -1) continue;
int i = inode*n_fields + ifield; // HERE : To periodic int i = inode*n_fields + ifield;
// int i = parent_tag[inode]*n_fields + ifield;
for (int jnode = 0; jnode< N_SF; ++jnode){ for (int jnode = 0; jnode< N_SF; ++jnode){
for (int jfield= 0; jfield< n_fields; ++jfield){ for (int jfield= 0; jfield< n_fields; ++jfield){
// int j = parent_tag[jnode]*n_fields+jfield
local_matrix[(inode*n_fields+ifield)*local_size+(jnode*n_fields+jfield)] = local_matrix[(inode*n_fields+ifield)*local_size+(jnode*n_fields+jfield)] =
(inode==jnode && jfield == forced_field) ? 1. : 0.; (inode==jnode && jfield == forced_field) ? 1. : 0.;
// local_matrix[i*local_size+j] = (inode==jnode && jfield == forced_field) ? 1. : 0.;
} }
} }
local_vector[inode+ifield*N_SF] = 0; local_vector[inode+ifield*N_SF] = 0;
// local_vector[parent_tag[inode]+ifield*N_SF] = 0;
} }
} }
} }
......
...@@ -562,7 +562,7 @@ void gmsh_mesh_read(GmshMesh *mesh, const char *file_name) ...@@ -562,7 +562,7 @@ void gmsh_mesh_read(GmshMesh *mesh, const char *file_name)
_gmsh_mesh_read_4(mesh, r); _gmsh_mesh_read_4(mesh, r);
else if (((int)version) == 2) else if (((int)version) == 2)
_gmsh_mesh_read_2(mesh, r); _gmsh_mesh_read_2(mesh, r);
else else
file_reader_error_at(r,"msh version is not supported"); file_reader_error_at(r,"msh version is not supported");
file_reader_delete(r); file_reader_delete(r);
} }
......
# MigFlow - Copyright (C) <2010-2018> # MigFlow - Copyright (C) <2010-2018>
# <Universite catholique de Louvain (UCL), Belgium # <Universite catholique de Louvain (UCL), Belgium
# Universite de Montpellier, France> # Universite de Montpellier, France>
# #
# List of the contributors to the development of MigFlow: see AUTHORS file. # List of the contributors to the development of MigFlow: see AUTHORS file.
# Description and complete License: see LICENSE file. # Description and complete License: see LICENSE file.
# #
# This program (MigFlow) is free software: # This program (MigFlow) is free software:
# you can redistribute it and/or modify it under the terms of the GNU Lesser General # 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 # Public License as published by the Free Software Foundation, either version
# 3 of the License, or (at your option) any later version. # 3 of the License, or (at your option) any later version.
# #
# This program is distributed in the hope that it will be useful, # This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of # but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details. # GNU Lesser General Public License for more details.
# #
# You should have received a copy of the GNU Lesser General Public License # 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, # along with this program (see COPYING and COPYING.LESSER files). If not,
# see <http://www.gnu.org/licenses/>. # see <http://www.gnu.org/licenses/>.
"""Model for Immersed Granular Flow -- Fluid user interface """Model for Immersed Granular Flow -- Fluid user interface
...@@ -24,7 +24,7 @@ ...@@ -24,7 +24,7 @@
Contact: jonathan.lambrechts@uclouvain.be Contact: jonathan.lambrechts@uclouvain.be
Webpage: www.migflow.be Webpage: www.migflow.be
MigFlow computes immersed granular flows using an unresolved FEM-DEM model. MigFlow computes immersed granular flows using an unresolved FEM-DEM model.
The continuous phase (mixture) is computed by solving averaged Navier-Stokes equations on unstructured meshes with the continuous finite element method. The continuous phase (mixture) is computed by solving averaged Navier-Stokes equations on unstructured meshes with the continuous finite element method.
""" """
...@@ -73,7 +73,7 @@ def _np2c(a,dtype=float) : ...@@ -73,7 +73,7 @@ def _np2c(a,dtype=float) :
r.tmp = tmp r.tmp = tmp
return r return r
class FluidProblem : class FluidProblem :
"""Create numerical representation of the fluid.""" """Create numerical representation of the fluid."""
...@@ -86,7 +86,7 @@ class FluidProblem : ...@@ -86,7 +86,7 @@ class FluidProblem :
viscosity -- list of fluid phases viscosities (this should be a vector whom dimension specify if there is one or two fluid) viscosity -- list of fluid phases viscosities (this should be a vector whom dimension specify if there is one or two fluid)
density -- list of fluid phases densities (this should be a vector whom dimension specify if there is one or two fluid) density -- list of fluid phases densities (this should be a vector whom dimension specify if there is one or two fluid)
sigma -- surface tension (only when there are two fluids) sigma -- surface tension (only when there are two fluids)
coeff_stab -- optional argument used as coefficient for extra diffusivity added to stabilise the advection of concentration (only for two fluid problem) coeff_stab -- optional argument used as coefficient for extra diffusivity added to stabilise the advection of concentration (only for two fluid problem)
petsc_solver_type -- optional argument to specify solver option (only when petsc is used) petsc_solver_type -- optional argument to specify solver option (only when petsc is used)
drag_in_stab -- state if the drag force is in the stabilisation term drag_in_stab -- state if the drag force is in the stabilisation term
drag_coefficient_factor -- factor multiplicating the drag coefficient drag_coefficient_factor -- factor multiplicating the drag coefficient
...@@ -323,7 +323,7 @@ class FluidProblem : ...@@ -323,7 +323,7 @@ class FluidProblem :
dt -- computation time step dt -- computation time step
""" """
torques = np.ndarray([self.n_particles,1],"d",order="C") torques = np.ndarray([self.n_particles,1],"d",order="C")
self._lib.fluid_problem_compute_node_particle_torque(self._fp, c_double(dt), c_void_p(torques.ctypes.data)) self._lib.fluid_problem_compute_node_particle_torque(self._fp, c_double(dt), c_void_p(torques.ctypes.data))
return torques return torques
...@@ -344,17 +344,22 @@ class FluidProblem : ...@@ -344,17 +344,22 @@ class FluidProblem :
rhs = np.zeros_like(solution_old) rhs = np.zeros_like(solution_old)
self._lib.fluid_problem_reset_boundary_force(self._fp) self._lib.fluid_problem_reset_boundary_force(self._fp)
self._lib.fluid_problem_compute_stab_parameters(self._fp,c_double(dt)) self._lib.fluid_problem_compute_stab_parameters(self._fp,c_double(dt))
# Create mapping
map = self.
if self.sys is None : if self.sys is None :
self.sys = LinearSystem(self.elements(),self.n_fields(),self.solver_options self.sys = LinearSystem(map[self.elements()],self.n_fields(),self.solver_options)
#self.sys = LinearSystem(self.elements(),self.n_fields(),self.solver_options
) )
sys = self.sys sys = self.sys
localv = np.ndarray([sys.localsize*self.n_elements()]) localv = np.ndarray([sys.localsize*self.n_elements()])
localm = np.ndarray([sys.localsize**2*self.n_elements()]) localm = np.ndarray([sys.localsize**2*self.n_elements()])
self._lib.fluid_problem_assemble_system(self._fp,_np2c(rhs),_np2c(solution_old),c_double(dt),_np2c(localv),_np2c(localm)) self._lib.fluid_problem_assemble_system(self._fp,_np2c(rhs),_np2c(solution_old),c_double(dt),_np2c(localv),_np2c(localm))
sys.local_to_global(localv,localm,rhs) sys.local_to_global(localv,localm,rhs)
# Periodic Conditions
if check_residual_norm > 0 : if check_residual_norm > 0 :
norm0 = np.linalg.norm(rhs) norm0 = np.linalg.norm(rhs)
self.solution()[:] -= sys.solve(rhs) self.solution()[:] -= sys.solve(rhs)[map]
#self.solution()[:] -= sys.solve(rhs)
if check_residual_norm > 0 : if check_residual_norm > 0 :
sys.zero_matrix() sys.zero_matrix()
rhs[:] = 0 rhs[:] = 0
...@@ -396,7 +401,7 @@ class FluidProblem : ...@@ -396,7 +401,7 @@ class FluidProblem :
Keyword arguments: Keyword arguments:
mass -- list of particles mass mass -- list of particles mass
volume -- list of particles volume volume -- list of particles volume
position -- list of particles centre positions position -- list of particles centre positions
velocity -- list of particles velocity velocity -- list of particles velocity
init -- if True set the old_porosity to zero (has to be used at the first call) init -- if True set the old_porosity to zero (has to be used at the first call)
reload -- optional arguments specifying if the simulation is reloaded or if it is the first iteration reload -- optional arguments specifying if the simulation is reloaded or if it is the first iteration
...@@ -441,7 +446,7 @@ class FluidProblem : ...@@ -441,7 +446,7 @@ class FluidProblem :
"""Return the number of mesh nodes.""" """Return the number of mesh nodes."""
self._lib.fluid_problem_n_elements.restype = c_int self._lib.fluid_problem_n_elements.restype = c_int
return self._lib.fluid_problem_n_elements(self._fp); return self._lib.fluid_problem_n_elements(self._fp);
def n_fields(self): def n_fields(self):
"""Return the number of fluid fields.""" """Return the number of fluid fields."""
self._lib.fluid_problem_n_fields.restype = c_int self._lib.fluid_problem_n_fields.restype = c_int
...@@ -451,7 +456,7 @@ class FluidProblem : ...@@ -451,7 +456,7 @@ class FluidProblem :
"""Return the number of mesh nodes.""" """Return the number of mesh nodes."""
self._lib.fluid_problem_n_nodes.restype = c_int self._lib.fluid_problem_n_nodes.restype = c_int
return self._lib.fluid_problem_n_nodes(self._fp); return self._lib.fluid_problem_n_nodes(self._fp);
def porosity(self): def porosity(self):
"""Return the porosity at nodes""" """Return the porosity at nodes"""
return self._get_matrix("porosity", self.n_nodes(), 1) return self._get_matrix("porosity", self.n_nodes(), 1)
...@@ -476,7 +481,7 @@ class FluidProblem : ...@@ -476,7 +481,7 @@ class FluidProblem :
def old_porosity(self): def old_porosity(self):
"""Return the porosity at previous time step""" """Return the porosity at previous time step"""
return self._get_matrix("old_porosity", self.n_nodes(), 1) return self._get_matrix("old_porosity", self.n_nodes(), 1)
def volume_flux(self, bnd_tag): def volume_flux(self, bnd_tag):
"""Compute the integral of the (outgoing) normal velocity through boundary with tag bnd_tag""" """Compute the integral of the (outgoing) normal velocity through boundary with tag bnd_tag"""
self._lib.fluid_problem_volume_flux.restype = c_double self._lib.fluid_problem_volume_flux.restype = c_double
...@@ -528,7 +533,7 @@ class FluidProblem : ...@@ -528,7 +533,7 @@ class FluidProblem :
phys_n_boundaries = c_int() phys_n_boundaries = c_int()
phys_boundaries = POINTER(c_int)() phys_boundaries = POINTER(c_int)()
f(self._fp,c_int(ibnd),byref(phys_name),byref(phys_dim),byref(phys_id),byref(phys_n_nodes),byref(phys_nodes),byref(phys_n_boundaries),byref(phys_boundaries)) f(self._fp,c_int(ibnd),byref(phys_name),byref(phys_dim),byref(phys_id),byref(phys_n_nodes),byref(phys_nodes),byref(phys_n_boundaries),byref(phys_boundaries))
nodes = np.ctypeslib.as_array(phys_nodes,shape=(phys_n_nodes.value,)) if phys_nodes else np.ndarray([0],np.int32) nodes = np.ctypeslib.as_array(phys_nodes,shape=(phys_n_nodes.value,)) if phys_nodes else np.ndarray([0],np.int32)
bnd = np.ctypeslib.as_array(phys_boundaries,shape=(phys_n_boundaries.value,2)) if phys_boundaries else np.ndarray([2,0],np.int32) bnd = np.ctypeslib.as_array(phys_boundaries,shape=(phys_n_boundaries.value,2)) if phys_boundaries else np.ndarray([2,0],np.int32)
return phys_name.value,phys_dim.value,phys_id.value,nodes,bnd return phys_name.value,phys_dim.value,phys_id.value,nodes,bnd
......
Supports Markdown
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