The idea is to create a Python class (e.g., in the input file) that derives from the C++ class `PythonHyperMaterial` and overwrites the following functions:

• `def PythonHyperMaterial.init(self, _no):`, with `_no` being the material number.
• `def computeStressPython(self, FTotal):`
  • Input: The total deformation gradient.
    • Sent from C++ to Python through a `std::vector` [F11, F12, F13, F21, F22, F23, F31, F32, F33]
  • Output: Cauchy stress, in the current configuration.
    • Sent from Python to C++ through a Python list [Sig11, Sig22, Sig33, Sig12, Sig13, Sig23]

The relation between stress and strain is defined by the user.

• Why does `PythonHyperMaterial` derive from `HyperMaterial`?

  Hyper (Cauchy) Materials are based on the total deformation gradient. In the next merge, I will (try to) add Material perturbation to evaluate (finite difference) the Material Tangent Operator for all Cauchy Hypermaterial, including `PythonHyperMaterial`.

  :!: Therefore, so far, global perturbation (of the internal forces) must be used :!:

  ```python
  elementsProps.put(STIFFMETHOD, STIFF_NUMERIC)
  ```

• Why doesn't `ComputeStressPython` (i.e., the user-defined material) use Matr3/… etc?

  This user-defined feature should be kept **as simple as possible** (KISS). Therefore, the user can rely on simple data types which they can easily transform, for instance, to `np.array` (see example). A student has already used this feature without any questions after receiving the example below. If the application requires an efficient implementation, adding the material directly in Metafor remains the preferred procedure (especially with its tangent moduli).

```python class UserMat(PythonHyperMaterial):

  def __init__(self, _no):
      PythonHyperMaterial.__init__(self, _no)

  def computeStressPython(self, FTotal):
      C1 = 0.5
      K = 1000

      F = np.array(FTotal).reshape(3, 3)
      J = np.linalg.det(FTotal)

      Fbar = J**(-1/3) * F
      Bbar = np.dot(Fbar.T, Fbar)

      cauchyStress = C1/J * Bbar + (-1/3 * np.trace(Bbar) + K*(J-1)) * np.identity(3)
      return [cauchyStress[0, 0], cauchyStress[1, 1], cauchyStress[2, 2], cauchyStress[0, 1], cauchyStress[0, 2], cauchyStress[1, 2]]

def getMetafor(_p):

  # ... some code ... #
  
  myNeoHookean = UserMat(1)
  domain.getMaterialSet().add(myNeoHookean)

  # ... some code ... #

  elementsProps = ElementProperties(Volume2DElement)
  elementsProps.put(MATERIAL, 1)
  elementsProps.put(STIFFMETHOD, STIFF_NUMERIC)

  # ... some code ... #
  
  return metafor