Commit aa0d3ee4 by Chernyshenko

### latest features in slice_func

parent 24f3848d
 ... ... @@ -4,25 +4,7 @@ using namespace INMOST; double func(double x, double y, double z, int n) { double Lx = 0., Rx = 20, Ly = -2.5, Ry = 17.5; if (x > Rx){ if (y > Ry) return -sqrt( (x-Rx)*(x-Rx) + (y-Ry)*(y-Ry) ); if (y < Ly) return -sqrt( (x-Rx)*(x-Rx) + (y-Ly)*(y-Ly) ); return -(x-Rx); } if (x < Lx){ if (y < Ly) return -sqrt( (x-Lx)*(x-Lx) + (y-Ly)*(y-Ly) ); if (y > Ry) return -sqrt( (x-Lx)*(x-Lx) + (y-Ry)*(y-Ry) ); 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) ); } double func2(double x, double y, double z, int n) { if( n == 0 ) return sqrt(x*x+y*y)-0.25; else if( n == 1 ) ... ... @@ -57,7 +39,7 @@ double func2(double x, double y, double z, int n) return fmin( fmin(x-Lx,Rx-x), fmin(y-Ly, Ry-y) ); */ double Lx = 0., Rx = 2, Ly = 0.0, Ry = 4.5; double Lx = 0., Rx = 20, Ly = -2.5, Ry = 17.5; if (x > Rx) { if (y > Ry) return -sqrt( (x-Rx)*(x-Rx) + (y-Ry)*(y-Ry) ); ... ... @@ -129,55 +111,6 @@ double search(double r0, double r1, double c0[3], double c1[3], double p[3],int return rp; } double search_zero(double r0, double r1, double c0[3], double c1[3], double p[3],int type) { int iters = 0; double rp = 1.0e+20; do { p[0] = 0.5*c1[0] + 0.5*c0[0]; p[1] = 0.5*c1[1] + 0.5*c0[1]; p[2] = 0.5*c1[2] + 0.5*c0[2]; rp = func(p[0],p[1],p[2],type); if( fabs(rp) < 1.0e-8 ) { if( fabs(r0) < 1.0e-8 ) { c0[0] = p[0]; c0[1] = p[1]; c0[2] = p[2]; r0 = rp; } else { c1[0] = p[0]; c1[1] = p[1]; c1[2] = p[2]; r1 = rp; } } else { if( fabs(r0) > 1.0e-8 ) { c0[0] = p[0]; c0[1] = p[1]; c0[2] = p[2]; r0 = rp; } else { c1[0] = p[0]; c1[1] = p[1]; c1[2] = p[2]; r1 = rp; } } iters++; } while( iters < 20 ); return rp; } double search_zero(double r0, double r1, double c0[3], double c1[3], double p[3],int type) ... ... @@ -252,1118 +185,21 @@ int main(int argc, char ** argv) //m.RemTopologyCheck(THROW_EXCEPTION); 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); TagReal level = m.CreateTag("level",DATA_REAL,NODE,NONE,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(); MarkerType mrk = m.CreateMarker(); int nslice = 0, nmark = 0; for(Mesh::iteratorEdge it = m.BeginEdge(); it != m.EndEdge(); ++it) if( !it->GetMarker(mrk) ) { double p[3],pc0[3],pc1[3]; Storage::real_array c0 = it->getBeg()->Coords(); 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); 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) || (fabs(r0*r1) < 1.0e-12 && ((fabs(r0) < 1.0e-6) ^ (fabs(r1) < 1.0e-6))) ) { pc0[0] = c0[0], pc0[1] = c0[1], pc0[2] = c0[2]; pc1[0] = c1[0], pc1[1] = c1[1], pc1[2] = c1[2]; if((fabs(r0) < 1.0e-6) ^ (fabs(r1) < 1.0e-6)) { search_zero(r0,r1,pc0,pc1,p,type); if( fabs(r0) < 1.0e-6 ) { material[it->getBeg()] = 2; it->getBeg()->SetMarker(slice); nmark++; } if( fabs(r1) < 1.0e-6 ) { material[it->getEnd()] = 2; it->getEnd()->SetMarker(slice); nmark++; } } else { 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; //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 < 5.0e-2*l ) { material[it->getBeg()] = 2; it->getBeg()->SetMarker(slice); nmark++; } else if( l1 < 5.0e-2*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 ) { material[it->getBeg()] = 2; it->getBeg()->SetMarker(slice); nmark++; } if( fabs(r1) < 1.0e-6 ) { material[it->getEnd()] = 2; it->getEnd()->SetMarker(slice); nmark++; } } } std::cout << "sliced edges: " << nslice << " marked nodes: " << nmark << std::endl; if( !Element::CheckConnectivity(&m) ) std::cout << "Connectivity is broken" << std::endl; { int mat[3] = {0,0,0}; int tot = 0; for(Mesh::iteratorNode it = m.BeginNode(); it != m.EndNode(); ++it) { mat[material[*it]]++; tot++; } std::cout << "node materials, 0: " << mat[0] << " 1: " << mat[1] << " 2: " << mat[2] << std::endl; } for(Mesh::iteratorEdge it = m.BeginEdge(); it != m.EndEdge(); ++it) { int mat[3] = {0,0,0}; mat[material[it->getBeg()]]++; mat[material[it->getEnd()]]++; if( !(mat[0] == 0 || mat[1] == 0) ) { 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 ) material[*it] = 0; else if( mat[1] != 0 ) material[*it] = 1; else material[*it] = 2; } nslice = 0; MarkerType unique = m.CreateMarker(); 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++; //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 " << 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; 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) || (fabs(r0*r1) < 1.0e-12 && ((fabs(r0) < 1.0e-6) ^ (fabs(r1) < 1.0e-6)))) { pc0[0] = c0[0], pc0[1] = c0[1], pc0[2] = c0[2]; pc1[0] = c1[0], pc1[1] = c1[1], pc1[2] = c1[2]; if((fabs(r0) < 1.0e-6) ^ (fabs(r1) < 1.0e-6)) { search_zero(r0,r1,pc0,pc1,p,type); } else { 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 < 5.0e-2*l ) //edge goes through centernode { if( !centernode.isValid() ) centernode = m.CreateNode(c0); cutnodes[q] = centernode; //std::cout << "selected centernode " << std::endl; } else if( l1 > 5.0e-2*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]; std::pair en = m.CreateEdge(edge_nodes); //if( en.second ) { Edge e = en.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]; std::pair en = m.CreateEdge(edge_nodes); //if( en.second ) { Edge e = en.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]; std::pair en = m.CreateEdge(edge_nodes); //if( en.second ) { Edge e = en.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() ) { int k = 0; for(int q = 0; q < split_edges.size(); ++q) { if( !split_edges[q].GetMarker(unique) ) { split_edges[q].SetMarker(unique); split_edges[k++] = split_edges[q]; } } split_edges.RemMarker(unique); split_edges.resize(k); //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(unique); 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( material[*it] != 2 ) std::cout << "Edge supposed to get material 2, but have " << material[*it] << std::endl; it->SetMarker(slice); nmark++; } std::cout << "sliced faces: " << nslice << " marked edges: " << nmark << std::endl; if( !Element::CheckConnectivity(&m) ) std::cout << "Connectivity is broken" << std::endl; { int mat[3] = {0,0,0}; int tot = 0; for(Mesh::iteratorEdge it = m.BeginEdge(); it != m.EndEdge(); ++it) #pragma omp parallel for for(int k = 0; k < m.NodeLastLocalID(); ++k) if(m.isValidNode(k)) { mat[material[*it]]++; tot++; Node it = m.NodeByLocalID(k); it->SetMarker(original); material[it] = 3; Storage::real_array c0 = it->Coords(); level[it] = func(c0[0],c0[1],c0[2],type); } std::cout << "edge materials, 0: " << mat[0] << " 1: " << mat[1] << " 2: " << mat[2] << std::endl; } for(Mesh::iteratorFace it = m.BeginFace(); it != m.EndFace(); ++it) { int mat[3] = {0,0,0}; ElementArray edges = it->getEdges(); for(int k = 0; k < (int)edges.size(); ++k) mat[material[edges[k]]]++; if( !(mat[0] == 0 || mat[1] == 0) ) { 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 ) 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(); 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.empty() ) //these should form a triangle { //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);