Loading the stored the position of plank, but it wrong

Asked by yang yi on 2020-05-29

Dear friend:

It is me again. Sorry. The new problem.

I design five plank(window), and they rotates around a axis. So I storage the position of them, and hope display the rotation process of them, but the reload location is wrong.

# (1) Define the plank:
temp = [
            Vector3(widthHydr * i, wind_y_0, wind_z_0),
            Vector3(widthHydr * (i + 1), wind_y_0, wind_z_0),
            Vector3(widthHydr * (i + 1), wind_y_1, wind_z_1),
            Vector3(widthHydr * i, wind_y_1, wind_z_1)
        ]
        positionWind.append(temp)

    Wind1 = pack.sweptPolylines2gtsSurface([positionWind[0]], capStart=True, capEnd=True)
    Wind2 = pack.sweptPolylines2gtsSurface([positionWind[1]], capStart=True, capEnd=True)
    Wind3 = pack.sweptPolylines2gtsSurface([positionWind[2]], capStart=True, capEnd=True)
    Wind4 = pack.sweptPolylines2gtsSurface([positionWind[3]], capStart=True, capEnd=True)
    Wind5 = pack.sweptPolylines2gtsSurface([positionWind[4]], capStart=True, capEnd=True)

     WindList = IDWind1 + IDWind2 +IDWind3 + IDWind4 + IDWind5

   # (2) Rotation the plank(window) in the O.engines

    RotationEngine(rotationAxis=(1, 0, 0), rotateAroundZero=True, ids=IDWind[0], zeroPoint=positionWind[0][2],
                   label='RotationW1'),
    RotationEngine(rotationAxis=(1, 0, 0), rotateAroundZero=True, ids=IDWind[1], zeroPoint=positionWind[1][2],
                   label='RotationW2'),
    RotationEngine(rotationAxis=(1, 0, 0), rotateAroundZero=True, ids=IDWind[2], zeroPoint=positionWind[2][2],
                   label='RotationW3'),
    RotationEngine(rotationAxis=(1, 0, 0), rotateAroundZero=True, ids=IDWind[3], zeroPoint=positionWind[3][2],
                   label='RotationW4'),
    RotationEngine(rotationAxis=(1, 0, 0), rotateAroundZero=True, ids=IDWind[4], zeroPoint=positionWind[4][2],
                   label='RotationW5'),

# (3) Control the speed and direction of the windows

def WindowsAction(IDWind):
    global WinAction, windPosition
    RotationW = [RotationW1, RotationW2, RotationW3, RotationW4, RotationW5]
    for nW in range(0, numWinds):
        ## action
        Pos_z = sum(O.bodies[i].state.pos[2] for i in IDWind[nW]) / len(IDWind[nW])
        if WinAction[nW] == 0:
            RotationW[nW].angularVelocity = -VelocityWindNegative
        else:
            RotationW[nW].angularVelocity = VelocityWindPositive

        NegtiveStop = (RotationW[nW].angularVelocity > 0) & (Pos_z <= windLowerBoundary)
        PostiveStop = (RotationW[nW].angularVelocity < 0) & (Pos_z >= windUpperBoundary)

        if NegtiveStop: windPosition[nW] = windPositionNegative
        elif PostiveStop: windPosition[nW] = windPositionPositive
        else: windPosition[nW] = 0

        if NegtiveStop | PostiveStop:
            RotationW[nW].angularVelocity = 0

# (4) save the location of windows
def SaveProcessLocation(Episode, iter):
    global WindList
    outputDir = "Output/location/"
    if not os.path.exists(outputDir):
        os.makedirs(outputDir)

    PositonWind = []

    for i in WindList:
        temp = [i, o.bodies[i].state.pos]
        PositonWind.append(temp)

    wind_name = outputDir + 'wind_' + str(Episode) + '_' + str(iter)

    np.save(wind_name, PositonWind)

#(5) display the location of the windows

def Relocation():
    global nEpisode, nIter, index
    path = 'Output/location/'
    wind_name = path + 'wind_' + str(nEpisode) + '_' + str(nIter) + '.npy'

    if nIter >= index[nEpisode, 1]:
        nIter = 0
        nEpisode += 1
    else:
        nIter += 1

    if nEpisode >= len(index):
        O.pause()

    locationWind = np.load(wind_name, allow_pickle=True)

    for i in range(0, len(locationWind)):
        n = locationWind[i][0]
        o.bodies[n].state.pos = locationWind[i][1]
        print(locationWind[i][1])

Problem: the each window constitute of 2 triangle. So if reload the location, the 2 triangles comes off. and the rotation angle does not look right.

Thank you very much

Question information

Language:
English Edit question
Status:
Answered
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
2020-05-29
Last reply:
2020-06-01
Jan Stránský (honzik) said : #1

Hello,

please try to provide a MWE [1], where:
- M stands for minimal. Do you have the same problem with one plank? then please provide a script with just one plank..
- W stands for working. I have tried your code with IndentationError: unexpected indent on line 8

> So I storage the position of them, and hope display the rotation process of them, but the reload location is wrong.
> So if reload the location ... the rotation angle does not look right.

store also the rotation (I have not found it in the script)

> the rotation angle does not look right.

as suggested before, use a numerical evidence rather than visual evidence (e.g. print rotation / vertices before loading and after loading)

cheers
Jan

[1] https://www.yade-dem.org/wiki/Howtoask

yang yi (yangyizaixian) said : #2

To Jan Stránský (honzik) :
Thank you very much. I attach the script as follow. There are two scripts: (1) the simulation to produce the position of 5 windows (2) display the process.

====================================================================================
script (1) Simulation script:

#!/usr/bin/env python
#  encoding: utf-8
from __future__ import print_function
import sys

sys.path.append('~/PycharmProjects/Yade20191229/')
# from yadeImport import *

import math
import numpy as np
import os
from yade.gridpfacet import *
o = Omega()
o.dt = 1e-12
outputDir = 'Output'
checkCounter = 0
numWinds = 5
WinAction = np.zeros((5,1), dtype=int)

## =========== Environment ============##
widthSpace = 6.8 ## width is the distance from the front coal wall to the back
widthHydr = 1.5
lengthSpace = widthHydr * numWinds ## length is the width of 5 windows
highHydr = 3.8
CheckThick = 1
highDummy = 3
colorWind = (0, 1, 1)
angleShield = 50 * 3.1415926 / 180
angleSwingPositive = 15 * 3.1415926 / 180
angleSwingNegtive = 40 * 3.1415926 / 180
lengthShield = 3
lengthTail = 2
windUpperBoundary = 0.9
windLowerBoundary = 0.5019
stateUpperBoundary = highHydr
stateLowerBoundary = windLowerBoundary

positionWind = []
windPosition = np.zeros(5)

shield_y_0 = lengthTail * math.cos(angleShield - angleSwingPositive)
shield_y_1 = shield_y_0 + lengthShield * math.cos(angleShield)
shield_z_0 = highHydr - lengthShield * math.sin(angleShield)
shield_z_1 = highHydr #
wind_y_0 = lengthTail * (math.cos(angleShield - angleSwingPositive) - math.cos(angleShield))
wind_y_1 = shield_y_0
wind_z_0 = highHydr - (lengthShield + lengthTail) * math.sin(angleShield)
wind_z_1 = highHydr - lengthShield * math.sin(angleShield)
HyMat = O.materials.append(FrictMat(frictionAngle=0.1, density=3000, young=2e9))
myGraviaty = -1200
def HydraulicSupport():
    for i in range(0, numWinds):

        temp = [
            Vector3(widthHydr * i, wind_y_0, wind_z_0),
            Vector3(widthHydr * (i + 1), wind_y_0, wind_z_0),
            Vector3(widthHydr * (i + 1), wind_y_1, wind_z_1),
            Vector3(widthHydr * i, wind_y_1, wind_z_1)
        ]
        positionWind.append(temp)

    Wind1 = pack.sweptPolylines2gtsSurface([positionWind[0]], capStart=True, capEnd=True)
    Wind2 = pack.sweptPolylines2gtsSurface([positionWind[1]], capStart=True, capEnd=True)
    Wind3 = pack.sweptPolylines2gtsSurface([positionWind[2]], capStart=True, capEnd=True)
    Wind4 = pack.sweptPolylines2gtsSurface([positionWind[3]], capStart=True, capEnd=True)
    Wind5 = pack.sweptPolylines2gtsSurface([positionWind[4]], capStart=True, capEnd=True)

    IDWind1 = O.bodies.append(pack.gtsSurface2Facets(Wind1, color=colorWind))
    IDWind2 = O.bodies.append(pack.gtsSurface2Facets(Wind2, color=colorWind))
    IDWind3 = O.bodies.append(pack.gtsSurface2Facets(Wind3, color=colorWind))
    IDWind4 = O.bodies.append(pack.gtsSurface2Facets(Wind4, color=colorWind))
    IDWind5 = O.bodies.append(pack.gtsSurface2Facets(Wind5, color=colorWind))

    IDWind = [IDWind1, IDWind2, IDWind3, IDWind4, IDWind5]
    WindList = IDWind1 + IDWind2 + IDWind3 + IDWind4 + IDWind5
    return IDWind, WindList

##---------------------------------------##
def WindowsAction(IDWind):
    global WinAction, windPosition, saveCounter, nEpisode
    RotationW = [RotationW1, RotationW2, RotationW3, RotationW4, RotationW5]

    for nW in range(0, numWinds):
        RotationW[nW].angularVelocity = -0.000001
    SaveProcessLocation(nEpisode, saveCounter)
    saveCounter += 1

def SaveProcessLocation(Episode, iter):
    outputDir = "Output/location/"
    if not os.path.exists(outputDir):
        os.makedirs(outputDir)
    Position = []
    for i in WindList:
        print(o.bodies[i].state.pos)
        temp = [i, o.bodies[i].state.pos]
        Position.append(temp)
    location_name = outputDir + 'wind_' + str(Episode)+ '_' + str(iter)
    np.save(location_name, Position)

def WindowsControl():
    global saveCounter,nEpisode
    nEpisode += 1
    saveCounter = 0

nEpisode = 0
saveCounter = 0
IDWind, WindList = HydraulicSupport()

O.engines = [
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()]),
    InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom(), Ig2_Wall_Sphere_ScGeom()],
        [Ip2_FrictMat_FrictMat_FrictPhys(),
         Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(xSectionWeibullScaleParameter=0.5,
        xSectionWeibullShapeParameter=0.5,
        weibullCutOffMin=0,
        weibullCutOffMax=10)],
        [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=False, recordMoments=False,label='interactionLaw'),
        Law2_ScGeom_FrictPhys_CundallStrack()]
        ),
    GlobalStiffnessTimeStepper(),

    NewtonIntegrator(gravity=(0, 0, myGraviaty), damping=0.5, label='down'),
    RotationEngine(rotationAxis=(1, 0, 0), rotateAroundZero=True, ids=IDWind[0], zeroPoint=positionWind[0][2],
                   label='RotationW1'),
    RotationEngine(rotationAxis=(1, 0, 0), rotateAroundZero=True, ids=IDWind[1], zeroPoint=positionWind[1][2],
                   label='RotationW2'),
    RotationEngine(rotationAxis=(1, 0, 0), rotateAroundZero=True, ids=IDWind[2], zeroPoint=positionWind[2][2],
                   label='RotationW3'),
    RotationEngine(rotationAxis=(1, 0, 0), rotateAroundZero=True, ids=IDWind[3], zeroPoint=positionWind[3][2],
                   label='RotationW4'),
    RotationEngine(rotationAxis=(1, 0, 0), rotateAroundZero=True, ids=IDWind[4], zeroPoint=positionWind[4][2],
                   label='RotationW5'),

    PyRunner(command="WindowsAction(IDWind)", iterPeriod=200000),
    PyRunner(command="WindowsControl()", iterPeriod=2000000),
]

