Planar relationship constraint

[AliceCatalano]

Good evening,
I am trying to implement something from some months but I tried every component related to it and still is not working, so I hoped I could find some suggestions here.
I have a simple cube mesh made soft and deformable through the Mass-Spring ForceField.
The idea is to define some flat areas of the surface of the cube, and constrain them in a planar relationship, and then change the orientation of the constrained plane to obtain different "slopes".
I tired PlaneConstraint, PlaneROI, ProjectToPlaneConstraint, yet I was not able to obtain the desired behavior, because I want the area enclosed in the 4 mesh indexes to be rigid, while in all the above cases it remained with the same deformable characteristics as the remaining mesh.

This is the scene
import Sofa
import os
import numpy as np
os.path.expanduser('~/Desktop/SOFA_v22.12.00_Linux/share/sofa/examples/Tutorials/Sofa_scene')
from imu1 import MyController

def createScene(root):
root.gravity = [0, 0, 0]
root.dt = 0.1
root.animate=True
root.addObject('OglSceneFrame', style="Arrows", alignment="TopRight")
root.addObject('MeshGmshLoader', name="Loaderk", filename="mesh/cube.msh")

# Required Plugins
root.addObject('RequiredPlugin', pluginName=["Sofa.Component.Constraint.Projective","Sofa.Component.Engine.Select","Sofa.Component.LinearSolver.Iterative","Sofa.Component.Mass",
                                        "Sofa.Component.ODESolver.Backward","Sofa.Component.SolidMechanics.Spring","Sofa.Component.StateContainer","Sofa.Component.Topology.Container.Grid",
                                        "Sofa.Component.Visual","Sofa.Component.MechanicalLoad", "Sofa.Component.IO.Mesh", "Sofa.Component.Topology.Container.Dynamic", "Sofa.Component.Mapping.Linear", "Sofa.GL.Component.Rendering3D"])


root.addObject('DefaultAnimationLoop')
root.addObject('VisualStyle', displayFlags="showVisualModels showBehavior")

square = root.addChild('Square')
square.addObject('EulerImplicitSolver', printLog=0, rayleighStiffness = 0.5, rayleighMass = 0.5, vdamping=0)
square.addObject('CGLinearSolver', name="CG Solver", iterations=2, tolerance=1e-09, threshold=1e-09)
square.addObject('TetrahedronSetTopologyContainer', name="topo", src="@../../Loaderk")
square.addObject('TetrahedronSetGeometryAlgorithms', template="Vec3d", name="GeomAlgo")
square.addObject('MechanicalObject', src='@topo', template="Vec3d", name="MechanicalModel", showObject=1, dx = 0, dy =0, dz = 0, showObjectScale=1, showIndices=1, showIndicesScale=0.01)
square.addObject('MeshSpringForceField', name="Springs", stiffness=200, damping=1)
square.addObject('DiagonalMass', name="Mass", massDensity=26)
square.addObject('FixedConstraint', indices=[0, 1, 2, 3, 4, 5, 6, 7])

"""the following coordinates are the coordinates of the 3 highlighted points on the mesh"""
square.addObject('PlaneROI', name='ProjectToPlane', plane=[[-0.018125, 0.05249999, 0.025, 0.018125, 0.05249999, 0.025, -0.018125, 0.0175, 0.025, 0.025]], position='@MechanicalModel.rest_position', drawBoxes=True)
square.addObject('ProjectToPlaneConstraint', name='ptpcN', normal=[0, 0, 0], indices='@ProjectToPlane.indices', origin= [0, 0.035, 0.025], drawSize=0.001, listening=1)

visual = square.addChild('Visual')
visual.addObject('OglModel', name="VisualModel", src="@../../Loaderk", color="red")
visual.addObject('BarycentricMapping',template="Vec3d,Vec3d", name="visual mapping", input="@../MechanicalModel", output="@VisualModel")
visual.addObject('IdentityMapping', input="@../MechanicalModel", output="@VisualModel")

#root.addObject(KeyPressedController(name = "control"))
root.addObject(MyController(name="MyController", node=root, printLog=0))

return root