From fb7edb963708d3724bef2c0615348ed3e38abc75 Mon Sep 17 00:00:00 2001
From: Laurent Van Migroet <l.vanmiegroet@ulg.ac.be>
Date: Fri, 22 Feb 2008 17:58:12 +0000
Subject: [PATCH] Added getTypeOfElements  and addThisTypeOfElement to
 GEntities

added
getTypeOfElements  from Gentities and Physical Groups in GModel
---
 Geo/GEdge.cpp   | 35 +++++++++++-------
 Geo/GEdge.h     |  7 +++-
 Geo/GEntity.cpp | 30 +++++++++++----
 Geo/GEntity.h   | 16 +++++---
 Geo/GFace.cpp   | 98 ++++++++++++++++++++++++++++---------------------
 Geo/GFace.h     | 23 +++++++-----
 Geo/GModel.cpp  | 60 +++++++++++++++++++++++++++---
 Geo/GModel.h    |  6 +++
 Geo/GRegion.cpp | 47 ++++++++++++++++++------
 Geo/GRegion.h   |  8 +++-
 10 files changed, 232 insertions(+), 98 deletions(-)

diff --git a/Geo/GEdge.cpp b/Geo/GEdge.cpp
index 5aad415051..4fcb2ccb48 100644
--- a/Geo/GEdge.cpp
+++ b/Geo/GEdge.cpp
@@ -1,4 +1,4 @@
-// $Id: GEdge.cpp,v 1.40 2008-02-21 12:11:12 geuzaine Exp $
+// $Id: GEdge.cpp,v 1.41 2008-02-22 17:58:12 miegroet Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -16,7 +16,7 @@
 // along with this program; if not, write to the Free Software
 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
 // USA.
-// 
+//
 // Please report all bugs and problems to <gmsh@geuz.org>.
 
 #include <algorithm>
@@ -41,21 +41,21 @@ GEdge::GEdge(GModel *model, int tag, GVertex *_v0, GVertex *_v1)
   resetMeshAttributes();
 }
 
-GEdge::~GEdge() 
+GEdge::~GEdge()
 {
   if(v0) v0->delEdge(this);
   if(v1 && v1 != v0) v1->delEdge(this);
 
-  for(unsigned int i = 0; i < mesh_vertices.size(); i++) 
+  for(unsigned int i = 0; i < mesh_vertices.size(); i++)
     delete mesh_vertices[i];
 
-  for(unsigned int i = 0; i < lines.size(); i++) 
+  for(unsigned int i = 0; i < lines.size(); i++)
     delete lines[i];
 }
 
