Commit 78310df9 authored by Kirill Terekhov's avatar Kirill Terekhov
Browse files

Update

parent 18394864
Pipeline #146 passed with stages
in 10 minutes and 37 seconds
......@@ -3,7 +3,8 @@ project(GridGen)
add_executable(FixFaults fix_faults.cpp)
add_executable(FixTiny fix_tiny.cpp)
add_executable(UniteFaces unite_faces.cpp)
add_executable(Dual dual.cpp)
add_executable(Dual dual.cpp)
add_executable(Tetra tetra.cpp)
add_executable(Slice slice.cpp)
add_executable(SliceFunc slice_func.cpp)
add_executable(SplitNonplanar split_nonplanar.cpp)
......@@ -66,7 +67,17 @@ if(USE_MPI)
set_target_properties(Dual PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
endif()
endif(USE_MPI)
install(TARGETS Dual EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
install(TARGETS Dual EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
target_link_libraries(Tetra inmost)
if(USE_MPI)
message("linking Tetra with MPI")
target_link_libraries(Tetra ${MPI_LIBRARIES})
if(MPI_LINK_FLAGS)
set_target_properties(Tetra PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
endif()
endif(USE_MPI)
install(TARGETS Tetra EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
target_link_libraries(Slice inmost)
......
#include "inmost.h"
#include <stdio.h>
#include <deque>
//#include "tetgen/tetgen.h"
using namespace INMOST;
typedef Storage::real real;
typedef Storage::real_array real_array;
const bool divide_cells = true;
real dot(real a[3], real b[3])
{
return a[0]*b[0] + a[1]*b[1] + a[2]*b[2];
}
void cross(real a[3], real b[3], real out[3])
{
out[0] = a[1]*b[2] - a[2]*b[1];
out[1] = a[2]*b[0] - a[0]*b[2];
out[2] = a[0]*b[1] - a[1]*b[0];
}
struct node
{
real coo[2];
node(const real _coo[2])
{
coo[0] = _coo[0];
coo[1] = _coo[1];
}
node(const node & b)
{
coo[0] = b.coo[0];
coo[1] = b.coo[1];
}
node & operator = (node const & b)
{
coo[0] = b.coo[0];
coo[1] = b.coo[1];
return *this;
}
};
struct segment
{
real beg[2];
real end[2];
segment(const real _beg[2], const real _end[2])
{
for(int k = 0; k < 2; ++k)
{
beg[k] = _beg[k];
end[k] = _end[k];
}
}
segment(const segment & b)
{
for(int k = 0; k < 2; ++k)
{
beg[k] = b.beg[k];
end[k] = b.end[k];
}
}
segment & operator = (segment const & b)
{
for(int k = 0; k < 2; ++k)
{
beg[k] = b.beg[k];
end[k] = b.end[k];
}
return *this;
}
};
//check whether two segments intersect
bool intersect_segments(const segment & a, const segment & b, bool print = false)
{
const real eps = 1.0e-9;
real div1 = (a.beg[0] - a.end[0])*(b.beg[1] - b.end[1]);
real div2 = (a.beg[1] - a.end[1])*(b.beg[0] - b.end[0]);
real div = div1 - div2;
real t, find[2];
if( print )
{
std::cout << "(" << a.beg[0] << "," << a.beg[1] << ",0) <-> (" << a.end[0] << "," << a.end[1] << ",0)" << std::endl;
std::cout << "(" << b.beg[0] << "," << b.beg[1] << ",0) <-> (" << b.end[0] << "," << b.end[1] << ",0)" << std::endl;
}
if( print ) std::cout << "div = " << div << " div1 " << div1 << " div2 " << div2 << " against " << eps*(std::abs(div1) + std::abs(div2)) << std::endl;
if (std::abs(div) <= eps*(std::abs(div1) + std::abs(div2))) //segments are parallel
{
//return true;
if( print ) std::cout << "parallel segments, div = " << div << std::endl;
if( std::abs(a.beg[0]-a.end[0]) < eps*(std::abs(a.beg[0])+std::abs(a.end[0])) ) //parallel to x
{
real a0 = a.beg[1], a1 = a.end[1];
real b0 = b.beg[1], b1 = b.end[1];
if( print ) std::cout << "parallel to x, projection on Oy: a [" << a0 << "," << a1 << "], b [" << b0 << "," << b1 << "]" << std::endl;
if( a0 > a1 ) std::swap(a0,a1);
if( b0 > b1 ) std::swap(b0,b1);
if( a1 > b1 ) //change a and b
{
std::swap(a0,b0);
std::swap(a1,b1);
}
if( a1 > b0+eps )
{
if( print ) std::cout << "intersection" << std::endl;
return true;
}
if( print ) std::cout << "no intersection" << std::endl;
return false;
}
//reconstruct k in y = kx+b
real ak = (a.end[1]-a.beg[1])/(a.end[0]-a.beg[0]);
real bk = (b.end[1]-b.beg[1])/(b.end[0]-b.beg[0]);
//reconstruct bs in y = kx+b
real ab = (a.end[0]*a.beg[1] - a.beg[0]*a.end[1])/(a.end[0]-a.beg[0]);
real bb = (b.end[0]*b.beg[1] - b.beg[0]*b.end[1])/(b.end[0]-b.beg[0]);
if( print )
{
std::cout << "parallel lines:" << std::endl;
std::cout << "a: y = " << ak << "*x + " << ab << std::endl;
std::cout << "b: y = " << bk << "*x + " << bb << std::endl;
std::cout << "kdiff: " << ak - bk << std::endl;
std::cout << "bdiff: " << ab - bb << " against " << eps*(std::abs(ab)+std::abs(bb)) << std::endl;
}
if( std::abs(ab-bb) <= eps*(std::abs(ab)+std::abs(bb)) )
{
//test segments
int ind = ak > 1.0 ? 1 : 0; //choose projection of segments
real a0 = a.beg[ind], a1 = a.end[ind];
real b0 = b.beg[ind], b1 = b.end[ind];
if( print ) std::cout << "projection on O" << (ind?"y":"x") << ": a [" << a0 << "," << a1 << "], b [" << b0 << "," << b1 << "]" << std::endl;
if( a0 > a1 ) std::swap(a0,a1);
if( b0 > b1 ) std::swap(b0,b1);
if( a1 > b1 ) //change a and b
{
std::swap(a0,b0);
std::swap(a1,b1);
}
if( a1 > b0+eps )
{
if( print ) std::cout << "intersection" << std::endl;
return true;
}
if( print ) std::cout << "no intersection" << std::endl;
return false;
}
//printf("divisor is zero\n");
return false;
}
find[0] = ((a.beg[0]*a.end[1] - a.beg[1]*a.end[0])*(b.beg[0] - b.end[0]) - (a.beg[0] - a.end[0])*(b.beg[0]*b.end[1] - b.beg[1]*b.end[0])) / div;
find[1] = ((a.beg[0]*a.end[1] - a.beg[1]*a.end[0])*(b.beg[1] - b.end[1]) - (a.beg[1] - a.end[1])*(b.beg[0]*b.end[1] - b.beg[1]*b.end[0])) / div;
if( print ) std::cout << "intersection: (" << find[0] << "," << find[1] << ")" << std::endl;
//probably some of these tests are redundant
bool inter = true;
if (std::abs(a.end[0] - a.beg[0]) > eps)
{
t = (find[0] - a.beg[0]) / (a.end[0] - a.beg[0]);
if( print ) std::cout << "coef a on Ox " << t << std::endl;
if (t < eps || t > 1.0 - eps) inter = false;
}
if (std::abs(a.end[1] - a.beg[1]) > eps)
{
t = (find[1] - a.beg[1]) / (a.end[1] - a.beg[1]);
if( print ) std::cout << "coef a on Oy " << t << std::endl;
if (t < eps || t > 1.0 - eps) inter = false;
}
if (std::abs(b.end[0] - b.beg[0]) > eps)
{
t = (find[0] - b.beg[0]) / (b.end[0] - b.beg[0]);
if( print ) std::cout << "coef b on Ox " << t << std::endl;
if (t < eps || t > 1.0 - eps) inter = false;
}
if (std::abs(b.end[1] - b.beg[1]) > eps)
{
t = (find[1] - b.beg[1]) / (b.end[1] - b.beg[1]);
if( print ) std::cout << "coef b on Oy " << t << std::endl;
if (t < eps || t > 1.0 - eps) inter = false;
}
if( inter )
{
if( print ) std::cout << "intersection" << std::endl;
return true;
}
if( print ) std::cout << "no intersection" << std::endl;
return false;
}
//check whether any of the segments from setb intersect any segments from seta
bool check_intersect(const std::vector<segment> & segmentsa, const std::vector<segment> & segmentsb, bool print = false)
{
if( print )
{
std::cout << "segments [" << segmentsa.size() << "]:" << std::endl;
for(int k = 0; k < segmentsa.size(); ++k)
std::cout << "(" << segmentsa[k].beg[0] << "," << segmentsa[k].beg[1] << ",0) <-> (" << segmentsa[k].end[0] << "," << segmentsa[k].end[1] << ",0)" << std::endl;
std::cout << "joints [" << segmentsb.size() << "]:" << std::endl;
for(int k = 0; k < segmentsb.size(); ++k)
std::cout << "(" << segmentsb[k].beg[0] << "," << segmentsb[k].beg[1] << ",0) <-> (" << segmentsb[k].end[0] << "," << segmentsb[k].end[1] << ",0)" << std::endl;
}
for(int i = 0; i < (int)segmentsb.size(); ++i)
{
for(int j = 0; j < (int)segmentsa.size(); ++j)
{
if( print ) std::cout << "check segment " << j << " and joint " << i << std::endl;
if( intersect_segments(segmentsa[j],segmentsb[i],print) )
{
if( print ) std::cout << "there is intersection of segment " << j << " and joint " << i << std::endl;
return true;
}
}
}
if( print ) std::cout << "there is no intersection" << std::endl;
return false;
}
//check whether any of the centers of setb drop outside of loop represented by seta
bool check_dropout(const std::vector<node> & loop, const std::vector<segment> & segmentsb)
{
for(int i = 0; i < (int)segmentsb.size(); ++i)
{
real cnt[2];
cnt[0] = (segmentsb[i].beg[0]+segmentsb[i].end[0])*0.5;
cnt[1] = (segmentsb[i].beg[1]+segmentsb[i].end[1])*0.5;
int cn = 0; // the crossing number counter
// loop through all edges of the polygon
for (int j = 0; j < loop.size()-1; j++) // edge from V[i] to V[i+1]
{
real px = cnt[0];
real py = cnt[1];
real v0x = loop[j].coo[0];
real v0y = loop[j].coo[1];
real v1x = loop[j+1].coo[0];
real v1y = loop[j+1].coo[1];
if ( ((v0y <= py) && (v1y > py)) // an upward crossing
|| ((v0y > py) && (v1y <= py)))// a downward crossing
{
// compute the actual edge-ray intersect x-coordinate
real vt = (py - v0y) / (v1y - v0y);
if (px < v0x + vt * (v1x - v0x)) // P.x < intersect
++cn; // a valid crossing of y=P.y right of P.x
}
}
if( !(cn&1) ) //out of polygon
return true;
}
return false;
}
bool Planarity(Face f, double tol)
{
real nrm[3], cnt[3], area = f->Area();
f->Centroid(cnt);
f->UnitNormal(nrm);
ElementArray<Node> nodes = f->getNodes();
for(ElementArray<Node>::iterator it = nodes.begin(); it != nodes.end(); ++it)
{
real_array c = it->Coords();
if( fabs((c[0]-cnt[0])*nrm[0] + (c[1]-cnt[1])*nrm[1] + (c[2]-cnt[2])*nrm[2]) > tol*area )
return false;
}
return true;
}
int main(int argc, char ** argv)
{
if( argc < 2 )
{
printf("Usage: %s input_mesh [output_mesh]\n",argv[0]);
return -1;
}
Mesh m;
m.SetFileOption("ECL_CURVILINEAR","FALSE");
m.SetFileOption("ECL_SPLIT_GLUED","TRUE");
m.Load(argv[1]);
Mesh::GeomParam table;
table[ORIENTATION] = FACE;
m.PrepareGeometricData(table);
//reorder nodes according to RCM
Tag order = m.CreateTag("ORDER",DATA_INTEGER,NODE,NONE,1);
Tag idnum = m.CreateTag("RCM_ID",DATA_INTEGER,NODE,NONE,1);
for(Mesh::iteratorNode it = m.BeginNode(); it != m.EndNode(); ++it)
{
it->Integer(order) = (int)it->BridgeAdjacencies2Node(EDGE).size();
it->Integer(idnum) = -1;
}
int id = 0;
while( id < m.NumberOfNodes() )
{
Node select = InvalidNode();
Mesh::iteratorNode kt;
for(Mesh::iteratorNode it = m.BeginNode(); it != m.EndNode(); ++it)
{
if( it->Integer(idnum) == -1 )
{
select = it->self();
kt = it;
break;
}
}
for(Mesh::iteratorNode it = kt; it != m.EndNode(); ++it)
{
if( it->Integer(idnum) == -1 && it->Integer(order) < select->Integer(order) )
select = it->self();
}
std::deque<Node> q;
q.push_back(select);
select->Integer(idnum) = id++;
while(!q.empty())
{
select = q.front();
q.pop_front();
ElementArray<Node> adjacent = select->BridgeAdjacencies2Node(EDGE);
std::sort(adjacent.begin(),adjacent.end(),Mesh::IntegerComparator(&m,order));
for(ElementArray<Node>::iterator it = adjacent.begin(); it != adjacent.end(); ++it)
{
if( it->Integer(idnum) == -1 )
{
it->Integer(idnum) = id++;
q.push_back(it->self());
}
}
}
}
ElementArray<Node> nodes_rcm(&m,m.NumberOfNodes());
for(Mesh::iteratorNode it = m.BeginNode(); it != m.EndNode(); ++it)
nodes_rcm[m.NumberOfNodes() - it->Integer(idnum) - 1] = it->self();
m.DeleteTag(order);
m.DeleteTag(idnum);
//choose nodes to divide each cell
Tag divnode = m.CreateTag("DIVIDING_NODE",DATA_REFERENCE,CELL|FACE,NONE,1);
Tag proj_coords = m.CreateTag("PROJECTED_COORDS",DATA_REAL,NODE,NONE,2);
Tag mat = m.CreateTag("MATERIAL",DATA_INTEGER,CELL,NONE,1);
MarkerType sharenode = m.CreateMarker();
MarkerType centernode = m.CreateMarker();
std::vector<segment> segments, joints;
std::vector<node> loop;
tiny_map<HandleType,real,16> hits;
real nrm[3], cnt[3], ncnt[3], scnt[3], pcnt[3], fcnt[3], ray[3], orthx[3], orthy[3], d, nd;
//for(Mesh::iteratorNode it = m.BeginNode(); it != m.EndNode(); ++it)
for(ElementArray<Node>::iterator it = nodes_rcm.begin(); it != nodes_rcm.end(); ++it)
{
ElementArray<Face> node_faces = it->getFaces();
ElementArray<Edge> node_edges = it->getEdges();
ElementArray<Cell> cells = it->getCells();
it->SetMarker(sharenode);
node_faces.SetMarker(sharenode);
for(ElementArray<Cell>::iterator jt = cells.begin(); jt != cells.end(); ++jt) if( jt->GetGeometricType() != Element::Tet ) //skip tets
{
//retrive shared faces
ElementArray<Face> faces = jt->getFaces();
ElementArray<Edge> edges = jt->getEdges();
//check that cell or some face is not already divided
bool fail = false;
if( jt->Reference(divnode) != InvalidHandle() ) fail = true;
for(ElementArray<Face>::iterator kt = faces.begin(); kt != faces.end() && !fail; ++kt)
if( kt->GetMarker(sharenode) && kt->GetGeometricType() != Element::Tri && ((kt->Reference(divnode) != InvalidHandle() && kt->Reference(divnode) != *it) /*|| !Planarity(kt->self(),0.0001)*/) ) fail = true;
if( !fail )
{
//resolve a wired degenerate case
//mark edges of faces not shared by current node
node_edges.SetMarker(sharenode);
for(ElementArray<Face>::iterator kt = faces.begin(); kt != faces.end(); ++kt) if( !kt->GetMarker(sharenode) )
kt->getEdges().SetMarker(sharenode);
//count that all edges are marked
int count = 0;
for(ElementArray<Edge>::iterator kt = edges.begin(); kt != edges.end(); ++kt) if( kt->GetMarker(sharenode) ) count++;
edges.RemMarker(sharenode);
node_edges.RemMarker(sharenode);
//not all edges were marked - degenerate case
if( count != (int)edges.size() ) fail = true;
//check that all the edges are visible from the node
for(ElementArray<Face>::iterator kt = faces.begin(); kt != faces.end() && !fail; ++kt) if( kt->GetMarker(sharenode) && kt->GetGeometricType() != Element::Tri ) //skip tris
{
int id = kt->LocalID();
if( kt->Reference(divnode) == *it ) continue; //already chosen
kt->UnitNormal(nrm);
kt->Centroid(cnt);
d = dot(cnt,nrm);
ElementArray<Node> face_nodes = kt->getNodes();
//project coordinates
for(ElementArray<Node>::iterator mt = face_nodes.begin(); mt != face_nodes.end(); ++mt)
{
mt->Centroid(ncnt);
nd = dot(ncnt,nrm)-d;
//project onto plane
pcnt[0] = ncnt[0] - nrm[0]*nd;
pcnt[1] = ncnt[1] - nrm[1]*nd;
pcnt[2] = ncnt[2] - nrm[2]*nd;
//relative to center
pcnt[0] -= cnt[0];
pcnt[1] -= cnt[1];
pcnt[2] -= cnt[2];
//choose the basis
if( mt == face_nodes.begin() )
{
orthx[0] = pcnt[0];
orthx[1] = pcnt[1];
orthx[2] = pcnt[2];
cross(orthx,nrm,orthy);
}
//compute and record the coords
real_array c = mt->RealArray(proj_coords);
c[0] = dot(orthx,pcnt);
c[1] = dot(orthy,pcnt);
}
//gather all the segments in projected coordinates
ElementArray<Edge> face_edges = kt->getEdges();
for(ElementArray<Edge>::iterator mt = face_edges.begin(); mt != face_edges.end(); ++mt)
segments.push_back(segment(mt->getBeg().RealArray(proj_coords).data(),mt->getEnd().RealArray(proj_coords).data()));
//check that the node sees centers of all segments
for(ElementArray<Edge>::iterator mt = face_edges.begin(); mt != face_edges.end(); ++mt)
{
if( mt->getBeg() != it->self() && mt->getEnd() != it->self() ) //skip adjacent to the node
{
scnt[0] = (mt->getBeg()->RealArray(proj_coords)[0]+mt->getEnd()->RealArray(proj_coords)[0])*0.5;
scnt[1] = (mt->getBeg()->RealArray(proj_coords)[1]+mt->getEnd()->RealArray(proj_coords)[1])*0.5;
joints.push_back(segment(it->RealArray(proj_coords).data(),scnt));
}
}
for(ElementArray<Node>::iterator mt = face_nodes.begin(); mt != face_nodes.end(); ++mt)
loop.push_back(node(mt->RealArray(proj_coords).data()));
loop.push_back(loop.front());
if( check_intersect(segments,joints) || check_dropout(loop,joints) )//,kt->LocalID() == 206357) )
fail = true;
segments.clear();
joints.clear();
loop.clear();
}
if( !fail )
{
//check that all the faces, other then shared with the node, are visible from the node
it->Centroid(cnt);
for(ElementArray<Face>::iterator kt = faces.begin(); kt != faces.end() && !fail; ++kt) if( !kt->GetMarker(sharenode) )
{
kt->Centroid(fcnt);
ray[0] = fcnt[0]-cnt[0];
ray[1] = fcnt[1]-cnt[1];
ray[2] = fcnt[2]-cnt[2];
it->CastRay(cnt,ray,hits);
int counter = 0;
for(tiny_map<HandleType,real,16>::iterator qt = hits.begin(); qt != hits.end(); ++qt)
if( qt->second > 1.0e-5 ) counter++;
if( counter > 1 ) //expect only one hit
fail = true;
hits.clear();
}
if( !fail )
{
//select node as the divisor for faces and cell
for(ElementArray<Face>::iterator kt = faces.begin(); kt != faces.end(); ++kt) if( kt->GetMarker(sharenode) && kt->GetGeometricType() != Element::Tri )
kt->Reference(divnode) = *it;
jt->Reference(divnode) = *it;
}
}
}
}
node_faces.RemMarker(sharenode);
it->RemMarker(sharenode);
}
//all faces that do not get divisor node should be attempted to be divided from their center
int nfcenter = 0;
int nffail = 0;
for(Mesh::iteratorFace it = m.BeginFace(); it != m.EndFace(); ++it) if( it->GetGeometricType() != Element::Tet && it->Reference(divnode) == InvalidHandle() )
{
bool success = false;
ElementArray<Node> face_nodes = it->getNodes();
it->UnitNormal(nrm);
it->Centroid(cnt);
d = dot(cnt,nrm);
//project coordinates
for(ElementArray<Node>::iterator mt = face_nodes.begin(); mt != face_nodes.end(); ++mt)
{
mt->Centroid(ncnt);
nd = dot(ncnt,nrm)-d;
//project onto plane
pcnt[0] = ncnt[0] - nrm[0]*nd;
pcnt[1] = ncnt[1] - nrm[1]*nd;
pcnt[2] = ncnt[2] - nrm[2]*nd;
//relative to center
pcnt[0] -= cnt[0];
pcnt[1] -= cnt[1];
pcnt[2] -= cnt[2];
//choose the basis
if( mt == face_nodes.begin() )
{
orthx[0] = pcnt[0];
orthx[1] = pcnt[1];
orthx[2] = pcnt[2];
cross(orthx,nrm,orthy);
}
//compute and record the coords
real_array c = mt->RealArray(proj_coords);
c[0] = dot(orthx,pcnt);
c[1] = dot(orthy,pcnt);
}
//gather all the segments in projected coordinates
ElementArray<Edge> face_edges = it->getEdges();
for(ElementArray<Edge>::iterator mt = face_edges.begin(); mt != face_edges.end(); ++mt)
segments.push_back(segment(mt->getBeg().RealArray(proj_coords).data(),mt->getEnd().RealArray(proj_coords).data()));
//if( Planarity(it->self(),0.0001) )
{
//form a loop out of face nodes
for(ElementArray<Node>::iterator mt = face_nodes.begin(); mt != face_nodes.end() && !success; ++mt)
loop.push_back(node(mt->RealArray(proj_coords).data()));
loop.push_back(loop.front());
//check that some face node is suitable
for(ElementArray<Node>::iterator mt = face_nodes.begin(); mt != face_nodes.end() && !success; ++mt)
{
int id = it->LocalID();
//check that the node sees centers of all segments
for(ElementArray<Edge>::iterator qt = face_edges.begin(); qt != face_edges.end(); ++qt)
{
if( qt->getBeg() != mt->self() && qt->getEnd() != mt->self() )
{
scnt[0] = (qt->getBeg()->RealArray(proj_coords)[0]+qt->getEnd()->RealArray(proj_coords)[0])*0.5;
scnt[1] = (qt->getBeg()->RealArray(proj_coords)[1]+qt->getEnd()->RealArray(proj_coords)[1])*0.5;
joints.push_back(segment(mt->RealArray(proj_coords).data(),scnt));
}
}
if( !check_intersect(segments,joints) && !check_dropout(loop,joints) )//,it->LocalID() == 206357) )
{
it->Reference(divnode) = mt->GetHandle();
success = true;
}
joints.clear();
}
}
//try face center
if( !success )
{