Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

[all] Clean due to doc #1398

Merged
merged 4 commits into from
Jun 3, 2020
Merged

[all] Clean due to doc #1398

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
7a3249b
Clean cpp code and release full example scene
hugtalbot May 27, 2020
1a37d6a
Add an example scene for the PointsFromIndices
hugtalbot May 27, 2020
50016d6
Fix compilation in test Newmark
hugtalbot May 29, 2020
d3282c7
Remove name of NewmarkImplicitSolver in msg
hugtalbot Jun 2, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
28 changes: 28 additions & 0 deletions examples/Components/engine/PointsFromIndices.scn
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,28 @@
<?xml version="1.0" ?>
<Node name="root" gravity="0 -9 1" dt="0.05">
<RequiredPlugin name="SofaOpenglVisual"/>
<VisualStyle displayFlags="showForceFields" />

<Node name="SquareGravity" gravity="0 -9.81 0">
<EulerImplicitSolver name="cg_odesolver" printLog="false" rayleighStiffness="0.1" rayleighMass="0.1" />
<CGLinearSolver iterations="25" name="linear solver" tolerance="1.0e-9" threshold="1.0e-9" />
<MeshGmshLoader name="loader" filename="mesh/square3.msh" createSubelements="true"/>
<MechanicalObject src="@loader" template="Vec3d" name="mecaObj" scale3d="10 10 10" restScale="1" />
<TriangleSetTopologyContainer src="@loader" name="Container" />
<TriangleSetTopologyModifier name="Modifier" />
<TriangleSetTopologyAlgorithms template="Vec3d" name="TopoAlgo" />
<TriangleSetGeometryAlgorithms template="Vec3d" name="GeomAlgo" />
<DiagonalMass template="Vec3d" name="default5" massDensity="0.15" />
<BoxROI template="Vec3d" name="FixedROI" box="2 9.5 -0.5 8 10.5 0.5" drawBoxes="1" position="@mecaObj.rest_position" computeTriangles="0" computeTetrahedra="0" computeEdges="0" />
<FixedConstraint template="Vec3d" name="default6" indices="@FixedROI.indices" />
<TriangularFEMForceField template="Vec3d" name="FEM" method="large" poissonRatio="0.3" youngModulus="60" />
<TriangularBendingSprings template="Vec3d" name="FEM-Bend" stiffness="300" damping="1" />

<PointsFromIndices name="PFI" position="@mecaObj.position" indices="10 20 30" />

<Node name="Selection" >
<MechanicalObject template="Vec3d" position="@../PFI.indices_position" name="SelectedDOFs" showIndices="1" showIndicesScale="0.2" />
</Node>

</Node>
</Node>
65 changes: 34 additions & 31 deletions examples/Components/solver/NewmarkImplicitSolver.scn
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,41 +14,44 @@
<UniformMass totalMass="15" />
<FixedConstraint indices="0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 268 269 270 271 343 345" />
<MeshSpringForceField name="Spring" tetrasStiffness="1870" tetrasDamping="0" />

<Node>
<MeshObjLoader name="meshLoader_3" filename="mesh/truthcylinder1.obj" handleSeams="1" />
<OglModel name="Visual" src="@meshLoader_3" color="yellow" dx="15" />
<BarycentricMapping input="@.." output="@Visual" />
</Node>
</Node>
<!--
<Node name="CoFEM">
<NewmarkImplicitSolver />
<CGLinearSolver iterations="100" />
<MeshGmshLoader name="loader" filename="mesh/truthcylinder1.msh"/>
<MeshTopology src="@loader" />
<MechanicalObject src="@loader" dx="30" />
<UniformMass totalmass="15" />
<FixedConstraint indices="0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 268 269 270 271 343 345" />
<TetrahedronFEMForceField name="FEM" youngModulus="1116" poissonRatio="0.49" method="polar" />
<Node>
<MeshObjLoader name="meshLoader_2" filename="mesh/truthcylinder1.obj" handleSeams="1" />
<OglModel name="Visual" src="@meshLoader_2" color="cyan" dx="30" />
<BarycentricMapping input="@.." output="@Visual" />
</Node>
</Node>
<Node name="LinearFEM">
<NewmarkImplicitSolver />
<CGLinearSolver iterations="100" />
<MeshGmshLoader name="loader" filename="mesh/truthcylinder1.msh"/>
<MeshTopology src="@loader" />
<MechanicalObject src="@loader" dx="45" />
<UniformMass totalmass="15" />
<FixedConstraint indices="0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 268 269 270 271 343 345" />
<TetrahedronFEMForceField name="FEM" youngModulus="1116" poissonRatio="0.49" method="small" />
<Node>
<MeshObjLoader name="meshLoader_1" filename="mesh/truthcylinder1.obj" handleSeams="1" />
<OglModel name="Visual" src="@meshLoader_1" color="red" dx="45" />
<BarycentricMapping input="@.." output="@Visual" />
</Node>
</Node>-->

