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

Feature/improve gmsh writer

Closes #17

See merge request tools/frackit!202

appl/example1/README.md

@@ -147,8 +147,8 @@ This network can then be written to disk, for example in [Gmsh][1] (.geo) file f
 
 ```cpp
 GmshWriter writer(network);
-writer.write("network", // filename of the .geo files (will add extension .geo automatically)
-             0.1);      // element size to be used on the network entities
+writer.setMeshSize(GmshWriter::GeometryTag::entity, 0.1); // mesh size at entities
+writer.write("network"); // filename of the .geo files (will add extension .geo automatically)
 ```
 
 Note that with the `EntityNetworkBuilder` class we have created a network that

appl/example1/example1.cc

@@ -133,8 +133,8 @@ int main(int argc, char** argv)
     // meshed by a two-dimensional surface mesh
     std::cout << "\n --- Writing .geo file ---\n" << std::endl;
     GmshWriter writer(network);
-    writer.write("network", // filename of the .geo files (will add extension .geo automatically)
-                 0.1);      // element size to be used
+    writer.setMeshSize(GmshWriter::GeometryTag::entity, 0.1); // mesh size at entities
+    writer.write("network"); // filename of the .geo files (will add extension .geo automatically)
     std::cout << "\n --- Finished writing .geo file ---\n" << std::endl;
 
     return 0;

appl/example1/example1.py

@@ -115,6 +115,7 @@ network = builder.build();
 print("\n --- Writing .geo file ---\n")
 from frackit.io import GmshWriter
 writer = GmshWriter(network);
-writer.write("network", # filename of the .geo files (will add extension .geo automatically)
-             0.1);      # element size to be used
+writer.setMeshSize(GmshWriter.GeometryTag.entity, 0.1)
+writer.write("network") # filename of the .geo files (will add extension .geo automatically)
+
 print("\n --- Finished writing .geo file ---\n")

appl/example2/README.md

@@ -143,9 +143,9 @@ which we can also pass to the `GmshWriter`:
 
 ```cpp
 GmshWriter writer(network);
-writer.write("network", // filename of the .geo files (will add extension .geo automatically)
-             0.1,       // element size to be used on the quadrilaterals
-             0.2);      // element size to be used in the domain away from the quadrilaterals
+writer.setMeshSize(GmshWriter::GeometryTag::entity, 0.1); // mesh size at entities
+writer.setMeshSize(GmshWriter::GeometryTag::subDomain, 0.2); // mesh size in the domain away from entities
+writer.write("network"); // filename of the .geo files (will add extension .geo automatically)
 ```
 
 As you can see, one can specify different mesh sizes to be used on the fracture

appl/example2/example2.cc

@@ -133,9 +133,9 @@ int main(int argc, char** argv)
     // to the two-dimensional quadrilaterals
     std::cout << "\n --- Writing .geo file ---\n" << std::endl;
     GmshWriter writer(network);
-    writer.write("network", // filename of the .geo files (will add extension .geo automatically)
-                 0.1,       // element size to be used on the quadrilaterals
-                 0.2);      // element size to be used in the domain away from the quadrilaterals
+    writer.setMeshSize(GmshWriter::GeometryTag::entity, 0.1); // mesh size at entities
+    writer.setMeshSize(GmshWriter::GeometryTag::subDomain, 0.2); // mesh size in the domain away from entities
+    writer.write("network"); // filename of the .geo files (will add extension .geo automatically)
     std::cout << "\n --- Finished writing .geo file ---\n" << std::endl;
 
     return 0;

appl/example2/example2.py

@@ -148,6 +148,6 @@ network = builder.build();
 print("\n --- Writing .geo file ---\n")
 from frackit.io import GmshWriter
 writer = GmshWriter(network);
-writer.write("network", # filename of the .geo files (will add extension .geo automatically)
-             0.1);      # element size to be used
+writer.setMeshSize(GmshWriter.GeometryTag.entity, 0.1)
+writer.write("network") # filename of the .geo files (will add extension .geo automatically)
 print("\n --- Finished writing .geo file ---\n")

appl/example3/example3.cc

@@ -235,9 +235,9 @@ int main(int argc, char** argv)
 
     // now we can build and write out the network in Gmsh file format
     GmshWriter gmshWriter(builder.build());
-    gmshWriter.write("contained_confined", // body of the filename to be used (will add .geo)
-                     2.5,                  // mesh size to be used on entities
-                     5.0);                 // mesh size to be used on domain boundaries
+    gmshWriter.setMeshSize(GmshWriter::GeometryTag::entity, 2.5);
+    gmshWriter.setMeshSize(GmshWriter::GeometryTag::subDomain, 5.0);
+    gmshWriter.write("contained_confined");
 
     // we can also not confine the network to its sub-domain,
     // simply by adding the sub-domains as non-confining
@@ -249,7 +249,9 @@ int main(int argc, char** argv)
     entitySets.exportEntitySets(builder, Id(2));
 
     gmshWriter = GmshWriter(builder.build());
-    gmshWriter.write("contained_unconfined", 2.5, 5.0);
+    gmshWriter.setMeshSize(GmshWriter::GeometryTag::entity, 2.5);
+    gmshWriter.setMeshSize(GmshWriter::GeometryTag::subDomain, 5.0);
+    gmshWriter.write("contained_unconfined");
 
     // We could also only write out the network, without the domain
     // For example, confining the network to the sub-domain...
@@ -258,14 +260,16 @@ int main(int argc, char** argv)
     uncontainedBuilder.addConfiningSubDomain(networkDomain, Id(2));
     entitySets.exportEntitySets(uncontainedBuilder, Id(2));
     gmshWriter = GmshWriter(uncontainedBuilder.build());
-    gmshWriter.write("uncontained_confined", 2.5);
+    gmshWriter.setMeshSize(GmshWriter::GeometryTag::entity, 2.5);
+    gmshWriter.write("uncontained_confined");
 
     // ... or not confining it
     std::cout << "Building and writing uncontained, unconfined network" << std::endl;
     uncontainedBuilder.clear();
     entitySets.exportEntitySets(uncontainedBuilder);
     gmshWriter = GmshWriter(uncontainedBuilder.build());
-    gmshWriter.write("uncontained_unconfined", 2.5);
+    gmshWriter.setMeshSize(GmshWriter::GeometryTag::entity, 2.5);
+    gmshWriter.write("uncontained_unconfined");
 
     return 0;
 }

appl/example3/example3.py

@@ -205,9 +205,9 @@ for setId in entitySets: builder.addSubDomainEntities(entitySets[setId], Id(2))
 # now we can build and write out the network in Gmsh file format
 from frackit.io import GmshWriter
 gmshWriter = GmshWriter(builder.build());
-gmshWriter.write("contained_confined", # body of the filename to be used (will add .geo)
-                 2.5,                  # mesh size to be used on entities
-                 5.0);                 # mesh size to be used on domain boundaries
+gmshWriter.setMeshSize(GmshWriter.GeometryTag.entity, 2.5)
+gmshWriter.setMeshSize(GmshWriter.GeometryTag.subDomain, 5.0)
+gmshWriter.write("contained_confined") # body of the filename to be used (will add .geo)
 
 # we can also not confine the network to its sub-domain,
 # simply by adding the sub-domains as non-confining
@@ -219,7 +219,9 @@ builder.addSubDomain(solids[2],     Id(3));
 for setId in entitySets: builder.addSubDomainEntities(entitySets[setId], Id(2))
 
 gmshWriter = GmshWriter(builder.build());
-gmshWriter.write("contained_unconfined", 2.5, 5.0);
+gmshWriter.setMeshSize(GmshWriter.GeometryTag.entity, 2.5)
+gmshWriter.setMeshSize(GmshWriter.GeometryTag.subDomain, 5.0)
+gmshWriter.write("contained_unconfined");
 
 # We could also only write out the network, without the domain
 # For example, confining the network to the sub-domain...
@@ -229,7 +231,8 @@ uncontainedBuilder.addConfiningSubDomain(networkDomain, Id(2));
 for setId in entitySets: uncontainedBuilder.addSubDomainEntities(entitySets[setId], Id(2))
 
 gmshWriter = GmshWriter(uncontainedBuilder.build());
-gmshWriter.write("uncontained_confined", 2.5);
+gmshWriter.setMeshSize(GmshWriter.GeometryTag.entity, 2.5)
+gmshWriter.write("uncontained_confined");
 
 # ... or not confining it
 print("Building and writing uncontained, unconfined network")
@@ -237,4 +240,5 @@ uncontainedBuilder.clear();
 for setId in entitySets: uncontainedBuilder.addSubDomainEntities(entitySets[setId], Id(2))
 
 gmshWriter = GmshWriter(uncontainedBuilder.build());
