Functions in FEM-DEM

Hello

> using FEM-DEM

please be more specific, there are several approaches / projects on this. ​

> if clump coupling is also supported.

see above

> the functions in YADE model are not recognized

please be more specific, e.g. if/what error you get etc.

> Do functions literally work in FEM-DEM coupling?

yes, they should

> If yes do I need to adjust something scripts other than the YADE mode?

what does "other mode" mean?

Cheers
Jan

Thanks so much, Jan

> please be more specific, there are several approaches / projects on this. ​
My problem falls in the soil-tool interaction category. The original DEM geometry is constructed outside this script and is then imported here. That DEM model is mostly composed of clumps of spheres of different shapes and sizes but there exist some standalone spheres as well.

Soil Geometry
I am interested to implement a volumetric coupling at the bottom of a 2.5 mm DEM model with their original clumps but am not sure if FEM-DEM coupling works for clumps. As a result, I considered a separate transition zone at the bottom of this 2.5 mm which is purely constructed from spherical particles. In a separate script, the clumps are first deposited on this transition and at the end of the deposition, the DEM geometry which now has two parts (spheres for the transition and clumps for the upper 2.5mm layer) are exported. That is why I have two ymport in my script.

This DEM is then coupled from the bottom to a coarse FEM mesh with constrained BCs (all 6 DOFs) at the bottom and sides.
To constrain the lateral movement of the DEM, I have used a facetBox (wallMask=51) which lacks the top and bottom boundaries. Therefore, the model owns two types of boundaries: nodes constraints of the FEM mesh, and the facetBox for the upper layer DEM.

I hope that I have been able to explain the geometry.

Motions
Upon the generation of the model geometry (soil) and application of boundary conditions, I try to model the interaction between a face wheel over the surface of the DEM terrain:
Using a combinedKinematicsEngine:
1- The facet wheel is first penetrated vertically (in the Y direction) up to a certain depth (1mm). In this phase, the Translation component (of the combinedKinematicsEngine) controls the movement of the wheel. The motion of the wheel is continuously monitored using the zero point of the Rotation Engine component which in the penetration phase has zero angular velocity but is NOT dead. I mean that the difference between the updating coordination of the zero point in direction [1] and the wheel radius, determines when to stop the penetration and start the 2 phase (wheel rotation).
2- Upon the arrival of the wheel to the desired penetration depth, the Translation component (of the combinedKinematicsEngine) is deactivated and the Rotation Engine component imposes the rotational movement of the wheel using an angular velocity which is calculated from the penetration velocity by means of a SLIDING RATIO.

As you can see, I need functions for enabling/disabling different motions of the wheel.

THE PREVIOUS SCRIPT HAD SOME ISSUES. THIS ONE IS THE UPDATED ONE:
THE MOST ANNOYING ERROR WHILE RUNNING IS THAT THE WheelPenetration() FUNCTION HAS NOT BEEN DEFINED. Therefore the PyRunner doesn't work. As a consequence, the next functions are not executed as well.

>what does "other mode" mean?

I meant model not mode. I mean if we add some functions in YADE script (DEM), do we need to change something in the coupling or FEM scripts?

BTW, functions () are not callable and this is my first and foremost issue.
Note that a pure DEM model with the same configuration works as a charm.

Thanks so much
######################################################################
######################################################################
import math
from math import sqrt
from libyade import yade
from yade import *
from yade import pack, plot
import numpy
from pyquaternion import Quaternion
from yade import pack, export, Vector3
from yade import ymport
import sys
import os
import os.path
import IPython
############################################
### INPUT PARAMETERS ###
############################################
## Wheel PROPERTIES
widthWheel=0.00125
Wheel_R=0.0025
Wheel_AdditionalNormalForce=-19.6

##WHELL MOVEMENT
slidingratio=0.3
velocity=0.01 #m/s THAT SHOULD BE 0.01 m/s
angularVelocity=velocity/((1-slidingratio)*Wheel_R) #RAD/S
gravityAcc=-9.81

deposFricDegree = 28.5 # INITIAL CONTACT FRICTION FOR SAMPLE PREPARATION
normalDamp=0.7 # NUMERICAL DAMPING
shearDamp=0.7
youngSoil=0.7e8# CONTACT STIFFNESS FOR SOIL
youngContainer=210e9 # CONTACT STIFFNESS FOR CONTAINER
poissonSoil=0.3 # POISSION'S RATIO FOR SOIL
poissionContainer=0.25 # POISSION'S RATIO FOR CONTAINER
densSoil=2650 # DENSITY FOR SOIL
densContainer=7850 # DENSITY FOR CONTAINER
numDamp=0.4
binSize=0.4*Wheel_R
activeDistance=0.01
iniDistanceWheelfromBoundary=0.004
differenceWCenterActiveEnd=activeDistance - iniDistanceWheelfromBoundary
initialPeneterationofWheel=0.001
iniWheelDisfromMaxY=0.0001