<Node name="CoFEM">
<NewmarkImplicitSolver />
<CGLinearSolver iterations="100" />
<MeshGmshLoader name="loader" filename="mesh/truthcylinder1.msh"/>
<MeshTopology src="@loader" />
<MechanicalObject src="@loader" dx="30" />
<UniformMass totalMass="15" />
<FixedConstraint indices="0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 268 269 270 271 343 345" />
<TetrahedronFEMForceField name="FEM" youngModulus="1116" poissonRatio="0.49" method="polar" />

<Node>
<MeshObjLoader name="meshLoader_2" filename="mesh/truthcylinder1.obj" handleSeams="1" />
<OglModel name="Visual" src="@meshLoader_2" color="cyan" dx="30" />
<BarycentricMapping input="@.." output="@Visual" />
</Node>
</Node>
<Node name="LinearFEM">
<NewmarkImplicitSolver />
<CGLinearSolver iterations="100" />
<MeshGmshLoader name="loader" filename="mesh/truthcylinder1.msh"/>
<MeshTopology src="@loader" />
<MechanicalObject src="@loader" dx="45" />
<UniformMass totalMass="15" />
<FixedConstraint indices="0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 268 269 270 271 343 345" />
<TetrahedronFEMForceField name="FEM" youngModulus="1116" poissonRatio="0.49" method="small" />

<Node>
<MeshObjLoader name="meshLoader_1" filename="mesh/truthcylinder1.obj" handleSeams="1" />
<OglModel name="Visual" src="@meshLoader_1" color="red" dx="45" />
<BarycentricMapping input="@.." output="@Visual" />
</Node>
</Node>
</Node>
105 changes: 31 additions & 74 deletions modules/SofaMiscSolver/NewmarkImplicitSolver.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -45,12 +45,11 @@ using namespace sofa::defaulttype;
using namespace core::behavior;

NewmarkImplicitSolver::NewmarkImplicitSolver()
: f_rayleighStiffness(initData(&f_rayleighStiffness,0.0,"rayleighStiffness","Rayleigh damping coefficient related to stiffness") )
, f_rayleighMass( initData(&f_rayleighMass,0.0,"rayleighMass","Rayleigh damping coefficient related to mass"))
, f_velocityDamping( initData(&f_velocityDamping,0.0,"vdamping","Velocity decay coefficient (no decay if null)") )
, f_verbose( initData(&f_verbose,false,"verbose","Dump system state at each iteration") )
, f_gamma( initData(&f_gamma, 0.5, "gamma", "Newmark scheme gamma coefficient"))
, f_beta( initData(&f_beta, 0.25, "beta", "Newmark scheme beta coefficient") )
: d_rayleighStiffness(initData(&d_rayleighStiffness,0.0,"rayleighStiffness","Rayleigh damping coefficient related to stiffness") )
, d_rayleighMass( initData(&d_rayleighMass,0.0,"rayleighMass","Rayleigh damping coefficient related to mass"))
, d_velocityDamping( initData(&d_velocityDamping,0.0,"vdamping","Velocity decay coefficient (no decay if null)") )
, d_gamma( initData(&d_gamma, 0.5, "gamma", "Newmark scheme gamma coefficient"))
, d_beta( initData(&d_beta, 0.25, "beta", "Newmark scheme beta coefficient") )
, d_threadSafeVisitor(initData(&d_threadSafeVisitor, false, "threadSafeVisitor", "If true, do not use realloc and free visitors in fwdInteractionForceField."))
{
cpt=0;
Expand All @@ -63,7 +62,6 @@ void NewmarkImplicitSolver::solve(const core::ExecParams* params, SReal dt, sofa
sofa::simulation::common::MechanicalOperations mop( params, this->getContext() );
MultiVecCoord pos(&vop, core::VecCoordId::position() );
MultiVecDeriv vel(&vop, core::VecDerivId::velocity() );
// MultiVecDeriv f(&vop, core::VecDerivId::force() );
MultiVecDeriv b(&vop);
MultiVecDeriv aResult(&vop);
MultiVecCoord newPos(&vop, xResult );
Expand All @@ -73,63 +71,36 @@ void NewmarkImplicitSolver::solve(const core::ExecParams* params, SReal dt, sofa
// dx is no longer allocated by default (but it will be deleted automatically by the mechanical objects)
MultiVecDeriv dx(&vop, core::VecDerivId::dx()); dx.realloc(&vop, !d_threadSafeVisitor.getValue(), true);


const SReal h = dt;
const double gamma = f_gamma.getValue();
const double beta = f_beta.getValue();
const double rM = f_rayleighMass.getValue();
const double rK = f_rayleighStiffness.getValue();
const bool verbose = f_verbose.getValue();

/* This integration scheme is based on the following equations:
*
* $x_{t+h} = x_t + h v_t + h^2/2 ( (1-2\beta) a_t + 2\beta a_{t+h} )$
* $v_{t+h} = v_t + h ( (1-\gamma) a_t + \gamma a_{t+h} )$
*
* Applied to a mechanical system where $ M a_t + (r_M M + r_K K) v_t + K x_t
= f_ext$, we need to solve the following system:
*
* $ M a_{t+h} + (r_M M + r_K K) v_{t+h} + K x_{t+h} = f_ext $
* $ M a_{t+h} + (r_M M + r_K K) ( v_t + h ( (1-\gamma) a_t + \gamma a_{t+h} ) ) + K ( x_t + h v_t + h^2/2 ( (1-2\beta) a_t + 2\beta a_{t+h} ) ) =
f_ext $
* $ ( M + h \gamma (r_M M + r_K K) + h^2 \beta K ) a_{t+h} = f_ext - (r_M
M + r_K K) ( v_t + h (1-\gamma) a_t ) - K ( x_t + h v_t + h^2/2 (1-2\beta) a_t )
$
* $ ( (1 + h \gamma r_M) M + (h^2 \beta + h \gamma r_K) K ) a_{t+h} =
f_ext - (r_M M + r_K K) v_t - K x_t - (r_M M + r_K K) ( h (1-\gamma) a_t ) - K (
h v_t + h^2/2 (1-2\beta) a_t) $
* $ ( (1 + h \gamma r_M) M + (h^2 \beta + h \gamma r_K) K ) a_{t+h} = Ma_t
- (r_M M + r_K K) ( h (1-\gamma) a_t ) - K ( h v_t + h^2/2 (1-2\beta) a_t ) $
*
* The current implementation first computes $a_t$ directly (as in the
explicit solvers), then solves the previous system to compute $a_{t+dt}$, and
finally computes the new position and velocity.
*/

//we need to initialize a_t and to store it as a vecId to be used in the resolution of this solver (using as well old xand v). Once we have a_{t+dt} we
//can update the new x and v.
const double gamma = d_gamma.getValue();
const double beta = d_beta.getValue();
const double rM = d_rayleighMass.getValue();
const double rK = d_rayleighStiffness.getValue();

// 1. Initialize a_t and to store it as a vecId to be used in the resolution of this solver (using as well old xand v)
// Once we have a_{t+dt} we can update the new x and v.

if (cpt == 0 || this->getContext()->getTime()==0.0)
{
vop.v_alloc(pID);

}

// Define a
MultiVecDeriv a(&vop, pID);
a.realloc(&vop, !d_threadSafeVisitor.getValue(), true);
if(cpt ==0)

if(cpt == 0)
{
a.clear();
mop.computeAcc(0,a,pos,vel);
}
cpt++;

if( verbose )
{
msg_info() << "NewmarkImplicitSolver, aPrevious = " << a;
msg_info() << "NewmarkImplicitSolver, xPrevious = " << pos;
msg_info() << "NewmarkImplicitSolver, vPrevious = " << vel;
}
msg_info() << "aPrevious = " << a;
msg_info() << "xPrevious = " << pos;
msg_info() << "vPrevious = " << vel;


// 2. Compute right hand term of equation on a_{t+h}

Expand All @@ -140,37 +111,27 @@ void NewmarkImplicitSolver::solve(const core::ExecParams* params, SReal dt, sofa
{
mop.propagateDx(a);

mop.addMBKdx(b, -h*(1-gamma)*rM, h*(1-gamma), h*(1-gamma)*rK + h*h*(1-2*beta)/2.0,true,true);
// b += (-h (1-\gamma)(r_M M + r_K K) - h^2/2 (1-2\beta) K ) a

mop.addMBKdx(b, -h*(1-gamma)*rM, h*(1-gamma), h*(1-gamma)*rK + h*h*(1-2*beta)/2.0,true,true);
}

mop.addMBKv(b, -rM, 1,rK+h);
// b += -h K v
mop.addMBKv(b, -rM, 1,rK+h);

if (verbose)
msg_info() << "NewmarkImplicitSolver, b = " << b;

msg_info() << "b = " << b;
mop.projectResponse(b); // b is projected to the constrained space
msg_info() << "projected b = " << b;

if (verbose)
msg_info() << "NewmarkImplicitSolver, projected b = " << b;

// 3. Solve system of equations on a_{t+h}

core::behavior::MultiMatrix<simulation::common::MechanicalOperations> matrix(&mop);

matrix = MechanicalMatrix::K * (-h*h*beta - h*rK*gamma) + MechanicalMatrix::B*(-h)*gamma + MechanicalMatrix::M * (1 + h*gamma*rM);

// if( verbose )
// msg_info()<<"NewmarkImplicitSolver, matrix = "<< MechanicalMatrix::K *(-h*h*beta + -h*rK*gamma) + MechanicalMatrix::M * (1 + h*gamma*rM) << " = " << matrix<<sendl;


msg_info()<<"matrix = "<< MechanicalMatrix::K *(-h*h*beta + -h*rK*gamma) + MechanicalMatrix::M * (1 + h*gamma*rM) << " = " << matrix<<sendl;
matrix.solve(aResult, b);
//mop.projectResponse(aResult);

if (verbose)
msg_info() << "NewmarkImplicitSolver, a1 = " << aResult;
msg_info() << "a1 = " << aResult;


// 4. Compute the new position and velocity
Expand Down Expand Up @@ -208,23 +169,19 @@ void NewmarkImplicitSolver::solve(const core::ExecParams* params, SReal dt, sofa
#endif

mop.addSeparateGravity(dt, newVel); // v += dt*g . Used if mass wants to add G separately from the other forces to v.
if (f_velocityDamping.getValue()!=0.0)
newVel *= exp(-h*f_velocityDamping.getValue());
if (d_velocityDamping.getValue()!=0.0)
newVel *= exp(-h*d_velocityDamping.getValue());

if( verbose )
{
msg_info() << "NewmarkImplicitSolver, final x = " << newPos;
msg_info() << "NewmarkImplicitSolver, final v = " << newVel;
}

msg_info() << "final x = " << newPos;
msg_info() << "final v = " << newVel;

// Increment
a.eq(aResult);
}

int NewmarkImplicitSolverClass = core::RegisterObject("Implicit time integratorusing Newmark scheme")
.add< NewmarkImplicitSolver >()
.addAlias("Newmark");
;
.add< NewmarkImplicitSolver >();

} // namespace odesolver

Expand Down
11 changes: 5 additions & 6 deletions modules/SofaMiscSolver/NewmarkImplicitSolver.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -62,13 +62,12 @@ class SOFA_MISC_SOLVER_API NewmarkImplicitSolver : public sofa::core::behavior::

public:
SOFA_CLASS(NewmarkImplicitSolver, sofa::core::behavior::OdeSolver);
Data<double> f_rayleighStiffness; ///< Rayleigh damping coefficient related to stiffness
Data<double> f_rayleighMass; ///< Rayleigh damping coefficient related to mass
Data<double> f_velocityDamping; ///< Velocity decay coefficient (no decay if null)
Data<bool> f_verbose; ///< Dump system state at each iteration
Data<double> d_rayleighStiffness; ///< Rayleigh damping coefficient related to stiffness
Data<double> d_rayleighMass; ///< Rayleigh damping coefficient related to mass
Data<double> d_velocityDamping; ///< Velocity decay coefficient (no decay if null)

Data<double> f_gamma; ///< Newmark scheme gamma coefficient
Data<double> f_beta; ///< Newmark scheme beta coefficient
Data<double> d_gamma; ///< Newmark scheme gamma coefficient
Data<double> d_beta; ///< Newmark scheme beta coefficient

Data<bool> d_threadSafeVisitor;

Expand Down
4 changes: 2 additions & 2 deletions modules/SofaMiscSolver/SofaMiscSolver_test/NewmarkImplicitSolverDynamic_test.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -88,8 +88,8 @@ struct NewmarkImplicitDynamic_test : public Elasticity_test<_DataTypes>

// Solver
NewmarkImplicitSolver::SPtr newmarkSolver = addNew<NewmarkImplicitSolver> (root);
newmarkSolver->f_rayleighStiffness.setValue(rk);
newmarkSolver->f_rayleighMass.setValue(rm);
newmarkSolver->d_rayleighStiffness.setValue(rk);
newmarkSolver->d_rayleighMass.setValue(rm);

CGLinearSolver::SPtr cgLinearSolver = addNew<CGLinearSolver> (root);
cgLinearSolver->f_maxIter=3000;
Expand Down