HMatrixUtils.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) OpenGeoSys Community (opengeosys.org)
// SPDX-License-Identifier: BSD-3-Clause
 
#pragma once
 
#include "NumLib/Fem/ShapeMatrixPolicy.h"
 
namespace ProcessLib
{
template <typename ShapeFunction, int DisplacementDim>
class HMatrixPolicyType
{
private:
    template <int N>
    using VectorType =
        typename ShapeMatrixPolicyType<ShapeFunction,
                                       DisplacementDim>::template VectorType<N>;

    template <int N, int M>
    using MatrixType =
        typename ShapeMatrixPolicyType<ShapeFunction, DisplacementDim>::
            template MatrixType<N, M>;
 
    // Dimensions of specific H-matrix for n-points and displacement dimension.
    static int const _number_of_dof = ShapeFunction::NPOINTS * DisplacementDim;
 
public:
    using StiffnessMatrixType = MatrixType<_number_of_dof, _number_of_dof>;
    using NodalForceVectorType = VectorType<_number_of_dof>;
 
    using HMatrixType = MatrixType<DisplacementDim, _number_of_dof>;
 
    using ConstitutiveMatrixType = MatrixType<DisplacementDim, DisplacementDim>;
    using ForceVectorType = VectorType<DisplacementDim>;
};

template <int DisplacementDim,
          int NPOINTS,
          typename N_Type,
          typename HMatrixType>
void computeHMatrix(N_Type const& N, HMatrixType& H)
{
    static_assert(1 < DisplacementDim && DisplacementDim <= 3,
                  "LinearHMatrix::computeHMatrix: DisplacementDim must be in "
                  "range (1,3].");
 
    H.setZero();
 
    for (int j = 0; j < DisplacementDim; j++)
    {
        H.block(j, j * NPOINTS, 1, NPOINTS) = N;
    }
}
 
}  // namespace ProcessLib