GeoData.h

Go to the documentation of this file.

#ifndef GEO_DATA_H_INCLUDED
#define GEO_DATA_H_INCLUDED
/*
 * Copyright 2011-2016 Universiteit Antwerpen
 *
 * Licensed under the EUPL, Version 1.1 or  as soon they will be approved by
 * the European Commission - subsequent versions of the EUPL (the "Licence");
 * You may not use this work except in compliance with the Licence.
 * You may obtain a copy of the Licence at: http://ec.europa.eu/idabc/eupl5
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the Licence is distributed on an "AS IS" basis,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the Licence for the specific language governing
 * permissions and limitations under the Licence.
 */
#include "AreaMoment.h"

#include "math/array3.h"

#include <array>
#include <iostream>
#include <tuple>

namespace SimPT_Sim {

struct GeoData
{
        GeoData()
                : m_area(0.0)
        {
        }

        GeoData(double area, const std::array<double, 3>& centroid, const AreaMoment& area_moment)
                : m_area(area), m_centroid(centroid), m_area_moment(area_moment)
        {
        }

        std::tuple<double, std::array<double, 3>, double> GetEllipseAxes() const
        {
                double const ea_xx   = m_area_moment.m_xx - (m_centroid[0] * m_centroid[0]) * m_area;
                double const ea_xy   = m_area_moment.m_xy + (m_centroid[0] * m_centroid[1]) * m_area;
                double const ea_yy   = m_area_moment.m_yy - (m_centroid[1] * m_centroid[1]) * m_area;
                double const rhs1    = 0.5 * (ea_xx + ea_yy);
                double const rhs2    = 0.5 * sqrt((ea_xx - ea_yy) * (ea_xx - ea_yy) + 4 * ea_xy * ea_xy);

                const double                 length { 4.0 * sqrt((rhs1 + rhs2) / m_area) };
                const double                 width  { 4.0 * sqrt((rhs1 - rhs2) / m_area) };
                const std::array<double, 3>  long_axis {{ -ea_xy, rhs1 + rhs2 - ea_xx, 0.0 }};

                return std::make_tuple(length, long_axis, width);
        }


        GeoData& operator+=(const GeoData& rhs)
        {
                m_area         += rhs.m_area;
                m_centroid     += rhs.m_centroid;
                m_area_moment  += rhs.m_area_moment;
                return *this;
        }

        GeoData& operator-=(const GeoData& rhs)
        {
                m_area         -= rhs.m_area;
                m_centroid     -= rhs.m_centroid;
                m_area_moment  -= rhs.m_area_moment;
                return *this;
        }

        double                  m_area;
        std::array<double, 3>   m_centroid;
        AreaMoment              m_area_moment;
};


inline GeoData ModifyForMove(const GeoData& old_geo,
                const std::array<double, 3>& nb1_p, const std::array<double, 3>& nb2_p,
                const std::array<double, 3>& now_p, const std::array<double, 3>& mov_p)
{
        const double del_area
                = 0.5 * ((mov_p[0] - now_p[0]) * (nb1_p[1] - nb2_p[1])
                       + (mov_p[1] - now_p[1]) * (nb2_p[0] - nb1_p[0]));
        const double new_area = old_geo.m_area - del_area;

        const double del_x
                = (nb1_p[0] + mov_p[0]) * (mov_p[0] * nb1_p[1] - nb1_p[0] * mov_p[1])
                + (mov_p[0] + nb2_p[0]) * (nb2_p[0] * mov_p[1] - mov_p[0] * nb2_p[1])
                - (nb1_p[0] + now_p[0]) * (now_p[0] * nb1_p[1] - nb1_p[0] * now_p[1])
                - (now_p[0] + nb2_p[0]) * (nb2_p[0] * now_p[1] - now_p[0] * nb2_p[1]);

        const double del_y
                = (nb1_p[1] + mov_p[1]) * (mov_p[0] * nb1_p[1] - nb1_p[0] * mov_p[1])
                + (mov_p[1] + nb2_p[1]) * (nb2_p[0] * mov_p[1] - mov_p[0] * nb2_p[1])
                - (nb1_p[1] + now_p[1]) * (now_p[0] * nb1_p[1] - nb1_p[0] * now_p[1])
                - (now_p[1] + nb2_p[1]) * (nb2_p[0] * now_p[1] - now_p[0] * nb2_p[1]);

        const double del_xx
                = (mov_p[0] * mov_p[0] + nb1_p[0] * mov_p[0] + nb1_p[0] * nb1_p[0])
                        * (mov_p[0] * nb1_p[1] - nb1_p[0] * mov_p[1])
                + (nb2_p[0] * nb2_p[0] + mov_p[0] * nb2_p[0] + mov_p[0] * mov_p[0])
                        * (nb2_p[0] * mov_p[1] - mov_p[0] * nb2_p[1])
                - (now_p[0] * now_p[0] + nb1_p[0] * now_p[0] + nb1_p[0] * nb1_p[0])
                        * (now_p[0] * nb1_p[1] - nb1_p[0] * now_p[1])
                - (nb2_p[0] * nb2_p[0] + now_p[0] * nb2_p[0] + now_p[0] * now_p[0])
                        * (nb2_p[0] * now_p[1] - now_p[0] * nb2_p[1]);

        const double del_xy
                = (nb1_p[0] * mov_p[1] - mov_p[0] * nb1_p[1])
                        * (mov_p[0] * (2 * mov_p[1] + nb1_p[1]) + nb1_p[0] * (mov_p[1] + 2.0 * nb1_p[1]))
                + (mov_p[0] * nb2_p[1] - nb2_p[0] * mov_p[1])
                        * (nb2_p[0] * (2 * nb2_p[1] + mov_p[1]) + mov_p[0] * (nb2_p[1]  + 2.0 * mov_p[1]))
                - (nb1_p[0] * now_p[1] - now_p[0] * nb1_p[1])
                        * (now_p[0] * (2 * now_p[1] + nb1_p[1]) + nb1_p[0] * (now_p[1] + 2.0 * nb1_p[1]))
                - (now_p[0] * nb2_p[1] - nb2_p[0] * now_p[1])
                        * (nb2_p[0] * (2 * nb2_p[1] + now_p[1]) + now_p[0] * (nb2_p[1] + 2.0 * now_p[1]));

        const double del_yy
                = (mov_p[0] * nb1_p[1] - nb1_p[0] * mov_p[1])
                        * (mov_p[1] * mov_p[1] + nb1_p[1] * mov_p[1] + nb1_p[1] * nb1_p[1])
                + (nb2_p[0] * mov_p[1] - mov_p[0] * nb2_p[1])
                        * (nb2_p[1] * nb2_p[1] + mov_p[1] * nb2_p[1] + mov_p[1] * mov_p[1])
                - (now_p[0] * nb1_p[1] - nb1_p[0] * now_p[1])
                        * (now_p[1] * now_p[1] + nb1_p[1] * now_p[1] + nb1_p[1] * nb1_p[1])
                - (nb2_p[0] * now_p[1] - now_p[0] * nb2_p[1])
                                * (nb2_p[1] * nb2_p[1] + now_p[1] * nb2_p[1] + now_p[1] * now_p[1]);

        return GeoData(
                new_area,
                (old_geo.m_centroid * 6.0 * old_geo.m_area - std::array<double, 3> {{del_x, del_y, 0.0 }}) / (6.0 * new_area),
                old_geo.m_area_moment - AreaMoment(del_xx/ 12.0, del_xy / 24.0, del_yy / 12.0)
                );
}

inline GeoData operator+(const GeoData& i1, const GeoData& i2)
{
        return GeoData(i1.m_area + i2.m_area,
                i1.m_centroid + i2.m_centroid, i1.m_area_moment + i2.m_area_moment);
}

inline GeoData operator-(const GeoData& i1, const GeoData& i2)
{
        return GeoData(i1.m_area - i2.m_area,
                i1.m_centroid - i2.m_centroid, i1.m_area_moment - i2.m_area_moment);
}

inline std::ostream& operator<<(std::ostream& os, const GeoData& c)
{
        os << "GeoData: {" << std::endl
                << "\t area: " << c.m_area << std::endl
                << "\t centroid: " << c.m_centroid << std::endl
                << "\t area moment of inertia: " << c.m_area_moment << std::endl
                << "}";

        return os;
}

}

#endif //end_of_include_guard