Commit 7b6cd069 authored by Kirill Terekhov's avatar Kirill Terekhov
Browse files

Extend SearchKDTree class to search all cells intersecting segment

parent 239df1cf
......@@ -16,7 +16,7 @@ add_executable(Sew sew.cpp)
add_executable(Scale scale.cpp)
add_executable(Move move.cpp)
add_executable(Convert convert.cpp)
add_executable(SameMeshDifference difference_same.cpp)
add_executable(SameMeshDifference difference_same.cpp)
add_executable(MatchSameMeshDifference difference_same_match.cpp)
add_executable(MeshDifference difference_map.cpp)
add_executable(Kmeans kmeans.cpp)
......@@ -34,6 +34,7 @@ add_executable(split_faces split_faces.cpp)
add_executable(glue_faces glue_faces.cpp)
add_executable(check_collapse check_collapse.cpp)
add_executable(test_fracture test_fracture.cpp)
add_executable(segment_data segment_data.cpp)
add_library(FractureLib fracture.cpp fracture.h)
target_link_libraries(test_fracture inmost)
......@@ -269,8 +270,8 @@ if(USE_MPI)
set_target_properties(SameMeshDifference PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
endif()
endif(USE_MPI)
install(TARGETS SameMeshDifference EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
install(TARGETS SameMeshDifference EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
target_link_libraries(MatchSameMeshDifference inmost)
if(USE_MPI)
message("linking MatchSameMeshDifference with MPI")
......@@ -405,6 +406,16 @@ if(USE_MPI)
endif(USE_MPI)
install(TARGETS sinusoidal_mesh EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
target_link_libraries(segment_data inmost)
if(USE_MPI)
message("linking segment_data with MPI")
target_link_libraries(segment_data ${MPI_LIBRARIES})
if(MPI_LINK_FLAGS)
set_target_properties(segment_data PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
endif()
endif(USE_MPI)
install(TARGETS segment_data EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
set_property(TARGET FractureLib PROPERTY PUBLIC_HEADER "${CMAKE_CURRENT_SOURCE_DIR}/fracture.h")
......
#include "inmost.h"
using namespace INMOST;
int main(int argc, char ** argv)
{
if (argc < 9)
{
std::cout << "Usage: " << argv[0] << " mesh x1 y1 z1 x2 y2 z2 data:component" << std::endl;
return -1;
}
Mesh m;
m.Load(argv[1]);
double p1[3], p2[3];
p1[0] = atof(argv[2]);
p1[1] = atof(argv[3]);
p1[2] = atof(argv[4]);
p2[0] = atof(argv[5]);
p2[1] = atof(argv[6]);
p2[2] = atof(argv[7]);
std::string odata = std::string(argv[8]);
size_t pos = odata.find(":");
int comp = -1; //output all components
if( pos != std::string::npos )
{
comp = atoi(odata.substr(pos+1).c_str());
odata = odata.substr(0,pos);
}
/*
std::cout << "output ";
if( comp == -1 )
std::cout << "all components";
else
std::cout << "component " << comp;
std::cout << " of " << odata << std::endl;
*/
if( !m.HaveTag(odata) )
{
std::cout << "data " << odata << " not found in mesh" << std::endl;
return -1;
}
Tag otag = m.GetTag(odata);
if( comp >= 0 && otag.GetSize() < comp )
{
std::cout << "data " << odata << " size is " << otag.GetSize() << " which is less then " << comp << std::endl;
return -1;
}
SearchKDTree tree(&m);
ElementArray<Cell> cells;
tree.IntersectSegment(cells,p1,p2);
if( cells.empty() )
{
std::cout << "no cells found" << std::endl;
return -1;
}
std::cout << "x; y; z; ";
if( comp != -1 )
std::cout << odata << "[" << comp << "]; ";
else if( otag.GetSize() != ENUMUNDEF )
for(unsigned k = 0; k < otag.GetSize(); ++k) std::cout << odata << "[" << k << "]; ";
else std::cout << odata << "; ";
std::cout << std::endl;
for(ElementArray<Cell>::iterator it = cells.begin(); it != cells.end(); ++it)
{
double cnt[3];
it->Centroid(cnt);
std::cout << cnt[0] << "; " << cnt[1] << "; " << cnt[2] << "; ";
if( otag.GetDataType() == DATA_REAL )
{
Storage::real_array oarr = it->RealArray(otag);
if( comp == -1 ) for(Storage::real_array::iterator jt = oarr.begin(); jt != oarr.end(); ++jt) std::cout << *jt << "; ";
else if( comp < oarr.size() ) std::cout << oarr[comp] << "; ";
else std::cout << "NAN; ";
}
else if( otag.GetDataType() == DATA_VARIABLE )
{
Storage::var_array oarr = it->VariableArray(otag);
if( comp == -1 ) for(Storage::var_array::iterator jt = oarr.begin(); jt != oarr.end(); ++jt) std::cout << get_value(*jt) << "; ";
else if( comp < oarr.size() ) std::cout << get_value(oarr[comp]) << "; ";
else std::cout << "NAN; ";
}
else if( otag.GetDataType() == DATA_INTEGER )
{
Storage::integer_array oarr = it->IntegerArray(otag);
if( comp == -1 ) for(Storage::integer_array::iterator jt = oarr.begin(); jt != oarr.end(); ++jt) std::cout << *jt << "; ";
else if( comp < oarr.size() ) std::cout << oarr[comp] << "; ";
else std::cout << "NAN; ";
}
else if( otag.GetDataType() == DATA_BULK )
{
Storage::bulk_array oarr = it->BulkArray(otag);
if( comp == -1 ) for(Storage::bulk_array::iterator jt = oarr.begin(); jt != oarr.end(); ++jt) std::cout << *jt << "; ";
else if( comp < oarr.size() ) std::cout << oarr[comp] << "; ";
else std::cout << "NAN; ";
}
std::cout << std::endl;
}
return 0;
}
......@@ -3509,11 +3509,18 @@ namespace INMOST
Cell SubSearchCell(const Storage::real p[3], bool print);
void clear_children();
inline int segment_bbox(const Storage::real p1[3], const Storage::real p2[3]);
inline int segment_tri(const Storage::real tri[3][3], const Storage::real p1[3], const Storage::real p2[3]);
inline bool segment_face(const Element & f, const Storage::real p1[3], const Storage::real p2[3]);
inline bool segment_cell(const Element & c, const Storage::real p1[3], const Storage::real p2[3]);
void sub_intersect_segment(ElementArray<Element> & hits, MarkerType mrk, const Storage::real p1[3], const Storage::real p2[3]);
public:
SearchKDTree(Mesh * m);
SearchKDTree(Mesh * m, HandleType * _set, unsigned set_size);
~SearchKDTree();
Cell SearchCell(const Storage::real * point, bool print = false);
void IntersectSegment(ElementArray<Cell> & cells, const Storage::real p1[3], const Storage::real p2[3]);
};
......
......@@ -6,6 +6,17 @@
namespace INMOST
{
inline static void crossproduct(const Storage::real vecin1[3],const Storage::real vecin2[3], Storage::real vecout[3])
{
vecout[0] = vecin1[1]*vecin2[2] - vecin1[2]*vecin2[1];
vecout[1] = vecin1[2]*vecin2[0] - vecin1[0]*vecin2[2];
vecout[2] = vecin1[0]*vecin2[1] - vecin1[1]*vecin2[0];
}
inline static Storage::real dotproduct(const Storage::real * vecin1, const Storage::real * vecin2)
{
return vecin1[0]*vecin2[0]+vecin1[1]*vecin2[1]+vecin1[2]*vecin2[2];
}
template<typename bbox_type>
inline int SearchKDTree::bbox_point(const Storage::real p[3], const bbox_type bbox[6], bool print)
{
......@@ -292,6 +303,179 @@ namespace INMOST
else set = NULL;
}
inline int SearchKDTree::segment_tri(const Storage::real tri[3][3], const Storage::real p1[3], const Storage::real p2[3])
{
const Storage::real eps = 1.0e-7;
Storage::real a[3],b[3],c[3],n[3], ray[3], d, k, m, l;
Storage::real dot00,dot01,dot02,dot11,dot12,invdenom, uq,vq;
ray[0] = p2[0]-p1[0];
ray[1] = p2[1]-p1[1];
ray[2] = p2[2]-p1[2];
l = sqrt(dotproduct(ray,ray));
ray[0] /= l;
ray[1] /= l;
ray[2] /= l;
a[0] = tri[0][0] - tri[2][0];
a[1] = tri[0][1] - tri[2][1];
a[2] = tri[0][2] - tri[2][2];
b[0] = tri[1][0] - tri[2][0];
b[1] = tri[1][1] - tri[2][1];
b[2] = tri[1][2] - tri[2][2];
crossproduct(a,b,n);
d = -dotproduct(n,tri[0]);
m = dotproduct(n,ray);
if( fabs(m) < 1.0e-25 )
return 0;
k = -(d + dotproduct(n,p1))/m;
if( k < 0 )
return 0;
if( k > l )
return 0;
c[0] = p1[0] + k*ray[0] - tri[2][0];
c[1] = p1[1] + k*ray[1] - tri[2][1];
c[2] = p1[2] + k*ray[2] - tri[2][2];
dot00 = dotproduct(a,a);
dot01 = dotproduct(a,b);
dot02 = dotproduct(a,c);
dot11 = dotproduct(b,b);
dot12 = dotproduct(b,c);
invdenom = (dot00*dot11 - dot01*dot01);
uq = (dot11*dot02-dot01*dot12);
vq = (dot00*dot12-dot01*dot02);
if( fabs(invdenom) < 1.0e-25 && fabs(uq) > 0.0 && fabs(vq) > 0.0 )
return 0;
uq = uq/invdenom;
vq = vq/invdenom;
if( uq >= -eps && vq >= -eps && 1.0-(uq+vq) >= -eps )
return 1;
return 0;
}
inline int SearchKDTree::segment_bbox(const Storage::real p1[3], const Storage::real p2[3])
{
Storage::real tnear = -1.0e20, tfar = 1.0e20;
Storage::real t1,t2,c;
for(int i = 0; i < 3; i++)
{
if( fabs(p2[i]-p1[i]) < 1.0e-15 )
{
if( p1[i] < bbox[i*2] || p1[i] > bbox[i*2+1] )
return 0;
}
else
{
t1 = (bbox[i*2+0] - p1[i])/(p2[i]-p1[i]);
t2 = (bbox[i*2+1] - p1[i])/(p2[i]-p1[i]);
if( t1 > t2 )
{
c = t1;
t1 = t2;
t2 = c;
}
if( t1 > tnear ) tnear = t1;
if( t2 < tfar ) tfar = t2;
if( tnear > tfar ) return 0;
if( tfar < 0 ) return 0;
if( tnear > 1 ) return 0;
}
}
return 1;
}
inline bool SearchKDTree::segment_face(const Element & f, const Storage::real p1[3], const Storage::real p2[3])
{
Mesh * m = f->GetMeshLink();
Storage::real tri[3][3] = {{0,0,0},{0,0,0},{0,0,0}};
m->GetGeometricData(f->GetHandle(),CENTROID,tri[2]);
ElementArray<Node> nodes = f->getNodes();
m->GetGeometricData(nodes[0]->GetHandle(),CENTROID,tri[1]);
for(ElementArray<Node>::size_type k = 0; k < nodes.size(); ++k)
{
memcpy(tri[0],tri[1],sizeof(Storage::real)*3);
//m->GetGeometricData(nodes[k]->GetHandle(),CENTROID,tri[0]);
m->GetGeometricData(nodes[(k+1)%nodes.size()]->GetHandle(),CENTROID,tri[1]);
if( segment_tri(tri,p1,p2) ) return true;
}
return false;
}
inline bool SearchKDTree::segment_cell(const Element & c, const Storage::real p1[3], const Storage::real p2[3])
{
ElementArray<Face> faces = c->getFaces();
for(ElementArray<Face>::iterator it = faces.begin(); it != faces.end(); ++it)
if( segment_face(it->self(),p1,p2) ) return true;
return false;
}
void SearchKDTree::IntersectSegment(ElementArray<Cell> & cells, const Storage::real p1[3], const Storage::real p2[3])
{
ElementArray<Element> temp(m);
MarkerType mrk = m->CreateMarker();
sub_intersect_segment(temp, mrk, p1, p2);
m->RemMarkerArray(temp.data(), static_cast<Storage::enumerator>(temp.size()), mrk);
for (ElementArray<Element>::iterator it = temp.begin(); it != temp.end(); ++it)
{
if (it->GetElementType() == CELL && !it->GetMarker(mrk))
{
cells.push_back(it->getAsCell());
it->SetMarker(mrk);
}
else if (it->GetElementType() == FACE)
{
ElementArray<Cell> f_cells = it->getCells();
for (ElementArray<Cell>::iterator kt = f_cells.begin(); kt != f_cells.end(); ++kt) if (!kt->GetMarker(mrk))
{
cells.push_back(kt->self());
kt->SetMarker(mrk);
}
}
}
m->RemMarkerArray(cells.data(), static_cast<Storage::enumerator>(cells.size()), mrk);
m->ReleaseMarker(mrk);
}
void SearchKDTree::sub_intersect_segment(ElementArray<Element> & hits, MarkerType mrk, const Storage::real p1[3], const Storage::real p2[3])
{
if( size == 1 )
{
if( !m->GetMarker(set[0].e,mrk) )
{
Storage::integer edim = Element::GetGeometricDimension(m->GetGeometricType(set[0].e));
if( edim == 2 )
{
if( segment_face(Element(m,set[0].e),p1,p2) )
{
hits.push_back(set[0].e);
m->SetMarker(set[0].e,mrk);
}
}
else if( edim == 3 )
{
if( segment_cell(Element(m,set[0].e),p1,p2) )
{
hits.push_back(set[0].e);
m->SetMarker(set[0].e,mrk);
}
}
else
{
std::cout << __FILE__ << ":" << __LINE__ << " kd-tree structure is not suitable to intersect edges with segments" << std::endl;
exit(-1);
}
}
}
else
{
assert(size > 1);
if( segment_bbox(p1,p2) )
{
children[0].sub_intersect_segment(hits,mrk,p1,p2);
children[1].sub_intersect_segment(hits,mrk,p1,p2);
}
}
}
SearchKDTree::~SearchKDTree()
{
if( size )
......
Markdown is supported
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