|
| 1 | +<?xml version="1.0"?> |
| 2 | +<!-- |
| 3 | + 2D cantilever beam under gravity — SOFA equivalent of beam-2d.edp (FreeFEM) |
| 4 | + Two formulations shown side by side: |
| 5 | +
|
| 6 | + In SOFA with the Vec1, Vec2, Vec3 templates. Depending on how these are set up, |
| 7 | + in order to generate a similar outcome as FreeFEM, we need to make some adjustments. |
| 8 | + Just like in many FE codes, we can solve 2D problems in 3D space. This is what using a Vec3 |
| 9 | + template does. The 'd' in Vec3d is the 'double' type of C++. |
| 10 | +
|
| 11 | + Depending on the template we choose, SOFA will do things differently: |
| 12 | +
|
| 13 | + In the Node with name: |
| 14 | + 1. beam_plane_stress (Vec2) y in [0, 2] |
| 15 | + Lame lambda = E*nu / ((1+nu)*(1-nu)) |
| 16 | + DOFs: 2 per node (u_x, u_y). Displacement in the x and y directions. |
| 17 | +
|
| 18 | + 2. beam_plane_strain (Vec3) y in [3, 5] |
| 19 | + Lame lambda = E*nu / ((1+nu)*(1-2*nu)) !!! Notice the 1-2*nu contrary to the above |
| 20 | + DOFs: 3 per node (u_x, u_y, u_z) |
| 21 | +
|
| 22 | + The FreeFEM file (beam-2d.edp) uses 3D Lame coefficients, so formulation 2 is the exact replicate. |
| 23 | + Formulation 1 shows the Vec2d approach. There will be a difference in the results between the |
| 24 | + two. |
| 25 | +
|
| 26 | + Young's modulus: E=21.5 |
| 27 | + Poisson's ratio: nu=0.29 |
| 28 | +--> |
| 29 | + |
| 30 | +<!-- Gravity is applied at the root node --> |
| 31 | +<Node name="root" gravity="0 -0.05 0" dt="1"> |
| 32 | + |
| 33 | + <Node name="plugins"> |
| 34 | + <RequiredPlugin name="Elasticity"/> |
| 35 | + <RequiredPlugin name="Sofa.Component.Constraint.Projective"/> |
| 36 | + <RequiredPlugin name="Sofa.Component.Engine.Select"/> |
| 37 | + <RequiredPlugin name="Sofa.Component.LinearSolver.Direct"/> |
| 38 | + <RequiredPlugin name="Sofa.Component.Mass"/> |
| 39 | + <RequiredPlugin name="Sofa.Component.ODESolver.Backward"/> |
| 40 | + <RequiredPlugin name="Sofa.Component.StateContainer"/> |
| 41 | + <RequiredPlugin name="Sofa.Component.Topology.Container.Dynamic"/> |
| 42 | + <RequiredPlugin name="Sofa.Component.Topology.Container.Grid"/> |
| 43 | + <RequiredPlugin name="Sofa.Component.Visual"/> |
| 44 | + <RequiredPlugin name="Sofa.GL.Component.Rendering3D"/> |
| 45 | + </Node> |
| 46 | + |
| 47 | + <DefaultAnimationLoop/> |
| 48 | + <VisualStyle displayFlags="showBehaviorModels showForceFields showWireframe"/> |
| 49 | + |
| 50 | + <!-- All template='' cases here are Vec2d--> |
| 51 | + <Node name="beam_plane_stress"> |
| 52 | + |
| 53 | + <NewtonRaphsonSolver name="newtonSolver" |
| 54 | + printLog="true" |
| 55 | + warnWhenLineSearchFails="true" |
| 56 | + maxNbIterationsNewton="1" |
| 57 | + maxNbIterationsLineSearch="1" |
| 58 | + lineSearchCoefficient="1" |
| 59 | + relativeSuccessiveStoppingThreshold="0" |
| 60 | + absoluteResidualStoppingThreshold="1e-7" |
| 61 | + absoluteEstimateDifferenceThreshold="1e-12" |
| 62 | + relativeInitialStoppingThreshold="1e-12" |
| 63 | + relativeEstimateDifferenceThreshold="0" |
| 64 | + /> |
| 65 | + <SparseLDLSolver name="linearSolver" template="CompressedRowSparseMatrixd"/> |
| 66 | + <StaticSolver name="staticSolver" newtonSolver="@newtonSolver" linearSolver="@linearSolver"/> |
| 67 | + |
| 68 | + <RegularGridTopology name="grid" min="0 0 0" max="10 2 0" n="41 11 1"/> |
| 69 | + <MechanicalObject template="Vec2d" name="dofs"/> |
| 70 | + |
| 71 | + <Node name="triangles"> |
| 72 | + <TriangleSetTopologyContainer name="topology" src="@../grid"/> |
| 73 | + |
| 74 | + <MeshMatrixMass template="Vec2d,Vec2d" totalMass="20" topology="@topology"/> |
| 75 | + <LinearSmallStrainFEMForceField name="FEM" template="Vec2d" |
| 76 | + youngModulus="21.5" poissonRatio="0.29" |
| 77 | + topology="@topology" |
| 78 | + /> |
| 79 | + </Node> |
| 80 | + |
| 81 | + <BoxROI name="fixed_roi" template="Vec2d" |
| 82 | + box="-0.01 -0.01 -0.01 0.01 2.01 0.01" drawBoxes="true"/> |
| 83 | + <FixedProjectiveConstraint template="Vec2d" indices="@fixed_roi.indices"/> |
| 84 | + |
| 85 | + </Node> |
| 86 | + |
| 87 | + <!-- All template='' cases here are Vec3d--> |
| 88 | + <Node name="beam_plane_strain"> |
| 89 | + |
| 90 | + <NewtonRaphsonSolver name="newtonSolver" |
| 91 | + printLog="true" |
| 92 | + warnWhenLineSearchFails="true" |
| 93 | + maxNbIterationsNewton="1" |
| 94 | + maxNbIterationsLineSearch="1" |
| 95 | + lineSearchCoefficient="1" |
| 96 | + relativeSuccessiveStoppingThreshold="0" |
| 97 | + absoluteResidualStoppingThreshold="1e-7" |
| 98 | + absoluteEstimateDifferenceThreshold="1e-12" |
| 99 | + relativeInitialStoppingThreshold="1e-12" |
| 100 | + relativeEstimateDifferenceThreshold="0" |
| 101 | + /> |
| 102 | + <SparseLDLSolver name="linearSolver" template="CompressedRowSparseMatrixd"/> |
| 103 | + <StaticSolver name="staticSolver" newtonSolver="@newtonSolver" linearSolver="@linearSolver"/> |
| 104 | + |
| 105 | + <RegularGridTopology name="grid" min="0 3 0" max="10 5 0" n="41 11 1"/> |
| 106 | + <MechanicalObject template="Vec3d" name="dofs"/> |
| 107 | + |
| 108 | + <Node name="triangles"> |
| 109 | + <TriangleSetTopologyContainer name="topology" src="@../grid"/> |
| 110 | + <TriangleSetTopologyModifier/> |
| 111 | + <TriangleSetGeometryAlgorithms template="Vec3d" drawTriangles="true"/> |
| 112 | + |
| 113 | + <MeshMatrixMass totalMass="20" topology="@topology"/> |
| 114 | + <LinearSmallStrainFEMForceField name="FEM" |
| 115 | + youngModulus="21.5" poissonRatio="0.29" |
| 116 | + topology="@topology" |
| 117 | + mstate="@../dofs"/> |
| 118 | + </Node> |
| 119 | + |
| 120 | + <!-- Clamp left edge: x=0, y in [3,5]. |
| 121 | + FixedProjectiveConstraint fixes all 3 DOFs (u_x, u_y, u_z) |
| 122 | + of the left-edge nodes, which is all that is needed. --> |
| 123 | + <BoxROI name="fixed_roi" template="Vec3d" |
| 124 | + box="-0.01 2.99 -1.0 0.01 5.01 1.0" drawBoxes="true"/> |
| 125 | + <FixedProjectiveConstraint template="Vec3d" indices="@fixed_roi.indices"/> |
| 126 | + |
| 127 | + </Node> |
| 128 | + |
| 129 | +</Node> |
0 commit comments