OOF2: The Manual

FluxProperty
	Property and Material Classes

Name

FluxProperty — A Property that contributes to a Flux

Synopses

C++ Synopsis

#include "engine/property.h"

class FluxProperty: , public PhysicalProperty {
  virtual void* begin_point(const FEMesh* mesh,
                            const Element element,
                            const Flux* flux,
                            const MasterPosition& position,
                            double time) const;
  virtual void end_point(const FEMesh* mesh,
                         const Element* element,
                         const Flux* flux,
                         const MasterPosition& position,
                         double time,
                         void* data) const;
  virtual void flux_matrix(const FEMesh* mesh,
                           const Element* element,
                           const ElementFuncNodeIterator& node,
                           const Flux* flux,
                           const MasterPosition& position,
                           double time,
                           void* data,
                           SmallSystem* linsys) const;
  virtual void flux_offset(const FEMesh* mesh,
                           const Element* element,
                           const Flux* flux,
                           const MasterPosition& position,
                           double time,
                           void* data,
                           SmallSystem* linsys);
  virtual void flux_value(const FEMesh* mesh,
                          const Element* element,
                          const Flux* flux,
                          const MasterPosition& position,
                          double time,
                          void* data,
                          SmallSystem* linsys) const;
  virtual void static_flux_value(const FEMesh* mesh,
                                 const Element* element,
                                 const Flux* flux,
                                 const MasterPosition& position,
                                 double time,
                                 void* data,
                                 SmallSystem* linsys) const;
}

Python Synopsis

PyFluxProperty is a swigged C++ class that is derived (in C++) from FluxProperty. It is used as a base class for FluxProperties that are defined in Python. When the methods listed below are called from C++, they use the Python C API to invoke methods in the Python derived class.

from ooflib.SWIG.engine import pypropertywrapper

class PyFluxProperty(FluxProperty):
  def begin_point(self, mesh, element, flux, position, time)
  def end_point(self, mesh, element, flux, position, time, data)
  def flux_matrix(self, mesh, element, node, flux, position, time, data, linsys)
  def flux_offset(self, mesh, element, flux, position, time, data, linsys)
  def flux_value(self, mesh, element, flux, position, time, data, linsys)
  def static_flux_value(self, mesh, element, flux, position, time, data, linsys)

Source Files

SRC/engine/property.C: C++ code
SRC/engine/property.h: C++ header
SRC/engine/property.swg: SWIG source code
SRC/engine/property.swg: Python code included in the SWIG output
SRC/engine/pypropertywrapper.C: C++ code for Properties implemented in Python
SRC/engine/pypropertywrapper.h: C++ header
SRC/engine/pypropertywrapper.swg: SWIG source code
SRC/engine/pypropertywrapper.spy: Python code inserted into the SWIG output.

Overview

A FluxProperty is Property that contributes to a Flux. For example, Elasticity contributes to Stress, and ThermalConductivity contributes to Heat_Flux. The principal job of the FluxProperty is to compute the coefficients of the Fields and their time derivatives in a linear expansion of a Flux, $\sigma$ in terms of a Field or Fields, $\varphi$ , and its spatial and time derivatives:

\sigma=\sigma^{0}({\bf x})+{\cal A}({\bf x})\varphi+{\cal K}({\bf x})\nabla% \varphi+{\cal D}({\bf x})\frac{\partial\varphi}{\partial t}+{\cal C}({\bf x})% \frac{\partial\nabla\varphi}{\partial t}

(187)

$\sigma^{0}$ is an field-independent offset, to first order. The coefficients are written as functions of position, ${\bf x}$ , but $\sigma^{0}$ , $\mathcal{A}$ and $\mathcal{K}$ may be nonlinear functions of the Fields and their derivatives. A future version of OOF2 may allow nonlinearities in $\mathcal{C}$ and $\mathcal{D}$ as well.

A new subclass of FluxProperty must redefine at least one of the following four virtual methods:

