OOF2: The Manual

void* begin_point(...) const

void *begin_point(const FEMesh* mesh,
                  const Element* element,
                  const Equation* equation,
                  const MasterPosition& position,
                  double time) const 

When building the finite element matrices and vectors, OOF2 loops over the Gauss points within each Element. At each Gauss point, for each active Equation, it calls EqnProperty methods for each Property that contributes the Equation, beginning with begin_point(), followed by time_deriv_matrices(), force_deriv_matrix() (for nonlinear properties), and finally force_value(). To avoid recomputing quantities that are required by more than one of these methods, computations can be done in begin_point() and their results stored in a data structure. begin_point() should return this structure (as a void* pointer in C++ or an object in Python). The data structure will be passed through to the other FluxProperty methods. If the data object was allocated in C++, end_point() should delete it.

The arguments to begin_point() are:

The default implementation of begin_point() simply returns nullptr in C++ and None in Python.

void end_point(...) const

void end_point(const FEMesh* mesh,
               const Element* element,
               const Equation* equation,
               const MasterPosition& position,
               double time,
               void *data) const 

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.

The arguments for end_point() are the same as the arguments for begin_point(), with the addition of a void* which points to the data returned by begin_point().

The default implementation of end_point() does nothing.

void time_deriv_matrices(...) const

void time_deriv_matrices(const FEMesh* mesh,
                         const Element* element,
                         const ElementFuncNodeIterator& node,
                         const Equation* equation,
                         const MasterPosition& position,
                         double time,
                         void *data,
                         SmallSystem* linsys) const 

EqnProperty classes that contribute to $𝐂$ or $𝐌$, the coefficients of the first and second time derivatives of the Fields in the time-dependent divergence equation (2.9), need to define time_deriv_matrices(). It is analogous to FluxProperty::flux_matrix(), but generally much simpler.

As in FluxProperty::flux_matrix(), the derivation begins by expanding the Field components in terms of the shape functions:^[77]

where $N_{\nu }(𝐱)$ is the shape function that is 1 at node $\nu$ and 0 at all other nodes, and

is the $k^{\mathrm{th}}$ component of the field at node $\nu$. OOF2 stores the values of ${\phi }_{k\nu }$ for all nodes in a vector, where $k\nu$ should be understood as a single index.

Because ${\phi }_{k\nu }$ is a vector, Equation (2.9) is a matrix equation, so the first derivative term is

where $i$ is an equation index. Using Eq. (201), this becomes

The derivation for $𝐌$ is identical, except that the second time derivative is used instead of the first:

The job of time_deriv_matrices() is to compute the terms in square brackets on the right hand sides of (204) and (205).

The arguments to time_deriv_matrices() are:

The default version of time_deriv_matrices() does nothing.

void IsotropicDampingProp::time_deriv_matrices(
					   const FEMesh* mesh,
					   const Element* element,
					   const ElementFuncNodeIterator& node, 
					   const Equation* eqn,
					   const MasterPosition& mpos,
					   double time,
					   void *data, 
					   SmallSystem* linsys)
  const
{
  if(*eqn != force_balance)  
     throw ErrProgrammingError("Unexpected equation.", __FILE__, __LINE__);
  double shapeFuncVal = node.shapefunction(mpos); 
  for(IndexP component : *eqn->components()) { 
    eqndata->damping_matrix_element(component, displacement, component, node) 
      += coeff*shapeFuncVal; 
  }
} 

void MassDensityProp::time_deriv_matrices(
				       const FEMesh* mesh,
				       const Element* element,
				       const ElementFuncNodeIterator& node,
				       const Equation *eqn,
				       const MasterPosition& mpos,
				       double time,
				       void*,
				       SmallSystem *eqdata) const
{
  double shapeFuncVal = node.shapefunction(mpos);
  for(IndexP eqncomp : *eqn->components()) {
    eqdata->mass_matrix_element(component, disp, component, node) 
      += rho_ * shapeFuncVal;
  }
} 

void force_value(...) const

void force_value(const FEMesh* mesh,
                 const Element* element,
                 const Equation* equation,
                 const MasterPosition& position,
                 double time,
                 void* data,
                 SmallSystem* linsys) const 

A Property that makes a direct contribution to the force term in Eq. (2.9) must redefine force_value(), which must use SmallSystem::forceVector() to set the value of the force at the given point in the given Element. The arguments to force_value() are^[83]:

The default version of force_value() does nothing.

void ForceDensity::force_value(
                       const FEMesh *mesh,
                       const Element *element,
                       const Equation *eqn,
                       const MasterPosition &x,
                       double time, void *data,
                       SmallSystem *linsys) const
{
  linsys->forceVector() += g; 
} 

void force_deriv_matrix(...) const

 void force_deriv_matrix(const FEMesh* mesh,
                         const Element* element,
                         const ElementFuncNodeIterator& node,
                         const Equation* equation,
                         const MasterPosition& position,
                         double time,
                         void* data, 
                         SmallSystem* linsys) const 

force_deriv_matrix() needs to be defined for any EqnProperty that defines a nonlinear force_value().^[84] It computes the derivatives of the force with respect to the Fields. That is, it computes a ???

where $f_i$ is a force component and ${\phi }_j$ is a Field component. The result is stored in the given SmallSystem by calling SmallSystem::force_deriv_matrix_element(), which handles all of the indexing.

The arguments to force_deriv_matrix() are the same as the arguments to time_deriv_matrices().

The default version of force_deriv_matrix() does nothing.

void NonlinearForceDensity::force_deriv_matrix(const FEMesh *mesh,
					       const Element *element,
					       const ElementFuncNodeIterator &node,
					       const Equation *eqn,
					       const MasterPosition &point,
					       double time, void*,
					       SmallSystem *linsys) const
{
  DoubleVec fieldVal = displacement->values(mesh, element, point); 
  Coord coord = element->from_master(point); 
  SmallMatrix forceDeriv = nonlin_force_density_deriv(coord, time, fieldVal); 
  
  for(IndexP eqncomp : *eqn->components()) { 
    for (IndexP fieldcomp : *displacement->components(ALL_INDICES)) { 

      linsys->force_deriv_matrix_element(eqncomp, displacement, fieldcomp, node)
        += forceDeriv(eqno, fieldno); 
    }
  }
} 

displacement = dynamic_cast<TwoVectorField*>(Field::getField("Displacement"));