=====================================================================================
script (2) Display the process

#!/usr/bin/env python
#  encoding: utf-8
from __future__ import print_function
import sys
sys.path.append('~/PycharmProjects/Yade20191229/')
# from yadeImport import *
import math
import numpy as np
import os
from yade.gridpfacet import *
o = Omega()
o.dt = 1e-12
outputDir = 'Output'
checkCounter = 0
numWinds = 5
WinAction = np.zeros((5,1), dtype=int)

## =========== Environment ============##
widthSpace = 6.8 ## width is the distance from the front coal wall to the back
widthHydr = 1.5
lengthSpace = widthHydr * numWinds ## length is the width of 5 windows
highHydr = 3.8
CheckThick = 1
highDummy = 3
colorWind = (0, 1, 1)
angleShield = 50 * 3.1415926 / 180
angleSwingPositive = 15 * 3.1415926 / 180
angleSwingNegtive = 40 * 3.1415926 / 180
lengthShield = 3
lengthTail = 2
windUpperBoundary = 0.9
windLowerBoundary = 0.5019
stateUpperBoundary = highHydr
stateLowerBoundary = windLowerBoundary
positionWind = []
windPosition = np.zeros(5)

shield_y_0 = lengthTail * math.cos(angleShield - angleSwingPositive)
shield_y_1 = shield_y_0 + lengthShield * math.cos(angleShield)
shield_z_0 = highHydr - lengthShield * math.sin(angleShield)
shield_z_1 = highHydr #
wind_y_0 = lengthTail * (math.cos(angleShield - angleSwingPositive) - math.cos(angleShield))
wind_y_1 = shield_y_0
wind_z_0 = highHydr - (lengthShield + lengthTail) * math.sin(angleShield)
wind_z_1 = highHydr - lengthShield * math.sin(angleShield)

HyMat = O.materials.append(FrictMat(frictionAngle=0.1, density=3000, young=2e9))
myGraviaty = -1200
def Relocation():
    global nEpisode, nIter, index
    path = 'Output/location/'
    wind_name = path + 'wind_' + str(nEpisode) + '_' + str(nIter) + '.npy'
    if nIter >= index[nEpisode, 1]:
        nIter = 0
        nEpisode += 1
    else:
        nIter += 1
    if nEpisode >= len(index):
        O.pause()

    locationWind = np.load(wind_name, allow_pickle=True)
    print(locationWind)

    for i in range(0, len(locationWind)):
        n = locationWind[i][0]
        o.bodies[n].state.pos = locationWind[i][1]

def HydraulicSupport():
    for i in range(0, numWinds):
        temp = [
            Vector3(widthHydr * i, wind_y_0, wind_z_0),
            Vector3(widthHydr * (i + 1), wind_y_0, wind_z_0),
            Vector3(widthHydr * (i + 1), wind_y_1, wind_z_1),
            Vector3(widthHydr * i, wind_y_1, wind_z_1)
        ]
        positionWind.append(temp)

    Wind1 = pack.sweptPolylines2gtsSurface([positionWind[0]], capStart=True, capEnd=True)
    Wind2 = pack.sweptPolylines2gtsSurface([positionWind[1]], capStart=True, capEnd=True)
    Wind3 = pack.sweptPolylines2gtsSurface([positionWind[2]], capStart=True, capEnd=True)
    Wind4 = pack.sweptPolylines2gtsSurface([positionWind[3]], capStart=True, capEnd=True)
    Wind5 = pack.sweptPolylines2gtsSurface([positionWind[4]], capStart=True, capEnd=True)

    IDWind1 = O.bodies.append(pack.gtsSurface2Facets(Wind1, color=colorWind))
    IDWind2 = O.bodies.append(pack.gtsSurface2Facets(Wind2, color=colorWind))
    IDWind3 = O.bodies.append(pack.gtsSurface2Facets(Wind3, color=colorWind))
    IDWind4 = O.bodies.append(pack.gtsSurface2Facets(Wind4, color=colorWind))
    IDWind5 = O.bodies.append(pack.gtsSurface2Facets(Wind5, color=colorWind))

    IDWind = [IDWind1, IDWind2, IDWind3, IDWind4, IDWind5]
    WindList = IDWind1 + IDWind2 + IDWind3 + IDWind4 + IDWind5
    return IDWind, WindList

IDWind, WindList = HydraulicSupport()
import os.path
rootdir = "Output/location/"
wind_locations = []
for parent, dirnames, filenames in os.walk(rootdir):
    for filename in filenames:
        if filename[0:4] == 'wind':
            wind_locations.append(filename)
episodeMax = 0

for i in wind_locations:
   n = int(i[5])
   if n >= episodeMax:
       episodeMax = n

index = np.zeros((episodeMax+1, 2))
for i in range(episodeMax+1):
    index[i,0] = i
    for name in wind_locations:
        episode = int(name[5])
        iter = int(name[7])
        if episode == i:
            if (iter) >= index[i, 1]:
                index[i, 1] = iter

nEpisode =0
nIter = 0

O.engines = [
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()]),
    InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom(), Ig2_Wall_Sphere_ScGeom()],
        [Ip2_FrictMat_FrictMat_FrictPhys(),
         Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(xSectionWeibullScaleParameter=0.5,
        xSectionWeibullShapeParameter=0.5,
        weibullCutOffMin=0,
        weibullCutOffMax=10)],
        [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=False, recordMoments=False,label='interactionLaw'),
        Law2_ScGeom_FrictPhys_CundallStrack()]
        ),
    GlobalStiffnessTimeStepper(),
    NewtonIntegrator(gravity=(0, 0, myGraviaty), damping=0.5, label='down'),
    PyRunner(command="Relocation()", iterPeriod=200000),

]

The position loaded in the second script is not right. That seems the stored location is wrong.
Thank you very much.

Jan Stránský (honzik) said : #3

Instant replay:
> M stands for minimal. Do you have the same problem with one plank? then please provide a script with just one plank..

please read [1]

A true MWE for this case would be:
- create ONE (or TWO) facet(s)
- move / rotate it a bit, a few iterations is enough
- save it
- load it
does it work ?
- no, then ask your question providing this true MWE
- yes, then add some other minimal part of the code, goto "does it work?"

> The position loaded in the second script is not right

how do you know?
what is "position"?
I have tested your code, and both saved and loaded printed values are the same..

> That seems the stored location is wrong.

seems? either it is or not.. and from the printed values they are correct..

The standard way to view / visualize the simulation is to export data for a dedicated software for visualization, like Paraview.
You can use Yade to visualize the results, of course, but then I don't understand most of the second script, which seems to be "too much simulation" and too similar to the first script to be the visualization script..

cheers
Jan

[1] https://www.yade-dem.org/wiki/Howtoask

yang yi (yangyizaixian) said : #4

To Jan Stránský:
Thank you very much. I modified the script as follow, there just a plank. In the first It rotates around a axis and the location stored. Then the location load by the second script and display in the 3D show.
The plank consist of two triangles. The 3D show of the first script that the two triangles rotate together, however in the display script, the triangles are separated. As you found the location is the same. So please help me how to modify that.

(1) Rotation script
#!/usr/bin/env python
# encoding: utf-8
from __future__ import print_function
import sys

sys.path.append('~/PycharmProjects/Yade20191229/')
# from yadeImport import *
import numpy as np
import os
from yade.gridpfacet import *
o = Omega()
o.dt = 1e-12
outputDir = 'Output'
positionWind = []

def HydraulicSupport():
    global positionWind
    temp = [Vector3(0, 0, 0),
            Vector3(0, 10, 0),
            Vector3(10, 10, 0),
            Vector3(10, 0, 0)]
    positionWind.append(temp)
    Wind1 = pack.sweptPolylines2gtsSurface([positionWind[0]], capStart=True, capEnd=True)
    IDWind = O.bodies.append(pack.gtsSurface2Facets(Wind1))
    return IDWind

def Ground():
    O.bodies.append(utils.wall(position=-10, sense=0, axis=2, material=-1))

##---------------------------------------##
def WindowsAction(IDWind):
    global WinAction, windPosition, saveCounter, nEpisode,RotationW1
    RotationW1.angularVelocity = -0.000001
    SaveProcessLocation(nEpisode, saveCounter)
    saveCounter += 1

def SaveProcessLocation(Episode, iter):
    outputDir = "Output/location/"
    if not os.path.exists(outputDir):
        os.makedirs(outputDir)
    Position = []

    for i in IDWind:
        print(o.bodies[i].state.pos)
        temp = [i, o.bodies[i].state.pos]
        Position.append(temp)
    location_name = outputDir + 'wind_' + str(Episode)+ '_' + str(iter)
    np.save(location_name, Position)

def WindowsControl():
    global saveCounter, nEpisode
    nEpisode += 1
    saveCounter = 0

nEpisode = 0
saveCounter = 0
IDWind = HydraulicSupport()
Ground()

O.engines = [
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()]),
    InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom(), Ig2_Wall_Sphere_ScGeom()],
        [Ip2_FrictMat_FrictMat_FrictPhys(),
         Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(xSectionWeibullScaleParameter=0.5,
        xSectionWeibullShapeParameter=0.5,
        weibullCutOffMin=0,
        weibullCutOffMax=10)],
        [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=False, recordMoments=False,label='interactionLaw'),
        Law2_ScGeom_FrictPhys_CundallStrack()]
        ),
    GlobalStiffnessTimeStepper(),
    NewtonIntegrator(gravity=(0, 0, 9.8), damping=0.5, label='down'),
    RotationEngine(rotationAxis=(1, 0, 0), rotateAroundZero=True, ids=IDWind, zeroPoint=positionWind[0][0],
                   label='RotationW1'),
    PyRunner(command="WindowsAction(IDWind)", iterPeriod=200000),
    PyRunner(command="WindowsControl()", iterPeriod=2000000),
]
# o.run(1200000000)

================================================================
(2) Display script
from __future__ import print_function
import sys
sys.path.append('~/PycharmProjects/Yade20191229/')
# from yadeImport import *
import numpy as np
import os
from yade.gridpfacet import *
o = Omega()
o.dt = 1e-12

positionWind=[]

def HydraulicSupport():
    global positionWind
    temp = [Vector3(0, 0, 0),
            Vector3(0, 10, 0),
            Vector3(10,10, 0),
            Vector3(10,0, 0)]
    positionWind.append(temp)
    Wind1 = pack.sweptPolylines2gtsSurface([positionWind[0]], capStart=True, capEnd=True)
    IDWind = O.bodies.append(pack.gtsSurface2Facets(Wind1))
    return IDWind

def Ground():
    O.bodies.append(utils.wall(position=-10, sense=0, axis=2, material=-1))

def Relocation():
    global nEpisode, nIter, index
    path = 'Output/location/'
    wind_name = path + 'wind_' + str(nEpisode) + '_' + str(nIter) + '.npy'

    if nIter >= index[nEpisode, 1]:
        nIter = 0
        nEpisode += 1
    else:
        nIter += 1
    if nEpisode >= len(index):
        O.pause()

    locationWind = np.load(wind_name, allow_pickle=True)
    print(nEpisode,nIter)
    print(locationWind)
    print('--------------')
    for i in range(0, len(locationWind)):
        n = locationWind[i][0]
        o.bodies[n].state.pos = locationWind[i][1]

import os.path
rootdir = "Output/location/"
wind_locations = []

