ConstraintDirichletBoundaryConditionLocalAssembler.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 "MeshLib/Elements/Elements.h"
#include "MeshLib/Elements/MapBulkElementPoint.h"
#include "MeshLib/Elements/Utils.h"
#include "NumLib/DOF/DOFTableUtil.h"
#include "NumLib/Fem/InitShapeMatrices.h"
#include "NumLib/Fem/Integration/GenericIntegrationMethod.h"
#include "NumLib/Fem/ShapeMatrixPolicy.h"
#include "ParameterLib/Parameter.h"
#include "ProcessLib/Process.h"
 
namespace ProcessLib
{
struct IntegrationPointData final
{
    IntegrationPointData(double const& detJ,
                         double const& integral_measure,
                         double const& integration_weight,
                         MathLib::Point3d&& bulk_element_point_)
        : detJ_times_integralMeasure_times_weight(detJ * integral_measure *
                                                  integration_weight),
          bulk_element_point(bulk_element_point_)
    {
    }
 
    double const detJ_times_integralMeasure_times_weight;
    MathLib::Point3d bulk_element_point;
 
    EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
};
 
class ConstraintDirichletBoundaryConditionLocalAssemblerInterface
{
public:
    virtual ~ConstraintDirichletBoundaryConditionLocalAssemblerInterface() =
        default;
 
    virtual double integrate(
        std::vector<GlobalVector*> const& x, double const t,
        std::function<Eigen::Vector3d(
            std::size_t const, MathLib::Point3d const&, double const,
            std::vector<GlobalVector*> const&)> const& getFlux) = 0;
};
 
template <typename ShapeFunction, int GlobalDim>
class ConstraintDirichletBoundaryConditionLocalAssembler final
    : public ConstraintDirichletBoundaryConditionLocalAssemblerInterface
{
protected:
    using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;
    using NodalMatrixType = typename ShapeMatricesType::NodalMatrixType;
    using NodalVectorType = typename ShapeMatricesType::NodalVectorType;
 
public:
    ConstraintDirichletBoundaryConditionLocalAssembler(
        MeshLib::Element const& surface_element,
        std::size_t /* local_matrix_size */,
        NumLib::GenericIntegrationMethod const& integration_method,
        bool const is_axially_symmetric, MeshLib::Mesh const& bulk_mesh,
        std::vector<std::pair<std::size_t, unsigned>> bulk_ids)
        : _surface_element(surface_element),
          _integration_method(integration_method),
          _bulk_element_id(bulk_ids[_surface_element.getID()].first),
          _surface_element_normal(MeshLib::calculateNormalizedSurfaceNormal(
              _surface_element, *(bulk_mesh.getElements()[_bulk_element_id])))
    {
        auto const shape_matrices =
            NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
                                      GlobalDim, NumLib::ShapeMatrixType::N_J>(
                _surface_element, is_axially_symmetric, _integration_method);
 
        auto const bulk_face_id = bulk_ids[_surface_element.getID()].second;
        auto const& bulk_element = *bulk_mesh.getElement(_bulk_element_id);
 
        auto const n_integration_points =
            _integration_method.getNumberOfPoints();
        _ip_data.reserve(n_integration_points);
 
        for (unsigned ip = 0; ip < n_integration_points; ++ip)
        {
            auto const& wp = _integration_method.getWeightedPoint(ip);
            auto bulk_element_point =
                MeshLib::getBulkElementPoint(bulk_element, bulk_face_id, wp);
            _ip_data.emplace_back(shape_matrices[ip].detJ,
                                  shape_matrices[ip].integralMeasure,
                                  wp.getWeight(),
                                  std::move(bulk_element_point));
        }
    }

    double integrate(
        std::vector<GlobalVector*> const& x, double const t,
        std::function<Eigen::Vector3d(
            std::size_t const, MathLib::Point3d const&, double const,
            std::vector<GlobalVector*> const&)> const& getFlux) override
    {
        auto const n_integration_points =
            _integration_method.getNumberOfPoints();
        // integrated_value +=
        //   int_{\Gamma_e} \alpha \cdot flux \cdot normal \d \Gamma
        double integrated_value = 0;
        for (unsigned ip = 0; ip < n_integration_points; ip++)
        {
            auto const bulk_flux = getFlux(
                _bulk_element_id, _ip_data[ip].bulk_element_point, t, x);
 
            // TODO find solution for 2d case
            double const bulk_grad_times_normal(
                Eigen::Map<Eigen::RowVectorXd const>(bulk_flux.data(),
                                                     bulk_flux.size())
                    .dot(_surface_element_normal));
 
            integrated_value +=
                bulk_grad_times_normal *
                _ip_data[ip].detJ_times_integralMeasure_times_weight;
        }
        return integrated_value;
    }
 
private:
    MeshLib::Element const& _surface_element;
 
    std::vector<IntegrationPointData> _ip_data;
 
    NumLib::GenericIntegrationMethod const& _integration_method;
    std::size_t const _bulk_element_id;
    Eigen::Vector3d const _surface_element_normal;
};
 
}  // namespace ProcessLib