all files for serial version of ALU3dGrid.

git-svn-id: https://dune.mathematik.uni-freiburg.de/svn/alugrid/trunk@148 0d966ed9-3843-0410-af09-ebfb50bd7c74

src/serial/alloca.h

+#include <malloc.h>

src/serial/gitter_impl.h

+// (c) bernhard schupp 1997 - 1998
+
+#ifndef GITTER_IMPL_H_INCLUDED
+#define GITTER_IMPL_H_INCLUDED
+
+#ifdef IBM_XLC
+  #define _ANSI_HEADER
+#endif
+
+#ifdef _ANSI_HEADER
+  using namespace std;
+  #include <fstream>
+  #include <vector>
+  #include <utility>
+#else
+  #include <fstream.h>
+  #include <vector.h>
+  #include <pair.h>
+#endif
+
+#include "gitter_sti.h"
+
+#include "mapp_tetra_3d.h"  
+
+#include "gitter_hexa_top.h"
+#include "gitter_tetra_top.h"
+
+//static ofstream logFile ("logfile");
+
+class GitterBasis : public virtual Gitter, public Gitter :: Geometric {
+  public :
+    class Objects {
+      public :
+        class VertexEmpty : public VertexGeo {
+          public :
+            inline VertexEmpty (int, double, double, double,
+                IndexManagerType &im) ;
+            inline VertexEmpty (int, double, double, double,
+                VertexGeo & ) ;
+     ~VertexEmpty () {}
+            virtual inline int ident () const ;
+        } ;
+
+        class VertexEmptyMacro : public VertexEmpty {
+          public :
+            inline VertexEmptyMacro (double, double, double, int,
+                IndexManagerType &im) ;
+     ~VertexEmptyMacro () {}
+            virtual inline int ident () const ;
+          private :
+            int _idn ;
+        } ;
+
+      // organizes the indices for boundary faces and the opposite vertices for
+      // ghost cells 
+      template <int pdimvx>
+      class Dune_hbndDefault 
+      {
+        protected:
+          enum { dimvx = pdimvx }; 
+          int    _index;
+          double _oppVx[dimvx][3];
+
+        public:
+          inline Dune_hbndDefault ();
+          inline void setOppPoint (int i, const double (&p)[3]);
+          inline const double (& oppositeVertex (int i) const) [3];
+
+          inline int getIndex () const; 
+          inline void setIndex ( int idx ); 
+          inline int dimVx () const; 
+
+        protected:  
+          inline void splitGhost () {} 
+          inline void setGhost (Gitter::helement_STI *) {}
+      };
+       
+        class Hbnd3Default : public hbndseg3_GEO 
+#ifdef _DUNE_USES_BSGRID_
+           , public Dune_hbndDefault<1>
+#endif
+        {
+          protected :
+           inline Hbnd3Default (myhface3_t *, int, ProjectVertex * ) ;
+           virtual ~Hbnd3Default () {}
+          public :
+            typedef hbndseg3_GEO :: bnd_t bnd_t;
+            virtual inline bnd_t bndtype () const ;
+            virtual int ghostLevel () const ; 
+            
+            // default implementation is doing nothing for these 3 methods
+            // these methods are overloades just on HbndPll
+            virtual helement_STI * getGhost () { return 0; }
+
+            // points inside ghosts 
+            void faceNormal( double * normal) const;
+        };
+        typedef Hbnd3Top < Hbnd3Default > hbndseg3_IMPL ;
+
+        class Hbnd4Default : public hbndseg4_GEO 
+#ifdef _DUNE_USES_BSGRID_
+           , public Dune_hbndDefault<4>
+#endif
+        {
+          protected :
+            inline Hbnd4Default (myhface4_t *, int, ProjectVertex *) ;
+            virtual ~Hbnd4Default () {}
+          public :
+            typedef hbndseg4_GEO :: bnd_t bnd_t;
+            virtual inline bnd_t bndtype () const ;
+            virtual int ghostLevel () const ;  
+
+            // default implementation is doing nothing for these 3 methods
+            // these methods are overloades just on HbndPll
+            virtual helement_STI * getGhost () { return 0; }
+
+            // points inside ghosts 
+            void faceNormal( double * normal) const;
+        };
+        typedef Hbnd4Top < Hbnd4Default > hbndseg4_IMPL ;
+
+        class Hedge1Empty : public hedge1_GEO {
+          protected :
+            typedef VertexEmpty innervertex_t ;
+            inline Hedge1Empty (myvertex_t *,myvertex_t *) ;
+     ~Hedge1Empty () {}
+           // Methode um einen Vertex zu verschieben; f"ur die Randanpassung
+           virtual inline void projectInnerVertex(const ProjectVertex &pv) ; 
+        } ;
+  
+        typedef Hedge1Top < Hedge1Empty > hedge1_IMPL ;
+
+        class Hface3Empty : public hface3_GEO {
+          protected :
+            typedef VertexEmpty   innervertex_t ;
+            typedef hedge1_IMPL   inneredge_t ;
+            inline Hface3Empty (myhedge1_t *,int, myhedge1_t *,int, myhedge1_t *,int) ;
+     ~Hface3Empty () {}
+           // Methode um einen Vertex zu verschieben; f"ur die Randanpassung
+           virtual inline void projectVertex(const ProjectVertex &pv) ; 
+  } ;
+        typedef Hface3Top < Hface3Empty > hface3_IMPL ;
+
+        class Hface4Empty : public hface4_GEO {
+          protected :
+            typedef VertexEmpty innervertex_t ;
+            typedef hedge1_IMPL     inneredge_t ;
+            inline Hface4Empty (myhedge1_t *,int, myhedge1_t *,int, myhedge1_t *,int,myhedge1_t *,int) ;
+     ~Hface4Empty () {}
+           // Methode um einen Vertex zu verschieben; f"ur die Randanpassung
+           virtual inline void projectVertex(const ProjectVertex &pv) ; 
+        } ;
+        typedef Hface4Top < Hface4Empty > hface4_IMPL ;
+
+        class TetraEmpty : public tetra_GEO {
+    protected :
+            typedef hface3_IMPL innerface_t ;
+            typedef hedge1_IMPL inneredge_t ;
+            typedef VertexEmpty innervertex_t ;
+            inline TetraEmpty (myhface3_t *,int,myhface3_t *,int,myhface3_t *,int,myhface3_t *,int) ;
+           ~TetraEmpty () {}
+    public :
+  } ;
+  typedef TetraTop < TetraEmpty > tetra_IMPL ;
+  
+  class Periodic3Empty : public periodic3_GEO {
+    protected :
+            typedef hface3_IMPL innerface_t ;
+            typedef hedge1_IMPL inneredge_t ;
+            typedef VertexEmpty innervertex_t ;
+      inline Periodic3Empty (myhface3_t *,int,myhface3_t *,int) ;
+     ~Periodic3Empty () {}
+    public:
+  } ;
+  typedef Periodic3Top < Periodic3Empty > periodic3_IMPL ;
+
+// Anfang - Neu am 23.5.02 (BS)
+  class Periodic4Empty : public periodic4_GEO {
+    protected :
+            typedef hface4_IMPL innerface_t ;
+            typedef hedge1_IMPL inneredge_t ;
+            typedef VertexEmpty innervertex_t ;
+      inline Periodic4Empty (myhface4_t *,int,myhface4_t *,int) ;
+     ~Periodic4Empty () {}
+    public:
+  } ;
+  typedef Periodic4Top < Periodic4Empty > periodic4_IMPL ;
+// Ende - Neu am 23.5.02 (BS)
+
+        class HexaEmpty : public hexa_GEO {
+          protected :
+            typedef hface4_IMPL innerface_t ;
+            typedef hedge1_IMPL inneredge_t ;
+            typedef VertexEmpty innervertex_t ;
+            inline HexaEmpty (myhface4_t *,int,myhface4_t *,int,myhface4_t *,int,myhface4_t *,int,myhface4_t *,int,myhface4_t *,int) ;
+           ~HexaEmpty () {}
+        } ;
+        typedef HexaTop < HexaEmpty > hexa_IMPL ;
+    } ;
+  public :
+    class MacroGitterBasis : public virtual BuilderIF {
+      protected :
+        virtual inline VertexGeo     * insert_vertex (double, double, double, int,int = 0) ;
+        virtual inline hedge1_GEO    * insert_hedge1 (VertexGeo *, VertexGeo *) ;
+  virtual inline hface3_GEO    * insert_hface3 (hedge1_GEO *(&)[3], int (&)[3]) ;
+        virtual inline hface4_GEO    * insert_hface4 (hedge1_GEO *(&)[4], int (&)[4]) ;
+  virtual inline hbndseg3_GEO  * insert_hbnd3 (hface3_GEO *, int, Gitter :: hbndseg_STI :: bnd_t) ;
+  // version with point , returns insert_hbnd3 here 
+  virtual inline hbndseg3_GEO  * insert_hbnd3 (hface3_GEO *, int, Gitter :: hbndseg_STI :: bnd_t, const double (&p)[3]) ;
+        virtual inline hbndseg4_GEO  * insert_hbnd4 (hface4_GEO *, int, Gitter :: hbndseg_STI :: bnd_t) ;
+  virtual inline tetra_GEO     * insert_tetra (hface3_GEO *(&)[4], int (&)[4]) ;
+  virtual inline periodic3_GEO * insert_periodic3 (hface3_GEO *(&)[2], int (&)[2]) ;
+// Anfang - Neu am 23.5.02 (BS)
+  virtual inline periodic4_GEO * insert_periodic4 (hface4_GEO *(&)[2], int (&)[2]) ;
+// Ende - Neu am 23.5.02 (BS)
+        virtual inline hexa_GEO      * insert_hexa (hface4_GEO *(&)[6], int (&)[6]) ;
+      public :
+        inline MacroGitterBasis (istream &) ;
+        inline MacroGitterBasis () ;
+        virtual ~MacroGitterBasis () {}
+
+        // return index manager mostly for restore and backup 
+        inline IndexManagerType & indexManager(int codim); 
+        
+      protected: 
+        // index provider, for every codim one , 4 is for boundary
+        IndexManagerType _indexmanager[ numOfIndexManager ];
+    } ;
+} ;
+
+class GitterBasisImpl : public GitterBasis {
+  MacroGitterBasis * _macrogitter ;
+  public:
+   //us fuer Globalmethode levelwalk
+  inline Makrogitter & container () ;
+  inline const Makrogitter & container () const ;
+  public :
+    inline IndexManagerType & indexManager(int codim);
+        
+    inline GitterBasisImpl () ;
+    inline GitterBasisImpl (istream &) ;
+    inline GitterBasisImpl (const char *) ;
+    inline ~GitterBasisImpl () ;
+} ;
+
+
+  //
+  //    #    #    #  #          #    #    #  ######
+  //    #    ##   #  #          #    ##   #  #
+  //    #    # #  #  #          #    # #  #  #####
+  //    #    #  # #  #          #    #  # #  #
+  //    #    #   ##  #          #    #   ##  #
+  //    #    #    #  ######     #    #    #  ######
+  //
+
+inline GitterBasis :: Objects :: VertexEmpty :: VertexEmpty (int l, double x, double y, double z, IndexManagerType & im)
+  : GitterBasis :: VertexGeo (l,x,y,z,im) {
+  return ;
+}
+
+inline GitterBasis :: Objects :: VertexEmpty :: VertexEmpty (int l, double x, double y, double z, VertexGeo & vx )
+  : GitterBasis :: VertexGeo (l,x,y,z,vx) {
+  return ;
+}
+
+inline int GitterBasis :: Objects :: VertexEmpty :: ident () const {
+  cerr << "**FEHLER (FATAL) vertex :: ident () nur f\"ur level-0 Vertices zul\"assig " << endl ;
+  return (abort (), -1) ;
+}
+
+inline GitterBasis :: Objects :: VertexEmptyMacro :: VertexEmptyMacro (double x,double y,double z,int i, IndexManagerType &im) 
+  : GitterBasis :: Objects :: VertexEmpty (0,x,y,z,im), _idn (i) {
+  return ;
+}
+
+inline int GitterBasis :: Objects :: VertexEmptyMacro :: ident () const {
+  return _idn ;
+}
+
+inline GitterBasis :: Objects :: Hedge1Empty :: Hedge1Empty (myvertex_t * a, myvertex_t * b) 
+  : Gitter :: Geometric :: hedge1_GEO (a,b) {
+  return ;
+}
+
+inline void GitterBasis :: Objects :: Hedge1Empty :: projectInnerVertex(const ProjectVertex &pv) {
+  if (innerVertex()) {
+    assert(!leaf());
+    innerVertex()->project(pv);
+  }
+}
+
+inline GitterBasis :: Objects :: Hface3Empty :: Hface3Empty (myhedge1_t *e0, int s0, 
+  myhedge1_t *e1, int s1, myhedge1_t *e2, int s2) : Gitter :: Geometric :: hface3_GEO (e0, s0, e1, s1, e2, s2) {
+  return ;
+}
+
+inline void GitterBasis :: Objects :: Hface3Empty :: projectVertex(const ProjectVertex &pv) {
+  assert(!leaf());
+  for (int e = 0; e < polygonlength; e++)
+    myhedge1(e)->projectInnerVertex(pv);
+  if (innerVertex())
+    innerVertex()->project(pv);
+}
+
+inline GitterBasis :: Objects :: Hface4Empty :: Hface4Empty (myhedge1_t *e0, int s0, 
+  myhedge1_t *e1, int s1, myhedge1_t *e2, int s2, myhedge1_t *e3, int s3)
+  : Gitter :: Geometric :: hface4_GEO (e0, s0, e1, s1, e2, s2, e3, s3) {
+  return ;
+}
+
+inline void GitterBasis :: Objects :: Hface4Empty :: projectVertex(const ProjectVertex &pv) {
+  for (int e = 0; e < polygonlength; e++)
+    myhedge1(e)->projectInnerVertex(pv);
+  if (innerVertex())
+    innerVertex()->project(pv);
+}
+
+
+// Dune_hbndDefault 
+template <int pdimvx> 
+inline GitterBasis :: Objects :: Dune_hbndDefault<pdimvx> :: Dune_hbndDefault () : _index (-1) 
+{
+  for(int i=0; i<dimvx; i++)
+    for(int j=0; j<3; j++)
+      _oppVx[i][j] = 0.0;
+}
+
+
+template <int pdimvx> 
+inline void GitterBasis :: Objects :: Dune_hbndDefault<pdimvx> :: setOppPoint (int i, const double (&p)[3]) 
+{ 
+  assert((i >= 0) && (i < dimvx));
+  _oppVx[i][0] = p[0];
+  _oppVx[i][1] = p[1];
+  _oppVx[i][2] = p[2];
+  return;
+}
+           
+template <int pdimvx> 
+inline const double (& GitterBasis :: Objects :: Dune_hbndDefault<pdimvx> ::oppositeVertex (int i) const) [3] 
+{
+  assert((i >= 0) && (i < dimvx));
+  return _oppVx[i];
+}
+
+template <int pdimvx> inline int 
+GitterBasis :: Objects :: Dune_hbndDefault<pdimvx> :: getIndex () const 
+{ 
+  assert( _index >= 0 );
+  return _index;
+}
+
+template <int pdimvx> inline void  
+GitterBasis :: Objects :: Dune_hbndDefault<pdimvx> :: setIndex ( int idx ) 
+{
+  _index = idx;
+  return ;
+}
+
+template <int pdimvx> inline int 
+GitterBasis :: Objects :: Dune_hbndDefault<pdimvx> :: dimVx () const { return dimvx; }
+
+// end of Dune_hbndDefault 
+
+
+inline GitterBasis :: Objects :: Hbnd3Default :: Hbnd3Default (myhface3_t * f, int i, ProjectVertex *ppv) : Gitter :: Geometric :: hbndseg3_GEO (f, i, ppv) 
+{
+  return ;
+}
+
+inline GitterBasis :: Objects ::Hbnd3Default :: bnd_t GitterBasis :: Objects :: Hbnd3Default :: bndtype () const {
+  return undefined ;
+}
+
+inline int GitterBasis :: Objects :: Hbnd3Default :: ghostLevel () const {
+  assert(false);
+  return level() ;
+}
+
+inline void GitterBasis :: Objects :: Hbnd3Default :: faceNormal( double * normal) const {
+
+  hface3_GEO * face = this->myhface3(0);
+  int tw = this->twist(0);
+  BSGridLinearSurfaceMapping 
+    LSM(face->myvertex( (tw < 0) ? 0 : 2 )->Point(),
+        face->myvertex( 1                )->Point(),
+        face->myvertex( (tw < 0) ? 2 : 0 )->Point()
+       );
+  LSM.normal(normal);
+  return ;
+}
+
+inline GitterBasis :: Objects :: Hbnd4Default :: Hbnd4Default (myhface4_t * f, int i, ProjectVertex *ppv) : Gitter :: Geometric :: hbndseg4_GEO (f, i,ppv) 
+{
+  return ;
+}
+
+inline GitterBasis :: Objects ::Hbnd4Default :: bnd_t GitterBasis :: Objects :: Hbnd4Default :: bndtype () const {
+  return undefined ;
+}
+
+inline int GitterBasis :: Objects :: Hbnd4Default :: ghostLevel () const {
+  return level() ;
+}
+
+inline void GitterBasis :: Objects :: Hbnd4Default :: faceNormal( double * normal) const {
+  cerr << "ERORR: Hbnd4Default :: faceNormal not implemented! " << __FILE__ << __LINE__ << endl;
+  return ;
+}
+
+inline GitterBasis :: Objects :: TetraEmpty :: TetraEmpty (myhface3_t * f0, int t0, myhface3_t * f1, int t1,
+  myhface3_t * f2, int t2, myhface3_t * f3, int t3) : Gitter :: Geometric :: Tetra (f0, t0, f1, t1, f2, t2, f3, t3) {
+  return ;
+}
+
+inline GitterBasis :: Objects :: Periodic3Empty :: Periodic3Empty (myhface3_t * f0, int t0, myhface3_t * f1, int t1) 
+  : Gitter :: Geometric :: Periodic3 (f0, t0, f1, t1) {
+  return ;
+}
+
+// Anfang - Neu am 23.5.02 (BS)
+inline GitterBasis :: Objects :: Periodic4Empty :: Periodic4Empty (myhface4_t * f0, int t0, myhface4_t * f1, int t1) 
+  : Gitter :: Geometric :: Periodic4 (f0, t0, f1, t1) {
+  return ;
+}
+// Ende - Neu am 23.5.02 (BS)
+
+inline GitterBasis :: Objects :: HexaEmpty :: HexaEmpty (myhface4_t * f0, int t0, myhface4_t * f1, int t1,
+  myhface4_t * f2, int t2, myhface4_t * f3, int t3, myhface4_t * f4, int t4, myhface4_t * f5, int t5)
+  : Gitter :: Geometric :: hexa_GEO (f0, t0, f1, t1, f2, t2, f3, t3, f4, t4, f5, t5) {
+  return ;
+}
+
+inline GitterBasisImpl :: GitterBasisImpl () : _macrogitter (0) {
+  _macrogitter = new MacroGitterBasis () ;
+  assert (_macrogitter) ;
+  notifyMacroGridChanges () ;
+  return ;
+}
+
+inline GitterBasisImpl :: GitterBasisImpl (istream & in) : _macrogitter (0) {
+  _macrogitter = new MacroGitterBasis (in) ;
+  assert (_macrogitter) ;
+  notifyMacroGridChanges () ;
+  return ;
+}
+
+inline GitterBasisImpl :: GitterBasisImpl (const char * file) : _macrogitter (0) {
+  ifstream in (file) ;
+  if (!in) {
+    cerr << "  GitterBasisImpl :: GitterBasisImpl (const char *) FEHLER (IGNORIERT) " ;
+    cerr << "beim \"Offnen der Datei " << (file ? file : "\"null\"" ) << endl ;
+    _macrogitter = new MacroGitterBasis () ;
+  } else {
+    _macrogitter = new MacroGitterBasis (in) ;
+  }
+  assert (_macrogitter) ;
+  notifyMacroGridChanges () ;
+  return ;
+}
+
+inline GitterBasisImpl :: ~GitterBasisImpl () {
+  delete _macrogitter ;
+  return ;
+}
+
+inline Gitter :: Makrogitter & GitterBasisImpl :: container () {
+  return * _macrogitter ;
+}
+
+inline const Gitter :: Makrogitter & GitterBasisImpl :: container () const {
+  return * _macrogitter ;
+}
+
+inline IndexManagerType & GitterBasisImpl :: indexManager (int codim) 
+{ 
+  return _macrogitter->indexManager(codim);
+}
+
+inline GitterBasis :: MacroGitterBasis :: MacroGitterBasis (istream & in) {
+  macrogridBuilder (in) ;
+  return ;
+}
+
+inline GitterBasis :: MacroGitterBasis :: MacroGitterBasis () {
+  return ;
+}
+
+inline GitterBasis :: VertexGeo * GitterBasis :: MacroGitterBasis :: insert_vertex (double x, double y, double z, int id,int) {
+  return new Objects :: VertexEmptyMacro (x, y, z, id, indexManager(3)) ;
+}
+
+inline GitterBasis :: hedge1_GEO * GitterBasis :: MacroGitterBasis :: insert_hedge1 (VertexGeo * a, GitterBasis :: VertexGeo * b) {
+  return new Objects :: hedge1_IMPL (0, a, b, indexManager(2) ) ;
+}
+
+inline GitterBasis :: hface3_GEO * GitterBasis :: MacroGitterBasis :: insert_hface3 (hedge1_GEO *(&e)[3], int (&s)[3]) {
+  return new Objects :: hface3_IMPL (0,e[0],s[0],e[1],s[1],e[2],s[2], indexManager(1) ) ;
+}
+
+inline GitterBasis :: hface4_GEO * GitterBasis :: MacroGitterBasis :: insert_hface4 (hedge1_GEO *(&e)[4], int (&s)[4]) {
+  return new Objects :: hface4_IMPL (0, e[0],s[0],e[1],s[1],e[2],s[2],e[3],s[3], indexManager(1) ) ;
+}
+
+inline GitterBasis :: tetra_GEO * GitterBasis :: MacroGitterBasis :: insert_tetra (hface3_GEO *(&f)[4], int (&t)[4]) {
+  return new Objects :: tetra_IMPL (0,f[0],t[0],f[1],t[1],f[2],t[2],f[3],t[3], indexManager(0) ) ;
+}
+
+inline GitterBasis :: periodic3_GEO * GitterBasis :: MacroGitterBasis :: insert_periodic3 (hface3_GEO *(&f)[2], int (&t)[2]) {
+  return new Objects :: periodic3_IMPL (0,f[0],t[0],f[1],t[1]) ;
+}
+
+// Anfang - Neu am 23.5.02 (BS)
+inline GitterBasis :: periodic4_GEO * GitterBasis :: MacroGitterBasis :: insert_periodic4 (hface4_GEO *(&f)[2], int (&t)[2]) {
+  return new Objects :: periodic4_IMPL (0,f[0],t[0],f[1],t[1]) ;
+}
+// Ende - Neu am 23.5.02 (BS)
+
+inline GitterBasis :: hexa_GEO * GitterBasis :: MacroGitterBasis :: insert_hexa (hface4_GEO *(&f)[6], int (&t)[6]) {
+  return new Objects :: hexa_IMPL (0,f[0],t[0],f[1],t[1],f[2],t[2],f[3],t[3],f[4],t[4],f[5],t[5], indexManager(0) ) ;
+}
+
+inline GitterBasis :: hbndseg3_GEO * GitterBasis :: MacroGitterBasis :: 
+insert_hbnd3 (hface3_GEO * f, int i, Gitter :: hbndseg_STI :: bnd_t b) {
+  return new Objects :: hbndseg3_IMPL (0,f,i,NULL,NULL ,b, indexManager(4) , 0 ) ;
+}
+
+inline GitterBasis :: hbndseg3_GEO * GitterBasis :: MacroGitterBasis :: 
+insert_hbnd3 (hface3_GEO * f, int i, Gitter :: hbndseg_STI :: bnd_t b, const double (&p)[3]) {
+  return insert_hbnd3(f,i,b); 
+}
+
+inline GitterBasis :: hbndseg4_GEO * GitterBasis :: MacroGitterBasis :: insert_hbnd4 (hface4_GEO * f, int i, Gitter :: hbndseg_STI :: bnd_t b) {
+  return new Objects :: hbndseg4_IMPL (0,f,i,NULL, b,indexManager(4));
+}
+
+inline IndexManagerType & GitterBasis :: MacroGitterBasis :: indexManager (int codim ) 
+{ 
+  assert((codim >= 0) && (codim < numOfIndexManager ));
+  return _indexmanager[codim];
+}
+
+#endif  //  GITTER_IMPL_H_INCLUDED

src/serial/gitter_mgb.h

+// (c) bernhard schupp 1997 - 1998
+
+#ifndef GITTER_MGB_H_INCLUDED
+#define GITTER_MGB_H_INCLUDED
+
+#ifdef IBM_XLC
+  #define _ANSI_HEADER
+#endif
+
+#include <assert.h>
+
+#ifdef _ANSI_HEADER
+  using namespace std;
+  #include <vector>
+  #include <functional>
+  #include <utility>
+  #include <map>
+  #include <algorithm>
+#else
+  #include <vector.h>
+  #include <function.h>
+  #include <pair.h>
+  #include <map.h>
+  #include <algo.h>
+#endif
+
+#include "key.h"
+#include "gitter_sti.h"
+
+static volatile char RCSId_gitter_mgb_h [] = "$Id$" ;
+
+template < class RandomAccessIterator > inline int cyclicReorder (RandomAccessIterator begin, RandomAccessIterator end) {
+  RandomAccessIterator middle = min_element (begin,end) ;
+  int pos = middle == begin ? 0 : (rotate (begin,middle,end), (end - middle)) ;
+  if (*(begin + 1) < *(end - 1)) return pos ;
+  else {
+    reverse (begin,end) ;
+    rotate (begin,end - 1,end) ;
+    return - pos - 1 ;
+  }
+}
+
+class MacroGridBuilder : protected Gitter :: Geometric {
+  protected :
+  class Hbnd3IntStorage
+  {
+    double _p[3]; 
+    hface3_GEO * _first;
+    int          _second;
+    bool _pInit; // true if p was initialized with a value 
+  public:  
+    // store point and face and twist  
+    Hbnd3IntStorage( hface3_GEO * f, int tw, const double (&p) [3] );
+    
+    // store face and twist and set point to default 
+    Hbnd3IntStorage( hface3_GEO * f, int tw ); 
+
+    // return reference to _p
+    const double (& getPoint () const )[3];
+
+    hface3_GEO * first  () const { return _first;  }
+    int          second () const { return _second; }
+  };
+
+  protected :
+    enum ElementRawID {TETRA_RAW=4, HEXA_RAW=8, PERIODIC3_RAW=33, PERIODIC4_RAW=44} ; 
+  protected :
+    typedef long    vertexKey_t ;
+    typedef pair < int, int >   edgeKey_t ;
+    typedef Key3 < int >  faceKey_t ;
+    typedef Key4 < int >  elementKey_t ;
+
+    typedef map < vertexKey_t , VertexGeo *, less < vertexKey_t > >     vertexMap_t ;
+    typedef map < edgeKey_t,    hedge1_GEO *, less < edgeKey_t > >      edgeMap_t ;
+    typedef map < faceKey_t,    void *, less < faceKey_t > >            faceMap_t ;
+    typedef map < faceKey_t,    Hbnd3IntStorage *, less < faceKey_t > > hbndintMap_t ;
+    typedef map < elementKey_t, void *, less < elementKey_t > >         elementMap_t ;
+  
+    vertexMap_t  _vertexMap ;
+    edgeMap_t    _edgeMap ;
+    
+    faceMap_t    _face4Map, _face3Map, _hbnd3Map, _hbnd4Map, _hbnd4Int ;
+    hbndintMap_t _hbnd3Int; // new type here, so we dont have to cast to void *
+    // todo here: same thing for hbnd4int 
+    
+    elementMap_t _hexaMap, _tetraMap, _periodic3Map, _periodic4Map ;
+    
+    inline BuilderIF & myBuilder () ;
+    inline const BuilderIF & myBuilder () const ;
+    void removeElement (const elementKey_t &) ;
+  public :
+    virtual pair < VertexGeo *, bool >     InsertUniqueVertex (double, double, double, int) ;
+    virtual pair < hedge1_GEO *, bool >    InsertUniqueHedge1 (int,int) ;
+    virtual pair < hface3_GEO *, bool >    InsertUniqueHface3 (int (&)[3]) ;
+    virtual pair < hface4_GEO *, bool >    InsertUniqueHface4 (int (&)[4]) ;
+    virtual pair < tetra_GEO *, bool >     InsertUniqueTetra (int (&)[4]) ;
+    virtual pair < periodic3_GEO *, bool > InsertUniquePeriodic3 (int (&)[6]) ;
+
+    virtual pair < periodic4_GEO *, bool > InsertUniquePeriodic4 (int (&)[8]) ;
+    virtual pair < hexa_GEO *, bool >      InsertUniqueHexa (int (&)[8]) ;
+    
+    virtual bool InsertUniqueHbnd3 (int (&)[3], Gitter :: hbndseg :: bnd_t) ;
+    virtual bool InsertUniqueHbnd4 (int (&)[4], Gitter :: hbndseg :: bnd_t) ;
+
+  public :
+    static bool debugOption (int) ;
+    static void generateRawHexaImage (istream &, ostream &) ;
+    static void generateRawTetraImage (istream &, ostream &) ;
+    static void cubeHexaGrid (int, ostream &) ;
+    MacroGridBuilder (BuilderIF &, bool init = true) ;
+    virtual ~MacroGridBuilder () ;
+    void inflateMacroGrid (istream &) ;
+    void backupMacroGrid (ostream &) ;
+
+    // former constructor 
+    void initialize ();
+    // former destructor 
+    void finalize ();
+  protected:  
+    bool _initialized;
+    bool _finalized;
+  private :
+    BuilderIF & _mgb ;
+} ;
+
+
+//
+//    #    #    #  #          #    #    #  ######
+//    #    ##   #  #          #    ##   #  #
+//    #    # #  #  #          #    # #  #  #####
+//    #    #  # #  #          #    #  # #  #
+//    #    #   ##  #          #    #   ##  #
+//    #    #    #  ######     #    #    #  ######
+//
+inline Gitter :: Geometric :: BuilderIF & MacroGridBuilder :: myBuilder () {
+  return _mgb ;
+}
+
+inline const Gitter :: Geometric :: BuilderIF & MacroGridBuilder :: myBuilder () const {
+  return _mgb ;
+}
+
+inline bool MacroGridBuilder :: debugOption (int level) {
+  return (getenv ("VERBOSE_MGB") ? ( atoi (getenv ("VERBOSE_MGB")) > level ? true : (level == 0)) : false) ;
+}
+
+inline MacroGridBuilder :: Hbnd3IntStorage :: 
+Hbnd3IntStorage( hface3_GEO * f, int tw, const double (&p) [3] )
+ : _first(f) , _second(tw) , _pInit(true)
+{
+  for(int i=0; i<3; i++) _p[i] = p[i];
+}
+    
+inline MacroGridBuilder :: Hbnd3IntStorage :: 
+Hbnd3IntStorage( hface3_GEO * f, int tw )
+ : _first(f) , _second(tw) , _pInit(false)
+{
+  for(int i=0; i<3; i++) _p[i] = -666.0;
+}
+
+inline const double (& MacroGridBuilder :: Hbnd3IntStorage :: getPoint () const )[3]
+{ 
+  assert(_pInit);
+  return _p;
+}
+
+#endif

src/serial/headers.h

+#ifndef __FAKE_HEADERS_INCLUDED__
+#define __FAKE_HEADERS_INCLUDED__
+
+#include "key.h"
+#include "serialize.h"
+#include "lock.h"
+#include "gitter_sti.h"
+
+#endif

src/serial/indexstack.h

+// (c) Robert Kloefkorn 2004 - 2005 
+#ifndef INDEXSTACK_H_INCLUDED
+#define INDEXSTACK_H_INCLUDED
+
+#ifdef IBM_XLC
+  #define _ANSI_HEADER
+#endif
+
+#include <assert.h>
+
+#ifdef _ANSI_HEADER
+  using namespace std;
+  #include <stack>
+#else
+  #include <stack.h>
+#endif
+
+#ifdef _DUNE_USES_BSGRID_
+
+template<class T, int length>
+class FiniteStack {
+public :
+  // Makes empty stack
+  FiniteStack () : _f(0) {}
+
+  // Returns true if the stack is empty
+  bool empty () const { return _f==0; }
+
+  // Returns true if the stack is full
+  bool full () const { return (_f >= length); }
+
+  // Puts a new object onto the stack
+  void push (const T& t) { _s[_f++] = t; }
+
+  // Removes and returns the uppermost object from the stack
+  T pop () { return _s[--_f]; }
+
+  // Returns the uppermost object on the stack
+  T top () const { return _s[_f-1]; }
+
+  // stacksize
+  int size () const { return _f; }
+
+private:
+   T   _s[length]; // the stack 
+   int _f;         // actual position in stack  
+};
+
+
+//******************************************************
+//
+//  IndexStack providing indices via getIndex and freeIndex
+//  indices that are freed, are put on a stack and get
+//
+//******************************************************
+template <class T, int length> 
+class IndexStack 
+{
+  typedef FiniteStack<T,length> StackType;
+  typedef stack < StackType * > StackListType;
+  
+  StackListType fullStackList_;
+  StackListType emptyStackList_;
+  
+  //typedef typename StackListType::Iterator DListIteratorType;
+  StackType * stack_; 
+
+  // current maxIndex 
+  int maxIndex_; 
+public:
+  //! Constructor, create new IndexStack
+  IndexStack(); 
+
+  //! Destructor, deleting all stacks 
+  inline ~IndexStack (); 
+
+  //! set index as maxIndex if index is bigger than maxIndex
+  void checkAndSetMax(T index) { if(index > maxIndex_) maxIndex_ = index;  }
+  
+  //! set index as maxIndex
+  void setMaxIndex(T index) { maxIndex_ = index; }
+
+  //! set index as maxIndex
+  int getMaxIndex() const { return maxIndex_;  }
+  
+  //! restore index from stack or create new index 
+  T getIndex (); 
+
+  //! store index on stack 
+  void freeIndex(T index);
+
+  //! test stack funtcionality
+  void test ();
+
+  // backup set to out stream 
+  void backupIndexSet ( ostream & os ); 
+
+  // restore from in stream 
+  void restoreIndexSet ( istream & is );
+private:
+  // no copy constructor allowed 
+  IndexStack( const IndexStack<T,length> & s) : maxIndex_ (0) , stack_(0) {}
+ 
+  // no assignment operator allowed 
+  IndexStack<T,length> & operator = ( const IndexStack<T,length> & s) 
+  {
+    cout << "IndexStack::operator = () not allowed! in: " __FILE__ << " line:" << __LINE__ << "\n";
+    abort();
+    return *this; 
+  }
+  
+  // clear all stored indices 
+  void clearStack ();
+};  // end class IndexStack 
+
+//****************************************************************
+// Inline implementation 
+// ***************************************************************
+template <class T, int length>
+inline IndexStack<T,length>::IndexStack()
+  : stack_ ( new StackType () ) , maxIndex_ (0) {} 
+  
+template <class T, int length>
+inline IndexStack<T,length>::~IndexStack () 
+{
+  if(stack_) 
+  {
+    delete stack_;
+    stack_ = new StackType();
+    assert(stack_);
+  }
+
+  while( !fullStackList_.empty() )
+  {
+    StackType * st = fullStackList_.top();
+    if(st) delete st; 
+    fullStackList_.pop();
+  }
+  while( !emptyStackList_.empty() )
+  {
+    StackType * st = emptyStackList_.top();
+    if(st) delete st; 
+    emptyStackList_.pop();
+  }
+}
+
+template <class T, int length>
+inline T IndexStack<T,length>::getIndex () 
+{
+  if((*stack_).empty()) 
+  {
+    if( fullStackList_.size() <= 0)
+    {
+      return maxIndex_++;
+    }
+    else 
+    {
+      emptyStackList_.push( stack_ );
+      stack_ = fullStackList_.top();
+      fullStackList_.pop();
+    }
+  }
+  return (*stack_).pop();
+}
+
+template <class T, int length>
+inline void IndexStack<T,length>::freeIndex ( T index ) 
+{
+  if((*stack_).full())
+  {
+    fullStackList_.push(  stack_ );
+    if(emptyStackList_.size() <= 0)
+    {
+      stack_ = new StackType (); 
+    }
+    else 
+    {
+      stack_ = emptyStackList_.top();
+      emptyStackList_.pop();
+    }
+  }
+  (*stack_).push(index); 
+}
+
+template <class T, int length>
+inline void IndexStack<T,length>::test () 
+{
+  T vec[2*length];
+
+  for(int i=0; i<2*length; i++)
+    vec[i] = getIndex();
+
+  for(int i=0; i<2*length; i++)
+    freeIndex(vec[i]);
+  
+  for(int i=0; i<2*length; i++)
+    vec[i] = getIndex();
+  
+  for(int i=0; i<2*length; i++)
+    printf(" index [%d] = %d \n",i,vec[i]);
+}
+
+template <class T, int length>
+inline void IndexStack<T,length>::backupIndexSet ( ostream & os ) 
+{
+  // holes are not stored at the moment 
+  os.write( ((const char *) &maxIndex_ ), sizeof(int) ) ;
+  return ;
+}
+
+template <class T, int length>
+inline void IndexStack<T,length>::restoreIndexSet ( istream & is )
+{
+  is.read ( ((char *) &maxIndex_), sizeof(int) );
+  clearStack ();
+
+  return ;
+}
+
+template <class T, int length>
+inline void IndexStack<T,length>::clearStack () 
+{
+  if(stack_) 
+  {
+    delete stack_;
+    stack_ = new StackType();
+    assert(stack_);
+  }
+
+  while( !fullStackList_.empty() )
+  {
+    StackType * st = fullStackList_.top();
+    if(st) delete st; 
+    fullStackList_.pop();
+  }
+  return;
+}
+
+
+#else 
+
+//*********************************************************************
+//
+// Dummy implementation for the index stack, if index is not used 
+//
+//*********************************************************************
+template <class T>
+struct DummyIndexStack 
+{
+  //! set index as maxIndex if index is bigger than maxIndex
+  inline void checkAndSetMax(T index)
+  {
+  }
+  
+  //! set index as maxIndex
+  inline void setMaxIndex(T index)
+  {
+  }
+  
+  //! restore index from stack or create new index 
+  T getIndex () 
+  {
+    return -1;
+  }
+
+  //! store index on stack 
+  inline void freeIndex(T index)
+  {
+  }
+
+  //! test stack funtcionality
+  inline void test ()
+  {
+  }
+};  // end class DummyIndexStack 
+
+#endif 
+
+#endif

src/serial/key.h

+// (c) bernhard schupp 1997 - 1998
+#ifndef KEY_H_INCLUDED
+#define KEY_H_INCLUDED
+
+template < class A > class Key3 {
+  protected :
+    A _a, _b, _c ;
+  public :
+    inline Key3 () ;
+    inline Key3 (const A &,const A &,const A &) ;
+    inline Key3 (const Key3 < A > &) ;
+    inline const Key3 < A > & operator = (const Key3 < A > &) ;
+    inline bool operator < (const Key3 < A > &) const ;
+} ;
+
+template < class A > class Key4 {
+  protected :
+    A _a, _b, _c, _d ;
+  public :
+    inline Key4 () ;
+    inline Key4 (const A &, const A &, const A &, const A &) ;
+    inline Key4 (const Key4 < A > &) ;
+    inline const Key4 < A > & operator = (const Key4 < A > &) ;
+    inline bool operator < (const Key4 < A > &) const ;
+} ;
+
+template < class A > inline Key3 < A > :: Key3 () : _a (-1), _b (-1), _c (-1) {
+  return ;
+}
+
+template < class A > inline Key3 < A > :: Key3 (const A & a, const A & b, const A & c) : _a (a), _b (b), _c (c) {
+  return ;
+}
+
+template < class A > inline Key3 < A > :: Key3 (const Key3 < A > & k) : _a (k._a), _b (k._b), _c (k._c) {
+  return ;
+}
+
+template < class A > inline const Key3 < A > & Key3 < A > :: operator = (const Key3 < A > & k) {
+  _a = k._a ;
+  _b = k._b ;
+  _c = k._c ;
+  return * this ;
+}
+
+template < class A > inline bool Key3 < A > :: operator < (const Key3 < A > & k) const {
+  return _a < k._a ? true : (_a == k._a ? (_b < k._b ? true : (_b == k._b ? (_c < k._c ? true : false) : false)) : false) ;
+}
+
+template < class A > inline Key4 < A > :: Key4 () : _a (), _b (), _c (), _d () {
+  return ;
+}
+
+template < class A > inline Key4 < A > :: Key4 (const A & a, const A & b, const A & c, const A & d) : _a(a), _b (b), _c (c), _d (d) {
+  return ;
+}
+
+template < class A > inline Key4 < A > :: Key4 (const Key4 < A > & k) : _a (k._a), _b (k._b), _c(k._c), _d (k._d) {
+  return ;
+}
+
+template < class A > inline const Key4 < A > & Key4 < A > :: operator = (const Key4 < A > & k) {
+  _a = k._a ;
+  _b = k._b ;
+  _c = k._c ;
+  _d = k._d ;
+  return * this ;
+}
+
+template < class A > inline bool Key4 < A >  :: operator < (const Key4 < A > & k) const {
+  return _a < k._a ? true : (_a == k._a ? (_b < k._b ? true : (_b == k._b ? (_c < k._c ? true :
+  (_c == k._c ? (_d < k._d ? true : false) : false)) : false)) : false) ;
+}
+
+#endif

src/serial/lock.h

+// (c) bernhard schupp, 1997 - 1998
+#ifndef LOCK_H_INCLUDED
+#define LOCK_H_INCLUDED
+
+#ifdef IBM_XLC
+  #define _ANSI_HEADER
+#endif
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <assert.h>
+
+#ifdef _ANSI_HEADER
+  using namespace std;
+  #include <iostream>
+#else
+  #include <iostream.h>
+#endif
+
+	// Einfache Klasse, die w"ahrend ihrer Lebnsdauer ein
+	// Lockfile mit einem vorgegebenen Namen (Pfad) h"alt.
+
+class FSLock {
+  char * _fname ;
+  public :
+    FSLock (const char * = "") ;
+   ~FSLock () ;
+} ;
+
+inline FSLock :: FSLock (const char * name) : _fname (0) {
+  _fname = new char [strlen(name) + 100] ;
+  assert (_fname) ;
+  sprintf (_fname, "%s.lock", name) ;
+  FILE * fp = fopen (_fname, "w") ;
+  if (fp == NULL) {
+    delete [] _fname ;
+    _fname = 0 ;
+    cerr << "**WARNUNG (IGNORIERT) Lockfile konnte nicht erzeugt werden" << endl ;
+  } else {
+    int test = fclose (fp) ;
+    assert (test == 0) ;
+  }
+  return ;
+}
+
+inline FSLock :: ~FSLock () {
+  if (_fname) {
+    int test = remove (_fname) ;
+    if (test != 0) {
+      cerr << "**WARNUNG (IGNORIERT) Lockfile konnte nicht gel\"oscht werden." << endl ;
+    }
+    delete [] _fname ;
+    _fname = 0 ;
+  }
+  return ;
+}
+
+#endif  // LOCK_H_INCLUDED

src/serial/mapp_cube_3d.cc

+	// (c) bernhard schupp 1997 - 1998
+
+	// $Source$
+	// $Revision$
+	// $Name$
+	// $State$
+	
+/* $Id$
+ * $Log$
+ * Revision 1.1  2005/03/23 14:57:55  robertk
+ * all files for serial version of ALU3dGrid.
+ *
+ * Revision 1.2  2004/10/25 16:38:11  robertk
+ * All header end with .h now. Like the original.
+ *
+ * In the .cc this changes are done.
+ *
+ * Revision 1.1  2004/10/15 09:48:37  robertk
+ * Inititial version. Some extenxions for Dune made. Schould be compatible
+ * with all other applications done so far.
+ *
+ * Revision 1.2  2001/12/10 13:57:23  wesenber
+ * RCS Log history and/or RCSId-variable added
+ *
+ ***/ 
+
+#include <math.h>
+#include <assert.h>
+#include <stdlib.h>
+
+#include "mapp_cube_3d.h"
+
+static volatile char RCSId_mapp_cube_3d_cc [] = "$Id$" ;
+
+const double TrilinearMapping :: _epsilon = 1.0e-8 ;
+const double QuadraturCube3Dbasis :: _p2 [4][3] = { { .816496580927726, .0,  .5773502691896258},
+                                                    { .0, .816496580927726, -.5773502691896258},
+                                                    { -.816496580927726, .0, .5773502691896258},
+                                                    { .0, -.816496580927726, -.5773502691896258} } ;
+const double QuadraturCube3Dbasis :: _p3 [8][3] = { {  .5773502691896258,  .5773502691896258,  .5773502691896258},
+                                                    { -.5773502691896258,  .5773502691896258,  .5773502691896258},
+                                                    {  .5773502691896258, -.5773502691896258,  .5773502691896258},
+                                                    { -.5773502691896258, -.5773502691896258,  .5773502691896258},
+                                                    {  .5773502691896258,  .5773502691896258, -.5773502691896258},
+                                                    { -.5773502691896258,  .5773502691896258, -.5773502691896258},
+                                                    {  .5773502691896258, -.5773502691896258, -.5773502691896258},
+                                                    { -.5773502691896258, -.5773502691896258, -.5773502691896258} } ;
+
+const double QuadraturCube2Dbasis :: _p1 [2] = { .0, .0 } ;
+const double QuadraturCube2Dbasis :: _p3 [4][2] = { { .5773502691896258,  .5773502691896258 },
+                                                    {-.5773502691896258,  .5773502691896258 },
+                                                    { .5773502691896258, -.5773502691896258 },
+                                                    {-.5773502691896258, -.5773502691896258 } } ;
+
+
+void TrilinearMapping :: linear(const double (&p)[3]) {
+  double x = .5 * (p[0] + 1.) ;
+  double y = .5 * (p[1] + 1.) ;
+  double z = .5 * (p[2] + 1.) ;
+  double t0 = .5 ;
+  double t3 = y * z ;
+  double t8 = x * z ;
+  double t13 = x * y ;
+  Df[2][0] = t0 * ( a[1][2] + y * a[4][2] + z * a[6][2] + t3 * a[7][2] ) ;
+  Df[2][1] = t0 * ( a[2][2] + x * a[4][2] + z * a[5][2] + t8 * a[7][2] ) ;
+  Df[1][2] = t0 * ( a[3][1] + y * a[5][1] + x * a[6][1] + t13 * a[7][1] ) ;
+  Df[2][2] = t0 * ( a[3][2] + y * a[5][2] + x * a[6][2] + t13 * a[7][2] ) ;
+  Df[0][0] = t0 * ( a[1][0] + y * a[4][0] + z * a[6][0] + t3 * a[7][0] ) ;
+  Df[0][2] = t0 * ( a[3][0] + y * a[5][0] + x * a[6][0] + t13 * a[7][0] ) ;
+  Df[1][0] = t0 * ( a[1][1] + y * a[4][1] + z * a[6][1] + t3 * a[7][1] ) ;
+  Df[0][1] = t0 * ( a[2][0] + x * a[4][0] + z * a[5][0] + t8 * a[7][0] ) ;
+  Df[1][1] = t0 * ( a[2][1] + x * a[4][1] + z * a[5][1] + t8 * a[7][1] ) ;
+
+}
+
+double TrilinearMapping :: det(const double (&point)[3]) {
+	//  Determinante der Abbildung f:[-1,1]^3 -> Hexaeder im Punkt point.
+  linear (point) ;
+  return (DetDf = Df[0][0] * Df[1][1] * Df[2][2] - Df[0][0] * Df[1][2] * Df[2][1] - 
+	          Df[1][0] * Df[0][1] * Df[2][2] + Df[1][0] * Df[0][2] * Df[2][1] + 
+	          Df[2][0] * Df[0][1] * Df[1][2] - Df[2][0] * Df[0][2] * Df[1][1]) ;
+}
+
+void TrilinearMapping :: inverse(const double (&p)[3]) {
+	//  Kramer - Regel, det() rechnet Df und DetDf neu aus.
+  double val = 1.0 / det(p) ;
+  Dfi[0][0] = ( Df[1][1] * Df[2][2] - Df[1][2] * Df[2][1] ) * val ;
+  Dfi[0][1] = ( Df[0][2] * Df[2][1] - Df[0][1] * Df[2][2] ) * val ;
+  Dfi[0][2] = ( Df[0][1] * Df[1][2] - Df[0][2] * Df[1][1] ) * val ;
+  Dfi[1][0] = ( Df[1][2] * Df[2][0] - Df[1][0] * Df[2][2] ) * val ;
+  Dfi[1][1] = ( Df[0][0] * Df[2][2] - Df[0][2] * Df[2][0] ) * val ;
+  Dfi[1][2] = ( Df[0][2] * Df[1][0] - Df[0][0] * Df[1][2] ) * val ;
+  Dfi[2][0] = ( Df[1][0] * Df[2][1] - Df[1][1] * Df[2][0] ) * val ;
+  Dfi[2][1] = ( Df[0][1] * Df[2][0] - Df[0][0] * Df[2][1] ) * val ;
+  Dfi[2][2] = ( Df[0][0] * Df[1][1] - Df[0][1] * Df[1][0] ) * val ;
+  return ;
+}
+
+void TrilinearMapping :: world2map (const double (&wld)[3], double (&map)[3]) {
+	//  Newton - Iteration zum Invertieren der Abbildung f.
+  double err = 10.0 * _epsilon ;
+#ifndef NDEBUG
+  int count = 0 ;
+#endif
+  map [0] = map [1] = map [2] = .0 ;
+  do {
+    double upd [3] ;
+    map2world (map, upd) ;
+    inverse (map) ;
+    double u0 = upd [0] - wld [0] ;
+    double u1 = upd [1] - wld [1] ;
+    double u2 = upd [2] - wld [2] ;
+    double c0 = Dfi [0][0] * u0 + Dfi [0][1] * u1 + Dfi [0][2] * u2 ;
+    double c1 = Dfi [1][0] * u0 + Dfi [1][1] * u1 + Dfi [1][2] * u2 ;
+    double c2 = Dfi [2][0] * u0 + Dfi [2][1] * u1 + Dfi [2][2] * u2 ;
+    map [0] -= c0 ;
+    map [1] -= c1 ;
+    map [2] -= c2 ;
+    err = fabs (c0) + fabs (c1) + fabs (c2) ;
+    assert (count ++ < 1000) ;
+  } while (err > _epsilon) ;
+  return ;
+}

src/serial/mapp_tetra_3d_ext.h

+#ifndef __MAPP_TETRA_3D_EXT_HH__
+#define __MAPP_TETRA_3D_EXT_HH__
+
+// BSGridVecType is defined in BSGrid interface to dune or in gitter_sti.hh 
+class BSGridLinearSurfaceMapping : public LinearSurfaceMapping 
+{
+public: 
+  inline BSGridLinearSurfaceMapping 
+    (const double (&)[3], const double (&)[3], const double (&)[3]);
+    
+  // same as method normal of LinearSurfaceMapping, just for Dune Vecs 
+  inline void normal(double * normal) const;
+};
+
+
+inline BSGridLinearSurfaceMapping :: BSGridLinearSurfaceMapping (const double (&x0)[3], 
+const double (&x1)[3], const double (&x2)[3])
+  : LinearSurfaceMapping (x0,x1,x2)
+{
+}
+
+inline void BSGridLinearSurfaceMapping :: normal (double * normal) const
+{
+  normal[0] = this->_n[0];
+  normal[1] = this->_n[1];
+  normal[2] = this->_n[2];
+  return ;
+}
+#endif