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

Improvements to high-order mesh generation: worked on parametric coordinates... 

Improvements to high-order mesh generation: worked on parametric coordinates for optimization (to be completed), changed fast curving algorithm to sweep over faces
parent 40d47785
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Showing
with 183 additions and 266 deletions
contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp
+ 47
− 112
View file @ d6537ff2
......@@ -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,36 +354,27 @@ 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();
......
contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
+ 28
− 31
View file @ d6537ff2
......@@ -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 ((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;
}
if (isFV) {
int iFV = addFreeVert(vert,iV,nPCV,toFix);
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]);
......
contrib/HighOrderMeshOptimizer/OptHomMesh.h
+ 4
− 4
View file @ d6537ff2
......@@ -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];
......
contrib/HighOrderMeshOptimizer/ParamCoord.cpp
+ 34
− 60
View file @ d6537ff2
......@@ -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]);
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());
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;
}
}
}
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();
}
contrib/HighOrderMeshOptimizer/ParamCoord.h
+ 45
− 29
View file @ d6537ff2
......@@ -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
contrib/HighOrderMeshOptimizer/SuperEl.cpp
+ 21
− 26
View file @ d6537ff2
......@@ -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,
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;
......
contrib/HighOrderMeshOptimizer/SuperEl.h
+ 4
− 4
View file @ d6537ff2
......@@ -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);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment