diff --git a/contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp b/contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp
index 0c3ab4cea37fd8edc436c811e41f6406f415a957..4d057c9d00173db394bf4e2ddc1080eadff73a41 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp
@@ -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)
+// Get the face index of an element given a MFace
+inline int getElementFace(const MFace &face, MElement *el)
{
- 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);
+ for (int iFace = 0; iFace < el->getNumFaces(); iFace++)
+ if (el->getFace(iFace) == face) return iFace;
+ return -1;
}
-// Get first domain element for a given boundary face
-MElement *getFirstEl(const std::map<MVertex*, std::vector<MElement*> > &vertex2elements,
- const MFace &baseFace)
+// Get the edge index of an element given a MEdge
+inline int getElementEdge(const MEdge &edge, 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;
-}
-
-
-
-// Get base vertices (in order of first element) for an edge
-bool getBaseVertices(const MEdge &baseEd, MElement *firstEl, std::vector<MVertex*> &baseVert)
-{
- 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);
+ foundCol = getColumn3D(face2el, p.maxNumLayers, baseFace, baseVert,
+ topPrimVert, blob);
}
- if (!foundVert) continue; // Skip bnd. el. if top vertices not found
- int order = firstEl->getPolynomialOrder();
+ 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();