second order mesh merge

libsrc/meshing/meshclass.cpp

@@ -1470,15 +1470,28 @@ namespace netgen
 
         // int si = sel.GetIndex();
 
-        for (int j = 0; j < sel.GetNP(); j++)
+        if (sel.GetNP() <= 4)
-          {
+          for (int j = 0; j < sel.GetNP(); j++)
-            INDEX_2 i2;
+            {
-            i2.I1() = sel.PNumMod(j+1);
+              INDEX_2 i2;
-            i2.I2() = sel.PNumMod(j+2);
+              i2.I1() = sel.PNumMod(j+1);
-            i2.Sort();
+              i2.I2() = sel.PNumMod(j+2);
-            if (sel.GetNP() <= 4)
+              i2.Sort();
               boundaryedges->Set (i2, 1);
+            }
+        else if (sel.GetType()==TRIG6)
+          {
+            for (int j = 0; j < 3; j++)
+              {
+                INDEX_2 i2;
+                i2.I1() = sel[j];
+                i2.I2() = sel[(j+1)%3];
+                i2.Sort();
+                boundaryedges->Set (i2, 1);
+              }
           }
+        else 
+          cerr << "illegal elemenet for buildboundaryedges" << endl;
       }
 
 
@@ -1825,6 +1838,7 @@ namespace netgen
           for (ii = 0; ii < row.Size(); ii++)
             {
               hel = SurfaceElement(row[ii]);
+              if (hel.GetType() == TRIG6) hel.SetType(TRIG);
               int ind = hel.GetIndex();	  
 
               if (GetFaceDescriptor(ind).DomainIn() && 
@@ -4347,11 +4361,9 @@ namespace netgen
     int i = 0;
 
     const int maxits = 30;
-
     while(delta > 1e-16 && i<maxits)
       {
         curvedelems->CalcElementTransformation(lam,element-1,x,Jac);
-
         rhs = p-x;
         Jac.Solve(rhs,deltalam);
 
@@ -4367,7 +4379,6 @@ namespace netgen
     if(i==maxits)
       return false;
 
-
     for(i=0; i<3; i++)
       lami[i] = lam(i);
 
@@ -4568,7 +4579,6 @@ namespace netgen
               ii = locels.Get(i);
             else
               ii = i;
-
             if(ii == ps_startelement) continue;
 
             if(indices != NULL && indices->Size() > 0)

libsrc/meshing/meshtype.hpp

@@ -354,11 +354,22 @@ namespace netgen
     ///
     int GetNV() const
     {
+      if (typ == TRIG || typ == TRIG6)
+        return 3;
+      else
+        {
+#ifdef DEBUG
+          if (typ != QUAD && typ != QUAD6 && typ != QUAD8)
+            PrintSysError ("element2d::GetNV not implemented for typ", typ)
+#endif
+          return 4;
+        }
+      /*
       switch (typ)
 	{
 	case TRIG:
 	case TRIG6: return 3;
-
+          
 	case QUAD:
 	case QUAD8:
 	case QUAD6: return 4;
@@ -369,6 +380,7 @@ namespace netgen
 	    ;
 	}
       return np;
+      */
     }
 
     ///

libsrc/meshing/secondorder.cpp

@@ -1,4 +1,4 @@
-#include <mystdlib.h>
+ #include <mystdlib.h>
 #include "meshing.hpp"
 
 
@@ -13,13 +13,64 @@ namespace netgen
   
   void Refinement :: MakeSecondOrder (Mesh & mesh)
   {
-    int nseg, nse, ne;
+    /*
+      Berlin, 2014: if we have curved surface elements, keep them !
+    */
 
     mesh.ComputeNVertices();
-    mesh.SetNP(mesh.GetNV());
+    // mesh.SetNP(mesh.GetNV());
-  
+    mesh.SetNP(mesh.GetNP());  // setup multilevel-table
+
     INDEX_2_HASHTABLE<PointIndex> between(mesh.GetNP() + 5);
 
+    for (SurfaceElementIndex sei = 0; sei < mesh.GetNSE(); sei++)
+      {
+	const Element2d & el = mesh[sei];
+
+	static int betw_trig[3][3] =
+	  { { 1, 2, 3 }, { 0, 2, 4 }, { 0, 1, 5 } };
+	static int betw_quad6[2][3] =
+	  { { 0, 1, 4 }, { 3, 2, 5 } };
+	static int betw_quad8[4][3] =
+	  { { 0, 1, 4 }, { 3, 2, 5 },
+	    { 0, 3, 6 }, { 1, 2, 7 } };
+      
+	int onp = 0;
+	int (*betw)[3] = NULL;
+	switch (el.GetType())
+	  {
+	  case TRIG6:
+	    {
+	      betw = betw_trig;
+	      onp = 3;
+	      break;
+	    }
+	  case QUAD6: 
+            {
+              betw = betw_quad6;
+              onp = 4;
+              break;
+            }
+	  case QUAD8: 
+            {
+              betw = betw_quad8;
+              onp = 4;
+              break;
+            }
+          default:
+            ;
+          }
+        
+        if (betw)
+          for (int j = 0; j < el.GetNP()-onp; j++)
+            {
+              int pi1 = el[betw[j][0]];
+              int pi2 = el[betw[j][1]];
+              INDEX_2 i2 = INDEX_2::Sort (pi1, pi2);
+              between.Set (i2, el[onp+j]);
+            }
+      }
+      
 
     bool thinlayers = 0;
     for (ElementIndex ei = 0; ei < mesh.GetNE(); ei++)
@@ -28,7 +79,7 @@ namespace netgen
 	thinlayers = 1;
     
 
-    nseg = mesh.GetNSeg();
+    int nseg = mesh.GetNSeg();
     for (SegmentIndex si = 0; si < nseg; si++)
       {
 	Segment & el = mesh.LineSegment(si);
@@ -54,8 +105,7 @@ namespace netgen
       }
 
     // refine surface elements
-    nse = mesh.GetNSE();
+    for (SurfaceElementIndex sei = 0; sei < mesh.GetNSE(); sei++)
-    for (SurfaceElementIndex sei = 0; sei < nse; sei++)
       {
 	int j;
 	const Element2d & el = mesh.SurfaceElement(sei);
@@ -149,8 +199,7 @@ namespace netgen
 
 
     // refine volume elements
-    ne = mesh.GetNE();
+    for (int i = 1; i <= mesh.GetNE(); i++)
-    for (int i = 1; i <= ne; i++)
       {
 	const Element & el = mesh.VolumeElement(i);
 	int onp = 0;

libsrc/stlgeom/stlgeomchart.cpp

@@ -22,6 +22,9 @@ void STLGeometry :: MakeAtlas(Mesh & mesh)
   int timer1 = NgProfiler::CreateTimer ("makeatlas");
   int timer2 = NgProfiler::CreateTimer ("makeatlas - part 2");
   int timer3 = NgProfiler::CreateTimer ("makeatlas - part 3");
+  int timer4 = NgProfiler::CreateTimer ("makeatlas - part 4");
+  int timer4a = NgProfiler::CreateTimer ("makeatlas - part 4a");
+  int timer5 = NgProfiler::CreateTimer ("makeatlas - part 5");
 
   PushStatusF("Make Atlas");
 
@@ -113,6 +116,8 @@ void STLGeometry :: MakeAtlas(Mesh & mesh)
       STLChart * chart = new STLChart(this);
       atlas.Append(chart);
 
+      // *testout << "Chart " << atlas.Size() << endl;
+
       //find unmarked trig
       int prelastunmarked = lastunmarked;
 
@@ -143,7 +148,8 @@ void STLGeometry :: MakeAtlas(Mesh & mesh)
 	  
       Vec<3> sn = GetTriangle(starttrig).Normal();
       chart->SetNormal (startp, sn);
-
+      
+      // *testout << "first trig " << starttrig << ", n = " << sn << endl;
 
       SetMarker(starttrig, chartnum);
       markedtrigcnt++;
@@ -181,14 +187,16 @@ void STLGeometry :: MakeAtlas(Mesh & mesh)
 		  for (int j = 1; j <= NONeighbourTrigs(i); j++)
 		    {
 		      int nt = NeighbourTrig(i,j);
+                      // *testout << "check trig " << nt << endl;
 		      int np1, np2;
 		      GetTriangle(i).GetNeighbourPoints(GetTriangle(nt),np1,np2);
 		      if (GetMarker(nt) == 0 && !IsEdge(np1,np2))
 			{
 			  Vec<3> n2 = GetTriangle(nt).Normal();
+                          // *testout << "acos = " << 180/M_PI*acos (n2*sn) << endl;
 			  if ( (n2 * sn) >= coschartangle )
 			    {
-			      
+                              // *testout << "good angle " << endl;
 			      accepted = 1;
 			      /*
 				//alter spiralentest, schnell, aber ungenau
@@ -229,19 +237,22 @@ void STLGeometry :: MakeAtlas(Mesh & mesh)
 				      accepted = chartbound.TestSeg(GetPoint(nnp1),
 								    GetPoint(nnp2),
 								    sn,sinchartangle,1 /*chartboundarydivisions*/ ,points, eps);
-
+                                      
+                                      // if (!accepted) *testout << "not acc due to testseg" << endl;
 
 				      Vec<3> n3 = GetTriangle(nnt).Normal();
 				      if ( (n3 * sn) >= coschartangle  &&
 					   IsSmoothEdge (nnp1, nnp2) )
 					accepted = 1;
 				    }
-				  if (!accepted) {break;}
+				  if (!accepted) 
+                                    break;
 				}
 			      
-
+                              
 			      if (accepted)
 				{
+                                  // *testout << "trig accepted" << endl;
 				  SetMarker(nt, chartnum); 
 				  changed = true;
 				  markedtrigcnt++;
@@ -318,6 +329,8 @@ void STLGeometry :: MakeAtlas(Mesh & mesh)
 		    {
 		      accepted = 1;
 		      
+                      NgProfiler::StartTimer (timer4);
+
 		      bool isdirtytrig = false;
 		      Vec<3> gn = GetTriangle(nt).GeomNormal(points);
 		      double gnlen = gn.Length();
@@ -339,11 +352,15 @@ void STLGeometry :: MakeAtlas(Mesh & mesh)
 				int nnp1, nnp2; 
 				GetTriangle(nt).GetNeighbourPoints(GetTriangle(nnt),nnp1,nnp2);
 				
+                                // NgProfiler::StartTimer (timer4a);
+
 				accepted = 
 				  chartbound.TestSeg(GetPoint(nnp1),GetPoint(nnp2),
 						     sn,sinouterchartangle, 0 /*chartboundarydivisions*/ ,points, eps);
 				
-				
+                                // NgProfiler::StopTimer (timer4a);				
+
+
 				Vec<3> n3 = GetTriangle(nnt).Normal();
 				if ( (n3 * sn) >= cosouterchartangle  &&
 				     IsSmoothEdge (nnp1, nnp2) )
@@ -352,7 +369,11 @@ void STLGeometry :: MakeAtlas(Mesh & mesh)
 			    if (!accepted) break;
 			  }
 		      
-		      
+                      NgProfiler::StopTimer (timer4);		      
+
+                      NgProfiler::RegionTimer reg5(timer5);		      
+
+
 		      // outer chart is only small environment of
 		      //    inner chart:
 		      

libsrc/stlgeom/stltool.cpp

@@ -894,10 +894,10 @@ void STLBoundary ::AddTriangle(const STLTriangle & t)
 int STLBoundary :: TestSeg(const Point<3>& p1, const Point<3> & p2, const Vec<3> & sn, 
 			   double sinchartangle, int divisions, Array<Point<3> >& points, double eps)
 {
-
   if (usechartnormal)
     return TestSegChartNV (p1, p2, sn);
 
+#ifdef NONE
   // for statistics
   {
     int i;
@@ -936,7 +936,7 @@ int STLBoundary :: TestSeg(const Point<3>& p1, const Point<3> & p2, const Vec<3>
 	*/
       }
   }
-
+#endif
 
   int i,j,k;
   Point<3> seg1p/*, seg2p*/;
@@ -1087,6 +1087,9 @@ int STLBoundary :: TestSeg(const Point<3>& p1, const Point<3> & p2, const Vec<3>
 int STLBoundary :: TestSegChartNV(const Point3d & p1, const Point3d& p2, 
 				  const Vec3d& sn)
 {
+  int timer = NgProfiler::CreateTimer ("TestSegChartNV");
+  NgProfiler::StartTimer (timer);
+
   int nseg = NOSegments();
 
   Point<2> p2d1 = chart->Project2d (p1);
@@ -1103,11 +1106,10 @@ int STLBoundary :: TestSegChartNV(const Point3d & p1, const Point3d& p2,
 
   for (int j = 1; j <= nseg; j++)
     {
-      if (!ok) continue;
       const STLBoundarySeg & seg = GetSegment(j);
 
-      if (!box2d.Intersect (seg.BoundingBox())) continue;
       if (seg.IsSmoothEdge()) continue;
+      if (!box2d.Intersect (seg.BoundingBox())) continue;
 
       const Point<2> & sp1 = seg.P2D1();
       const Point<2> & sp2 = seg.P2D2();
@@ -1115,14 +1117,18 @@ int STLBoundary :: TestSegChartNV(const Point3d & p1, const Point3d& p2,
       Line2d l2 (sp1, sp2);
       double lam1, lam2;
       
-      int err =
+      int err = CrossPointBarycentric (l1, l2, lam1, lam2);
-	CrossPointBarycentric (l1, l2, lam1, lam2);
 
       if (!err && lam1 > eps && lam1 < 1-eps &&
 	  lam2 > eps && lam2 < 1-eps)
-	ok = false;
+        {
+          ok = false;
+          break;
+        }
     }
 
+  NgProfiler::StopTimer (timer);
+
   return ok;
 }