Development of a library for constraint types in the dockSIM research code

J. Boungard, J. Wackerfuß

For each constraint, a constraint method requires the following three quantities: the constraints c, the constraint Jacobian G and the constraint Hessian_Hi of all constraints. This has three implications for the implementation in the framework of the finite element code.

First, this allows the re-use of existing functions for the computation of constraint quantities. Therefore, a new constraint method can be easily added to existing FE-codes where other constraint methods have already been implemented.

Second, a library concept for the constraints can be employed, similarly to the element and material libraries usually exist in FEM codes. This allows for an easy extension of existing code with respect to new (arbitrary) constraints.

Third, the calculation of the three constraint quantities can be performed in a similar manner to the global quantities in the finite element method, the residual vector and the tangential stiffness matrix. For each constraint, the specific constraint value _ci, the specific constraint Jacobian_Gi and the specific constraint Hessian_Hi can be computed independently of other constraints, similarly to the element residual vector and element stiffness matrix in the FEM. Afterwards, these quantities are assembled in the same manner as in the FEM. Here, the assembly is performed over the number of constraints instead of the number of elements. The assembly of the constraint vector and the constraint Jacobian is even more simple compared to the FEM because no summation has to be performed (no overlapping within each column and row). For the assembly of all three quantities, exiting functions of an FEM code can be used, simplifying the implementation.

In most applications, the constraint Jacobian as well as the Hessian are sparse. This allows for the use of efficient storage techniques like the compressed sparse row format (CSR). To be able to allocate the matrices in CSR format, its structure of zero and non-zero entries has to be determined. This can be done in a pre-processing step in which the connectivity of the dofs by constraints is assessed. For this purpose, two information with respect to each constraint i have to be assessed and stored in the data structure of a finite element software:

First, the list of dofs _DOFSi, that are involved in the constraint equation of the specific constraint i, is required. Based on this list, the Jacobian is allocated. The i-th row of the Jacobian corresponds to the i-th constraint. The number of non-zero entries per row is equal to the number of affected dofs of constraint i, i.e. the number of elements in _DOFSi. The column indices of the non-zero entries are equal to the incides of the affected dofs, i.e. the values of elements in _DOFSi.
Second, a list of dof pairs of the constraint _PAIRSi is required. In this list, every dof affected by the specific constraint j is paired with every dof of this constraint k, resulting in (j,k)-pairs. Based on this list, the Hessian is allocated. For this purpose, all pair lists of the constraint are combined and sorted. Afterwards, duplicates are eliminated. Each pair (j,k) in the resulting total pair list _PAIRSi corresponds to a non-zero value in the Hessian for the j-th row and the k-th column.

Based on the previous remarks, the specific constraint library in dockSIM was designed. In this library concept, templates of functions are given. For each constraint, the corresponding functions have to be programmed. The specific functions of the constraints are encapsulated so that the programmer of new constraint solely operates with the constraint library and its templates. No additional changes in other parts of the software are required. The constraint library uses the following four functions:

processConstraint_<constraintTypeName>: This function checks the input of the specific constraint given by the user of the finite element software. If the input is valid, all relevant parameters of the constraint are stored in the data structure. If the input is invalid, the user gets constraint specific error messages.
addConstraintDofsToLists_<constraintTypeName>: This function generates the dof list _DOFSi and the dof pair list _PAIRSi of the specific constraint.
getConstraintVectorAndJacobian_<constraintTypeName>: This function computes the constraint value _ci as wells as the constraint Jacobian_Gi of the specific constraint and performs the assembly in the global constraint vector c and global constraint Jacobian G.
getProductStrategySpecificVectorAndHessian_<constraintTypeName>: This function computes the effective constraint Hessian_Hi of the specific constraint. This functions also performs the assembly into the global effective constraint Hessian Ĥ. This function is not required for linear constraints. For some constraint methods, the effective constraint Hessian_Hi depends on the constraint vector c and the constraint Jacobian G. Therefore, it is not possible to compute all three constraint quantities at once.