ProcessLib::ForwardDifferencesJacobianAssembler Class Reference

Detailed Description

Assembles the Jacobian matrix using forward differences.

Definition at line 11 of file ForwardDifferencesJacobianAssembler.h.

#include <ForwardDifferencesJacobianAssembler.h>

Inheritance diagram for ProcessLib::ForwardDifferencesJacobianAssembler:

Collaboration diagram for ProcessLib::ForwardDifferencesJacobianAssembler:

Public Member Functions
	ForwardDifferencesJacobianAssembler (std::vector< double > &&absolute_epsilons)
void	assembleWithJacobian (std::size_t const mesh_item_id, LocalAssemblerInterface &local_assembler, double const t, double const dt, std::vector< double > const &local_x_data, std::vector< double > const &local_x_prev_data, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data) override
std::unique_ptr< AbstractJacobianAssembler >	copy () const override
Public Member Functions inherited from ProcessLib::NumericalJacobianAssembler
	NumericalJacobianAssembler (std::vector< double > &&absolute_epsilons)
bool	needsPicardAssembly () const override
void	checkPerturbationSize (int const max_non_deformation_dofs_per_node) const override
void	setNonDeformationComponentIDs (std::vector< int > const &non_deformation_component_ids) override
void	setNonDeformationComponentIDsNoSizeCheck (std::vector< int > const &non_deformation_component_ids) override
Public Member Functions inherited from ProcessLib::AbstractJacobianAssembler
virtual void	assembleWithJacobianForStaggeredScheme (LocalAssemblerInterface &, double const, double const, Eigen::VectorXd const &, Eigen::VectorXd const &, int const, std::vector< double > &, std::vector< double > &)
virtual	~AbstractJacobianAssembler ()=default
virtual void	preIteration (const unsigned)

Private Attributes
std::vector< double >	_local_M_data
std::vector< double >	_local_K_data
std::vector< double >	_local_b_data
std::vector< double >	_local_x_perturbed_data

Additional Inherited Members
Protected Member Functions inherited from ProcessLib::NumericalJacobianAssembler
auto	getVariableComponentEpsilonsView () const
Protected Attributes inherited from ProcessLib::NumericalJacobianAssembler
std::vector< double > const	absolute_epsilons_
std::vector< int >	non_deformation_component_ids_

Constructor & Destructor Documentation

ForwardDifferencesJacobianAssembler()

ProcessLib::ForwardDifferencesJacobianAssembler::ForwardDifferencesJacobianAssembler ( std::vector< double > && absolute_epsilons )

explicit

Constructs a new instance.

Parameters

absolute_epsilons

perturbations of the variable components used for evaluating the finite differences.

Note

The size of absolute_epsilons defines the "number of components" of the local solution vector (This is not the number of elements of the vector!). Therefore the size of the local solution vector must be divisible by the size of absolute_epsilons. This is the only consistency check performed. It is not checked whether said "number of components" is sensible. E.g., one could pass one epsilon per node, which would be valid but would not make sense at all.

Definition at line 12 of file ForwardDifferencesJacobianAssembler.cpp.

    : NumericalJacobianAssembler(std::move(absolute_epsilons))
{
}

References ProcessLib::NumericalJacobianAssembler::NumericalJacobianAssembler().

Member Function Documentation

assembleWithJacobian()

void ProcessLib::ForwardDifferencesJacobianAssembler::assembleWithJacobian ( std::size_t const mesh_item_id, LocalAssemblerInterface & local_assembler, double const t, double const dt, std::vector< double > const & local_x_data, std::vector< double > const & local_x_prev_data, std::vector< double > & local_b_data, std::vector< double > & local_Jac_data )

overridevirtual

Assembles the Jacobian, the matrices \(M\) and \(K\), and the vector \(b\). For the assembly the assemble() method of the given local_assembler is called several times and the Jacobian is built from finite differences. The number of calls of the assemble() method is \(2N+1\) if \(N\) is the size of local_x.

Attention

It is assumed that the local vectors and matrices are ordered by component.

Implements ProcessLib::AbstractJacobianAssembler.

Definition at line 18 of file ForwardDifferencesJacobianAssembler.cpp.

{
    std::vector<double> local_M_data(local_Jac_data.size());
    std::vector<double> local_K_data(local_Jac_data.size());
 
    auto const num_r_c =
        static_cast<Eigen::MatrixXd::Index>(local_x_data.size());
 
    auto const x = MathLib::toVector<Eigen::VectorXd>(local_x_data, num_r_c);
    auto const x_prev =
        MathLib::toVector<Eigen::VectorXd>(local_x_prev_data, num_r_c);
 
    Eigen::VectorXd const local_xdot = (x - x_prev) / dt;
 
    auto local_Jac =
        MathLib::createZeroedMatrix(local_Jac_data, num_r_c, num_r_c);
 
    assert(this->non_deformation_component_ids_.size() > 0);
 
    auto const num_dofs_per_component =
        local_x_data.size() / this->non_deformation_component_ids_.size();
 
    // Assemble with unperturbed local x to get M0, K0, and b0 used in the
    // finite differences below.
    local_assembler.assemble(t, dt, local_x_data, local_x_prev_data,
                             local_M_data, local_K_data, local_b_data);
 
    auto const nved = local_assembler.getNumberOfVectorElementsForDeformation();
 
    // Residual  res := M xdot + K x - b
    // Computing Jac := dres/dx
    //                = M dxdot/dx + dM/dx xdot + K dx/dx + dK/dx x - db/dx
    //                  with dxdot/dx = 1/dt and dx/dx = 1
    //                  (Note: dM/dx and dK/dx actually have the second and
    //                  third index transposed.)
    // The loop computes the dM/dx, dK/dx and db/dx terms, the rest is computed
    // afterwards. The loop skips the entries corresponding to the deformation
    // part of the solution vector if a vector segment size is given by nved.
    // This is to avoid recomputing the analytic block of the deformation
    // process.
    auto const num_purterbated_colums = num_r_c - nved;
    auto const perturbations = this->getVariableComponentEpsilonsView();
    for (Eigen::MatrixXd::Index i = 0; i < num_purterbated_colums; ++i)
    {
        // assume that local_x_data is ordered by component.
        auto const component = i / num_dofs_per_component;
        auto const eps = perturbations[component];
 
        // Assemble with perturbed local x.
        _local_x_perturbed_data = local_x_data;
        _local_x_perturbed_data[i] = local_x_data[i] + eps;
 
        local_assembler.assemble(t, dt, _local_x_perturbed_data,
                                 local_x_prev_data, _local_M_data,
                                 _local_K_data, _local_b_data);
 
        if (!local_M_data.empty() && !_local_M_data.empty())
        {
            auto const local_M_0 =
                MathLib::toMatrix(local_M_data, num_r_c, num_r_c);
            auto const local_M_p =
                MathLib::toMatrix(_local_M_data, num_r_c, num_r_c);
            local_Jac.col(i).noalias() +=
                // dM/dxi * x_dot
                (local_M_p - local_M_0) * local_xdot / eps;
            _local_M_data.clear();
        }
        if (!local_K_data.empty() && !_local_K_data.empty())
        {
            auto const local_K_0 =
                MathLib::toMatrix(local_K_data, num_r_c, num_r_c);
            auto const local_K_p =
                MathLib::toMatrix(_local_K_data, num_r_c, num_r_c);
 
            local_Jac.col(i).noalias() +=
                // dK/dxi * x
                (local_K_p - local_K_0) * x / eps;
            _local_K_data.clear();
        }
        if (!local_b_data.empty() && !_local_b_data.empty())
        {
            auto const local_b_0 =
                MathLib::toVector<Eigen::VectorXd>(local_b_data, num_r_c);
            auto const local_b_p =
                MathLib::toVector<Eigen::VectorXd>(_local_b_data, num_r_c);
            local_Jac.col(i).noalias() -= (local_b_p - local_b_0) / eps;
            _local_b_data.clear();
        }
    }
 
    // Assemble with unperturbed local x, i.e. compute M dxdot/dx + K dx/dx =
    // M/dt + K
    // Compute remaining terms of the Jacobian.
    if (!local_M_data.empty())
    {
        auto local_M = MathLib::toMatrix(local_M_data, num_r_c, num_r_c);
        local_Jac.noalias() += local_M / dt;
    }
    if (!local_K_data.empty())
    {
        auto local_K = MathLib::toMatrix(local_K_data, num_r_c, num_r_c);
        local_Jac.noalias() += local_K;
    }
 
    // Move the M and K contributions to the residuum for evaluation of nodal
    // forces, flow rates, and the like. Cleaning up the M's and K's storage so
    // it is not accounted for twice.
    auto b = [&]()
    {
        if (!local_b_data.empty())
        {
            return MathLib::toVector<Eigen::VectorXd>(local_b_data, num_r_c);
        }
        return MathLib::createZeroedVector<Eigen::VectorXd>(local_b_data,
                                                            num_r_c);
    }();
 
    if (!local_M_data.empty())
    {
        auto M = MathLib::toMatrix(local_M_data, num_r_c, num_r_c);
        b -= M * (x - x_prev) / dt;
        local_M_data.clear();
    }
    if (!local_K_data.empty())
    {
        auto K = MathLib::toMatrix(local_K_data, num_r_c, num_r_c);
 
        // Note: The deformation segment of \c b is already computed as
        // int{B^T sigma}dA, which is identical to K_uu * u. Therefore the
        // corresponding K block is set to zero.
        if (nved != 0)
        {
            auto const dm_start_index = num_purterbated_colums;
            auto const dm_size = nved;
            K.block(dm_start_index, dm_start_index, dm_size, dm_size).setZero();
        }
 
        b -= K * x;
        local_K_data.clear();
    }
}

References _local_b_data, _local_K_data, _local_M_data, _local_x_perturbed_data, ProcessLib::LocalAssemblerInterface::assemble(), MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), ProcessLib::LocalAssemblerInterface::getNumberOfVectorElementsForDeformation(), ProcessLib::NumericalJacobianAssembler::getVariableComponentEpsilonsView(), ProcessLib::NumericalJacobianAssembler::non_deformation_component_ids_, MathLib::toMatrix(), and MathLib::toVector().

copy()

std::unique_ptr< AbstractJacobianAssembler > ProcessLib::ForwardDifferencesJacobianAssembler::copy ( ) const

overridevirtual

Implements ProcessLib::AbstractJacobianAssembler.

Definition at line 166 of file ForwardDifferencesJacobianAssembler.cpp.

{
    return std::make_unique<ForwardDifferencesJacobianAssembler>(*this);
}

Member Data Documentation

_local_b_data

std::vector<double> ProcessLib::ForwardDifferencesJacobianAssembler::_local_b_data

private

Definition at line 55 of file ForwardDifferencesJacobianAssembler.h.

Referenced by assembleWithJacobian().

_local_K_data

std::vector<double> ProcessLib::ForwardDifferencesJacobianAssembler::_local_K_data

private

Definition at line 54 of file ForwardDifferencesJacobianAssembler.h.

Referenced by assembleWithJacobian().

_local_M_data

std::vector<double> ProcessLib::ForwardDifferencesJacobianAssembler::_local_M_data

private

Definition at line 53 of file ForwardDifferencesJacobianAssembler.h.

Referenced by assembleWithJacobian().

_local_x_perturbed_data

std::vector<double> ProcessLib::ForwardDifferencesJacobianAssembler::_local_x_perturbed_data

private

Definition at line 56 of file ForwardDifferencesJacobianAssembler.h.

Referenced by assembleWithJacobian().

---

The documentation for this class was generated from the following files:

• ForwardDifferencesJacobianAssembler.h
• ForwardDifferencesJacobianAssembler.cpp