parallel topology + curved elements

libsrc/csg/csgeom.cpp

@@ -1497,4 +1497,80 @@ namespace netgen
 		 boundingbox.PMin()(2));
     return max2 (maxs, -mins) * 1.1;
   }
+
+
+
+  class CSGeometryRegister : public GeometryRegister
+  {
+  public:
+    virtual NetgenGeometry * Load (string filename) const;
+    virtual NetgenGeometry * LoadFromMeshFile (istream & ist) const;
+    // virtual VisualScene * GetVisualScene (const NetgenGeometry * geom) const;
+  };
+
+  extern CSGeometry * ParseCSG (istream & istr);
+
+  NetgenGeometry *  CSGeometryRegister :: Load (string filename) const
+  {
+    const char * cfilename = filename.c_str();
+    if (strcmp (&cfilename[strlen(cfilename)-3], "geo") == 0)
+      {
+	PrintMessage (1, "Load CSG geometry file ", cfilename);
+
+	ifstream infile(cfilename);
+
+	CSGeometry * hgeom = ParseCSG (infile);
+	if (!hgeom)
+	  throw NgException ("geo-file should start with 'algebraic3d'");
+
+	hgeom -> FindIdenticSurfaces(1e-8 * hgeom->MaxSize()); 
+	return hgeom;
+      }
+
+    if (strcmp (&cfilename[strlen(cfilename)-3], "ngg") == 0)
+      {
+	PrintMessage (1, "Load new CSG geometry file ", cfilename);
+
+	ifstream infile(cfilename);
+	CSGeometry * hgeom = new CSGeometry("");
+	hgeom -> Load (infile);
+
+	return hgeom;
+      }
+
+
+    
+    return NULL;
+  }
+
+  NetgenGeometry * CSGeometryRegister :: LoadFromMeshFile (istream & ist) const 
+  {
+    string auxstring;
+    if (ist.good())
+      {
+	ist >> auxstring;
+	if (auxstring == "csgsurfaces")
+	  {
+	    CSGeometry * geometry = new CSGeometry ("");
+	    geometry -> LoadSurfaces(ist);
+	    return geometry;
+	  }
+	// else
+	// ist.putback (auxstring);
+      }
+    return NULL;
+  }
+
+
+  class CSGInit
+  {
+  public:
+    CSGInit()
+    {
+      geometryregister.Append (new CSGeometryRegister);
+    }
+  };
+
+  CSGInit csginit;
+
 }

libsrc/csg/csgpkg.cpp

@@ -546,6 +546,7 @@ namespace netgen
   }
 
 
+  /*
   class CSGeometryRegister : public GeometryRegister
   {
   public:
@@ -617,8 +618,27 @@ namespace netgen
       }
     return NULL;
   }
+  */
 
 
+
+  class CSGeometryVisRegister : public GeometryRegister
+  {
+  public:
+    virtual NetgenGeometry * Load (string filename) const { return NULL; }
+    virtual VisualScene * GetVisualScene (const NetgenGeometry * geom) const;
+  };
+
+  VisualScene * CSGeometryVisRegister :: GetVisualScene (const NetgenGeometry * geom) const
+  {
+    CSGeometry * geometry = dynamic_cast<CSGeometry*> (ng_geometry);
+    if (geometry)
+      {
+	vsgeom.SetGeometry (geometry);
+	return &vsgeom;
+      }
+    return NULL;
+  }
 }
 
 
@@ -626,10 +646,10 @@ using namespace netgen;
 
 int Ng_CSG_Init (Tcl_Interp * interp)
 {
-  geometryregister.Append (new CSGeometryRegister);
+  geometryregister.Append (new CSGeometryVisRegister);
+  if (interp == NULL) return TCL_OK;
   
 
-
   Tcl_CreateCommand (interp, "Ng_ParseGeometry", Ng_ParseGeometry,
 		     (ClientData)NULL,
 		     (Tcl_CmdDeleteProc*) NULL);

libsrc/general/dynamicmem.hpp

@@ -27,8 +27,8 @@ protected:
   void Alloc (size_t s);
   void ReAlloc (size_t s);
   void Free ();
-  char * Ptr() { return ptr; }
+  // char * Ptr() { return ptr; }
-  const char * Ptr() const { return ptr; }
+  char * Ptr() const { return ptr; }
   void Swap (BaseDynamicMem & m2);
 public:
   void SetName (const char * aname);
@@ -64,22 +64,24 @@ public:
     BaseDynamicMem::Free ();
   }
 
+  /*
   const T * Ptr() const
   {
     return reinterpret_cast<const T*> (BaseDynamicMem::Ptr());
   }
-
+  */
-  T * Ptr()
+  const T * Ptr()
   {
     return reinterpret_cast<T*> (BaseDynamicMem::Ptr());
   }
 
+  /*
   operator const T* () const
   {
     return reinterpret_cast<const T*> (BaseDynamicMem::Ptr());
   }
-
+  */
-  operator T* () 
+  operator T* () const
   {
     return reinterpret_cast<T*> (BaseDynamicMem::Ptr());
   }

libsrc/general/mpi_interface.hpp

@@ -186,8 +186,8 @@ namespace netgen
     receive-table entries will be set
    */
   template <typename T>
-  inline void MyMPI_ExchangeTable (TABLE<T> & send_verts, 
+  inline void MyMPI_ExchangeTable (TABLE<T> & send_data, 
-				   TABLE<T> & recv_verts, int tag,
+				   TABLE<T> & recv_data, int tag,
 				   MPI_Comm comm = MPI_COMM_WORLD)
   {
     int ntasks, rank;
@@ -197,7 +197,7 @@ namespace netgen
     Array<MPI_Request> requests;
     for (int dest = 0; dest < ntasks; dest++)
       if (dest != rank)
-	requests.Append (MyMPI_ISend (send_verts[dest], dest, tag, comm));
+	requests.Append (MyMPI_ISend (send_data[dest], dest, tag, comm));
 
     for (int i = 0; i < ntasks-1; i++)
       {
@@ -205,9 +205,10 @@ namespace netgen
 	MPI_Probe (MPI_ANY_SOURCE, tag, comm, &status);
 	int size, src = status.MPI_SOURCE;
 	MPI_Get_count (&status, MPI_INT, &size);
-	recv_verts.SetEntrySize (src, size, sizeof(int));
+	recv_data.SetEntrySize (src, size, sizeof(int));
-	requests.Append (MyMPI_IRecv (recv_verts[src], src, tag, comm));
+	requests.Append (MyMPI_IRecv (recv_data[src], src, tag, comm));
       }
+    MPI_Barrier (comm);
     MPI_Waitall (requests.Size(), &requests[0], MPI_STATUS_IGNORE);
   }
 

libsrc/general/template.hpp

@@ -188,6 +188,14 @@ public:
 };
 
 
+inline INDEX_2 Sort (const INDEX_2 & i2)
+{
+  INDEX_2 tmp = i2;
+  tmp.Sort();
+  return tmp;
+}
+
+
 ///
 class INDEX_3
 {

libsrc/include/mystdlib.h

@@ -1,7 +1,7 @@
 #ifndef FILE_MYSTDLIB
 #define FILE_MYSTDLIB
 
-#ifndef WIN32
+#ifdef HAVE_CONFIG_H
 #include <config.h>
 #endif
 

libsrc/interface/nginterface.cpp

@@ -119,9 +119,9 @@ namespace netgen
 
 
 #ifdef OPENGL 
-  extern VisualSceneSolution vssolution;
+  // extern VisualSceneSolution vssolution;
 #endif
-  extern CSGeometry * ParseCSG (istream & istr);
+  // extern CSGeometry * ParseCSG (istream & istr);
 
 #ifdef SOCKETS
   extern AutoPtr<ClientSocket> clientsocket;
@@ -161,8 +161,8 @@ void Ng_LoadGeometry (const char * filename)
       return;
     }
 
-  if (id == 0)
+  // if (id == 0)
-    cerr << "cannot load geometry '" << filename << "'" << endl;
+  cerr << "cannot load geometry '" << filename << "'" << ", id = " << id << endl;
 }                          
 
 
@@ -1182,41 +1182,12 @@ void Ng_HPRefinement (int levels, double parameter, bool setorders,
 void Ng_HighOrder (int order, bool rational)
 {
   NgLock meshlock (mesh->MajorMutex(), true);
-  /*
-    Refinement * ref;
 
-    if (stlgeometry)
+  mesh -> GetCurvedElements().BuildCurvedElements 
-    ref = new RefinementSTLGeometry (*stlgeometry);
+    (&const_cast<Refinement&> (ng_geometry -> GetRefinement()),
-    #ifdef OCCGEOMETRY
+     order, rational);
-    else if (occgeometry)
-    ref = new OCCRefinementSurfaces (*occgeometry);
-    #endif
-    #ifdef ACIS
-    else if (acisgeometry)
-    {
-    ref = new ACISRefinementSurfaces (*acisgeometry);
-    }
-    #endif
-    else if (geometry2d)
-    ref = new Refinement2d (*geometry2d);
-    else
-    {
-    ref = new RefinementSurfaces (*geometry);
-    }
-  */
-  // cout << "parameter 1: " << argv[1] << " (conversion to int = " << atoi(argv[1]) << ")" << endl;
- 
 
-  mesh -> GetCurvedElements().BuildCurvedElements (&const_cast<Refinement&> (ng_geometry -> GetRefinement()),
-						   order, rational);
   mesh -> SetNextMajorTimeStamp();
-
-  /*
-    if(mesh)
-    mesh -> GetCurvedElements().BuildCurvedElements (ref, order, rational);
-  */
-
-  // delete ref;
 }
 
 
@@ -1690,13 +1661,18 @@ int Ng_GetNVertexElements (int vnr)
 
 void Ng_GetVertexElements (int vnr, int * els)
 {
-  FlatArray<int> ia(0,0);
   if (mesh->GetDimension() == 3)
-    ia = mesh->GetTopology().GetVertexElements(vnr);
+    {
+      FlatArray<ElementIndex> ia = mesh->GetTopology().GetVertexElements(vnr);
+      for (int i = 0; i < ia.Size(); i++)
+	els[i] = ia[i]+1;
+    }
   else
-    ia = mesh->GetTopology().GetVertexSurfaceElements(vnr);
+    {
-  for (int i = 0; i < ia.Size(); i++)
+      FlatArray<int> ia = mesh->GetTopology().GetVertexSurfaceElements(vnr);
-    els[i] = ia[i];
+      for (int i = 0; i < ia.Size(); i++)
+	els[i] = ia[i];
+    }
 }
 
 
@@ -1963,11 +1939,11 @@ int Ng_IsRunning()
 int Ng_GetVertex_Elements( int vnr, int* elems )
 {
   const MeshTopology& topology = mesh->GetTopology();
-  ArrayMem<int,4> indexArray;
+  ArrayMem<ElementIndex,4> indexArray;
   topology.GetVertexElements( vnr, indexArray );
   
   for( int i=0; i<indexArray.Size(); i++ )
-    elems[i] = indexArray[i];
+    elems[i] = indexArray[i]+1;
   
   return indexArray.Size();
 }
@@ -1989,7 +1965,7 @@ int Ng_GetVertex_SurfaceElements( int vnr, int* elems )
 int Ng_GetVertex_NElements( int vnr )
 {
   const MeshTopology& topology = mesh->GetTopology();
-  ArrayMem<int,4> indexArray;
+  ArrayMem<ElementIndex,4> indexArray;
   topology.GetVertexElements( vnr, indexArray );
   
   return indexArray.Size();

libsrc/interface/writediffpack.cpp

@@ -28,6 +28,8 @@ void WriteDiffPackFormat (const Mesh & mesh,
 
   ofstream outfile(filename.c_str());
 
+
+
   if (mesh.GetDimension() == 3)
 
     {
@@ -80,6 +82,17 @@ void WriteDiffPackFormat (const Mesh & mesh,
 	"   - the boundary indicators that are set (ON) if any.\n"
 	"#\n";
 
+
+      // setup point-to-surfaceelement table 
+      TABLE<SurfaceElementIndex, PointIndex::BASE> point2sel(np);
+      for (SurfaceElementIndex sei = 0; sei < nse; sei++)
+	{
+	  const Element2d & el = mesh[sei];
+	  for (int j = 0; j < el.GetNP(); j++)
+	    point2sel.Add (el[j], sei);
+	}
+
+
       for (i = 1; i <= np; i++)
         {
           const Point3d & p = mesh.Point(i);
@@ -96,13 +109,17 @@ void WriteDiffPackFormat (const Mesh & mesh,
 	  if(mesh[PointIndex(i)].Type() != INNERPOINT) 
 	    {
 	      BCsinpoint.DeleteAll();
+	      /*
 	      for (j = 1; j <= nse; j++) 
+	      */
+	      FlatArray<SurfaceElementIndex> sels = point2sel[i];
+	      for (int jj = 0; jj < sels.Size(); jj++)
 		{
-		  for (k = 1; k <= mesh.SurfaceElement(j).GetNP(); k++) 
+		  for (k = 1; k <= mesh[sels[jj]].GetNP(); k++) 
 		    {
-		      if(mesh.SurfaceElement(j).PNum(k)==i) 
+		      if(mesh[sels[jj]].PNum(k)==i) 
 			{
-			  int BC=mesh.GetFaceDescriptor(mesh.SurfaceElement(j).GetIndex()).BCProperty();
+			  int BC=mesh.GetFaceDescriptor(mesh[sels[jj]].GetIndex()).BCProperty();
 			  int nbcsp=BCsinpoint.Size();
 			  int found = 0;
 			  for (l = 1; l <= nbcsp; l++)

libsrc/meshing/boundarylayer.cpp

@@ -546,7 +546,7 @@ namespace netgen
          
          for (i = 1; i <= np; i++)
          {
-            Array<int> vertelems;
+            Array<ElementIndex> vertelems;
 
             if(bndnodes.Test(i))
             {

libsrc/meshing/meshclass.cpp

@@ -177,7 +177,7 @@ namespace netgen
     return pi;
   }
 
-
+  /*
 #ifdef PARALLEL
   PointIndex Mesh :: AddPoint (const Point3d & p, bool isghost,  int layer)
   { 
@@ -210,7 +210,7 @@ namespace netgen
   }
 
 #endif
-
+  */
 
 
   SegmentIndex Mesh :: AddSegment (const Segment & s)
@@ -462,26 +462,14 @@ namespace netgen
 
     for (ElementIndex ei = 0; ei < GetNE(); ei++)
       {
-        outfile.width(8);
         outfile << (*this)[ei].GetIndex();
-        outfile.width(8);
+        outfile << " " << (*this)[ei].GetNP();
-        outfile << (*this)[ei].GetNP();
 
         Element el = (*this)[ei];
-        if (inverttets)
+        if (inverttets) el.Invert();
-          el.Invert();
-
-        /*
-          for (j = 0; j < el.GetNP(); j++)
-          for (int k = 0; k < el.GetNP()-1; k++)
-          if (el[k] > el[k+1]) swap (el[k], el[k+1]);
-        */
 
         for (j = 0; j < el.GetNP(); j++)
-          {
+	  outfile << " " << el[j];
-            outfile.width(8);
-            outfile << el[j];
-          }
         outfile << "\n";
       }
 
@@ -920,7 +908,7 @@ namespace netgen
                 if (inverttets)
                   el.Invert();
 
-                AddVolumeElement (el);
+		AddVolumeElement (el);
               }
           }
 
@@ -1193,11 +1181,12 @@ namespace netgen
 
 
     CalcSurfacesOfNode ();
-    //  BuildConnectedNodes ();
+ 
-    topology -> Update();
+    if (ntasks == 1) // sequential run only
-    clusters -> Update();
+      {
-
+	topology -> Update();
-    // Distribute();
+	clusters -> Update();
+      }
 
     SetNextMajorTimeStamp();
     //  PrintMemInfo (cout);

libsrc/meshing/meshclass.hpp

@@ -216,10 +216,12 @@ namespace netgen
 
     DLL_HEADER PointIndex AddPoint (const Point3d & p, int layer = 1);
     DLL_HEADER PointIndex AddPoint (const Point3d & p, int layer, POINTTYPE type);
+    /*
 #ifdef PARALLEL
     PointIndex AddPoint (const Point3d & p, bool aisghost, int layer = 1);
     PointIndex AddPoint (const Point3d & p, bool aisghost, int layer, POINTTYPE type);
 #endif
+    */
     int GetNP () const { return points.Size(); }
 
     MeshPoint & Point(int i) { return points.Elem(i); }
@@ -659,13 +661,6 @@ namespace netgen
   
 
 
-    /*
-   /// build connected nodes along prism stack
-   void BuildConnectedNodes ();
-   void ConnectToNodeRec (int node, int tonode, 
-   const TABLE<int> & conto);
-    */
-
     bool PureTrigMesh (int faceindex = 0) const;
     bool PureTetMesh () const;
 

libsrc/meshing/meshtype.hpp

@@ -814,6 +814,11 @@ namespace netgen
     ///
     int meshdocval;
 
+#ifdef PARALLEL
+    /// number of partition for parallel computation 
+    int partitionNumber;
+#endif
+
   private:
     string* bcname;
 
@@ -863,6 +868,13 @@ namespace netgen
   
     PointIndex & operator[] (int i) { return pnums[i]; }
     const PointIndex & operator[] (int i) const { return pnums[i]; }
+
+#ifdef PARALLEL
+    int GetPartition () const { return partitionNumber; }
+    void SetPartition (int nr) { partitionNumber = nr; }; 
+#else
+    int GetPartition () const { return 0; }
+#endif
   };
 
   ostream & operator<<(ostream  & s, const Segment & seg);

libsrc/meshing/parallelmesh.cpp

@@ -3,19 +3,21 @@
 #include <meshing.hpp>
 #include "paralleltop.hpp"
 
-#define METIS4
+// #define METIS4
 
 
 #ifdef METIS
 namespace metis {
   extern "C" {
 
-#ifdef METIS4
-#include <metis.h>
-    typedef idxtype idx_t;   
-#else
 #include <metis.h>
+
+#if METIS_VER_MAJOR >= 5
+#define METIS5
     typedef idx_t idxtype;   
+#else
+#define METIS4
+    typedef idxtype idx_t;   
 #endif
   } 
 }
@@ -23,19 +25,14 @@ namespace metis {
 using namespace metis;
 #endif
 
-
-
 namespace netgen
 {
 
   template <>
-  inline MPI_Datatype MyGetMPIType<PointIndex> ( ) 
+  inline MPI_Datatype MyGetMPIType<PointIndex> ( )
   { return MPI_INT; }
 
 
-
-
-
   void Mesh :: SendRecvMesh ()
   {
     if (id == 0)
@@ -64,9 +61,6 @@ namespace netgen
 
 
 
-
-
-
     {
       // distribute number of local elements
       if (id == 0)
@@ -94,9 +88,7 @@ namespace netgen
     else
       ReceiveParallelMesh();
 
-
     paralleltop -> UpdateCoarseGrid();
-
   }
 
 
@@ -131,6 +123,15 @@ namespace netgen
     for (SurfaceElementIndex ei = 0; ei < GetNSE(); ei++)
       sels_of_proc.Add ( (*this)[ei].GetPartition(), ei);
 
+    Array<int> num_segs_on_proc(ntasks);
+    num_segs_on_proc = 0;
+    for (SegmentIndex ei = 0; ei < GetNSeg(); ei++)
+      num_segs_on_proc[(*this)[ei].GetPartition()]++;
+
+    TABLE<SegmentIndex> segs_of_proc (num_segs_on_proc);
+    for (SegmentIndex ei = 0; ei < GetNSeg(); ei++)
+      segs_of_proc.Add ( (*this)[ei].GetPartition(), ei);
+
 
 
 
@@ -142,16 +143,10 @@ namespace netgen
     num_procs_on_vert = 0;
     vert_flag = -1;
 
-    Array<int> nelloc (ntasks);
-    nelloc = 0;
-    Array<int> nselloc (ntasks);
-    nselloc = 0;
-
 
     for (int dest = 1; dest < ntasks; dest++)
       {
 	FlatArray<ElementIndex> els = els_of_proc[dest];
-
 	for (int hi = 0; hi < els.Size(); hi++)
 	  {
 	    const Element & el = (*this) [ els[hi] ];
@@ -165,16 +160,13 @@ namespace netgen
 		    num_verts_on_proc[dest]++;
 		    num_procs_on_vert[epi]++;
 
-		    paralleltop -> SetDistantPNum (dest, epi); // , num_verts_on_proc[dest]);
+		    paralleltop -> SetDistantPNum (dest, epi); 
 		  }
 	      }
-	    nelloc[dest] ++;
-	    // paralleltop -> SetDistantEl ( dest, els[hi]+1, nelloc[dest] );
 	  }
 
 
 	FlatArray<SurfaceElementIndex> sels = sels_of_proc[dest];
-
 	for (int hi = 0; hi < sels.Size(); hi++)
 	  {
 	    const Element2d & el = (*this) [ sels[hi] ];
@@ -188,11 +180,28 @@ namespace netgen
 		    num_verts_on_proc[dest]++;
 		    num_procs_on_vert[epi]++;
 
-		    paralleltop -> SetDistantPNum ( dest, epi); //  num_verts_on_proc[dest]);
+		    paralleltop -> SetDistantPNum (dest, epi);
+		  }
+	      }
+	  }
+
+	FlatArray<SegmentIndex> segs = segs_of_proc[dest];
+	for (int hi = 0; hi < segs.Size(); hi++)
+	  {
+	    const Segment & el = (*this) [segs[hi]];
+	    for (int i = 0; i < 2; i++)
+	      {
+		PointIndex epi = el[i]; 
+		if (vert_flag[epi] < dest)
+		  {
+		    vert_flag[epi] = dest;
+
+		    num_verts_on_proc[dest]++;
+		    num_procs_on_vert[epi]++;
+
+		    paralleltop -> SetDistantPNum (dest, epi);
 		  }
 	      }
-	    nselloc[dest] ++;
-	    // paralleltop -> SetDistantSurfEl ( dest, sels[hi]+1, nselloc[dest] );
 	  }
       }
 
@@ -233,6 +242,21 @@ namespace netgen
 	      }
 	  }
 
+	FlatArray<SegmentIndex> segs = segs_of_proc[dest];
+	for (int hi = 0; hi < segs.Size(); hi++)
+	  {
+	    const Segment & el = (*this) [ segs[hi] ];
+	    for (int i = 0; i < 2; i++)
+	      {
+		PointIndex epi = el[i];
+		if (vert_flag[epi] < dest)
+		  {
+		    vert_flag[epi] = dest;
+		    procs_of_vert.Add (epi, dest);
+		  }
+	      }
+	  }
+
       }
 
     for (int vert = 1; vert <= GetNP(); vert++ )
@@ -383,8 +407,53 @@ namespace netgen
     for (int dest = 1; dest < ntasks; dest++)
       sendrequests.Append (MyMPI_ISend(selbuf[dest], dest, MPI_TAG_MESH+4));
 
+
     PrintMessage ( 3, "Sending Edge Segments");
 
+    
+    Array <int> nloc(ntasks); 
+    nloc = 0;
+    bufsize = 1;
+      
+    for (int i = 1; i <= GetNSeg(); i++ )
+      {
+	const Segment & seg = LineSegment (i);
+	int dest = seg.GetPartition();
+	nloc[dest] ++;
+	bufsize[dest] += 14;
+      }
+    
+    TABLE<double> segm_buf(bufsize);
+    
+    /*
+    for (int dest = 1; dest < ntasks; dest++ )
+      segm_buf.Add (dest, nloc[dest]);
+    */
+    for (int i = 1; i <= GetNSeg(); i++ )
+      {
+	const Segment & seg = LineSegment (i);
+	int dest = seg.GetPartition();
+	
+	segm_buf.Add (dest, i);
+	segm_buf.Add (dest, seg.si);
+	segm_buf.Add (dest, seg.pnums[0]);
+	segm_buf.Add (dest, seg.pnums[1]);
+	segm_buf.Add (dest, seg.geominfo[0].trignum);
+	segm_buf.Add (dest, seg.geominfo[1].trignum);
+	segm_buf.Add (dest, seg.surfnr1);
+	segm_buf.Add (dest, seg.surfnr2);
+	segm_buf.Add (dest, seg.edgenr);
+	segm_buf.Add (dest, seg.epgeominfo[0].dist);
+	segm_buf.Add (dest, seg.epgeominfo[1].edgenr);
+	segm_buf.Add (dest, seg.epgeominfo[1].dist);
+	segm_buf.Add (dest, seg.singedge_right);
+	segm_buf.Add (dest, seg.singedge_left);
+      }
+    
+    for (int dest = 1; dest < ntasks; dest++)
+      sendrequests.Append (MyMPI_ISend(segm_buf[dest], dest, MPI_TAG_MESH+5));
+
+    /*
     Array <int> nlocseg(ntasks), segi(ntasks);
     for ( int i = 0; i < ntasks; i++)
       {
@@ -452,6 +521,9 @@ namespace netgen
 
     for ( int dest = 1; dest < ntasks; dest++)
       sendrequests.Append (MyMPI_ISend(segmbuf[dest], dest, MPI_TAG_MESH+5));
+    */
+
+    PrintMessage ( 3, "now wait ...");
 
     MPI_Waitall (sendrequests.Size(), &sendrequests[0], MPI_STATUS_IGNORE);
 
@@ -592,6 +664,7 @@ namespace netgen
       int segi = 1;
       int nsegloc = int ( segmbuf.Size() / 14 ) ;
       paralleltop -> SetNSegm ( nsegloc );
+
       while ( ii < segmbuf.Size() )
 	{
 	  globsegi = int (segmbuf[ii++]);
@@ -681,10 +754,71 @@ namespace netgen
   }
   
 
+#ifdef METIS5
+  void Mesh :: ParallelMetis ( )  
+  {
+    PrintMessage (3, "call metis 5 ...");
+
+    int timer = NgProfiler::CreateTimer ("Mesh::Partition");
+    NgProfiler::RegionTimer reg(timer);
+
+    idx_t ne = GetNE() + GetNSE() + GetNSeg();
+    idx_t nn = GetNP();
+
+    Array<idx_t> eptr, eind;
+    for (int i = 0; i < GetNE(); i++)
+      {
+	eptr.Append (eind.Size());
+	const Element & el = VolumeElement(i+1);
+	for (int j = 0; j < el.GetNP(); j++)
+	  eind.Append (el[j]-1);
+      }
+    for (int i = 0; i < GetNSE(); i++)
+      {
+	eptr.Append (eind.Size());
+	const Element2d & el = SurfaceElement(i+1);
+	for (int j = 0; j < el.GetNP(); j++)
+	  eind.Append (el[j]-1);
+      }
+    for (int i = 0; i < GetNSeg(); i++)
+      {
+	eptr.Append (eind.Size());
+	const Segment & el = LineSegment(i+1);
+	eind.Append (el[0]);
+	eind.Append (el[1]);
+      }
+    eptr.Append (eind.Size());
+    Array<idx_t> epart(ne), npart(nn);
+
+    int nparts = ntasks-1;
+    int edgecut;
+
+    /*
+    int ncommon = 3;
+    METIS_PartMeshDual (&ne, &nn, &eptr[0], &eind[0], NULL, NULL, &ncommon, &nparts,
+			NULL, NULL,
+			&edgecut, &epart[0], &npart[0]);
+    */
+    METIS_PartMeshNodal (&ne, &nn, &eptr[0], &eind[0], NULL, NULL, &nparts,
+			 NULL, NULL,
+			 &edgecut, &epart[0], &npart[0]);
+
+    PrintMessage (3, "metis complete");
+    // cout << "done" << endl;
+
+    for (int i = 0; i < GetNE(); i++)
+      VolumeElement(i+1).SetPartition(epart[i] + 1);
+    for (int i = 0; i < GetNSE(); i++)
+      SurfaceElement(i+1).SetPartition(epart[i+GetNE()] + 1);
+    for (int i = 0; i < GetNSeg(); i++)
+      LineSegment(i+1).SetPartition(epart[i+GetNE()+GetNSE()] + 1);
+  }
+#endif
 
 
 
-#ifdef METIS
+
+#ifdef METIS4
   void Mesh :: ParallelMetis ( )  
   {
     int timer = NgProfiler::CreateTimer ("Mesh::Partition");
@@ -783,7 +917,7 @@ namespace netgen
 			    &edgecut, &epart[0], &npart[0]);
 #else
 	cout << "call metis(5)_PartMeshDual ... " << endl;
-	idx_t options[METIS_NOPTIONS];
+	// idx_t options[METIS_NOPTIONS];
 	
 	Array<idx_t> eptr(ne+1);
 	for (int j = 0; j < ne+1; j++)

libsrc/meshing/paralleltop.cpp

@@ -113,7 +113,7 @@ namespace netgen
   void ParallelMeshTopology :: UpdateCoarseGridGlobal ()
   {
     if (id == 0)
-      PrintMessage ( 3, "UPDATE GLOBAL COARSEGRID STARTS" );      // JS
+      PrintMessage ( 3, "UPDATE GLOBAL COARSEGRID STARTS" );      
 
     int timer = NgProfiler::CreateTimer ("UpdateCoarseGridGlobal");
     NgProfiler::RegionTimer reg(timer);
@@ -151,6 +151,7 @@ namespace netgen
 
 	    for ( int i = 0; i < edges.Size(); i++ )
 	      sendarray.Append (edges[i]);
+	    sendarray.Append (topology.GetSurfaceElementFace (el));
 	  }
 
 	Array<MPI_Request> sendrequests;
@@ -188,9 +189,9 @@ namespace netgen
 	    topology.GetSurfaceElementEdges (surfel, edges);
 	    for (int i = 0; i  < edges.Size(); i++)
 	      SetLoc2Glob_Edge (edges[i], recvarray[ii++]);
+	    int face = topology.GetSurfaceElementFace (surfel);
+	    SetLoc2Glob_Face ( face, recvarray[ii++]);
 	  }
-
-
       }
     
     is_updated = true;
@@ -213,13 +214,7 @@ namespace netgen
       PrintMessage (1, "UPDATE COARSE GRID PARALLEL TOPOLOGY ");
 
 
-    // find exchange edges - first send exchangeedges locnum, v1, v2
+    // UpdateCoarseGridGlobal();
-    // receive distant distnum, v1, v2
-    // find matching
-    const MeshTopology & topology = mesh.GetTopology();
-
-    UpdateCoarseGridGlobal();
-
 
 
     
@@ -246,11 +241,12 @@ namespace netgen
     static int timerf = NgProfiler::CreateTimer ("UpdateCoarseGrid - ex faces");
 
 
-    Array<int,1> glob2loc;
     Array<int> cnt_send(ntasks-1);
 
     NgProfiler::StartTimer (timere);
 
+
+    const MeshTopology & topology = mesh.GetTopology();
     int nfa = topology . GetNFaces();
     int ned = topology . GetNEdges();
     
@@ -278,13 +274,11 @@ namespace netgen
 	gv2e.Set (es, edge);
 
 	for (int dest = 1; dest < ntasks; dest++)
-	  {
+	  if (IsExchangeVert (dest, v1) && IsExchangeVert (dest, v2))
-	    if (IsExchangeVert (dest, v1) && IsExchangeVert (dest, v2))
+	    {
-	      {
+	      send_edges.Add (dest-1, es[0]);
-		send_edges.Add (dest-1, es[0]);
+	      send_edges.Add (dest-1, es[1]);
-		send_edges.Add (dest-1, es[1]);
+	    }
-	      }
-	  }
       }
 
     TABLE<int> recv_edges(ntasks-1);
@@ -305,13 +299,69 @@ namespace netgen
  
     NgProfiler::StopTimer (timere);
 
-
     MPI_Barrier (MPI_LocalComm);
 
+
     if (mesh.GetDimension() == 3)
       {
 	NgProfiler::StartTimer (timerf);
+
+	// exchange faces
+	cnt_send = 0;
+	Array<int> verts;
+	for (int face = 1; face <= nfa; face++)
+	  {
+	    topology.GetFaceVertices (face, verts);
+	    for (int dest = 1; dest < ntasks; dest++)
+	      if (IsExchangeVert (dest, verts[0]) && 
+		  IsExchangeVert (dest, verts[1]) &&
+		  IsExchangeVert (dest, verts[2]))
+		cnt_send[dest-1]+=3;
+	  }
+    
+	TABLE<int> send_faces(cnt_send);
+	INDEX_3_HASHTABLE<int> gv2f(2*nfa);
+
+	for (int face = 1; face <= nfa; face++)
+	  {
+	    topology.GetFaceVertices (face, verts);
+	    INDEX_3 fs (GetGlobalPNum(verts[0]),
+			GetGlobalPNum(verts[1]),
+			GetGlobalPNum(verts[2]));
+	    fs.Sort();
+
+	    gv2f.Set (fs, face);
+
+	    for (int dest = 1; dest < ntasks; dest++)
+	      if (IsExchangeVert (dest, verts[0]) && 
+		  IsExchangeVert (dest, verts[1]) &&
+		  IsExchangeVert (dest, verts[2]))
+		{
+		  send_faces.Add (dest-1, fs[0]);
+		  send_faces.Add (dest-1, fs[1]);
+		  send_faces.Add (dest-1, fs[2]);
+		}
+	  }
 	
+	TABLE<int> recv_faces(ntasks-1);
+	MyMPI_ExchangeTable (send_faces, recv_faces, MPI_TAG_MESH+9, MPI_LocalComm);
+	
+	for (int sender = 1; sender < ntasks; sender ++)
+	  if (id != sender)
+	    {
+	      FlatArray<int> recvarray = recv_faces[sender-1];
+	      for (int ii = 0; ii < recvarray.Size(); ii+=3)
+		{ 
+		  INDEX_3 gv123 (recvarray[ii],recvarray[ii+1],recvarray[ii+2]);
+		  if (gv2f.Used (gv123))
+		    SetDistantFaceNum (sender, gv2f.Get(gv123));
+		}
+	    }
+
+
+	/*
+	  Array<int,1> glob2loc;
+
 	int maxface = 0;
 	for (int face = 1; face <= nfa; face++)
 	  maxface = max (maxface, GetGlobalFaceNum (face));
@@ -351,7 +401,7 @@ namespace netgen
 		    send_faces.Add (dest-1, face);
 		  }
 	      }
-      }
+	  }
 	TABLE<int> recv_faces(ntasks-1);
 	MyMPI_ExchangeTable (send_faces, recv_faces, MPI_TAG_MESH+8, MPI_LocalComm);
 	
@@ -373,7 +423,7 @@ namespace netgen
 		    }
 		}
 	    } 
-	
+	*/	
 	
 	NgProfiler::StopTimer (timerf);
       }

libsrc/meshing/paralleltop.hpp

@@ -76,6 +76,11 @@ namespace netgen
       for ( int i = 0; i < loc2distface[locfacenum-1].Size(); i++ )
 	distfacenums[i] = loc2distface[locfacenum-1][i];
     } 
+
+    void GetDistantFaceNums (int locfacenum, Array<int> & distfacenums ) const
+    {
+      distfacenums = loc2distface[locfacenum-1];
+    }
     
     void GetDistantEdgeNums (int locedgenum, int * distedgenums ) const
     {
@@ -83,15 +88,14 @@ namespace netgen
 	distedgenums[i] = loc2distedge[locedgenum-1][i];
     } 
 
+    void GetDistantEdgeNums (int locedgenum, Array<int> & distedgenums ) const
+    {
+      distedgenums = loc2distedge[locedgenum-1];
+    } 
+
     bool IsExchangeVert (int dest, int vnum) const
     {
       return loc2distvert[vnum-1].Contains (dest);
-      /*
-      FlatArray<int> exchange = loc2distvert[vnum-1];
-      for (int i = 0; i < exchange.Size(); i ++)
-	if (exchange[i] == dest) return true;
-      return false;
-      */
     }
   };
  

