Merge branch 'feature/gmsh-writer' into 'master'

Feature/gmsh writer

See merge request DennisGlaeser/frackit!15

frackit/CMakeLists.txt

@@ -3,6 +3,7 @@ add_subdirectory(distance)
 add_subdirectory(entitynetwork)
 add_subdirectory(geometry)
 add_subdirectory(intersection)
+add_subdirectory(io)
 add_subdirectory(magnitude)
 add_subdirectory(occ)
 add_subdirectory(precision)

frackit/io/CMakeLists.txt

+install(FILES
+brepwriter.hh
+gmshwriter.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/frackit/io)

frackit/io/brepwriter.hh

+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 2 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief Class that writes entity networks into the .brep file format.
+ */
+#ifndef FRACKIT_BREP_WRITER_HH
+#define FRACKIT_BREP_WRITER_HH
+
+#include <string>
+#include <stdexcept>
+#include <unordered_map>
+
+#include <BRepTools.hxx>
+#include <BRep_Builder.hxx>
+
+#include <TopTools_IndexedMapOfShape.hxx>
+#include <TopTools_DataMapOfShapeInteger.hxx>
+#include <TopTools_ListOfShape.hxx>
+
+#include <TopExp_Explorer.hxx>
+#include <TopoDS_Compound.hxx>
+#include <TopoDS_Shape.hxx>
+#include <TopoDS.hxx>
+
+#include <frackit/occ/breputilities.hh>
+#include <frackit/entitynetwork/entitynetwork.hh>
+#include <frackit/entitynetwork/containedentitynetwork.hh>
+
+namespace Frackit {
+
+/*!
+ * \brief Class that writes entity networks into the .brep file format.
+ *        This creates a single TopoDS_Compound shape in which each sub-shape
+ *        is uniquely defined.
+ */
+class BRepWriter
+{
+public:
+
+    /*!
+     * \brief Construction from an entity network.
+     */
+    BRepWriter(const EntityNetwork& network)
+    {
+        if (network.entityDimension() != 2)
+            throw std::runtime_error("BRepWriter only implemented for 2d networks so far");
+
+        makeEntityMap_(network);
+        findEntityIntersectionShapes_();
+        makeSubShapeMaps_(network);
+        makeCompound_();
+    }
+
+    /*!
+     * \brief Construction from a contained entity network.
+     */
+    BRepWriter(const ContainedEntityNetwork& network)
+    {
+        if (network.entityDimension() != 2)
+            throw std::runtime_error("BRepWriter only implemented for 2d networks so far");
+        if (network.domainDimension() != network.entityDimension() + 1)
+            throw std::runtime_error("BRepWriter expects entityDim = domainDim - 1");
+
+        makeEntityMap_(network);
+        findEntityIntersectionShapes_();
+        makeSubShapeMaps_(network);
+        makeCompound_();
+    }
+
+    /*!
+     * \brief Write the entity network to disk.
+     */
+    void write(const std::string& fileName)
+    {
+        const std::string brepFile = fileName + ".brep";
+        BRepTools::Write(compound_, brepFile.c_str());
+    }
+
+protected:
+    /*!
+     * \brief Returns the map that maps a primary entity index to the list
+     *        of fragment indices that were inserted to the compound for
+     *        that primary entity index. The indices refer to the indices
+     *        within the compound.
+     */
+    const std::unordered_map<std::size_t, std::vector<std::size_t>>& entityToFragmentsMap() const
+    { return entityToFragmentsMap_; }
+
+    /*!
+     * \brief Returns the map that maps a primary entity intersectionindex to the
+     *        list of intersection fragments that were inserted to the compound for
+     *        that intersection. The indices refer to the indices within the compound.
+     */
+    const std::unordered_map<std::size_t, std::vector<std::size_t>>& entityIntersectionsToFragmentsMap() const
+    { return entityIntersectionsToFragmentsMap_; }
+
+    /*!
+     * \brief Returns the map that maps to each primary sub-domain index the list
+     *        of domain fragment indices that were inserted to the compound for
+     *        that sub-domain. The indices refer to the indices within the compound.
+     */
+    const std::unordered_map<std::size_t, std::vector<std::size_t>>& domainToFragmentsMap() const
+    { return domainToFragmentsMap_; }
+
+private:
+    /*!
+     * \brief Copies all entities of a network into a single map.
+     * \note This overload is for contained entity networks.
+     */
+    void makeEntityMap_(const ContainedEntityNetwork& network)
+    {
+        std::size_t subNetworkOffset = 0;
+        for (auto idx : network.subDomainIndices())
+        {
+            // the map maps to the entity index within the sub-domain network
+            // Therefore, we add the offset in order to ensure that we don't add
+            // entities from another network to the same primary entity index
+            const auto& map = network.subDomainEntityFragmentsIndexMap(idx);
+            for (TopTools_DataMapIteratorOfDataMapOfShapeInteger it(map); it.More(); it.Next())
+            {
+                allEntities_.Append(it.Key());
+                allEntitiesIndexMap_.Bind(it.Key(), it.Value() + subNetworkOffset);
+            }
+
+            // the next sub-domain
+            subNetworkOffset += map.Extent();
+        }
+    }
+
+    /*!
+     * \brief Copies all entities of a network into a single map.
+     * \note This overload is for uncontained entity networks.
+     */
+    void makeEntityMap_(const EntityNetwork& network)
+    {
+        const auto& map = network.entityFragmentsIndexMap();
+        for (TopTools_DataMapIteratorOfDataMapOfShapeInteger it(map); it.More(); it.Next())
+        {
+            allEntities_.Append(it.Key());
+            allEntitiesIndexMap_.Bind(it.Key(), it.Value());
+        }
+    }
+
+    /*!
+     * \brief Determins the shapes that describe intersections of entities.
+     */
+    void findEntityIntersectionShapes_()
+    {
+        // TODO: Start second iterator from ++it, somehow. A first implementation
+        //       of that led to some intersection edges not being meshed at the end
+        for (TopTools_ListIteratorOfListOfShape it(allEntities_); it.More(); it.Next())
+        {
+            for (TopTools_ListIteratorOfListOfShape it2(allEntities_); it2.More(); it2.Next())
+            {
+                if (it.Value().IsSame(it2.Value()))
+                    continue;
+
+                const auto edges = OCCUtilities::getOverlapEdges(TopoDS::Face(it.Value()),
+                                                                 TopoDS::Face(it2.Value()));
+
+                std::vector<std::size_t> handled;
+                for (unsigned int eIdx = 0; eIdx < edges.size(); ++eIdx)
+                    for (const auto& isList : entityIntersections_)
+                        if (isList.Contains(edges[eIdx]))
+                            handled.push_back(eIdx);
+
+                if (handled.size() != edges.size())
+                {
+                    entityIntersections_.resize(entityIntersections_.size()+1);
+                    for (unsigned int eIdx = 0; eIdx < edges.size(); ++eIdx)
+                        if (!std::count(handled.begin(), handled.end(), eIdx))
+                            entityIntersections_.back().Append(edges[eIdx]);
+                }
+            }
+        }
+    }
+
+    /*!
+     * \brief Fill the sub-shape maps with all shapes that describe the network.
+     * \note This overload is for contained entity networks.
+     */
+    void makeSubShapeMaps_(const ContainedEntityNetwork& network)
+    {
+        for (auto idx : network.subDomainIndices())
+        {
+            const auto& sdFragments = network.subDomainFragments(idx);
+            for (TopTools_ListIteratorOfListOfShape it(sdFragments); it.More(); it.Next())
+                addSolids_(it.Value(), idx);
+        }
+
+        for (auto idx : network.subDomainIndices())
+        {
+            const auto& sdEntityFragments = network.subDomainEntityFragments(idx);
+            for (TopTools_ListIteratorOfListOfShape it(sdEntityFragments); it.More(); it.Next())
+                addFaces_(it.Value());
+        }
+    }
+
+    /*!
+     * \brief Fill the sub-shape maps with all shapes that describe the network.
+     * \note This overload is for uncontained entity networks.
+     */
+    void makeSubShapeMaps_(const EntityNetwork& network)
+    {
+        const auto& sdEntityFragments = network.entityFragments();
+        for (TopTools_ListIteratorOfListOfShape it(sdEntityFragments); it.More(); it.Next())
+            addFaces_(it.Value());
+    }
+
+    /*!
+     * \brief Construct the single compound describing the network.
+     */
+    void makeCompound_()
+    {
+        BRep_Builder b;
+        b.MakeCompound(compound_);
+        for(std::size_t i = 1; i <= vmap_.Extent(); i++) b.Add(compound_, vmap_(i));
+        for(std::size_t i = 1; i <= emap_.Extent(); i++) b.Add(compound_, emap_(i));
+        for(std::size_t i = 1; i <= wmap_.Extent(); i++) b.Add(compound_, wmap_(i));
+        for(std::size_t i = 1; i <= fmap_.Extent(); i++) b.Add(compound_, fmap_(i));
+        for(std::size_t i = 1; i <= shmap_.Extent(); i++) b.Add(compound_, shmap_(i));
+        for(std::size_t i = 1; i <= somap_.Extent(); i++) b.Add(compound_, somap_(i));
+    }
+
+    /*!
+     * \brief Add solids to the sub-shape maps
+     * \param shape The shape of which the solids are to be extracted
+     * \param subDomainIdx The index of the sub-domain these solids are embedded in
+     */
+    void addSolids_(const TopoDS_Shape& shape, std::size_t subDomainIdx)
+    {
+        for(TopExp_Explorer solidExp(shape, TopAbs_SOLID); solidExp.More(); solidExp.Next())
+        {
+            TopoDS_Solid solid = TopoDS::Solid(solidExp.Current());
+            if(somap_.FindIndex(solid) < 1) // solid not yet in map
+            {
+                const auto curEntityIdx = somap_.Add(solid);
+                domainToFragmentsMap_[subDomainIdx].push_back(curEntityIdx);
+
+                for(TopExp_Explorer shellExp(solid, TopAbs_SHELL); shellExp.More(); shellExp.Next())
+                {
+                    TopoDS_Shell shell = TopoDS::Shell(shellExp.Current());
+                    if(shmap_.FindIndex(shell) < 1)
+                    {
+                        shmap_.Add(shell);
+                        addFaces_(shellExp.Current());
+                    }
+                }
+            }
+        }
+    }
+
+    /*!
+     * \brief Add faces to the sub-shape maps
+     * \param shape The shape of which the faces are to be extracted
+     * \todo TODO: For support of 1d networks (in 2d space) we would have
+     *             to add an overload taking the sub-domain index which is used
+     *             in that case and which would not search for the entity index.
+     */
+    void addFaces_(const TopoDS_Shape& shape)
+    {
+        for(TopExp_Explorer faceExp(shape, TopAbs_FACE); faceExp.More(); faceExp.Next())
+        {
+            TopoDS_Face face = TopoDS::Face(faceExp.Current());
+            const auto faceIdx = fmap_.FindIndex(face);
+
+            if (faceIdx < 1)
+            {
+                const auto curEntityIdx = fmap_.Add(face);
+
+                int entityIdx;
+                if (allEntitiesIndexMap_.Find(faceExp.Current(), entityIdx))
+                    entityToFragmentsMap_[entityIdx].push_back(curEntityIdx);
+
+                for (TopExp_Explorer wireExp(face.Oriented(TopAbs_FORWARD), TopAbs_WIRE); wireExp.More(); wireExp.Next())
+                {
+                    TopoDS_Wire wire = TopoDS::Wire(wireExp.Current());
+                    if(wmap_.FindIndex(wire) < 1)
+                    {
+                        wmap_.Add(wire);
+                        addEdges_(wireExp.Current());
+                    }
+                }
+            }
+        }
+    }
+
+    /*!
+     * \brief Add edges to the sub-shape maps
+     * \param shape The shape of which the edges are to be extracted
+     * \todo TODO: For support of 1d networks (in 2d space) we would have
+     *             to distinguish here and search for the entity index instead
+     *             of the intersection index.
+     */
+    void addEdges_(const TopoDS_Shape& shape)
+    {
+        for (TopExp_Explorer edgeExp(shape, TopAbs_EDGE); edgeExp.More(); edgeExp.Next())
+        {
+            TopoDS_Edge edge = TopoDS::Edge(edgeExp.Current());
+            const auto edgeIdx = emap_.FindIndex(edge);
+
+            if(edgeIdx < 1)
+            {
+                const auto curEntityIdx = emap_.Add(edge);
+
+                // start from index 1 in the map
+                for (std::size_t isIdx = 0; isIdx < entityIntersections_.size(); ++isIdx)
+                    if (entityIntersections_[isIdx].Contains(edge))
+                        entityIntersectionsToFragmentsMap_[isIdx+1].push_back(curEntityIdx);
+
+                addVertices_(edgeExp.Current());
+            }
+        }
+    }
+
+    /*!
+     * \brief Add vertices to the sub-shape maps
+     * \param shape The shape of which the vertices are to be extracted
+     * \todo TODO: For support of 1d networks (in 2d space) we would have
+     *            to  check if the vertices describe entity intersections.
+     */
+    void addVertices_(const TopoDS_Shape& shape)
+    {
+        for (TopExp_Explorer vertExp(shape, TopAbs_VERTEX); vertExp.More(); vertExp.Next())
+        {
+            TopoDS_Vertex vertex = TopoDS::Vertex(vertExp.Current());
+            if(vmap_.FindIndex(vertex) < 1)
+                vmap_.Add(vertex);
+        }
+    }
+
+    // containers to store the primary entities
+    TopTools_ListOfShape allEntities_;
+    TopTools_DataMapOfShapeInteger allEntitiesIndexMap_;
+    std::vector<TopTools_ListOfShape> entityIntersections_;
+
+    // sub-shape maps
+    TopTools_IndexedMapOfShape vmap_, emap_, wmap_, fmap_, shmap_, somap_;
+
+    // single compound describing the overall network
+    TopoDS_Compound compound_;
+
+    // index maps from primary to fragment indices
+    std::unordered_map<std::size_t, std::vector<std::size_t>> entityToFragmentsMap_;
+    std::unordered_map<std::size_t, std::vector<std::size_t>> entityIntersectionsToFragmentsMap_;
+    std::unordered_map<std::size_t, std::vector<std::size_t>> domainToFragmentsMap_;
+};
+
+} // end namespace Frackit
+
+#endif // FRACKIT_BREP_WRITER_HH

frackit/io/gmshwriter.hh

+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 2 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief Class that writes entity networks into Gmsh .geo file format.
+ *        Effectively, this will write a .brep file and a .geo file in
+ *        which the .brep file is included. Then, physical definitions
+ *        are added to all entity and sub-domain fragments.
+ */
+#ifndef FRACKIT_GMSH_WRITER_HH
+#define FRACKIT_GMSH_WRITER_HH
+
+#include <string>
+#include <fstream>
+
+#include "brepwriter.hh"
+
+namespace Frackit {
+
+/*!
+ * \brief Class that writes entity networks into Gmsh .geo file format.
+ *        Effectively, this will write a .brep file and a .geo file in
+ *        which the .brep file is included. Then, physical definitions
+ *        are added to all entity and sub-domain fragments.
+ */
+class GmshWriter : public BRepWriter
+{
+
+public:
+
+    /*!
+     * \brief Constructor from an entity network.
+     */
+    template<class EntityNetwork>
+    GmshWriter(const EntityNetwork& network)
+    : BRepWriter(network)
+    , entityDimension_(network.entityDimension())
+    {}
+
+    /*!
+     * \brief Write the .geo file.
+     * \param fileName The body of the file name to be used.
+     * \param meshSizeAtEntities Characteristic mesh size at entities.
+     * \param meshSizeAtBoundary Characteristic mesh size at domain boundaries.
+     * \note Two files are created: _fileName_.brep and _fileName_.geo
+     */
+    template<class ctype>
+    void write(const std::string& fileName,
+               ctype meshSizeAtEntities,
+               ctype meshSizeAtBoundary)
+    {
+        BRepWriter::write(fileName);
+
+        std::ofstream geoFile(fileName + ".geo");
+        geoFile << "Merge \"";
+        geoFile << fileName + ".brep";
+        geoFile << "\";\n";
+
+        if (!this->domainToFragmentsMap().empty())
+        {
+            geoFile << "\n// Physical domain definitions\n";
+            auto it = this->domainToFragmentsMap().begin();
+            for (; it != this->domainToFragmentsMap().end(); ++it)
+            {
+                geoFile << "Physical ";
+                if (entityDimension_ == 2) geoFile << "Volume(";
+                else if (entityDimension_ == 1) geoFile << "Surface(";
+                geoFile << it->first;
+                geoFile << ") = {";
+                geoFile << it->second[0];
+                for (unsigned int j = 1; j < it->second.size(); ++j)
+                    geoFile << ", " << it->second[j];
+                geoFile << "};\n";
+
+                // mesh size specification
+                geoFile << "Characteristic Length{ PointsOf{";
+                if (entityDimension_ == 2) geoFile << "Volume{";
+                else if (entityDimension_ == 1) geoFile << "Surface{";
+                geoFile << it->second[0];
+                for (unsigned int j = 1; j < it->second.size(); ++j)
+                    geoFile << ", " << it->second[j];
+                geoFile << "};} } = ";
+                geoFile << meshSizeAtBoundary << ";\n";
+            }
+        }
+
+        if (!this->entityToFragmentsMap().empty())
+        {
+            geoFile << "\n// Physical entity definitions\n";
+            auto it = this->entityToFragmentsMap().begin();
+            for (; it != this->entityToFragmentsMap().end(); ++it)
+            {
+                geoFile << "Physical ";
+                if (entityDimension_ == 2) geoFile << "Surface(";
+                else if (entityDimension_ == 1) geoFile << "Line(";
+                geoFile << it->first;
+                geoFile << ") = {";
+                geoFile << it->second[0];
+                for (unsigned int j = 1; j < it->second.size(); ++j)
+                    geoFile << ", " << it->second[j];
+                geoFile << "};\n";
+
+                // mesh size specification
+                geoFile << "Characteristic Length{ PointsOf{";
+                if (entityDimension_ == 2) geoFile << "Surface{";
+                else if (entityDimension_ == 1) geoFile << "Line{";
+                geoFile << it->second[0];
+                for (unsigned int j = 1; j < it->second.size(); ++j)
+                    geoFile << ", " << it->second[j];
+                geoFile << "};} } = ";
+                geoFile << meshSizeAtEntities << ";\n";
+            }
+        }
+
+        if (!this->entityIntersectionsToFragmentsMap().empty())
+        {
+            geoFile << "\n// Physical entity intersections definitions\n";
+            auto it = this->entityIntersectionsToFragmentsMap().begin();
+            for (; it != this->entityIntersectionsToFragmentsMap().end(); ++it)
+            {
+                geoFile << "Physical ";
+                if (entityDimension_ == 2) geoFile << "Line(";
+                else if (entityDimension_ == 1) geoFile << "Point(";
+                geoFile << it->first;
+                geoFile << ") = {";
+                geoFile << it->second[0];
+                for (unsigned int j = 1; j < it->second.size(); ++j)
+                    geoFile << ", " << it->second[j];
+                geoFile << "};\n";
+            }
+        }
+    }
+
+    /*!
+     * \brief Write the .geo file.
+     * \param fileName The body of the file name to be used.
+     * \param meshSizeAtEntities Characteristic mesh size at entities.
+     * \note This uses the provided mesh size on both entities and boundaries.
+     * \note Two files are created: _fileName_.brep and _fileName_.geo.
+     */
+    template<class ctype>
+    void write(const std::string& fileName, ctype meshSizeAtEntities)
+    { write(fileName, meshSizeAtEntities, meshSizeAtEntities); }
+
+private:
+    int entityDimension_;
+};
+
+} // end namespace Frackit
+
+#endif // FRACKIT_GMSH_WRITER_HH

