Coordinates of Mesh-Nodes

Hi Oguz,

ESyS-Particle currently does not have a method to obtain the coordinates of mesh nodes during simulations.

In the past, when I have wanted to rotate a mesh, I have written a runnable that initially reads the mesh file to store the initial coordinates and ids of mesh nodes. Then, during the Runnable.run() call, I update the positions of the nodes and inform the simulation to move the nodes via a LsmMpi.moveSingleMeshNodeBy(..) call.

This probably isn't the most elegant solution but it works...

Cheers,

Dion

Hi Dion,
Is it possible, that you can send me that Runnable? So i have an example how to do that. Only if it's okay for you :)

Greetings
Oguz

i tried here a runnable:

from __future__ import division
from esys.lsm import *
from esys.lsm.util import *
import math
import numpy as np
from numpy import matrix,array

class kinematic (Runnable):
   def __init__ (self,
                 LsmMpi=None,
                 meshName=None,
                 startTime=0,
                 rampTime = 200):
      """
      Subroutine to initialise the Runnable and store parameter values.
      """
      Runnable.__init__(self)
      self.sim = LsmMpi
      self.meshName = meshName
      self.dt = self.sim.getTimeStepSize()
      self.rampTime = rampTime
      self.startTime = startTime
      self.Nt = 0

   def run (self):
  alpha=((2*math.pi)/360)*45*self.Nt
  R_z = np.matrix( [[math.cos(alpha),-math.sin(alpha),0],[math.sin(alpha),math.cos(alpha),0],[0,0,1]] )
  #node coordinates
  node0 = np.array( ((-5),(0),(0)) )
  node1 = np.array( ((-5),(0),(5)) )
  node2 = np.array( ((0),(0),(-5)) )
  node3 = np.array( ((0),(0),(0)) )
  node4 = np.array( ((5.0),(0),(-5)) )
  node5 = np.array( ((5.0),(0),(5)) )
  rotate0 = np.dot(node0,R_z)
  rotate1 = np.dot(node1,R_z)
  rotate2 = np.dot(node2,R_z)
  rotate3 = np.dot(node3,R_z)
  rotate4 = np.dot(node4,R_z)
  rotate5 = np.dot(node5,R_z)
       #moving nodes
  self.sim.moveSingleMeshNodeBy(
   self.meshName,
   nodeId = 0,
   delta = Vec3(rotate0[0,0]-node0[0],rotate0[0,1]-node0[1],rotate0[0,2]-node0[2])
  )

  self.sim.moveSingleMeshNodeBy(
   self.meshName,
   nodeId = 1,
   delta = Vec3(rotate1[0,0]-node1[0],rotate1[0,1]-node1[1],rotate1[0,2]-node1[2])
  )

  self.sim.moveSingleMeshNodeBy(
   self.meshName,
   nodeId = 2,
   delta = Vec3(rotate2[0,0]-node2[0],rotate2[0,1]-node2[1],rotate2[0,2]-node2[2])
  )

  self.sim.moveSingleMeshNodeBy(
   self.meshName,
   nodeId = 3,
   delta = Vec3(rotate3[0,0]-node3[0],rotate3[0,1]-node3[1],rotate3[0,2]-node3[2])
  )

  self.sim.moveSingleMeshNodeBy(
   self.meshName,
   nodeId = 4,
   delta = Vec3(rotate4[0,0]-node4[0],rotate4[0,1]-node4[1],rotate4[0,2]-node4[2])
  )

  self.sim.moveSingleMeshNodeBy(
   self.meshName,
   nodeId = 5,
   delta = Vec3(rotate5[0,0]-node5[0],rotate5[0,1]-node5[1],rotate5[0,2]-node5[2])
   )
          self.Nt += .00075

where R_z is the rotation matrix in Z direction, and alpha the angular ,which changes every step. My problem is, that the mesh behaves very strange wenn I visualize it. It gets infinital long and does not rotate correctly. The mesh is the example from the Tutorial.
Here is the part from the main code:

sim.readMesh(
    fileName = "trianglemesh.lsm",
    meshName = "mesh"
)

movemesh = kinematic(
 LsmMpi = sim,
 meshName = "mesh",
    startTime = 0,
    rampTime = 0
)
sim.addPreTimeStepRunnable (movemesh)

Did i miss something? maybe it's a dump question :D
Cheers Oguz

Hi Oguz,

I think the problem might be that you are moving each node by an amount equal to the difference between the current desired position of the nodes (rotate*) and the initial locations of the nodes (node*). LsmMpi.moveSingleMeshNodeBy(..) expects you provide an incremental displacement at each calling timestep i.e. the difference between the desired node positions (rotate*) and the "current" node positions (not the initial positions).

Perhaps your problem will be solved by setting node* = rotate* after the six moveSingleMeshNodeBy(..) calls?

Cheers,

Dion

Hi,
First of all thanks for your help! I tried your approach but nothing changed.
There is some problem with the "delta" in "self.sim.moveSingleMeshNodeBy(...)". If i choose "delta = Vec3(0,0,0)" nothing happens, like its supposed to be. But when I choose "delta = Vec3(2,2,2)", an constent value, the mesh splatters and disappers like before.

Greetings
Oguz

Hi Oguz,

To get your runnable to work correctly you must do two things:

1) set the rotation angle to be a constant = the angle by which you wish to rotate the mesh in a single timestep. Currently alpha increases with time and this is the main reason your mesh explodes. The script is progressively rotating the mesh by larger and larger angles each timestep.

2) ensure that you set node* = rotate* after the sim.moveSingleMeshNode(..) calls

