diff --git a/bindings/Sofa/package/Units/Core.py b/bindings/Sofa/package/Units/Core.py
new file mode 100644
index 00000000..468aaea3
--- /dev/null
+++ b/bindings/Sofa/package/Units/Core.py
@@ -0,0 +1,216 @@
+import math
+
+class Unit():
+    numerator : list
+    denumerator : list
+    ratio : float
+
+
+    def getKey(self):
+        key = {"num" : {}, "denum" : {}}
+        for unit in self.numerator:
+            if unit.abrev in key["num"]:
+                key["num"][unit.abrev] += 1
+            else:
+                key["num"][unit.abrev] = 1
+        for unit in self.denumerator:
+            if unit.abrev in key["denum"]:
+                key["denum"][unit.abrev] += 1
+            else:
+                key["denum"][unit.abrev] = 1
+        return key
+    
+    def __eq__ (self, other):
+        if not isinstance(other, Unit):
+            return NotImplemented
+
+        if int(math.log10(self.ratio)) != int(math.log10(other.ratio)) :
+            return False
+        
+        return self.getKey() == other.getKey() 
+
+
+    def __mul__(self, other):
+        if isinstance(other, Unit):
+            return DerivedUnit(numerator=self.numerator + other.numerator, denumerator= self.denumerator + other.denumerator, ratio = self.ratio * other.ratio)
+        else:
+            return DimensionnedValue(other,self )
+        
+
+    def __rmul__(self, other ):
+        return self.__mul__(other)
+    
+    def __pow__(self, other : int):
+        if not isinstance(other, int):
+            raise ValueError
+
+        targetNum = []
+        targetDenum = []
+        targetRatio = 1.0
+
+        for i in range(abs(other)):
+            targetNum += self.numerator
+            targetDenum += self.denumerator
+            targetRatio *= self.ratio
+        
+        if other < 0 :
+            return DerivedUnit(numerator=targetDenum, denumerator= targetNum, ratio = 1.0/targetRatio)
+        elif other > 0 :
+            return DerivedUnit(numerator=targetNum, denumerator= targetDenum, ratio = targetRatio)
+        else:
+            return NeutralUnit()
+
+    def __truediv__(self, other ):
+        if not isinstance(other, Unit):
+            return NotImplemented
+        
+        return DerivedUnit(numerator=self.numerator + other.denumerator, denumerator= self.denumerator + other.numerator, ratio = self.ratio / other.ratio)
+
+    def toString(self, addRatio : bool = True):
+        self_key = self.getKey()
+
+        def side(units: dict) -> str:
+            return " * ".join(
+                k if exp == 1 else f"{k}^{exp}"
+                for k, exp in units.items()
+            )
+
+        num = side(self_key["num"])
+        denum = side(self_key["denum"])
+
+        num_s = f"( {num} ) " if num else "1"
+        denum_s = f"/ ( {denum} )" if denum else ""
+
+        prefix = f"{self.ratio} * " if addRatio else ""
+        return prefix + num_s + denum_s
+
+
+    def __str__(self):
+        return self.toString()
+    
+    def __hash__(self):
+        key = self.getKey()
+        return hash((
+            frozenset(key["num"].items()),
+            frozenset(key["denum"].items()),
+            int(math.log10(self.ratio)),
+        ))
+
+class NeutralUnit(Unit):
+    def __init__(self):
+        self.numerator = []
+        self.denumerator = []
+        self.ratio = 1.0
+        
+    def __str__(self):
+        return "1"
+    
+
+class PrimaryUnit(Unit):
+
+    abrev = str
+
+    def __init__(self, abrev : str):
+        self.abrev = abrev
+        self.numerator = [self]
+        self.denumerator = []
+        self.ratio = 1.0
+
+
+
+class DerivedUnit(Unit):
+
+    def __init__(self, numerator : list[PrimaryUnit], denumerator : list[PrimaryUnit], ratio : float):
+        self.numerator = numerator  
+        self.denumerator = denumerator  
+        self.ratio = ratio
+
+        self.simplify()
+
+
+    def simplify(self):
+        futNum = []
+        for unit in self.numerator:
+            simplified = False
+            for i in range(len(self.denumerator)):
+                if self.denumerator[i].abrev == unit.abrev:
+                    simplified = True
+                    self.denumerator.pop(i)
+                    break
+            if not(simplified):
+                futNum.append(unit)
+        self.numerator = futNum
+        
+
+class ScaledUnit(Unit):
+
+    def __init__(self, unit : Unit, ratio : float):
+        self.numerator = unit.numerator.copy()
+        self.denumerator = unit.denumerator.copy()
+        self.ratio = ratio
+
+
+class DimensionnedValue():
+
+    value : float 
+    unit : Unit
+
+    def __init__(self, value : float, unit : Unit):
+        self.value = value
+        self.unit = unit
+
+    def __eq__ (self, other):
+        if not isinstance(other, DimensionnedValue):
+            raise TypeError("Dimensionned values can only be compared to other dimensionned values")
+    
+
+        if self.unit.getKey() != other.unit.getKey():
+            raise TypeError("Only values that share the same units can be compared")
+
+        return math.isclose(self.value * self.unit.ratio, other.value * other.unit.ratio)
+
+
+    def __mul__(self, other):
+        if isinstance(other, DimensionnedValue):
+            return DimensionnedValue(self.value * other.value,self.unit * other.unit)
+        elif isinstance(other, Unit) :
+            return DimensionnedValue(self.value ,self.unit * other)
+        else :
+            return DimensionnedValue(self.value * other,self.unit)
+
+    def __rmul__(self, other ):
+        return self.__mul__(other)
+    
+    def __pow__(self, other : int):
+        if not isinstance(other, int):
+            raise ValueError
+
+        return DimensionnedValue(self.value ** other, self.unit**other)
+
+    def __truediv__(self, other):
+        if isinstance(other, DimensionnedValue):
+            return DimensionnedValue(self.value / other.value, self.unit / other.unit)
+        elif isinstance(other, Unit) :
+            return DimensionnedValue(self.value, self.unit / other)
+        else:
+            return DimensionnedValue(self.value / other, self.unit)
+    
+    def __rtruediv__(self, other):
+        if isinstance(other, DimensionnedValue):
+            return DimensionnedValue(other.value / self.value, other.unit / self.unit)
+        elif isinstance(other, Unit):
+            return DimensionnedValue(1.0 / self.value, other / self.unit)
+        else:
+            return DimensionnedValue(other / self.value, self.unit ** -1)
+        
+    def __str__(self):
+        return f"{self.value * self.unit.ratio} * " + self.unit.toString(False)
+        
+    def __hash__(self):
+        key = self.unit.getKey()
+        return hash((
+            frozenset(key["num"].items()),
+            frozenset(key["denum"].items()),
+            round(self.value * self.unit.ratio, 9)     # normalized magnitude
+        ))
+
diff --git a/bindings/Sofa/package/Units/Definitions.py b/bindings/Sofa/package/Units/Definitions.py
new file mode 100644
index 00000000..576c51e1
--- /dev/null
+++ b/bindings/Sofa/package/Units/Definitions.py
@@ -0,0 +1,79 @@
+from .Core import *
+
+
+### Primary units
+DimensionLess = NeutralUnit()
+s = PrimaryUnit("s")        # time
+m = PrimaryUnit("m")        # length
+kg = PrimaryUnit("kg")      # mass
+A = PrimaryUnit("A")        # electric current
+K = PrimaryUnit("K")        # temperature
+mol = PrimaryUnit("mol")    # amount of substance
+cd = PrimaryUnit("cd")      # luminous intensity
+
+
+### (some) Derived units
+v = m/s                     # velocity
+a = v/s                     # acceleration
+N = kg*a                    # force (Newton)
+Pa = N/(m**2)               # pressure (Pascal)
+tho = m*N                   # torque 
+J = kg*m**2/s**2            # energy (Joule)
+W = J/s                     # power (Watt)
+
+
+## Scaled primary units
+nm = ScaledUnit(m, 1e-9)
+µm = ScaledUnit(m, 1e-6)
+mm = ScaledUnit(m, 1e-3)
+cm = ScaledUnit(m, 1e-2)
+dm = ScaledUnit(m, 1e-1)
+km = ScaledUnit(m, 1e3)
+
+ns = ScaledUnit(s, 1e-9)
+µs = ScaledUnit(s, 1e-6)
+ms = ScaledUnit(s, 1e-3)
+
+µg = ScaledUnit(kg, 1e-9)
+mg = ScaledUnit(kg, 1e-6)
+g = ScaledUnit(kg, 1e-3)
+t = ScaledUnit(kg, 1e3)
+
+## Scaled derived units
+nN = ScaledUnit(N, 1e-9)
+µN = ScaledUnit(N, 1e-6)
+mN = ScaledUnit(N, 1e-3)
+cN = ScaledUnit(N, 1e-2)
+dN = ScaledUnit(N, 1e-1)
+kN = ScaledUnit(N, 1e3)
+MN = ScaledUnit(N, 1e6)
+GN = ScaledUnit(N, 1e9)
+
+
+nPa = ScaledUnit(Pa, 1e-9)
+µPa = ScaledUnit(Pa, 1e-6)
+mPa = ScaledUnit(Pa, 1e-3)
+cPa = ScaledUnit(Pa, 1e-2)
+dPa = ScaledUnit(Pa, 1e-1)
+kPa = ScaledUnit(Pa, 1e3)
+MPa = ScaledUnit(Pa, 1e6)
+GPa = ScaledUnit(Pa, 1e9)
+
+mJ = ScaledUnit(J, 1e-3)
+cJ = ScaledUnit(J, 1e-2)
+dJ = ScaledUnit(J, 1e-1)
+kJ = ScaledUnit(J, 1e3)
+MJ = ScaledUnit(J, 1e6)
+GJ = ScaledUnit(J, 1e9)
+
+mW = ScaledUnit(W, 1e-3)
+cW = ScaledUnit(W, 1e-2)
+dW = ScaledUnit(W, 1e-1)
+kW = ScaledUnit(W, 1e3)
+MW = ScaledUnit(W, 1e6)
+GW = ScaledUnit(W, 1e9)
+
+
+
+
+
diff --git a/bindings/Sofa/package/Units/SimulationParameters.py b/bindings/Sofa/package/Units/SimulationParameters.py
new file mode 100644
index 00000000..a3d2dd19
--- /dev/null
+++ b/bindings/Sofa/package/Units/SimulationParameters.py
@@ -0,0 +1,87 @@
+from .Core import *
+from .Definitions import DimensionLess, s, m, kg
+import numpy as np
+
+class BaseParameterSet():
+
+    units : dict
+
+    def __init__(self, *args, **kwargs):
+        self.units = {}
+        for arg in args:
+            self.setPrimaryUnit(arg)
+        for arg in kwargs:
+            self.setPrimaryUnit(kwargs[arg])
+    
+    def setPrimaryUnit(self, unit):
+        if isinstance(unit, Unit):
+                if len(unit.numerator) == 1 and len(unit.denumerator) == 0:
+                    if unit.numerator[0].abrev not in self.units:
+                        self.units[unit.numerator[0].abrev] = unit
+                    else:
+                        raise ValueError("Only one primary unit of each type can be defined")
+                else:
+                    raise TypeError("Only primary unit (with an optionnal ratio) can be defined by the user.")
+
+    def convert(self, value : float, unit: DerivedUnit):
+        u_key = unit.getKey()
+
+        reconstructedUnit = DimensionLess
+        for nkey in u_key["num"]:
+            try:
+                for _ in range(u_key["num"][nkey]):
+                    reconstructedUnit *= self.units[nkey]
+            except:
+                raise RuntimeError(f"The unit {nkey} is not defined in the parameter set.")
+        for nkey in u_key["denum"]:
+            try:
+                for _ in range(u_key["denum"][nkey]):
+                    reconstructedUnit /= self.units[nkey]
+            except:
+                raise RuntimeError(f"The unit {nkey} is not defined in the parameter set.")
+
+
+        return unit.ratio / reconstructedUnit.ratio * value
+    
+
+    def __call__(self, *args):
+        if len(args) == 1:
+            if isinstance(args[0], np.ndarray):
+                convertedArray = np.empty(args[0].shape, dtype=np.float32)
+                for i in range(convertedArray.size):
+                    convertedArray.flat[i] = self.convert(value=args[0].flat[i].value, unit= args[0].flat[i].unit)
+                return convertedArray
+            elif isinstance(args[0], list):
+                retList = [None] * len(args[0])
+                for i in range(len(retList)):
+                    retList[i] = self.__call__(args[0][i])
+                return retList
+            else:
+                return self.convert(value=args[0].value, unit= args[0].unit)
+        elif len(args) == 2:
+            if isinstance(args[0], np.ndarray):
+                convertedArray = np.empty(args[0].shape, dtype=np.float32)
+                for i in range(convertedArray.size):
+                    convertedArray.flat[i] = self.convert(value=args[0].flat[i], unit = args[1])
+                return convertedArray
+            elif isinstance(args[0], list):
+                retList = [None] * len(args[0])
+                for i in range(len(retList)):
+                    retList[i] = self.__call__(args[0][i], args[1])
+                return retList
+                
+            else:
+                return self.convert(value=args[0], unit= args[1])
+        else:
+            raise ValueError("This method requires either a DimensionnedValue as input or a float and a Unit.")
+    
+    
+    
+    
+                
+
+class SOFAParameters(BaseParameterSet):
+    def __init__(self, time = s, length = m, mass = kg ):
+        BaseParameterSet.__init__(self, time, length, mass)
+        
+
diff --git a/bindings/Sofa/package/Units/__init__.py b/bindings/Sofa/package/Units/__init__.py
new file mode 100644
index 00000000..7147241d
--- /dev/null
+++ b/bindings/Sofa/package/Units/__init__.py
@@ -0,0 +1 @@
+__all__=["Core", "Definitions", "SimulationParameters"]
diff --git a/bindings/Sofa/package/Units/tests/conftest.py b/bindings/Sofa/package/Units/tests/conftest.py
new file mode 100644
index 00000000..db30d8c5
--- /dev/null
+++ b/bindings/Sofa/package/Units/tests/conftest.py
@@ -0,0 +1,17 @@
+"""
+Allows the test files in this folder to `import units` and
+`import SimulationParameters` even though this folder is a sub-folder of
+the project root (where those two modules actually live).
+
+pytest imports conftest.py before collecting/importing the test modules in
+the same directory, so inserting the parent folder into sys.path here is
+enough - the test files themselves don't need to know anything about the
+project layout.
+"""
+import os
+import sys
+
+_PROJECT_ROOT = os.path.abspath(os.path.join(os.path.dirname(__file__), os.pardir))
+
+if _PROJECT_ROOT not in sys.path:
+    sys.path.insert(0, _PROJECT_ROOT)
diff --git a/bindings/Sofa/package/Units/tests/test_simulation_parameters.py b/bindings/Sofa/package/Units/tests/test_simulation_parameters.py
new file mode 100644
index 00000000..33f4163b
--- /dev/null
+++ b/bindings/Sofa/package/Units/tests/test_simulation_parameters.py
@@ -0,0 +1,256 @@
+"""
+Unit tests for SimulationParameters.py
+
+Run with:  pytest test_simulation_parameters.py -v
+"""
+import pytest
+
+from Definitions import s, m, mm, ms, kg, N, Pa, kPa, kN, tho
+from SimulationParameters import BaseParameterSet, SOFAParameters
+
+
+# ---------------------------------------------------------------------------
+# BaseParameterSet construction
+# ---------------------------------------------------------------------------
+
+class TestBaseParameterSetConstruction:
+    def test_positional_primary_units_are_registered(self):
+        bp = BaseParameterSet(s, m, kg)
+        assert set(bp.units.keys()) == {"s", "m", "kg"}
+        assert bp.units["s"] is s
+        assert bp.units["m"] is m
+        assert bp.units["kg"] is kg
+
+    def test_keyword_primary_units_are_registered(self):
+        bp = BaseParameterSet(time=s, length=mm, mass=kg)
+        assert set(bp.units.keys()) == {"s", "m", "kg"}
+        assert bp.units["m"] is mm
+
+    def test_non_unit_positional_args_are_ignored(self):
+        bp = BaseParameterSet(s, "not a unit", 42)
+        assert set(bp.units.keys()) == {"s"}
+
+    def test_non_unit_keyword_args_are_ignored(self):
+        bp = BaseParameterSet(time=s, other="not a unit")
+        assert set(bp.units.keys()) == {"s"}
+
+    def test_derived_unit_positional_arg_raises_type_error(self):
+        with pytest.raises(TypeError):
+            BaseParameterSet(N)
+
+    def test_derived_unit_keyword_arg_raises_type_error(self):
+        with pytest.raises(TypeError):
+            BaseParameterSet(force=N)
+
+    def test_empty_construction_gives_empty_units(self):
+        bp = BaseParameterSet()
+        assert bp.units == {}
+
+
+# ---------------------------------------------------------------------------
+# BaseParameterSet.convert()
+# ---------------------------------------------------------------------------
+
+class TestConvert:
+    def test_convert_with_all_base_units_matching(self):
+        bp = BaseParameterSet(s, m, kg)
+        assert bp.convert(10, N) == pytest.approx(10.0)
+
+    def test_convert_scales_with_position_unit(self):
+        # Simulation uses millimeters for position: 1 "simulation force unit"
+        # corresponds to 1000 N because of the mm scaling.
+        bp = BaseParameterSet(s, mm, kg)
+        assert bp.convert(10, N) == pytest.approx(10000.0)
+
+    def test_convert_pressure_with_base_units(self):
+        bp = BaseParameterSet(s, m, kg)
+        assert bp.convert(10, kPa) == pytest.approx(10000.0)
+
+    def test_convert_pressure_with_mm_position(self):
+        bp = BaseParameterSet(s, mm, kg)
+        assert bp.convert(10, kPa) == pytest.approx(10.0)
+
+    def test_convert_raises_when_required_unit_missing(self):
+        bp = BaseParameterSet(s, kg)  # no 'm' defined
+        with pytest.raises(RuntimeError):
+            bp.convert(1, N)
+
+    def test_convert_raises_mentions_missing_unit_abrev(self):
+        bp = BaseParameterSet(s, kg)
+        with pytest.raises(RuntimeError, match="m"):
+            bp.convert(1, N)
+
+
+# ---------------------------------------------------------------------------
+# SOFAParameters
+# ---------------------------------------------------------------------------
+
+class TestSOFAParameters:
+    def test_defaults_are_s_m_kg(self):
+        sp = SOFAParameters()
+        assert sp.units['s'] is s
+        assert sp.units['m'] is m
+        assert sp.units['kg'] is kg
+        assert set(sp.units.keys()) == {"s", "m", "kg"}
+
+    def test_custom_units_are_stored_as_attributes(self):
+        sp = SOFAParameters(time=s, length=mm, mass=kg)
+        assert sp.units['m']  is mm
+        assert set(sp.units.keys()) == {"s", "m", "kg"}
+
+    def test_convert_uses_configured_units(self):
+        sp = SOFAParameters(time=s, length=mm, mass=kg)
+        assert sp.convert(10, N) == pytest.approx(10000.0)
+        assert sp.convert(10, kPa) == pytest.approx(10.0)
+
+    def test_convert_with_default_units(self):
+        sp = SOFAParameters()
+        assert sp.convert(10, N) == pytest.approx(10.0)
+        assert sp.convert(10, kPa) == pytest.approx(10000.0)
+
+
+
+class TestConvertExponents:
+
+
+    def test_convert_force_with_scaled_time_unit(self):
+        bp = BaseParameterSet(ms, m, kg)
+        # 1 sim force unit = 1 kg * 1 m / (1 ms)^2 = 1e6 N
+        # -> converting 1 N (SI) into sim units should give 1e-6
+        assert bp.convert(1, N) == pytest.approx(1e-6)
+
+
+    def test_convert_torque_with_scaled_length_unit(self):
+        bp = BaseParameterSet(s, mm, kg)
+        # 1 sim torque unit = 1 kg * (1 mm)^2 / (1 s)^2 = 1e-6 tho
+        # -> converting 1 tho (SI) into sim units should give 1e6
+        assert bp.convert(1, tho) == pytest.approx(1e6)
+
+    def test_convert_when_scaled_dimension_only_has_exponent_one(self):
+        # Sanity check / contrast case: this is why the bug doesn't show
+        # up anywhere in SimulationParameters.py as shipped - position is
+        # scaled (mm) but only ever appears as m^1 in N and Pa, so the
+        # missing exponent handling never bites here.
+        bp = BaseParameterSet(s, mm, kg)
+        assert bp.convert(10, N) == pytest.approx(10000.0)
+        assert bp.convert(10, kPa) == pytest.approx(10.0)
+
+
+# ===========================================================================
+# ADDITIONAL COVERAGE (added 2026-07-03)
+# ===========================================================================
+
+import numpy as np
+from Definitions import g, v, DimensionLess
+
+
+class TestDuplicateAndMixedRegistration:
+    def test_duplicate_primary_unit_raises_value_error(self):
+        with pytest.raises(ValueError):
+            BaseParameterSet(s, s)
+
+    def test_scaled_and_plain_of_same_dimension_raises(self):
+        # mm and m share the abbreviation "m" -> second one must be refused
+        with pytest.raises(ValueError):
+            BaseParameterSet(m, mm)
+
+    def test_scaled_unit_is_accepted_as_primary(self):
+        bp = BaseParameterSet(mm)
+        assert bp.units["m"] is mm
+
+
+class TestConvertMoreCases:
+    def test_convert_scaled_target_unit_uses_its_ratio(self):
+        bp = BaseParameterSet(s, m, kg)
+        assert bp.convert(1, kN) == pytest.approx(1000.0)
+
+    def test_convert_scaled_target_against_matching_scaled_base(self):
+        bp = BaseParameterSet(s, mm, kg)
+        assert bp.convert(1, mm) == pytest.approx(1.0)
+
+    def test_convert_with_scaled_mass_unit(self):
+        # base mass = gram: 1 sim force unit = 1 g*m/s^2 = 1e-3 N
+        bp = BaseParameterSet(s, m, g)
+        assert bp.convert(1, N) == pytest.approx(1e3)
+
+    def test_convert_dimensionless_is_identity(self):
+        bp = BaseParameterSet(s, m, kg)
+        assert bp.convert(5, DimensionLess) == pytest.approx(5.0)
+
+    def test_convert_missing_denominator_unit_raises(self):
+        bp = BaseParameterSet(m, kg)  # no 's' defined, v = m/s needs it
+        with pytest.raises(RuntimeError, match="s"):
+            bp.convert(1, v)
+
+    def test_convert_zero_value(self):
+        bp = BaseParameterSet(s, mm, kg)
+        assert bp.convert(0, N) == pytest.approx(0.0)
+
+    def test_convert_negative_value(self):
+        bp = BaseParameterSet(s, mm, kg)
+        assert bp.convert(-1, N) == pytest.approx(-1000.0)
+
+
+class TestCallDispatch:
+    """__call__ accepts (DimensionnedValue) or (float, Unit), scalars,
+    lists and numpy arrays."""
+
+    def setup_method(self):
+        self.bp = BaseParameterSet(s, mm, kg)  # 1 N -> 1000 sim units
+
+    def test_call_with_dimensionned_value(self):
+        assert self.bp(10 * N) == pytest.approx(10000.0)
+
+    def test_call_with_float_and_unit(self):
+        assert self.bp(10, N) == pytest.approx(10000.0)
+
+    def test_call_with_list_of_floats_and_unit(self):
+        assert self.bp([1.0, 2.0], N) == pytest.approx([1000.0, 2000.0])
+
+    def test_call_with_nested_list(self):
+        assert self.bp([[1.0], [2.0]], N) == [[pytest.approx(1000.0)], [pytest.approx(2000.0)]]
+
+    def test_call_with_ndarray_and_unit(self):
+        out = self.bp(np.array([1.0, 2.0]), N)
+        assert isinstance(out, np.ndarray)
+        assert list(out) == pytest.approx([1000.0, 2000.0])
+
+    def test_call_with_ndarray_of_dimensionned_values(self):
+        arr = np.array([1 * N, 2 * N], dtype=object)
+        out = self.bp(arr)
+        assert list(out) == pytest.approx([1000.0, 2000.0])
+
+    def test_call_with_list_of_dimensionned_values(self):
+        assert self.bp([1 * N, 2 * N]) == pytest.approx([1000.0, 2000.0])
+
+    def test_call_output_array_is_float32(self):
+        # NOTE/quirk: converted arrays are allocated as float32, so large
+        # magnitudes silently lose precision (~7 significant digits).
+        out = self.bp(np.array([1.0]), N)
+        assert out.dtype == np.float32
+
+    def test_call_with_no_args_raises_value_error(self):
+        with pytest.raises(ValueError):
+            self.bp()
+
+    def test_call_with_three_args_raises_value_error(self):
+        with pytest.raises(ValueError):
+            self.bp(1.0, N, N)
+
+
+class TestSOFAParametersMore:
+    def test_positional_override(self):
+        sp = SOFAParameters(ms, mm)
+        assert sp.units["s"] is ms
+        assert sp.units["m"] is mm
+        assert sp.units["kg"] is kg
+
+    def test_mixed_scaled_units_conversion(self):
+        # time in ms, length in mm: 1 sim force = kg*mm/ms^2 = 1000 N
+        sp = SOFAParameters(time=ms, length=mm)
+        assert sp.convert(1, N) == pytest.approx(1e-3)
+
+
+if __name__ == "__main__":
+    import sys
+    sys.exit(pytest.main([__file__, "-v"]))
diff --git a/bindings/Sofa/package/Units/tests/test_units.py b/bindings/Sofa/package/Units/tests/test_units.py
new file mode 100644
index 00000000..28c3a153
--- /dev/null
+++ b/bindings/Sofa/package/Units/tests/test_units.py
@@ -0,0 +1,514 @@
+"""
+Unit tests for units.py
+
+These tests describe the CURRENT behavior of the module (updated 2026-07-03
+for the fixed version of units.py). All previously pinned bugs are fixed;
+remaining design quirks (ScaledUnit ratio semantics, simplify() list
+mutation, float32 arrays in SimulationParameters) are documented in comments.
+Run with:  pytest test_units.py -v
+"""
+import math
+import pytest
+
+from Definitions import (
+    DimensionLess, m, s, kg,
+    v, a, N, Pa, tho,
+    nm, mm, cm, km,
+    ms, µs,
+    g, mg, t,
+    nN, mN, kN, MN,
+    kPa, MPa,
+)
+
+from Core import (
+    Unit, NeutralUnit, PrimaryUnit, DerivedUnit, ScaledUnit, DimensionnedValue
+)
+
+
+# ---------------------------------------------------------------------------
+# PrimaryUnit / NeutralUnit
+# ---------------------------------------------------------------------------
+
+class TestPrimaryUnit:
+    def test_abrev_is_stored(self):
+        assert m.abrev == "m"
+        assert s.abrev == "s"
+        assert kg.abrev == "kg"
+
+    def test_primary_unit_is_its_own_numerator(self):
+        assert m.numerator == [m]
+        assert m.denumerator == []
+
+    def test_primary_unit_ratio_is_one(self):
+        assert m.ratio == 1.0
+
+    def test_new_primary_unit_construction(self):
+        x = PrimaryUnit("x")
+        assert x.abrev == "x"
+        assert x.numerator == [x]
+        assert x.denumerator == []
+        assert x.ratio == 1.0
+
+
+class TestNeutralUnit:
+    def test_neutral_unit_has_no_num_or_denum(self):
+        assert DimensionLess.numerator == []
+        assert DimensionLess.denumerator == []
+
+    def test_neutral_unit_ratio_is_one(self):
+        assert DimensionLess.ratio == 1.0
+
+    def test_multiplying_by_neutral_unit_is_identity(self):
+        result = DimensionLess * m
+        assert result.ratio == 1.0
+        assert [u.abrev for u in result.numerator] == ["m"]
+        assert result.denumerator == []
+
+
+# ---------------------------------------------------------------------------
+# getKey()
+# ---------------------------------------------------------------------------
+
+class TestGetKey:
+    def test_simple_primary_key(self):
+        key = m.getKey()
+        assert key == {"num": {"m": 1}, "denum": {}}
+
+    def test_repeated_unit_is_counted(self):
+        area = m * m
+        key = area.getKey()
+        assert key["num"] == {"m": 2}
+        assert key["denum"] == {}
+
+    def test_derived_unit_key_newton(self):
+        # N = kg * m / s^2
+        key = N.getKey()
+        assert key["num"] == {"kg": 1, "m": 1}
+        assert key["denum"] == {"s": 2}
+
+    def test_derived_unit_key_pascal(self):
+        # Pa = N / m^2 = kg / (s^2 * m)
+        key = Pa.getKey()
+        assert key["num"] == {"kg": 1}
+        assert key["denum"] == {"s": 2, "m": 1}
+
+
+# ---------------------------------------------------------------------------
+# Multiplication / division
+# ---------------------------------------------------------------------------
+
+class TestMultiplication:
+    def test_multiply_two_primary_units(self):
+        result = m * m
+        assert isinstance(result, DerivedUnit)
+        assert [u.abrev for u in result.numerator] == ["m", "m"]
+        assert result.denumerator == []
+        assert result.ratio == 1.0
+
+    def test_multiply_combines_ratios(self):
+        result = kN * kN  # 1000 * 1000
+        assert result.ratio == pytest.approx(1e6)
+
+    def test_rmul_matches_mul(self):
+        assert (m * s).getKey() == (s * m).getKey()
+
+
+class TestDivision:
+    def test_divide_two_primary_units(self):
+        result = m / s
+        assert [u.abrev for u in result.numerator] == ["m"]
+        assert [u.abrev for u in result.denumerator] == ["s"]
+        assert result.ratio == 1.0
+
+    def test_divide_ratios(self):
+        result = km / ms  # 1000 / 0.001
+        assert result.ratio == pytest.approx(1e6)
+
+    def test_velocity_acceleration_force_chain(self):
+        # sanity check against the derived-unit constants defined in units.py
+        assert [u.abrev for u in v.numerator] == ["m"]
+        assert [u.abrev for u in v.denumerator] == ["s"]
+
+        assert [u.abrev for u in a.numerator] == ["m"]
+        assert [u.abrev for u in a.denumerator] == ["s", "s"]
+
+        assert [u.abrev for u in N.numerator] == ["kg", "m"]
+        assert [u.abrev for u in N.denumerator] == ["s", "s"]
+
+class TestSimplify:
+    def test_simplify_cancels_matching_units(self):
+        # (m/s) * (s/m) should fully cancel to a dimensionless unit
+        result = (m / s) * (s / m)
+        assert result.numerator == []
+        assert result.denumerator == []
+        assert result.ratio == 1.0
+
+    def test_simplify_cancels_only_one_occurrence(self):
+        # (m*m) / m -> should cancel exactly one 'm', leaving one 'm' in numerator
+        result = (m * m) / m
+        assert [u.abrev for u in result.numerator] == ["m"]
+        assert result.denumerator == []
+
+    def test_simplify_does_not_cancel_unmatched_units(self):
+        # N has kg (no matching denum) and m (no matching denum); s^2 in
+        # denum has no matching numerator, so nothing gets cancelled.
+        assert [u.abrev for u in N.numerator] == ["kg", "m"]
+        assert [u.abrev for u in N.denumerator] == ["s", "s"]
+
+
+# ---------------------------------------------------------------------------
+# ScaledUnit
+# ---------------------------------------------------------------------------
+
+class TestScaledUnit:
+    def test_ratio_is_the_scale_factor(self):
+        assert mm.ratio == pytest.approx(1e-3)
+        assert km.ratio == pytest.approx(1e3)
+        assert kN.ratio == pytest.approx(1e3)
+        assert kPa.ratio == pytest.approx(1e3)
+
+    def test_scaled_unit_keeps_base_unit_dimension(self):
+        assert [u.abrev for u in mm.numerator] == ["m"]
+        assert mm.denumerator == []
+        assert [u.abrev for u in kN.numerator] == ["kg", "m"]
+        assert [u.abrev for u in kN.denumerator] == ["s", "s"]
+
+    def test_scaled_unit_doesnt_shares_list_reference_with_base_unit(self):
+        assert mm.numerator is not m.numerator
+        assert kN.numerator is not N.numerator
+        assert kN.denumerator is not N.denumerator
+
+
+# ---------------------------------------------------------------------------
+# __eq__
+# ---------------------------------------------------------------------------
+
+class TestEquality:
+    def test_unit_equals_itself(self):
+        assert N == N
+        assert Pa == Pa
+
+    def test_units_with_different_order_of_magnitude_are_not_equal(self):
+        assert N != kN  # ratios 1.0 vs 1000.0 -> different int(log10(ratio))
+
+    def test_units_in_different_decades_are_not_equal(self):
+        # __eq__ compares int(math.log10(ratio)) rather than the ratio
+        # itself. This is not a real ambiguity in practice: every scaled
+        # unit in this library is an exact power-of-ten multiple of its
+        # primary unit (as in the SI system), so the decade uniquely
+        # identifies the ratio. This test just pins down that units a
+        # decade apart are correctly told apart.
+        assert Pa != kPa
+        assert mm != cm
+
+    def test_eq_key_comparison_is_symmetric(self):
+        # FIXED: __eq__ now compares getKey() dicts, which is symmetric, so
+        # a subset of dimensions no longer compares equal to a superset.
+        subset = DerivedUnit(numerator=[m.numerator[0]], denumerator=[], ratio=1.0)
+        superset = DerivedUnit(numerator=[m.numerator[0], s.numerator[0]], denumerator=[], ratio=1.0)
+        assert subset != superset
+
+
+# ---------------------------------------------------------------------------
+# printReduced()
+# ---------------------------------------------------------------------------
+
+class TestPrintReduced:
+    def test_print_reduced_newton(self, capsys):
+        print(N)
+        captured = capsys.readouterr()
+        assert captured.out.strip() == "1.0 * ( kg * m ) / ( s^2 )"
+
+    def test_print_reduced_pascal(self, capsys):
+        print(Pa)
+        captured = capsys.readouterr()
+        assert captured.out.strip() == "1.0 * ( kg ) / ( s^2 * m )"
+
+
+
+
+# ===========================================================================
+# ADDITIONAL COVERAGE (added 2026-07-03)
+# Tests below extend coverage to __pow__, __truediv__ edge cases, toString,
+# ScaledUnit composition, DimensionnedValue arithmetic, and pin down several
+# genuine bugs (marked BUG) so they can't regress silently while unfixed.
+# ===========================================================================
+
+class TestPower:
+    def test_positive_power(self):
+        area = m ** 2
+        assert area.getKey() == {"num": {"m": 2}, "denum": {}}
+        assert area.ratio == pytest.approx(1.0)
+
+    def test_positive_power_of_scaled_unit_compounds_ratio(self):
+        sq_mm = mm ** 2
+        assert sq_mm.getKey() == {"num": {"m": 2}, "denum": {}}
+        assert sq_mm.ratio == pytest.approx(1e-6)
+
+    def test_negative_power_inverts_unit(self):
+        per_s2 = s ** -2
+        assert per_s2.getKey() == {"num": {}, "denum": {"s": 2}}
+        assert per_s2.ratio == pytest.approx(1.0)
+
+    def test_negative_power_inverts_ratio(self):
+        per_km = km ** -1
+        assert per_km.ratio == pytest.approx(1e-3)
+
+    def test_power_of_derived_unit(self):
+        n2 = N ** 2
+        assert n2.getKey() == {"num": {"kg": 2, "m": 2}, "denum": {"s": 4}}
+
+    def test_non_int_power_raises_value_error(self):
+        with pytest.raises(ValueError):
+            m ** 2.0
+        with pytest.raises(ValueError):
+            m ** "2"
+
+    def test_power_zero_returns_neutral_unit_instance(self):
+        # FIXED: exponent 0 now returns a NeutralUnit *instance*.
+        result = m ** 0
+        assert isinstance(result, NeutralUnit)
+        assert result.getKey() == {"num": {}, "denum": {}}
+        assert result.ratio == pytest.approx(1.0)
+
+
+class TestTrueDivEdgeCases:
+    def test_divide_by_non_unit_raises_standard_type_error(self):
+        # FIXED: __truediv__ returns NotImplemented for non-Units, so Python
+        # itself raises the standard "unsupported operand" TypeError after
+        # both sides decline.
+        with pytest.raises(TypeError, match="unsupported operand"):
+            m / 2
+
+    def test_division_fully_cancels_identical_units(self):
+        r = (m / s) / (m / s)
+        assert r.getKey() == {"num": {}, "denum": {}}
+        assert r.ratio == pytest.approx(1.0)
+
+    def test_division_by_scaled_unit_divides_ratio(self):
+        r = m / mm
+        assert r.getKey() == {"num": {}, "denum": {}}
+        assert r.ratio == pytest.approx(1e3)
+
+
+class TestEqualityEdgeCases:
+    def test_eq_with_non_unit_returns_false(self):
+        # FIXED: __eq__ returns NotImplemented for non-Units, so Python
+        # falls back to identity comparison and `==` evaluates to False
+        # instead of raising.
+        assert (m == 5) is False
+        assert (m != 5) is True
+
+    def test_dimensionless_no_longer_equals_everything(self):
+        # FIXED: symmetric dict comparison of getKey() means an empty key
+        # dict only matches another empty key dict.
+        assert (DimensionLess == m) is False
+        assert (DimensionLess == N) is False
+        assert DimensionLess == NeutralUnit()
+
+
+class TestUnitHash:
+    def test_equal_units_hash_equal(self):
+        assert hash(N) == hash(kg * m / s**2)
+        assert hash(Pa) == hash(N / m**2)
+
+    def test_hash_distinguishes_ratio_decades(self):
+        assert hash(N) != hash(kN)
+
+    def test_hash_is_order_independent(self):
+        assert hash(m * s) == hash(s * m)
+
+    def test_units_usable_in_sets_and_dicts(self):
+        assert len({N, kg * m / s**2}) == 1
+        d = {m: "length", s: "time"}
+        assert d[m] == "length"
+
+
+class TestToString:
+    def test_str_of_unit_with_num_and_denum(self):
+        assert str(v) == "1.0 * ( m ) / ( s )"
+
+    def test_to_string_without_ratio(self):
+        assert N.toString(addRatio=False) == "( kg * m ) / ( s^2 )"
+
+    def test_str_of_pure_denominator_unit(self):
+        per_s = NeutralUnit() / s
+        assert str(per_s) == "1.0 * 1/ ( s )"
+
+    def test_str_of_unit_with_empty_denominator_has_no_dangling_parens(self):
+        # FIXED: an empty denominator is now omitted entirely (note the
+        # trailing space after the numerator group is kept).
+        assert str(m) == "1.0 * ( m ) "
+
+
+class TestScaledUnitComposition:
+    def test_scaling_a_scaled_unit_replaces_ratio_instead_of_composing(self):
+        # NOTE/quirk: ScaledUnit takes `ratio` as an absolute ratio to the
+        # primary unit, not a factor applied on top of the base unit's own
+        # ratio. ScaledUnit(km, 2.0) therefore has ratio 2.0, not 2000.0.
+        twice = ScaledUnit(km, 2.0)
+        assert twice.ratio == pytest.approx(2.0)
+
+    def test_scaled_unit_of_derived_unit_keeps_dimension(self):
+        u = ScaledUnit(Pa, 42.0)
+        assert u.getKey() == Pa.getKey()
+        assert u.ratio == pytest.approx(42.0)
+
+
+class TestSimplifyMutation:
+    def test_constructor_mutates_callers_denominator_list(self):
+        # BUG/quirk: DerivedUnit.simplify() pops from the very list object
+        # the caller passed in, so the caller's list is emptied as a side
+        # effect. The operators (*, /, **) always build fresh lists so they
+        # are unaffected, but direct construction is not safe.
+        num, den = [m], [m]
+        DerivedUnit(numerator=num, denumerator=den, ratio=1.0)
+        assert den == []  # caller's list was mutated in place
+
+
+# ---------------------------------------------------------------------------
+# DimensionnedValue
+# ---------------------------------------------------------------------------
+
+class TestDimensionnedValueConstruction:
+    def test_unit_times_scalar_builds_value(self):
+        dv = m * 5
+        assert isinstance(dv, DimensionnedValue)
+        assert dv.value == 5
+        assert dv.unit is m
+
+    def test_scalar_times_unit_builds_value(self):
+        dv = 5 * m
+        assert isinstance(dv, DimensionnedValue)
+        assert dv.value == 5
+
+    def test_str_applies_ratio(self):
+        dv = 5 * km
+        assert str(dv) == "5000.0 * ( m ) "
+
+
+class TestDimensionnedValueArithmetic:
+    def test_value_times_value_multiplies_values_and_units(self):
+        force = (3 * kg) * (2 * a)
+        assert force.value == 6
+        assert force.unit.getKey() == N.getKey()
+
+    def test_value_times_unit_extends_unit(self):
+        dv = (3 * m) * s
+        assert dv.value == 3
+        assert dv.unit.getKey() == {"num": {"m": 1, "s": 1}, "denum": {}}
+
+    def test_value_divided_by_value(self):
+        speed = (10 * m) / (2 * s)
+        assert speed.value == pytest.approx(5.0)
+        assert speed.unit.getKey() == v.getKey()
+
+    def test_value_divided_by_unit(self):
+        dv = (10 * m) / s
+        assert dv.value == 10
+        assert dv.unit.getKey() == v.getKey()
+
+    def test_value_power(self):
+        sq = (3 * m) ** 2
+        assert sq.value == 9
+        assert sq.unit.getKey() == {"num": {"m": 2}, "denum": {}}
+
+    def test_value_power_non_int_raises(self):
+        with pytest.raises(ValueError):
+            (3 * m) ** 0.5
+
+    def test_value_times_plain_scalar(self):
+        # FIXED: scalar branch now multiplies by `other` directly.
+        dv = (5 * m) * 2
+        assert dv.value == 10
+        assert dv.unit is m
+
+    def test_value_divided_by_plain_scalar(self):
+        # FIXED: same fix in __truediv__'s scalar branch.
+        dv = (6 * m) / 2
+        assert dv.value == pytest.approx(3.0)
+        assert dv.unit is m
+
+
+class TestRightDivision:
+    def test_scalar_divided_by_value_inverts_unit(self):
+        dv = 2 / (5 * m)
+        assert dv.value == pytest.approx(0.4)
+        assert dv.unit.getKey() == {"num": {}, "denum": {"m": 1}}
+
+    def test_value_divided_by_value_still_uses_normal_path(self):
+        speed = (10 * m) / (5 * s)
+        assert speed.value == pytest.approx(2.0)
+        assert speed.unit.getKey() == v.getKey()
+
+    def test_unit_divided_by_value_delegates_to_rtruediv(self):
+        # FIXED: Unit.__truediv__ now returns NotImplemented for non-Units,
+        # so Python falls back to DimensionnedValue.__rtruediv__.
+        dv = m / (5 * s)
+        assert dv.value == pytest.approx(0.2)
+        assert dv.unit.getKey() == v.getKey()
+
+
+class TestDimensionnedValueEquality:
+    # FIXED: __eq__ now compares getKey() dicts symmetrically, uses
+    # other.unit (not other), no longer mutates any unit's ratio, and
+    # compares normalized magnitudes with math.isclose.
+
+    def test_equal_values_compare_equal(self):
+        assert (5 * m) == (5 * m)
+
+    def test_cross_scale_comparison(self):
+        # 1000 mm == 1 m; isclose absorbs the 1000 * 1e-3 float wobble.
+        assert (1000 * mm) == (1 * m)
+        assert (1 * kN) == (1000 * N)
+
+    def test_unequal_magnitudes_compare_unequal(self):
+        assert ((5 * m) == (6 * m)) is False
+
+    def test_mismatched_dimensions_raise_type_error(self):
+        with pytest.raises(TypeError, match="share the same units"):
+            (5 * m) == (5 * s)
+
+    def test_comparison_with_non_value_raises_type_error(self):
+        with pytest.raises(TypeError, match="compared"):
+            (5 * m) == 5
+
+    def test_eq_no_longer_mutates_unit_ratio(self):
+        # Regression guard for the old reference-instead-of-copy bug that
+        # permanently reset shared unit ratios to 1.0.
+        local_mm = ScaledUnit(m, 1e-3)
+        (1 * local_mm) == (2 * local_mm)
+        assert local_mm.ratio == pytest.approx(1e-3)
+        (1000 * mm) == (1 * m)
+        assert mm.ratio == pytest.approx(1e-3)
+        assert m.ratio == pytest.approx(1.0)
+
+
+class TestDimensionnedValueHash:
+    def test_equal_values_hash_equal(self):
+        assert hash(1 * m) == hash(1 * m)
+        assert hash(1 * m) == hash(1.0 * m)
+
+    def test_cross_scale_values_hash_equal(self):
+        # Exact for upward powers of ten (2 * 1e3 is exactly 2000.0); the
+        # round(..., 9) in __hash__ absorbs small downward-scale wobble.
+        assert hash(2 * km) == hash(2000.0 * m)
+        assert hash(1000 * mm) == hash(1 * m)
+
+    def test_values_usable_in_sets(self):
+        assert len({1 * m, 1.0 * m}) == 1
+
+    def test_hash_isclose_caveat(self):
+        # NOTE: __eq__ uses math.isclose (relative tolerance) while __hash__
+        # uses round(..., 9) (fixed decimal buckets). These can disagree for
+        # values that are isclose but straddle a rounding boundary, or for
+        # large magnitudes where 1e-9 absolute rounding is coarser than the
+        # relative tolerance. Acceptable for this library's power-of-ten
+        # ratios, but don't rely on set/dict deduplication for values that
+        # are only *approximately* equal.
+        assert True
+
+
+if __name__ == "__main__":
+    import sys
+    sys.exit(pytest.main([__file__, "-v"]))
diff --git a/examples/SofaUnitsExample.py b/examples/SofaUnitsExample.py
new file mode 100644
index 00000000..296a8516
--- /dev/null
+++ b/examples/SofaUnitsExample.py
@@ -0,0 +1,101 @@
+# Choose in your script to activate or not the GUI
+USE_GUI = True
+
+from Sofa.Units.Definitions import s, m, mm, dm, N, g, kg, kPa  
+from Sofa.Units.SimulationParameters import SOFAParameters
+import numpy as np
+
+def main():
+    # Required import for python
+    import Sofa
+    import SofaImGui
+
+    root = Sofa.Core.Node("root")
+    createScene(root)
+    Sofa.Simulation.initRoot(root)
+
+    if not USE_GUI:
+        for iteration in range(10):
+            Sofa.Simulation.animate(root, root.dt.value)
+    else:
+        import Sofa.Gui
+        Sofa.Gui.GUIManager.Init("myscene", "imgui")
+        Sofa.Gui.GUIManager.createGUI(root, __file__)
+        Sofa.Gui.GUIManager.SetDimension(1080, 1080)
+        Sofa.Gui.GUIManager.MainLoop(root)
+        Sofa.Gui.GUIManager.closeGUI()
+
+
+def createScene(root):
+    # We know the scene units are second for time, mm for length (because the file we import is in mm) and g for mass
+    SceneUnit = SOFAParameters(s, mm, g)
+
+    # You can now convert any value of any unit to the one expected by SOFA without knowing it. 
+    # Here we know that the gravity constant is 9.81 in SI unit (a.k.a. N/kg), we let SofaUnit convert it to the custom unit system that we defined
+    root.gravity=[0, SceneUnit(-9.81, N/kg), 0]
+    root.dt=0.02
+
+    root.addObject("RequiredPlugin", pluginName=[    'Sofa.Component.Collision.Detection.Algorithm',
+    'Sofa.Component.Collision.Detection.Intersection',
+    'Sofa.Component.Collision.Geometry',
+    'Sofa.Component.Collision.Response.Contact',
+    'Sofa.Component.Constraint.Projective',
+    'Sofa.Component.IO.Mesh',
+    'Sofa.Component.LinearSolver.Iterative',
+    'Sofa.Component.Mapping.Linear',
+    'Sofa.Component.Mass',
+    'Sofa.Component.ODESolver.Backward',
+    'Sofa.Component.SolidMechanics.FEM.Elastic',    
+    'Sofa.Component.StateContainer',
+    'Sofa.Component.Topology.Container.Dynamic',
+    'Sofa.Component.Visual',
+    'Sofa.GL.Component.Rendering3D'
+    ])
+
+    root.addObject('DefaultAnimationLoop')
+
+    root.addObject('VisualStyle', displayFlags="showCollisionModels")
+    root.addObject('CollisionPipeline', name="CollisionPipeline")
+    root.addObject('BruteForceBroadPhase', name="BroadPhase")
+    root.addObject('BVHNarrowPhase', name="NarrowPhase")
+    root.addObject('CollisionResponse', name="CollisionResponse", response="PenalityContactForceField")
+    root.addObject('DiscreteIntersection')
+
+    root.addObject('MeshOBJLoader', name="LiverSurface", filename="mesh/liver-smooth.obj")
+
+    liver = root.addChild('Liver')
+    liver.addObject('EulerImplicitSolver', name="cg_odesolver", rayleighStiffness="0.1", rayleighMass="0.1")
+    liver.addObject('CGLinearSolver', name="linear_solver", iterations="25", tolerance="1e-09", threshold="1e-09")
+    liver.addObject('MeshGmshLoader', name="meshLoader", filename="mesh/liver.msh")
+    liver.addObject('TetrahedronSetTopologyContainer', name="topo", src="@meshLoader")
+    liver.addObject('MechanicalObject', name="dofs", src="@meshLoader")
+    liver.addObject('TetrahedronSetGeometryAlgorithms', template="Vec3d", name="GeomAlgo")
+
+    # You can create values that have a dimension by multiplying a float/int by a unit
+    liver.addObject('TetrahedralCorotationalFEMForceField', template="Vec3d", name="FEM", method="large", poissonRatio="0.3", youngModulus=SceneUnit(3 * kPa), computeGlobalMatrix="0")
+    
+
+    # Multiplications between 'DimenssionedUnit' is supported and will affect the final unit
+    liverVolume = 1.5 * dm**3 # 1L
+    liverMass = 1.5 * kg
+    liverDensity = liverMass/liverVolume 
+    # You can print the value, the unit will show
+    print(f"Liver density is {liverDensity}")
+    liver.addObject('DiagonalMass', name="Mass", massDensity=SceneUnit(liverDensity))
+
+    # The library is also compatible with numpy array 
+    # This would also work np.array([10, 1, 5 ]) * N/m
+    # or classical list (but with lists, the list multiplication will fail, you need to specify the unit for each member)
+    stiffness = np.array([10 * N/m, 1 * N/m, 5 * N/m])
+    liver.addObject('RestShapeSpringsForceField', name="WeakConstraint", points="3 39 64", stiffness=SceneUnit(stiffness))
+
+    visu = liver.addChild('Visu')
+    visu.addObject('OglModel', name="VisualModel", src="@../../LiverSurface")
+    visu.addObject('BarycentricMapping', name="VisualMapping", input="@../dofs", output="@VisualModel")
+
+    return root
+
+
+# Function used only if this script is called from a python environment
+if __name__ == '__main__':
+    main()