add element to tetra mesh sim problem

[chc2]

i am changing a tetra mesh by modifying existing tetras and adding new tetras. I am using the recommended methods like using topologicalModifiers. in terms of the mesh it seems like its working.

the problem is the sim behaviour after i have changed the topology. it almost seems as if the tetras have no mass or gravity.

this video demonstrates the problem

https://www.youtube.com/watch?v=yoyBpAV36_Y

i made another video where i start a fresh new sim that starts with the modified mesh. as you can see in the video it seems to behave as expected.

https://www.youtube.com/watch?v=1XVaM_68LEQ

i checked the new tetras and they seem to have a force, mass, velocity, etc.

do i need to reinitialize the sim or something?

the sofa version is pretty old and was retrieved before sofa was moved to github.

[Hugo]

hi chc2,

Could you tell us which physics components you are using ? what ForceField and what Mass ?
Because all components in SOFA do not support topological changes.

What is your SOFA version btw ?

Cheers,

Hugo

[chc2]

<Node          name=”tetramesh”  gravity=”0 -10 0″  dt=”0.01″  time=”0″  animate=”0″  >
                <MechanicalObject template=”Vec3d” name=”DOF”  />
                <TetrahedronSetTopologyContainer name=”Container”  nbPoints=”21″  position=”@DOF.position”  />
                <EulerImplicitSolver name=”solver” />
                <CGLinearSolver template=”GraphScattered” name=”linearSolver”  iterations=”2000″  tolerance=”1e-15″  threshold=”1e-15″ />
                <UniformMass template=”Vec3d” name=”mass”  mass=”1″  totalmass=”21″ />
                <TetrahedralCorotationalFEMForceField template=”Vec3d” name=”FEM”  tetrahedronInfo=””  poissonRatio=”0.1″  youngModulus=”1e+09″  updateStiffnessMatrix=”1″ />
                <TetrahedronSetGeometryAlgorithms template=”Vec3d” name=”tetgeoalgo” />
                <TetrahedronSetTopologyAlgorithms template=”Vec3d” name=”tettopoalgo” />
                <TetrahedronSetTopologyModifier name=”tettopomod” />
                <TetrahedronModel name=”currentTetrahedronModel”  selfCollision=”1″  contactStiffness=”1″  contactFriction=”0.01″  contactRestitution=”0″  color=”1 0.5 0 1″ />
                <FixedConstraint template=”Vec3d” indices=”20 11 12 13″ />
                <UniformMass template=”Vec3d” name=”mass2″  mass=”10000″  totalmass=”210000″  localRange=”11 13″ />
                <UniformMass template=”Vec3d” name=”mass3″  mass=”10000″  totalmass=”210000″  localRange=”20 20″ />
                <Node          name=”tet2tri”  gravity=”0 -10 0″  dt=”0.01″  time=”0″  animate=”0″  >
                        <TriangleSetTopologyContainer name=”Container”  nbPoints=”21″   />
                        <TriangleSetGeometryAlgorithms template=”Vec3d” name=”trigeoalgo” />
                        <TriangleSetTopologyAlgorithms template=”Vec3d” name=”tritopoalgo” />
                        <TriangleSetTopologyModifier name=”tritopomod” />
                        <Tetra2TriangleTopologicalMapping name=”tet2trimap”  input=”@../Container”  output=”@Container” />
                        <TTriangleModel template=”Vec3d” name=”currentTTriangleModel”  selfCollision=”1″  contactStiffness=”1.5″  contactFriction=”0.01″  contactRestitution=”0″  color=”1 0.5 0 1″  bothSide=”1″ />
                </Node>
        </Node>

[chc2]

i forgot to mention that i am using sofa as a library in c++ code.

the part of scn file that i posted is written out for the frame from the c++ code.

i’m not sure what version number it is. this was before sofa was moved to github and before some big fork that was mentioned.

just to make it cleaner i am using:

MechanicalObject
TetrahedronSetTopologyContainer
EulerImplicitSolver
CGLinearSolver
UniformMass
TetrahedralCorotationalFEMForceField
TetrahedronSetGeometryAlgorithms
TetrahedronSetTopologyAlgorithms
TetrahedronSetTopologyModifier
TetrahedronModel
FixedConstraint template=”Vec3d” indices=”20 11 12 13″
UniformMass template=”Vec3d” name=”mass2″ mass=”10000″ totalmass=”210000″ localRange=”11 13″
UniformMass template=”Vec3d” name=”mass3″ mass=”10000″ totalmass=”210000″ localRange=”20 20″

[Hugo]

Hi,

What is the component responsible for the topological change ? TetrahedronModel inheriting from topologicalModifiers ?

Also, why having three UniformMass ?

Did you try the same simulation without the subnode “tet2tri” ?

Cheers,

Hugo

[chc2]