test/CMakeLists.txt

@@ -2,5 +2,6 @@ add_subdirectory(distance)
 add_subdirectory(entitynetwork)
 add_subdirectory(geometry)
 add_subdirectory(intersection)
+add_subdirectory(io)
 add_subdirectory(magnitude)
 add_subdirectory(sampling)

test/io/CMakeLists.txt

+# test write of networks to gmsh output format
+frackit_add_test(NAME test_write_gmsh
+                 SOURCES test_write_gmsh.cc
+                 LABELS io gmsh)
+
+set(CMAKE_BUILD_TYPE Debug)

test/io/layer.brep

+DBRep_DrawableShape
+
+CASCADE Topology V1, (c) Matra-Datavision
+Locations 6
+1
+              1               0               0               0 
+              0               1               0             100 
+              0               0               1               0 
+1
+              1               0               0               0 
+              0               1               0               0 
+              0               0               1              25 
+2  1 1 2 1 0
+2  1 -1 0
+2  2 -1 0
+2  2 -1 1 -1 0
+Curve2ds 12
+1 0 0 1 0 
+1 0 -100 1 0 
+1 0 0 1 0 
+1 0 -25 1 0 
+1 0 0 1 0 
+1 0 -25 1 0 
+1 0 0 0 -1 
+1 1 0 0 -1 
+1 0 0 0 -1 
+1 0 0 0 -1 
+1 1 0 0 -1 
+1 1 0 0 -1 
+Curves 7
+6 0 3 0 0 0  50 0 25  100 0 0  150 0 0  
+1 0 0 0 0 1 0 
+1 150 0 0 0 1 0 
+1 0 0 0 0 0 1 
+1 0 100 0 0 0 1 
+1 150 0 0 0 0 1 
+1 150 100 0 0 0 1 
+Polygon3D 0
+PolygonOnTriangulations 0
+Surfaces 5
+6 -0 -1 -0 
+6 0 3 0 0 0  50 0 25  100 0 0  150 0 0  
+6 -0 -0 -1 
+6 0 3 0 0 0  50 0 25  100 0 0  150 0 0  
+6 -0 -0 -1 
+6 0 3 0 100 0  50 100 25  100 100 0  150 100 0  
+1 0 0 0 -1 0 0 0 1 0 0 0 -1 
+1 150 0 0 -1 0 0 0 1 0 0 0 -1 
+Triangulations 0
+
+TShapes 24