ug4
UserData

An IElemDisc object can assemble the local stiffness/mass matrix and/or the local defect/rhs. In order to accomplish this task user defined data is often needed for several integration points on the element. On this page, we describe the idea of the UserData import, that is used in ug.

While this page is concentrated on the technical usage, a more detailed and mathematical introduction to the import/export strategy can be found on Import/Export Concept.

Data Import

Each IElemDisc may have several DataImports, describing some user defined functions or data given from outside of the IElemDisc. To give an idea, this may be a Diffusion Tensor, a Velocity field, some reaction constants, but also data computed by other IElemDisc, such as their local solution values at those points or some computed velocity field etc. This data must be known at several points on each element. We call this points IntegrationPoints, as they are usually linked to some quadrature procedure. Please note, that the position of the integration points is uniquely and only determined by the requesting IElemDisc. Other (data exporting) IElemDiscs or the user defined Data must not need any information on how and why these integration points have been chosen. This is important to ensure, that IElemDiscs remain independent from each other.

The IElemDisc tells the DataImport, at which integration points it wants to evaluate the data. This is done in two fashions:

  • local integration points: A vector of all ips expressed in reference element coordinates. This is usually needed for data exported by other IElemDiscs
  • global integration points: A vector if all ips expressed in global (i.e. physical world) coordinates. This is usually needed for user-defined data

Please note, that only a pointer to the vector of coordinates together with the number of integration points is stored when requesting user data. Some reason for that:

  • Storing positions in a vector-like structure gives good cache hit rates
  • Computing the data for all ips with only one virtual call
  • The identity of two requested ip series is compared using the pointer to the integration points. Iff two pointers are equal, the computation of the second request is skipped, since not needed. This comparison and scheduling of evaluation is done before the loop over all elements starts.

Lets show the usage with the example of a Diffusion Tensor Import. The Import is realized by class member

DataImport<MathMatrix<dim,dim>, dim, algebra_type> m_Diff;
dim
static const int dim

This indicates, that the computed/imported data has the C++-type MathMatrix<dim,dim>, the dimension of the global coordinates is dim. This Import must be registered, in order to include it to the evaluation process, by calling (e.g. in the constructor):

// register imports
register_import(m_Diff);

The request of a local integration point series must be done in prep_elem_loop(const ReferenceObjectID roid, const int si) , like e.g.

template< template <class TElem, int TWorldDim> class TFVGeom,
typename TDomain,
typename TAlgebra>
template<typename TElem>
bool FVConvectionDiffusion<TFVGeom, TDomain, TAlgebra>::
prep_elem_loop(const ReferenceObjectID roid, const int si)
{
// get reference element type
typedef typename reference_element_traits<TElem>::reference_element_type
ref_elem_type;
// get Finite-Volume Geometry
TFVGeom<TElem, dim>& geo = FVGeometryProvider::get_geom<TFVGeom, TElem, dim>();
// set local positions
m_Diff.template set_local_ips<ref_elem_type::dim>(geo.scvf_local_ips(),
geo.num_scvf_ips());
}
ug::ReferenceObjectID
ReferenceObjectID
these ids are used to identify the shape of a geometric object.
Definition: grid_base_objects.h:74

The request of the global integration points takes place in prep_elem(), after the global ips have been updated for the element:

template< template <class TElem, int TWorldDim> class TFVGeom,
typename TDomain,
typename TAlgebra>
template<typename TElem >
bool FVConvectionDiffusion<TFVGeom, TDomain, TAlgebra>::
prep_elem(TElem* elem, const LocalVector& u,
const LocalIndices& glob_ind)
{
// Load corners of this element
for(size_t i = 0; i < m_vCornerCoords.size(); ++i)
{
Vertex* vert = elem->vertex(i);
m_vCornerCoords[i] = m_aaPos[vert];
}
// Update Geometry for this element
TFVGeom<TElem, dim>& geo = FVGeometryProvider::get_geom<TFVGeom, TElem,dim>();
try{
geo.update(elem, m_pDomain->subset_handler(), vCornerCoords))
}
UG_CATCH_THROW("FVConvectionDiffusion::prep_elem:"
" Cannot update Finite Volume Geometry.");
// set global positions for rhs
m_Diff.set_global_ips(geo.scvf_global_ips(), geo.num_scvf_ips());
}
UG_CATCH_THROW
#define UG_CATCH_THROW(msg)
Definition: error.h:64

Then, in each assemble_...() function the data can be accessed at the integration points like this

