Skip to content
Snippets Groups Projects
Commit 3b4e8b21 authored by Peter Bastian's avatar Peter Bastian
Browse files

everything implemented now 

[[Imported from SVN: r9]]
parent b7448b48
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
grid/sgrid/numbering.cc
+ 279
− 257
View file @ 3b4e8b21
// NOTE: The current revision of this file was left untouched when the DUNE source files were reindented!
// NOTE: It contained invalid syntax that could not be processed by uncrustify.
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef __NUMBERING_CC__
#define __NUMBERING_CC__
namespace Dune {
template<int dim>
inline void LexOrder<dim>::init (Tupel<int,dim>& _N)
{
// store argument
N=_N;
// build P array
P[0] = 1;
for (int i=1; i<=dim; i++) P[i] = P[i-1]*N[i-1];
}
template<int dim>
inline int LexOrder<dim>::tupels ()
{
return P[dim];
}
template<int dim>
inline int LexOrder<dim>::n (Tupel<int,dim>& z)
{
int r=0;
for (int i=0; i<dim; i++) r += z[i]*P[i];
return r;
}
template<int dim>
inline Tupel<int,dim> LexOrder<dim>::z (int n)
{
Tupel<int,dim> z;
for (int i=0; i<dim; i++)
{
z[i] = n%N[i];
n = n/N[i];
}
return z;
}
//************************************************************************
template<int dim>
inline void JoinOrder<dim>::init (Tupel<int,dim>& _N)
{
// store argument
N=_N;
// build P array
offset[0] = 0;
for (int i=1; i<=dim; i++) offset[i] = offset[i-1]+N[i-1];
}
template<int dim>
inline int JoinOrder<dim>::size ()
{
return offset[dim];
}
template<int dim>
inline int JoinOrder<dim>::n (int subset, int index)
{
return index+offset[subset];
}
template<int dim>
inline int JoinOrder<dim>::index (int n)
{
for (int i=0; i<dim; i++)
if (N[i]!=0)
{
if (n<N[i]) return n;
n -= N[i];
}
return n; // oops, error !
}
template<int dim>
inline int JoinOrder<dim>::subset (int n)
{
for (int i=0; i<dim; i++)
if (N[i]!=0)
{
if (n<N[i]) return i;
n -= N[i];
}
return 0; // oops, error !
}
//************************************************************************
template<int dim>
CubeMapper<dim>::CubeMapper (Tupel<int,dim>& _N)
{
// store argument
N=_N;
// preprocess binary partitions
for (int i=0; i<=dim; i++) ne[i] = 0;
for (int b=0; b<power2(dim); b++) // loop over all binary partitions
{
int mask=1;
Tupel<int,dim> t;
for (int i=0; i<dim; i++)
{
if (b&mask)
{
// bit i is even
t[i] = N[i]+1;
}
else
{
// bit i is odd
t[i] = N[i];
}
mask *= 2;
}
lex[b].init(t); // set up lex ordering of tupels
nb[b] = lex[b].tupels();
cb[b] = ones(b);
ne[cb[b]] += nb[b];
}
// preprocess lex ordering for each codimension
for (int c=0; c<=dim; c++)
{
Tupel<int,1<<dim> t;
for (int b=0; b<power2(dim); b++) // loop over all binary partitions
if (ones(b)==c)
{
// partition b is part of codim c
t[b] = nb[b];
}
else
{
// partition does not contribute to codim c
t[b] = 0;
}
join[c].init(t); // set join mapper
}
cout << "CubeMapper [" ;
for (int i=0; i<dim; i++)
cout << N[i] << " ";
cout << "]" << endl;
for (int i=0; i<=dim; i++)
cout << " " << ne[i] << " elements of codim " << i
<< " in dimension " << dim << endl;
}
template<int dim>
inline int CubeMapper<dim>::elements (int codim)
{
return ne[codim];
}
template<int dim>
inline int CubeMapper<dim>::codim (Tupel<int,dim>& z)
{
int r=0;
for (int i=0; i<dim; i++)
if (z[i]%2==0) r++; // count even components
return r;
}
template<int dim>
inline int CubeMapper<dim>::n (Tupel<int,dim>& z)
{
int p = partition(z); // get partition
Tupel<int,dim> r=reduce(z); // get reduced coordinate
// treat easy cases first
if (p==0 || p==power2(dim)-1)
{
// all components are either odd or even
return lex[p].n(r);
}
// general case
return join[ones(p)].n(p,lex[p].n(r));
}
template<int dim>
inline Tupel<int,dim> CubeMapper<dim>::z (int i, int codim)
{
Tupel<int,dim> r;
// easy cases first
if (codim==0)
{
r = lex[0].z(i);
return expand(r,0);
}
if (codim==dim)
{
r = lex[power2(dim)-1].z(i);
return expand(r,power2(dim)-1);
}
// general case
int p = join[codim].subset(i);
int n = join[codim].index(i);
r = lex[p].z(n);
return expand(r,p);
}
template<int dim>
inline int CubeMapper<dim>::ones (int b)
{
int r = 0;
int mask = 1;
for (int i=0; i<dim; i++)
{
if (b&mask) r++;
mask *=2 ;
}
return r;
}
template<int dim>
inline int CubeMapper<dim>::partition (Tupel<int,dim>& z)
{
int r = 0;
int mask = 1;
for (int i=0; i<dim; i++)
{
if (z[i]%2==0) r += mask;
mask *=2 ;
}
return r;
}
template<int dim>
inline Tupel<int,dim> CubeMapper<dim>::reduce (Tupel<int,dim>& z)
{
Tupel<int,dim> r;
for (int i=0; i<dim; i++)
if (z[i]%2==0)
r[i] = z[i]/2; // even component
else
r[i] = (z[i]-1)/2; // odd component
return r;
}
template<int dim>
inline Tupel<int,dim> CubeMapper<dim>::expand (Tupel<int,dim>& r, int b)
{
Tupel<int,dim> z;
int mask = 1;
for (int i=0; i<dim; i++)
if (b&mask)
z[i] = r[i]*2; // even component
else
z[i] = 2*r[i]+1); // odd component
return r;
}
template<int dim>
inline void LexOrder<dim>::init (Tupel<int,dim>& _N)
{
// store argument
N=_N;
// build P array
P[0] = 1;
for (int i=1; i<=dim; i++) P[i] = P[i-1]*N[i-1];
}
template<int dim>
inline int LexOrder<dim>::tupels ()
{
return P[dim];
}
template<int dim>
inline int LexOrder<dim>::n (Tupel<int,dim>& z)
{
int r=0;
for (int i=0; i<dim; i++) r += z[i]*P[i];
return r;
}
template<int dim>
inline Tupel<int,dim> LexOrder<dim>::z (int n)
{
Tupel<int,dim> z;
for (int i=0; i<dim; i++)
{
z[i] = n%N[i];
n = n/N[i];
}
return z;
}
//************************************************************************
template<int dim>
inline void JoinOrder<dim>::init (Tupel<int,dim>& _N)
{
// store argument
N=_N;
// build P array
offset[0] = 0;
for (int i=1; i<=dim; i++) offset[i] = offset[i-1]+N[i-1];
}
template<int dim>
inline int JoinOrder<dim>::size ()
{
return offset[dim];
}
template<int dim>
inline int JoinOrder<dim>::n (int subset, int index)
{
return index+offset[subset];
}
template<int dim>
inline int JoinOrder<dim>::index (int n)
{
for (int i=0; i<dim; i++)
if (N[i]!=0)
{
if (n<N[i]) return n;
n -= N[i];
}
return n; // oops, error !
}
template<int dim>
inline int JoinOrder<dim>::subset (int n)
{
for (int i=0; i<dim; i++)
if (N[i]!=0)
{
if (n<N[i]) return i;
n -= N[i];
}
return 0; // oops, error !
}
//************************************************************************
template<int dim>
CubeMapper<dim>::CubeMapper (Tupel<int,dim> _N)
{
make(_N);
}
template<int dim>
CubeMapper<dim>::CubeMapper ()
{
Tupel<int,dim> M;
// make mesh of single cube
for (int i=0; i<dim; i++) M[i]=1;
make(M);
}
template<int dim>
void CubeMapper<dim>::make (Tupel<int,dim>& _N)
{
// store argument
N=_N;
// preprocess binary partitions
for (int i=0; i<=dim; i++) ne[i] = 0;
for (int b=0; b<power2(dim); b++) // loop over all binary partitions
{
int mask=1;
Tupel<int,dim> t;
for (int i=0; i<dim; i++)
{
if (b&mask)
{
// bit i is even
t[i] = N[i]+1;
}
else
{
// bit i is odd
t[i] = N[i];
}
mask *= 2;
}
lex[b].init(t); // set up lex ordering of tupels
nb[b] = lex[b].tupels();
cb[b] = ones(b);
ne[cb[b]] += nb[b];
}
// preprocess lex ordering for each codimension
for (int c=0; c<=dim; c++)
{
Tupel<int,1<<dim> t;
for (int b=0; b<power2(dim); b++) // loop over all binary partitions
if (ones(b)==c)
{
// partition b is part of codim c
t[b] = nb[b];
}
else
{
// partition does not contribute to codim c
t[b] = 0;
}
join[c].init(t); // set join mapper
}
}
template<int dim>
inline void CubeMapper<dim>::print (std::ostream& ss, int indent)
{
for (int k=0; k<indent; k++) ss << " ";
ss << "CubeMapper [" ;
for (int i=0; i<dim; i++)
cout << N[i] << " ";
ss << "]" << endl;
for (int i=0; i<=dim; i++)
{
for (int k=0; k<indent; k++) ss << " ";
ss << " " << ne[i] << " elements of codim " << i
<< " in dimension " << dim << endl;
}
}
template<int dim>
inline int CubeMapper<dim>::elements (int codim)
{
return ne[codim];
}
template<int dim>
inline int CubeMapper<dim>::codim (Tupel<int,dim>& z)
{
int r=0;
for (int i=0; i<dim; i++)
if (z[i]%2==0) r++; // count even components
return r;
}
template<int dim>
inline int CubeMapper<dim>::n (Tupel<int,dim>& z)
{
int p = partition(z); // get partition
Tupel<int,dim> r=compress(z); // get compressd coordinate
// treat easy cases first
if (p==0 || p==power2(dim)-1)
{
// all components are either odd or even
return lex[p].n(r);
}
// general case
return join[ones(p)].n(p,lex[p].n(r));
}
template<int dim>
inline Tupel<int,dim> CubeMapper<dim>::z (int i, int codim)
{
Tupel<int,dim> r;
// easy cases first
if (codim==0)
{
r = lex[0].z(i);
return expand(r,0);
}
if (codim==dim)
{
r = lex[power2(dim)-1].z(i);
return expand(r,power2(dim)-1);
}
// general case
int p = join[codim].subset(i);
int n = join[codim].index(i);
r = lex[p].z(n);
return expand(r,p);
}
template<int dim>
inline int CubeMapper<dim>::ones (int b)
{
int r = 0;
int mask = 1;
for (int i=0; i<dim; i++)
{
if (b&mask) r++;
mask *=2 ;
}
return r;
}
template<int dim>
inline int CubeMapper<dim>::partition (Tupel<int,dim>& z)
{
int r = 0;
int mask = 1;
for (int i=0; i<dim; i++)
{
if (z[i]%2==0) r += mask;
mask *=2 ;
}
return r;
}
template<int dim>
inline Tupel<int,dim> CubeMapper<dim>::compress (Tupel<int,dim>& z)
{
Tupel<int,dim> r;
for (int i=0; i<dim; i++)
if (z[i]%2==0)
r[i] = z[i]/2; // even component
else
r[i] = (z[i]-1)/2; // odd component
return r;
}
template<int dim>
inline Tupel<int,dim> CubeMapper<dim>::expand (Tupel<int,dim>& r, int b)
{
Tupel<int,dim> z;
for (int i=0; i<dim; i++)
if (b&(1<<i))
z[i] = r[i]*2; // even component
else
z[i] = 2*r[i]+1; // odd component
return z;
}
}
#endif
grid/sgrid/numbering.hh
+ 93
− 5
View file @ 3b4e8b21
......@@ -12,6 +12,57 @@ namespace Dune {
//! make empty tupel
Tupel() {}
//! Tupel from components
Tupel(T t0)
{
for (int i=0; i<d; i++) x[i] = t0;
}
//! Tupel from components
Tupel(T t0, T t1)
{
x[0] = t0;
x[1] = t1;
}
//! Tupel from components
Tupel(T t0, T t1, T t2)
{
x[0] = t0;
x[1] = t1;
x[2] = t2;
}
//! Tupel from components
Tupel(T t0, T t1, T t2, T t3)
{
x[0] = t0;
x[1] = t1;
x[2] = t2;
x[3] = t3;
}
//! Tupel from components
Tupel(T t0, T t1, T t2, T t3, T t4)
{
x[0] = t0;
x[1] = t1;
x[2] = t2;
x[3] = t3;
x[4] = t4;
}
//! Tupel from components
Tupel(T t0, T t1, T t2, T t3, T t4, T t5)
{
x[0] = t0;
x[1] = t1;
x[2] = t2;
x[3] = t3;
x[4] = t4;
x[5] = t5;
}
//! read/write components
int& operator[] (int i) {return x[i];}
......@@ -64,12 +115,40 @@ namespace Dune {
int offset[dim+1]; // P[i] = Sum_{i=0}^{i} N[i];
};
//! associate a unique consecutive index to all entities of a given codimension in a d-dimension cube
/*! The CubeMapper assigns an id to all entities of all codimensions of a structured mesh
with arbitrary number of elements (codim 0 entities) in each direction. The ids
are unique and consecutive within each codimension.
The idea is as follows: Consider a structured mesh in \f$d\f$ dimensions with \f$N\f$ elements per
direction. This mesh has \f$N^d\f$ elements in total. Now imagine refined mesh where each element
is halfened in every coordinate direction. This refined mesh has \f$(2N+1)^d\f$ vertices (entities
of codimension \f$d\f$). Each vertex of the refined mesh now corresponds to a grid entity of the
original mesh. Moreover, a vertex in the refined mesh can be identified by integer coordintes \f$z\f$
where \f$z_i\in\{0,\ldots,2N\}, 0\leq i < d\f$. Let \f$c(z)\f$ be the number of even components in \f$z\f$.
Then, \f$c(z)\f$ is the codimension of the mesh entity with coordinate \f$z\f$. E.~g.~
entities of codimension 0 have odd coordinates, all entities of codim \f$d\f$ have \f$d\f$
even coordinates.
In order to number all entities of one codimension consecutively we observe that the
refined mesh can be subdivided into \f$2^d\f$subsets. Subset number \f$b\f$ with
binary representation \f$(b_{d-1},\ldots,b_0)\f$ corresponds to all \f$z\in [0,2N]^d\f$ where
\f$z_i\f$ is even if \f$b_i\f$ is 1 and \f$z_i\f$ is odd if \f$b_i\f$ is 0. The entities of codimension
\f$c\f$ now consist of \f$\left ( \begin{array}{cc}d\\c\end{array} \right)\f$ of those
subsets. Within the subsets the numbering is lexicographic and then the corrsponding
subsets are numbered consecutively.
*/
template<int dim>
class CubeMapper {
public:
//! construct with number of elements (of codim 0) in each direction
CubeMapper (Tupel<int,dim>& _N);
CubeMapper (Tupel<int,dim> _N);
//! make cube of single element
CubeMapper ();
//! (re)initialize with number of elements (of codim 0) in each direction
void make (Tupel<int,dim>& _N);
//! get number of elements in each codimension
int elements (int codim);
......@@ -85,6 +164,18 @@ namespace Dune {
//! compute coordinates from number and codimension
Tupel<int,dim> z (int i, int codim);
//! compress from expanded coordinates to grid for a single partition number
Tupel<int,dim> compress (Tupel<int,dim>& z);
//! expand with respect to partition number
Tupel<int,dim> expand (Tupel<int,dim>& r, int b);
//! There are \f$2^d\f$ possibilities of having even/odd coordinates. The binary representation is called partition number
int partition (Tupel<int,dim>& z);
//! print internal data
void print (std::ostream& ss, int indent);
private:
Tupel<int,dim> N; // number of elements per direction
int ne[dim+1]; // number of elements per codimension
......@@ -95,9 +186,6 @@ namespace Dune {
int power2 (int i) {return 1<<i;}
int ones (int b); // count number of bits set in binary rep of b
int partition (Tupel<int,dim>& z); // get binary representation of partition
Tupel<int,dim> reduce (Tupel<int,dim>& z);
Tupel<int,dim> expand (Tupel<int,dim>& r, int b);
};
} // end namespace
......
grid/sgrid/sgrid.cc
+ 596
− 240
View file @ 3b4e8b21
This diff is collapsed.
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