hex7 elements

libsrc/meshing/classifyhpel.hpp

@@ -1698,6 +1698,49 @@ HPREF_ELEMENT_TYPE ClassifyHex(HPRefElement & el, INDEX_2_HASHTABLE<int> & edges
 
 }
 
+
+
+
+
+
+HPREF_ELEMENT_TYPE ClassifyHex7 (HPRefElement & el, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HASHTABLE<int> & edgepoint_dom, 
+                                 NgBitArray & cornerpoint, NgBitArray & edgepoint, INDEX_3_HASHTABLE<int> & faces, INDEX_2_HASHTABLE<int> & face_edges, 
+                                 INDEX_2_HASHTABLE<int> & surf_edges, NgArray<int, PointIndex::BASE> & facepoint)
+{
+  HPREF_ELEMENT_TYPE type = HP_NONE;
+  
+  // no singular
+  // singular bottom
+  // singular top
+  
+  // indices of bot,top-faces combinations
+  int index[6][2] = {{0,1},{1,0},{2,4},{4,2},{3,5},{5,3}}; 
+  int p[8]; 
+  const ELEMENT_FACE * elfaces  = MeshTopology::GetFaces1 (HEX);
+  const ELEMENT_EDGE * eledges = MeshTopology::GetEdges1 (HEX);
+
+  INDEX_3 fbot = { el.pnums[0], el.pnums[1], el.pnums[2] };
+  INDEX_3 ftop = { el.pnums[4], el.pnums[5], el.pnums[6] };
+  fbot.Sort();
+  ftop.Sort();
+  
+  bool singbot = faces.Used(fbot);
+  bool singtop = faces.Used(ftop);
+
+  if (singbot)
+    el.type =  HP_HEX7_1FA;
+  else if (singtop)
+    el.type = HP_HEX7_1FB;
+  else
+    el.type = HP_HEX7;
+
+  return el.type;
+}
+
+
+
+
+
 HPREF_ELEMENT_TYPE ClassifySegm(HPRefElement & hpel, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HASHTABLE<int> & edgepoint_dom, 
                                 NgBitArray & cornerpoint, NgBitArray & edgepoint, INDEX_3_HASHTABLE<int> & faces, INDEX_2_HASHTABLE<int> & face_edges, 
                                 INDEX_2_HASHTABLE<int> & surf_edges, NgArray<int, PointIndex::BASE> & facepoint)

libsrc/meshing/hpref_hex.hpp

@@ -54,6 +54,8 @@ HPRef_Struct refhex_1f_0e_0v =
 
 
 
+
+
 // HP_HEX_1FA_1FB  ... face (1 - 4 - 3 -2) and face (1-2-6-5) singular 
 int refhex_1fa_1fb_0e_0v_splitedges[][3] =
   {
@@ -234,3 +236,95 @@ HPRef_Struct refhex_1fa_1fb_0e_0v =
     };
 
 
