diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp index 74eefdf12ae54017b9d254b585a82dcfc579e2f6..67c7cfee75632ad9b5466a4edf87b9804931e13e 100644 --- a/Mesh/HighOrder.cpp +++ b/Mesh/HighOrder.cpp @@ -21,6 +21,7 @@ #include "Context.h" #include "fullMatrix.h" #include "BasisFactory.h" +#include "MVertexRTree.h" // --------- Functions that help optimizing placement of points on geometry ----------- @@ -1307,6 +1308,73 @@ static void updateHighOrderVertices(GEntity *e, e->deleteVertexArrays(); } +static void updatePeriodicEdgesAndFaces(GModel *m) +{ + for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); ++it) { + GEdge *slave = *it; + GEdge *master = dynamic_cast<GEdge*>(slave->meshMaster()); + if (master == slave) continue; + + std::map<MVertex*,MVertex*> &v2v = slave->correspondingVertices; + v2v.clear(); + + MVertexRTree rtree = MVertexRTree(1e-8 * CTX::instance()->lc); + for (unsigned int i = 0; i < slave->getNumMeshVertices(); ++i) + rtree.insert(slave->getMeshVertex(i)); + + std::vector<double> tfo = slave->affineTransform; + + for (unsigned int i = 0; i < master->getNumMeshVertices(); ++i) { + MVertex *vs = master->getMeshVertex(i); + double ps[4] = {vs->x(), vs->y(), vs->z(), 1.}; + double res[4] = {0., 0., 0., 0.}; + int idx = 0; + for(int i = 0; i < 4; i++) + for(int j = 0; j < 4; j++) + res[i] += tfo[idx++] * ps[j]; + + SPoint3 p3 (res[0], res[1], res[2]); + double u = slave->parFromPoint(p3); + GPoint gp = slave->point(u); + MVertex *vt = rtree.find(gp.x(), gp.y(), gp.z()); + if (!vt) Msg::Error("Couldn't find a vertex for updating periodicity"); + else v2v[vt] = vs; + } + } + + for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it) { + GFace *slave = *it; + GFace *master = dynamic_cast<GFace*>(slave->meshMaster()); + if (master == slave) continue; + + std::map<MVertex*,MVertex*> &v2v = slave->correspondingVertices; + v2v.clear(); + + MVertexRTree rtree = MVertexRTree(1e-8 * CTX::instance()->lc); + for (unsigned int i = 0; i < slave->getNumMeshVertices(); ++i) + rtree.insert(slave->getMeshVertex(i)); + + std::vector<double> tfo = slave->affineTransform; + + for (unsigned int i = 0; i < master->getNumMeshVertices(); ++i) { + MVertex *vs = master->getMeshVertex(i); + double ps[4] = {vs->x(), vs->y(), vs->z(), 1.}; + double res[4] = {0., 0., 0., 0.}; + int idx = 0; + for(int i = 0; i < 4; i++) + for(int j = 0; j < 4; j++) + res[i] += tfo[idx++] * ps[j]; + + SPoint3 p3 (res[0], res[1], res[2]); + SPoint2 p2 = slave->parFromPoint(p3); + GPoint gp = slave->point(p2); + MVertex *vt = rtree.find(gp.x(), gp.y(), gp.z()); + if (!vt) Msg::Error("Couldn't find a vertex for updating periodicity"); + else v2v[vt] = vs; + } + } +} + void SetOrder1(GModel *m, bool onlyVisible) { m->destroyMeshCaches(); @@ -1334,6 +1402,8 @@ void SetOrder1(GModel *m, bool onlyVisible) updateHighOrderVertices(*it, newHOVert, onlyVisible); for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it) updateHighOrderVertices(*it, newHOVert, onlyVisible); + + updatePeriodicEdgesAndFaces(m); } void checkHighOrderTriangles(const char* cc, GModel *m, @@ -1500,6 +1570,8 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete, bool onlyVisi for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it) updateHighOrderVertices(*it, newHOVert[*it], onlyVisible); + updatePeriodicEdgesAndFaces(m); + double t2 = Cpu(); std::vector<MElement*> bad; @@ -1559,6 +1631,8 @@ void SetHighOrderComplete(GModel *m, bool onlyVisible) (*it)->mesh_vertices.clear(); (*it)->mesh_vertices.insert((*it)->mesh_vertices.begin(), newV.begin(), newV.end()); } + + updatePeriodicEdgesAndFaces(m); } void SetHighOrderInComplete (GModel *m, bool onlyVisible) @@ -1609,4 +1683,6 @@ void SetHighOrderInComplete (GModel *m, bool onlyVisible) } (*it)->mesh_vertices = newV; } + + updatePeriodicEdgesAndFaces(m); } diff --git a/contrib/HighOrderMeshOptimizer/CMakeLists.txt b/contrib/HighOrderMeshOptimizer/CMakeLists.txt index 559fc502ff3c1a08884307801c1e347a1a6bf293..6f7c6f2a4572e845fc8341341f21a0f8422eaac0 100644 --- a/contrib/HighOrderMeshOptimizer/CMakeLists.txt +++ b/contrib/HighOrderMeshOptimizer/CMakeLists.txt @@ -13,6 +13,7 @@ set(SRC SuperEl.cpp OptHomElastic.cpp OptHomFastCurving.cpp + OptHomPeriodicity.cpp CADDistances.cpp ) diff --git a/contrib/HighOrderMeshOptimizer/OptHomPeriodicity.cpp b/contrib/HighOrderMeshOptimizer/OptHomPeriodicity.cpp new file mode 100644 index 0000000000000000000000000000000000000000..e6cb5ab968b2354cc72615f5760f8d9956f94d00 --- /dev/null +++ b/contrib/HighOrderMeshOptimizer/OptHomPeriodicity.cpp @@ -0,0 +1,186 @@ +// Copyright (C) 2013 ULg-UCL +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, copy, +// modify, merge, publish, distribute, and/or sell copies of the +// Software, and to permit persons to whom the Software is furnished +// to do so, provided that the above copyright notice(s) and this +// permission notice appear in all copies of the Software and that +// both the above copyright notice(s) and this permission notice +// appear in supporting documentation. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF +// MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT OF THIRD PARTY RIGHTS. IN NO EVENT SHALL THE +// COPYRIGHT HOLDER OR HOLDERS INCLUDED IN THIS NOTICE BE LIABLE FOR +// ANY CLAIM, OR ANY SPECIAL INDIRECT OR CONSEQUENTIAL DAMAGES, OR ANY +// DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, +// WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS +// ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE +// OF THIS SOFTWARE. +// +// Please report all bugs and problems to the public mailing list +// <gmsh@onelab.info>. +// +// Contributors: Amaury Johnen + +#include "OptHomPeriodicity.h" +#include "GEntity.h" +#include "GFace.h" +#include "GEdge.h" +#include "MVertex.h" +#include "fullMatrix.h" + +OptHomPeriodicity::OptHomPeriodicity(std::vector<GEntity*> &entities) +{ + for (unsigned int i = 0; i < entities.size(); ++i) { + // MVertex on GVertex cannot move + if (entities[i]->dim() == 0) continue; + + GEntity *master = entities[i]->meshMaster(); + if (master != entities[i]) { + _master2slave.insert(std::make_pair(master, entities[i])); + } + } +} + +void OptHomPeriodicity::fixPeriodicity() +{ + _relocateMasterVertices(); + _copyBackMasterVertices(); +} + +void OptHomPeriodicity::_relocateMasterVertices() +{ + std::multimap<GEntity*, GEntity*>::iterator it; + for (it = _master2slave.begin(); it != _master2slave.end(); ++it) { + if (it->first->dim() == 2) { + GFace *master = dynamic_cast<GFace*>(it->first); + GFace *slave = dynamic_cast<GFace*>(it->second); + std::vector<double> tfo = _inverse(slave->affineTransform); + + std::map<MVertex*, MVertex*> &vertS2M = slave->correspondingVertices; + std::map<MVertex*, MVertex*>::iterator vit; + for (vit = vertS2M.begin(); vit != vertS2M.end(); ++vit) { + GPoint p = _transform(vit->first, master, tfo); + MFaceVertex *v = dynamic_cast<MFaceVertex*>(vit->second); + SPoint3 p3 ((v->x()+p.x())/2, (v->y()+p.y())/2, (v->z()+p.z())/2); + SPoint2 p2 = master->parFromPoint(p3); + GPoint gp = master->point(p2); + v->setXYZ(gp.x(), gp.y(), gp.z()); + v->setParameter(0, gp.u()); + v->setParameter(1, gp.v()); + } + } + else { + GEdge *master = dynamic_cast<GEdge*>(it->first); + int numSlave = _master2slave.count(master); + + for (int i = 0; i < numSlave; ++i) { + if (i > 0) ++it; + GEntity *slave = it->second; + std::vector<double> tfo = _inverse(slave->affineTransform); + + std::map<MVertex*, MVertex*> &vertS2M = slave->correspondingVertices; + std::map<MVertex*, MVertex*>::iterator vit; + for (vit = vertS2M.begin(); vit != vertS2M.end(); ++vit) { + GPoint p = _transform(vit->first, master, tfo); + MVertex *v = vit->second; + v->setXYZ(v->x()+p.x(), v->y()+p.y(), v->z()+p.z()); + } + } + + double coeff = 1./(1+numSlave); + for (unsigned int k = 0; k < master->getNumMeshVertices(); ++k) { + MEdgeVertex *v = dynamic_cast<MEdgeVertex*>(master->getMeshVertex(k)); + SPoint3 p3 (v->x()*coeff, v->y()*coeff, v->z()*coeff); + double u = master->parFromPoint(p3); + GPoint gp = master->point(u); + v->setXYZ(gp.x(), gp.y(), gp.z()); + v->setParameter(0, gp.u()); + } + } + } +} + +void OptHomPeriodicity::_copyBackMasterVertices() +{ + std::multimap<GEntity*, GEntity*>::iterator it; + for (it = _master2slave.begin(); it != _master2slave.end(); ++it) { + if (it->first->dim() == 2) { + GFace *master = dynamic_cast<GFace*>(it->first); + GFace *slave = dynamic_cast<GFace*>(it->second); + std::vector<double> tfo = slave->affineTransform; + + std::map<MVertex*, MVertex*> &vertS2M = slave->correspondingVertices; + std::map<MVertex*, MVertex*>::iterator vit; + for (vit = vertS2M.begin(); vit != vertS2M.end(); ++vit) { + GPoint p = _transform(vit->second, slave, tfo); + MFaceVertex *v = dynamic_cast<MFaceVertex*>(vit->first); + v->setXYZ(p.x(), p.y(), p.z()); + v->setParameter(0, p.u()); + v->setParameter(1, p.v()); + } + } + else { + GEdge *master = dynamic_cast<GEdge*>(it->first); + GEdge *slave = dynamic_cast<GEdge*>(it->second); + std::vector<double> tfo = slave->affineTransform; + + std::map<MVertex*, MVertex*> &vertS2M = slave->correspondingVertices; + std::map<MVertex*, MVertex*>::iterator vit; + for (vit = vertS2M.begin(); vit != vertS2M.end(); ++vit) { + GPoint p = _transform(vit->second, slave, tfo); + MEdgeVertex *v = dynamic_cast<MEdgeVertex*>(vit->first); + v->setXYZ(p.x(), p.y(), p.z()); + v->setParameter(0, p.u()); + } + } + } +} + +GPoint OptHomPeriodicity::_transform(MVertex *vsource, GEntity *target, std::vector<double> &tfo) +{ + double ps[4] = {vsource->x(), vsource->y(), vsource->z(), 1.}; + double res[4] = {0., 0., 0., 0.}; + int idx = 0; + for(int i = 0; i < 4; i++) + for(int j = 0; j < 4; j++) + res[i] += tfo[idx++] * ps[j]; + + SPoint3 p3 (res[0], res[1], res[2]); + if (target->dim() == 2) { + SPoint2 p2 = dynamic_cast<GFace*>(target)->parFromPoint(p3); + return dynamic_cast<GFace*>(target)->point(p2); + } + else if (target->dim() == 1) { + double u = dynamic_cast<GEdge*>(target)->parFromPoint(p3); + return dynamic_cast<GEdge*>(target)->point(u); + } + else { + Msg::Error("Expected a face or an edge for computing " + "periodicity transformation."); + return GPoint(); + } +} + +std::vector<double> OptHomPeriodicity::_inverse(std::vector<double> &tfo) +{ + // Note that the last row of tfo must be (0 0 0 1)... + std::vector<double> result(16); + fullMatrix<double> mat(4, 4), inv; + int idx = 0; + for(int i = 0; i < 4; i++) + for(int j = 0; j < 4; j++) + mat(i, j) = tfo[idx++]; + mat.invert(inv); + + idx = 0; + for(int i = 0; i < 4; i++) + for(int j = 0; j < 4; j++) + result[idx++] = inv(i, j); + return result; +} diff --git a/contrib/HighOrderMeshOptimizer/OptHomPeriodicity.h b/contrib/HighOrderMeshOptimizer/OptHomPeriodicity.h new file mode 100644 index 0000000000000000000000000000000000000000..919c40869e34f83fd587e5459c5ccc94d8e03403 --- /dev/null +++ b/contrib/HighOrderMeshOptimizer/OptHomPeriodicity.h @@ -0,0 +1,59 @@ +// Copyright (C) 2013 ULg-UCL +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, copy, +// modify, merge, publish, distribute, and/or sell copies of the +// Software, and to permit persons to whom the Software is furnished +// to do so, provided that the above copyright notice(s) and this +// permission notice appear in all copies of the Software and that +// both the above copyright notice(s) and this permission notice +// appear in supporting documentation. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF +// MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT OF THIRD PARTY RIGHTS. IN NO EVENT SHALL THE +// COPYRIGHT HOLDER OR HOLDERS INCLUDED IN THIS NOTICE BE LIABLE FOR +// ANY CLAIM, OR ANY SPECIAL INDIRECT OR CONSEQUENTIAL DAMAGES, OR ANY +// DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, +// WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS +// ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE +// OF THIS SOFTWARE. +// +// Please report all bugs and problems to the public mailing list +// <gmsh@onelab.info>. +// +// Contributors: Amaury Johnen + +#ifndef _OPTHOMPERIODICITY_H_ +#define _OPTHOMPERIODICITY_H_ + +#include <map> +#include <vector> + +class GEntity; +class GFace; +class MVertex; +class GPoint; + +class OptHomPeriodicity { +private: + std::multimap<GEntity*, GEntity*> _master2slave; + //std::map<GEntity*, std::map<MVertex*, MVertex*> > _ent2vv; + +public: + OptHomPeriodicity(std::vector<GEntity*>&); + + void fixPeriodicity(); + +private: + void _relocateMasterVertices(); + void _copyBackMasterVertices(); + + static GPoint _transform(MVertex*, GEntity*, std::vector<double>&); + static std::vector<double> _inverse(std::vector<double>&); +}; + +#endif diff --git a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp index dfe512dae405ad9df3194a71a4895a8b74fa3609..7999af6218105707840cec91d7e5daf714a35ab2 100644 --- a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp +++ b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp @@ -34,6 +34,7 @@ #include "GmshConfig.h" #include "OptHOM.h" #include "OptHomRun.h" +#include "OptHomPeriodicity.h" #include "GModel.h" #include "Gmsh.h" #include "MTriangle.h" @@ -118,12 +119,12 @@ void exportMeshToDassault(GModel *gm, const std::string &fn, int dim) for (GModel::fiter itf = gm->firstFace(); itf != gm->lastFace(); ++itf){ std::vector<MTriangle*> &tris = (*itf)->triangles; for (size_t i=0;i<tris.size();i++){ - MTriangle *t = tris[i]; - fprintf(f,"%d ", count++); - for (int j=0;j<t->getNumVertices();j++){ - fprintf(f,"%d ", t->getVertex(j)->getIndex()); - } - fprintf(f,"\n"); + MTriangle *t = tris[i]; + fprintf(f,"%d ", count++); + for (int j=0;j<t->getNumVertices();j++){ + fprintf(f,"%d ", t->getVertex(j)->getIndex()); + } + fprintf(f,"\n"); } } int ne = 0; @@ -136,12 +137,12 @@ void exportMeshToDassault(GModel *gm, const std::string &fn, int dim) for (GModel::eiter ite = gm->firstEdge(); ite != gm->lastEdge(); ++ite){ std::vector<MLine*> &l = (*ite)->lines; for (size_t i=0;i<l.size();i++){ - MLine *t = l[i]; - fprintf(f,"%d ", count++); - for (int j=0;j<t->getNumVertices();j++){ - fprintf(f,"%d ", t->getVertex(j)->getIndex()); - } - fprintf(f,"%d \n",(*ite)->tag()); + MLine *t = l[i]; + fprintf(f,"%d ", count++); + for (int j=0;j<t->getNumVertices();j++){ + fprintf(f,"%d ", t->getVertex(j)->getIndex()); + } + fprintf(f,"%d \n",(*ite)->tag()); } } } @@ -437,6 +438,7 @@ static void optimizeConnectedBlobs(const vertElVecMap &vertex2elements, std::vector<elSetVertSetPair> toOptimize = getConnectedBlobs(vertex2elements, badasses, p.nbLayers, p.distanceFactor, weakMerge, p.optPrimSurfMesh); + p.numBlobs = static_cast<int>(toOptimize.size()); //#pragma omp parallel for schedule(dynamic, 1) for (int i = 0; i < toOptimize.size(); ++i) { @@ -507,6 +509,7 @@ static void optimizeOneByOne const int initNumBadElts = badElts.size(); Msg::Info("%d badasses, starting to iterate...", initNumBadElts); + p.numBlobs = initNumBadElts; elElSetMap element2elements; // Element to element connectivity, built progressively @@ -603,10 +606,10 @@ double ComputeDistanceToGeometry (GEntity *ge , int distanceDefinition, double t if (ge->dim() == el->getDim()){ const double DISTE =computeBndDist(el,distanceDefinition, tolerance); if (DISTE != 0.0){ - NUM++; - // if(distanceDefinition == 1)printf("%d %12.5E\n",iEl,DISTE); - maxd = std::max(maxd,DISTE); - sum += DISTE; + NUM++; + // if(distanceDefinition == 1)printf("%d %12.5E\n",iEl,DISTE); + maxd = std::max(maxd,DISTE); + sum += DISTE; } } } @@ -630,7 +633,6 @@ void HighOrderMeshOptimizer(GModel *gm, OptHomParameters &p) std::map<MVertex*, std::vector<MElement *> > vertex2elements; std::map<MElement*,GEntity*> element2entity; - elSet badasses; double maxdist = 0.; // TODO: To be cleaned? std::map<MElement*,double> distances; @@ -643,33 +645,37 @@ void HighOrderMeshOptimizer(GModel *gm, OptHomParameters &p) entity->tag()); calcVertex2Elements(p.dim,entity,vertex2elements); if (p.optPrimSurfMesh) calcElement2Entity(entity,element2entity); - for (int iEl = 0; iEl < entity->getNumMeshElements();iEl++) { // Detect bad elements - double jmin, jmax; - MElement *el = entity->getMeshElement(iEl); - if (el->getDim() == p.dim) { - // FIXME TEST - // badasses.insert(el); + } + + OptHomPeriodicity periodicity = OptHomPeriodicity(entities); + do { + // Detect bad elements + elSet badasses; + for (int iEnt = 0; iEnt < entities.size(); ++iEnt) { + GEntity* &entity = entities[iEnt]; + if (entity->dim() != p.dim || (p.onlyVisible && !entity->getVisibility())) continue; + for (int iEl = 0; iEl < entity->getNumMeshElements();iEl++) { + MElement *el = entity->getMeshElement(iEl); if (p.optCAD) { - // const double DISTE =computeBndDist(el,2,fabs(p.discrTolerance)); - const double DISTE =distances[el]; - // if (DISTE > 0)printf("El %d dist %12.5E vs %12.5E\n",iEl,DISTE,p.optCADDistMax); + const double DISTE =distances[el]; maxdist = std::max(DISTE, maxdist); if (DISTE > p.optCADDistMax) badasses.insert(el); } + double jmin, jmax; el->scaledJacRange(jmin, jmax, p.optPrimSurfMesh ? entity : 0); if (p.BARRIER_MIN_METRIC > 0) jmax = jmin; if (jmin < p.BARRIER_MIN || jmax > p.BARRIER_MAX) badasses.insert(el); } } - } - printf("maxdist = %g badasses size = %lu\n", maxdist, badasses.size()); - - if (p.strategy == 0) - optimizeConnectedBlobs(vertex2elements, element2entity, badasses, p, samples, false); - else if (p.strategy == 2) - optimizeConnectedBlobs(vertex2elements, element2entity, badasses, p, samples, true); - else - optimizeOneByOne(vertex2elements, element2entity, badasses, p, samples); + printf("maxdist = %g badasses size = %lu\n", maxdist, badasses.size()); + if (p.strategy == 0) + optimizeConnectedBlobs(vertex2elements, element2entity, badasses, p, samples, false); + else if (p.strategy == 2) + optimizeConnectedBlobs(vertex2elements, element2entity, badasses, p, samples, true); + else + optimizeOneByOne(vertex2elements, element2entity, badasses, p, samples); + if (p.numBlobs) periodicity.fixPeriodicity(); + } while (p.numBlobs); if (p.SUCCESS == 1) Msg::Info("Optimization succeeded"); diff --git a/contrib/HighOrderMeshOptimizer/OptHomRun.h b/contrib/HighOrderMeshOptimizer/OptHomRun.h index 88b162e25051920439d8f42892b1e69e4f59f4c4..8c1b2b61a120a18e1cfe114f3ee417b94a4f98b1 100644 --- a/contrib/HighOrderMeshOptimizer/OptHomRun.h +++ b/contrib/HighOrderMeshOptimizer/OptHomRun.h @@ -59,7 +59,8 @@ struct OptHomParameters { double discrTolerance; // OUTPUT ------> - int SUCCESS ; // 0 --> success , 1 --> Not converged + int SUCCESS; // 0 --> success , 1 --> Not converged + int numBlobs; double minJac, maxJac; // after optimization, range of jacobians double CPU; // Time for optimization