curved elements access functions

2d meshing cleaning (a bit)

libsrc/meshing/curvedelems.cpp

@@ -2743,11 +2743,61 @@ namespace netgen
     ;
   }
 
+
+  template <int DIM_SPACE>
+  void CurvedElements :: 
+  CalcMultiPointSegmentTransformation (SegmentIndex elnr, int n,
+				       const double * xi, size_t sxi,
+				       double * x, size_t sx,
+				       double * dxdxi, size_t sdxdxi)
+  {
+    for (int ip = 0; ip < n; ip++)
+      {
+        Point<3> xg;
+        Vec<3> dx;
+
+        // mesh->GetCurvedElements().
+	CalcSegmentTransformation (xi[ip*sxi], elnr, xg, dx);
+      
+        if (x)
+          for (int i = 0; i < DIM_SPACE; i++)
+            x[ip*sx+i] = xg(i);
+	  
+        if (dxdxi)
+          for (int i=0; i<DIM_SPACE; i++)
+            dxdxi[ip*sdxdxi+i] = dx(i);
+      }
+  }
+
+  template void CurvedElements :: 
+  CalcMultiPointSegmentTransformation<2> (SegmentIndex elnr, int npts,
+                                          const double * xi, size_t sxi,
+                                          double * x, size_t sx,
+                                          double * dxdxi, size_t sdxdxi);
+
+  template void CurvedElements :: 
+  CalcMultiPointSegmentTransformation<3> (SegmentIndex elnr, int npts,
+                                          const double * xi, size_t sxi,
+                                          double * x, size_t sx,
+                                          double * dxdxi, size_t sdxdxi);
+
+
+
   void CurvedElements :: 
   CalcMultiPointSurfaceTransformation (Array< Point<2> > * xi, SurfaceElementIndex elnr,
 				       Array< Point<3> > * x,
 				       Array< Mat<3,2> > * dxdxi)
   {
+    double * px = (x) ? &(*x)[0](0) : NULL;
+    double * pdxdxi = (dxdxi) ? &(*dxdxi)[0](0) : NULL;
+
+    CalcMultiPointSurfaceTransformation <3> (elnr, xi->Size(),
+					     &(*xi)[0](0), 2, 
+					     px, 3,
+					     pdxdxi, 6);
+					    
+    return;
+#ifdef OLD
     if (mesh.coarsemesh)
       {
 	const HPRefElement & hpref_el =
@@ -2876,6 +2926,7 @@ namespace netgen
             (*dxdxi)[ip] = ds;
 	  }
       }
+#endif
   }
 
 
@@ -2918,8 +2969,7 @@ namespace netgen
                                                           &coarse_xi[0](0), &coarse_xi[1](0)-&coarse_xi[0](0), 
                                                           x, sx, dxdxi, sdxdxi);
 
-        Mat<2,2> trans;
+        // Mat<3,2> dxdxic;
-        Mat<3,2> dxdxic;
 	if (dxdxi)
 	  {
 	    MatrixFixWidth<2> dlami(4);
@@ -2930,6 +2980,7 @@ namespace netgen
                 Point<2> hxi(xi[pi*sxi], xi[pi*sxi+1]);
 		mesh[elnr].GetDShapeNew ( hxi, dlami);	  
 		
+		Mat<2,2> trans;
 		trans = 0;
 		for (int k = 0; k < 2; k++)
 		  for (int l = 0; l < 2; l++)
@@ -3052,12 +3103,22 @@ namespace netgen
 
 
 
-
   void CurvedElements :: 
   CalcMultiPointElementTransformation (Array< Point<3> > * xi, ElementIndex elnr,
 				       Array< Point<3> > * x,
 				       Array< Mat<3,3> > * dxdxi)
   {
+    double * px = (x) ? &(*x)[0](0) : NULL;
+    double * pdxdxi = (dxdxi) ? &(*dxdxi)[0](0) : NULL;
+
+    CalcMultiPointElementTransformation (elnr, xi->Size(),
+					 &(*xi)[0](0), 3, 
+					 px, 3,
+					 pdxdxi, 9);
+    
+    return;
+#ifdef OLD
+
     if (mesh.coarsemesh)
       {
 	const HPRefElement & hpref_el =
@@ -3193,11 +3254,11 @@ namespace netgen
 	      (*dxdxi)[ip] = ds;
 	    }
       }
+#endif
   }
 
 
 
-
   void  CurvedElements :: 
   CalcMultiPointElementTransformation (ElementIndex elnr, int n,
                                        const double * xi, size_t sxi,

libsrc/meshing/curvedelems.hpp

@@ -119,6 +119,13 @@ public:
 					    Array<Point<3> > * x,
 					    Array<Vec<3> > * dxdxi);
 
+  template <int DIM_SPACE>
+  void CalcMultiPointSegmentTransformation (SegmentIndex elnr, int n,
+                                            const double * xi, size_t sxi,
+                                            double * x, size_t sx,
+                                            double * dxdxi, size_t sdxdxi);
+
+
   void CalcMultiPointSurfaceTransformation (Array< Point<2> > * xi, SurfaceElementIndex elnr,
 					    Array< Point<3> > * x,
 					    Array< Mat<3,2> > * dxdxi);
@@ -129,7 +136,6 @@ public:
                                             double * x, size_t sx,
                                             double * dxdxi, size_t sdxdxi);
 
-
   void CalcMultiPointElementTransformation (Array< Point<3> > * xi, ElementIndex elnr,
 					    Array< Point<3> > * x,
 					    Array< Mat<3,3> > * dxdxi);

libsrc/meshing/netrule2.cpp

@@ -13,24 +13,14 @@ netrule :: netrule ()
 
 netrule ::  ~netrule()
 {
-  // if(name != NULL) 
   delete [] name;
-  for(int i=0; i<oldutofreearea_i.Size(); i++)
+  for(int i = 0; i < oldutofreearea_i.Size(); i++)
     delete oldutofreearea_i[i];
+  for(int i = 0; i < freezone_i.Size(); i++)
+    delete freezone_i[i];
 }
 
 
-/*
-void netrule :: GetFreeArea (Array<Point2d> & afreearea)
-  {
-  int i;
-
-  afreearea.SetSize (freearea.Size());
-  for (i = 1; i <= freearea.Size(); i++)
-    afreearea[i] = freearea[i];
-  }
-*/
-
 
 void netrule :: SetFreeZoneTransformation (const Vector & devp, int tolclass)
 {
@@ -41,37 +31,46 @@ void netrule :: SetFreeZoneTransformation (const Vector & devp, int tolclass)
 
   int vs = oldutofreearea.Height();
   FlatVector devfree(vs, mem1);
-  FlatVector devfree1(vs, mem2);
-  FlatVector devfree2(vs, mem3);
-
-  if (tolclass <= oldutofreearea_i.Size())
-    {
-      oldutofreearea_i[tolclass-1] -> Mult (devp, devfree);
-    }
-  else
-    {
-      oldutofreearea.Mult (devp, devfree1);
-      oldutofreearealimit.Mult (devp, devfree2);
-      devfree.Set2 (lam1, devfree1, lam2, devfree2);
-    }
-
 
   int fzs = freezone.Size();
   transfreezone.SetSize (fzs);
 
+  if (tolclass <= oldutofreearea_i.Size())
+    {
+      oldutofreearea_i[tolclass-1] -> Mult (devp, devfree);
+
+      Array<Point2d> & fzi = *freezone_i[tolclass-1];
+      for (int i = 0; i < fzs; i++)
+	{
+	  transfreezone[i].X() = fzi[i].X() + devfree[2*i];
+	  transfreezone[i].Y() = fzi[i].Y() + devfree[2*i+1];
+	}
+    }
+  else
+    {
+      FlatVector devfree1(vs, mem2);
+      FlatVector devfree2(vs, mem3);
+
+      oldutofreearea.Mult (devp, devfree1);
+      oldutofreearealimit.Mult (devp, devfree2);
+      devfree.Set2 (lam1, devfree1, lam2, devfree2);
+
+      for (int i = 0; i < fzs; i++)
+	{
+	  transfreezone[i].X() = lam1 * freezone[i].X() + lam2 * freezonelimit[i].X() + devfree[2*i];
+	  transfreezone[i].Y() = lam1 * freezone[i].Y() + lam2 * freezonelimit[i].Y() + devfree[2*i+1];
+	}
+    }
+
+
   if (fzs > 0)
     {
-      transfreezone[0].X() = lam1 * freezone[0].X() + lam2 * freezonelimit[0].X() + devfree[0];
-      transfreezone[0].Y() = lam1 * freezone[0].Y() + lam2 * freezonelimit[0].Y() + devfree[1];
       fzmaxx = fzminx = transfreezone[0].X();
       fzmaxy = fzminy = transfreezone[0].Y();
     }
 
   for (int i = 1; i < fzs; i++)
     {
-      transfreezone[i].X() = lam1 * freezone[i].X() + lam2 * freezonelimit[i].X() + devfree[2*i];
-      transfreezone[i].Y() = lam1 * freezone[i].Y() + lam2 * freezonelimit[i].Y() + devfree[2*i+1];
-
       if (transfreezone[i].X() > fzmaxx) fzmaxx = transfreezone[i].X();
       if (transfreezone[i].X() < fzminx) fzminx = transfreezone[i].X();
       if (transfreezone[i].Y() > fzmaxy) fzmaxy = transfreezone[i].Y();

libsrc/meshing/parser2.cpp

@@ -407,8 +407,6 @@ void netrule :: LoadRule (istream & ist)
     }
   while (!ist.eof() && strcmp (buf, "endrule") != 0);
 
-  //(*testout) << "loadr1" << endl;
-
   oldutonewu.SetSize (2 * (points.Size() - noldp), 2 * noldp);
   oldutofreearea.SetSize (2 * freezone.Size(), 2 * noldp);
   oldutofreearealimit.SetSize (2 * freezone.Size(), 2 * noldp);
@@ -428,8 +426,6 @@ void netrule :: LoadRule (istream & ist)
   freesetinequ.SetSize (freezone.Size());
 
 
-  //(*testout) << "loadr2" << endl;
-
   {
     char ok;
     int minn;
@@ -441,7 +437,6 @@ void netrule :: LoadRule (istream & ist)
     for (j = 1; j <= 2; j++)
       pnearness.Elem(GetPointNr (1, j)) = 0;
 
-    //(*testout) << "loadr3" << endl;
     do
       {
 	ok = 1;
@@ -461,7 +456,6 @@ void netrule :: LoadRule (istream & ist)
 	  }
       }
     while (!ok);
-    //(*testout) << "loadr4" << endl;
 
     lnearness.SetSize (noldl);
 
@@ -472,19 +466,25 @@ void netrule :: LoadRule (istream & ist)
 	  lnearness.Elem(i) += pnearness.Get(GetPointNr (i, j));
       }
   }
