Triangulate IP Cells

CMakeLists.txt

add_unit_test(TriangulateIPCellsTest
    mechanics/integrationtypes/IntegrationTypeBase.cpp
    mechanics/integrationtypes/IntegrationTypeEnum.cpp
    mechanics/integrationtypes/IntegrationType2D3NGauss3Ip.cpp
    mechanics/integrationtypes/IntegrationType2D4NGauss1Ip.cpp
    mechanics/integrationtypes/IntegrationType2D4NGauss4Ip.cpp
    mechanics/integrationtypes/IntegrationType2D4NGauss9Ip.cpp)
target_link_libraries(TriangulateIPCellsTest -lgmp -lCGAL)

TriangulateIPCells.h

#include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <iterator>
#include "mechanics/integrationtypes/IntegrationTypeBase.h"

typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::Point_2 Point_2;
typedef K::Iso_rectangle_2 Iso_rectangle_2;
typedef K::Segment_2 Segment_2;
typedef K::Ray_2 Ray_2;
typedef K::Line_2 Line_2;
typedef CGAL::Delaunay_triangulation_2<K> Delaunay_triangulation_2;
typedef K::Triangle_2 Triangle_2;

// A class to recover Voronoi diagram from stream.
// Rays, lines and segments are cropped to a rectangle
// so that only segments are stored
struct Cropped_voronoi_from_delaunay {
    std::list<Segment_2> m_cropped_vd;
    Iso_rectangle_2 m_bbox;

    Cropped_voronoi_from_delaunay(const Iso_rectangle_2& bbox) : m_bbox(bbox) {}

    template <class RSL>
    void crop_and_extract_segment(const RSL& rsl) {
        CGAL::Object obj = CGAL::intersection(rsl, m_bbox);
        const Segment_2* s = CGAL::object_cast<Segment_2>(&obj);
        if (s) m_cropped_vd.push_back(*s);
    }

    void operator<<(const Ray_2& ray) { crop_and_extract_segment(ray); }
    void operator<<(const Line_2& line) { crop_and_extract_segment(line); }
    void operator<<(const Segment_2& seg) { crop_and_extract_segment(seg); }
};

struct Cropped_voronoi_from_delaunay_tri {
    std::list<Segment_2> m_cropped_vd;
    Triangle_2 m_bbox;

    Cropped_voronoi_from_delaunay_tri(const Triangle_2& bbox) : m_bbox(bbox) {}

    template <class RSL>
    void crop_and_extract_segment(const RSL& rsl) {
        CGAL::Object obj = CGAL::intersection(rsl, m_bbox);
        const Segment_2* s = CGAL::object_cast<Segment_2>(&obj);
        if (s) m_cropped_vd.push_back(*s);
    }

    void operator<<(const Ray_2& ray) { crop_and_extract_segment(ray); }
    void operator<<(const Line_2& line) { crop_and_extract_segment(line); }
    void operator<<(const Segment_2& seg) { crop_and_extract_segment(seg); }
};

void TriangulateIPCells(NuTo::IntegrationTypeBase& integration)
{
    int nPoints = integration.GetNumIntegrationPoints();

    std::vector<Point_2> points;
    for (int i = 0; i < nPoints; ++i)
    {
        auto point = integration.GetLocalIntegrationPointCoordinates(i);
        points.push_back(Point_2(point(0), point(1)));
    }

    Delaunay_triangulation_2 dt2;
    // insert points into the triangulation
    dt2.insert(points.begin(), points.end());
    // construct a rectangle
    Iso_rectangle_2 bbox(-1, -1, 1, 1);
    Cropped_voronoi_from_delaunay vor(bbox);
    // extract the cropped Voronoi diagram
    dt2.draw_dual(vor);
    // print the cropped Voronoi diagram as segments
    std::copy(vor.m_cropped_vd.begin(), vor.m_cropped_vd.end(),
              std::ostream_iterator<Segment_2>(std::cout, "\n"));
}


void TriangulateIP3NCells(NuTo::IntegrationTypeBase& integration)
{
    int nPoints = integration.GetNumIntegrationPoints();

    std::vector<Point_2> points;
    for (int i = 0; i < nPoints; ++i)
    {
        auto point = integration.GetLocalIntegrationPointCoordinates(i);
        points.push_back(Point_2(point(0), point(1)));
    }

    Delaunay_triangulation_2 dt2;
    // insert points into the triangulation
    dt2.insert(points.begin(), points.end());
    // construct a rectangle
    Triangle_2 bbox(Point_2(0, 0), Point_2(1, 0), Point_2(0, 1));
    Cropped_voronoi_from_delaunay_tri vor(bbox);
    // extract the cropped Voronoi diagram
    dt2.draw_dual(vor);
    // print the cropped Voronoi diagram as segments
    std::copy(vor.m_cropped_vd.begin(), vor.m_cropped_vd.end(),
              std::ostream_iterator<Segment_2>(std::cout, "\n"));
}

TriangulateIPCellsTest.cpp

#include "BoostUnitTest.h"

#include "visualize/TriangulateIPCells.h"
#include "mechanics/integrationtypes/IntegrationType2D3NGauss3Ip.h"
#include "mechanics/integrationtypes/IntegrationType2D4NGauss1Ip.h"
#include "mechanics/integrationtypes/IntegrationType2D4NGauss4Ip.h"
#include "mechanics/integrationtypes/IntegrationType2D4NGauss9Ip.h"

BOOST_AUTO_TEST_CASE(one)
{
    std::cout << "1 Ips:\n";
    NuTo::IntegrationType2D4NGauss1Ip integration;
    TriangulateIPCells(integration);
}

BOOST_AUTO_TEST_CASE(four)
{
    std::cout << "4 Ips:\n";
    NuTo::IntegrationType2D4NGauss4Ip integration;
    TriangulateIPCells(integration);
}

BOOST_AUTO_TEST_CASE(nine)
{
    std::cout << "9 Ips:\n";
    NuTo::IntegrationType2D4NGauss9Ip integration;
    TriangulateIPCells(integration);
}

BOOST_AUTO_TEST_CASE(triangle_three)
{
    std::cout << "Tri 3 Ips:\n";
    NuTo::IntegrationType2D3NGauss3Ip integration;
    TriangulateIP3NCells(integration);
}