(hopefully) fixed point in quad for the last time

libsrc/meshing/meshclass.cpp

@@ -43,8 +43,8 @@ namespace netgen
     bcnames.SetSize(0);
     cd2names.SetSize(0);
 
+    this->comm = netgen :: ng_comm;
 #ifdef PARALLEL
-    this->comm = MPI_COMM_WORLD;
     paralleltop = new ParallelMeshTopology (*this);
 #endif
   }
@@ -83,6 +83,10 @@ namespace netgen
 #endif
   }
 
+  void Mesh :: SetCommunicator(MPI_Comm acomm)
+  {
+    this->comm = acomm;
+  }
 
   Mesh & Mesh :: operator= (const Mesh & mesh2)
   {
@@ -314,7 +318,7 @@ namespace netgen
     SurfaceElementIndex si = surfelements.Size();
     surfelements.Append (el); 
 
-    if (el.index > facedecoding.Size())
+    if (el.index<=0 || el.index > facedecoding.Size())
       cerr << "has no facedecoding: fd.size = " << facedecoding.Size() << ", ind = " << el.index << endl;
 
     surfelements.Last().next = facedecoding[el.index-1].firstelement;
@@ -1173,7 +1177,7 @@ namespace netgen
 		for (i = 1; i<= GetNSeg(); i++)
 		  {
 		    Segment & seg = LineSegment(i);
-		    if ( seg.cd2i <= n )
+		    if ( seg.edgenr <= n )
 		      seg.SetBCName (cd2names[seg.edgenr-1]);
 		    else
 		      seg.SetBCName(0);
@@ -1304,7 +1308,7 @@ namespace netgen
   }
 
 
-  void Mesh :: DoArchive (ngstd::Archive & archive)
+  void Mesh :: DoArchive (Archive & archive)
   {
     archive & dimension;
     archive & points;
@@ -1316,25 +1320,14 @@ namespace netgen
 
     archive & *ident;
 
-
+    archive.Shallow(geometry);
-    // archive geometry
+    archive & *curvedelems;
-    if (archive.Output())
-      {
-        ostringstream ost;
-        if (geometry)
-          geometry -> SaveToMeshFile (ost);
-        archive << ost.str();
-      }
-    else
-      {
-        string str;
-        archive & str;
-        istringstream ist(str);
-        geometry = geometryregister.LoadFromMeshFile (ist);        
-      }
     
     if (archive.Input())
       {
+	int rank = MyMPI_GetId(GetCommunicator());
+	int ntasks = MyMPI_GetNTasks(GetCommunicator());
+	
         RebuildSurfaceElementLists();
         
         CalcSurfacesOfNode ();
@@ -1687,6 +1680,7 @@ namespace netgen
 
   void Mesh :: CalcSurfacesOfNode ()
   {
+    static Timer t("Mesh::CalcSurfacesOfNode"); RegionTimer reg (t);    
     // surfacesonnode.SetSize (GetNP());
     TABLE<int,PointIndex::BASE> surfacesonnode(GetNP());
 
@@ -1871,8 +1865,7 @@ namespace netgen
 
   void Mesh :: FindOpenElements (int dom)
   {
-    static int timer = NgProfiler::CreateTimer ("Mesh::FindOpenElements");
+    static Timer t("Mesh::FindOpenElements"); RegionTimer reg (t);
-    NgProfiler::RegionTimer reg (timer);
 
     int np = GetNP();
     int ne = GetNE();
@@ -2724,6 +2717,8 @@ namespace netgen
 
   void Mesh :: CalcLocalH (double grading) 
   {
+    static Timer t("Mesh::CalcLocalH"); RegionTimer reg(t);
+    
     if (!lochfunc)
       {
         Point3d pmin, pmax;
@@ -3249,6 +3244,8 @@ namespace netgen
 
   void Mesh :: Compress ()
   {
+    static Timer t("Mesh::Compress"); RegionTimer reg(t);
+    
     Array<PointIndex,PointIndex::BASE,PointIndex> op2np(GetNP());
     Array<MeshPoint> hpoints;
     BitArrayChar<PointIndex::BASE> pused(GetNP());
@@ -3388,7 +3385,7 @@ namespace netgen
 
     for (int i = 0; i < lockedpoints.Size(); i++)
       lockedpoints[i] = op2np[lockedpoints[i]];
-
+    /*
     for (int i = 0; i < facedecoding.Size(); i++)
       facedecoding[i].firstelement = -1;
     for (int i = surfelements.Size()-1; i >= 0; i--)
@@ -3397,7 +3394,8 @@ namespace netgen
         surfelements[i].next = facedecoding[ind-1].firstelement;
         facedecoding[ind-1].firstelement = i;
       }
-
+    */
+    RebuildSurfaceElementLists ();
     CalcSurfacesOfNode();
 
 
@@ -3515,6 +3513,8 @@ namespace netgen
 
   int Mesh :: CheckOverlappingBoundary () 
   {
+    static Timer t("Mesh::CheckOverlappingBoundary"); RegionTimer reg(t);
+    
     int i, j, k;
 
     Point3d pmin, pmax;
@@ -4327,10 +4327,12 @@ namespace netgen
           elementsearchtree = NULL;
           
           int ne = (dimension == 2) ? GetNSE() : GetNE();
+          if (dimension == 3 && !GetNE() && GetNSE())
+            ne = GetNSE();
 
           if (ne) 
             {
-              if (dimension == 2)
+              if (dimension == 2 || (dimension == 3 && !GetNE()) )
                 {
                   Box<3> box (Box<3>::EMPTY_BOX);
                   for (SurfaceElementIndex sei = 0; sei < ne; sei++)
@@ -4394,6 +4396,10 @@ namespace netgen
     return 0;
   }
 
+  bool ValidBarCoord(double lami[3], double eps=1e-12)
+  {
+    return (lami[0]<=1.+eps && lami[0]>=0.-eps && lami[1]<=1.+eps && lami[1]>=0.-eps && lami[2]<=1.+eps && lami[2]>=0.-eps );
+  }
 
   bool Mesh :: PointContainedIn2DElement(const Point3d & p,
                                          double lami[3],
@@ -4424,16 +4430,167 @@ namespace netgen
         Vec3d c = p4 - a; 
         Vec3d d = p3 - a - b - c;
 
+        /*cout << "p = " << p << endl;
+        cout << "p1 = " << p1 << endl;
+        cout << "p2 = " << p2 << endl;
+        cout << "p3 = " << p3 << endl;
+        cout << "p4 = " << p4 << endl;
 
+        cout << "a = " << a << endl;
+        cout << "b = " << b << endl;
+        cout << "c = " << c << endl;
+        cout << "d = " << d << endl;*/
+
+
+        Vec3d pa = p-a;
         double dxb = d.X()*b.Y()-d.Y()*b.X();
         double dxc = d.X()*c.Y()-d.Y()*c.X();
+        double bxc = b.X()*c.Y()-b.Y()*c.X();
+        double bxpa = b.X()*pa.Y()-b.Y()*pa.X();
+        double cxpa = c.X()*pa.Y()-c.Y()*pa.X();
+        double dxpa = d.X()*pa.Y()-d.Y()*pa.X();
+
+        /*cout << "dxb = " << dxb << endl;
+        cout << "dxc = " << dxc << endl;
+        cout << "bxc = " << bxc << endl;
+        cout << "bxpa = " << bxpa << endl;
+        cout << "cxpa = " << cxpa << endl;
+        cout << "dxpa = " << dxpa << endl;*/
+
+        /*
+          P = a + b x + c y + d x y
+          1) P1 = a1 + b1 x + c1 y + d1 x y
+          2) P2 = a2 + b2 x + c2 y + d2 x y
+          
+          -> det(x,d) = det(a,d) + det(b,d) x + det(c,d) y
+            -> x = 1/det(b,d) *( det(P-a,d)-det(c,d) y )
+            -> y = 1/det(c,d) *( det(P-a,d)-det(b,d) x )
+          
+          -> x = (P1 - a1 - c1 y)/(b1 + d1 y)
+            -> det(c,d) y**2 + [det(d,P-a) + det(c,b)] y + det(b,P-a) = 0
+          ( same if we express x = (P2 - a2 - c2 y)/(b2 + d2 y) )
+
+          -> y = (P1 - a1 - b1 x)/(c1 + d1 x)
+            -> det(b,d) x**2 + [det(d,P-a) + det(b,c)] x + det(c,P-a) = 0
+          ( same if we express y = (P2 - a2 - b2 x)/(c2 + d2 x)
+         */
+
+        lami[2]=0.; 
+        double eps = 1.E-12;
+        double c1,c2,r;
+
+        //First check if point is "exactly" a vertex point
+        Vec3d d1 = p-p1;
+        Vec3d d2 = p-p2;
+        Vec3d d3 = p-p3;
+        Vec3d d4 = p-p4;
+
+        //cout << " d1 = " << d1 << ", d2 = " << d2 << ", d3 = " << d3 << ", d4 = " << d4 << endl;
+        
+        if (d1.Length2() < sqr(eps)*d2.Length2() && d1.Length2() < sqr(eps)*d3.Length2() && d1.Length2() < sqr(eps)*d4.Length2())
+          {
+            lami[0] = lami[1] = 0.;
+            return true;
+          }
+        else if (d2.Length2() < sqr(eps)*d1.Length2() && d2.Length2() < sqr(eps)*d3.Length2() && d2.Length2() < sqr(eps)*d4.Length2())
+          {
+            lami[0] = 1.;
+            lami[1] = 0.;
+            return true;
+          }
+        else if (d3.Length2() < sqr(eps)*d1.Length2() && d3.Length2() < sqr(eps)*d2.Length2() && d3.Length2() < sqr(eps)*d4.Length2())
+          {
+            lami[0] = lami[1] = 1.;
+            return true;
+          }
+        else if (d4.Length2() < sqr(eps)*d1.Length2() && d4.Length2() < sqr(eps)*d2.Length2() && d4.Length2() < sqr(eps)*d3.Length2())
+          {
+            lami[0] = 0.;
+            lami[1] = 1.;
+            return true;
+          }//if d is nearly 0: solve resulting linear system
+        else if (d.Length2() < sqr(eps)*b.Length2() && d.Length2() < sqr(eps)*c.Length2())
+          {
+            Vec2d sol;
+            SolveLinearSystemLS (b, c, p-a, sol);
+            lami[0] = sol.X();
+            lami[1] = sol.Y();
+	    return ValidBarCoord(lami, eps);
+          }// if dxc is nearly 0: solve resulting linear equation for y and compute x
+        else if (fabs(dxc) < sqr(eps))
+          {
+            lami[1] = -bxpa/(dxpa-bxc);
+            lami[0] = (dxpa-dxc*lami[1])/dxb;
+            return ValidBarCoord(lami, eps);
+          }// if dxb is nearly 0: solve resulting linear equation for x and compute y
+        else if (fabs(dxb) < sqr(eps))
+          {
+            lami[0] = -cxpa/(dxpa+bxc);
+            lami[1] = (dxpa-dxb*lami[0])/dxc;
+            return ValidBarCoord(lami, eps);
+          }//if dxb >= dxc: solve quadratic equation in y and compute x
+        else if (fabs(dxb) >= fabs(dxc))
+          {
+            c1 = (bxc-dxpa)/dxc;
+            c2 = -bxpa/dxc;
+            r = c1*c1/4.0-c2;
+
+            //quadratic equation has only 1 (unstable) solution
+            if (fabs(r) < eps) //not eps^2!
+              {
+                lami[1] = -c1/2;
+                lami[0] = (dxpa-dxc*lami[1])/dxb;
+                return ValidBarCoord(lami, eps);
+              }
+            if (r < 0) return false;
+
+            lami[1] = -c1/2+sqrt(r);
+            lami[0] = (dxpa-dxc*lami[1])/dxb;
+
+            if (ValidBarCoord(lami, eps))
+                return true;
+            else
+              {
+                lami[1] = -c1/2-sqrt(r);
+                lami[0] = (dxpa-dxc*lami[1])/dxb;
+                return ValidBarCoord(lami, eps);
+              }
+          }//if dxc > dxb: solve quadratic equation in x and compute y
+        else
+          {
+            c1 = (-bxc-dxpa)/dxb;
+            c2 = -cxpa/dxb;
+            r = c1*c1/4.0-c2;
+
+            //quadratic equation has only 1 (unstable) solution
+            if (fabs(r) < eps) //not eps^2!
+              {
+                lami[0] = -c1/2;
+                lami[1] = (dxpa-dxb*lami[0])/dxc;
+                return ValidBarCoord(lami, eps);
+              }
+            if (r < 0) return false;
+
+            lami[0] = -c1/2+sqrt(r);
+            lami[1] = (dxpa-dxb*lami[0])/dxc;
+
+            if (ValidBarCoord(lami, eps))
+                return true;
+            else
+              {
+                lami[0] = -c1/2-sqrt(r);
+                lami[1] = (dxpa-dxb*lami[0])/dxc;
+                return ValidBarCoord(lami, eps);
+              }
+          }
+        
+        /*
         double dxa = d.X()*a.Y()-d.Y()*a.X(); 
         double dxp = d.X()*p.Y()-d.Y()*p.X();
 	
 	
         double c0,c1,c2; // ,rt; 
-        lami[2]=0.; 
+        
-        double eps = 1.E-12;
 
 	Vec3d dp13 = p3-p1;
 	Vec3d dp24 = p4-p2;
@@ -4453,12 +4610,12 @@ namespace netgen
 	    if(lami[1]<=1.+eps && lami[1]>=0.-eps && lami[0]<=1.+eps && lami[0]>=0.-eps)
 	      return true;
 	    
-            /*
+            
-            lami[0]=(c.Y()*(p.X()-a.X())-c.X()*(p.Y()-a.Y()))/
+              //lami[0]=(c.Y()*(p.X()-a.X())-c.X()*(p.Y()-a.Y()))/
-              (b.X()*c.Y() -b.Y()*c.X()); 
+              //(b.X()*c.Y() -b.Y()*c.X()); 
-            lami[1]=(-b.Y()*(p.X()-a.X())+b.X()*(p.Y()-a.Y()))/
+            //lami[1]=(-b.Y()*(p.X()-a.X())+b.X()*(p.Y()-a.Y()))/
-              (b.X()*c.Y() -b.Y()*c.X()); 
+             // (b.X()*c.Y() -b.Y()*c.X()); 
-            */
+            
           } 
         else
           if(fabs(dxb) <= eps*fabs(dxc))
@@ -4552,7 +4709,7 @@ namespace netgen
 
 		if(lami[1]<=1.+eps && lami[1]>=0.-eps && lami[0]<=1.+eps && lami[0]>=0.-eps)
 		  return true;
-	      }
+                  }*/
 
       
         //cout << "lam0,1 = " << lami[0] << ", " << lami[1] << endl;
@@ -4712,7 +4869,7 @@ namespace netgen
                   lam(2) > -eps && 
                   lam(0) + lam(1) + lam(2) < 1+eps);
       }
-    else if (el.GetType() == PRISM)
+    else if (el.GetType() == PRISM || el.GetType() == PRISM15)
       {
         retval = (lam(0) > -eps &&
                   lam(1) > -eps &&
@@ -4720,7 +4877,7 @@ namespace netgen
                   lam(2) < 1+eps &&
                   lam(0) + lam(1) < 1+eps);
       }
-    else if (el.GetType() == PYRAMID)
+    else if (el.GetType() == PYRAMID || el.GetType() == PYRAMID13)
       {
         retval = (lam(0) > -eps &&
                   lam(1) > -eps &&
@@ -4728,7 +4885,7 @@ namespace netgen
                   lam(0) + lam(2) < 1+eps &&
                   lam(1) + lam(2) < 1+eps);
       }
-    else if (el.GetType() == HEX)
+    else if (el.GetType() == HEX || el.GetType() == HEX20)
       {
         retval = (lam(0) > -eps && lam(0) < 1+eps &&
                   lam(1) > -eps && lam(1) < 1+eps &&
@@ -4831,10 +4988,7 @@ namespace netgen
                                  bool build_searchtree,
                                  const bool allowindex) const
   {
-    const double pointtol = 1e-12;
+    if (dimension == 2 || (dimension==3 && !GetNE() && GetNSE()))
-    netgen::Point<3> pmin = p - Vec<3> (pointtol, pointtol, pointtol);
-    netgen::Point<3> pmax = p + Vec<3> (pointtol, pointtol, pointtol);
-    if (dimension == 2)
       {
         int ne;
         int ps_startelement = 0;  // disable global buffering
@@ -4847,7 +5001,10 @@ namespace netgen
           {
             // update if necessary:
             const_cast<Mesh&>(*this).BuildElementSearchTree ();
-            elementsearchtree->GetIntersecting (pmin, pmax, locels);
+            // double tol = elementsearchtree->Tolerance();
+            // netgen::Point<3> pmin = p - Vec<3> (tol, tol, tol);
+            // netgen::Point<3> pmax = p + Vec<3> (tol, tol, tol);
+            elementsearchtree->GetIntersecting (p, p, locels);
             ne = locels.Size();
           }
         else
@@ -4890,7 +5047,10 @@ namespace netgen
           {
             // update if necessary:
             const_cast<Mesh&>(*this).BuildElementSearchTree ();
-            elementsearchtree->GetIntersecting (pmin, pmax, locels);
+            // double tol = elementsearchtree->Tolerance();            
+            // netgen::Point<3> pmin = p - Vec<3> (tol, tol, tol);
+            // netgen::Point<3> pmax = p + Vec<3> (tol, tol, tol);
+            elementsearchtree->GetIntersecting (p, p, locels);
             ne = locels.Size();
           }
         else
@@ -4988,6 +5148,14 @@ namespace netgen
         //(*testout) << "p " << p << endl;
         //(*testout) << "velement " << velement << endl;
 
+        if (!GetNE() && GetNSE() )
+          {
+            lami[0] = vlam[0];
+            lami[1] = vlam[1];
+            lami[2] = vlam[2];
+            return velement;
+          }
+        
         Array<int> faces;
         topology.GetElementFaces(velement,faces);
 
@@ -6030,7 +6198,7 @@ namespace netgen
   void Mesh :: SetUserData(const char * id, Array<int> & data)
   {
     if(userdata_int.Used(id))
-      delete userdata_int.Get(id);
+      delete userdata_int[id];
 
     Array<int> * newdata = new Array<int>(data);
 
@@ -6040,10 +6208,10 @@ namespace netgen
   {
     if(userdata_int.Used(id))
       {
-        if(data.Size() < (*userdata_int.Get(id)).Size()+shift)
+        if(data.Size() < (*userdata_int[id]).Size()+shift)
-          data.SetSize((*userdata_int.Get(id)).Size()+shift);
+          data.SetSize((*userdata_int[id]).Size()+shift);
-        for(int i=0; i<(*userdata_int.Get(id)).Size(); i++)
+        for(int i=0; i<(*userdata_int[id]).Size(); i++)
-          data[i+shift] = (*userdata_int.Get(id))[i];
+          data[i+shift] = (*userdata_int[id])[i];
         return true;
       }
     else
@@ -6055,7 +6223,7 @@ namespace netgen
   void Mesh :: SetUserData(const char * id, Array<double> & data)
   {
     if(userdata_double.Used(id))
-      delete userdata_double.Get(id);
+      delete userdata_double[id];
 
     Array<double> * newdata = new Array<double>(data);
 
@@ -6065,10 +6233,10 @@ namespace netgen
   {
     if(userdata_double.Used(id))
       {
-        if(data.Size() < (*userdata_double.Get(id)).Size()+shift)
+        if(data.Size() < (*userdata_double[id]).Size()+shift)
-          data.SetSize((*userdata_double.Get(id)).Size()+shift);
+          data.SetSize((*userdata_double[id]).Size()+shift);
-        for(int i=0; i<(*userdata_double.Get(id)).Size(); i++)
+        for(int i=0; i<(*userdata_double[id]).Size(); i++)
-          data[i+shift] = (*userdata_double.Get(id))[i];
+          data[i+shift] = (*userdata_double[id])[i];
         return true;
       }
     else