-  //(*testout) << "loadr5" << endl;
 
   oldutofreearea_i.SetSize (10);
+  freezone_i.SetSize (10);
+
   for (i = 0; i < oldutofreearea_i.Size(); i++)
     {
+      double lam1 = 1.0/(i+1);
+
       oldutofreearea_i[i] = new DenseMatrix (oldutofreearea.Height(), oldutofreearea.Width());
       DenseMatrix & mati = *oldutofreearea_i[i];
       for (j = 0; j < oldutofreearea.Height(); j++)
 	for (int k = 0; k < oldutofreearea.Width(); k++)
-	  mati(j,k) = 1.0 / (i+1) * oldutofreearea(j,k) + (1 - 1.0/(i+1)) * oldutofreearealimit(j,k);
+	  mati(j,k) = lam1 * oldutofreearea(j,k) + (1 - lam1) * oldutofreearealimit(j,k);
-    }
 
-  //(*testout) << "loadr6" << endl;
+      freezone_i[i] = new Array<Point2d> (freezone.Size());
+      Array<Point2d> & fzi = *freezone_i[i];
+      for (int j = 0; j < freezone.Size(); j++)
+	fzi[j] = freezonelimit[j] + lam1 * (freezone[j] - freezonelimit[j]);
+    }
 }
 
 

libsrc/meshing/ruler2.hpp

@@ -29,6 +29,8 @@ private:
   ///
   Array<Point2d> freezone, freezonelimit;
   ///
+  Array<Array<Point2d>*> freezone_i;
+  ///
   Array<Point2d> transfreezone;
 
   ///

ng/nginterface_v2.cpp

@@ -216,6 +216,8 @@ namespace netgen
                                                       double * x, size_t sx,
                                                       double * dxdxi, size_t sdxdxi)
   {
+    mesh->GetCurvedElements().CalcMultiPointSegmentTransformation<2> (elnr, npts, xi, sxi, x, sx, dxdxi, sdxdxi);
+    /*
     for (int ip = 0; ip < npts; ip++)
       {
         Point<3> xg;
@@ -231,6 +233,7 @@ namespace netgen
           for (int i=0; i<2; i++)
             dxdxi[ip*sdxdxi+i] = dx(i);
       }
+    */
   }
 
   template <>