+
+
+
+
+// HP_HEX  ... no refinement
+int refhex7_splitedges[][3] =
+  {
+      { 0, 0, 0 }
+  };
+HPREF_ELEMENT_TYPE refhex7_newelstypes[] =
+  {
+    HP_HEX7,
+    HP_NONE,
+  };
+int refhex7_newels[][8] =
+  {
+    { 1, 2, 3, 4, 5, 6, 7 }
+  };
+HPRef_Struct refhex7 =
+  {
+    HP_HEX7,
+    refhex7_splitedges, 
+    0, 0,
+    refhex7_newelstypes, 
+    refhex7_newels
+  };
+
+// HP_HEX_1FA  ... face (1 - 4 - 3 -2) singular 
+int refhex7_1fa_splitedges[][3] =
+  {
+    { 1, 5, 8 },
+    { 2, 6, 9 },
+    { 3, 7, 10 },
+    { 4, 7, 11 }, 
+    { 0, 0, 0 }
+  };
+HPREF_ELEMENT_TYPE refhex7_1fa_newelstypes[] =
+  {
+    HP_HEX_1F_0E_0V,
+    HP_HEX7,
+    HP_NONE,
+  };
+int  refhex7_1fa_newels[][8] =
+  {
+    { 1, 2, 3, 4, 8, 9, 10, 11 },
+    { 8, 9, 10, 11, 5, 6, 7 }
+ }; 
+HPRef_Struct refhex7_1fa =
+  {
+    HP_HEX7,
+    refhex7_1fa_splitedges, 
+    0, 0,
+    refhex7_1fa_newelstypes, 
+    refhex7_1fa_newels
+  };
+
+
+
+
+// HP_HEX_1FB  ... face (5,6,7) singular 
+int refhex7_1fb_splitedges[][3] =
+  {
+    { 5, 1, 8 },
+    { 6, 2, 9 },
+    { 7, 3, 10 },
+    { 7, 4, 11 }, 
+    { 0, 0, 0 }
+  };
+HPREF_ELEMENT_TYPE refhex7_1fb_newelstypes[] =
+  {
+    HP_HEX,
+    HP_HEX7_1FB,
+    HP_NONE,
+  };
+int  refhex7_1fb_newels[][8] =
+  {
+    { 1, 2, 3, 4, 8, 9, 10, 11 },
+    { 8, 9, 10, 11, 5, 6, 7 }
+ }; 
+HPRef_Struct refhex7_1fb =
+  {
+    HP_HEX7,
+    refhex7_1fb_splitedges, 
+    0, 0,
+    refhex7_1fb_newelstypes, 
+    refhex7_1fb_newels
+  };
+
+
+
+
+

libsrc/meshing/hpref_tet.hpp

@@ -3007,14 +3007,16 @@ int reftet_1f_0e_1va_splitedges[][3] =
 };
 HPREF_ELEMENT_TYPE reftet_1f_0e_1va_newelstypes[] =
 {
-  HP_HEX_1F_0E_0V,
+  // HP_HEX_1F_0E_0V,
+  HP_HEX7_1FA,
   HP_TET_1F_0E_1VA,
   HP_TET,
   HP_NONE,
 };
 int reftet_1f_0e_1va_newels[][8] =
 {
-  { 3, 6, 7, 4, 8, 5, 5, 9 },
+  // { 3, 6, 7, 4, 8, 5, 5, 9 },
+  { 4, 3, 6, 7, 9, 8, 5 },  
   { 5, 2, 6, 7 },
   { 5, 9, 8, 1 },
 };
@@ -3247,6 +3249,7 @@ HPREF_ELEMENT_TYPE reftet_1f_1e_2va_newelstypes[] =
   HP_TET_1F_0E_1VA,
   HP_TET_1E_1VA,  
   HP_TET_1E_1VA,
+  HP_HEX7_1FB,
   HP_NONE,  
 };
 int reftet_1f_1e_2va_newels[][8] =
@@ -3257,7 +3260,9 @@ int reftet_1f_1e_2va_newels[][8] =
   { 10, 3, 12, 11 },
   { 14, 2, 8, 9 },
   { 2, 7, 15, 14 },
-  { 2, 9, 14, 15 } 
+  { 2, 9, 14, 15 },
+  // { 13, 10, 14, 15, 4, 12, 9 }
+  { 10, 13, 15, 14, 12, 4, 9 }
 };
 HPRef_Struct reftet_1f_1e_2va =
 {

libsrc/meshing/hprefinement.cpp

@@ -564,6 +564,13 @@ namespace netgen
       case HP_HEX_1FA_1FB_0E_0V: 
 	hps = &refhex_1fa_1fb_0e_0v; break; 
 
+      case HP_HEX7:
+	hps = &refhex7; break;
+      case HP_HEX7_1FA:
+	hps = &refhex7_1fa; break;
+      case HP_HEX7_1FB:
+	hps = &refhex7_1fb; break;
+       
 
 
       default:
@@ -740,6 +747,7 @@ namespace netgen
 	  case HP_PYRAMID: oldnp = 5; break;
 	  case HP_PRISM: oldnp = 6; break;
 	  case HP_HEX: oldnp = 8; break;
+	  case HP_HEX7: oldnp = 7; break;
             
           default:
             cerr << "HPRefElement: illegal type (3) " << hprs->geom << endl;
@@ -753,6 +761,7 @@ namespace netgen
 	    el.type == HP_TET ||
 	    el.type == HP_PRISM ||
 	    el.type == HP_HEX || 
+	    el.type == HP_HEX7 || 
 	    el.type == HP_PYRAMID)
 	  newlevel = el.levelx;
 
