Skip to content
Snippets Groups Projects
Commit 7929e086 authored by Adrian Burri's avatar Adrian Burri
Browse files

Added support for hexahedra - works serially, not yet tested with adaptation 

[[Imported from SVN: r2038]]
parent ad6b7d73
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 1358 additions and 605 deletions
grid/alu3dgrid.hh
+ 295
− 42
View file @ 7929e086
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef __DUNE_ALU3DGRID_HH__
#define __DUNE_ALU3DGRID_HH__
#ifndef DUNE_ALU3DGRID_HH
#define DUNE_ALU3DGRID_HH
#include <vector>
#include <assert.h>
#include <cassert>
#include <dune/common/misc.hh>
#include <dune/common/fmatrix.hh>
......@@ -13,6 +13,7 @@
#include "common/defaultindexsets.hh"
#include "alu3dgrid/alu3dinclude.hh"
#include "alu3dgrid/alu3dmappings.hh"
#include <dune/common/exceptions.hh>
#include <dune/common/stdstreams.hh>
......@@ -27,7 +28,71 @@ namespace Dune
//#undef DUNE_THROW
//#define DUNE_THROW(e,m) assert(false);
enum ALU3dGridElementType { tetra = 4, hexa = 7 };
enum ALU3dGridElementType { tetra = 4, hexa = 7, mixed, error };
//ALU3dGridElementType convertGeometryType2ALU3dGridElementType(GeometryType);
//GeometryType convertALU3dGridElementType2GeometryType(ALU3dGridElementType);
template <ALU3dGridElementType elType>
struct ALU3dImplTraits;
template <>
struct ALU3dImplTraits<tetra> {
typedef ALU3DSPACE GEOFace3Type GEOFaceType;
typedef ALU3DSPACE GEOEdgeT GEOEdgeType;
typedef ALU3DSPACE GEOVertexT GEOVertexType;
typedef ALU3DSPACE IMPLTetraElementType IMPLElementType;
typedef ALU3DSPACE GEOTetraElementType GEOElementType;
typedef ALU3DSPACE HasFace3Type HasFaceType;
typedef ALU3DSPACE BNDFace3Type BNDFaceType;
typedef ALU3DSPACE ImplBndFace3Type ImplBndFaceType;
typedef ALU3DSPACE BNDFace3Type PLLBndFaceType;
// refinement and coarsening enum for tetrahedons
enum { refine_element_t =
ALU3DSPACE GitterType::Geometric::TetraRule::iso8 };
enum { coarse_element_t =
ALU3DSPACE GitterType::Geometric::TetraRule::crs };
typedef std::pair<GEOFaceType*, int> NeighbourFaceType;
typedef std::pair<HasFaceType*, int> NeighbourPairType;
typedef std::pair<PLLBndFaceType*, int> GhostPairType;
static int alu2duneFace(int index) { return index; }
static int dune2aluFace(int index) { return index; }
};
template <>
struct ALU3dImplTraits<hexa> {
typedef ALU3DSPACE GEOFace4Type GEOFaceType;
typedef ALU3DSPACE GEOEdgeT GEOEdgeType;
typedef ALU3DSPACE GEOVertexT GEOVertexType;
typedef ALU3DSPACE IMPLHexaElementType IMPLElementType;
typedef ALU3DSPACE GEOHexaElementType GEOElementType;
typedef ALU3DSPACE HasFace4Type HasFaceType;
typedef ALU3DSPACE BNDFace4Type BNDFaceType;
typedef ALU3DSPACE ImplBndFace4Type ImplBndFaceType;
typedef ALU3DSPACE BNDFace4Type PLLBndFaceType;
// refinement and coarsening enum for hexahedrons
enum { refine_element_t = ALU3DSPACE GitterType::Geometric::HexaRule::iso8 };
enum { coarse_element_t = ALU3DSPACE GitterType::Geometric::HexaRule::crs };
typedef std::pair<GEOFaceType*, int> NeighbourFaceType;
typedef std::pair<HasFaceType*, int> NeighbourPairType;
typedef std::pair<PLLBndFaceType*, int> GhostPairType;
//inline static int alu2duneVertex(int);
//inline static int dune2aluVertex(int);
static int alu2duneFace(int index) { return alu2duneFace_[index]; }
static int dune2aluFace(int index) { return dune2aluFace_[index]; }
//inline static int alu2duneQuad(int);
//inline static int dune2aluQuad(int);
private:
static const int alu2duneFace_[6];
static const int dune2aluFace_[6];
};
// i.e. double or float
typedef double alu3d_ctype;
......@@ -41,10 +106,12 @@ namespace Dune
template<class GridImp> class ALU3dGridHierarchicIterator;
template<class GridImp> class ALU3dGridIntersectionIterator;
template<class GridImp> class ALU3dGridLeafIterator;
template<int dim, int dimworld> class ALU3dGrid;
template<int dim, int dimworld,
ALU3dGridElementType elType> class ALU3dGrid;
// the hierarhic index set
template<int dim, int dimworld> class ALU3dGridHierarchicIndexSet;
template<int dim, int dimworld,
ALU3dGridElementType elType> class ALU3dGridHierarchicIndexSet;
// singleton holding reference elements
template<int dim,class GridImp> struct ALU3dGridReferenceGeometry;
......@@ -75,6 +142,9 @@ namespace Dune
GridImp, ALU3dGridGeometry>
{
typedef Geometry<mydim, coorddim, GridImp, ALU3dGridGeometry> GeometryType;
typedef ALU3dImplTraits<GridImp::elementType>::PLLBndFaceType PLLBndFaceType;
friend class ALU3dGridIntersectionIterator<GridImp>;
public:
ALU3dGridMakeableGeometry(bool makeRefelem=false) :
GeometryType (ALU3dGridGeometry<mydim, coorddim,GridImp>(makeRefelem)) {}
......@@ -90,8 +160,13 @@ namespace Dune
return this->realGeometry.buildGeom(item);
}
bool buildGeom(const ALU3DSPACE HFaceType& item,
int twist, int faceIdx) {
return this->realGeometry.buildGeom(item, twist, faceIdx);
}
// call buildGhost of realGeometry
bool buildGhost(const ALU3DSPACE PLLBndFaceType & ghost)
bool buildGhost(const PLLBndFaceType & ghost)
{
return this->realGeometry.buildGhost(ghost);
}
......@@ -111,12 +186,39 @@ namespace Dune
};
//! ALU3dGridGeometry
template<int mydim, int cdim, class GridImp>
class ALU3dGridGeometry :
public GeometryDefault <mydim,cdim,GridImp,ALU3dGridGeometry>
// calculates m^p at compile-time
template <int m, int p>
struct POWER_M_P
{
// power stores m^p
enum { power = (m * POWER_M_P<m,p-1>::power ) };
};
// end of recursion via specialization
template <int m>
struct POWER_M_P< m , 0>
{
// m^0 = 1
enum { power = 1 };
};
//! Empty definition, needs to be specialized for element type
template <int mydim, int cdim, class GridImp>
class ALU3dGridGeometry :
public GeometryDefault <mydim,cdim,GridImp,ALU3dGridGeometry> {};
template <int mydim, int cdim>
class ALU3dGridGeometry<mydim, cdim, const ALU3dGrid<3, 3, tetra> > :
public GeometryDefault<mydim, cdim, const ALU3dGrid<3, 3, tetra>,
ALU3dGridGeometry> {
typedef const ALU3dGrid<3, 3, tetra> GridImp;
friend class ALU3dGridBoundaryEntity<GridImp>;
typedef ALU3dImplTraits<tetra>::IMPLElementType IMPLElementType;
typedef ALU3dImplTraits<tetra>::PLLBndFaceType PLLBndFaceType;
typedef ALU3dImplTraits<tetra>::GEOFaceType GEOFaceType;
typedef ALU3dImplTraits<tetra>::GEOEdgeType GEOEdgeType;
typedef ALU3dImplTraits<tetra>::GEOVertexType GEOVertexType;
//! know dimension of barycentric coordinates
enum { dimbary=mydim+1};
public:
......@@ -160,13 +262,13 @@ namespace Dune
//! Methods that not belong to the Interface, but have to be public
//***********************************************************************
//! generate the geometry for out of given ALU3dGridElement
bool buildGeom(const ALU3DSPACE IMPLElementType & item);
bool buildGeom(const ALU3DSPACE HFaceType & item);
bool buildGeom(const IMPLElementType & item);
bool buildGeom(const ALU3DSPACE HFaceType & item, int twist, int faceIdx);
bool buildGeom(const ALU3DSPACE HEdgeType & item);
bool buildGeom(const ALU3DSPACE VertexType & item);
//! build ghost out of internal boundary segment
bool buildGhost(const ALU3DSPACE PLLBndFaceType & ghost);
bool buildGhost(const PLLBndFaceType & ghost);
//! print internal data
//! no interface method
......@@ -183,7 +285,7 @@ namespace Dune
void calcElMatrix () const;
//! the vertex coordinates
mutable FieldMatrix<alu3d_ctype,mydim+1,cdim> coord_;
mutable FieldMatrix<alu3d_ctype, POWER_M_P<2,mydim>::power, cdim> coord_;
//! is true if Jinv_, A and detDF_ is calced
mutable bool builtinverse_;
......@@ -202,6 +304,105 @@ namespace Dune
mutable FieldVector<alu3d_ctype,cdim> tmpU_; //! temporary memory
};
template <int mydim, int cdim>
class ALU3dGridGeometry<mydim, cdim, const ALU3dGrid<3, 3, hexa> > :
public GeometryDefault<mydim, cdim, const ALU3dGrid<3, 3, hexa>,
ALU3dGridGeometry> {
typedef const ALU3dGrid<3, 3, hexa> GridImp;
friend class ALU3dGridBoundaryEntity<GridImp>;
friend class ALU3dGridIntersectionIterator<GridImp>;
typedef ALU3dImplTraits<hexa>::IMPLElementType IMPLElementType;
typedef ALU3dImplTraits<hexa>::PLLBndFaceType PLLBndFaceType;
typedef ALU3dImplTraits<hexa>::GEOFaceType GEOFaceType;
typedef ALU3dImplTraits<hexa>::GEOEdgeType GEOEdgeType;
typedef ALU3dImplTraits<hexa>::GEOVertexType GEOVertexType;
public:
//! for makeRefGeometry == true a Geometry with the coordinates of the
//! reference element is made
ALU3dGridGeometry(bool makeRefGeometry=false);
//! Destructor
~ALU3dGridGeometry();
//! return the element type identifier
//! line , triangle or tetrahedron, depends on dim
GeometryType type () const;
//! return the number of corners of this element. Corners are numbered 0..n-1
int corners () const;
//! access to coordinates of corners. Index is the number of the corner
const FieldVector<alu3d_ctype, cdim>& operator[] (int i) const;
//! return reference element corresponding to this element. If this is
//! a reference element then self is returned.
static const Dune::Geometry<mydim,mydim,GridImp,Dune::ALU3dGridGeometry> & refelem ();
//! maps a local coordinate within reference element to
//! global coordinate in element
FieldVector<alu3d_ctype, cdim> global (const FieldVector<alu3d_ctype, mydim>& local) const;
//! maps a global coordinate within the element to a
//! local coordinate in its reference element
FieldVector<alu3d_ctype, mydim> local (const FieldVector<alu3d_ctype, cdim>& global) const;
//! returns true if the point in local coordinates is inside reference
//! element
bool checkInside(const FieldVector<alu3d_ctype, mydim>& local) const;
//! A(l) , see grid.hh
alu3d_ctype integrationElement (const FieldVector<alu3d_ctype, mydim>& local) const;
//! can only be called for dim=dimworld! (Trivially true, since there is no
//! other specialization...)
const FieldMatrix<alu3d_ctype,mydim,mydim>& jacobianInverse (const FieldVector<alu3d_ctype, cdim>& local) const;
//***********************************************************************
//! Methods that not belong to the Interface, but have to be public
//***********************************************************************
//! generate the geometry out of a given ALU3dGridElement
bool buildGeom(const IMPLElementType & item);
bool buildGeom(const ALU3DSPACE HFaceType & item, int twist, int faceIdx);
bool buildGeom(const ALU3DSPACE HEdgeType & item);
bool buildGeom(const ALU3DSPACE VertexType & item);
//! build ghost out of internal boundary segment
bool buildGhost(const PLLBndFaceType & ghost);
//! print internal data
//! no interface method
void print (std::ostream& ss) const;
// for changing the coordinates of one element
FieldVector<alu3d_ctype, cdim> & getCoordVec (int i);
private:
//! calculates local index of ALU3dGrid Face using the twist of the face regarding the element prototype
int faceTwist(int val, int idx) const;
//! the vertex coordinates
mutable FieldMatrix<alu3d_ctype, POWER_M_P<2,mydim>::power, cdim> coord_;
//mutable FieldVector<alu3d_ctype, mydim> tmp1_;
mutable FieldVector<alu3d_ctype, cdim> tmp2_;
TrilinearMapping* triMap_;
BilinearSurfaceMapping* biMap_;
mutable FieldMatrix<alu3d_ctype, 3, 3> jInv_;
const static int alu2duneVol[8];
const static int dune2aluVol[8];
const static int alu2duneFace[6];
const static int dune2aluFace[6];
const static int alu2duneQuad[4];
const static int dune2aluQuad[4];
const static int dune2aluFaceVertex[6][4];
};
//**********************************************************************
//
......@@ -215,6 +416,10 @@ namespace Dune
{
// typedef typename GridImp::template codim<codim>::Entity EntityType;
friend class ALU3dGridEntity<codim, dim, GridImp>;
typedef ALU3dImplTraits<GridImp::elementType>::PLLBndFaceType PLLBndFaceType;
typedef ALU3dImplTraits<GridImp::elementType>::IMPLElementType IMPLElementType;
public:
// Constructor creating the realEntity
......@@ -237,7 +442,7 @@ namespace Dune
}
//! set original element pointer to fake entity
void setGhost(ALU3DSPACE PLLBndFaceType &ghost)
void setGhost(PLLBndFaceType &ghost)
{
this->realEntity.setGhost(ghost);
}
......@@ -275,11 +480,11 @@ namespace Dune
{
enum { dimworld = GridImp::dimensionworld };
friend class ALU3dGrid < dim , dimworld >;
friend class ALU3dGrid < dim , dimworld, GridImp::elementType >;
friend class ALU3dGridEntity < 0, dim, GridImp >;
friend class ALU3dGridLevelIterator < cd, All_Partition, GridImp >;
friend class ALU3dGridHierarchicIndexSet<dim,dimworld>;
friend class ALU3dGridHierarchicIndexSet<dim,dimworld,GridImp::elementType>;
public:
typedef typename ALU3DSPACE ALUHElementType<cd>::ElementType BSElementType;
......@@ -376,8 +581,15 @@ namespace Dune
: public EntityDefault<0,dim,GridImp,ALU3dGridEntity>
{
enum { dimworld = GridImp::dimensionworld };
typedef ALU3dImplTraits<GridImp::elementType>::GEOElementType GEOElementType;
typedef ALU3dImplTraits<GridImp::elementType>::PLLBndFaceType PLLBndFaceType;
friend class ALU3dGrid < dim , dimworld >;
enum { refine_element_t =
ALU3dImplTraits<GridImp::elementType>::refine_element_t };
enum { coarse_element_t =
ALU3dImplTraits<GridImp::elementType>::coarse_element_t };
friend class ALU3dGrid < dim , dimworld, GridImp::elementType>;
friend class ALU3dGridIntersectionIterator < GridImp >;
friend class ALU3dGridHierarchicIterator < const GridImp >;
friend class ALU3dGridHierarchicIterator < GridImp >;
......@@ -387,7 +599,7 @@ namespace Dune
friend class ALU3dGridLevelIterator <3,All_Partition,GridImp>;
friend class ALU3dGridLeafIterator <GridImp>;
friend class ALU3dGridHierarchicIndexSet<dim,dimworld>;
friend class ALU3dGridHierarchicIndexSet<dim,dimworld,GridImp::elementType>;
public:
typedef typename GridImp::template codim<0>::Geometry Geometry;
......@@ -496,7 +708,7 @@ namespace Dune
void setGhost(ALU3DSPACE HElementType &ghost);
//! set original element pointer to fake entity
void setGhost(ALU3DSPACE PLLBndFaceType &ghost);
void setGhost(PLLBndFaceType &ghost);
//! set actual walk level
void reset ( int l );
......@@ -509,6 +721,7 @@ namespace Dune
void setEntity ( const ALU3dGridEntity<0,dim,GridImp> & org );
private:
typedef ALU3dImplTraits<GridImp::elementType>::IMPLElementType IMPLElementType;
//! index is unique within the grid hierachie and per codim
int getIndex () const;
......@@ -517,10 +730,10 @@ namespace Dune
const GridImp & grid_;
// the current element of grid
ALU3DSPACE IMPLElementType *item_;
IMPLElementType *item_;
// the current ghost, if element is ghost
ALU3DSPACE PLLBndFaceType * ghost_;
PLLBndFaceType * ghost_;
mutable bool isGhost_; //! true if entity is ghost entity
//! the cuurent geometry
......@@ -735,6 +948,12 @@ namespace Dune
enum { dim = GridImp::dimension };
enum { dimworld = GridImp::dimensionworld };
typedef ALU3dImplTraits<GridImp::elementType>::GEOElementType GEOElementType;
typedef ALU3dImplTraits<GridImp::elementType>::GEOFaceType GEOFaceType;
typedef ALU3dImplTraits<GridImp::elementType>::NeighbourPairType NeighbourPairType;
typedef ALU3dImplTraits<GridImp::elementType>::PLLBndFaceType PLLBndFaceType;
typedef ALU3dImplTraits<GridImp::elementType>::BNDFaceType BNDFaceType;
friend class ALU3dGridEntity<0,dim,GridImp>;
public:
......@@ -749,7 +968,8 @@ namespace Dune
//! The default Constructor , level tells on which level we want
//! neighbours
ALU3dGridIntersectionIterator(const GridImp & grid, ALU3DSPACE HElementType *el,
ALU3dGridIntersectionIterator(const GridImp & grid,
ALU3DSPACE HElementType *el,
int wLevel,bool end=false);
//! The copy constructor
......@@ -812,9 +1032,36 @@ namespace Dune
NormalType integrationOuterNormal (const FieldVector<alu3d_ctype, dim-1>& local) const;
private:
typedef ALU3dImplTraits<GridImp::elementType>::NeighbourFaceType NeighbourFaceType;
// calculate normal
void calculateNormal(const FieldVector<alu3d_ctype, dim-1>& local,
NormalType& result) const;
// calculate normal seen from neighbor
void calculateNormalNeighbor(const FieldVector<alu3d_ctype, dim-1>& local,
NormalType& result) const;
// get the face corresponding to the index
ALU3dImplTraits<tetra>::GEOFaceType& getFace(int, Int2Type<tetra>) const;
ALU3dImplTraits<hexa>::GEOFaceType& getFace(int, Int2Type<hexa>) const;
// if neighbour exists , do setup of new neighbour
void setNeighbor () const ;
// \brief get neighboring face \index of the element where the iteration started
// Index conversion from Dune to Alu3dGrid reference element is done in this
// method
// \param index Local face index in the Dune reference element
NeighbourPairType getNeighPair (int index) const;
// Index conversion from Dune to Alu3dGrid reference element is done in this
// method
// \param index Local face index in the Dune reference element
NeighbourFaceType getNeighFace (int index) const;
// check wether we at ghost boundary, only parallel grid
void checkGhost () const ;
......@@ -829,36 +1076,38 @@ namespace Dune
//! the grid
//const GridImp & grid_;
const int nFaces_;
int walkLevel_;
//EntityImp * entity_; //! neighbour entity
// current element from which we started the intersection iterator
mutable ALU3DSPACE GEOElementType *item_;
mutable GEOElementType *item_;
//! current neighbour
mutable ALU3DSPACE GEOElementType *neigh_;
mutable GEOElementType *neigh_;
//! current ghost if we have one
mutable ALU3DSPACE PLLBndFaceType *ghost_;
mutable PLLBndFaceType *ghost_;
mutable int index_; //! internal index of intersection
mutable int numberInNeigh_; //! index of intersection in neighbor
mutable bool theSituation_; //! true if the "situation" occurs :-)
mutable bool daOtherSituation_; //! true if the "da other situation" occurs :-)
//! see bsgrid.cc for descritption
//! see alugrid.cc for description
mutable bool isBoundary_; //! true if intersection is with boundary
mutable bool isGhost_; //! true if intersection is with internal boundary (only parallel grid)
mutable bool needSetup_; //! true if setup is needed
mutable bool needNormal_; //! true if normal has to be calculated
// pair holding pointer to face and twist
mutable ALU3DSPACE NeighbourFaceType neighpair_;
mutable NeighbourFaceType neighpair_;
mutable bool initInterGl_; //! true if interSelfGlobal_ was initialized
mutable bool twist_; //! true if orientation is such that normal points into hexahedron
GeometryImp * interSelfGlobal_; //! intersection_self_global
MakeableBndEntityImp * bndEntity_; //! boundaryEntity
......@@ -883,7 +1132,7 @@ namespace Dune
friend class ALU3dGridEntity<2,dim,GridImp>;
friend class ALU3dGridEntity<1,dim,GridImp>;
friend class ALU3dGridEntity<0,dim,GridImp>;
friend class ALU3dGrid < dim , dimworld >;
friend class ALU3dGrid < dim , dimworld, GridImp::elementType >;
public:
typedef typename GridImp::template codim<cd>::Entity Entity;
......@@ -991,14 +1240,14 @@ namespace Dune
@}
*/
template <int dim, int dimworld>
class ALU3dGrid : public GridDefault < dim, dimworld, alu3d_ctype,ALU3dGrid<dim,dimworld> >
template <int dim, int dimworld, ALU3dGridElementType elType>
class ALU3dGrid : public GridDefault < dim, dimworld, alu3d_ctype,ALU3dGrid<dim,dimworld, elType> >
{
//CompileTimeChecker<dim == 3> ALU3dGrid_only_implemented_for_3dp;
//CompileTimeChecker<dimworld == 3> ALU3dGrid_only_implemented_for_3dw;
//CompileTimeChecker< (eltype == ALU3DSPACE tetra_t) || (eltype == ALU3DSPACE hexa_t ) > ALU3dGrid_only_implemented_for_tetra_or_hexa;
typedef ALU3dGrid<dim,dimworld> MyType;
typedef ALU3dGrid<dim,dimworld,elType> MyType;
friend class ALU3dGridEntity <0,dim,MyType>;
friend class ALU3dGridEntity <0,dim,const MyType>;
friend class ALU3dGridIntersectionIterator<MyType>;
......@@ -1015,7 +1264,8 @@ namespace Dune
// The Interface Methods
//**********************************************************
public:
enum { myElementType = tetra };
static const ALU3dGridElementType elementType = elType;
typedef GridTraits<dim,dimworld, MyType ,
ALU3dGridGeometry,ALU3dGridEntity,
ALU3dGridBoundaryEntity,
......@@ -1037,7 +1287,7 @@ namespace Dune
//typedef typename std::pair < ObjectStreamType * , ALU3dGridEntity<0,dim,dimworld> * >
// DataCollectorParamType;
typedef ALU3dGridHierarchicIndexSet<dim,dimworld> HierarchicIndexSetType;
typedef ALU3dGridHierarchicIndexSet<dim,dimworld,elType> HierarchicIndexSetType;
typedef DefaultLevelIndexSet<MyType> LevelIndexSetType;
typedef typename Traits::LeafIterator LeafIterator;
......@@ -1114,6 +1364,8 @@ namespace Dune
const HierarchicIndexSetType & hierarchicIndexSet () const { return hIndexSet_; }
const LevelIndexSetType & levelIndexSet () const
{
// * This is pure evil when adapting
assert(false);
if(!levelIndexSet_) levelIndexSet_ = new LevelIndexSetType (*this);
return *levelIndexSet_;
}
......@@ -1179,7 +1431,7 @@ namespace Dune
bool mark( int refCount , const typename Traits::template codim<0>::Entity & en );
template <int cd>
ALU3dGridEntity<cd,dim,const ALU3dGrid<dim,dimworld> >&
ALU3dGridEntity<cd,dim,const MyType >&
getRealEntity(typename Traits::template codim<cd>::Entity& entity)
{
return entity.realEntity;
......@@ -1187,7 +1439,7 @@ namespace Dune
//private:
template <int cd>
const ALU3dGridEntity<cd,dim,const ALU3dGrid<dim,dimworld> >&
const ALU3dGridEntity<cd,dim,const MyType >&
getRealEntity(const typename Traits::template codim<cd>::Entity& entity) const
{
return entity.realEntity;
......@@ -1195,10 +1447,11 @@ namespace Dune
private:
//! Copy constructor should not be used
ALU3dGrid( const ALU3dGrid<dim,dimworld> & g);
ALU3dGrid( const MyType & g);
//! assignment operator should not be used
ALU3dGrid<dim,dimworld> & operator = (const ALU3dGrid<dim,dimworld> & g);
ALU3dGrid<dim,dimworld,elType>&
operator = (const MyType & g);
// reset size and global size
void calcExtras();
......@@ -1256,13 +1509,13 @@ namespace Dune
mutable BndProvider bndProvider_;
//mutable VertexProvider vertexProvider_;
}; // end Class ALU3dGridGrid
}; // end class ALU3dGridGrid
//! hierarchic index set of ALU3dGrid
template <int dim, int dimworld>
template <int dim, int dimworld, ALU3dGridElementType elType>
class ALU3dGridHierarchicIndexSet
{
typedef ALU3dGrid<dim,dimworld> GridType;
typedef ALU3dGrid<dim,dimworld,elType> GridType;
enum { numCodim = 4 };
public:
......
grid/alu3dgrid/Makefile.am
+ 6
− 1
View file @ 7929e086
AM_CPPFLAGS = $(ALU3DGRID_CPPFLAGS)
noinst_LTLIBRARIES = libalu3d.la
libalu3d_la_SOURCES = alu3dmappings.cc
alu3dgriddir = $(includedir)/dune/grid/alu3dgrid/
alu3dgrid_HEADERS = alu3dgrid.cc alu3dinclude.hh alu3dgridcomm.hh \
datahandle.hh leafwalk.hh alumemory.hh myautoptr.hh alu3dgeometry.cc \
alu3dgeometry.hh alu3dmappings.cc alu3dmappings.hh
include $(top_srcdir)/am/sourcescheck
grid/alu3dgrid/alu3dgeometry.cc
+ 1
− 215
View file @ 7929e086
......@@ -3,221 +3,6 @@
#ifndef DUNE_ALU3DGEOMETRY_CC
#define DUNE_ALU3DGEOMETRY_CC
const double TrilinearMapping :: _epsilon = 1.0e-8 ;
inline TrilinearMapping ::
TrilinearMapping (const coord_t& x0, const coord_t& x1,
const coord_t& x2, const coord_t& x3,
const coord_t& x4, const coord_t& x5,
const coord_t& x6, const coord_t& x7)
: p0(x0), p1(x1), p2(x2), p3(x3), p4(x4), p5(x5), p6(x6), p7(x7) {
a [0][0] = p0 [0] ;
a [0][1] = p0 [1] ;
a [0][2] = p0 [2] ;
a [1][0] = p1 [0] - p0 [0] ;
a [1][1] = p1 [1] - p0 [1] ;
a [1][2] = p1 [2] - p0 [2] ;
a [2][0] = p2 [0] - p0 [0] ;
a [2][1] = p2 [1] - p0 [1] ;
a [2][2] = p2 [2] - p0 [2] ;
a [3][0] = p4 [0] - p0 [0] ;
a [3][1] = p4 [1] - p0 [1] ;
a [3][2] = p4 [2] - p0 [2] ;
a [4][0] = p3 [0] - p2 [0] - a [1][0] ;
a [4][1] = p3 [1] - p2 [1] - a [1][1] ;
a [4][2] = p3 [2] - p2 [2] - a [1][2] ;
a [5][0] = p6 [0] - p4 [0] - a [2][0] ;
a [5][1] = p6 [1] - p4 [1] - a [2][1] ;
a [5][2] = p6 [2] - p4 [2] - a [2][2] ;
a [6][0] = p5 [0] - p1 [0] - a [3][0] ;
a [6][1] = p5 [1] - p1 [1] - a [3][1] ;
a [6][2] = p5 [2] - p1 [2] - a [3][2] ;
a [7][0] = p7 [0] - p5 [0] + p4 [0] - p6 [0] - p3 [0] + p1 [0] + a [2][0] ;
a [7][1] = p7 [1] - p5 [1] + p4 [1] - p6 [1] - p3 [1] + p1 [1] + a [2][1] ;
a [7][2] = p7 [2] - p5 [2] + p4 [2] - p6 [2] - p3 [2] + p1 [2] + a [2][2] ;
return ;
}
inline TrilinearMapping :: TrilinearMapping (const TrilinearMapping & map)
: p0(map.p0), p1(map.p1), p2(map.p2), p3(map.p3),
p4(map.p4), p5(map.p5), p6(map.p6), p7(map.p7) {
for (int i = 0 ; i < 8 ; i ++)
for (int j = 0 ; j < 3 ; j ++)
a [i][j] = map.a [i][j] ;
return ;
}
FieldMatrix<alu3d_ctype, 3, 3>
TrilinearMapping::jacobianInverse(const coord_t& p) {
inverse (p);
return Dfi;
}
inline void TrilinearMapping ::
map2world(const coord_t& p, coord_t& world) const {
alu3d_ctype x = p [0];
alu3d_ctype y = p [1];
alu3d_ctype z = p [2];
alu3d_ctype t3 = y * z ;
alu3d_ctype t8 = x * z ;
alu3d_ctype t13 = x * y ;
alu3d_ctype t123 = x * t3 ;
world [0] = a [0][0] + a [1][0] * x + a [2][0] * y + a [3][0] * z + a [4][0] * t13 + a [5][0] * t3 + a [6][0] * t8 + a [7][0] * t123 ;
world [1] = a [0][1] + a [1][1] * x + a [2][1] * y + a [3][1] * z + a [4][1] * t13 + a [5][1] * t3 + a [6][1] * t8 + a [7][1] * t123 ;
world [2] = a [0][2] + a [1][2] * x + a [2][2] * y + a [3][2] * z + a [4][2] * t13 + a [5][2] * t3 + a [6][2] * t8 + a [7][2] * t123 ;
return ;
}
inline void TrilinearMapping ::
map2world(const alu3d_ctype x1, const alu3d_ctype x2,
const alu3d_ctype x3, coord_t& world ) const {
coord_t map ;
map [0] = x1 ;
map [1] = x2 ;
map [2] = x3 ;
map2world (map, world) ;
return ;
}
void TrilinearMapping :: linear(const coord_t& p ) {
alu3d_ctype x = p[0];
alu3d_ctype y = p[1];
alu3d_ctype z = p[2];
// ? get rid of that alu3d_ctype t0 = .5 ;
alu3d_ctype t3 = y * z ;
alu3d_ctype t8 = x * z ;
alu3d_ctype t13 = x * y ;
Df[2][0] = a[1][2] + y * a[4][2] + z * a[6][2] + t3 * a[7][2] ;
Df[2][1] = a[2][2] + x * a[4][2] + z * a[5][2] + t8 * a[7][2] ;
Df[1][2] = a[3][1] + y * a[5][1] + x * a[6][1] + t13 * a[7][1] ;
Df[2][2] = a[3][2] + y * a[5][2] + x * a[6][2] + t13 * a[7][2] ;
Df[0][0] = a[1][0] + y * a[4][0] + z * a[6][0] + t3 * a[7][0] ;
Df[0][2] = a[3][0] + y * a[5][0] + x * a[6][0] + t13 * a[7][0] ;
Df[1][0] = a[1][1] + y * a[4][1] + z * a[6][1] + t3 * a[7][1] ;
Df[0][1] = a[2][0] + x * a[4][0] + z * a[5][0] + t8 * a[7][0] ;
Df[1][1] = a[2][1] + x * a[4][1] + z * a[5][1] + t8 * a[7][1] ;
}
alu3d_ctype TrilinearMapping :: det(const coord_t& point ) {
// Determinante der Abbildung f:[-1,1]^3 -> Hexaeder im Punkt point.
linear (point) ;
return (DetDf = Df.determinant());
}
void TrilinearMapping :: inverse(const coord_t& p ) {
// Kramer - Regel, det() rechnet Df und DetDf neu aus.
alu3d_ctype 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 coord_t& wld , coord_t& map ) {
// Newton - Iteration zum Invertieren der Abbildung f.
alu3d_ctype err = 10.0 * _epsilon ;
#ifndef NDEBUG
int count = 0 ;
#endif
map [0] = map [1] = map [2] = .0 ;
do {
coord_t upd ;
map2world (map, upd) ;
inverse (map) ;
alu3d_ctype u0 = upd [0] - wld [0] ;
alu3d_ctype u1 = upd [1] - wld [1] ;
alu3d_ctype u2 = upd [2] - wld [2] ;
alu3d_ctype c0 = Dfi [0][0] * u0 + Dfi [0][1] * u1 + Dfi [0][2] * u2 ;
alu3d_ctype c1 = Dfi [1][0] * u0 + Dfi [1][1] * u1 + Dfi [1][2] * u2 ;
alu3d_ctype 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 ;
}
//- Bilinear surface mapping
inline BilinearSurfaceMapping ::
BilinearSurfaceMapping (const coord3_t& x0, const coord3_t& x1,
const coord3_t& x2, const coord3_t& x3)
: _p0 (x0), _p1 (x1), _p2 (x2), _p3 (x3) {
_b [0][0] = _p0 [0] ;
_b [0][1] = _p0 [1] ;
_b [0][2] = _p0 [2] ;
_b [1][0] = _p1 [0] - _p0 [0] ;
_b [1][1] = _p1 [1] - _p0 [1] ;
_b [1][2] = _p1 [2] - _p0 [2] ;
_b [2][0] = _p2 [0] - _p0 [0] ;
_b [2][1] = _p2 [1] - _p0 [1] ;
_b [2][2] = _p2 [2] - _p0 [2] ;
_b [3][0] = _p2 [0] - _p3 [0] - _b [1][0] ;
_b [3][1] = _p2 [1] - _p3 [1] - _b [1][1] ;
_b [3][2] = _p2 [2] - _p3 [2] - _b [1][2] ;
_n [0][0] = _b [1][1] * _b [2][2] - _b [1][2] * _b [2][1] ;
_n [0][1] = _b [1][2] * _b [2][0] - _b [1][0] * _b [2][2] ;
_n [0][2] = _b [1][0] * _b [2][1] - _b [1][1] * _b [2][0] ;
_n [1][0] = _b [1][1] * _b [3][2] - _b [1][2] * _b [3][1] ;
_n [1][1] = _b [1][2] * _b [3][0] - _b [1][0] * _b [3][2] ;
_n [1][2] = _b [1][0] * _b [3][1] - _b [1][1] * _b [3][0] ;
_n [2][0] = _b [3][1] * _b [2][2] - _b [3][2] * _b [2][1] ;
_n [2][1] = _b [3][2] * _b [2][0] - _b [3][0] * _b [2][2] ;
_n [2][2] = _b [3][0] * _b [2][1] - _b [3][1] * _b [2][0] ;
return ;
}
inline BilinearSurfaceMapping ::
BilinearSurfaceMapping (const BilinearSurfaceMapping & m)
: _p0(m._p0), _p1(m._p1), _p2(m._p2), _p3(m._p3) {
{for (int i = 0 ; i < 4 ; i ++)
for (int j = 0 ; j < 3 ; j ++ )
_b [i][j] = m._b [i][j] ;}
{for (int i = 0 ; i < 3 ; i ++)
for (int j = 0 ; j < 3 ; j ++ )
_n [i][j] = m._n [i][j] ;}
return ;
}
inline void BilinearSurfaceMapping ::
map2world (const coord2_t& map, coord3_t& wld) const {
double x = map [0];
double y = map [1];
double xy = x * y ;
wld[0] = _b [0][0] + x * _b [1][0] + y * _b [2][0] + xy * _b [3][0] ;
wld[1] = _b [0][1] + x * _b [1][1] + y * _b [2][1] + xy * _b [3][1] ;
wld[2] = _b [0][2] + x * _b [1][2] + y * _b [2][2] + xy * _b [3][2] ;
return ;
}
inline void BilinearSurfaceMapping ::
map2world (double x, double y, coord3_t& w) const {
coord2_t p ;
p [0] = x ;
p [1] = y ;
map2world (p,w) ;
return ;
}
inline void BilinearSurfaceMapping ::
normal (const coord2_t& map, coord3_t& normal) const {
double x = map [0];
double y = map [1];
normal [0] = -( _n [0][0] + _n [1][0] * x + _n [2][0] * y) ;
normal [1] = -( _n [0][1] + _n [1][1] * x + _n [2][1] * y) ;
normal [2] = -( _n [0][2] + _n [1][2] * x + _n [2][2] * y) ;
return ;
}
//- Tetra specialization
template<int mydim, int cdim>
inline ALU3dGridGeometry<mydim,cdim,const ALU3dGrid<3, 3, tetra> >::
......@@ -739,6 +524,7 @@ template<>
inline alu3d_ctype
ALU3dGridGeometry<2, 3, const ALU3dGrid<3, 3, hexa> >::
integrationElement (const FieldVector<alu3d_ctype, 2>& local) const {
// * is this right
biMap_->normal(local, tmp2_);
return tmp2_.two_norm();
}
......
grid/alu3dgrid/alu3dgeometry.hh
+ 1
− 57
View file @ 7929e086
......@@ -3,6 +3,7 @@
#ifndef DUNE_ALU3DGEOMETRY_HH
#define DUNE_ALU3DGEOMETRY_HH
#include "alu3dmappings.hh"
// * temporary, should be moved to misc.hh or something
// calculates m^p at compile-time
template <int m, int p>
......@@ -20,63 +21,6 @@ struct POWER_M_P< m , 0>
enum { power = 1 };
};
//! A trilinear mapping from the Dune reference hexahedron into the physical
//! space (same as in mapp_cube_3d.h, but for a different reference hexahedron)
class TrilinearMapping {
typedef FieldVector<alu3d_ctype, 3> coord_t;
typedef FieldMatrix<alu3d_ctype, 3, 3> mat_t;
static const alu3d_ctype _epsilon ;
const coord_t& p0;
const coord_t& p1;
const coord_t& p2;
const coord_t& p3;
const coord_t& p4;
const coord_t& p5;
const coord_t& p6;
const coord_t& p7;
double a [8][3] ;
mat_t Df;
mat_t Dfi;
alu3d_ctype DetDf ;
void linear (const coord_t&) ;
void inverse (const coord_t&) ;
public:
inline TrilinearMapping (const coord_t&, const coord_t&,
const coord_t&, const coord_t&,
const coord_t&, const coord_t&,
const coord_t&, const coord_t&);
inline TrilinearMapping (const TrilinearMapping &) ;
~TrilinearMapping () {}
alu3d_ctype det (const coord_t&) ;
mat_t jacobianInverse(const coord_t&);
inline void map2world (const coord_t&, coord_t&) const ;
inline void map2world (const double , const double , const double ,
coord_t&) const ;
void world2map (const coord_t&, coord_t&) ;
};
//! A bilinear surface mapping
class BilinearSurfaceMapping {
typedef FieldVector<alu3d_ctype, 3> coord3_t;
typedef FieldVector<alu3d_ctype, 2> coord2_t;
const coord3_t& _p0;
const coord3_t& _p1;
const coord3_t& _p2;
const coord3_t& _p3;
double _b [4][3] ;
double _n [3][3] ;
public:
inline BilinearSurfaceMapping (const coord3_t&, const coord3_t&,
const coord3_t&, const coord3_t&) ;
inline BilinearSurfaceMapping (const BilinearSurfaceMapping &) ;
~BilinearSurfaceMapping () {}
inline void map2world(const coord2_t&, coord3_t&) const ;
inline void map2world(double x, double y, coord3_t&) const ;
inline void normal(const coord2_t&, coord3_t&) const ;
} ;
//! Empty definition, needs to be specialized for element type
template <int mydim, int cdim, class GridImp>
class ALU3dGridGeometry :
......
grid/alu3dgrid/alu3dgrid.cc
+ 939
− 236
View file @ 7929e086
This diff is collapsed.
grid/alu3dgrid/alu3dinclude.hh
+ 27
− 28
View file @ 7929e086
......@@ -39,12 +39,12 @@ namespace ALU3dGridSpace {
typedef GitterDunePll GitterType;
typedef GitterDunePll GitterImplType;
typedef Hbnd3PllInternal < GitterType :: Objects :: Hbnd3Default,
BndsegPllBaseXClosure < GitterType :: hbndseg3_GEO > ,
BndsegPllBaseXMacroClosure < GitterType :: hbndseg3_GEO > > :: micro_t MicroType;
typedef Hbnd3PllInternal<GitterType::Objects::Hbnd3Default,
BndsegPllBaseXClosure<GitterType::hbndseg3_GEO>,
BndsegPllBaseXMacroClosure<GitterType::hbndseg3_GEO> > :: micro_t MicroType;
// value for boundary to other processes
static const int ProcessorBoundary_t = GitterImplType:: hbndseg_STI :: closure;
static const int ProcessorBoundary_t = GitterImplType::hbndseg_STI::closure;
#else
typedef GatherScatter GatherScatterType;
......@@ -52,7 +52,7 @@ namespace ALU3dGridSpace {
// the header
typedef Gitter GitterType;
typedef GitterDuneImpl GitterImplType;
typedef GitterType::hface_STI PLLFaceType; // Interface Face
typedef GitterType::hface_STI PLLFaceType; // Interface Face
#endif
......@@ -60,28 +60,27 @@ namespace ALU3dGridSpace {
#include "leafwalk.hh"
// typedefs of Element types
typedef GitterType::helement_STI HElementType; // Interface Element
typedef GitterType::hface_STI HFaceType; // Interface Face
typedef GitterType::hedge_STI HEdgeType; // Interface Edge
typedef GitterType::vertex_STI VertexType; // Interface Vertex
typedef GitterType::Geometric::hface3_GEO GEOFaceType; // real Face
typedef GitterType::Geometric::hedge1_GEO GEOEdgeType; // real Face
typedef GitterType::Geometric::VertexGeo GEOVertexType; // real Face
typedef GitterImplType::Objects::tetra_IMPL IMPLElementType; // impl Element
typedef GitterType::Geometric::tetra_GEO GEOElementType; // real Element
typedef GitterType::Geometric::hasFace3 HasFace3Type; // has Face with 3 polygons
typedef GitterImplType::Objects::Hbnd3Default BNDFaceType; // boundary segment
typedef GitterImplType::Objects::hbndseg3_IMPL ImplBndFaceType; // boundary segment
typedef BNDFaceType PLLBndFaceType;
// refinement and coarsening enum for tetrahedons
enum { refine_element_t = GitterType::Geometric::TetraRule::iso8 };
enum { coarse_element_t = GitterType::Geometric::TetraRule::crs };
typedef pair < GEOFaceType * , int > NeighbourFaceType;
typedef pair < HasFace3Type* , int > NeighbourPairType;
typedef pair < PLLBndFaceType * , int > GhostPairType;
typedef GitterType::helement_STI HElementType; // Interface Element
typedef GitterType::hface_STI HFaceType; // Interface Face
typedef GitterType::hedge_STI HEdgeType; // Interface Edge
typedef GitterType::vertex_STI VertexType; // Interface Vertex
typedef GitterType::Geometric::hface3_GEO GEOFace3Type; // Tetra Face
typedef GitterType::Geometric::hface4_GEO GEOFace4Type; // Hexa Face
typedef GitterType::Geometric::hedge1_GEO GEOEdgeT; // * stays real Face
typedef GitterType::Geometric::VertexGeo GEOVertexT; // * stays real Face
typedef GitterImplType::Objects::tetra_IMPL IMPLTetraElementType; //impl Element
typedef GitterImplType::Objects::hexa_IMPL IMPLHexaElementType;
typedef GitterType::Geometric::tetra_GEO GEOTetraElementType; // real Element
typedef GitterType::Geometric::hexa_GEO GEOHexaElementType;
typedef GitterType::Geometric::hasFace3 HasFace3Type; // has Face with 3 polygons
typedef GitterType::Geometric::hasFace4 HasFace4Type;
typedef GitterImplType::Objects::Hbnd3Default BNDFace3Type; // boundary segment
typedef GitterImplType::Objects::Hbnd4Default BNDFace4Type;
typedef GitterImplType::Objects::hbndseg3_IMPL ImplBndFace3Type; // boundary segment
typedef GitterImplType::Objects::hbndseg4_IMPL ImplBndFace4Type;
// * end new
//*************************************************************
// definition of original LeafIterators of ALU3dGrid
......@@ -90,7 +89,7 @@ namespace ALU3dGridSpace {
template <int codim>
struct BSMacroIterator
{
typedef AccessIterator < GitterType::helement_STI > :: Handle IteratorType;
typedef AccessIterator<GitterType::helement_STI>::Handle IteratorType;
};
//******************************************************************
......
grid/alu3dgrid/alu3dmappings.cc
+ 82
− 25
View file @ 7929e086
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#include "alu3dmappings.hh"
#include "../alu3dgrid.hh"
namespace Dune {
const int ALU3dImplTraits<hexa>::
alu2duneFace_[6] = {4, 5, 1, 3, 0, 2};
//{4, 5, 2, 1, 3, 0};
const int ALU3dImplTraits<hexa>::
dune2aluFace_[6] = {4, 2, 5, 3, 0, 1};
//{5, 3, 2, 4, 0, 1};
const double TrilinearMapping :: _epsilon = 1.0e-8 ;
void TrilinearMapping :: world2map (const coord_t& wld , coord_t& map ) {
// 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 {
coord_t upd ;
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 ;
}
//- Bilinear surface mapping
/*
BilinearSurfaceMapping ::
BilinearSurfaceMapping (const coord3_t& x0, const coord3_t& x1,
const coord3_t& x2, const coord3_t& x3)
: _p0 (x0), _p1 (x1), _p2 (x2), _p3 (x3) {
_b [0][0] = _p0 [0] ;
_b [0][1] = _p0 [1] ;
_b [0][2] = _p0 [2] ;
_b [1][0] = _p1 [0] - _p0 [0] ;
_b [1][1] = _p1 [1] - _p0 [1] ;
_b [1][2] = _p1 [2] - _p0 [2] ;
_b [2][0] = _p2 [0] - _p0 [0] ;
_b [2][1] = _p2 [1] - _p0 [1] ;
_b [2][2] = _p2 [2] - _p0 [2] ;
_b [3][0] = _p3 [0] - _p2 [0] - _b [1][0] ;
_b [3][1] = _p3 [1] - _p2 [1] - _b [1][1] ;
_b [3][2] = _p3 [2] - _p2 [2] - _b [1][2] ;
_n [0][0] = _b [1][1] * _b [2][2] - _b [1][2] * _b [2][1] ;
_n [0][1] = _b [1][2] * _b [2][0] - _b [1][0] * _b [2][2] ;
_n [0][2] = _b [1][0] * _b [2][1] - _b [1][1] * _b [2][0] ;
_n [1][0] = _b [1][1] * _b [3][2] - _b [1][2] * _b [3][1] ;
_n [1][1] = _b [1][2] * _b [3][0] - _b [1][0] * _b [3][2] ;
_n [1][2] = _b [1][0] * _b [3][1] - _b [1][1] * _b [3][0] ;
_n [2][0] = _b [3][1] * _b [2][2] - _b [3][2] * _b [2][1] ;
_n [2][1] = _b [3][2] * _b [2][0] - _b [3][0] * _b [2][2] ;
_n [2][2] = _b [3][0] * _b [2][1] - _b [3][1] * _b [2][0] ;
return ;
}
BilinearSurfaceMapping ::
BilinearSurfaceMapping (const BilinearSurfaceMapping & m)
: _p0(m._p0), _p1(m._p1), _p2(m._p2), _p3(m._p3) {
{for (int i = 0 ; i < 4 ; i ++)
for (int j = 0 ; j < 3 ; j ++ )
_b [i][j] = m._b [i][j] ;
}
{for (int i = 0 ; i < 3 ; i ++)
for (int j = 0 ; j < 3 ; j ++ )
_n [i][j] = m._n [i][j] ;
}
return ;
}
void BilinearSurfaceMapping ::
map2world (const coord2_t& map, coord3_t& wld) const {
double x = map [0];
double y = map [1];
double xy = x * y ;
wld[0] = _b [0][0] + x * _b [1][0] + y * _b [2][0] + xy * _b [3][0] ;
wld[1] = _b [0][1] + x * _b [1][1] + y * _b [2][1] + xy * _b [3][1] ;
wld[2] = _b [0][2] + x * _b [1][2] + y * _b [2][2] + xy * _b [3][2] ;
return ;
}
void BilinearSurfaceMapping ::
map2world (double x, double y, coord3_t& w) const {
coord2_t p ;
p [0] = x ;
p [1] = y ;
map2world (p,w) ;
return ;
}
void BilinearSurfaceMapping ::
normal (const coord2_t& map, coord3_t& normal) const {
double x = map [0];
double y = map [1];
normal [0] = -( _n [0][0] + _n [1][0] * x + _n [2][0] * y) ;
normal [1] = -( _n [0][1] + _n [1][1] * x + _n [2][1] * y) ;
normal [2] = -( _n [0][2] + _n [1][2] * x + _n [2][2] * y) ;
return ;
}
*/
} // end namespace Dune
grid/alu3dgrid/alu3dmappings.hh
+ 1
− 1
View file @ 7929e086
......@@ -4,7 +4,6 @@
#define DUNE_ALU3DMAPPINGS_HH
//- Dune includes
//#include <dune/grid/alu3dgrid.hh>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
......@@ -66,6 +65,7 @@ namespace Dune {
void normal(const coord2_t&, coord3_t&) const ;
} ;
} // end namespace Dune
#endif
grid/alu3dgrid/datahandle.hh
+ 6
− 0
View file @ 7929e086
......@@ -3,6 +3,12 @@
#ifndef __DUNE_ALU3DGRID_DATAHANDLE_HH__
#define __DUNE_ALU3DGRID_DATAHANDLE_HH__
#include <iostream>
using std::endl;
using std::cout;
using std::flush;
namespace ALU3dGridSpace {
//! the corresponding interface class is defined in bsinclude.hh
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment