diff --git a/Applications/FileIO/CMakeLists.txt b/Applications/FileIO/CMakeLists.txt
index c83d1009486b8296e85b58f1cc4f84436a5a06c1..91b4dfc6a13b055d0c8f5940b2248b0699f2175e 100644
--- a/Applications/FileIO/CMakeLists.txt
+++ b/Applications/FileIO/CMakeLists.txt
@@ -6,6 +6,10 @@ if(NOT Shapelib_FOUND)
list(REMOVE_ITEM SOURCES SHPInterface.h SHPInterface.cpp)
endif()
+# GO2OGS
+GET_SOURCE_FILES(SOURCES_GO2OGS GocadIO)
+set(SOURCES ${SOURCES} ${SOURCES_GO2OGS})
+
if(Qt5XmlPatterns_FOUND)
APPEND_SOURCE_FILES(SOURCES XmlIO/Qt)
APPEND_SOURCE_FILES(SOURCES FEFLOW)
diff --git a/Applications/FileIO/GocadIO/CMakeLists.txt b/Applications/FileIO/GocadIO/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..7ca90dceddc1db1c45a524b9f51cf6db2fc3ea36
--- /dev/null
+++ b/Applications/FileIO/GocadIO/CMakeLists.txt
@@ -0,0 +1,6 @@
+include_directories(
+ ${CMAKE_SOURCE_DIR}/BaseLib
+ ${CMAKE_SOURCE_DIR}/FileIO
+ ${CMAKE_SOURCE_DIR}/GeoLib
+ ${CMAKE_SOURCE_DIR}/MeshLib
+)
diff --git a/Applications/FileIO/GocadIO/GenerateFaceSetMeshes.cpp b/Applications/FileIO/GocadIO/GenerateFaceSetMeshes.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d0c068710aef03ea084660a6b565a98a340c6c72
--- /dev/null
+++ b/Applications/FileIO/GocadIO/GenerateFaceSetMeshes.cpp
@@ -0,0 +1,40 @@
+/**
+ *
+ * @copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/LICENSE.txt
+ */
+
+#include "GenerateFaceSetMeshes.h"
+#include "MeshLib/Elements/Quad.h"
+#include "MeshLib/IO/writeMeshToFile.h"
+
+namespace FileIO
+{
+namespace Gocad
+{
+
+void generateFaceSets(GocadSGridReader const& reader, std::string const& path)
+{
+ for (std::size_t l(0); l < 128; l++)
+ {
+ std::unique_ptr<MeshLib::Mesh> face_set(reader.getFaceSetMesh(l));
+
+ if (!face_set)
+ {
+ continue;
+ }
+ INFO("Face set mesh created. #nodes: %d, #elements: %d",
+ face_set->getNumberOfNodes(),
+ face_set->getNumberOfElements());
+
+ std::string const mesh_out_fname(path + face_set->getName() + ".vtu");
+ INFO("Writing face set mesh to '%s'.", mesh_out_fname.c_str());
+ MeshLib::IO::writeMeshToFile(*face_set, mesh_out_fname);
+ }
+}
+
+} // namespace Gocad
+} // namespace FileIO
diff --git a/Applications/FileIO/GocadIO/GenerateFaceSetMeshes.h b/Applications/FileIO/GocadIO/GenerateFaceSetMeshes.h
new file mode 100644
index 0000000000000000000000000000000000000000..ec197c85661e31a7a8d0158c351341ac783664cf
--- /dev/null
+++ b/Applications/FileIO/GocadIO/GenerateFaceSetMeshes.h
@@ -0,0 +1,27 @@
+
+/**
+ *
+ * @copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/LICENSE.txt
+ */
+
+#pragma once
+
+#include <string>
+
+#include "MeshLib/Mesh.h"
+
+#include "GocadSGridReader.h"
+
+namespace FileIO
+{
+namespace Gocad
+{
+
+void generateFaceSets(GocadSGridReader const& reader, std::string const& path);
+
+} // namespace Gocad
+} // namespace FileIO
diff --git a/Applications/FileIO/GocadIO/GocadSGridReader.cpp b/Applications/FileIO/GocadIO/GocadSGridReader.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..008db648a3274d5b54cbcef325d856a3a98773e3
--- /dev/null
+++ b/Applications/FileIO/GocadIO/GocadSGridReader.cpp
@@ -0,0 +1,856 @@
+/**
+ *
+ * @copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/LICENSE.txt
+ */
+
+#include "GocadSGridReader.h"
+
+#include <algorithm>
+#include <cstddef>
+#include <fstream>
+#include <iostream>
+#include <sstream>
+
+#include <boost/tokenizer.hpp>
+
+#include "MeshLib/Elements/Hex.h"
+#include "MeshLib/Elements/Quad.h"
+#include "MeshLib/Mesh.h"
+
+namespace FileIO
+{
+namespace Gocad
+{
+using Bitset = boost::dynamic_bitset<>;
+
+GocadSGridReader::GocadSGridReader(std::string const& fname)
+ : _fname(fname),
+ _path(BaseLib::extractPath(fname)),
+ _n_face_sets(0),
+ _double_precision_binary(false),
+ _bin_pnts_in_double_precision(false)
+{
+ // check if file exists
+ std::ifstream in(_fname.c_str());
+ if (!in)
+ {
+ ERR("Could not open '%s'.", _fname.c_str());
+ in.close();
+ return;
+ }
+
+ bool pnts_read(false);
+
+ // read information in the stratigraphic grid file
+ std::string line;
+ while (std::getline(in, line))
+ {
+ if (line.empty()) // skip empty lines
+ {
+ continue;
+ }
+ if (line.back() == '\r') // check dos line ending
+ {
+ line.pop_back();
+ }
+ if (line.substr(0, 8) == "HEADER {")
+ {
+ parseHeader(in);
+ }
+ if (line.substr(0, 32) == "GOCAD_ORIGINAL_COORDINATE_SYSTEM")
+ {
+ _coordinate_system.parse(in);
+ continue;
+ }
+ if (line.substr(0, 7) == "AXIS_N ")
+ {
+ parseDims(line);
+ }
+ else if (line.substr(0, 12) == "POINTS_FILE ")
+ {
+ parseFileName(line, _pnts_fname);
+ }
+ else if (line.substr(0, 9) == "PROPERTY ")
+ {
+ _property_meta_data_vecs.push_back(
+ Gocad::parseGocadPropertyMetaData(line, in, _path));
+ }
+ else if (line.substr(0, 35) == "BINARY_POINTS_IN_DOUBLE_PRECISION 1")
+ {
+ _bin_pnts_in_double_precision = true;
+ }
+ else if (line.substr(0, 11) == "FLAGS_FILE ")
+ {
+ parseFileName(line, _flags_fname);
+ }
+ else if (line.substr(0, 18) == "REGION_FLAGS_FILE ")
+ {
+ parseFileName(line, _region_flags_fname);
+ }
+ else if (line.substr(0, 7) == "REGION " ||
+ line.substr(0, 13) == "MODEL_REGION ")
+ {
+ regions.push_back(Gocad::parseRegion(line));
+ }
+ else if (line.substr(0, 12) == "MODEL_LAYER ")
+ {
+ layers.push_back(Gocad::parseLayer(line, regions));
+ }
+ else if (line.substr(0, 24) == "REGION_FLAGS_BIT_LENGTH ")
+ {
+ std::istringstream iss(line);
+ std::istream_iterator<std::string> it(iss);
+ it++;
+ std::size_t bit_length = std::atoi(it->c_str());
+ if (regions.size() != bit_length)
+ {
+ ERR("%d regions read but %d expected.\n", regions.size(),
+ bit_length);
+ throw std::runtime_error(
+ "Number of read regions differs from expected.\n");
+ }
+ }
+ else if (line.substr(0, 9) == "FACE_SET ")
+ {
+ // first read the points
+ if (!pnts_read)
+ {
+ readNodesBinary();
+ pnts_read = true;
+ }
+ parseFaceSet(line, in);
+ }
+ }
+
+#ifndef NDEBUG
+ std::stringstream regions_ss;
+ regions_ss << regions.size() << " regions read:\n";
+ std::copy(regions.begin(), regions.end(),
+ std::ostream_iterator<Gocad::Region>(regions_ss, "\t"));
+ DBUG("%s", regions_ss.str().c_str());
+
+ std::stringstream layers_ss;
+ layers_ss << layers.size() << " layers read:\n";
+ std::copy(layers.begin(), layers.end(),
+ std::ostream_iterator<Gocad::Layer>(layers_ss, "\n"));
+ DBUG("%s", layers_ss.str().c_str());
+
+ std::stringstream properties_ss;
+ properties_ss << "meta data for " << _property_meta_data_vecs.size()
+ << " properties read:\n";
+ std::copy(_property_meta_data_vecs.begin(), _property_meta_data_vecs.end(),
+ std::ostream_iterator<Gocad::Property>(properties_ss, "\n"));
+ DBUG("%s", properties_ss.str().c_str());
+#endif
+
+ // if not done already read the points
+ if (!pnts_read)
+ {
+ readNodesBinary();
+ pnts_read = true;
+ }
+ readElementPropertiesBinary();
+ std::vector<Bitset> region_flags = readRegionFlagsBinary();
+ mapRegionFlagsToCellProperties(region_flags);
+
+ readSplitInformation();
+
+ in.close();
+}
+
+GocadSGridReader::~GocadSGridReader()
+{
+ for (auto node : _nodes)
+ {
+ delete node;
+ }
+ for (auto node : _split_nodes)
+ {
+ delete node;
+ }
+}
+
+std::unique_ptr<MeshLib::Mesh> GocadSGridReader::getMesh() const
+{
+ std::vector<MeshLib::Node*> nodes;
+ std::transform(_nodes.cbegin(), _nodes.cend(), std::back_inserter(nodes),
+ [](MeshLib::Node const* const node) {
+ return new MeshLib::Node(*node);
+ });
+
+ std::vector<MeshLib::Element*> elements(createElements(nodes));
+ applySplitInformation(nodes, elements);
+
+ DBUG("Creating mesh from Gocad SGrid.");
+ std::unique_ptr<MeshLib::Mesh> mesh(new MeshLib::Mesh(
+ BaseLib::extractBaseNameWithoutExtension(_fname), nodes, elements));
+ addGocadPropertiesToMesh(*mesh);
+ DBUG("Mesh created.");
+
+ return mesh;
+}
+
+void GocadSGridReader::addGocadPropertiesToMesh(MeshLib::Mesh& mesh) const
+{
+ std::vector<std::string> const& prop_names(getPropertyNames());
+ for (auto const& name : prop_names)
+ {
+ boost::optional<Gocad::Property const&> prop(getProperty(name));
+ if (!prop)
+ {
+ continue;
+ }
+
+ DBUG("Adding Gocad property '%s' with %d values.", name.c_str(),
+ prop->_property_data.size());
+
+ auto pv = MeshLib::getOrCreateMeshProperty<double>(
+ mesh, name, MeshLib::MeshItemType::Cell, 1);
+ if (pv == nullptr)
+ {
+ ERR("Could not create mesh property '%s'.", name.c_str());
+ continue;
+ }
+
+ pv->resize(prop->_property_data.size());
+ std::copy(prop->_property_data.cbegin(), prop->_property_data.cend(),
+ pv->begin());
+ }
+}
+
+void GocadSGridReader::parseHeader(std::istream& in)
+{
+ std::string line;
+ while (std::getline(in, line))
+ {
+ if (line.front() == '}')
+ {
+ return;
+ }
+ if (line.substr(0, 27) == "double_precision_binary: on")
+ {
+ _double_precision_binary = true;
+ }
+ }
+ if (_double_precision_binary)
+ {
+ DBUG(
+ "GocadSGridReader::parseHeader(): _double_precision_binary == "
+ "true.");
+ }
+}
+
+void GocadSGridReader::parseDims(std::string const& line)
+{
+ std::size_t x_dim(0), y_dim(0), z_dim(0);
+ boost::tokenizer<> tok(line);
+ auto it(tok.begin());
+ it++; // overread token "AXIS"
+ it++; // overread "N"
+ std::stringstream ssx(*(it),
+ std::stringstream::in | std::stringstream::out);
+ ssx >> x_dim;
+ it++;
+ std::stringstream ssy(*it, std::stringstream::in | std::stringstream::out);
+ ssy >> y_dim;
+ it++;
+ std::stringstream ssz(*it, std::stringstream::in | std::stringstream::out);
+ ssz >> z_dim;
+ _index_calculator = Gocad::IndexCalculator(x_dim, y_dim, z_dim);
+ DBUG("x_dim = %d, y_dim = %d, z_dim = %d => #nodes = %d, #cells = %d",
+ x_dim, y_dim, z_dim, _index_calculator._n_nodes,
+ _index_calculator._n_cells);
+}
+
+void GocadSGridReader::parseFileName(std::string const& line,
+ std::string& result_string) const
+{
+ boost::char_separator<char> sep(" \r");
+ boost::tokenizer<boost::char_separator<char>> tok(line, sep);
+ auto it(tok.begin());
+ ++it; // overread POINTS_FILE or FLAGS_FILE or REGION_FLAGS_FILE
+ result_string = _path + *it;
+}
+
+/**
+ * @param line input/output
+ * @param in input stream containing the face set
+ */
+void GocadSGridReader::parseFaceSet(std::string& line, std::istream& in)
+{
+ // create and initialize a Gocad::Property object for storing face set data
+ Gocad::Property face_set_property;
+ face_set_property._property_id = _n_face_sets;
+ face_set_property._property_name = "FaceSet";
+ face_set_property._property_class_name = "FaceSetData";
+ face_set_property._property_unit = "unitless";
+ face_set_property._property_data_type = "double";
+ face_set_property._property_data_fname = "";
+ face_set_property._property_no_data_value = -1.0;
+ face_set_property._property_data.resize(_index_calculator._n_cells);
+ std::fill(face_set_property._property_data.begin(),
+ face_set_property._property_data.end(),
+ face_set_property._property_no_data_value);
+
+ std::istringstream iss(line);
+ std::istream_iterator<std::string> it(iss);
+ // Check first word is FACE_SET
+ if (*it != std::string("FACE_SET"))
+ {
+ ERR("Expected FACE_SET keyword but '%s' found.", it->c_str());
+ throw std::runtime_error(
+ "In GocadSGridReader::parseFaceSet() expected FACE_SET keyword not "
+ "found.");
+ }
+ ++it;
+ face_set_property._property_name += *it;
+ ++it;
+ std::size_t const n_of_face_set_ids(
+ static_cast<std::size_t>(std::atoi(it->c_str())));
+ std::size_t face_set_id_cnt(0);
+
+ while (std::getline(in, line) && face_set_id_cnt < n_of_face_set_ids)
+ {
+ if (line.back() == '\r')
+ {
+ line.pop_back();
+ }
+ boost::char_separator<char> sep("\t ");
+ boost::tokenizer<boost::char_separator<char>> tokens(line, sep);
+
+ for (auto tok_it = tokens.begin(); tok_it != tokens.end();)
+ {
+ std::size_t const id(
+ static_cast<std::size_t>(std::atoi(tok_it->c_str())));
+ tok_it++;
+ std::size_t const face_indicator(
+ static_cast<std::size_t>(std::atoi(tok_it->c_str())));
+ tok_it++;
+
+ if (id >= _index_calculator._n_nodes)
+ {
+ ERR("Face set id %d is greater than the number of nodes (%d).",
+ id, _index_calculator._n_nodes);
+ }
+ else
+ {
+ static_cast<GocadNode*>(_nodes[id])
+ ->setFaceSet(_n_face_sets, face_indicator);
+ std::array<std::size_t, 3> const c(
+ _index_calculator.getCoordsForID(id));
+ if (c[0] >= _index_calculator._x_dim - 1)
+ ERR("****** i coord %d to big for id %d.", c[0], id);
+ if (c[1] >= _index_calculator._y_dim - 1)
+ ERR("****** j coord %d to big for id %d.", c[1], id);
+ if (c[2] >= _index_calculator._z_dim - 1)
+ ERR("****** k coord %d to big for id %d.", c[2], id);
+ std::size_t const cell_id(
+ _index_calculator.getCellIdx(c[0], c[1], c[2]));
+ face_set_property._property_data[cell_id] = face_indicator;
+ }
+ face_set_id_cnt++;
+ }
+ }
+
+ if (face_set_id_cnt != n_of_face_set_ids)
+ {
+ ERR("Expected %d number of face set ids, read %d.", n_of_face_set_ids,
+ face_set_id_cnt);
+ throw std::runtime_error(
+ "Expected number of face set points does not match number of read "
+ "points.");
+ }
+ _n_face_sets++;
+
+ // pre condition: split nodes are read already
+ for (auto split_node : _split_nodes)
+ {
+ std::size_t const id(_index_calculator(split_node->getGridCoords()));
+ split_node->transmitFaceDirections(*_nodes[id]);
+ }
+
+ _property_meta_data_vecs.push_back(face_set_property);
+}
+
+// Reads given number of bits (rounded up to next byte) into a bitset.
+// Used for reading region information which can be represented by some
+// number of bits.
+Bitset readBits(std::ifstream& in, const std::size_t bits)
+{
+ using block_t = Bitset::block_type;
+ std::size_t const bytes = static_cast<std::size_t>(std::ceil(bits / 8.));
+ std::size_t const blocks =
+ bytes + 1 < sizeof(block_t) ? 1 : (bytes + 1) / sizeof(block_t);
+
+ std::vector<block_t> data;
+ data.resize(blocks);
+ std::fill_n(data.data(), blocks, 0);
+ in.read(reinterpret_cast<char*>(data.data()), bytes);
+
+ return Bitset(data.begin(), data.end());
+}
+
+void GocadSGridReader::readNodesBinary()
+{
+ std::ifstream in(_pnts_fname.c_str(), std::ios::binary);
+ if (!in)
+ {
+ ERR("Could not open points file '%s'.", _pnts_fname.c_str());
+ throw std::runtime_error("Could not open points file.");
+ }
+
+ std::size_t const n = _index_calculator._n_nodes;
+ _nodes.resize(n);
+
+ double coords[3];
+
+ std::size_t k = 0;
+ while (in && k < n * 3)
+ {
+ if (_bin_pnts_in_double_precision)
+ {
+ coords[k % 3] =
+ BaseLib::swapEndianness(BaseLib::readBinaryValue<double>(in));
+ }
+ else
+ {
+ coords[k % 3] =
+ BaseLib::swapEndianness(BaseLib::readBinaryValue<float>(in));
+ }
+ if ((k + 1) % 3 == 0)
+ {
+ const std::size_t layer_transition_idx(
+ _index_calculator.getCoordsForID(k / 3)[2]);
+ _nodes[k / 3] =
+ new GocadNode(coords, k / 3, layer_transition_idx);
+ }
+ k++;
+ }
+ if (k != n * 3 && !in.eof())
+ ERR("Read different number of points. Expected %d floats, got %d.\n",
+ n * 3, k);
+}
+
+void GocadSGridReader::mapRegionFlagsToCellProperties(
+ std::vector<Bitset> const& rf)
+{
+ DBUG(
+ "GocadSGridReader::mapRegionFlagsToCellProperties region_flags.size: "
+ "%d",
+ rf.size());
+
+ Gocad::Property region_flags;
+ region_flags._property_id = 0;
+ region_flags._property_name = "RegionFlags";
+ region_flags._property_class_name = "RegionFlags";
+ region_flags._property_unit = "unitless";
+ region_flags._property_data_type = "int";
+ region_flags._property_data_fname = "";
+ region_flags._property_no_data_value = -1;
+ std::size_t const n = _index_calculator._n_cells;
+ region_flags._property_data.resize(n);
+ std::fill(region_flags._property_data.begin(),
+ region_flags._property_data.end(), -1);
+
+ // region flags are stored in each node ijk and give the region index for
+ // the ijk-th cell.
+ for (std::size_t i(0); i < _index_calculator._x_dim - 1; i++)
+ {
+ for (std::size_t j(0); j < _index_calculator._y_dim - 1; j++)
+ {
+ for (std::size_t k(0); k < _index_calculator._z_dim - 1; k++)
+ {
+ std::size_t const cell_id(
+ _index_calculator.getCellIdx(i, j, k));
+ std::size_t const node_id(_index_calculator({i, j, k}));
+ for (auto& region : regions)
+ {
+ if (rf[node_id].test(region.bit))
+ {
+ region_flags._property_data[cell_id] += region.bit + 1;
+ }
+ }
+ }
+ }
+ }
+
+ _property_meta_data_vecs.push_back(region_flags);
+}
+
+void GocadSGridReader::readElementPropertiesBinary()
+{
+ for (auto prop_it(_property_meta_data_vecs.begin());
+ prop_it != _property_meta_data_vecs.end();
+ prop_it++)
+ {
+ std::string const& fname(prop_it->_property_data_fname);
+ if (prop_it->_property_data_fname.empty())
+ {
+ WARN("Empty filename for property %s.",
+ prop_it->_property_name.c_str());
+ continue;
+ }
+ std::vector<float> float_properties =
+ BaseLib::readBinaryArray<float>(fname, _index_calculator._n_cells);
+ DBUG(
+ "GocadSGridReader::readElementPropertiesBinary(): Read %d float "
+ "properties from binary file.",
+ _index_calculator._n_cells);
+
+ std::transform(float_properties.cbegin(), float_properties.cend(),
+ float_properties.begin(), [](float const& val) {
+ return BaseLib::swapEndianness(val);
+ });
+
+ prop_it->_property_data.resize(float_properties.size());
+ std::copy(float_properties.begin(), float_properties.end(),
+ prop_it->_property_data.begin());
+ if (prop_it->_property_data.empty())
+ {
+ ERR("Reading of element properties file '%s' failed.",
+ fname.c_str());
+ }
+ }
+}
+
+std::vector<int> GocadSGridReader::readFlagsBinary() const
+{
+ std::vector<int> result;
+ if (!_double_precision_binary)
+ {
+ result = BaseLib::readBinaryArray<int32_t>(_flags_fname,
+ _index_calculator._n_nodes);
+ std::for_each(result.begin(), result.end(),
+ [](int32_t& val) { BaseLib::swapEndianness(val); });
+ }
+ else
+ {
+ result = BaseLib::readBinaryArray<int>(_flags_fname,
+ _index_calculator._n_nodes);
+ std::for_each(result.begin(), result.end(),
+ [](int& val) { BaseLib::swapEndianness(val); });
+ }
+
+ if (result.empty())
+ ERR("Reading of flags file '%s' failed.", _flags_fname.c_str());
+
+ return result;
+}
+
+std::vector<Bitset> GocadSGridReader::readRegionFlagsBinary() const
+{
+ std::vector<Bitset> result;
+
+ std::ifstream in(_region_flags_fname.c_str());
+ if (!in)
+ {
+ ERR("readRegionFlagsBinary(): Could not open file '%s' for input.\n",
+ _region_flags_fname.c_str());
+ in.close();
+ return result;
+ }
+
+ std::size_t const n = _index_calculator._n_nodes;
+ result.resize(n);
+
+ std::size_t k = 0;
+ while (in && k < n)
+ {
+ result[k++] = readBits(in, regions.size());
+ }
+ if (k != n && !in.eof())
+ ERR("Read different number of values. Expected %d, got %d.\n", n, k);
+
+ return result;
+}
+std::vector<MeshLib::Element*> GocadSGridReader::createElements(
+ std::vector<MeshLib::Node*> const& nodes) const
+{
+ std::vector<MeshLib::Element*> elements;
+ elements.resize(_index_calculator._n_cells);
+ std::array<MeshLib::Node*, 8> element_nodes{};
+ std::size_t cnt(0);
+ for (std::size_t k(0); k < _index_calculator._z_dim - 1; k++)
+ {
+ for (std::size_t j(0); j < _index_calculator._y_dim - 1; j++)
+ {
+ for (std::size_t i(0); i < _index_calculator._x_dim - 1; i++)
+ {
+ element_nodes[0] = nodes[_index_calculator({i, j, k})];
+ element_nodes[1] = nodes[_index_calculator({i + 1, j, k})];
+ element_nodes[2] = nodes[_index_calculator({i + 1, j + 1, k})];
+ element_nodes[3] = nodes[_index_calculator({i, j + 1, k})];
+ element_nodes[4] = nodes[_index_calculator({i, j, k + 1})];
+ element_nodes[5] = nodes[_index_calculator({i + 1, j, k + 1})];
+ element_nodes[6] =
+ nodes[_index_calculator({i + 1, j + 1, k + 1})];
+ element_nodes[7] = nodes[_index_calculator({i, j + 1, k + 1})];
+ elements[cnt] = new MeshLib::Hex(
+ element_nodes, _index_calculator.getCellIdx(i, j, k));
+ cnt++;
+ }
+ }
+ }
+ return elements;
+}
+
+void GocadSGridReader::readSplitInformation()
+{
+ std::ifstream in(_fname.c_str());
+ if (!in)
+ {
+ ERR("Could not open '%s'.", _fname.c_str());
+ in.close();
+ return;
+ }
+
+ // read split information from the stratigraphic grid file
+ std::string line;
+ std::stringstream ss;
+ while (std::getline(in, line))
+ {
+ std::size_t pos(line.find("SPLIT "));
+ if (pos != std::string::npos)
+ {
+ ss << line.substr(pos + 6, line.size() - (pos + 6));
+ // read position in grid
+ std::array<std::size_t, 3> grid_coords{};
+ ss >> grid_coords[0];
+ ss >> grid_coords[1];
+ ss >> grid_coords[2];
+ // read coordinates for the split node
+ double coords[3];
+ ss >> coords[0];
+ ss >> coords[1];
+ ss >> coords[2];
+ // read the id
+ std::size_t id;
+ ss >> id;
+ // read the affected cells
+ std::bitset<8> affected_cells{};
+ for (std::size_t ac = 0; ac < 8; ++ac)
+ {
+ char bit;
+ ss >> bit;
+ affected_cells[ac] = bit == '0' ? false : true;
+ }
+ const std::size_t layer_transition_index(
+ _nodes[id]->getLayerTransitionIndex());
+ _split_nodes.push_back(new GocadSplitNode(coords, id, grid_coords,
+ affected_cells,
+ layer_transition_index));
+ }
+ }
+}
+
+void GocadSGridReader::applySplitInformation(
+ std::vector<MeshLib::Node*>& nodes,
+ std::vector<MeshLib::Element*>& elements) const
+{
+ for (auto split_node : _split_nodes)
+ {
+ std::size_t const new_node_pos(nodes.size());
+ nodes.push_back(
+ new MeshLib::Node(split_node->getCoords(), new_node_pos));
+
+ // get grid coordinates
+ std::array<std::size_t, 3> const& gc(split_node->getGridCoords());
+ // get affected cells
+ auto const& affected_cells(split_node->getAffectedCells());
+ // get mesh node to substitute in elements
+ MeshLib::Node const* const node2sub(nodes[_index_calculator(gc)]);
+
+ if (affected_cells[0] && gc[0] < _index_calculator._x_dim - 1 &&
+ gc[1] < _index_calculator._y_dim - 1 &&
+ gc[2] < _index_calculator._z_dim - 1)
+ {
+ const std::size_t idx(
+ _index_calculator.getCellIdx(gc[0], gc[1], gc[2]));
+ modifyElement(elements[idx], node2sub, nodes[new_node_pos]);
+ }
+ if (affected_cells[1] && gc[0] > 0 &&
+ gc[1] < _index_calculator._y_dim - 1 &&
+ gc[2] < _index_calculator._z_dim - 1)
+ {
+ const std::size_t idx(
+ _index_calculator.getCellIdx(gc[0] - 1, gc[1], gc[2]));
+ modifyElement(elements[idx], node2sub, nodes[new_node_pos]);
+ }
+ if (affected_cells[2] && gc[1] > 0 &&
+ gc[0] < _index_calculator._x_dim - 1 &&
+ gc[2] < _index_calculator._z_dim - 1)
+ {
+ const std::size_t idx(
+ _index_calculator.getCellIdx(gc[0], gc[1] - 1, gc[2]));
+ modifyElement(elements[idx], node2sub, nodes[new_node_pos]);
+ }
+ if (affected_cells[3] && gc[0] > 0 && gc[1] > 0 &&
+ gc[2] < _index_calculator._z_dim - 1)
+ {
+ const std::size_t idx(
+ _index_calculator.getCellIdx(gc[0] - 1, gc[1] - 1, gc[2]));
+ modifyElement(elements[idx], node2sub, nodes[new_node_pos]);
+ }
+ if (affected_cells[4] && gc[2] > 0 &&
+ gc[0] < _index_calculator._x_dim - 1 &&
+ gc[1] < _index_calculator._y_dim - 1)
+ {
+ const std::size_t idx(
+ _index_calculator.getCellIdx(gc[0], gc[1], gc[2] - 1));
+ modifyElement(elements[idx], node2sub, nodes[new_node_pos]);
+ }
+ if (affected_cells[5] && gc[0] > 0 && gc[2] > 0 &&
+ gc[1] < _index_calculator._y_dim - 1)
+ {
+ const std::size_t idx(
+ _index_calculator.getCellIdx(gc[0] - 1, gc[1], gc[2] - 1));
+ modifyElement(elements[idx], node2sub, nodes[new_node_pos]);
+ }
+ if (affected_cells[6] && gc[1] > 0 && gc[2] > 0 &&
+ gc[0] < _index_calculator._x_dim - 1)
+ {
+ const std::size_t idx(
+ _index_calculator.getCellIdx(gc[0], gc[1] - 1, gc[2] - 1));
+ modifyElement(elements[idx], node2sub, nodes[new_node_pos]);
+ }
+ if (affected_cells[7] && gc[0] > 0 && gc[1] > 0 && gc[2] > 0)
+ {
+ const std::size_t idx(
+ _index_calculator.getCellIdx(gc[0] - 1, gc[1] - 1, gc[2] - 1));
+ modifyElement(elements[idx], node2sub, nodes[new_node_pos]);
+ }
+ }
+}
+
+void GocadSGridReader::modifyElement(MeshLib::Element* hex,
+ MeshLib::Node const* node2sub,
+ MeshLib::Node* substitute_node) const
+{
+ // get the node pointers of the cell
+ MeshLib::Node* const* hex_nodes(hex->getNodes());
+ // search for the position the split node will be set to
+ MeshLib::Node* const* node_pos(
+ std::find(hex_nodes, hex_nodes + 8, node2sub));
+ // set the split node instead of the node2sub
+ if (node_pos != hex_nodes + 8)
+ {
+ const_cast<MeshLib::Node**>(
+ hex_nodes)[std::distance(hex_nodes, node_pos)] = substitute_node;
+ }
+}
+
+boost::optional<Gocad::Property const&> GocadSGridReader::getProperty(
+ std::string const& name) const
+{
+ auto const it(std::find_if(_property_meta_data_vecs.begin(),
+ _property_meta_data_vecs.end(),
+ [&name](Gocad::Property const& p) {
+ return p._property_name == name;
+ }));
+ if (it == _property_meta_data_vecs.end())
+ {
+ return boost::optional<Gocad::Property const&>();
+ }
+ return boost::optional<Gocad::Property const&>(*it);
+}
+
+std::vector<std::string> GocadSGridReader::getPropertyNames() const
+{
+ std::vector<std::string> names;
+ std::transform(_property_meta_data_vecs.begin(),
+ _property_meta_data_vecs.end(),
+ std::back_inserter(names),
+ [](Gocad::Property const& p) { return p._property_name; });
+ return names;
+}
+
+std::unique_ptr<MeshLib::Mesh> GocadSGridReader::getFaceSetMesh(
+ std::size_t const face_set_number) const
+{
+ std::vector<MeshLib::Node*> face_set_nodes;
+ std::vector<MeshLib::Element*> face_set_elements;
+
+ for (auto const node : _nodes)
+ {
+ if (node->isMemberOfFaceSet(face_set_number))
+ {
+ addFaceSetQuad(node, face_set_number, face_set_nodes,
+ face_set_elements);
+ }
+ }
+
+ if (face_set_nodes.empty())
+ {
+ return nullptr;
+ }
+
+ for (auto const node : _split_nodes)
+ {
+ if (node->isMemberOfFaceSet(face_set_number))
+ {
+ if (node->getAffectedCells()[0])
+ {
+ addFaceSetQuad(node, face_set_number, face_set_nodes,
+ face_set_elements);
+ }
+ }
+ }
+
+ std::string const mesh_name("FaceSet-" + std::to_string(face_set_number));
+ return std::make_unique<MeshLib::Mesh>(mesh_name, face_set_nodes,
+ face_set_elements);
+}
+
+void GocadSGridReader::addFaceSetQuad(
+ GocadNode* face_set_node, std::size_t face_set_number,
+ std::vector<MeshLib::Node*>& face_set_nodes,
+ std::vector<MeshLib::Element*>& face_set_elements) const
+{
+ std::array<MeshLib::Node*, 4> quad_nodes{};
+ quad_nodes[0] = new GocadNode(*face_set_node);
+ const std::size_t id(face_set_node->getID());
+ std::array<std::size_t, 3> const c(_index_calculator.getCoordsForID(id));
+
+ const FaceDirection dir(face_set_node->getFaceDirection(face_set_number));
+ switch (dir)
+ {
+ case FaceDirection::U:
+ quad_nodes[1] = new GocadNode(
+ *_nodes[_index_calculator({c[0], c[1] + 1, c[2]})]);
+ quad_nodes[2] = new GocadNode(
+ *_nodes[_index_calculator({c[0], c[1] + 1, c[2] + 1})]);
+ quad_nodes[3] = new GocadNode(
+ *_nodes[_index_calculator({c[0], c[1], c[2] + 1})]);
+ break;
+ case FaceDirection::V:
+ quad_nodes[1] = new GocadNode(
+ *_nodes[_index_calculator({c[0] + 1, c[1], c[2]})]);
+ quad_nodes[2] = new GocadNode(
+ *_nodes[_index_calculator({c[0] + 1, c[1], c[2] + 1})]);
+ quad_nodes[3] = new GocadNode(
+ *_nodes[_index_calculator({c[0], c[1], c[2] + 1})]);
+ break;
+ case FaceDirection::W:
+ quad_nodes[1] = new GocadNode(
+ *_nodes[_index_calculator({c[0] + 1, c[1], c[2]})]);
+ quad_nodes[2] = new GocadNode(
+ *_nodes[_index_calculator({c[0] + 1, c[1] + 1, c[2]})]);
+ quad_nodes[3] = new GocadNode(
+ *_nodes[_index_calculator({c[0], c[1] + 1, c[2]})]);
+ break;
+ default:
+ ERR("Could not create face for node with id %d.", id);
+ }
+ for (auto quad_node : quad_nodes)
+ {
+ face_set_nodes.push_back(quad_node);
+ }
+ face_set_elements.push_back(new MeshLib::Quad(quad_nodes));
+}
+
+} // end namespace Gocad
+} // end namespace FileIO
diff --git a/Applications/FileIO/GocadIO/GocadSGridReader.h b/Applications/FileIO/GocadIO/GocadSGridReader.h
new file mode 100644
index 0000000000000000000000000000000000000000..1c21b8a84cac782bf0d154837b871c942cbb51e6
--- /dev/null
+++ b/Applications/FileIO/GocadIO/GocadSGridReader.h
@@ -0,0 +1,115 @@
+/**
+ *
+ * @copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/LICENSE.txt
+ */
+
+#pragma once
+
+#include <cstddef>
+#include <limits>
+#include <numeric>
+#include <string>
+#include <vector>
+
+#include <boost/dynamic_bitset.hpp>
+
+#include <logog/include/logog.hpp>
+
+#include "MeshLib/Elements/Element.h"
+
+#include "CoordinateSystem.h"
+#include "GocadNode.h"
+#include "IndexCalculator.h"
+#include "Layer.h"
+#include "Property.h"
+#include "Region.h"
+
+namespace FileIO
+{
+namespace Gocad
+{
+class GocadSGridReader final
+{
+public:
+ /**
+ * Constructor takes as argument the Gocad .sg text file.
+ * @param fname file name
+ */
+ explicit GocadSGridReader(std::string const& fname);
+ ~GocadSGridReader();
+
+ GocadSGridReader() = delete;
+ GocadSGridReader(GocadSGridReader&& src) = delete;
+ GocadSGridReader(GocadSGridReader const& src) = delete;
+ GocadSGridReader& operator=(GocadSGridReader&& rhs) = delete;
+ GocadSGridReader& operator=(GocadSGridReader const& rhs) = delete;
+
+ std::unique_ptr<MeshLib::Mesh> getMesh() const;
+ std::unique_ptr<MeshLib::Mesh> getFaceSetMesh(
+ std::size_t const face_set_number) const;
+
+ std::vector<std::string> getPropertyNames() const;
+
+private:
+ using Bitset = boost::dynamic_bitset<>;
+
+ void parseDims(std::string const& line);
+ void parseFileName(std::string const& line,
+ std::string& result_string) const;
+ void parseHeader(std::istream& in);
+ void parseFaceSet(std::string& line, std::istream& in);
+
+ void readNodesBinary();
+ std::vector<int> readFlagsBinary() const;
+ std::vector<Bitset> readRegionFlagsBinary() const;
+ void readElementPropertiesBinary();
+ void mapRegionFlagsToCellProperties(std::vector<Bitset> const& rf);
+
+ std::vector<MeshLib::Element*> createElements(
+ std::vector<MeshLib::Node*> const& nodes) const;
+
+ // split handling
+ void readSplitInformation();
+ void applySplitInformation(std::vector<MeshLib::Node*>& nodes,
+ std::vector<MeshLib::Element*>& elements) const;
+ void modifyElement(MeshLib::Element* hex, MeshLib::Node const* node2sub,
+ MeshLib::Node* substitute_node) const;
+
+ void addFaceSetQuad(
+ GocadNode* face_set_node, std::size_t face_set_number,
+ std::vector<MeshLib::Node*>& face_set_nodes,
+ std::vector<MeshLib::Element*>& face_set_elements) const;
+
+ boost::optional<Gocad::Property const&> getProperty(
+ std::string const& name) const;
+ void addGocadPropertiesToMesh(MeshLib::Mesh& mesh) const;
+
+ std::string const& _fname;
+ std::string const _path;
+ // data read from sg file
+ Gocad::IndexCalculator _index_calculator;
+ Gocad::CoordinateSystem _coordinate_system;
+ std::string _pnts_fname;
+ std::string _flags_fname;
+ std::string _region_flags_fname;
+
+ std::vector<Gocad::Region> regions;
+ std::vector<Gocad::Layer> layers;
+ std::size_t _n_face_sets;
+
+ bool _double_precision_binary;
+ bool _bin_pnts_in_double_precision;
+
+ // data read from binary points file
+ std::vector<GocadNode*> _nodes;
+ std::vector<GocadSplitNode*> _split_nodes;
+ // properties
+ std::vector<Gocad::Property> _property_meta_data_vecs;
+}; // end class GocadSGridReader
+
+} // end namespace Gocad
+} // end namespace FileIO
diff --git a/Applications/Utils/FileConverter/CMakeLists.txt b/Applications/Utils/FileConverter/CMakeLists.txt
index 34f7437ba3ba830dd7489b917c11213b238f5306..8db49249a1b3908de206a4fd0417e26890d23e8e 100644
--- a/Applications/Utils/FileConverter/CMakeLists.txt
+++ b/Applications/Utils/FileConverter/CMakeLists.txt
@@ -43,11 +43,19 @@ add_executable(TecPlotTools TecPlotTools.cpp)
set_target_properties(TecPlotTools PROPERTIES FOLDER Utilities)
target_link_libraries(TecPlotTools GeoLib MeshLib)
+add_executable(GocadSGridReader GocadSGridReaderMain.cpp)
+set_target_properties(GocadSGridReader PROPERTIES FOLDER Utilities)
+target_link_libraries(GocadSGridReader
+ GeoLib
+ MeshLib
+ ApplicationsFileIO
+ ${Boost_LIBRARIES}
+)
####################
### Installation ###
####################
-install(TARGETS generateMatPropsFromMatID GMSH2OGS OGS2VTK VTK2OGS VTK2TIN
+install(TARGETS generateMatPropsFromMatID GMSH2OGS OGS2VTK VTK2OGS VTK2TIN GocadSGridReader
RUNTIME DESTINATION bin COMPONENT ogs_converter)
if(Qt5XmlPatterns_FOUND)
diff --git a/Applications/Utils/FileConverter/GocadSGridReaderMain.cpp b/Applications/Utils/FileConverter/GocadSGridReaderMain.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..3fbf80f92eae00354097f15a950f33e31cbadf59
--- /dev/null
+++ b/Applications/Utils/FileConverter/GocadSGridReaderMain.cpp
@@ -0,0 +1,95 @@
+/**
+ *
+ * @copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/LICENSE.txt
+ */
+
+#include <fstream>
+#include <sstream>
+#include <string>
+
+#include <logog/include/logog.hpp>
+#include <tclap/CmdLine.h>
+
+#include "Applications/ApplicationsLib/LogogSetup.h"
+#include "Applications/FileIO/GocadIO/GenerateFaceSetMeshes.h"
+#include "Applications/FileIO/GocadIO/GocadSGridReader.h"
+#include "BaseLib/BuildInfo.h"
+#include "BaseLib/FileTools.h"
+#include "MeshLib/Elements/Element.h"
+#include "MeshLib/IO/writeMeshToFile.h"
+#include "MeshLib/Mesh.h"
+#include "MeshLib/Node.h"
+
+int main(int argc, char* argv[])
+{
+ ApplicationsLib::LogogSetup logog_setup;
+
+ TCLAP::CmdLine cmd(
+ "Reads a Gocad stratigraphic grid file (file ending sg) and writes the "
+ "data in the vtk unstructured grid file format.\n\n OpenGeoSys-6 "
+ "software, version " +
+ BaseLib::BuildInfo::git_describe +
+ ".\n"
+ "Copyright (c) 2012-2019, OpenGeoSys Community "
+ "(http://www.opengeosys.org)",
+ ' ', BaseLib::BuildInfo::git_describe);
+
+ TCLAP::ValueArg<bool> face_set_arg(
+ "f", "generate-face-sets",
+ "generate face sets; default false, i.e. do not generate face sets",
+ false, false, "true/false");
+ cmd.add(face_set_arg);
+
+ TCLAP::ValueArg<std::string> mesh_output_arg(
+ "o", "output-mesh", "vtk unstructured grid file name", true, "",
+ "file_name.vtu");
+ cmd.add(mesh_output_arg);
+
+ TCLAP::ValueArg<std::string> sg_file_arg(
+ "s", "sg", "Gocad stratigraphic grid file name", true, "",
+ "file_name.sg");
+ cmd.add(sg_file_arg);
+
+ cmd.parse(argc, argv);
+
+ // read the Gocad SGrid
+ INFO("Start reading Gocad SGrid.");
+ FileIO::Gocad::GocadSGridReader reader(sg_file_arg.getValue());
+ INFO("End reading Gocad SGrid.");
+
+ if (face_set_arg.getValue())
+ {
+ INFO("Generating a mesh for every face set.");
+ FileIO::Gocad::generateFaceSets(
+ reader, BaseLib::extractPath(sg_file_arg.getValue()));
+ }
+
+ std::unique_ptr<MeshLib::Mesh> mesh(reader.getMesh());
+
+ INFO("The mesh contains the following PropertyVectors:");
+ auto property_names(mesh->getProperties().getPropertyVectorNames());
+ for (auto const& property_name : property_names)
+ {
+ INFO("- %s (#values: %d)", property_name.c_str(),
+ mesh->getProperties()
+ .getPropertyVector<double>(property_name)
+ ->size());
+ auto bounds(
+ std::minmax_element(mesh->getProperties()
+ .getPropertyVector<double>(property_name)
+ ->cbegin(),
+ mesh->getProperties()
+ .getPropertyVector<double>(property_name)
+ ->cend()));
+ INFO("\tvalues in range [%e, %e].", *(bounds.first), *(bounds.second));
+ }
+
+ INFO("Writing mesh to '%s'.", mesh_output_arg.getValue().c_str());
+ MeshLib::IO::writeMeshToFile(*mesh, mesh_output_arg.getValue());
+
+ return 0;
+}