-gmshWriter.write("uncontained_unconfined", 2.5);
+gmshWriter.setMeshSize(GmshWriter.GeometryTag.entity, 2.5)
+gmshWriter.write("uncontained_unconfined");

appl/example4/README.md

@@ -262,6 +262,7 @@ Gmsh file format:
 
 ```cpp
 // now we can build and write out the network in Gmsh file format
-GmshWriter gmshWriter(builder.build());
-gmshWriter.write("network", /*sizeAtEntities*/0.5, /*sizeInDomain*/5.0);
+gmshWriter.setMeshSize(GmshWriter::GeometryTag::entity, 0.5);
+gmshWriter.setMeshSize(GmshWriter::GeometryTag::subDomain, 5.0);
+gmshWriter.write("network");
 ```

appl/example4/example4.cc

@@ -439,7 +439,9 @@ int main(int argc, char** argv)
 
     // now we can build and write out the network in Gmsh file format
     GmshWriter gmshWriter(builder.build());
-    gmshWriter.write("network", /*sizeAtEntities*/0.5, /*sizeInDomain*/5.0);
+    gmshWriter.setMeshSize(GmshWriter::GeometryTag::entity, 0.5);
+    gmshWriter.setMeshSize(GmshWriter::GeometryTag::subDomain, 5.0);
+    gmshWriter.write("network");
 
     // print time it took to generate the network
     const auto end = std::chrono::steady_clock::now();

appl/example4/example4.py

@@ -363,7 +363,9 @@ for id in entitySets: builder.addSubDomainEntities(entitySets[id], Id(1))
 # now we can build and write out the network in Gmsh file format
 from frackit.io import GmshWriter
 gmshWriter = GmshWriter(builder.build());
-gmshWriter.write("network", 0.5, 5.0)
+gmshWriter.setMeshSize(GmshWriter.GeometryTag.entity, 0.5)
+gmshWriter.setMeshSize(GmshWriter.GeometryTag.subDomain, 5.0)
+gmshWriter.write("network")
 
 stopTime = process_time()
 print("Overall CPU time was {:.2f} seconds.".format(stopTime-startTime))

frackit/io/CMakeLists.txt

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

frackit/io/brepwriter.hh

@@ -76,7 +76,7 @@ public:
     /*!
      * \brief Write the entity network to disk.
      */
