Constraints and time integration schemes

This is merely intended as a reminder that we have to tackle that, and maybe as a space for @joergfunger to articulate his ideas.

Currently, there are multiple people using Boost Odeint for time stepping, and pass an Eigen::VectorXd as the independent DoF values. Then you have to use the DofType, Dof numbering and constraint information to create a full DofVector, and then merge that into the nodes.

a) It would probably beneficial if not everyone did that themselves. b) When/where/how to apply the constraints?

For b), the idea is to provide a ConstraintApplicator (hopefully with a less stupid name), that has for example

    SparseMatrix ApplyBCs(SparseMatrix, Constraints::Constraints);

where you could put in an unconstrained matrix and get back the constrained matrix. Whether that happens via the C matrix, or Lagrange parameters, or penalty could then be decided via switching to another ConstraintApplicator. The ConstraintApplicator could either be an abstract base class or a concept.

I also prefer to completely remove constraints as much as possible from all classes such as time integration schemes or the build gradient/hessian stuff. Currently, I feel that we have not found yet an intuitive class design and what information is passed between the classes.

Regarding the constraints, in what situation do I need the constraints? In my opinion, the only occurance is the DoStep-routine of a time integration scheme. So I would pass the constraints as a parameter to these functions. This has the additional side effect of easily being able to change constraints (not only the rhs) between time steps.
Where should we use only the independent dofs and where the full dof vector? This question is very important in order to define where we should use a full dof vector and where directly an Eigen::Vector. The constraint scheme is then the one that has to perform the conversion between these two spaces. In my opinion, the time integration scheme (boost::odeint, ...) gets as input a parameter vector and the gradients (hessians), and perform one step. At this point, we do not need to have any information regarding the original dof type of each dof, and we also only require to use the independent dofs. Even in the case of different scalings between different dof types, this should be handled in a preconditioning step. As a consequence, all these vectors and matrices (vector of input dofs, gradient and hessian) should directly be an Eigen::Matrix.
Naming and what should not be covered ConstraintApplicator is imo not fully correct. The class performs a linear mapping between all the nodal values (function space of dimension n = number of dofs) and the function space of dimension j with j<n that we solve for. The standard case is that we reduce the dimension by fixing some dependent dofs, but in general, there is no reason why the class should not handle any other linear mapping. A straightforward example is the POD approach (that will have to be handled in the next project starting in August 2018), where linear combinations of global basis functions computed from a SVD of some training samples are used to represent the nodal values. This would mean that any other approaches to enforce constraint equations such as the penalty approach or the Lagrange multiplier approach should not be handled by this approach. This is because in the penalty approach there is no reduction of the function space, we just add additional "regularization terms to our global potential that enforce the constraint equations. As such, this should conceptually be implemented by just adding some contributions (via the cell interface) to the global matrices. The Langrange approach introduces additional dofs (the Lagrange multiplier), which actually increases the dimension of the problem. These Lagrange multipliers are part of the solution vector. As a consequence, they should also be implemented differently.
Where should we apply constraints then? As a consequence, I suggest that we apply constrains to the (not yet existing) mesh nodes, and at this point decide, whether we want to do Penalty, (Nonlinear) Lagrange or a linear mapping d = C d_j + b(t). I would currently only implement the last one, Peter will certainly have to implement the Langrange approach for the contact models. In the current implementation (without this mesh), we
What is the functionality that a LinearConstraintApplicator class should have?
- The class should perform the mapping between the reduced space and the full space, with the time t as additional parameter.
- The class should tranform gradients/hessians from the full space to the reduced space, either inplace or by returning the tranformed vector/matrix.

@potto1 and I have discussed, where and how the constraint application should be modified. Mygeneral idea is to have a class LinearConstraintApplicator (naming is not fixed yet) that does all the conversions.

class LinearConstraintApplicator
{
public:
    LinearConstraintApplicator(const std::vector<DofType>& dofTypes, const DofContainer<int>& mNumTotalDofs,
                               const Constraints& constraints)
        : mDofTypes(dofTypes)
        , mDofContainer<int> mNumTotalDofs
        , mConstraints(constraints);

    //! @brief transforms a hessian matrix from the full system by computing C^T K C to the reduced system
    //! @param Mat matrix in the original space
    //! @return hessian in the reduced space
    Eigen::SparseMatrix<double> HessianToReducedBasis(const Eigen::SparseMatrix<double>& m); 

    //! @brief transforms a gradient from the full system by computing C^T v to the reduced system
    //! @param v vector in the original space
    //! @return hessian in the reduced space
    Eigen::VectorXd GradientToReducedBasis(Eigen::VectorXd& v);

    //! @brief transforms a reduced solution vector back to the full system
    //! @param uj vector of independent dofs
    //! @param t time to compute the dependent dofs based on the constraint equations
    //! @return C*uj+b(t)
    Eigen::VectorXd ToFull(const Eigen::VectorXd& uj, double t);

    //! @brief transforms a full solution vector back to a dofVector and replaces the corresponding entries
    //! @param u vector of all active dofs
    //! @param dovV dof vector of all dofs
    Eigen::VectorXd ToDofVector(const Eigen::VectorXd& uj, DofVector& dofV);

    //! @brief transforms a delta of the reduced solution vector back to the full system (assuming t=const)
    //! @param uj vector of independent dofs
    //! @return C*duj
    Eigen::VectorXd ToDeltaFull(const Eigen::VectorXd& uj);

    //! @brief reduces the dof solution vector to the vector of independent dofs
    //! @param dofV actual state of the system
    //! @param t time in order to verify consistent initial conditions (constraint equations are fulfilled)
    //! @return eigenvector with only the independent dofs
    Eigen::VectorXd ToReducedBasis(const DofVector& dofV, t);

private:
    // these mdofTypes and the mNumTotalDofs are only tmp, there should be a reference to the function space that is
    // constraint by the LinearConstraintapplicator
    //! @brief doftypes used in the solution procedure
    std::vector<DofType> mDofTypes;
    //! @brief number of total dofs for each dof type in the solution procedure
    DofContainer<int> mNumTotalDofs;
    //! @brief the linear transformation matrix to map the reduced space to the full space d_tot = C d_red + b(t)
    Eigen::Sparse<double> mCmat;
    //! @brief constraints
    Constraints mConstraints;
};

This class would be passed as parameter to the DoStep function of the time integration scheme and would perform all the conversions. mDofs and mConstraints would be removed from the Quasistatic solver (or Newmark). This class is essentially to map the DofVectors/Matrices to the reduced variables (of type eigen) including constraints.

IMO, this is the right way to go. The function to integrate, its derivatives and the constraints are passed to compute the next step. The alternative to compute the constraints in the function and its derivatives is less efficient from the implementation point of view. The Lagrange multipliers are additional dofs and will then be treated as other dofs, too.

nutofem / nuto

Constraints and time integration schemes #253