From 0990df3c54f4c9aace189141405ae7304ae11b76 Mon Sep 17 00:00:00 2001 From: "Morteza H. Siboni" Date: Tue, 5 Jun 2018 14:16:23 -0400 Subject: [PATCH 1/4] Added maRegionCollapse This will be used for removing flat regions near the model boundary. --- ma/CMakeLists.txt | 1 + ma/maRegionCollapse.cc | 384 +++++++++++++++++++++++++++++++++++++++++ ma/maRegionCollapse.h | 45 +++++ 3 files changed, 430 insertions(+) create mode 100644 ma/maRegionCollapse.cc create mode 100644 ma/maRegionCollapse.h diff --git a/ma/CMakeLists.txt b/ma/CMakeLists.txt index 7b2421d6c..fa51053ba 100644 --- a/ma/CMakeLists.txt +++ b/ma/CMakeLists.txt @@ -20,6 +20,7 @@ set(SOURCES maSize.cc maOperator.cc maCollapse.cc + maRegionCollapse.cc maMatchedCollapse.cc maLayerCollapse.cc maMatch.cc diff --git a/ma/maRegionCollapse.cc b/ma/maRegionCollapse.cc new file mode 100644 index 000000000..38e42bd52 --- /dev/null +++ b/ma/maRegionCollapse.cc @@ -0,0 +1,384 @@ +/****************************************************************************** + + Copyright 2013 Scientific Computation Research Center, + Rensselaer Polytechnic Institute. All rights reserved. + + The LICENSE file included with this distribution describes the terms + of the SCOREC Non-Commercial License this program is distributed under. + +*******************************************************************************/ +#include "maRegionCollapse.h" +#include "maAdapt.h" +#include "maShape.h" +#include +#include + +namespace ma { + +void RegionCollapse::Init(Adapt* a, double fa) +{ + adapt = a; + region = 0; + reclassifyEdge = 0; + numBdryFaces = 0; + faces[0] = faces[1] = faces[2] = faces[3] = 0; + flatAngle = fa; +} + +bool RegionCollapse::requestLocality(apf::CavityOp* o) +{ +/* get vertices again since this is sometimes used + before setVerts can work */ + Entity* v[4]; + Mesh* m = adapt->mesh; + m->getDownward(region,0,v); + return o->requestLocality(v,4); +} + +bool RegionCollapse::setRegion(Entity* r) +{ + if (getFlag(adapt,r,DONT_COLLAPSE)) + return false; + region = r; + return true; +} + +static int relativeDirection(Vector n1, Vector n2, double flatAngle) +{ + double cosFlatAngleSquare; + double dot,ll1,ll2; + + cosFlatAngleSquare = cos(flatAngle*0.017453293); // PI/180 + cosFlatAngleSquare = cosFlatAngleSquare * cosFlatAngleSquare; + + dot = n1 * n2; + + ll1 = n1 * n1; + ll2 = n2 * n2; + + if(dot*dot/(ll1*ll2) < cosFlatAngleSquare) return 0; + if(dot > 0) return 1; + return -1; +} + +static Entity* getTetVertOppositeFace(Mesh* m, Entity* r, Entity* f) +{ + Entity* rfs[4]; + m->getDownward(r, 2, rfs); + int index = findIn(rfs, 4, f); + PCU_ALWAYS_ASSERT(index > -1); + + Entity* rvs[4]; + Entity* fvs[3]; + + m->getDownward(r, 0, rvs); + m->getDownward(f, 0, fvs); + + for (int i = 0; i < 4; i++) { + int index = findIn(fvs, 3, rvs[i]); + if (index == -1) + return rvs[i]; + } +} + +Vector faceNormal(Mesh* mesh, Entity* face, Entity* rgn) +{ + Vector xyz[4]; + Vector v01, v02, v03; + Entity* fvs[3]; + mesh->getDownward(face, 0, fvs); + for (int i = 0; i < 3; i++) { + xyz[i] = getPosition(mesh, fvs[i]); + } + + + Entity* vertex = getTetVertOppositeFace(mesh, rgn, face); + xyz[3] = getPosition(mesh, vertex); + + v01 = xyz[1] - xyz[0]; + v02 = xyz[2] - xyz[0]; + v03 = xyz[3] - xyz[0]; + + Vector normal = cross(v01, v02); + if( (normal * v03) > 0 ) + normal = normal * (-1); + + return normal; +} + +bool RegionCollapse::checkGeom() +{ + Vector normal1, normal2; + Mesh* mesh = adapt->mesh; + switch (numBdryFaces) { + case 1: + normal1 = faceNormal(mesh, faces[0], region); + for (int i = 1; i < 4; i++) { + normal2 = faceNormal(mesh, faces[i], region); + if (relativeDirection(normal1, normal2, flatAngle) != -1) + return false; + } + break; + case 2: + normal1 = faceNormal(mesh, faces[0], region); + normal2 = faceNormal(mesh, faces[1], region); + if (relativeDirection(normal1, normal2, flatAngle) != 1) + return false; + normal1 = faceNormal(mesh, faces[2], region); + normal2 = faceNormal(mesh, faces[3], region); + if (relativeDirection(normal1, normal2, flatAngle) != 1) + return false; + break; + case 3: + normal1 = faceNormal(mesh, faces[3], region); + for (int i = 0; i < 3; i++) { + normal2 = faceNormal(mesh, faces[i], region); + if (!relativeDirection(normal1, normal2, flatAngle)) + return false; + } + break; + default: + return false; + } + return true; +} + + +bool RegionCollapse::checkTopo() +{ + + Mesh* mesh = adapt->mesh; + + Model* modelFace = 0; + Model* uniqueModelFace = 0; + numBdryFaces = 0; + + Entity* fs[4]; + mesh->getDownward(region, 2, fs); + int m = 0; + for (int i = 0; i < 4; i++) { + Entity* face = fs[i]; + + int faceModelType = mesh->getModelType(mesh->toModel(face)); + if (faceModelType != 2) { + faces[3-m] = face; + ++m; + continue; + } + + if (mesh->countUpward(face) == 2) { + faces[3-m] = face; + ++m; + continue; + } + + faces[numBdryFaces] = face; + ++numBdryFaces; + + modelFace = mesh->toModel(face); + if (uniqueModelFace) { + if (uniqueModelFace != modelFace) + return false; + } + else + uniqueModelFace = modelFace; + } + + Entity* oppositeEdge = 0; + switch (numBdryFaces) { + case 0: + return false; + case 1: + { + Entity* vert = getTetVertOppositeFace(mesh, region, faces[0]); + if (mesh->getModelType(mesh->toModel(vert)) != 3) + return false; + break; + } + case 2: + { + Entity* f2es[3]; + Entity* f3es[3]; + mesh->getDownward(faces[2], 1, f2es); + mesh->getDownward(faces[3], 1, f3es); + int j; + for (j = 0; j < 3; j++) { + int index = findIn(f3es, 3, f2es[j]); + if (index > -1) + break; + } + reclassifyEdge = f2es[j]; // this is the edge being reclassified + + Entity* f0es[0]; + Entity* f1es[1]; + mesh->getDownward(faces[0], 1, f0es); + mesh->getDownward(faces[1], 1, f1es); + for (j = 0; j < 3; j++) { + int index = findIn(f1es, 3, f0es[j]); + if (index > -1) + break; + } + oppositeEdge = f0es[j]; // this is the edge being deleted + + if (mesh->getModelType(mesh->toModel(reclassifyEdge)) != 3) + return false; + if (mesh->getModelType(mesh->toModel(oppositeEdge)) == 1) + return false; + break; + } + case 3: + { + Entity* f3es[3]; + for (int j = 0; j < 2; j++) { + Entity* fes[3]; + mesh->getDownward(faces[j], 1, fes); + for (int k = 0; k < 3; k++) { + int index = findIn(f3es, 3, fes[k]); + if (index > -1) continue; + if (mesh->getModelType(mesh->toModel(fes[k])) != 2) + return false; + } + } + Entity* vert = getTetVertOppositeFace(mesh, region, vert); + if (mesh->getModelType(mesh->toModel(vert)) != 2) + return false; + if (mesh->countUpward(vert) != 3) + return false; + break; + } + default: + return false; + } + + return true; +} + +static void processNewBdryVert(Mesh* m, Entity* v) +{ + Model* c = m->toModel(v); + int modelType = m->getModelType(c); + PCU_ALWAYS_ASSERT(modelType == 1 || modelType == 2); + + Vector coords = getPosition(m, v); + Vector newCoords; + Vector newParams; + m->getClosestPoint(c, coords, newCoords, newParams); + m->setParam(v, newParams); + // TODO: this must be done via snapping + /* m->setPoint(v, 0, newCoords); */ +} + +void RegionCollapse::apply() +{ + Mesh* mesh = adapt->mesh; + Model* modelFace = mesh->toModel(faces[0]); + + switch (numBdryFaces) { + case 1: + { + // reclassify faces + for (int i = 1; i < 4; i++) { + mesh->setModelEntity(faces[i], modelFace); + // TODO: do we need to change direction of faces? + } + // reclassify edges + Entity* f0es[3]; + mesh->getDownward(faces[0], 1, f0es); + Entity* fes[3]; + mesh->getDownward(faces[1], 1, fes); + for (int i = 0; i < 3; i++) { + Entity* edge = fes[i]; + int index = findIn(f0es, 3, edge); + if (index > -1) continue; + mesh->setModelEntity(edge, modelFace); + } + mesh->getDownward(faces[2], 1, fes); + for (int i = 0; i < 3; i++) { + Entity* edge = fes[i]; + int index = findIn(f0es, 3, edge); + if (index > -1) continue; + mesh->setModelEntity(edge, modelFace); + } + Entity* vert = getTetVertOppositeFace(mesh, region, faces[0]); + mesh->setModelEntity(vert, modelFace); + processNewBdryVert(mesh, vert); + + mesh->destroy(region); + mesh->destroy(faces[0]); + break; + } + case 2: + { + Entity* edgeToRemove; + Entity* f0es[3]; + Entity* f1es[3]; + mesh->getDownward(faces[0], 1, f0es); + mesh->getDownward(faces[1], 1, f1es); + for (int i = 0; i < 3; i++) { + edgeToRemove = f0es[i]; + int index = findIn(f1es, 3, edgeToRemove); + if (index > -1) break; + } + + // reclassify faces + for (int i = 2; i < 4; i++) { + mesh->setModelEntity(faces[i], modelFace); + // TODO: do we need to change direction of faces? + } + + mesh->setModelEntity(reclassifyEdge, modelFace); + + mesh->destroy(region); + mesh->destroy(faces[0]); + mesh->destroy(faces[1]); + mesh->destroy(edgeToRemove); + break; + } + case 3: + { + mesh->setModelEntity(faces[3], modelFace); + // TODO: do we need to change direction of faces? + + Entity* vert = getTetVertOppositeFace(mesh, region, faces[3]); + int n = mesh->countUpward(vert); + PCU_ALWAYS_ASSERT(n == 3); + Entity* edges[3]; + for (int i = 0; i < n; i++) { + edges[i] = mesh->getUpward(vert, i); + } + + + mesh->destroy(region); + mesh->destroy(faces[0]); + mesh->destroy(faces[1]); + mesh->destroy(faces[2]); + mesh->destroy(edges[0]); + mesh->destroy(edges[1]); + mesh->destroy(edges[2]); + mesh->destroy(vert); + break; + } + } +} + +void RegionCollapse::unmark() +{ + clearFlag(adapt,region,COLLAPSE); +} + + +bool setupRegionCollapse(RegionCollapse& rcollapse, Entity* region) +{ + Adapt* adapter = rcollapse.adapt; + PCU_ALWAYS_ASSERT(adapter->mesh->getType(region) == apf::Mesh::TET); + if ( ! rcollapse.setRegion(region)) + return false; + if ( ! rcollapse.checkGeom()) + return false; + if ( ! rcollapse.checkTopo()) + return false; + + return true; +} + +} diff --git a/ma/maRegionCollapse.h b/ma/maRegionCollapse.h new file mode 100644 index 000000000..aadcc7f21 --- /dev/null +++ b/ma/maRegionCollapse.h @@ -0,0 +1,45 @@ +/****************************************************************************** + + Copyright 2013 Scientific Computation Research Center, + Rensselaer Polytechnic Institute. All rights reserved. + + The LICENSE file included with this distribution describes the terms + of the SCOREC Non-Commercial License this program is distributed under. + +*******************************************************************************/ +#ifndef MA_REGION_COLLAPSE_H +#define MA_REGION_COLLAPSE_H + +#include "maAdapt.h" + +namespace apf { +class CavityOp; +} + +namespace ma { + +class Adapt; + +class RegionCollapse +{ + public: + void Init(Adapt* a, double fa); + bool requestLocality(apf::CavityOp* o); + bool setRegion(Entity* r); + bool checkGeom(); + bool checkTopo(); + void apply(); + void unmark(); + Adapt* adapt; + Entity* region; + Entity* reclassifyEdge; + int numBdryFaces; + Entity* faces[4]; + double flatAngle; +}; + +bool setupRegionCollapse(RegionCollapse& rcollapse, Entity* region); + +} + +#endif From 3419ad04d9e2b6070efc1483d7a2347a427c45f8 Mon Sep 17 00:00:00 2001 From: "Morteza H. Siboni" Date: Tue, 5 Jun 2018 14:18:35 -0400 Subject: [PATCH 2/4] Updated pkg_tribits.cmake to include maRegionCollapse.cc --- ma/pkg_tribits.cmake | 1 + 1 file changed, 1 insertion(+) diff --git a/ma/pkg_tribits.cmake b/ma/pkg_tribits.cmake index 4ae30f688..cbcc94195 100644 --- a/ma/pkg_tribits.cmake +++ b/ma/pkg_tribits.cmake @@ -17,6 +17,7 @@ set(SOURCES maSize.cc maOperator.cc maCollapse.cc + maRegionCollapse.cc maMatchedCollapse.cc maLayerCollapse.cc maMatch.cc From 233afd5f1c1380b7832b266714c18075b46ca494 Mon Sep 17 00:00:00 2001 From: "Morteza H. Siboni" Date: Tue, 12 Jun 2018 08:41:25 -0400 Subject: [PATCH 3/4] bug fixed in maRegionCollapse.cc --- ma/maRegionCollapse.cc | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/ma/maRegionCollapse.cc b/ma/maRegionCollapse.cc index 38e42bd52..4408e0be2 100644 --- a/ma/maRegionCollapse.cc +++ b/ma/maRegionCollapse.cc @@ -79,6 +79,7 @@ static Entity* getTetVertOppositeFace(Mesh* m, Entity* r, Entity* f) if (index == -1) return rvs[i]; } + return 0; } Vector faceNormal(Mesh* mesh, Entity* face, Entity* rgn) @@ -209,8 +210,8 @@ bool RegionCollapse::checkTopo() } reclassifyEdge = f2es[j]; // this is the edge being reclassified - Entity* f0es[0]; - Entity* f1es[1]; + Entity* f0es[3]; + Entity* f1es[3]; mesh->getDownward(faces[0], 1, f0es); mesh->getDownward(faces[1], 1, f1es); for (j = 0; j < 3; j++) { @@ -239,7 +240,7 @@ bool RegionCollapse::checkTopo() return false; } } - Entity* vert = getTetVertOppositeFace(mesh, region, vert); + Entity* vert = getTetVertOppositeFace(mesh, region, faces[3]); if (mesh->getModelType(mesh->toModel(vert)) != 2) return false; if (mesh->countUpward(vert) != 3) From f3ce4091ff72dd9d6a6df88edf0ff0369cc8a146 Mon Sep 17 00:00:00 2001 From: "Morteza H. Siboni" Date: Tue, 12 Jun 2018 08:42:00 -0400 Subject: [PATCH 4/4] fixes the hangup in shape correction When adapting meshes on 100+ processors FaceVertFixer in shape correction used to cause hangup. --- ma/maShape.cc | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/ma/maShape.cc b/ma/maShape.cc index 7e715b281..9526b63fa 100644 --- a/ma/maShape.cc +++ b/ma/maShape.cc @@ -294,6 +294,10 @@ class FaceVertFixer : public TetFixerBase edges[0] = 0; edges[1] = 0; edges[2] = 0; + verts[0] = 0; + verts[1] = 0; + verts[2] = 0; + verts[3] = 0; face = 0; oppVert = 0; tet = 0; @@ -308,17 +312,18 @@ class FaceVertFixer : public TetFixerBase are too close, the key edges are those that bound face v(0,1,2) */ apf::findTriDown(mesh,v,edges); - face = apf::findUpward(mesh, apf::Mesh::TRIANGLE, edges); tet = apf::findElement(mesh, apf::Mesh::TET, v); oppVert = v[3]; + verts[0] = v[0]; + verts[1] = v[1]; + verts[2] = v[2]; + verts[3] = v[3]; } virtual bool requestLocality(apf::CavityOp* o) { - /* Request locality for edges (for swaps) and v[3] (for face-split- - collapse, the face is already handled by the edges) */ - bool edgesLocalized = o->requestLocality(edges,3); - bool oppVertLocalized = o->requestLocality(&oppVert,1); - return edgesLocalized && oppVertLocalized; + /* by requesting locality for all the verts we can be sure + * that all the desired entities for this operator are local */ + return o->requestLocality(verts,4); } virtual bool run() { @@ -328,6 +333,7 @@ class FaceVertFixer : public TetFixerBase ++nes; return true; } + face = apf::findUpward(mesh, apf::Mesh::TRIANGLE, edges); if (faceSplitCollapse.run(face, tet)) { ++nfsc; @@ -339,6 +345,7 @@ class FaceVertFixer : public TetFixerBase private: Mesh* mesh; Entity* edges[3]; + Entity* verts[4]; Entity *face, *oppVert; Entity* tet; FaceSplitCollapse faceSplitCollapse;