OOF2: The Manual

Property
Prev  8.6.3. Property- and Material-related Classes  Next

Name

Property — Base class for material properties

Synopses

C++ Synopsis

#include "engine/property.h"
	  
class Property {
  Property(const std::string& name,
           PyObject* registration);
  virtual void cross_reference(Material* material);
  virtual void precompute();
  virtual void begin_element(const CSubProblem* subproblem,
                             const Element* element);
  virtual void begin_point(const FEMesh* mesh,
                           const Element* element,
                           const Flux* flux,
                           const MasterPosition& masterpos);
  virtual void end_point(const FEMesh* mesh,
                         const Element* element,
                         const Flux* flux,
                         const MasterPosition& masterpos);
  virtual void end_element(const CSubProblem* subproblem,
                           const Element* element);
  virtual void post_process(CSubProblem* mesh,
                            const Element* element) const;
  virtual int integration_order(const CSubProblem* mesh,
                                const Element* element) const;
  virtual void fluxmatrix(const FEMesh* mesh,
                          const Element* element,
                          const ElementFuncNodeIterator& node,
                          const Flux* flux,
                          const FluxData* fluxdata,
                          const MasterPosition& position) const;
  virtual void fluxrhs(const FEMesh* mesh,
                       const Element* element,
                       const Flux* flux,
                       FluxData* fluxdata,
                       MasterPosition& position) const;
  virtual void output(const FEMesh* mesh,
                      const Element* element,
                      const PropertyOutput* p_output,
                      const MasterPosition& position,
                      OutputVal* value) const;
  virtual bool is_symmetric(const CSubProblem* subproblem) const;
}

Python Synopsis

from oof2.SWIG.engine.pypropertywrapper import PyPropertyWrapper
	  
class PyPropertyWrapper(Property):
  def __init__(self, referent, registration, name)
  def integration_order(self, subproblem, element)
  def cross_reference(self, material)
  def precompute(self)
  def begin_element(self, subproblem, element)
  def begin_point(self, mesh, element, flux, master_position)
  def end_point(self, mesh, element, flux, master_position)
  def end_element(self, subproblem, element)
  def post_process(self, subproblem, element)
  def fluxmatrix(self, mesh, element, node, flux, fluxdata, master_position)
  def fluxrhs(self, mesh, element, flux, fluxdata, master_position)
  def output(self, mesh, element, p_output, position)
  def is_symmetric(self, subproblem)

Overview

Property classes can be constructed either in C++ or Python. Python classes are somewhat simpler to program and build, but C++ classes are significantly more efficient, computationally.

In C++, new Properties are constructed by creating subclasses of the Property class. In Python, they're constructed by creating subclasses of the PyPropertyWrapper class. (PyPropertyWrapper is a swig wrapper class around a C++ class of the same name, which is itself derived from Property.) In both cases, the base classes contain trivial implementations of all methods except integration_order, so it is safe to omit from the derived classes functions which are irrelevant to a particular physical property. The integration_order is always relevant to the physical property to which it belongs.

In order for both C++ and Python Properties to work properly and to be included in the GUI, every subclass must have an associated PropertyRegistration object. PropertyRegistrations must always be created in Python, even for C++ Property subclasses.

Constructors

C++

C++ Property subclasses must have a constructor of the form 

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

and they must invoke the base class constructor like this:[57] 

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

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 PyPropertyWrapper subclass must look like this: 

def __init__(self, registration, name, ...):
    PyPropertyWrapper.__init__(self, registration, name)

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.

Methods

The following methods can be defined in C++ Property and Python PyPropertyWrapper subclasses. Most are optional, and can be omitted if the default behavior is satisfactory. Unless otherwise noted, the default behavior is to do nothing. The exception is integration_order, which must be provided explicitly for all properties. A base-class version does exist, but its default behavior is to indicate an error.

All of the functions can be accessed from either C++ or Python. If there are non-obvious differences between the C++ and Python invocations, they are noted below. Minor and standard punctuation differences are not noted. Function prototypes are specified in the C++ format, because it is more informative.

