Skip to content
Snippets Groups Projects
Commit a4e17957 authored by Thomas Toulorge's avatar Thomas Toulorge
Browse files

Updated fast high-order BL mesh curving algorithm

parent 19db7760
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
contrib/HighOrderMeshOptimizer/CMakeLists.txt
+ 1
− 1
View file @ a4e17957
......@@ -10,7 +10,7 @@ set(SRC
OptHomIntegralBoundaryDist.cpp
OptHomCADDist.cpp
ParamCoord.cpp
SuperEl.cpp
MetaEl.cpp
OptHomElastic.cpp
OptHomFastCurving.cpp
OptHomPeriodicity.cpp
......
contrib/HighOrderMeshOptimizer/SuperEl.cpp → contrib/HighOrderMeshOptimizer/MetaEl.cpp
+ 290
− 0
View file @ a4e17957
......@@ -36,87 +36,119 @@
#include "MPrism.h"
#include "MHexahedron.h"
#include "BasisFactory.h"
#include "SuperEl.h"
#include "MetaEl.h"
std::map<int,SuperEl::superInfoType> SuperEl::_superInfo;
SuperEl::superInfoType::superInfoType(int type, int order)
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)
{
int iBaseFace = 0, iTopFace = 0;
typedef std::vector<int>::iterator VIntIter;
int iBaseFace = 0, iTopFace = 0, nbEdVert = 0;
switch (type) {
case TYPE_QUA:
iBaseFace = 0; iTopFace = 2;
iBaseFace = 0; iTopFace = 2; nbEdVert = 4 + 4*(order-1);
break;
case TYPE_PRI:
iBaseFace = 0; iTopFace = 1;
iBaseFace = 0; iTopFace = 1; nbEdVert = 6 + 9*(order-1);
break;
case TYPE_HEX:
iBaseFace = 0; iTopFace = 5;
iBaseFace = 0; iTopFace = 5; nbEdVert = 8 + 12*(order-1);
break;
default:
Msg::Error("SuperEl not implemented for element of type %d", type);
nV = 0;
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){
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);
// 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);
}
}
SuperEl::SuperEl(int type, int order, const std::vector<MVertex*> &baseVert,
const std::vector<MVertex*> &topPrimVert) : _superEl(0), _superEl0(0)
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,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;
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.nV == 0) return;
if (sInfo.nbVert == 0) return;
// References for easier writing
const int &nV = sInfo.nV;
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)
_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]);
_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:
_superEl0 = new MQuadrangle(_superVert);
_metaEl0 = new MQuadrangle(_metaVert);
break;
case TYPE_PRI:
_superEl0 = new MPrism(_superVert);
_metaEl0 = new MPrism(_metaVert);
break;
case TYPE_HEX:
_superEl0 = new MHexahedron(_superVert);
_metaEl0 = new MHexahedron(_metaVert);
break;
default:
Msg::Error("SuperEl not implemented for element of type %d", type);
......@@ -124,85 +156,108 @@ SuperEl::SuperEl(int type, int order, const std::vector<MVertex*> &baseVert,
}
// Add HO vertices in top face
for (int i=topPrimVert.size(); i<topInd.size(); ++i) {
for (int iV = topPrimVert.size(); iV < topInd.size(); iV++) {
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());
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 not in base & top faces
for (int i=0; i<otherInd.size(); ++i) {
// Add vertices on edges (excluding base and top faces)
for (int iV = 0; iV < edgeInd.size(); iV++) {
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());
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:
_superEl = new MQuadrangleN(_superVert, order);
_metaEl = new MQuadrangleN(_metaVert, order);
break;
case TYPE_PRI:
_superEl = new MPrismN(_superVert, order);
_metaEl = new MPrismN(_metaVert, order);
break;
case TYPE_HEX:
_superEl = new MHexahedronN(_superVert, order);
_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 (_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());
}
}
// // 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());
// }
// }
}
SuperEl::~SuperEl()
MetaEl::~MetaEl()
{
for (int i = 0; i < _superVert.size(); i++) delete _superVert[i];
_superVert.clear();
delete _superEl;
if(_superEl0) delete _superEl0;
for (int i = 0; i < _metaVert.size(); i++) delete _metaVert[i];
_metaVert.clear();
if (_metaEl) delete _metaEl;
if (_metaEl0) delete _metaEl0;
}
bool SuperEl::isPointIn(const SPoint3 p) const
bool MetaEl::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]);
_metaEl0->xyz2uvw(xyz, uvw);
return _metaEl0->isInside(uvw[0], uvw[1], uvw[2]);
}
bool SuperEl::straightToCurved(double *xyzS, double *xyzC) const
bool MetaEl::straightToCurved(double *xyzS, double *xyzC) const
{
double uvw[3];
_superEl0->xyz2uvw(xyzS,uvw);
if (!_superEl0->isInside(uvw[0],uvw[1],uvw[2])) return false;
_metaEl0->xyz2uvw(xyzS, uvw);
if (!_metaEl0->isInside(uvw[0], uvw[1], uvw[2])) return false;
SPoint3 pC;
_superEl->pnt(uvw[0],uvw[1],uvw[2],pC);
_metaEl->pnt(uvw[0], uvw[1], uvw[2], pC);
xyzC[0] = pC[0];
xyzC[1] = pC[1];
xyzC[2] = pC[2];
......@@ -210,20 +265,22 @@ bool SuperEl::straightToCurved(double *xyzS, double *xyzC) const
return true;
}
std::string SuperEl::printPOS()
std::string MetaEl::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();
_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 << _superVert[i]->x() << "," << _superVert[i]->y() << "," << _superVert[i]->z();
oss << _metaVert[i]->x() << "," << _metaVert[i]->y() << ","
<< _metaVert[i]->z();
}
oss << "){0.";
for(int i = 1; i < n; i++) oss << ",0.";
......
contrib/HighOrderMeshOptimizer/SuperEl.h → contrib/HighOrderMeshOptimizer/MetaEl.h
+ 20
− 18
View file @ a4e17957
......@@ -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_
contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp
+ 77
− 85
View file @ a4e17957
......@@ -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();
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment