implemented jacobianInverseTransposed for 1d-in-2d and 2d-triangles-in-3d

[[Imported from SVN: r4188]]

grid/uggrid/uggridgeometry.cc

@@ -305,6 +305,14 @@ integrationElement (const Dune::FieldVector<typename GridImp::ctype, 1>& local)
   return diff.two_norm();
 }
 
+template <class GridImp>
+const Dune::FieldMatrix<typename GridImp::ctype,1,1>& Dune::UGGridGeometry<1,2,GridImp>::
+jacobianInverseTransposed(const Dune::FieldVector<typename GridImp::ctype, 1>& local) const
+{
+  jacobianInverseTransposed_[0][0] = 1 / (coord_[0]-coord_[1]).two_norm();
+  return jacobianInverseTransposed_;
+}
+
 template <class GridImp>
 inline typename GridImp::ctype Dune::UGGridGeometry<2,3,GridImp>::
 integrationElement (const Dune::FieldVector<typename GridImp::ctype, 2>& local) const
@@ -314,3 +322,39 @@ integrationElement (const Dune::FieldVector<typename GridImp::ctype, 2>& local)
   /** \todo Maybe there should be a test for this */
   return UG_NS<3>::SurfaceElement(corners(), (const double(*)[3])&coord_,&local[0]);
 }
+
+template <class GridImp>
+const Dune::FieldMatrix<typename GridImp::ctype,2,2>& Dune::UGGridGeometry<2,3,GridImp>::
+jacobianInverseTransposed(const Dune::FieldVector<typename GridImp::ctype, 2>& local) const
+{
+  // I don't really know how to implement this for quadrilateral faces,
+  // especially since they may be nonplanar.
+  if (!type().isTriangle())
+    DUNE_THROW(NotImplemented, "jacobianInverse only implemented for triangular faces!");
+
+  // The spatial triangle is first mapped isometrically onto the plane.  We map
+  // the first vertex onto the origin, the second one on the positive x-axis,
+  // and the third one such that is has positive y-coordinate.  Then we call
+  // the UG-routine for planar triangles.  This is certainly not the most elegant
+  // way, but the first one that comes to my mind.
+
+  double l0 = (coord_[2]-coord_[1]).two_norm();
+  double l1 = (coord_[2]-coord_[0]).two_norm();
+  double l2 = (coord_[1]-coord_[0]).two_norm();
+
+  double q0 = (l2*l2 - l0*l0 + l1*l1) / (2*l2);
+  double h  = sqrt(l1*l1 - q0*q0);
+
+  FieldVector<double,2> p0(0);
+  FieldVector<double,2> p1(0);    p1[0] = l2;
+  FieldVector<double,2> p2(0);    p2[0] = q0;   p2[1] = h;
+
+  // Check that this was really an isometry
+  assert( fabs(p2.two_norm()      - l1)      < 1e-6 );
+  assert( fabs((p2-p1).two_norm() - l0) < 1e-6 );
+
+  double* cornerCoords[3] = {&p0[0], &p1[0], &p2[0]};
+
+  UG_NS<2>::Transformation(3, cornerCoords, local, jacobianInverseTransposed_);
+  return jacobianInverseTransposed_;
+}

grid/uggrid/uggridgeometry.hh

@@ -334,7 +334,7 @@ namespace Dune {
     mutable FixedArray<FieldVector<UGCtype, 3>, 4> coord_;
 
     //! The jacobian inverse
-    mutable FieldMatrix<UGCtype,3,3> jac_inverse_;
+    mutable FieldMatrix<UGCtype,2,2> jacobianInverseTransposed_;
 
   };
 
@@ -416,6 +416,9 @@ namespace Dune {
     //! the vertex coordinates
     FixedArray<FieldVector<UGCtype, 2>, 2> coord_;
 
+    //! The jacobian inverse
+    mutable FieldMatrix<UGCtype,1,1> jacobianInverseTransposed_;
+
   };
 
 }  // namespace Dune