yade-dev team mailing list archive

[Vasileios Angelidakis <b7063391@xxxxxxxxxxxxxxx>]

Hi,

I have started working on the "PotentialBlock" code in YADE for the generation of polyhedra using the Potential Particles approach. I want to use these particles in a periodic cell, but it seems the PotentialBlock class is not compatible with periodic boundaries. Would be grateful to get any advice on whether this is the case and where I should focus to implement it myself? (FYI I am still a rookie in C++ development).

In the following lines I paste a minimal working script to demonstrate a potential block falling through the periodic cell.
Visualisation of the simulation is available only in a VTK format (not in qt).

Cheers,

Vasileios



Vasileios Angelidakis

Post-Graduate Researcher in Geotechnical Engineering

School of Engineering, Newcastle University



The script:

from yade import pack
import math

# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# Clear the directory where the VTK results will be written of all files (but not subdirectories).
# The directory is cleared, so files from previous runs do not interfere with files from new runs.
# If the directory does not exist, it is created.

import os, shutil
folder = './vtk_ele'
try:
    os.mkdir(folder[2:])
except:
    for the_file in os.listdir(folder):
        file_path = os.path.join(folder, the_file)
        try:
        if os.path.isfile(file_path):
            os.unlink(file_path)
    #        elif os.path.isdir(file_path): shutil.rmtree(file_path) #uncomment to also delete the subdirectories
        except Exception as e:
        print(e)
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# Engines
O.engines=[
    ForceResetter(),
    InsertionSortCollider([PotentialBlock2AABB()],verletDist=0.01),
    InteractionLoop(
        [Ig2_PB_PB_ScGeom()],
        [Ip2_FrictMat_FrictMat_KnKsPBPhys(kn_i=1e8, ks_i=1e8, Knormal = 1e8, Kshear = 1e8, useFaceProperties=False, calJointLength=False, twoDimension=False, unitWidth2D=1.0, viscousDamping=0.05)],
        [Law2_SCG_KnKsPBPhys_KnKsPBLaw(label='law',neverErase=False)]
    ),
    NewtonIntegrator(damping=0.0,exactAsphericalRot=False,gravity=[0,-10,0]),
    VTKRecorder(fileName=folder+'/vtkPeriodicCell-VTKRecorder-',recorders=['pericell'],iterPeriod=200)
]

# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# Basic dimensions used in the model
distanceToCentre= 0.5
meanSize = 1.0
heightOfBase = 5.0*meanSize

# Material definition
powderDensity = 2000
O.materials.append(FrictMat(young=150e6,poisson=.4,frictionAngle=radians(0.0),density=powderDensity,label='frictionless'))

# Creation of a spheres packing, following a desired PSD
sp=pack.SpherePack()

mn,mx=(Vector3(0,0,0),
       Vector3(4*heightOfBase, 4*heightOfBase, 4*heightOfBase))

sphereRad = sqrt(3.0)*0.5*meanSize
sp.makeCloud(mn,mx,sphereRad,0,10,periodic=True)

# Replacement of the spheres with cuboids
count= 0
rPP=0.05*meanSize
for s in sp:
    b=Body()
    dynamic=True
    wire=False
    color=[0,0,255.0]
    highlight=False
    b.shape=PotentialBlock(
k=0.2, r=0.05*meanSize, R=1.02*sphereRad,
a=[1.0,-1.0,0.0,0.0,0.0,0.0],
b=[0.0,0.0,1.0,-1.0,0.0,0.0],
c=[0.0,0.0,0.0,0.0,1.0,-1.0],
d=[distanceToCentre-rPP,distanceToCentre-rPP,distanceToCentre-rPP,distanceToCentre-rPP,distanceToCentre-rPP,distanceToCentre-rPP], isBoundary=False, color=color, wire=wire, highlight=highlight,
minAabb=Vector3(1.0*sphereRad,1.0*sphereRad,1.0*sphereRad),
maxAabb=Vector3(1.0*sphereRad,1.0*sphereRad,1.0*sphereRad),
maxAabbRotated=Vector3(1.0*sphereRad,1.0*sphereRad,1.0*sphereRad),
minAabbRotated=Vector3(1.0*sphereRad,1.0*sphereRad,1.0*sphereRad),
AabbMinMax=True,fixedNormal=False)

    length=meanSize
    V= 1.0
    geomInert=(2./5.)*powderDensity*V*sphereRad**2
    utils._commonBodySetup(b,V,Vector3(geomInert,geomInert,geomInert), material='frictionless',pos=s[0],  dynamic=True, fixed=False)
    b.state.pos = s[0] #s[0] stores center
    b.state.ori = Quaternion((random.random(),random.random(),random.random()),random.random()) #s[2]
    b.state.mass =V*powderDensity
    O.bodies.append(b)
    b.dynamic = True
    count =count+1

# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# Set up Periodic Boundaries
O.periodic = True
O.cell.setBox(4*heightOfBase, 4*heightOfBase, 4*heightOfBase)

# PotentialBlockVTKRecorder recorder
sampleQuality=50 # Increase to improve visual quality of the particles in Paraview
O.engines=O.engines+[PotentialBlockVTKRecorder(fileName=folder+'/eleTest_', iterPeriod=200, twoDimension=False, sampleX=sampleQuality,  sampleY=sampleQuality,  sampleZ=sampleQuality, maxDimension=0.2, label='PBvtkRecorder')]

# Time step
O.dt = 0.2*sqrt(0.3*O.bodies[0].state.mass/1.0e9)
import yade.timing
O.timingEnabled = True
yade.timing.reset()