Commit a26f6451 by Kirill Terekhov

Updates

Update geometric algorithms, improve SplitNonplanar tool, introduce
CollapseDegenerate tool (currently only deletes cells).
parent 00f3e3ce
......@@ -7,6 +7,7 @@ add_executable(Dual dual.cpp)
add_executable(Slice slice.cpp)
add_executable(SliceFunc slice_func.cpp)
add_executable(SplitNonplanar split_nonplanar.cpp)
add_executable(CollapseDegenerate collapse_degenerate.cpp)
target_link_libraries(FixFaults inmost)
if(USE_MPI)
......@@ -83,4 +84,15 @@ endif(USE_MPI)
install(TARGETS SplitNonplanar EXPORT inmost-targets RUNTIME DESTINATION bin)
target_link_libraries(CollapseDegenerate inmost)
if(USE_MPI)
message("linking CollapseDegenerate with MPI")
target_link_libraries(CollapseDegenerate ${MPI_LIBRARIES})
if(MPI_LINK_FLAGS)
set_target_properties(CollapseDegenerate PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
endif()
endif(USE_MPI)
install(TARGETS CollapseDegenerate EXPORT inmost-targets RUNTIME DESTINATION bin)
#include <stdio.h>
#include "inmost.h"
using namespace INMOST;
typedef Storage::real real;
typedef Storage::real_array real_array;
int main(int argc, char ** argv)
{
if( argc < 2 )
{
printf("Usage: %s input_mesh [output_mesh]\n",argv[0]);
return -1;
}
Mesh A("");
A.Load(argv[1]);
Mesh::GeomParam table;
table[ORIENTATION] = FACE;
table[BARYCENTER] = FACE | CELL | EDGE;
A.PrepareGeometricData(table);
rMatrix n(3,1), x(3,1), E(3,3), I(3,3), cx(3,1);
I = rMatrix::Unit(3);
double vol;
int collapsed = 0;
MarkerType rev = A.CreatePrivateMarker();
for(Mesh::iteratorCell it = A.BeginCell(); it != A.EndCell(); ++it)
{
ElementArray<Face> faces = it->getFaces();
// sum(x*n^T - (n^T*x)/3.0*I) = 0
double vol = 0,vol0 = 0;
E.Zero();
it->Centroid(cx.data());
A.FacesOrientation(faces,rev);
for(ElementArray<Face>::iterator jt = faces.begin(); jt != faces.end(); ++jt)
{
jt->Barycenter(x.data());
jt->Normal(n.data());
n*= jt->GetPrivateMarker(rev) ? -1.0:1.0;
E += (x-cx)*n.Transpose();
vol += n.DotProduct(x-cx)/3.0;
vol0 += n.DotProduct(cx)/3.0;
}
if( vol < 0.0 )
{
vol = -vol;
E = -E;
vol0 = -vol0;
}
faces.RemPrivateMarker(rev);
E = E-vol*I;
if( E.FrobeniusNorm() > 1.0e-4 || vol0 > 1.0e-5)
{
std::cout << "Collapsing cell " << it->LocalID() << " err " << E.FrobeniusNorm() << std::endl;
std::cout << "Volume " << it->Volume() << " vol " << vol << " vol0 " << vol0 << std::endl;
std::cout << "E:" << std::endl;
E.Print();
//replace with node
//it->Centroid(x.data());
it->Delete();
collapsed ++;
}
}
A.ReleasePrivateMarker(rev);
std::cout << "total collapsed: " << collapsed << std::endl;
for(ElementType et = FACE; et > NONE; et = PrevElementType(et) )
for(Mesh::iteratorElement it = A.BeginElement(et); it != A.EndElement(); ++it)
if( it->nbAdjElements(CELL) == 0 ) it->Delete();
if( A.HaveTag("GRIDNAME") )
{
Storage::bulk_array nameA = A.self().BulkArray(A.GetTag("GRIDNAME"));
std::string ins = "_collapse_degenerate";
nameA.insert(nameA.end(),ins.c_str(),ins.c_str()+ins.size());
}
if( argc > 2 )
{
std::cout << "Save to " << argv[2] << std::endl;
A.Save(argv[2]);
}
else
{
std::cout << "Save to out.vtk" << std::endl;
A.Save("out.vtk");
}
return 0;
}
\ No newline at end of file
......@@ -275,8 +275,10 @@ int main(int argc, char ** argv)
A.PrepareGeometricData(table);
Tag proj_coords = A.CreateTag("PROJECTED_COORDS",DATA_REAL,NODE,NONE,2);
MarkerType myedges = A.CreateMarker();
real nrm[3], cnt[3], ncnt[3], scnt[3], pcnt[3], fcnt[3], ray[3], orthx[3], orthy[3], d, nd, c[2];
std::cout << "Start splitting faces" << std::endl;
int nsplit = 0, had_faces = A.NumberOfFaces();
for(Mesh::iteratorFace it = A.BeginFace(); it != A.EndFace(); ++it)
{
if( !it->Planarity() )
......@@ -312,6 +314,7 @@ int main(int argc, char ** argv)
//gather all the segments in projected coordinates
ElementArray<Edge> face_edges = it->getEdges();
face_edges.SetMarker(myedges);
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()));
//form a loop out of face nodes
......@@ -347,17 +350,82 @@ int main(int argc, char ** argv)
{
edge_nodes[1] = kt->self();
std::pair<Edge,bool> edge = A.CreateEdge(edge_nodes);
if( edge.second ) //this is actually a new edge
//if( edge.second ) //this is actually a new edge
if( !edge.first.GetMarker(myedges) )
new_edges.push_back(edge.first);
}
}
int had_id = it->LocalID();
int had_nodes = it->nbAdjElements(NODE);
ElementArray<Face> split_faces = Face::SplitFace(it->self(),new_edges,0);
for(int q = 0; q < split_faces.size(); ++q)
if( !split_faces[q].Planarity() )
{
std::cout << __FILE__ << ":" << __LINE__ << "Face " << split_faces[q].LocalID() << " non-planar after split, nodes " << split_faces[q].nbAdjElements(NODE) << " original face " << had_id << " with " << had_nodes << " nodes, split by " << new_edges.size() << " edges, resulted in " << split_faces.size() << " new faces" << std::endl;
std::cout << "original nodes:";
for(ElementArray<Node>::iterator kt = face_nodes.begin(); kt != face_nodes.end(); ++kt)
std::cout << " " << kt->LocalID();
std::cout << std::endl;
std::cout << "split node " << jt->LocalID() << std::endl;
std::cout << "split edges:";
for(ElementArray<Edge>::iterator kt = new_edges.begin(); kt != new_edges.end(); ++kt)
std::cout << " " << kt->LocalID() << " (" << kt->getBeg().LocalID() << "," << kt->getEnd().LocalID() << ")";
std::cout << std::endl;
}
nsplit++;
success = true;
}
}
face_edges.RemMarker(myedges);
if( !success )
{
std::cout << "Face " << it->LocalID() << ", nodes " << it->nbAdjElements(NODE) << ", testing center" << std::endl;
std::vector<segment> joints;
//check that the node sees centers of all segments
for(int k = 0; k < (int)segments.size(); ++k)
{
const real cnt0[2] = {0.0,0.0};
scnt[0] = (segments[k].beg[0]+segments[k].end[0])*0.5;
scnt[1] = (segments[k].beg[1]+segments[k].end[1])*0.5;
joints.push_back(segment(cnt0,scnt));
}
if( check_intersect(segments,joints) )
{
std::cout << "Center of face " << it->LocalID() << " is not suitable to divide it" << std::endl;
}
else
{
//found a suitable node
ElementArray<Edge> new_edges(&A);
ElementArray<Node> edge_nodes(&A,2);
edge_nodes[0] = A.CreateNode(cnt);
for(ElementArray<Node>::iterator kt = face_nodes.begin(); kt != face_nodes.end(); ++kt)
{
edge_nodes[1] = kt->self();
std::pair<Edge,bool> edge = A.CreateEdge(edge_nodes);
new_edges.push_back(edge.first);
}
ElementArray<Face> split_faces = Face::SplitFace(it->self(),new_edges,0);
int had_id = it->LocalID();
int had_nodes = it->nbAdjElements(NODE);
for(int q = 0; q < split_faces.size(); ++q)
if( !split_faces[q].Planarity() )
std::cout << __FILE__ << ":" << __LINE__ << "Face " << split_faces[q].LocalID() << " non-planar after split, nodes " << split_faces[q].nbAdjElements(NODE) << " original face " << had_id << " with " << had_nodes << " nodes, split by " << new_edges.size() << " edges, resulted in " << split_faces.size() << " new faces" << std::endl;
nsplit++;
success = true;
}
}
if( !success ) std::cout << "Cannot split face " << it->LocalID() << std::endl;
}
}
std::cout <<"total split: " << nsplit << " / " << had_faces << std::endl;
A.ReleaseMarker(myedges);
for(Mesh::iteratorFace it = A.BeginFace(); it != A.EndFace(); ++it)
if( !it->Planarity() )
std::cout << it->LocalID() << " is still non-planar, nodes " << it->nbAdjElements(NODE) << std::endl;
if( A.HaveTag("GRIDNAME") )
{
Storage::bulk_array nameA = A.self().BulkArray(A.GetTag("GRIDNAME"));
......
......@@ -1200,8 +1200,8 @@ namespace INMOST
char have_perm = 0;
std::cout << std::scientific;
bool perform_splitting = false;
bool project_perm = true;
int split_degenerate = 1;
bool project_perm = false;
int split_degenerate = 0;
bool check_topology = false;
int verbosity = 0;
for (INMOST_DATA_ENUM_TYPE k = 0; k < file_options.size(); ++k)
......
......@@ -914,7 +914,8 @@ namespace INMOST
real vol = 0, a, at;
real x[3] = {0,0,0}, n[3] = {0,0,0}, n0[3] = {0,0,0}, s, ss;
real l1[3] = {0,0,0}, l2[3] = {0,0,0};
real nt[3] = {0,0,0};
real nt[3] = {0,0,0}, cx[3] = {0,0,0};
me.Centroid(cx);
for(unsigned j = 0; j < faces.size(); j++)
{
//compute normal to face
......@@ -948,7 +949,7 @@ namespace INMOST
n[q] += nt[q];
//same as faces[j].Centroid(x)
for(int q = 0; q < 3; ++q)
x[q] += at*(v0[q]+v1[q]+v2[q])/3.0;
x[q] += at*((v0[q]-cx[q])+(v1[q]-cx[q])+(v2[q]-cx[q]))/3.0;
a += at;
}
for(int q = 0; q < 3; ++q) x[q] /= a;
......@@ -1020,7 +1021,7 @@ namespace INMOST
{
*ret = 0;
real nt[3] = {0,0,0}, l1[3] = {0,0,0}, l2[3] = {0,0,0};
real c[3] = {0,0,0}, n0[3] = {0,0,0}, ss;
real c[3] = {0,0,0}, n0[3] = {0,0,0}, ss, cx[3] = {0,0,0};
real_array v0 = nodes[0].Coords();
real_array v1 = nodes[1].Coords();
real_array v2 = nodes[2].Coords();
......@@ -1028,6 +1029,7 @@ namespace INMOST
vec_diff(v2,v0,l2,mdim);
vec_cross_product(l1,l2,n0);
real a = 0, at;
Element(this,e).Centroid(cx);
for(int i = 1; i < (int)nodes.size()-1; i++)
{
real_array v1 = nodes[i].Coords();
......@@ -1039,10 +1041,10 @@ namespace INMOST
ss /= fabs(ss);
at = sqrt(vec_dot_product(nt,nt,3))*0.5*ss;
for(int q = 0; q < mdim; ++q)
c[q] += at*(v0[q]+v1[q]+v2[q])/3.0;
c[q] += at*((v0[q]-cx[q])+(v1[q]-cx[q])+(v2[q]-cx[q]))/3.0;
a += at;
}
for(int q = 0; q < mdim; ++q) ret[q] = c[q]/a;
for(int q = 0; q < mdim; ++q) ret[q] = c[q]/a+cx[q];
}
//std::cout << ret[0] << " " << ret[1] << " " << ret[2] << std::endl;
}
......@@ -1062,6 +1064,8 @@ namespace INMOST
real c[3] = {0,0,0}, n[3] = {0,0,0};
real n0[3] = {0,0,0}, ss;
real l1[3] = {0,0,0}, l2[3] = {0,0,0};
real cx[3] = {0,0,0};
me.Centroid(cx);
for(unsigned j = 0; j < faces.size(); j++)
{
//compute normal to face
......@@ -1096,13 +1100,13 @@ namespace INMOST
n[q] += nt[q];
//same as faces[j].Centroid(x)
for(int q = 0; q < 3; ++q)
x[q] += at*(v0[q]+v1[q]+v2[q])/3.0;
x[q] += at*((v0[q]-cx[q])+(v1[q]-cx[q])+(v2[q]-cx[q]))/3.0;
a += at;
//second-order midpoint formula
for(int q = 0; q < 3; ++q)
c[q] += s*nt[q]*(pow(v0[q]+v1[q],2)+pow(v0[q]+v2[q],2)+pow(v1[q]+v2[q],2))/24.0;
c[q] += s*nt[q]*(pow((v0[q]-cx[q])+(v1[q]-cx[q]),2)+pow((v0[q]-cx[q])+(v2[q]-cx[q]),2)+pow((v1[q]-cx[q])+(v2[q]-cx[q]),2))/24.0;
}
for(int q = 0; q < 3; ++q) x[q] /= a;
for(int q = 0; q < 3; ++q) x[q] = x[q]/a;
vol += s*vec_dot_product(x,n,3);
}
if( ornt )
......@@ -1118,7 +1122,7 @@ namespace INMOST
}
vol /= 3.0;
for(int q = 0; q < mdim; ++q)
ret[q] = c[q]/(vol);
ret[q] = c[q]/(vol) + cx[q];
//std::cout << ret[0] << " " << ret[1] << " " << ret[2] << std::endl;
}
}
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment