plane_approximation_utility.h

Go to the documentation of this file.

//    |  /           |
//    ' /   __| _` | __|  _ \   __|
//    . \  |   (   | |   (   |\__ `
//   _|\_\_|  \__,_|\__|\___/ ____/
//                   Multi-Physics
//
//  License:         BSD License
//                   Kratos default license: kratos/license.txt
//
//  Main authors:    Ruben Zorrilla
//
 
#if !defined(KRATOS_PLANE_APPROXIMATION_UTILITY )
#define  KRATOS_PLANE_APPROXIMATION_UTILITY
 
/* System includes */
 
/* External includes */
 
/* Project includes */
#include "utilities/math_utils.h"
 
namespace Kratos
{
 
 
 
 
 
 
template <unsigned int TDim>
class PlaneApproximationUtility
{
public:
 
 
    KRATOS_CLASS_POINTER_DEFINITION( PlaneApproximationUtility );
 
 
 
    PlaneApproximationUtility(){}
 
    virtual ~PlaneApproximationUtility(){}
 
 
    static void ComputePlaneApproximation(
        const std::vector< array_1d<double,3> > &rPointsCoords,
        array_1d<double,3> &rPlaneBasePointCoords,
        array_1d<double,3> &rPlaneNormal)
    {
        ComputeBasePoint(rPointsCoords,rPlaneBasePointCoords);
        ComputePlaneNormal(rPointsCoords, rPlaneBasePointCoords, rPlaneNormal);
    }
 
 
 
 
 
private:
 
 
 
 
 
    static void ComputeBasePoint(
        const std::vector< array_1d<double,3> > &rPointsCoords,
        array_1d<double,3> &rBasePointCoords)
    {
        const unsigned int n_points = rPointsCoords.size();
        KRATOS_ERROR_IF(n_points == 0) << "No base point can be computed. Points container is empty." << std::endl;
 
        noalias(rBasePointCoords) = ZeroVector(3);
        for (unsigned int j = 0; j < n_points; ++j){
            rBasePointCoords += rPointsCoords[j];
        }
        rBasePointCoords /= n_points;
    }
 
    static void SetMatrixA(
        const std::vector< array_1d< double,3 > > &rPointsCoords,
        const array_1d<double,3> &rPlaneBasePointCoords,
        BoundedMatrix<double,TDim,TDim> &rA)
    {
        noalias(rA) = ZeroMatrix(TDim,TDim);
        const unsigned int n_points = rPointsCoords.size();
 
        for (unsigned int i = 0; i < TDim; ++i){
            const double base_i = rPlaneBasePointCoords(i);
            for (unsigned int j = 0; j < TDim; ++j){
                const double base_j = rPlaneBasePointCoords(j);
                for (unsigned int m = 0; m < n_points; ++m){
                    const double pt_m_i = rPointsCoords[m](i);
                    const double pt_m_j = rPointsCoords[m](j);
                    rA(i,j) += (base_i - pt_m_i)*(base_j - pt_m_j);
                }
            }
        }
    }
 
    static void ComputePlaneNormal(
        const std::vector< array_1d<double,3> > &rPointsCoords,
        const array_1d<double,3> &rPlaneBasePointCoords,
        array_1d<double,3> &rPlaneNormal)
    {
        // Solve the A matrix eigenvalue problem
        BoundedMatrix<double, TDim, TDim> a_mat, eigenval_mat, eigenvector_mat;
        SetMatrixA(rPointsCoords, rPlaneBasePointCoords, a_mat);
        bool converged = MathUtils<double>::GaussSeidelEigenSystem(a_mat, eigenvector_mat, eigenval_mat);
        KRATOS_ERROR_IF(!converged) << "Plane normal can't be computed. Eigenvalue problem did not converge." << std::endl;
 
        // Find the minimum eigenvalue
        double min_eigval = std::numeric_limits<double>::max();
        unsigned int min_eigval_id = 0;
        for (unsigned int i = 0; i < TDim; ++i){
            if (eigenval_mat(i,i) < min_eigval){
                min_eigval_id = i;
                min_eigval = eigenval_mat(i,i);
            }
        }
 
        // Set as plane normal the eigenvector associated to the minimum eigenvalue
        noalias(rPlaneNormal) = ZeroVector(3);
        for (unsigned int i = 0; i < TDim; ++i){
            rPlaneNormal(i) = eigenvector_mat(i, min_eigval_id);
        }
    }
 
 
 
 
 
}; /* Class PlaneApproximationUtility */
 
 
 
 
}  /* namespace Kratos.*/
 
#endif /* KRATOS_PLANE_APPROXIMATION_UTILITY  defined */