Skip to content
Snippets Groups Projects
Commit ca3d8536 authored by Karsten Rink's avatar Karsten Rink
Browse files

added direct source terms with surface integration and shortened direct value method

parent 3b174211
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Gui/DataView/DirectConditionGenerator.cpp
+ 25
− 111
View file @ ca3d8536
......@@ -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;
}
......
MeshLib/MeshSurfaceExtraction.cpp
+ 6
− 5
View file @ ca3d8536
......@@ -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;
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment