diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index a7dd1906ce8793f24ade025e12c9a6537c78c36c..7a0eab8af78e71ac2bf2020e8a385854278b21a2 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -1135,6 +1135,9 @@ StringXNumber MeshOptions_Number[] = {
{ F|O, "MshFilePartitioned" , opt_mesh_msh_file_partitioned , 0. ,
"Split MSH file by mesh partition" },
+ { F|O, "PartitionByExtrusion" , opt_mesh_partition_by_extrusion, 0. ,
+ "Special partitioner that annotates all all extruded elements to the same node as the source element" },
+
{ F, "NbHexahedra" , opt_mesh_nb_hexahedra , 0. ,
"Number of hexahedra in the current mesh (read-only)" },
{ F, "NbNodes" , opt_mesh_nb_nodes , 0. ,
diff --git a/Common/Options.cpp b/Common/Options.cpp
index f1077ee93827502e99760f2a756bdab6fa607ec3..a94e6f4856569800df16e3b803c94c048ec71d8d 100644
--- a/Common/Options.cpp
+++ b/Common/Options.cpp
@@ -5811,6 +5811,14 @@ double opt_mesh_partition_partitioner(OPT_ARGS_NUM)
}
return CTX::instance()->partitionOptions.partitioner;
}
+double opt_mesh_partition_by_extrusion(OPT_ARGS_NUM)
+{
+ if(action & GMSH_SET) {
+ const int ival = (int)val;
+ CTX::instance()->partitionOptions.partitionByExtrusion = ival;
+ }
+ return CTX::instance()->partitionOptions.partitionByExtrusion;
+}
double opt_mesh_partition_num(OPT_ARGS_NUM)
{
diff --git a/Common/Options.h b/Common/Options.h
index 5a7cbc189f33cfd2f31ad040147b2082816548ed..58eccd97a9a4a3b6247c2fb6e57c87caaed45107 100644
--- a/Common/Options.h
+++ b/Common/Options.h
@@ -563,6 +563,7 @@ double opt_mesh_partition_chaco_terminal_propogation(OPT_ARGS_NUM);
double opt_mesh_partition_metis_algorithm(OPT_ARGS_NUM);
double opt_mesh_partition_metis_edge_matching(OPT_ARGS_NUM);
double opt_mesh_partition_metis_refine_algorithm(OPT_ARGS_NUM);
+double opt_mesh_partition_by_extrusion(OPT_ARGS_NUM);
double opt_mesh_clip(OPT_ARGS_NUM);
double opt_solver_listen(OPT_ARGS_NUM);
double opt_solver_plugins(OPT_ARGS_NUM);
diff --git a/Geo/ExtrudeParams.cpp b/Geo/ExtrudeParams.cpp
index 8841de0dd7b4a9e5ff762b783ea32d401024bf8b..ea12b7d9990d0de2bf7e6e489aa4a79b46cea1e2 100644
--- a/Geo/ExtrudeParams.cpp
+++ b/Geo/ExtrudeParams.cpp
@@ -21,7 +21,7 @@ static void Projette(double p[3], double mat[3][3])
p[2] = Z;
}
-ExtrudeParams::ExtrudeParams(int ModeEx)
+ExtrudeParams::ExtrudeParams(int ModeEx) : elementMap(this)
{
geo.Mode = ModeEx;
geo.Source = -1;
@@ -121,3 +121,76 @@ double ExtrudeParams::u(int iLayer, int iElemLayer)
double t = (double)iElemLayer / (double)mesh.NbElmLayer[iLayer];
return t0 + t * (t1 - t0);
}
+
+ExtrudeParams::ExtrusionElementMap::ExtrusionElementMap(ExtrudeParams*const parent)
+{
+ _parent = parent;
+}
+
+std::vector< MElement* >* ExtrudeParams::ExtrusionElementMap::getExtrudedElems(MElement* source)
+{
+ std::map<MElement*,std::vector<MElement*> >::iterator it = _extrudedElements.find(source);
+ if(it != _extrudedElements.end())
+ return &(it->second);
+ return NULL;
+}
+void ExtrudeParams::ExtrusionElementMap::clear()
+{
+ _extrudedElements.clear();
+}
+
+bool ExtrudeParams::ExtrusionElementMap::empty()
+{
+ return _extrudedElements.empty();
+}
+
+void ExtrudeParams::ExtrusionElementMap::addExtrudedElem(MElement* source, MElement* extrudedElem)
+{
+ std::map<MElement*,std::vector<MElement*> >::iterator it = _extrudedElements.find(source);
+ if(it != _extrudedElements.end())
+ it->second.push_back(extrudedElem);
+ else {
+ std::vector<MElement*>* vec = new std::vector<MElement*>();
+ int totalNbElems = 0;
+ for (int i = 0;i<_parent->mesh.NbLayer;i++)
+ totalNbElems += _parent->mesh.NbElmLayer[i];
+ vec->reserve(totalNbElems);
+ vec->push_back(extrudedElem);
+ _extrudedElements.insert(std::pair<MElement*,std::vector<MElement*> >(source,*vec));
+ }
+ SPoint3 np = extrudedElem->barycenter(), sp = source->barycenter();
+// if(extrudedElem->getDim() == 3)
+// printf("Inst (%.3f,%.3f,%.3f) -> (%.3f,%.3f,%.3f) %p->%p\n",sp.x(),sp.y(),sp.z(),np.x(),np.y(),np.z(),source,extrudedElem);
+}
+
+// Propagates the partition information from the source elements to the extruded elements.
+// Increments the partition sizes if partitionSizes is given as argument.
+// For efficient looping, this routine is within ExtrusionElementMap class.
+void ExtrudeParams::ExtrusionElementMap::propagatePartitionInformation(std::vector< int >* partitionSizes)
+{
+ if (_extrudedElements.empty())
+ Msg::Error("No extrusion information found!");
+ std::map<MElement*,std::vector<MElement*> >::iterator columnit;
+ for(columnit=_extrudedElements.begin(); columnit != _extrudedElements.end(); columnit++) {
+ MElement* sourceElem = (*columnit).first;
+ if(!sourceElem) {
+ Msg::Warning("No source found!");
+ continue;
+ }
+ std::vector<MElement*> extrudedElements = (*columnit).second;
+// if(!extrudedElements) {
+// // Msg::Warning("No element vector found!");
+// continue;
+// }
+ for(int iE = 0;iE < extrudedElements.size();iE++) {
+ if(extrudedElements[iE]) {
+ extrudedElements[iE]->setPartition(sourceElem->getPartition());
+ if (partitionSizes)
+ ++(*partitionSizes)[sourceElem->getPartition()-1];
+// SPoint3 np = extrudedElements[iE]->barycenter(), sp = sourceElem->barycenter();
+// if(extrudedElements[iE]->getDim() == 3)
+// printf("Read (%.3f,%.3f,%.3f) %d -> (%.3f,%.3f,%.3f) %d %p->%p\n",sp.x(),sp.y(),sp.z(),extrudedElements[iE]->getPartition() ,np.x(),np.y(),np.z(),sourceElem->getPartition(),sourceElem,extrudedElements[iE]);
+ }
+ }
+ }
+}
diff --git a/Geo/ExtrudeParams.h b/Geo/ExtrudeParams.h
index 532b99c84b4a3025e440bd3925f1192ed85d00d0..ce6c344f0437f97121a0a4d27bde4f46f1cd2cff 100644
--- a/Geo/ExtrudeParams.h
+++ b/Geo/ExtrudeParams.h
@@ -10,6 +10,7 @@
#include <map>
#include <string>
#include "SmoothData.h"
+#include "MElement.h"
// geo.Mode
#define EXTRUDED_ENTITY 1
@@ -23,6 +24,19 @@
class ExtrudeParams{
public :
+ class ExtrusionElementMap {
+ private:
+ ExtrudeParams* _parent;
+ // maps source element to all extruded elements
+ std::map<MElement*,std::vector<MElement*> > _extrudedElements;
+ public:
+ ExtrusionElementMap(ExtrudeParams* const parent);
+ std::vector<MElement*>* getExtrudedElems(MElement* source);
+ void addExtrudedElem(MElement* source,MElement* extrudedElem);
+ void clear();
+ bool empty();
+ void propagatePartitionInformation(std::vector<int>* partitionSizes = NULL);
+ } elementMap;
static smooth_data *normals;
ExtrudeParams(int Mode = EXTRUDED_ENTITY);
void fill(int type,
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 421d93724aa22dfef1115f37d0f9383c417d51a4..3585a71a1e34d3384ce7d4997669fcfbf421187a 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -315,7 +315,7 @@ void MElement::xyz2uvw(double xyz[3], double uvw[3])
uvw[0] = uvw[1] = uvw[2] = 0.;
int iter = 1, maxiter = 20;
double error = 1., tol = 1.e-6;
-
+
while (error > tol && iter < maxiter){
double jac[3][3];
if(!getJacobian(uvw[0], uvw[1], uvw[2], jac)) break;
diff --git a/Geo/MPrism.cpp b/Geo/MPrism.cpp
index 6039f44b8dbd52881c1ccca093c4c7321fa3574b..29de1872482827a7e120ed038081956818cfa0ae 100644
--- a/Geo/MPrism.cpp
+++ b/Geo/MPrism.cpp
@@ -53,6 +53,20 @@ const polynomialBasis* MPrism::getFunctionSpace(int o) const
return 0;
}
+double MPrism::getInnerRadius()
+{
+ double dist[3], k = 0.;
+ double triEdges[3] = {0,1,3};
+ for (int i = 0; i < 3; i++){
+ MEdge e = getEdge(triEdges[i]);
+ dist[i] = e.getVertex(0)->distance(e.getVertex(1));
+ k += 0.5 * dist[i];
+ }
+ double radTri = sqrt(k * (k - dist[0]) * (k - dist[1]) * (k - dist[2])) / k;
+ double radVert = 0.5*getVertex(0)->distance(getVertex(3));
+ return std::min(radTri,radVert);
+}
+
void MPrism::getFaceInfo(const MFace &face, int &ithFace, int &sign, int &rot) const
{
for (ithFace = 0; ithFace < 5; ithFace++){
diff --git a/Geo/MPrism.h b/Geo/MPrism.h
index 68713e0920897958be6732cb0ebe32e4d54cd7b8..1efd463a12dc71f453883f202d1b48e663543aef 100644
--- a/Geo/MPrism.h
+++ b/Geo/MPrism.h
@@ -62,6 +62,7 @@ class MPrism : public MElement {
~MPrism(){}
virtual int getDim() const { return 3; }
virtual int getNumVertices() const { return 6; }
+ virtual double getInnerRadius();
virtual MVertex *getVertex(int num){ return _v[num]; }
virtual MVertex *getVertexMED(int num)
{
diff --git a/Mesh/meshGEdgeExtruded.cpp b/Mesh/meshGEdgeExtruded.cpp
index 7bac2545a2852fd5ce00e3d2f120d9c2d7d23873..c06de9aaf21a90deca6f06e69c8fae9cdd43b25d 100644
--- a/Mesh/meshGEdgeExtruded.cpp
+++ b/Mesh/meshGEdgeExtruded.cpp
@@ -54,13 +54,13 @@ int MeshExtrudedCurve(GEdge *ge)
if(!ep || !ep->mesh.ExtrudeMesh)
return 0;
+ GEdge *from = ge->model()->getEdgeByTag(std::abs(ep->geo.Source));
if(ep->geo.Mode == EXTRUDED_ENTITY) {
// curve is extruded from a point
extrudeMesh(ge->getBeginVertex(), ge);
}
else {
// curve is a copy of another curve (the "top" of the extrusion)
- GEdge *from = ge->model()->getEdgeByTag(std::abs(ep->geo.Source));
if(!from){
Msg::Error("Unknown source curve %d for extrusion", ep->geo.Source);
return 0;
@@ -74,8 +74,14 @@ int MeshExtrudedCurve(GEdge *ge)
ge->getBeginVertex()->mesh_vertices[0] : ge->mesh_vertices[i - 1];
MVertex *v1 = (i == ge->mesh_vertices.size()) ?
ge->getEndVertex()->mesh_vertices[0] : ge->mesh_vertices[i];
- ge->lines.push_back(new MLine(v0, v1));
+ MLine* newElem = new MLine(v0, v1);
+ ge->lines.push_back(newElem);
+ if(ep->geo.Mode == COPIED_ENTITY) {
+ // Extrusion information is only stored for copied edges
+ // (the source of extruded edge is a vertex, not an element)
+ MElement* sourceElem = from->getMeshElement(i);
+ ep->elementMap.addExtrudedElem(sourceElem,newElem);
+ }
}
-
return 1;
}
diff --git a/Mesh/meshGFaceExtruded.cpp b/Mesh/meshGFaceExtruded.cpp
index 780b408e4d2bbe38745db8779eec8c38f498f7eb..8e7e254c0e2072537e4309b8a647fc4254f9ca22 100644
--- a/Mesh/meshGFaceExtruded.cpp
+++ b/Mesh/meshGFaceExtruded.cpp
@@ -12,34 +12,45 @@
#include "Context.h"
#include "GmshMessage.h"
+static void addTriangle(MVertex* v1, MVertex* v2, MVertex* v3, GFace *to, MElement* source) {
+ MTriangle* newTri = new MTriangle(v1, v2, v3);
+ to->triangles.push_back(newTri);
+ to->meshAttributes.extrude->elementMap.addExtrudedElem((MElement*)source,(MElement*)newTri);
+}
+
+static void addQuadrangle(MVertex* v1, MVertex* v2, MVertex* v3, MVertex* v4, GFace *to, MElement* source) {
+ MQuadrangle* newQuad = new MQuadrangle(v1, v2, v3, v4);
+ to->quadrangles.push_back(newQuad);
+ to->meshAttributes.extrude->elementMap.addExtrudedElem((MElement*)source,(MElement*)newQuad);
+}
+
static void createQuaTri(std::vector<MVertex*> &v, GFace *to,
- std::set<std::pair<MVertex*, MVertex*> > *constrainedEdges)
+ std::set<std::pair<MVertex*, MVertex*> > *constrainedEdges,MLine* source)
{
ExtrudeParams *ep = to->meshAttributes.extrude;
-
if(v[0] == v[1] || v[1] == v[3])
- to->triangles.push_back(new MTriangle(v[0], v[3], v[2]));
+ addTriangle(v[0], v[3], v[2],to,source);
else if(v[0] == v[2] || v[2] == v[3])
- to->triangles.push_back(new MTriangle(v[0], v[1], v[3]));
+ addTriangle(v[0], v[1], v[3],to,source);
else if(v[0] == v[3] || v[1] == v[2])
Msg::Error("Uncoherent extruded quadrangle in surface %d", to->tag());
else{
if(ep->mesh.Recombine){
- to->quadrangles.push_back(new MQuadrangle(v[0], v[1], v[3], v[2]));
+ addQuadrangle(v[0], v[1], v[3], v[2],to,source);
}
else if(!constrainedEdges){
- to->triangles.push_back(new MTriangle(v[0], v[1], v[3]));
- to->triangles.push_back(new MTriangle(v[0], v[3], v[2]));
+ addTriangle(v[0], v[1], v[3],to,source);
+ addTriangle(v[0], v[3], v[2],to,source);
}
else{
std::pair<MVertex*, MVertex*> p(std::min(v[1], v[2]), std::max(v[1], v[2]));
if(constrainedEdges->count(p)){
- to->triangles.push_back(new MTriangle(v[2], v[1], v[0]));
- to->triangles.push_back(new MTriangle(v[2], v[3], v[1]));
+ addTriangle(v[2], v[1], v[0],to,source);
+ addTriangle(v[2], v[3], v[1],to,source);
}
else{
- to->triangles.push_back(new MTriangle(v[2], v[3], v[0]));
- to->triangles.push_back(new MTriangle(v[0], v[3], v[1]));
+ addTriangle(v[2], v[3], v[0],to,source);
+ addTriangle(v[0], v[3], v[1],to,source);
}
}
}
@@ -100,7 +111,7 @@ static void extrudeMesh(GEdge *from, GFace *to,
}
verts.push_back(*itp);
}
- createQuaTri(verts, to, constrainedEdges);
+ createQuaTri(verts, to, constrainedEdges,from->lines[i]);
}
}
}
@@ -143,7 +154,7 @@ static void copyMesh(GFace *from, GFace *to,
}
verts.push_back(*itp);
}
- to->triangles.push_back(new MTriangle(verts));
+ addTriangle(verts[0],verts[1],verts[2],to,from->triangles[i]);
}
for(unsigned int i = 0; i < from->quadrangles.size(); i++){
std::vector<MVertex*> verts;
@@ -164,7 +175,7 @@ static void copyMesh(GFace *from, GFace *to,
}
verts.push_back(*itp);
}
- to->quadrangles.push_back(new MQuadrangle(verts));
+ addQuadrangle(verts[0],verts[1],verts[2],verts[3],to,from->quadrangles[i]);
}
}
diff --git a/Mesh/meshGRegionExtruded.cpp b/Mesh/meshGRegionExtruded.cpp
index 3aaea62fd7d9b43b8e756f336f5eb4085a0d51e0..72dcff796ede15fc1b4facec63b912f29f2283e5 100644
--- a/Mesh/meshGRegionExtruded.cpp
+++ b/Mesh/meshGRegionExtruded.cpp
@@ -17,7 +17,30 @@
#include "Context.h"
#include "GmshMessage.h"
-static void createPriPyrTet(std::vector<MVertex*> &v, GRegion *to)
+static void addTetrahedron(MVertex* v1, MVertex* v2, MVertex* v3, MVertex* v4, GRegion *to, MElement* source) {
+ MTetrahedron* newElem = new MTetrahedron(v1, v2, v3, v4);
+ to->tetrahedra.push_back(newElem);
+ to->meshAttributes.extrude->elementMap.addExtrudedElem((MElement*)source,(MElement*)newElem);
+}
+
+static void addPyramid(MVertex* v1, MVertex* v2, MVertex* v3, MVertex* v4, MVertex* v5, GRegion *to, MElement* source) {
+ MPyramid* newElem = new MPyramid(v1, v2, v3, v4, v5);
+ to->pyramids.push_back(newElem);
+ to->meshAttributes.extrude->elementMap.addExtrudedElem((MElement*)source,(MElement*)newElem);
+}
+
+static void addPrism(MVertex* v1, MVertex* v2, MVertex* v3, MVertex* v4, MVertex* v5, MVertex* v6, GRegion *to, MElement* source) {
+ MPrism* newElem = new MPrism(v1, v2, v3, v4, v5, v6);
+ to->prisms.push_back(newElem);
+ to->meshAttributes.extrude->elementMap.addExtrudedElem((MElement*)source,(MElement*)newElem);
+}
+static void addHexahedron(MVertex* v1, MVertex* v2, MVertex* v3, MVertex* v4, MVertex* v5, MVertex* v6, MVertex* v7, MVertex* v8, GRegion *to, MElement* source) {
+ MHexahedron* newElem = new MHexahedron(v1, v2, v3, v4, v5, v6, v7, v8);
+ to->hexahedra.push_back(newElem);
+ to->meshAttributes.extrude->elementMap.addExtrudedElem((MElement*)source,(MElement*)newElem);
+}
+
+static void createPriPyrTet(std::vector<MVertex*> &v, GRegion *to, MElement* source)
{
int dup[3];
int j = 0;
@@ -27,27 +50,27 @@ static void createPriPyrTet(std::vector<MVertex*> &v, GRegion *to)
if(j == 2) {
if(dup[0] == 0 && dup[1] == 1)
- to->tetrahedra.push_back(new MTetrahedron(v[0], v[1], v[2], v[5]));
+ addTetrahedron(v[0], v[1], v[2], v[5],to,source);
else if(dup[0] == 1 && dup[1] == 2)
- to->tetrahedra.push_back(new MTetrahedron(v[0], v[1], v[2], v[3]));
+ addTetrahedron(v[0], v[1], v[2], v[3],to,source);
else
- to->tetrahedra.push_back(new MTetrahedron(v[0], v[1], v[2], v[4]));
+ addTetrahedron(v[0], v[1], v[2], v[4],to,source);
}
else if(j == 1) {
if(dup[0] == 0)
- to->pyramids.push_back(new MPyramid(v[1], v[4], v[5], v[2], v[0]));
+ addPyramid(v[1], v[4], v[5], v[2], v[0],to,source);
else if(dup[0] == 1)
- to->pyramids.push_back(new MPyramid(v[0], v[2], v[5], v[3], v[1]));
+ addPyramid(v[0], v[2], v[5], v[3], v[1],to,source);
else
- to->pyramids.push_back(new MPyramid(v[0], v[1], v[4], v[3], v[2]));
+ addPyramid(v[0], v[1], v[4], v[3], v[2],to,source);
}
else {
- to->prisms.push_back(new MPrism(v));
+ addPrism(v[0], v[1], v[2], v[3], v[4], v[5],to,source);
if(j) Msg::Error("Degenerated prism in extrusion of volume %d", to->tag());
}
}
-static void createHexPri(std::vector<MVertex*> &v, GRegion *to)
+static void createHexPri(std::vector<MVertex*> &v, GRegion *to, MElement* source)
{
int dup[4];
int j = 0;
@@ -57,26 +80,26 @@ static void createHexPri(std::vector<MVertex*> &v, GRegion *to)
if(j == 2) {
if(dup[0] == 0 && dup[1] == 1)
- to->prisms.push_back(new MPrism(v[0], v[3], v[7], v[1], v[2], v[6]));
+ addPrism(v[0], v[3], v[7], v[1], v[2], v[6],to,source);
else if(dup[0] == 1 && dup[1] == 2)
- to->prisms.push_back(new MPrism(v[0], v[1], v[4], v[3], v[2], v[7]));
+ addPrism(v[0], v[1], v[4], v[3], v[2], v[7],to,source);
else if(dup[0] == 2 && dup[1] == 3)
- to->prisms.push_back(new MPrism(v[0], v[3], v[4], v[1], v[2], v[5]));
+ addPrism(v[0], v[3], v[4], v[1], v[2], v[5],to,source);
else if(dup[0] == 0 && dup[1] == 3)
- to->prisms.push_back(new MPrism(v[0], v[1], v[5], v[3], v[2], v[6]));
+ addPrism(v[0], v[1], v[5], v[3], v[2], v[6],to,source);
else
Msg::Error("Uncoherent hexahedron in extrusion of volume %d", to->tag());
}
else {
- to->hexahedra.push_back(new MHexahedron(v));
+ addHexahedron(v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7],to,source);
if(j) Msg::Error("Degenerated hexahedron in extrusion of volume %d", to->tag());
}
}
-static void createTet(MVertex *v1, MVertex *v2, MVertex *v3, MVertex *v4, GRegion *to)
+static void createTet(MVertex *v1, MVertex *v2, MVertex *v3, MVertex *v4, GRegion *to, MElement* source)
{
if(v1 != v2 && v1 != v3 && v1 != v4 && v2 != v3 && v2 != v4 && v3 != v4)
- to->tetrahedra.push_back(new MTetrahedron(v1, v2, v3, v4));
+ addTetrahedron(v1, v2, v3, v4,to,source);
}
static int getExtrudedVertices(MElement *ele, ExtrudeParams *ep, int j, int k,
@@ -140,8 +163,9 @@ static void extrudeMesh(GFace *from, GRegion *to,
for(int j = 0; j < ep->mesh.NbLayer; j++) {
for(int k = 0; k < ep->mesh.NbElmLayer[j]; k++) {
std::vector<MVertex*> verts;
- if(getExtrudedVertices(from->triangles[i], ep, j, k, pos, verts) == 6)
- createPriPyrTet(verts, to);
+ if(getExtrudedVertices(from->triangles[i], ep, j, k, pos, verts) == 6) {
+ createPriPyrTet(verts, to,from->triangles[i]);
+ }
}
}
}
@@ -155,7 +179,7 @@ static void extrudeMesh(GFace *from, GRegion *to,
for(int k = 0; k < ep->mesh.NbElmLayer[j]; k++) {
std::vector<MVertex*> verts;
if(getExtrudedVertices(from->quadrangles[i], ep, j, k, pos, verts) == 8)
- createHexPri(verts, to);
+ createHexPri(verts, to, from->quadrangles[i]);
}
}
}
@@ -360,51 +384,52 @@ static void phase3(GRegion *gr,
if(!from) return;
for(unsigned int i = 0; i < from->triangles.size(); i++){
+ MTriangle* tri = from->triangles[i];
for(int j = 0; j < ep->mesh.NbLayer; j++) {
for(int k = 0; k < ep->mesh.NbElmLayer[j]; k++) {
std::vector<MVertex*> v;
- if(getExtrudedVertices(from->triangles[i], ep, j, k, pos, v) == 6){
+ if(getExtrudedVertices(tri, ep, j, k, pos, v) == 6){
if(edgeExists(v[3], v[1], edges) &&
edgeExists(v[4], v[2], edges) &&
edgeExists(v[3], v[2], edges)) {
- createTet(v[0], v[1], v[2], v[3], gr);
- createTet(v[3], v[4], v[5], v[2], gr);
- createTet(v[1], v[3], v[4], v[2], gr);
+ createTet(v[0], v[1], v[2], v[3], gr, tri);
+ createTet(v[3], v[4], v[5], v[2], gr, tri);
+ createTet(v[1], v[3], v[4], v[2], gr, tri);
}
else if(edgeExists(v[3], v[1], edges) &&
edgeExists(v[1], v[5], edges) &&
edgeExists(v[3], v[2], edges)) {
- createTet(v[0], v[1], v[2], v[3], gr);
- createTet(v[3], v[4], v[5], v[1], gr);
- createTet(v[3], v[1], v[5], v[2], gr);
+ createTet(v[0], v[1], v[2], v[3], gr, tri);
+ createTet(v[3], v[4], v[5], v[1], gr, tri);
+ createTet(v[3], v[1], v[5], v[2], gr, tri);
}
else if(edgeExists(v[3], v[1], edges) &&
edgeExists(v[1], v[5], edges) &&
edgeExists(v[5], v[0], edges)) {
- createTet(v[0], v[1], v[2], v[5], gr);
- createTet(v[3], v[4], v[5], v[1], gr);
- createTet(v[1], v[3], v[5], v[0], gr);
+ createTet(v[0], v[1], v[2], v[5], gr, tri);
+ createTet(v[3], v[4], v[5], v[1], gr, tri);
+ createTet(v[1], v[3], v[5], v[0], gr, tri);
}
else if(edgeExists(v[4], v[0], edges) &&
edgeExists(v[4], v[2], edges) &&
edgeExists(v[3], v[2], edges)) {
- createTet(v[0], v[1], v[2], v[4], gr);
- createTet(v[3], v[4], v[5], v[2], gr);
- createTet(v[0], v[3], v[4], v[2], gr);
+ createTet(v[0], v[1], v[2], v[4], gr, tri);
+ createTet(v[3], v[4], v[5], v[2], gr, tri);
+ createTet(v[0], v[3], v[4], v[2], gr, tri);
}
else if(edgeExists(v[4], v[0], edges) &&
edgeExists(v[4], v[2], edges) &&
edgeExists(v[5], v[0], edges)) {
- createTet(v[0], v[1], v[2], v[4], gr);
- createTet(v[3], v[4], v[5], v[0], gr);
- createTet(v[0], v[2], v[4], v[5], gr);
+ createTet(v[0], v[1], v[2], v[4], gr, tri);
+ createTet(v[3], v[4], v[5], v[0], gr, tri);
+ createTet(v[0], v[2], v[4], v[5], gr, tri);
}
else if(edgeExists(v[4], v[0], edges) &&
edgeExists(v[1], v[5], edges) &&
edgeExists(v[5], v[0], edges)) {
- createTet(v[0], v[1], v[2], v[5], gr);
- createTet(v[3], v[4], v[5], v[0], gr);
- createTet(v[0], v[1], v[4], v[5], gr);
+ createTet(v[0], v[1], v[2], v[5], gr, tri);
+ createTet(v[3], v[4], v[5], v[0], gr, tri);
+ createTet(v[0], v[1], v[4], v[5], gr, tri);
}
}
}
@@ -429,7 +454,6 @@ int SubdivideExtrudedMesh(GModel *m)
}
if(regions.empty()) return 0;
-
Msg::Info("Subdividing extruded mesh");
// create edges on lateral sides of "prisms"
@@ -465,36 +489,59 @@ int SubdivideExtrudedMesh(GModel *m)
for(unsigned int i = 0; i < gr->pyramids.size(); i++)
delete gr->pyramids[i];
gr->pyramids.clear();
+ gr->meshAttributes.extrude->elementMap.clear(); // reconstruct extrusion info
phase3(gr, pos, edges);
- }
- // re-Extrude bounding surfaces using edges as constraint
- std::set<GFace*> faces;
- for(unsigned int i = 0; i < regions.size(); i++){
- std::list<GFace*> f = regions[i]->faces();
- faces.insert(f.begin(), f.end());
- }
- for(std::set<GFace*>::iterator it = faces.begin(); it != faces.end(); it++){
- ExtrudeParams *ep = (*it)->meshAttributes.extrude;
- if(ep && ep->mesh.ExtrudeMesh && ep->geo.Mode == EXTRUDED_ENTITY &&
- !ep->mesh.Recombine){
- GFace *gf = *it;
- Msg::Info("Remeshing surface %d", gf->tag());
- for(unsigned int i = 0; i < gf->triangles.size(); i++)
- delete gf->triangles[i];
- gf->triangles.clear();
- for(unsigned int i = 0; i < gf->quadrangles.size(); i++)
- delete gf->quadrangles[i];
- gf->quadrangles.clear();
- MeshExtrudedSurface(gf, &edges);
+ // re-Extrude bounding surfaces using edges as constraint
+ std::list<GFace*> faces = gr->faces();
+ for(std::list<GFace*>::iterator it = faces.begin(); it != faces.end(); it++){
+ ExtrudeParams *ep = (*it)->meshAttributes.extrude;
+ if(ep && ep->mesh.ExtrudeMesh && ep->geo.Mode == EXTRUDED_ENTITY &&
+ !ep->mesh.Recombine){
+ GFace *gf = *it;
+ Msg::Info("Remeshing surface %d", gf->tag());
+ for(unsigned int i = 0; i < gf->triangles.size(); i++)
+ delete gf->triangles[i];
+ gf->triangles.clear();
+ for(unsigned int i = 0; i < gf->quadrangles.size(); i++)
+ delete gf->quadrangles[i];
+ gf->quadrangles.clear();
+ ep->elementMap.clear(); // reconstruct extrusion info
+ MeshExtrudedSurface(gf, &edges);
+ }
}
}
+
+// std::set<GFace*> faces;
+// for(unsigned int i = 0; i < regions.size(); i++){
+// std::list<GFace*> f = regions[i]->faces();
+// faces.insert(f.begin(), f.end());
+// }
+// for(std::set<GFace*>::iterator it = faces.begin(); it != faces.end(); it++){
+// ExtrudeParams *ep = (*it)->meshAttributes.extrude;
+// if(ep && ep->mesh.ExtrudeMesh && ep->geo.Mode == EXTRUDED_ENTITY &&
+// !ep->mesh.Recombine){
+// GFace *gf = *it;
+// Msg::Info("Remeshing surface %d", gf->tag());
+// for(unsigned int i = 0; i < gf->triangles.size(); i++)
+// delete gf->triangles[i];
+// gf->triangles.clear();
+// for(unsigned int i = 0; i < gf->quadrangles.size(); i++)
+// delete gf->quadrangles[i];
+// gf->quadrangles.clear();
+// ep->elementMap.clear(); // reconstruct extrusion info
+// MeshExtrudedSurface(gf, &edges);
+// }
+// }
+
// carve holes if any
+ // TODO: update extrusion information
for(unsigned int i = 0; i < regions.size(); i++){
GRegion *gr = regions[i];
ExtrudeParams *ep = gr->meshAttributes.extrude;
if(ep->mesh.Holes.size()){
+ ep->elementMap.clear();
std::map<int, std::pair<double, std::vector<int> > >::iterator it;
for(it = ep->mesh.Holes.begin(); it != ep->mesh.Holes.end(); it++)
carveHole(gr, it->first, it->second.first, it->second.second);
diff --git a/Mesh/meshPartition.cpp b/Mesh/meshPartition.cpp
index e50c3ad8981d5992dc3ac331841ef4046dba0f94..5f84c66b7bcb09442dbd90464541aa16307f03df 100644
--- a/Mesh/meshPartition.cpp
+++ b/Mesh/meshPartition.cpp
@@ -28,6 +28,7 @@
#include "discreteFace.h"
#include "GFaceCompound.h"
#include "Context.h"
+#include "ExtrudeParams.h"
//--Prototype for Chaco interface
@@ -84,14 +85,14 @@ struct BoElemGr;
typedef std::vector<BoElemGr> BoElemGrVec;
int MakeGraph(GModel *const model, Graph &graph,
- BoElemGrVec *const boElemGrVec = 0);
+ meshPartitionOptions &options, BoElemGrVec *const boElemGrVec = 0);
template <unsigned DIM, typename EntIter, typename EntIterBE>
void MakeGraphDIM(const EntIter begin, const EntIter end,
const EntIterBE beginBE, const EntIterBE endBE,
Graph &graph, BoElemGrVec *const boElemGrVec);
-
+
/*******************************************************************************
*
* Routine partitionMesh
@@ -120,7 +121,7 @@ int RenumberMesh(GModel *const model, meshPartitionOptions &options, std::vector
BoElemGrVec boElemGrVec;
int ier;
Msg::StatusBar(1, true, "Building graph...");
- ier = MakeGraph(model, graph, &boElemGrVec);
+ ier = MakeGraph(model, graph, options, &boElemGrVec);
Msg::StatusBar(1, true, "Renumbering graph...");
if(!ier) ier = RenumberGraph(graph, options);
if(ier) {
@@ -262,21 +263,74 @@ int PartitionMesh(GModel *const model, meshPartitionOptions &options)
BoElemGrVec boElemGrVec;
int ier;
Msg::StatusBar(1, true, "Building graph...");
- ier = MakeGraph(model, graph, &boElemGrVec);
+ ier = MakeGraph(model, graph, options, &boElemGrVec);
Msg::StatusBar(1, true, "Partitioning graph...");
if(!ier) ier = PartitionGraph(graph, options);
if(ier) {
Msg::StatusBar(1, false, "Mesh");
return 1;
}
-
+
// Count partition sizes and assign partitions to internal elements
std::vector<int> ssize(options.num_partitions, 0);
const int n = graph.getNumVertex();
for(int i = 0; i != n; ++i) {
- ++ssize[graph.partition[i] - 1];
+ if (!options.partitionByExtrusion)
+ ++ssize[graph.partition[i] - 1];
graph.element[i]->setPartition(graph.partition[i]);
}
+ // Assign partitions to boundary elements
+ const int nb = boElemGrVec.size();
+ for(int i = 0; i != nb; ++i) {
+ boElemGrVec[i].elem->setPartition(graph.partition[boElemGrVec[i].grVert]);
+ }
+ // propagate partition information to extruded elements
+ if (options.partitionByExtrusion) {
+ unsigned numElem[5];
+ const int meshDim = model->getNumMeshElements(numElem);
+ switch(meshDim) {
+ case 2:
+ // assing partitions for 2D faces and count number of elems in each partition
+ for (GModel:: fiter fit = model->firstFace(); fit != model->lastFace(); fit++) {
+ ExtrudeParams* epFace = (*fit)->meshAttributes.extrude;
+ GEdge *sourceEdge = model->getEdgeByTag(std::abs(epFace->geo.Source));
+ epFace->elementMap.propagatePartitionInformation(&ssize);
+ // loop over extruded edges to set partitions
+ std::list<GEdge*> regionEdges = (*fit)->edges();
+ for (std::list<GEdge*>::iterator eit = regionEdges.begin(); eit != regionEdges.end(); eit++) {
+ ExtrudeParams* epEdge = (*eit)->meshAttributes.extrude;
+ if((*eit)!=sourceEdge && epEdge && !epEdge->elementMap.empty())
+ epEdge->elementMap.propagatePartitionInformation();
+ }
+ }
+ break;
+ case 3:
+ // assing partitions for 3D volumes and count number of elems in each partition
+ for (GModel:: riter rit = model->firstRegion(); rit != model->lastRegion(); rit++) {
+ ExtrudeParams* epRegion = (*rit)->meshAttributes.extrude;
+ GFace *sourceFace = model->getFaceByTag(std::abs(epRegion->geo.Source));
+ epRegion->elementMap.propagatePartitionInformation(&ssize);
+ // loop over extruded faces to set partitions
+ std::list<GFace*> regionFaces = (*rit)->faces();
+ for (std::list<GFace*>::iterator fit = regionFaces.begin(); fit != regionFaces.end(); fit++) {
+ if ((*fit) == sourceFace)
+ continue;
+ ExtrudeParams* epFace = (*fit)->meshAttributes.extrude;
+ epFace->elementMap.propagatePartitionInformation();
+ // loop over bounding edges to set partitions
+ std::list<GEdge*> edges = (*fit)->edges();
+ GEdge *sourceEdge = model->getEdgeByTag(std::abs(epFace->geo.Source));
+ for (std::list<GEdge*>::iterator eit = edges.begin(); eit != edges.end(); eit++) {
+ ExtrudeParams* epEdge = (*eit)->meshAttributes.extrude;
+ if((*eit)!=sourceEdge && epEdge && epEdge->geo.Mode == COPIED_ENTITY)
+ // extruded edges have no source elements (extruded from vertex)
+ epEdge->elementMap.propagatePartitionInformation();
+ }
+ }
+ } // TODO: generalize to 2D too
+ break;
+ }
+ }
int sMin = graph.getNumVertex();
int sMax = 0;
for(int i = 0; i != options.num_partitions; ++i) {
@@ -286,12 +340,6 @@ int PartitionMesh(GModel *const model, meshPartitionOptions &options)
model->setMinPartitionSize(sMin);
model->setMaxPartitionSize(sMax);
- // Assign partitions to boundary elements
- const int nb = boElemGrVec.size();
- for(int i = 0; i != nb; ++i) {
- boElemGrVec[i].elem->setPartition(graph.partition[boElemGrVec[i].grVert]);
- }
-
model->recomputeMeshPartitions();
if (options.createPartitionBoundaries || options.createGhostCells)
@@ -500,7 +548,7 @@ int RenumberGraph(Graph &graph, meshPartitionOptions &options)
*
******************************************************************************/
-int MakeGraph(GModel *const model, Graph &graph, BoElemGrVec *const boElemGrVec)
+int MakeGraph(GModel *const model, Graph &graph, meshPartitionOptions &options, BoElemGrVec *const boElemGrVec)
{
enum {
@@ -516,15 +564,15 @@ int MakeGraph(GModel *const model, Graph &graph, BoElemGrVec *const boElemGrVec)
ElemTypePolyg = 2
};
- int ier = 0;
+ int ier = 0;
// switch(partitionScope) {
// case PartitionEntireDomain:
-
+ if(!options.partitionByExtrusion){
/*--------------------------------------------------------------------*
* Make a graph for the entire domain
*--------------------------------------------------------------------*/
- {
+
//--Get the dimension of the mesh and count the numbers of elements
@@ -584,6 +632,79 @@ int MakeGraph(GModel *const model, Graph &graph, BoElemGrVec *const boElemGrVec)
break;
}
}
+ else {
+/*--------------------------------------------------------------------*
+ * Create graph for extrusion source entity (propagate group ids later)
+ *--------------------------------------------------------------------*/
+ unsigned tmp[5];
+ const int meshDim = model->getNumMeshElements(tmp);
+ if(meshDim < 2) {
+ Msg::Error("No mesh elements were found");
+ return 1;
+ }
+ std::vector<const GEntity*> sourceEntities;
+ std::set<const GEntity*> bannedEntities;
+ unsigned numElem[5] = {0,0,0,0,0};
+ switch(meshDim) {
+ case 2:
+ for(GModel::fiter it = model->firstFace(); it != model->lastFace(); ++it){
+ ExtrudeParams* ep = (*it)->meshAttributes.extrude;
+ const GEdge *edge = model->getEdgeByTag(std::abs(ep->geo.Source));
+ if(!edge){
+ Msg::Error("Unknown source edge %d for extrusion", ep->geo.Source);
+ continue;
+ }
+ // mark face for partitioning if not associated to already processed face
+ if (bannedEntities.find(edge) == bannedEntities.end()) {
+ sourceEntities.push_back(edge);
+ edge->getNumMeshElements(numElem);
+ }
+ // do not allow partitioning of the extruded edges
+ std::list<GEdge*> edges = (*it)->edges();
+ for(std::list<GEdge*>::iterator eit = edges.begin(); eit != edges.end(); eit++) {
+ bannedEntities.insert((*eit));
+ }
+ }
+ break;
+ case 3:
+ for(GModel::riter it = model->firstRegion(); it != model->lastRegion(); ++it){
+ ExtrudeParams* ep = (*it)->meshAttributes.extrude;
+ const GFace *face = model->getFaceByTag(std::abs(ep->geo.Source));
+ if(!face){
+ Msg::Error("Unknown source surface %d for extrusion", ep->geo.Source);
+ continue;
+ }
+ // mark face for partitioning if not associated to already processed region
+ if (bannedEntities.find(face) == bannedEntities.end()) {
+ sourceEntities.push_back(face);
+ face->getNumMeshElements(numElem);
+ }
+ // do not allow partitioning of the extruded faces
+ std::list<GFace*> faces = (*it)->faces();
+ for(std::list<GFace*>::iterator fit = faces.begin(); fit != faces.end(); fit++) {
+ bannedEntities.insert((*fit));
+ }
+ }
+ break;
+ }
+ // build a graph only on the source entities
+ try {
+ // Allocate memory for the graph
+ const int numGrVert = numElem[ElemTypeTri] + numElem[ElemTypeQuad]
+ + numElem[ElemTypePolyg];
+ const int maxGrEdge = (numElem[ElemTypeTri]*3 + numElem[ElemTypeQuad]*4
+ + numElem[ElemTypePolyg]*4)/2;
+ graph.allocate(numGrVert, maxGrEdge);
+ // Make the graph
+ MakeGraphDIM<2>(sourceEntities.begin(), sourceEntities.end(),
+ model->firstEdge(), model->lastEdge(), graph,
+ boElemGrVec);
+ }
+ catch (...){
+ Msg::Error("Exception thrown during graph generation");
+ ier = 2;
+ }
+ }
// break;
// case PartitionByPhysical:
diff --git a/Mesh/meshPartitionOptions.h b/Mesh/meshPartitionOptions.h
index d4f73da384dcf2217ab964a63d6d7394fe9a7550..2867ab086c199ae2f31b4778e8a87a064e9f0ede 100644
--- a/Mesh/meshPartitionOptions.h
+++ b/Mesh/meshPartitionOptions.h
@@ -69,6 +69,8 @@ struct meshPartitionOptions
// 3 - Random boundary refinement (with
// minimization of connectivity
// along sub-domains)
+ int partitionByExtrusion; // if true, all extruded elements belong
+ // to the same partition as the source element
//--NODAL WEIGHT
std::vector<int> nodalWeights;
@@ -109,6 +111,7 @@ struct meshPartitionOptions
refine_algorithm = 3;
createPartitionBoundaries = true;
createGhostCells = true;
+ partitionByExtrusion =false;
}
};
diff --git a/benchmarks/extrude/partitionByExtrusion_spirale.geo b/benchmarks/extrude/partitionByExtrusion_spirale.geo
new file mode 100644
index 0000000000000000000000000000000000000000..401625dbd92c90af69489b991aed89a75a140829
--- /dev/null
+++ b/benchmarks/extrude/partitionByExtrusion_spirale.geo
@@ -0,0 +1,115 @@
+// Extrusion and rotation. Christophe Geuzaine, Feb 2002.
+
+lc = 0.005;
+lc2 = 4*lc;
+
+Point(1)={0.09,0,0,lc};
+Point(2)={0.11,0,0,lc};
+Point(3)={0,0.09,0,lc};
+Point(4)={0,0.11,0,lc};
+Point(5)={0.2,0.09,0,lc};
+Point(6)={0.2,0.11,0,lc};
+Point(7)={0.09,0.2,0,lc};
+Point(8)={0.11,0.2,0,lc};
+Point(9)={0.09,0.09,0,lc2};
+Point(10)={0.09,0.11,0,lc2};
+Point(11)={0.11,0.11,0,lc2};
+Point(12)={0.11,0.09,0,lc2};
+
+Line(1) = {1,2};
+Line(2) = {2,12};
+Line(3) = {12,5};
+Line(4) = {5,6};
+Line(5) = {6,11};
+Line(6) = {11,8};
+Line(7) = {8,7};
+Line(8) = {7,10};
+Line(9) = {10,4};
+Line(10) = {4,3};
+Line(11) = {3,9};
+Line(12) = {9,1};
+
+Line Loop(13) = {2,3,4,5,6,7,8,9,10,11,12,1};
+Plane Surface(14) = {13};
+
+turns = 1/3;
+zz = 0.01;
+cc = 0.1;
+
+Extrude Surface {news-1, {0,0,2*zz}, {0,0,1} , {cc,cc,0} , 0}
+ { Layers {2}; Recombine; };
+
+For j In {1:3}
+ Extrude Surface {news-1, {0,0,zz}, {0,0,1} , {cc,cc,0} , Pi/(40*(4-j))}
+ { Layers {1}; Recombine; };
+EndFor
+
+For j In {1:turns*8}
+ Extrude Surface {news-1, {0,0,10*zz}, {0,0,1} , {cc,cc,0} , Pi/4}
+ { Layers {10}; Recombine; };
+EndFor
+
+For j In {1:3}
+ Extrude Surface {news-1, {0,0,zz}, {0,0,1} , {cc,cc,0} , Pi/(40*j)}
+ { Layers {1}; Recombine; };
+EndFor
+
+Extrude Surface {news-1, {0,0,2*zz}, {0,0,1} , {cc,cc,0} , 0}
+ { Layers {2}; Recombine; };
+
+p = newp;
+
+Point(p+1000)={0.1, -0.02, 0,lc};
+Point(p+1001)={0.1, -0.01, 0,lc};
+Point(p+1002)={0.1, -0.02, 0.01,lc};
+Point(p+1003)={0.1, -0.03, 0,lc};
+Point(p+1004)={0.1, -0.02, -0.01,lc};
+
+Circle(1635) = {p+1001,p+1000,p+1002};
+Circle(1636) = {p+1002,p+1000,p+1003};
+Circle(1637) = {p+1003,p+1000,p+1004};
+Circle(1638) = {p+1004,p+1000,p+1001};
+
+Line Loop(1639) = {1635,1636,1637,1638};
+Plane Surface(1640) = {1639};
+
+turns = 2;
+zz = 0.01/4;
+
+For j In {1:8*turns}
+ Extrude Surface {news-1, {0,0,10*zz}, {0,0,1} , {cc,cc,0} , Pi/4}
+ { Layers {10}; Recombine; };
+EndFor
+
+Point(p+4000)={0.1, -0.02 + 0.24, 0, lc};
+Point(p+4001)={0.1, -0.01 + 0.24, 0, lc};
+Point(p+4002)={0.1, -0.02 + 0.24, 0.01, lc};
+Point(p+4003)={0.1, -0.03 + 0.24, 0, lc};
+Point(p+4004)={0.1, -0.02 + 0.24, -0.01, lc};
+Circle(4635) = {p+4001,p+4000,p+4002};
+Circle(4636) = {p+4002,p+4000,p+4003};
+Circle(4637) = {p+4003,p+4000,p+4004};
+Circle(4638) = {p+4004,p+4000,p+4001};
+
+Line Loop(4639) = {4635,4636,4637,4638};
+Plane Surface(4640) = {4639};
+
+For j In {1:8*turns}
+ Extrude Surface {news-1, {0,0,10*zz}, {0,0,1} , {cc,cc,0} , Pi/4}
+ { Layers {10}; Recombine; };
+EndFor
+
+// options
+Mesh.PartitionByExtrusion=1; // forces partitioning by columns
+Mesh.NbPartitions=20;
+Mesh.ColorCarousel = 3;
+Mesh.SurfaceFaces= 1;
+Mesh.VolumeEdges = 0;
+Mesh.VolumeFaces = 1;
+Mesh.Color.Zero = {255,0,0};
+Mesh.Color.One = {0,0,255};
+Mesh.Color.Two = {0,255,0};
+Mesh.Color.Three = {255,255,0};
+Mesh.Color.Four = {255,0,255};
+Mesh.Color.Six = {255,128,0};
+Mesh.Color.Seven = {64,128,32};