the component responsible for topological change should be TetrahedronSetTopologyModifier

i have 3 uniformmass just so i can set the nodes at the end to be heavy so the bar deforms.

it seems different without tet2tri node. it seems tetras are stretching. maybe it doesnt have the stiffness applied to them?

https://www.youtube.com/watch?v=N87LdmiUD9Y

[Hugo]

Hi chc2,

I am sorry but I don’t get how the scene you pasted above could work: you have no mesh (i.e. no dofs) loaded and the order of the components in the scene looks strange. Is this the entire scene ?

Moreover, the TetrahedronSetTopologyModifier implements some functions for topology but, at my best knowledge, is not responsible for the cutting. For cutting, you need to have some component triggering the topological change, like the TopologicalChangeProcessor in all scenes located in examples/Components/topology/TopologicalModifiers/..

Thanks for the video, it helps in the understanding btw.
Cheers,

Hugo

[chc2]

** just a side note, the forum keeps preventing me from posting. it forbids me then i have to go to cleantalk and delete my IP, but it still doesnt work. do you know what is wrong with the forum? **

that scn file i posted was only to show what components i was using and i didn’t want to make a large post with so many points and tetra information.

i am using sofa as a c++ library.

so i created a small demo that adds a tetra to a tetra mesh. i modified sofa/applications/tutorials/oneTetrahedron.cpp

this creates a *.obj file for each frame. at frame 2 it creates a new tetra. you can see the new point of the tetra not behaving as expected.

you need to change YOURPATH to a path on your system.

this was compiled on ubuntu 16.04 LTS (virtual machine in this case) using sofa 16.12

/******************************************************************************
*       SOFA, Simulation Open-Framework Architecture, development version     *
*                (c) 2006-2017 INRIA, USTL, UJF, CNRS, MGH                    *
*                                                                             *
* This program is free software; you can redistribute it and/or modify it     *
* under the terms of the GNU General Public License as published by the Free  *
* Software Foundation; either version 2 of the License, or (at your option)   *
* any later version.                                                          *
*                                                                             *
* This program is distributed in the hope that it will be useful, but WITHOUT *
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or       *
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for    *
* more details.                                                               *
*                                                                             *
* You should have received a copy of the GNU General Public License along     *
* with this program. If not, see <http://www.gnu.org/licenses/>.              *
*******************************************************************************
* Authors: The SOFA Team and external contributors (see Authors.txt)          *
*                                                                             *
* Contact information: contact@sofa-framework.org                             *
******************************************************************************/

#include <SofaComponentCommon/initComponentCommon.h>
#include <SofaComponentBase/initComponentBase.h>
#include <SofaComponentGeneral/initComponentGeneral.h>
#include <SofaComponentAdvanced/initComponentAdvanced.h>
#include <SofaComponentMisc/initComponentMisc.h>

#include <SofaGraphComponent/Gravity.h>
#include <SofaImplicitOdeSolver/EulerImplicitSolver.h>
#include <SofaBaseLinearSolver/CGLinearSolver.h>
#include <sofa/component/typedef/Sofa_typedef.h>
//#include <SofaOpenglVisual/OglModel.h>
#include <SofaBaseVisual/VisualStyle.h>

#include <sofa/core/objectmodel/Context.h>
#include <sofa/core/VecId.h>

#include <sofa/gui/GUIManager.h>
#include <sofa/gui/Main.h>

#include <sofa/helper/ArgumentParser.h>
#include <sofa/helper/system/FileRepository.h>
#include <sofa/helper/system/glut.h>

#include <SofaSimulationTree/init.h>
#include <SofaSimulationTree/GNode.h>
#include <SofaSimulationTree/TreeSimulation.h>

#include <SofaBaseTopology/TetrahedronSetTopologyContainer.h>
#include <SofaBaseTopology/TetrahedronSetTopologyModifier.h>

#include <iostream>
#include <fstream>

using namespace sofa::simulation::tree;
using namespace sofa;
using sofa::simulation::Node;
using sofa::component::odesolver::EulerImplicitSolver;
typedef component::linearsolver::CGLinearSolver<component::linearsolver::GraphScatteredMatrix, component::linearsolver::GraphScatteredVector> CGLinearSolver;
using sofa::component::topology::MeshTopology;
//using sofa::component::visualmodel::OglModel;
using sofa::core::objectmodel::Data;
using sofa::helper::ReadAccessor;
using sofa::helper::WriteAccessor;
using sofa::core::VecId;
using sofa::core::objectmodel::New;

