MeshGeometry.cpp

Go to the documentation of this file.

/*
 * 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 "MeshGeometry.h"

#include "bio/Cell.h"
#include "bio/Node.h"
#include "bio/Mesh.h"
#include "math/ChainHull.h"

#include <memory>
#include <string>
#include <vector>

using namespace std;
using namespace SimPT_Sim::Util;

namespace SimPT_Sim {

tuple<array<double, 3>, array<double, 3>> MeshGeometry::BoundingBox(shared_ptr<Mesh> mesh)
{
        const auto& nodes = mesh->GetBoundaryPolygon()->GetNodes();
        array<double, 3>  lower_left = **(nodes.begin());
        array<double, 3>  upper_right(lower_left);

        for (auto const& node : nodes) {
                array<double, 3> a(*node);
                lower_left[0]  = min(lower_left[0],  a[0]);
                lower_left[1]  = min(lower_left[1],  a[1]);
                lower_left[2]  = min(lower_left[2],  a[2]);
                upper_right[0] = max(upper_right[0], a[0]);
                upper_right[1] = max(upper_right[1], a[1]);
                upper_right[2] = max(upper_right[2], a[2]);
        }
        return make_tuple(lower_left, upper_right);
}

tuple<double, double, double> MeshGeometry::Compactness(const Mesh* mesh)
{
        // Calculate compactness using the convex hull of the cells
        // We use Andrew's Monotone Chain Algorithm (see hull.cpp)

        // Step 1: Prepare data for 2D hull code: get boundary
        std::vector<std::array<double, 2> > boundary;
        const Cell*  bp = mesh->GetBoundaryPolygon();
        for (auto const& node : bp->GetNodes()) {
                boundary.push_back({ {(*node)[0], (*node)[1]} });
        }

        // Step 2: call 2D Hull code
        std::vector<std::array<double, 2> > hull = ChainHull::Compute(boundary);

        // Step 3: calculate area and circumference of convex hull
        double hull_area = 0.;
        double hull_circumference = 0.;
        for (unsigned int i = 0; i < hull.size() - 1; i++) {
                hull_area += hull[i][0] * hull[i + 1][1] - hull[i + 1][0] * hull[i][1];
                double s_dx = (hull[i + 1][0] - hull[i][0]);
                double s_dy = (hull[i + 1][1] - hull[i][1]);
                double l = sqrt(s_dx * s_dx + s_dy * s_dy);
                hull_circumference += l;
        }
        hull_area /= 2.;

        // Step 4: get area of boundary polygon
        const double boundary_pol_area = bp->GetArea();

        return make_tuple(boundary_pol_area / hull_area, hull_area, hull_circumference);
}

} // namespace