Commit 96d808c1 authored by Tomáš Jakubec's avatar Tomáš Jakubec
Browse files

UnstructuredMesh tests

parent 0304c8b2
Loading
Loading
Loading
Loading
+1 −0
Original line number Diff line number Diff line
@@ -61,6 +61,7 @@ DISTFILES += \
    ../src/UnitTests/Traits/ArithmeticTraitsTest.cpp \
    ../src/UnitTests/Traits/TraitsTest.cpp \
    ../src/UnitTests/UnstructuredMesh/MeshDataContainerTest.cpp \
    ../src/UnitTests/UnstructuredMesh/UnstructuredMeshTest.cpp \
    ../README.md \
    ../src/GTMesh/Debug/README.md \
    ../src/GTMesh/Traits/TraitsAlgorithm/README.md \
+76 −115
Original line number Diff line number Diff line
@@ -390,11 +390,20 @@ void testMesh2D() {

    DBGMSG("2D cells distances");

    auto distances = ComputeCellsDistance(mesh);
    auto distances = computeCellsDistance(mesh);
    for(auto& edge : mesh.getEdges()){
        DBGVAR(edge.getIndex(),distances.at(edge));
    }

    DBGVAR((mesh.connections<2,0>().getDataByPos<0>()),
           (mesh.connections<2,1>().getDataByPos<0>()),
           (mesh.connections<1,2>().getDataByPos<0>()),
           (mesh.connections<0,2>().getDataByPos<0>()));

    DBGVAR((mesh.neighborhood<2,0>().getDataByPos<0>()),
           (mesh.neighborhood<2,1,0>().getDataByPos<0>()),
           (mesh.neighborhood<1,2,0>().getDataByPos<0>()),
           (mesh.neighborhood<0,2>().getDataByPos<0>()));
}


