Commit c3a618bf authored by Stenger, Florian's avatar Stenger, Florian
Browse files

bugfix: singularity-backgrid

parent 96fa0112
...@@ -7,7 +7,6 @@ ...@@ -7,7 +7,6 @@
#include <dune/curvedsurfacegrid/surfacedistance/vertexmap.hh> #include <dune/curvedsurfacegrid/surfacedistance/vertexmap.hh>
#include <dune/curvedsurfacegrid/surfacedistance/vtureader.hh> #include <dune/curvedsurfacegrid/surfacedistance/vtureader.hh>
#include <dune/grid/geometrygrid/coordfunction.hh> #include <dune/grid/geometrygrid/coordfunction.hh>
#include <cfloat>
#define COLLECTOR_BENCHMARK 0 #define COLLECTOR_BENCHMARK 0
...@@ -227,6 +226,18 @@ namespace Dune ...@@ -227,6 +226,18 @@ namespace Dune
#if COLLECTOR_BENCHMARK != 0 #if COLLECTOR_BENCHMARK != 0
void collectorBenchmark(bool fresh_cache){ void collectorBenchmark(bool fresh_cache){
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
//! $debug: output vertex-list
/*std::ofstream out("vertex_list.txt", std::ios_base::out | std::ios_base::binary | std::ios_base::trunc);
int v_cnt = inputs.size();
out.write(reinterpret_cast<char*>(&v_cnt), 4);
for(auto &input : inputs){
out.write(reinterpret_cast<char*>(&input[0]), 8);
out.write(reinterpret_cast<char*>(&input[1]), 8);
out.write(reinterpret_cast<char*>(&input[2]), 8);
}
out.close();*/
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
if(cached){ if(cached){
if(fresh_cache) cache->clear(); if(fresh_cache) cache->clear();
int uncached_cnt = 0; int uncached_cnt = 0;
......
...@@ -6,6 +6,7 @@ ...@@ -6,6 +6,7 @@
#include <vector> #include <vector>
#include <list> #include <list>
#include <stdint.h> #include <stdint.h>
#include <cfloat>
namespace Dune{ namespace Curved{ namespace Dune{ namespace Curved{
...@@ -237,8 +238,8 @@ struct Mesh{ ...@@ -237,8 +238,8 @@ struct Mesh{
//determine_coordinate_extrema //determine_coordinate_extrema
void determine_coordinate_extrema(Vertex &min_value, Vertex &max_value) const{ void determine_coordinate_extrema(Vertex &min_value, Vertex &max_value) const{
for(int p = 0; p < dim_of_world; ++p){ for(int p = 0; p < dim_of_world; ++p){
max_value[p] = -9e100; max_value[p] = -DBL_MAX;
min_value[p] = 9e100; min_value[p] = DBL_MAX;
} }
for(int i = 0, i_end = vertices.size(); i < i_end; ++i){ for(int i = 0, i_end = vertices.size(); i < i_end; ++i){
for(int p = 0; p < dim_of_world; ++p){ for(int p = 0; p < dim_of_world; ++p){
......
...@@ -2,6 +2,7 @@ ...@@ -2,6 +2,7 @@
#define DUNE_CURVED_SURFACE_GRID_VERTEXMAP_HH #define DUNE_CURVED_SURFACE_GRID_VERTEXMAP_HH
#include <dune/curvedsurfacegrid/surfacedistance/datastructures.hh> #include <dune/curvedsurfacegrid/surfacedistance/datastructures.hh>
#include <cmath>
#define VERTEXMAP_DEBUG_MODE 0 #define VERTEXMAP_DEBUG_MODE 0
...@@ -56,7 +57,7 @@ template<typename T> int VertexMap<T>::partition_threshold = 64; ...@@ -56,7 +57,7 @@ template<typename T> int VertexMap<T>::partition_threshold = 64;
template<typename T> int VertexMap<T>::depth_threshold_1d = 19; template<typename T> int VertexMap<T>::depth_threshold_1d = 19;
template<typename T> int VertexMap<T>::depth_threshold_2d = 10; template<typename T> int VertexMap<T>::depth_threshold_2d = 10;
template<typename T> int VertexMap<T>::depth_threshold_3d = 7; template<typename T> int VertexMap<T>::depth_threshold_3d = 7;
template<typename T> double VertexMap<T>::config_epsilon = epsilon; template<typename T> double VertexMap<T>::config_epsilon = Curved::epsilon;
//--------------------< destructor > //--------------------< destructor >
template<typename T> VertexMap<T>::~VertexMap(){ //$ACL_exclude template<typename T> VertexMap<T>::~VertexMap(){ //$ACL_exclude
...@@ -184,7 +185,7 @@ template<typename T> bool VertexMap<T>::find_1d(const Vertex &v, T &result){ ...@@ -184,7 +185,7 @@ template<typename T> bool VertexMap<T>::find_1d(const Vertex &v, T &result){
for(typename std::vector<std::pair<Vertex, T> >::iterator it = entries.begin(), for(typename std::vector<std::pair<Vertex, T> >::iterator it = entries.begin(),
end = entries.end(); it != end; ++it){ end = entries.end(); it != end; ++it){
Vertex &t = it->first; Vertex &t = it->first;
if(abs(t.c[0] - v.c[0]) <= config_epsilon){ if(fabs(t.c[0] - v.c[0]) <= config_epsilon){
result = it->second; result = it->second;
return true; return true;
} }
...@@ -202,7 +203,7 @@ template<typename T> bool VertexMap<T>::find_2d(const Vertex &v, T &result){ ...@@ -202,7 +203,7 @@ template<typename T> bool VertexMap<T>::find_2d(const Vertex &v, T &result){
for(typename std::vector<std::pair<Vertex, T> >::iterator it = entries.begin(), for(typename std::vector<std::pair<Vertex, T> >::iterator it = entries.begin(),
end = entries.end(); it != end; ++it){ end = entries.end(); it != end; ++it){
Vertex &t = it->first; Vertex &t = it->first;
if(abs(t.c[0] - v.c[0]) <= config_epsilon && abs(t.c[1] - v.c[1]) <= config_epsilon){ if(fabs(t.c[0] - v.c[0]) <= config_epsilon && fabs(t.c[1] - v.c[1]) <= config_epsilon){
result = it->second; result = it->second;
return true; return true;
} }
...@@ -226,8 +227,8 @@ template<typename T> bool VertexMap<T>::find_3d(const Vertex &v, T &result){ ...@@ -226,8 +227,8 @@ template<typename T> bool VertexMap<T>::find_3d(const Vertex &v, T &result){
for(typename std::vector<std::pair<Vertex, T> >::iterator it = entries.begin(), for(typename std::vector<std::pair<Vertex, T> >::iterator it = entries.begin(),
end = entries.end(); it != end; ++it){ end = entries.end(); it != end; ++it){
Vertex &t = it->first; Vertex &t = it->first;
if(abs(t.c[0] - v.c[0]) <= config_epsilon && abs(t.c[1] - v.c[1]) <= config_epsilon if(fabs(t.c[0] - v.c[0]) <= config_epsilon && fabs(t.c[1] - v.c[1]) <= config_epsilon
&& abs(t.c[2] - v.c[2]) <= config_epsilon){ && fabs(t.c[2] - v.c[2]) <= config_epsilon){
result = it->second; result = it->second;
return true; return true;
} }
......
...@@ -621,7 +621,7 @@ void Dune::Curved::read_vtu_file(std::string file, Mesh &m, std::vector<DataVect ...@@ -621,7 +621,7 @@ void Dune::Curved::read_vtu_file(std::string file, Mesh &m, std::vector<DataVect
vvr->initialize(); vvr->initialize();
if(datatypes_point_data[e] == "Float32"){ if(datatypes_point_data[e] == "Float32"){
unsigned entry_count = vvr->total_size / 4; unsigned entry_count = vvr->total_size / 4;
if(int(entry_count / component_counts_point_data[e]) != m.vertices.size()){ if(entry_count / component_counts_point_data[e] != m.vertices.size()){
error("Number of vertices does not match length of appended data for vertex-vector " error("Number of vertices does not match length of appended data for vertex-vector "
+ ft(e) + "."); + ft(e) + ".");
} }
...@@ -632,7 +632,7 @@ void Dune::Curved::read_vtu_file(std::string file, Mesh &m, std::vector<DataVect ...@@ -632,7 +632,7 @@ void Dune::Curved::read_vtu_file(std::string file, Mesh &m, std::vector<DataVect
} }
}else if(datatypes_point_data[e] == "Float64"){ }else if(datatypes_point_data[e] == "Float64"){
unsigned entry_count = vvr->total_size / 8; unsigned entry_count = vvr->total_size / 8;
if(int(entry_count / component_counts_point_data[e]) != m.vertices.size()){ if(entry_count / component_counts_point_data[e] != m.vertices.size()){
error("Number of vertices does not match length of appended data for vertex-vector " error("Number of vertices does not match length of appended data for vertex-vector "
+ ft(e) + "."); + ft(e) + ".");
} }
...@@ -643,7 +643,7 @@ void Dune::Curved::read_vtu_file(std::string file, Mesh &m, std::vector<DataVect ...@@ -643,7 +643,7 @@ void Dune::Curved::read_vtu_file(std::string file, Mesh &m, std::vector<DataVect
} }
}else if(datatypes_point_data[e] == "Int64"){ }else if(datatypes_point_data[e] == "Int64"){
unsigned entry_count = vvr->total_size / 8; unsigned entry_count = vvr->total_size / 8;
if(int(entry_count / component_counts_point_data[e]) != m.vertices.size()){ if(entry_count / component_counts_point_data[e] != m.vertices.size()){
error("Number of vertices does not match length of appended data for vertex-vector " error("Number of vertices does not match length of appended data for vertex-vector "
+ ft(e) + "."); + ft(e) + ".");
} }
...@@ -664,7 +664,7 @@ void Dune::Curved::read_vtu_file(std::string file, Mesh &m, std::vector<DataVect ...@@ -664,7 +664,7 @@ void Dune::Curved::read_vtu_file(std::string file, Mesh &m, std::vector<DataVect
vvr->initialize(); vvr->initialize();
if(datatypes_cell_data[e] == "Float32"){ if(datatypes_cell_data[e] == "Float32"){
unsigned entry_count = vvr->total_size / 4; unsigned entry_count = vvr->total_size / 4;
if(int(entry_count / component_counts_cell_data[e]) != m.elements.size()){ if(entry_count / component_counts_cell_data[e] != m.elements.size()){
error("Number of elements does not match length of appended data for element-vector " error("Number of elements does not match length of appended data for element-vector "
+ ft(e) + "."); + ft(e) + ".");
} }
...@@ -675,7 +675,7 @@ void Dune::Curved::read_vtu_file(std::string file, Mesh &m, std::vector<DataVect ...@@ -675,7 +675,7 @@ void Dune::Curved::read_vtu_file(std::string file, Mesh &m, std::vector<DataVect
} }
}else if(datatypes_cell_data[e] == "Float64"){ }else if(datatypes_cell_data[e] == "Float64"){
unsigned entry_count = vvr->total_size / 8; unsigned entry_count = vvr->total_size / 8;
if(int(entry_count / component_counts_cell_data[e]) != m.elements.size()){ if(entry_count / component_counts_cell_data[e] != m.elements.size()){
error("Number of elements does not match length of appended data for element-vector " error("Number of elements does not match length of appended data for element-vector "
+ ft(e) + "."); + ft(e) + ".");
} }
...@@ -686,7 +686,7 @@ void Dune::Curved::read_vtu_file(std::string file, Mesh &m, std::vector<DataVect ...@@ -686,7 +686,7 @@ void Dune::Curved::read_vtu_file(std::string file, Mesh &m, std::vector<DataVect
} }
}else if(datatypes_cell_data[e] == "Int64"){ }else if(datatypes_cell_data[e] == "Int64"){
unsigned entry_count = vvr->total_size / 8; unsigned entry_count = vvr->total_size / 8;
if(int(entry_count / component_counts_cell_data[e]) != m.elements.size()){ if(entry_count / component_counts_cell_data[e] != m.elements.size()){
error("Number of elements does not match length of appended data for element-vector " error("Number of elements does not match length of appended data for element-vector "
+ ft(e) + "."); + ft(e) + ".");
} }
......
...@@ -2,7 +2,7 @@ ...@@ -2,7 +2,7 @@
// vi: set et ts=4 sw=2 sts=2: // vi: set et ts=4 sw=2 sts=2:
#ifdef HAVE_CONFIG_H #ifdef HAVE_CONFIG_H
# include "config.h" #include "config.h"
#endif #endif
#include <iostream> #include <iostream>
#include <dune/common/parallel/mpihelper.hh> // An initializer of MPI #include <dune/common/parallel/mpihelper.hh> // An initializer of MPI
...@@ -45,7 +45,7 @@ int main(int argc, char** argv) ...@@ -45,7 +45,7 @@ int main(int argc, char** argv)
//case1: projection to actual sphere //case1: projection to actual sphere
{ {
SphereCoordFunction sphereCF(/*center=*/{0.0, 0.0, 0.0}, /*radius=*/1.0); SphereCoordFunction sphereCF(/*center=*/{0.0, 0.0, 0.0}, /*radius=*/1.0);
CurvedSurfaceGrid grid(*hostGrid, sphereCF, std::integral_constant<int, ip_order>{}); CurvedSurfaceGrid grid(*hostGrid, std::ref(sphereCF), std::integral_constant<int, ip_order>{});
auto gridView = grid.leafGridView(); auto gridView = grid.leafGridView();
using GridView = decltype(gridView); using GridView = decltype(gridView);
...@@ -53,6 +53,10 @@ int main(int argc, char** argv) ...@@ -53,6 +53,10 @@ int main(int argc, char** argv)
vtkWriter.write("sphere1_real_elements"); vtkWriter.write("sphere1_real_elements");
SubsamplingVTKWriter<GridView> ssvtkWriter(gridView, refinementIntervals(ip_order)); SubsamplingVTKWriter<GridView> ssvtkWriter(gridView, refinementIntervals(ip_order));
ssvtkWriter.write("sphere1_fake_elements" + Curved::ft(ip_order, 2)); ssvtkWriter.write("sphere1_fake_elements" + Curved::ft(ip_order, 2));
#if COLLECTOR_BENCHMARK != 0
std::cout << "benchmark-result for SphereCoordFunction:" << std::endl;
sphereCF.collectorBenchmark(true);
#endif
} }
//case2: projection to a surface-mesh of a sphere //case2: projection to a surface-mesh of a sphere
...@@ -66,6 +70,10 @@ int main(int argc, char** argv) ...@@ -66,6 +70,10 @@ int main(int argc, char** argv)
vtkWriter.write("sphere2_real_elements"); vtkWriter.write("sphere2_real_elements");
SubsamplingVTKWriter<GridView> ssvtkWriter(gridView, refinementIntervals(ip_order)); SubsamplingVTKWriter<GridView> ssvtkWriter(gridView, refinementIntervals(ip_order));
ssvtkWriter.write("sphere2_fake_elements" + Curved::ft(ip_order, 2)); ssvtkWriter.write("sphere2_fake_elements" + Curved::ft(ip_order, 2));
#if COLLECTOR_BENCHMARK != 0
std::cout << "benchmark-result for SurfaceDistanceCoordFunction:" << std::endl;
surfdistCF.collectorBenchmark(true);
#endif
} }
return 0; return 0;
......
Supports Markdown
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