############################################
### DEFINING VARIABLES AND MATERIALS ###
############################################
SoilId=FrictMat(young=youngSoil,poisson=poissonSoil,frictionAngle=radians(deposFricDegree),density=densSoil)
O.materials.append(SoilId)

ContainerId=FrictMat(young=youngContainer,poisson=poissionContainer,frictionAngle=radians(0),density=densContainer)
O.materials.append(ContainerId)

###################################
##### CREATING GEOMETRIES #####
###################################
## CREATE WALLS AROUND THE PACKING
radius=8.e-5
XcenterofContainer = 0.0025
YcenterofContainer = 0.0054383384
minY=0.00241406
ZcenterofContainer = 10.0e-3
Container=yade.geom.facetBox((XcenterofContainer,YcenterofContainer,ZcenterofContainer), (XcenterofContainer,(YcenterofContainer-minY),ZcenterofContainer),wallMask=51,material=EOId)
O.bodies.append(Container)
## IMPORTING THE ALREADY COMPACTED BIN
ymport.textClumps("/home/ngoudarzi/Desktop/dem-fem-coupling-master/examples/vol5/DEMPClumps.txt",shift=Vector3(0,0,0),material=SoilId)
O.bodies.append(ymport.text("/home/ngoudarzi/Desktop/dem-fem-coupling-master/examples/vol5/DEMSpheres.txt",shift=Vector3(0,0,0),material=SoilId,color=Vector3(0.6,0.6,0.6)))
minX=min([b.state.pos[0]-b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])
maxX=max([b.state.pos[0]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])
minY=min([b.state.pos[1]-b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])
maxY=max([b.state.pos[1]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])
minZ=min([b.state.pos[2]-b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])
maxZ=max([b.state.pos[2]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])
print ("minX:",minX,"maxX:",maxX,"minY:",minY,"maxY:",maxY,"minZ:",minZ,"maxZ:",maxZ)

XcenterofWheel = widthWheel/2+minX
YcenterofWheel = Wheel_R+maxY+iniWheelDisfromMaxY
ZcenterofWheel = iniDistanceWheelfromBoundary+minZ

Wheel1 = yade.geom.facetCylinder((XcenterofWheel,YcenterofWheel,ZcenterofWheel),radius=Wheel_R,height=widthWheel,orientation=Quaternion((0,1,0),-pi/2),wallMask=7,segmentsNumber=200,material=ContainerId,wire=True,color=Vector3(255,255,255))
O.bodies.append(Wheel1)
idsr = [w.id for w in Wheel1]
facets = [b for b in O.bodies if isinstance(b.shape,Facet)] # list of facets in simulation
for b in facets:
  O.forces.setPermF(b.id,(0,Wheel_AdditionalNormalForce/800,0))

############################
### DEFINING ENGINES ###
############################
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_Facet_Aabb()], label="collider"),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_MindlinPhys(betan=normalDamp,betas=shearDamp,label='ContactModel')],
[Law2_ScGeom_MindlinPhys_Mindlin(label='Mindlin')]
  ),
  ## WE WILL USE THE GLOBAL STIFFNESS OF EACH BODY TO DETERMINE AN OPTIMAL TIMESTEP (SEE HTTPS://YADE-DEM.ORG/W/IMAGES/1/1B/CHAREYRE&VILLARD2005_LICENSED.PDF)
  #GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.25),
  CombinedKinematicEngine(ids=idsr,label='combEngine') + TranslationEngine(translationAxis=(0,-1,0),velocity=velocity) +\
  RotationEngine(rotationAxis=(1,0,0), angularVelocity=0, rotateAroundZero=True, zeroPoint=(XcenterofWheel,YcenterofWheel,ZcenterofWheel)),
  NewtonIntegrator(damping=numDamp,gravity=(0,gravityAcc,0)),
  PyRunner(iterPeriod=1, command="WheelPenetration()" ,label='checker'), #######THIS IS NOT CALLED
  PyRunner(iterPeriod=1000,command='history()',label='recorder'),
  ]
O.dt = 1e-7
O.step()
# Get TranslationEngine and RotationEngine from CombinedKinematicEngine
transEngine, rotEngine = combEngine.comb[0], combEngine.comb[1]