-    void write(const std::string& fileName)
+    void write(const std::string& fileName) const
     {
         const std::string brepFile = fileName + ".brep";
         BRepTools::Write(compound_, brepFile.c_str());

frackit/io/gmshbackend.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 3 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
+ * \ingroup IO
+ * \brief Backend for Gmsh Syntax, depending on the desired .geo file format version.
+ */
+#ifndef FRACKIT_GMSH_FORMAT_BACKEND_HH
+#define FRACKIT_GMSH_FORMAT_BACKEND_HH
+
+#include <array>
+#include <string>
+#include <limits>
+#include <cassert>
+#include <type_traits>
+
+namespace Frackit::Gmsh {
+
+/*!
+ * \ingroup IO
+ * \brief Backend for writing Gmsh Syntax to files,
+ *        depending on the desired .geo file format version.
+ */
+class GmshFormatBackend
+{
+
+public:
+
+    /*!
+     * \brief Default constructor.
+     */
+    GmshFormatBackend()
+    {
+        setVersion(4, 0);
+    }
+
+    /*!
+     * \brief Constructor defining the gmsh version
+     *        for which a .geo file is to be written.
+     * \param majVersion major gmsh version
+     * \param minVersion minor gmsh version
+     */
+    GmshFormatBackend(int maj, int min)
+    {
+        setVersion(maj, min);
+    }
+
+    /*!
+     * \brief Update version settings.
+     */
+    void setVersion(int maj, int min)
+    {
+        // we use two-digits minor version
+        if (min < 10) min *= 10;
+
+        assert(min < 100);
+        gmshVersion_ = maj*100 + min;
+        assert(gmshVersion_ < 1000);
+    }
+
+    /*!
+     * \brief Write a physical entity definition.
+     * \param geoFile The geoFile to write to
+     * \param entityDimension Entity dimension
+     * \param physicalId The id to give the physical entity
+     * \param elementaryIds Ids of the elementary entities composing the physical entity
+     */
+    template<class IdList>
+    void writePhysicalGeometry(std::ofstream& geoFile,
+                               int entityDimension,
+                               std::size_t physicalId,
+                               const IdList& elementaryIds) const
+    {
+        static_assert(std::is_integral_v<typename IdList::value_type>,
+                      "Ids must be integral types");
+
+        const auto& geomKey = getGeometryKeyword_(entityDimension);
+        geoFile << "Physical " << geomKey << "(" << physicalId << ") = ";
+        writeIdList_(geoFile, elementaryIds);
+        geoFile << ";\n";
+    }
+
+    /*!
+     * \brief Write variable definition.
+     * \param geoFile The geoFile to write to
+     * \param varName variable name
+     * \param varValue variable value
+     */
+    template<class ValueType>
+    void writeVariableDefinition(std::ofstream& geoFile,
+                                 const std::string& varName,
+                                 const ValueType& varValue) const
+    {
+        geoFile << varName << " = " << varValue << ";\n";
+    }
+
+    /*!
+     * \brief Write mesh size definition.
+     * \param geoFile The geoFile to write to
+     * \param entityDimension Dimension of the entity on which to set the mesh size.
+     * \param elementaryIds List of elementary entity ids (must correspond
+     *                      to entities with dimension entityDimension)
+     * \param size string containing the size or size variable name.
+     */
+    template<class IdList>
+    void writeMeshSizeDefinition(std::ofstream& geoFile,
+                                 int entityDimension,
+                                 const IdList& idList,
+                                 const std::string& size) const
+    {
+        const auto& meshSizeKey = getMeshSizeKeyword_();
+        const auto& geomKey = getGeometryKeyword_(entityDimension);
+        geoFile << meshSizeKey << "{ PointsOf{" << geomKey;
+        writeIdList_(geoFile, idList);
+        geoFile << ";} } = " << size << ";\n";
+    }
+
+private:
+    template<class IdList>
+    void writeIdList_(std::ofstream& geoFile, const IdList& idList) const
+    {
+        geoFile << "{";
+        geoFile << idList[0];
+        for (unsigned int j = 1; j < idList.size(); ++j)
+            geoFile << ", " << idList[j];
+        geoFile << "}";
+    }
+
+    std::string getGeometryKeyword_(unsigned int geomDim) const
+    {
+        assert(geomDim <= 3);
+        static const std::array<std::string, 4> keywords = {
+            "Point", "Curve", "Surface", "Volume"
+        };
+
+        return keywords[geomDim];
+    }
+
+    std::string getMeshSizeKeyword_() const
+    {
+        if (gmshVersion_ < 470)
+            return "Characteristic Length";
+        return "MeshSize";
+    }
+
+    int gmshVersion_;
+};
+
+namespace Detail {
+
+class MeshSizeVariableNameHelper
+{
+    static constexpr double maxDouble = std::numeric_limits<double>::max();
+    std::array<std::size_t, 4> counts = {0, 0, 0, 0};
+    std::array<double, 4> sizes = {maxDouble, maxDouble, maxDouble, maxDouble};
+    std::array<std::string, 4> names = {"subDomainMeshSize_",
+                                        "entityMeshSize_",
+                                        "intersectionMeshSize_",
+                                        "pointMeshSize_"};
+
+public:
+    void reset()
+    {
+        for (auto& c : counts) c = 0;
+        for (auto& s : sizes) s = maxDouble;
+    }
+
+    bool prepareNew(unsigned int codim, double size)
+    {
+        if (sizes[codim] > size)
+        {
+            sizes[codim] = size;
+            counts[codim]++;
+
+            const auto pos = names[codim].find_last_of("_");
+            const auto suffix = std::to_string(counts[codim]);
+            const auto suffixLength = suffix.size();
+            names[codim].resize(pos+1 + suffixLength);
+            names[codim].replace(pos+1, suffixLength, suffix);
+
+            return true;
+        }
+
+        return false;
+    }
+
+    const std::string& get(unsigned int codim)
+    { return names[codim]; }
+};
+
+} // end namespace Detail
+
+} // end namespace Frackit::Gmsh
+
+#endif // FRACKIT_GMSH_BACKEND_HH

frackit/python/io/gmshwriter.hh

@@ -19,9 +19,12 @@
 #ifndef FRACKIT_PYTHON_IO_GMSH_WRITER_HH
 #define FRACKIT_PYTHON_IO_GMSH_WRITER_HH
 
+#include <vector>
 #include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
 
 #include <frackit/io/gmshwriter.hh>
+#include <frackit/common/id.hh>
 #include <frackit/entitynetwork/entitynetwork.hh>
 #include <frackit/python/entitynetwork/traits.hh>
 
@@ -36,8 +39,70 @@ void registerGmshWriter(py::module& module)
 
     py::class_<GmshWriter> cls(module, "GmshWriter");
     cls.def(py::init<const EntityNetwork&>());
+    cls.def(py::init<const EntityNetwork&, int, int>());
 
     using namespace py::literals;
+
+    // overloading with py::overload_cast seemed to fail when having mixed template/non-template functions
+    cls.def("write",
+            static_cast<void (GmshWriter::*)(const std::string&) const>(&GmshWriter::write),
+            "fileName"_a, "write the entity network to a .geo file");
+
+    // register geometry tags for setter functions
+    py::enum_<GmshWriter::GeometryTag>(cls, "GeometryTag")
+    .value("subDomain", GmshWriter::GeometryTag::subDomain)
+    .value("entity", GmshWriter::GeometryTag::entity)
+    .value("intersection", GmshWriter::GeometryTag::intersection)
+    .value("junction", GmshWriter::GeometryTag::junction);
+
+    // register geometry tag getter function for given dimension
+    cls.def("getGeometryTag", &GmshWriter::getGeometryTag,
+            "returns the geometry tag associated with the geometries of the given dimension");
+
+    // physical status setter functions
+    cls.def("setGmshVersion", &GmshWriter::setGmshVersion, "set gmsh format version");
+
+    // we define overloads for vector<Id> which is automatically converted
+    // by pybind11 such that it can be used with lists from within Python
+    using IdList = std::vector<Id>;
+
+    cls.def("setNonPhysical",
+            static_cast<void (GmshWriter::*)(GmshWriter::GeometryTag)>(&GmshWriter::setNonPhysical),
+            "set all geometries with the given type tag to be non-physical");
+    cls.def("setNonPhysical",
+            static_cast<void (GmshWriter::*)(GmshWriter::GeometryTag, const Id&)>(&GmshWriter::setNonPhysical),
+            "set the geometry with the given type tag and identifier to be non-physical");
+    cls.def("setNonPhysical",
+            py::overload_cast<GmshWriter::GeometryTag, const IdList&>(&GmshWriter::template setNonPhysical<IdList>),
+            "set all geometries with a specific type tag and with the ids of the given list to be non-physical");
+    cls.def("setPhysical",
+            static_cast<void (GmshWriter::*)(GmshWriter::GeometryTag)>(&GmshWriter::setPhysical),
+            "set all geometries with the given type tag to be physical");
+    cls.def("setPhysical",
+            static_cast<void (GmshWriter::*)(GmshWriter::GeometryTag, const Id&)>(&GmshWriter::setPhysical),
+            "set the geometry with the given type tag and identifier to be physical");
+    cls.def("setPhysical",
+            py::overload_cast<GmshWriter::GeometryTag, const IdList&>(&GmshWriter::template setPhysical<IdList>),
+            "set all geometries with a specific type tag and with the ids of the given list to be physical");
+
+    // mesh size setters
+    cls.def("setMeshSize",
+            static_cast<void (GmshWriter::*)(GmshWriter::GeometryTag, ctype)>(&GmshWriter::template setMeshSize<ctype>),
+            "set the mesh size for all geometries with the given type tag");
+    cls.def("setMeshSize",
+            static_cast<void (GmshWriter::*)(GmshWriter::GeometryTag, const Id&, ctype)>(&GmshWriter::template setMeshSize<ctype>),
+            "set the mesh size for the geometry with the given type tag and identifier");
+    cls.def("setMeshSize",
+            py::overload_cast<GmshWriter::GeometryTag, const IdList&, ctype>(&GmshWriter::template setMeshSize<IdList, ctype>),
+            "set the mesh size for all geometries with a specific type tag and with the ids of the given list");
+    cls.def("resetMeshSizes",
+            py::overload_cast<>(&GmshWriter::resetMeshSizes),
+            "Reset all mesh size settings");
+    cls.def("resetMeshSizes",
+            py::overload_cast<GmshWriter::GeometryTag>(&GmshWriter::resetMeshSizes),
+            "Reset the mesh size settings for geometries with the given type tag");
+
+    // deprecated (TODO: remove afer 1.2 release)
     cls.def("write",
             py::overload_cast<const std::string&, ctype, ctype>(&GmshWriter::template write<ctype>),
             "fileName"_a, "meshSizeAtEntities"_a, "meshSizeAtBoundary"_a,

test/io/test_write_gmsh.cc

 #include <iostream>
-#include <random>
 #include <fstream>
-#include <unordered_map>
-
-#include <BRepTools.hxx>
-#include <BRep_Builder.hxx>
-#include <TopTools_IndexedMapOfShape.hxx>
-#include <TopTools_DataMapOfShapeInteger.hxx>
-#include <TopoDS_Compound.hxx>
-#include <TopoDS_Solid.hxx>
-#include <TopoDS.hxx>
 
+#include <frackit/common/id.hh>
 #include <frackit/geometry/cylinder.hh>
 #include <frackit/geometry/disk.hh>
 
-#include <frackit/magnitude/magnitude.hh>
-#include <frackit/magnitude/containedmagnitude.hh>
+#include <frackit/entitynetwork/networkbuilder.hh>
+#include <frackit/io/gmshwriter.hh>
 
-#include <frackit/common/id.hh>
-#include <frackit/occ/breputilities.hh>
-#include <frackit/precision/precision.hh>
+void testGeoFile(const std::string& fileName,
+                 const std::string& testString,
+                 std::size_t numOccurrences,
+                 bool atBeginning = false)
+{
+    std::ifstream geoFile(fileName + ".geo");
+    if (geoFile.fail()) throw std::runtime_error("Could not open geo file");
+
+    std::size_t foundCount = 0;
+    while (!geoFile.eof())
+    {
+        std::string line;
+        std::getline(geoFile, line);
+        const auto pos = line.find(testString);
+        if (pos != std::string::npos)
+            if (!atBeginning || pos == 0)
+                foundCount++;
+    }
+
+    if (foundCount != numOccurrences)
+    {
+        std::cout << "Found " << testString << " " << foundCount
+                  << " times instead of " << numOccurrences << std::endl;
+        throw std::runtime_error("Could not find given number of occurrences");
+    }
+
+    std::cout << "Successfully found " << testString << " " << foundCount << " times." << std::endl;
+}
 
-#include <frackit/sampling/makeuniformpointsampler.hh>
-#include <frackit/sampling/disksampler.hh>
 
-#include <frackit/entitynetwork/networkbuilder.hh>
-#include <frackit/io/gmshwriter.hh>
+template<class Call>
+void testException(const Call& call)
+{
+    bool passed = false;
+    try { call(); }
+    catch (...) { passed = true; }
+
+    if (!passed) throw std::runtime_error("No exception caught!");
+}
 
 //! create 3 disk-shaped fractures and write geo (gmsh file format) file
 int main(int argc, char** argv)
@@ -46,31 +67,148 @@ int main(int argc, char** argv)
     Domain domain(0.5, 1.0);
     Domain domain2(Disk(Point(0.0, 0.0, 1.0), e1, e2, 1.0, 1.0), 1.0);
 
-    // make and write non-contained network
-    EntityNetworkBuilder<ctype> builder;
-    builder.addEntity(Disk(Point(0.0, 0.0, 0.1), e1, e2, 1.0, 1.0));
-    builder.addEntity(Disk(Point(0.0, 0.0, 0.5), e1, e3, 2.0, 2.0));
-    builder.addEntity(Disk(Point(0.0, 0.0, 0.75), e1, e2, 1.0, 1.0));
-    builder.addEntity(Disk(Point(0.0, 0.0, 1.1), e1, e2, 1.0, 1.0));
-    builder.addEntity(Disk(Point(0.0, 0.0, 1.5), e1, e3, 2.0, 2.0));
-    builder.addEntity(Disk(Point(0.0, 0.0, 1.75), e1, e2, 1.0, 1.0));
+    // disks to be added
+    std::vector<Disk> disks; disks.reserve(7);
+    disks.emplace_back(Disk(Point(0.0, 0.0, 0.1), e1, e2, 1.0, 1.0));
+    disks.emplace_back(Disk(Point(0.0, 0.1, 0.5), e1, e3, 2.0, 2.0));
+    disks.emplace_back(Disk(Point(0.0, 0.0, 0.75), e1, e2, 1.0, 1.0));
+    disks.emplace_back(Disk(Point(0.0, 0.0, 0.5), e2, e3, 1.0, 1.0));