int main(int argc, char** argv)
{

	sofa::simulation::setSimulation(new sofa::simulation::tree::TreeSimulation());
    sofa::simulation::Node::SPtr root = sofa::simulation::getSimulation()->createNewGraph("root");
    root->setGravity( Coord3(0,-10,0) );

    simulation::Node::SPtr nodeTetra = root->createChild("tetramesh");

    MechanicalObject3::SPtr DOF = sofa::core::objectmodel::New<MechanicalObject3>();
    nodeTetra->addObject(DOF);
    DOF->resize(21);
    DOF->setName("DOF");

    sofa::component::topology::TetrahedronSetTopologyContainer::SPtr topology = sofa::core::objectmodel::New<sofa::component::topology::TetrahedronSetTopologyContainer>();

    WriteAccessor<Data<VecCoord3> > xx = *DOF->write(VecId::position());

    xx[0] = Coord3(0.000000,  0.000000,  0.000000);
    xx[1] = Coord3(0.200000,  0.000000,  0.000000);
    xx[2] = Coord3(0.000000,  0.200000,  0.000000);
    xx[3] = Coord3(0.200000,  0.200000,  0.000000);
    xx[4] = Coord3(0.000000,  0.000000,  1.200000);

    xx[5] = Coord3(0.200000,  0.000000,  1.200000);
    xx[6] = Coord3(0.000000,  0.200000,  1.200000);
    xx[7] = Coord3(0.200000,  0.200000,  1.200000);
    xx[8] = Coord3(0.200000,  0.200000,  0.965324);
    xx[9] = Coord3(0.200000,  0.200000,  0.595933);

    xx[10] = Coord3(0.200000,  0.200000,  0.246792);
    xx[11] = Coord3(0.000000,  0.200000,  0.958366);
    xx[12] = Coord3(0.000000,  0.200000,  0.611182);
    xx[13] = Coord3(0.000000,  0.200000,  0.247715);
    xx[14] = Coord3(0.200000,  0.000000,  0.958366);

    xx[15] = Coord3(0.200000,  0.000000,  0.611182);
    xx[16] = Coord3(0.200000,  0.000000,  0.247715);
    xx[17] = Coord3(0.000000,  0.000000,  1.006005);
    xx[18] = Coord3(0.000000,  0.000000,  0.746736);
    xx[19] = Coord3(0.000000,  0.000000,  0.444605);

    xx[20] = Coord3(0.000000,  0.000000,  0.206183);

    topology->addTetra(17,8,18,14);
    topology->addTetra(12,8,15,18);
    topology->addTetra(18,11,12,8);
    topology->addTetra(20,10,13,2);
    topology->addTetra(6,5,17,4);
    topology->addTetra(5,8,6,17);
    topology->addTetra(20,3,1,10);
    topology->addTetra(12,9,19,15);
    topology->addTetra(15,8,12,9);
    topology->addTetra(18,12,19,15);

    topology->addTetra(9,16,19,15);
    topology->addTetra(18,8,17,11);
    topology->addTetra(20,16,19,10);
    topology->addTetra(20,3,2,0);
    topology->addTetra(16,20,1,10);
    topology->addTetra(18,14,8,15);
    topology->addTetra(19,9,13,10);
    topology->addTetra(14,5,17,8);
    topology->addTetra(6,8,5,7);
    topology->addTetra(20,3,10,2);

    topology->addTetra(12,9,13,19);
    topology->addTetra(11,8,17,6);
    topology->addTetra(3,20,1,0);
    topology->addTetra(20,19,13,10);
    topology->addTetra(19,9,10,16);

	EulerImplicitSolver::SPtr solver = sofa::core::objectmodel::New<EulerImplicitSolver>();
    solver->setName("solver");
    nodeTetra->addObject(solver);

    CGLinearSolver::SPtr linearSolver = New<CGLinearSolver>();
    linearSolver->setName("linearSolver");
    linearSolver->f_maxIter.setValue(2000);
    linearSolver->f_tolerance.setValue(1e-15);
    linearSolver->f_smallDenominatorThreshold.setValue(1e-15);
    nodeTetra->addObject(linearSolver);

    topology->setName("Container");
    topology->setSrc("@"+DOF->getName(), DOF.get());
    nodeTetra->addObject(topology);

	sofa::component::topology::TetrahedronSetTopologyModifier::SPtr tettopomod = sofa::core::objectmodel::New<sofa::component::topology::TetrahedronSetTopologyModifier>();
	tettopomod->setName("tettopomod");
	nodeTetra->addObject(tettopomod);

	TetrahedronSetTopologyAlgorithms3::SPtr tettopoalgo = sofa::core::objectmodel::New<TetrahedronSetTopologyAlgorithms3>();
	tettopoalgo->setName("tettopoalgo");
	nodeTetra->addObject(tettopoalgo);

	TetrahedronSetGeometryAlgorithms3::SPtr tetgeoalgo = sofa::core::objectmodel::New<TetrahedronSetGeometryAlgorithms3>();
	tetgeoalgo->setName("tetgeoalgo");
	nodeTetra->addObject(tetgeoalgo);

    UniformMass3::SPtr mass = sofa::core::objectmodel::New<UniformMass3>();
    nodeTetra->addObject(mass);
    mass->setMass(1);
    mass->setName("mass");

	FixedConstraint3::SPtr fixed1 = sofa::core::objectmodel::New<FixedConstraint3>();
	fixed1->removeConstraint(0);

	fixed1->addConstraint(0);
	fixed1->addConstraint(1);
	fixed1->addConstraint(2);
	fixed1->addConstraint(3);

	nodeTetra->addObject(fixed1);

	TetrahedralCorotationalFEMForceField3::SPtr fem = sofa::core::objectmodel::New<TetrahedralCorotationalFEMForceField3>();
	fem->setName("FEM");
	fem->setSrc("@"+DOF->getName(), DOF.get());
	fem->setMethod(2);
	fem->setUpdateStiffnessMatrix(true);
	fem->setComputeGlobalMatrix(true);
	fem->setYoungModulus(5e+5);
	fem->setPoissonRatio(0.1);
	nodeTetra->addObject(fem);

	sofa::simulation::getSimulation()->init(root.get());

	std::ofstream myfile;
    std::string f1 = "YOURPATH";
    std::string f2 = ".obj";
    string filename = "";

	for(int i=0;i<10;++i)
	{
		sofa::simulation::getSimulation()->animate(root.get(),0.01);

		if(i==2)
		{

			sofa::core::topology::BaseMeshTopology::Tetra tetra = topology->getTetra(18);

            sofa::component::topology::TetrahedronSetTopologyModifier* topoMod;
            topology->getContext()->get(topoMod);

            //add point
            topoMod->addPoints(1,true);

            MechanicalObject3::WriteVecCoord x_pos = DOF->writePositions();
            MechanicalObject3::WriteVecCoord x_restpos = DOF->writeRestPositions();

            int idx = DOF->getSize();

            x_pos[idx-1] = Coord3( topology->getPX(tetra[0]),topology->getPY(tetra[0])+0.2,topology->getPZ(tetra[0]) );
            x_restpos[idx-1] = Coord3( 0.2,1.0,1.2 );

            //add tetra
			helper::vector<core::topology::Topology::Tetrahedron > vitems;
			vitems.resize(1);
			vitems[0][0] = 6;
			vitems[0][1] = 8;
			vitems[0][2] = 7;
			vitems[0][3] = idx-1;
			topoMod->addTetrahedra(vitems);

            topoMod->notifyEndingEvent();
            topoMod->propagateTopologicalChanges();

            sofa::simulation::tree::getSimulation()->init( root.get() );

		}

		std::ostringstream filestring;
		filestring << f1 << i+1 << f2;
		filename = filestring.str();
		myfile.open(filename.c_str());

		for(int a=0;a<DOF->getSize();++a)
		{
			myfile << "v " << topology->getPX(a) << " " << topology->getPY(a) << " " << topology->getPZ(a) << "\n";
		}

		myfile.close();

	}
	
    return 0;
}