libsrc/meshing/topology.cpp

@@ -120,12 +120,12 @@ namespace netgen
 	  cnt[el[j]]++;
       }
 
-    vert2element = new TABLE<int,PointIndex::BASE> (cnt);
+    vert2element = new TABLE<ElementIndex,PointIndex::BASE> (cnt);
     for (ElementIndex ei = 0; ei < ne; ei++)
       {
 	const Element & el = mesh[ei];
 	for (int j = 0; j < el.GetNV(); j++)
-	  vert2element->AddSave (el[j], ei+1);
+	  vert2element->AddSave (el[j], ei);
       }
 
     cnt = 0;
@@ -160,15 +160,25 @@ namespace netgen
 	vert2segment->AddSave (seg[1], i);
       }
 
-
     if (buildedges)
       {
 	static int timer1 = NgProfiler::CreateTimer ("topology::buildedges");
-	NgProfiler::RegionTimer reg1 (timer1);
+	static int t1a = NgProfiler::CreateTimer ("topology - edges 1");
+	static int t1b = NgProfiler::CreateTimer ("topology - edges 2");
+	static int t1b1 = NgProfiler::CreateTimer ("topology - edges 21");
+	static int t1b11 = NgProfiler::CreateTimer ("topology - edges 211");
+	static int t1b2 = NgProfiler::CreateTimer ("topology - edges 22");
+	static int t1b3 = NgProfiler::CreateTimer ("topology - edges 23");
+	static int t1c = NgProfiler::CreateTimer ("topology - edges 3");
 
+
+	NgProfiler::RegionTimer reg1 (timer1);
+	
 	if (id == 0)
 	  PrintMessage (5, "Update edges ");
       
+	NgProfiler::StartTimer (t1a);
+
 	edges.SetSize(ne);
 	surfedges.SetSize(nse); 
 	segedges.SetSize(nseg);
@@ -192,44 +202,26 @@ namespace netgen
 	cnt = 0;
 	for (int i = mesh.mlbetweennodes.Begin(); i < mesh.mlbetweennodes.End(); i++)
 	  {
-	    /*  JS, Oct 2009
+	    INDEX_2 parents = Sort (mesh.mlbetweennodes[i]);
-		int pa[2];
-		pa[0] = mesh.mlbetweennodes[i].I1();
-		pa[1] = mesh.mlbetweennodes[i].I2();
-		if (pa[0] > pa[1]) Swap (pa[0], pa[1]);
-		if (pa[0] > 0)
-		cnt.Elem(pa[0])++;
-	    */
-	    INDEX_2 parents = mesh.mlbetweennodes[i];
-	    parents.Sort();
 	    if (parents[0] >= PointIndex::BASE) cnt[parents[0]]++;
 	  }
-
 	TABLE<int,PointIndex::BASE> vert2vertcoarse (cnt);
 	for (int i = mesh.mlbetweennodes.Begin(); i < mesh.mlbetweennodes.End(); i++)
 	  {
-	    /*
+	    INDEX_2 parents = Sort (mesh.mlbetweennodes[i]);
-	      int pa[2];
-	      pa[0] = mesh.mlbetweennodes[i].I1();
-	      pa[1] = mesh.mlbetweennodes[i].I2();
-	      if (pa[0] > pa[1]) swap (pa[0], pa[1]);
-	      if (pa[0] > 0)
-	      vert2vertcoarse.AddSave1 (pa[0], pa[1]);
-	    */
-	    INDEX_2 parents = mesh.mlbetweennodes[i];
-	    parents.Sort();
 	    if (parents[0] > PointIndex::BASE) vert2vertcoarse.AddSave (parents[0], parents[1]);
 	  }
 
-
 	Array<int,PointIndex::BASE> edgenr(nv);
 	Array<int,PointIndex::BASE> edgeflag(nv);
 	Array<int> vertex2;
 
-	edgeflag = 0;
+	edgeflag = PointIndex::BASE-1;
 
 	ned = edge2vert.Size();
-	Array<INDEX_3> missing;
+
+	NgProfiler::StopTimer (t1a);
+	NgProfiler::StartTimer (t1b);
 
 	for (int i = PointIndex::BASE; i < nv+PointIndex::BASE; i++)
 	  {
@@ -243,34 +235,30 @@ namespace netgen
 		edgenr[i2] = ednr;
 	      }
 
-	    for (int j = 0; j < vert2vertcoarse[i].Size(); j++)      // fix by Markus
+	    for (int j = 0; j < vert2vertcoarse[i].Size(); j++)    
 	      {
 		int v2 = vert2vertcoarse[i][j];
 		if (edgeflag[v2] < i)
 		  {
-		    *testout << "do we really need that ??????????????" << endl;
-		    // ned++;
-		    // edgenr[v2] = ned;
 		    edgeflag[v2] = i;
 		    vertex2.Append (v2);
-		    // missing.Append (INDEX_3(i,v2,ned));
 		  }
 	      }
 
-	    for (int j = 0; j < (*vert2element)[i].Size(); j++)
+	    NgProfiler::StartTimer (t1b1);
-	      {
-		int elnr = (*vert2element)[i][j];
-		const Element & el = mesh.VolumeElement (elnr);
 
+	    FlatArray<ElementIndex> v2els = (*vert2element)[i];
+	    for (int j = 0; j < v2els.Size(); j++)
+	      {
+		const Element & el = mesh[v2els[j]];
 		int neledges = GetNEdges (el.GetType());
 		const ELEMENT_EDGE * eledges = GetEdges0 (el.GetType());
-	  
 		for (int k = 0; k < neledges; k++)
 		  {
 		    INDEX_2 edge(el[eledges[k][0]], el[eledges[k][1]]);
 		    edge.Sort();
 		    if (edge.I1() != i) continue;
-	     
+		    
 		    if (edgeflag[edge.I2()] < i)
 		      {
 			vertex2.Append (edge.I2());
@@ -279,6 +267,8 @@ namespace netgen
 		  }
 	      }
 
+	    NgProfiler::StopTimer (t1b1);	    
+
 	    for (int j = 0; j < (*vert2surfelement)[i].Size(); j++)
 	      {
 		int elnr = (*vert2surfelement)[i][j];
@@ -317,27 +307,23 @@ namespace netgen
 		  }   
 	      }
 
-	    
+	    NgProfiler::StartTimer (t1b2);	    
+
 	    QuickSort (vertex2);
 	    for (int j = 0; j < vertex2.Size(); j++)
-	      edgenr[vertex2[j]] = ++ned;
-
-	    for (int j = 0; j < vert2vertcoarse[i].Size(); j++)      // fix by Markus
 	      {
-		int v2 = vert2vertcoarse[i][j];
+		edgenr[vertex2[j]] = ++ned;
-		if (edgeflag[v2] < i)
+		edge2vert.Append (INDEX_2 (i, vertex2[j]));
-		  {
-		    // ned++;
-		    // edgenr[v2] = ned;
-		    // edgeflag[v2] = i;
-		    missing.Append (INDEX_3(i,v2,edgenr[v2]));
-		  }
 	      }
 
+	    NgProfiler::StopTimer (t1b2);	    
+	    NgProfiler::StartTimer (t1b3);	    
+
+
 	    for (int j = 0; j < (*vert2element)[i].Size(); j++)
 	      {
-		int elnr = (*vert2element)[i][j];
+		ElementIndex elnr = (*vert2element)[i][j];
-		const Element & el = mesh.VolumeElement (elnr);
+		const Element & el = mesh[elnr];
 
 		int neledges = GetNEdges (el.GetType());
 		const ELEMENT_EDGE * eledges = GetEdges0 (el.GetType());
@@ -349,12 +335,11 @@ namespace netgen
 		    int edgedir = (edge.I1() > edge.I2());
 		    if (edgedir) swap (edge.I1(), edge.I2());
 	     
-		    if (edge.I1() != i)
+		    if (edge.I1() != i) continue;
-		      continue;
 
 		    int edgenum = edgenr[edge.I2()];
 		    if (edgedir) edgenum *= -1;
-		    edges.Elem(elnr)[k] = edgenum;
+		    edges[elnr][k] = edgenum;
 		  }
 	      }
 
@@ -397,100 +382,16 @@ namespace netgen
 		if (edgedir) edgenum *= -1;
 		segedges.Elem(elnr) = edgenum;
 	      }
+	    NgProfiler::StopTimer (t1b3);	    
 	  }
 
-	edge2vert.SetSize (ned);
+	NgProfiler::StopTimer (t1b);
-	for (int i = 1; i <= ne; i++)
-	  {
-	    const Element & el = mesh.VolumeElement (i);
-      
-	    int neledges = GetNEdges (el.GetType());
-	    const ELEMENT_EDGE * eledges = GetEdges0 (el.GetType());
-	  
-	    for (int k = 0; k < neledges; k++)
-	      {
-		INDEX_2 edge(el[eledges[k][0]], el[eledges[k][1]]);
-	  
-		int edgedir = (edge.I1() > edge.I2());
-		if (edgedir) swap (edge.I1(), edge.I2());
-
-		int edgenum = abs (edges.Elem(i)[k]);
-
-		edge2vert.Elem(edgenum)[0] = edge.I1();
-		edge2vert.Elem(edgenum)[1] = edge.I2();
-	      }
-	  }
-
-	/*
-	*testout << "edge 2 vert:" << endl;
-	for (int i = 0; i < edge2vert.Size(); i++)
-	  *testout << edge2vert[i][0] << " " << edge2vert[i][1] << endl;
-	  */
-
-	for (int i = 1; i <= nse; i++)
-	  {
-	    const Element2d & el = mesh.SurfaceElement (i);
-      
-	    int neledges = GetNEdges (el.GetType());
-	    const ELEMENT_EDGE * eledges = GetEdges0 (el.GetType());
-	  
-	    for (int k = 0; k < neledges; k++)
-	      {
-		INDEX_2 edge(el[eledges[k][0]], el[eledges[k][1]]);
-	  
-		int edgedir = (edge.I1() > edge.I2());
-		if (edgedir) swap (edge.I1(), edge.I2());
-
-		int edgenum = abs (surfedges.Elem(i)[k]);
-
-		edge2vert.Elem(edgenum)[0] = edge.I1();
-		edge2vert.Elem(edgenum)[1] = edge.I2();
-	      }
-	  }
-
-	for (int i = 1; i <= nseg; i++)
-	  {
-	    const Segment & el = mesh.LineSegment (i);
-      
-	    INDEX_2 edge(el[0], el[1]);
-	    int edgedir = (edge.I1() > edge.I2());
-	    if (edgedir) swap (edge.I1(), edge.I2());
-	  
-	    int edgenum = abs (segedges.Elem(i));
-	  
-	    edge2vert.Elem(edgenum)[0] = edge.I1();
-	    edge2vert.Elem(edgenum)[1] = edge.I2();
-	  }
-
-	for (int i = 1; i <= missing.Size(); i++)
-	  {
-	    INDEX_3 i3 = missing.Get(i);
-	    edge2vert.Elem(i3.I3())[0] = i3.I1();
-	    edge2vert.Elem(i3.I3())[1] = i3.I2();
-	  }
-	
-      
-	/*
-	  (*testout) << "edge table:" << endl;
-	  (*testout) << "edge2vert:" << endl;
-	  for (int i = 1; i <= edge2vert.Size(); i++)
-	  (*testout) << "edge " << i << ", v1,2 = " << edge2vert.Elem(i)[0] << ", " << edge2vert.Elem(i)[1] << endl;
-	  (*testout) << "surfedges:" << endl;
-	  for (int i = 1; i <= surfedges.Size(); i++)
-	  (*testout) << "el " << i << ", edges = " 
-	  << surfedges.Elem(i)[0] << ", "
-	  << surfedges.Elem(i)[1] << ", "
-	  << surfedges.Elem(i)[2] << endl;
-	*/ 
       }
 
 
-    //  cout << "build edges done" << endl;
-
-
 
     // generate faces
-    if (buildfaces) //  && mesh.GetDimension() == 3)
+    if (buildfaces) 
       {
 	static int timer2 = NgProfiler::CreateTimer ("topology::buildfaces");
 	NgProfiler::RegionTimer reg2 (timer2);
@@ -804,8 +705,8 @@ namespace netgen
 
 		for (int j = 0; j < (*vert2element)[v].Size(); j++)
 		  {
-		    int elnr = (*vert2element)[v][j];
+		    ElementIndex elnr = (*vert2element)[v][j];
-		    const Element & el = mesh.VolumeElement (elnr);
+		    const Element & el = mesh[elnr];
 	  
 		    int nelfaces = GetNFaces (el.GetType());
 		    const ELEMENT_FACE * elfaces = GetFaces1 (el.GetType());
@@ -840,7 +741,7 @@ namespace netgen
 			      INDEX_4 hface(face.I1(),face.I2(),face.I3(),0);
 			      face2vert.Append (hface);
 			    }
-			  faces.Elem(elnr)[j] = 8*(facenum-1)+facedir+1;
+			  faces[elnr][j] = 8*(facenum-1)+facedir+1;
 			}
 		
 		      else
@@ -893,7 +794,7 @@ namespace netgen
 			      face2vert.Append (hface);
 			    }
 		      