void cross_reference(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, here is the cross_reference routine from the AnisotropicThermalExpansion class, which needs to access the Elasticity (for the elastic modulus) and the Orientation: 

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

The ThermalExpansion and AnisotropicThermalExpansion classes store the elasticity and orientation variables: 

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

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

The default base class cross_reference function does nothing.

void precompute()

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.

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.

void begin_element(const CSubProblem*, const const Element*)

When building the finite element stiffness matrix and right hand side vector, OOF2 loops over the elements of the mesh and computes each element's contribution by calling each Property's fluxmatrix and fluxrhs methods. 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.

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.

void begin_point(const FEMesh*, const const Element*, const Flux*, const MasterPosition&)

Within each element, OOF2 loops over gausspoints. Property objects can use this begin_point hook to perform expensive computations that are required at an evaluation point, and cache the results. There is a matching end_point function which can be used to clear the cached data.

void end_point(const FEMesh*, const Element*, const Flux*, const MasterPosition&)

end_point is just like begin_point, except that it's called after OOF2 is done computing the contribution at a point. The end_point method can be used to clean up any temporary objects allocated by begin_point.

void post_process(CSubProblem*, const Element*)

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.

int integration_order(const CSubProblem *subproblem, const Element *element) const

integration_order returns the polynomial order (or degree) of the part of the flux matrix (computed by fluxmatrix or fluxrhs) that this Property is responsible for. All Properties that define fluxmatrix or fluxrhs should define integration_order.

Because fluxmatrix and fluxrhs use the finite element shape functions and their derivatives, integration_order must find out the polynomial degree of the shapefunctions. The current mesh Element is passed in as an argument. The shapefunction's degree can be found by calling Element::shapefun_degree, and it's derivative's degree can be found by calling Element::dshapefun_degree.

[Warning] Warning

Do not compute the shapefunction's derivative's degree by subtracting 1 from the shapefunction's degree. For purposes of Gaussian integration, the degree of the shapefunction is sometimes less than its actual polynomial degree. For example, the linear quadrilateral shapefunction \((1-x)(1-y)/4\) can be integrated exactly with a single Gauss point at (0,0), although its polynomial degree is 2. For this function, both shapefun_degree and dshapefun_degree return 1.

As an example, here is the integration_order method for the Elasticity class in OOF2. The fluxmatrix routine always adds a (constant) modulus times a shape function derivative to the flux matrix, but when the displacement field has out-of-plane components, there are terms proportional to the shape functions as well. integration_order must return the largest relevant degree, so it has to check for out-of-plane fields: 

int Elasticity::integration_order(const CSubProblem *subproblem, const Element *el) const {
  if(displacement->in_plane(subproblem))
    return el->dshapefun_degree();
  return el->shapefun_degree();
}

void fluxmatrix(...)

OOF2 Properties that represent terms in a constitutive equation must define the function Property::fluxmatrix. fluxmatrix computes a matrix that, when multiplied by a vector of degrees of freedom, produces a vector containing the components of a Flux. Explaining this properly requires a brief review of the finite element machinery.

Let \(u_n(r)\) be the nth component of a field \(u\) at position \(r\). If the field values are known at nodes \(\nu\) at positions \(r_\nu\), then the finite element shape functions \(N_\nu(r)\) can be used to approximate the field at any point:

\[u_n({\bf r}) = \sum_\nu N_\nu({\bf r}) u_{n\nu}\] (8.2)

where \(N_\nu({\bf r})\) is the shape function that is 1 at node \(\nu\) and 0 at all other nodes, and \(u_{n\nu}\equiv u_n({\bf r_\nu})\).

A constitutive relation connects a Field, \(u\), and a Flux, \(\sigma\), through a modulus. For now, let's assume that it can be represented as a linear operator, \(M\):

\[\sigma({\bf r}) = M({\bf r})\cdot u({\bf r})\] (8.3)

Inserting (8.2) into (8.3) shows that \(M\) can be written as a matrix, which we call the flux matrix. The columns of \(M\) correspond to degrees of freedom \(u_{n\nu}\), and its rows to components of the flux \(\sigma\).

For example, for elasticity

\begin{align*} \sigma_{ij} &= C_{ijkl} \epsilon_{kl} \\ &= \frac12 C_{ijkl} \left(\frac{\partial u_l}{\partial r_k} + \frac{\partial u_k}{\partial r_l}\right) \\ &= C_{ijkl} \frac{\partial u_l}{\partial r_k} \\ &= C_{ijkl} \frac\partial{\partial r_k} \left(u_{l\nu}N_\nu(r)\right) \\ &= \left(C_{ijkl} \frac{\partial N_\nu(r)}{\partial r_k}\right) u_{l\nu} \end{align*} \\ (8.4)

from which we get

\[ M_{ij,l\nu} = C_{ijkl}\frac{\partial N_\nu(r)}{\partial r_k} \] (8.5)

Repeated indices are summed in all of the above equations. In particular, in (8.5), k runs only over the in-plane components x and y.

Actually, (8.5) isn't quite correct, because we ignored the out-of-plane components. The out-of-plane part of the displacement field is the vector of derivatives \( \partial u_z/\partial x, \partial u_z/\partial y, \partial u_z/\partial z \), which is just \( 2\epsilon_{xz}, 2\epsilon_{yz}, \epsilon_{zz} \) (using the fact that \(\partial u_x/\partial z\) and \(\partial u_y/\partial z\) must both be zero). The part of the flux \(\sigma\) due to the out-of-plane strains is

\[ \sigma_{ij} = \frac12\left[C_{ijxz}\frac{\partial u_z}{\partial x} +C_{ijyz}\frac{\partial u_z}{\partial y}\right] + C_{ijzz}\frac{\partial u_z}{\partial z} \] (8.6)

Since \(\partial u_z/\partial r_k\) is the \(k\) component of a field (an out-of-plane field, in the OOF2 sense — it's computed at nodes just like in-plane fields are), it can be expanded in terms of shape functions and the field values at nodes, as in (8.2). From this we see that the extra terms in \(M\) are:

\[ M'_{ij,k\nu} = \frac12 C_{ijzk}(1+\delta_{kz})N_\nu \] (8.7)

(Note: \(k\) is not summed in (8.7).)

On each call, fluxmatrix must make contributions to \(M\) for all of the components of a given flux (\(\sigma\)), all of the components of the relevant fields (\(u\) and its out-of-plane part), and one given node (\(\nu\)).

It is important to note that it is not necessary for the fluxmatrix routine or its author to know anything about the following:

  • The topology or number of nodes in the element.

  • How nodal degrees of freedom are mapped into columns of the flux matrix.

  • How flux components are mapped into rows of the flux matrix.

  • What equations are being solved.

The arguments to fluxmatrix are:

const FEMesh *mesh

The finite element mesh that's being solved. fluxmatrix probably doesn't have to use the mesh object directly, but it does have to pass it through to other functions.

const Element *element

The finite element under consideration. This shouldn't be explicitly needed except in cases in which the material parameters depend on physical space coordinates or in which history-dependent fields are stored in the element.

const ElementFuncNodeIterator &node

This is the node \(\nu\) referred to in the discussion above. Its passed in in the guise of an iterator, which can iterate over all of the nodes of the element, although it should be thought of here simply as a way of accessing the node's indices and shape functions.[58] The node's shape function and its derivatives can be obtained from ElementFuncNodeIterator::shapefunction and ElementFuncNodeIterator::dshapefunction.

const Flux *flux

The Flux \(\sigma\) referred to in the discussion above. Properties that contribute to more than one Flux need to check this variable to know which flux they're computing now.

FluxData *fluxdata

The FluxData class stores the actual flux matrix \(M\) and other flux data. (This data can't be stored directly in the Flux object because one Flux is shared among many FEMeshes. The FluxData object is local to this computation.) The actual flux matrix \(M\) can be accessed only through the method FluxData::matrix_element. This function handles the mapping from node, field, and flux component indices to actual row and column indices.

const MasterPosition &x

The flux matrix is evaluated at a physical coordinate \(r\) (often a Gauss integration point, but not always). x is the point in the Element's master coordinate space[59] corresponding to r.

Here is the fluxmatrix routine for the Elasticity class, from SRC/engine/property/elasticity/elasticity.C in the OOF2 source code: 

void Elasticity::fluxmatrix(const FEMesh *mesh, const Element *element,
			    const ElementFuncNodeIterator &nu,
			    const Flux *flux, FluxData *fluxdata,
			    const MasterPosition &x) const 
{
  if(*flux != *stress_flux) { 1
    throw ErrProgrammingError("Unexpected flux", __FILE__, __LINE__); 2
  }

  const Cijkl modulus = cijkl(mesh, element, x); 3
  double sf = nu.shapefunction(x); 4
  double dsf0 = nu.dshapefunction(0, x); 5
  double dsf1 = nu.dshapefunction(1, x); 

  for(SymTensorIterator ij; !ij.end(); ++ij) { 6
    for(IteratorP ell=displacement->iterator(); !ell.end(); ++ell) { 7

      SymTensorIndex ell0(0, ell.integer()); 8
      SymTensorIndex ell1(1, ell.integer());
      fluxdata->matrix_element(mesh, ij, displacement, ell, nu) += 9
	modulus(ij, ell0)*dsf0 + modulus(ij, ell1)*dsf1; 10
    }

    if(!displacement->in_plane(mesh)) { 11
      Field *oop = displacement->out_of_plane(); 12
      for(IteratorP ell=oop->iterator(ALL_INDICES); !ell.end(); ++ell) { 13
	double diag_factor = ( ell.integer()==2 ? 1.0 : 0.5); 14
	fluxdata->matrix_element(mesh, ij, oop, ell, nu) += 15
	  modulus(ij, SymTensorIndex(2, ell.integer())) * sf * diag_factor; 16
      }
    }
  }
}

1

We don't really need to do this. It's a sanity check to make sure that we got the Flux we wanted. stress_flux was set in the Elasticity constructor by calling Flux::getFlux("Stress") and dynamically casting the result to a SymmetricTensorFlux*. If a Property made contributions to more than one Flux, then a statement like this would be used to determine which flux it was being asked to compute.

2

ErrProgrammingError is declared in SRC/common/ooferror.h.

3

cijkl is a virtual function defined in the Elasticity subclasses. It takes FEMesh*, Element*, and MasterPosition& arguments so that a subclass can define a position-dependent modulus, or something even more bizarre.

4

Here is where the shape functions are evaluated. This is \(N_\nu\). Values are precomputed and cached at the Gauss integration points, so the evaluation here is fast.

5

These are the derivatives of the shape functions. The first argument is the component of the gradient. 0 is x and 1 is y. These values are also cached at Gauss points.

6

This is the loop over components of the flux. We know that the flux (stress) is a symmetric tensor, so we use a SymTensorIterator to loop over the components. If we didn't know the type, we would have instead written 

for(IteratorP ij = flux->iterator(); !ij.end(); ++ij) ...

The loop here extends over all components of the stress, without regard to whether or not the stress is in-plane. That's because the out-of-plane components of the flux matrix are used to construct the plane-stress constraint equation.

7

This is the loop over the in-plane components of the displacement field. In this case it's been written in the generic form. displacement is a TwoVectorField, with only x and y components, so we don't have to specify the planarity of the iterator. The associated out-of-plane field, containing the out-of-plane strains, will be treated separately.

8

This simply creates a SymTensorIndex object for addressing a component of the elastic modulus. \(C_{ijkl}\) requires two SymTensorIndex or SymTensorIterator objects to extract a single component.

9

This retrieves a C++ reference to the matrix element that couples flux component ij to the component ell of the field displacement at node nu. Note that we use += instead of = because a previous call to fluxmatrix may already have addressed this matrix element.

10

The two terms in this line are the implied sum over \(k\) in (8.5).

11

If the displacement field has no out-of-plane part, we don't need to compute the out-of-plane part's contributions to the stress.

12

This retrieves the Field for the out-of-plane part of the displacement.

13

This loops over the components of the out-of-plane field. Since this field has three components (\(\partial u_z/\partial x\), \(\partial u_z/\partial y\), and \(\partial u_z/\partial z\)), the iterator needs to be told to loop over all of them.

14

This takes care of the numerical factor in (8.7). All IteratorP objects can be converted to an integer with their integer function, which provides a handy way of checking their value if you know the meaning of the integer. In this case, 2 means z.

15

This line is just like 9, but it refers to the out-of-plane field, oop.

16

This is the right hand side of (8.7).

void Property::fluxrhs(...) const

OOF2 Properties that represent external (generalized) forces or otherwise contribute to the right hand side of a divergence equation must define Property::fluxrhs. fluxrhs is similar to Property::fluxmatrix in its role and its arguments, but is generally simpler to implement. Each fluxrhs implementation must compute a quantity at a given point within an element, but does not have to concern itself with nodes and shapefunctions.

There are two kinds of contributions to the right hand side of the divergence equation. Body forces contribute directly to the right hand side of the divergence equation (2.8), but offsets contribute a field-independent value to the flux on the left hand side of (2.8). (Field-dependent contributions are made by Property::fluxmatrix.)

For an example, consider linear thermal expansion. The flux (stress) is

\[ \sigma_{ij} = C_{ijkl} \left(\epsilon_{kl} - \alpha_{kl}(T-T_0)\right) \] (8.8)

\(T_0\) is the temperature at which the stress-free strain vanishes, and makes a field independent offset, \(C_{ijkl}\alpha_{kl}T_0\), to the flux. On the other hand, gravitational forces do not contribute to the flux, but appear as body forces:

\[\nabla\cdot\sigma = -g\hat{\mathrm{\bf y}}\] (8.9)

The arguments to Property::fluxrhs are:

const FEMesh *mesh

The finite element mesh that's being solved. fluxrhs probably doesn't have to use the mesh object directly, but it might need to pass it through to other functions.

const Element *element

The finite element under consideration. This shouldn't be explicitly needed except in cases in which the material parameters depend on physical space coordinates or in which history-dependent fields are stored in the element.

const Flux *flux

The Flux under consideration.

FluxData *fluxdata

The fluxdata object stores the results computed by fluxrhs. It's the same object that was passed to fluxmatrix. FluxData::offset_element accumulates offsets, and FluxData::rhs_element accumulates body forces.

const MasterPosition &x

Each call to fluxrhs evaluates the rhs or flux offset at a given point within the given Element. x is the master space[59] coordinate of the point.

Here is the fluxrhs function from the ForceDensity class, from SRC/engine/property/forcedensity/forcedensity.C in the OOF2 source code. It provides an example of a body force: 

void ForceDensity::fluxrhs(const FEMesh *mesh, const Element *element,
			      const Flux *flux, FluxData *fluxdata,
			      const MasterPosition &x) const 
{
  if(*flux != *stress_flux) { 1
    throw ErrProgrammingError("Unexpected flux", __FILE__, __LINE__); 2
  }
  
  fluxdata->rhs_element(0) -= gx; 3
  fluxdata->rhs_element(1) -= gy;
}

1

As in fluxmatrix, above, this is a sanity check.

2

ErrProgrammingError is declared in SRC/common/ooferror.h.

3

FluxData::rhs_element(i) returns a reference to the ith component of the divergence of the flux. Some other Property may already have made a contribution to it, so here its value is changed with -= instead of =. All of the components of the divergence should be addressed.

[Note] Note

As far as FluxData is concerned, the divergence of a symmetric 3×3 tensor flux, such as stress, has only 2 components! This is because the out-of-plane forces are irrelevant.[60] The function Flux::divergence_dim can be used to find the number of components of the divergence of a flux.

Here is the fluxrhs method from the ThermalExpansion property, which is an example of a flux offset. It computes the field-independent term

\(C_{ijkl}\alpha_{kl}T_0\) (8.10)

The original version can be found in SRC/engine/property/thermalexpansion/thermalexpansion.C in the OOF2 source code. 

void ThermalExpansion::fluxrhs(const FEMesh *mesh, const Element *element,
				   const Flux *flux, FluxData *fluxdata,
				   const MasterPosition &x) const {
  
  if(*flux!=*stress_flux) {
    throw ErrProgrammingError("Unexpected flux." __FILE__, __LINE__);
  }
  const Cijkl modulus = elasticity->cijkl(mesh, element, x); 1

  for(SymTensorIterator ij; !ij.end(); ++ij) { 2
    double &offset_el = fluxdata->offset_element(mesh, ij); 3
    for(SymTensorIterator kl; !kl.end(); ++kl) { 4
      if(kl.diagonal()) { 5
	offset_el += modulus(ij,kl)*expansiontensor[kl]*tzero_;
      }
      else {
	offset_el += 2.0*modulus(ij,kl)*expansiontensor[kl]*tzero_;
      }
    }
  }
}

1

This retrieves the elastic modulus \(C_{ijkl}\) from the Material's Elasticity property. The variable elasticity was set by this Property's cross_reference function, like this: 

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

2

This is the loop over components of the flux, just like item 6 in the fluxmatrix example.

3

FluxData::offset_element(i) returns a reference to a component of the flux offset. The variable i must be the appropriate kind of FieldIndex or FieldIterator object.

4

This is the loop over the indices \(kl\) in Eq. (8.10).

5

The next few lines compute terms in the sum, accumulating them in the FluxData object. expansiontensor is a SymmMatrix3 object stored in the ThermalExpansion object, containing the matrix \(\alpha_{kl}\). The loop 4 addresses only the independent symmetric tensor entries \(kl\), but the sum in (8.10) includes the other half of the tensor as well. Therefore, we need to include a factor of 2 for the off-diagonal terms. tzero_ is \(T_0\), which is set when the ThermalExpansion object is constructed.

void output(..)

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.5). 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 a given OutputVal object.

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

const FEMesh *mesh

The mesh on which the output is being computed. The output function will probably not need to use this variable directly, but must pass it through to other functions.

const Element *element

The element of the mesh containing the point at which output values are desired.

const PropertyOutput *output

The PropertyOutput object being computed. The object is created by the OOF2 menu system and, depending on the type of output, may contain Python arguments specifying exactly what's to be computed. For example, a strain output will specify whether it's computing the total, elastic, thermal or other variety of strain. See the example below.

const MasterPosition & pos

The position in the element's master coordinate space[59] at which the output is to be computed.

OutputVal *data

The object in which the computed data should be stored. In C++, the object must be first cast to an appropriate derived type. The new value should be added to data, to ensure that values computed by other Properties are retained.

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/property/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. 

void ThermalExpansion::output(const FEMesh *mesh,
				  const Element *element,
				  const PropertyOutput *output,
				  const MasterPosition &pos,
				  OutputVal *data)
  const 
{
  const std::string &outputname = output->name(); 1
  if(outputname == "Strain") {
    // The parameter is a Python StrainType instance.  Extract its name.
    std::string stype = output->getRegisteredParamName("type"); 2
    SymmMatrix3 *sdata = dynamic_cast<SymmMatrix3*>(data); 3
    // Compute alpha*T with T interpolated to position pos
    const OutputValue tfield = element->outputField(mesh, *temperature, pos); 4
    const ScalarOutputVal *tval = 
      dynamic_cast<const ScalarOutputVal*>(tfield.valuePtr()); 5
    double t = tval->value(); 6
    
    if(stype == "Thermal") 7
      *sdata += expansiontensor*(t-tzero_);
    else if(stype == "Elastic")
      *sdata -= expansiontensor*(t-tzero_);
  } // strain output ends here

  if(outputname == "Energy") {
    // The parameter is a Python Enum instance.  Extract its name.
    std::string etype = output->getEnumParam("etype"); 8
    if(etype == "Total" || etype == "Elastic") { 9
      ScalarOutputVal *edata = 
	dynamic_cast<ScalarOutputVal*>(data); 10
      SymmMatrix3 thermalstrain;
      const Cijkl modulus = elasticity->cijkl(mesh, element, pos); 11
      const OutputValue tfield = element->outputField(mesh, *temperature, pos); 12
      const ScalarOutputVal *tval = 
	dynamic_cast<const ScalarOutputVal*>(tfield.valuePtr());
      double t = tval->value();
      thermalstrain = expansiontensor*(t-tzero_);
      SymmMatrix3 thermalstress(modulus*thermalstrain);
      SymmMatrix3 strain;
      findGeometricStrain(mesh, element, pos, strain); 13
      double e = 0;
      for(int i=0; i<3; i++) { 14
	e += thermalstress(i,i)*strain(i,i);
	int j = (i+1)%3;
	e += 2*thermalstress(i,j)*strain(i,j);
      }
      *edata += -0.5*e;
    }
  } //energy output ends here
}

1

The name of the output indicates what type of quantity is to be computed.

2

This call to PropertyOutput::getRegisteredParamName retrieves the PropertyOutput's type parameter. Parameters for PropertyOutputs are defined in their PropertyOutputRegistrations. The PropertyOutput class provides a number of functions for retrieving parameters of different kinds (integer, float, string, etc.). In this case, type is a RegisteredParameter, so we use getRegisteredParamName to retrieve its name.

3

Because we know that this is a Strain output, we can cast the data pointer from its base class, OutputVal, to its derived class, SymmMatrix3, which holds a 3×3 symmetric matrix output value.

4

ThermalExpansion's contribution to the strain depends on the value of the Temperature Field. The variable temperature was set during the ThermalExpansion constructor, like this: 

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

This line calls Element::outputField to evaluate the temperature at the given position, pos. The return value is an OutputValue, which is a generic wrapper around different kinds of OutputVals. The generic wrapper allows different types of outputs to be handled by the same code elsewhere in OOF2. Here, however, it's a bit of a nuisance since we know that the temperature is a scalar field.

5

This line peels off the wrapper around the ScalarOutputVal contained within the OutputValue retrieved in 4. For different types of Field, make sure to use the correct OutputVal subclass.

6

This line extracts the actual temperature from the ScalarOutputVal. Whew![61]

7

These lines add the thermal expansion contribution to the strain, with the appropriate sign (depending on the type of strain being computed). expansiontensor is a SymmMatrix3 object which can be added directly to *sdata, which is the SymmMatrix3 cast of data.

8

The etype parameter for the Energy PropertyOutput is an EnumParameter instance, representing the EnergyType OOF2 Enum. The function PropertyOutput::getEnumParam retrieves the value of the parameter, and returns its name.

9

There are many different types of Energy output, but this Property only contributes to these two.

10

Energy is a scalar, so the OutputVal that was passed in is really a ScalarOutputVal. It must be cast to the derived class in order to be used.

11

See 1 in the fluxrhs discussion.

12

See 4, 5, and 6.

13

findGeometricStrain is defined in SRC/engine/cstrain.C and declared in SRC/engine/cstrain.h. It computes \(\epsilon_{ij}=0.5(\partial u_i/\partial x_j + \partial u_j/\partial x_i)\) where \(u\) is the displacement field.

14

This loop computes the contribution of thermal expansion to the elastic energy.

bool is_symmetric(const CSubProblem *mesh) const

is_symmetric indicates whether or not the finite element stiffness matrix constructed from the Property can be made symmetric by using the relationships established by ConjugatePair calls.

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 thermal expansion Property couples the temperature (with no derivatives) to the gradient of the displacement. Therefore, in the ThermalExpansion class, is_symmetric is defined like this: 

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

[57] Sorry about the inconsistent order of the arguments.

[58] It's not possible to use a Node instead of an ElementFuncNodeIterator here. Nodes are unable to evaluate shape functions because they don't know which Element is being computed. ElementFuncNodeIterators know which Element they're looping over.

[59] Elements are first defined in a master coordinate space, where geometry is easy, and then mapped to their actual positions in physical space. The master quadrilateral is a square of side 2 centered on the origin. The master triangle is a right isosceles triangle with vertices (0,0), (1,0), and (0,1). A master space coordinate can be converted to a physical point by Element::from_master.

[60] If the flux is in-plane, then \(\sigma_{iz}=0\) and the out-of-plane force must be zero too. If the flux is not in-plane, then \(\sigma_{iz}\) will be found during the solution process, and the out-of-plane forces can be computed. In both cases it's meaningless to specify the forces ahead of time.

[61] This method of obtaining a field value may seem baroque, especially for a scalar field, but it's necessary for preserving generality. We may implement a short-cut for scalar fields in the future.

Prev  Up  Next
Planarity  Home  PropertyRegistration