From a4e1795797a0144745387078a72be497876aa1fe Mon Sep 17 00:00:00 2001
From: Thomas Toulorge <thomas.toulorge@mines-paristech.fr>
Date: Tue, 7 Mar 2017 23:50:17 +0000
Subject: [PATCH] Updated fast high-order BL mesh curving algorithm
---
contrib/HighOrderMeshOptimizer/CMakeLists.txt | 2 +-
contrib/HighOrderMeshOptimizer/MetaEl.cpp | 290 ++++++++++++++++++
.../{SuperEl.h => MetaEl.h} | 38 +--
.../OptHomFastCurving.cpp | 162 +++++-----
contrib/HighOrderMeshOptimizer/SuperEl.cpp | 233 --------------
5 files changed, 388 insertions(+), 337 deletions(-)
create mode 100644 contrib/HighOrderMeshOptimizer/MetaEl.cpp
rename contrib/HighOrderMeshOptimizer/{SuperEl.h => MetaEl.h} (70%)
delete mode 100644 contrib/HighOrderMeshOptimizer/SuperEl.cpp
diff --git a/contrib/HighOrderMeshOptimizer/CMakeLists.txt b/contrib/HighOrderMeshOptimizer/CMakeLists.txt
index 79c7deb4be..4f4e6d6dd6 100644
--- a/contrib/HighOrderMeshOptimizer/CMakeLists.txt
+++ b/contrib/HighOrderMeshOptimizer/CMakeLists.txt
@@ -10,7 +10,7 @@ set(SRC
OptHomIntegralBoundaryDist.cpp
OptHomCADDist.cpp
ParamCoord.cpp
- SuperEl.cpp
+ MetaEl.cpp
OptHomElastic.cpp
OptHomFastCurving.cpp
OptHomPeriodicity.cpp
diff --git a/contrib/HighOrderMeshOptimizer/MetaEl.cpp b/contrib/HighOrderMeshOptimizer/MetaEl.cpp
new file mode 100644
index 0000000000..8a3c3ea7a3
--- /dev/null
+++ b/contrib/HighOrderMeshOptimizer/MetaEl.cpp
@@ -0,0 +1,290 @@
+// 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>.
+//
+// Contributor: Thomas Toulorge
+
+#include <iostream>
+#include <sstream>
+#include "GmshMessage.h"
+#include "GEdge.h"
+#include "MLine.h"
+#include "MQuadrangle.h"
+#include "MPrism.h"
+#include "MHexahedron.h"
+#include "BasisFactory.h"
+#include "MetaEl.h"
+
+
+std::map<int, MetaEl::metaInfoType> MetaEl::_metaInfo;
+
+
+//namespace {
+//
+//
+//void getEdgeIndQuad(int order, std::vector<int> edInd) {
+// edInd.clear();
+// const int nbEdVert = 4*(order-1);
+// edInd.resize(nbEdVert);
+// for (int iV = 0; iV < nbEdVert; iV++) edInd[iV] = iV;
+//}
+//
+//
+//
+//}
+
+
+MetaEl::metaInfoType::metaInfoType(int type, int order)
+{
+ typedef std::vector<int>::iterator VIntIter;
+
+ int iBaseFace = 0, iTopFace = 0, nbEdVert = 0;
+ switch (type) {
+ case TYPE_QUA:
+ iBaseFace = 0; iTopFace = 2; nbEdVert = 4 + 4*(order-1);
+ break;
+ case TYPE_PRI:
+ iBaseFace = 0; iTopFace = 1; nbEdVert = 6 + 9*(order-1);
+ break;
+ case TYPE_HEX:
+ iBaseFace = 0; iTopFace = 5; nbEdVert = 8 + 12*(order-1);
+ break;
+ default:
+ Msg::Error("MetaEl not implemented for element of type %d", type);
+ nbVert = 0;
+ return;
+ }
+
+ // Get HO nodal basis
+// const int tag = ElementType::getTag(type, order, true); // Get tag for incomplete element type
+ const int tag = ElementType::getTag(type, order); // Get tag for complete element type
+ if (!tag) return;
+ const nodalBasis *basis = BasisFactory::getNodalBasis(tag);
+ nbVert = basis->getNumShapeFunctions();
+ points = basis->points;
+ std::set<int> foundVerts;
+
+ // Vertices on base and top faces
+ baseInd = basis->getClosure(basis->getClosureId(iBaseFace, 1, 0));
+ topInd = basis->getClosure(basis->getClosureId(iTopFace, 0, 0));
+ foundVerts.insert(baseInd.begin(), baseInd.end());
+ foundVerts.insert(topInd.begin(), topInd.end());
+
+ // Vertices on edges (excluding base and top faces)
+ for (int iV = 0; iV < nbEdVert; iV++)
+ if (foundVerts.find(iV) == foundVerts.end()) {
+ edgeInd.push_back(iV);
+ foundVerts.insert(iV);
+ }
+
+ // Remaining face and interior vertices
+ for(int iV = 0; iV < nbVert; iV++)
+ if (foundVerts.find(iV) == foundVerts.end()) {
+ otherInd.push_back(iV);
+ foundVerts.insert(iV);
+ }
+}
+
+
+MetaEl::MetaEl(int type, int order, const std::vector<MVertex*> &baseVert,
+ const std::vector<MVertex*> &topPrimVert) : _metaEl(0), _metaEl0(0)
+{
+
+ // Get useful info on meta-element type if not already there
+ std::map<int, MetaEl::metaInfoType>::iterator itSInfo = _metaInfo.find(type);
+ if (itSInfo == _metaInfo.end()) {
+ const metaInfoType mInfo(type, order);
+ itSInfo = _metaInfo.insert(std::pair<int,metaInfoType>(type, mInfo)).first;
+ }
+ MetaEl::metaInfoType &sInfo = itSInfo->second;
+
+ // Exit if unknown type
+ if (sInfo.nbVert == 0) return;
+
+ // References for easier writing
+ const int &nbVert = sInfo.nbVert;
+ const fullMatrix<double> &points = sInfo.points;
+ const std::vector<int> &baseInd = sInfo.baseInd;
+ const std::vector<int> &topInd = sInfo.topInd;
+ const std::vector<int> &edgeInd = sInfo.edgeInd;
+ const std::vector<int> &otherInd = sInfo.otherInd;
+
+ // Add copies of vertices in base & top faces (only first-order vertices for top face)
+ _metaVert.resize(nbVert);
+ for (int iV = 0; iV < baseInd.size(); iV++)
+ _metaVert[baseInd[iV]] = new MVertex(*baseVert[iV]);
+ for (int iV = 0; iV < topPrimVert.size(); iV++)
+ _metaVert[topInd[iV]] = new MVertex(*topPrimVert[iV]);
+
+ // Create first-order meta-element
+ switch (type) {
+ case TYPE_QUA:
+ _metaEl0 = new MQuadrangle(_metaVert);
+ break;
+ case TYPE_PRI:
+ _metaEl0 = new MPrism(_metaVert);
+ break;
+ case TYPE_HEX:
+ _metaEl0 = new MHexahedron(_metaVert);
+ break;
+ default:
+ Msg::Error("SuperEl not implemented for element of type %d", type);
+ return;
+ }
+
+ // Add HO vertices in top face
+ for (int iV = topPrimVert.size(); iV < topInd.size(); iV++) {
+ SPoint3 p;
+ const int ind = topInd[iV];
+ _metaEl0->pnt(points(ind, 0), points(ind, 1), points(ind, 2), p);
+ _metaVert[ind] = new MVertex(p.x(), p.y(), p.z());
+ }
+
+ // Add vertices on edges (excluding base and top faces)
+ for (int iV = 0; iV < edgeInd.size(); iV++) {
+ SPoint3 p;
+ const int ind = edgeInd[iV];
+ _metaEl0->pnt(points(ind, 0), points(ind, 1), points(ind, 2), p);
+ _metaVert[ind] = new MVertex(p.x(), p.y(), p.z());
+ }
+
+ // Get incomplete high-order meta-element basis
+ const int tag = ElementType::getTag(type, order, true); // Get tag for incomplete element type
+ if (!tag) return;
+ const nodalBasis *incBasis = BasisFactory::getNodalBasis(tag);
+
+ // Add remaining face and interior points
+ for (int iV = 0; iV < otherInd.size(); iV++) {
+ const int ind = otherInd[iV];
+ double shapeFunc[1256];
+ incBasis->f(points(ind, 0), points(ind, 1), points(ind, 2), shapeFunc);
+ double x = 0., y = 0., z = 0.;
+ for (int iSF = 0; iSF < incBasis->getNumShapeFunctions(); iSF++) {
+ x += shapeFunc[iSF] * _metaVert[iSF]->x();
+ y += shapeFunc[iSF] * _metaVert[iSF]->y();
+ z += shapeFunc[iSF] * _metaVert[iSF]->z();
+ }
+ _metaVert[ind] = new MVertex(x, y, z);
+ }
+
+ // Create high-order meta-element
+ switch (type) {
+ case TYPE_QUA:
+ _metaEl = new MQuadrangleN(_metaVert, order);
+ break;
+ case TYPE_PRI:
+ _metaEl = new MPrismN(_metaVert, order);
+ break;
+ case TYPE_HEX:
+ _metaEl = new MHexahedronN(_metaVert, order);
+ break;
+ }
+
+// // Scale meta-element if not valid
+// // TODO: Scale depending on target Jmin?
+// // TODO: Could be improved by using complete meta-element and use optimization
+// for (int iter = 0; iter < 10; iter++) {
+// if (_metaEl->distoShapeMeasure() > 0) break;
+// if (iter == 0)
+// Msg::Warning("A meta-element has a negative distortion, trying to scale");
+// for (int i = 0; i < topPrimVert.size(); ++i) { // Move top primary vert.
+// MVertex *&vb = _metaVert[baseInd[i]], *&vt = _metaVert[topInd[i]];
+// SPoint3 pb = vb->point(), pt = vt->point();
+// pt += SPoint3(0.25*(pt-pb));
+// vt->setXYZ(pt.x(), pt.y(), pt.z());
+// }
+// for (int i=topPrimVert.size(); i<topInd.size(); ++i) { // Recompute HO vert. in top face
+// SPoint3 p;
+// const int ind = topInd[i];
+// _metaEl0->pnt(points(ind,0), points(ind,1), points(ind,2),p);
+// _metaVert[ind]->setXYZ(p.x(), p.y(), p.z());
+// }
+// for (int i=0; i<otherInd.size(); ++i) { // Recompute vert. not in base & top faces
+// SPoint3 p;
+// const int ind = otherInd[i];
+// _metaEl0->pnt(points(ind,0), points(ind,1), points(ind,2),p);
+// _metaVert[ind]->setXYZ(p.x(), p.y(), p.z());
+// }
+// }
+
+}
+
+
+MetaEl::~MetaEl()
+{
+ for (int i = 0; i < _metaVert.size(); i++) delete _metaVert[i];
+ _metaVert.clear();
+ if (_metaEl) delete _metaEl;
+ if (_metaEl0) delete _metaEl0;
+}
+
+
+bool MetaEl::isPointIn(const SPoint3 p) const
+{
+ double xyz[3] = {p.x(), p.y(), p.z()}, uvw[3];
+ _metaEl0->xyz2uvw(xyz, uvw);
+ return _metaEl0->isInside(uvw[0], uvw[1], uvw[2]);
+}
+
+
+bool MetaEl::straightToCurved(double *xyzS, double *xyzC) const
+{
+ double uvw[3];
+ _metaEl0->xyz2uvw(xyzS, uvw);
+ if (!_metaEl0->isInside(uvw[0], uvw[1], uvw[2])) return false;
+
+ SPoint3 pC;
+ _metaEl->pnt(uvw[0], uvw[1], uvw[2], pC);
+ xyzC[0] = pC[0];
+ xyzC[1] = pC[1];
+ xyzC[2] = pC[2];
+
+ return true;
+}
+
+
+std::string MetaEl::printPOS()
+{
+ std::vector<MVertex*> verts;
+ _metaEl->getVertices(verts);
+ std::string posStr = _metaEl0->getStringForPOS();
+ int n = (posStr[posStr.size()-1] == '2' ?
+ _metaEl : _metaEl0)->getNumVertices();
+
+ std::ostringstream oss;
+
+ oss << posStr << "(";
+ for(int i = 0; i < n; i++){
+ if(i) oss << ",";
+ oss << _metaVert[i]->x() << "," << _metaVert[i]->y() << ","
+ << _metaVert[i]->z();
+ }
+ oss << "){0.";
+ for(int i = 1; i < n; i++) oss << ",0.";
+ oss << "};\n";
+
+ return oss.str();
+}
diff --git a/contrib/HighOrderMeshOptimizer/SuperEl.h b/contrib/HighOrderMeshOptimizer/MetaEl.h
similarity index 70%
rename from contrib/HighOrderMeshOptimizer/SuperEl.h
rename to contrib/HighOrderMeshOptimizer/MetaEl.h
index 180c013b12..37ad9bbb12 100644
--- a/contrib/HighOrderMeshOptimizer/SuperEl.h
+++ b/contrib/HighOrderMeshOptimizer/MetaEl.h
@@ -27,41 +27,43 @@
//
// Contributor: Thomas Toulorge
-#ifndef _SUPEREL_H_
-#define _SUPEREL_H_
+#ifndef _METAEL_H_
+#define _METAEL_H_
#include <string>
#include "MElement.h"
-class SuperEl
+class MetaEl
{
public:
- SuperEl(int type, int order, const std::vector<MVertex*> &baseVert,
+ MetaEl(int type, int order, const std::vector<MVertex*> &baseVert,
const std::vector<MVertex*> &topPrimVert);
- ~SuperEl();
- bool isOK() const { return _superEl; }
+ ~MetaEl();
+ bool isOK() const { return _metaEl; }
bool isPointIn(const SPoint3 p) const;
bool straightToCurved(double *xyzS, double *xyzC) const;
std::string printPOS();
void printCoord()
{
std::cout << "DBGTT: superEl -> ";
- for(int i = 0; i < _superVert.size(); i++){
- std::cout << "v" << i << " = (" << _superVert[i]->x() << ","
- << _superVert[i]->y() << "," << _superVert[i]->z() << ")";
- if (i == _superVert.size()-1) std::cout << "\n"; else std::cout << ", ";
+ for(int i = 0; i < _metaVert.size(); i++){
+ std::cout << "v" << i << " = (" << _metaVert[i]->x() << ","
+ << _metaVert[i]->y() << "," << _metaVert[i]->z() << ")";
+ if (i == _metaVert.size()-1) std::cout << "\n"; else std::cout << ", ";
}
}
+ MElement *getMElement() { return _metaEl; }
+
private:
- struct superInfoType {
- int nV;
+ struct metaInfoType {
+ int nbVert;
fullMatrix<double> points;
- std::vector<int> baseInd, topInd, otherInd;
- superInfoType(int type, int order);
+ std::vector<int> baseInd, topInd, edgeInd, otherInd;
+ metaInfoType(int type, int order);
};
- static std::map<int,superInfoType> _superInfo;
- std::vector<MVertex*> _superVert;
- MElement *_superEl, *_superEl0;
+ static std::map<int, metaInfoType> _metaInfo;
+ std::vector<MVertex*> _metaVert;
+ MElement *_metaEl, *_metaEl0;
};
-#endif
+#endif // _METAEL_H_
diff --git a/contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp b/contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp
index f71e14d542..0c3ab4cea3 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp
@@ -44,10 +44,10 @@
#include "MLine.h"
#include "OS.h"
#include <stack>
-#include "SuperEl.h"
#include "SVector3.h"
#include "BasisFactory.h"
#include "Field.h"
+#include "MetaEl.h"
@@ -300,31 +300,23 @@ bool getBaseVertices(const MFace &baseFace, MElement *firstEl, std::vector<MVert
template<class bndType>
bool getTopPrimVertices(std::map<MVertex*, std::vector<MElement *> > &vertex2elements,
const std::map<MVertex*, SVector3> &normVert,
- const BoundaryLayerColumns *blc, const bndType &baseBnd,
- int maxNumLayers, const std::vector<MVertex*> &baseVert,
+ const bndType &baseBnd, int maxNumLayers,
+ const std::vector<MVertex*> &baseVert,
std::vector<MVertex*> &topPrimVert)
{
int nbVert = baseBnd.getNumVertices();
topPrimVert = std::vector<MVertex*>(nbVert);
for(int i = 0; i < nbVert; i++) {
- if (blc && (blc->size() > 1)) { // Is there BL data?
- return false;
- // const BoundaryLayerData &c = blc->getColumn(baseVert[i], baseBnd);
- // if (c._column.size() == 0) return false; // Give up element if column not found
- // topPrimVert[i] = *(c._column.end()-2);
- }
- else { // No BL data, look for columns of vertices
- std::map<MVertex*, SVector3>::const_iterator itNorm = normVert.find(baseVert[i]);
- if (itNorm == normVert.end()) {
- Msg::Error("Normal to vertex not found in getTopPrimVertices");
- itNorm = normVert.begin();
- }
- std::vector<MVertex*> col;
- bool colFound = findColumn(vertex2elements, baseVert[i], itNorm->second,
- baseBnd, col, maxNumLayers);
- if (!colFound) return false; // Give up element if column not found
- topPrimVert[i] = *(col.end()-2);
+ std::map<MVertex*, SVector3>::const_iterator itNorm = normVert.find(baseVert[i]);
+ if (itNorm == normVert.end()) {
+ Msg::Error("Normal to vertex not found in getTopPrimVertices");
+ itNorm = normVert.begin();
}
+ std::vector<MVertex*> col;
+ bool colFound = findColumn(vertex2elements, baseVert[i], itNorm->second,
+ baseBnd, col, maxNumLayers);
+ if (!colFound) return false; // Give up element if column not found
+ topPrimVert[i] = *(col.end()-2);
}
return true;
@@ -335,7 +327,7 @@ bool getTopPrimVertices(std::map<MVertex*, std::vector<MElement *> > &vertex2ele
// Get blob of elements encompassed by a given meta-element
// FIXME: Implement 3D
void get2DBLBlob(const std::map<MVertex*, std::vector<MElement*> > &vertex2elements,
- MElement *firstEl, const SuperEl *supEl, std::set<MElement*> &blob)
+ MElement *firstEl, const MetaEl *supEl, std::set<MElement*> &blob)
{
static const int maxDepth = 100;
typedef std::list<MElement*> elLType;
@@ -420,10 +412,52 @@ inline void insertIfCurved(MElement *el, std::list<MElement*> &bndEl)
}
+class DbgOutput {
+public:
+ void addMetaEl(MetaEl &mEl) {
+ MElement *elt = mEl.getMElement();
+ elType_.push_back(elt->getTypeForMSH());
+ nbVertEl_.push_back(elt->getNumVertices());
+ for (int iV = 0; iV < elt->getNumVertices(); iV++)
+ point_.push_back(elt->getVertex(iV)->point());
+ }
+ void write(std::string fNameBase, int tag) {
+ std::ostringstream oss;
+ oss << fNameBase << "_" << tag << ".msh";
+ std::string fName = oss.str();
+ FILE *fp = fopen(fName.c_str(), "w");
+ fprintf(fp, "$MeshFormat\n2.2 0 8\n$EndMeshFormat\n");
+ fprintf(fp, "$Nodes\n");
+ fprintf(fp, "%d\n", point_.size());
+ for (int iV = 0; iV < point_.size(); iV++)
+ fprintf(fp, "%i %g %g %g\n", iV+1, point_[iV].x(),
+ point_[iV].y(), point_[iV].z());
+ fprintf(fp, "$EndNodes\n");
+ fprintf(fp, "$Elements\n");
+ fprintf(fp, "%d\n", elType_.size());
+ int iV = 0;
+ for (int iEl = 0; iEl < elType_.size(); iEl++) {
+ fprintf(fp, "%i %i 2 0 0 ", iEl+1, elType_[iEl]);
+ for (int iVEl = 1; iVEl <= nbVertEl_[iEl]; iVEl++)
+ fprintf(fp, " %i", iV+iVEl);
+ fprintf(fp, "\n");
+ iV += nbVertEl_[iEl];
+ }
+ fprintf(fp, "$EndElements\n");
+ fclose(fp);
+ }
+
+private:
+ std::vector<int> elType_;
+ std::vector<int> nbVertEl_;
+ std::vector<SPoint3> point_;
+};
+
+
// Curve elements adjacent to a boundary model entity
void curveMeshFromBnd(std::map<MVertex*, std::vector<MElement *> > &vertex2elements,
- const std::map<MVertex*, SVector3> &normVert, BoundaryLayerColumns *blc,
+ const std::map<MVertex*, SVector3> &normVert,
GEntity *bndEnt, FastCurvingParameters &p)
{
// Build list of bnd. elements to consider
@@ -443,28 +477,25 @@ void curveMeshFromBnd(std::map<MVertex*, std::vector<MElement *> > &vertex2eleme
else
Msg::Error("Cannot treat model entity %i of dim %i", bndEnt->tag(), bndEnt->dim());
- std::ostringstream oss;
- oss << "meta-elements_" << bndEnt->tag() << ".pos";
- std::ofstream of(oss.str().c_str());
- of << "View \" \"{\n";
+ DbgOutput dbgOut;
std::set<MVertex*> movedVert;
for(std::list<MElement*>::iterator itBE = bndEl.begin();
itBE != bndEl.end(); itBE++) { // Loop over bnd. elements
const int bndType = (*itBE)->getType();
- int seType;
+ int metaElType;
bool foundVert;
MElement *firstEl = 0;
std::vector<MVertex*> baseVert, topPrimVert;
if (bndType == TYPE_LIN) { // 1D boundary?
MVertex *vb0 = (*itBE)->getVertex(0);
MVertex *vb1 = (*itBE)->getVertex(1);
- seType = TYPE_QUA;
+ metaElType = TYPE_QUA;
MEdge baseEd(vb0, vb1);
firstEl = getFirstEl(vertex2elements, baseEd);
foundVert = getBaseVertices(baseEd, firstEl, baseVert);
if (!foundVert) continue; // Skip bnd. el. if base vertices not found
- foundVert = getTopPrimVertices(vertex2elements, normVert, blc, baseEd,
+ foundVert = getTopPrimVertices(vertex2elements, normVert, baseEd,
p.maxNumLayers, baseVert, topPrimVert);
}
else { // 2D boundary
@@ -474,11 +505,11 @@ void curveMeshFromBnd(std::map<MVertex*, std::vector<MElement *> > &vertex2eleme
MVertex *vb3;
if (bndType == TYPE_QUA) {
vb3 = (*itBE)->getVertex(3);
- seType = TYPE_HEX;
+ metaElType = TYPE_HEX;
}
else {
vb3 = 0;
- seType = TYPE_PRI;
+ metaElType = TYPE_PRI;
}
MFace baseFace(vb0, vb1, vb2, vb3);
firstEl = getFirstEl(vertex2elements, baseFace);
@@ -489,23 +520,23 @@ void curveMeshFromBnd(std::map<MVertex*, std::vector<MElement *> > &vertex2eleme
}
foundVert = getBaseVertices(baseFace, firstEl, baseVert);
if (!foundVert) continue; // Skip bnd. el. if base vertices not found
- foundVert = getTopPrimVertices(vertex2elements, normVert, blc,
+ foundVert = getTopPrimVertices(vertex2elements, normVert,
baseFace, p.maxNumLayers,
baseVert, topPrimVert);
}
if (!foundVert) continue; // Skip bnd. el. if top vertices not found
int order = firstEl->getPolynomialOrder();
- SuperEl se(seType, order, baseVert, topPrimVert);
- of << se.printPOS();
+ MetaEl metaEl(metaElType, order, baseVert, topPrimVert);
+ dbgOut.addMetaEl(metaEl);
std::set<MElement*> blob;
- get2DBLBlob(vertex2elements, firstEl, &se, blob); // TODO: Implement for 3D
+ get2DBLBlob(vertex2elements, firstEl, &metaEl, blob); // TODO: Implement for 3D
for (std::set<MElement*>::iterator itE = blob.begin(); itE != blob.end(); ++itE) {
makeStraight(*itE, movedVert);
for (int i = (*itE)->getNumPrimaryVertices(); i < (*itE)->getNumVertices(); ++i) { // Loop over HO vert. of each el. in blob
MVertex* vert = (*itE)->getVertex(i);
if (movedVert.find(vert) == movedVert.end()) { // Try to move vert. not already moved
double xyzS[3] = {vert->x(), vert->y(), vert->z()}, xyzC[3];
- if (se.straightToCurved(xyzS,xyzC)) {
+ if (metaEl.straightToCurved(xyzS,xyzC)) {
vert->setXYZ(xyzC[0], xyzC[1], xyzC[2]);
movedVert.insert(vert);
}
@@ -514,8 +545,7 @@ void curveMeshFromBnd(std::map<MVertex*, std::vector<MElement *> > &vertex2eleme
}
}
- of << "};\n";
- of.close();
+ dbgOut.write("meta-elements", bndEnt->tag());
}
@@ -539,64 +569,26 @@ void HighOrderMeshFastCurving(GModel *gm, FastCurvingParameters &p)
for (int iEnt = 0; iEnt < allEntities.size(); ++iEnt)
calcVertex2Elements(p.dim, allEntities[iEnt], vertex2elements);
- // Get BL field (if any)
- BoundaryLayerField *blf = getBLField(gm);
-
// Build multimap of each geometric entity to its boundaries
std::multimap<GEntity*,GEntity*> entities;
- if (blf) { // BF field?
- for (int iEnt = 0; iEnt < allEntities.size(); ++iEnt) {
- GEntity* &entity = allEntities[iEnt];
- if (entity->dim() == p.dim && (!p.onlyVisible || entity->getVisibility())) // Consider only "domain" entities
- if (p.dim == 2) { // "Domain" face?
- std::list<GEdge*> edges = entity->edges();
- for (std::list<GEdge*>::iterator itEd = edges.begin(); itEd != edges.end(); itEd++) // Loop over model boundary edges
- if (blf->isEdgeBL((*itEd)->tag())) // Already skip model edge if no BL there
- entities.insert(std::pair<GEntity*,GEntity*>(entity, *itEd));
- }
- // else if (p.dim == 3) { // "Domain" region?
- // std::list<GFace*> faces = entity->faces();
- // for (std::list<GFace*>::iterator itF = faces.begin(); itF != faces.end(); itF++) // Loop over model boundary faces
- // if (blf->isFaceBL((*itF)->tag())) // Already skip model face if no BL there
- // entities.insert(std::pair<GEntity*,GEntity*>(entity, *itF));
- // }
- }
- }
- else { // No BL field
- for (int iEnt = 0; iEnt < allEntities.size(); ++iEnt) {
- GEntity *dummy = 0;
- GEntity* &entity = allEntities[iEnt];
- if (entity->dim() == p.dim-1 && (!p.onlyVisible || entity->getVisibility())) // Consider boundary entities
- entities.insert(std::pair<GEntity*,GEntity*>(dummy, entity));
- }
+ for (int iEnt = 0; iEnt < allEntities.size(); ++iEnt) {
+ GEntity *dummy = 0;
+ GEntity* &entity = allEntities[iEnt];
+ if (entity->dim() == p.dim-1 && (!p.onlyVisible || entity->getVisibility())) // Consider boundary entities
+ entities.insert(std::pair<GEntity*,GEntity*>(dummy, entity));
}
- // Build normals if necessary
+ // Build normals to base vertices
std::map<GEntity*, std::map<MVertex*, SVector3> > normVertEnt; // Normal to each vertex for each geom. entity
- if (!blf) {
- Msg::Warning("Boundary layer data not found, trying to detect columns");
- buildNormals(vertex2elements, entities, p, normVertEnt);
- }
+ buildNormals(vertex2elements, entities, p, normVertEnt);
// Loop over geometric entities
for (std::multimap<GEntity*,GEntity*>::iterator itBE = entities.begin();
itBE != entities.end(); itBE++) {
GEntity *domEnt = itBE->first, *bndEnt = itBE->second;
- BoundaryLayerColumns *blc = 0;
- if (blf) {
- Msg::Info("Curving elements for entity %d bounding entity %d...",
- bndEnt->tag(), domEnt->tag());
- if (p.dim == 2)
- blc = domEnt->cast2Face()->getColumns();
- else if (p.dim == 3)
- blc = domEnt->cast2Region()->getColumns();
- else
- Msg::Error("Fast curving implemented only in dim. 2 and 3");
- }
- else
- Msg::Info("Curving elements for boundary entity %d...", bndEnt->tag());
+ Msg::Info("Curving elements for boundary entity %d...", bndEnt->tag());
std::map<MVertex*, SVector3> &normVert = normVertEnt[bndEnt];
- curveMeshFromBnd(vertex2elements, normVert, blc, bndEnt, p);
+ curveMeshFromBnd(vertex2elements, normVert, bndEnt, p);
}
double t2 = Cpu();
diff --git a/contrib/HighOrderMeshOptimizer/SuperEl.cpp b/contrib/HighOrderMeshOptimizer/SuperEl.cpp
deleted file mode 100644
index 845c548d93..0000000000
--- a/contrib/HighOrderMeshOptimizer/SuperEl.cpp
+++ /dev/null
@@ -1,233 +0,0 @@
-// 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>.
-//
-// Contributor: Thomas Toulorge
-
-#include <iostream>
-#include <sstream>
-#include "GmshMessage.h"
-#include "GEdge.h"
-#include "MLine.h"
-#include "MQuadrangle.h"
-#include "MPrism.h"
-#include "MHexahedron.h"
-#include "BasisFactory.h"
-#include "SuperEl.h"
-
-std::map<int,SuperEl::superInfoType> SuperEl::_superInfo;
-
-SuperEl::superInfoType::superInfoType(int type, int order)
-{
- int iBaseFace = 0, iTopFace = 0;
- switch (type) {
- case TYPE_QUA:
- iBaseFace = 0; iTopFace = 2;
- break;
- case TYPE_PRI:
- iBaseFace = 0; iTopFace = 1;
- break;
- case TYPE_HEX:
- iBaseFace = 0; iTopFace = 5;
- break;
- default:
- Msg::Error("SuperEl not implemented for element of type %d", type);
- nV = 0;
- return;
- }
-
- // Get HO nodal basis
- const int tag = ElementType::getTag(type, order, true); // Get tag for incomplete element type
-// const int tag = ElementType::getTag(type, order); // Get tag for complete element type
- if(tag){
- const nodalBasis *basis = BasisFactory::getNodalBasis(tag);
-
- nV = basis->getNumShapeFunctions();
- // _superInfo[type].nV1 = basis->getNumShapeFunctions();
- points = basis->points;
-
- baseInd = basis->getClosure(basis->getClosureId(iBaseFace,1,0));
- topInd = basis->getClosure(basis->getClosureId(iTopFace,0,0));
- otherInd.reserve(nV-baseInd.size()-topInd.size());
- for(int i=0; i<nV; ++i) {
- const std::vector<int>::iterator inBaseFace = find(baseInd.begin(),baseInd.end(),i);
- const std::vector<int>::iterator inTopFace = find(topInd.begin(),topInd.end(),i);
- if (inBaseFace == baseInd.end() && inTopFace == topInd.end()) otherInd.push_back(i);
- }
- }
-}
-
-
-
-SuperEl::SuperEl(int type, int order, const std::vector<MVertex*> &baseVert,
- const std::vector<MVertex*> &topPrimVert) : _superEl(0), _superEl0(0)
-{
-
- // Get useful info on meta-element type if not already there
- std::map<int,SuperEl::superInfoType>::iterator itSInfo = _superInfo.find(type);
- if (itSInfo == _superInfo.end())
- itSInfo = _superInfo.insert(std::pair<int,superInfoType>(type,superInfoType(type, order))).first;
- SuperEl::superInfoType &sInfo = itSInfo->second;
-
- // Exit if unknown type
- if (sInfo.nV == 0) return;
-
- // References for easier writing
- const int &nV = sInfo.nV;
- const fullMatrix<double> &points = sInfo.points;
- const std::vector<int> &baseInd = sInfo.baseInd;
- const std::vector<int> &topInd = sInfo.topInd;
- const std::vector<int> &otherInd = sInfo.otherInd;
-
- // Add copies of vertices in base & top faces (only first-order vertices for top face)
- _superVert.resize(nV);
- for (int i=0; i<baseInd.size(); ++i) _superVert[baseInd[i]] = new MVertex(*baseVert[i]);
- for (int i=0; i<topPrimVert.size(); ++i) _superVert[topInd[i]] = new MVertex(*topPrimVert[i]);
-
- // Create first-order meta-element
- switch (type) {
- case TYPE_QUA:
- _superEl0 = new MQuadrangle(_superVert);
- break;
- case TYPE_PRI:
- _superEl0 = new MPrism(_superVert);
- break;
- case TYPE_HEX:
- _superEl0 = new MHexahedron(_superVert);
- break;
- default:
- Msg::Error("SuperEl not implemented for element of type %d", type);
- return;
- }
-
- // Add HO vertices in top face
- for (int i=topPrimVert.size(); i<topInd.size(); ++i) {
- SPoint3 p;
- const int ind = topInd[i];
- _superEl0->pnt(points(ind,0),points(ind,1),points(ind,2),p);
- _superVert[ind] = new MVertex(p.x(),p.y(),p.z());
- }
-
- // Add vertices not in base & top faces
- for (int i=0; i<otherInd.size(); ++i) {
- SPoint3 p;
- const int ind = otherInd[i];
- _superEl0->pnt(points(ind,0),points(ind,1),points(ind,2),p);
- _superVert[ind] = new MVertex(p.x(),p.y(),p.z());
- }
-
- // Create high-order meta-element
- switch (type) {
- case TYPE_QUA:
- _superEl = new MQuadrangleN(_superVert, order);
- break;
- case TYPE_PRI:
- _superEl = new MPrismN(_superVert, order);
- break;
- case TYPE_HEX:
- _superEl = new MHexahedronN(_superVert, order);
- break;
- }
-
- // Scale meta-element if not valid
- // TODO: Scale depending on target Jmin?
- // TODO: Could be improved by using complete meta-element and use optimization
- for (int iter = 0; iter < 10; iter++) {
- if (_superEl->distoShapeMeasure() > 0) break;
- if (iter == 0) Msg::Warning("A meta-element has a negative distortion, trying to scale");
- for (int i=0; i<topPrimVert.size(); ++i) { // Move top primary vert.
- MVertex *&vb = _superVert[baseInd[i]], *&vt = _superVert[topInd[i]];
- SPoint3 pb = vb->point(), pt = vt->point();
- pt += SPoint3(0.25*(pt-pb));
- vt->setXYZ(pt.x(), pt.y(), pt.z());
- }
- for (int i=topPrimVert.size(); i<topInd.size(); ++i) { // Recompute HO vert. in top face
- SPoint3 p;
- const int ind = topInd[i];
- _superEl0->pnt(points(ind,0),points(ind,1),points(ind,2),p);
- _superVert[ind]->setXYZ(p.x(),p.y(),p.z());
- }
- for (int i=0; i<otherInd.size(); ++i) { // Recompute vert. not in base & top faces
- SPoint3 p;
- const int ind = otherInd[i];
- _superEl0->pnt(points(ind,0),points(ind,1),points(ind,2),p);
- _superVert[ind]->setXYZ(p.x(),p.y(),p.z());
- }
- }
-
-}
-
-SuperEl::~SuperEl()
-{
- for (int i = 0; i < _superVert.size(); i++) delete _superVert[i];
- _superVert.clear();
- delete _superEl;
- if(_superEl0) delete _superEl0;
-}
-
-bool SuperEl::isPointIn(const SPoint3 p) const
-{
- double xyz[3] = {p.x(), p.y(), p.z()}, uvw[3];
- _superEl0->xyz2uvw(xyz,uvw);
- return _superEl0->isInside(uvw[0],uvw[1],uvw[2]);
-}
-
-bool SuperEl::straightToCurved(double *xyzS, double *xyzC) const
-{
- double uvw[3];
- _superEl0->xyz2uvw(xyzS,uvw);
- if (!_superEl0->isInside(uvw[0],uvw[1],uvw[2])) return false;
-
- SPoint3 pC;
- _superEl->pnt(uvw[0],uvw[1],uvw[2],pC);
- xyzC[0] = pC[0];
- xyzC[1] = pC[1];
- xyzC[2] = pC[2];
-
- return true;
-}
-
-std::string SuperEl::printPOS()
-{
- std::vector<MVertex*> verts;
- _superEl->getVertices(verts);
-// std::string posStr = _superEl->getStringForPOS();
- std::string posStr = _superEl0->getStringForPOS();
- int n = (posStr[posStr.size()-1]=='2' ? _superEl : _superEl0)->getNumVertices();
-
- std::ostringstream oss;
-
- oss << posStr << "(";
- for(int i = 0; i < n; i++){
- if(i) oss << ",";
- oss << _superVert[i]->x() << "," << _superVert[i]->y() << "," << _superVert[i]->z();
- }
- oss << "){0.";
- for(int i = 1; i < n; i++) oss << ",0.";
- oss << "};\n";
-
- return oss.str();
-}
--
GitLab