diff --git a/examples/Freefem/edp/3d-cube.edp b/examples/Freefem/edp/3d-cube.edp
new file mode 100644
index 0000000..74916a8
--- /dev/null
+++ b/examples/Freefem/edp/3d-cube.edp
@@ -0,0 +1,78 @@
+load "msh3"
+load "medit"
+load "iovtk"
+
+
+// Pb elasticite dans un cube  
+mesh3 Th = cube(10, 10, 10);
+
+// Vérifier les labels du cube
+cout << "Labels bords : " << endl;
+int[int] labs = labels(Th);
+for(int i = 0; i < labs.n; i++)
+    cout << "  label " << labs[i] << endl;
+
+real E       = 210e9;
+real sigma   = 0.3;
+real rho     = 7800.;
+real gravity = -9.81 * rho;
+real mu      = E / (2.*(1.+sigma));
+real lambda  = E*sigma / ((1.+sigma)*(1.-2.*sigma));
+
+// Zmin/zmax
+real zmin = 1e100, zmax = -1e100;
+for (int i = 0; i < Th.nv; i++) {
+    zmin = min(zmin, Th(i).z);
+    zmax = max(zmax, Th(i).z);
+}
+cout << "zmin=" << zmin << "  zmax=" << zmax << endl;
+
+fespace Vh(Th, [P1,P1,P1]);
+Vh [ux,uy,uz],[vx,vy,vz];
+
+macro Epsilon(u1,u2,u3) [dx(u1),dy(u2),dz(u3),(dy(u3)+dz(u2))/sqrt(2.),(dz(u1)+dx(u3))/sqrt(2.),(dx(u2)+dy(u1))/sqrt(2.)] // EOM
+macro Div(u1,u2,u3) (dx(u1)+dy(u2)+dz(u3)) // EOM
+ 
+
+problem Elasticite3D([ux,uy,uz],[vx,vy,vz], solver=sparsesolver) =
+    int3d(Th)(
+        lambda * Div(vx,vy,vz) * Div(ux,uy,uz)
+      + 2.*mu * (Epsilon(vx,vy,vz)' * Epsilon(ux,uy,uz))
+    )
+    - int3d(Th)(gravity * vz)
+    + on(1, ux=0, uy=0, uz=0);  // encastrement 
+
+Elasticite3D;
+
+cout << "Deplacement max uz = " << uz[].linfty << " m" << endl;
+cout << "Theorique uz_max   = " << rho*9.81*(zmax-zmin)^2/(2.*E) << " m" << endl;
+
+fespace Wh(Th, P1);
+Wh normu = sqrt(ux^2 + uy^2 + uz^2);
+medit("Norme", Th, normu, wait=1);
+medit("uz", Th, uz, wait=1);
+
+real umax = max(ux[].linfty, max(uy[].linfty, uz[].linfty));
+real coef = 0.1*(zmax-zmin)/umax;
+cout << "Coefficient amplification = " << coef << endl;
+mesh3 Thd = movemesh3(Th, transfo=[x+coef*ux, y+coef*uy, z+coef*uz]);
+//savemesh(Thd, "cube_deforme.mesh");
+medit("Cube deforme", Thd, wait=1);
+
+
+
+
+
+
+
+// Exporter la solution en format VTK
+int[int] order = [1, 1, 1]; // Ordre P1
+savevtk("cube_deplacement.vtu", Th, [ux, uy, uz], dataname="Deplacement", order=order);
+savevtk("cube_norme.vtu", Th, normu, dataname="Norme", order=order);
+savevtk("cube_uz.vtu", Th, uz, dataname="Uz", order=order);
+
+// Exporter le maillage déformé
+mesh3 Thd = movemesh3(Th, transfo=[x + 1e6*ux, y + 1e6*uy, z + 1e6*uz]);
+savevtk("cube_deforme.vtu", Thd, [ux, uy, uz], dataname="Deplacement", order=order);
+
+cout << "Fichiers VTK exportés. Ouvrez-les avec ParaView." << endl;
\ No newline at end of file
diff --git a/examples/Freefem/edp/barre-1d-solution-man.edp b/examples/Freefem/edp/barre-1d-solution-man.edp
new file mode 100644
index 0000000..406c39f
--- /dev/null
+++ b/examples/Freefem/edp/barre-1d-solution-man.edp
@@ -0,0 +1,44 @@
+// Solution manufacturée normalisée
+// On résout : -u'' = pi^2*sin(pi*x)  EA simplifie =1
+real L = 1.0;
+real E = 1.0;   
+real A = 1.0;   
+real pi = 3.14159265358979;
+
+int[int] nvals = [2, 4, 8, 16, 32, 64];
+
+ofstream file("convergence_grossier.txt");
+file << "# n\th\terrL2\terrH1" << endl;
+
+for (int idx = 0; idx < nvals.n; idx++) {
+    int n = nvals[idx];
+    real h = L / n;
+
+    meshL Th = segment(n, [x * L]);
+    fespace Vh(Th, P1);
+    Vh u, v;
+
+    problem Traction(u, v) =
+        int1d(Th)( dx(u) * dx(v) )
+        - int1d(Th)( pi * pi * sin(pi * x) * v )
+        + on(1, u = 0)
+        + on(2, u = 0);
+
+    Traction;
+
+    Vh uexact = sin(pi * x);
+
+    // Erreur sur les noeuds uniquement
+real errL2 = 0, errH1 = 0;
+
+errL2 = sqrt(int1d(Th)( (u - sin(pi*x))^2 ));
+errH1 = sqrt(int1d(Th)( (dx(u) - pi*cos(pi*x))^2 ));
+
+
+    cout << "n = " << n
+         << ",  h = " << h
+         << ",  errL2 = " << errL2
+         << ",  errH1 = " << errH1 << endl;
+
+    file << n << "\t" << h << "\t" << errL2 << "\t" << errH1 << endl;
+}
\ No newline at end of file
diff --git a/examples/Freefem/edp/barre_exacte.edp b/examples/Freefem/edp/barre_exacte.edp
new file mode 100644
index 0000000..55adbef
--- /dev/null
+++ b/examples/Freefem/edp/barre_exacte.edp
@@ -0,0 +1,42 @@
+// === Élasticité linéaire - BARRE en traction ===
+// Version pour étude de convergence
+
+real L = 1.0;
+real E = 1e6;
+real A = 0.001;
+real F = 1000;
+
+// Différentes valeurs de n pour l'étude de convergence
+
+int[int] nvals = [2, 4, 8, 16, 32, 64];
+ofstream file("convergence_data.txt");
+file << "# n\th\terrL2\terrH1" << endl;
+
+for (int idx = 0; idx < nvals.n; idx++) {
+    int n = nvals[idx];
+    real h = L/n;
+    
+    meshL Th = segment(n, [x*L]);
+    fespace Vh(Th, P1);
+    Vh u, v;
+    
+    problem Traction(u, v) =
+        int1d(Th)( E * A * dx(u) * dx(v) )
+        - int0d(Th, 2)( F * v )
+        + on(1, u=0);
+    
+    Traction;
+    
+    Vh uexact = (F/(E*A)) * x;
+    
+    // Erreur L2
+    real errL2 = sqrt(int1d(Th)((u - uexact)^2));
+    
+    // Erreur H1 (semi-norme)
+    real errH1 = sqrt(int1d(Th)((dx(u) - dx(uexact))^2));
+    
+    cout << "n = " << n << ", h = " << h << ", errL2 = " << errL2 << ", errH1 = " << errH1 << endl;
+    
+    // Sauvegarde
+    file << n << "\t " << h << "\t " << errL2 << "\t " << errH1 << endl;
+}
\ No newline at end of file
diff --git a/examples/Freefem/edp/beam-2d.edp b/examples/Freefem/edp/beam-2d.edp
new file mode 100644
index 0000000..2f9269e
--- /dev/null
+++ b/examples/Freefem/edp/beam-2d.edp
@@ -0,0 +1,104 @@
+// Beam under its own weight 
+
+// Parameters
+real E = 21.5;
+real sigma = 0.29;
+real gravity = -0.05;
+int n=140;
+int m=35;
+
+
+// Maillage non structure " avec courbes parametriques "
+border a(t=2, 0){x=0; y=t ;label=1;}      // gauche cote de la poutre à encastre 
+border b(t=0, 10){x=t; y=0 ;label=2;}     // bas (libre)
+border c(t=0, 2){x=10; y=t ;label=3;}     // droite (libre)
+border d(t=0, 10){x=10-t; y=2; label=4;}  // haut (libre)
+mesh th = buildmesh(b(n) + c(m) + d(n) + a(m)); // un maillage equilibre previlige celui là 
+// n=100 unite de long pour les cotes haut et bas, et n=25 unite de haut pour les cote gauche et droite 
+//mesh th = buildmesh(b(n) + c(n) + d(n) + a(n)); // evite celui là le maillage est raffine uniquement aux extrimite a et c aka gauche et droite 
+// Fespace element finis espace P1 ;
+plot (th,cmm="maillage-poutre", wait=1);
+fespace Vh(th, [P1, P1]);
+Vh [uu, vv], [w, s];
+
+// Macro
+real sqrt2 = sqrt(2.);
+macro epsilon(u1, u2) [dx(u1), dy(u2), (dy(u1) + dx(u2))/sqrt2] // EOM
+macro div(u, v) (dx(u) + dy(v)) // EOM
+
+// Coefficients de Lamé
+real mu = E/(2*(1 + sigma));
+real lambda = E*sigma/((1 + sigma)*(1 - 2*sigma));
+
+// Problème
+solve Elasticity ([uu, vv], [w, s], solver=sparsesolver)
+    = int2d(th)(
+          lambda*div(w, s)*div(uu, vv)
+        + 2.*mu*(epsilon(w, s)' * epsilon(uu, vv))
+    )
+    - int2d(th)(gravity*s)
+    + on(1, uu=0, vv=0);  
+
+// Résultats
+cout << "Max displacement = " << uu[].linfty << ", " << vv[].linfty << endl;
+
+// Visualisation
+plot([uu, vv], wait=1, cmm="Deplacement");
+
+
+mesh th1 = movemesh(th, [x + uu, y + vv]);
+plot( th1, wait=1, cmm="Maillage initial  et deformation");
+
+// Espace pour les contraintes
+fespace Wh(th, P1);
+Wh sigmaxx, sigmayy, sigmaxy, sigmavm;
+
+// Calcul des contraintes
+sigmaxx = lambda * (dx(uu) + dy(vv)) + 2*mu*dx(uu);
+sigmayy = lambda * (dx(uu) + dy(vv)) + 2*mu*dy(vv);
+sigmaxy = mu*(dy(uu) + dx(vv));
+
+// Contrainte de von Mises pour plasticité
+sigmavm = sqrt(sigmaxx^2 + sigmayy^2 - sigmaxx*sigmayy + 3*sigmaxy^2);
+
+// Visualisation
+plot(sigmaxx, wait=1, fill=1, value=1, cmm="Contrainte sigma_xx (MPa)");
+plot(sigmavm, wait=1, fill=1, value=1, cmm="Contrainte von Mises (MPa)");
+
+// Valeurs extrêmes
+cout << "sigma_xx min = " << sigmaxx[].min << ", max = " << sigmaxx[].max << endl;
+cout << "sigma_vm max = " << sigmavm[].max << endl;
+
+
+
+//  Calcul des réactions par intégration sur le bord encastré (label=1)
+real reactionx = int1d(th, 1)( (sigmaxx * N.x + sigmaxy * N.y) );
+real reactiony = int1d(th, 1)( (sigmaxy * N.x + sigmayy * N.y) );
+
+//  Poids total 
+real poidstotal = abs(gravity) * int2d(th)(1);
+
+//  Affichage des résultats
+
+cout << "Reaction horizontale a l'encastrement : " << reactionx << endl;
+cout << "Reaction verticale a l'encastrement   : " << reactiony << endl;
+cout << "Poids total de la poutre               : " << poidstotal << endl;
+cout << "Rapport reaction verticale / poids      : " << reactiony / poidstotal << endl;
+
+// Test d'équilibre
+real erreurequilibre = abs(reactiony/poidstotal - 1) * 100;
+cout << "Erreur d'equilibre : " << erreurequilibre << " %" << endl;
+
+if (erreurequilibre < 5)
+    cout << "equilibre verifie " << endl;
+else
+    cout << " Attention : desequilibre detecte !" << endl;
+
+//  Vérification du moment à l'encastrement
+real moment = int1d(th, 1)( (y - 1) * (sigmaxy * N.x + sigmayy * N.y) - (x) * (sigmaxx * N.x + sigmaxy * N.y) );
+cout << "Moment a l'encastrement : " << moment << endl;
+
+
+
+
+
diff --git a/examples/Freefem/edp/elasticity_3D.edp b/examples/Freefem/edp/elasticity_3D.edp
new file mode 100644
index 0000000..902cc75
--- /dev/null
+++ b/examples/Freefem/edp/elasticity_3D.edp
@@ -0,0 +1,105 @@
+load "msh3"
+load "medit"
+load "iovtk"
+ // remarque : la comparaison theorique qu'on a obtenue n'est pas vraiment valide 
+ // car on compare une geometrie tres difficle à une poutre, mais on reste tjrs dans le cadre du respect de la loi d la physique 
+mesh3 Th("C:/Program Files (x86)/FreeFem++/examples/3d/dodecaedre01.mesh");
+
+// Calculer zmin
+real zmin = 1e100, zmax = -1e100;
+for (int i = 0; i < Th.nv; i++) {
+    zmin = min(zmin, Th(i).z);
+    zmax = max(zmax, Th(i).z);
+}
+cout << "zmin=" << zmin << "  zmax=" << zmax << endl;
+
+// Changer le label des faces du bas (proches de zmin)
+Th = change(Th, flabel=(abs(z - zmin) < 0.1) ? 1 : 0);
+
+// Vérifier les labels modifiés
+cout << "Nouveaux labels  :" << endl;
+for (int i = 0; i < Th.nbe; i++) {
+    if (Th.be(i).label == 1)
+        cout << "Face " << i << " label = " << Th.be(i).label << endl;
+}
+
+// Vérifier tous les labels présents
+int[int] labs = labels(Th);
+cout << "Labels presents dans le maillage :" << endl;
+for(int i = 0; i < labs.n; i++)
+    cout << "  label " << labs[i] << endl;
+
+// Paramètres physiques
+real E       = 210e9;
+real sigma   = 0.3;
+real rho     = 7800.;
+real gravity = -9.81 * rho;
+real mu      = E / (2.*(1.+sigma));
+real lambda  = E*sigma / ((1.+sigma)*(1.-2.*sigma));
+
+// Espace éléments finis
+fespace Vh(Th, [P1,P1,P1]);
+Vh [ux,uy,uz],[vx,vy,vz];
+
+// Macros (attention aux parenthèses superflues)
+macro Epsilon(u1,u2,u3) [dx(u1), dy(u2), dz(u3), 
+                         (dy(u3)+dz(u2))/sqrt(2.), 
+                         (dz(u1)+dx(u3))/sqrt(2.), 
+                         (dx(u2)+dy(u1))/sqrt(2.)] // EOM
+macro Div(u1,u2,u3) (dx(u1) + dy(u2) + dz(u3)) // EOM
+
+// Problème d'élasticité  on fixe les faces label=1
+problem Elasticite3D([ux,uy,uz],[vx,vy,vz], solver=sparsesolver) =
+    int3d(Th)(
+        lambda * Div(vx,vy,vz) * Div(ux,uy,uz)
+      + 2.*mu * (Epsilon(vx,vy,vz)' * Epsilon(ux,uy,uz))
+    )
+    - int3d(Th)(gravity * vz)
+    + on(1, ux=0, uy=0, uz=0);  
+
+// Résolution
+Elasticite3D;
+
+// Résultats
+cout << "Deplacement max |ux| = " << ux[].linfty << " m" << endl;
+cout << "Deplacement max |uy| = " << uy[].linfty << " m" << endl;
+cout << "Deplacement max |uz| = " << uz[].linfty << " m" << endl;
+
+// Création d'un espace scalaire P1
+fespace Wh(Th, P1);
+Wh normu;
+
+
+normu = sqrt(ux*ux + uy*uy + uz*uz);
+
+cout << "Deplacement max total = " << normu[].linfty << " m" << endl;
+
+
+
+cout << "Theorique (poutre)    = " << rho*9.81*(zmax-zmin)^2/(2.*E) << " m" << endl;
+
+// Visualisation
+medit("Norme du deplacement", Th, normu, wait=1);
+medit("Deplacement vertical uz", Th, uz, wait=1);
+
+// Maillage déformé 
+real umax = max(ux[].linfty, max(uy[].linfty, uz[].linfty));
+real coef = 0.1*(zmax-zmin)/umax;
+mesh3 Thd = movemesh3(Th, transfo=[x+coef*ux, y+coef*uy, z+coef*uz]);
+medit("Dodecaedre deforme", Thd, wait=1);
+
+
+
+
+
+
+// Exporter la solution en format VTK
+//int[int] order = [1, 1, 1]; 
+//savevtk("cube_deplacement.vtu", Th, [ux, uy, uz], dataname="Deplacement", order=order);
+//savevtk("cube_norme.vtu", Th, normu, dataname="Norme", order=order);
+//savevtk("cube_uz.vtu", Th, uz, dataname="Uz", order=order);
+
+// Exporter le maillage déformé
+//mesh3 Thd = movemesh3(Th, transfo=[x + 1e6*ux, y + 1e6*uy, z + 1e6*uz]);
+//savevtk("cube_deforme.vtu", Thd, [ux, uy, uz], dataname="Deplacement", order=order);
+
diff --git a/examples/Freefem/sofa/barre-traction.scn b/examples/Freefem/sofa/barre-traction.scn
new file mode 100644
index 0000000..1ba754f
--- /dev/null
+++ b/examples/Freefem/sofa/barre-traction.scn
@@ -0,0 +1,104 @@
+<?xml version="1.0"?>
+<!--
+    Bar under axial traction — SOFA equivalent to FreeFEM file
+
+    The bar is expanded from x=0 to L=1
+    The traction force F=1000N
+    Young's Modulus E = 1e6
+    Poisson's Ratio rho = 0
+
+    Exact solution: u(x) = (F / EA) * x  
+    This means that u(L) = 1.0
+
+    Important To Note:
+        The Elasticity plugin's Edge element integrates stiffness over element length 
+        only. It does not take into account any assumed cross section
+        The axial stiffness per unit length is (lambda + 2*mu). 
+        In the FreeFEM script, the cross-section A = 0.001.
+        Setting E_eff = E * A = 1e6 * 0.001 = 1000 recovers the same EA stiffness
+        as FreeFEM. Traction force is unchanged.
+
+    Expected Result:
+        u(L) = F * L / E_eff = 1000 * 1 / 1000 = 1.0
+        Displacement profile: linear  u_x(x) = x
+-->
+<Node name="root" gravity="0 0 0" dt="1">
+
+    <Node name="plugins">
+        <RequiredPlugin name="Elasticity"/>
+        <RequiredPlugin name="Sofa.Component.Constraint.Projective"/>  <!-- FixedProjectiveConstraint -->
+        <RequiredPlugin name="Sofa.Component.Engine.Select"/>           <!-- BoxROI -->
+        <RequiredPlugin name="Sofa.Component.LinearSolver.Direct"/>     <!-- SparseLDLSolver -->
+        <RequiredPlugin name="Sofa.Component.MechanicalLoad"/>          <!-- ConstantForceField -->
+        <RequiredPlugin name="Sofa.Component.ODESolver.Backward"/>      <!-- StaticSolver -->
+        <RequiredPlugin name="Sofa.Component.StateContainer"/>          <!-- MechanicalObject -->
+        <RequiredPlugin name="Sofa.Component.Topology.Container.Dynamic"/> <!-- EdgeSet* -->
+        <RequiredPlugin name="Sofa.Component.Topology.Container.Grid"/> <!-- RegularGridTopology -->
+        <RequiredPlugin name="Sofa.Component.Visual"/>                  <!-- VisualStyle, LineAxis -->
+        <RequiredPlugin name="Sofa.GL.Component.Rendering3D"/>          <!-- OglSceneFrame -->
+    </Node>
+
+    <DefaultAnimationLoop/>
+    <VisualStyle displayFlags="showBehaviorModels showForceFields showWireframe"/>
+
+    <Node name="bar">
+
+        <!-- This is where you define the  components that represent your -->
+        <!-- integration schemes for equillibrium problems-->
+        <!-- Even though we are solving a linear problem, the NewtonRaphsonSolver is employed-->
+        <!-- It will just do 1 iteration -->
+        <NewtonRaphsonSolver name="newtonSolver"
+            printLog="true"
+            warnWhenLineSearchFails="true"
+            maxNbIterationsNewton="1"
+            maxNbIterationsLineSearch="1"
+            lineSearchCoefficient="1"
+            relativeSuccessiveStoppingThreshold="0"
+            absoluteResidualStoppingThreshold="1e-7"
+            absoluteEstimateDifferenceThreshold="1e-12"
+            relativeInitialStoppingThreshold="1e-12"
+            relativeEstimateDifferenceThreshold="0"
+        />
+
+        <!-- Linear solver defined first-->
+        <!-- Equivalent to FreeFEM's direct linear solver-->
+        <SparseLDLSolver name="linearSolver" template="CompressedRowSparseMatrixd"/>
+
+        <!-- StaticSolver to find equilibrium-->
+        <StaticSolver name="staticSolver" newtonSolver="@newtonSolver" linearSolver="@linearSolver"/>
+
+        <!-- 8 equal elements: 9 nodes from x=0 to x=1              -->
+        <!-- Notice the 9 1 1, meaning 9 nodes in the x-axis,       -->
+        <!-- only 1 in each other direction of the grid i.e. a line -->
+        <RegularGridTopology name="grid" min="0 0 0" max="1 0 0" n="9 1 1"/>
+        <!-- The MechanicalObject is reponsible for storing information about mechanics-->
+        <MechanicalObject template="Vec3d" name="dofs" showObject="true" showObjectScale="0.02"/>
+
+        <Node name="edges">
+            <!-- The Topology must be defined explicitly; Important -->
+            <EdgeSetTopologyContainer name="topology"
+                position="@../dofs.position"
+                edges="@../grid.edges"/>
+
+            <!--
+                youngModulus = E * A = 1e6 * 0.001 = 1000  (effective axial stiffness)
+                poissonRatio = 0  =>  lambda=0, lambda+2*mu = E_eff  (no 3D coupling factor)
+            -->
+            <LinearSmallStrainFEMForceField name="FEM"
+                youngModulus="1000" poissonRatio="0.0"
+                topology="@topology"/>
+        </Node>
+
+        <!-- The BoxROI components help us select DOF indices in order to apply BCs-->
+        <BoxROI name="fixed_roi" box="-0.01 -1 -1  0.01 1 1" drawBoxes="false"/>
+        <!-- Dirichlet BC: u(x=0) = 0, that's the "Fixed" in ProjectiveConstraint-->
+        <FixedProjectiveConstraint indices="@fixed_roi.indices"/>
+
+        <!-- Neumann BC: point traction F=1000 at x=L=1 -->
+        <BoxROI name="load_roi" box="0.99 -1 -1  1.01 1 1" drawBoxes="false"/>
+        <!-- The "ForceField" in the name means that we apply a force on certain DOFs-->
+        <ConstantForceField indices="@load_roi.indices" forces="1000 0 0" showArrowSize="1e-4"/>
+
+    </Node>
+
+</Node>
diff --git a/examples/Freefem/sofa/beam-2d.scn b/examples/Freefem/sofa/beam-2d.scn
new file mode 100644
index 0000000..f08d903
--- /dev/null
+++ b/examples/Freefem/sofa/beam-2d.scn
@@ -0,0 +1,129 @@
+<?xml version="1.0"?>
+<!--
+    2D cantilever beam under gravity — SOFA equivalent of beam-2d.edp (FreeFEM)
+    Two formulations shown side by side:
+
+    In SOFA with the Vec1, Vec2, Vec3 templates. Depending on how these are set up, 
+    in order to generate a similar outcome as FreeFEM, we need to make some adjustments.
+    Just like in many FE codes, we can solve 2D problems in 3D space. This is what using a Vec3
+    template does. The 'd' in Vec3d is the 'double' type of C++.
+
+    Depending on the template we choose, SOFA will do things differently:
+
+    In the Node with name:
+    1. beam_plane_stress (Vec2)  y in [0, 2]
+        Lame lambda = E*nu / ((1+nu)*(1-nu))
+        DOFs: 2 per node (u_x, u_y). Displacement in the x and y directions.
+
+    2. beam_plane_strain (Vec3)  y in [3, 5]
+        Lame lambda = E*nu / ((1+nu)*(1-2*nu)) !!! Notice the 1-2*nu contrary to the above
+        DOFs: 3 per node (u_x, u_y, u_z)
+
+    The FreeFEM file (beam-2d.edp) uses 3D Lame coefficients, so formulation 2 is the exact replicate.
+    Formulation 1 shows the Vec2d approach. There will be a difference in the results between the
+    two.
+
+    Young's modulus: E=21.5 
+    Poisson's ratio: nu=0.29
+-->
+
+<!-- Gravity is applied at the root node -->
+<Node name="root" gravity="0 -0.05 0" dt="1">
+
+    <Node name="plugins">
+        <RequiredPlugin name="Elasticity"/>
+        <RequiredPlugin name="Sofa.Component.Constraint.Projective"/>
+        <RequiredPlugin name="Sofa.Component.Engine.Select"/>
+        <RequiredPlugin name="Sofa.Component.LinearSolver.Direct"/>
+        <RequiredPlugin name="Sofa.Component.Mass"/>
+        <RequiredPlugin name="Sofa.Component.ODESolver.Backward"/>
+        <RequiredPlugin name="Sofa.Component.StateContainer"/>
+        <RequiredPlugin name="Sofa.Component.Topology.Container.Dynamic"/>
+        <RequiredPlugin name="Sofa.Component.Topology.Container.Grid"/>
+        <RequiredPlugin name="Sofa.Component.Visual"/>
+        <RequiredPlugin name="Sofa.GL.Component.Rendering3D"/>
+    </Node>
+
+    <DefaultAnimationLoop/>
+    <VisualStyle displayFlags="showBehaviorModels showForceFields showWireframe"/>
+
+    <!-- All template='' cases here are Vec2d-->
+    <Node name="beam_plane_stress">
+
+        <NewtonRaphsonSolver name="newtonSolver"
+            printLog="true"
+            warnWhenLineSearchFails="true"
+            maxNbIterationsNewton="1"
+            maxNbIterationsLineSearch="1"
+            lineSearchCoefficient="1"
+            relativeSuccessiveStoppingThreshold="0"
+            absoluteResidualStoppingThreshold="1e-7"
+            absoluteEstimateDifferenceThreshold="1e-12"
+            relativeInitialStoppingThreshold="1e-12"
+            relativeEstimateDifferenceThreshold="0"
+        />
+        <SparseLDLSolver name="linearSolver" template="CompressedRowSparseMatrixd"/>
+        <StaticSolver name="staticSolver" newtonSolver="@newtonSolver" linearSolver="@linearSolver"/>
+
+        <RegularGridTopology name="grid" min="0 0 0" max="10 2 0" n="41 11 1"/>
+        <MechanicalObject template="Vec2d" name="dofs"/>
+
+        <Node name="triangles">
+            <TriangleSetTopologyContainer name="topology" src="@../grid"/>
+
+            <MeshMatrixMass template="Vec2d,Vec2d" totalMass="20" topology="@topology"/>
+            <LinearSmallStrainFEMForceField name="FEM" template="Vec2d"
+                youngModulus="21.5" poissonRatio="0.29"
+                topology="@topology"
+            />
+        </Node>
+
+        <BoxROI name="fixed_roi" template="Vec2d"
+            box="-0.01 -0.01 -0.01  0.01 2.01 0.01" drawBoxes="true"/>
+        <FixedProjectiveConstraint template="Vec2d" indices="@fixed_roi.indices"/>
+
+    </Node>
+
+    <!-- All template='' cases here are Vec3d-->
+    <Node name="beam_plane_strain">
+
+        <NewtonRaphsonSolver name="newtonSolver"
+            printLog="true"
+            warnWhenLineSearchFails="true"
+            maxNbIterationsNewton="1"
+            maxNbIterationsLineSearch="1"
+            lineSearchCoefficient="1"
+            relativeSuccessiveStoppingThreshold="0"
+            absoluteResidualStoppingThreshold="1e-7"
+            absoluteEstimateDifferenceThreshold="1e-12"
+            relativeInitialStoppingThreshold="1e-12"
+            relativeEstimateDifferenceThreshold="0"
+        />
+        <SparseLDLSolver name="linearSolver" template="CompressedRowSparseMatrixd"/>
+        <StaticSolver name="staticSolver" newtonSolver="@newtonSolver" linearSolver="@linearSolver"/>
+
+        <RegularGridTopology name="grid" min="0 3 0" max="10 5 0" n="41 11 1"/>
+        <MechanicalObject template="Vec3d" name="dofs"/>
+
+        <Node name="triangles">
+            <TriangleSetTopologyContainer name="topology" src="@../grid"/>
+            <TriangleSetTopologyModifier/>
+            <TriangleSetGeometryAlgorithms template="Vec3d" drawTriangles="true"/>
+
+            <MeshMatrixMass totalMass="20" topology="@topology"/>
+            <LinearSmallStrainFEMForceField name="FEM"
+                youngModulus="21.5" poissonRatio="0.29"
+                topology="@topology"
+                mstate="@../dofs"/>
+        </Node>
+
+        <!-- Clamp left edge: x=0, y in [3,5].
+             FixedProjectiveConstraint fixes all 3 DOFs (u_x, u_y, u_z)
+             of the left-edge nodes, which is all that is needed. -->
+        <BoxROI name="fixed_roi" template="Vec3d"
+            box="-0.01 2.99 -1.0  0.01 5.01 1.0" drawBoxes="true"/>
+        <FixedProjectiveConstraint template="Vec3d" indices="@fixed_roi.indices"/>
+
+    </Node>
+
+</Node>