Possible bug of a two sided wall using mesh wall?

Hi Feng,

TriMeshes are not designed to work in 2D simulations. There is a 2D equivalent of the TriMesh called a LinMesh. This is a 2D surface consisting of line-segments. Unfortunately I have not used this so am not so sure how it works or the expected mesh-file format. I've had a look at the source-code and here's my best guess.

Try the following format for the 2D meshfile (although it may not be entirely correct):
2D-Nodes 3
0 0 1 10.0 10.0
1 1 1 10.0 20.0
2 2 1 20.0 20.0
Line2 2
0 0 0 1
1 0 1 2

For the nodes list, each node line is [nodeID] [dummy] [tag] [X] [Y]
For the lines list, each line is [lineID] [tag] [nodeID1] [nodeID2]

Looks like you read the 2D meshfile using:

LsmMpi.readMesh2D (
   fileName = "myMesh.msh",
   meshName = "my2Dmesh",
   tag = 1
)

and define elastic interactions via:

LsmMpi.createInteractionGroup (
   NRotElasticLinMeshPrms (
       name = my2Dmesh_repell",
       meshName = "my2Dmesh",
       normalK = 12345.0
   )
)

I suspect 2D linear meshes are also one-way so you will need to put two of them back-to-back for a two-way wall. I've no idea the correct order for nodes so you may have to experiment with that.

Good luck and let me know how you go with this.

Cheers,

Dion.

Hi, Dr. Dion Weatherley:

After playing with the 2D mesh for 2 days I think the 2D piece-wise mesh are not working properly as expected, (even if we change the order of the nodes to clockwise/counterclockwise), I attached my scripts, mesh files and screenshots below:

https://sites.google.com/site/fengchenmaple/test

basically this seems resemble the type of problem described in the following thread:

https://answers.launchpad.net/esys-particle/+question/98840

My experience is without the mesh wall, the simulation runs quickly and accurately and finish around 200sec(50000steps) on my desktop, however, due to the strange oscillating and trembling behavior of the particles, the simulation are very slow (need around 4000sec), and definitely incorrect results. I also tried switch to the latest svn version (svn-1350), with the same result.

Below are what I learned by looking at the source code, in order to use the 2D mesh, we need to use the following in the python script:

sim.readMesh2D (

   fileName = "test_mesh2d.msh",

   meshName = mshName,

   tag = 1

)

sim.createInteractionGroup (

   NRotElasticMesh2DPrms ( #However the source code says the Mesh2DPrms is deprecated...

       name = "mesh2d_repell",

       meshName = mshName,

       normalK = kn,

   )

)

If we need to use NRotElasticLinMeshPrms, we need to add:

.def(

    "createInteractionGroup",

    &LsmMpiPy::createNRotElasticLinMeshInteractGrp,

    (arg("prms"))

)

and

void LsmMpiPy::createNRotElasticLinMeshInteractGrp(

      const NRotElasticLinMeshPrmsPy &prms

    )

{

      getLatticeMaster().addMesh2DIG(prms); //I copied the same thing from NRotElasticMesh2DPrmsPy

}

in LsmMpiPy.cpp

and declaration of

"class NRotElasticLinMeshPrmsPy"

void createNRotElasticLinMeshInteractGrp(const NRotElasticLinMeshPrmsPy &prms);

in LsmMpiPy.h

Having done all these, I tried both "NRotElasticLinMeshPrms" and "NRotElasticMesh2DPrms", however it seems neither is working properly, especially the particles at the back of the active side, is there any method I can trace this problem or I am still not using the mesh wall in a correct way?

Thanks for any suggestions!

Feng

Hi Feng,

I think you may be correct that the problem is similar to the notorious "Hanging spheres" bug for triangle meshes. We did not check if the same bug existed for the 2D mesh interaction at the time. Looking at the source code, it would appear that the bug still exists for the 2D meshes.

Since you have a development version already, try adding an extra check of the separation between corners and particles in ECorner2DInteraction.cpp. ECorner2DInteraction::calcForces() should look like this:

void ECorner2DInteraction::calcForces()
{
  Vec3 ppos=m_p->getPos();
  if(m_corner->isValidContact(ppos)){ // if no contact to adjacent edges or triangles
    double sep=m_corner->sep(ppos);
    if (sep < m_p->getRad()) {
       Vec3 force=m_k*(m_p->getRad()-sep)*m_corner->getDirectionFromPoint(ppos);
       m_p->applyForce(force,ppos);
    }
  }
}

Hopefully the bug fix above will solve your oscillating particles problem. I suspect the explosion of the model may be a different problem. Your two meshes are very close to one-another (only 0.1 units apart) If a particle overlaps one of these meshes by more than 0.1 units, it will experience an incorrect force from the back side of the other mesh. You may need to either separate your meshes by a greater distance or make your mesh stiffnesses much larger to avoid this.

Thanks for information on the other changes you made. These need to be incorporated into the development version. All of these are updates to the Python interface to the NRotElasticMesh2DPrms interaction but none of them change the way forces are calculated. Consequently, you should get the same behaviour whether you use NRotElasticMesh2DPrms or NRotElasticLinMeshPrms.

Cheers,

Dion.

Hi, Dr. Dion Weatherley:

After making the ECorner2DInteraction::calcForces() change, the problem was solved, thanks for the fix. As it might be useful for later users with 2D mesh, I also attached a sample figure showing the active side of the 2D mesh (line), which is still critical to make simulation correct. My experience with the 2D piece-wise mesh direction is somewhat similar to the right hand rule but with only two fingers (do not guarantee to be correct, only from my experience):

http://i55.tinypic.com/8y6id0.jpg

Below is a 2-sided 2D mesh wall mesh file example (vertical, 0.2 unit in thickness, also work with thickness 0.1):

2D-Nodes 4
0 0 1 20.0 10.0
1 1 1 20.0 20.0
2 2 1 20.2 20.0
3 3 1 20.2 10.0
Line2 3
0 1 0 1
1 1 2 3
3 1 3 0

To make this complete, in order to use it in the python script, we need to use:

mshName="2dshplWall"

sim.readMesh2D (
   fileName = "test_mesh2d.msh",
   meshName = mshName,
   tag = 1
)

sim.createInteractionGroup (
   NRotElasticMesh2DPrms (
       name = "mesh2d_repell",
       meshName = mshName,
       normalK = kn,
   )
)

This should work for both svn and the stable 2.0 version.

Feng

Also, just a quick thought, how to make the 2D mesh visible with the dump2vtk or the dump2pov command, because currently the 2D mesh seems not recorded in the checkpointer file, although there is writeCheckPoints function in the 2D Mesh implementation, i.e. how to invoke the writeCheckPoint function for the 2D mesh?

Feng

Hi Feng,

Great to hear your simulation is working now! Thanks for the additional information about use of 2D meshes and the "right-hand rule". That's a really nice illustration.

I'll ask one of the developers to add the bug-fix to the development trunk next week and to have a look into how to record 2D mesh information in the checkpoint files, dump2vtk etc. We will need to determine whether the current 2D mesh writeCheckPoints function is complete or a stub that requires update before invoking it during checkPointing.

Have fun!

Dion.

Hi, Dr. Dion Weatherley:

Thanks for your response, a quick modification to the "right hand rule" figure, the mid-finger also needs to be used and always pointed to the negative Z direction ("-Z"), otherwise the two finger rule will be useless because the active side (thumb finger) can be either side :-))), see the updated figure:

http://i54.tinypic.com/2008g1j.jpg

Feng

Hi Feng:

Thanks for confirming the code changes needed for your simulation to work correctly, and for picking up on the additional code required in Python/esys/lsm/LsmMpiPy.cpp for NRotElasticLinMeshPrms, which was an oversight on my part. I have committed the changes to the repository and checked your script's results using your updated mesh. The particles no longer explode from the wall or pass through the wall, although you may wish to verify this with the new development version of ESyS-Particle. Now I will get cracking on the problem of outputting mesh information to checkpoint files.

Thank you for all your valuable input.

Vince