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

Simplified fast curving (only 2D at the moment)

parent d2aea603
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp
+ 146
− 315
View file @ dfd86002
......@@ -55,304 +55,54 @@ namespace {
// Compute vertex -> element connectivity
void calcVertex2Elements(int dim, GEntity *entity,
std::map<MVertex*, std::vector<MElement *> > &vertex2elements)
{
for (size_t i = 0; i < entity->getNumMeshElements(); ++i) {
MElement *element = entity->getMeshElement(i);
if (element->getDim() == dim)
for (int j = 0; j < element->getNumPrimaryVertices(); ++j)
vertex2elements[element->getVertex(j)].push_back(element);
}
}
// Given a starting vertex and a direction, find the next vertex in a column
MVertex *findVertFromDir(const std::map<MVertex*, std::vector<MElement*> > &vertex2elements,
MVertex *v0, MVertex *vExclude, const SVector3 &dir, double spLimit)
{
const std::vector<MElement*> &elts = (*vertex2elements.find(v0)).second; // All elements connected to vertex
MVertex *vFound = 0;
double spMax = 0.;
for (std::vector<MElement*>::const_iterator itEl = elts.begin(); itEl != elts.end(); ++itEl) // Loop over all elements
for (int i=0; i<(*itEl)->getNumEdges(); i++) { // Loop over all edges of each element
MEdge edge = (*itEl)->getEdge(i);
MVertex *vA = edge.getVertex(0), *vB = edge.getVertex(1), *vTmp;
if (v0 == vA) // Skip if edge unconnected to vertex...
if (vB == vExclude) continue; // ... or if connecting to preceding vertex
else vTmp = vB;
else if (v0 == vB)
if (vA == vExclude) continue;
else vTmp = vA;
else continue;
SVector3 tanVec = edge.tangent();
double sp = fabs(dot(tanVec, dir));
if ((sp > spMax) && (sp > spLimit)){ // Retain the edge giving max. dot product with normal
spMax = sp;
vFound = vTmp;
}
}
return vFound;
}
// Find a column of vertices starting with a given vertex
template<class bndType>
bool findColumn(const std::map<MVertex*, std::vector<MElement*> > &vertex2elements,
MVertex *baseVert, const SVector3 &norm,
const bndType &baseEd, std::vector<MVertex*> &col, int maxNumLayers)
{
static const double spMin = 0.866025404; // Threshold in cosine between 1st edge and normal
static const double spTol = 0.999; // Threshold in cosine between an edge and 1st edge
MVertex *vert = findVertFromDir(vertex2elements, baseVert, 0, norm, spMin);
if (!vert) { // If 1st normal edge not found...
vert = findVertFromDir(vertex2elements, baseVert, 0, baseEd.normal(), spMin); // ... tries normal to base edge/face
if (!vert) {
Msg::Warning("Could not find normal edge for vertex %i ", baseVert->getNum());
return false;
}
}
SVector3 firstEdgeDir(baseVert->point(), vert->point());
firstEdgeDir.normalize();
SVector3 dir(baseVert->point(), vert->point());
col.push_back(vert);
MVertex *oldVert = baseVert, *newVert;
for (int iLayer = 0; iLayer < maxNumLayers; iLayer++) {
newVert = findVertFromDir(vertex2elements, vert, oldVert, firstEdgeDir, spTol);
col.push_back(newVert); // Will add null pointer at the end, which is OK
if (!newVert) break;
oldVert = vert;
vert = newVert;
}
return (col.size() > 1);
}
// Add normal to edge/face to data structure for vertices
void addNormVertEl(MElement *el, const SVector3 &norm,
std::map<MVertex*, SVector3> &normVert)
{
for(unsigned int iV = 0; iV<el->getNumVertices(); iV++) {
MVertex *vert = el->getVertex(iV);
std::pair<std::map<MVertex*, SVector3>::iterator, bool> insNorm = // If nothing for vert...
normVert.insert(std::pair<MVertex*, SVector3>(vert, SVector3(0.))); // ... create zero normal
insNorm.first->second += norm; // Cumulate normals
}
}
typedef std::map<MEdge, std::vector<MElement*>, Less_Edge> MEdgeVecMEltMap;
typedef std::map<MFace, std::vector<MElement*>, Less_Face> MFaceVecMEltMap;
// Compute normals to boundary edges/faces and store them per vertex
void buildNormals(const std::map<MVertex*, std::vector<MElement*> > &vertex2elements,
const std::multimap<GEntity*,GEntity*> &entities, const FastCurvingParameters &p,
std::map<GEntity*, std::map<MVertex*, SVector3> > &normVertEnt)
// Compute edge -> element connectivity (for 2D elements)
void calcEdge2Elements(GEntity *entity, MEdgeVecMEltMap &ed2el)
{
normVertEnt.clear();
for (std::multimap<GEntity*,GEntity*>::const_iterator itBE = entities.begin();
itBE != entities.end(); itBE++) { // Loop over model entities
GEntity* const &bndEnt = itBE->second;
std::map<MVertex*, SVector3> &normVert = normVertEnt[bndEnt];
if (bndEnt->dim() == 1) {
GEdge *ge = bndEnt->cast2Edge();
for(unsigned int iEl = 0; iEl<ge->lines.size(); iEl++) {
MLine* const &line = ge->lines[iEl];
SVector3 norm = line->getEdge(0).normal();
addNormVertEl(line, norm, normVert);
}
}
else if (bndEnt->dim() == 2) {
GFace *gf = bndEnt->cast2Face();
for(unsigned int iEl = 0; iEl<gf->triangles.size(); iEl++) {
MTriangle* const &tri = gf->triangles[iEl];
SVector3 norm = tri->getFace(0).normal();
addNormVertEl(tri, norm, normVert);
}
for(unsigned int iEl = 0; iEl<gf->quadrangles.size(); iEl++) {
MQuadrangle* const &quad = gf->quadrangles[iEl];
SVector3 norm = quad->getFace(0).normal();
addNormVertEl(quad, norm, normVert);
for (size_t iEl = 0; iEl < entity->getNumMeshElements(); iEl++) {
MElement *elt = entity->getMeshElement(iEl);
if (elt->getDim() == 2)
for (int iEdge = 0; iEdge < elt->getNumEdges(); iEdge++) {
ed2el[elt->getEdge(iEdge)].push_back(elt);
}
}
for (std::map<MVertex*, SVector3> ::iterator itNV = // Re-normalize for each geom. entity
normVert.begin(); itNV != normVert.end(); itNV++)
itNV->second.normalize();
}
}
// Get first domain element for a given boundary edge
MElement *getFirstEl(const std::map<MVertex*, std::vector<MElement*> > &vertex2elements,
const MEdge &baseEd)
// Compute face -> element connectivity (for 3D elements)
void calcFace2Elements(GEntity *entity, MFaceVecMEltMap &face2el)
{
MVertex *vb0 = baseEd.getVertex(0), *vb1 = baseEd.getVertex(1);
const std::vector<MElement*> &nb0 = vertex2elements.find(vb0)->second;
const std::vector<MElement*> &nb1 = vertex2elements.find(vb1)->second;
MElement *firstEl = 0;
for (int iEl=0; iEl<nb0.size(); iEl++) {
if (find(nb1.begin(), nb1.end(), nb0[iEl]) != nb1.end()) {
firstEl = nb0[iEl];
break;
for (size_t iEl = 0; iEl < entity->getNumMeshElements(); iEl++) {
MElement *elt = entity->getMeshElement(iEl);
if (elt->getDim() == 3)
for (int iFace = 0; iFace < elt->getNumFaces(); iFace++)
face2el[elt->getFace(iFace)].push_back(elt);
}
}
return firstEl;
}
// Get intersection of two vectors of MElement*
void intersect(const std::vector<MElement*> &vi1, const std::vector<MElement*> &vi2,
std::vector<MElement*> &vo)
{
vo.clear();
for (std::vector<MElement*>::const_iterator it = vi1.begin(); it != vi1.end(); it++)
if (std::find(vi2.begin(), vi2.end(), *it) != vi2.end())
vo.push_back(*it);
}
// Get first domain element for a given boundary face
MElement *getFirstEl(const std::map<MVertex*, std::vector<MElement*> > &vertex2elements,
const MFace &baseFace)
// Get the face index of an element given a MFace
inline int getElementFace(const MFace &face, MElement *el)
{
MVertex *vb0 = baseFace.getVertex(0), *vb1 = baseFace.getVertex(1);
// Get elements common to vertices 0 and 1
const std::vector<MElement*> &nb0 = vertex2elements.find(vb0)->second;
const std::vector<MElement*> &nb1 = vertex2elements.find(vb1)->second;
std::vector<MElement*> nb01;
intersect(nb0, nb1, nb01);
if (nb01.empty()) return 0;
if (nb01.size() == 1) return nb01.front();
// Get elements common to vertices 0, 1 and 2
MVertex *vb2 = baseFace.getVertex(2);
const std::vector<MElement*> &nb2 = vertex2elements.find(vb2)->second;
std::vector<MElement*> nb012;
intersect(nb01, nb2, nb012);
if (nb012.empty()) return 0;
if (nb012.size() == 1) return nb012.front();
// If quad, get elements common to all 4 vertices
if (baseFace.getNumVertices() == 4) {
MVertex *vb3 = baseFace.getVertex(3);
const std::vector<MElement*> &nb3 = vertex2elements.find(vb3)->second;
std::vector<MElement*> nb0123;
intersect(nb012, nb3, nb0123);
if (nb0123.empty()) return 0;
if (nb0123.size() == 1) return nb0123.front();
}
// Too many elements, return error
return 0;
for (int iFace = 0; iFace < el->getNumFaces(); iFace++)
if (el->getFace(iFace) == face) return iFace;
return -1;
}
// Get base vertices (in order of first element) for an edge
bool getBaseVertices(const MEdge &baseEd, MElement *firstEl, std::vector<MVertex*> &baseVert)
// Get the edge index of an element given a MEdge
inline int getElementEdge(const MEdge &edge, MElement *el)
{
baseVert.clear();
for (int iEd = 0; iEd < firstEl->getNumEdges(); iEd++)
if (firstEl->getEdge(iEd) == baseEd)
firstEl->getEdgeVertices(iEd, baseVert);
if (baseVert.empty()) {
Msg::Error("Base edge vertices not found in getBaseVertices");
return false;
}
else
return true;
}
// Get base vertices (in order of first element) for a face
bool getBaseVertices(const MFace &baseFace, MElement *firstEl, std::vector<MVertex*> &baseVert)
{
baseVert.clear();
for (int iFace = 0; iFace < firstEl->getNumFaces(); iFace++)
if (firstEl->getFace(iFace) == baseFace)
firstEl->getFaceVertices(iFace, baseVert);
if (baseVert.empty()) {
Msg::Error("Base face vertices not found in getBaseVertices");
return false;
}
else
return true;
}
// Top primary vertices (in order corresponding to the primary base vertices)
template<class bndType>
bool getTopPrimVertices(std::map<MVertex*, std::vector<MElement *> > &vertex2elements,
const std::map<MVertex*, SVector3> &normVert,
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++) {
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;
}
// 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 MetaEl *supEl, std::set<MElement*> &blob)
{
static const int maxDepth = 100;
typedef std::list<MElement*> elLType;
typedef std::vector<MElement*> elVType;
// Build blob by looping over layers of elements (assuming that if an el is too far, its neighbours are too far as well)
blob.clear();
elLType currentLayer, lastLayer;
blob.insert(firstEl);
lastLayer.push_back(firstEl);
for (int d = 0; d < maxDepth; ++d) { // Loop over layers
currentLayer.clear();
for (elLType::iterator it = lastLayer.begin(); it != lastLayer.end(); ++it) { // Loop over elements of last layer
for (int i = 0; i < (*it)->getNumPrimaryVertices(); ++i) {
const elVType &neighbours = vertex2elements.find((*it)->getVertex(i))->second;
for (elVType::const_iterator itN = neighbours.begin(); itN != neighbours.end(); ++itN) { // Loop over neighbours of last layer
SPoint3 p = (*itN)->barycenter(true);
if (supEl->isPointIn(p)) // Add neighbour to blob if inside test element
if (blob.insert(*itN).second) currentLayer.push_back(*itN); // Add to current layer if new in blob
}
}
}
if (currentLayer.empty()) break; // Finished if no new element in blob
lastLayer = currentLayer;
}
for (int iEdge = 0; iEdge < el->getNumEdges(); iEdge++)
if (el->getEdge(iEdge) == edge) return iEdge;
return -1;
}
......@@ -383,7 +133,8 @@ inline void insertIfCurved(MElement *el, std::list<MElement*> &bndEl)
const int dim = el->getDim();
// Compute unit normal to straight edge/face
SVector3 normal = (dim == 1) ? el->getEdge(0).normal() : el->getFace(0).normal();
SVector3 normal = (dim == 1) ? el->getEdge(0).normal() :
el->getFace(0).normal();
// Get functional spaces
const nodalBasis *lagBasis = el->getFunctionSpace();
......@@ -412,6 +163,102 @@ inline void insertIfCurved(MElement *el, std::list<MElement*> &bndEl)
}
void getOppositeEdge(MElement *el, const MEdge &elBaseEd, MEdge &elTopEd,
double &edLenMin, double &edLenMax)
{
static const double BIG = 1e300;
const int iElBaseEd = getElementEdge(elBaseEd, el);
edLenMin = BIG;
edLenMax = -BIG;
MEdge elMaxEd;
// Look for largest edge except base one
for (int iElEd = 0; iElEd < el->getNumEdges(); iElEd++) {
if (iElEd != iElBaseEd) {
MEdge edTest = el->getEdge(iElEd);
double edLenTest = edTest.length();
if (edLenTest < edLenMin) edLenMin = edLenTest;
if (edLenTest > edLenMax) {
edLenMax = edLenTest;
elMaxEd = edTest;
}
}
}
// Build top edge from max edge (with right orientation)
bool sameOrient = false;
if (el->getType() == TYPE_TRI) { // Triangle: look for common vertex
sameOrient = (elBaseEd.getVertex(0) == elMaxEd.getVertex(0));
}
else { // Quad: look for edge between base and top vert.
MEdge sideEdTest(elBaseEd.getVertex(0), elMaxEd.getVertex(0));
for (int iElEd = 0; iElEd < el->getNumEdges(); iElEd++)
if (el->getEdge(iElEd) == sideEdTest) {
sameOrient = true;
break;
}
}
if (sameOrient) elTopEd = elMaxEd;
else elTopEd = MEdge(elMaxEd.getVertex(1), elMaxEd.getVertex(0));
}
bool getColumn2D(MEdgeVecMEltMap &ed2el,
int maxNumLayers, const MEdge &baseEd,
std::vector<MVertex*> &baseVert,
std::vector<MVertex*> &topPrimVert,
std::vector<MElement*> &blob)
{
static const double tol = 2.;
// Get first element and base vertices
std::vector<MElement*> firstElts = ed2el[baseEd];
if (firstElts.size() != 1) {
Msg::Error("%i edge(s) found for edge (%i, %i)", firstElts.size(),
baseEd.getVertex(0)->getNum(), baseEd.getVertex(1)->getNum());
return false;
}
MElement *el = firstElts[0];
const int iFirstElEd = getElementEdge(baseEd, el);
el->getEdgeVertices(iFirstElEd, baseVert);
MEdge elBaseEd(baseVert[0], baseVert[1]);
MEdge elTopEd;
// Sweep column upwards by choosing largest edges in each element
for (int iLayer = 0; iLayer < maxNumLayers; iLayer++) {
double edLenMin, edLenMax;
getOppositeEdge(el, elBaseEd, elTopEd, edLenMin, edLenMax);
if (edLenMax < edLenMin*tol) break;
blob.push_back(el);
std::vector<MElement*> newElts = ed2el[elTopEd];
el = (newElts[0] == el) ? newElts[1] : newElts[0];
elBaseEd = elTopEd;
}
// TODO: improve by taking all vertices?
topPrimVert.resize(2);
topPrimVert[0] = elBaseEd.getVertex(0);
topPrimVert[1] = elBaseEd.getVertex(1);
return true;
}
// TODO: Implement 3D
bool getColumn3D(MFaceVecMEltMap &face2el,
int maxNumLayers, const MFace &baseFace,
std::vector<MVertex*> &baseVert,
std::vector<MVertex*> &topPrimVert,
std::vector<MElement*> &blob)
{
return true;
}
class DbgOutput {
public:
void addMetaEl(MetaEl &mEl) {
......@@ -455,9 +302,7 @@ private:
// Curve elements adjacent to a boundary model entity
void curveMeshFromBnd(std::map<MVertex*, std::vector<MElement *> > &vertex2elements,
const std::map<MVertex*, SVector3> &normVert,
void curveMeshFromBnd(MEdgeVecMEltMap &ed2el, MFaceVecMEltMap &face2el,
GEntity *bndEnt, FastCurvingParameters &p)
{
// Build list of bnd. elements to consider
......@@ -475,7 +320,8 @@ void curveMeshFromBnd(std::map<MVertex*, std::vector<MElement *> > &vertex2eleme
insertIfCurved(gFace->quadrangles[i], bndEl);
}
else
Msg::Error("Cannot treat model entity %i of dim %i", bndEnt->tag(), bndEnt->dim());
Msg::Error("Cannot treat model entity %i of dim %i", bndEnt->tag(),
bndEnt->dim());
DbgOutput dbgOut;
......@@ -484,19 +330,16 @@ void curveMeshFromBnd(std::map<MVertex*, std::vector<MElement *> > &vertex2eleme
itBE != bndEl.end(); itBE++) { // Loop over bnd. elements
const int bndType = (*itBE)->getType();
int metaElType;
bool foundVert;
MElement *firstEl = 0;
bool foundCol;
std::vector<MVertex*> baseVert, topPrimVert;
std::vector<MElement*> blob;
if (bndType == TYPE_LIN) { // 1D boundary?
MVertex *vb0 = (*itBE)->getVertex(0);
MVertex *vb1 = (*itBE)->getVertex(1);
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, baseEd,
p.maxNumLayers, baseVert, topPrimVert);
foundCol = getColumn2D(ed2el, p.maxNumLayers, baseEd, baseVert,
topPrimVert, blob);
}
else { // 2D boundary
MVertex *vb0 = (*itBE)->getVertex(0);
......@@ -512,28 +355,19 @@ void curveMeshFromBnd(std::map<MVertex*, std::vector<MElement *> > &vertex2eleme
metaElType = TYPE_PRI;
}
MFace baseFace(vb0, vb1, vb2, vb3);
firstEl = getFirstEl(vertex2elements, baseFace);
if (!firstEl) {
Msg::Error("Did not find first domain element for boundary element %i",
(*itBE)->getNum());
continue;
}
foundVert = getBaseVertices(baseFace, firstEl, baseVert);
if (!foundVert) continue; // Skip bnd. el. if base vertices not found
foundVert = getTopPrimVertices(vertex2elements, normVert,
baseFace, p.maxNumLayers,
baseVert, topPrimVert);
}
if (!foundVert) continue; // Skip bnd. el. if top vertices not found
int order = firstEl->getPolynomialOrder();
foundCol = getColumn3D(face2el, p.maxNumLayers, baseFace, baseVert,
topPrimVert, blob);
}
if (!foundCol) continue; // Skip bnd. el. if top vertices not found
int order = blob[0]->getPolynomialOrder();
MetaEl metaEl(metaElType, order, baseVert, topPrimVert);
dbgOut.addMetaEl(metaEl);
std::set<MElement*> blob;
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);
for (int iEl = 0; iEl < blob.size(); iEl++) {
MElement *&elt = blob[iEl];
makeStraight(elt, movedVert);
for (int iV = elt->getNumPrimaryVertices();
iV < elt->getNumVertices(); ++iV) { // Loop over HO vert. of each el. in blob
MVertex* vert = elt->getVertex(iV);
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 (metaEl.straightToCurved(xyzS,xyzC)) {
......@@ -565,9 +399,11 @@ void HighOrderMeshFastCurving(GModel *gm, FastCurvingParameters &p)
// Compute vert. -> elt. connectivity
Msg::Info("Computing connectivity...");
std::map<MVertex*, std::vector<MElement *> > vertex2elements;
MEdgeVecMEltMap ed2el;
MFaceVecMEltMap face2el;
for (int iEnt = 0; iEnt < allEntities.size(); ++iEnt)
calcVertex2Elements(p.dim, allEntities[iEnt], vertex2elements);
if (p.dim == 2) calcEdge2Elements(allEntities[iEnt], ed2el);
else calcFace2Elements(allEntities[iEnt], face2el);
// Build multimap of each geometric entity to its boundaries
std::multimap<GEntity*,GEntity*> entities;
......@@ -578,17 +414,12 @@ void HighOrderMeshFastCurving(GModel *gm, FastCurvingParameters &p)
entities.insert(std::pair<GEntity*,GEntity*>(dummy, entity));
}
// Build normals to base vertices
std::map<GEntity*, std::map<MVertex*, SVector3> > normVertEnt; // Normal to each vertex for each geom. entity
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;
Msg::Info("Curving elements for boundary entity %d...", bndEnt->tag());
std::map<MVertex*, SVector3> &normVert = normVertEnt[bndEnt];
curveMeshFromBnd(vertex2elements, normVert, bndEnt, p);
curveMeshFromBnd(ed2el, face2el, bndEnt, p);
}
double t2 = Cpu();
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment