OOF2: The Manual

C++

All non-abstract C++ Property subclasses must have a constructor of the form

SubProperty::SubProperty(const std::string& name, PyObject *registration, ...) 

and they must invoke the base class constructor like this:

FluxProperty(const std::string& name, PyObject *registration) // or EqnProperty or AuxiliaryProperty

The derived class should treat name and registration as opaque variables, used only by the base class. They just need to be passed on through. The ... in the subclass constructor refers to the Property's parameters, which must appear in the constructor argument list in the same order that they appear in the corresponding PropertyRegistration.

C++ Property classes should not have a copy constructor. The base class copy constructor is private, to prohibit unauthorized duplication, by any means, electronic or otherwise.

Python

The Python constructor for a PyFluxProperty, PyEqnProperty, or PyAuxProperty subclass must look like this:

def __init__(self, name, registration, ...):
    PyFluxProperty.__init__(self, name, registration) # or PyEqnProperty or PyAuxProperty 

As in the C++ constructor, the name and registration parameters should be passed through to the base class unchanged.

The ... in the constructor arguments refers to any parameters that the Property might have. These must have the same names and be in the same order as the parameters in the associated PropertyRegistration. If the Property has no parameters, the derived class constructor does not need to be defined.

void cross_reference(Material* material)

cross_reference is called whenever a new Property is added to a Material. It gives each Property a chance to locate other Properties in the same Material from which it may need information. For example, many Properties of an anisotropic material will need to know the Material's Orientation Property. The argument to cross_reference is the Material object to which the Property belongs.

Material::fetchProperty should be used to look for other Properties. It will raise an exception if the sought-for Property does not exist. This exception should not be caught by the Property — it's handled in the Material.

For example, an anisotropic thermal expansion Property needs to access its Material's Elasticity and Orientation. It stores pointers to them in the class:

class ThermalExpansion : public FluxProperty {
protected:
  Elasticity *elasticity;
  ...
};

class AnisotropicThermalExpansion : public ThermalExpansion {
private:
  OrientationProp *orientation;
  ...
};    

and obtains their values by defining cross_reference() like this:

void AnisotropicThermalExpansion::cross_reference(Material *mat) {
  elasticity = dynamic_cast<Elasticity*>(mat->fetchProperty("Elasticity"));
  orientation = dynamic_cast<OrientationProp*>(mat->fetchProperty("Orientation"));
} 

In Python, it would look like this:

from ooflib.SWIG.engine import pypropertywrapper

class ThermalExpansion(pypropertywrapper.PyFluxProperty):
    def cross_reference(self, material):
        self.elasticity = material.fetchProperty("Elasticity")

class AnisotropicThermalExpansion(ThermalExpansion):
    def cross_reference(self, material):
        self.orientation = material.fetchProperty("Orientation") 

bool constant_in_space() const

constant_in_space() should return false if the Property explicitly depends on position. The default implementation returns true.

constant_in_space() can be used by other Properties to decide whether or not it's possible for them to precompute something. For example, an anisotropic elasticity Property can precompute its rotated stiffness modulus only if its orientation has no spatial variation.

void precompute(FEMesh* mesh)

precompute() is called when a Mesh is being solved. It's called before any other Property methods, other than cross_reference(). It should compute any quantities that the Property needs that depend only upon Property data. It can use Property data from other Properties, having located those other Properties during the cross_reference() step, but it should not assume that precompute() has been called for those Properties yet.

