- Python Interface for flexible and fast input scripting
- Pyck Support for pre-preocessing
- 3D Material Point Method (MPM) Simulation
- Explicit/Leap Frog time integration
- OpenMP Parallel Processing
- VTP/VTI/CSV output formats (some require additional libraries - see below), can be viewed in ParaView
- Current Physics Supported
- Getting the Code
- Running a Simulation
- Code Features
- Creating Visualisations in Houdini
- Graphyt API Documentation
- Solid Mechanics: small strain theory
- Plasticity: Drucker-Prager, Perfect
- Damage: Grady-Kipp
- Fluid Mechanics: Newtonian fluids
- Fluid-Solid Interaction
- Multi-body Contact
- Rigid Bodies (moving & fixed)
- Thermal conductivity
Binaries can be downloaded here:
If you are interested in developing the code, please contact us.
- Make an output directory "output" in the same folder as the python scripts to run
- Create a python script of the model
- Type "python nameOfScript.py" in the terminal
- Interpolation is based on simple GIMP implementation
- Boundary Conditions can be applied to particles (vel, stress, force(accel)) and nodes (vel, force). Currently only one BC type per particle/node
- Domain is periodic by default, for finite simulations, nodal BCs can be applied within 2 cells of the domain extent defined in the scripts
- Output data using VTI file format
- Open files in Paraview: - Create a 'threshold' - Select filter: 'Cell to Point Data' - Create a contour - Sava Data as ".stl" - Import stl's into Houdini
Python modules are required for a graphyt script. Aside from the graphyt module, for more complex geometries, pyck (detailed below) may be used, for visualization, the VTKWriter module (also detailed below) should be included.
import graphyt
import pyck # if using pyck for geometry
import VTKWriter # for writing output filesTo initialize the background grid within which MPM simulations occur, the domain and resolution need to be defined:
L = [Lx,Ly,Lz] # Lx,y,z are the lengths of the grid
in each dimension in meters
cellsize = dx # This is the distance between the nodes
(only uniform resolution supported)
psep = dx / 2 # This is the distance between particles
(optimal value is dx/2)
The nodes and material points' properties are controlled and contained in the graphyt.Nodes and graphyt.MaterialPoints python object. They require some input values during creation as shown below.
nodes = graphyt.Nodes([Lx,Ly,Lz], cellsize)
matpoints = graphyt.MaterialPoints(nodes)These objects are used for the majority of the rest of the simulation process.
A number of material models are available in graphyt. To define a material, the format follows the same pattern:
# Create the material object and set the material-type
material_1 = graphyt.MaterialModels()
material_1.setModel(1)
# 1 - Elastic
# 2 - Newtonian Fluid
# 3 - Elastic + Perfect Plastic
# 4 - Elastic + Drucker Prager Plasticity# Elastic
material_1.setBulkMod(val)
material_1.setShearMod(val)# Newtonian Fluid
material_1.setBulkMod(val)
material_1.setViscosity(val)# Elastic + Perfect Plastic
material_1.setBulkMod(val)
material_1.setShearMod(val)
material_1.setYieldStress(val)# Elastic + Drucker Prager
material_1.setBulkMod(val)
material_1.setShearMod(val)
material_1.setCohesion(val)
material_1.setCriticalStrain(val)
material_1.setFrictionAngle(val)
material_1.setDilatancyAngle(val)Damage models can be added to both elastic and elastic-plastic solids.
# Grady-Kipp
material_1.setIsDamage(boolean)
material_1.setDamageModel(val)
material_1.setCrackSpeed(val)
material_1.setDamageM(val)
material_1.setDamageK(val)# Set Material Density
material_1.setDensity(val)
# Add the material to a materials array, more than one material can be defined
materials = [material_1]##PYCK *required objectID, materialID
*required objectID, materialID
Initializing particle properties. These commands will set particle properties for the time step they are applied. For example, if set at the beginning of the simulation they are initial conditions. Once set these values are free to change with the dynamics of the system. For persistent values the user can either set these values for each timestep or it may be more appropriate to use the Boundary Conditions (see the Particle Boundary Conditions section).
setVel(particle_index, vx, vy, vz)
#For the given particle_index, that particle's velocity-vector is set to: [vx,vy,vz]
setDamage(particle_index, damage_value)
#For the given particle_index, that particle's damage value is set to: damage_value
setStress(particle_index,sxx,syy,szz,sxy,syz,sxz)
#For the given particle_index, that particle's stress-tensor is set to:
#[sxx, sxy, sxz
# sxy, syy, syz
# sxz, syz, szz]
setTemp(particle_index, Temperature)
#For the given particle_index, that particle's temperature is set to: Temperature
holdParticle(particle_index)
#For the given particle_index, this particle's position is not updated, it stays fixed in place. Note for rigid particles, see the Materials section.
#### Setting Particle Boundary Conditions*required
node.setBC(nodeID,graphyt.BCTypes.grid_type,[vals])
*required
is3D,tmax, graphyt.Parameters(max,is3D)
*required
set arrays, VTPWriter.VTPWriter(L,cellsize,numParticles)
*required
graphyt.Solver(nodes, matpoints, parameters, materials) solver.iterate_CPU(t, step,parameters)
# Import the necessary modules
import graphyt
import pyck
import VTKWriter
import math
## GEOMETRY ##
# Set the resolution and the pack
cellsize = 0.02
psep = cellsize/2.0
# Background Grid
L = [1.0,1.0,0.0]
nodes = graphyt.Nodes(L, cellsize)
print("+++++ Grid INFO+++++ No. Nodes:", nodes.numNodes)
# Material Points
matpoints = graphyt.MaterialPoints(nodes)
matpoints.psep = psep
# PYCK
# create materials
water = graphyt.MaterialModels()
water.setBulkMod(1e6)
water.setShearMod(0.0)
water.setDensity(1000)
water.setModel(2)
materials = [water]
# create packers and pack
cubic = pyck.CubicPacker(L,psep,[0,0,0])
pack = pyck.StructuredPack(cubic)
# create shapes
liquidLowerLeft = [2.75*cellsize,2.75*cellsize,0.0]
liquidUpperRight = [0.5*L[0],0.75*L[1]-2*cellsize,0.0]
liquid = pyck.Cuboid(1,liquidLowerLeft,liquidUpperRight)
# Add shapes to pack
pack.AddShape(liquid)
pack.Process()
model = pyck.Model(pack)
# Set the objectID's for all objects
objID = model.CreateIntField("objectID",1)
model.SetIntField(objID,1,0)
# Set the material ID for all objects
matID = model.CreateIntField("material",1)
model.SetIntField(matID,1,0)
# Add the geometry from pyck to the matpoints
matIDfield = model.GetIntField(matID)
objIDfield = model.GetIntField(objID)
numParticles = pack.GetNumParticles()
pyckPos = pack.GetPositions()
matpoints.addPyckGeom(pyckPos,numParticles)
matpoints.addPyckIntField(matIDfield)
matpoints.addPyckIntField(objIDfield)
print("+++++ MAT POINT INFO+++++ No. Particles:", matpoints.numParticles)
# /////////-----------BOUNDARY CONDITIONS--------------///////
# //// MATERIAL POINT BOUNDARY CONDITIONS /////
# //// NODAL BOUNDARY CONDITIONS /////
#// Wall boundaries for Grid Nodes
wallthickness = 2.0*nodes.cellsize
for n in range(nodes.numNodes):
posX = nodes.getPosX(n)
posY = nodes.getPosY(n)
posZ = nodes.getPosZ(n)
if (((posX>=L[0]-wallthickness) or (posX <= wallthickness)) or ((posY>=L[1]-wallthickness) or (posY <= wallthickness))):
nodes.setBC(n,graphyt.BCTypes.grid_vx_vy,[0,0])
print("+++++ MODEL INFORMATION+++++ Boundary Conditions Applied")
# /////////// Parameters Set /////////////// #
tmax = 5 # //(time in seconds)
is3D = False
parameters = graphyt.Parameters(tmax,is3D)
parameters.setGravity(0, -10, 0)
parameters.setCFL(0.0001)
parameters.setDampingCoef(0.5)
print("+++++ MODEL INFORMATION+++++ Parameters Set")
# // //+++++++++++++++++++++++++++++++++++++++++++++++//
# // //+++++++++++++ MODEL RUN +++++++++++++++++++++++//
# // //+++++++++++++++++++++++++++++++++++++++++++++++//
#Now Mat points and nodes are set and ready to go, create a solver to run an MPM analysis on them
print("+++++ SOLVER INFORMATION+++++ Starting Solver")
solver = graphyt.Solver(nodes, matpoints, parameters, materials)
matpntPos = matpoints.getPosArr()
matpntMat = matpoints.getMatArr()
matpntMass = matpoints.getMassArr()
matpntVel = matpoints.getVelArr()
matpntStress = matpoints.getStressArr()
Le = [1+int(L[0]/cellsize), 1+int(L[1]/cellsize), 1+int(L[2]/cellsize)]
vtp = VTKWriter.VTPWriter(Le,cellsize,numParticles)
vtp.AddPositions("Position", matpntPos,3, 2)
vtp.AddArray("Material", matpntMat,1, 2)
vtp.AddArray("Mass", matpntMass,1, 2)
vtp.AddArray("Velocity",matpntVel,3, 2)
vtp.AddArray("Stress",matpntStress,6, 2)
# Iterate over nsteps number of timesteps
for step in range(int(solver.tmax/solver.dt)): # ; t<solver.tmax; t+=solver.dt):
t = solver.dt * step
solver.iterate_CPU(t, step, parameters)
vtpfname='output/dambreak2D'+str(step)+'.vtp'
if step % 50000 == 0:
vtp.Write(vtpfname)
print("Total Simulation Duration (s)")