def WheelPenetration():
  transEngine.velocity = 7*velocity
  rotEngine.zeroPoint += Vector3(0,-1,0)*transEngine.velocity*O.dt
  if (rotEngine.zeroPoint[1]-Wheel_R)<=(maxY-initialPeneterationofWheel):
    checker.command='WheelRolling()'

def WheelRolling():
  transEngine.translationAxis=(0,0,1)
  transEngine.velocity = velocity
  rotEngine.angularVelocity = angularVelocity
  rotEngine.zeroPoint += Vector3(0,0,1)*velocity*O.dt
  if rotEngine.zeroPoint[2]>=maxZ-1.05*Wheel_R:
    O.pause()

def history():
  global Fx,Fy,Fz,Dx,Dy,Dz
  Fx=0
  Fy=0
  Fz=0
  for b in facets:
    Fx+=O.forces.f(b.id,sync=True)[0]
    Fy+=O.forces.f(b.id,sync=True)[1]
    Fz+=O.forces.f(b.id,sync=True)[2]
  Dx=rotEngine.zeroPoint[0]
  Dy=rotEngine.zeroPoint[1]
  Dz=rotEngine.zeroPoint[2]
  yade.plot.addData({'i':O.iter,'Fx':Fx,'Fy':Fy,'Fz':Fz,'Dx':Dx,'Dy':Dy,'Dz':Dz,})
  ## In that case we can still save the data to a text file at the end of the simulation, with:
  plot.saveDataTxt('./Wheel_OutputData_5by6by20_n057_Relaxed_FreeSurface.txt')
def vtkExport(i):
  from yade import export
  name = '/tmp/vol1_yade'
  export.VTKExporter(name,i).exportSpheres()
  export.VTKExporter(name,i).exportFacets()

Hi Nima,

thanks for more information.

> FEM-DEM

what does this exactly mean? [1]? [2]? something other 3rd party? Something yours?

What is the difference in running pure DEM and DEM-FEM?

> FUNCTION HAS NOT BEEN DEFINED

Ok, it might be some scope problem e.g. if you do not use Yade executable but pure python. In this case, "saving" the functions as "global" might help, using "builtin(s)" module.
What Python version do you use?

### python 2
import __builtin__
...
def someFunction(...):
    ...
__builtin__.someFunction = someFunction
###

### python 3
import builtins
...
def someFunction(...):
    ...
builtins.someFunction = someFunction
###

The same approach can be used for other functions or variables.
Putting them into "builtin", they are "visible" also outside their scope.

Have a try. Be careful not to overwrite some important names (like 'list' or 'dict'), but seeing the functions names it should not be the case here.

Cheers
Jan

[1] https://github.com/stranskyjan/dem-fem-coupling
[2] https://yade-dem.org/doc/FEMxDEM.html

Hi and thanks Jan,

Indeed, I'm using [1]. However, I didn't include other dependencies including oofem file in my previous posts.
The pure DEM (with moving wheel as facet) works by YADE executable. And yes, my FEM-DEM uses Python purely which might be the source of the issue.

My python is 2.7 and my Ubuntu 18.04 since I kept everything old for being able to install FEM-DEM after struggling a lot for installation on Ubuntu 20.04 with no avails. This is my other question: Can FEM-DEM be installed on Ubunto 20.04 with newer versions of Python?

Thanks for your proposed solution. I'll give it a try. My question here is: Can I use __builtin__.someFunction inside engines when I implement the Pyrunner. For example, does this work?
import __builtin__
.
.
.
 PyRunner(iterPeriod=1, command="__builtin__.WheelPenetration()" ,label='checker')

or do I need to declare it somewhere else?

Thanks so much

> This is my other question: Can FEM-DEM be installed on Ubunto 20.04 with newer versions of Python?

No.
IMO Ubuntu itself should not be a problem, but Python needs to be version 2 (as the versions of Yade and Oofem used uses only Python 2)
You can open an issue on the project sites if you have some terminal outputs e.g.

theoretically it is possible, but I do not have time for the update :-)
It would need new version of both Yade and Oofem, selecting a specific versions, create the patches to make it work together, .......

> Can I use __builtin__.someFunction inside engines when I implement the Pyrunner. For example, does this work?
> import __builtin__
> PyRunner(iterPeriod=1, command="__builtin__.WheelPenetration()" ,label='checker')
> or do I need to declare it somewhere else?

Yes. But putting it into __builtin__, also this should work (not tested) without specifying __builtin__ explicitly:
PyRunner(iterPeriod=1, command="WheelPenetration()" ,label='checker')

Cheers
Jan