diff --git a/CMakeLists.txt b/CMakeLists.txt
index d908ffcd09107d4b9e99e41dbcc45f19b5d6f7ee..039716e734ac529447e7afbfa15b97df855246b6 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -5,6 +5,7 @@ set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH}
${CMAKE_SOURCE_DIR}/ogs/scripts/cmake/cmake
${CMAKE_SOURCE_DIR}/ogs/scripts/cmake)
include(CompilerSetup)
+include(ProjectSetup)
add_subdirectory(${CMAKE_SOURCE_DIR}/ogs)
diff --git a/ErtData2Mesh/ErtData2Mesh.cpp b/ErtData2Mesh/ErtData2Mesh.cpp
index ed33b4d6eec15bf5ccb984d46f53f46de0f11db7..76593d8c679c322db8e0b9c69c18e721f9b85458 100644
--- a/ErtData2Mesh/ErtData2Mesh.cpp
+++ b/ErtData2Mesh/ErtData2Mesh.cpp
@@ -26,9 +26,11 @@
// FileIO
#include "FileIO/VtkIO/VtuInterface.h"
#include "FileIO/CsvInterface.h"
+#include "FileIO/AsciiRasterInterface.h"
// GeoLib
#include "GeoLib/Point.h"
+#include "GeoLib/Raster.h"
// MeshLib
#include "MeshLib/Mesh.h"
@@ -40,6 +42,20 @@
#include "MeshLib/MeshEditing/MeshRevision.h"
+std::vector<double> getElevationCorrectionValues(GeoLib::Raster const& dem, std::vector<MeshLib::Node*> &nodes)
+{
+ std::size_t const n_sfc_nodes (nodes.size());
+ std::vector<double> elevation_correction;
+ for (std::size_t i=0; i<n_sfc_nodes; ++i)
+ {
+ double a = (*nodes[i])[2];
+ double b = dem.getValueAtPoint(*nodes[i]);
+ elevation_correction.push_back((*nodes[i])[2] - dem.getValueAtPoint(*nodes[i]));
+ (*nodes[i])[2] -= elevation_correction[i];
+ }
+ return elevation_correction;
+}
+
int main (int argc, char* argv[])
{
LOGOG_INITIALIZE();
@@ -58,53 +74,107 @@ int main (int argc, char* argv[])
"CSV-file containing ERT information.",
true, "", "name of the csv input file");
cmd.add(csv_in);
+ TCLAP::ValueArg<std::string> dem_in("s", "DEM-file",
+ "Surface DEM for mapping ERT data", false,
+ "", "file name of the Surface DEM");
+ cmd.add(dem_in);
cmd.parse(argc, argv);
+ std::vector<std::size_t> ids;
+ int const e = FileIO::CsvInterface::readColumn<std::size_t>(csv_in.getValue(), '\t', ids, "x2/m");
+ std::size_t const n_nodes_per_layer (ids.back()+1);
std::vector<GeoLib::Point*> points1;
std::vector<GeoLib::Point*> points2;
- int const r1 = FileIO::CsvInterface::readPoints("hang.txt", '\t', points1, "E1", "N1", "H1");
- int const r2 = FileIO::CsvInterface::readPoints("hang.txt", '\t', points2, "E2", "N2", "H2");
+ int const r1 = FileIO::CsvInterface::readPoints(csv_in.getValue(), '\t', points1, "E1", "N1", "H1");
+ int const r2 = FileIO::CsvInterface::readPoints(csv_in.getValue(), '\t', points2, "E2", "N2", "H2");
std::vector<double> z1;
std::vector<double> z2;
- int const c1 = FileIO::CsvInterface::readColumn<double>("hang.txt", '\t', z1, "z1/m");
- int const c2 = FileIO::CsvInterface::readColumn<double>("hang.txt", '\t', z2, "z2/m");
+ int const c1 = FileIO::CsvInterface::readColumn<double>(csv_in.getValue(), '\t', z1, "z1/m");
+ int const c2 = FileIO::CsvInterface::readColumn<double>(csv_in.getValue(), '\t', z2, "z2/m");
- if (r1 < 0 || r2 < 0 || c1 < 0 || c2 < 0)
+ if (e!= 0 || r1 != 0 || r2 != 0 || c1 != 0 || c2 != 0)
{
ERR("Error reading data from file");
return 1;
}
+ std::size_t const n_quads (points1.size());
+ for (std::size_t i=1; i<n_quads; ++i)
+ {
+ if ((*points1[i-1])[0] == (*points2[i])[0] ||
+ (*points1[i-1])[1] == (*points2[i])[1] ||
+ (*points1[i-1])[2] == (*points2[i])[2])
+ {
+ ERR ("Error in ERT file.");
+ return 1;
+ }
+ }
std::vector<MeshLib::Node*> nodes;
- std::size_t const n_quads (points1.size());
std::vector<MeshLib::Element*> quads;
quads.reserve(n_quads);
std::vector<int> materials;
materials.reserve(n_quads);
- int mat_id(0);
- for (std::size_t i=0; i<n_quads; ++i)
+
+ nodes.push_back(new MeshLib::Node((*points1[0])[0], (*points1[0])[1], (*points1[0])[2] - z1[0], nodes.size() ));
+ for (std::size_t i=0; i<n_nodes_per_layer; ++i)
+ nodes.push_back(new MeshLib::Node((*points2[i])[0], (*points2[i])[1], (*points2[i])[2] - z1[i], nodes.size() ));
+
+ std::vector<double> elevation_correction (nodes.size(), 0.0);
+ if (dem_in.isSet())
{
- std::array<MeshLib::Node*, 4> quad_nodes;
- quad_nodes[0] = new MeshLib::Node((*points1[i])[0], (*points1[i])[1], (*points1[i])[2] - z1[i]);
- quad_nodes[1] = new MeshLib::Node((*points1[i])[0], (*points1[i])[1], (*points1[i])[2] - z2[i]);
- quad_nodes[2] = new MeshLib::Node((*points2[i])[0], (*points2[i])[1], (*points2[i])[2] - z2[i]);
- quad_nodes[3] = new MeshLib::Node((*points2[i])[0], (*points2[i])[1], (*points2[i])[2] - z1[i]);
- if (i>0 && (z1[i] != z1[i-1]))
- mat_id++;
- for (std::size_t j=0; j<4; ++j)
- nodes.push_back(quad_nodes[j]);
- quads.push_back(new MeshLib::Quad(quad_nodes));
- materials.push_back(mat_id);
+ GeoLib::Raster* dem = FileIO::AsciiRasterInterface::readRaster(dem_in.getValue());
+ elevation_correction = getElevationCorrectionValues(*dem, nodes);
+ delete dem;
}
-
+ std::size_t const n_layers ((n_quads / n_nodes_per_layer));
+ for (std::size_t mat_id=0; mat_id<n_layers; ++mat_id)
+ {
+ std::size_t const base_idx (mat_id * n_nodes_per_layer);
+ double z = (*points1[base_idx])[2] - elevation_correction[0] - z2[base_idx];
+ nodes.push_back(new MeshLib::Node((*points1[base_idx])[0], (*points1[base_idx])[1], z, nodes.size() ));
+ for (std::size_t i=0; i<n_nodes_per_layer; ++i)
+ {
+ std::size_t const idx (base_idx + i);
+ std::array<MeshLib::Node*, 4> quad_nodes;
+ quad_nodes[0] = nodes[base_idx + mat_id + i];
+ quad_nodes[1] = nodes[nodes.size()-1];
+ z = (*points2[idx])[2] - elevation_correction[i+1] - z2[base_idx+i];
+ nodes.push_back(new MeshLib::Node((*points2[idx])[0], (*points2[idx])[1], z, nodes.size() ));
+ quad_nodes[2] = nodes.back();
+ quad_nodes[3] = nodes[base_idx + mat_id + i+1];
+ quads.push_back(new MeshLib::Quad(quad_nodes));
+ materials.push_back(mat_id);
+ }
+ }
+
MeshLib::Mesh mesh("ERT Mesh", nodes, quads);
+
boost::optional<std::vector<int>&> mat_prop (mesh.getProperties().createNewPropertyVector<int>("MaterialIDs", MeshLib::MeshItemType::Cell));
mat_prop->insert(mat_prop->end(), materials.begin(), materials.end());
- /*
- MeshLib::MeshRevision rev(mesh);
- MeshLib::Mesh* final_mesh (rev.collapseNodes("ERT Mesh", std::numeric_limits<double>::epsilon()));
+ boost::optional<std::vector<double>&> resistance (mesh.getProperties().createNewPropertyVector<double>("Resistance", MeshLib::MeshItemType::Cell));
+ int const errors_res = FileIO::CsvInterface::readColumn<double>(csv_in.getValue(), '\t', *resistance, "rho/Ohmm ");
+ if (errors_res != 0 || resistance->size() != n_quads)
+ WARN ("Erros reading resistance values.");
+
+ boost::optional<std::vector<double>&> coverage (mesh.getProperties().createNewPropertyVector<double>("Coverage", MeshLib::MeshItemType::Cell));
+ int const errors_cov = FileIO::CsvInterface::readColumn<double>(csv_in.getValue(), '\t', *coverage, "coverage");
+ if (errors_cov != 0 || coverage->size() != n_quads)
+ WARN ("Erros reading coverage values.");
+
+ /* TODO writing arrays with tuple size > 1 needs fixing
+ double const max_cov = *std::max_element(coverage->cbegin(), coverage->cend());
+ boost::optional<std::vector<double>&> conduct (mesh.getProperties().createNewPropertyVector<double>("Conductivity", MeshLib::MeshItemType::Cell, 1));
+ if (errors_res == 0 && errors_cov == 0)
+ {
+ for (std::size_t i=0; i<n_quads; ++i)
+ {
+ conduct->push_back(1.0 / (*resistance)[i]);
+ conduct->push_back((*coverage)[i] / max_cov);
+ }
+ }
*/
+
INFO ("Writing result...");
FileIO::VtuInterface vtu(&mesh);
vtu.writeToFile(mesh_out.getValue());