diff --git a/crv/CMakeLists.txt b/crv/CMakeLists.txt index 992589d23..80ffb6ecc 100644 --- a/crv/CMakeLists.txt +++ b/crv/CMakeLists.txt @@ -1,6 +1,7 @@ # Package sources set(SOURCES crv.cc + crvDBG.cc crvAdapt.cc crvBernstein.cc crvBezier.cc @@ -20,11 +21,20 @@ set(SOURCES crvTables.cc crvQuality.cc crvVtk.cc + LBFGS.cc + crvOptimizations.cc + crvObjectiveFunctions.cc ) # Package headers set(HEADERS crv.h + crvAdapt.h + crvBezier.h + crvQuality.h + crvOptimizations.h + crvObjectiveFunctions.h + LBFGS.h ) # Add the crv library diff --git a/crv/LBFGS.cc b/crv/LBFGS.cc new file mode 100644 index 000000000..9946a339f --- /dev/null +++ b/crv/LBFGS.cc @@ -0,0 +1,189 @@ +#include "LBFGS.h" +#include "crvObjectiveFunctions.h" +#include +#include "apfMatrix.h" + +namespace crv{ + +//void LBFGS::setInitialValue(std::vector x) +//{ +// x0 = x; +//} + +std::vector LBFGS::getCurrentX() +{ + return currentX; +} + +static double dotP(const std::vector &v1, const std::vector &v2) +{ + double sum = 0.0; + for (std::size_t i = 0; i < v1.size(); i++) { + sum += v1[i]*v2[i]; + } + return sum; +} + +double LBFGS::getFvalueAfter() +{ + return fValAfter; +} + +double LBFGS::lineSearch(std::vector &xold, std::vector &g, std::vector &direction, double stpmax) +{ + double alpha = 1.0e-4, tolOndeltaX = 1.0e-8; + int itrs = 2000; + double a, alam, alam2 = 0.0, alamin, b, disc, f2 = 0.0; + double rhs1, rhs2, slope = 0.0, sum = 0.0, temp, test, tmplam; + int n = xold.size(); + + std::vector xnew(xold.size(), 0.0); + sum = sqrt(dotP(direction,direction)); + if (sum > stpmax) { + for (int i = 0; i < n; i++) { + direction[i] = direction[i]*stpmax/sum; + } + } +// for (int i = 0; i < n; i++) slope += g[i]*p[i]; + slope = dotP(g,direction); +//if (slope >= 0.0) std::cout<<"slope positive"< test) test = temp; + } + alamin = tolOndeltaX/test; + alam = 1.0; + + double fold = objFunc->getValue(xold); + + for (int k = 0; k < itrs; k++) { + for (int i =0; i < n; i++) xnew[i] = xold[i] + alam*direction[i]; + + double fnew = objFunc->getValue(xnew); + if (alam < alamin || fnew <= fold + alpha*alam*slope) { + return alam; + } + else { + if (alam == 1.0) + tmplam = -slope/(2.0* (fnew-fold-slope)); + else { + rhs1 = fnew - fold - alam*slope; + rhs2 = f2 - fold - alam2*slope; + a = (rhs1/(alam*alam) - rhs2/(alam2*alam2))/(alam - alam2); + b = (-alam2*rhs1/(alam*alam) + alam*rhs2/(alam2*alam2))/(alam-alam2); + if (a == 0.0) tmplam = -slope/(2.0*b); + else { + disc = b*b - 3.0*a*slope; + if (disc < 0.0) tmplam = 0.5*alam; + else if (b <= 0.0) tmplam = (-b + sqrt(disc))/(3.0*a); + else tmplam = -slope/(b + sqrt(disc)); + } + if (tmplam > 0.5*alam) tmplam = 0.5*alam; + } + } + alam2 = alam; + alam = std::max(tmplam, 0.1*alam); + } + return alam; +} + +void LBFGS::moveArrayToLeft(std::vector a[], int r) +{ + for (int i = 0; i < r-1; i++) + a[i] = a[i+1]; +} + +bool LBFGS::run() +{ + std::vector p(x0.size(), 0.0); + std::vector xs[r]; + std::vector gs[r]; + std::vector gdiffs[r]; + double gammas[r]; + + for(int i = 0; i < r; ++i) + { + xs[i] = std::vector(x0.size(), 0.0); + gs[i] = std::vector(x0.size(), 0.0); + gdiffs[i] = std::vector(x0.size(), 0.0); + gammas[i] = 0.0; + } + + xs[0] = x0; + gs[0] = objFunc->getGrad(x0); + + for (std::size_t i = 0; i < xs[0].size(); i++) p[i] = -gs[0][i]; + + for (int k = 0; k < iter; k++) { + int I = 0; + int J = 0; + if (k+1 < r) { + I = k+1; + J = k; + } + else { + I = r-1; + J = I-1; + moveArrayToLeft(xs, r); + moveArrayToLeft(gs, r); + } + double stpmax = (std::max(sqrt(dotP(p,p)), double(objFunc->getSpaceDim()))); + double lambda = lineSearch(xs[J], gs[J], p, stpmax); + + for (std::size_t j = 0; j < xs[I].size(); j++) + xs[I][j] = xs[J][j] + lambda * p[j]; + + gs[I] = objFunc->getGrad(xs[I]); + + if ( I > 0) { + for (std::size_t jj = 0; jj < gs[I].size(); jj++) + gdiffs[I-1][jj] = gs[I][jj] - gs[I-1][jj]; + } + + if ((dotP(gs[I],gs[I]) < tol) || (dotP(gdiffs[I-1], gdiffs[I-1]) < tol)) { + currentX = xs[I]; + fValAfter = objFunc->getValue(xs[I]); + //std::cout<<"number of LBFGS iterations: "<= 0; --i) { + std::vector s(xs[I].size(), 0.0); + std::vector y(xs[I].size(), 0.0); + for( std::size_t j = 0; j < xs[i+1].size(); j++) { + s[j] = xs[i+1][j] - xs[i][j]; + y[j] = gs[i+1][j] - gs[i][j]; + } + double rho = 1 / dotP(s, y); + gammas[i] = rho * dotP(s, p); + for(std::size_t j = 0; j < p.size(); j++) + p[j] = p[j] - gammas[i]*y[j]; + } + + for (int i = 0; i <= I-1; i++) { + std::vector s(xs[I].size(), 0.0); + std::vector y(xs[I].size(), 0.0); + for( std::size_t j = 0; j < xs[i+1].size(); j++) { + s[j] = xs[i+1][j] - xs[i][j]; + y[j] = gs[i+1][j] - gs[i][j]; + } + double rho = 1 / dotP(s, y); + double phi = rho* dotP(y, p); + for(std::size_t j = 0; j < p.size(); j++) + p[j] = p[j] + s[j]* (gammas[i] - phi); + } + + if (dotP(p, gs[I]) > 0) + for (std::size_t j = 0; j < p.size(); j++) p[j] = -p[j]; + + } + return false; + +} + +} // end of namespace + diff --git a/crv/LBFGS.h b/crv/LBFGS.h new file mode 100644 index 000000000..198e82964 --- /dev/null +++ b/crv/LBFGS.h @@ -0,0 +1,47 @@ + +#ifndef LBFGS_H +#define LBFGS_H + +#include +#include +#include + +namespace crv{ + +// forward declare ObjFunction so the type is known to the following classes +class ObjFunction; + +class LBFGS +{ +public: + LBFGS(double inTol, int inIter, const std::vector &x, ObjFunction *inObjFunc): + tol(inTol), iter(inIter), x0(x), objFunc(inObjFunc) + { + //setInitialValue(x0); + } + + ~LBFGS() {} + +public: + //void setInitialValue(std::vector x); + std::vector getCurrentX(); + double getFvalueAfter(); + double lineSearch(std::vector &xold, std::vector &g, std::vector &direction, double stpmax); + void moveArrayToLeft(std::vector a[], int r); + bool run(); + +public: + double tol; + int iter; + std::vector x0; + ObjFunction *objFunc; + std::vector currentX; + double fValAfter; + +private: + int r = 50; +}; + +} + +#endif diff --git a/crv/crv.h b/crv/crv.h index 60e8c2c1a..7d0ce1224 100644 --- a/crv/crv.h +++ b/crv/crv.h @@ -22,6 +22,14 @@ * \brief the curving functions are contained in this namespace */ namespace crv { + +/** \brief choices of objective function for reshape*/ +enum { + NIJK, // based on Jacobian coefficients + DETJ, // based on Jacobian Determinant Values at Integration Points + DETJNIJK // based on TODO : to be decided +}; + /** \brief actually 1 greater than max order */ static unsigned const MAX_ORDER = 19; diff --git a/crv/crvAdapt.cc b/crv/crvAdapt.cc index 3b6baec15..c3128169f 100644 --- a/crv/crvAdapt.cc +++ b/crv/crvAdapt.cc @@ -163,9 +163,10 @@ ma::Input* configureShapeCorrection( static int fixInvalidElements(crv::Adapt* a) { - a->input->shouldForceAdaptation = true; + a->input->shouldForceAdaptation = false; int count = crv::fixLargeBoundaryAngles(a) - + crv::fixInvalidEdges(a); + + crv::fixInvalidFaces(a) + + crv::fixInvalidEdges(a); int originalCount = count; int prev_count; int i = 0; @@ -174,12 +175,14 @@ static int fixInvalidElements(crv::Adapt* a) break; prev_count = count; count = crv::fixLargeBoundaryAngles(a) - + crv::fixInvalidEdges(a); + + crv::fixInvalidFaces(a) + + crv::fixInvalidEdges(a); ++i; } while(count < prev_count); crv::fixLargeBoundaryAngles(a); ma::clearFlagFromDimension(a,ma::COLLAPSE | ma::BAD_QUALITY,1); + ma::clearFlagFromDimension(a, ma::SNAP, 2); a->input->shouldForceAdaptation = false; return originalCount - count; } @@ -220,7 +223,6 @@ void adapt(ma::Input* in) fixCrvElementShapes(a); } - allowSplitCollapseOutsideLayer(a); if (in->maximumIterations > 0) { fixInvalidElements(a); diff --git a/crv/crvCurveMesh.cc b/crv/crvCurveMesh.cc index 5cc39bbfb..aa1a9ea0b 100644 --- a/crv/crvCurveMesh.cc +++ b/crv/crvCurveMesh.cc @@ -9,8 +9,9 @@ #include "crvBezier.h" #include "crvShape.h" #include "crvSnap.h" - +#include "crvMath.h" #include +#include namespace crv { @@ -193,9 +194,11 @@ bool BezierCurver::run() convertInterpolatingToBezier(); - if( m_mesh->getDimension() >= 2 && m_order == 2){ - ma::Input* shapeFixer = configureShapeCorrection(m_mesh); - crv::adapt(shapeFixer); + if(m_mesh->getDimension() >= 2) { + if (m_order == 2 || m_order == 3) { + ma::Input* shapeFixer = configureShapeCorrection(m_mesh); + crv::adapt(shapeFixer); + } } m_mesh->acceptChanges(); diff --git a/crv/crvDBG.cc b/crv/crvDBG.cc new file mode 100644 index 000000000..5a6ed8a10 --- /dev/null +++ b/crv/crvDBG.cc @@ -0,0 +1,337 @@ +#include "crv.h" +#include "crvQuality.h" +#include "crvDBG.h" +#include + + +namespace crv_dbg +{ + +void printTetNumber(apf::Mesh2* m, apf::MeshEntity* e, const char* numberingName) +{ + apf::Numbering* n = m->findNumbering(numberingName); + if (!n) return; + printf("TET:: %d\n", apf::getNumber(n, e, 0, 0)); +} + +void printEdgeCavityInvalidities(apf::Mesh2* m, apf::MeshEntity* e, apf::Numbering* n) +{ + if (n) + printf("entity %d, ", apf::getNumber(n, e, 0, 0)); + else + printf("entity, "); + printf("on model %d\n", m->getModelTag(m->toModel(e))); + apf::Adjacent upTets; + m->getAdjacent(e, 3, upTets); + for (std::size_t i = 0; i < upTets.getSize(); i++) { + std::vector ai = crv::getAllInvalidities(m, upTets[i]); + for (std::size_t j = 0; j < ai.size(); j++) { + printf("%d ", ai[j]); + } + printf("\n"); + } +} + +void visualizeCavityMesh(apf::Mesh2* m, apf::MeshEntity* ent, + const char* prefix, apf::Numbering* n, int resolution) +{ + int num = 0; + if (n) + num = apf::getNumber(n, ent, 0, 0); + + int dim = m->getDimension(); + apf::Adjacent e; + m->getAdjacent(ent, dim, e); + + + gmi_model* g = gmi_load(".null"); + apf::Mesh2* outMesh = apf::makeEmptyMdsMesh(g, dim, false); + + apf::MeshEntity* newEnt[99]; + + for (std::size_t ii = 0; ii < e.getSize(); ii++) { + // Verts + apf::MeshEntity* vs[12]; + apf::MeshEntity* newVs[12]; + int nv = m->getDownward(e[ii], 0, vs); + for(int i = 0; i < nv; ++i) + { + apf::Vector3 p; + apf::Vector3 param(0., 0., 0.); + m->getPoint(vs[i], 0, p); + newVs[i] = outMesh->createVertex(0, p, param); + } + + // Edges + apf::MeshEntity* es[12]; + apf::MeshEntity* newEs[12]; + int ne = m->getDownward(e[ii], 1, es); + for(int i = 0; i < ne; ++i) + { + apf::MeshEntity* evs[2]; + apf::MeshEntity* new_evs[2]; + m->getDownward(es[i], 0, evs); + for (int j = 0; j < 2; j++) { + new_evs[j] = newVs[apf::findIn(vs, nv, evs[j])]; + } + + newEs[i] = outMesh->createEntity(apf::Mesh::EDGE, 0, new_evs); + } + + // Faces + apf::MeshEntity* fs[12]; + apf::MeshEntity* newFs[12]; + int nf = m->getDownward(e[ii], 2, fs); + for(int i = 0; i < nf; ++i) + { + apf::MeshEntity* fes[3]; + apf::MeshEntity* new_fes[3]; + m->getDownward(fs[i], 1, fes); + for (int j = 0; j < 3; j++) { + new_fes[j] = newEs[apf::findIn(es, ne, fes[j])]; + } + newFs[i] = outMesh->createEntity(apf::Mesh::TRIANGLE, 0, new_fes); + } + + // Regions + apf::MeshEntity* tet = 0; + if (dim == 3) { + tet = outMesh->createEntity(apf::Mesh::TET, 0, newFs); + } + + newEnt[ii] = dim == 2 ? newFs[0] : tet; + + PCU_ALWAYS_ASSERT(m->getType(e[ii]) == outMesh->getType(newEnt[ii])); + outMesh->acceptChanges(); + } + apf::changeMeshShape(outMesh, crv::getBezier(3), true); + outMesh->acceptChanges(); + + for (std::size_t ii = 1; ii < e.getSize(); ii++) { + for (int d = 1; d <= dim; d++) + { + if (!m->getShape()->hasNodesIn(d)) continue; + apf::MeshEntity* eds[12]; + int counter = m->getDownward(e[ii], d, eds); + apf::MeshEntity* new_eds[12]; + outMesh->getDownward(newEnt[ii], d, new_eds); + int non = outMesh->getShape()->countNodesOn(apf::Mesh::simplexTypes[d]); + for(int n = 0; n < counter; ++n) { + for(int i = 0; i < non; ++i) { + apf::Vector3 p; + m->getPoint(eds[n], i, p); + outMesh->setPoint(new_eds[n], i, p); + } + } + } + outMesh->acceptChanges(); + } + + std::stringstream ss; + ss << prefix<< num; + crv::writeCurvedVtuFiles(outMesh, apf::Mesh::TET, resolution, ss.str().c_str()); + crv::writeCurvedVtuFiles(outMesh, apf::Mesh::TRIANGLE, resolution, ss.str().c_str()); + crv::writeCurvedWireFrame(outMesh, resolution, ss.str().c_str()); + + outMesh->destroyNative(); + apf::destroyMesh(outMesh); +} + + +void visualizeIndividualCavityEntities(apf::Mesh2* m, apf::MeshEntity* ent, + const char* prefix, apf::Numbering* n, int resolution) +{ + int num = 0; + if (n) + num = apf::getNumber(n, ent, 0, 0); + + int dim = m->getDimension(); + apf::Adjacent e; + m->getAdjacent(ent, dim, e); + + + gmi_model* g = gmi_load(".null"); + + int nat = e.getSize(); + apf::Mesh2* outMesh[nat]; + for (int i = 0; i < nat; i++) { + outMesh[i] = apf::makeEmptyMdsMesh(g, dim, false); + } + + apf::MeshEntity* newEnt[nat]; + + for (int ii = 0; ii < nat; ii++) { + // Verts + apf::MeshEntity* vs[12]; + apf::MeshEntity* newVs[12]; + int nv = m->getDownward(e[ii], 0, vs); + for(int i = 0; i < nv; ++i) + { + apf::Vector3 p; + apf::Vector3 param(0., 0., 0.); + m->getPoint(vs[i], 0, p); + newVs[i] = outMesh[ii]->createVertex(0, p, param); + } + + // Edges + apf::MeshEntity* es[12]; + apf::MeshEntity* newEs[12]; + int ne = m->getDownward(e[ii], 1, es); + for(int i = 0; i < ne; ++i) + { + apf::MeshEntity* evs[2]; + apf::MeshEntity* new_evs[2]; + m->getDownward(es[i], 0, evs); + for (int j = 0; j < 2; j++) { + new_evs[j] = newVs[apf::findIn(vs, nv, evs[j])]; + } + + newEs[i] = outMesh[ii]->createEntity(apf::Mesh::EDGE, 0, new_evs); + } + // Faces + apf::MeshEntity* fs[12]; + apf::MeshEntity* newFs[12]; + int nf = m->getDownward(e[ii], 2, fs); + for(int i = 0; i < nf; ++i) + { + apf::MeshEntity* fes[3]; + apf::MeshEntity* new_fes[3]; + m->getDownward(fs[i], 1, fes); + for (int j = 0; j < 3; j++) { + new_fes[j] = newEs[apf::findIn(es, ne, fes[j])]; + } + newFs[i] = outMesh[ii]->createEntity(apf::Mesh::TRIANGLE, 0, new_fes); + } + + // Regions + apf::MeshEntity* tet = 0; + if (dim == 3) { + tet = outMesh[ii]->createEntity(apf::Mesh::TET, 0, newFs); + } + + newEnt[ii] = dim == 2 ? newFs[0] : tet; + + PCU_ALWAYS_ASSERT(m->getType(e[ii]) == outMesh[ii]->getType(newEnt[ii])); + outMesh[ii]->acceptChanges(); + + apf::changeMeshShape(outMesh[ii], crv::getBezier(m->getShape()->getOrder()), true); + outMesh[ii]->acceptChanges(); + + + for (int d = 1; d <= dim; d++) + { + if (!m->getShape()->hasNodesIn(d)) continue; + apf::MeshEntity* eds[12]; + int counter = m->getDownward(e[ii], d, eds); + apf::MeshEntity* new_eds[12]; + outMesh[ii]->getDownward(newEnt[ii], d, new_eds); + int non = outMesh[ii]->getShape()->countNodesOn(apf::Mesh::simplexTypes[d]); + for(int n = 0; n < counter; ++n) { + for(int i = 0; i < non; ++i) { + apf::Vector3 p; + m->getPoint(eds[n], i, p); + outMesh[ii]->setPoint(new_eds[n], i, p); + } + } + } + outMesh[ii]->acceptChanges(); + + std::stringstream ss; + ss << prefix<< num << "_TET_"<< ii; + crv::writeCurvedVtuFiles(outMesh[ii], apf::Mesh::TET, resolution, ss.str().c_str()); + crv::writeCurvedVtuFiles(outMesh[ii], apf::Mesh::TRIANGLE, resolution, ss.str().c_str()); + crv::writeCurvedWireFrame(outMesh[ii], resolution, ss.str().c_str()); + + outMesh[ii]->destroyNative(); + apf::destroyMesh(outMesh[ii]); + } +} + +void visualizeTetFaces(apf::Mesh2* m, apf::MeshEntity* e, const char* prefix, int resolution) +{ + if (m->getType(e) != apf::Mesh::TET) + return; + int dim = 3; + + gmi_model* g = gmi_load(".null"); + + apf::Mesh2* outMesh[4]; + for (int i = 0; i < 4; i++) { + outMesh[i] = apf::makeEmptyMdsMesh(g, 2, false); + } + + apf::MeshEntity* face[4]; + apf::MeshEntity* newface[4]; + int nf = m->getDownward(e, 2, face); + + for (int i = 0; i < nf; i++) { + + //Verts + apf::MeshEntity* vs[3]; + apf::MeshEntity* newVs[3]; + int nv = m->getDownward(face[i], 0, vs); + for (int j = 0; j < nv; j++) { + apf::Vector3 p; + apf::Vector3 param(0., 0., 0.); + m->getPoint(vs[j], 0, p); + newVs[j] = outMesh[i]->createVertex(0, p, param); + } + + //Edges + apf::MeshEntity* es[3]; + apf::MeshEntity* newEs[3]; + int ne = m->getDownward(face[i], 1, es); + for (int j = 0; j < ne; j++) { + apf::MeshEntity* evs[2]; + apf::MeshEntity* newEvs[2]; + m->getDownward(es[j], 0, evs); + for (int k = 0; k < 2; k++) { + int kk = apf::findIn(vs,nv, evs[k]); + newEvs[k] = newVs[kk]; + } + newEs[j] = outMesh[i]->createEntity(apf::Mesh::EDGE, 0, newEvs); + } + + //Faces + apf::MeshEntity* fes[3]; + apf::MeshEntity* newFes[3]; + m->getDownward(face[i], 1, fes); + for (int j = 0; j < 3; j++) + newFes[j] = newEs[apf::findIn(es, ne, fes[j])]; + + newface[i] = outMesh[i]->createEntity(apf::Mesh::TRIANGLE, 0, newFes); + + PCU_ALWAYS_ASSERT(m->getType(face[i]) == outMesh[i]->getType(newface[i])); + + outMesh[i]->acceptChanges(); + + apf::changeMeshShape(outMesh[i], crv::getBezier(m->getShape()->getOrder()), true); + outMesh[i]->acceptChanges(); + + for (int d = 1; d < dim; d++) { + if (!m->getShape()->hasNodesIn(d)) continue; + apf::MeshEntity* ent[12]; + int counter = m->getDownward(face[i], d, ent); + apf::MeshEntity* newent[12]; + outMesh[i]->getDownward(newface[i], d, newent); + int non = outMesh[i]->getShape()->countNodesOn(apf::Mesh::simplexTypes[d]); + for (int n = 0; n < counter; n++) { + for (int j = 0; j < non; j++) { + apf::Vector3 p; + m->getPoint(ent[n], j, p); + outMesh[i]->setPoint(newent[n], j, p); + } + } + } + + outMesh[i]->acceptChanges(); + + std::stringstream ss; + ss << prefix << "_Face_"<< i; + crv::writeCurvedVtuFiles(outMesh[i], apf::Mesh::TRIANGLE, resolution, ss.str().c_str()); + crv::writeCurvedWireFrame(outMesh[i], resolution, ss.str().c_str()); + } +} + + +} diff --git a/crv/crvDBG.h b/crv/crvDBG.h new file mode 100644 index 000000000..e5254ba18 --- /dev/null +++ b/crv/crvDBG.h @@ -0,0 +1,24 @@ +#ifndef CRVDBG_H +#define CRVDBG_H + +#include +#include +#include +#include + +namespace crv_dbg +{ +void printTetNumber(apf::Mesh2* m, apf::MeshEntity* e, const char* numberingName); +void printCavityInvalidities(apf::Mesh2* m, apf::MeshEntity* e, + apf::Numbering* n = 0); + +void visualizeCavityMesh(apf::Mesh2* m, apf::MeshEntity* ent, + const char* prefix, apf::Numbering* n = 0, int resolution = 8); + +void visualizeIndividualCavityEntities(apf::Mesh2* m, apf::MeshEntity* ent, + const char* prefix, apf::Numbering* n = 0, int resolution = 8); + +void visualizeTetFaces(apf::Mesh2* m, apf::MeshEntity* e, + const char* prefix, int resolution = 8); +} +#endif diff --git a/crv/crvObjectiveFunctions.cc b/crv/crvObjectiveFunctions.cc new file mode 100644 index 000000000..61ea410be --- /dev/null +++ b/crv/crvObjectiveFunctions.cc @@ -0,0 +1,1000 @@ +#include +#include "crv.h" +#include "crvSnap.h" +#include "crvObjectiveFunctions.h" +#include "crvBezier.h" +#include "crvQuality.h" +#include "crvMath.h" + + +static double getAr(apf::Vector3 p0, apf::Vector3 p1, apf::Vector3 p2) +{ + p1 = p1 - p0; + p2 = p2 - p0; + double area = (apf::cross(p1, p2)).getLength(); + return (area/2.0); +} + +namespace crv +{ + +int InternalEdgeReshapeObjFunc :: getSpaceDim() +{ + return P*d; +} + +double InternalEdgeReshapeObjFunc :: getValue(const vector &x) +{ + setNodes(x); + double sum = 0.0; + if (d == 2) { + apf::Adjacent adjF; + apf::NewArray allNodes; + mesh->getAdjacent(edge, 2, adjF); + for (std::size_t i = 0; i < adjF.getSize(); i++) { + sum = sum + f(mesh, adjF[i], size); + } + restoreInitialNodes(); + } + if (d == 3) { + apf::Adjacent adjT; + apf::NewArray allNodes; + mesh->getAdjacent(edge, 3, adjT); + for (std::size_t i = 0; i < adjT.getSize(); i++) { + sum = sum + f(mesh, adjT[i], size); + } + // TODO: In the original code this line is commented. Why? + restoreInitialNodes(); + } + return sum; +} + + +/* vector InternalEdgeReshapeObjFunc :: getGrad(const vector &_x) */ +/* { */ +/* vector x; */ +/* for (int i = 0; i < _x.size(); i++) */ +/* x.push_back(_x[i]); */ + +/* // TODO: Why these values? */ +/* double eps = 1.0e-5; */ +/* double delta = 1.0; */ +/* double h = eps; */ + +/* vector g; */ +/* double xmx = x[0]; */ +/* double xmn = x[0]; */ +/* double ymx = x[1]; */ +/* double ymn = x[1]; */ +/* double zmx = x[2]; */ +/* double zmn = x[2]; */ +/* double df = 0.0, dff = 0.0, dbb = 0.0; */ + +/* for (std::size_t i = 0; i < x.size(); i+=3) { */ +/* if (x[i] >= xmx) xmx = x[i]; */ +/* if (x[i] <= xmn) xmn = x[i]; */ +/* } */ + +/* for (std::size_t i = 1; i < x.size(); i+=3) { */ +/* if (x[i] >= ymx) ymx = x[i]; */ +/* if (x[i] <= ymn) ymn = x[i]; */ +/* } */ + +/* for (std::size_t i = 2; i < x.size(); i+=3) { */ +/* if (x[i] >= zmx) zmx = x[i]; */ +/* if (x[i] <= zmn) zmn = x[i]; */ +/* } */ + +/* double delx = std::abs(xmx - xmn); */ +/* double dely = std::abs(ymx - ymn); */ +/* double delz = std::abs(zmx - zmn); */ + +/* for (std::size_t i = 0; i < x.size(); i++) { */ +/* if (i % 3 == 0) delta = delx; */ +/* if (i % 3 == 1) delta = dely; */ +/* if (i % 3 == 2) delta = delz; */ + +/* h = eps * delta; */ + +/* x[i] = x[i] + h; */ +/* double ff = getValue(x); */ +/* x[i] = x[i] - h; */ + +/* x[i] = x[i] - h; */ +/* double fb = getValue(x); */ +/* x[i] = x[i] + h; */ + +/* df = (ff - fb)/(2.0 * h); */ +/* g.push_back(df); */ +/* } */ +/* return g; */ +/* } */ + +vector InternalEdgeReshapeObjFunc :: getInitialGuess() +{ + return convertNodeVectorToX(ien, ifn, itn); +} + +void InternalEdgeReshapeObjFunc :: setNodes(const vector &x) +{ + std::vector en;// (ien.begin(), ien.end()); + std::vector fn;// (ifn.begin(), ifn.end()); + std::vector tn;// (itn.begin(), itn.end()); + std::vector nod = convertXtoNodeVector(x); + //blendTris(en, fn); + //from all internal node vector to distict vector of nodes + + int nEN = mesh->getShape()->countNodesOn(mesh->getType(edge)); + for (int i = 0; i getShape()->countNodesOn(mesh->TRIANGLE); + apf::Adjacent adjF; + mesh->getAdjacent(edge, 2, adjF); + + for (std::size_t i = nEN; i < nEN + nFN*adjF.getSize(); i++) + fn.push_back(nod[i]); + + if (d > 2 && P > 3) { + //int nTN = mesh->getShape()->countNodesOn(mesh->TET); + for (std::size_t i = nEN+nFN; i getShape()->countNodesOn(mesh->getType(edge)); + + for (int i = 0; i < numENodes; i++) { + mesh->getPoint(edge, i, intEdgeX); + ien.push_back(intEdgeX); + } +} + +void InternalEdgeReshapeObjFunc :: getInitFaceN() +{ + apf::Adjacent adjF; + apf::Vector3 intFaceX; + mesh->getAdjacent(edge, 2, adjF); + int numFNodes = mesh->getShape()->countNodesOn(mesh->TRIANGLE); + for (std::size_t i = 0; i < adjF.getSize(); i++) { + for (int j = 0; j < numFNodes; j++) { + mesh->getPoint(adjF[i], j, intFaceX); + ifn.push_back(intFaceX); + } + } +} + +void InternalEdgeReshapeObjFunc :: getInitTetN() +{ + apf::Adjacent adjT; + apf::Vector3 intTetX; + mesh->getAdjacent(edge, 3, adjT); + int numTNodes = mesh->getShape()->countNodesOn(mesh->TET); + for (std::size_t i = 0; i < adjT.getSize(); i++) { + for (int j = 0; j < numTNodes; j++) { + mesh->getPoint(adjT[i], j, intTetX); + itn.push_back(intTetX); + } + } +} + +vector InternalEdgeReshapeObjFunc :: convertNodeVectorToX( + const vector &en, + const vector &fn, + const vector &tn) +{ + vector x0; + for (int i = 0; i < P-1; i++) + for (int j = 0; j < 3; j++) + x0.push_back(en[i][j]); + + for (std::size_t i = 0; i < fn.size(); i++) + for (int j = 0; j < 3; j++) + x0.push_back(fn[i][j]); + + // TODO: Why Is There a Need For if conditions? + // Isn't tn.size() == 0 when there are no tets in the mesh? + // or when the order of the mesh is less than 4 + if (d > 2 && P > 3) { + for (std::size_t i = 0; i < tn.size(); i++) + for (int j = 0; j < 3; j++) + x0.push_back(tn[i][j]); + } + return x0; +} + +vector InternalEdgeReshapeObjFunc :: convertXtoNodeVector(const vector &x) +{ + std::vector a; + apf::Vector3 v; + if (d == 2) { + std::size_t num = x.size()/d; // TODO: this seems unsafe + //int numENodes = mesh->getShape()->countNodesOn(mesh->getType(edge)); + //check later for 2D case:: x should not include z coordinate in optim search + for (std::size_t i = 0; i < num; i++) { + v = {x[d*i], x[d*i + 1], 0.0}; + a.push_back(v); + } + } + if (d == 3) { + std::size_t num = x.size()/d; // TODO: this seems unsafe + //int numENodes = mesh->getShape()->countNodesOn(mesh->getType(edge)); + for (std::size_t i = 0; i < num; i++) { + v = {x[d*i], x[d*i + 1], x[d*i + 2]}; + a.push_back(v); + } + } + return a; +} + +void InternalEdgeReshapeObjFunc :: blendTris( + const vector &egn, + vector &faceNodes) +{ + apf::Vector3 xi; + apf::Adjacent adjF; + apf::Adjacent adjE; + std::vector cien (ien.begin(),ien.end()); + + mesh->getAdjacent(edge, 2, adjF); + for (std::size_t i = 0; i < adjF.getSize(); i++) { + mesh->getAdjacent(adjF[i], 1, adjE); + int numFNodes = mesh->getShape()->countNodesOn(mesh->getType(adjF[i])); + for (std::size_t j = 0; j < adjE.getSize(); j++) { + if (adjE[j] == edge) { + int jj = 1; + if ( j == 0) + jj = 2; + else if ( j == 1) + jj = 0; + else + jj = 1; + + for (int k = 0; k < numFNodes; k++) { + getBezierNodeXi(mesh->TRIANGLE, P, k, xi); + for (std::size_t ii = 0; ii < egn.size(); ii++) { + double factor = + 0.5 * (xi[j]/(1-xi[jj])) * binomial(P, ii+1) * intpow(1-xi[jj], ii+1) * intpow(xi[jj], P-ii-1) + + 0.5 * (xi[jj]/(1-xi[j])) * binomial(P, ii+1) * intpow(xi[j], ii+1) * intpow(1-xi[j], P-ii-1); + faceNodes[numFNodes*i+k] = faceNodes[numFNodes*i+k] + (egn[ii] - cien[ii]) * factor; + } + } + } + } + } +} + +void InternalEdgeReshapeObjFunc :: updateNodes( + const vector &ed, + const vector &fa, + const vector &te) +{ + int numENodes = mesh->getShape()->countNodesOn(mesh->getType(edge)); + for (int i = 0; i < numENodes; i++) + mesh->setPoint(edge, i, ed[i]); + + apf::Adjacent adjF; + mesh->getAdjacent(edge, 2, adjF); + for (std::size_t i = 0; i < adjF.getSize(); i++) { + int numFNodes = mesh->getShape()->countNodesOn(mesh->getType(adjF[i])); + for (int j = 0; j < numFNodes; j++) + mesh->setPoint(adjF[i], j, fa[numFNodes*i + j]); + } + + if (d == 3 && P > 3) { + apf::Adjacent adjT; + mesh->getAdjacent(edge, 3, adjT); + for (std::size_t i = 0; i < adjT.getSize(); i++) { + int numTNodes = mesh->getShape()->countNodesOn(mesh->getType(adjT[i])); + for (int j =0; j < numTNodes; j++) + mesh->setPoint(adjT[i], j, te[numTNodes*i + j]); + } + } +} + + + + + +// Boundary Edge Objective Functions Impl + +int BoundaryEdgeReshapeObjFunc :: getSpaceDim() +{ + return P*d; +} + +double BoundaryEdgeReshapeObjFunc :: getValue(const vector &x) +{ + setNodes(x); + + double sum = 0.0; + if (d == 2) { + apf::Adjacent adjF; + apf::NewArray allNodes; + mesh->getAdjacent(edge, 2, adjF); + for (std::size_t i = 0; i < adjF.getSize(); i++) { + sum = sum + f(mesh, adjF[i], size); + } + restoreInitialNodes(); + } + + if (d == 3) { + apf::Adjacent adjT; + apf::NewArray allNodes; + mesh->getAdjacent(edge, 3, adjT); + for (std::size_t i = 0; i < adjT.getSize(); i++) { + sum = sum + f(mesh, adjT[i], size); + } + + double ad = 0.0; + double xr = 1.0; + double xir = 1.0; + int nEN = mesh->getShape()->countNodesOn(mesh->getType(edge)); + double beta = 0.0; + std::vector xs; + xs.clear(); + std::vector xis; + xis.clear(); + xis.push_back(apf::Vector3(-1.0, 0.0, 0.0)); + for (int i = 0; i getType(edge), P, i, currentXi); + xis.push_back(currentXi); + } + xis.push_back(apf::Vector3(+1.0, 0.0, 0.0)); + + apf::MeshElement* me = apf::createMeshElement(mesh, edge); + for (std::size_t i = 0; i < xis.size(); i++) { + apf::Vector3 scord; + apf::mapLocalToGlobal(me, xis[i], scord); + xs.push_back(scord); + } + apf::destroyMeshElement(me); + + for (std::size_t i = 1; i < xs.size()-1; i++) { + + xr = (xs[i] - xs[0]).getLength() / + (xs[xs.size()-1] - xs[0]).getLength(); + xir = (xis[i] - xis[0]).getLength() / + (xis[xs.size()-1] - xis[0]).getLength(); + ad = (1.0*xr/xir - 1.0); //(alpha*alpha); + beta = beta + ad*ad; + //sum = sum + ad*ad; + } + + //sum = sum*(1 + beta); + // apf::destroyElement(Edl); + + apf::Adjacent adjF; + mesh->getAdjacent(edge, 2, adjF); + + std::vector xfs; + std::vector xifs; + double arPhys[3] = {1.0, 1.0, 1.0}; + double arParnt[3] = {1.0, 1.0, 1.0}; + double adf; + double gamma = 0.0; + + for (std::size_t i = 0; i < adjF.getSize(); i++) { + xfs.clear(); + xifs.clear(); + adf = 0.0; + + if (mesh->getModelType(mesh->toModel(adjF[i])) == 2) { + int nFN = mesh->getShape()->countNodesOn(mesh->getType(adjF[i])); + apf::Vector3 faceXi; + for (int j = 0; j < nFN; j++) { + getBezierNodeXi(mesh->getType(adjF[i]), P, j, faceXi); + xifs.push_back(faceXi); + } + xifs.push_back(apf::Vector3(0.0, 0.0, 0.0)); + xifs.push_back(apf::Vector3(1.0, 0.0, 0.0)); + xifs.push_back(apf::Vector3(0.0, 1.0, 0.0)); + + apf::MeshElement* mef = apf::createMeshElement(mesh, adjF[i]); + for (size_t k = 0; k < xifs.size(); k++) { + apf::Vector3 fcord; + apf::mapLocalToGlobal(mef, xifs[k], fcord); + xfs.push_back(fcord); + } + apf::destroyMeshElement(mef); + + double triPhys = getAr(xfs[xifs.size()-1], xfs[xifs.size()-2], xfs[xifs.size()-3]); + double triParnt = getAr(xifs[xifs.size()-1], xifs[xifs.size()-2], xifs[xifs.size()-3]); + + for (int j = 0; j < nFN; j++) { + for (int k = 0; k < 2; k++) { + arPhys[k] = getAr(xfs[j], xfs[xifs.size()-(3-k)], xfs[xifs.size()-(2-k)]); + arParnt[k] = getAr(xifs[j], xifs[xifs.size()-(3-k)], xifs[xifs.size()-(2-k)]); + adf = (1.0*arPhys[k]*triParnt/(arParnt[k]*triPhys) - 1); + gamma = gamma + adf*adf; + } + } + } + } + // double gamma = 0.0; + sum = sum*(1 + beta + 0.3*gamma); + restoreInitialNodes(); + } + return sum; +} + +/* vector BoundaryEdgeReshapeObjFunc :: getGrad(const vector &_x) */ +/* { */ +/* vector x; */ +/* for (int i = 0; i < _x.size(); i++) { */ +/* x[i] = _x[i]; */ +/* } */ + +/* //double fold = getValue(x); */ +/* double eps = 1.0e-5; */ +/* double h = eps; */ +/* vector g; */ +/* for (std::size_t i = 0; i < x.size(); i++) { */ +/* if (std::abs(x[i]) > eps) */ +/* h = eps * std::abs(x[i]); */ +/* else */ +/* h = eps; */ + +/* x[i] = x[i] + h; */ +/* double ff = getValue(x); */ +/* x[i] = x[i] - h; */ + +/* x[i] = x[i] - h; */ +/* double fb = getValue(x); */ +/* x[i] = x[i] + h; */ +/* double df = (ff - fb)/(2.0 * h); */ +/* g.push_back(df); */ +/* } */ +/* return g; */ +/* } */ + +vector BoundaryEdgeReshapeObjFunc :: getInitialGuess() +{ + return getParamCoords(); +} + +void BoundaryEdgeReshapeObjFunc :: setNodes(const vector &x) +{ + vector en; + vector fn; + vector tn; + vector nod = convertParamCoordsToNodeVector(x); + + int nEN = mesh->getShape()->countNodesOn(mesh->getType(edge)); + for (int i = 0; i < nEN; i++) + en.push_back(nod[i]); + + apf::Adjacent adjF; + mesh->getAdjacent(edge, 2, adjF); + int nFN = mesh->getShape()->countNodesOn(mesh->TRIANGLE); + + for (std::size_t i = 0; i < adjF.getSize(); i++) { + for (int j = 0; j < nFN; j++) + fn.push_back(nod[nEN+i*nFN+j]); + } + + if (d > 2 && P > 3) { + for (std::size_t i = nEN+nFN; i getShape()->countNodesOn(mesh->getType(edge)); + + for (int i = 0; i < numENodes; i++) { + mesh->getPoint(edge, i, intEdgeX); + ien.push_back(intEdgeX); + } +} + +void BoundaryEdgeReshapeObjFunc :: getInitFaceN() +{ + apf::Adjacent adjF; + apf::Vector3 intFaceX; + apf::Vector3 ipFN; + mesh->getAdjacent(edge, 2, adjF); + int numFNodes = mesh->getShape()->countNodesOn(mesh->TRIANGLE); + for (std::size_t i = 0; i < adjF.getSize(); i++) { + for (int j = 0; j < numFNodes; j++) { + mesh->getPoint(adjF[i], j, intFaceX); + ifn.push_back(intFaceX); + //ipFN = getInterpolatingPointOnFace(mesh, adjF[i], P, j); + //itpfn.push_back(ipFN); + } + } +} + +void BoundaryEdgeReshapeObjFunc :: getInitTetN() +{ + apf::Adjacent adjT; + apf::Vector3 intTetX; + mesh->getAdjacent(edge, 3, adjT); + int numTNodes = mesh->getShape()->countNodesOn(mesh->TET); + for (std::size_t i = 0; i < adjT.getSize(); i++) { + for (int j = 0; j < numTNodes; j++) { + mesh->getPoint(adjT[i], j, intTetX); + itn.push_back(intTetX); + } + } +} + +vector BoundaryEdgeReshapeObjFunc :: getParamCoords() +{ + apf::Vector3 xi; + apf::Vector3 param; + std::vector xp; + + int numENodes = mesh->getShape()->countNodesOn(mesh->getType(edge)); + for (int i = 0; i < numENodes; i++) { + getBezierNodeXi(mesh->getType(edge), P, i, xi); + transferParametricOnEdgeSplit(mesh, edge, 0.5*(xi[0]+1.0), param); + for (int j = 0; j < 3; j++) { + xp.push_back(param[j]); + } + } + + apf::Vector3 xif; + apf::Vector3 paramf; + apf::Adjacent adjF; + mesh->getAdjacent(edge, 2, adjF); + for (std::size_t i = 0; i < adjF.getSize(); i++) { + int numFNodes = mesh->getShape()->countNodesOn(mesh->getType(adjF[i])); + if (mesh->getModelType(mesh->toModel(adjF[i])) == 2) { + for (int j = 0; j < numFNodes; j++) { + getBezierNodeXi(mesh->getType(adjF[i]), P, j, xif); + transferParametricOnTriSplit(mesh, adjF[i], xif, paramf); + for (int k = 0; k < 3; k++) { + xp.push_back(paramf[k]); + } + } + } + else { + apf::Vector3 intFN; + for (int j = 0; j < numFNodes; j++) { + mesh->getPoint(adjF[i], j, intFN); + for (int k = 0; k < 3; k++) + xp.push_back(intFN[k]); + } + } + } + return xp; +} + +vector BoundaryEdgeReshapeObjFunc :: convertParamCoordsToNodeVector( + const vector &x) +{ + apf::ModelEntity* me = mesh->toModel(edge); + std::vector edn; + std::vector vn = convertXtoNodeVector(x); + int nENodes = mesh->getShape()->countNodesOn(mesh->EDGE); + for (int i = 0; i < nENodes; i++) { + apf::Vector3 coorde; + mesh->snapToModel(me, vn[i], coorde); + edn.push_back(coorde); + } + + apf::Adjacent adjF; + mesh->getAdjacent(edge, 2, adjF); + int numFNTotal = 0; + for (std::size_t i = 0; i < adjF.getSize(); i++) { + int nFNodes = mesh->getShape()->countNodesOn(mesh->getType(adjF[i])); + + if (mesh->getModelType(mesh->toModel(adjF[i])) == 2) { + apf::ModelEntity* mef = mesh->toModel(adjF[i]); + for (int j = 0; j < nFNodes; j++) { + apf::Vector3 coordf; + mesh->snapToModel(mef, vn[nENodes+numFNTotal], coordf); + edn.push_back(coordf); + numFNTotal++; + } + } + else { + for (int j = 0; j < nFNodes; j++) { + edn.push_back(vn[nENodes+numFNTotal]); + numFNTotal++; + } + } + } + return edn; +} + +vector BoundaryEdgeReshapeObjFunc :: convertXtoNodeVector( + const std::vector &x) +{ + std::vector a; + apf::Vector3 v; + + if (d == 3) { + std::size_t num = x.size()/d; + //int numENodes = mesh->getShape()->countNodesOn(mesh->getType(edge)); + for (std::size_t i = 0; i < num; i++) { + v = {x[d*i], x[d*i + 1], x[d*i + 2]}; + a.push_back(v); + } + } + return a; +} + + + +void BoundaryEdgeReshapeObjFunc :: blendTris( + const vector &egn, + vector &faceNodes) +{ + apf::Vector3 xi; + apf::Adjacent adjF; + apf::Adjacent adjE; + + mesh->getAdjacent(edge, 2, adjF); + + for (std::size_t i = 0; i < adjF.getSize(); i++) { + if (mesh->getModelType(mesh->toModel(adjF[i])) == 2) { + mesh->getAdjacent(adjF[i], 1, adjE); + int numFNodes = mesh->getShape()->countNodesOn(mesh->getType(adjF[i])); + + for (std::size_t j = 0; j < adjE.getSize(); j++) { + if (adjE[j] == edge) { + int jj = 1; + + if ( j == 0) + jj = 2; + else if ( j == 1) + jj = 0; + else + jj = 1; + + for (int k = 0; k < numFNodes; k++) { + getBezierNodeXi(mesh->TRIANGLE, P, k, xi); + for (std::size_t ii = 0; ii < egn.size(); ii++) { + double factor = + 0.5 * (xi[j]/(1-xi[jj])) * binomial(P, ii+1) * intpow(1-xi[jj], ii+1) * intpow(xi[jj], P-ii-1) + + 0.5 * (xi[jj]/(1-xi[j])) * binomial(P, ii+1) * intpow(xi[j], ii+1) * intpow(1-xi[j], P-ii-1); + faceNodes[numFNodes*i+k] = faceNodes[numFNodes*i+k] + (egn[ii] - ien[ii])*factor; + } + } + } + } + } + } +} + + +vector BoundaryEdgeReshapeObjFunc :: getFaceControlPointsFromInterpolatingPoints( + apf::MeshEntity* face, + const vector &faceInterpolatingP) +{ + vector faceControlP; + apf::Vector3 xi; + apf::NewArray allCntrlP; + apf::Element* Fa = apf::createElement(mesh->getCoordinateField(), face); + apf::getVectorNodes(Fa, allCntrlP); + + int n = mesh->getShape()->countNodesOn(mesh->getType(face)); + mth::Matrix A(n, n); + mth::Matrix Ainv(n, n); + apf::NewArray rhs(n); + int j = 0; + + for (int i = 0; i < n; i++) { + getBezierNodeXi(apf::Mesh::TRIANGLE, P, i, xi); + rhs[i].zero(); + for (int ii = 0; ii < P+1; ii++) { + for (int jj = 0; jj < P+1-ii; jj++) { + if (ii == 0 || jj == 0 || (ii+jj == P)) { + double bFactor = trinomial(P, ii, jj) * Bijk(ii, jj, P-ii-jj, 1.-xi[0]-xi[1], xi[0], xi[1]); + rhs[i] += allCntrlP[getTriNodeIndex(P, ii, jj)] * bFactor; + } + else { + j = getTriNodeIndex(P, ii, jj) - 3*P; //3P is the total number of nodes on all edges + A(i, j) = trinomial(P, ii, jj) * Bijk(ii, jj, P-ii-jj, 1.-xi[0]-xi[1], xi[0], xi[1]); + } + } + } + rhs[i] = faceInterpolatingP[i] - rhs[i]; + } + apf::destroyElement(Fa); + + if (n > 1) + invertMatrixWithPLU(n, A, Ainv); + else + Ainv(0,0) = 1./A(0,0); + + for (int i = 0; i < n; i++) { + apf::Vector3 fcp(0., 0., 0.); + for (int j = 0; j < n; j++) + fcp += rhs[j]*Ainv(i, j); + faceControlP.push_back(fcp); + } + return faceControlP; +} + +void BoundaryEdgeReshapeObjFunc :: updateNodes( + const vector &ed, + const vector &fa, + const vector &te, + bool isInitialX) +{ + if (!isInitialX) { + apf::NewArray eIntpCords; + apf::Element* Ed = apf::createElement(mesh->getCoordinateField(), edge); + apf::getVectorNodes(Ed, eIntpCords); + apf::NewArray trsCoff; + + int nEN = mesh->getShape()->countNodesOn(mesh->getType(edge)); + apf::NewArray contP(nEN); + + for (int i = 0; i < nEN; i++) { + for (int j = 0; j < 3; j++) + eIntpCords[2+i][j] = ed[i][j]; + } + + getBezierTransformationCoefficients(P, mesh->getType(edge), trsCoff); + crv::convertInterpolationPoints(nEN+2, nEN, eIntpCords, trsCoff, contP); + + for (int i = 0; i < nEN; i++) + mesh->setPoint(edge, i, contP[i]); + + apf::destroyElement(Ed); + } + else { + int nEN = mesh->getShape()->countNodesOn(mesh->getType(edge)); + for (int i = 0; i < nEN; i++) + mesh->setPoint(edge, i, ed[i]); + } + + apf::Adjacent adjF; + mesh->getAdjacent(edge, 2, adjF); + std::vector fNd; + + if (!isInitialX) { + for (std::size_t i = 0; i < adjF.getSize(); i++) { + int numFNodes = mesh->getShape()->countNodesOn(mesh->getType(adjF[i])); + fNd.clear(); + + for (int j = 0; j < numFNodes; j++) + fNd.push_back(fa[numFNodes*i + j]); + + if (mesh->getModelType(mesh->toModel(adjF[i])) == 2) { + std::vector fCN = getFaceControlPointsFromInterpolatingPoints(adjF[i], fNd); + + for (int j = 0; j < numFNodes; j++) + mesh->setPoint(adjF[i], j, fCN[j]); + } + else { + for (int j = 0; j < numFNodes; j++) + mesh->setPoint(adjF[i], j, fa[numFNodes*i + j]); + } + } + } + else { + for (std::size_t i = 0; i < adjF.getSize(); i++) { + int numFNodes = mesh->getShape()->countNodesOn(mesh->getType(adjF[i])); + for (int j = 0; j < numFNodes; j++) + mesh->setPoint(adjF[i], j, fa[numFNodes*i+j]); + } + } + + if (d == 3 && P > 3) { + apf::Adjacent adjT; + mesh->getAdjacent(edge, 3, adjT); + for (std::size_t i = 0; i < adjT.getSize(); i++) { + int numTNodes = mesh->getShape()->countNodesOn(mesh->getType(adjT[i])); + for (int j =0; j < numTNodes; j++) + mesh->setPoint(adjT[i], j, te[numTNodes*i + j]); + } + } + + apf::synchronize(mesh->getCoordinateField()); +} + + + + + +// Face Reshape Objective Function +int FaceReshapeObjFunc :: getSpaceDim() +{ + return P*d; +} + +double FaceReshapeObjFunc :: getValue(const std::vector &x) +{ + setNodes(x); + + double sum = 0.0; + + if (d == 3) { + apf::Adjacent adjT; + apf::NewArray allNodes; + mesh->getAdjacent(face, 3, adjT); + for (std::size_t i = 0; i < adjT.getSize(); i++) { + sum = sum + f(mesh, adjT[i], size); + } + restoreInitialNodes(); + } + return sum; +} + +/* vector FaceReshapeObjFunc :: getGrad(const vector &_x) */ +/* { */ +/* vector x; */ +/* for (int i = 0; i < _x.size(); i++) { */ +/* x[i] = _x[i]; */ +/* } */ + +/* //double fold = getValue(x); */ +/* // TODO : Why the following value */ +/* double eps = 1.0e-4; */ +/* double h = eps; */ +/* std::vector g; */ +/* double xmx = x[0]; */ +/* double xmn = x[0]; */ +/* double ymx = x[1]; */ +/* double ymn = x[1]; */ +/* double zmx = x[2]; */ +/* double zmn = x[2]; */ +/* double df = 0.0, dff = 0.0, dbb = 0.0; */ + +/* for (std::size_t i = 0; i < x.size(); i+=3) { */ +/* if (x[i] >= xmx) xmx = x[i]; */ +/* if (x[i] <= xmn) xmn = x[i]; */ +/* } */ + +/* for (std::size_t i = 1; i < x.size(); i+=3) { */ +/* if (x[i] >= ymx) ymx = x[i]; */ +/* if (x[i] <= ymn) ymn = x[i]; */ +/* } */ + +/* for (std::size_t i = 2; i < x.size(); i+=3) { */ +/* if (x[i] >= zmx) zmx = x[i]; */ +/* if (x[i] <= zmn) zmn = x[i]; */ +/* } */ + +/* double delx = std::abs(xmx - xmn); */ +/* double dely = std::abs(ymx - ymn); */ +/* double delz = std::abs(zmx - zmn); */ +/* double delta = 1.0; */ + +/* for (std::size_t i = 0; i < x.size(); i++) { */ +/* if (i % 3 == 0) delta = delx; */ +/* if (i % 3 == 1) delta = dely; */ +/* if (i % 3 == 2) delta = delz; */ + +/* if (delta < eps) */ +/* h = eps * std::abs(x[i]); */ +/* else */ +/* h = eps * delta; */ + +/* x[i] = x[i] + h; */ +/* double ff = getValue(x); */ +/* x[i] = x[i] - h; */ + +/* x[i] = x[i] - h; */ +/* double fb = getValue(x); */ +/* x[i] = x[i] + h; */ + +/* df = (ff - fb)/(2.0 * h); */ +/* g.push_back(df); */ +/* } */ +/* return g; */ +/* } */ + +vector FaceReshapeObjFunc :: getInitialGuess() +{ + return convertNodeVectorToX(ifn); +} + +void FaceReshapeObjFunc :: setNodes(const vector &x) +{ + vector fn;// (ifn.begin(), ifn.end()); + vector tn;// (itn.begin(), itn.end()); + vector nod = convertXtoNodeVector(x); + //blendTris(en, fn); + //from all internal node vector to distict vector of nodes + + int nFN = mesh->getShape()->countNodesOn(mesh->getType(face)); + for (int i = 0; i 2 && P > 3) { + //int nTN = mesh->getShape()->countNodesOn(mesh->TET); + for (std::size_t i = nFN; i getShape()->countNodesOn(mesh->TRIANGLE); + for (int j = 0; j < numFNodes; j++) { + mesh->getPoint(face, j, intFaceX); + ifn.push_back(intFaceX); + } +} + +void FaceReshapeObjFunc :: getInitTetN() +{ + apf::Adjacent adjT; + apf::Vector3 intTetX; + mesh->getAdjacent(face, 3, adjT); + int numTNodes = mesh->getShape()->countNodesOn(mesh->TET); + for (std::size_t i = 0; i < adjT.getSize(); i++) { + for (int j = 0; j < numTNodes; j++) { + mesh->getPoint(adjT[i], j, intTetX); + itn.push_back(intTetX); + } + } +} + +vector FaceReshapeObjFunc :: convertNodeVectorToX(const vector &fn) +{ + vector x0; + + for (std::size_t i = 0; i < fn.size(); i++) + for (int j = 0; j < 3; j++) + x0.push_back(fn[i][j]); + + return x0; +} + +vector FaceReshapeObjFunc :: convertXtoNodeVector(const vector &x) +{ + std::vector a; + apf::Vector3 v; + std::size_t num = x.size()/d; + //int numENodes = mesh->getShape()->countNodesOn(mesh->getType(edge)); + for (std::size_t i = 0; i < num; i++) { + v = {x[d*i], x[d*i + 1], x[d*i + 2]}; + a.push_back(v); + } + return a; +} + +void FaceReshapeObjFunc :: updateNodes( + const vector &fa, + const vector &te) +{ + int numFNodes = mesh->getShape()->countNodesOn(mesh->getType(face)); + for (int j = 0; j < numFNodes; j++) + mesh->setPoint(face, j, fa[j]); + + if (d == 3 && P > 3) { + apf::Adjacent adjT; + mesh->getAdjacent(face, 3, adjT); + for (std::size_t i = 0; i < adjT.getSize(); i++) { + int numTNodes = mesh->getShape()->countNodesOn(mesh->getType(adjT[i])); + for (int j =0; j < numTNodes; j++) + mesh->setPoint(adjT[i], j, te[numTNodes*i + j]); + } + } +} + +} diff --git a/crv/crvObjectiveFunctions.h b/crv/crvObjectiveFunctions.h new file mode 100644 index 000000000..86a551d11 --- /dev/null +++ b/crv/crvObjectiveFunctions.h @@ -0,0 +1,275 @@ +#ifndef CRVOBJECTIVEFUNCTIONS_H +#define CRVOBJECTIVEFUNCTIONS_H + +#include +#include + +#include +#include +#include +#include +#include +#include +#include "crvAdapt.h" + +using namespace std; + +static double getLinearVolPhys(apf::Mesh2* m, apf::MeshEntity* e) +{ + apf::MeshEntity* vs[12]; + int n = m->getDownward(e, 0, vs); + apf::Vector3 coords[12]; + for (int i = 0; i < n; i++) { + m->getPoint(vs[i], 0, coords[i]); + } + + if (m->getType(e) == apf::Mesh::TRIANGLE) + { + return apf::cross(coords[1]-coords[0], coords[2]-coords[0]).getLength()/2.; + } + else if (m->getType(e) == apf::Mesh::TET) + { + apf::Matrix3x3 J; + J[0] = coords[1] - coords[0]; + J[1] = coords[2] - coords[0]; + J[2] = coords[3] - coords[0]; + return apf::getDeterminant(J) / 6.; + } + else + PCU_ALWAYS_ASSERT_VERBOSE(0, + "Not implemented for entities of type other than tri or tet!"); + return 0.; +} + +namespace crv +{ + +class ObjFunction +{ + public: + ObjFunction(){}; + virtual ~ObjFunction(){}; + virtual int getSpaceDim() = 0; + virtual double getTol() = 0; + virtual double getValue(const vector &x) = 0; + /* virtual std::vector getInitialGuess() = 0; */ + /* virtual void setNodes(const vector &x) = 0; */ + /* virtual void restoreInitialNodes() = 0; */ + // TODO :: can we do this once for all the objective functions? + vector getGrad(const vector &_x) + { + double h; + vector x = _x; + double eps = this->getTol(); + vector g; + for (size_t i = 0; i < x.size(); i++) { + h = abs(x[i]) > eps ? eps * abs(x[i]) : eps; + + // forward diff + x[i] += h; + double ff = this->getValue(x); + x[i] -= h; + + // backward diff + x[i] -= h; + double fb = this->getValue(x); + x[i] += h; + + g.push_back( (ff - fb) / 2./ h ); + } + return(g); + } +}; + + +// TODO: can these be reused +/* void getInitEdgeN(); */ +/* void getInitFaceN(); */ +/* void getInitTetN(); */ + +// Internal Edge Objective Functions +class InternalEdgeReshapeObjFunc : public ObjFunction +{ + public: + InternalEdgeReshapeObjFunc( + crv::Adapt* a, + apf::MeshEntity* e, + apf::MeshEntity* t, + double (*fPtr)(apf::Mesh2*, apf::MeshEntity*, ma::SizeField*)) : + adapt(a), edge(e), tet(t), f(fPtr) + { + mesh = adapt->mesh; + size = adapt->sizeField; + P = mesh->getShape()->getOrder(); + d = mesh->getDimension(); + getSpaceDim(); + getInitEdgeN(); + getInitFaceN(); + getInitTetN(); + tol = cbrt(getLinearVolPhys(mesh, t)) * 1.e-3; + } + ~InternalEdgeReshapeObjFunc(){} + int getSpaceDim(); + double getTol() {return tol;} + double getValue(const vector &x); + /* vector getGrad(const vector &x); */ + vector getInitialGuess(); + void setNodes(const vector &x); + void restoreInitialNodes(); + int P; //order + int d; //dimension + private: + void getInitEdgeN(); + void getInitFaceN(); + void getInitTetN(); + vector convertNodeVectorToX( + const vector &en, + const vector &fn, + const vector &tn); + vector convertXtoNodeVector(const vector &x); + void blendTris( + const vector &edgeNodes, + vector &faceNodes); + /* void blendTets( */ + /* const vector &edgeNodes, */ + /* const vector faceNodes, */ + /* vector &tetNodes); */ + void updateNodes( + const vector &ed, + const vector &fa, + const vector &te); + protected: + crv::Adapt* adapt; + apf::MeshEntity* edge; + apf::MeshEntity* tet; + double (*f)(apf::Mesh2*, apf::MeshEntity*, ma::SizeField*); + apf::Mesh2* mesh; + ma::SizeField* size; + vector ien; + vector ifn; + vector itn; + double tol; +}; + + +// Boundary Edge Objective Function +class BoundaryEdgeReshapeObjFunc : public ObjFunction +{ + public: + BoundaryEdgeReshapeObjFunc( + crv::Adapt* a, + apf::MeshEntity* e, + apf::MeshEntity* t, + double (*fPtr)(apf::Mesh2*, apf::MeshEntity*, ma::SizeField*)) : + adapt(a), edge(e), tet(t), f(fPtr) + { + mesh = adapt->mesh; + size = adapt->sizeField; + P = mesh->getShape()->getOrder(); + d = mesh->getDimension(); + getSpaceDim(); + getInitEdgeN(); + getInitFaceN(); + getInitTetN(); + tol = cbrt(getLinearVolPhys(mesh, t)) * 1.e-3; + } + ~BoundaryEdgeReshapeObjFunc(){} + int getSpaceDim(); + double getTol() {return tol;} + double getValue(const vector &x); + /* vector getGrad(const vector &x); */ + vector getInitialGuess(); + void setNodes(const vector &x); + void restoreInitialNodes(); + int P; //order + int d; //dimension + private: + void getInitEdgeN(); + void getInitFaceN(); + void getInitTetN(); + vector getParamCoords(); + vector convertParamCoordsToNodeVector( + const vector &x); + vector convertXtoNodeVector(const vector &x); + void blendTris( + const vector &edgeNodes, + vector &faceNodes); + /* void blendTets( */ + /* const vector &edgeNodes, */ + /* const vector faceNodes, */ + /* vector &tetNodes); */ + vector getFaceControlPointsFromInterpolatingPoints( + apf::MeshEntity* face, + const vector &faceInterpolatingP); + void updateNodes( + const vector &ed, + const vector &fa, + const vector &te, + bool isInitialX); + protected: + crv::Adapt* adapt; + apf::MeshEntity* edge; + apf::MeshEntity* tet; + double (*f)(apf::Mesh2*, apf::MeshEntity*, ma::SizeField*); + apf::Mesh2* mesh; + ma::SizeField* size; + vector ien; + vector ifn; + vector itpfn; + vector itn; + double tol; +}; + +class FaceReshapeObjFunc : public ObjFunction +{ + public: + FaceReshapeObjFunc( + crv::Adapt* a, + apf::MeshEntity* f, + apf::MeshEntity* t, + double (*fPtr)(apf::Mesh2*, apf::MeshEntity*, ma::SizeField*)) : + adapt(a), face(f), tet(t), f(fPtr) + { + mesh = adapt->mesh; + size = adapt->sizeField; + P = mesh->getShape()->getOrder(); + d = mesh->getDimension(); + getSpaceDim(); + getInitFaceN(); + getInitTetN(); + tol = cbrt(getLinearVolPhys(mesh, t)) * 1.e-3; + } + ~FaceReshapeObjFunc() {} + int getSpaceDim(); + double getTol() {return tol;} + double getValue(const vector &x); + /* vector getGrad(const vector &_x); */ + vector getInitialGuess(); + void setNodes(const vector &x); + void restoreInitialNodes(); + int P; //order + int d; //dimension + private: + void getInitFaceN(); + void getInitTetN(); + vector convertNodeVectorToX(const vector &fn); + vector convertXtoNodeVector(const vector &x); + //void blendTets(const std::vector &edgeNodes, const std::vector faceNodes, std::vector &tetNodes); + void updateNodes( + const vector &fa, + const vector &te); + protected: + crv::Adapt* adapt; + apf::MeshEntity* face; + apf::MeshEntity* tet; + double (*f)(apf::Mesh2*, apf::MeshEntity*, ma::SizeField*); + apf::Mesh2* mesh; + ma::SizeField* size; + vector ifn; + vector itn; + double tol; +}; + +} + +#endif diff --git a/crv/crvOptimizations.cc b/crv/crvOptimizations.cc new file mode 100644 index 000000000..cc0a93287 --- /dev/null +++ b/crv/crvOptimizations.cc @@ -0,0 +1,597 @@ +#include +#include +#include +#include +#include +#include "crvOptimizations.h" +#include "crv.h" +#include "crvQuality.h" +#include "crvBezier.h" +#include "crvMath.h" +#include "crvTables.h" +#include "crvDBG.h" +#include + +/* static int global_counter = 0; */ +/* static apf::MeshEntity* tetra[100]; */ +/* static int number = 0; */ + + + +static void printInvalidities(apf::Mesh2* m, apf::Adjacent e, apf::MeshEntity* edge) +{ + return; + int nat = e.getSize(); + apf::Numbering* n = m->findNumbering("debug_num_edge"); + PCU_ALWAYS_ASSERT(n); + int num = apf::getNumber(n, edge, 0, 0); + int tag = m->getModelTag(m->toModel(edge)); + std::cout<<"at edge "<< num <<" tag: "< ai = crv::getAllInvalidities(m, e[i]); + for (std::size_t j = 0; j < ai.size(); j++) { + printf("%d ", ai[j]); + } + printf("\n"); + } +} + +static double computeFValNIJKL(apf::Mesh2* m, apf::MeshEntity* e, ma::SizeField* s = 0) +{ + s = 0; + PCU_ALWAYS_ASSERT_VERBOSE(s == 0, "Not implemented for non-zero sizefield!"); + + int d = m->getDimension(); + int P = m->getShape()->getOrder(); + + apf::NewArray nodes; + apf::Element* el = apf::createElement(m->getCoordinateField(), e); + apf::getVectorNodes(el, nodes); + apf::destroyElement(el); + + double volm = getLinearVolPhys(m, e); + + int weight = 1; + double sumf = 0; + if (d == 3) { + for (int I = 0; I <= d*(P-1); I++) { + for (int J = 0; J <= d*(P-1); J++) { + for (int K = 0; K <= d*(P-1); K++) { + for (int L = 0; L <= d*(P-1); L++) { + if ((I == J && J == K && I == 0) || + (J == K && K == L && J == 0) || + (I == K && K == L && I == 0) || + (I == J && J == L && I == 0)) + weight = 4; + else if ((I == J && I == 0) || + (I == K && I == 0) || + (I == L && I == 0) || + (J == K && J == 0) || + (J == L && J == 0) || + (K == L && K == 0)) + weight = 2; + else + weight = 1; + if (I + J + K + L == d*(P-1)) { + double f = crv::Nijkl(nodes,P,I,J,K)/(6.0*volm) - 1.0; + sumf = sumf + weight*f*f; + } + } + } + } + } + } + + if (d == 2) { + for (int I = 0; I <= d*(P-1); I++) { + for (int J = 0; J <= d*(P-1); J++) { + for (int K = 0; K <= d*(P-1); K++) { + if ((I == J && I == 0) || + (J == K && J == 0) || + (I == K && I == 0)) + weight = 2; + else + weight = 1; + if (I + J + K == d*(P-1)) { + double f = crv::Nijk(nodes,P,I,J)/(4.0*volm) - 1.0; + sumf = sumf + weight*f*f; + } + } + } + } + } + return sumf; +} + +static double computeDetW(apf::MeshElement* me, const apf::Vector3 &xi) +{ + apf::Vector3 xyz; + apf::mapLocalToGlobal(me, xi, xyz); + + //Metric size field + //r1 = 0.0325, r2 = 0.025, h = 0.2 + double tsq = 0.0075 * 0.0075; + double hsq = 0.1 * 0.1; + double ssq = 3.141 * 3.141 * 0.02875 * 0.02875/4.0; + + double rsq = (xyz[0] * xyz[0] + xyz[1] * xyz[1]); + + double xsq = xyz[0]*xyz[0]; + double ysq = xyz[1]*xyz[1]; + + double M11 = (xsq/tsq + ysq/ssq)/rsq; + double M12 = (xyz[0]*xyz[1])*(1./tsq - 1./ssq)/rsq; + double M22 = (ysq/tsq + xsq/ssq)/rsq; + double M33 = 1./hsq; + + apf::Matrix3x3 M(M11, M12, 0., M12, M22, 0., 0., 0., M33); + double detM = apf::getDeterminant(M); + + return (sqrt(detM)); + +} + +static double computeFValDetJNIJKL(apf::Mesh2* m, + apf::MeshEntity* e, ma::SizeField* s = 0) +{ + s = 0; + PCU_ALWAYS_ASSERT_VERBOSE(s == 0, "Not implemented for non-zero sizefield!"); + + int d = m->getDimension(); + int P = m->getShape()->getOrder(); + + apf::NewArray nodes; + apf::Element* el = apf::createElement(m->getCoordinateField(), e); + apf::getVectorNodes(el, nodes); + apf::destroyElement(el); + + apf::MeshElement* me = apf::createMeshElement(m, e); + double volm = getLinearVolPhys(m, e); + + int n = crv::getNumControlPoints(apf::Mesh::TET, d*(P-1)); + apf::NewArray xi; + xi.allocate(n); + + crv::collectNodeXi(apf::Mesh::TET, apf::Mesh::TET, d*(P-1), + crv::elem_vert_xi[apf::Mesh::TET], xi); + + int weight = 1; + double sumf = 0; + if (d == 3) { + for (int I = 0; I <= d*(P-1); I++) { + for (int J = 0; J <= d*(P-1); J++) { + for (int K = 0; K <= d*(P-1); K++) { + for (int L = 0; L <= d*(P-1); L++) { + if ((I == J && J == K && I == 0) || + (J == K && K == L && J == 0) || + (I == K && K == L && I == 0) || + (I == J && J == L && I == 0)) + weight = 4; + else if ((I == J && I == 0) || + (I == K && I == 0) || + (I == L && I == 0) || + (J == K && J == 0) || + (J == L && J == 0) || + (K == L && K == 0)) + weight = 2; + else + weight = 1; + if (I + J + K + L == d*(P-1)) { + double detW = computeDetW(me, xi[crv::getTetNodeIndex(d*(P-1),I,J,K)]); + double f = detW * crv::Nijkl(nodes,P,I,J,K)/(6.0) - 1.0; + sumf = sumf + weight*f*f; + } + } + } + } + } + + apf::destroyMeshElement(me); + } + + if (d == 2) { + for (int I = 0; I <= d*(P-1); I++) { + for (int J = 0; J <= d*(P-1); J++) { + for (int K = 0; K <= d*(P-1); K++) { + if ((I == J && I == 0) || + (J == K && J == 0) || + (I == K && I == 0)) + weight = 2; + else + weight = 1; + if (I + J + K == d*(P-1)) { + double f = crv::Nijk(nodes,P,I,J)/(4.0*volm) - 1.0; + sumf = sumf + weight*f*f; + } + } + } + } + } + return sumf; +} +static double computeFValDetJ(apf::Mesh2* m, apf::MeshEntity* e, ma::SizeField* s) +{ + int order = m->getShape()->getOrder(); + apf::MeshElement* me = apf::createMeshElement(m, e); + apf::Matrix3x3 J; + apf::Matrix3x3 T; + apf::Matrix3x3 Jm; + + double jDet, sum = 0.; + double jDet1, jDet2; + + double volm = getLinearVolPhys(m, e); + for (int i = 0; i < apf::countIntPoints(me, order) ; i++) { + apf::Vector3 qp; + double w = apf::getIntWeight(me, order, i); + apf::getIntPoint(me, order, i, qp); + + apf::getJacobian(me, qp, J); + s->getTransform(me, qp, T); + Jm = J*T; // Jacobian in metric space + jDet = apf::getDeterminant(Jm); + jDet = apf::getDeterminant(J) * computeDetW(me, qp); + + sum += w * (jDet/6.0 - 1.) * (jDet/6.0 - 1.); + apf::Vector3 xyz; + apf::mapLocalToGlobal(me, qp, xyz); + std::cout<<"Xi: "<< qp <<" XYZ: "<getTransform(me, v[i], T); + //Jm = J*T; + jDet1 = apf::getDeterminant(J); + jDet2 = computeDetW(me, v[i]); + jDet = jDet1 * jDet2; + sum += (jDet1/(6.0*volm) - 1.0) * (jDet1/(6.0*volm) - 1.0); + } + apf::destroyMeshElement(me); + + return sum; +} + + +namespace crv{ + +void CrvInternalEdgeOptim :: setMaxIter(int n) +{ + iter = n; +} + +void CrvInternalEdgeOptim :: setTol(double tolerance) +{ + tol = tolerance; +} + +bool CrvInternalEdgeOptim :: run(int &invaliditySize) +{ + bool res = false; + std::vector sizeHolder; + apf::Adjacent adj; + mesh->getAdjacent(edge, 3, adj); + int thisTetSize = 0; + InternalEdgeReshapeObjFunc* objF = 0; + LBFGS* l; + + for (std::size_t i = 0; i < adj.getSize(); i++) { + std::vector ai = crv::getAllInvalidities(mesh, adj[i]); + if (adj[i] == tet) thisTetSize = ai.size(); + sizeHolder.push_back(ai.size()); + + } + + switch (mode) { + case NIJK: + objF = new InternalEdgeReshapeObjFunc(adapt, edge, tet, computeFValNIJKL); + break; + case DETJ: + objF = new InternalEdgeReshapeObjFunc(adapt, edge, tet, computeFValDetJ); + break; + case DETJNIJK: + objF = new InternalEdgeReshapeObjFunc(adapt, edge, tet, computeFValDetJNIJKL); + break; + default: + break; + } + std::vector x0 = objF->getInitialGuess(); + //double f0 = objF->getValue(x0); + //std::cout<< "fval at x0 " << f0<getDownward(adj[i], 1, ed); + int edgeIndex = apf::findIn(ed, 6, edge); + printf("reshape tried on %d edge, TET %d; ", edgeIndex, thisTETnum); + //crv_dbg::printTetNumber(mesh, adj[i], "debug_num_tet"); + } + + bool hasDecreased = false; + invaliditySize = 0; + + if (l->run() && thisTetSize > 0) { + finalX = l->currentX; + fval = l->fValAfter; + objF->setNodes(finalX); + + + for (std::size_t i = 0; i < adj.getSize(); i++) { + std::vector aiNew = crv::getAllInvalidities(mesh, adj[i]); + invaliditySize = invaliditySize + aiNew.size(); + hasDecreased = hasDecreased || (aiNew.size() > (std::size_t)sizeHolder[i]); + } + + if (hasDecreased == false ) { + res = true; + } + else { + //makeIndividualTetsFromFacesOrEdges(mesh, adj_array, edge, "after_cavity_indv_tet_of_edge_", adj.getSize()); + objF->restoreInitialNodes(); + /* printInvalidities(mesh, adj_array, edge, adj.getSize()); */ + std::cout<<"Size DID NOT decrease"<currentX; + //objF->setNodes(finalX); + //makeMultipleEntityMesh(mesh, adj_array, edge, "after_cavity_of_edge_", adj.getSize()); + if (thisTetSize == 0) { + std::cout<<" No Optimization tried"< sizeHolder; + apf::Adjacent adj; + mesh->getAdjacent(edge, 3, adj); + int thisTetSize = 0; + BoundaryEdgeReshapeObjFunc* objF = 0; + LBFGS* l; + + for (std::size_t i = 0; i < adj.getSize(); i++) { + std::vector ai = crv::getAllInvalidities(mesh, adj[i]); + if (adj[i] == tet) thisTetSize = ai.size(); + sizeHolder.push_back(ai.size()); + } + + switch (mode) { + case NIJK: + objF = new BoundaryEdgeReshapeObjFunc(adapt, edge, tet, computeFValNIJKL); + break; + case DETJ: + objF = new BoundaryEdgeReshapeObjFunc(adapt, edge, tet, computeFValDetJ); + break; + case DETJNIJK: + objF = new BoundaryEdgeReshapeObjFunc(adapt, edge, tet, computeFValDetJNIJKL); + break; + default: + break; + } + std::vector x0 = objF->getInitialGuess(); + + + + //double f0 = objF->getValue(x0); + //std::cout<< "fval at x0 " << f0<getDownward(adj[i], 1, ed); + int edgeIndex = apf::findIn(ed, 6, edge); + printf("reshape tried on %d Medge, TET %d ", edgeIndex, thisTETnum); + /* printTetNumber(mesh, adj[i]); */ + } + + bool hasDecreased = false; + invaliditySize = 0; + + if (l->run() && thisTetSize > 0) { + finalX = l->currentX; + fval = l->fValAfter; + objF->setNodes(finalX); + + for (std::size_t i = 0; i < adj.getSize(); i++) { + std::vector aiNew = crv::getAllInvalidities(mesh, adj[i]); + invaliditySize = invaliditySize + aiNew.size(); + hasDecreased = hasDecreased || (aiNew.size() > (std::size_t) sizeHolder[i]); + } + + + if (hasDecreased == false) { + printInvalidities(mesh, adj, edge); + std::cout<<"--------------------------------------"<restoreInitialNodes(); + printInvalidities(mesh, adj, edge); + std::cout<<"size did not decrease"<currentX; + //objF->setNodes(finalX); + + if (thisTetSize == 0) { + std::cout<<"No Optimization tried"< sizeHolder; + apf::Adjacent adj; + mesh->getAdjacent(face, 3, adj); + int thisTetSize = 0; + FaceReshapeObjFunc* objF = 0; + LBFGS* l; + for (std::size_t i = 0; i < adj.getSize(); i++) { + //std::vector ai = crv::getAllInvalidities(mesh, adj[i]); + //if (adj[i] == tet) thisTetSize = ai.size(); + //sizeHolder.push_back(ai.size()); + } + + std::vector ai = crv::getAllInvalidities(mesh, tet); + invaliditySize = ai.size(); + //makeMultipleEntityMesh(mesh, adj_array, face, "before_cavity_of_face_", adj.getSize()); + //makeIndividualTetsFromFacesOrEdges(mesh, adj_array, face, "before_cavity_indv_tet_of_face_", adj.getSize()); + /* printTetNumber(mesh, tet); */ + printInvalidities(mesh, adj, face); + switch (mode) { + case NIJK: + objF = new FaceReshapeObjFunc(adapt, face, tet, computeFValNIJKL); + break; + case DETJ: + objF = new FaceReshapeObjFunc(adapt, face, tet, computeFValDetJ); + break; + case DETJNIJK: + objF = new FaceReshapeObjFunc(adapt, face, tet, computeFValDetJNIJKL); + break; + default: + break; + } + std::vector x0 = objF->getInitialGuess(); + //double f0 = objF->getValue(x0); + //std::cout<< "fval at x0 " << f0<getDownward(adj[i], 2, fc); + int faceIndex = apf::findIn(fc, 4, face); + printf("reshape tried on %d face, TET %d ",faceIndex, thisTETnum); + /* printTetNumber(mesh, adj[i]); */ + } + + //bool hasDecreased = false; + //invaliditySize = 0; + + //if (l->run() && thisTetSize > 0) { + if (l->run() && invaliditySize > 0) { + finalX = l->currentX; + fval = l->fValAfter; + objF->setNodes(finalX); + + std::vector aiNew = crv::getAllInvalidities(mesh, tet); + /* + for (int i = 0; i < adj.getSize(); i++) { + std::vector aiNew = crv::getAllInvalidities(mesh, adj[i]); + invaliditySize = invaliditySize + aiNew.size(); + hasDecreased = hasDecreased || (aiNew.size() > sizeHolder[i]); + } +*/ + // if (hasDecreased == false ) { + if (aiNew.size() <= (std::size_t) invaliditySize) { + invaliditySize = aiNew.size(); + //makeMultipleEntityMesh(mesh, adj_array, face, "after_cavity_of_face_", adj.getSize()); + //makeIndividualTetsFromFacesOrEdges(mesh, adj_array, face, "after_cavity_indv_tet_of_face_", adj.getSize()); + printInvalidities(mesh, adj, face); + std::cout<<"----------------------------------------------------"<restoreInitialNodes(); + printInvalidities(mesh, adj, face); + std::cout<<"Size didNOT decrease"< +#include +#include +#include "apf.h" +#include "apfMesh.h" +#include "apfMesh2.h" +#include "apfShape.h" + +#include "LBFGS.h" +#include "crvObjectiveFunctions.h" + + +namespace crv { + +class CrvEntityOptim +{ + public: + CrvEntityOptim() {} + virtual ~CrvEntityOptim() {} + virtual void setMaxIter(int n) = 0; + virtual void setTol(double tolerance) = 0; + virtual bool run(int &invaliditySize) = 0; + protected: +}; + +class CrvInternalEdgeOptim : public CrvEntityOptim +{ + public: + CrvInternalEdgeOptim( + crv::Adapt* a, + apf::MeshEntity* e, + apf::MeshEntity* t, + int m) : + adapt(a), edge(e), tet(t), mode(m) + { + mesh = adapt->mesh; + } + ~CrvInternalEdgeOptim() + { + } + public: + void setMaxIter(int n); + void setTol(double tolerance); + bool run(int &invaliditySize); + public: + crv::Adapt* adapt; + apf::MeshEntity* edge; + apf::MeshEntity* tet; + int mode; + apf::Mesh2* mesh; + int iter; + double tol; + std::vector finalX; + double fval; +}; + +class CrvBoundaryEdgeOptim : public CrvEntityOptim +{ + public: + CrvBoundaryEdgeOptim( + crv::Adapt* a, + apf::MeshEntity* e, + apf::MeshEntity* t, + int m) : + adapt(a), edge(e), tet(t), mode(m) + { + mesh = adapt->mesh; + } + ~CrvBoundaryEdgeOptim() + { + } + + public: + void setMaxIter(int n); + void setTol(double tolerance); + bool run(int &invaliditySize); + public: + crv::Adapt* adapt; + apf::MeshEntity* edge; + apf::MeshEntity* tet; + int mode; + apf::Mesh2* mesh; + int iter; + double tol; + std::vector finalX; + double fval; +}; + +class CrvFaceOptim : public CrvEntityOptim +{ + public: + CrvFaceOptim( + crv::Adapt* a, + apf::MeshEntity* f, + apf::MeshEntity* t, + int m) : + adapt(a), face(f), tet(t), mode(m) + { + mesh = adapt->mesh; + } + ~CrvFaceOptim() + { + } + + public: + void setMaxIter(int n); + void setTol(double tolerance); + bool run(int &invaliditySize); + public: + crv::Adapt* adapt; + apf::MeshEntity* face; + apf::MeshEntity* tet; + int mode; + apf::Mesh2* mesh; + int iter; + double tol; + std::vector finalX; + double fval; +}; + +} + +#endif diff --git a/crv/crvQuality.cc b/crv/crvQuality.cc index 04119933a..2eefda236 100644 --- a/crv/crvQuality.cc +++ b/crv/crvQuality.cc @@ -10,6 +10,7 @@ #include "crvMath.h" #include "crvTables.h" #include "crvQuality.h" +#include namespace crv { @@ -141,7 +142,7 @@ static double getTetPartialJacobianDet(apf::NewArray& nodes, * 2D planar meshes. * */ -static double Nijk(apf::NewArray& nodes, +double Nijk(apf::NewArray& nodes, int d, int I, int J) { double sum = 0.; @@ -157,7 +158,7 @@ static double Nijk(apf::NewArray& nodes, return sum*d*d/CD; } -static double Nijkl(apf::NewArray& nodes, +double Nijkl(apf::NewArray& nodes, int d, int I, int J, int K) { double sum = 0.; @@ -474,6 +475,7 @@ int Quality3D::checkValidity(apf::MeshEntity* e) int validityTag = computeJacDetNodes(e,nodes,true); if (validityTag > 1) return validityTag; + else return 1; // TODO: Make sure this is actually needed here! // check verts apf::Downward verts; mesh->getDownward(e,0,verts); @@ -614,6 +616,63 @@ double computeTetJacobianDetFromBezierFormulation(apf::Mesh* m, return detJ; } +std::vector getAllInvalidities(apf::Mesh* mesh,apf::MeshEntity* e) +{ + int order = mesh->getShape()->getOrder(); + int n = getNumControlPoints(apf::Mesh::TET, 3*(order-1)); + + apf::NewArray xi; + xi.allocate(n); + + collectNodeXi(apf::Mesh::TET, apf::Mesh::TET, 3*(order-1), + elem_vert_xi[apf::Mesh::TET], xi); + std::vector ai; + apf::NewArray interNodes(n); + apf::MeshElement* me = apf::createMeshElement(mesh,e); + + for (int i = 0; i < 4; ++i){ + interNodes[i] = apf::getDV(me,xi[i]); + if(interNodes[i] < 1e-10){ + ai.push_back(i+2); + } + } + + for (int edge = 0; edge < 6; ++edge){ + for (int i = 0; i < 3*(order-1)-1; ++i){ + int index = 4+edge*(3*(order-1)-1)+i; + interNodes[index] = apf::getDV(me,xi[index]); + if(interNodes[index] < 1e-10){ + ai.push_back(edge+8); + break; + } + } + } + + for (int face = 0; face < 4; ++face){ + for (int i = 0; i < (3*order-4)*(3*order-5)/2; ++i){ + int index = 18*order-20+face*(3*order-4)*(3*order-5)/2+i; + interNodes[index] = apf::getDV(me,xi[index]); + if(interNodes[index] < 1e-10){ + ai.push_back(face+14); + break; + } + } + } + + for (int i = 0; i < (3*order-4)*(3*order-5)*(3*order-6)/6; ++i){ + int index = 18*order*order-36*order+20+i; + interNodes[index] = apf::getDV(me,xi[index]); + if(interNodes[index] < 1e-10){ + ai.push_back(20); + break; + } + } + + apf::destroyMeshElement(me); + + return ai; +} + int Quality3D::computeJacDetNodes(apf::MeshEntity* e, apf::NewArray& nodes, bool validity) { diff --git a/crv/crvQuality.h b/crv/crvQuality.h index e06ad2d1a..7e6f78731 100644 --- a/crv/crvQuality.h +++ b/crv/crvQuality.h @@ -38,6 +38,14 @@ extern const SubdivisionFunction subdivideBezierJacobianDet[apf::Mesh::TYPES]; void elevateBezierJacobianDet(int type, int P, int r, apf::NewArray& nodes, apf::NewArray& elevatedNodes); + +double Nijk(apf::NewArray& nodes, int d, int I, int J); +double Nijkl(apf::NewArray& nodes, int d, int I, int J, int K); + +std::vector getAllInvalidities( + apf::Mesh* mesh, + apf::MeshEntity* e); + } #endif diff --git a/crv/crvShape.cc b/crv/crvShape.cc index 33f421292..4f6d2a9f4 100644 --- a/crv/crvShape.cc +++ b/crv/crvShape.cc @@ -7,11 +7,6 @@ #include #include -#include "crv.h" -#include "crvAdapt.h" -#include "crvShape.h" -#include "crvTables.h" -#include "crvShapeFixer.h" #include #include #include @@ -19,6 +14,13 @@ #include #include #include +#include "crv.h" +#include "crvAdapt.h" +#include "crvQuality.h" +#include "crvShape.h" +#include "crvTables.h" +#include "crvShapeFixer.h" +#include "crvOptimizations.h" /* This is similar to maShape.cc, conceptually, but different enough * that some duplicate code makes sense */ @@ -52,6 +54,361 @@ static bool hasTwoEntitiesOnBoundary(apf::Mesh* m, apf::MeshEntity* e, int dimen /* Mark Edges based on the invalidity code the element has been * tagged with. */ + +static std::vector getEdgeSequenceFromInvalidVertex(ma::Mesh* mesh, ma::Entity* e, int index) +{ + apf::MeshEntity* edges[6]; + mesh->getDownward(e, 1, edges); + + apf::MeshEntity* f[4]; + int nf = mesh->getDownward(e, 2, f); + apf::MeshEntity* vf[3]; + apf::MeshEntity* vt[4]; + mesh->getDownward(e, 0, vt); + + std::vector a; + std::vector b; + + // get the vertex local number for face from the invalid vertex + for (int i = 0; i < nf; i++) { + mesh->getDownward(f[i], 0, vf); + int j = apf::findIn(vf, 3, vt[index]); + + if (j != -1) { + a.push_back(i); + b.push_back(j); + } + } + + //getJacobianDeterminant on each face at + double c[3]; //contains jacobian value + int cc[3]; //face index + apf::Vector3 xi; + apf::Matrix3x3 Jac; + for(size_t k = 0; k < a.size(); k++) { + if (b[k] == 0) xi = {0, 0, 0}; + else if (b[k] == 1) xi = {1, 0, 0}; + else xi = {0, 1, 0}; + + cc[k] = k; + + apf::MeshElement* me = apf::createMeshElement(mesh, f[a[k]]); + apf::getJacobian(me, xi, Jac); + c[k] = apf::getJacobianDeterminant(Jac, 2); + + apf::destroyMeshElement(me); + } + + // sort min to max jacobian determinant + for (int i = 0; i < 3; i++) { + for (int j = i-1; j >= 0; --j) { + double k = c[j+1]; + int kk = cc[j+1]; + if ( k < c[j]) { + c[j+1] = c[j]; + c[j] = k; + cc[j+1] = cc[j]; + cc[j] = kk; + } + } + } + + apf::MeshEntity* ef[3]; + apf::MeshEntity* ve[2]; + std::vector aa; + //std::vector bb; + + // need to flag all adjacent edges + // hence only 2 faces would do + for (int i = 0; i < 2; i++) { + int nef = mesh->getDownward(f[a[cc[i]]], 1, ef); + for (int j = 0; j < nef; j++) { + mesh->getDownward(ef[j], 0, ve); + if (apf::findIn(ve, 2, vt[index]) != -1) { + int k = apf::findIn(edges, 6, ef[j]); + if (k != -1) { + if (aa.size() >= 2) { + bool contains = false; + for (std::size_t ii = 0; ii < aa.size(); ii++) + contains = contains || (k == aa[ii]); + if (contains == false) + aa.push_back(k); + } + else { + aa.push_back(k); + //bb.push_back(mesh->getModelType(mesh->toModel(ef[j]))); + } + } + } + } + } + + // aa has the index of the edges ordered + // min to max of adj face jacobian + + return aa; + +} + +static void sortDesWrtFreqVector(std::vector &f, std::vector &a) +{ + for (std::size_t i = 0; i < f.size(); i++) { + for (int j = i-1; j >= 0; j--) { + int ko = f[j+1]; + int ke = a[j+1]; + if (ko > f[j]) { + f[j+1] = f[j]; + a[j+1] = a[j]; + f[j] = ko; + a[j] = ke; + } + } + } +} + +static std::vector sortEdgeIndexByFrequency(std::vector &all) +{ + std::vector freq, unqEntries; + freq.push_back(1); + unqEntries.push_back(all[0]); + + for (std::size_t i = 1; i < all.size(); i++) { + bool milila = false; + + for (std::size_t j = 0; j < unqEntries.size(); j++) { + milila = milila || (all[i] == unqEntries[j]); + + if (all[i] == unqEntries[j]) { + freq[j] = freq[j] + 1; + } + } + if (milila == false) { + unqEntries.push_back(all[i]); + freq.push_back(1); + } + } + + sortDesWrtFreqVector(freq, unqEntries); + + return unqEntries; +} + +static std::vector numOfUniqueFaces(std::vector &ai) +{ + std::vector aiEdgeNVtx; + std::vector aiOnlyFace; + std::vector aiOnlyEdge; + std::vector aiOnlyVtx; + std::vector aiOnlyTet; + + std::vector aiUniqueFace; + std::vector FFE; + + for (size_t i = 0; i < ai.size(); i++) { + if (ai[i] < 8) aiOnlyVtx.push_back(ai[i]); + else if (ai[i] > 7 && ai[i] < 14) aiOnlyEdge.push_back(ai[i]); + else if (ai[i] > 13 && ai[i] < 20) aiOnlyFace.push_back(ai[i]); + else aiOnlyTet.push_back(ai[i]); + + } + + int faceToEdgeInd[4][4] = {{-1, 0, 1, 2}, + {0, -1, 4, 3}, + {1, 4, -1, 5}, + {2, 3, 5, -1}}; + + int edgeToFaceInd[6][2] = {{0, 1}, + {0, 2}, + {0, 3}, + {1, 3}, + {1, 2}, + {2, 3}}; + + int edgeToVtx[6][2] = {{0, 1}, + {1, 2}, + {0, 2}, + {0, 3}, + {1, 3}, + {2, 3}}; + + for (std::size_t i = 0; i < aiOnlyFace.size(); i++) { + for (std::size_t j = 0; j < aiOnlyFace.size(); j++) { + if ((i != j) && (j > i)) { + for (std::size_t k = 0; k < aiOnlyEdge.size(); k++) { + if (faceToEdgeInd[aiOnlyFace[i]-14][aiOnlyFace[j]-14] + 8 == aiOnlyEdge[k]) { + FFE.push_back(aiOnlyEdge[k]); + FFE.push_back(aiOnlyFace[j]); + FFE.push_back(aiOnlyFace[i]); + break; + } + } + break; + } + } + } + + bool hasDownVtx, alreadyIn; + + for (size_t j = 0; j < aiOnlyEdge.size(); j++) { + hasDownVtx = false; + for (int jj = 0; jj < 2; jj++) { + //hasDownVtx = false; + for (size_t ii = 0; ii < aiOnlyVtx.size(); ii++) { + hasDownVtx = hasDownVtx || (edgeToVtx[aiOnlyEdge[j]-8][jj] == aiOnlyVtx[ii] - 2); + } + } + + if (hasDownVtx == false) { + for (int i = 0; i < 2; i++) { + alreadyIn = false; + for (size_t k = 0; k < aiOnlyFace.size(); k++) { + alreadyIn = alreadyIn || (edgeToFaceInd[aiOnlyEdge[j]-8][i] == aiOnlyFace[k]-14); + } + if (alreadyIn == false) { + aiOnlyFace.push_back(edgeToFaceInd[aiOnlyEdge[j]-8][i] + 14); + } + } + } + else { + for (int i = 0; i < 2; i++) { + for (size_t k = 0; k < aiOnlyFace.size(); k++) { + if (edgeToFaceInd[aiOnlyEdge[j]-8][i] == aiOnlyFace[k]-14) { + //aiOnlyFace.erase(aiOnlyFace.begin()+k); + } + } + } + } + } + + for (size_t j = 0; j < aiOnlyFace.size(); j++) { + aiUniqueFace.push_back(aiOnlyFace[j]); + } + + ai.clear(); + + std::vector forFaceMarking; + + // TODO: to be cleaned up + //[number of unique faces, {unique face index}, {pairs of face face coomon edge invalids}] + + if (aiUniqueFace.size() > 0) { + forFaceMarking.push_back(aiUniqueFace.size()); + for (std::size_t i = 0; i < aiOnlyVtx.size(); i++) + ai.push_back(aiOnlyVtx[i]); + for (std::size_t i = 0; i < aiOnlyEdge.size(); i++) + ai.push_back(aiOnlyEdge[i]); + for (std::size_t i = 0; i < aiOnlyFace.size(); i++) + ai.push_back(aiOnlyFace[i]); + for (std::size_t i = 0; i < aiOnlyTet.size(); i++) + ai.push_back(aiOnlyTet[i]); + for (std::size_t i = 0; i < aiUniqueFace.size(); i++) { + forFaceMarking.push_back(aiUniqueFace[i]); + } + for (std::size_t i = 0; i < FFE.size(); i++) + forFaceMarking.push_back(FFE[i]); + + } + else if (aiUniqueFace.size() == 0) { + forFaceMarking.push_back(0); + if (FFE.size() > 0) { + for (std::size_t i = 0; i < FFE.size(); i++) + forFaceMarking.push_back(FFE[i]); + } + for (std::size_t i = 0; i < aiOnlyVtx.size(); i++) + ai.push_back(aiOnlyVtx[i]); + for (std::size_t i = 0; i < aiOnlyEdge.size(); i++) + ai.push_back(aiOnlyEdge[i]); + for (std::size_t i = 0; i < aiOnlyFace.size(); i++) + ai.push_back(aiOnlyFace[i]); + } + return forFaceMarking; +} + +static int markAllEdges(ma::Mesh* m, ma::Entity* e, + std::vector ai, ma::Entity* edges[6]) +{ + std::vector bb; + apf::MeshEntity* edf[3]; + apf::MeshEntity* faces[4]; + ma::Downward ed; + m->getDownward(e,1,ed); + + //std::vector ai = inspectInvalidies(aiall); + std::vector faceInvalid = numOfUniqueFaces(ai); + + if (ai.size() == 0) return 0; + if (ai.size() == (std::size_t) faceInvalid[0]) return 0; + + for (size_t ii = 0; ii < ai.size(); ii++) { + + int dim = (ai[ii]-2)/6; + int index = (ai[ii]-2) % 6; + + switch (dim) { + case 0: + { + //ma::Downward ed; + //m->getDownward(e,1,ed); + + std::vector aa = getEdgeSequenceFromInvalidVertex(m, e, index); + PCU_ALWAYS_ASSERT(index < 4); + for (size_t i = 0; i < aa.size(); i++) + bb.push_back(aa[i]); + + break; + } + + case 1: + { + //ma::Downward ed; + //m->getDownward(e,1,ed); + bb.push_back(index); + bb.push_back(index); + break; + } + //break; + case 2: + { + // if we have an invalid face, operate on its edges + //ma::Downward edf, faces; + m->getDownward(e,2,faces); + m->getDownward(faces[index],1,edf); + + // TODO : This loop seems to be doing nothing + for (int i = 0; i < 3; i++) { + int j = apf::findIn(ed, 6, edf[i]); + j++; + //if (j != -1) + //bb.push_back(j); + } + break; + } + //break; + case 3: + { + m->getDownward(e,1,edges); + return 6; + } + default: + fail("invalid quality tag in markEdges\n"); + break; + } + } + + int n = 0; + //std::vector allinvEdges = sortEdgeIndexByType(m, e, bb); + std::vector allinvEdges = sortEdgeIndexByFrequency(bb); + + for (size_t i = 0; i < allinvEdges.size(); i++) { + if (m->getModelType(m->toModel(ed[allinvEdges[i]])) != 1) { + n++; + edges[n-1] = ed[allinvEdges[i]]; + } + } + + return n; +} + static int markEdges(ma::Mesh* m, ma::Entity* e, int tag, ma::Entity* edges[6]) { @@ -61,7 +418,6 @@ static int markEdges(ma::Mesh* m, ma::Entity* e, int tag, int index = (tag-2) % 6; int n = 0; int md = m->getDimension(); - switch (dim) { case 0: { @@ -69,6 +425,7 @@ static int markEdges(ma::Mesh* m, ma::Entity* e, int tag, ma::Downward ed; m->getDownward(e,1,ed); n = md; + if(md == 2){ edges[0] = ed[index]; edges[1] = ed[(index+2) % 3]; @@ -78,8 +435,8 @@ static int markEdges(ma::Mesh* m, ma::Entity* e, int tag, edges[1] = ed[tetVertEdges[index][1]]; edges[2] = ed[tetVertEdges[index][2]]; } + break; } - break; case 1: { // if we have a single invalid edge, operate on it @@ -87,8 +444,8 @@ static int markEdges(ma::Mesh* m, ma::Entity* e, int tag, m->getDownward(e,1,ed); edges[0] = ed[index]; n = 1; + break; } - break; case 2: { // if we have an invalid face, operate on its edges @@ -99,12 +456,14 @@ static int markEdges(ma::Mesh* m, ma::Entity* e, int tag, edges[0] = ed[0]; edges[1] = ed[1]; edges[2] = ed[2]; + break; } - break; case 3: + { m->getDownward(e,1,edges); n = 6; break; + } default: fail("invalid quality tag in markEdges\n"); break; @@ -113,6 +472,31 @@ static int markEdges(ma::Mesh* m, ma::Entity* e, int tag, return n; } +static int markUniqueFaces(ma::Mesh* m, ma::Entity* e, std::vector ai, + ma::Entity* faces[4]) +{ + + if (ai.size() == 0) return 0; + std::vector faceInvalid = numOfUniqueFaces(ai); + int n = 0; + + ma::Downward fc; + m->getDownward(e, 2, fc); + ma::Downward ed; + m->getDownward(e, 1, ed); + if (faceInvalid[0] > 0) { + for (int i = 1; i < faceInvalid[0]+1; i++) { + int index = (faceInvalid[i] -2) % 6; + + if (m->getModelType(m->toModel(fc[index])) == 3) { + faces[n] = fc[index]; + n++; + } + } + } + return n; +} + static apf::Vector3 getTriOrientation(apf::Mesh2* m, apf::MeshEntity* tri) { PCU_ALWAYS_ASSERT(m->getType(tri) == apf::Mesh::TRIANGLE); @@ -158,10 +542,22 @@ class EdgeSwapper : public ma::Operator virtual int getTargetDimension() {return md;} virtual bool shouldApply(ma::Entity* e) { + std::vector ai = crv::getAllInvalidities(mesh, e); + mesh->getDownward(e, 1, edges); + int niv = ai.size(); + ne = 0; + if (niv != 0) { + ne = markAllEdges(mesh, e, ai, edges); + simplex = e; + } + return (ne > 0); + /* int tag = crv::getTag(adapter,e); ne = markEdges(mesh,e,tag,edges); + simplex = e; return (ne > 0); + */ } virtual bool requestLocality(apf::CavityOp* o) { @@ -172,7 +568,7 @@ class EdgeSwapper : public ma::Operator for (int i = 0; i < ne; ++i){ if (edgeSwap->run(edges[i])){ ns++; - crv::clearTag(adapter,simplex); + // crv::clearTag(adapter,simplex); ma::clearFlag(adapter,edges[i],ma::COLLAPSE | ma::BAD_QUALITY); break; } @@ -211,6 +607,7 @@ class EdgeReshaper : public ma::Operator virtual bool shouldApply(ma::Entity* e) { int tag = crv::getTag(adapter,e); + ne = markEdges(mesh,e,tag,edges); simplex = e; return (ne > 0); @@ -488,6 +885,160 @@ class EdgeReshaper : public ma::Operator int nr; }; +class FaceOptimizer : public ma::Operator +{ +public: + FaceOptimizer(Adapt* a) { + adapter = a; + mesh = a->mesh; + faces[0] = faces[1] = faces[2] = faces[3] = 0; + simplex = 0; + md = mesh->getDimension(); + numf = 0; + ns = 0; + nf = 0; + } + ~FaceOptimizer() { + } + virtual int getTargetDimension() {return md;} + virtual bool shouldApply(ma::Entity* e) { + + std::vector ai = crv::getAllInvalidities(mesh, e); + int niv = ai.size(); + mesh->getDownward(e, 2, faces); + numf = 0; + if (niv != 0) { + numf = markUniqueFaces(mesh, e, ai, faces); + simplex = e; + } + return (numf > 0); + } + + virtual bool requestLocality(apf::CavityOp* o) + { + return o->requestLocality(faces, numf); + } + + virtual void apply(){ + int invaliditySize = 0; + //mesh->getDownward(simplex, 2, faces); + if (numf == 0) return; + for (int i = 0; i < numf; i++ ) { + if (mesh->getModelType(mesh->toModel(faces[i])) != 3) continue; + //CrvFaceOptim *cfo = new CrvFaceOptim(adapter, faces[i], simplex, DETJ); + CrvFaceOptim *cfo = new CrvFaceOptim(adapter, faces[i], simplex, NIJK); + cfo->setMaxIter(100); + cfo->setTol(1e-8); + if (cfo->run(invaliditySize)) { + if (invaliditySize < 1) { + ns++; + delete cfo; + break; + } + } + else { + nf++; + if (invaliditySize < 1) { + delete cfo; + break; + } + } + delete cfo; + } + } +protected: + Adapt* adapter; + ma::Mesh* mesh; +public: + ma::Entity* faces[4]; + ma::Entity* simplex; + int numf; + int md; + int ns; + int nf; +}; + +class EdgeOptimizer : public ma::Operator +{ +public: + EdgeOptimizer(Adapt* a) { + adapter = a; + mesh = a->mesh; + edges[0] = edges[1] = edges[2] = edges[3] = edges[4] = edges[5] = 0; + simplex = 0; + md = mesh->getDimension(); + ns = 0; + nf = 0; + ne = 0; + } + ~EdgeOptimizer() { + } + virtual int getTargetDimension() {return md;} + virtual bool shouldApply(ma::Entity* e) { + std::vector ai = crv::getAllInvalidities(mesh, e); + mesh->getDownward(e, 1, edges); + int niv = ai.size(); + ne = 0; + if (niv != 0) { + ne = markAllEdges(mesh, e, ai, edges); + simplex = e; + } + return (ne > 0); + } + + virtual bool requestLocality(apf::CavityOp* o) + { + return o->requestLocality(edges, ne); + } + + virtual void apply() { + if (ne == 0) return; + for (int i = 0; i < ne; i++ ) { + int invaliditySize = 0; + CrvEntityOptim* ceo; + + if (mesh->getModelType(mesh->toModel(edges[i])) == 3) { + //ceo = new CrvInternalEdgeOptim(adapter, edges[i], simplex, DETJ); + ceo = new CrvInternalEdgeOptim(adapter, edges[i], simplex, NIJK); + } + else { + //ceo = new CrvBoundaryEdgeOptim(adapter, edges[i], simplex, DETJ); + ceo = new CrvBoundaryEdgeOptim(adapter, edges[i], simplex, NIJK); + } + + ceo->setMaxIter(100); + ceo->setTol(1e-8); + + if (ceo->run(invaliditySize)) { + ns++; + if (invaliditySize < 1) { + break; + delete ceo; + } + } + else { + nf++; + if (invaliditySize < 1) { + delete ceo; + break; + } + } + delete ceo; + } + } +protected: + Adapt* adapter; + ma::Mesh* mesh; + ma::Entity* edge; +public: + int ne; + int md; + ma::Entity* edges[6]; + ma::Entity* simplex; + int ns; + int nf; +}; + static bool isCornerTriAngleLargeMetric(crv::Adapt *a, ma::Entity* tri, int index) { @@ -670,32 +1221,76 @@ static int markEdgesOppLargeAnglesTet(Adapt* a) * and then use the results to mark edges, etc for * fixing */ -static int markEdgesToFix(Adapt* a, int flag) + +/* static int markEdgesToFix(Adapt* a, int flag) */ +/* { */ +/* // do a invalidity check first */ +/* int invalid = markInvalidEntities(a); */ +/* if ( !invalid ) */ +/* return 0; */ +/* int count = 0; */ + +/* ma::Mesh* m = a->mesh; */ +/* ma::Entity* e; */ + +/* // markEdges could have upto 6 edges marked!!! */ +/* ma::Entity* edges[6]; */ +/* ma::Iterator* it = m->begin(m->getDimension()); */ + +/* while ((e = m->iterate(it))) */ +/* { */ +/* //int tag = crv::getTag(a,e); */ +/* //int n = markEdges(m,e,tag,edges); */ + +/* std::vector ai = crv::getAllInvalidities(m, e); */ +/* int niv = ai.size(); */ + +/* if (niv != 0) { */ +/* int n = markAllEdges(m, e, ai, edges); */ +/* for (int i = 0; i < n; ++i){ */ +/* ma::Entity* edge = edges[i]; */ +/* PCU_ALWAYS_ASSERT(edge); */ + +/* if (edge && !ma::getFlag(a,edge,flag)) { */ +/* ma::setFlag(a,edge,flag); */ +/* if (a->mesh->isOwned(edge)) */ +/* ++count; */ +/* } */ +/* } */ +/* } */ +/* } */ +/* m->end(it); */ + +/* return PCU_Add_Long(count); */ +/* } */ + +/* +static int markFacesToFix(Adapt* a, int flag) { - // do a invalidity check first int invalid = markInvalidEntities(a); - if ( !invalid ) + if (!invalid) return 0; int count = 0; ma::Mesh* m = a->mesh; ma::Entity* e; - // markEdges could have upto 6 edges marked!!! - ma::Entity* edges[6]; + ma::Entity* faces[4]; ma::Iterator* it = m->begin(m->getDimension()); while ((e = m->iterate(it))) { - int tag = crv::getTag(a,e); - int n = markEdges(m,e,tag,edges); - for (int i = 0; i < n; ++i){ - ma::Entity* edge = edges[i]; - PCU_ALWAYS_ASSERT(edge); - if (edge && !ma::getFlag(a,edge,flag)) - { - ma::setFlag(a,edge,flag); - if (a->mesh->isOwned(edge)) - ++count; + std::vector ai = crv::getAllInvalidities(m, e); +// int tag = crv::getTag(a, e); +// int n = markFaces(m, e, tag, faces); + int n = markUniqueFaces(m, e, ai, faces); + for (int i = 0; i < n; i++) { + ma::Entity* face = faces[i]; + PCU_ALWAYS_ASSERT(face); + if (face && !ma::getFlag(a, face, flag)) { + if (m->getModelType(m->toModel(face)) != 2) + ma::setFlag(a, face, flag); + if (a->mesh->isOwned(face)) + count++; } } } @@ -703,7 +1298,7 @@ static int markEdgesToFix(Adapt* a, int flag) return PCU_Add_Long(count); } - +*/ int fixLargeBoundaryAngles(Adapt* a) { double t0 = PCU_Time(); @@ -720,7 +1315,7 @@ int fixLargeBoundaryAngles(Adapt* a) return 0; } -static void collapseInvalidEdges(Adapt* a) +static int collapseInvalidEdges(Adapt* a) { double t0 = PCU_Time(); ma::Mesh* m = a->mesh; @@ -737,9 +1332,10 @@ static void collapseInvalidEdges(Adapt* a) double t1 = PCU_Time(); ma::print("Collapsed %d bad edges " "in %f seconds",successCount, t1-t0); + return successCount; } -static void swapInvalidEdges(Adapt* a) +static int swapInvalidEdges(Adapt* a) { double t0 = PCU_Time(); EdgeSwapper es(a); @@ -747,9 +1343,10 @@ static void swapInvalidEdges(Adapt* a) double t1 = PCU_Time(); ma::print("Swapped %d bad edges " "in %f seconds",es.ns, t1-t0); + return es.ns; } -static void repositionInvalidEdges(Adapt* a) +static int repositionInvalidEdges(Adapt* a) { double t0 = PCU_Time(); EdgeReshaper es(a); @@ -757,23 +1354,50 @@ static void repositionInvalidEdges(Adapt* a) double t1 = PCU_Time(); ma::print("Repositioned %d bad edges " "in %f seconds",es.nr, t1-t0); + return es.nr; +} + +static int optimizeInvalidEdges(Adapt* a) +{ + double t0 = PCU_Time(); + EdgeOptimizer eo(a); + ma::applyOperator(a,&eo); + double t1 = PCU_Time(); + ma::print("Optimized %d bad edges, failed %d edges " + "in %f seconds",eo.ns, eo.nf, t1-t0); + return eo.ns; +} + +static int optimizeInvalidFaces(Adapt* a) +{ + double t0 = PCU_Time(); + FaceOptimizer fo(a); + ma::applyOperator(a, &fo); + double t1 = PCU_Time(); + ma::print("Optimized %d bad faces, failed %d faces " + "in %f seconds", fo.ns, fo.nf, t1-t0); + return fo.ns; +} + +int fixInvalidFaces(Adapt* a) +{ + return optimizeInvalidFaces(a); } int fixInvalidEdges(Adapt* a) { - int count = markEdgesToFix(a,ma::BAD_QUALITY | ma::COLLAPSE ); - if (! count){ - return 0; - } + int count = 0; if(a->mesh->getShape()->getOrder() == 2) - repositionInvalidEdges(a); - collapseInvalidEdges(a); - swapInvalidEdges(a); + count += repositionInvalidEdges(a); + else + count += optimizeInvalidEdges(a); + + count += collapseInvalidEdges(a); + count += swapInvalidEdges(a); return count; } - struct IsBadCrvQuality : public ma::Predicate { IsBadCrvQuality(Adapt* a_):a(a_) @@ -847,6 +1471,7 @@ void fixCrvElementShapes(Adapt* a) /* PCU_Add_Ints(&numEdgeRemoved,1); */ /* if (PCU_Comm_Self() == 0) */ /* lion_oprint(1,"==> %d edges removal operations succeeded.\n", numEdgeRemoved); */ + //fixInvalidFaces(a); count = markCrvBadQuality(a); ++i; } while(count < prev_count && i < 6); // the second conditions is to make sure this does not take long diff --git a/crv/crvShape.h b/crv/crvShape.h index fcd3ed09f..2536a13a5 100644 --- a/crv/crvShape.h +++ b/crv/crvShape.h @@ -59,6 +59,10 @@ int fixLargeBoundaryAngles(Adapt* a); try and collapse or swap it away */ int fixInvalidEdges(Adapt* a); +/** \brief If a face is flagged as invalid + optimize the face nodes */ +int fixInvalidFaces(Adapt* a); + /** \brief attempts to fix the shape of the elements in a same manner as ma::fixElementShape */ void fixCrvElementShapes(Adapt* a); diff --git a/crv/crvTables.h b/crv/crvTables.h index 357301193..b31e1cc2a 100644 --- a/crv/crvTables.h +++ b/crv/crvTables.h @@ -43,14 +43,13 @@ static int const tet_tri_edges[4][3] = corresponding to tet_tri_edges, 0 -> it is correctly oriented, 1 -> it is flipped canonically */ static bool const flip_tet_tri_edges[4][3] = -{{0,0,0},{0,0,1},{0,0,1},{1,0,1}}; +{{1,0,0},{0,0,1},{0,0,1},{1,0,1}}; /** \brief edge indices connected to a vertex of a tet, this does not comment on their orientation wrt to the vertex \details ordered as XJ Luo's thesis */ static int const tetVertEdges[4][3] = {{3,0,2},{0,4,1},{1,5,2},{3,5,4}}; - /** \brief edge indices connected to a vertex of a tri, this does not comment on their orientation wrt to the vertex \details [refer to the comment for tets] */ diff --git a/crv/crvVtk.cc b/crv/crvVtk.cc index 489e1f3ae..0964b2e6f 100644 --- a/crv/crvVtk.cc +++ b/crv/crvVtk.cc @@ -311,6 +311,7 @@ static void writeEdgeJacobianDet(std::ostream& file, apf::Mesh* m, int n) file << "\n"; } + static void writeTriJacobianDet(std::ostream& file, apf::Mesh* m, int n) { file << "