FluxProperty::flux_matrix() computes the linearization of the Flux at a given point in an Element, as a function of the relevant Fields, their gradients, and possibly their time derivatives. That is, it determines the coefficients $\mathcal{A}$ , $\mathcal{K}$ , $\mathcal{D}$ , and $\mathcal{C}$ in Eqn. (187).
FluxProperty::flux_offset() computes the value of the Flux at a given point in an Element when the Fields and their gradients are zero. This is $\sigma^{0}$ in (187).
FluxProperty::flux_value() computes the actual value of the Flux at a given point in an Element. If the Property is linear in the Fields, flux_value() does not need to be redefined, since the Flux value can be deduced from flux_matrix() and flux_offset().
static_flux_value() computes the part of the Flux value that depends only on the Fields and not on their time derivatives, if the Flux depends on the Field derivatives.

The functions listed above all have a SmallSystem* argument, which represents a small part of the LinearizedSystem that is currently being built. SmallSystem methods handle all of the issues related to matrix indexing, element order, and so forth. The derived FluxProperty methods (and the developer) do not need to know anything about what equations are being solved or how they are discretized by the finite element machinery. Details on all of the functions are provided below.

In addition, two virtual utility functions, begin_point() and end_point() can be redefined to facilitate passing data between the other routines.

Constructor

The C++ and Python constructors for FluxProperty are the same as the base class constructors.

Methods

`void* begin_point(...) const`