-			  faces.Elem(elnr)[j] = 8*(facenum-1)+facedir+1;
+			  faces[elnr][j] = 8*(facenum-1)+facedir+1;
 			}
 		  }
 
@@ -1121,20 +1022,20 @@ namespace netgen
 			    (*testout) << mesh[(PointIndex)face2vert[i].I(j+1)] << " ";
 			(*testout) << endl;
 
-			FlatArray<int> vertels = GetVertexElements (face2vert[i].I(1));
+			FlatArray<ElementIndex> vertels = GetVertexElements (face2vert[i].I(1));
 			for (int k = 0; k < vertels.Size(); k++)
 			  {
 			    int elfaces[10], orient[10];
-			    int nf = GetElementFaces (vertels[k], elfaces, orient);
+			    int nf = GetElementFaces (vertels[k]+1, elfaces, orient);
 			    for (int l = 0; l < nf; l++)
 			      if (elfaces[l] == i)
 				{
-				  (*testout) << "is face of element " << vertels[k] << endl;
+				  // (*testout) << "is face of element " << vertels[k] << endl;
 			    
 				  if (mesh.coarsemesh && mesh.hpelements->Size() == mesh.GetNE() )
 				    {
 				      const HPRefElement & hpref_el =
-					(*mesh.hpelements) [ mesh.VolumeElement (vertels[k]).hp_elnr];
+					(*mesh.hpelements) [ mesh[vertels[k]].hp_elnr];
 				      (*testout) << "coarse eleme = " << hpref_el.coarse_elnr << endl;
 				    }
 
@@ -1408,21 +1309,18 @@ namespace netgen
 
   void MeshTopology :: GetSurfaceElementEdges (int elnr, Array<int> & eledges) const
   {
-    int i;
+    int ned = GetNEdges (mesh.SurfaceElement(elnr).GetType());
-    if (mesh.GetDimension()==3 || 1)
+    eledges.SetSize (ned);
-      {
+    for (int i = 0; i < ned; i++)
-	int ned = GetNEdges (mesh.SurfaceElement(elnr).GetType());
+      eledges[i] = abs (surfedges.Get(elnr)[i]);
-	eledges.SetSize (ned);
+  }
-	for (i = 1; i <= ned; i++)
+
-	  eledges.Elem(i) = abs (surfedges.Get(elnr)[i-1]);
+  void MeshTopology :: GetEdges (SurfaceElementIndex elnr, Array<int> & eledges) const
-      }
+  {
-    else
+    int ned = GetNEdges (mesh[elnr].GetType());
-      {
+    eledges.SetSize (ned);
-	cout << "surfeledge(" << elnr << ") = " << flush;
+    for (int i = 0; i < ned; i++)
-	eledges.SetSize(1); 
+      eledges[i] = abs (surfedges[elnr][i])-1;
-	eledges.Elem(1) = abs (segedges.Get(elnr));
-	cout << eledges.Elem(1) << endl;
-      }
   }
 
   int MeshTopology :: GetSurfaceElementFace (int elnr) const
@@ -1430,13 +1328,19 @@ namespace netgen
     return (surffaces.Get(elnr)-1) / 8 + 1;  
   }
 
+  int MeshTopology :: GetFace (SurfaceElementIndex elnr) const
+  {
+    return (surffaces[elnr]-1) / 8;  
+  }
+
+
   void MeshTopology :: 
   GetSurfaceElementEdgeOrientations (int elnr, Array<int> & eorient) const
   {
     int ned = GetNEdges (mesh.SurfaceElement(elnr).GetType());
     eorient.SetSize (ned);
-    for (int i = 1; i <= ned; i++)
+    for (int i = 0; i < ned; i++)
-      eorient.Elem(i) = (surfedges.Get(elnr)[i-1] > 0) ? 1 : -1;
+      eorient[i] = (surfedges.Get(elnr)[i] > 0) ? 1 : -1;
   }
 
   int MeshTopology :: GetSurfaceElementFaceOrientation (int elnr) const
@@ -1481,10 +1385,9 @@ namespace netgen
   void MeshTopology :: GetFaceVertices (int fnr, Array<int> & vertices) const
   {
     vertices.SetSize(4);
-    int i;
+    for (int i = 0; i < 4; i++)
-    for (i = 1; i <= 4; i++)
+      vertices[i] = face2vert.Get(fnr)[i];
-      vertices.Elem(i) = face2vert.Get(fnr)[i-1];
+    if (vertices[3] == 0)
-    if (vertices.Elem(4) == 0)
       vertices.SetSize(3);
   }
 
@@ -1501,6 +1404,12 @@ namespace netgen
     v2 = edge2vert.Get(ednr)[1];
   }
 