For example, if a Property has an anisotropic modulus, the Property's precompute() method should rotate the modulus to lab coordinates using the Material's Orientation (but only if the Orientation is constant in space.

void clear_mesh_data(FEMesh* mesh, void* data) const

clear_mesh_data() is called when data stored by set_mesh_data() is about to be overwritten, or if the FEMesh is about to be deleted. In the example given in set_mesh_data() above, clear_mesh_data() would simply delete the allocated pointer:

void MyProperty::clear_mesh_data(FEMesh* mesh, void* data) const {
    delete static_cast<MyMeshData*>(data);
} 

void begin_element(const CSubProblem*, const Element*)

When building the finite element stiffness matrix and other numerical components, OOF2 loops over the elements of the mesh and computes each element's contribution by calling its Properties' computational methods, e.g, FluxProperty::flux_matrix(), or EqnProperty::force_value(). Before calling these functions, however, it first calls begin_element(), passing the current Element as an argument. This allows the Property to precompute any Element dependent properties.

The base class implementation of begin_element() does nothing.

void end_element(const CSubProblem*, const Element*)

end_element is just like begin_element, except that it's called after OOF2 is done computing an Element's contribution. end_element() can be used to clean up any temporary objects allocated by begin_element().

The base class implementation of end_element() does nothing.

void post_process(CSubProblem*, const Element*) const

After solving the finite element equations, post_process is called once on each Element in a CSubProblem. It can be used for any sort of post-processing, such as adjusting material parameters, or generating output files.

Do not confuse post_process() with end_element(). end_element is called during finite element matrix construction, before the equations are solved.

bool is_symmetric_[KCM](const CSubProblem* mesh) const

is_symmetric_X indicates whether or not the finite element matrix constructed from the Property can be made symmetric by using the relationships established by ConjugatePair calls. There are three functions:

In most cases is_symmetric will simply return true, if the matrix can be symmetrized, or false, if it can't. Some cases are more complicated, however. For example, thermal expansion makes the matrix asymmetric if the temperature field is an active field, because the Property couples the temperature (with no derivatives) to the gradient of the displacement. Therefore, in the ThermalExpansion class, is_symmetric_K is defined like this:

bool ThermalExpansion::is_symmetric_K(const CSubProblem* subproblem) const {
  Equation *forcebalance = Equation::getEquation("Force_Balance");
  return !(forcebalance->is_active(subproblem) &&
	   temperature->is_defined(subproblem) &&
	   temperature->is_active(subproblem));
} 

The base class implementations of all three of these functions return true.

void output(FEMesh* mesh, const Element* element, const PropertyOutput* output, const MasterPosition& pos, OutputVal *data)

A Property's output function is called when quantities that depend on the Property are being computed, usually (but not necessarily) after a solution has been obtained on a mesh. Many different kinds of PropertyOutputs can be defined (see Section 8.7). Each Property class's PropertyRegistration indicates which kinds of PropertyOutput the Property can compute. The output function must determine which PropertyOutput is being computed, get the output's parameter values if necessary, compute its value at a given point in the mesh, and store the value in the given OutputVal object.

The arguments to Property::output in C++ are:

In Python, the arguments are the same, except that there's no data argument. Instead, the function returns an OutputVal (of the appropriate type) containing the Property's contribution to the output quantity.

Here is an example of a fairly complicated output function, from SRC/engine/properties/thermalexpansion/thermalexpansion.C in the OOF2 source code. It computes one of two types of output, Strain and Energy, each of which has two variants. It can also return the material parameters.

void ThermalExpansion::output(FEMesh *mesh,
                              const Element *element,
                              const PropertyOutput *output,
                              const MasterPosition &pos,
                              OutputVal *data)
  const 
{
  const std::string &outputname = output->name(); 
  if(outputname == "Strain") {
    // The parameter is a Python StrainType instance.  Extract its name.
    std::string stype = output->getRegisteredParamName("type"); 
    SymmMatrix3 *sdata = dynamic_cast<SymmMatrix3*>(data); 
    // Compute Temperature at position pos.
    double t = temperature->value(mesh, element, pos);  

    // Compute alpha*T and add it to the output.
    if(stype == "Thermal") 
      *sdata += expansiontensor(mesh, element, pos)*(t-T0);  
    else if(stype == "Elastic")
      *sdata -= expansiontensor(mesh, element, pos)*(t-T0);
  } // end outputname == "Strain"

  if(outputname == "Energy") { 
    // The parameter is a Python Enum instance.  Extract its name.
    std::string etype = output->getEnumParam("etype"); 
    if(etype == "Total" || etype == "Elastic") { 
      ScalarOutputVal *edata = dynamic_cast<ScalarOutputVal*>(data); 
      SymmMatrix3 thermalstrain;
      const Cijkl modulus = elasticity->cijkl(mesh, element, pos); 
      double t = temperature->value(mesh, element, pos);
      thermalstrain = expansiontensor(mesh, element, pos)*(t-T0);
      SymmMatrix3 thermalstress(modulus*thermalstrain);
      SymmMatrix3 strain;
      findGeometricStrain(mesh, element, pos, strain); 
      double e = 0;
      for(SpaceIndex i=0; i<3; i++) { 
        e += thermalstress(i,i)*(-strain(i,i) + 0.5*thermalstrain(i,i));
        SpaceIndex j = (i+1)%3;
        e += 2*thermalstress(i,j)*(-strain(i,j) + 0.5*thermalstrain(i,j));
      }
      *edata += e;
    }
  } // end outputname == "Energy"

  if(outputname == "Material Constants:Couplings:Thermal Expansion T0") {  
    ScalarOutputVal *sdata = dynamic_cast<ScalarOutputVal*>(data); 
    *sdata = T0;
  }
} 

temperature = dynamic_cast<ScalarField*>(Field::getField("Temperature")); 

PropertyRegistration(
  name="Couplings:ThermalExpansion:Isotropic",
  params = [
       parameter.FloatParameter("T0", 0.0, tip="Reference temperature"),
       parameter.FloatParameter("alpha", 1.0, tip="Expansion coefficient")],
  ...) 

void ThermalExpansion::cross_reference(Material *mat) {
  elasticity = dynamic_cast<Elasticity*>(mat->fetchProperty("Elasticity"));
} 

As an example in Python, here is output from SRC/engine/properties/elasticity/pyelasticity.py, which defines an elasticity Property in Python (for illustrative purposes only):

def output(self, mesh, element, output, pos):
    if output.name() == "Energy":
        etype = output.getEnumParam("etype")
        if etype in ("Total", "Elastic"):
            mod = self.modulus()
            # strain is a SymmMatrix3.  modulus is a cijkl.Cijkl
            strain = cstrain.getGeometricStrain(mesh, element, pos, False)
            stress = mod*strain # another SymmMatrix3.
            return outputval.ScalarOutputVal(0.5*stress.contract(strain)) 

Note that a new OutputVal is created and returned, containing the Property's contribution to the output quantity. Because energy is a scalar, the OutputVal is a ScalarOutputVal.

void set_mesh_data(const FEMesh* mesh, void *data) const

set_mesh_data() allows a Property to store arbitrary data in a FEMesh. It is meant to be called from within precompute() to cache mesh-dependent data that will be needed later. Mesh-dependent data cannot be stored directly in a Property because Properties can be shared between FEMeshes. The data will be accessible only to the Property that stored it.

For example, a Property called MyProperty could store its data in a class called MyMeshData and cache it in the mesh like this:

void MyProperty::precompute(FEMesh *mesh) {
    MyMeshData* data = new MyMeshData(mesh);  // some expensive computation is hidden here
    set_mesh_data(mesh, data);
} 

The Python implementation of set_mesh_data is completely independent of the C++ implementation. C++ Properties store data in C++, and Python Properties store data in Python.

Before data stored by set_mesh_data is overwritten or deleted, OOF2 will call clear_mesh_data. A Property subclass should define clear_mesh_data if it needs to do any cleaning up before the data is deleted.

void* get_mesh_data(const FEMesh* mesh) const

get_mesh_data() retrieves the data stored by set_mesh_data(). For example, if MyProperty in the example above is a FluxProperty, its flux_matrix method could retrieve the data like this:

void MyProperty::flux_matrix(const FEMesh* mesh, ...) const {
    MyMeshData* data = static_cast<MyMeshData*>(get_mesh_data(mesh));
    ...
}