for parent, dirnames, filenames in os.walk(rootdir):
    for filename in filenames:
        if filename[0:4] == 'wind':
            wind_locations.append(filename)

episodeMax = 0
endNumber = 100000000
for i in wind_locations:
   n = int(i[5])
   if n >= episodeMax:
       episodeMax = n

index = np.zeros((episodeMax+1, 2))
for i in range(episodeMax+1):
    index[i,0] = i
    for name in wind_locations:
        episode = int(name[5])
        iter = int(name[7])
        if episode == i:
            if (iter) >= index[i, 1]:
                index[i, 1] = iter

IDWind = HydraulicSupport()
Ground()
nEpisode =0
nIter = 0

O.engines = [
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()]),
    InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom(), Ig2_Wall_Sphere_ScGeom()],
        [Ip2_FrictMat_FrictMat_FrictPhys(),
         Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(xSectionWeibullScaleParameter=0.5,
        xSectionWeibullShapeParameter=0.5,
        weibullCutOffMin=0,
        weibullCutOffMax=10)],
        [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=False, recordMoments=False,label='interactionLaw'),
        Law2_ScGeom_FrictPhys_CundallStrack()]
        ),
    GlobalStiffnessTimeStepper(),

    NewtonIntegrator(gravity=(0, 0, 9.8), damping=0.5, label='down'),
    PyRunner(command="Relocation()", iterPeriod=100000),
]
# o.run(1200000000, True)

Thank you very much

Jan Stránský (honzik) said : #5

Sorry, but this is still too far from MWE. If you want help from us, you have to do some work on your side first.

If your problem is saving/loading a facet, then the MWE should be about saving and loading one facet:

### 1st script
f = facet(((0,0,0),(1,0,0),(0,1,0)))
O.bodies.append(f)
f.state.vel = (.1,0,0)
f.state.angVel = (0,0,.1)
O.run(10,True)
numpy.save("f",[[f.id,f.state.pos]])
###

### 2nd script
f = facet(((0,0,0),(1,0,0),(0,1,0)))
O.bodies.append(f)
[[i,pos]] = numpy.load("f.npy",allow_pickle=True)
f.state.pos = pos
###

is this the same problem you are facing?

cheers
Jan

yang yi (yangyizaixian) said : #6

To Jan Stránský :
Thank you very much. I am sorry that the script can not meet MWE. My script is similar with yours. The facet rotate and restore the location , then display the rotation process by reload the location. How ever my facet is a quadrangle, and consist with two triangles.
      I used the " pack.sweptPolylines2gtsSurface([positionWind[0]], capStart=True, capEnd=True)" to create the plank, but I do not know why the quadrangle consists with two triangles. That is source of my problem. In the rotation process, the quadrangle moving with the two triangles sticks together in the 3D show. However, If I read the location from the storage, in the 3D show, the two triangles separated.

    I desire your help. so I modified the script as follow, I hope it can meet the MWE

(1) Rotation process:

positionWind = []
temp = [Vector3(0, 0, 0),
        Vector3(0, 10, 0),
        Vector3(10, 10, 0),
        Vector3(10, 0, 0)]
positionWind.append(temp)
Wind1 = pack.sweptPolylines2gtsSurface([positionWind[0]], capStart=True, capEnd=True)
IDWind = O.bodies.append(pack.gtsSurface2Facets(Wind1))
saveCounter = 0
##---------------------------------------##
def WindowsAction_SaveLoaction(IDWind):
    global WinAction, windPosition, saveCounter, RotationW1
    RotationW1.angularVelocity = -0.000001
    outputDir = "Output/location/"
    if not os.path.exists(outputDir):
        os.makedirs(outputDir)
    Position = []
    for i in IDWind:
        temp = [i, o.bodies[i].state.pos]
        Position.append(temp)
    location_name = outputDir + 'wind' + '_' + str(saveCounter)
    np.save(location_name, Position)
    saveCounter += 1

