From 36e517fbb1188757d4722046431225ca2b388989 Mon Sep 17 00:00:00 2001
From: Joachim Schoeberl <joachim.schoeberl@tuwien.ac.at>
Date: Fri, 26 Feb 2010 17:06:00 +0000
Subject: [PATCH] fixes for curved elements

---
 libsrc/meshing/curvedelems.cpp | 10 ++---
 libsrc/meshing/meshfunc.cpp    | 76 ++++++++++++++++------------------
 2 files changed, 40 insertions(+), 46 deletions(-)

diff --git a/libsrc/meshing/curvedelems.cpp b/libsrc/meshing/curvedelems.cpp
index b8a30470..fe6c0e27 100644
--- a/libsrc/meshing/curvedelems.cpp
+++ b/libsrc/meshing/curvedelems.cpp
@@ -2080,7 +2080,7 @@ namespace netgen
 	  double lami[6] = { xi(0), xi(1), 1-xi(0)-xi(1), xi(0), xi(1), 1-xi(0)-xi(1) };
 	  double lamiz[6] = { 1-xi(2), 1-xi(2), 1-xi(2), xi(2), xi(2), xi(2) };
 	  for (int i = 0; i < 6; i++)
-	    shapes(i) = lami[i%3] * ( (i < 3) ? (1-xi(2)) : xi(2) );
+	    shapes(i) = lami[i] * lamiz[i]; 
 	  for (int i = 6; i < info.ndof; i++)
 	    shapes(i) = 0;
 
@@ -2094,10 +2094,8 @@ namespace netgen
 	      int eorder = edgeorder[info.edgenrs[i]];
 	      if (eorder >= 2)
 		{
-		  int vi1 = (edges[i][0]-1), vi2 = (edges[i][1]-1);
+		  int vi1 = edges[i][0]-1, vi2 = edges[i][1]-1;
 		  if (el[vi1] > el[vi2]) swap (vi1, vi2);
-		  vi1 = vi1 % 3;
-		  vi2 = vi2 % 3;
 
 		  CalcScaledEdgeShape (eorder, lami[vi1]-lami[vi2], lami[vi1]+lami[vi2], &shapes(ii));
 		  double facz = (i < 3) ? (1-xi(2)) : xi(2);
@@ -2113,12 +2111,12 @@ namespace netgen
 	      int eorder = edgeorder[info.edgenrs[i]];
 	      if (eorder >= 2)
 		{
-		  int vi1 = (edges[i][0]-1), vi2 = (edges[i][1]-1);
+		  int vi1 = edges[i][0]-1, vi2 = edges[i][1]-1;
 		  if (el[vi1] > el[vi2]) swap (vi1, vi2);
 
 		  double bubz = lamiz[vi1]*lamiz[vi2];
 		  double polyz = lamiz[vi1] - lamiz[vi2];
-		  double bubxy = lami[(vi1)%3];
+		  double bubxy = lami[vi1];
 
 		  for (int j = 0; j < eorder-1; j++)
 		    {
diff --git a/libsrc/meshing/meshfunc.cpp b/libsrc/meshing/meshfunc.cpp
index f2d6fd8e..3c6a7b65 100644
--- a/libsrc/meshing/meshfunc.cpp
+++ b/libsrc/meshing/meshfunc.cpp
@@ -65,50 +65,46 @@ namespace netgen
      double globmaxh = mp.maxh;
 
      for (int k = 1; k <= mesh3d.GetNDomains(); k++)
-     {
-        if (multithread.terminate)
+       {
+	 if (multithread.terminate)
            break;
-
-        PrintMessage (2, "");
-        PrintMessage (1, "Meshing subdomain ", k, " of ", mesh3d.GetNDomains());
-        (*testout) << "Meshing subdomain " << k << endl;
-
-        mp.maxh = min2 (globmaxh, mesh3d.MaxHDomain(k));
-
-        mesh3d.CalcSurfacesOfNode();
-        mesh3d.FindOpenElements(k);
-
-        if (!mesh3d.GetNOpenElements())
+	 
+	 PrintMessage (2, "");
+	 PrintMessage (1, "Meshing subdomain ", k, " of ", mesh3d.GetNDomains());
+	 (*testout) << "Meshing subdomain " << k << endl;
+	 
+	 mp.maxh = min2 (globmaxh, mesh3d.MaxHDomain(k));
+	 
+	 mesh3d.CalcSurfacesOfNode();
+	 mesh3d.FindOpenElements(k);
+	 
+	 if (!mesh3d.GetNOpenElements())
            continue;
+	 
+	 
 
-
-
-        Box<3> domain_bbox( Box<3>::EMPTY_BOX ); 
-        /*
-        Point<3> (1e10, 1e10, 1e10),
-        Point<3> (-1e10, -1e10, -1e10));
-        */
-
-        for (SurfaceElementIndex sei = 0; sei < mesh3d.GetNSE(); sei++)
-        {
-           const Element2d & el = mesh3d[sei];
-           if (el.IsDeleted() ) continue;
-
-           if (mesh3d.GetFaceDescriptor(el.GetIndex()).DomainIn() == k ||
-              mesh3d.GetFaceDescriptor(el.GetIndex()).DomainOut() == k)
-
-              for (int j = 0; j < el.GetNP(); j++)
-                 domain_bbox.Add (mesh3d[el[j]]);
-        }
-        domain_bbox.Increase (0.01 * domain_bbox.Diam());
-
-
+	 Box<3> domain_bbox( Box<3>::EMPTY_BOX ); 
+	 
+	 for (SurfaceElementIndex sei = 0; sei < mesh3d.GetNSE(); sei++)
+	   {
+	     const Element2d & el = mesh3d[sei];
+	     if (el.IsDeleted() ) continue;
+	     
+	     if (mesh3d.GetFaceDescriptor(el.GetIndex()).DomainIn() == k ||
+		 mesh3d.GetFaceDescriptor(el.GetIndex()).DomainOut() == k)
+	       
+	       for (int j = 0; j < el.GetNP(); j++)
+		 domain_bbox.Add (mesh3d[el[j]]);
+	   }
+	 domain_bbox.Increase (0.01 * domain_bbox.Diam());
+	 
+	
         for (int qstep = 1; qstep <= 3; qstep++)
-        {
-           if (mesh3d.HasOpenQuads())
-           {
-              string rulefile = ngdir;
-
+	  {
+	    if (mesh3d.HasOpenQuads())
+	      {
+		string rulefile = ngdir;
+		
               const char ** rulep = NULL;
               switch (qstep)
               {