diff --git a/Gui/DataView/DirectConditionGenerator.cpp b/Gui/DataView/DirectConditionGenerator.cpp
index 40e244c9d6cef914d89ea354b35241db6cf0cd08..45be168dd3db9f7f5cc071be38e59d0cd8bb02cd 100644
--- a/Gui/DataView/DirectConditionGenerator.cpp
+++ b/Gui/DataView/DirectConditionGenerator.cpp
@@ -23,6 +23,7 @@
#include "MeshSurfaceExtraction.h"
#include "PointWithID.h"
#include "Mesh.h"
+#include "Node.h"
#include <cmath>
#include <limits>
@@ -31,42 +32,23 @@ const std::vector< std::pair<size_t,double> >& DirectConditionGenerator::directT
{
if (_direct_values.empty())
{
- GeoLib::Raster* raster(GeoLib::Raster::getRasterFromASCFile(filename));
+ GeoLib::Raster* raster (GeoLib::Raster::readRaster(filename));
if (! raster) {
- ERR("Error in DirectConditionGenerator::directWithSurfaceIntegration() - could not load raster file.");
+ ERR("Error in DirectConditionGenerator::directToSurfaceNodes() - could not load raster file.");
return _direct_values;
}
- double origin_x(raster->getOrigin()[0]);
- double origin_y(raster->getOrigin()[1]);
- double delta(raster->getRasterPixelDistance());
- double no_data(raster->getNoDataValue());
- unsigned imgwidth(raster->getNCols()), imgheight(raster->getNRows());
- double const*const img(raster->begin());
-
- const MathLib::Vector3 dir(0,0,1);
+ const MathLib::Vector3 dir(0,0,-1);
const std::vector<GeoLib::PointWithID*> surface_nodes(MeshLib::MeshSurfaceExtraction::getSurfaceNodes(mesh, dir) );
const size_t nNodes(surface_nodes.size());
+ const double no_data (raster->getNoDataValue());
_direct_values.reserve(nNodes);
for (size_t i=0; i<nNodes; i++)
{
- const double* coords (surface_nodes[i]->getCoords());
-
- if (coords[0]>=origin_x && coords[0]<(origin_x+(delta*imgwidth)) && coords[1]>=origin_y && coords[1]<(origin_y+(delta*imgheight)))
- {
- int cell_x = static_cast<int>(floor((coords[0] - origin_x)/delta));
- int cell_y = static_cast<int>(floor((coords[1] - origin_y)/delta));
-
- // if node outside of raster use raster boundary values
- cell_x = (cell_x < 0) ? 0 : ((cell_x > static_cast<int>(imgwidth )) ? (imgwidth-1) : cell_x);
- cell_y = (cell_y < 0) ? 0 : ((cell_y > static_cast<int>(imgheight)) ? (imgheight-1) : cell_y);
-
- size_t index = cell_y*imgwidth+cell_x;
- if (fabs(img[index] - no_data) > std::numeric_limits<float>::epsilon())
- _direct_values.push_back( std::pair<size_t, double>(surface_nodes[i]->getID(),img[index]) );
- }
+ double val (raster->getValueAtPoint(*surface_nodes[i]));
+ val = (fabs(val-no_data) < std::numeric_limits<double>::epsilon()) ? 0 : val;
+ _direct_values.push_back (std::pair<size_t, double>(surface_nodes[i]->getID(), val));
}
-
delete raster;
}
else
@@ -78,102 +60,34 @@ const std::vector< std::pair<size_t,double> >& DirectConditionGenerator::directT
const std::vector< std::pair<size_t,double> >& DirectConditionGenerator::directWithSurfaceIntegration(MeshLib::Mesh &mesh, const std::string &filename, double scaling)
{
- (void)mesh;
- (void)filename;
- (void)scaling;
-/* TODO6
- double no_data_value (-9999); // TODO: get this from asc-reader!
-
if (_direct_values.empty())
{
- //mesh.MarkInterface_mHM_Hydro_3D(); // mark element faces on the surface
- //----
- const double dir[3] = {0,0,1};
- MeshLib::Mesh* sfc_mesh (MeshLib::MeshSurfaceExtraction::getMeshSurface(mesh, dir));
- std::vector<double> node_area_vec (sfc_mesh->getNNodes());
- MeshLib::MeshSurfaceExtraction::getSurfaceAreaForNodes(sfc_mesh, node_area_vec);
- //----
- double origin_x(0), origin_y(0), delta(0);
- size_t imgwidth(0), imgheight(0);
- double node_val[8] = {0,0,0,0,0,0,0,0}; // maximum possible number of nodes per face (just in case ...)
-
- FiniteElement::CElement* fem ( new FiniteElement::CElement(mesh.GetCoordinateFlag()) );
-
- float* img = 0;
- if (filename.substr(filename.length()-3,3).compare("asc") == 0)
- img = VtkRaster::loadDataFromASC(filename, origin_x, origin_y, imgwidth, imgheight, delta);
- else if (filename.substr(filename.length()-3,3).compare("grd") == 0)
- img = VtkRaster::loadDataFromSurfer(filename, origin_x, origin_y, imgwidth, imgheight, delta);
-
- if (img == 0)
- {
- std::cout << "Error in DirectConditionGenerator::directWithSurfaceIntegration() - could not load vtk raster." << std::endl;
+ GeoLib::Raster* raster (GeoLib::Raster::readRaster(filename));
+ if (!raster) {
+ ERR("Error in DirectConditionGenerator::directWithSurfaceIntegration() - could not load raster file.");
return _direct_values;
}
- const size_t nNodes(mesh.nod_vector.size());
- std::vector<double> val(nNodes, 0.0);
- for(size_t i = 0; i < nNodes; i++)
- mesh.nod_vector[i]->SetMark(false);
-
- // copied from CFEMesh::Precipitation2NeumannBC() by WW
- size_t nFaces = mesh.face_vector.size();
- for(size_t i=0; i<nFaces; i++)
- {
- MeshLib::CElem* elem = mesh.face_vector[i];
- if (!elem->GetMark())
- continue;
-
- // if face is on the surface of the mesh
- size_t nElemNodes = elem->GetNodesNumber(false);
- for(size_t k=0; k<nElemNodes; k++)
- node_val[k] = 0.0;
-
- // get values from the raster for all nodes of the face
- for(size_t k=0; k<nElemNodes; k++)
- {
- double const* const pnt_k (elem->GetNode(k)->getData());
- int cell_x = static_cast<int>(floor((pnt_k[0] - origin_x) / delta));
- int cell_y = static_cast<int>(floor((pnt_k[1] - origin_y) / delta));
-
- // if node outside of raster use raster boundary values
- cell_x = (cell_x < 0) ? 0 : ((static_cast<size_t>(cell_x) > imgwidth) ? (imgwidth-1) : cell_x);
- cell_y = (cell_y < 0) ? 0 : ((static_cast<size_t>(cell_y) > imgheight) ? (imgheight-1) : cell_y);
-
- node_val[k] = img[ 2 * (cell_y * imgwidth + cell_x) ];
- if (fabs(node_val[k] - no_data_value) < std::numeric_limits<double>::epsilon())
- node_val[k] = 0.;
- }
-
- // get area of the surface element face
- elem->ComputeVolume();
-
- // do the actual surface integration
- fem->setOrder(mesh.getOrder() + 1);
- fem->ConfigElement(elem);
- fem->FaceIntegration(node_val);
-
- // add up the integrated values (nodes get values added for all faces they are part of)
- for(size_t k=0; k<elem->GetNodesNumber(false); k++)
- {
- MeshLib::CNode* node = elem->GetNode(k);
- node->SetMark(true);
- val[node->GetIndex()] += node_val[k];
- }
- }
-
+ const MathLib::Vector3 dir(0,0,-1);
+ MeshLib::Mesh* sfc_mesh (MeshLib::MeshSurfaceExtraction::getMeshSurface(mesh, dir, true));
+ std::vector<double> node_area_vec;
+ MeshLib::MeshSurfaceExtraction::getSurfaceAreaForNodes(*sfc_mesh, node_area_vec);
+ const std::vector<MeshLib::Node*> &surface_nodes (sfc_mesh->getNodes());
+ const size_t nNodes(sfc_mesh->getNNodes());
+ const double no_data (raster->getNoDataValue());
_direct_values.reserve(nNodes);
- for(size_t k=0; k<nNodes; k++)
+ for (size_t i=0; i<nNodes; ++i)
{
- if (!mesh.nod_vector[k]->GetMark())
- continue;
- // Assuming the unit of precipitation is mm/day
- _direct_values.push_back( std::pair<size_t, double>(k, val[k] / scaling) );
+ double val (raster->getValueAtPoint(*surface_nodes[i]));
+ val = (fabs(val-no_data) < std::numeric_limits<double>::epsilon()) ? 0 : (val*node_area_vec[i]*scaling);
+ _direct_values.push_back (std::pair<size_t, double>(surface_nodes[i]->getID(), val));
}
+
+ delete raster;
}
else
std::cout << "Error in DirectConditionGenerator::directWithSurfaceIntegration() - Data vector contains outdated values..." << std::endl;
-*/
+
return _direct_values;
}
diff --git a/MeshLib/MeshSurfaceExtraction.cpp b/MeshLib/MeshSurfaceExtraction.cpp
index d9af34468011223b0622c4c0dcf87d433518d0bf..7165d20665cc927ee7799994163d01428318be41 100644
--- a/MeshLib/MeshSurfaceExtraction.cpp
+++ b/MeshLib/MeshSurfaceExtraction.cpp
@@ -37,6 +37,7 @@ void MeshSurfaceExtraction::getSurfaceAreaForNodes(const MeshLib::Mesh &mesh, st
// for each node, a vector containing all the element idget every element
std::vector<MeshLib::Node*> nodes = mesh.getNodes();
const size_t nNodes ( mesh.getNNodes() );
+ node_area_vec.reserve(nNodes);
for (size_t n=0; n<nNodes; ++n)
{
double node_area (0);
@@ -101,14 +102,14 @@ MeshLib::Mesh* MeshSurfaceExtraction::getMeshSurface(const MeshLib::Mesh &mesh,
std::vector<std::size_t> id_map;
if (keep3dMeshIds)
{
- id_map.reserve(mesh.getNNodes());
- std::size_t idx(0);
- std::generate(id_map.begin(), id_map.end(), [&sfc_nodes, &idx](){ return sfc_nodes[idx++]->getID(); });
+ id_map.reserve(sfc_nodes.size());
+ for (auto node = sfc_nodes.cbegin(); node != sfc_nodes.cend(); ++node)
+ id_map.push_back((*node)->getID());
}
MeshLib::Mesh* result (new Mesh("SurfaceMesh", sfc_nodes, new_elements));
if (keep3dMeshIds)
- for (MeshLib::Node* node : sfc_nodes)
- node->setID(id_map[node->getID()]);
+ for (auto node = sfc_nodes.begin(); node != sfc_nodes.end(); ++node)
+ (*node)->setID(id_map[(*node)->getID()]);
return result;
}