size_t ip = 0;
if(!m_Diff.zero_data())
{
// loop Sub Control Volume Faces (SCVF)
for(size_t i = 0; i < geo.num_scvf(); ++i)
{
// get current SCVF
const typename TFVGeom<TElem, dim>::SCVF& scvf = geo.scvf(i);
// loop integration point of SCVF
for(size_t i = 0; i < scvf.num_ip(); ++i, ++ip)
{
// loop shape functions
for(size_t j = 0; j < scvf.num_sh(); ++j)
{
// Compute Diffusion Tensor times Gradient
MathMatrix<dim, dim> ipDiffTensor = m_Diff[ip];
}
}
}
}

UserData

The base class for Data, that can be plugged into the Data Import is given by the class UserData. In ug, we distinguish four kinds of UserData:

  • default data (usually zero or identity data)
  • constant data
  • data depending on \( (\vec{x}, t) \), i.e. on the global space coordinate and the time.
  • data produced by other IElemDiscs.

Default data is used, iff no data has been set to the Data Import. The the IElemDisc may decide to use a predefined default value. This is usually zero, sometimes identity.

The base class for the last three data, is the class

template <typename TData>
class UserData

This class can provide the data for the templated type. It knows:

  • The number of integration points
  • The Data at the integration points

Note, that we have to serve a variety of possible setting. This ranges between to extremes:

  • Data is used by several IElemDiscs at exactly the same points. Thus, we want to compute the data only once.
  • Data is needed at different points. Thus, several data sets must be produced using the same evaluation function, but at different points.

In order to implement this, every DataImport is connected to exactly one UserData. Then, it requests the needed ips using register_local_ip_series( ... ) and gets a Series ID. This id is like a ticket for the Import. Each time it requests the data at its points, the Import must use the ID to indicate which series it has requested. Internally, requests to identical ip points, will produce identical series IDs and the data at those points will only be computed once, while several Imports have access to the data.

Data Export

Each IElemDisc can provide as many DataExports as it wants. In general, an Export can be regarded as the offer to compute data at integration points (and its derivative w.r.t. the local unknowns in case of Jacobian). Please note: if no request is made to use this offer, no computation or call to an export takes place, resulting in a zero overhead.

In order to implement the a DataExport, a class defines a member of this type.

DataExport<MathVector<dim>, dim, algebra_type> m_DarcyVelExport;

Now, a evaluation function must be written in the same class, with the signature:

template <typename TElem>
bool
IElemDisc<TAlgebra>::
compute_darcy_export(const LocalVector& u, bool compDeriv);

A DataExport uses the local integration points to compute the data.

On creation the IElemDisc must register all element types it can evalute the data for, by

template <typename TElem>
void register_all_assemble_functions(int id)
{
typedef DensityDrivenFlow T;
m_DarcyVelExport.reg_export_fct(id, this, &T::template compute_darcy_export<TElem>);
}

EvaluationOrder

The order of evaluation is as follows.

  1. Setup (before loop over all elements)
    1. Get Imports from every IElemDisc.
    2. Sort the connected data by its behaviour: constant, user-defined and data-export data. Each type is scheduled for computation at a certain point.
    3. Remember imports that are connected to a non-zero-derivative data. This imports are scheduled for the evaluation of the linearized defect.
    4. TODO: Skip Exports, that are already scheduled.
    5. set_roid(...): Set geometric object type for each IElemDisc to indicate which evaluation function must be used
    6. prep_elem_loop(...): Call preparation of element loop for each IElemDisc. Local positions are set for DataImports here.
    7. For all DataImports: Set the geometric object type to indicate which evaluation function must be used
    8. For all DataExports: Set the geometric object type to indicate which evaluation function must be used
    9. Compute Constant User Data , since they only depend on number of integration points
  2. Loop all elements. For each element do:
    1. prep_elem(...): Prepare element for each IElemDisc. Here, the new global ips are computed and set for each IDataImport
    2. Compute all user-defined data
    3. Compute all exports
    4. For Jacobian: Compute linearized defect w.r.t. to Imports
    5. assemble_....(...): Call assemble function of each IElemDisc
    6. For Jacobian: Assemble inter-IElemDisc coupling by multiplying linearized defect with Data derivative for each import.
  3. Post process (after loop over all elements)
    1. fsh_elem_loop(): Finish element loop for each IElemDisc

The main point is, that first the ip positions are set and then the UserData is computed before the local assemble_... functions are called for each IElemDisc. Therefore, the computed user data can be used in these functions.