I suspect the reason delta=Vec3(2,2,2) splatters the mesh is because this displacement increment is too large. Consequently the simulation becomes numerically unstable. To avoid this, displacement increments should always be smaller than the verletDist...ideally considerably smaller. If a mesh node moves by more than the verletDist, the neighbour table must be re-built (just as if a particle had moved by more than the verletDist).

Pasted below is my modified version of your Runnable (rotate.py) and a small simulation script (sim.py) to test out the mesh rotations. Having run the simulation and generated snapshot files, run the following to generate VTK files of the rotating mesh:
dump2vtk -i snapshot -o snaps_ -mesh -t 0 1001 1
These can then be visualised in Paraview.

Have fun!

Cheers,

Dion

#rotate.py:
from __future__ import division
from esys.lsm import *
from esys.lsm.util import *
import math
import numpy as np
from numpy import matrix,array

class kinematic (Runnable):
   def __init__ (self,
                 LsmMpi=None,
                 meshName=None,
                 startTime=0,
                 rampTime = 200):
      """
      Subroutine to initialise the Runnable and store parameter values.
      """
      Runnable.__init__(self)
      self.sim = LsmMpi
      self.meshName = meshName
      self.rampTime = rampTime
      self.startTime = startTime
      self.Nt = 0
      self.Nt = .00075

   def run (self):
      alpha=((2*math.pi)/360)*45*self.Nt
      R_z = np.matrix( [[math.cos(alpha),-math.sin(alpha),0],[math.sin(alpha),math.cos(alpha),0],[0,0,1]] )
      #node coordinates
      node0 = np.array( ((-5),(0),(0)) )
      node1 = np.array( ((-5),(0),(5)) )
      node2 = np.array( ((0),(0),(-5)) )
      node3 = np.array( ((0),(0),(0)) )
      node4 = np.array( ((5.0),(0),(-5)) )
      node5 = np.array( ((5.0),(0),(5)) )
      rotate0 = np.dot(node0,R_z)
      rotate1 = np.dot(node1,R_z)
      rotate2 = np.dot(node2,R_z)
      rotate3 = np.dot(node3,R_z)
      rotate4 = np.dot(node4,R_z)
      rotate5 = np.dot(node5,R_z)
           #moving nodes
      self.sim.moveSingleMeshNodeBy(
       self.meshName,
       nodeId = 0,
       delta = Vec3(rotate0[0,0]-node0[0],rotate0[0,1]-node0[1],rotate0[0,2]-node0[2])
      )

      self.sim.moveSingleMeshNodeBy(
       self.meshName,
       nodeId = 1,
       delta = Vec3(rotate1[0,0]-node1[0],rotate1[0,1]-node1[1],rotate1[0,2]-node1[2])
      )

      self.sim.moveSingleMeshNodeBy(
       self.meshName,
       nodeId = 2,
       delta = Vec3(rotate2[0,0]-node2[0],rotate2[0,1]-node2[1],rotate2[0,2]-node2[2])
      )

      self.sim.moveSingleMeshNodeBy(
       self.meshName,
       nodeId = 3,
       delta = Vec3(rotate3[0,0]-node3[0],rotate3[0,1]-node3[1],rotate3[0,2]-node3[2])
      )

      self.sim.moveSingleMeshNodeBy(
       self.meshName,
       nodeId = 4,
       delta = Vec3(rotate4[0,0]-node4[0],rotate4[0,1]-node4[1],rotate4[0,2]-node4[2])
      )

      self.sim.moveSingleMeshNodeBy(
       self.meshName,
       nodeId = 5,
       delta = Vec3(rotate5[0,0]-node5[0],rotate5[0,1]-node5[1],rotate5[0,2]-node5[2])
      )
      """
      """
      #self.Nt += .00075
      #print '#',self.Nt
      node0[0] = rotate0[0,0]
      node1[0] = rotate1[0,0]
      node2[0] = rotate2[0,0]
      node3[0] = rotate3[0,0]
      node4[0] = rotate4[0,0]
      node5[0] = rotate5[0,0]

      node0[1] = rotate0[0,1]
      node1[1] = rotate1[0,1]
      node2[1] = rotate2[0,1]
      node3[1] = rotate3[0,1]
      node4[1] = rotate4[0,1]
      node5[1] = rotate5[0,1]
      """
      print rotate0[0,0],rotate0[0,1],rotate0[0,2]
      print rotate1[0,0],rotate1[0,1],rotate1[0,2]
      print rotate2[0,0],rotate2[0,1],rotate2[0,2]
      print rotate3[0,0],rotate3[0,1],rotate3[0,2]
      print rotate4[0,0],rotate4[0,1],rotate4[0,2]
      print rotate5[0,0],rotate5[0,1],rotate5[0,2]
      """

if __name__=="__main__":
   runner = kinematic()
   runner.run()

#sim.py:
from esys.lsm import *
from esys.lsm.util import *
from rotate import kinematic

sim = LsmMpi(1,[1,1,1])
sim.initNeighbourSearch (
   particleType = "NRotSphere",
   gridSpacing = 2.5000,
   verletDist = 0.1000
)

sim.setNumTimeSteps (1000)
sim.setTimeStepSize (1.0000e-04)

domain = BoundingBox(Vec3(-20,-20,-20), Vec3(20,20,20))
sim.setSpatialDomain(domain)

sim.readMesh(
   fileName = "mesh.lsm",
   meshName = "mesh_wall"
)

runner = kinematic(LsmMpi = sim, meshName="mesh_wall")
sim.addPreTimeStepRunnable(runner)

sim.createCheckPointer (
   CheckPointPrms (
      fileNamePrefix = "snapshot",
      beginTimeStep = 0,
      endTimeStep = 1000,
      timeStepIncr = 1
   )
)

sim.run()