i just load the obj files into a 3d application and you can see the point shoot off into space.

here is a different example:

https://www.youtube.com/watch?v=0MMVekCvaAg

i was just trying to think what could be wrong. is there some corresponding rest mesh that i have to update? do i need to update the stiffness? maybe the order of the vertices in the tetra?

do you know what could be the issue?

[Hugo]

Hmmm I am no expert of topological changes..
you can have a look at the following topics:

• Changing topology during runtime
• Proper way of removing points in a Topological Mapping

Hope this helps.

Hugo

[epernod]

Hi chc2,

did you solved your problem? are you mesh still messed up?

Erik

[chc2]

i havent resolved the issue yet. from what i can tell it seems that the actual mesh itself is fine and consistent. maybe it might be some initialization issue in regards to forces?

i am not simply deleting or adding tetras. i am actually reassigning vertices to existing teras and also assigning previously created vertices to new tetras. maybe this can be an issue?

thanks for any input or feedback you may have.

[jnbrunet]

I will make a long shot guess here, but since you are using a dynamic simulation scheme (time integration), by only setting the position of your new vertex, you could be experiencing the consequence of a huge jump in velocity.

I can see two ways of testing this:

1. Manually set the velocity of your new vertex with the value of a nearby point.
2. Use the ancestors and coef parameters of the addPoints function. You choose a set of nearby points of your new vertex as the ancestors and choose some coefficients for each one. I think there is also a way to simply give a tetrahedron as the ancestor and the nearby points will be the element’s vertices and the coefs will be set automatically with the barycentric coordinates. Then, the position and the velocity of your new vertex should be initialized automatically (I think!).

Like I said, this is a wild guess of mine, I may be wrong. But I think it’s worth a try.

Keep us updated.

[Author]

Posts