Extending OOFEM in Python

The python interface allows user to define custom derived classes in Python and inject them inside oofem. We first illustate the concept on defining user defined time function (called MyFunc) in Python. The class should be derived from corresponding oofem base class (Function). It is necessary to implement at least methods to evaluate the function value (evaluateAtTime) and its derivatives (evaluateVelocityAtTime and evaluateAccelerationAtTime).

import oofempy
class MyFunc(oofempy.Function):
    def __init__(self, num, domain):
        oofempy.Function.__init__(self, num, domain)
    def evaluateAtTime (self, t):
        return 2.5*t
    def evaluateVelocityAtTime(self, t):
        return 0.0;
    def evaluateAccelerationAtTime(self, t):
        return 0.0;

After the definition, we can directly use the MyFunc class in the problem setup:

problem = oofempy.linearStatic(nSteps=1, outFile="test3.out")
domain = oofempy.domain(1, 1, problem, oofempy.domainType._1dTrussMode, tstep_all=True, dofman_all=True, element_all=True)
problem.setDomain(1, domain, True)
ltf1 = oofempy.peakFunction(1, domain, t=1, f_t=1)
ltf2 = MyFunc(2, domain) # use custom ltf here
ltfs = (ltf1, ltf2)
# boundary conditions
bc1   = oofempy.boundaryCondition(    1, domain, loadTimeFunction=1,    prescribedValue=0.0)
n2    = oofempy.nodalLoad(            2, domain, loadTimeFunction=2,    components=(1.,), dofs=(1,))
bcs = (bc1, n2)
#nodes
n1 = oofempy.node(1, domain, coords=(0.,  0., 0. ), bc=(bc1,))
n2 = oofempy.node(2, domain, coords=(2.4, 0., 0. ), load=(n2,))
nodes = (n1, n2)
# material and cross section
mat = oofempy.isoLE(1, domain, d=1., E=30.e6, n=0.2, tAlpha=1.2e-5)
cs  = oofempy.simpleCS(1, domain, area=0.5, Iy=0.0, beamShearCoeff=1.e18, thick=0.5)
# elements
e1 = oofempy.truss1d(1, domain, nodes=(1,n2),  mat=1,   crossSect=1)
elems = (e1,)
# add eveything to domain (resize container first to save some time, but it is not necessary 0 see ltfs)
domain.setup(nodes, elems, (mat,), (cs,), bcs, ltfs, (,))

problem.checkProblemConsistency()
problem.init()
problem.postInitialize()
problem.solveYourself()
problem.terminateAnalysis()
print("\nProblem solved")

Follow https://github.com/oofem/oofem/blob/master/bindings/python/tests/test_3.py for a full example.

Custom element

This section illustrates how to implement custom element in Python. In this case, we are implementing element for structural analysis, so it should be derived from cooresponding oofem base class, which is StructuralElement in this case. There are various options how to implement an element. The simplest option is to directly define how to evaluate element stifness matrix and internal force vector. Additional functuos to determine degrees of freedom required by the element anre also to be provided. Other possibilities exists, for example one can rely on default implementation of StructuralElement to evalaute stifness matrix, but then methods delivering strain-displacement and constituve matrix have to be provided together with element integration rules. For details, please consult oofem manuals.

import oofempy
class MyElement(oofempy.StructuralElement):
def __init__(self, num, domain):
    oofempy.StructuralElement.__init__(self, num, domain)
    self.setNumberOfDofManagers(2)
def computeStiffnessMatrix(self, answer, rMode, tStep):
    answer.resize(2,2);
    answer[0,0] = 1.0;
    answer[0,1] = -1.0;
    answer[1,0] = 1.0;
    answer[1,1] = -1.0;
def computeBmatrixAt(self, gp, answer, lowerIndx, upperIndx):
    answer.resize(1,2);
    answer[0,0]=-1.0;
    answer[0,1]=1.0;
