Commit bab94ecf by Kirill Terekhov

Update to slice by function grid tool

parent 1b73e719
 ... ... @@ -2,25 +2,61 @@ using namespace INMOST; double func(double x, double y, double z, int n) { if( n == 0 ) return sqrt(x*x+y*y)-0.25; else if( n == 1 ) return -(sqrt((x-0.5)*(x-0.5)+(y-0.5)*(y-0.5)+(z-0.5)*(z-0.5))-0.5); else if( n == 2 ) { if( x > 0.5 ) return -(sqrt((x-0.7)*(x-0.7)+(y-0.2)*(y-0.2)+(z-0.4)*(z-0.4))-0.4); else return sqrt((x-0.3)*(x-0.3)+(y-0.4)*(y-0.4)+(z-0.5)*(z-0.5))-0.3; } else if( n == 3 ) return y-(4*x*x*x*x-2*x*x+0.5)+z*z*z; else if( n == 4 ) return y-(10*x*x*x*x-8*x*x*x-5*x*x+0.2+4*x)+z*z*z; else if( n == 5 ) { double Lx = 0.2, Rx = 0.4, Ly = 0.1, Ry = 0.3; if (x > Rx){ if (y > Ry) return -sqrt( (x-Rx)*(x-Rx) + (y-Ry)*(y-Ry) ); else return -(x-Rx); } if (x < Lx){ if (y < Ly) return -sqrt( (x-Lx)*(x-Lx) + (y-Ly)*(y-Ly) ); else return (x - Lx); } if (y > Ry) return Ry - y; if (y < Ly) return y - Ly; return fmin( fmin(x-Lx,Rx-x), fmin(y-Ly, Ry-y) ); } return 1; } void search(double r0, double r1, double c0[3], double c1[3], double p[3],int type) double search(double r0, double r1, double c0[3], double c1[3], double p[3],int type, bool binary = false) { double rp; double rp = 1.0e20, rp_min = 1.0e20, p_min[3]; int iters = 0; do { p[0] = (r0*c1[0] - r1*c0[0])/(r0-r1); p[1] = (r0*c1[1] - r1*c0[1])/(r0-r1); p[2] = (r0*c1[2] - r1*c0[2])/(r0-r1); double m1 = r0/(r0-r1), m0; if( m1 < 0.0 || m1 > 1.0 || binary ) m1 = 0.5; m0 = 1.0-m1; p[0] = m1*c1[0] + m0*c0[0]; p[1] = m1*c1[1] + m0*c0[1]; p[2] = m1*c1[2] + m0*c0[2]; rp = func(p[0],p[1],p[2],type); if( rp_min > rp ) { p_min[0] = p[0]; p_min[1] = p[1]; p_min[2] = p[2]; rp_min = rp; } if( rp*r0 < 0.0 ) { c1[0] = p[0]; ... ... @@ -35,7 +71,18 @@ void search(double r0, double r1, double c0[3], double c1[3], double p[3],int ty c0[2] = p[2]; r0 = rp; } } while( fabs(rp) > 1.0e-9 ); iters++; } while( fabs(rp) > 1.0e-9 && iters < 100 ); if( rp > rp_min ) { p[0] = p_min[0]; p[1] = p_min[1]; p[2] = p_min[2]; rp = rp_min; } return rp; } int main(int argc, char ** argv) ... ... @@ -58,13 +105,14 @@ int main(int argc, char ** argv) m.Load(argv[1]); m.SetTopologyCheck(NEED_TEST_CLOSURE|PROHIBIT_MULTILINE|PROHIBIT_MULTIPOLYGON|GRID_CONFORMITY|DEGENERATE_EDGE|DEGENERATE_FACE|DEGENERATE_CELL | FACE_EDGES_ORDER); //m.RemTopologyCheck(THROW_EXCEPTION); Tag material = m.CreateTag("MATERIAL",DATA_INTEGER,CELL|FACE|EDGE|NODE,NONE,1); TagInteger material = m.CreateTag("MATERIAL",DATA_INTEGER,CELL|FACE|EDGE|NODE,NONE,1); Tag sliced = m.CreateTag("SLICED",DATA_BULK,FACE|EDGE|NODE,FACE|EDGE|NODE,1); std::cout << "Cells: " << m.NumberOfCells() << std::endl; std::cout << "Faces: " << m.NumberOfFaces() << std::endl; MarkerType original = m.CreateMarker(); for(Mesh::iteratorNode it = m.BeginNode(); it != m.EndNode(); ++it) it->SetMarker(original); for(Mesh::iteratorEdge it = m.BeginEdge(); it != m.EndEdge(); ++it) it->SetMarker(original); MarkerType slice = m.CreateMarker(); ... ... @@ -77,46 +125,79 @@ int main(int argc, char ** argv) Storage::real_array c1 = it->getEnd()->Coords(); double r0 = func(c0[0],c0[1],c0[2],type); double r1 = func(c1[0],c1[1],c1[2],type); it->getBeg()->Integer(material) = (r0 <= 0)? 0 : 1; it->getEnd()->Integer(material) = (r1 <= 0)? 0 : 1; material[it->getBeg()] = (r0 <= 0)? 0 : 1; material[it->getEnd()] = (r1 <= 0)? 0 : 1; //std::cout << "r0 " << r0 << " r1 " << r1 << std::endl; if( r0*r1 < -1.0e-12 ) { pc0[0] = c0[0], pc0[1] = c0[1], pc0[2] = c0[2]; pc1[0] = c1[0], pc1[1] = c1[1], pc1[2] = c1[2]; search(r0,r1,pc0,pc1,p,type); double rp = search(r0,r1,pc0,pc1,p,type); if( rp > 1.0e-3 ) //cannot find intersection { rp = search(r0,r1,pc0,pc1,p,type,true); p[0] = c0[0]*0.5+c1[0]*0.5; p[1] = c0[1]*0.5+c1[1]*0.5; p[2] = c0[2]*0.5+c1[2]*0.5; } //p[0] = (r0*c1[0] - r1*c0[0])/(r0-r1); //p[1] = (r0*c1[1] - r1*c0[1])/(r0-r1); //p[2] = (r0*c1[2] - r1*c0[2])/(r0-r1); //std::cout << "p " << p[0] << " " << p[1] << " " << p[2] << std::endl; Node n = m.CreateNode(p); n->Integer(material) = 2; n.Bulk(sliced) = 1; n.SetMarker(slice); bool was_sliced = it->HaveData(sliced) ? true : false; ElementArray ret = Edge::SplitEdge(it->self(),ElementArray(&m,1,n.GetHandle()),0); ret.SetMarker(mrk); if( was_sliced ) for(int q = 0; q < ret.size(); ++q) ret[q]->Bulk(sliced) = 1; nslice++; //distance to the corners double l0 = 0, l1 = 0, l; for(int r = 0; r < 3; ++r) { l0 += (p[r]-c0[r])*(p[r]-c0[r]); l1 += (p[r]-c1[r])*(p[r]-c1[r]); } l0 = sqrt(l0); l1 = sqrt(l1); l = l0+l1; if( l0 < 1.0e-5*l ) { material[it->getBeg()] = 2; it->getBeg()->SetMarker(slice); nmark++; } else if( l1 < 1.0e-5*l ) { material[it->getEnd()] = 2; it->getEnd()->SetMarker(slice); nmark++; } else { Node n = m.CreateNode(p); material[n] = 2; n.Bulk(sliced) = 1; n.SetMarker(slice); bool was_sliced = it->HaveData(sliced) ? true : false; ElementArray ret = Edge::SplitEdge(it->self(),ElementArray(&m,1,n.GetHandle()),0); ret.SetMarker(mrk); if( was_sliced ) for(int q = 0; q < ret.size(); ++q) ret[q]->Bulk(sliced) = 1; nslice++; } } else { if( fabs(r0) < 1.0e-6 ) { it->getBeg()->Integer(material) = 2; material[it->getBeg()] = 2; it->getBeg()->SetMarker(slice); nmark++; } if( fabs(r1) < 1.0e-6 ) { it->getEnd()->Integer(material) = 2; material[it->getEnd()] = 2; it->getEnd()->SetMarker(slice); nmark++; } } } m.ReleaseMarker(mrk,EDGE); std::cout << "sliced edges: " << nslice << " marked nodes: " << nmark << std::endl; ... ... @@ -128,7 +209,7 @@ int main(int argc, char ** argv) int tot = 0; for(Mesh::iteratorNode it = m.BeginNode(); it != m.EndNode(); ++it) { mat[it->Integer(material)]++; mat[material[*it]]++; tot++; } std::cout << "node materials, 0: " << mat[0] << " 1: " << mat[1] << " 2: " << mat[2] << std::endl; ... ... @@ -139,52 +220,228 @@ int main(int argc, char ** argv) for(Mesh::iteratorEdge it = m.BeginEdge(); it != m.EndEdge(); ++it) { int mat[3] = {0,0,0}; mat[it->getBeg()->Integer(material)]++; mat[it->getEnd()->Integer(material)]++; mat[material[it->getBeg()]]++; mat[material[it->getEnd()]]++; if( !(mat[0] == 0 || mat[1] == 0) ) { it->Integer(material) = 2; material[*it] = 2; std::cout << "oops, materials for edge nodes were not split, 0: " << mat[0] << " ,1: " << mat[1] << " ,2: " << mat[2] << std::endl; } else if( mat[0] != 0 ) it->Integer(material) = 0; else if( mat[1] != 0 ) it->Integer(material) = 1; else it->Integer(material) = 2; else if( mat[0] != 0 ) material[*it] = 0; else if( mat[1] != 0 ) material[*it] = 1; else material[*it] = 2; } nslice = 0; for(Mesh::iteratorFace it = m.BeginFace(); it != m.EndFace(); ++it) for(Mesh::iteratorFace it = m.BeginFace(); it != m.EndFace(); ++it) if( !it->GetMarker(mrk) ) { ElementArray nodes = it->getNodes(slice); //those nodes should be ordered so that each pair forms an edge if( nodes.size() > 1 ) // if there is 1, then only one vertex touches the plane { int nsliced = 0; for(int q = 0; q < (int)nodes.size(); ++q) if( !nodes[q].GetMarker(original) ) nsliced++; //int nsliced = 0; //for(int q = 0; q < (int)nodes.size(); ++q) if( !nodes[q].GetMarker(original) ) nsliced++; //if there is more then two, then original face is non-convex if( nsliced && nodes.size() > 2 ) std::cout << "Looks like face " << it->LocalID() << " is nonconvex, there is " << nodes.size() << " nodes, out of them " << nsliced << " new cuts on edge" << std::endl; else if( nodes.size() == 2 ) //if( nsliced && nodes.size() > 2 ) std::cout << "Looks like face " << it->LocalID() << " is nonconvex, there is " << it->nbAdjElements(NODE) << " nodes, out of them " << nsliced << " new cuts on face" << " slice " << slice << " original " << original << std::endl; //else if( nodes.size() == 2 ) { Edge e = m.CreateEdge(nodes).first; e.Integer(material) = 2; //on slice material[e] = 2; //on slice e.Bulk(sliced) = 1; e.SetMarker(slice); bool was_sliced = it->HaveData(sliced) ? true : false; ElementArray ret = Face::SplitFace(it->self(),ElementArray(&m,1,e.GetHandle()),0); ret.SetMarker(mrk); if( was_sliced ) for(int q = 0; q < ret.size(); ++q) ret[q]->Bulk(sliced) = 1; nslice++; } else if( it->nbAdjElements(NODE,slice) != it->nbAdjElements(NODE) ) //not entire face is sliced { //std::cout << "Sliced nodes " << nodes.size() << " total nodes: "; nodes = it->getNodes(); //std::cout << nodes.size() << std::endl; for(int q = 0; q < (int)nodes.size(); ++q) { Storage::real_array c0 = nodes[q].Coords(); double r0 = func(c0[0],c0[1],c0[2],type); //std::cout << "NODE:" << nodes[q].LocalID() << " " << (nodes[q]->GetMarker(slice)?"":"not ") << "slice " << (nodes[q]->GetMarker(original)?"":"not ") << "original r="< fedges = it->getEdges(); for(int q = 0; q < (int)fedges.size(); ++q) { //std::cout << "EDGE:" << fedges[q].LocalID() << " " << fedges[q].getBeg().LocalID() << "<->" << fedges[q].getEnd().LocalID() << " " << (fedges[q]->GetMarker(slice)?"":"not ") << "slice material " << material[fedges[q]] << std::endl; } double c0[3],c1[3],pc0[3],pc1[3],p[3]; it->Centroid(c0); double r0 = func(c0[0],c0[1],c0[2],type); int material0 = (r0 <= 0)? 0 : 1; bool s0 = false; if( fabs(r0) < 1.0e-6 ) s0 = true; //std::cout << "Centernode " << (s0?"sliced":"") << " r=" << r0 << std::endl; Node centernode = InvalidNode(); ElementArray split_edges(&m); ElementArray cutnodes(&m,nodes.size()), edge_nodes(&m,2); //calculate nodes that cut along the edges connecting centernode for(int q = 0; q < (int)nodes.size(); ++q) if( !nodes[q].GetMarker(slice) ) { nodes[q].Centroid(c1); double r1 = func(c1[0],c1[1],c1[2],type); //std::cout << "NODE:" << nodes[q].LocalID() << " r0 " << r0 << " r1 " << r1 << " r0*r1 " << r0*r1 << std::endl; if( r0*r1 < -1.0e-12 ) { pc0[0] = c0[0], pc0[1] = c0[1], pc0[2] = c0[2]; pc1[0] = c1[0], pc1[1] = c1[1], pc1[2] = c1[2]; double rp = search(r0,r1,pc0,pc1,p,type); if( rp > 1.0e-3 ) { //std::cout << "inaccurate search " << rp; rp = search(r0,r1,pc0,pc1,p,type,true); //std::cout << " binary " << rp << std::endl; p[0] = c0[0]*0.5+c1[0]*0.5; p[1] = c0[1]*0.5+c1[1]*0.5; p[2] = c0[2]*0.5+c1[2]*0.5; } //distance to the center node double l0 = 0, l1 = 0, l; for(int r = 0; r < 3; ++r) { l0 += (p[r]-c0[r])*(p[r]-c0[r]); l1 += (p[r]-c1[r])*(p[r]-c1[r]); } l0 = sqrt(l0); l1 = sqrt(l1); l = l0+l1; if( l0 < 1.0e-3*l ) //edge goes through centernode { if( !centernode.isValid() ) centernode = m.CreateNode(c0); cutnodes[q] = centernode; //std::cout << "selected centernode " << std::endl; } else if( l1 > 1.0e-3*l ) { cutnodes[q] = m.CreateNode(p); //std::cout << "created new node " << std::endl; } else { material[nodes[q]] = 2; nodes[q].SetMarker(slice); } //else std::cout << "use old node " << std::endl; } else if( s0 ) { if( !centernode.isValid() ) centernode = m.CreateNode(c0); cutnodes[q] = centernode; } } for(int q = 0; q < (int)cutnodes.size(); ++q) if( cutnodes[q].isValid() ) { //std::cout << "New cut node " << cutnodes[q].LocalID() << " at " << q << " on line with " << nodes[q].LocalID() << " r=" << func(cutnodes[q].Coords()[0],cutnodes[q].Coords()[1],cutnodes[q].Coords()[2],type) << std::endl; Node n = cutnodes[q]; material[n] = 2; n.Bulk(sliced) = 1; n.SetMarker(slice); } //go over triangles and figure out the path of the cut for(int q = 0; q < (int)nodes.size(); ++q) { int i1 = q; int i2 = (q+1)%nodes.size(); Node n1 = nodes[i1]; Node n2 = nodes[i2]; bool s1 = n1->GetMarker(slice); bool s2 = n2->GetMarker(slice); if( s1 && s2 ) { //cut passing through edge //skip this } else if( s2 ) { //check cut on oposite edge if( cutnodes[i1].isValid() ) { edge_nodes[0] = n2; edge_nodes[1] = cutnodes[i1]; Edge e = m.CreateEdge(edge_nodes).first; material[e] = 2; //on slice e.Bulk(sliced) = 1; e.SetMarker(slice); split_edges.push_back(e); } } else if( s1 ) { //check cut on oposite edge if( cutnodes[i2].isValid() ) { edge_nodes[0] = n1; edge_nodes[1] = cutnodes[i2]; Edge e = m.CreateEdge(edge_nodes).first; material[e] = 2; //on slice e.Bulk(sliced) = 1; e.SetMarker(slice); split_edges.push_back(e); } } else if( cutnodes[i1].isValid() && cutnodes[i2].isValid() ) { edge_nodes[0] = cutnodes[i1]; edge_nodes[1] = cutnodes[i2]; Edge e = m.CreateEdge(edge_nodes).first; material[e] = 2; //on slice e.Bulk(sliced) = 1; e.SetMarker(slice); split_edges.push_back(e); } } //split face with multiple edges if( !split_edges.empty() ) { //std::cout << "Split edges " << split_edges.size() << std::endl; //for(int q = 0; q < (int)split_edges.size(); ++q) // std::cout << split_edges[q].getBeg().LocalID() << "<->" << split_edges[q].getEnd().LocalID() << std::endl; bool was_sliced = it->HaveData(sliced) ? true : false; ElementArray ret = Face::SplitFace(it->self(),split_edges,0); ret.SetMarker(mrk); //std::cout << "New faces: " << ret.size() << ":"; for(int q = 0; q < ret.size(); ++q) { int mat[3] = {0,0,0}; ElementArray fe = ret[q].getEdges(); for( int l = 0; l < fe.size(); ++l) mat[material[fe[l]]]++; //std::cout << " " << ret[q].LocalID() << "[" << mat[0] << "," << mat[1] << "," << mat[2] << "]"; } //std::cout << std::endl; if( was_sliced ) for(int q = 0; q < ret.size(); ++q) ret[q]->Bulk(sliced) = 1; nslice++; } //else std::cout << "No split edges " << std::endl; } //else std::cout << "No new edges on face " << it->LocalID() << std::endl; } //else std::cout << "Only one adjacent slice node, face " << it->LocalID() << std::endl; } m.ReleaseMarker(original,NODE); nmark = 0; for(Mesh::iteratorEdge it = m.BeginEdge(); it != m.EndEdge(); ++it) if( !it->GetMarker(slice) && it->getBeg()->GetMarker(slice) && it->getEnd()->GetMarker(slice) ) { if( it->Integer(material) != 2 ) std::cout << "Edge supposed to get material 2, but have " << it->Integer(material) << std::endl; if( material[*it] != 2 ) std::cout << "Edge supposed to get material 2, but have " << material[*it] << std::endl; it->SetMarker(slice); nmark++; } ... ... @@ -199,7 +456,7 @@ int main(int argc, char ** argv) int tot = 0; for(Mesh::iteratorEdge it = m.BeginEdge(); it != m.EndEdge(); ++it) { mat[it->Integer(material)]++; mat[material[*it]]++; tot++; } std::cout << "edge materials, 0: " << mat[0] << " 1: " << mat[1] << " 2: " << mat[2] << std::endl; ... ... @@ -210,49 +467,517 @@ int main(int argc, char ** argv) int mat[3] = {0,0,0}; ElementArray edges = it->getEdges(); for(int k = 0; k < (int)edges.size(); ++k) mat[edges[k]->Integer(material)]++; mat[material[edges[k]]]++; if( !(mat[0] == 0 || mat[1] == 0) ) { it->Integer(material) = 2; std::cout << "oops, materials for face edges were not split, 0: " << mat[0] << " ,1: " << mat[1] << " ,2: " << mat[2] << std::endl; material[*it] = 2; std::cout << "oops, materials for face " << it->LocalID() << " edges were not split, 0: " << mat[0] << " ,1: " << mat[1] << " ,2: " << mat[2] << std::endl; } else if( mat[0] != 0 ) it->Integer(material) = 0; else if( mat[1] != 0 ) it->Integer(material) = 1; else it->Integer(material) = 2; else if( mat[0] != 0 ) material[*it] = 0; else if( mat[1] != 0 ) material[*it] = 1; else material[*it] = 2; } nslice = 0; TagInteger indx = m.CreateTag("TEMP_INDX",DATA_INTEGER,NODE|EDGE,NONE,1); MarkerType visited = m.CreateMarker(); for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it) MarkerType cmrk = m.CreateMarker(); MarkerType isolate = m.CreateMarker(); for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it) if( !it->GetMarker(mrk) ) { ElementArray edges = it->getEdges(slice); if( edges.size() >= 3 ) //these should form a triangle if( !edges.empty() ) //these should form a triangle { //order edges ElementArray order_edges(&m); order_edges.push_back(edges[0]); order_edges.SetMarker(visited); while(order_edges.size() != edges.size() ) //check edges form a simple loops, each node should be visited twice ElementArray split_faces(&m); std::map visit_count; for(ElementArray::iterator jt = edges.begin(); jt != edges.end(); ++jt) { visit_count[jt->getBeg()]++; visit_count[jt->getEnd()]++; } bool simple = true; for(std::map::iterator jt = visit_count.begin(); jt != visit_count.end(); ++jt) simple &= (jt->second == 2); if( simple ) { //gather loop by loop and create faces int loop_cnt = 0; //order edges ElementArray order_edges(&m); Node last; int nvisited = 0; bool found = false; for(int k = 0; k < edges.size(); ++k) if( !edges[k]->GetMarker(original) ) { nvisited++; last = InvalidNode(); order_edges.push_back(edges[k]); order_edges.SetMarker(visited); found = true; break; } if( !found ) continue; while(nvisited != edges.size() ) { found = false; for(int k = 0; k < edges.size(); ++k) if( !edges[k]->GetMarker(visited) ) { bool match = false; if( last.isValid() ) { if( edges[k]->getBeg() == last ) { match = true; last = edges[k]->getEnd(); } else if( edges[k]->getEnd() == last ) { match = true; last = edges[k]->getBeg(); } } else { if( edges[k]->getBeg() == order_edges.back()->getBeg() || edges[k]->getBeg() == order_edges.back()->getEnd() ) { match = true; last = edges[k]->getEnd(); } else if( edges[k]->getEnd() == order_edges.back()->getBeg() || edges[k]->getEnd() == order_edges.back()->getEnd() ) { match = true; last = edges[k]->getBeg(); } } if( match ) { nvisited++; order_edges.push_back(edges[k]); order_edges.back().SetMarker(visited); found = true; } } if( !found || nvisited == edges.size() ) { if( !order_edges.empty() ) { //std::cout << "New loop " << ++loop_cnt << ": " << order_edges.size() << " total edges " << edges.size() << std::endl; //for(int k = 0; k < order_edges.size(); ++k) // std::cout << order_edges[k].getBeg().LocalID() << "<->" << order_edges[k].getEnd().LocalID() << std::endl; //std::cout << "All:" << std::endl; //for(int k = 0; k < edges.size(); ++k) // std::cout << edges[k].getBeg().LocalID() << "<->" << edges[k].getEnd().LocalID() << " " << (edges[k].GetMarker(visited) ? "visited":"") << " " << (edges[k].GetMarker(original) ? "original":"")<< std::endl; if( order_edges.size() > 2 ) { std::pair f = m.CreateFace(order_edges); material[f.first] = 2; f.first->Bulk(sliced) = 1; split_faces.push_back(f.first); } order_edges.clear(); found = false; for(int k = 0; k < edges.size(); ++k) if( !edges[k]->GetMarker(visited) && !edges[k]->GetMarker(original)) { nvisited++; last = InvalidNode(); order_edges.push_back(edges[k]); order_edges.back().SetMarker(visited); found = true; break; } if( !found ) break; } else break; } } //if( !order_edges.empty() ) std::cout << "size: " << order_edges.size() << std::endl; edges.RemMarker(visited); } else { //edges do not breakup into simple loops, have to run complicated algorithm //split into pyramids joining faces ElementArray cfaces = it->getFaces(); ElementArray cedges = it->getEdges(); ElementArray cnodes = it->getNodes(); ElementArray cutnodes(&m,cnodes.size()); ElementArray cutedges(&m,cedges.size()); ElementArray edge_nodes(&m,2); for(int k = 0; k < (int)cnodes.size(); ++k) indx[cnodes[k]] = k; for(int k = 0; k < (int)cedges.size(); ++k) indx[cedges[k]] = k; double c0[3],c1[3],pc0[3],pc1[3],p[3]; it->Centroid(c0); double r0 = func(c0[0],c0[1],c0[2],type); int material0 = (r0 <= 0)? 0 : 1; bool s0 = false; if( fabs(r0) < 1.0e-6 ) s0 = true; //std::cout << "Number of cut edges: " << edges.size() << std::endl; //for(std::map::iterator jt = visit_count.begin(); jt != visit_count.end(); ++jt) // std::cout << "NODE:" << jt->first.LocalID() << " visited " << jt->second << " times" << std::endl; //std::cout << "Centernode " << (s0?"sliced":"") << " r=" << r0 << std::endl; Node centernode = InvalidNode(); //calculate nodes that cut along the edges connecting centernode for(int q = 0; q < (int)cnodes.size(); ++q) { if( !cnodes[q].GetMarker(slice) ) { cnodes[q].Centroid(c1); double r1 = func(c1[0],c1[1],c1[2],type); //std::cout << "NODE:" << cnodes[q].LocalID() << " r0 " << r0 << " r1 " << r1 << " r0*r1 " << r0*r1 << " " << (cnodes[q].GetMarker(slice)?"":"not ") << "sliced" << std::endl; if( r0*r1 < -1.0e-12 ) { pc0[0] = c0[0], pc0[1] = c0[1], pc0[2] = c0[2]; pc1[0] = c1[0], pc1[1] = c1[1], pc1[2] = c1[2]; double rp = search(r0,r1,pc0,pc1,p,type); if( rp > 1.0e-3 ) { //std::cout << "inaccurate search " << rp; rp = search(r0,r1,pc0,pc1,p,type,true); //std::cout << " binary " << rp << std::endl; p[0] = c0[0]*0.5+c1[0]*0.5; p[1] = c0[1]*0.5+c1[1]*0.5; p[2] = c0[2]*0.5+c1[2]*0.5;