+  void MeshTopology :: GetEdgeVertices (int ednr, PointIndex & v1, PointIndex & v2) const
+  {
+    v1 = edge2vert.Get(ednr)[0];
+    v2 = edge2vert.Get(ednr)[1];
+  }
+
 
   void MeshTopology :: GetFaceEdges (int fnr, Array<int> & fedges, bool withorientation) const
   {
@@ -1525,12 +1434,12 @@ namespace netgen
   
 
     //  GetVertexElements (pi[0], els);
-    FlatArray<int> els= GetVertexElements (pi[0]);
+    FlatArray<ElementIndex> els = GetVertexElements (pi[0]);
 
     // find one element having all vertices of the face
     for (int i = 0; i < els.Size(); i++)
       {
-	const Element & el = mesh.VolumeElement(els[i]);
+	const Element & el = mesh[els[i]];
 	int nref_faces = GetNFaces (el.GetType());
 	const ELEMENT_FACE * ref_faces = GetFaces1 (el.GetType());
 	int nfa_ref_edges = GetNEdges (GetFaceType(fnr));
@@ -1556,7 +1465,7 @@ namespace netgen
 	  {
 	    const ELEMENT_EDGE * fa_ref_edges = GetEdges1 (GetFaceType(fnr)); 
 	    fedges.SetSize(nfa_ref_edges);
-	    GetElementEdges (els[i], eledges);
+	    GetElementEdges (els[i]+1, eledges);
 	  
 	    for (int j = 0; j < eledges.Size(); j++)
 	      {
@@ -1600,6 +1509,13 @@ namespace netgen
 	    return;
 	  }
       }   
+
+    int surfel = GetFace2SurfaceElement(fnr);
+    if (surfel != 0)
+      {
+	GetSurfaceElementEdges (surfel, fedges);
+	return;
+      }
   }
 
 
@@ -1609,24 +1525,23 @@ namespace netgen
   }
 
 
-  void MeshTopology :: GetVertexElements (int vnr, Array<int> & elements) const
+  void MeshTopology :: GetVertexElements (int vnr, Array<ElementIndex> & elements) const
   {
     if (vert2element)
       {
-	int i; 
 	int ne = vert2element->EntrySize(vnr);
 	elements.SetSize(ne);
-	for (i = 1; i <= ne; i++)
+	for (int i = 1; i <= ne; i++)
 	  elements.Elem(i) = vert2element->Get(vnr, i);
       }
   }
 
 
-  FlatArray<int> MeshTopology :: GetVertexElements (int vnr) const
+  FlatArray<ElementIndex> MeshTopology :: GetVertexElements (int vnr) const
   {
     if (vert2element)
       return (*vert2element)[vnr];
-    return FlatArray<int> (0,0);
+    return FlatArray<ElementIndex> (0,0);
   }
 
   FlatArray<int> MeshTopology :: GetVertexSurfaceElements (int vnr) const
