SurfaceFluxLocalAssembler.h

Go to the documentation of this file.

 
#pragma once
 
#include "MeshLib/Elements/Elements.h"
#include "MeshLib/Elements/FaceRule.h"
#include "MeshLib/Elements/MapBulkElementPoint.h"
#include "NumLib/DOF/DOFTableUtil.h"
#include "NumLib/Fem/InitShapeMatrices.h"
#include "NumLib/Fem/ShapeMatrixPolicy.h"
 
namespace ProcessLib
{
class SurfaceFluxLocalAssemblerInterface
{
public:
    virtual ~SurfaceFluxLocalAssemblerInterface() = default;
 
    virtual void integrate(
        std::size_t const element_id,
        std::vector<GlobalVector*> const& x,
        MeshLib::PropertyVector<double>& specific_flux,
        double const t,
        MeshLib::Mesh const& bulk_mesh,
        std::function<Eigen::Vector3d(std::size_t const,
                                      MathLib::Point3d const&, double const,
                                      std::vector<GlobalVector*> const&)>) = 0;
};
 
template <typename ShapeFunction, typename IntegrationMethod, int GlobalDim>
class SurfaceFluxLocalAssembler final
    : public SurfaceFluxLocalAssemblerInterface
{
protected:
    using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;
    using NodalMatrixType = typename ShapeMatricesType::NodalMatrixType;
    using NodalVectorType = typename ShapeMatricesType::NodalVectorType;
 
public:
    SurfaceFluxLocalAssembler(
        MeshLib::Element const& surface_element,
        std::size_t /*const local_matrix_size*/,
        bool const is_axially_symmetric,
        unsigned const integration_order,
        MeshLib::PropertyVector<std::size_t> const& bulk_element_ids,
        MeshLib::PropertyVector<std::size_t> const& bulk_face_ids)
        : _surface_element(surface_element),
          _integration_method(integration_order),
          _bulk_element_id(bulk_element_ids[surface_element.getID()]),
          _bulk_face_id(bulk_face_ids[surface_element.getID()])
    {
        auto const n_integration_points =
            _integration_method.getNumberOfPoints();
 
        auto const shape_matrices =
            NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
                                      GlobalDim, NumLib::ShapeMatrixType::N_J>(
                _surface_element, is_axially_symmetric, _integration_method);
        for (std::size_t ip = 0; ip < n_integration_points; ++ip)
        {
            _detJ_times_integralMeasure.push_back(
                shape_matrices[ip].detJ * shape_matrices[ip].integralMeasure);
        }
    }
 
    void integrate(std::size_t const element_id,
                   std::vector<GlobalVector*> const& x,
                   MeshLib::PropertyVector<double>& specific_flux,
                   double const t,
                   MeshLib::Mesh const& bulk_mesh,
                   std::function<Eigen::Vector3d(
                       std::size_t const, MathLib::Point3d const&, double const,
                       std::vector<GlobalVector*> const&)>
                       getFlux) override
    {
        auto get_surface_normal =
            [this, &bulk_mesh](
                MeshLib::Element const& surface_element) -> Eigen::Vector3d {
            Eigen::Vector3d surface_element_normal;
            if (surface_element.getGeomType() == MeshLib::MeshElemType::LINE)
            {
                auto const bulk_normal = MeshLib::FaceRule::getSurfaceNormal(
                    *bulk_mesh.getElements()[_bulk_element_id]);
                auto const l0 = Eigen::Map<Eigen::Vector3d const>(
                    _surface_element.getNode(0)->getCoords());
                auto const l1 = Eigen::Map<Eigen::Vector3d const>(
                    _surface_element.getNode(1)->getCoords());
                Eigen::Vector3d const line = l1 - l0;
                surface_element_normal = line.cross(bulk_normal);
            }
            else
            {
                surface_element_normal =
                    MeshLib::FaceRule::getSurfaceNormal(surface_element);
            }
            surface_element_normal.normalize();
            // At the moment (2016-09-28) the surface normal is not oriented
            // according to the right hand rule. Thus for an intuitive flux
            // output, i.e., inflow has positive sign, outflow has negative
            // sign, the normal must not be multiplied by -1.
            return surface_element_normal;
        };
        auto const surface_element_normal =
            get_surface_normal(_surface_element);
 
        double element_area = 0.0;
        std::size_t const n_integration_points =
            _integration_method.getNumberOfPoints();
        // specific_flux[id_of_element] +=
        //   int_{\Gamma_e} \alpha \cdot flux \cdot normal \d \Gamma /
        //   element_area
        for (std::size_t ip(0); ip < n_integration_points; ip++)
        {
            auto const& wp = _integration_method.getWeightedPoint(ip);
 
            auto const bulk_element_point = MeshLib::getBulkElementPoint(
                bulk_mesh, _bulk_element_id, _bulk_face_id, wp);
            auto const bulk_flux =
                getFlux(_bulk_element_id, bulk_element_point, t, x);
            for (int component_id(0);
                 component_id < specific_flux.getNumberOfGlobalComponents();
                 ++component_id)
            {
                // 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));
 
                specific_flux.getComponent(element_id, component_id) +=
                    bulk_grad_times_normal * _detJ_times_integralMeasure[ip] *
                    wp.getWeight();
            }
            element_area += _detJ_times_integralMeasure[ip] * wp.getWeight();
        }
        for (int component_id(0);
             component_id < specific_flux.getNumberOfGlobalComponents();
             ++component_id)
        {
            specific_flux.getComponent(element_id, component_id) /=
                element_area;
        }
    }
 
private:
    MeshLib::Element const& _surface_element;
 
    std::vector<double> _detJ_times_integralMeasure;
 
    IntegrationMethod const _integration_method;
    std::size_t const _bulk_element_id;
    std::size_t const _bulk_face_id;
};
 
}  // namespace ProcessLib