@@ -859,6 +868,7 @@ namespace netgen
                 newel.type == HP_TET ||
                 newel.type == HP_PRISM ||
                 newel.type == HP_HEX || 
+                newel.type == HP_HEX7 || 
                 newel.type == HP_PYRAMID)
               newel.levelx = newel.levely = newel.levelz = newlevel;
             else
@@ -870,6 +880,7 @@ namespace netgen
 	      case HP_QUAD:    newel.np=4; break;
 	      case HP_TRIG:    newel.np=3; break;
 	      case HP_HEX:     newel.np=8; break; 
+	      case HP_HEX7:    newel.np=7; break; 
 	      case HP_PRISM:   newel.np=6; break;
 	      case HP_TET:     newel.np=4; break; 
 	      case HP_PYRAMID: newel.np=5; break; 
@@ -1456,6 +1467,7 @@ namespace netgen
 		    break; 
 		  } 
 		case HP_HEX:
+		case HP_HEX7:
 		case HP_TET:
 		case HP_PRISM:
 		case HP_PYRAMID:
@@ -1516,6 +1528,9 @@ namespace netgen
 		ord_dir[2] = 2; 
 		ned = 9; 
 		break;
+	      case HEX7:
+                // ????
+                break;
 	      case HEX: 
 		/* cout << " HEX " ; 
 		for(int k=0;k<8;k++) cout << el[k] << "\t" ; 
@@ -1922,6 +1937,12 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HAS
 					face_edges, surf_edges, facepoint); 	    	    
 	 
 	    
+	      break;
+	    }
+	  case HP_HEX7:
+	    { 
+	      hpel.type = ClassifyHex7(hpel, edges, edgepoint_dom, cornerpoint, edgepoint, faces,
+                                       face_edges, surf_edges, facepoint); 	    	    
 	      break; 
 	    } 
 	  case HP_HEX:

libsrc/meshing/hprefinement.hpp

@@ -246,7 +246,11 @@ enum HPREF_ELEMENT_TYPE {
   HP_HEX_1E_0V,
   HP_HEX_3E_0V,
   HP_HEX_1F_0E_0V,
-  HP_HEX_1FA_1FB_0E_0V
+  HP_HEX_1FA_1FB_0E_0V,
+
+  HP_HEX7 = 3100,
+  HP_HEX7_1FA,      // singular quad face 1,2,3,4
+  HP_HEX7_1FB       // singular trig face 5,6,7
 };
 
 

libsrc/meshing/meshtype.cpp

@@ -1057,6 +1057,7 @@ namespace netgen
       case 4: typ = TET; break;
       case 5: typ = PYRAMID; break;
       case 6: typ = PRISM; break;
+      case 7: typ = HEX7; break;
       case 8: typ = HEX; break;
       case 10: typ = TET10; break;
       case 13: typ = PYRAMID13; break;
@@ -1139,6 +1140,7 @@ namespace netgen
       case 4: typ = TET; break;
       case 5: typ = PYRAMID; break;
       case 6: typ = PRISM; break;
+      case 7: typ = HEX7; break;
       case 8: typ = HEX; break;
       case 10: typ = TET10; break;
       case 13: typ = PYRAMID13; break;
@@ -1160,6 +1162,7 @@ namespace netgen
       case TET: np = 4; break;
       case PYRAMID: np = 5; break;
       case PRISM: np = 6; break;
+      case HEX7: np = 7; break;
       case HEX: np = 8; break;
       case TET10: np = 10; break;
       case PYRAMID13: np = 13; break;
@@ -1256,6 +1259,16 @@ namespace netgen
         { 4, 3, 1, 4, 6 }
       };
     
+    static const int hex7faces[][5] =
+      {
+        { 4, 4, 3, 2, 1 },
+        { 4, 3, 7, 6, 2 },
+        { 3, 7, 5, 6 },
+        { 4, 7, 4, 1, 5 },
+        { 4, 1, 2, 6, 5 },
+        { 3, 3, 4, 7 }
+      };
+    
     static const int hexfaces[][5] =
       {
         { 4, 4, 3, 2, 1 },
@@ -1266,7 +1279,6 @@ namespace netgen
         { 4, 3, 4, 8, 7 }
       };
 
-
     switch (np)
       {
       case 4: // tet
@@ -1301,6 +1313,14 @@ namespace netgen
             face.PNum(j) = PNum(prismfaces[i-1][j]);
           break;
         }
+      case 7: // hex7
+        {
+          //	face.SetNP(prismfaces[i-1][0]);
+          face.SetType ( ((i == 3) || (i==6)) ? TRIG : QUAD);
+          for (int j = 1; j <= face.GetNP(); j++)
+            face.PNum(j) = PNum(hex7faces[i-1][j]);
+          break;
+        }
       case 8:
         {
           face.SetType(QUAD);
@@ -1684,6 +1704,20 @@ namespace netgen
         { 6, 1, 4 }
       };
   
+    static int hex7trigs[][3] = 
+      {
+        { 1, 3, 2 },
+        { 1, 4, 3 }, 
+        { 5, 6, 7 },
+        { 1, 2, 6 },
+        { 1, 6, 5 },
+        { 2, 3, 7 },
+        { 2, 7, 6 },
+        { 3, 4, 7 },
+        { 4, 1, 7 },
+        { 1, 5, 7 }
+      };
+    
     static int hextrigs[][3] = 
       {
         { 1, 3, 2 },
@@ -1700,6 +1734,8 @@ namespace netgen
         { 1, 5, 8 }
       };
 
+
+    
     int j;
 
     int nf;
@@ -1732,6 +1768,12 @@ namespace netgen
           fp = tet10trigs;
           break;
         }
+      case HEX7:
+        {
+          nf = 10;
+          fp = hex7trigs;
+          break;
+        }
       case HEX:
         {
           nf = 12;

libsrc/meshing/meshtype.hpp

@@ -29,7 +29,7 @@ namespace netgen
     TRIG = 10, QUAD=11, TRIG6 = 12, QUAD6 = 13, QUAD8 = 14,
     TET = 20, TET10 = 21, 
     PYRAMID = 22, PRISM = 23, PRISM12 = 24, PRISM15 = 27, PYRAMID13 = 28,
-    HEX = 25, HEX20 = 26
+    HEX = 25, HEX20 = 26, HEX7 = 29
   };
 
   /*
@@ -796,7 +796,7 @@ namespace netgen
     ///
     uint8_t GetNV() const
     {
-      __assume(typ >= TET && typ <= PYRAMID13);
+      // __assume(typ >= TET && typ <= PYRAMID13);
       switch (typ)
 	{
         case TET: 
@@ -809,6 +809,8 @@ namespace netgen
 	case PYRAMID:
         case PYRAMID13:
 	  return 5;
+	case HEX7:
+	  return 7;
 	case HEX:
 	case HEX20:
 	  return 8;
@@ -915,6 +917,7 @@ namespace netgen
 	case PRISM:
         case PRISM15:
 	case PRISM12: return 5;
+        case HEX7: return 6;
         case HEX: case HEX20:
           return 6;
 	default:

libsrc/meshing/topology.hpp

@@ -262,6 +262,9 @@ inline short int MeshTopology :: GetNVertices (ELEMENT_TYPE et)
     case PRISM15:
       return 6;
 
+    case HEX7:
+      return 7;
+      
     case HEX:
     case HEX20:
       return 8;
@@ -311,6 +314,9 @@ inline short int MeshTopology :: GetNPoints (ELEMENT_TYPE et)
     case PRISM15:
       return 15;
 
+    case HEX7:
+      return 7;
+
     case HEX:
       return 8;
 
@@ -326,7 +332,7 @@ inline short int MeshTopology :: GetNPoints (ELEMENT_TYPE et)
 
 inline short int MeshTopology :: GetNEdges (ELEMENT_TYPE et)
 {
-  __assume(et >= SEGMENT && et <= PYRAMID13);
+  // __assume(et >= SEGMENT && et <= PYRAMID13);
   switch (et)
     {
     case SEGMENT:
@@ -355,11 +361,12 @@ inline short int MeshTopology :: GetNEdges (ELEMENT_TYPE et)
     case PRISM15:
       return 9;
 
+    case HEX7:
+      return 11;
+      
     case HEX:
     case HEX20:
       return 12;
-    default:
-      return 0;
       // default:
       // cerr << "Ng_ME_GetNEdges, illegal element type " << et << endl;
     }
@@ -369,7 +376,7 @@ inline short int MeshTopology :: GetNEdges (ELEMENT_TYPE et)
 
 inline short int MeshTopology :: GetNFaces (ELEMENT_TYPE et)
 {
-  __assume(et >= SEGMENT && et <= PYRAMID13);
+  // __assume(et >= SEGMENT && et <= PYRAMID13);
   switch (et)
     {
     case SEGMENT:
@@ -400,6 +407,7 @@ inline short int MeshTopology :: GetNFaces (ELEMENT_TYPE et)
 
     case HEX:
     case HEX20:
+    case HEX7:      
       return 6;
 
     default:
@@ -460,6 +468,21 @@ const ELEMENT_EDGE * MeshTopology :: GetEdges1 (ELEMENT_TYPE et)
       { 3, 5 },
       { 4, 5 }};
 
+  static ELEMENT_EDGE hex7_edges[11] =
+    {
+      { 1, 2 },
+      { 3, 4 },
+      { 4, 1 },
+      { 2, 3 },
+      { 5, 6 },
+      { 7, 5 },
+      { 6, 7 },
+      { 1, 5 },
+      { 2, 6 },
+      { 3, 7 },
+      { 4, 7 },
+    };
+
   static ELEMENT_EDGE hex_edges[12] =
     {
       { 1, 2 },
@@ -476,6 +499,7 @@ const ELEMENT_EDGE * MeshTopology :: GetEdges1 (ELEMENT_TYPE et)
       { 4, 8 },
     };
 
+  
   switch (et)
     {
     case SEGMENT:
@@ -504,6 +528,9 @@ const ELEMENT_EDGE * MeshTopology :: GetEdges1 (ELEMENT_TYPE et)
     case PRISM15:
       return prism_edges;
 
+    case HEX7:
+      return hex7_edges;
+      
     case HEX:
     case HEX20:
       return hex_edges;
@@ -561,6 +588,21 @@ const ELEMENT_EDGE * MeshTopology :: GetEdges0 (ELEMENT_TYPE et)
       { 2, 4 },
       { 3, 4 }};
 
+  static ELEMENT_EDGE hex7_edges[11] =
+    {
+      { 0, 1 },
+      { 2, 3 },
+      { 3, 0 },
+      { 1, 2 },
+      { 4, 5 },
+      { 6, 4 },
+      { 5, 6 },
+      { 0, 4 },
+      { 1, 5 },
+      { 2, 6 },
+      { 3, 6 },
+    };
+
   static ELEMENT_EDGE hex_edges[12] =
     {
       { 0, 1 },
@@ -577,6 +619,7 @@ const ELEMENT_EDGE * MeshTopology :: GetEdges0 (ELEMENT_TYPE et)
       { 3, 7 },
     };
 
+  
   switch (et)
     {
     case SEGMENT:
@@ -605,6 +648,9 @@ const ELEMENT_EDGE * MeshTopology :: GetEdges0 (ELEMENT_TYPE et)
     case PRISM15:
       return prism_edges;
 
+    case HEX7:
+      return hex7_edges;
+      
     case HEX:
     case HEX20:
       return hex_edges;
@@ -661,6 +707,21 @@ FlatArray<ELEMENT_EDGE> MeshTopology :: GetEdges (ELEMENT_TYPE et)
       { 2, 4 },
       { 3, 4 }};
 
+  static ELEMENT_EDGE hex7_edges[11] =
+    {
+      { 0, 1 },
+      { 2, 3 },
+      { 3, 0 },
+      { 1, 2 },
+      { 4, 5 },
+      { 6, 4 },
+      { 5, 6 },
+      { 0, 4 },
+      { 1, 5 },
+      { 2, 6 },
+      { 3, 6 },
+    };
+
   static ELEMENT_EDGE hex_edges[12] =
     {
       { 0, 1 },
@@ -705,6 +766,9 @@ FlatArray<ELEMENT_EDGE> MeshTopology :: GetEdges (ELEMENT_TYPE et)
     case PRISM15:
       return { 9, prism_edges };
 
+    case HEX7:
+      return { 11, hex7_edges };
+
     case HEX:
     case HEX20:
       return { 12, hex_edges };
@@ -752,6 +816,17 @@ inline const ELEMENT_FACE * MeshTopology :: GetFaces1 (ELEMENT_TYPE et)
       { 1, 4, 3, 2 } 
     };
 
+  static const ELEMENT_FACE hex7_faces[6] =
+    {
+      { 1, 4, 3, 2 },
+      { 5, 6, 7, 0  },
+      { 1, 2, 6, 5 },
+      { 2, 3, 7, 6 },
+      { 3, 4, 7, 0 },
+      { 4, 1, 5, 7 }
+    };
+
+  
   static const ELEMENT_FACE hex_faces[6] =
     {
       { 1, 4, 3, 2 },
@@ -792,6 +867,9 @@ inline const ELEMENT_FACE * MeshTopology :: GetFaces1 (ELEMENT_TYPE et)
     case SEGMENT:
     case SEGMENT3:
 
+    case HEX7:
+      return hex7_faces;
+    
     case HEX:
     case HEX20:
       return hex_faces;
@@ -837,6 +915,16 @@ inline const ELEMENT_FACE * MeshTopology :: GetFaces0 (ELEMENT_TYPE et)
       { 0, 3, 2, 1 } 
     };
 
+  static const ELEMENT_FACE hex7_faces[6] =
+    {
+      { 0, 3, 2, 1 },
+      { 4, 5, 6, -1},
+      { 0, 1, 5, 4 },
+      { 1, 2, 6, 5 },
+      { 2, 3, 6, -1},
+      { 3, 0, 4, 6 }
+    };
+
   static const ELEMENT_FACE hex_faces[6] =
     {
       { 0, 3, 2, 1 },
@@ -877,6 +965,9 @@ inline const ELEMENT_FACE * MeshTopology :: GetFaces0 (ELEMENT_TYPE et)
     case SEGMENT:
     case SEGMENT3:
 
+    case HEX7:
+      return hex7_faces;
+
     case HEX:
     case HEX20:
       return hex_faces;

libsrc/visualization/vsmesh.cpp

@@ -2672,6 +2672,119 @@ namespace netgen
 	      }
 	  }
       }
+
+
+    for (ElementIndex ei = 0; ei < mesh->GetNE(); ei++)
+      {
+	const Element & el = (*mesh)[ei];
+	if (el.GetType() == HEX7 && !el.IsDeleted())
+	  {
+            /*
+            CurvedElements & curv = mesh->GetCurvedElements();
+            if (curv.IsHighOrder()) 
+	      {
+		const ELEMENT_FACE * faces = MeshTopology :: GetFaces1 (HEX);
+		const Point3d * vertices = MeshTopology :: GetVertices (HEX);
+
+		Point<3> grid[11][11];
+		Point<3> fpts[4];
+		int order = subdivisions+1;
+
+		for (int quad = 0; quad<6; quad++)
+		  {
+		    for (int j = 0; j < 4; j++)
+		      fpts[j] = vertices[faces[quad][j]-1];
+
+		    static Point<3> c(0.5, 0.5, 0.5);
+		    if (vispar.shrink < 1)
+		      for (int j = 0; j < 4; j++)
+                        fpts[j] += (1-vispar.shrink) * (c-fpts[j]);
+
+		    for (int ix = 0; ix <= order; ix++)
+		      for (int iy = 0; iy <= order; iy++)
+			{
+			  double lami[4] =
+			    { (1-double(ix)/order) * (1-double(iy)/order),
+			      (  double(ix)/order) * (1-double(iy)/order),
+			      (  double(ix)/order) * (  double(iy)/order),
+			      (1-double(ix)/order) * (  double(iy)/order) };
+
+			  Point<3> xl;
+			  for (int l = 0; l < 3; l++)
+			    xl(l) = lami[0] * fpts[0](l) + lami[1] * fpts[1](l) +
+			      lami[2] * fpts[2](l) + lami[3] * fpts[3](l);
+
+			  curv.CalcElementTransformation (xl, ei, grid[ix][iy]);
+			}
+
+		    for (int j = 0; j <= order; j++)
+		      ToBernstein (order, &grid[j][0], &grid[0][1]-&grid[0][0]);
+		    for (int j = 0; j <= order; j++)
+		      ToBernstein (order, &grid[0][j], &grid[1][0]-&grid[0][0]);
+
+		    glMap2d(GL_MAP2_VERTEX_3,
+			    0.0, 1.0, &grid[0][1](0)-&grid[0][0](0), order+1,
+			    0.0, 1.0, &grid[1][0](0)-&grid[0][0](0), order+1,
+			    &grid[0][0](0));
+		    glEnable(GL_MAP2_VERTEX_3);
+		    glEnable(GL_AUTO_NORMAL);
+
+		    glMapGrid2f(8, 0.0, 1.0, 8, 0.0, 1.0);
+		    glEvalMesh2(GL_FILL, 0, 8, 0, 8);
+
+		    glDisable (GL_AUTO_NORMAL);
+		    glDisable (GL_MAP2_VERTEX_3);
+		  }
+	      }
+            else
+            */
+	      {
+		Point3d c(0,0,0);
+		if (vispar.shrink < 1)
+		  {
+		    for (int j = 1; j <= 7; j++)
+		      {
+			Point3d p = mesh->Point(el.PNum(j));
+			c.X() += p.X();
+			c.Y() += p.Y();
+			c.Z() += p.Z();
+		      }
+		    c.X() /= 7;
+		    c.Y() /= 7;
+		    c.Z() /= 7;
+		  }
+
+		glBegin (GL_TRIANGLES);
+
+		el.GetSurfaceTriangles (faces);
+		for (int j = 1; j <= faces.Size(); j++)
+		  {
+		    Element2d & face = faces.Elem(j);
+		    Point<3> lp1 = mesh->Point (el.PNum(face.PNum(1)));
+		    Point<3> lp2 = mesh->Point (el.PNum(face.PNum(2)));
+		    Point<3> lp3 = mesh->Point (el.PNum(face.PNum(3)));
+		    Vec<3> n = Cross (lp3-lp1, lp2-lp1);
+		    n.Normalize();
+		    glNormal3dv (n);
+
+		    if (vispar.shrink < 1)
+		      {
+			lp1 = c + vispar.shrink * (lp1 - c);
+			lp2 = c + vispar.shrink * (lp2 - c);
+			lp3 = c + vispar.shrink * (lp3 - c);
+		      }
+
+		    glVertex3dv (lp1);
+		    glVertex3dv (lp2);
+		    glVertex3dv (lp3);
+		  }
+
+		glEnd();
+	      }
+	  }
+      }
+
+    
     glEndList ();
   }