def giveInternalForcesVector(self, answer, tStep, useUpdatedGpRecord):
    u = oofempy.FloatArray()
    k = oofempy.FloatMatrix(2,2)
    self.computeVectorOf(oofempy.ValueModeType.VM_Total, tStep, u)
    self.computeStiffnessMatrix(k, oofempy.CharType.StiffnessMatrix, tStep)
    answer.beProductOf (k, u)
def giveNumberOfDofs(self):
    return 2;
def computeNumberOfDofs(self):
    return 2;
def giveDofManDofIDMask(self, inode, answer):
    print ("giveDofManDofIDMask for %d"%(inode,))
    print (answer)
    answer.resize(1)
    answer[0] = oofempy.DofIDItem.D_u
    #answer.pY()
    print(answer)
def giveClassName(self):
    return "MyElement"
def giveInputRecordName(self):
    return "MyElement"

The element can again added into domain and used from oofem

# nodes
n1 = oofempy.node(1, domain, coords=(0.,  0., 0. ), bc=(bc1,))
n2 = oofempy.node(2, domain, coords=(2.4, 0., 0. ), load=(n2,))
nodes = (n1, n2)
# elements
e1 = MyElement(1, domain) # additinal entries should go to to the custom element constructor
e1.setDofManagers((1,2))
ir = oofempy.OOFEMTXTInputRecord()
ir.setRecordString ("nodes 2 1 2 mat 1 crosssect 1")
# pass input record to elem
e1.initializeFrom(ir)
elems = (e1,)
...

You can follow https://github.com/oofem/oofem/blob/master/bindings/python/tests/test_3.py for a complete illustration.

Custom material

This section illustrates how to implement custom constitutive model in Python. First we have to define Python class implementing the model derived from corresponding oofem class. In this case, we are going to implement constututive model for structural analysis, so we derive our class from StructuralMaterial. The presented implementation is mininimalistic one, we overload or define just methods to support 1d stress strain state, by overriding give1dStiffnessMatrix and giveRealStressVector_1d methods, that simply returns constitutive matrix and evaluate 1d stress from given strain. The more elaborate implementation would be necessary to support sevaral stress-strain modes, nonlinear materials (need to create custom status to track history variables). For more details, please refer to oofem programmer’s manual. Each matteril model must define its material status contatining its internal state variables. Even if material models does not need to track internal variables, it should provide status. In such case, it is sufficient to create instance of StructuralMaterialStatus. In oofem, the method responsible for status creation is CreateStatus, that would normally be overriden in Python as well. However, doe to the issue in Pybind11 (https://github.com/pybind/pybind11/issues/1962) this is not yet possible. The workaroud is to override giveStatus, as illustrate in the following example.

class MyMaterial(oofempy.StructuralMaterial):
  def __init__(self, num, domain):
      oofempy.StructuralMaterial.__init__(self, num, domain)
      self.k = 1.5;
  def giveClassName(self):
      return "MyMaterial"
  def giveInputRecordName(self):
      return "MyElement"
  # Overloading this method is not yet possible in pybind11
  # see https://github.com/pybind/pybind11/issues/1962
  # However a workaround is to override giveStatus, see below
  def CreateStatus(self, gp):
      return oofempy.StructuralMaterialStatus (gp)
  def give1dStressStiffMtrx(self, answer, mode, gp, tStep):
      answer.resize(1,1)
      answer[0,0] = self.k
      return
  def giveRealStressVector_1d (self, answer, gp, reducedStrain, tStep):
      answer.resize(1)
      answer[0] = self.k * reducedStrain[0]
      status = self.giveStatus(gp)
      status.letTempStrainVectorBe(reducedStrain);
      status.letTempStressVectorBe(answer);
      return
  def giveStatus (self, gp):
      print ("getStatus")
      if (gp.giveMaterialStatus() is None):
          print ("getStatus creating")
          status = oofempy.StructuralMaterialStatus (gp)
          gp.setMaterialStatus(status)
      return gp.giveMaterialStatus()

You can follow https://github.com/oofem/oofem/blob/master/bindings/python/tests/test_4.py for a complete illustration.