From aa0d3ee4e05c4232c3940d2821ece7881f85e6b3 Mon Sep 17 00:00:00 2001 From: Chernyshenko Date: Thu, 7 Jun 2018 15:11:40 +0300 Subject: [PATCH] latest features in slice_func --- Examples/GridTools/slice_func.cpp | 1306 +---------------------------- 1 file changed, 29 insertions(+), 1277 deletions(-) diff --git a/Examples/GridTools/slice_func.cpp b/Examples/GridTools/slice_func.cpp index 44dbd80..142179f 100644 --- a/Examples/GridTools/slice_func.cpp +++ b/Examples/GridTools/slice_func.cpp @@ -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); - 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) || (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 " << cutnodes[q].LocalID() << std::endl; - } - else if( l1 > 5.0e-2*l ) - { - cutnodes[q] = m.CreateNode(p); - //std::cout << "created new node " << cutnodes[q].LocalID() << std::endl; - } - else - { - material[cnodes[q]] = 2; - cnodes[q].SetMarker(slice); - //std::cout << "use old node " << std::endl; - } - } - else if( s0 ) - { - if( !centernode.isValid() ) centernode = m.CreateNode(c0); - cutnodes[q] = centernode; - } - } - //else cutnodes[q] = cnodes[q]; - } - - 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 " << cnodes[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); - } - - for(int q = 0; q < (int)cnodes.size(); ++q) if( cnodes[q].GetMarker(slice) ) cutnodes[q] = cnodes[q]; - - //now find cutedges - - for(int q = 0; q < (int)cedges.size(); ++q ) - { - if( !cedges[q].GetMarker(slice ) ) - { - Node n1 = cutnodes[indx[cedges[q].getBeg()]]; - Node n2 = cutnodes[indx[cedges[q].getEnd()]]; - if( n1.isValid() && n2.isValid() && n1 != n2 ) - { - edge_nodes[0] = n1; - edge_nodes[1] = n2; - cutedges[q] = m.CreateEdge(edge_nodes).first; - } - } - } - - for(int q = 0; q < (int)cutedges.size(); ++q) if( cutedges[q].isValid() ) - { - //std::cout << "New cut edge " << cutedges[q].LocalID() << " at " << q << " on plane with " << cedges[q].LocalID() << " " << cutedges[q].getBeg().LocalID() << "<->" << cutedges[q].getEnd().LocalID() << " original " << cedges[q].getBeg().LocalID() << "<->" << cedges[q].getEnd().LocalID() << std::endl; - Edge n = cutedges[q]; - material[n] = 2; - n.Bulk(sliced) = 1; - n.SetMarker(slice); - } - - for(int q = 0; q < (int)cedges.size(); ++q) if( cedges[q].GetMarker(slice) ) cutedges[q] = cedges[q]; - - ElementArray split_edges(&m); - //run over pyramids and collect faces, they should be already ordered into loops - //although we still have to check they form closed loop - std::map vstcnt; - for(int q = 0; q < (int)cfaces.size(); ++q) - { - split_edges.clear(); - ElementArray fedges = cfaces[q].getEdges(); - for(int r = 0; r < (int)fedges.size(); ++r) - { - if( cutedges[indx[fedges[r]]].isValid() ) - split_edges.push_back(cutedges[indx[fedges[r]]]); - } - - if( split_edges.size() > 2 ) //at least a triangle - { - visit_count.clear(); - for(int r = 0; r < (int)split_edges.size(); ++r) - { - visit_count[split_edges[r].getBeg()]++; - visit_count[split_edges[r].getEnd()]++; - } - bool simple = true; - for(std::map::iterator jt = visit_count.begin(); jt != visit_count.end(); ++jt) - simple &= (jt->second == 2); - - if( simple ) - { - for(int r = 0; r < split_edges.size(); ++r) vstcnt[split_edges[r]]++; - Face f = m.CreateFace(split_edges).first; - //std::cout << "Created face " << f.LocalID() << " with " << split_edges.size() << " edges: "; - //for(int r = 0; r < split_edges.size(); ++r) std::cout << "EDGE:" << split_edges[r].LocalID() << " "; - //std::cout << " original FACE:" << cfaces[q].LocalID() << std::endl; - material[f] = 2; - f.Bulk(sliced) = 1; - split_faces.push_back(f); - } - else - { - //std::cout << "Face not created, node visits:"; - //for(std::map::iterator jt = visit_count.begin(); jt != visit_count.end(); ++jt) - // std::cout << " NODE:" << jt->first.LocalID() << " " << jt->second; - //std::cout << " edges: "; - //for(int r = 0; r < split_edges.size(); ++r) std::cout << "EDGE:" << split_edges[r].LocalID() << " "; - //std::cout << " original FACE:" << cfaces[q].LocalID() << std::endl; - } - } - else if( !split_edges.empty() ) - { - //std::cout << "Face not created, edges:"; - //for(int r = 0; r < (int)split_edges.size(); ++r) - // std::cout << " EDGE:" << split_edges[r].LocalID(); - //std::cout << " original FACE:" << cfaces[q].LocalID() << std::endl; - } - } - - //resulting surface may touch the element with one or more edge - //making it impossible to separate in conformal way - // - // in this case we have to isolate this edge with additional faces - // - // check if the edge is counted twice and appears on original element - visit_count.clear(); - cedges.SetMarker(cmrk); - - bool isolate_alogirthm = false; - for(std::map::iterator jt = vstcnt.begin(); jt != vstcnt.end(); ++jt) - if( jt->first.GetMarker(cmrk) && jt->second > 1 ) - { - isolate_alogirthm = true; - jt->first.SetMarker(isolate); - } - - if( isolate_alogirthm ) - //if( false ) - { - //std::cout << "Isolation algorithm" << std::endl; - //on pyramids, build faces that start from cutedges and end up with edge on original edge - for(int q = 0; q < (int)cfaces.size(); ++q) - { - split_edges.clear(); - - ElementArray fedges = cfaces[q].getEdges(); - bool have_isolate = false; - for(int r = 0; r < (int)fedges.size(); ++r) have_isolate |= fedges[r].GetMarker(isolate); - if( have_isolate ) - { - for(int r = 0; r < (int)fedges.size(); ++r) if( !fedges[r].GetMarker(isolate) && !fedges[r].GetMarker(slice) && cutedges[indx[fedges[r]]].isValid() ) - { - split_edges.clear(); - int cnt = 0; - Edge e; - Face f; - int mat = material[fedges[r]]; - e = fedges[r]; - Node e1n1 = e.getBeg(); - Node e1n2 = e.getEnd(); - Node e2n1 = cutnodes[indx[e.getBeg()]]; - Node e2n2 = cutnodes[indx[e.getEnd()]]; - - //std::cout << "Cell edge: " << e.LocalID() << " " << e.getBeg().LocalID() << "<->" << e.getEnd().LocalID() << std::endl; - //std::cout << "Cut edge: " << e.LocalID() << " " << cutnodes[indx[e.getBeg()]].LocalID() << "<->" << cutnodes[indx[e.getEnd()]].LocalID() << std::endl; - split_edges.push_back(e); - - if( e2n1.isValid() && e2n1 != e1n1 ) - { - edge_nodes[0] = e1n1; - edge_nodes[1] = e2n1; - e = m.CreateEdge(edge_nodes).first; - material[e] = mat; - e.Bulk(sliced) = 1; - //std::cout << "New edge: " << e.LocalID() << " " << e.getBeg().LocalID() << "<->" << e.getEnd().LocalID() << std::endl; - split_edges.push_back(e); - cnt++; - } - - e = cutedges[indx[fedges[r]]]; - //std::cout << "Cut edge: " << e.LocalID() << " " << e.getBeg().LocalID() << "<->" << e.getEnd().LocalID() << std::endl; - split_edges.push_back(e); - - if( e2n2.isValid() && e2n2 != e1n2 ) - { - edge_nodes[0] = e1n2; - edge_nodes[1] = e2n2; - e = m.CreateEdge(edge_nodes).first; - material[e] = mat; - e.Bulk(sliced) = 1; - //std::cout << "New edge: " << e.LocalID() << " " << e.getBeg().LocalID() << "<->" << e.getEnd().LocalID() << std::endl; - split_edges.push_back(e); - cnt++; - } - - if( cnt ) - { - f = m.CreateFace(split_edges).first; - //std::cout << "New face " << f.LocalID() << std::endl; - material[f] = mat; - f.Bulk(sliced) = 1; - split_faces.push_back(f); - } - } - } - } - cedges.RemMarker(isolate); - } - - cedges.RemMarker(cmrk); - } - - if( !split_faces.empty() ) - { - int lid = it->LocalID(); - ElementArray ret = Cell::SplitCell(it->self(),split_faces,0); - - - if(!simple ) - { - //std::cout << (simple?"simple":"complex") << " algorithm, split cell " << lid << " with " << split_faces.size() << " faces "; - //std::cout << " result in " << ret.size() << " cells:"; - //for(int q = 0; q < (int)ret.size(); ++q) std::cout << " " << ret[q].LocalID(); - //std::cout << std::endl; - } - - ret.SetMarker(mrk); - nslice++; - - if( false )//ret.size() == 1 && !simple ) - { - ElementArray cedges = ret[0].getEdges(); - - - cedges.SetMarker(cmrk); - - std::cout << "Original cut edges: " << std::endl; - for(int k = 0; k < edges.size(); ++k) std::cout << edges[k].LocalID() << " "; - std::cout << std::endl; - - std::cout << "Original cell edges: " << std::endl; - for(int k = 0; k < cedges.size(); ++k) std::cout << cedges[k].LocalID() << " "; - std::cout << std::endl; - - std::cout << "Cut faces edges: " << std::endl; - std::map vstcnt; - for(int q = 0; q < split_faces.size(); ++q) - { - ElementArray sedges = split_faces[q].getEdges(); - std::cout << "Face " << split_faces[q].LocalID() << ":"; - for(int r = 0; r < sedges.size(); ++r) - { - vstcnt[sedges[r]]++; - std::cout << " " << sedges[r].LocalID() << "(" << (sedges[r].GetMarker(cmrk)?"s":"i") << ")"; - } - std::cout << std::endl; - } - std::cout << "Visit counts:"; - for(std::map::iterator jt = vstcnt.begin(); jt != vstcnt.end(); ++jt) - std::cout << " EDGE:" << jt->first.LocalID() << "(" << (jt->first.GetMarker(cmrk)?"s":"i") << ") " << jt->second; - std::cout < cfaces = ret[0].getFaces(); - for(int k = 0; k < cfaces.size(); ++k) - { - ElementArray sedges = cfaces[k].getEdges(); - std::cout << "Face " << cfaces[k].LocalID() << ":"; - for(int r = 0; r < sedges.size(); ++r) - { - vstcnt[sedges[r]]++; - std::cout << " " << sedges[r].LocalID() << "(" << (sedges[r].GetMarker(cmrk)?"s":"i") << ")"; - } - std::cout << std::endl; - } - std::cout << "Visit counts including cell's faces:"; - for(std::map::iterator jt = vstcnt.begin(); jt != vstcnt.end(); ++jt) - std::cout << " EDGE:" << jt->first.LocalID() << "(" << (jt->first.GetMarker(cmrk)?"s":"i") << ") " << jt->second; - std::cout <(" << c2[0] << "," << c2[1] << "," << c2[2] << ")" << std::endl; - } - std::cout << "Faces:" << std::endl; - for(int q = 0; q < split_faces.size(); ++q) - { - ElementArray sedges = split_faces[q].getEdges(); - for(int r = 0; r < sedges.size(); ++r) - { - Storage::real_array c1 = sedges[r].getBeg().Coords(); - Storage::real_array c2 = sedges[r].getEnd().Coords(); - std::cout << "(" << c1[0] << "," << c1[1] << "," << c1[2] << ")<->(" << c2[0] << "," << c2[1] << "," << c2[2] << ")" << std::endl; - } - } - - cedges.RemMarker(cmrk); - - } - } - } - } - m.ReleaseMarker(isolate); - m.ReleaseMarker(cmrk); - m.ReleaseMarker(visited); - m.DeleteTag(indx); - - std::cout << "sliced cells: " << nslice << std::endl; - - if( !Element::CheckConnectivity(&m) ) - std::cout << "Connectivity is broken" << std::endl; - - { - int mat[3] = {0,0,0}; - int tot = 0; - for(Mesh::iteratorFace it = m.BeginFace(); it != m.EndFace(); ++it) - { - mat[material[*it]]++; - tot++; - } - std::cout << "face materials, 0: " << mat[0] << " 1: " << mat[1] << " 2: " << mat[2] << std::endl; - } - - for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it) - { - int mat[3] = {0,0,0}; - ElementArray faces = it->getFaces(); - for(int k = 0; k < (int)faces.size(); ++k) - mat[material[faces[k]]]++; - if( !(mat[0] == 0 || mat[1] == 0) ) - { - material[*it] = 2; - std::cout << "oops, materials for cell " << it->LocalID() << " faces were not split, 0: " << mat[0] << " ,1: " << mat[1] << " ,2: " << mat[2] << " slice edges " << it->getEdges(slice).size() << std::endl; - } - else if( mat[0] != 0 ) material[*it] = 0; - else if( mat[1] != 0 ) material[*it] = 1; - else - { - //double cnt[3]; - //it->Centroid(cnt); - //double v = func(cnt[0],cnt[1],cnt[2],type); - //it->Integer(material) = (v <= 0.0 ? 0 : 1); - material[*it] = 2; - std::cout << "oops cannot determine material for cell, all faces have type 2, set to " << material[*it]; - } - } - - - { - int mat[3] = {0,0,0}; - int tot = 0; - for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it) - { - mat[material[*it]]++; - tot++; - } - std::cout << "cell materials, 0: " << mat[0] << " 1: " << mat[1] << " 2: " << mat[2] << std::endl; - } - - Tag collapse = m.CreateTag("COLLAPSE",DATA_INTEGER,CELL,NONE,1); - Tag el_mat = m.CreateTag("ELLIPSE_MATRIX",DATA_REAL,CELL,CELL,9); - Tag el_cnt = m.CreateTag("ELLIPSE_CENTER",DATA_REAL,CELL,CELL,3); - Tag level = m.CreateTag("LEVEL",DATA_REAL,CELL|FACE|EDGE|NODE,NONE,1); - - int ncollapsed = 0; - for(Mesh::iteratorElement it = m.BeginElement(NODE|EDGE|FACE|CELL); it != m.EndElement(); ++it) - { - double cnt[3]; - it->Centroid(cnt); - it->Real(level) = func(cnt[0],cnt[1],cnt[2],type); - } - - for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it) - { - double vv, v, nvv, cnt[3]; - v = it->Volume(); - vv = v; - nvv = 1; - ElementArray cells = it->NeighbouringCells(); - for(ElementArray::iterator jt = cells.begin(); jt != cells.end(); ++jt) - { - vv += jt->Volume(); - nvv++; - } - vv /= nvv; - if( v < vv*0.01 ) - { - if( false ) - { - std::cout << "looks like cell " << it->LocalID() << " is too small, volume " << it->Volume() << std::endl; - for(ElementArray::iterator jt = cells.begin(); jt != cells.end(); ++jt) - { - std::cout << "adjacent cell " << jt->LocalID() << " volume " << jt->Volume() << std::endl; - } - } - it->Integer(collapse) = 1; - - /* - it->Centroid(cnt); - vMatrix A(3,3), b(3,1); - A(0,0) = unknown(1,0); - A(1,1) = unknown(1,1); - A(2,2) = unknown(1,2); - A(1,2) = A(2,1) = unknown(0,3); - A(1,3) = A(3,1) = unknown(0,4); - A(2,3) = A(3,2) = unknown(0,5); - b(0,0) = unknown(cnt[0],6); - b(1,0) = unknown(cnt[1],7); - b(2,0) = unknown(cnt[2],8); - - ElementArray nodes = it->getNodes(); - */ - - - //find biggest face - /* - ElementArray faces = it->getFaces(); - Cell n = InvalidCell(); - double area = 0; - for(int k = 0; k < (int)faces.size(); ++k) - { - Cell nn = it->Neighbour(faces[k]); - if( nn.isValid() ) - { - if( faces[k].Area() > area ) - { - n = nn; - area = faces[k].Area(); - } - } - } - //merge cells - if( n.isValid() ) - { - ElementArray unite(&m,2); - unite[0] = it->self(); - unite[1] = n; - Cell::UniteCells(unite,0); - } - */ - } - - } - - for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it) - if( material[*it] == 0 ) - it->Delete(); - - for(ElementType etype = FACE; etype >= NODE; etype = PrevElementType(etype) ) - for(Mesh::iteratorElement it = m.BeginElement(etype); it != m.EndElement(); ++it) - if( it->nbAdjElements(CELL) == 0 ) it->Delete(); - - - m.ReleaseMarker(slice,NODE|EDGE); - m.ReleaseMarker(original,NODE|EDGE); - m.ReleaseMarker(mrk,EDGE|FACE|NODE); - m.Save(grid_out); - return 0; -} -int main_old(int argc, char ** argv) -{ - if( argc < 2 ) - { - std::cout << "Usage: " << argv[0] << " mesh [mesh_out=grid.pmf] [type=perforated_strip]" << std::endl; - return -1; - } - - std::string grid_out = "grid.pmf"; - - int type = 0; - if( argc > 2 ) grid_out = std::string(argv[2]); - if( argc > 3 ) type = atoi(argv[3]); - - - - Mesh m; - 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); - 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(); @@ -1376,8 +212,8 @@ int main_old(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); - material[it->getBeg()] = (r0 <= 0)? 0 : 1; - material[it->getEnd()] = (r1 <= 0)? 0 : 1; + if (material[it->getBeg()] == 3) material[it->getBeg()] = (r0 <= 0)? 0 : 1; + if (material[it->getEnd()] == 3) 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))) ) { @@ -1496,6 +332,7 @@ int main_old(int argc, char ** argv) { 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; @@ -1533,24 +370,14 @@ int main_old(int argc, char ** argv) //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; + + bool s0 = false; if( fabs(r0) < 1.0e-6 ) s0 = true; //std::cout << "Centernode " << (s0?"sliced":"") << " r=" << r0 << std::endl; @@ -1805,7 +632,7 @@ int main_old(int argc, char ** argv) if( simple ) { //gather loop by loop and create faces - int loop_cnt = 0; + //int loop_cnt = 0; //order edges ElementArray order_edges(&m); @@ -1922,7 +749,7 @@ int main_old(int argc, char ** argv) 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; @@ -1974,22 +801,26 @@ int main_old(int argc, char ** argv) l0 = sqrt(l0); l1 = sqrt(l1); l = l0+l1; - if( l0 < 5.0e-2*l ) //edge goes through centernode + if( l0 < 1.0e-3*l ) //edge goes through centernode { if( !centernode.isValid() ) centernode = m.CreateNode(c0); cutnodes[q] = centernode; //std::cout << "selected centernode " << cutnodes[q].LocalID() << std::endl; } - else if( l1 > 5.0e-2*l ) + else if( l1 > 1.0e-3*l ) { cutnodes[q] = m.CreateNode(p); //std::cout << "created new node " << cutnodes[q].LocalID() << std::endl; } - else + else if( !cnodes[q].GetMarker(slice)) { material[cnodes[q]] = 2; cnodes[q].SetMarker(slice); + ElementArray nedges = cnodes[q]->getEdges(); + for(int r = 0; r < nedges.size();++r) + if( material[nedges[r]->getBeg()] == 2 && material[nedges[r]->getEnd()] == 2) + material[nedges[r]] = 2; //std::cout << "use old node " << std::endl; } } @@ -2192,7 +1023,7 @@ int main_old(int argc, char ** argv) if( !split_faces.empty() ) { - int lid = it->LocalID(); + //int lid = it->LocalID(); ElementArray ret = Cell::SplitCell(it->self(),split_faces,0); @@ -2342,90 +1173,7 @@ int main_old(int argc, char ** argv) std::cout << "cell materials, 0: " << mat[0] << " 1: " << mat[1] << " 2: " << mat[2] << std::endl; } - Tag collapse = m.CreateTag("COLLAPSE",DATA_INTEGER,CELL,NONE,1); - Tag el_mat = m.CreateTag("ELLIPSE_MATRIX",DATA_REAL,CELL,CELL,9); - Tag el_cnt = m.CreateTag("ELLIPSE_CENTER",DATA_REAL,CELL,CELL,3); - Tag level = m.CreateTag("LEVEL",DATA_REAL,CELL|FACE|EDGE|NODE,NONE,1); - int ncollapsed = 0; - for(Mesh::iteratorElement it = m.BeginElement(NODE|EDGE|FACE|CELL); it != m.EndElement(); ++it) - { - double cnt[3]; - it->Centroid(cnt); - it->Real(level) = func(cnt[0],cnt[1],cnt[2],type); - } - - for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it) - { - double vv, v, nvv, cnt[3]; - v = it->Volume(); - vv = v; - nvv = 1; - ElementArray cells = it->NeighbouringCells(); - for(ElementArray::iterator jt = cells.begin(); jt != cells.end(); ++jt) - { - vv += jt->Volume(); - nvv++; - } - vv /= nvv; - if( v < vv*0.01 ) - { - if( false ) - { - std::cout << "looks like cell " << it->LocalID() << " is too small, volume " << it->Volume() << std::endl; - for(ElementArray::iterator jt = cells.begin(); jt != cells.end(); ++jt) - { - std::cout << "adjacent cell " << jt->LocalID() << " volume " << jt->Volume() << std::endl; - } - } - it->Integer(collapse) = 1; - - /* - it->Centroid(cnt); - vMatrix A(3,3), b(3,1); - A(0,0) = unknown(1,0); - A(1,1) = unknown(1,1); - A(2,2) = unknown(1,2); - A(1,2) = A(2,1) = unknown(0,3); - A(1,3) = A(3,1) = unknown(0,4); - A(2,3) = A(3,2) = unknown(0,5); - b(0,0) = unknown(cnt[0],6); - b(1,0) = unknown(cnt[1],7); - b(2,0) = unknown(cnt[2],8); - - ElementArray nodes = it->getNodes(); - */ - - - //find biggest face - /* - ElementArray faces = it->getFaces(); - Cell n = InvalidCell(); - double area = 0; - for(int k = 0; k < (int)faces.size(); ++k) - { - Cell nn = it->Neighbour(faces[k]); - if( nn.isValid() ) - { - if( faces[k].Area() > area ) - { - n = nn; - area = faces[k].Area(); - } - } - } - //merge cells - if( n.isValid() ) - { - ElementArray unite(&m,2); - unite[0] = it->self(); - unite[1] = n; - Cell::UniteCells(unite,0); - } - */ - } - - } /* for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it) if( material[*it] == 0 || it->Integer(collapse) ) @@ -2439,6 +1187,10 @@ int main_old(int argc, char ** argv) m.ReleaseMarker(slice,NODE|EDGE); m.ReleaseMarker(original,NODE|EDGE); m.ReleaseMarker(mrk,EDGE|FACE|NODE); + + m.DeleteTag(material); + m.DeleteTag(level); + m.Save(grid_out); return 0; } -- 2.26.2