diff --git a/contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp b/contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp
index a758a22276eaeb9846ef1a8e57dca4004f80afe1..ea07f088b60f3c83f58f5549b86228c5c746f776 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp
@@ -173,7 +173,7 @@ SVector3 getNormalEdge(MVertex *vert, const SVector3 &n,
// Detect whether edge/face is curved, and give normal
bool isCurvedAndNormal(int type, int order, const std::vector<MVertex*> &faceVert,
- const SPoint3 pBar, SVector3 &normal, double &maxDist) {
+ SVector3 &normal, double &maxDist) {
// static const double eps = 1.e-10;
static const double eps = 1.e-6;
@@ -192,24 +192,22 @@ bool isCurvedAndNormal(int type, int order, const std::vector<MVertex*> &faceVer
sxyz[i] += sxyz[j] * f[j];
}
- // Compute (non-oriented) unit normal to straight edge/face and its scale [length]
+ // Compute unit normal to straight edge/face and its scale [length]
double scale;
const SPoint3 &p0 = sxyz[0], &p1 = sxyz[1];
if (type == TYPE_LIN) {
- normal = SVector3(p0.y()-p1.y(),p1.x()-p0.x(),0.);
+// normal = SVector3(p0.y()-p1.y(),p1.x()-p0.x(),0.);
+ normal = SVector3(p1.y()-p0.y(),p0.x()-p1.x(),0.);
scale = normal.normalize();
}
else {
const SPoint3 &p2 = sxyz[2];
SVector3 p01(p0,p1), p02(p0,p2);
- normal = crossprod(p01,p02);
+// normal = crossprod(p01,p02);
+ normal = crossprod(p02,p01);
scale = sqrt(normal.normalize());
}
- // Orient normal with help of pBar
- SVector3 p0C(p0,pBar);
- if (dot(normal,p0C) < 0.) normal *= -1.; // Make straight normal in-going
-
// Calc max. normal dist. from straight to HO points
maxDist = 0.;
for (int iV = nV1; iV < nV; iV++) {
@@ -226,67 +224,6 @@ bool isCurvedAndNormal(int type, int order, const std::vector<MVertex*> &faceVer
-bool getCurvedFace(MElement* badEl,
- const std::map<MVertex*, std::vector<MElement*> > &vertex2elements,
- std::vector<MVertex*> &baseVert, std::vector<SVector3> &normals,
- int &baseType, double &baseMaxDist) {
-
- bool found = false;
- const int dim = badEl->getDim();
- const int order = badEl->getPolynomialOrder();
- const int numFaces = (dim == 2) ? badEl->getNumEdges() : badEl->getNumFaces();
-
- for (int iFace=0; iFace<numFaces; iFace++) { // Loop on edges/faces
-
- // Get face info (type, vertices)
- std::vector<MVertex*> bv; // Vertices of edge/face
- int numPrimVert, type;
- if (dim == 2) {
- badEl->getEdgeVertices(iFace,bv);
- type = TYPE_LIN;
- numPrimVert = 2;
- }
- else {
- badEl->getFaceVertices(iFace,bv);
- MFace face = badEl->getFace(iFace);
- numPrimVert = face.getNumVertices();
- type = (numPrimVert == 3) ? TYPE_TRI : TYPE_QUA;
- }
-
- // Skip face if at least 1 vertex not on boundary
- bool inDom = false;
- for (int iV=0; iV<bv.size(); ++iV)
- if (bv[iV]->onWhat()->dim() == dim) inDom = true;
- if (inDom) continue;
-
- // std::cout << "DBGTT: Checking edge/face " << iFace << " in el. " << badEl->getNum() << ": "
-// << numPrimVert << " vert., type " << type << "\n";
- // If face curved, mark it and compute normals to its primary vertices
- SVector3 n; // Normal to straight edge/face
- double maxDist; // TOFIX: Max. normal. dist. to straight in curved face
- const SPoint3 pBar = badEl->barycenter(true); // Bary. of el. to help orienting normal
- if (isCurvedAndNormal(type,order,bv,pBar,n,maxDist)) { // Check whether edge/face is curved and compute straight normal
-// std::cout << "DBGTT: Curved edge/face in el. " << badEl->getNum() << "\n";
- if (found) { // If more than 1 curved edge/face in bad el., drop it
- std::cout << "DBGTT: More than 1 curved edge/face detected in el. " << badEl->getNum() << "\n";
- return false;
- }
- for (int iV=0; iV<numPrimVert; iV++) // Compute normals to prim. vert. of edge/face
- normals.push_back(getNormalEdge(bv[iV],n,vertex2elements));
- baseVert = bv;
- baseType = type;
- baseMaxDist = maxDist;
- found = true;
- }
-
- }
-
- return found;
-
-}
-
-
-
void makeStraight(MElement *el, const std::set<MVertex*> &movedVert) {
const nodalBasis *nb = el->getFunctionSpace();
@@ -343,27 +280,34 @@ std::set<MElement*> getSuperElBlob(MElement *el, const std::map<MVertex*,
-void curveMesh(std::map<MVertex*, std::vector<MElement *> > &vertex2elements,
- std::set<MElement*> &badElements, FastCurvingParameters &p, int samples)
+void curveMeshFromFaces(std::map<MVertex*, std::vector<MElement *> > &vertex2elements,
+ std::set<MElement*> &faceElements, FastCurvingParameters &p)
{
- const int nbBadElts = badElements.size();
- std::vector<MElement*> badElts;
+ const int nbFaceElts = faceElements.size();
+ std::vector<MElement*> faceElts;
std::vector<SuperEl*> superElts;
- badElts.reserve(nbBadElts);
- superElts.reserve(nbBadElts);
+ faceElts.reserve(nbFaceElts);
+ superElts.reserve(nbFaceElts);
std::ofstream of("dum.pos");
of << "View \" \"{\n";
- for (std::set<MElement*>::const_iterator itBE = badElements.begin(); itBE != badElements.end(); ++itBE) {
- std::vector<MVertex*> faceBE;
- std::vector<SVector3> normalsBE;
- int typeBE;
- double maxDistBE;
- if (getCurvedFace(*itBE,vertex2elements,faceBE,normalsBE,typeBE,maxDistBE)) {
- badElts.push_back(*itBE);
- superElts.push_back(new SuperEl(*itBE,maxDistBE*p.distanceFactor,typeBE,faceBE,normalsBE));
+ for (std::set<MElement*>::const_iterator itFE = faceElements.begin(); itFE != faceElements.end(); ++itFE) {
+ const int dim = (*itFE)->getDim();
+ const int order = (*itFE)->getPolynomialOrder();
+ const int numPrimVert = (*itFE)->getNumPrimaryVertices();
+ const int type = (*itFE)->getType();
+ std::vector<MVertex*> faceVert;
+ (*itFE)->getVertices(faceVert);
+ double maxDist;
+ SVector3 faceNormal;
+ if (isCurvedAndNormal(type,order,faceVert,faceNormal,maxDist)) {
+ std::vector<SVector3> baseNormal;
+ for (int iV=0; iV<numPrimVert; iV++) // Compute normals to prim. vert. of edge/face
+ baseNormal.push_back(getNormalEdge(faceVert[iV],faceNormal,vertex2elements));
+ faceElts.push_back(*itFE);
+ superElts.push_back(new SuperEl(order,maxDist*p.distanceFactor,type,faceVert,baseNormal));
of << superElts.back()->printPOS();
}
}
@@ -372,12 +316,12 @@ void curveMesh(std::map<MVertex*, std::vector<MElement *> > &vertex2elements,
of.close();
std::set<MVertex*> movedVert;
- for (int iBE=0; iBE<badElts.size(); ++iBE) {
- std::set<MElement*> blob = getSuperElBlob(badElts[iBE], vertex2elements, superElts[iBE]);
-// std::cout << "DBGTT: Blob of bad el. " << badElts[iBE]->getNum() << " contains elts.";
+ for (int iFE=0; iFE<faceElts.size(); ++iFE) {
+ std::set<MElement*> blob = getSuperElBlob(faceElts[iFE], vertex2elements, superElts[iFE]);
+// std::cout << "DBGTT: Blob of bad el. " << faceElts[iBE]->getNum() << " contains elts.";
// for (std::set<MElement*>::iterator itE = blob.begin(); itE != blob.end(); ++itE) std::cout << " " << (*itE)->getNum();
// std::cout << "\n";
- makeStraight(badElts[iBE],movedVert); // Make bad. el. straight
+// makeStraight(faceElts[iFE],movedVert); // Make bad. el. straight
for (std::set<MElement*>::iterator itE = blob.begin(); itE != blob.end(); ++itE) {
makeStraight(*itE,movedVert);
for (int i = 0; i < (*itE)->getNumVertices(); ++i) { // For each vert. of each el. in blob
@@ -385,14 +329,14 @@ void curveMesh(std::map<MVertex*, std::vector<MElement *> > &vertex2elements,
MVertex* vert = (*itE)->getVertex(i);
if (movedVert.find(vert) == movedVert.end()) { // If vert. not already moved
double xyzS[3] = {vert->x(), vert->y(), vert->z()}, xyzC[3];
- if (superElts[iBE]->straightToCurved(xyzS,xyzC)) {
+ if (superElts[iFE]->straightToCurved(xyzS,xyzC)) {
// std::cout << "DBGTT: moving vertex " << vert->getNum() << " from (" << xyzS[0] << "," << xyzS[1] << "," << xyzS[2] << ") to (" << xyzC[0] << "," << xyzC[1] << "," << xyzC[2] << ")\n";
vert->setXYZ(xyzC[0],xyzC[1],xyzC[2]);
movedVert.insert(vert);
}
-// else std::cout << "DBGTT: Failed to move vertex " << vert->getNum() << " with bad. el " << badElts[iBE]->getNum() << "\n";
+// else std::cout << "DBGTT: Failed to move vertex " << vert->getNum() << " with bad. el " << faceElts[iBE]->getNum() << "\n";
}
-// else std::cout << "DBGTT: Already moved vertex " << vert->getNum() << " with bad. el " << badElts[iBE]->getNum() << "\n";
+// else std::cout << "DBGTT: Already moved vertex " << vert->getNum() << " with bad. el " << faceElts[iBE]->getNum() << "\n";
}
}
}
@@ -410,37 +354,28 @@ void HighOrderMeshFastCurving(GModel *gm, FastCurvingParameters &p)
double t1 = Cpu();
- int samples = 30;
-
- double tf1 = Cpu();
-
Msg::StatusBar(true, "Optimizing high order mesh...");
std::vector<GEntity*> entities;
gm->getEntities(entities);
+ const int bndDim = p.dim-1;
+ // Compute vert. -> elt. connectivity
+ Msg::Info("Computing connectivity...");
std::map<MVertex*, std::vector<MElement *> > vertex2elements;
- std::set<MElement*> badasses;
+ for (int iEnt = 0; iEnt < entities.size(); ++iEnt)
+ calcVertex2Elements(p.dim,entities[iEnt],vertex2elements);
+
+ // Loop over geometric entities
for (int iEnt = 0; iEnt < entities.size(); ++iEnt) {
GEntity* &entity = entities[iEnt];
- if (entity->dim() != p.dim || (p.onlyVisible && !entity->getVisibility())) continue;
- Msg::Info("Computing connectivity and bad elements for entity %d...",
- entity->tag());
- calcVertex2Elements(p.dim,entity,vertex2elements); // Compute vert. -> elt. connectivity
-// std::cout << "DBGTT: num. el. = " << entity->getNumMeshElements() << "\n";
- for (int iEl = 0; iEl < entity->getNumMeshElements();iEl++) { // Detect bad elements
- double jmin, jmax;
- MElement *el = entity->getMeshElement(iEl);
- if (el->getDim() == p.dim) {
- el->scaledJacRange(jmin, jmax);
-// std::cout << "El. " << iEl << ", jmin = " << jmin << ", jmax = " << jmax << "\n";
-// Msg::Info("El %i, jmin = %g, jmax = %g\n",iEl,jmin,jmax);
- if (jmin < p.BARRIER_MIN || jmax > p.BARRIER_MAX) badasses.insert(el);
- }
- }
+ if (entity->dim() != bndDim || (p.onlyVisible && !entity->getVisibility())) continue;
+ Msg::Info("Curving elements for entity %d...",entity->tag());
+ std::set<MElement*> faceElements;
+ for (int iEl = 0; iEl < entity->getNumMeshElements();iEl++)
+ faceElements.insert(entity->getMeshElement(iEl));
+ curveMeshFromFaces(vertex2elements, faceElements, p);
}
- curveMesh(vertex2elements, badasses, p, samples);
-
double t2 = Cpu();
Msg::StatusBar(true, "Done curving high order mesh (%g s)", t2-t1);
diff --git a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
index 4b19087dde1ddadf92570447ad13aca58b5d052d..4b81c780dd9933cc4c4079c19cc478995e5e8d3d 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
@@ -43,16 +43,6 @@ Mesh::Mesh(const std::map<MElement*,GEntity*> &element2entity,
_dim = (*els.begin())->getDim();
- if (fixBndNodes) {
- if (_dim == 2) _pc = new ParamCoordPhys2D;
- else _pc = new ParamCoordPhys3D;
- Msg::Debug("METHOD: Fixing boundary nodes and using physical coordinates");
- }
- else {
- _pc = new ParamCoordParent;
- Msg::Debug("METHOD: Freeing boundary nodes and using parent parametric coordinates");
- }
-
// Initialize elements, vertices, free vertices and element->vertices
// connectivity
const int nElements = els.size();
@@ -64,6 +54,7 @@ Mesh::Mesh(const std::map<MElement*,GEntity*> &element2entity,
_nNodEl.resize(nElements);
_indPCEl.resize(nElements);
int iEl = 0;
+ bool nonGeoMove = false;
for(std::set<MElement*>::const_iterator it = els.begin();
it != els.end(); ++it, ++iEl) {
_el[iEl] = *it;
@@ -72,31 +63,36 @@ Mesh::Mesh(const std::map<MElement*,GEntity*> &element2entity,
_nNodEl[iEl] = jac->getNumMapNodes();
for (int iVEl = 0; iVEl < jac->getNumMapNodes(); iVEl++) {
MVertex *vert = _el[iEl]->getVertex(iVEl);
+ GEntity *ge = vert->onWhat();
+ const int vDim = ge->dim();
+ const bool hasParam = ge->haveParametrization();
int iV = addVert(vert);
_el2V[iEl].push_back(iV);
- const int nPCV = _pc->nCoord(vert);
- bool isFV = false;
- if (nPCV > 0) {
- if (fixBndNodes)
- isFV = (vert->onWhat()->dim() == _dim) && (toFix.find(vert) == toFix.end());
- else
- isFV = (vert->onWhat()->dim() >= 1) && (toFix.find(vert) == toFix.end());
- }
- if (isFV) {
- int iFV = addFreeVert(vert,iV,nPCV,toFix);
+ if ((vDim > 0) && (toFix.find(vert) == toFix.end()) && (!fixBndNodes || vDim == _dim)) { // Free vertex?
+ ParamCoord *param;
+ if (vDim == 3) param = new ParamCoordPhys3D();
+ else if (hasParam) param = new ParamCoordParent(vert);
+ else {
+ if (vDim == 2) param = new ParamCoordPhys2D(); //todo: make 2d local surf. param
+ else param = new ParamCoordLocalLine(vert);
+ nonGeoMove = true;
+ }
+ int iFV = addFreeVert(vert,iV,vDim,param,toFix);
_el2FV[iEl].push_back(iFV);
- for (int i=_startPCFV[iFV]; i<_startPCFV[iFV]+nPCV; i++)
- _indPCEl[iEl].push_back(i);
+ for (int i=_startPCFV[iFV]; i<_startPCFV[iFV]+vDim; i++) _indPCEl[iEl].push_back(i);
}
else _el2FV[iEl].push_back(-1);
}
}
+ if (nonGeoMove) Msg::Info("WARNING: Some vertices will be moved along local lines "
+ "or planes, they may not remain on the exact geometry");
+
// Initial coordinates
_ixyz.resize(nVert());
for (int iV = 0; iV < nVert(); iV++) _ixyz[iV] = _vert[iV]->point();
_iuvw.resize(nFV());
- for (int iFV = 0; iFV < nFV(); iFV++) _iuvw[iFV] = _pc->getUvw(_freeVert[iFV]);
+ for (int iFV = 0; iFV < nFV(); iFV++) _iuvw[iFV] = _paramFV[iFV]->getUvw(_freeVert[iFV]);
// Set current coordinates
_xyz = _ixyz;
@@ -167,7 +163,7 @@ int Mesh::addVert(MVertex* vert)
}
int Mesh::addFreeVert(MVertex* vert, const int iV, const int nPCV,
- std::set<MVertex*> &toFix)
+ ParamCoord *param, std::set<MVertex*> &toFix)
{
std::vector<MVertex*>::iterator itVert = find(_freeVert.begin(),
_freeVert.end(),vert);
@@ -175,6 +171,7 @@ int Mesh::addFreeVert(MVertex* vert, const int iV, const int nPCV,
const int iStart = (_startPCFV.size() == 0)? 0 : _startPCFV.back()+_nPCFV.back();
const bool forcedV = (vert->onWhat()->dim() < 2) || (toFix.find(vert) != toFix.end());
_freeVert.push_back(vert);
+ _paramFV.push_back(param);
_fv2V.push_back(iV);
_startPCFV.push_back(iStart);
_nPCFV.push_back(nPCV);
@@ -205,7 +202,7 @@ void Mesh::updateMesh(const double *it)
uvwV[0] = *it; it++;
if (_nPCFV[iFV] >= 2) { uvwV[1] = *it; it++; }
if (_nPCFV[iFV] == 3) { uvwV[2] = *it; it++; }
- _xyz[iV] = _pc->uvw2Xyz(_freeVert[iFV],uvwV);
+ _xyz[iV] = _paramFV[iFV]->uvw2Xyz(uvwV);
}
}
@@ -251,7 +248,7 @@ void Mesh::updateGEntityPositions()
for (int iV = 0; iV < nVert(); iV++)
_vert[iV]->setXYZ(_xyz[iV].x(),_xyz[iV].y(),_xyz[iV].z());
for (int iFV = 0; iFV < nFV(); iFV++)
- _pc->exportParamCoord(_freeVert[iFV], _uvw[iFV]);
+ _paramFV[iFV]->exportParamCoord(_uvw[iFV]);
}
void Mesh::metricMinAndGradients(int iEl, std::vector<double> &lambda,
@@ -299,7 +296,7 @@ void Mesh::metricMinAndGradients(int iEl, std::vector<double> &lambda,
gXyzV [l] = SPoint3(gradLambdaB(l,i+0*numMapNodes),
gradLambdaB(l,i+1*numMapNodes),/*BDB(l,i+2*nbNod)*/ 0.);
}
- _pc->gXyz2gUvw(_freeVert[iFVi],_uvw[iFVi],gXyzV,gUvwV);
+ _paramFV[iFVi]->gXyz2gUvw(_uvw[iFVi],gXyzV,gUvwV);
for (int l = 0; l < numJacNodes; l++) {
gradLambda[indGSJ(iEl,l,iPC)] = gUvwV[l][0];
if (_nPCFV[iFVi] >= 2) gradLambda[indGSJ(iEl,l,iPC+1)] = gUvwV[l][1];
@@ -354,7 +351,7 @@ void Mesh::scaledJacAndGradients(int iEl, std::vector<double> &sJ,
for (int l = 0; l < numJacNodes; l++)
gXyzV [l] = SPoint3(BDB(l,i+0*numMapNodes), BDB(l,i+1*numMapNodes),
BDB(l,i+2*numMapNodes));
- _pc->gXyz2gUvw(_freeVert[iFVi],_uvw[iFVi],gXyzV,gUvwV);
+ _paramFV[iFVi]->gXyz2gUvw(_uvw[iFVi],gXyzV,gUvwV);
for (int l = 0; l < numJacNodes; l++) {
gSJ[indGSJ(iEl,l,iPC)] = gUvwV[l][0];
if (_nPCFV[iFVi] >= 2) gSJ[indGSJ(iEl,l,iPC+1)] = gUvwV[l][1];
@@ -371,9 +368,9 @@ void Mesh::pcScale(int iFV, std::vector<double> &scale)
// Calc. derivative of x, y & z w.r.t. parametric coordinates
const SPoint3 dX(1.,0.,0.), dY(0.,1.,0.), dZ(0.,0.,1.);
SPoint3 gX, gY, gZ;
- _pc->gXyz2gUvw(_freeVert[iFV],_uvw[iFV],dX,gX);
- _pc->gXyz2gUvw(_freeVert[iFV],_uvw[iFV],dY,gY);
- _pc->gXyz2gUvw(_freeVert[iFV],_uvw[iFV],dZ,gZ);
+ _paramFV[iFV]->gXyz2gUvw(_uvw[iFV],dX,gX);
+ _paramFV[iFV]->gXyz2gUvw(_uvw[iFV],dY,gY);
+ _paramFV[iFV]->gXyz2gUvw(_uvw[iFV],dZ,gZ);
// Scale = inverse norm. of vector (dx/du, dy/du, dz/du)
scale[0] = 1./sqrt(gX[0]*gX[0]+gY[0]*gY[0]+gZ[0]*gZ[0]);
diff --git a/contrib/HighOrderMeshOptimizer/OptHomMesh.h b/contrib/HighOrderMeshOptimizer/OptHomMesh.h
index 481cc1f96e24e07adbe549430964f078eb9aeb41..ea54b5b2cc41668f3e9bc0652ea990fc3df6d9e0 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomMesh.h
+++ b/contrib/HighOrderMeshOptimizer/OptHomMesh.h
@@ -105,12 +105,12 @@ private:
std::vector<int> _nBezEl, _nNodEl;
// Index of parametric coord. for an el.
std::vector<std::vector<int> > _indPCEl;
-
- ParamCoord *_pc;
+ // Parametrization for a free vertex
+ std::vector<ParamCoord*> _paramFV;
int addVert(MVertex* vert);
int addFreeVert(MVertex* vert, const int iV, const int nPCV,
- std::set<MVertex*> &toFix);
+ ParamCoord *param, std::set<MVertex*> &toFix);
void calcScaledNormalEl2D(const std::map<MElement*,GEntity*> &element2entity, int iEl);
static inline int indJB2DBase(int nNod, int l, int i, int j)
{
@@ -140,7 +140,7 @@ void Mesh::gradDistSq(int iFV, std::vector<double> &gDSq)
{
SPoint3 gXyz = (_xyz[_fv2V[iFV]]-_ixyz[_fv2V[iFV]])*2.;
SPoint3 gUvw;
- _pc->gXyz2gUvw(_freeVert[iFV],_uvw[iFV],gXyz,gUvw);
+ _paramFV[iFV]->gXyz2gUvw(_uvw[iFV],gXyz,gUvw);
gDSq[0] = gUvw[0];
if (_nPCFV[iFV] >= 2) gDSq[1] = gUvw[1];
diff --git a/contrib/HighOrderMeshOptimizer/ParamCoord.cpp b/contrib/HighOrderMeshOptimizer/ParamCoord.cpp
index ff511388a6a03624cc149ae756e6c27ed0e941c5..1361610e1dc03175c3ee71ce2f87b6706671c0f8 100644
--- a/contrib/HighOrderMeshOptimizer/ParamCoord.cpp
+++ b/contrib/HighOrderMeshOptimizer/ParamCoord.cpp
@@ -28,9 +28,11 @@
// Contributor(s): Thomas Toulorge, Jonathan Lambrechts
#include <iostream>
+#include <algorithm>
#include "GFace.h"
#include "GEdge.h"
#include "MVertex.h"
+#include "MLine.h"
#include "ParamCoord.h"
#include "GmshMessage.h"
@@ -44,81 +46,53 @@ SPoint3 ParamCoordParent::getUvw(MVertex* vert)
switch (ge->dim()) {
case 1: {
SPoint3 p(0.,0.,0.);
- reparamMeshVertexOnEdge(vert,static_cast<GEdge*>(ge),p[0]);
+ reparamMeshVertexOnEdge(vert,static_cast<GEdge*>(ge),p[0]); // Overkill if vert. well classified and parametrized
return p;
break;
}
case 2: {
SPoint2 p;
- reparamMeshVertexOnFace(vert,static_cast<GFace*>(ge),p);
+ reparamMeshVertexOnFace(vert,static_cast<GFace*>(ge),p); // Overkill if vert. well classified and parametrized
return SPoint3(p[0],p[1],0.);
break;
}
- case 3: {
- return vert->point();
- break;
- }
}
}
-SPoint3 ParamCoordParent::uvw2Xyz(MVertex* vert, const SPoint3 &uvw)
+SPoint3 ParamCoordParent::uvw2Xyz(const SPoint3 &uvw)
{
- GEntity *ge = vert->onWhat();
-
- switch (ge->dim()) {
- case 1: {
- GPoint gp = static_cast<GEdge*>(ge)->point(uvw[0]);
- return SPoint3(gp.x(),gp.y(),gp.z());
- break;
- }
- case 2: {
- GPoint gp = static_cast<GFace*>(ge)->point(uvw[0],uvw[1]);
- return SPoint3(gp.x(),gp.y(),gp.z());
- break;
- }
- case 3: {
- return uvw;
- break;
- }
- }
+ GEntity *ge = _vert->onWhat();
+ GPoint gp = (ge->dim() == 1) ? static_cast<GEdge*>(ge)->point(uvw[0]) :
+ static_cast<GFace*>(ge)->point(uvw[0],uvw[1]);
+ return SPoint3(gp.x(),gp.y(),gp.z());
}
-void ParamCoordParent::gXyz2gUvw(MVertex* vert, const SPoint3 &uvw,
- const SPoint3 &gXyz, SPoint3 &gUvw)
+void ParamCoordParent::gXyz2gUvw(const SPoint3 &uvw, const SPoint3 &gXyz, SPoint3 &gUvw)
{
- GEntity *ge = vert->onWhat();
- switch (ge->dim()) {
- case 1: {
+ GEntity *ge = _vert->onWhat();
+
+ if (ge->dim() == 1) {
SVector3 der = static_cast<GEdge*>(ge)->firstDer(uvw[0]);
gUvw[0] = gXyz.x() * der.x() + gXyz.y() * der.y() + gXyz.z() * der.z();
- break;
}
- case 2: {
+ else {
Pair<SVector3,SVector3> der = static_cast<GFace*>(ge)->firstDer(SPoint2(uvw[0],uvw[1]));
gUvw[0] = gXyz.x() * der.first().x() + gXyz.y() * der.first().y() +
gXyz.z() * der.first().z();
gUvw[1] = gXyz.x() * der.second().x() + gXyz.y() * der.second().y() +
gXyz.z() * der.second().z();
- break;
- }
- case 3: {
- gUvw = gXyz;
- break;
- }
}
}
-void ParamCoordParent::gXyz2gUvw(MVertex* vert, const SPoint3 &uvw,
- const std::vector<SPoint3> &gXyz,
- std::vector<SPoint3> &gUvw)
+void ParamCoordParent::gXyz2gUvw(const SPoint3 &uvw,
+ const std::vector<SPoint3> &gXyz, std::vector<SPoint3> &gUvw)
{
- GEntity *ge = vert->onWhat();
+ GEntity *ge = _vert->onWhat();
- switch (ge->dim()) {
- case 1: {
+ if (ge->dim() == 1) {
SVector3 der = static_cast<GEdge*>(ge)->firstDer(uvw[0]);
std::vector<SPoint3>::iterator itUvw = gUvw.begin();
for (std::vector<SPoint3>::const_iterator itXyz=gXyz.begin();
@@ -126,9 +100,8 @@ void ParamCoordParent::gXyz2gUvw(MVertex* vert, const SPoint3 &uvw,
(*itUvw)[0] = itXyz->x() * der.x() + itXyz->y() * der.y() + itXyz->z() * der.z();
itUvw++;
}
- break;
}
- case 2: {
+ else {
Pair<SVector3,SVector3> der = static_cast<GFace*>(ge)->firstDer
(SPoint2(uvw[0],uvw[1]));
std::vector<SPoint3>::iterator itUvw=gUvw.begin();
@@ -140,24 +113,25 @@ void ParamCoordParent::gXyz2gUvw(MVertex* vert, const SPoint3 &uvw,
itXyz->z() * der.second().z();
itUvw++;
}
- break;
- }
- case 3: {
- std::vector<SPoint3>::iterator itUvw=gUvw.begin();
- for (std::vector<SPoint3>::const_iterator itXyz=gXyz.begin();
- itXyz != gXyz.end(); itXyz++) {
- *itUvw = *itXyz;
- itUvw++;
- }
- break;
- }
}
}
-void ParamCoordParent::exportParamCoord(MVertex *v, const SPoint3 &uvw)
+ParamCoordLocalLine::ParamCoordLocalLine(MVertex* v) : dir(0.)
{
- for (int d = 0; d < v->onWhat()->dim(); ++d) {
- v->setParameter(d, uvw[d]);
+
+ GEdge *ge = static_cast<GEdge*>(v->onWhat());
+
+ for (std::vector<MLine*>::iterator it = ge->lines.begin(); it != ge->lines.end(); it++) {
+ std::vector<MVertex*> lVerts;
+ (*it)->getVertices(lVerts);
+ if (std::find(lVerts.begin(),lVerts.end(),v) == lVerts.end()) {
+ SVector3 tan(lVerts[0]->point(),lVerts[1]->point());
+ tan.normalize();
+ dir += tan;
+ }
}
+
+ dir.normalize();
+
}
diff --git a/contrib/HighOrderMeshOptimizer/ParamCoord.h b/contrib/HighOrderMeshOptimizer/ParamCoord.h
index 9b50e50280a609c0dd36b8cb078074538d32276f..3c509cb83ed47896967bee3ab28465c9a70de036 100644
--- a/contrib/HighOrderMeshOptimizer/ParamCoord.h
+++ b/contrib/HighOrderMeshOptimizer/ParamCoord.h
@@ -33,18 +33,16 @@
class ParamCoord
{
public:
- virtual void exportParamCoord(MVertex *v, const SPoint3 &uvw) {};
- // Number of parametric coordinates for vertex
- virtual int nCoord(MVertex* vert) = 0;
+ // Set param. coord. of MVertex if applicable
+ virtual void exportParamCoord(const SPoint3 &uvw) {}
// Get parametric coordinates of vertex
- virtual SPoint3 getUvw(MVertex* vert) = 0;
+ virtual SPoint3 getUvw(MVertex* v) = 0;
// Calculate physical coordinates from parametric coordinates of vertex
- virtual SPoint3 uvw2Xyz(MVertex* vert, const SPoint3 &uvw) = 0;
+ virtual SPoint3 uvw2Xyz(const SPoint3 &uvw) = 0;
// Calculate derivatives w.r.t parametric coordinates
- virtual void gXyz2gUvw(MVertex* vert, const SPoint3 &uvw,
- const SPoint3 &gXyz, SPoint3 &gUvw) = 0;
+ virtual void gXyz2gUvw(const SPoint3 &uvw, const SPoint3 &gXyz, SPoint3 &gUvw) = 0;
// Calculate derivatives w.r.t parametric coordinates
- virtual void gXyz2gUvw(MVertex* vert, const SPoint3 &uvw,
+ virtual void gXyz2gUvw(const SPoint3 &uvw,
const std::vector<SPoint3> &gXyz, std::vector<SPoint3> &gUvw) = 0;
virtual ~ParamCoord() {}
};
@@ -52,13 +50,11 @@ public:
class ParamCoordPhys3D : public ParamCoord
{
public:
- int nCoord(MVertex* vert) { return 3; }
- virtual SPoint3 getUvw(MVertex* vert) { return vert->point(); }
- virtual SPoint3 uvw2Xyz(MVertex* vert, const SPoint3 &uvw) { return uvw; }
- virtual void gXyz2gUvw(MVertex* vert, const SPoint3 &uvw, const SPoint3 &gXyz,
+ SPoint3 getUvw(MVertex* v) { return v->point(); }
+ SPoint3 uvw2Xyz(const SPoint3 &uvw) { return uvw; }
+ void gXyz2gUvw(const SPoint3 &uvw, const SPoint3 &gXyz,
SPoint3 &gUvw){ gUvw = gXyz; }
- virtual void gXyz2gUvw(MVertex* vert, const SPoint3 &uvw,
- const std::vector<SPoint3> &gXyz, std::vector<SPoint3> &gUvw)
+ void gXyz2gUvw(const SPoint3 &uvw, const std::vector<SPoint3> &gXyz, std::vector<SPoint3> &gUvw)
{
std::vector<SPoint3>::iterator itUvw=gUvw.begin();
for (std::vector<SPoint3>::const_iterator itXyz=gXyz.begin(); itXyz != gXyz.end();
@@ -72,13 +68,10 @@ public:
class ParamCoordPhys2D : public ParamCoord
{
public:
- int nCoord(MVertex* vert) { return 2; }
- virtual SPoint3 getUvw(MVertex* vert) { return vert->point(); }
- virtual SPoint3 uvw2Xyz(MVertex* vert, const SPoint3 &uvw) { return uvw; }
- virtual void gXyz2gUvw(MVertex* vert, const SPoint3 &uvw, const SPoint3 &gXyz,
- SPoint3 &gUvw) { gUvw = gXyz; }
- virtual void gXyz2gUvw(MVertex* vert, const SPoint3 &uvw,
- const std::vector<SPoint3> &gXyz, std::vector<SPoint3> &gUvw)
+ SPoint3 getUvw(MVertex* v) { return v->point(); }
+ SPoint3 uvw2Xyz(const SPoint3 &uvw) { return uvw; }
+ void gXyz2gUvw(const SPoint3 &uvw, const SPoint3 &gXyz, SPoint3 &gUvw) { gUvw = gXyz; }
+ void gXyz2gUvw(const SPoint3 &uvw, const std::vector<SPoint3> &gXyz, std::vector<SPoint3> &gUvw)
{
std::vector<SPoint3>::iterator itUvw=gUvw.begin();
for (std::vector<SPoint3>::const_iterator itXyz=gXyz.begin(); itXyz != gXyz.end();
@@ -92,14 +85,37 @@ public:
class ParamCoordParent : public ParamCoord
{
public:
- int nCoord(MVertex* vert) { return vert->onWhat()->haveParametrization() ? vert->onWhat()->dim() : 0; }
- virtual void exportParamCoord(MVertex *v, const SPoint3 &uvw);
- virtual SPoint3 getUvw(MVertex* vert);
- virtual SPoint3 uvw2Xyz(MVertex* vert, const SPoint3 &uvw);
- virtual void gXyz2gUvw(MVertex* vert, const SPoint3 &uvw, const SPoint3 &gXyz,
- SPoint3 &gUvw);
- virtual void gXyz2gUvw(MVertex* vert, const SPoint3 &uvw,
- const std::vector<SPoint3> &gXyz, std::vector<SPoint3> &gUvw);
+ ParamCoordParent(MVertex* v) : _vert(v) {}
+ void exportParamCoord(const SPoint3 &uvw) {
+ for (int d = 0; d < _vert->onWhat()->dim(); ++d) _vert->setParameter(d, uvw[d]);
+ }
+ SPoint3 getUvw(MVertex* v);
+ SPoint3 uvw2Xyz(const SPoint3 &uvw);
+ void gXyz2gUvw(const SPoint3 &uvw, const SPoint3 &gXyz, SPoint3 &gUvw);
+ void gXyz2gUvw(const SPoint3 &uvw, const std::vector<SPoint3> &gXyz, std::vector<SPoint3> &gUvw);
+protected:
+ MVertex *_vert;
+};
+
+class ParamCoordLocalLine : public ParamCoord
+{
+public:
+ ParamCoordLocalLine(MVertex* v);
+ SPoint3 getUvw(MVertex* v) { return SPoint3(0.,0.,0.); }
+ SPoint3 uvw2Xyz(const SPoint3 &uvw) { return SPoint3(uvw[0]*dir[0],uvw[0]*dir[1],uvw[0]*dir[2]); }
+ void gXyz2gUvw(const SPoint3 &uvw, const SPoint3 &gXyz, SPoint3 &gUvw) {
+ gUvw[0] = gXyz.x()*dir[0] + gXyz.y()*dir[1] + gXyz.z()*dir[2];
+ }
+ void gXyz2gUvw(const SPoint3 &uvw, const std::vector<SPoint3> &gXyz, std::vector<SPoint3> &gUvw) {
+ std::vector<SPoint3>::iterator itUvw = gUvw.begin();
+ for (std::vector<SPoint3>::const_iterator itXyz=gXyz.begin();
+ itXyz != gXyz.end(); itXyz++) {
+ (*itUvw)[0] = itXyz->x()*dir[0] + itXyz->y()*dir[1] + itXyz->z()*dir[2];
+ itUvw++;
+ }
+ }
+protected:
+ SVector3 dir;
};
#endif
diff --git a/contrib/HighOrderMeshOptimizer/SuperEl.cpp b/contrib/HighOrderMeshOptimizer/SuperEl.cpp
index 76a4bece92f38d73df27007a0ef4a5dce07dabc2..d564553c8d3af23efb1939d8303c74936ade0ac4 100644
--- a/contrib/HighOrderMeshOptimizer/SuperEl.cpp
+++ b/contrib/HighOrderMeshOptimizer/SuperEl.cpp
@@ -38,20 +38,20 @@
-SuperEl::SuperEl(MElement *badEl, double dist, int type, const std::vector<MVertex*> &baseVert,
+SuperEl::SuperEl(int order, double dist, int type, const std::vector<MVertex*> &baseVert,
const std::vector<SVector3> &normals) {
// std::cout << "DBGTT: badEl = " << badEl->getNum() << ", _superEl = " << _superEl << std::endl;
switch (type) {
case TYPE_LIN:
- createSuperElQuad(badEl, dist, baseVert, normals[0], normals[1]);
+ createSuperElQuad(order, dist, baseVert, normals[0], normals[1]);
break;
case TYPE_TRI:
- createSuperElPrism(badEl, dist, baseVert, normals[0], normals[1], normals[2]);
+ createSuperElPrism(order, dist, baseVert, normals[0], normals[1], normals[2]);
break;
case TYPE_QUA:
- createSuperElHex(badEl, dist, baseVert, normals[0], normals[1], normals[2], normals[3]);
+ createSuperElHex(order, dist, baseVert, normals[0], normals[1], normals[2], normals[3]);
break;
default:
std::cout << "ERROR: SuperEl not implemented for element of type " << type << std::endl;
@@ -65,12 +65,12 @@ SuperEl::SuperEl(MElement *badEl, double dist, int type, const std::vector<MVert
-void SuperEl::createSuperElQuad(MElement *badEl, double dist, const std::vector<MVertex*> &baseVert,
+void SuperEl::createSuperElQuad(int order, double dist, const std::vector<MVertex*> &baseVert,
const SVector3 &n0, const SVector3 &n1) {
if (baseVert.empty()) {
- std::cout << "ERROR: Base edge not given for SuperEl linked to bad. el. " << badEl->getNum() << "\n";
+ std::cout << "ERROR: Base edge not given for SuperEl\n";
_superEl = 0;
return;
}
@@ -90,8 +90,7 @@ void SuperEl::createSuperElQuad(MElement *badEl, double dist, const std::vector<
// Get basis functions
MQuadrangle quad1(v0,v1,v2,v3);
- const int p = badEl->getPolynomialOrder();
- const nodalBasis *quadBasis = quad1.getFunctionSpace(p);
+ const nodalBasis *quadBasis = quad1.getFunctionSpace(order);
const int Nv = quadBasis->getNumShapeFunctions();
// Store/create all vertices for super-element
@@ -100,8 +99,8 @@ void SuperEl::createSuperElQuad(MElement *badEl, double dist, const std::vector<
_superVert[1] = v1;
_superVert[2] = v2;
_superVert[3] = v3;
- for (int i=2; i<p+1; ++i) _superVert[2+i] = new MVertex(*baseVert[i]); // HO vertices on base edge
- for (int i=3+p; i<Nv; ++i) { // Additional HO vertices
+ for (int i=2; i<order+1; ++i) _superVert[2+i] = new MVertex(*baseVert[i]); // HO vertices on base edge
+ for (int i=3+order; i<Nv; ++i) { // Additional HO vertices
SPoint3 p;
quad1.pnt(quadBasis->points(i,0),quadBasis->points(i,1),0.,p);
_superVert[i] = new MVertex(p.x(),p.y(),0.);
@@ -109,7 +108,7 @@ void SuperEl::createSuperElQuad(MElement *badEl, double dist, const std::vector<
// << "," << _superVert[i]->y()<< "," << _superVert[i]->z() << ")" << std::endl;
}
- _superEl = new MQuadrangleN(_superVert,p);
+ _superEl = new MQuadrangleN(_superVert,order);
_superEl0 = new MQuadrangle(_superVert[0],_superVert[1],_superVert[2],_superVert[3]);
// std::cout << "Created SuperEl at address " << _superEl << std::endl;
@@ -117,14 +116,13 @@ void SuperEl::createSuperElQuad(MElement *badEl, double dist, const std::vector<
-void SuperEl::createSuperElPrism(MElement *badEl, double dist, const std::vector<MVertex*> &baseVert,
- const SVector3 &n0, const SVector3 &n1, const SVector3 &n2) {
+void SuperEl::createSuperElPrism(int order, double dist, const std::vector<MVertex*> &baseVert,
+ const SVector3 &n0, const SVector3 &n1, const SVector3 &n2) {
if (baseVert.empty()) {
- std::cout << "ERROR: Base edge not given for SuperEl linked to bad. el. " << badEl->getNum() << "\n";
+ std::cout << "ERROR: Base edge not given for SuperEl\n";
_superEl = 0;
- _superEl0 = 0;
return;
}
@@ -141,12 +139,11 @@ void SuperEl::createSuperElPrism(MElement *badEl, double dist, const std::vector
// Get basis functions
MPrism prism1(v0,v1,v2,v3,v4,v5);
- const int p = badEl->getPolynomialOrder();
- const nodalBasis *prismBasis = prism1.getFunctionSpace(p);
+ const nodalBasis *prismBasis = prism1.getFunctionSpace(order);
const int Nv = prismBasis->getNumShapeFunctions(); // Number of vertices in HO prism
// Store/create all vertices for super-element
- if (p == 2) {
+ if (order == 2) {
_superVert.resize(18);
_superVert[0] = v0; // First-order vertices
_superVert[1] = v1;
@@ -167,7 +164,7 @@ void SuperEl::createSuperElPrism(MElement *badEl, double dist, const std::vector
_superEl = new MPrism18(_superVert);
}
else {
- std::cout << "ERROR: Prismatic superEl. of order " << p << " not supported for bad. el. " << badEl->getNum() << "\n";
+ std::cout << "ERROR: Prismatic superEl. of order " << order << " not supported\n";
_superEl = 0;
_superEl0 = 0;
return;
@@ -180,14 +177,13 @@ void SuperEl::createSuperElPrism(MElement *badEl, double dist, const std::vector
-void SuperEl::createSuperElHex(MElement *badEl, double dist, const std::vector<MVertex*> &baseVert,
+void SuperEl::createSuperElHex(int order, double dist, const std::vector<MVertex*> &baseVert,
const SVector3 &n0, const SVector3 &n1, const SVector3 &n2, const SVector3 &n3) {
if (baseVert.empty()) {
- std::cout << "ERROR: Base edge not given for SuperEl linked to bad. el. " << badEl->getNum() << "\n";
+ std::cout << "ERROR: Base edge not given for SuperEl\n";
_superEl = 0;
- _superEl0 = 0;
return;
}
@@ -207,12 +203,11 @@ void SuperEl::createSuperElHex(MElement *badEl, double dist, const std::vector<M
// Get basis functions
MHexahedron *hex1 = new MHexahedron(v0,v1,v2,v3,v4,v5,v6,v7);
- const int p = badEl->getPolynomialOrder();
- const nodalBasis *prismBasis = hex1->getFunctionSpace(p);
+ const nodalBasis *prismBasis = hex1->getFunctionSpace(order);
const int Nv = prismBasis->getNumShapeFunctions(); // Number of vertices in HO hex
// Store/create all vertices for super-element
- if (p == 2) {
+ if (order == 2) {
_superVert.resize(27);
_superVert[0] = v0; // First-order vertices
_superVert[1] = v1;
@@ -237,7 +232,7 @@ void SuperEl::createSuperElHex(MElement *badEl, double dist, const std::vector<M
_superEl = new MHexahedron27(_superVert);
}
else {
- std::cout << "ERROR: Hex. superEl. of order " << p << " not supported for bad. el. " << badEl->getNum() << "\n";
+ std::cout << "ERROR: Hex. superEl. of order " << order << " not supported\n";
_superEl = 0;
_superEl0 = 0;
return;
diff --git a/contrib/HighOrderMeshOptimizer/SuperEl.h b/contrib/HighOrderMeshOptimizer/SuperEl.h
index 595394011fd9242d1c6f5066ea438165e86ecdfc..6d129bd0adadbf31360b9d768e51d29500a7ebc4 100644
--- a/contrib/HighOrderMeshOptimizer/SuperEl.h
+++ b/contrib/HighOrderMeshOptimizer/SuperEl.h
@@ -39,7 +39,7 @@ class SuperEl
{
public:
- SuperEl(MElement *badEl, double dist, int type, const std::vector<MVertex*> &baseVert,
+ SuperEl(int order, double dist, int type, const std::vector<MVertex*> &baseVert,
const std::vector<SVector3> &normals);
~SuperEl() { _superVert.clear(); delete _superEl, _superEl0; }
@@ -63,11 +63,11 @@ private:
std::vector<MVertex*> _superVert;
MElement *_superEl, *_superEl0;
- void createSuperElQuad(MElement *badEl, double dist, const std::vector<MVertex*> &baseVert,
+ void createSuperElQuad(int order, double dist, const std::vector<MVertex*> &baseVert,
const SVector3 &n0, const SVector3 &n1);
- void createSuperElPrism(MElement *badEl, double dist, const std::vector<MVertex*> &baseVert,
+ void createSuperElPrism(int order, double dist, const std::vector<MVertex*> &baseVert,
const SVector3 &n0, const SVector3 &n1, const SVector3 &n2);
- void createSuperElHex(MElement *badEl, double dist, const std::vector<MVertex*> &baseVert,
+ void createSuperElHex(int order, double dist, const std::vector<MVertex*> &baseVert,
const SVector3 &n0, const SVector3 &n1, const SVector3 &n2, const SVector3 &n3);