Skip to content
GitLab
Menu
Projects
Groups
Snippets
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
fluidparticles
MigFlow
Commits
e4b05067
Commit
e4b05067
authored
Jun 25, 2020
by
Nathan Coppin
Browse files
Updated docu
parent
b089bbb3
Pipeline
#7779
failed with stages
in 2 minutes and 51 seconds
Changes
2
Pipelines
1
Expand all
Hide whitespace changes
Inline
Side-by-side
python/fluid.py
View file @
e4b05067
This diff is collapsed.
Click to expand it.
python/scontact.py
View file @
e4b05067
...
...
@@ -60,7 +60,7 @@ def _np2c(a,dtype=np.float64) :
class
ParticleProblem
:
"""Create the numerical structure containing all the physical particles that appear in the problem"""
"""Create
s
the numerical structure containing all the physical particles that appear in the problem"""
def
_get_array
(
self
,
fName
,
dtype
)
:
f
=
getattr
(
self
.
_lib
,
"particleProblem"
+
fName
)
...
...
@@ -93,14 +93,14 @@ class ParticleProblem :
return
np
.
ctypeslib
.
as_array
(
buf
)
def
__init__
(
self
,
dim
,
friction_enabled
=
False
,
rotation_enabled
=
False
)
:
"""Build the particles contanier structure.
"""Build
s
the particles contanier structure.
Keyword arguments:
dim --
d
imension of the problem
dim --
D
imension of the problem
Raises:
ValueError --
i
f dimension differs from 2 or 3
NameError --
i
f C builder cannot be found
ValueError --
I
f dimension differs from 2 or 3
NameError --
I
f C builder cannot be found
"""
self
.
_dim
=
dim
self
.
_friction_enabled
=
friction_enabled
...
...
@@ -133,72 +133,86 @@ class ParticleProblem :
self
.
_triangletype
=
np
.
dtype
(
bndtype
+
[(
'p'
,
np
.
float64
,(
3
,
dim
))])
def
__del__
(
self
):
"""Delete the particle solver structure."""
"""Delete
s
the particle solver structure."""
if
(
self
.
_p
is
not
None
)
:
self
.
_lib
.
particleProblemDelete
(
self
.
_p
)
def
volume
(
self
):
"""Return the list of particle volumes."""
"""Return
s
the list of particle volumes."""
if
self
.
_dim
==
2
:
return
np
.
pi
*
self
.
_get_matrix
(
"Particle"
,
2
)[:,
[
0
]]
**
2
else
:
return
4.
/
3.
*
np
.
pi
*
self
.
_get_matrix
(
"Particle"
,
2
)[:,
[
0
]]
**
3
def
r
(
self
)
:
"""Return the list of particle radii."""
"""Return
s
the list of particle radii."""
return
self
.
_get_matrix
(
"Particle"
,
2
)[:,
0
][:,
None
]
def
mass
(
self
):
"""Return the list of particle masses."""
"""Return
s
the list of particle masses."""
return
self
.
_get_matrix
(
"Particle"
,
2
)[:,
1
][:,
None
]
def
position
(
self
):
"""Return the list of particle centre position vectors."""
"""Return
s
the list of particle centre position vectors."""
return
self
.
_get_matrix
(
"Position"
,
self
.
_dim
)
def
velocity
(
self
):
"""Return the list of particle velocity vectors."""
"""Return
s
the list of particle velocity vectors."""
return
self
.
_get_matrix
(
"Velocity"
,
self
.
_dim
)
def
omega
(
self
):
"""Return the list of particle angular velocity."""
"""Return
s
the list of particle angular velocity."""
assert
self
.
_friction_enabled
and
self
.
_rotation_enabled
return
(
self
.
_get_matrix
(
"Omega"
,
1
))
if
self
.
_dim
==
2
else
self
.
_get_matrix
(
"Omega"
,
3
)
def
particle_material
(
self
):
"""Returns the list of particle materials."""
return
self
.
_get_idx
(
"ParticleMaterial"
)
def
disk_tag
(
self
):
"""Returns the list of boundary disk tag numbers."""
return
self
.
_get_idx
(
"DiskTag"
)
def
disks
(
self
)
:
"""Returns the list of boundary disks."""
return
self
.
_get_array
(
"Disk"
,
self
.
_disktype
)
def
forced_flag
(
self
):
"""Returns the list of particle flags indicating whether they are fixed or not."""
return
self
.
_get_idx
(
"ForcedFlag"
)
def
get_tag_id
(
self
,
tag
=
"default"
):
"""Returns the number associated to a string tag."""
return
self
.
_lib
.
particleProblemGetMaterialTagId
(
self
.
_p
,
tag
.
encode
())
def
n_particles
(
self
)
:
"""Returns the number of particles."""
self
.
_lib
.
particleProblemNParticle
.
restype
=
c_size_t
return
self
.
_lib
.
particleProblemNParticle
(
self
.
_p
)
def
mu
(
self
)
:
return
self
.
_get_matrix
(
"Mu"
,
1
)
"""Returns the matrix of the friction coefficients between materials."""
return
self
.
_get_matrix
(
"Mu"
,
self
.
_lib
.
particleProblemNMaterial
(
self
.
_p
))
def
segments
(
self
):
"""Returns the list of boundary segments."""
return
self
.
_get_array
(
"Segment"
,
self
.
_segmenttype
)
def
triangles
(
self
):
"""Returns the list of boundary triangles."""
if
self
.
_dim
==
2
:
return
np
.
ndarray
((
0
),
self
.
_triangletype
)
return
self
.
_get_array
(
"Triangle"
,
self
.
_triangletype
)
def
contact_forces
(
self
):
"""Returns the contact forces on each grain."""
return
(
self
.
_get_matrix
(
"ContactForces"
,
self
.
_dim
))
def
get_boundary_impulsion
(
self
,
tag
=
"default"
)
:
def
get_boundary_forces
(
self
,
tag
=
"default"
)
:
"""Returns the net force acting on a boundary because of the contacts.
Keyword arguments:
tag -- Name of the boundary
"""
self
.
_lib
.
particleProblemGetTagId
.
restype
=
c_size_t
tag
=
self
.
_lib
.
particleProblemGetTagId
(
self
.
_p
,
tag
.
encode
())
def
compute_fn_ft
(
contact_name
,
objects
,
tag
)
:
...
...
@@ -219,10 +233,12 @@ class ParticleProblem :
f
+=
compute_fn_ft
(
"particle_triangle"
,
self
.
triangles
(),
tag
)
return
f
def
disks
(
self
)
:
return
self
.
_get_array
(
"Disk"
,
self
.
_disktype
)
def
set_boundary_velocity
(
self
,
tag
,
v
)
:
"""Sets the velocity of a boundary to a given value.
Keyword arguments:
tag -- Name of the boundary
v -- Velocity vector to be given to the boundary
"""
self
.
_lib
.
particleProblemGetTagId
.
restype
=
c_size_t
tag
=
self
.
_lib
.
particleProblemGetTagId
(
self
.
_p
,
tag
.
encode
())
disks
=
self
.
disks
()
...
...
@@ -238,7 +254,7 @@ class ParticleProblem :
Only in 2D
Keyword arguments:
tag -- Name of the boundary
vt --
v
elocity in the tangent direction
vt --
V
elocity in the tangent direction
"""
assert
self
.
_friction_enabled
assert
self
.
_dim
==
2
...
...
@@ -254,27 +270,29 @@ class ParticleProblem :
""" Sets the friction coefficient between two materials.
Keyword arguments:
mu --
v
alue of the friction coefficient
mat0 --
f
irst material
mat1 --
s
econd material
mu --
V
alue of the friction coefficient
mat0 --
F
irst material
mat1 --
S
econd material
"""
assert
self
.
_friction_enabled
self
.
_lib
.
particleProblemSetFrictionCoefficient
(
self
.
_p
,
c_double
(
mu
),
mat0
.
encode
(),
mat1
.
encode
())
def
iterate
(
self
,
dt
,
forces
,
tol
=
1e-8
,
force_motion
=
0
)
:
"""Compute iteratively the set of velocities that respect the non-interpenetration constrain.
"""Compute
s
iteratively the set of velocities that respect the non-interpenetration constrain.
Returns 1 if the computation converged, 0 otherwise.
Keyword arguments:
dt -- numerical time step
forces -- list of vectors containing the forces applied on the particles
tol -- optional argument defining the interpenetration tolerance to stop the NLGS iterations of the NSCD
dt -- Numerical time step
forces -- List of vectors containing the forces applied on the particles
tol -- Optional argument defining the interpenetration tolerance to stop the NLGS iterations of the NSCD
force_motion -- Optional argument allowing to move the grains if convergence has not been reached when set to 1
"""
self
.
_lib
.
particleProblemIterate
.
restype
=
c_int
return
self
.
_lib
.
particleProblemIterate
(
self
.
_p
,
c_double
(
np
.
min
(
self
.
r
())
if
self
.
r
().
size
!=
0
else
tol
),
c_double
(
dt
),
c_double
(
tol
),
c_int
(
-
1
),
c_int
(
force_motion
),
_np2c
(
forces
))
def
get_contacts
(
self
,
contact_type
)
:
"""Gives the contact forces. Attention : during the resolution fo the contacts,
the considered quantities are impulses, while forces are written and read."""
"""Gives the contact forces. Warning : during the resolution fo the contacts,
the considered quantities are impulses, while forces are written and read.
"""
ctype
=
{
"particle_particle"
:
0
,
"particle_disk"
:
1
,
"particle_segment"
:
2
,
"particle_triangle"
:
3
}[
contact_type
]
n
=
self
.
_lib
.
particleProblemContactN
(
self
.
_p
,
c_int
(
ctype
))
basis
=
np
.
ndarray
((
n
,
self
.
_dim
*
self
.
_dim
),
dtype
=
np
.
float64
,
order
=
"C"
)
...
...
@@ -286,14 +304,23 @@ class ParticleProblem :
return
o
,
v
,
basis
def
computeStressTensor
(
self
,
nodes
,
radius
):
"""Compute the stres tensor of the granular material"""
"""Computes the stress tensor of the granular material
Keyword arguments:
nodes -- Array of nodes at which to compute the stress tensor
r -- Radius within which the contact forces will be averaged
"""
n_nodes
=
len
(
nodes
[:,
0
])
s
=
np
.
ndarray
((
n_nodes
,
self
.
_dim
**
2
))
self
.
_lib
.
particleProblemComputeStressTensor
(
self
.
_p
,
_np2c
(
nodes
[:,:
self
.
_dim
]),
c_double
(
radius
),
c_int
(
n_nodes
),
_np2c
(
s
))
return
s
def
write_vtk
(
self
,
odir
,
i
,
t
)
:
"""Write output files for post-visualization."""
"""Writes output files for post-visualization.
Keyword arguments:
odir -- Directory in which to write the file
i -- Number of the fiel to write
t -- Time at which the simulation is
"""
v
=
self
.
velocity
()
if
self
.
_rotation_enabled
:
...
...
@@ -342,7 +369,12 @@ class ParticleProblem :
VTK
.
write_multipart
(
odir
+
"/particle_problem"
,[
"particles_%05d.vtp"
%
i
,
"boundaries_%05d.vtu"
%
i
,
"contacts_%05d.vtp"
%
i
],
i
,
t
)
def
read_vtk
(
self
,
dirname
,
i
,
contact_factor
=
1
)
:
"""Read an existing output file to set the particle data."""
"""Reads an existing output file to set the particle data.
Keyword arguments:
dirname -- Path to the directory of the file to read
i -- Number of the file to read
contact_factor -- Factor that determines how to take the read contacts into account
"""
self
.
_lib
.
particleProblemClearBoundaries
(
self
.
_p
)
x
,
_
,
d
,
cdata
,
fdata
=
VTK
.
read
(
dirname
+
"/particles_%05d.vtp"
%
i
)
matnames
=
VTK
.
string_array_decode
(
fdata
[
"MaterialNames"
])
...
...
@@ -405,41 +437,67 @@ class ParticleProblem :
self
.
_lib
.
particleProblemSetContacts
(
self
.
_p
,
c_int
(
t
.
shape
[
0
]),
_np2c
(
o
,
np
.
uint64
),
_np2c
(
v
),
_np2c
(
basis
),
_np2c
(
t
,
np
.
int32
))
def
add_boundary_disk
(
self
,
x0
,
r
,
tag
,
material
=
"default"
)
:
"""Adds a boundary disk.
Keyword arguments:
x0 -- Tuple of the coordinates of the centre
r -- Disk radius
tag -- Disk tag
material -- Disk material
"""
self
.
_lib
.
particleProblemAddBoundaryDisk
.
restype
=
c_size_t
return
self
.
_lib
.
particleProblemAddBoundaryDisk
(
self
.
_p
,
self
.
_coord_type
(
*
x0
),
c_double
(
r
),
tag
.
encode
(),
material
.
encode
())
def
add_boundary_segment
(
self
,
x0
,
x1
,
tag
,
material
=
"default"
)
:
"""Adds a boundary segment.
Keyword arguments:
x0 -- Tuple of the coordinates of the first endpoint
x1 -- Tuple of the coordinates of the second endpoint
tag -- Segment tag
material -- Segment material
"""
self
.
_lib
.
particleProblemAddBoundarySegment
.
restype
=
c_size_t
return
self
.
_lib
.
particleProblemAddBoundarySegment
(
self
.
_p
,
self
.
_coord_type
(
*
x0
),
self
.
_coord_type
(
*
x1
),
tag
.
encode
(),
material
.
encode
())
def
add_boundary_triangle
(
self
,
x0
,
x1
,
x2
,
tag
,
material
=
"default"
)
:
"""Adds a boundary triangle.
Keyword arguments:
x0 -- Tuple of the coordinates of the first summit
x1 -- Tuple of the coordinates of the second summit
x2 -- Tuple of the coordinates of the third summit
tag -- Triangle tag
material -- Triangle material
"""
if
self
.
_dim
!=
3
:
raise
NameError
(
"Triangle boundaries only available in 3D"
)
self
.
_lib
.
particleProblemAddBoundaryTriangle
.
restype
=
c_size_t
return
self
.
_lib
.
particleProblemAddBoundaryTriangle
(
self
.
_p
,
self
.
_coord_type
(
*
x0
),
self
.
_coord_type
(
*
x1
),
self
.
_coord_type
(
*
x2
),
tag
.
encode
(),
material
.
encode
())
def
add_particle
(
self
,
x
,
r
,
m
,
tag
=
"default"
,
forced
=
0
):
"""Add a particle in the particle container.
"""Add
s
a particle in the particle container.
Keyword arguments:
x -- tuple to set the centre particle position
r -- particle radius
m -- particle mass
tag -- particle tag
x -- Tuple to set the centre particle position
r -- Particle radius
m -- Particle mass
tag -- Particle material
forced -- Flag indicating whether the particle is forced or not
"""
self
.
_lib
.
particleProblemAddParticle
(
self
.
_p
,
self
.
_coord_type
(
*
x
),
c_double
(
r
),
c_double
(
m
),
tag
.
encode
(),
c_int
(
forced
))
def
add_particles
(
self
,
x
,
r
,
m
,
v
=
None
,
omega
=
None
,
tag
=
"default"
,
forced
=
None
,
contact_forces
=
None
):
"""Add particles in the particle container.
"""Add
s
particles in the particle container.
Keyword arguments:
x --
a
rray with the particles's centers position
r --
a
rray with the particles's radii
m --
a
rray with the particles's masses
v --
a
rray with the particles's velocities
omega --
a
rray with the particles's angular velocities
tag --
p
article
tag
forced --
a
rray of flags
for
particle forc
ing
contact_forces --
a
rray of contact forces between the particles
x --
A
rray with the particles's centers position
r --
A
rray with the particles's radii
m --
A
rray with the particles's masses
v --
A
rray with the particles's velocities
omega --
A
rray with the particles's angular velocities
tag --
P
article
material
forced --
A
rray of flags
indicating whether the
particle
s are
forc
ed or not
contact_forces --
A
rray of contact forces between the particles
"""
n_particles
=
len
(
m
)
tags
=
np
.
ones
(
n_particles
)
*
self
.
get_tag_id
(
tag
)
...
...
@@ -454,41 +512,41 @@ class ParticleProblem :
self
.
_lib
.
particleProblemAddParticles
(
self
.
_p
,
c_int
(
n_particles
),
_np2c
(
x
),
_np2c
(
m
),
_np2c
(
r
),
_np2c
(
v
),
_np2c
(
omega
),
_np2c
(
tags
,
dtype
=
np
.
int32
),
_np2c
(
forced
,
dtype
=
np
.
int32
),
_np2c
(
contact_forces
))
def
set_use_queue
(
self
,
use_queue
=
1
):
"""Enable the use of the queue if 1 and disables it if 0."""
"""Enable
s
the use of the queue if 1 and disables it if 0."""
self
.
_lib
.
particleProblemSetUseQueue
(
self
.
_p
,
c_int
(
use_queue
))
def
get_tagnames
(
self
):
"""Return the names of the boundaries """
"""Return
s
the names of the boundaries """
self
.
_lib
.
particleProblemGetTagName
.
restype
=
c_char_p
ntagname
=
self
.
_lib
.
particleProblemNTag
(
self
.
_p
)
tagnames
=
list
([
self
.
_lib
.
particleProblemGetTagName
(
self
.
_p
,
i
)
for
i
in
range
(
ntagname
)])
return
tagnames
def
clear_boundaries
(
self
):
"""Remove the boundaries."""
"""Remove
s
the boundaries."""
self
.
_lib
.
particleProblemClearBoundaries
(
self
.
_p
)
def
_save_contacts
(
self
):
"""Save the contacts from the current time step."""
"""Save
s
the contacts from the current time step."""
self
.
_lib
.
particleProblemSaveContacts
(
self
.
_p
)
def
_restore_contacts
(
self
):
"""Restore the saved contacts from the previous time step."""
"""Restore
s
the saved contacts from the previous time step."""
self
.
_lib
.
particleProblemRestoreContacts
(
self
.
_p
)
def
remove_particles_flag
(
self
,
flag
)
:
"""Remove particles based on given flag."""
"""Remove
s
particles based on given flag."""
if
flag
.
shape
!=
(
self
.
n_particles
(),)
:
raise
NameError
(
"size of flag array should be the number of particles"
)
self
.
_lib
.
particleProblemRemoveParticles
(
self
.
_p
,
_np2c
(
flag
,
np
.
int32
))
def
load_msh_boundaries
(
self
,
filename
,
tags
,
shift
=
None
,
material
=
"default"
)
:
"""Load the numerical domain and set the physical boundaries the particles cannot go through.
"""Load
s
the numerical domain and set the physical boundaries the particles cannot go through.
Keyword arguments:
filename --
n
ame of the msh file
tags --
t
ags of physical boundary defined in the msh file
shift --
o
ptional argument to shift the numerical domain
material --
m
aterial of the boundary
filename --
N
ame of the msh file
tags --
T
ags of physical boundary defined in the msh file
shift --
O
ptional argument to shift the numerical domain
material --
M
aterial of the boundary
"""
if
shift
==
None
:
shift
=
[
0
]
*
self
.
_dim
...
...
Write
Preview
Markdown
is supported
0%
Try again
or
attach a new file
.
Attach a file
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment