ElementCoordinatesMappingLocal.cpp

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) OpenGeoSys Community (opengeosys.org)
// SPDX-License-Identifier: BSD-3-Clause
 
#include "ElementCoordinatesMappingLocal.h"
 
#include <cassert>
#include <limits>
 
#include "GeoLib/AnalyticalGeometry.h"
#include "MathLib/MathTools.h"
#include "MathLib/Point3d.h"
#include "MeshLib/Elements/Element.h"
#include "MeshLib/Node.h"
 
namespace detail
{
void rotateToLocal(const MeshLib::RotationMatrix& matR2local,
                   std::vector<MathLib::Point3d>& points)
{
    std::transform(points.begin(), points.end(), points.begin(),
                   [&matR2local](auto const& pnt) { return matR2local * pnt; });
}

MeshLib::RotationMatrix getRotationMatrixToGlobal(
    const unsigned element_dimension,
    const unsigned global_dim,
    const std::vector<MathLib::Point3d>& points)
{
    Eigen::Matrix3d matR;
    // compute R in x=R*x' where x are original coordinates and x' are local
    // coordinates
    if (element_dimension == 1)
    {
        Eigen::Vector3d const xx =
            (points[1].asEigenVector3d() - points[0].asEigenVector3d())
                .normalized();
        if (global_dim == 2)
        {
            matR = GeoLib::compute2DRotationMatrixToX(xx);
        }
        else
        {
            matR = GeoLib::compute3DRotationMatrixToX(xx);
        }
        matR.transposeInPlace();
    }
    else if (global_dim == 3 && element_dimension == 2)
    {
        // get plane normal
        auto const [plane_normal, d] = GeoLib::getNewellPlane(points);
        // compute a rotation matrix to XY
        matR = GeoLib::computeRotationMatrixToXY(plane_normal);
        // set a transposed matrix
        matR.transposeInPlace();
    }
    return matR;
}
}  // namespace detail
 
namespace MeshLib
{
ElementCoordinatesMappingLocal::ElementCoordinatesMappingLocal(
    const Element& e, const unsigned global_dim)
    : _global_dim(global_dim), _matR2global(Eigen::Matrix3d::Identity())
{
    unsigned const space_dim = std::max(e.space_dimension_, _global_dim);
    assert(e.getDimension() <= space_dim);
    _points.reserve(e.getNumberOfNodes());
    for (unsigned i = 0; i < e.getNumberOfNodes(); i++)
    {
        _points.emplace_back(*(e.getNode(i)));
    }
 
    auto const element_dim = e.getDimension();
 
    if (element_dim != space_dim)
    {
        _matR2global =
            detail::getRotationMatrixToGlobal(element_dim, space_dim, _points);
        detail::rotateToLocal(_matR2global.transpose(), _points);
    }
}
 
}  // namespace MeshLib