void* begin_point(const FEMesh* mesh,
                  const Element* element,
                  const Flux* flux,
                  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, and for each active Flux, it calls FluxProperty methods for each Property that contributes to the Flux, beginning with begin_point(), followed by flux_matrix(), flux_offset(), and (sometimes) flux_value(), and ending with end_point(). 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 data returned by begin_point() is not stored or transmitted to any parts of OOF2 other than the FluxProperty methods discussed here. If no data needs to be transmitted, begin_point() should return nullptr in C++ or None in Python.

The arguments to begin_point() are:

const FEMesh* mesh: The FEMesh currently being solved.
const Element* element: The Element containing the current point.
const Flux* flux: The Flux that this Property computes. If the Property computes more than one Flux, begin_point will be called once for each Flux at each point.
const MasterPosition& point: The position within the Element where the Property is being evaluated. The position is given in master coordinates, which can be converted to physical coordinates by Element::from_master() if necessary.
double time: The time at which the Property is being evaluated.

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 Flux* flux,
               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 flux_matrix(...) const`

void flux_matrix(const FEMesh* mesh,
                 const Element* element,
                 const ElementFuncNodeIterator& node,
                 const Flux* flux,
                 const MasterPosition& position,
                 double time,
                 void* data,
                 SmallSystem* linsys) const

FluxProperty subclasses that represent the $\mathcal{A}$ , $\mathcal{K}$ , $\mathcal{D}$ or $\mathcal{C}$ terms in the constitutive equation (187) must redefine the virtual function FluxProperty::flux_matrix(). flux_matrix() computes a matrix that, when multiplied by a vector of degrees of freedom and/or their time derivatives, produces a vector containing the components of a Flux. Explaining this properly requires a brief review of the finite element machinery.

Let $\varphi_{k}({\bf x})$ be the $k^{\mathrm{th}}$ component of a field $\varphi$ at position ${\bf x}$ . If the field values are known at nodes $\nu$ , which are at positions ${\bf x}_{\nu}$ , then the finite element shape functions $N_{\nu}({\bf x})$ can be used to approximate the field at ${\bf x}$ :^[77]

\varphi_{k}({\bf x})=N_{\nu}({\bf x})\varphi_{k\nu}

(188)

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

\varphi_{k\nu}\equiv\varphi_{k}({\bf x_{\nu}})

(189)

is the $k^{\mathrm{th}}$ component of the field at node $\nu$ . Similarly, the gradient of the field is

\nabla\varphi=\frac{\partial\varphi_{k}({\bf x})}{\partial{\bf x}}=\frac{% \partial N_{\nu}({\bf x})}{\partial{\bf x}}\varphi_{k\nu}.

(190)

Because OOF2 does its computations in 2D, gradients include only x and y derivatives, unless explicitly stated otherwise.

Temporarily ignoring the $\partial\varphi/\partial t$ terms, a constitutive relation connects a Field $\varphi$ and its gradient, $\nabla\varphi$ , to a Flux $\sigma$ , through moduli $\mathcal{A}$ and $\mathcal{K}$ , and an offset, $\sigma^{0}$ :

\sigma=\sigma^{0}+\mathcal{A}\cdot\varphi+\mathcal{K}\cdot\nabla\varphi.

(191)

This is the linearization of a more general relation. $\mathcal{A}$ , $\mathcal{K}$ and $\sigma^{0}$ may all be functions of ${\bf x}$ , $\varphi$ , $\nabla\varphi$ , or time. ( $\sigma^{0}$ must not have any linear dependence on $\varphi$ or $\nabla\varphi$ , though, because then it would be part of $\mathcal{A}$ or $\mathcal{K}$ .)

(191) can be written as a matrix equation:

\sigma_{i}=\sigma^{0}_{i}+\mathcal{A}_{ik}\cdot\varphi_{k}+\mathcal{K}_{ijk}% \cdot\nabla_{j}\varphi_{k}

(192)

(repeated indices are summed). The index $i$ is a flux component, $j$ is a space component, and $k$ is a field component. Inserting (188) and (190) into (192)

	$\displaystyle\sigma_{i}$	$\displaystyle=\mathcal{A}_{ik}N_{\nu}({\bf x})\varphi_{k\nu}+\mathcal{K}_{ijk}% \frac{\partial N_{\nu}}{\partial{\bf x}_{j}}\varphi_{k\nu}+\sigma^{0}_{i}$
		$\displaystyle=\left[\mathcal{A}_{ik}N_{\nu}({\bf x})+\mathcal{K}_{ijk}\frac{% \partial N_{\nu}}{\partial{\bf x}_{j}}\right]\varphi_{k\nu}+\sigma^{0}_{i}$

(193)

shows that the vector of flux values $\sigma$ at ${\bf x}$ is related to the vector of field values $\varphi$ at the nodes by a matrix we call the flux matrix:

{\bf K}_{ik\nu}=\mathcal{A}_{ik}({\bf x})N_{\nu}+\mathcal{K}_{ijk}\frac{% \partial N_{\nu}}{\partial{\bf x}_{j}}

(194)

and an offset $\sigma^{0}$ . The columns of $\bf K$ correspond to field values $\varphi_{k\nu}$ and the rows to flux components $\sigma_{i}({\bf x})$ . $\bf K$ is a part of what is generally called the stiffness matrix.

Computing $\bf K$ (194) is the purpose of FluxProperty::flux_matrix(). On each call, flux_matrix() must make contributions to $\bf K$ for all of the components of a given flux, all of the components of the relevant fields $\varphi$ (and their out-of-plane parts), and one given node $\nu$ .

The SmallSystem object that is passed in to flux_matrix is a container for the various matrices that Properties compute. flux_matrix() contributes to $\bf K$ by calling SmallSystem::stiffness_matrix_element(), as shown in the example below.

Similarly, if the Flux depends on the time derivative of the Field, $\partial\varphi/\partial t$ , or its gradient, $\partial\nabla\varphi/\partial t$ , through moduli $\mathcal{D}$ and $\mathcal{C}$ respectively, then flux_matrix() must also compute a damping matrix:

{\bf C}_{ik\nu}=\mathcal{D}_{ik}({\bf x})N_{\nu}+\mathcal{C}_{ijk}({\bf x})% \frac{\partial N_{\nu}}{\partial{\bf x}_{j}}

(195)

It stores its results by calling SmallSystem::damping_matrix_element(). OOF2 currently does not allow nonlinearities in the time derivative terms, so $\mathcal{D}$ and $\mathcal{C}$ cannot depend on the Fields, but they can be functions of position.

Note that it is not necessary for the flux_matrix() 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 matrix.
How flux components are mapped into rows of the matrix.
What equations are being solved.

The arguments to flux_matrix() are:

const FEMesh* mesh: The finite element mesh that's being solved. flux_matrix() probably doesn't have to use the mesh object directly, but it does 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 if history-dependent fields are stored in the element.
const ElementFuncNodeIterator& node: This is the node $\nu$ referred to in the discussion above. It's passed in in the form 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.^[78] The node's shape function and its derivatives can be obtained from ElementFuncNodeIterator::shapefunction() and ElementFuncNodeIterator::dshapefunction().
const Flux* flux: The Flux $\sigma$ . Properties that contribute to more than one Flux need to check this variable to know which flux they're computing now.
const MasterPosition& position: This is the position ${\bf x}$ in (193). position is a point in element's Element's master coordinate space. It is often a Gauss integration point, but not always. Master coordinates can be converted to physical coordinates by calling Element::from_master().
double time: The time at which the flux matrix is being evaluated.
void* data: The data object, if any, returned by begin_point().
SmallSystem* linsys: The SmallSystem class stores the portion of the the actual flux matrix $\bf K$ and damping matrix $\bf C$ that are currently being computed, along with other flux data.^[79]

Example of `flux_matrix()`

For example, for thermal conductivity the field $\varphi$ is the temperature T, and the heat flux is the vector ${\bf J}$ :

{\bf J}_{i}({\bf x})=-\kappa_{ij}\nabla_{j}T({\bf x})-\kappa_{iz}T_{z}

(196)

The z component in (196) has been separated out because, as explained in Section 2.5.5, the z derivatives of Fields are treated differently from the x and y derivatives in OOF2. $T_{z}\equiv\partial T/\partial z$ is treated as a separate scalar Field, and is only present if Temperature is not in-plane. Because gradients of $T_{z}$ don't appear in (196), it is treated as part of the $\mathcal{A}$ term in (194).

Expanding the fields and derivatives in terms of the shape functions via (188) and (190),

	$\displaystyle{\bf J}_{i}$	$\displaystyle=-\kappa_{ij}\frac{\partial}{\partial{\bf x}_{j}}(T_{\nu}N_{\nu}(% {\bf x}))-\kappa_{iz}T_{z\nu}N_{\nu}({\bf x})$
		$\displaystyle=-\left(\kappa_{ij}\frac{\partial N_{\nu}({\bf x})}{\partial{\bf x% }_{j}}\,\right)T_{\nu}-\left(\kappa_{iz}N_{\nu}({\bf x})\right)T_{z\nu}$

(197)

from which we get

{\bf K}_{i\nu}=-\kappa_{ij}\frac{\partial N_{\nu}({\bf x})}{\partial{\bf x}_{j}}

(198)

for $T$ and

{\bf K}_{z\nu}=-\kappa_{iz}N_{\nu}({\bf x})

(199)

for $T_{z}$ . Note that there is no $k$ index because Temperature is a scalar.

Here is the fluxmatrix routine for the HeatConductivity class, from SRC/engine/properties/heatconductivity/heatconductivity.C in the OOF2 source code:

void HeatConductivity::flux_matrix(const FEMesh  *mesh,
                                   const Element *el,
                                   const ElementFuncNodeIterator &nu,
                                   const Flux *flux,
                                   const MasterPosition &position,
                                   double time,
                                   void* data,
                                   SmallSystem *linsys) const
{
  if (*flux != *heat_flux) { 
    throw ErrProgrammingError("Unexpected flux", __FILE__, __LINE__); 
  }

  double sf   = nu.shapefunction(position); 
  double dsf0 = nu.dshapefunction(0, position); 
  double dsf1 = nu.dshapefunction(1, position); 

  const SymmMatrix3 cond(conductivitytensor(mesh, el, position)); 

  for(IndexP i : *flux->components(ALL_INDICES)) { 
     
    // in-plane temperature gradient contributions
    linsys->stiffness_matrix_element(i, temperature, nu) -=  
        cond(i, 0) * dsf0 +  
        cond(i, 1) * dsf1;   

    // out-of-plane temperature gradient contribution
    if(!temperature->in_plane(mesh)) 
      
      linsys->stiffness_matrix_element(i, temperature->out_of_plane(), nu) 
              -= cond(i, 2) * sf; 
  }
}

	We don't really need to do this. It's a sanity check to make sure that we got the `Flux` we wanted. `heat_flux` is a protected data member of `HeatConductivity` and was set in the constructor by calling `Flux::getFlux("Heat_Flux")` and dynamically casting the result to a `VectorFlux*`. If a Property makes contributions to more than one `Flux` then a statement like this must used to determine which flux it is being asked to compute.
	`ErrProgrammingError` is declared in `SRC/common/ooferror.h`.
	Here is where the shape functions are evaluated. This is $N_{\nu}$ . Shape function values are precomputed and cached at the Gauss integration points, so the evaluation here is fast.
	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.
	`conductivity_tensor()` is a virtual function defined in the `HeatConductivity` subclasses. It takes `FEMesh`, `Element`, and `MasterPosition&` arguments so that a subclass can define a position-dependent modulus, or something even more bizarre.
	Note that this loop is over all components of the `Flux`. When solving a plane-flux problem, the out-of-plane coefficients in $\bf K$ are used to generate the constraint equations. Looping over the components of a `Field` or `Flux` is discussed in Section 8.4.
	If the field weren't a scalar, there would be a loop over its components here.
	`temperature` is a variable stored in the `HeatConductivity` class and initialized by calling `Field::getField("Displacement")` in the `HeatConductivity` constructor.
	`SmallSystem::stiffness_matrix_element()` retrieves a C++ reference to the matrix element that couples flux component `i` to the `temperature` field at node `nu`. Note that we use `-=` instead of `=` because a previous call to `fluxmatrix` may already have addressed this matrix element.
	Here the `IndexP`, `i`, is implicitly converted to an integer when passed to `SymmMatrix3::operator(int, int)`.
	This is (198), summing over $j$ .
	This checks to see if the temperature field has an out-of-plane part. If it does, its contribution to the stress is computed.
	This retrieves a reference to the matrix element coupling the out-of-plane part of the Temperature `Field` to flux component $i$ .
	This is (199).

`void flux_offset(...) const`

void flux_offset(const FEMesh* mesh,
                 const Element* element,
                 const Flux* flux,
                 const MasterPosition& position,
                 double time,
                 void* data,
                 SmallSystem *linysys) const

FluxProperty::flux_offset() computes the $\sigma^{0}$ term in (187). It is the value of the linearized Flux when the Fields are zero. It is similar in purpose to FluxProperty::flux_matrix(), but generally simpler to implement.

Because $\sigma^{0}$ is not proportional to a Field or its gradient, flux_offset() does not have to worry about expansions and shape functions. It just computes the value of $\sigma^{0}$ and stores it in the given SmallSystem.^[80]

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

\sigma_{ij}=C_{ijkl}\left(\epsilon_{kl}-\alpha_{kl}(T-T_{0})\right).

(200)

$\alpha_{kl}$ is the thermal expansion coefficient and $T_{0}$ is the temperature at which the thermal stress vanishes. The $C_{ijkl}\left(\epsilon_{kl}-\alpha_{kl}\,T\right)$ term in (200) is handled by ThermalExpansion::flux_matrix() because of its field dependencies, but $C_{ijkl}\alpha_{kl}\,T_{0}$ makes a field independent contribution to the flux, and needs to be handled by flux_offset().

The arguments to FluxProperty::flux_offset() are:

const FEMesh* mesh: The finite element mesh that's being solved. flux_offset() 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 Flux* flux: The Flux $\sigma$ . Properties that contribute to more than one Flux need to check this variable to know which flux they're computing now.
const MasterPosition& position: The flux offset is evaluated at a physical coordinate ${\bf x}$ (often a Gauss integration point, but not always). position is the point in the Element's master coordinate space.
double time: The time at which the flux offset is being evaluated.
void* data: The data object, if any, returned by begin_point().
SmallSystem* linsys: The SmallSystem class stores the portion of the finite element matrices and vectors that are currently being computed. The results of the flux_offset() calculation are stored here.

Example of `flux_offset`

Here is the flux_offset routine for the ThermalExpansion class, from SRC/engine/properties/thermalexpansion/thermalexpansion.C:

void ThermalExpansion::flux_offset(const FEMesh* mesh,
                   const Element* element,
                   const Flux* flux,
                   const MasterPosition& x,
                   double time,
                   void* data, 
                   SmallSystem* linsys)
  const
{

  if(*flux!=*stress_flux) { 
    throw ErrProgrammingError("Unexpected flux.", __FILE__, __LINE__);
  }
  const Cijkl modulus = elasticity->cijkl(mesh, element, x); 
  SymmMatrix3 expten = expansiontensor(mesh, element, x); 
  linsys->offsetVector() += modulus*expten*T0; 
}

	This is the object returned by `begin_point()`.
	This is optional here, to check that we got the `Flux` we were expecting. If a `FluxProperty` contributes to more than one `Flux`, then a statement like this can be used to determine which one is being computed. The variable `stress_flux` was set by the `ThermalExpansion` constructor by calling `Flux::get_Flux("Stress")` and dynamically casting the result to a `SymmetricTensorFlux*`.
	This extracts the elastic modulus from the `Material`'s `Elasticity` `Property`. `elasticity` is a variable set in `ThermalExpansion::cross_reference()` by calling `Material::fetchProperty("Elasticity")`.
	This is $\alpha_{kl}$ in (200). `expansiontensor()` is a virtual function in the `ThermalExpansion` class. The various isotropic and anisotropic subclasses define `expansiontensor()` in different ways.
	`SmallSystem::offsetVector()` returns a reference to flux offset vector $\sigma^{0}$ . This line increments it according to (200). It is important to increment it instead of setting it, because some other `Property` may already have made a contribution to the same vector. $\sigma^{0}$ is stored as a vector even though the Stress is a tensor. It's treated as a list of values in an order determined by the `integer()` method of the `SymTensorIndex` class.

The example above uses C++ vectors and tensors. An equivalent way of writing it explicitly loops over their components:

void ThermalExpansion::flux_offset(const FEMesh* mesh,
                   const Element* element,
                   const Flux* flux,
                   const MasterPosition& x,
                   double time,
                   void* data,
                   SmallSystem* linsys) const {

  if(*flux != *stress_flux) { 
    throw ErrProgrammingError("Unexpected flux.", __FILE__, __LINE__);
  }
  const Cijkl modulus = elasticity->cijkl(mesh, element, x); 
  SymmMatrix3 expten = expansiontensor(mesh, element, x);

  for(IndexP ij : *flux->components(ALL_INDICES)) { 
    double &offset_el = linsys->offset_vector_element(ij); 
    for(SymTensorIndex kl : symTensorIJComponents) { 
      if(kl.diagonal()) { 
	offset_el += modulus(ij,kl)*expten[kl]*T0;
      }
      else {
	offset_el += 2.0*modulus(ij,kl)*expten[kl]*T0;
      }
    }
  }
}

	This loop is over all components of the stress, even if the problem is being solved in plane stress. The out-of-plane components are used in the constraint equations. Looping over components is discussed in Section 8.4.
	`offset_el` is a reference to the `ij` component of $\sigma^{0}$ in the `SmallSystem` object. This value will be incorporated into the global vectors constructed by `SmallSystem`.
	This loops over all components of the stress. The line could have been written for(IndexP kl : flux->components(ALL_INDICES)) { ... } but the generic `IndexP` would have to be cast into a `SymTensorIndex` in order to call `SymTensorIndex::diagonal()` on the next line.^[81] Using `symTensorIJComponents` (note the lower-case “s”!) avoids this complication.
	The loop over the components of the symmetric tensor flux hits each independent component once, but when multiplying tensors, the off-diagonal components occur twice, leading to a factor of two in the off-diagonal terms. This complication can be avoided by using the tensor multiplication operators directly, as in the simpler example above.

`void flux_value(...) const`

void flux_value(const FEMesh* mesh,
                const Element* element,
                const Flux* flux,
                const MasterPosition& position,
                double time,
                void* data,
                SmallSystem *linysys) const

flux_value() computes the value of a given Flux, at a given point in an Element, as in (187). It is used when solving nonlinear equations and when computing flux-dependent Outputs.

The default base class version of flux_value() uses flux_property(), flux_offset(), and the Field values to evaluate (187). A subclass of FluxProperty should redefine flux_value() if for some reason it is faster or more accurate to compute the Flux directly than it is to use (187). If the Flux is a nonlinear function of the Fields, it should always define flux_value().

The arguments to flux_value() are:

const FEMesh* mesh: The FEMesh currently being solved or evaluated.
const Element* element: The Element containing the current point.
const Flux* flux: The Flux that is being computed.
const MasterPosition& position: The position within the Element where the Property is being evaluated. The position is given in master coordinates, which can be converted to physical coordinates by Element::from_master() if necessary.
double time: The time at which the Flux is being computed.
void* data: The data object, if any, returned by begin_point().
SmallSystem* linsys: The SmallSystem class stores the portion of the finite element matrices and vectors that are currently being computed. The results of the flux_value() calculation are stored here.

Example of `flux_value()`

Here is the flux_value() routine for the thermal conductivity property, from SRC/engine/properties/heatconductivity/heatconductivity.C in the OOF2 source code. (Note that it could have been omitted, since in this case the heat flux can be computed from flux_matrix alone.)

void HeatConductivity::flux_value(
            const FEMesh *mesh,
            const Element *element,
            const Flux *flux,
            const MasterPosition &position,
            double time,
            void* data,
            SmallSystem *fluxdata)
  const
{
  DoubleVec fieldGradient = temperature->gradient(mesh, element, pt); 
  const SymmMatrix3 cond(conductivitytensor(mesh, element, pt)); 
  linsys->fluxVector() -= cond*fieldGradient; 
}

	`temperature` is a `ScalarField`*; stored in the `HeatConductivity` class. It was initialized by calling Field::getField() in the `HeatConductivity` constructor. `Field::gradient()` returns an object containing all the components of a `Field`.
	`conductivitytensor()` is a virtual function. The various isotropic and anisotropic subclasses of `HeatConductivity` define it differently.
	`SmallSystem::fluxVector()` returns a reference to the `Flux` components stored in the `SmallSystem`. Like the offset described in `flux_offset()`, it's always stored as a vector, independent of the type of the `Flux`.

`void static_flux_value(...) const`

void static_flux_value(const FEMesh* mesh,
                       const Element* element,
                       const Flux* flux,
                       const MasterPosition& position,
                       double time,
                       void* data,
                       SmallSystem *linysys) const

static_flux_value() computes the part of the Flux that depends only on the Fields and not on their time derivatives. It only needs to be defined if the Flux depends on the time derivatives of the Fields and is nonlinear. Its arguments are the same as the arguments to flux_value().

In the base class, static_flux_value() simply calls flux_value().

^[77] We are using the convention that repeated indices, such as $\nu$ in (188) are implicitly summed: $x_{i}y_{i}\equiv\sum_{i}x_{i}y_{i}$

^[78] 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.

^[79]This data can't be stored directly in the Flux object because one Flux is shared among many FEMeshes. The SmallSystem object is local to this computation.

^[80] The values of the flux offset are stored in the SmallSystem as a vector, for all types of Flux. This would be completely invisible to the user, if the name of the SmallSystem method for accessing the offset weren't offsetVector().

^[81] It would look like

for(IndexP kl : *flux->components(ALL_INDICES)) {
  const SymTensorIndex* kayell = dynamic_cast<const SymTensorIndex*>(kl.fieldindex());
  if(kayell->diagonal()) {
    offset_el -= modulus(ij,kl)*expten[kl]*T0;
  }
  else {
    offset_el -= 2.0*modulus(ij,kl)*expten[kl]*T0;
  }
}


PhysicalProperty		EqnProperty

OOF2: The Manual

Name

Synopses

C++ Synopsis

Python Synopsis

Source Files

Overview

Constructor

Methods

void* begin_point(...) const

void end_point(...) const

void flux_matrix(...) const

Example of flux_matrix()

void flux_offset(...) const

Example of flux_offset

void flux_value(...) const

Example of flux_value()

void static_flux_value(...) const