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()

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:
For:
Assignee:
No assignee Edit question
Last query:
2020-05-29
2020-06-01
 Jan Stránský (honzik) said on 2020-05-29: #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.

cheers
Jan

 yang yi (yangyizaixian) said on 2020-05-29: #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

import math
import numpy as np
import os
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
import math
import numpy as np
import os
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()

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 on 2020-05-31: #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..

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
does it work ?
- 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

 yang yi (yangyizaixian) said on 2020-06-01: #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

import numpy as np
import os
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
import numpy as np
import os
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()

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 on 2020-06-01: #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.

### 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)
f.state.pos = pos
###

is this the same problem you are facing?

cheers
Jan

 yang yi (yangyizaixian) said on 2020-06-01: #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),
]

==========================================

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'
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 on 2020-06-01: #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)
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