@@ -411,10 +420,7 @@ void testMesh2DLoadAndWrite(){
    DBGVAR(mesh.getVertices().size(), mesh.getVertices().at(4),mesh.getCells().size());


    DBGMSG("mesh apply test");
    Impl::MeshRun<2, 2, 0, 2,false, true>::run(mesh,size_t(4), size_t(4), [](size_t ori, size_t i){
        DBGVAR(ori,i);
    });
    DBGVAR((mesh.connections<2,0>().getDataByPos<0>()[4]));

    mesh.initializeCenters();
    auto normals = mesh.computeFaceNormals();
@@ -718,140 +724,95 @@ struct doubleEdgeNormalData {
};
MAKE_CUSTOM_TRAIT(doubleEdgeNormalData, "firstEdgeNormal", std::make_pair(&doubleEdgeNormalData::get, &doubleEdgeNormalData::set));

template<unsigned int Dim>
struct CellData {
    unsigned int color;
    Vertex<Dim, double> center;
};

MAKE_ATTRIBUTE_TEMPLATE_TRAIT((CellData<Dim>), (unsigned int Dim), center, color);


void testMeshRefine() {
    UnstructuredMesh<3, size_t, double, 6> mesh;
    twoPrisms(mesh);
    mesh.initializeCenters();

    auto col = mesh.coloring<3,1>();

    MeshDataContainer<CellData<3>, 3> cellData(mesh);
    for (auto& cell : mesh.getCells()){
        cellData[cell].color = col[cell];
        cellData[cell].center = cell.getCenter();
    }

    MeshDataContainer<MeshNativeType<3>::ElementType,3> types(mesh, MeshNativeType<3>::WEDGE);
    VTKMeshWriter<3, size_t, double> writer;
    ofstream out3D;
    std::ofstream out3D;
    out3D.open("mesh_refine_0.vtk");
    DBGVAR(bool(out3D));
    writer.writeHeader(out3D, "test data");
    writer.writeToStream(out3D, mesh, types);

    auto colours = MeshColoring<3,0>::color(mesh);

    out3D << "CELL_DATA " << mesh.getCells().size() << endl;
    out3D << "SCALARS cell_wrt_vertex_colour double 1\nLOOKUP_TABLE default" << endl;
    for (auto colour : colours.getDataByPos<0>()) {
        out3D << colour << ' ';
    }
    VTKMeshDataWriter<3>::writeToStream(out3D, cellData, writer);
    out3D.close();

    MeshDataContainer<MeshNativeType<3>::ElementType,3> types1(mesh, MeshNativeType<3>::POLYHEDRON);

    VTKMeshWriter<3, size_t, double> writer1;
    out3D.open("mesh_refine_1.vtk");
    writer1.writeHeader(out3D, "test data");
    writer1.writeToStream(out3D, mesh, types1);
    auto colours1 = MeshColoring<3,0>::color(mesh);

    MeshDataContainer<colorData, 3> cd(mesh);
    auto normals = mesh.computeFaceNormals();

    for(auto& cell : mesh.getCells()){
        cd.at(cell).color = colours1.at(cell);
        cd.at(cell).firstEdgeNormal = normals.getDataByDim<2>().at(mesh.getFaces().at(cell.getBoundaryElementIndex()).getNextBElem(cell.getIndex()));
    }
    DBGVAR(cd.getDataByDim<3>());

    VTKMeshDataWriter<3> dataWriter;

    cd.getDataByPos<0>()[0] = cd.getDataByPos<0>()[0] + cd.getDataByPos<0>()[1];
    dataWriter.writeToStream(out3D, cd, writer1);


    //out3D << "CELL_DATA " << writer1.cellVert.getDataByPos<0>().size() << endl;
    out3D << "SCALARS cell_wrt_vertex_colour double 1\nLOOKUP_TABLE default" << endl;
    size_t realIndex = 0;
    for (size_t i = 0; i < writer1.cellVert.getDataByPos<0>().size(); i++) {
        auto iterator = writer1.backwardCellIndexMapping.find(i);
        if (iterator == writer1.backwardCellIndexMapping.end()){
            out3D << colours1.getDataByPos<0>().at(realIndex) << ' ';
            realIndex++;
        } else {
            out3D << colours1.getDataByPos<0>().at(iterator->second) << ' ';
            realIndex = iterator->second;
        }
    }
    out3D.close();
    mesh.clear();


    auto reader_ptr = mesh.load("mesh_refine_1.vtk");
    DBGVAR(reader_ptr->getCellTypes().getDataByPos<0>());
    ifstream in3D;
    in3D.open("mesh_refine_1.vtk", std::ios::binary);
    VTKMeshReader<3> reader;
    //reader.loadFromStream(in3D, mesh);


    MeshDataContainer<colorData, 3> cd1(mesh);
    VTKMeshDataReader<3, size_t>::readFromStream(in3D, cd1);

    DBGVAR(cd1.getDataByDim<3>());

    MeshDataContainer<std::tuple<doubleColorData, doubleEdgeNormalData>, 3,3> cd2(mesh);
    VTKMeshDataReader<3, size_t>::readFromStream(in3D, cd2);
    std::ifstream ifst("mesh_refine_0.vtk", ios::binary | ios::in);
    DBGVAR(bool(ifst));
    reader.loadFromStream(ifst, mesh);
    MeshDataContainer<CellData<3>, 3> cellDataLoad(mesh);
    VTKMeshDataReader<3, size_t>::readFromStream(ifst, cellDataLoad);
    ifst.close();

    DBGVAR(cd2.getDataByPos<0>(), cd2.getDataByPos<1>());
    mesh.initializeCenters();
    // compute centers
    auto centers = computeCenters<METHOD_DEFAULT>(mesh);
    std::vector<Vertex<3, double>> expectCenter = { {0.333333, 0.333333, 0.5}, {0.666667, 0.666667, 0.5} };
    DBGVAR((centers.getDataByDim<3>()), expectCenter);

    in3D.close();
    // measures test
    auto measures = mesh.computeElementMeasures();
    std::vector<double> expectEdgeM = { 1, 1.41421, 1, 1, 1, 1, 1, 1.41421, 1, 1, 1, 1, 1, 1 };
    DBGVAR(measures.getDataByDim<1>(), expectEdgeM);
    std::vector<double> expectFaceM = {  0.5, 1, 1.41421, 1, 0.5, 0.5, 1, 1, 0.5 };
    DBGVAR(measures.getDataByDim<2>(), expectFaceM);
    std::vector<double> expectCellM = {  0.5, 0.5 };
    DBGVAR(measures.getDataByDim<3>(), expectCellM);

    for (auto& cell : mesh.getCells()){
        DBGVAR(cellDataLoad[cell].color, cellData[cell].color);
        DBGVAR(cellDataLoad[cell].center, cellData[cell].center);
    }


    mesh.initializeCenters();

    MeshDataContainer<MeshNativeType<3>::ElementType,3> types2(mesh, MeshNativeType<3>::POLYHEDRON);
    out3D.open("mesh_refine_2.vtk");
    MeshDataContainer<MeshNativeType<3>::ElementType,3> types1(mesh, MeshNativeType<3>::POLYHEDRON);
    VTKMeshWriter<3, size_t, double> writer1;
    out3D.open("mesh_refine_1.vtk");
    DBGVAR(bool(out3D));
    writer1.writeHeader(out3D, "test data");
    writer1.writeToStream(out3D, mesh, types2);

    auto colours2 = MeshColoring<3,0>::color(mesh);

    out3D << "CELL_DATA " << writer1.cellVert.getDataByPos<0>().size() << endl;
    out3D << "SCALARS cell_wrt_vertex_colour double 1\nLOOKUP_TABLE default" << endl;
    realIndex = 0;
    for (size_t i = 0; i < writer1.cellVert.getDataByPos<0>().size(); i++) {
        auto iterator = writer1.backwardCellIndexMapping.find(i);
        if (iterator == writer1.backwardCellIndexMapping.end()){
            out3D << colours2.getDataByPos<0>().at(realIndex) << ' ';
            realIndex++;
        } else {
            out3D << colours2.getDataByPos<0>().at(iterator->second) << ' ';
            realIndex = iterator->second;
        }
    }
    writer1.writeToStream(out3D, mesh, types1);
    VTKMeshDataWriter<3>::writeToStream(out3D, cellData, writer1);
    out3D.close();

    UnstructuredMesh<3, size_t, double, 6> meshRefined;

    in3D.open("mesh_refine_2.vtk");
    reader.loadFromStream(in3D, mesh);
    in3D.close();

    mesh.initializeCenters();
    ifst.open("mesh_refine_1.vtk", ios::binary | ios::in);
    DBGVAR(bool(ifst));
    reader.loadFromStream(ifst, meshRefined);
    DBGVAR(meshRefined.getCells().size());
    MeshDataContainer<CellData<3>, 3> cellDataLoadRefined(meshRefined);
    VTKMeshDataReader<3, size_t>::readFromStream(ifst, cellDataLoadRefined);
    ifst.close();

    MeshDataContainer<MeshNativeType<3>::ElementType,3> types3(mesh, MeshNativeType<3>::POLYHEDRON);
    out3D.open("mesh_refine_3.vtk");
    writer1.writeHeader(out3D, "test data");
    writer1.writeToStream(out3D, mesh, types3);
    auto colours3 = MeshColoring<3,0>::color(mesh);

    out3D << "CELL_DATA " << writer1.cellVert.getDataByPos<0>().size() << endl;
    out3D << "SCALARS cell_wrt_vertex_colour double 1\nLOOKUP_TABLE default" << endl;
    realIndex = 0;
    for (size_t i = 0; i < writer1.cellVert.getDataByPos<0>().size(); i++) {
        auto iterator = writer1.backwardCellIndexMapping.find(i);
        if (iterator == writer1.backwardCellIndexMapping.end()){
            out3D << colours3.getDataByPos<0>().at(realIndex) << ' ';
            realIndex++;
        } else {
            out3D << colours3.getDataByPos<0>().at(iterator->second) << ' ';
            realIndex = iterator->second;
        }
    for (auto& cell : meshRefined.getCells()){
        DBGVAR(cellDataLoadRefined[cell].color, cellData.getDataByPos<0>()[writer1.backwardCellIndexMapping[cell.getIndex()]].color);
        DBGVAR(cellDataLoadRefined[cell].center, cellData.getDataByPos<0>()[writer1.backwardCellIndexMapping[cell.getIndex()]].center);
    }
    out3D.close();

}


@@ -1091,11 +1052,11 @@ int main()
    //testMesh2DLoadAndWrite();
    //testMesh3D();
    //test3DMeshDeformedPrisms();
    //testMeshRefine();
    testMeshDataContainer();
    testMeshRefine();
    //testMeshDataContainer();
    //UnstructuredMesh<5, size_t, double, 6,5,4> m;
    //m.ComputeElementMeasures();
    test3DMeshLoad();
    //test3DMeshLoad();

    //testFPMA_poly();

+1 −1
Original line number Diff line number Diff line
@@ -15,7 +15,7 @@


template <unsigned int Dimension,typename IndexType, typename Real, unsigned int ...Reserve>
MeshDataContainer<Real, Dimension-1> ComputeCellsDistance(const MeshElements<Dimension, IndexType, Real, Reserve...>& mesh){
MeshDataContainer<Real, Dimension-1> computeCellsDistance(const MeshElements<Dimension, IndexType, Real, Reserve...>& mesh){

    MeshDataContainer<Real, Dimension-1> distances(mesh);

+16 −16
Original line number Diff line number Diff line
@@ -15,7 +15,7 @@ namespace Impl {
template <unsigned int CurrentDimension, unsigned int MeshDimension, ComputationMethod Method = METHOD_DEFAULT>
struct _ComputeMeasures{
    template <typename IndexType, typename Real, unsigned int ...Reserve>
    static void compute(MakeMeshDataContainer_t<Real, make_custom_integer_sequence_t<unsigned int, 1, MeshDimension>>&,MeshElements<MeshDimension, IndexType, Real, Reserve...>&){
    static void compute(MakeMeshDataContainer_t<Real, make_custom_integer_sequence_t<unsigned int, 1, MeshDimension>>&, const MeshElements<MeshDimension, IndexType, Real, Reserve...>&){
        static_assert (MeshDimension <= 3,"The measure computation of mesh of dimension higher than 3 is not implemented yet.");
    }
};
@@ -24,7 +24,7 @@ struct _ComputeMeasures{
template <unsigned int MeshDimension, ComputationMethod Method>
struct _ComputeMeasures<1, MeshDimension, Method>{
    template <typename IndexType, typename Real, unsigned int ...Reserve>
    static void compute(MakeMeshDataContainer_t<Real, make_custom_integer_sequence_t<unsigned int, 1, MeshDimension>>& measures,MeshElements<MeshDimension, IndexType, Real, Reserve...>& mesh){
    static void compute(MakeMeshDataContainer_t<Real, make_custom_integer_sequence_t<unsigned int, 1, MeshDimension>>& measures, const MeshElements<MeshDimension, IndexType, Real, Reserve...>& mesh){

        auto& edgeLengths = measures.template getDataByDim<1>();

@@ -42,7 +42,7 @@ struct _ComputeMeasures<1, MeshDimension, Method>{
template <ComputationMethod Method>
struct _ComputeMeasures<3, 3, Method>{
    template <typename IndexType, typename Real, unsigned int ...Reserve>
    static void compute(MakeMeshDataContainer_t<Real, make_custom_integer_sequence_t<unsigned int, 1, 3>>& measures,MeshElements<3, IndexType, Real, Reserve...>& mesh){
    static void compute(MakeMeshDataContainer_t<Real, make_custom_integer_sequence_t<unsigned int, 1, 3>>& measures, const MeshElements<3, IndexType, Real, Reserve...>& mesh){

        auto& cellMeasures = measures.template getDataByDim<3>();

@@ -58,9 +58,9 @@ struct _ComputeMeasures<3, 3, Method>{
                IndexType vAIndex = mesh.getEdges().at(mesh.getFaces().at(tmpFace).getSubelements()[0]).getVertexAIndex();
                IndexType vBIndex = mesh.getEdges().at(mesh.getFaces().at(tmpFace).getSubelements()[0]).getVertexBIndex();

                Vertex<3,Real>& a = mesh.getVertices().at(vAIndex);
                Vertex<3,Real>& b = mesh.getVertices().at(vBIndex);
                Vertex<3,Real>& c = mesh.getFaces().at(tmpFace).getCenter();
                const Vertex<3,Real>& a = mesh.getVertices().at(vAIndex);
                const Vertex<3,Real>& b = mesh.getVertices().at(vBIndex);
                const Vertex<3,Real>& c = mesh.getFaces().at(tmpFace).getCenter();

                std::array<Vertex<3,Real>, 3> gsVecs = {a - b, a - c, a - cellCenter};
                std::array<Real, 3> gsNorms = {};
@@ -83,18 +83,18 @@ struct _ComputeMeasures<3, 3, Method>{
template <ComputationMethod Method>
struct _ComputeMeasures<2, 2, Method>{
    template <typename IndexType, typename Real, unsigned int ...Reserve>
    static void compute(MakeMeshDataContainer_t<Real, make_custom_integer_sequence_t<unsigned int, 1, 2>>& measures,MeshElements<2, IndexType, Real, Reserve...>& mesh){
    static void compute(MakeMeshDataContainer_t<Real, make_custom_integer_sequence_t<unsigned int, 1, 2>>& measures, const MeshElements<2, IndexType, Real, Reserve...>& mesh){

        auto& surfaceMeasures = measures.template getDataByDim<2>();

        for (IndexType cellIndex = 0; cellIndex < mesh.getCells().size(); cellIndex++) {
            typename MeshElements<2, IndexType, Real, Reserve...>::template ElementType<2>& cell = mesh.getCells().at(cellIndex);
            auto& cell = mesh.getCells().at(cellIndex);
            IndexType tmpEdge = cell.getBoundaryElementIndex();
            Real measure = Real();
            Vertex<2,Real>& cellCenter = cell.getCenter();
            const Vertex<2,Real>& cellCenter = cell.getCenter();
            do {
                Vertex<2,Real>& a = mesh.getVertices().at(mesh.getEdges().at(tmpEdge).getVertexAIndex());
                Vertex<2,Real>& b = mesh.getVertices().at(mesh.getEdges().at(tmpEdge).getVertexBIndex());
                const Vertex<2,Real>& a = mesh.getVertices().at(mesh.getEdges().at(tmpEdge).getVertexAIndex());
                const Vertex<2,Real>& b = mesh.getVertices().at(mesh.getEdges().at(tmpEdge).getVertexBIndex());
                double tmp = (cellCenter[0] - a[0]) * (b[1] - a[1]);
                tmp -= (cellCenter[1] - a[1]) * (b[0] - a[0]);
                measure += 0.5 * fabs(tmp);
@@ -112,7 +112,7 @@ struct _ComputeMeasures<2, 2, Method>{
template <ComputationMethod Method>
struct _ComputeMeasures<2, 3, Method>{
    template <typename IndexType, typename Real, unsigned int ...Reserve>
    static void compute(MakeMeshDataContainer_t<Real, make_custom_integer_sequence_t<unsigned int, 1, 3>>& measures,MeshElements<3, IndexType, Real, Reserve...>& mesh){
    static void compute(MakeMeshDataContainer_t<Real, make_custom_integer_sequence_t<unsigned int, 1, 3>>& measures, const MeshElements<3, IndexType, Real, Reserve...>& mesh){

        auto& surfaceMeasures = measures.template getDataByDim<2>();

@@ -123,8 +123,8 @@ struct _ComputeMeasures<2, 3, Method>{
            const Vertex<3,Real>& faceCenter = face.getCenter();
            for(auto sube : face.getSubelements()){
                const auto& edge = mesh.getEdges().at(sube);
                Vertex<3,Real>& a = mesh.getVertices().at(edge.getVertexAIndex());
                Vertex<3,Real>& b = mesh.getVertices().at(edge.getVertexBIndex());
                const Vertex<3,Real>& a = mesh.getVertices().at(edge.getVertexAIndex());
                const Vertex<3,Real>& b = mesh.getVertices().at(edge.getVertexBIndex());

                std::array<Vertex<3,Real>, 2> gsVecs = {b - a, faceCenter - a};
                std::array<Real, 2> gsNorms = {};
@@ -149,7 +149,7 @@ struct _ComputeMeasures<2, 3, Method>{
template <>
struct _ComputeMeasures<3, 3, METHOD_TESSELLATED>{
    template <typename IndexType, typename Real, unsigned int ...Reserve>
    static void compute(MakeMeshDataContainer_t<Real, make_custom_integer_sequence_t<unsigned int, 1, 3>>& measures,MeshElements<3, IndexType, Real, Reserve...>& mesh){
    static void compute(MakeMeshDataContainer_t<Real, make_custom_integer_sequence_t<unsigned int, 1, 3>>& measures, const MeshElements<3, IndexType, Real, Reserve...>& mesh){

        auto& cellMeasures = measures.template getDataByDim<3>();

@@ -191,7 +191,7 @@ struct _ComputeMeasures<3, 3, METHOD_TESSELLATED>{


template <ComputationMethod Method, unsigned int MeshDimension,typename IndexType, typename Real, unsigned int ...Reserve>
MakeMeshDataContainer_t<Real, make_custom_integer_sequence_t<unsigned int, 1, MeshDimension>> computeMeasures(MeshElements<MeshDimension, IndexType, Real, Reserve...>& mesh){
MakeMeshDataContainer_t<Real, make_custom_integer_sequence_t<unsigned int, 1, MeshDimension>> computeMeasures(const MeshElements<MeshDimension, IndexType, Real, Reserve...>& mesh){
    MakeMeshDataContainer_t<Real, make_custom_integer_sequence_t<unsigned int, 1, MeshDimension>> measures(mesh);

    Impl::_ComputeMeasures<1, MeshDimension, Method>::compute(measures, mesh);
+5 −5
Original line number Diff line number Diff line
@@ -15,7 +15,7 @@ template<unsigned int ColoredDim, unsigned int ConnectingDim, ColoringMethod Met
struct _MeshColoring {
    template<unsigned int MeshDimension, typename IndexType, typename Real, unsigned int ...Reserve>
    static MeshDataContainer<unsigned int, ColoredDim> color(
            MeshElements<MeshDimension, IndexType, Real, Reserve...>& mesh
            const MeshElements<MeshDimension, IndexType, Real, Reserve...>& mesh
            ) {
        // resulting container of colours
        MeshDataContainer<unsigned int, ColoredDim> result(mesh);
@@ -70,7 +70,7 @@ template<unsigned int ColoredDim, unsigned int ConnectingDim>
struct _MeshColoring <ColoredDim, ConnectingDim, METHOD_GREEDY, true> {
    template<unsigned int MeshDimension, typename IndexType, typename Real, unsigned int ...Reserve>
    static MeshDataContainer<unsigned int, ColoredDim> color(
            MeshElements<MeshDimension, IndexType, Real, Reserve...>& mesh
            const MeshElements<MeshDimension, IndexType, Real, Reserve...>& mesh
            ) {
        MeshDataContainer<unsigned int, ColoredDim> result(mesh);

@@ -120,7 +120,7 @@ template<unsigned int ColoredDim, unsigned int ConnectingDim, bool Descend>
struct _MeshColoring <ColoredDim, ConnectingDim, METHOD_RANDOM, Descend> {
    template<unsigned int MeshDimension, typename IndexType, typename Real, unsigned int ...Reserve>
    static MeshDataContainer<unsigned int, ColoredDim> color(
            MeshElements<MeshDimension, IndexType, Real, Reserve...>& mesh,
            const MeshElements<MeshDimension, IndexType, Real, Reserve...>& mesh,
            unsigned int seed = 1562315
            ) {
        // calculating the reserve
@@ -195,7 +195,7 @@ struct MeshColoring{

template<unsigned int MeshDimension, typename IndexType, typename Real, unsigned int ...Reserve>
static MeshDataContainer<unsigned int, ColoredDim> color(
            MeshElements<MeshDimension, IndexType, Real, Reserve...>& mesh
            const MeshElements<MeshDimension, IndexType, Real, Reserve...>& mesh
        ){
    return Impl::_MeshColoring<ColoredDim, ConnectingDim, Method>::color(mesh);
}
@@ -207,7 +207,7 @@ struct MeshColoring<ColoredDim, ConnectingDim, METHOD_RANDOM>{

template<unsigned int MeshDimension, typename IndexType, typename Real, unsigned int ...Reserve>
static MeshDataContainer<unsigned int, ColoredDim> color(
            MeshElements<MeshDimension, IndexType, Real, Reserve...>& mesh,
            const MeshElements<MeshDimension, IndexType, Real, Reserve...>& mesh,
            unsigned int seed = 1562315
        ){
    return Impl::_MeshColoring<ColoredDim, ConnectingDim, METHOD_RANDOM>::color(mesh, seed);
Loading