O.engines = [
    ForceResetter(),
    GlobalStiffnessTimeStepper(),
    NewtonIntegrator(gravity=(0, 0, 9.8), damping=0.5, label='down'),
    RotationEngine(rotationAxis=(1, 0, 0), rotateAroundZero=True, ids=IDWind, zeroPoint=positionWind[0][0],
                   label='RotationW1'),
    PyRunner(command="WindowsAction_SaveLoaction(IDWind)", iterPeriod=200000),
]

==========================================
(2) Read the location

positionWind = []
temp = [Vector3(0, 0, 0),
        Vector3(0, 10, 0),
        Vector3(10, 10, 0),
        Vector3(10, 0, 0)]
positionWind.append(temp)
Wind1 = pack.sweptPolylines2gtsSurface([positionWind[0]], capStart=True, capEnd=True)
IDWind = O.bodies.append(pack.gtsSurface2Facets(Wind1))

saveCounter = 0
import os.path
rootdir = "Output/location/"
wind_locations = []
for parent, dirnames, filenames in os.walk(rootdir):
    for filename in filenames:
        if filename[0:4] == 'wind':
            wind_locations.append(filename)
fileNum = len(filenames)
saveCounter=0

def Display():
    global fileNum,saveCounter
    path = 'Output/location/'
    wind_name = path + 'wind_' + str(saveCounter) + '.npy'
    locationWind = np.load(wind_name, allow_pickle=True)
    print(locationWind)
    print('--------------')
    for i in range(0, len(locationWind)):
        n = locationWind[i][0]
        o.bodies[n].state.pos = locationWind[i][1]
    saveCounter +=1

O.engines = [
    ForceResetter(),
    GlobalStiffnessTimeStepper(),
    NewtonIntegrator(gravity=(0, 0, 9.8), damping=0.5, label='down'),
    PyRunner(command="Display()", iterPeriod=200000),
]

Thank you very much

Jan Stránský (honzik) said : #7

> my facet is a quadrangle, and consist with two triangles.
> I used the " pack.sweptPolylines2gtsSurface([positionWind[0]], capStart=True, capEnd=True)" to create the plank, but I do not know why the quadrangle consists with two triangles. That is source of my problem.

sweptPolylines2gtsSurface - gts means GNU **Triangulated** Surface [1], therefore the result are triangles
Yade facet = triangle. Yade currently does not support more general planar shapes.
You can use box [2] with zero or "small" thick if your aim is a rectangle.

> My script is similar with yours.

Ok, now we are targeting the source of the problem :-)
Yade Body has a State, which has position (pos), but also an orientation (ori). Depending on the shape, to be able to restore the position AND orientation, you should save state.pos, but also state.ori.

> I modified the script as follow, I hope it can meet the MWE

getting closer, thanks :-) but I will update my MWE:
### script 1
f = facet(((0,0,0),(1,0,0),(0,1,0)))
O.bodies.append(f)
f.state.vel = (.1,0,0)
f.state.angVel = (0,0,.1)
O.run(10,True)
data = []
for b in O.bodies:
   axis,angle = b.state.ori.toAxisAngle()
   data.append([b.id,b.state.pos,axis,angle])
numpy.save("f",data)
### script 2
f = facet(((0,0,0),(1,0,0),(0,1,0)))
O.bodies.append(f)
print(numpy.load("f.npy",allow_pickle=True))
[[i,pos,axis,angle]] = numpy.load("f.npy",allow_pickle=True)
f.state.pos = pos
f.state.ori = Quaternion(axis,angle)
###

a few notes:
- as pointed above, consider an external program for visualization
- instead of pos or ori, you can save vertices of facets, which is somehow more general than just pos and ori and depending on the geometry to be hard-coded (similarly to other shapes e.g. you switch to box)

cheers
Jan

[1] http://gts.sourceforge.net/
[2] https://yade-dem.org/doc/yade.utils.html#yade.utils.box

Can you help with this problem?

Provide an answer of your own, or ask yang yi for more information if necessary.

To post a message you must log in.