diff --git a/common/misc.hh b/common/misc.hh
index d621f0a8532787ffe13519d5e8c24ff234529c0c..6775461ba159269454e9d19701f7fd35cc787ec0 100644
--- a/common/misc.hh
+++ b/common/misc.hh
@@ -6,13 +6,14 @@
#include <iostream>
#include <math.h>
-template <bool flag>
-class CompileTimeChecker { };
+//! Check condition at compilation time
+template <bool flag> class CompileTimeChecker;
+//! it exists only an implementation for true so the compiler throws an
+//! error if the condition is false
template <> class CompileTimeChecker<true> { };
-
namespace Dune {
/** @addtogroup Common
diff --git a/fem/lagrangebase.hh b/fem/lagrangebase.hh
index b96ee76ee2f3ccfd9ab60eb545dfbdc084cafbaa..301c39e1963be1565fc1b6439fc10458be0792f9 100644
--- a/fem/lagrangebase.hh
+++ b/fem/lagrangebase.hh
@@ -97,9 +97,6 @@ namespace Dune {
}
};
-#if 0
-
- // not in use
//*****************************************************************
//
@@ -114,20 +111,6 @@ namespace Dune {
//! (0,0) (1,0)
//
//*****************************************************************
-#endif
- //*****************************************************************
- //
- //! (0,1)
- //! 2|\ coordinates and local node numbers
- //! | \
- // //! | \
- //! | \
- // //! | \
- // //! | \
- // //! 0|______\1
- //! (0,0) (1,0)
- //
- //*****************************************************************
template<class FunctionSpaceType>
class LagrangeBaseFunction < FunctionSpaceType , triangle , 1 >
: public BaseFunctionInterface<FunctionSpaceType>
@@ -136,7 +119,6 @@ namespace Dune {
RangeField factor[3];
public:
-#if 0
LagrangeBaseFunction ( FunctionSpaceType & f , int baseNum )
: BaseFunctionInterface<FunctionSpaceType> (f)
{
@@ -156,42 +138,19 @@ namespace Dune {
factor[i] = 0.0;
}
}
-#endif
-
- // this is the version with with phi(x,y) = x as base function 1
- LagrangeBaseFunction ( FunctionSpaceType & f , int baseNum )
- : BaseFunctionInterface<FunctionSpaceType> (f)
- {
- if(baseNum == 0)
- { // 1 - x - y
- factor[0] = 1.0;
- factor[1] = -1.0;
- factor[2] = -1.0;
- }
- else
- {
- factor[2] = 0.0;
- for(int i=1; i<3; i++) // x , y
- if(baseNum == i)
- factor[i] = 1.0;
- else
- factor[i] = 0.0;
- }
- }
-
virtual void evaluate ( const Vec<0, deriType> &diffVariable,
const Domain & x, Range & phi) const
{
- phi = factor[0];
- for(int i=1; i<3; i++)
- phi += factor[i] * x.get(i-1);
+ phi = factor[2];
+ for(int i=0; i<2; i++)
+ phi += factor[i] * x.get(i);
}
virtual void evaluate ( const Vec<1, deriType> &diffVariable,
const Domain & x, Range & phi) const
{
// x or y ==> 1 or 2
- int num = diffVariable.get(0)+1;
+ int num = diffVariable.get(0);
phi = factor[num];
}
@@ -201,7 +160,6 @@ namespace Dune {
// function is linear, therfore
phi = 0.0 ;
}
-
};
//*****************************************************************
diff --git a/grid/albertgrid.hh b/grid/albertgrid.hh
index 41e47d36b0801ba869af4fd61a76b34175e0a8ad..4c99f816a4124a907a1f9b79ed9cc9299330d40a 100644
--- a/grid/albertgrid.hh
+++ b/grid/albertgrid.hh
@@ -3,8 +3,8 @@
#ifndef DUNE_ALBERTGRID_HH
#define DUNE_ALBERTGRID_HH
-#include "../common/matvec.hh"
#include "../common/misc.hh"
+#include "../common/matvec.hh"
#include "../common/array.hh"
#include "common/grid.hh"
@@ -49,6 +49,18 @@ namespace Albert
#include <albert.h>
#include "albertgrid/albertextra.hh"
+#ifdef ABS
+#undef ABS
+#endif
+
+#ifdef MIN
+#undef MIN
+#endif
+
+#ifdef MAX
+#undef MAX
+#endif
+
#ifndef __ALBERTNAME__
} // end extern "C"
#endif
@@ -162,13 +174,6 @@ namespace Albert
//! local coordinate in its reference element
Vec<dim,albertCtype>& local (const Vec<dimworld,albertCtype>& global);
- template <int dimbary>
- Vec<dimbary,albertCtype>& localB (const Vec<dimworld,albertCtype>& global)
- {
- localBary_ = localBary(global);
- return localBary_;
- }
-
//! returns true if the point is in the current element
bool pointIsInside(const Vec<dimworld,albertCtype> &point);
@@ -201,23 +206,32 @@ namespace Albert
//! can only be called for dim=dimworld!
Mat<dim,dim>& Jacobian_inverse (const Vec<dim,albertCtype>& local);
+ //***********************************************************************
+ // Methods that not belong to the Interface, but have to be public
+ //***********************************************************************
+ bool builtGeom(ALBERT EL_INFO *elInfo, unsigned char face,
+ unsigned char edge, unsigned char vertex);
+ void initGeom();
+
+ private:
//! print internal data
void print (std::ostream& ss, int indent);
+#if 0
//! return unit outer normal of this element, work for Faces nad Edges
Vec<dimworld,albertCtype>& unit_outer_normal ();
//! return outer normal of this element, work for Faces nad Edges
Vec<dimworld,albertCtype>& outer_normal ();
+#endif
- //***********************************************************************
- // Methods that not belong to the Interface, but have to be public
- //***********************************************************************
- bool builtGeom(ALBERT EL_INFO *elInfo, unsigned char face,
- unsigned char edge, unsigned char vertex);
- void initGeom();
+ template <int dimbary>
+ Vec<dimbary,albertCtype>& localB (const Vec<dimworld,albertCtype>& global)
+ {
+ localBary_ = localBary(global);
+ return localBary_;
+ }
- private:
//! built the reference element
void makeRefElemCoords();
@@ -267,10 +281,16 @@ namespace Albert
}
+ //! the vertex coordinates
Mat<dimworld,dim+1,albertCtype> coord_;
+ //! storage for global coords
Vec<dimworld,albertCtype> globalCoord_;
+
+ //! storage for local coords
Vec<dim,albertCtype> localCoord_;
+
+ //! storage for barycentric coords
Vec<dimbary,albertCtype> localBary_;
@@ -287,11 +307,12 @@ namespace Albert
//! Which Edge of the Face of the Element 0...dim-1
unsigned char vertex_;
- Mat<dim,dim,albertCtype> Jinv_; //!< storage for inverse of jacobian
- albertCtype volume_;
+ //! is true if Jinv_ and volume_ is calced
bool builtinverse;
+ Mat<dim,dim,albertCtype> Jinv_; //!< storage for inverse of jacobian
+ albertCtype volume_; //!< storage of element volume
- Vec<dimworld,albertCtype> outerNormal_;
+ //Vec<dimworld,albertCtype> outerNormal_;
};
@@ -451,6 +472,7 @@ namespace Albert
/*! Intra-element access to entities of codimension cc > codim. Return number of entities
with codimension cc.
*/
+
#if 0
//!< Default codim 1 Faces and codim == dim Vertices
template <int cc> int count () { return ((dim-codim)+1); };
@@ -581,7 +603,8 @@ namespace Albert
//!< Default codim 1 Faces and codim == dim Vertices
template<int cc> int count () { return (dim+1); };
- // specialization only for codim == 2 , edges
+ // specialization only for codim == 2 , edges,
+ // a tetrahedron has always 6 edges
template<> int count<2 << dim> () { return 6; };
@@ -617,14 +640,14 @@ namespace Albert
//! Assumes that meshes are nested.
AlbertGridLevelIterator<0,dim,dimworld> father ();
- /*! Location of this element relative to the reference element element of the father.
- This is sufficient to interpolate all dofs in conforming case.
- Nonconforming may require access to neighbors of father and
- computations with local coordinates.
- On the fly case is somewhat inefficient since dofs are visited several times.
- If we store interpolation matrices, this is tolerable. We assume that on-the-fly
- implementation of numerical algorithms is only done for simple discretizations.
- Assumes that meshes are nested.
+ /*! Location of this element relative to the reference element
+ of the father. This is sufficient to interpolate all
+ dofs in conforming case. Nonconforming may require access to
+ neighbors of father and computations with local coordinates.
+ On the fly case is somewhat inefficient since dofs are visited
+ several times. If we store interpolation matrices, this is tolerable.
+ We assume that on-the-fly implementation of numerical algorithms
+ is only done for simple discretizations. Assumes that meshes are nested.
*/
AlbertGridElement<dim,dim>& father_relative_local ();
@@ -710,11 +733,11 @@ namespace Albert
//! know your own dimension of world
enum { dimensionworld=dimworld };
- // the normal Constructor
+ //! the normal Constructor
AlbertGridHierarchicIterator(AlbertGrid<dim,dimworld> &grid,
ALBERT TRAVERSE_STACK *travStack, int travLevel);
- // the default Constructor
+ //! the default Constructor
AlbertGridHierarchicIterator(AlbertGrid<dim,dimworld> &grid);
//! prefix increment
@@ -818,13 +841,6 @@ namespace Albert
//! return unit outer normal, if you know it is constant use this function instead
Vec<dimworld,albertCtype>& unit_outer_normal ();
- //! return outer normal, this should be dependent on local
- //! coordinates for higher order boundary
- Vec<dimworld,albertCtype>& outer_normal (Vec<dim-1,albertCtype>& local);
-
- //! return unit outer normal, if you know it is constant use this function instead
- Vec<dimworld,albertCtype>& outer_normal ();
-
//! intersection of codimension 1 of this neighbor with element where
//! iteration started.
//! Here returned element is in LOCAL coordinates of the element
@@ -850,6 +866,13 @@ namespace Albert
int number_in_neighbor ();
private:
+ //! return outer normal, this should be dependent on local
+ //! coordinates for higher order boundary
+ Vec<dimworld,albertCtype>& outer_normal (Vec<dim-1,albertCtype>& local);
+
+ //! return unit outer normal, if you know it is constant use this function instead
+ Vec<dimworld,albertCtype>& outer_normal ();
+
//! setup the virtual entity
void setupVirtualEntity(int neighbor);
@@ -969,7 +992,7 @@ namespace Albert
};
template <> ALBERT EL_INFO *
- goNextEntity<1>(ALBERT TRAVERSE_STACK *stack,ALBERT EL_INFO *elinfo_old)
+ goNextEntity<1 << dim>(ALBERT TRAVERSE_STACK *stack,ALBERT EL_INFO *elinfo_old)
{
return goNextFace(stack,elinfo_old);
};
@@ -1035,6 +1058,13 @@ namespace Albert
//friend class AlbertGridEntity <1 << dim-1 ,dim,dimworld>;
friend class AlbertGridEntity <dim,dim,dimworld>;
+ //! AlbertGrid is only implemented for 2 and 3 dimension
+ //! for 1d use SGrid or SimpleGrid
+ CompileTimeChecker<dimworld != 1> Do_not_use_AlbertGrid_for_1d_Grids;
+
+ //! At this time AlbertGrid only supports the dim == dimworld case
+ CompileTimeChecker<dimworld == dim> Only_dim_equal_dimworld_is_implemented;
+
//**********************************************************
// The Interface Methods
//**********************************************************
@@ -1057,12 +1087,7 @@ namespace Albert
//! Return maximum level defined in this grid. Levels are numbered
//! 0 ... maxlevel with 0 the coarsest level.
- int maxlevel();
-
- int maxlevel() const
- {
- return const_cast<AlbertGrid<dim,dimworld> *> (this)->maxlevel();
- };
+ int maxlevel() const;
//! Iterator to first entity of given codim on level
template<int codim>
@@ -1074,11 +1099,7 @@ namespace Albert
//! number of grid entities per level and codim
int size (int level, int codim) ;
-
- int size (int level, int codim) const
- {
- return const_cast<AlbertGrid<dim,dimworld> *> (this)->size(level,codim);
- }
+ int size (int level, int codim) const;
//**********************************************************
// End of Interface Methods
diff --git a/grid/albertgrid/albertgrid.cc b/grid/albertgrid/albertgrid.cc
index f8f1ad4b63eecd8bb8fcbd5cc52ceb5ea9b39180..c83e62cf48e5e9acf742d2f8bea84db342023428 100644
--- a/grid/albertgrid/albertgrid.cc
+++ b/grid/albertgrid/albertgrid.cc
@@ -491,12 +491,15 @@ namespace Dune
inline Vec<dimworld,albertCtype>& AlbertGridElement<dim,dimworld>::
global(const Vec<dim>& local)
{
- // Umrechnen von localen Koordinaten zu baryzentrischen Koordinaten
- Vec<dim+1> tmp (1.0); // Wichtig, da tmp(0) = 1 - tmp(1)- ... -tmp(dim+1)
+ // only the dim first components are needed, because we don't use
+ // barycentric coordinates
+ Vec<dim+1,albertCtype> tmp;
for(int i=0; i<dim; i++)
- tmp(0) -= local.read(i);
- for(int i=1; i<dim+1; i++)
- tmp(i) = local.read(i-1);
+ tmp(i) = local.read(i);
+
+ tmp(dim) = 1.0;
+ for(int i=0; i<dim; i++)
+ tmp(dim) -= local.read(i);
// globale Koordinaten ausrechnen
globalCoord_ = globalBary(tmp);
@@ -532,12 +535,14 @@ namespace Dune
inline Vec<dim>& AlbertGridElement<dim,dimworld>::
local(const Vec<dimworld>& global)
{
+ ALBERT REAL lambda[dim+1];
+
Vec<dim+1,albertCtype> tmp = localBary(global);
// Umrechnen von baryzentrischen localen Koordinaten nach
// localen Koordinaten,
for(int i=0; i<dim; i++)
- localCoord_(i) = tmp(i+1);
+ localCoord_(i) = tmp(i);
return localCoord_;
}
@@ -564,7 +569,7 @@ namespace Dune
ALBERT REAL *v = NULL;
/*
- * wir haben das gleichungssystem zu loesen:
+ * we got to solve the problem :
*/
/*
* ( q1x q2x ) (lambda1) = (qx)
@@ -573,7 +578,7 @@ namespace Dune
* ( q1y q2y ) (lambda2) = (qy)
*/
/*
- * mit qi=pi-p3, q=global-p3
+ * with qi=pi-p3, q=global-p3
*/
v = static_cast<ALBERT REAL *> (elInfo_->coord[0]);
@@ -616,11 +621,11 @@ namespace Dune
Vec<dim+1,albertCtype> lambda;
int j, k;
- //! wir haben das gleichungssystem zu loesen:
+ //! we got to solve the problem :
//! ( q1x q2x q3x) (lambda1) (qx)
//! ( q1y q2y q3y) (lambda2) = (qy)
//! ( q1z q2z q3z) (lambda3) (qz)
- //! mit qi=pi-p3, q=xy-p3
+ //! with qi=pi-p3, q=xy-p3
for (int j = 0; j < dimworld; j++)
{
@@ -728,6 +733,7 @@ namespace Dune
return ret;
}
+#if 0
template< int dim, int dimworld>
inline Vec<dimworld,albertCtype>& AlbertGridElement<dim,dimworld>::
unit_outer_normal()
@@ -785,6 +791,7 @@ namespace Dune
return outerNormal_;
}
+#endif
#if 0
//************************************************************************
@@ -1107,94 +1114,6 @@ namespace Dune
edgeEntity_ = NULL;
}
-#if 0
- //****************************************
- //
- // specialization of count and entity
- //
- //****************************************
- // specialization for codim == 2 , edges
- template <>
- inline int AlbertGridEntity<0,3,3>::count<2> ()
- {
- // number of edges of a tetrahedron
- return 6;
- }
-
- //****************************
- // the vertex entitys
- template <>
- inline AlbertGridLevelIterator<2,2,2>& AlbertGridEntity <0,2,2>::
- entity<2> ( int i)
- {
- vxEntity_.virtualEntity_.setElInfo(elInfo_,elInfo_->el->dof[i][0],0,0,i);
- return vxEntity_;
- }
- template <>
- inline AlbertGridLevelIterator<3,3,3>& AlbertGridEntity <0,3,3>::
- entity<3> ( int i)
- {
- vxEntity_.virtualEntity_.setElInfo(elInfo_,elInfo_->el->dof[i][0],0,0,i);
- return vxEntity_;
- }
-
- //*********************
- // the face access
- template <>
- inline AlbertGridLevelIterator<1,2,2>& AlbertGridEntity <0,2,2>::
- entity<1> ( int i)
- {
- std::cout << "Not checked yet! \n";
- if(!faceEntity_)
- {
- faceEntity_ = new AlbertGridLevelIterator<1,2,2>
- ( grid_ , elInfo_,index(),i,0,0 );
- return (*faceEntity_);
- }
-
- faceEntity_->virtualEntity_.setElInfo(elInfo_,index(),i,0,0);
- return (*faceEntity_);
- }
-
- template <>
- inline AlbertGridLevelIterator<1,3,3>& AlbertGridEntity <0,3,3>::
- entity<1> ( int i)
- {
- std::cout << "Not checked yet! \n";
- if(!faceEntity_)
- {
- faceEntity_ = new AlbertGridLevelIterator<1,3,3>
- ( grid_ , elInfo_,index(),i,0,0 );
- return (*faceEntity_);
- }
-
- faceEntity_->virtualEntity_.setElInfo(elInfo_,index(),i,0,0);
- return (*faceEntity_);
- }
- //*********************
- // the edge access, only for dim == 3
- template <>
- inline AlbertGridLevelIterator<2,3,3>& AlbertGridEntity <0,3,3>::
- entity<2> ( int i)
- {
- printf("Entity::entity<codim = %d>: Warning elNum may be not correct! \n",2);
- if(!edgeEntity_)
- {
- edgeEntity_ = new AlbertGridLevelIterator<2,3,3>
- ( grid_ , NULL , 0,0,0,0);
- }
-
- if(i < 3)
- { // 0,1,2
- edgeEntity_->virtualEntity_.setElInfo(elInfo_,index(),0,i,0);
- }
- else
- { // 3,4,5
- edgeEntity_->virtualEntity_.setElInfo(elInfo_,index(),i-2,1,0);
- }
- return (*edgeEntity_);
- }
-#endif
template <int dim, int dimworld>
inline AlbertGridLevelIterator<1,dim,dimworld>& AlbertGridEntity <0,dim,dimworld>::
faceEntity ( int i)
@@ -1214,6 +1133,10 @@ namespace Dune
inline AlbertGridLevelIterator<2,dim,dimworld>& AlbertGridEntity <0,dim,dimworld>::
edgeEntity ( int i)
{
+ // check whether dim == 3 or not, because this method is only for
+ // dim == 3 and codim == 2
+ CompileTimeChecker<dim == 3> edgeEntity_only_for_dim_equal_3;
+
printf("Entity::entity<codim = %d>: Warning elNum may be not correct! \n",2);
if(!edgeEntity_)
{
@@ -1291,15 +1214,16 @@ namespace Dune
inline AlbertGridLevelIterator<0,dim,dimworld>
AlbertGridEntity < 0, dim ,dimworld >::father()
{
- ALBERT TRAVERSE_STACK travStack;
- initTraverseStack(&travStack);
-
- travStack = (*travStack_);
-
- travStack.stack_used--;
+ ALBERT EL_INFO * fatherInfo = NULL;
+ // if level <= 0 return father = this entity
+ if(elInfo_->level <= 0)
+ fatherInfo = &travStack_->elinfo_stack[travStack_->stack_used-1];
+ else
+ fatherInfo = elInfo_;
+ // new LevelIterator with EL_INFO one level above
AlbertGridLevelIterator <0,dim,dimworld>
- it(travStack.elinfo_stack+travStack.stack_used);
+ it(grid_,fatherInfo,grid_.indexOnLevel<0>(fatherInfo->el->index,fatherInfo->level),0,0,0);
return it;
}
@@ -1307,7 +1231,15 @@ namespace Dune
inline AlbertGridElement<dim,dim>&
AlbertGridEntity < 0, dim ,dimworld >::father_relative_local()
{
- std::cout << "\nfather_realtive_local not implemented yet! \n";
+ //AlbertGridLevelIterator<0,dim,dimworld> daddy = father();
+ AlbertGridElement<dim,dimworld> daddy = (*father()).geometry();
+
+ fatherReLocal_.initGeom();
+ // compute the local coordinates in father refelem
+ for(int i=0; i<fatherReLocal_.corners(); i++)
+ fatherReLocal_[i] = daddy.local(geometry()[i]);
+
+ std::cout << "\nfather_realtive_local not tested yet! \n";
return fatherReLocal_;
}
@@ -2536,11 +2468,17 @@ namespace Dune
template < int dim, int dimworld >
- inline int AlbertGrid < dim, dimworld >::maxlevel()
+ inline int AlbertGrid < dim, dimworld >::maxlevel() const
{
return maxlevel_;
}
+ template < int dim, int dimworld >
+ inline int AlbertGrid < dim, dimworld >::size (int level, int codim) const
+ {
+ return const_cast<AlbertGrid<dim,dimworld> *> (this)->size(level,codim);
+ }
+
template < int dim, int dimworld >
inline int AlbertGrid < dim, dimworld >::size (int level, int codim)
{