Multi-mesh reimplemented for linear elements in 3d.

AMDiS/src/ElInfo2d.cc

@@ -32,6 +32,33 @@ namespace AMDiS {
   mtl::dense2D<double> ElInfo2d::mat_d1_right(mat_d1_right_val);
 
 
+
+
+  double ElInfo2d::mat_d2_val[6][6] = {{1.0, 0.0, 0.0, 0.0, 0.0, 0.0}, 
+				       {0.0, 1.0, 0.0, 0.0, 0.0, 0.0}, 
+				       {0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
+				       {0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
+				       {0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
+				       {0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
+  mtl::dense2D<double> ElInfo2d::mat_d2(mat_d2_val);
+
+  double ElInfo2d::mat_d2_left_val[6][6] = {{0.0, 1.0, 0.0, 0.375, -0.125, 0.0}, 
+					    {0.0, 0.0, 0.0, -0.125, -0.125, 0.0}, 
+					    {1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
+					    {0.0, 0.0, 0.0, 0.0, 0.5, 0.0},
+					    {0.0, 0.0, 0.0, 0.0, 0.5, 1.0},
+					    {0.0, 0.0, 1.0, 0.75, 0.25, 0.0}};
+  mtl::dense2D<double> ElInfo2d::mat_d2_left(mat_d2_left_val);
+
+  double ElInfo2d::mat_d2_right_val[6][6] = {{0.0, 0.0, 0.0, -0.125, -0.125, 0.0}, 
+					     {1.0, 0.0, 0.0, -0.125, 0.375, 0.0}, 
+					     {0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
+					     {0.0, 0.0, 0.0, 0.5, 0.0, 1.0},
+					     {0.0, 0.0, 0.0, 0.5, 0.0, 0.0},
+					     {0.0, 0.0, 1.0, 0.25, 0.75, 0.0}};
+  mtl::dense2D<double> ElInfo2d::mat_d2_right(mat_d2_right_val);
+
+
   ElInfo2d::ElInfo2d(Mesh *aMesh) 
     : ElInfo(aMesh) 
   {
@@ -714,6 +741,24 @@ namespace AMDiS {
 	  subElemMatrices[1][refinementPath] = mat;  
 	}
 	break;
+      case 2:
+	{
+	  dense2D<double> mat(mat_d2);
+	  dense2D<double> tmpMat(num_rows(mat), num_rows(mat));
+	  
+	  for (int i = 0; i < refinementPathLength; i++) {
+	    if (refinementPath & (1 << i)) {
+	      tmpMat = mat_d2_right * mat;
+	      mat = tmpMat;
+	    } else  {
+	      tmpMat = mat_d2_left * mat;
+	      mat = tmpMat;
+	    }
+	  }
+
+	  subElemMatrices[2][refinementPath] = mat;  
+	}
+	break;
       default:
 	ERROR_EXIT("Not supported for basis function degree: %d\n", degree);
       }

AMDiS/src/ElInfo2d.h

@@ -67,16 +67,22 @@ namespace AMDiS {
     WorldVector<double> *e1, *e2, *normal;
 
     static double mat_d1_val[3][3];
-
     static mtl::dense2D<double> mat_d1;
 
     static double mat_d1_left_val[3][3];
-
     static mtl::dense2D<double> mat_d1_left;
 
     static double mat_d1_right_val[3][3];
-
     static mtl::dense2D<double> mat_d1_right;
+
+    static double mat_d2_val[6][6];
+    static mtl::dense2D<double> mat_d2;
+
+    static double mat_d2_left_val[6][6];
+    static mtl::dense2D<double> mat_d2_left;
+
+    static double mat_d2_right_val[6][6];
+    static mtl::dense2D<double> mat_d2_right;
   };
 
 }

AMDiS/src/ElInfo3d.cc

@@ -14,6 +14,25 @@
 
 namespace AMDiS {
 
+  double ElInfo3d::mat_d1_left_val[4][4] = {{1.0, 0.0, 0.0, 0.5}, 
+					    {0.0, 0.0, 0.0, 0.5},
+					    {0.0, 1.0, 0.0, 0.0},
+					    {0.0, 0.0, 1.0, 0.0}};
+  mtl::dense2D<double> ElInfo3d::mat_d1_left(mat_d1_left_val);
+
+  
+  double ElInfo3d::mat_d1_l0_right_val[4][4] = {{0.0, 0.0, 0.0, 0.5}, 
+						{1.0, 0.0, 0.0, 0.5},
+						{0.0, 0.0, 1.0, 0.0},
+						{0.0, 1.0, 0.0, 0.0}};
+  double ElInfo3d::mat_d1_l12_right_val[4][4] = {{0.0, 0.0, 0.0, 0.5}, 
+						 {1.0, 0.0, 0.0, 0.5},
+						 {0.0, 1.0, 0.0, 0.0},
+						 {0.0, 0.0, 1.0, 0.0}};
+  mtl::dense2D<double> ElInfo3d::mat_d1_l0_right(mat_d1_l0_right_val);
+  mtl::dense2D<double> ElInfo3d::mat_d1_l12_right(mat_d1_l12_right_val);
+
+
   void ElInfo3d::fillMacroInfo(const MacroElement * mel)
   {
     FUNCNAME("ElInfo3d::fillMacroInfo()");
@@ -104,6 +123,7 @@ namespace AMDiS {
     }
   }
 
+
   double ElInfo3d::calcGrdLambda(DimVec<WorldVector<double> >& grd_lam)
   {
     FUNCNAME("ElInfo3d::calcGrdLambda()");
@@ -169,6 +189,7 @@ namespace AMDiS {
     return adet;
   }
 
+
   const int ElInfo3d::worldToCoord(const WorldVector<double>& xy,
 				   DimVec<double>* lambda) const
   {
@@ -343,8 +364,6 @@ namespace AMDiS {
   }
 
 
-
-
   void ElInfo3d::fillElInfo(int ichild, const ElInfo *elInfoOld)
   {
     FUNCNAME("ElInfo3d::fillElInfo()");
@@ -595,7 +614,35 @@ namespace AMDiS {
   {
     FUNCNAME("ElInfo3d::getSubElemCoordsMat()");
 
-    ERROR_EXIT("Not yet implemented!\n");
+    using namespace mtl;
+
+    if (subElemMatrices[degree].count(refinementPath) == 0) {
+      switch (degree) {
+      case 1:
+	{
+	  dense2D<double> mat(4, 4), tmpMat(4, 4);
+	  mat = 1;
+	  
+	  for (int i = 0; i < refinementPathLength; i++) {
+	    if (refinementPath & (1 << i)) {
+	      if ((level + i) % 3 == 0)
+		tmpMat = mat_d1_l0_right * mat;
+	      else
+		tmpMat = mat_d1_l12_right * mat;
+
+	      mat = tmpMat;
+	    } else  {
+	      tmpMat = mat_d1_left * mat;
+	      mat = tmpMat;
+	    }
+	  }
+
+	  subElemMatrices[1][refinementPath] = mat;  
+	}
+      default:
+	ERROR_EXIT("Not supported for basis function degree: %d\n", degree);
+      }
+    }
 
     return subElemMatrices[degree][refinementPath];
   }

AMDiS/src/ElInfo3d.h

@@ -93,6 +93,14 @@ namespace AMDiS {
 
     /// Tmp vectors used for calculations in calcGrdLambda and getNormal().
     std::vector< WorldVector<double> > tmpWorldVecs;
+
+    static double mat_d1_left_val[4][4];
+    static mtl::dense2D<double> mat_d1_left;
+
+    static double mat_d1_l0_right_val[4][4];
+    static double mat_d1_l12_right_val[4][4];
+    static mtl::dense2D<double> mat_d1_l0_right;
+    static mtl::dense2D<double> mat_d1_l12_right;
   };
 
 }