@@ -1653,13 +1568,14 @@ namespace netgen
 
   int MeshTopology :: GetVerticesEdge ( int v1, int v2 ) const
   {
-    Array<int> elements_v1, elementedges;
+    Array<ElementIndex> elements_v1;
+    Array<int> elementedges;
     GetVertexElements ( v1, elements_v1);
     int edv1, edv2;
 
     for ( int i = 0; i < elements_v1.Size(); i++ )
       {
-	GetElementEdges( elements_v1[i], elementedges );
+	GetElementEdges( elements_v1[i]+1, elementedges );
 	for ( int ed = 0; ed < elementedges.Size(); ed ++)
 	  {
 	    GetEdgeVertices( elementedges[ed], edv1, edv2 );
@@ -1678,14 +1594,14 @@ namespace netgen
   {
     int v1, v2;
     GetEdgeVertices ( GetSegmentEdge (segnr), v1, v2 );
-    Array<int> volels1, volels2;
+    Array<ElementIndex> volels1, volels2;
     GetVertexElements ( v1, volels1 );
     GetVertexElements ( v2, volels2 );
     volels.SetSize(0);
 
     for ( int eli1=1; eli1 <= volels1.Size(); eli1++)
       if ( volels2.Contains( volels1.Elem(eli1) ) )
-	volels.Append ( volels1.Elem(eli1) );
+	volels.Append ( volels1.Elem(eli1)+1 );
   }
 
   void MeshTopology :: 

libsrc/meshing/topology.hpp

@@ -28,7 +28,7 @@ class MeshTopology
   MoveableArray<int> surffaces;
   MoveableArray<INDEX_2> surf2volelement;
   MoveableArray<int> face2surfel;
-  TABLE<int,PointIndex::BASE> *vert2element;
+  TABLE<ElementIndex,PointIndex::BASE> *vert2element;
   TABLE<int,PointIndex::BASE> *vert2surfelement;
   TABLE<int,PointIndex::BASE> *vert2segment;
   int timestamp;
@@ -68,6 +68,7 @@ public:
   
   int GetSegmentEdge (int segnr) const { return abs(segedges[segnr-1]); }
   int GetSegmentEdgeOrientation (int segnr) const { return sgn(segedges[segnr-1]); }
+  int GetEdge (SegmentIndex segnr) const { return abs(segedges[segnr])-1; }
 
   void GetSegmentEdge (int segnr, int & enr, int & orient) const
   {
@@ -86,6 +87,7 @@ public:
   void GetFaceVertices (int fnr, Array<int> & vertices) const;
   void GetFaceVertices (int fnr, int * vertices) const;
   void GetEdgeVertices (int enr, int & v1, int & v2) const;
+  void GetEdgeVertices (int enr, PointIndex & v1, PointIndex & v2) const;
   const int * GetEdgeVerticesPtr (int enr) const { return &edge2vert[enr][0]; }
   const int * GetFaceVerticesPtr (int fnr) const { return &face2vert[fnr][0]; }
   void GetFaceEdges (int fnr, Array<int> & edges, bool withorientation = false) const;
@@ -96,6 +98,8 @@ public:
   int GetSurfaceElementFace (int elnr) const;
   void GetSurfaceElementEdgeOrientations (int elnr, Array<int> & eorient) const;
   int GetSurfaceElementFaceOrientation (int elnr) const;
+  void GetEdges (SurfaceElementIndex elnr, Array<int> & edges) const;
+  int GetFace (SurfaceElementIndex elnr) const;
 
   int GetSurfaceElementEdges (int elnr, int * edges, int * orient) const;
 
@@ -115,8 +119,8 @@ public:
 
   int GetFace2SurfaceElement (int fnr) const { return face2surfel[fnr-1]; }
   
-  void GetVertexElements (int vnr, Array<int> & elements) const;
+  void GetVertexElements (int vnr, Array<ElementIndex> & elements) const;
-  FlatArray<int> GetVertexElements (int vnr) const;
+  FlatArray<ElementIndex> GetVertexElements (int vnr) const;
 
   void GetVertexSurfaceElements( int vnr, Array<int>& elements ) const;
   FlatArray<int> GetVertexSurfaceElements (int vnr) const;

libsrc/occ/Makefile.am

@@ -5,7 +5,7 @@ Partition_Inter2d.ixx Partition_Loop2d.ixx Partition_Loop.ixx \
 Partition_Inter3d.ixx Partition_Loop3d.ixx Partition_Spliter.ixx \
 Partition_Inter2d.jxx Partition_Loop2d.jxx Partition_Loop.jxx \
 Partition_Inter3d.jxx Partition_Loop3d.jxx Partition_Spliter.jxx \
-utilities.h
+utilities.h vsocc.hpp
 
 
 AM_CPPFLAGS = -I$(top_srcdir)/libsrc/include  $(OCCFLAGS) $(TCL_INCLUDES)

libsrc/occ/Partition_Loop.cxx

@@ -178,7 +178,7 @@ static Standard_Boolean  SelectEdge(const TopoDS_Face&    F,
       }
     }
 
-    Standard_Real anglemax = - PI;
+    Standard_Real anglemax = - M_PI;
     TopoDS_Edge   SelectedEdge;	
     for ( itl.Initialize(LE); itl.More(); itl.Next()) {
       const TopoDS_Edge& E = TopoDS::Edge(itl.Value());

libsrc/visualization/vssolution.cpp

@@ -645,17 +645,8 @@ namespace netgen
         subdivision_timestamp > surfellinetimestamp ||
         (deform && solutiontimestamp > surfellinetimestamp) || 
         zoomall)
-      {
+      {      
-        if (linelist)
+        DrawSurfaceElementLines();      
-          glDeleteLists (linelist, 1);
-      
-        linelist = glGenLists (1);
-        glNewList (linelist, GL_COMPILE);
-      
-        DrawSurfaceElementLines();
-
-        glEndList ();
-      
         surfellinetimestamp = max2 (solutiontimestamp, mesh->GetTimeStamp());
       }
 
@@ -1394,28 +1385,51 @@ namespace netgen
 
   void  VisualSceneSolution :: DrawSurfaceElementLines ()
   {
+#ifdef PARALLELGL
+    if (id == 0 && ntasks > 1)
+      {
+	InitParallelGL();
+
+	par_surfellists.SetSize (ntasks);
+
+	MyMPI_SendCmd ("redraw");
+	MyMPI_SendCmd ("solsurfellinelist");
+
+	for ( int dest = 1; dest < ntasks; dest++ )
+	  MyMPI_Recv (par_surfellists[dest], dest, MPI_TAG_VIS);
+
+	if (linelist)
+	  glDeleteLists (linelist, 1);
+
+	linelist = glGenLists (1);
+	glNewList (linelist, GL_COMPILE);
+	
+	for ( int dest = 1; dest < ntasks; dest++ )
+	  glCallList (par_surfellists[dest]);
+	
+	glEndList();
+	return;
+      }
+#endif
+
+    if (linelist)
+      glDeleteLists (linelist, 1);
+    
+    linelist = glGenLists (1);
+    glNewList (linelist, GL_COMPILE);
+
     glLineWidth (1.0f);
-    // glNormal3d (1, 0, 0);
 
     int nse = mesh->GetNSE();
-
     CurvedElements & curv = mesh->GetCurvedElements();
 
     int n = 1 << subdivisions;
     ArrayMem<Point<2>, 65> ptsloc(n+1);
     ArrayMem<Point<3>, 65> ptsglob(n+1);
 
-    /*
-    double trigpts[3][2]  = { { 0, 0 }, { 1, 0 }, { 0, 1} };
-    double trigvecs[3][2] = { { 1, 0 }, { -1,1 }, { 0, -1} };
-    */
     double trigpts[3][2]  = { { 0, 0 }, { 0, 1 }, { 1, 0} };
     double trigvecs[3][2] = { { 1, 0 }, { 0, -1 }, { -1, 1} };
-    
+
-    /*
-    double quadpts[4][2] = { { 0, 0 }, { 1, 0 }, { 1, 1 }, { 0, 1} };
-    double quadvecs[4][2] = { { 1, 0 }, { 0, 1 }, { -1, 0}, { 0, -1} };
-    */
     double quadpts[4][2]  = { { 0, 0 },  { 1, 1 }, { 0, 1}, { 1, 0 } };
     double quadvecs[4][2] = { { 1, 0 },  { -1, 0}, { 0, -1}, { 0, 1 } };
 
@@ -1423,20 +1437,7 @@ namespace netgen
       {
         Element2d & el = (*mesh)[sei];
 
-        // bool curved = curv.IsSurfaceElementCurved (sei);
-
         int nv = (el.GetType() == TRIG || el.GetType() == TRIG6) ? 3 : 4;
-        /*
-        Point<3> p1, p2, p3, p4;
-        if (!curved)
-	{
-            p1 = (*mesh)[el[0]];
-            p2 = (*mesh)[el[1]];
-            p3 = (*mesh)[el[2]];
-            if (nv == 4) p4 = (*mesh)[el[3]];
-          }
-	*/
-
         for (int k = 0; k < nv; k++)
           {
             Point<2> p0;
@@ -1454,44 +1455,29 @@ namespace netgen
 
             glBegin (GL_LINE_STRIP);
 
-            // if (curved)
+	    for (int ix = 0; ix <= n; ix++)
-              {
+	      ptsloc[ix] = p0 + (double(ix) / n) * vtau;
-                for (int ix = 0; ix <= n; ix++)
+	    
-                  ptsloc[ix] = p0 + (double(ix) / n) * vtau;
+	    curv.CalcMultiPointSurfaceTransformation (&ptsloc, sei, &ptsglob, 0);
-                    
+	    
-                curv.CalcMultiPointSurfaceTransformation (&ptsloc, sei, &ptsglob, 0);
+	    for (int ix = 0; ix <= n; ix++)
-
+	      {
-                for (int ix = 0; ix <= n; ix++)
+		if (deform)
-                  {
+		  ptsglob[ix] += GetSurfDeformation (sei, k, ptsloc[ix](0), ptsloc[ix](1));
-                    if (deform)
+		glVertex3dv (ptsglob[ix]);
-                      ptsglob[ix] += GetSurfDeformation (sei, k, ptsloc[ix](0), ptsloc[ix](1));
+	      }
-                    glVertex3dv (ptsglob[ix]);
+	    
-                  }
-              }
-	      /*
-            else
-              {
-                for (int ix = 0; ix <= n; ix++)
-                  {
-                    Point<2> p = p0 + (double(ix) / n) * vtau;
-
-                    Point<3> pnt;
-
-                    if (nv == 3)
-                      pnt = p3 + p(0) * (p1-p3) + p(1) * (p2-p3);
-                    else
-                      pnt = p1 + p(0) * (p2-p1) + p(1) * (p4-p1) + p(0)*p(1) * ( (p1-p2)+(p3-p4) );
-                
-                    if (deform)
-                      pnt += GetSurfDeformation (sei, p(0), p(1) );
-                    
-                    glVertex3dv (pnt);
-                  }
-              }
-	      */
             glEnd ();
           }
       }
+    glEndList ();
+
+
+#ifdef PARALLELGL
+    glFinish();
+    if (id > 0)
+      MyMPI_Send (linelist, 0, MPI_TAG_VIS);
+#endif
   }