-void GEdge::resetMeshAttributes() 
-{ 
-  meshAttributes.Method = LIBRE; 
+void GEdge::resetMeshAttributes()
+{
+  meshAttributes.Method = LIBRE;
   meshAttributes.coeffTransfinite = 0.;
   meshAttributes.nbPointsTransfinite = 0;
   meshAttributes.typeTransfinite = 0;
@@ -64,12 +64,12 @@ void GEdge::resetMeshAttributes()
 }
 
 void GEdge::addFace(GFace *e)
-{ 
-  l_faces.push_back(e);  
+{
+  l_faces.push_back(e);
 }
 
 void GEdge::delFace(GFace *e)
-{ 
+{
   l_faces.erase(std::find(l_faces.begin(), l_faces.end(), e));
 }
 
@@ -129,7 +129,7 @@ int GEdge::containsParam(double pt) const
   return (pt >= rg.low() && pt <= rg.high());
 }
 
-SVector3 GEdge::secondDer(double par) const 
+SVector3 GEdge::secondDer(double par) const
 {
   // use central differences
   const double eps = 1.e-3;
@@ -146,7 +146,7 @@ SPoint2 GEdge::reparamOnFace(GFace *face, double epar,int dir) const
   return face->parFromPoint(sp3);
 }
 
-double GEdge::curvature(double par) const 
+double GEdge::curvature(double par) const
 {
   double eps1 = 1.e-5;
   double eps2 = 1.e-5;
@@ -190,3 +190,12 @@ double GEdge::length(const double &u0, const double &u1, const int nbQuadPoints)
   return L;
 #endif
 }
+
+void GEdge::getTypeOfElements(std::vector<int> &groups)
+{
+	for(unsigned int j = 0; j < lines.size(); j++)
+	{
+		int type = lines[j]->getTypeForMSH();
+		this->addThisTypeOfElement(type,groups);
+	}
+}
\ No newline at end of file
diff --git a/Geo/GEdge.h b/Geo/GEdge.h
index 487bb1591b..799322d997 100644
--- a/Geo/GEdge.h
+++ b/Geo/GEdge.h
@@ -17,7 +17,7 @@
 // along with this program; if not, write to the Free Software
 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
 // USA.
-// 
+//
 // Please report all bugs and problems to <gmsh@geuz.org>.
 
 #include "GEntity.h"
@@ -76,7 +76,7 @@ class GEdge : public GEntity {
   // Get second derivative of edge at the given parameter.
   // Default implentation using central differences
   virtual SVector3 secondDer(double par) const;
-  virtual double curvature(double par) const;  
+  virtual double curvature(double par) const;
 
   // Reparmaterize the point onto the given face.
   virtual SPoint2 reparamOnFace(GFace *face, double epar,int dir) const;
@@ -113,6 +113,9 @@ class GEdge : public GEntity {
   // True if start == end and no more than 2 segments
   bool isMeshDegenerated() const{ return (v0 == v1 && mesh_vertices.size() < 2); }
 
+// Returns each type of element in the gentity
+  void getTypeOfElements(std::vector<int> &groups);
+
   struct {
     char   Method;
     double coeffTransfinite;
diff --git a/Geo/GEntity.cpp b/Geo/GEntity.cpp
index 040e2c269a..89288341dc 100644
--- a/Geo/GEntity.cpp
+++ b/Geo/GEntity.cpp
@@ -1,4 +1,4 @@
-// $Id: GEntity.cpp,v 1.18 2008-02-17 08:47:58 geuzaine Exp $
+// $Id: GEntity.cpp,v 1.19 2008-02-22 17:58:12 miegroet Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -16,7 +16,7 @@
 // along with this program; if not, write to the Free Software
 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
 // USA.
-// 
+//
 // Please report all bugs and problems to <gmsh@geuz.org>.
 
 #include "GEntity.h"
@@ -49,15 +49,15 @@ void GEntity::deleteVertexArrays()
 }
 
 char GEntity::getVisibility()
-{ 
+{
   if(CTX.hide_unselected && !CTX.pick_elements && !getSelection() &&
-     geomType() != ProjectionFace) 
+     geomType() != ProjectionFace)
     return false;
-  return _visible; 
+  return _visible;
 }
 
 bool GEntity::useColor()
-{ 
+{
   int r = CTX.UNPACK_RED(_color);
   int g = CTX.UNPACK_GREEN(_color);
   int b = CTX.UNPACK_BLUE(_color);
@@ -75,7 +75,7 @@ std::string GEntity::getInfoString()
   std::string out = getTypeString() + tmp;
 
   std::string info = getAdditionalInfoString();
-  if(info.size()) 
+  if(info.size())
     out += " " + info;
 
   if(physicals.size()){
@@ -90,3 +90,19 @@ std::string GEntity::getInfoString()
 
   return out;
 }
+
+
+inline
+void GEntity::addThisTypeOfElement(int type,std::vector<int> &groups)
+{
+	unsigned int size=groups.size();
+	if(size==0)
+			groups.push_back(type);
+	for(unsigned int i=0;i<size;i++)
+	{
+		if(type==groups[i])
+			break;
+		else
+			groups.push_back(type);
+	}
+}
diff --git a/Geo/GEntity.h b/Geo/GEntity.h
index 14474b77ba..bb5c43dc5c 100644
--- a/Geo/GEntity.h
+++ b/Geo/GEntity.h
@@ -17,12 +17,13 @@
 // along with this program; if not, write to the Free Software
 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
 // USA.
-// 
+//
 // Please report all bugs and problems to <gmsh@geuz.org>.
 
 #include <list>
 #include <vector>
 #include <string>
+#include <map>
 #include "Range.h"
 #include "SPoint3.h"
 #include "SBoundingBox3d.h"
@@ -40,19 +41,19 @@ class GEntity {
  private:
   // All entities are owned by a GModel
   GModel *_model;
-  
+
   // The tag (the number) of this entity
   int _tag;
-  
+
   // The visibility and the selection flag
   char _visible, _selection;
-  
+
   // Flag storing if all mesh elements are visible
   char _allElementsVisible;
 
   // The color of the entity (ignored if set to transparent blue)
   unsigned int _color;
-  
+
  public:
 
   // All known native model types
@@ -149,7 +150,7 @@ class GEntity {
 
   void deleteVertexArrays();
 
-  // Spatial dimension of the entity 
+  // Spatial dimension of the entity
   virtual int dim() const { throw; }
 
   // Regions that bound this entity or that this entity bounds.
@@ -244,6 +245,9 @@ class GEntity {
   // Vertex arrays to draw the mesh efficiently
   VertexArray *va_lines, *va_triangles;
 
+  // Returns all type of element in the GEntity
+  virtual void getTypeOfElements(std::vector<int>  &)const { throw; }
+  virtual void addThisTypeOfElement(int type,std::vector<int> &groups);
 };
 
 class GEntityLessThan {
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 8aa5c49c73..6ce8a0b2d8 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -1,4 +1,4 @@
-// $Id: GFace.cpp,v 1.53 2008-02-21 12:11:12 geuzaine Exp $
+// $Id: GFace.cpp,v 1.54 2008-02-22 17:58:12 miegroet Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -16,7 +16,7 @@
 // along with this program; if not, write to the Free Software
 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
 // USA.
-// 
+//
 // Please report all bugs and problems to <gmsh@geuz.org>.
 
 #include "GModel.h"
@@ -44,7 +44,7 @@ void dsvdcmp(double **a, int m, int n, double w[], double **v);
 
 extern Context_T CTX;
 
-GFace::GFace(GModel *model, int tag) 
+GFace::GFace(GModel *model, int tag)
   : GEntity(model, tag), r1(0), r2(0), va_geom_triangles(0)
 {
   meshStatistics.status = GFace::PENDING;
@@ -52,7 +52,7 @@ GFace::GFace(GModel *model, int tag)
 }
 
 GFace::~GFace()
-{ 
+{
   std::list<GEdge*>::iterator it = l_edges.begin();
 
   while (it != l_edges.end()){
@@ -60,13 +60,13 @@ GFace::~GFace()
     ++it;
   }
 
-  for(unsigned int i = 0; i < mesh_vertices.size(); i++) 
+  for(unsigned int i = 0; i < mesh_vertices.size(); i++)
     delete mesh_vertices[i];
 
-  for(unsigned int i = 0; i < triangles.size(); i++) 
+  for(unsigned int i = 0; i < triangles.size(); i++)
     delete triangles[i];
 
-  for(unsigned int i = 0; i < quadrangles.size(); i++) 
+  for(unsigned int i = 0; i < quadrangles.size(); i++)
     delete quadrangles[i];
 
   if(va_geom_triangles)
@@ -174,7 +174,7 @@ void GFace::computeMeanPlane()
   std::list<GVertex*> verts = vertices();
   std::list<GVertex*>::const_iterator itv = verts.begin();
   for(; itv != verts.end(); itv++){
-    const GVertex *v = *itv; 
+    const GVertex *v = *itv;
     pts.push_back(SPoint3(v->x(), v->y(), v->z()));
   }
 
@@ -183,7 +183,7 @@ void GFace::computeMeanPlane()
     std::list<GEdge*> edg = edges();
     std::list<GEdge*>::const_iterator ite = edg.begin();
     for(; ite != edg.end(); ite++){
-      const GEdge *e = *ite; 
+      const GEdge *e = *ite;
       Range<double> b = e->parBounds(0);
       GPoint p1 = e->point(b.low() + 0.333 * (b.high() - b.low()));
       pts.push_back(SPoint3(p1.x(), p1.y(), p1.z()));
@@ -305,7 +305,7 @@ void GFace::computeMeanPlane(const std::vector<SPoint3> &points)
     norme(t2);
     // prodve(t1, t2, res2);
     // Warning: the rest of the code assumes res = t2 x t1, not t1 x t2 (WTF?)
-    prodve(t2, t1, res2); 
+    prodve(t2, t1, res2);
     norme(res2);
     prodve(res, res2, c);
     prosca(res, res2, &cosc);
@@ -350,27 +350,27 @@ end:
   meanPlane.d = res[3];
 
   meanPlane.x = meanPlane.y = meanPlane.z = 0.;
-  if(fabs(meanPlane.a) >= fabs(meanPlane.b) && 
+  if(fabs(meanPlane.a) >= fabs(meanPlane.b) &&
      fabs(meanPlane.a) >= fabs(meanPlane.c) ){
     meanPlane.x = meanPlane.d / meanPlane.a;
   }
-  else if(fabs(meanPlane.b) >= fabs(meanPlane.a) && 
+  else if(fabs(meanPlane.b) >= fabs(meanPlane.a) &&
 	  fabs(meanPlane.b) >= fabs(meanPlane.c)){
     meanPlane.y = meanPlane.d / meanPlane.b;
   }
   else{
     meanPlane.z = meanPlane.d / meanPlane.c;
   }
-  
+
   Msg(DEBUG1, "Surface: %d", tag());
   Msg(DEBUG2, "SVD    : %g,%g,%g (min=%d)", svd[0], svd[1], svd[2], min);
-  Msg(DEBUG2, "Plane  : (%g x + %g y + %g z = %g)", 
+  Msg(DEBUG2, "Plane  : (%g x + %g y + %g z = %g)",
       meanPlane.a, meanPlane.b, meanPlane.c, meanPlane.d);
-  Msg(DEBUG2, "Normal : (%g , %g , %g )", 
+  Msg(DEBUG2, "Normal : (%g , %g , %g )",
       meanPlane.a, meanPlane.b, meanPlane.c);
   Msg(DEBUG3, "t1     : (%g , %g , %g )", t1[0], t1[1], t1[2]);
   Msg(DEBUG3, "t2     : (%g , %g , %g )", t2[0], t2[1], t2[2]);
-  Msg(DEBUG3, "pt     : (%g , %g , %g )", 
+  Msg(DEBUG3, "pt     : (%g , %g , %g )",
       meanPlane.x, meanPlane.y, meanPlane.z);
 
   //check coherence for plane surfaces
@@ -380,8 +380,8 @@ end:
     std::list<GVertex*> verts = vertices();
     std::list<GVertex*>::const_iterator itv = verts.begin();
     for(; itv != verts.end(); itv++){
-      const GVertex *v = *itv; 
-      double d = meanPlane.a * v->x() + meanPlane.b * v->y() + 
+      const GVertex *v = *itv;
+      double d = meanPlane.a * v->x() + meanPlane.b * v->y() +
 	meanPlane.c * v->z() - meanPlane.d;
       if(fabs(d) > lc * 1.e-3) {
 	Msg(GERROR1, "Plane surface %d (%gx+%gy+%gz+%g=0) is not plane!",
@@ -395,7 +395,7 @@ end:
 #endif
 }
 
-void GFace::getMeanPlaneData(double VX[3], double VY[3], 
+void GFace::getMeanPlaneData(double VX[3], double VY[3],
 			     double &x, double &y, double &z) const
 {
   VX[0] = meanPlane.plan[0][0];
@@ -404,9 +404,9 @@ void GFace::getMeanPlaneData(double VX[3], double VY[3],
   VY[0] = meanPlane.plan[1][0];
   VY[1] = meanPlane.plan[1][1];
   VY[2] = meanPlane.plan[1][2];
-  x = meanPlane.x;  
-  y = meanPlane.y;  
-  z = meanPlane.z;  
+  x = meanPlane.x;
+  y = meanPlane.y;
+  z = meanPlane.z;
 }
 
 void GFace::getMeanPlaneData(double plan[3][3]) const
@@ -445,14 +445,14 @@ double GFace::curvature (const SPoint2 &param) const
   SVector3 dndu = 500 * (n2 - n1);
   SVector3 dndv = 500 * (n4 - n3);
 
-  double c = fabs(dot(dndu, du) +  dot(dndv, dv)) / detJ; 
+  double c = fabs(dot(dndu, du) +  dot(dndv, dv)) / detJ;
 
   // Msg(INFO, "c = %g detJ %g", c, detJ);
 
   return c;
 }
 
-void GFace::XYZtoUV(const double X, const double Y, const double Z, 
+void GFace::XYZtoUV(const double X, const double Y, const double Z,
 		    double &U, double &V, const double relax,
 		    const bool onSurface) const
 {
@@ -467,7 +467,7 @@ void GFace::XYZtoUV(const double X, const double Y, const double Z,
   double umin, umax, vmin, vmax;
   double initu[NumInitGuess] = {0.487, 0.6, 0.4, 0.7, 0.3, 0.8, 0.2, 1., 0., 0., 1.};
   double initv[NumInitGuess] = {0.487, 0.6, 0.4, 0.7, 0.3, 0.8, 0.2, 0., 1., 0., 1.};
-  
+
   Range<double> ru = parBounds(0);
   Range<double> rv = parBounds(1);
   umin = ru.low();
@@ -477,7 +477,7 @@ void GFace::XYZtoUV(const double X, const double Y, const double Z,
 
   for(int i = 0; i < NumInitGuess; i++){
     for(int j = 0; j < NumInitGuess; j++){
-    
+
       U = initu[i];
       V = initv[j];
       err = 1.0;
@@ -489,7 +489,7 @@ void GFace::XYZtoUV(const double X, const double Y, const double Z,
 
       while(err > Precision && iter < MaxIter) {
 	P = point(U,V);
-	Pair<SVector3,SVector3> der = firstDer(SPoint2(U,V)); 
+	Pair<SVector3,SVector3> der = firstDer(SPoint2(U,V));
 	mat[0][0] = der.left().x();
 	mat[0][1] = der.left().y();
 	mat[0][2] = der.left().z();
@@ -500,38 +500,38 @@ void GFace::XYZtoUV(const double X, const double Y, const double Z,
 	mat[2][1] = 0.;
 	mat[2][2] = 0.;
 	invert_singular_matrix3x3(mat, jac);
-	
+
 	Unew = U + relax *
 	  (jac[0][0] * (X - P.x()) + jac[1][0] * (Y - P.y()) +
 	   jac[2][0] * (Z - P.z()));
-	Vnew = V + relax * 
+	Vnew = V + relax *
 	  (jac[0][1] * (X - P.x()) + jac[1][1] * (Y - P.y()) +
 	   jac[2][1] * (Z - P.z()));
-	
+
 	err = DSQR(Unew - U) + DSQR(Vnew - V);
 	err2 = sqrt(DSQR(X - P.x()) + DSQR(Y - P.y()) + DSQR(Z - P.z()));
 	iter++;
 	U = Unew;
 	V = Vnew;
       }
-      
-      if(iter < MaxIter && err <= Precision && 
-	 Unew <= umax && Vnew <= vmax && 
+
+      if(iter < MaxIter && err <= Precision &&
+	 Unew <= umax && Vnew <= vmax &&
 	 Unew >= umin && Vnew >= vmin){
 	if (onSurface && err2 > 1.e-4 * CTX.lc)
-	  Msg(WARNING,"Converged for i=%d j=%d (err=%g iter=%d) BUT xyz error = %g", 
+	  Msg(WARNING,"Converged for i=%d j=%d (err=%g iter=%d) BUT xyz error = %g",
 	      i, j, err, iter, err2);
-	return;	
+	return;
       }
     }
   }
-  
+
   if(relax < 1.e-6)
     Msg(GERROR, "Could not converge: surface mesh will be wrong");
   else {
     Msg(INFO, "point %g %g %g : Relaxation factor = %g",X,Y,Z, 0.75 * relax);
     XYZtoUV(X, Y, Z, U, V, 0.75 * relax);
-  }  
+  }
 #endif
 }
 
@@ -546,10 +546,10 @@ int GFace::containsParam(const SPoint2 &pt) const
 {
   Range<double> uu = parBounds(0);
   Range<double> vv = parBounds(1);
-  if((pt.x() >= uu.low() && pt.x() <= uu.high()) && 
+  if((pt.x() >= uu.low() && pt.x() <= uu.high()) &&
      (pt.y() >= vv.low() && pt.y() <= vv.high()))
     return 1;
-  else 
+  else
     return 0;
 }
 
@@ -566,7 +566,7 @@ bool GFace::buildRepresentationCross()
   SBoundingBox3d bb;
   for(std::list<GEdge*>::iterator it = ed.begin(); it != ed.end(); it++){
     GEdge *ge = *it;
-    if(ge->geomType() == GEntity::DiscreteCurve || 
+    if(ge->geomType() == GEntity::DiscreteCurve ||
        ge->geomType() == GEntity::BoundaryLayerCurve){
       // don't try again
       cross.clear();
@@ -710,7 +710,7 @@ double GFace::length(const SPoint2 &pt1, const SPoint2 &pt2, int nbQuadPoints)
   for (int i = 0; i < nbQuadPoints; i++){
     const double ti = 0.5 * (1. + t[i]);
     SPoint2 pi = geodesic(pt1, pt2, ti);
-    Pair<SVector3, SVector3> der2 = firstDer(pi);    
+    Pair<SVector3, SVector3> der2 = firstDer(pi);
     SVector3 der = der2.left() * (pt2.x() - pt1.x()) + der2.right() * (pt2.y() - pt1.y());
     const double d = norm(der);
     L += d * w[i] ;
@@ -718,3 +718,19 @@ double GFace::length(const SPoint2 &pt1, const SPoint2 &pt2, int nbQuadPoints)
   return L;
 #endif
 }
+
+void GFace::getTypeOfElements(std::vector<int> &groups)
+{
+ 	for(unsigned int j = 0; j < triangles.size(); j++)
+	{
+		int type = triangles[j]->getTypeForMSH();
+		addThisTypeOfElement(type,groups);
+	}
+ 	for(unsigned int j = 0; j < quadrangles.size(); j++)
+	{
+		int type = quadrangles[j]->getTypeForMSH();
+		addThisTypeOfElement(type,groups);
+	}
+}
+
+
diff --git a/Geo/GFace.h b/Geo/GFace.h
index 0135f3acec..f4e1bef707 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -17,7 +17,7 @@
 // along with this program; if not, write to the Free Software
 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
 // USA.
-// 
+//
 // Please report all bugs and problems to <gmsh@geuz.org>.
 
 #include "GEntity.h"
@@ -45,10 +45,10 @@ struct surface_params
 
 class GRegion;
 
-// A model face. 
-class GFace : public GEntity 
+// A model face.
+class GFace : public GEntity
 {
- protected: 
+ protected:
   // edge loops, will replace what follows list of al the edges of the
   // face
   std::list<GEdge *> l_edges;
@@ -87,14 +87,14 @@ class GFace : public GEntity
   virtual void setVisibility(char val, bool recursive=false);
 
   // compute the parameters UV from a point XYZ
-  void XYZtoUV(const double X, const double Y, const double Z, 
-	       double &U, double &V, const double relax, 
+  void XYZtoUV(const double X, const double Y, const double Z,
+	       double &U, double &V, const double relax,
 	       const bool onSurface=true) const;
 
   // The bounding box
-  virtual SBoundingBox3d bounds() const; 
+  virtual SBoundingBox3d bounds() const;
 
-  // retrieve surface params 
+  // retrieve surface params
   virtual surface_params getSurfaceParams() const { throw; }
 
   // Return the point on the face corresponding to the given parameter.
@@ -133,7 +133,7 @@ class GFace : public GEntity
 
   // Return the curvature i.e. the divergence of the normal
   virtual double curvature(const SPoint2 &param) const;
-  
+
   // Recompute the mesh partitions defined on this face.
   void recomputeMeshPartitions();
 
@@ -158,7 +158,7 @@ class GFace : public GEntity
   void computeMeanPlane();
 
   // Get the mean plane info
-  void getMeanPlaneData(double VX[3], double VY[3], 
+  void getMeanPlaneData(double VX[3], double VY[3],
 			double &x, double &y, double &z) const;
   void getMeanPlaneData(double plan[3][3]) const;
 
@@ -206,6 +206,9 @@ class GFace : public GEntity
 
   std::vector<MTriangle*> triangles;
   std::vector<MQuadrangle*> quadrangles;
+
+  // Returns each type of element in the gentity
+  void getTypeOfElements(std::vector<int> &groups);
 };
 
 #endif
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 984eda8ba8..b42671dc07 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -1,4 +1,4 @@
-// $Id: GModel.cpp,v 1.62 2008-02-22 07:49:38 geuzaine Exp $
+// $Id: GModel.cpp,v 1.63 2008-02-22 17:58:12 miegroet Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -171,14 +171,14 @@ void GModel::remove(GVertex *v)
 void GModel::snapVertices()
 {
   viter vit = firstVertex();
-  
-  double tol = CTX.geom.tolerance; 	       
-  
+
+  double tol = CTX.geom.tolerance;
+
   while (vit != lastVertex()){
     std::list<GEdge*> edges = (*vit)->edges();
     for (std::list<GEdge*>::iterator it = edges.begin(); it != edges.end(); ++it){
       Range<double> parb = (*it)->parBounds(0);
-      double t;	
+      double t;
       if ((*it)->getBeginVertex() == *vit){
 	t = parb.low();
       }
@@ -404,7 +404,7 @@ void GModel::buildMeshVertexCache()
   _vertexVectorCache.clear();
   _vertexMapCache.clear();
   bool dense = (getNumMeshVertices() == MVertex::getGlobalNumber());
- 
+
   if(dense){
     _vertexVectorCache.resize(MVertex::getGlobalNumber());
     for(viter it = firstVertex(); it != lastVertex(); ++it)
@@ -720,3 +720,51 @@ void GModel::checkMeshCoherence()
     MElementLessThanLexicographic::tolerance = old_tol;
   }
 }
+
+//maybe UseLess as we can get each vector directly through the GEntities
+void GModel::getTypeOfElements(std::map<GEntity*, std::vector<int> > groups[3])
+{
+  //for(viter it = firstVertex(); it != lastVertex(); ++it)
+  //	(*it)->getTypeOfElements(groups[0][(*it)]);
+
+  for(eiter it = firstEdge(); it != lastEdge(); ++it)
+  	(*it)->getTypeOfElements(groups[0][(*it)]);
+
+
+  for(fiter it = firstFace(); it != lastFace(); ++it)
+	(*it)->getTypeOfElements(groups[1][(*it)]);
+
+  for(riter it = firstRegion(); it != lastRegion(); ++it)
+	(*it)->getTypeOfElements(groups[2][(*it)]);
+
+}
+
+void GModel::getTypeOfElements(std::map<int, std::vector<int> > groups[3])
+{
+	//for(unsigned int i=0;i<(*it)->physicals.size();i++)
+	//	(*it)->getTypeOfElements(groups[0][i]);
+
+	for(eiter it = firstEdge(); it != lastEdge(); ++it){
+		if((*it)->physicals.size()){
+		(*it)->getTypeOfElements(groups[0][(*it)->physicals[0]]);
+		for(unsigned int i=1;i<(*it)->physicals.size();i++)
+			groups[0][(*it)->physicals[i]]=groups[0][(*it)->physicals[i-1]];
+		}
+	}
+
+	for(fiter it = firstFace(); it != lastFace(); ++it){
+		if((*it)->physicals.size()){
+		(*it)->getTypeOfElements(groups[1][(*it)->physicals[0]]);
+		for(unsigned int i=1;i<(*it)->physicals.size();i++)
+			groups[1][(*it)->physicals[i]]=groups[1][(*it)->physicals[i-1]];
+		}
+	}
+	for(riter it = firstRegion(); it != lastRegion(); ++it){
+		if((*it)->physicals.size()){
+		(*it)->getTypeOfElements(groups[2][(*it)->physicals[0]]);
+		for(unsigned int i=1;i<(*it)->physicals.size();i++)
+			groups[2][(*it)->physicals[i]]=groups[2][(*it)->physicals[i-1]];
+		}
+	}
+
+}
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 3c20afcf26..8154579d5c 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -121,6 +121,12 @@ class GModel
   // Returns all physical groups (one map per dimension: 0-D to 3-D)
   void getPhysicalGroups(std::map<int, std::vector<GEntity*> > groups[4]);
 
+  // Returns all type of element in each GEntity (one map per dimension: 1-D to 3-D)
+  void getTypeOfElements(std::map<GEntity*, std::vector<int> > groups[3]);
+
+  // Returns all type of element in each physical group (one map per dimension: 1-D to 3-D)
+  void getTypeOfElements(std::map<int, std::vector<int> > groups[3]);
+
   // Deletes physical groups in the model
   void deletePhysicalGroups();
   void deletePhysicalGroup(int dim, int num);
diff --git a/Geo/GRegion.cpp b/Geo/GRegion.cpp
index fd403b58e1..b4a63c2c5e 100644
--- a/Geo/GRegion.cpp
+++ b/Geo/GRegion.cpp
@@ -1,4 +1,4 @@
-// $Id: GRegion.cpp,v 1.22 2008-02-17 08:47:58 geuzaine Exp $
+// $Id: GRegion.cpp,v 1.23 2008-02-22 17:58:12 miegroet Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -16,7 +16,7 @@
 // along with this program; if not, write to the Free Software
 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
 // USA.
-// 
+//
 // Please report all bugs and problems to <gmsh@geuz.org>.
 
 #include "GModel.h"
@@ -36,32 +36,32 @@ GRegion::GRegion(GModel *model, int tag) : GEntity (model, tag)
 }
 
 GRegion::~GRegion()
-{ 
+{
   std::list<GFace*>::iterator it = l_faces.begin();
   while(it != l_faces.end()){
     (*it)->delRegion(this);
     ++it;
   }
 
-  for(unsigned int i = 0; i < mesh_vertices.size(); i++) 
+  for(unsigned int i = 0; i < mesh_vertices.size(); i++)
     delete mesh_vertices[i];
 
-  for(unsigned int i = 0; i < tetrahedra.size(); i++) 
+  for(unsigned int i = 0; i < tetrahedra.size(); i++)
     delete tetrahedra[i];
 
-  for(unsigned int i = 0; i < hexahedra.size(); i++) 
+  for(unsigned int i = 0; i < hexahedra.size(); i++)
     delete hexahedra[i];
 
-  for(unsigned int i = 0; i < prisms.size(); i++) 
+  for(unsigned int i = 0; i < prisms.size(); i++)
     delete prisms[i];
 
-  for(unsigned int i = 0; i < pyramids.size(); i++) 
+  for(unsigned int i = 0; i < pyramids.size(); i++)
     delete pyramids[i];
 }
 
 void GRegion::resetMeshAttributes()
 {
-  meshAttributes.Method = LIBRE; 
+  meshAttributes.Method = LIBRE;
   meshAttributes.extrude = 0;
 }
 
@@ -71,7 +71,7 @@ SBoundingBox3d GRegion::bounds() const
   if(geomType() != DiscreteVolume){
     std::list<GFace*>::const_iterator it = l_faces.begin();
     for(; it != l_faces.end(); it++)
-      res += (*it)->bounds();  
+      res += (*it)->bounds();
   }
   else{
     for(unsigned int i = 0; i < mesh_vertices.size(); i++)
@@ -168,7 +168,7 @@ bool GRegion::edgeConnected(GRegion *r) const
 {
   std::list<GEdge*> e = edges();
   std::list<GEdge*> e2 = r->edges();
-  
+
   std::list<GEdge*>::const_iterator it = e.begin();
   while(it != e.end()){
     if(std::find(e2.begin(), e2.end(), *it) != e2.end())
@@ -177,3 +177,28 @@ bool GRegion::edgeConnected(GRegion *r) const
   }
   return false;
 }
+
+
+
+void GRegion::getTypeOfElements(std::vector<int>  &groups)
+{
+  for(unsigned int i = 0; i < tetrahedra.size(); i++){
+    int type = tetrahedra[i]->getTypeForMSH();
+	addThisTypeOfElement(type,groups);
+  }
+
+  for(unsigned int i = 0; i < hexahedra.size(); i++){
+    int type = hexahedra[i]->getTypeForMSH();
+	addThisTypeOfElement(type,groups);
+  }
+
+  for(unsigned int i = 0; i < prisms.size(); i++){
+    int type = prisms[i]->getTypeForMSH();
+  	addThisTypeOfElement(type,groups);
+  }
+
+  for(unsigned int i = 0; i < pyramids.size(); i++)  {
+    int type = pyramids[i]->getTypeForMSH();
+  	addThisTypeOfElement(type,groups);
+	}
+}
\ No newline at end of file
diff --git a/Geo/GRegion.h b/Geo/GRegion.h
index 547889ef73..c9a3d86207 100644
--- a/Geo/GRegion.h
+++ b/Geo/GRegion.h
@@ -17,7 +17,7 @@
 // along with this program; if not, write to the Free Software
 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
 // USA.
-// 
+//
 // Please report all bugs and problems to <gmsh@geuz.org>.
 
 #include "GEntity.h"
@@ -45,7 +45,7 @@ class GRegion : public GEntity {
   void set(std::list<GFace*> &f) { l_faces= f; }
 
   // The bounding box
-  virtual SBoundingBox3d bounds() const; 
+  virtual SBoundingBox3d bounds() const;
 
   // Checks if the region is connected to another region by an edge
   bool edgeConnected(GRegion *r) const;
@@ -74,6 +74,10 @@ class GRegion : public GEntity {
   // indices
   std::vector<std::vector<std::vector<MVertex*> > > transfinite_vertices;
 
+// Returns each type of element in the gentity
+  void getTypeOfElements(std::vector<int> &groups);
+
+
   std::vector<MTetrahedron*> tetrahedra;
   std::vector<MHexahedron*> hexahedra;
   std::vector<MPrism*> prisms;
-- 
GitLab