From 3478a00b4569f1ca2ab34407ae33d3f0936afa94 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@dune-project.org>
Date: Mon, 18 Oct 2004 08:52:34 +0000
Subject: [PATCH] FieldVector replaces Vec
[[Imported from SVN: r880]]
---
common/fixedvector.hh | 180 -----
common/function.hh | 2 +-
common/functionspace.hh | 6 +-
common/fvector.hh | 941 ++++++++++++++++++++++
common/matvec.hh | 27 +-
fem/common/basefunctions.hh | 18 +-
fem/common/fastbase.cc | 11 +-
fem/common/fastbase.hh | 8 +-
fem/dgspace/monomialbase.cc | 12 +-
fem/dgspace/monomialbase.hh | 14 +-
fem/discfuncarray.hh | 15 +-
fem/discfuncarray/discfuncarray.cc | 4 +-
fem/fastquad.hh | 18 +-
fem/fastquad/gaussquadimp.cc | 8 +-
fem/fastquad/quadtetratri.hh | 4 +-
fem/feoperator.hh | 2 +
fem/gaussquadrature.hh | 4 +-
fem/gaussquadrature/gaussquadrature.cc | 12 +-
fem/lagrangebase.hh | 22 +-
fem/lagrangebase/lagrangebasefunctions.hh | 64 +-
fem/lagrangebase/lagrangespace.cc | 16 +-
fem/norms/l2norm.hh | 2 +-
fem/operator/discreteoperator.cc | 8 +-
fem/operator/laplace.hh | 6 +-
fem/operator/massmatrix.hh | 5 +-
fem/transfer/mgtransfer.cc | 68 +-
grid/albertgrid.hh | 42 +-
grid/albertgrid/albertgrid.cc | 100 +--
grid/bsgrid.hh | 28 +-
grid/bsgrid/bsgrid.cc | 28 +-
grid/bsgrid/bsinclude.hh | 2 +-
grid/common/grid.cc | 74 +-
grid/common/grid.hh | 58 +-
grid/sgrid.hh | 57 +-
grid/sgrid/numbering.hh | 6 +-
grid/sgrid/sgrid.cc | 61 +-
grid/simplegrid.hh | 54 +-
grid/uggrid/ugfunctions.hh | 19 +-
grid/uggrid/uggridelement.cc | 26 +-
grid/uggrid/uggridelement.hh | 194 ++++-
grid/uggrid/uggridentity.cc | 60 +-
grid/uggrid/uggridentity.hh | 4 +-
grid/uggrid/uggridhieriterator.cc | 8 +
grid/uggrid/uggridleveliterator.cc | 8 +-
grid/uggrid/ugintersectionit.cc | 16 +-
grid/uggrid/ugintersectionit.hh | 10 +-
grid/yaspgrid.hh | 90 ++-
grid/yaspgrid/grids.hh | 30 +-
48 files changed, 1735 insertions(+), 717 deletions(-)
delete mode 100644 common/fixedvector.hh
create mode 100644 common/fvector.hh
diff --git a/common/fixedvector.hh b/common/fixedvector.hh
deleted file mode 100644
index 11cd9f0c7..000000000
--- a/common/fixedvector.hh
+++ /dev/null
@@ -1,180 +0,0 @@
-// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
-// vi: set et ts=4 sw=2 sts=2:
-#ifndef __DUNE_FIXEDVECTOR_HH__
-#define __DUNE_FIXEDVECTOR_HH__
-
-#include <iostream>
-#include <math.h>
-
-#include "misc.hh"
-#include "fixedarray.hh"
-
-namespace Dune {
- /** @defgroup Common Dune Common Module
-
- This module contains classes for general usage in dune, such as e.g.
- (small) dense matrices and vectors or containers.
-
- @{
- */
-
-
- //************************************************************************
- /*!
- Generic vector class for short vectors in d dimensions. Used e.g.
- for global or local coordinates.
- */
- template<int dim, class T = double>
- class Vec : public FixedArray<T, dim> {
- public:
-
- //! Constructor making uninitialized vector
- Vec() {}
-
- //! Constructor making vector from built-in array
- Vec (T* y) {for(int i=0; i<dim; i++) this->a[i]=y[i];}
-
- //! Constructor making vector with one coordinate set, others zeroed
- Vec (int k, T t)
- {
- for (int i=0; i<dim; i++) this->a[i] = 0;
- this->a[k] = t;
- }
-
- //! Constructor making vector with identical coordinates
- Vec (T t)
- {
- for (int i=0; i<dim; i++) this->a[i] = t;
- }
-
- //! operator () for read/write access to element of the vector
- T& operator() (int i) {return this->a[i];}
-
- //! operator () for read access to element of the vector
- const T& operator() (int i) const {return this->a[i];}
-
- //! operator+ adds two vectors
- Vec<dim,T>& operator+= (const Vec<dim,T>& b)
- {
- for (int i=0; i<dim; i++) this->a[i] += b.a[i];
- return *this;
- }
-
- //! Vector addition
- Vec<dim,T> operator+ (const Vec<dim,T>& b) const
- {
- Vec<dim,T> z = *this;
- return (z+=b);
- }
-
- //! operator- binary minus
- Vec<dim,T>& operator-= (const Vec<dim,T>& b)
- {
- for (int i=0; i<dim; i++) this->a[i] -= b.a[i];
- return *this;
- }
-
- //! Vector subtraction
- Vec<dim,T> operator- (const Vec<dim,T>& b) const
- {
- Vec<dim,T> z = *this;
- return (z-=b);
- }
-
- //! scalar product of two vectors with operator*
- T operator* (const Vec<dim,T>& b) const
- {
- T s=0; for (int i=0; i<dim; i++) s += this->a[i]*b.a[i];return s;
- }
-
- //! multiplication of vector with scalar
- Vec<dim,T> operator* (T k) const
- {
- Vec<dim,T> z; for (int i=0; i<dim; i++) z.a[i] = k*this->a[i];return z;
- }
-
- //! multiplication assignment with scalar
- Vec<dim,T>& operator*= (T k)
- {
- for (int i=0; i<dim; i++) this->a[i] *= k;
- return *this;
- }
-
- //! operator ==
- bool operator== (const Vec<dim,T> & b) const
- {
- for (int i=0; i<dim; i++) if (this->a[i] != b.a[i]) return false;
- return true;
- }
-
- //! operator !=
- bool operator!= (const Vec<dim,T> & b) const
- {
- return ! operator== (b);
- }
-
- //! 1 norm
- T norm1 () const
- {
- T s=0.0;
- for (int i=0; i<dim; i++) s += ABS(this->a[i]);
- return s;
- }
-
- //! 2 norm
- T norm2 () const
- {
- T s=0.0;
- for (int i=0; i<dim; i++) s += this->a[i]*this->a[i];
- return sqrt(s);
- }
-
- //! Infinity norm
- T norminfty () const
- {
- T s=0.0;
- for (int i=0; i<dim; i++)
- if (ABS(this->a[i])>s) s = ABS(this->a[i]);
- return s;
- }
-
- //! Euclidean distance of two vectors
- T distance (const Vec<dim,T>& b) const
- {
- T s=0.0;
- for (int i=0; i<dim; i++) s += (this->a[i]-b.a[i])*(this->a[i]-b.a[i]);
- return sqrt(s);
- }
- };
-
- //! multiplication of scalar with vector
- template<int dim, class T>
- inline Vec<dim,T> operator* (T k, Vec<dim,T> b)
- {
- Vec<dim,T> z; for (int i=0; i<dim; i++) z(i) = k*b(i);return z;
- }
-
- //! unary minus
- template<int dim, class T>
- inline Vec<dim,T> operator- (Vec<dim,T> b)
- {
- Vec<dim,T> z; for (int i=0; i<dim; i++) z(i) = -b(i);return z;
- }
-
- // Ausgabe
- template <int dim, class T>
- std::ostream& operator<< (std::ostream& s, const Vec<dim,T>& a)
- {
- s << "Vec<" << dim << "> = [";
- for (int i=0; i<dim; i++)
- s << " " << a(i);
- s << " ]";
- return s;
- }
-
- /** @} */
-
-} // end namespace Dune
-
-
-#endif
diff --git a/common/function.hh b/common/function.hh
index ca2bef8b6..f101a42d0 100644
--- a/common/function.hh
+++ b/common/function.hh
@@ -51,7 +51,7 @@ namespace Dune {
template <int derivation>
//! ???
- void evaluate ( const Vec<derivation,deriType> &diffVariable,
+ void evaluate ( const FieldVector<deriType, derivation> &diffVariable,
const Domain & , Range &) const {};
//! ???
diff --git a/common/functionspace.hh b/common/functionspace.hh
index 3ea73493a..bbe2f20c9 100644
--- a/common/functionspace.hh
+++ b/common/functionspace.hh
@@ -28,12 +28,12 @@ namespace Dune {
/** \todo Please doc me! */
typedef Mat < n, m, RangeField> JacobianRange;
/** \todo Please doc me! */
- typedef Vec < m, Mat< n, n, RangeField> > HessianRange ;
+ typedef FieldVector<Mat< n, n, RangeField>, m> HessianRange ;
/** \todo Please doc me! */
- typedef Vec<n, DomainField> Domain;
+ typedef FieldVector<DomainField, n> Domain;
/** \todo Please doc me! */
- typedef Vec<m, RangeField> Range;
+ typedef FieldVector<RangeField, m> Range;
/** \todo Please doc me! */
enum { DimDomain = n, DimRange = m};
diff --git a/common/fvector.hh b/common/fvector.hh
new file mode 100644
index 000000000..41c9d383c
--- /dev/null
+++ b/common/fvector.hh
@@ -0,0 +1,941 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef __DUNE_FVECTOR_HH__
+#define __DUNE_FVECTOR_HH__
+
+#include <math.h>
+#include <complex>
+
+#ifdef DUNE_ISTL_WITH_CHECKING
+#include "istlexception.hh"
+#endif
+
+/*! \file __FILE__
+
+ This file implements a vector constructed from a given type
+ representing a field and a compile-time given size.
+ */
+
+namespace Dune {
+
+ /** @defgroup ISTL Iterative Solvers Template Library
+ @addtogroup ISTL
+ @{
+ */
+
+ // forward declaration of template
+ template<class K, int n> class FieldVector;
+
+ // template meta program for assignment from scalar
+ template<int I>
+ struct fvmeta_assignscalar {
+ template<class K, int n>
+ static K& assignscalar (FieldVector<K,n>& x, const K& k)
+ {
+ x[I] = fvmeta_assignscalar<I-1>::assignscalar(x,k);
+ return x[I];
+ }
+ };
+ template<>
+ struct fvmeta_assignscalar<0> {
+ template<class K, int n>
+ static K& assignscalar (FieldVector<K,n>& x, const K& k)
+ {
+ x[0] = k;
+ return x[0];
+ }
+ };
+
+ // template meta program for operator+=
+ template<int I>
+ struct fvmeta_plusequal {
+ template<class K, int n>
+ static void plusequal (FieldVector<K,n>& x, const FieldVector<K,n>& y)
+ {
+ fvmeta_plusequal<I-1>::plusequal(x,y);
+ x[I] += y[I];
+ }
+ template<class K, int n>
+ static void plusequal (FieldVector<K,n>& x, const K& y)
+ {
+ fvmeta_plusequal<I-1>::plusequal(x,y);
+ x[I] += y;
+ }
+ };
+ template<>
+ struct fvmeta_plusequal<0> {
+ template<class K, int n>
+ static void plusequal (FieldVector<K,n>& x, const FieldVector<K,n>& y)
+ {
+ x[0] += y[0];
+ }
+ template<class K, int n>
+ static void plusequal (FieldVector<K,n>& x, const K& y)
+ {
+ x[0] += y;
+ }
+ };
+
+ // template meta program for operator-=
+ template<int I>
+ struct fvmeta_minusequal {
+ template<class K, int n>
+ static void minusequal (FieldVector<K,n>& x, const FieldVector<K,n>& y)
+ {
+ fvmeta_minusequal<I-1>::minusequal(x,y);
+ x[I] -= y[I];
+ }
+ template<class K, int n>
+ static void minusequal (FieldVector<K,n>& x, const K& y)
+ {
+ fvmeta_minusequal<I-1>::minusequal(x,y);
+ x[I] -= y;
+ }
+ };
+ template<>
+ struct fvmeta_minusequal<0> {
+ template<class K, int n>
+ static void minusequal (FieldVector<K,n>& x, const FieldVector<K,n>& y)
+ {
+ x[0] -= y[0];
+ }
+ template<class K, int n>
+ static void minusequal (FieldVector<K,n>& x, const K& y)
+ {
+ x[0] -= y;
+ }
+ };
+
+ // template meta program for operator*=
+ template<int I>
+ struct fvmeta_multequal {
+ template<class K, int n>
+ static void multequal (FieldVector<K,n>& x, const K& k)
+ {
+ fvmeta_multequal<I-1>::multequal(x,k);
+ x[I] *= k;
+ }
+ };
+ template<>
+ struct fvmeta_multequal<0> {
+ template<class K, int n>
+ static void multequal (FieldVector<K,n>& x, const K& k)
+ {
+ x[0] *= k;
+ }
+ };
+
+ // template meta program for operator/=
+ template<int I>
+ struct fvmeta_divequal {
+ template<class K, int n>
+ static void divequal (FieldVector<K,n>& x, const K& k)
+ {
+ fvmeta_divequal<I-1>::divequal(x,k);
+ x[I] /= k;
+ }
+ };
+ template<>
+ struct fvmeta_divequal<0> {
+ template<class K, int n>
+ static void divequal (FieldVector<K,n>& x, const K& k)
+ {
+ x[0] /= k;
+ }
+ };
+
+ // template meta program for axpy
+ template<int I>
+ struct fvmeta_axpy {
+ template<class K, int n>
+ static void axpy (FieldVector<K,n>& x, const K& a, const FieldVector<K,n>& y)
+ {
+ fvmeta_axpy<I-1>::axpy(x,a,y);
+ x[I] += a*y[I];
+ }
+ };
+ template<>
+ struct fvmeta_axpy<0> {
+ template<class K, int n>
+ static void axpy (FieldVector<K,n>& x, const K& a, const FieldVector<K,n>& y)
+ {
+ x[0] += a*y[0];
+ }
+ };
+
+ // template meta program for dot
+ template<int I>
+ struct fvmeta_dot {
+ template<class K, int n>
+ static K dot (const FieldVector<K,n>& x, const FieldVector<K,n>& y)
+ {
+ return x[I]*y[I] + fvmeta_dot<I-1>::dot(x,y);
+ }
+ };
+ template<>
+ struct fvmeta_dot<0> {
+ template<class K, int n>
+ static K dot (const FieldVector<K,n>& x, const FieldVector<K,n>& y)
+ {
+ return x[0]*y[0];
+ }
+ };
+
+ // some abs functions needed for the norms
+ template<class K>
+ inline double fvmeta_abs (const K& k)
+ {
+ return (k >= 0) ? k : -k;
+ }
+
+ template<class K>
+ inline double fvmeta_abs (const std::complex<K>& c)
+ {
+ return sqrt(c.real()*c.real() + c.imag()*c.imag());
+ }
+
+ template<class K>
+ inline double fvmeta_absreal (const K& k)
+ {
+ return fvmeta_abs(k);
+ }
+
+ template<class K>
+ inline double fvmeta_absreal (const std::complex<K>& c)
+ {
+ return fvmeta_abs(c.real()) + fvmeta_abs(c.imag());
+ }
+
+ template<class K>
+ inline double fvmeta_abs2 (const K& k)
+ {
+ return k*k;
+ }
+
+ template<class K>
+ inline double fvmeta_abs2 (const std::complex<K>& c)
+ {
+ return c.real()*c.real() + c.imag()*c.imag();
+ }
+
+ // template meta program for one_norm
+ template<int I>
+ struct fvmeta_one_norm {
+ template<class K, int n>
+ static double one_norm (const FieldVector<K,n>& x)
+ {
+ return fvmeta_abs(x[I]) + fvmeta_one_norm<I-1>::one_norm(x);
+ }
+ };
+ template<>
+ struct fvmeta_one_norm<0> {
+ template<class K, int n>
+ static double one_norm (const FieldVector<K,n>& x)
+ {
+ return fvmeta_abs(x[0]);
+ }
+ };
+
+ // template meta program for one_norm_real
+ template<int I>
+ struct fvmeta_one_norm_real {
+ template<class K, int n>
+ static double one_norm_real (const FieldVector<K,n>& x)
+ {
+ return fvmeta_absreal(x[I]) + fvmeta_one_norm_real<I-1>::one_norm_real(x);
+ }
+ };
+ template<>
+ struct fvmeta_one_norm_real<0> {
+ template<class K, int n>
+ static double one_norm_real (const FieldVector<K,n>& x)
+ {
+ return fvmeta_absreal(x[0]);
+ }
+ };
+
+ // template meta program for two_norm squared
+ template<int I>
+ struct fvmeta_two_norm2 {
+ template<class K, int n>
+ static double two_norm2 (const FieldVector<K,n>& x)
+ {
+ return fvmeta_abs2(x[I]) + fvmeta_two_norm2<I-1>::two_norm2(x);
+ }
+ };
+ template<>
+ struct fvmeta_two_norm2<0> {
+ template<class K, int n>
+ static double two_norm2 (const FieldVector<K,n>& x)
+ {
+ return fvmeta_abs2(x[0]);
+ }
+ };
+
+ // template meta program for infinity norm
+ template<int I>
+ struct fvmeta_infinity_norm {
+ template<class K, int n>
+ static double infinity_norm (const FieldVector<K,n>& x)
+ {
+ return std::max(fvmeta_abs(x[I]),fvmeta_infinity_norm<I-1>::infinity_norm(x));
+ }
+ };
+ template<>
+ struct fvmeta_infinity_norm<0> {
+ template<class K, int n>
+ static double infinity_norm (const FieldVector<K,n>& x)
+ {
+ return fvmeta_abs(x[0]);
+ }
+ };
+
+ // template meta program for simplified infinity norm
+ template<int I>
+ struct fvmeta_infinity_norm_real {
+ template<class K, int n>
+ static double infinity_norm_real (const FieldVector<K,n>& x)
+ {
+ return std::max(fvmeta_absreal(x[I]),fvmeta_infinity_norm_real<I-1>::infinity_norm_real(x));
+ }
+ };
+ template<>
+ struct fvmeta_infinity_norm_real<0> {
+ template<class K, int n>
+ static double infinity_norm_real (const FieldVector<K,n>& x)
+ {
+ return fvmeta_absreal(x[0]);
+ }
+ };
+
+ /**! Construct a vector space out of a tensor product of fields.
+ K is the field type (use float, double, complex, etc) and n
+ is the number of components.
+
+ It is generally assumed that K is a numerical type compatible with double
+ (E.g. norms are always computed in double precision).
+
+ Implementation of all members uses template meta programs where appropriate
+ */
+ template<class K, int SIZE>
+ class FieldVector
+ {
+ //! The actual number of elements that gets allocated.
+ //! It's always at least 1.
+ /** \todo Wie genau gehen wir mit dem Fall N==0 um? So? */
+ enum { n = (SIZE > 0) ? SIZE : 1 };
+
+ public:
+ // standard constructor and everything is sufficient ...
+
+
+ //===== type definitions and constants
+
+ //! export the type representing the field
+ typedef K field_type;
+
+ //! export the type representing the components
+ typedef K block_type;
+
+ //! We are at the leaf of the block recursion
+ enum {blocklevel = 1};
+
+ //! export size
+ enum {size = n};
+
+ //! Constructor making uninitialized vector
+ FieldVector() {}
+
+ //! Constructor making vector with identical coordinates
+ explicit FieldVector (const K& t)
+ {
+ for (int i=0; i<n; i++) p[i] = t;
+ }
+
+ //===== assignment from scalar
+ FieldVector& operator= (const K& k)
+ {
+ fvmeta_assignscalar<n-1>::assignscalar(*this,k);
+ return *this;
+ }
+
+
+ //===== access to components
+
+ //! random access
+ K& operator[] (int i)
+ {
+#ifdef DUNE_ISTL_WITH_CHECKING
+ if (i<0 || i>=n) DUNE_THROW(ISTLError,"index out of range");
+#endif
+ return p[i];
+ }
+
+ //! same for read only access
+ const K& operator[] (int i) const
+ {
+#ifdef DUNE_ISTL_WITH_CHECKING
+ if (i<0 || i>=n) DUNE_THROW(ISTLError,"index out of range");
+#endif
+ return p[i];
+ }
+
+ // forward declaration
+ class ConstIterator;
+
+ //! Iterator class for sequential access
+ class Iterator
+ {
+ public:
+ //! constructor
+ Iterator (K* _p, int _i)
+ {
+ p = _p;
+ i = _i;
+ }
+
+ //! empty constructor, use with care
+ Iterator ()
+ { }
+
+ //! prefix increment
+ Iterator& operator++()
+ {
+ ++i;
+ return *this;
+ }
+
+ //! prefix decrement
+ Iterator& operator--()
+ {
+ --i;
+ return *this;
+ }
+
+ //! equality
+ bool operator== (const Iterator& it) const
+ {
+ return (p+i)==(it.p+it.i);
+ }
+
+ //! inequality
+ bool operator!= (const Iterator& it) const
+ {
+ return (p+i)!=(it.p+it.i);
+ }
+
+ //! dereferencing
+ K& operator* ()
+ {
+ return p[i];
+ }
+
+ //! arrow
+ K* operator-> ()
+ {
+ return p+i;
+ }
+
+ //! return index
+ int index ()
+ {
+ return i;
+ }
+
+ friend class ConstIterator;
+
+ private:
+ K* p;
+ int i;
+ };
+
+ //! begin iterator
+ Iterator begin ()
+ {
+ return Iterator(p,0);
+ }
+
+ //! end iterator
+ Iterator end ()
+ {
+ return Iterator(p,n);
+ }
+
+ //! begin iterator
+ Iterator rbegin ()
+ {
+ return Iterator(p,n-1);
+ }
+
+ //! end iterator
+ Iterator rend ()
+ {
+ return Iterator(p,-1);
+ }
+
+ //! return iterator to given element or end()
+ Iterator find (int i)
+ {
+ if (i>=0 && i<n)
+ return Iterator(p,i);
+ else
+ return Iterator(p,n);
+ }
+
+ //! ConstIterator class for sequential access
+ class ConstIterator
+ {
+ public:
+ //! constructor
+ ConstIterator ()
+ {
+ p = 0;
+ i = 0;
+ }
+
+ //! constructor from pointer
+ ConstIterator (const K* _p, int _i) : p(_p), i(_i)
+ { }
+
+ //! constructor from non_const iterator
+ ConstIterator (const Iterator& it) : p(it.p), i(it.i)
+ { }
+
+ //! prefix increment
+ ConstIterator& operator++()
+ {
+ ++i;
+ return *this;
+ }
+
+ //! prefix decrement
+ ConstIterator& operator--()
+ {
+ --i;
+ return *this;
+ }
+
+ //! equality
+ bool operator== (const ConstIterator& it) const
+ {
+ return (p+i)==(it.p+it.i);
+ }
+
+ //! inequality
+ bool operator!= (const ConstIterator& it) const
+ {
+ return (p+i)!=(it.p+it.i);
+ }
+
+ //! equality
+ bool operator== (const Iterator& it) const
+ {
+ return (p+i)==(it.p+it.i);
+ }
+
+ //! inequality
+ bool operator!= (const Iterator& it) const
+ {
+ return (p+i)!=(it.p+it.i);
+ }
+
+ //! dereferencing
+ const K& operator* () const
+ {
+ return p[i];
+ }
+
+ //! arrow
+ const K* operator-> () const
+ {
+ return p+i;
+ }
+
+ // return index corresponding to pointer
+ int index () const
+ {
+ return i;
+ }
+
+ friend class Iterator;
+
+ private:
+ const K* p;
+ int i;
+ };
+
+ //! begin ConstIterator
+ ConstIterator begin () const
+ {
+ return ConstIterator(p,0);
+ }
+
+ //! end ConstIterator
+ ConstIterator end () const
+ {
+ return ConstIterator(p,n);
+ }
+
+ //! begin ConstIterator
+ ConstIterator rbegin () const
+ {
+ return ConstIterator(p,n-1);
+ }
+
+ //! end ConstIterator
+ ConstIterator rend () const
+ {
+ return ConstIterator(p,-1);
+ }
+
+ //! return iterator to given element or end()
+ ConstIterator find (int i) const
+ {
+ if (i>=0 && i<n)
+ return ConstIterator(p,i);
+ else
+ return ConstIterator(p,n);
+ }
+
+ //===== vector space arithmetic
+
+ //! vector space addition
+ FieldVector& operator+= (const FieldVector& y)
+ {
+ fvmeta_plusequal<n-1>::plusequal(*this,y);
+ return *this;
+ }
+
+ //! vector space subtraction
+ FieldVector& operator-= (const FieldVector& y)
+ {
+ fvmeta_minusequal<n-1>::minusequal(*this,y);
+ return *this;
+ }
+
+ //! Binary vector addition
+ FieldVector<K, size> operator+ (const FieldVector<K, size>& b) const
+ {
+ FieldVector<K, size> z = *this;
+ return (z+=b);
+ }
+
+ //! Binary vector subtraction
+ FieldVector<K, size> operator- (const FieldVector<K, size>& b) const
+ {
+ FieldVector<K, size> z = *this;
+ return (z-=b);
+ }
+
+ //! vector space add scalar to all comps
+ FieldVector& operator+= (const K& k)
+ {
+ fvmeta_plusequal<n-1>::plusequal(*this,k);
+ return *this;
+ }
+
+ //! vector space subtract scalar from all comps
+ FieldVector& operator-= (const K& k)
+ {
+ fvmeta_minusequal<n-1>::minusequal(*this,k);
+ return *this;
+ }
+
+ //! vector space multiplication with scalar
+ FieldVector& operator*= (const K& k)
+ {
+ fvmeta_multequal<n-1>::multequal(*this,k);
+ return *this;
+ }
+
+ //! vector space division by scalar
+ FieldVector& operator/= (const K& k)
+ {
+ fvmeta_divequal<n-1>::divequal(*this,k);
+ return *this;
+ }
+
+ //! vector space axpy operation
+ FieldVector& axpy (const K& a, const FieldVector& y)
+ {
+ fvmeta_axpy<n-1>::axpy(*this,a,y);
+ return *this;
+ }
+
+
+ //===== Euclidean scalar product
+
+ //! scalar product
+ const K operator* (const FieldVector& y) const
+ {
+ return fvmeta_dot<n-1>::dot(*this,y);
+ }
+
+
+ //===== norms
+
+ //! one norm (sum over absolute values of entries)
+ double one_norm () const
+ {
+ return fvmeta_one_norm<n-1>::one_norm(*this);
+ }
+
+ //! simplified one norm (uses Manhattan norm for complex values)
+ double one_norm_real () const
+ {
+ return fvmeta_one_norm_real<n-1>::one_norm_real(*this);
+ }
+
+ //! two norm sqrt(sum over squared values of entries)
+ double two_norm () const
+ {
+ return sqrt(fvmeta_two_norm2<n-1>::two_norm2(*this));
+ }
+
+ //! sqare of two norm (sum over squared values of entries), need for block recursion
+ double two_norm2 () const
+ {
+ return fvmeta_two_norm2<n-1>::two_norm2(*this);
+ }
+
+ //! infinity norm (maximum of absolute values of entries)
+ double infinity_norm () const
+ {
+ return fvmeta_infinity_norm<n-1>::infinity_norm(*this);
+ }
+
+ //! simplified infinity norm (uses Manhattan norm for complex values)
+ double infinity_norm_real () const
+ {
+ return fvmeta_infinity_norm_real<n-1>::infinity_norm_real(*this);
+ }
+
+
+ //===== sizes
+
+ //! number of blocks in the vector (are of size 1 here)
+ int N () const
+ {
+ return n;
+ }
+
+ //! dimension of the vector space
+ int dim () const
+ {
+ return n;
+ }
+
+ //===== conversion operator
+
+ operator K () {return *p;}
+
+ void print (std::ostream& s) const
+ {
+ for (int i=0; i<n; i++)
+ if (i>0)
+ s << " " << p[i];
+ else
+ s << p[i];
+ }
+
+ private:
+ // the data, very simply a built in array
+ K p[n];
+ };
+
+
+ // Ausgabe
+ template<class K, int n>
+ std::ostream& operator<< (std::ostream& s, const FieldVector<K,n>& v)
+ {
+ v.print(s);
+ return s;
+ }
+
+
+ // forward declarations
+ template<class K> class K1Vector;
+ template<class K> class K11Matrix;
+
+ /**! Vectors containing only one component
+ */
+ template<class K>
+ class K1Vector
+ {
+ public:
+ friend class K11Matrix<K>;
+
+ //===== type definitions and constants
+
+ //! export the type representing the field
+ typedef K field_type;
+
+ //! export the type representing the components
+ typedef K block_type;
+
+ //! We are at the leaf of the block recursion
+ enum {blocklevel = 1};
+
+ //! export size
+ enum {size = 1};
+
+ //===== construction
+
+ K1Vector ()
+ { }
+
+ K1Vector (const K& k)
+ {
+ p = k;
+ }
+
+ //===== assignment from basetype
+ K1Vector& operator= (const K& k)
+ {
+ p = k;
+ return *this;
+ }
+
+ //===== access to data
+
+ //! random access
+ K& operator() (void)
+ {
+ return p;
+ }
+
+ //! same for read only access
+ const K& operator() (void) const
+ {
+ return p;
+ }
+
+ //===== vector space arithmetic
+
+ //! vector space addition
+ K1Vector& operator+= (const K1Vector& y)
+ {
+ p += y.p;
+ return *this;
+ }
+
+ //! vector space subtraction
+ K1Vector& operator-= (const K1Vector& y)
+ {
+ p -= y.p;
+ return *this;
+ }
+
+ //! vector space add scalar to each comp
+ K1Vector& operator+= (const K& k)
+ {
+ p += k;
+ return *this;
+ }
+
+ //! vector space subtract scalar from each comp
+ K1Vector& operator-= (const K& k)
+ {
+ p -= k;
+ return *this;
+ }
+
+ //! vector space multiplication with scalar
+ K1Vector& operator*= (const K& k)
+ {
+ p *= k;
+ return *this;
+ }
+
+ //! vector space division by scalar
+ K1Vector& operator/= (const K& k)
+ {
+ p /= k;
+ return *this;
+ }
+
+ //! vector space axpy operation
+ K1Vector& axpy (const K& a, const K1Vector& y)
+ {
+ p += a*y.p;
+ return *this;
+ }
+
+
+ //===== Euclidean scalar product
+
+ //! scalar product
+ const K operator* (const K1Vector& y) const
+ {
+ return p*y.p;
+ }
+
+
+ //===== norms
+
+ //! one norm (sum over absolute values of entries)
+ double one_norm () const
+ {
+ return fvmeta_abs(p);
+ }
+
+ //! simplified one norm (uses Manhattan norm for complex values)
+ double one_norm_real () const
+ {
+ return fvmeta_abs_real(p);
+ }
+
+ //! two norm sqrt(sum over squared values of entries)
+ double two_norm () const
+ {
+ return sqrt(fvmeta_abs2(p));
+ }
+
+ //! sqare of two norm (sum over squared values of entries), need for block recursion
+ double two_norm2 () const
+ {
+ return fvmeta_abs2(p);
+ }
+
+ //! infinity norm (maximum of absolute values of entries)
+ double infinity_norm () const
+ {
+ return fvmeta_abs(p);
+ }
+
+ //! simplified infinity norm (uses Manhattan norm for complex values)
+ double infinity_norm_real () const
+ {
+ return fvmeta_abs_real(p);
+ }
+
+
+ //===== sizes
+
+ //! number of blocks in the vector (are of size 1 here)
+ int N () const
+ {
+ return 1;
+ }
+
+ //! dimension of the vector space
+ int dim () const
+ {
+ return 1;
+ }
+
+ //===== conversion operator
+
+ operator K () {return p;}
+
+ private:
+ // the data
+ K p;
+ };
+
+
+ /** @} end documentation */
+
+} // end namespace
+
+#endif
diff --git a/common/matvec.hh b/common/matvec.hh
index 91bd7f758..e740b1ac4 100644
--- a/common/matvec.hh
+++ b/common/matvec.hh
@@ -7,7 +7,7 @@
#include <math.h>
#include "misc.hh"
-#include "fixedvector.hh"
+#include "fvector.hh"
#include <dune/common/exceptions.hh>
namespace Dune {
@@ -43,23 +43,24 @@ namespace Dune {
}
//! operator () for read/write access to element in matrix
- T& operator() (int i, int j) {return a[j](i);}
+ T& operator() (int i, int j) {return a[j][i];}
//! operator () for read/write access to element in matrix
- const T& operator() (int i, int j) const {return a[j](i);}
+ const T& operator() (int i, int j) const {return a[j][i];}
//! operator () for read/write access to column vector
- Vec<n,T>& operator() (int j) {return a[j];}
+ FieldVector<T,n>& operator() (int j) {return a[j];}
//! operator [] for read/write access to column vector
- Vec<n,T>& operator[] (int j) {return a[j];}
+ FieldVector<T,n>& operator[] (int j) {return a[j];}
//! matrix/vector multiplication
- Vec<n,T> operator* (const Vec<m,T>& x)
+ FieldVector<T,n> operator* (const FieldVector<T,m>& x)
{
- Vec<n,T> z(0.0);
+ FieldVector<T,n> z(0.0);
for (int j=0; j<m; j++)
- for (int i=0; i<n; i++) z(i) += a[j](i) * x[j];
+ for (int i=0; i<n; i++)
+ z[i] += a[j][i] * x[j];
return z;
}
@@ -73,12 +74,12 @@ namespace Dune {
}
//! matrix/vector multiplication with transpose of matrix
- Vec<m,T> mult_t (const Vec<n,T>& x) const
+ FieldVector<T,m> mult_t (const FieldVector<T,n>& x) const
{
- Vec<m,T> z(0.0);
+ FieldVector<T,m> z(0.0);
const Mat<n,m,T> &matrix = (*this);
for (int j=0; j<n; j++)
- for (int i=0; i<m; i++) z(i) += matrix(j,i) * x[j];
+ for (int i=0; i<m; i++) z[i] += matrix(j,i) * x[j];
return z;
}
@@ -126,7 +127,7 @@ namespace Dune {
}
//! scalar product of two vectors; one of them stored in matrixform
- T operator* (const Vec<n,T>& b) const
+ T operator* (const FieldVector<T,n>& b) const
{
if (m != 1)
DUNE_THROW(MathError,
@@ -138,7 +139,7 @@ namespace Dune {
private:
//! built-in array to hold the data
- Vec<n,T> a[m];
+ FieldVector<T,n> a[m];
};
namespace HelpMat {
diff --git a/fem/common/basefunctions.hh b/fem/common/basefunctions.hh
index dbc0607ff..fe416f0a7 100644
--- a/fem/common/basefunctions.hh
+++ b/fem/common/basefunctions.hh
@@ -31,7 +31,7 @@ namespace Dune {
template <int dim>
struct DiffVariable
{
- typedef Vec<dim, deriType> Type;
+ typedef FieldVector<deriType, dim> Type;
};
//*************************************************************************
@@ -64,16 +64,16 @@ namespace Dune {
//! methods with template parameter "length of Vec".
//! Though the evaluate Methods can be spezialized for each
//! differentiation order
- virtual void evaluate ( const Vec<0,deriType> &diffVariable,
+ virtual void evaluate ( const FieldVector<deriType, 0> &diffVariable,
const Domain & , Range &) const {}; // = 0 ?
//! diffVariable contain the component of the gradient which is delivered.
//! for example gradient of the basefunction x component ==>
//! diffVariable(0) == 0, y component ==> diffVariable(0) == 1 ...
- virtual void evaluate ( const Vec<1,deriType> &diffVariable,
+ virtual void evaluate ( const FieldVector<deriType, 1> &diffVariable,
const Domain & , Range &) const {}; // = 0 ?
- virtual void evaluate ( const Vec<2,deriType> &diffVariable,
+ virtual void evaluate ( const FieldVector<deriType, 2> &diffVariable,
const Domain & , Range &) const {}; // = 0 ?
private:
@@ -146,14 +146,14 @@ namespace Dune {
};
template <int diffOrd>
- void evaluate ( int baseFunct, const Vec<diffOrd,deriType> &diffVariable, const
+ void evaluate ( int baseFunct, const FieldVector<deriType, diffOrd> &diffVariable, const
Domain & x, Range & phi ) const {
std::cout << "BaseFunctionSetInterface::evaluate \n";
asImp().evaluate( baseFunct, diffVariable, x, phi );
}
template <int diffOrd, class QuadratureType >
- void evaluate ( int baseFunct, const Vec<diffOrd,deriType> &diffVariable, QuadratureType & quad, int quadPoint, Range & phi ) const {
+ void evaluate ( int baseFunct, const FieldVector<deriType, diffOrd> &diffVariable, QuadratureType & quad, int quadPoint, Range & phi ) const {
asImp().evaluate( baseFunct, diffVariable, quad, quadPoint, phi );
}
protected:
@@ -241,19 +241,19 @@ namespace Dune {
{
asImp().evaluate( baseFunct, jacobianDiffVar_[i] , quad, quadPoint, tmp_ );
for(int j=0; j<dimRow; j++)
- phi(i,j) = tmp_(j);
+ phi(i,j) = tmp_[j];
}
return;
}
private:
//! just diffVariable for evaluation of the functions
- const Vec<0,deriType> diffVariable_;
+ const FieldVector<deriType, 0> diffVariable_;
//! temporary Range vec
mutable Range tmp_;
- Vec<1,deriType> jacobianDiffVar_[dimCol];
+ FieldVector<deriType, 1> jacobianDiffVar_[dimCol];
//! Barton-Nackman trick
BaseFunctionSetImp &asImp() { return static_cast<BaseFunctionSetImp&>(*this); }
diff --git a/fem/common/fastbase.cc b/fem/common/fastbase.cc
index 010cde264..51b93e2f6 100644
--- a/fem/common/fastbase.cc
+++ b/fem/common/fastbase.cc
@@ -21,14 +21,14 @@ FastBaseFunctionSet( FunctionSpaceType & fuspace, int numOfBaseFct )
template <class FunctionSpaceType> template <int diffOrd>
void FastBaseFunctionSet<FunctionSpaceType >::
-evaluate( int baseFunct, const Vec<diffOrd, deriType> &diffVariable, const Domain & x, Range & phi ) const
+evaluate( int baseFunct, const FieldVector<deriType, diffOrd> &diffVariable, const Domain & x, Range & phi ) const
{
getBaseFunction( baseFunct ).evaluate( diffVariable, x, phi );
}
template <class FunctionSpaceType> template <int diffOrd, class QuadratureType>
void FastBaseFunctionSet<FunctionSpaceType >::
-evaluate( int baseFunct, const Vec<diffOrd, deriType> &diffVariable, QuadratureType & quad,
+evaluate( int baseFunct, const FieldVector<deriType, diffOrd> &diffVariable, QuadratureType & quad,
int quadPoint, Range & phi ) const
{
// check if cache is for the used quadrature
@@ -68,13 +68,14 @@ evaluateInit( const QuadratureType &quad )
// go through all possible combinations of 0,1,2
if(diffOrd == 2)
{
- diffVariable(0) = i % DimDomain;
- diffVariable(1) = diffCount;
+ diffVariable[0] = i % DimDomain;
+ diffVariable[1] = diffCount;
diffCount++;
if (diffCount >= DimDomain) diffCount = 0;
}
else
- for(int j=0; j<diffOrd; j++) diffVariable(j) = i % DimDomain;
+ for(int j=0; j<diffOrd; j++)
+ diffVariable[j] = i % DimDomain;
for ( int baseFunc = 0; baseFunc < this->getNumberOfBaseFunctions(); baseFunc++ )
{
diff --git a/fem/common/fastbase.hh b/fem/common/fastbase.hh
index eef0bf23c..6b3506301 100644
--- a/fem/common/fastbase.hh
+++ b/fem/common/fastbase.hh
@@ -69,7 +69,7 @@ namespace Dune {
//! evaluate base function baseFunct with the given diffVariable and a
//! point x and range phi
template <int diffOrd>
- void evaluate ( int baseFunct, const Vec<diffOrd, deriType> &diffVariable,
+ void evaluate ( int baseFunct, const FieldVector<deriType, diffOrd> &diffVariable,
const Domain & x, Range & phi ) const;
//! evaluate base fucntion baseFunct at a given quadrature point
@@ -77,7 +77,7 @@ namespace Dune {
//! qaudrature has changed an the values at the quadrature have to be
//! calulated again
template <int diffOrd, class QuadratureType>
- void evaluate ( int baseFunct, const Vec<diffOrd, deriType> &diffVariable,
+ void evaluate ( int baseFunct, const FieldVector<deriType, diffOrd> &diffVariable,
QuadratureType & quad, int quadPoint, Range & phi ) const;
//! get a reference of the base function baseFunct
@@ -114,12 +114,12 @@ namespace Dune {
//! method to navigate through the vector vecEvaluate, which holds
//! precalculated values
template <int diffOrd>
- int index( int baseFunct, const Vec<diffOrd, deriType> &diffVariable,
+ int index( int baseFunct, const FieldVector<deriType, diffOrd> &diffVariable,
int quadPt, int numQuadPoints ) const
{
int n = 0;
for ( int i = 0; i < diffOrd; i++ )
- n = diffVariable(i) + i * DimDomain;
+ n = diffVariable[i] + i * DimDomain;
return numQuadPoints*(getNumberOfBaseFunctions()*n + baseFunct) + quadPt;
};
diff --git a/fem/dgspace/monomialbase.cc b/fem/dgspace/monomialbase.cc
index 44fc7be3b..cc29ff7dc 100644
--- a/fem/dgspace/monomialbase.cc
+++ b/fem/dgspace/monomialbase.cc
@@ -38,7 +38,7 @@ MonomialBaseFunctionSet( FunctionSpaceType & fuspace, int polOrder )
template <class FunctionSpaceType>
void MonomialBaseFunctionSet<FunctionSpaceType >::
-real_evaluate( int baseFunct, const Vec<0, deriType> &diffVariable,
+real_evaluate( int baseFunct, const FieldVector<deriType, 0> &diffVariable,
const Domain & x, Range & phi ) const
{
phi = power(x[0], Phi_[baseFunct][0]) * power(x[1], Phi_[baseFunct][1]);
@@ -46,7 +46,7 @@ real_evaluate( int baseFunct, const Vec<0, deriType> &diffVariable,
template <class FunctionSpaceType>
void MonomialBaseFunctionSet<FunctionSpaceType >::
-real_evaluate( int baseFunct, const Vec<1, deriType> &diffVariable,
+real_evaluate( int baseFunct, const FieldVector<deriType, 1> &diffVariable,
const Domain & x, Range & phi ) const
{
if (diffVariable(0) == 0) {
@@ -61,7 +61,7 @@ real_evaluate( int baseFunct, const Vec<1, deriType> &diffVariable,
template <class FunctionSpaceType>
void MonomialBaseFunctionSet<FunctionSpaceType >::
-real_evaluate( int baseFunct, const Vec<2, deriType> &diffVariable,
+real_evaluate( int baseFunct, const FieldVector<deriType, 2> &diffVariable,
const Domain & x, Range & phi ) const
{
if (diffVariable(0) == 0) {
@@ -94,7 +94,7 @@ real_evaluate( int baseFunct, const Vec<2, deriType> &diffVariable,
template <class FunctionSpaceType> template <int diffOrd>
void MonomialBaseFunctionSet<FunctionSpaceType >::
-evaluate( int baseFunct, const Vec<diffOrd, deriType> &diffVariable,
+evaluate( int baseFunct, const FieldVector<deriType, diffOrd> &diffVariable,
const Domain & x, Range & phi ) const
{
assert(baseFunct < numOfBaseFct_);
@@ -103,7 +103,7 @@ evaluate( int baseFunct, const Vec<diffOrd, deriType> &diffVariable,
template <class FunctionSpaceType> template <int diffOrd, class QuadratureType>
void MonomialBaseFunctionSet<FunctionSpaceType >::
-evaluate( int baseFunct, const Vec<diffOrd, deriType> &diffVariable,
+evaluate( int baseFunct, const FieldVector<deriType, diffOrd> &diffVariable,
QuadratureType & quad, int quadPoint, Range & phi ) const
{
Domain x = quad.point(quadPoint);
@@ -112,7 +112,7 @@ evaluate( int baseFunct, const Vec<diffOrd, deriType> &diffVariable,
template <class FunctionSpaceType>
void MonomialBaseFunctionSet<FunctionSpaceType >::
-print (std::ostream& s, const Dune::Vec<2,int> & pol) const
+print (std::ostream& s, const Dune::FieldVector<int, 2> & pol) const
{
if (pol[0] > 0 && pol[1] > 0) {
s << "x^" << pol[0]
diff --git a/fem/dgspace/monomialbase.hh b/fem/dgspace/monomialbase.hh
index ea31d34c0..bd8a0bf59 100644
--- a/fem/dgspace/monomialbase.hh
+++ b/fem/dgspace/monomialbase.hh
@@ -47,7 +47,7 @@ namespace Dune {
//! point x and range phi
template <int diffOrd>
inline
- void evaluate ( int baseFunct, const Vec<diffOrd, deriType> &diffVariable,
+ void evaluate ( int baseFunct, const FieldVector<deriType, diffOrd> &diffVariable,
const Domain & x, Range & phi ) const;
//! evaluate base fucntion baseFunct at a given quadrature point
@@ -56,7 +56,7 @@ namespace Dune {
//! calulated again
template <int diffOrd, class QuadratureType>
inline
- void evaluate ( int baseFunct, const Vec<diffOrd, deriType> &diffVariable,
+ void evaluate ( int baseFunct, const FieldVector<deriType, diffOrd> &diffVariable,
QuadratureType & quad, int quadPoint, Range & phi ) const;
void print (std::ostream& s, int baseFunct) const {
@@ -79,7 +79,7 @@ namespace Dune {
int numOfBaseFct_;
//! vector which holds the base function pointers
- Dune::SimpleVector<Dune::Vec<DimDomain,int> > Phi_ ;
+ SimpleVector<FieldVector<int, DimDomain> > Phi_ ;
double power(double x, int p) const {
if (p <= 0)
@@ -88,16 +88,16 @@ namespace Dune {
}
inline
- void real_evaluate ( int baseFunct, const Vec<0, deriType> &diffVariable,
+ void real_evaluate ( int baseFunct, const FieldVector<deriType, 0> &diffVariable,
const Domain & x, Range & phi ) const;
inline
- void real_evaluate ( int baseFunct, const Vec<1, deriType> &diffVariable,
+ void real_evaluate ( int baseFunct, const FieldVector<deriType, 1> &diffVariable,
const Domain & x, Range & phi ) const;
inline
- void real_evaluate ( int baseFunct, const Vec<2, deriType> &diffVariable,
+ void real_evaluate ( int baseFunct, const FieldVector<deriType, 2> &diffVariable,
const Domain & x, Range & phi ) const;
- void print (std::ostream& s, const Dune::Vec<2,int> & pol) const;
+ void print (std::ostream& s, const FieldVector<int, 2> & pol) const;
}; // end class MonomialBaseFunctionSet
diff --git a/fem/discfuncarray.hh b/fem/discfuncarray.hh
index 73e478337..b038bb3ed 100644
--- a/fem/discfuncarray.hh
+++ b/fem/discfuncarray.hh
@@ -1,7 +1,7 @@
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
-#ifndef __DUNE_DISFUNCARRAY_HH__
-#define __DUNE_DISFUNCARRAY_HH__
+#ifndef __DUNE_DISCFUNCARRAY_HH__
+#define __DUNE_DISCFUNCARRAY_HH__
#include <dune/common/array.hh>
@@ -44,14 +44,21 @@ namespace Dune {
enum { myId_ = 0};
public:
+
+ //! ???
typedef typename DiscreteFunctionSpaceType::RangeField RangeFieldType;
+ //! ???
typedef typename DiscreteFunctionSpaceType::GridType GridType;
+ //! ???
typedef DiscFuncArray <DiscreteFunctionSpaceType> DiscreteFunctionType;
+ //! ???
typedef LocalFunctionArray < DiscreteFunctionSpaceType > LocalFunctionType;
+ //! ???
typedef DofIteratorArray < typename DiscreteFunctionSpaceType::RangeField > DofIteratorType;
+ //! ???
typedef DiscreteFunctionSpaceType FunctionSpaceType;
//! Constructor make empty DiscFuncArray
@@ -109,6 +116,9 @@ namespace Dune {
//! return if allLevels are used
bool allLevels () { return allLevels_; }
+ //! Return the name of the discrete function
+ const std::string& name() const {return name_;}
+
//! set all dofs to zero
void clearLevel( int level );
@@ -323,6 +333,7 @@ namespace Dune {
// --DofIteratorArray
//
//***********************************************************************
+ /** \brief ??? */
template < class DofType >
class DofIteratorArray : public
DofIteratorDefault < DofType , DofIteratorArray < DofType > >
diff --git a/fem/discfuncarray/discfuncarray.cc b/fem/discfuncarray/discfuncarray.cc
index 5772611ef..c5f76f2c2 100644
--- a/fem/discfuncarray/discfuncarray.cc
+++ b/fem/discfuncarray/discfuncarray.cc
@@ -657,14 +657,14 @@ namespace Dune
if(eval)
{
for(int l=0; l<dimrange; l++)
- ret(l) += (* (values_[i])) * tmp_(l);
+ ret[l] += (* (values_[i])) * tmp_[l];
}
}
}
else
{
for(int l=0; l<dimrange; l++)
- ret(l) = (* (values_[ l ]));
+ ret[l] = (* (values_[ l ]));
}
}
diff --git a/fem/fastquad.hh b/fem/fastquad.hh
index f474d2ce2..47c820dd5 100644
--- a/fem/fastquad.hh
+++ b/fem/fastquad.hh
@@ -36,7 +36,7 @@ namespace Dune {
FastQuad < RangeFieldType, DomainType, poly_order > >
{
private:
- enum { dim = DomainType::dimension };
+ enum { dim = DomainType::size };
//! my Type
typedef FastQuad < RangeFieldType , DomainType , poly_order > FastQuadType;
@@ -89,13 +89,13 @@ namespace Dune {
//! return weight for point i
const RangeFieldType& weight ( int i) const
{
- return (weights_(i));
+ return weights_[i];
}
//! return point i local coordinates
const DomainType& point (int i) const
{
- return (points_(i));
+ return points_[i];
}
private:
@@ -112,8 +112,8 @@ namespace Dune {
for(int i=0; i<numberOfQuadPoints_; i++)
{
- points_(i) = QuadInitializer::getPoint(i);
- weights_(i) = QuadInitializer::getWeight(i);
+ points_[i] = QuadInitializer::getPoint(i);
+ weights_[i] = QuadInitializer::getWeight(i);
}
int myType = (int) ElType;
@@ -130,8 +130,8 @@ namespace Dune {
int order_;
//! Vecs with constant length holding the weights and points
- Vec < maxQuadPoints , RangeFieldType > weights_;
- Vec < maxQuadPoints , DomainType > points_;
+ FieldVector<RangeFieldType, maxQuadPoints> weights_;
+ FieldVector<DomainType, maxQuadPoints> points_;
}; // end class FastQuadrature
@@ -241,8 +241,8 @@ namespace Dune {
int order_;
//! Vecs with constant length holding the weights and points
- Vec < maxQuadPoints , RangeFieldType > weights_;
- Vec < maxQuadPoints , DomainType > points_;
+ FieldVector<RangeFieldType, maxQuadPoints> weights_;
+ FieldVector<DomainType, maxQuadPoints> points_;
}; // end class FastQuadrature
diff --git a/fem/fastquad/gaussquadimp.cc b/fem/fastquad/gaussquadimp.cc
index fd9edf14a..f3d80e222 100644
--- a/fem/fastquad/gaussquadimp.cc
+++ b/fem/fastquad/gaussquadimp.cc
@@ -174,18 +174,18 @@ namespace Dune {
for (int i=0; i<n; ++i)
{
// compute integer coordinates of Gauss point from number
- Vec<dim,int> x (0);
+ FieldVector<int, dim> x (0);
int z = i;
for (int k=0; k<dim; ++k)
{
- x(k) = z%m;
+ x[k] = z%m;
z = z/m;
}
weight[i] = 1.0;
for (int k=0; k<dim; k++) {
- local[i](k) = G[x(k)];
- weight[i] *= W[x(k)];
+ local[i][k] = G[x[k]];
+ weight[i] *= W[x[k]];
}
}
}
diff --git a/fem/fastquad/quadtetratri.hh b/fem/fastquad/quadtetratri.hh
index 5fa736aee..50ed9a2bb 100644
--- a/fem/fastquad/quadtetratri.hh
+++ b/fem/fastquad/quadtetratri.hh
@@ -64,7 +64,7 @@ namespace Dune {
QUADRATURE * quad = UG_Quadratures::GetQuadrature(dim,numberOfCorners,polOrd);
Domain tmp;
for(int j=0; j<dim; j++)
- tmp(j) = quad->local[i][j];
+ tmp[j] = quad->local[i][j];
return tmp;
}
@@ -177,7 +177,7 @@ namespace Dune {
QUADRATURE * quad = UG_Quadratures::GetQuadrature(dim,numberOfCorners,polOrd);
Domain tmp;
for(int j=0; j<dim; j++)
- tmp(j) = quad->local[i][j];
+ tmp[j] = quad->local[i][j];
return tmp;
}
diff --git a/fem/feoperator.hh b/fem/feoperator.hh
index fedd97d61..6992724f8 100644
--- a/fem/feoperator.hh
+++ b/fem/feoperator.hh
@@ -127,6 +127,7 @@ namespace Dune {
}
}
+#if 0
{
// eliminate the Dirichlet rows and columns
typedef typename GridType::template Traits<0>::Entity EntityType;
@@ -194,6 +195,7 @@ namespace Dune {
}
}
}
+#endif
matrix_assembled_ = true;
}
diff --git a/fem/gaussquadrature.hh b/fem/gaussquadrature.hh
index cd778f3fc..5cd0e0654 100644
--- a/fem/gaussquadrature.hh
+++ b/fem/gaussquadrature.hh
@@ -45,14 +45,14 @@ namespace Dune {
int nip ();
//! return local coordinates of integration point i
- Vec<dim,ct>& ip (int i);
+ FieldVector<ct, dim>& ip (int i);
//! return weight associated with integration point i
double w (int i);
private:
int n;
- Vec<dim,ct> *local;
+ FieldVector<ct, dim> *local;
double *weight;
int power (int y, int d);
diff --git a/fem/gaussquadrature/gaussquadrature.cc b/fem/gaussquadrature/gaussquadrature.cc
index 655459681..4ccbbf291 100644
--- a/fem/gaussquadrature/gaussquadrature.cc
+++ b/fem/gaussquadrature/gaussquadrature.cc
@@ -180,7 +180,7 @@ namespace Dune {
// allocate memory for the quadrature points and weights
try {
- local = new Vec<dim,ct>[n];
+ local = new FieldVector<ct, dim>[n];
}
catch (std::bad_alloc) {
std::cout << "failed to allocate memory for quadrature points" << std::endl;
@@ -196,7 +196,7 @@ namespace Dune {
for (int i=0; i<n; ++i)
{
// compute integer coordinates of Gauss point from number
- Vec<dim,short> x = 0;
+ FieldVector<short, dim> x = 0;
int z = i;
for (int k=0; k<dim; ++k)
{
@@ -209,7 +209,7 @@ namespace Dune {
local[i](k) = G[x(k)];
weight[i] *= W[x(k)];
}
- // std::cout << i << " " << x << " " << local[i] << " " << w[i] << std::endl;
+ // std::cout << i << " " << x << " " << local[i] << " " << w[i] << std::endl;
}
}
@@ -221,7 +221,7 @@ namespace Dune {
// allocate memory for the quadrature points and weights
try {
- local = new Vec<dim,ct>[n];
+ local = new FieldVector<ct, dim>[n];
}
catch (std::bad_alloc) {
std::cout << "failed to allocate memory for quadrature points" << std::endl;
@@ -266,7 +266,7 @@ namespace Dune {
// allocate memory for the quadrature points and weights
try {
- local = new Vec<dim,ct>[n];
+ local = new FieldVector<ct, dim>[n];
}
catch (std::bad_alloc) {
std::cout << "failed to allocate memory for quadrature points" << std::endl;
@@ -297,7 +297,7 @@ namespace Dune {
}
template<int dim, class ct>
- inline Vec<dim,ct>& GaussQuadrature<dim,ct>::ip (int i)
+ inline FieldVector<ct, dim>& GaussQuadrature<dim,ct>::ip (int i)
{
return local[i];
}
diff --git a/fem/lagrangebase.hh b/fem/lagrangebase.hh
index 11cbf736d..8d07a8eeb 100644
--- a/fem/lagrangebase.hh
+++ b/fem/lagrangebase.hh
@@ -141,7 +141,7 @@ namespace Dune {
//! the corresponding vector of base function sets
//! length is different element types
- Vec< numOfDiffBase_ , FastBaseFunctionSetType * > baseFuncSet_;
+ FieldVector<FastBaseFunctionSetType*, numOfDiffBase_ > baseFuncSet_;
protected:
//! DofManager manages the memory
@@ -243,17 +243,17 @@ namespace Dune {
}
/** \todo Please doc me! */
- virtual void evaluate ( const Vec<0, deriType> &diffVariable,
+ virtual void evaluate ( const FieldVector<deriType, 0> &diffVariable,
const Domain & x, Range & phi) const
{ // q(x) = (x - point_ ) * 1/(2|T|) mit |T|=0.5
phi = (x - point_);
}
/** \todo Please doc me! */
- virtual void evaluate ( const Vec<1, deriType> &diffVariable,
+ virtual void evaluate ( const FieldVector<deriType, 1> &diffVariable,
const Domain & x, Range & phi) const
{
- int comp = diffVariable(0);
+ int comp = diffVariable[0];
phi = 0.0;
phi(comp) = 1.0;
}
@@ -302,7 +302,7 @@ namespace Dune {
/** \todo Please doc me! */
- virtual void evaluate ( const Vec<0, deriType> &diffVariable,
+ virtual void evaluate ( const FieldVector<deriType, 0> &diffVariable,
const Domain & x, Range & phi) const
{
phi = factor[2];
@@ -311,11 +311,11 @@ namespace Dune {
}
/** \todo Please doc me! */
- virtual void evaluate ( const Vec<1, deriType> &diffVariable,
+ virtual void evaluate ( const FieldVector<deriType, 1> &diffVariable,
const Domain & x, Range & phi) const
{
// x or y ==> 1 or 2
- int num = diffVariable(0);
+ int num = diffVariable[0];
phi = factor[num];
}
@@ -408,7 +408,7 @@ namespace Dune {
}
private:
//! Vector with all base functions corresponding to the base function set
- Vec < numOfBaseFct , RaviartThomasBaseFunctionType *> baseFuncList_;
+ FieldVector <RaviartThomasBaseFunctionType*, numOfBaseFct> baseFuncList_;
};
//*******************************************************************
@@ -463,7 +463,7 @@ namespace Dune {
}
private:
//! Vector with all base functions corresponding to the base function set
- Vec < numOfBaseFct , EdgeBaseFunctionType *> baseFuncList_;
+ FieldVector<EdgeBaseFunctionType*, numOfBaseFct> baseFuncList_;
};
@@ -651,7 +651,7 @@ namespace Dune {
//! the corresponding vector of base function sets
//! lenght is diffrent element types
- Vec< numOfDiffBase_ , FastBaseFunctionSetType * > baseFuncSet_;
+ FieldVector< FastBaseFunctionSetType*, numOfDiffBase_ > baseFuncSet_;
//! make base function set depending on ElementType and polynomial order
template <ElementType ElType, int pO >
@@ -872,7 +872,7 @@ namespace Dune {
//! the corresponding vector of base function sets
//! lenght is diffrent element types
- Vec< numOfDiffBase_ , FastBaseFunctionSetType * > baseFuncSet_;
+ FieldVector<FastBaseFunctionSetType*, numOfDiffBase_ > baseFuncSet_;
//! make base function set depending on ElementType and polynomial order
template <ElementType ElType, int pO >
diff --git a/fem/lagrangebase/lagrangebasefunctions.hh b/fem/lagrangebase/lagrangebasefunctions.hh
index 5e54685af..4feba6e5d 100644
--- a/fem/lagrangebase/lagrangebasefunctions.hh
+++ b/fem/lagrangebase/lagrangebasefunctions.hh
@@ -34,20 +34,20 @@ namespace Dune {
assert((baseNum_ >= 0) || (baseNum_ < DimRange));
};
- virtual void evaluate ( const Vec<0, deriType> &diffVariable,
+ virtual void evaluate ( const FieldVector<deriType, 0> &diffVariable,
const Domain & x, Range & phi) const
{
phi = 0.0;
phi(baseNum_) = 1.0;
}
- virtual void evaluate ( const Vec<1, deriType> &diffVariable,
+ virtual void evaluate ( const FieldVector<deriType, 1> &diffVariable,
const Domain & x, Range & phi) const
{
phi = 0.0;
}
- virtual void evaluate ( const Vec<2,deriType> &diffVariable,
+ virtual void evaluate ( const FieldVector<deriType, 2> &diffVariable,
const Domain & x, Range & phi) const
{
phi = 0.0 ;
@@ -88,23 +88,23 @@ namespace Dune {
}
//! evaluate the function
- virtual void evaluate ( const Vec<0, deriType> &diffVariable,
+ virtual void evaluate ( const FieldVector<deriType, 0> &diffVariable,
const Domain & x, Range & phi) const
{
phi = factor[1];
- phi += factor[0] * x(0);
+ phi += factor[0] * x[0];
}
//! evaluate first derivative
- virtual void evaluate ( const Vec<1, deriType> &diffVariable,
+ virtual void evaluate ( const FieldVector<deriType, 1> &diffVariable,
const Domain & x, Range & phi) const
{
- int num = diffVariable(0);
+ int num = diffVariable[0];
phi = factor[num];
}
//! evaluate second derivative
- virtual void evaluate ( const Vec<2,deriType> &diffVariable,
+ virtual void evaluate ( const FieldVector<deriType, 2> &diffVariable,
const Domain & x, Range & phi) const
{
// function is linear, therefore
@@ -163,19 +163,19 @@ namespace Dune {
}
}
//! evaluate the function
- virtual void evaluate ( const Vec<0, deriType> &diffVariable,
+ virtual void evaluate ( const FieldVector<deriType, 0> &diffVariable,
const Domain & x, Range & phi) const
{
phi = factor[0];
for(int i=1; i<3; i++)
- phi += factor[i] * x(i-1);
+ phi += factor[i] * x[i-1];
}
- virtual void evaluate ( const Vec<1, deriType> &diffVariable,
+ virtual void evaluate ( const FieldVector<deriType, 1> &diffVariable,
const Domain & x, Range & phi) const
{
// x or y ==> 1 or 2
- int num = diffVariable(0);
+ int num = diffVariable[0];
assert( (num >= 0) && ( num <= 1));
phi = factor[num+1];
}
@@ -225,20 +225,20 @@ namespace Dune {
}
//! evaluate function
- virtual void evaluate ( const Vec<0, deriType> &diffVariable,
+ virtual void evaluate ( const FieldVector<deriType, 0> &diffVariable,
const Domain & x, Range & phi) const
{
phi = factor[0];
for(int i=1; i<4; i++)
- phi += factor[i]*x(i-1);
+ phi += factor[i]*x[i-1];
}
//! first Derivative
- virtual void evaluate ( const Vec<1, deriType> &diffVariable,
+ virtual void evaluate ( const FieldVector<deriType, 1> &diffVariable,
const Domain & x, Range & phi) const
{
// num = 0 ==> derivative respect to x
- int num = diffVariable(0);
+ int num = diffVariable[0];
phi = factor[num+1];
}
@@ -290,23 +290,23 @@ namespace Dune {
};
//! evaluate the basefunction on point x
- virtual void evaluate ( const Vec<0,deriType> &diffVariable,
+ virtual void evaluate ( const FieldVector<deriType, 0> &diffVariable,
const Domain & x, Range & phi) const
{
// dim == 2, tested
phi = 1.0;
for(int i=0; i<dim; i++)
- phi *= (factor[i][0] + (factor[i][1] * x(i)));
+ phi *= (factor[i][0] + (factor[i][1] * x[i]));
}
//! derivative with respect to x or y
//! diffVariable(0) == 0 ==> x
//! diffVariable(0) == 1 ==> y
//! diffVariable(0) == 2 ==> z, and so on
- virtual void evaluate ( const Vec<1,deriType> &diffVariable,
+ virtual void evaluate ( const FieldVector<deriType, 1> &diffVariable,
const Domain & x, Range & phi) const
{
- int num = diffVariable(0);
+ int num = diffVariable[0];
assert( (num >= 0) && ( num <= 1));
phi = 1.0;
for(int i=0; i<dim; i++)
@@ -314,17 +314,17 @@ namespace Dune {
if(num == i)
phi *= factor[num][1];
else
- phi *= (factor[i][0] + factor[i][1] * x(i));
+ phi *= (factor[i][0] + factor[i][1] * x[i]);
}
return;
}
//! evaluate second derivative
- virtual void evaluate ( const Vec<2,deriType> &diffVariable,
+ virtual void evaluate ( const FieldVector<deriType, 2> &diffVariable,
const Domain & x, Range & phi) const
{
// which means derivative xx or yy
- if(diffVariable(0) == diffVariable(1))
+ if(diffVariable[0] == diffVariable[1])
{
phi = 0.0;
return;
@@ -400,36 +400,36 @@ namespace Dune {
};
//! evaluate the basefunction on point x
- virtual void evaluate ( const Vec<0,deriType> &diffVariable,
+ virtual void evaluate ( const FieldVector<deriType, 0> &diffVariable,
const Domain & x, Range & phi) const
{
// dim == 3, tested
phi = 1.0;
for(int i=0; i<dim; i++)
- phi *= (factor[i][0] + (factor[i][1] * x(i)));
+ phi *= (factor[i][0] + (factor[i][1] * x[i]));
}
//! derivative with respect to x or y or z
//! diffVariable(0) == 0 ==> x
//! diffVariable(0) == 1 ==> y
//! diffVariable(0) == 2 ==> z, and so on
- virtual void evaluate ( const Vec<1,deriType> &diffVariable,
+ virtual void evaluate ( const FieldVector<deriType, 1> &diffVariable,
const Domain & x, Range & phi) const
{
- int num = diffVariable(0);
+ int num = diffVariable[0];
phi = 1.0;
for(int i=0; i<dim; i++)
{
if(num == i)
phi *= factor[num][1];
else
- phi *= (factor[i][0] + factor[i][1] * x(i));
+ phi *= (factor[i][0] + factor[i][1] * x[i]);
}
return;
}
//! evaluate second derivative
- virtual void evaluate ( const Vec<2,deriType> &diffVariable,
+ virtual void evaluate ( const FieldVector<deriType, 2> &diffVariable,
const Domain & x, Range & phi) const
{
std::cout << "BaseFunction for hexahedron, evaluate 2nd derivative not implemented! \n";
@@ -517,8 +517,8 @@ namespace Dune {
{
for(int k=0; k<dimrange; k++)
{
- baseFuncList_(i*dimrange + k) = new LagrangeBaseFunctionType ( fuSpace, i ) ;
- this->setBaseFunctionPointer ( i*dimrange + k , baseFuncList_ (i*dimrange + k) );
+ baseFuncList_[i*dimrange + k] = new LagrangeBaseFunctionType ( fuSpace, i ) ;
+ this->setBaseFunctionPointer ( i*dimrange + k , baseFuncList_[i*dimrange + k] );
}
}
this->setNumOfDiffFct ( numOfDifferentFuncs );
@@ -542,7 +542,7 @@ namespace Dune {
}
private:
//! Vector with all base functions corresponding to the base function set
- Vec < numOfBaseFct , LagrangeBaseFunctionType *> baseFuncList_;
+ FieldVector <LagrangeBaseFunctionType*, numOfBaseFct> baseFuncList_;
};
} // end namespace Dune
diff --git a/fem/lagrangebase/lagrangespace.cc b/fem/lagrangebase/lagrangespace.cc
index b0ceb711c..cdb076db5 100644
--- a/fem/lagrangebase/lagrangespace.cc
+++ b/fem/lagrangebase/lagrangespace.cc
@@ -16,7 +16,7 @@ namespace Dune {
//std::cout << "Constructor of LagrangeDiscreteFunctionSpace! \n";
for(int i=0; i<numOfDiffBase_; i++)
- baseFuncSet_(i) = 0;
+ baseFuncSet_[i] = 0;
{
// search the macro grid for diffrent element types
@@ -25,15 +25,15 @@ namespace Dune {
for(LevelIteratorType it = g.template lbegin<0>(0); it != endit; ++it)
{
ElementType type = (*it).geometry().type(); // Hack
- if(baseFuncSet_( type ) == 0 )
- baseFuncSet_ ( type ) = setBaseFuncSetPointer(*it);
+ if(baseFuncSet_[type] == 0 )
+ baseFuncSet_[type] = setBaseFuncSetPointer(*it);
}
}
for(int i=0; i<numOfDiffBase_; i++)
{
- if(baseFuncSet_(i))
- maxNumBase_ = MAX(baseFuncSet_(i)->getNumberOfBaseFunctions(),maxNumBase_);
+ if(baseFuncSet_[i])
+ maxNumBase_ = MAX(baseFuncSet_[i]->getNumberOfBaseFunctions(),maxNumBase_);
}
}
@@ -42,8 +42,8 @@ namespace Dune {
~LagrangeDiscreteFunctionSpace ()
{
for(int i=0; i<numOfDiffBase_; i++)
- if (baseFuncSet_(i) != 0)
- delete baseFuncSet_(i);
+ if (baseFuncSet_[i] != 0)
+ delete baseFuncSet_[i];
}
template< class FunctionSpaceType, class GridType,int polOrd, class DofManagerType >
@@ -62,7 +62,7 @@ namespace Dune {
getBaseFunctionSet (EntityType &en) const
{
ElementType type = en.geometry().type();
- return (*baseFuncSet_( type ));
+ return *baseFuncSet_[type];
}
template< class FunctionSpaceType, class GridType,int polOrd, class DofManagerType >
diff --git a/fem/norms/l2norm.hh b/fem/norms/l2norm.hh
index c1065cb83..8348f3b8c 100644
--- a/fem/norms/l2norm.hh
+++ b/fem/norms/l2norm.hh
@@ -48,7 +48,7 @@ namespace Dune {
for(int qP = 0; qP < quad.nop(); qP++) {
lf.evaluate((*it),quad,qP,phi);
- sum += det * quad.weight(qP) * SQR(phi(0));
+ sum += det * quad.weight(qP) * SQR(phi[0]);
}
}
diff --git a/fem/operator/discreteoperator.cc b/fem/operator/discreteoperator.cc
index 9e8589833..8d9fada5d 100644
--- a/fem/operator/discreteoperator.cc
+++ b/fem/operator/discreteoperator.cc
@@ -150,7 +150,7 @@ namespace Dune
enum { dimdef = 2};
- Vec<GRID::dimension> tmp(1.0);
+ FieldVector<double, GRID::dimension> tmp(1.0);
GRID &grid = (*func.getGrid());
FuncSpace &funcSpace = (*func.getFuncSpace());
@@ -168,7 +168,7 @@ namespace Dune
for(int p=0; p<vx; p++)
{
row = funcSpace.mapIndex((*it),p);
- Vec<FuncSpace::dimdef,double> vec1 =
+ FieldVector<double, FuncSpace::dimdef> vec1 =
//grad[p];
(funcSpace.getLocalBaseFunc(p)->evalFirstDrv (tmp))(0);
vec1 = inv*vec1;
@@ -182,7 +182,7 @@ namespace Dune
val = vol;
col = funcSpace.mapIndex((*it),oth);
- Vec<FuncSpace::dimdef,double> vec2 =
+ FieldVector<double, FuncSpace::dimdef> vec2 =
//grad[oth];
(funcSpace.getLocalBaseFunc(oth)->evalFirstDrv (tmp))(0);
@@ -328,7 +328,7 @@ namespace Dune
}
SpaceDiscr& fake = (*fv_);
- Vec<1> tmp(0.0);
+ FieldVector<double, 1> tmp(0.0);
TimeGrid::LevelIterator endit = timegrid_->lend<0>(0);
for(TimeGrid::LevelIterator it = timegrid_->lbegin<0>(0); it != endit; ++it)
diff --git a/fem/operator/laplace.hh b/fem/operator/laplace.hh
index b0c591b39..f957a1400 100644
--- a/fem/operator/laplace.hh
+++ b/fem/operator/laplace.hh
@@ -173,15 +173,15 @@ namespace Dune
if(stiffTensor_) {
stiffTensor_->evaluate(entity.geometry().global(quad.point(pt)),ret);
- ret(0) *= quad.weight( pt );
+ ret[0] *= quad.weight( pt );
for(i=0; i<matSize; i++)
for (j=0; j<=i; j++ )
- mat(i,j) += ( mygrad[i](0) * mygrad[j](0) ) * ret(0);
+ mat(i,j) += ( mygrad[i][0] * mygrad[j][0] ) * ret[0];
}
else{
for(i=0; i<matSize; i++)
for (j=0; j<=i; j++ )
- mat(i,j) += ( mygrad[i](0) * mygrad[j](0) ) * quad.weight( pt );
+ mat(i,j) += ( mygrad[i][0] * mygrad[j][0] ) * quad.weight( pt );
}
diff --git a/fem/operator/massmatrix.hh b/fem/operator/massmatrix.hh
index 172b720c5..1507f91ad 100644
--- a/fem/operator/massmatrix.hh
+++ b/fem/operator/massmatrix.hh
@@ -61,7 +61,7 @@ namespace Dune {
typedef typename FunctionSpaceType::BaseFunctionSetType BaseFunctionSetType;
double val = 0.;
- Vec<GridType::dimension> tmp(1.0);
+ FieldVector<double, GridType::dimension> tmp(1.0);
const BaseFunctionSetType & baseSet
= this->functionSpace_.getBaseFunctionSet( entity );
@@ -90,7 +90,8 @@ namespace Dune {
const BaseFunctionSetType & baseSet
= this->functionSpace_.getBaseFunctionSet( entity );
- static Vec<GridType::dimension> tmp(1.0);
+ /** \todo What's the correct type here? */
+ static FieldVector<double, GridType::dimension> tmp(1.0);
const double vol = entity.geometry().integration_element(tmp);
static RangeType v[4];
diff --git a/fem/transfer/mgtransfer.cc b/fem/transfer/mgtransfer.cc
index ca63c39d8..9460bde11 100644
--- a/fem/transfer/mgtransfer.cc
+++ b/fem/transfer/mgtransfer.cc
@@ -67,7 +67,7 @@ void Dune::MGTransfer<DiscFuncType>::setup(const FunctionSpaceType& fS, int cL,
// Evaluate coarse grid base function at the location of the fine grid dof
// first determine local fine grid dof position
- Vec<GridType::dimension, double> local = cIt->geometry().local(fIt->geometry()[j]);
+ FieldVector<double, GridType::dimension> local = cIt->geometry().local(fIt->geometry()[j]);
//std::cout << "finePos" << fIt->geometry()[j] << "\n";
//std::cout << "locally" << local << "\n";
@@ -75,7 +75,7 @@ void Dune::MGTransfer<DiscFuncType>::setup(const FunctionSpaceType& fS, int cL,
// Evaluate coarse grid base function
typename FunctionSpaceType::Range value;
- Vec<0, int> diffVariable;
+ FieldVector<int, 0> diffVariable;
coarseBaseSet.evaluate(i, diffVariable, local, value);
//std::cout << "value: " << value << "\n";
@@ -154,5 +154,67 @@ template<class DiscFuncType>
Dune::SparseRowMatrix<double> Dune::MGTransfer<DiscFuncType>::
galerkinRestrict(const Dune::SparseRowMatrix<double>& fineMat) const
{
- return matrix_.applyFromLeftAndRightTo(fineMat);
+ typedef typename SparseRowMatrix<double>::ColumnIterator ColumnIterator;
+ //return matrix_.applyFromLeftAndRightTo(fineMat);
+
+ // Hack: We need the transposed matrix
+
+ SparseRowMatrix<double> transpose(matrix_.cols(), matrix_.rows(), matrix_.numNonZeros());
+
+ for (int row=0; row<matrix_.rows(); row++) {
+
+ /** \todo Remove casts */
+ ColumnIterator cIt = const_cast<SparseRowMatrix<double>&>(matrix_).template rbegin(row);
+ ColumnIterator cEndIt = const_cast<SparseRowMatrix<double>&>(matrix_).template rend(row);
+
+ for(; cIt!=cEndIt; ++cIt)
+ transpose.set(cIt.col(), row, *cIt);
+
+ }
+
+ SparseRowMatrix<double> result(matrix_.rows() ,matrix_.rows() , fineMat.numNonZeros());
+
+ //for (v=FIRSTVECTOR(FineGrid); v!=NULL; v=SUCCVC(v))
+ for (int row=0; row<fineMat.rows(); row++) {
+
+ // for (m=VSTART(v); m!=NULL; m=MNEXT(m)) {
+ // w = MDEST(m);
+ ColumnIterator cIt = const_cast<SparseRowMatrix<double>&>(fineMat).template rbegin(row);
+ ColumnIterator cEndIt = const_cast<SparseRowMatrix<double>&>(fineMat).template rend(row);
+
+ for(; cIt!=cEndIt; ++cIt) {
+
+ //mvalue = MVALUE(m,mc);
+ double mvalue = *cIt;
+
+ ColumnIterator tciIt = transpose.template rbegin(row);
+ ColumnIterator tciEndIt = transpose.template rend(row);
+
+ //for (im=VISTART(v); im!= NULL; im = NEXT(im)) {
+ for (; tciIt!=tciEndIt; ++tciIt) {
+
+ //fac = mvalue*MVALUE(im,0);
+ double fac = mvalue* (*tciIt);
+
+ ColumnIterator tcjIt = transpose.template rbegin(cIt.col());
+ ColumnIterator tcjEndIt = transpose.template rend(cIt.col());
+
+ //for (jm=VISTART(v); jm!= NULL; jm = NEXT(jm)) {
+ for (; tcjIt!=tcjEndIt; ++tcjIt) {
+
+ // jv = MDEST(jm);
+ // cm = GetMatrix(iv,jv);
+ // MVALUE(cm,mc) +=
+ // fac * MVALUE(jm,0);
+ result.add( tciIt.col(), tcjIt.col(), fac* (*tcjIt));
+ }
+
+ }
+
+ }
+
+
+ }
+
+ return result;
}
diff --git a/grid/albertgrid.hh b/grid/albertgrid.hh
index 1dcc6f13f..ff08f238b 100644
--- a/grid/albertgrid.hh
+++ b/grid/albertgrid.hh
@@ -1,7 +1,7 @@
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
-#ifndef DUNE_ALBERTGRID_HH
-#define DUNE_ALBERTGRID_HH
+#ifndef __DUNE_ALBERTGRID_HH__
+#define __DUNE_ALBERTGRID_HH__
#include <iostream>
#include <fstream>
@@ -190,7 +190,7 @@ namespace Dune
int corners ();
//! access to coordinates of corners. Index is the number of the corner
- Vec<dimworld,albertCtype>& operator[] (int i);
+ FieldVector<albertCtype, dimworld>& operator[] (int i);
/*! return reference element corresponding to this element. If this is
a reference element then self is returned.
@@ -199,14 +199,14 @@ namespace Dune
//! maps a local coordinate within reference element to
//! global coordinate in element
- Vec<dimworld,albertCtype> global (const Vec<dim,albertCtype>& local);
+ FieldVector<albertCtype, dimworld> global (const FieldVector<albertCtype, dim>& local);
//! maps a global coordinate within the element to a
//! local coordinate in its reference element
- Vec<dim,albertCtype> local (const Vec<dimworld,albertCtype>& global);
+ FieldVector<albertCtype, dim> local (const FieldVector<albertCtype, dimworld>& global);
//! returns true if the point in local coordinates is inside reference element
- bool checkInside(const Vec<dim,albertCtype>& local);
+ bool checkInside(const FieldVector<albertCtype, dim>& local);
/*!
Copy from sgrid.hh:
@@ -232,13 +232,13 @@ namespace Dune
*/
// A(l)
- albertCtype integration_element (const Vec<dim,albertCtype>& local);
+ albertCtype integration_element (const FieldVector<albertCtype, dim>& local);
//! can only be called for dim=dimworld!
//! Note that if both methods are called on the same element, then
//! call Jacobian_inverse first because integration element is calculated
//! during calculation of the Jacobian_inverse
- Mat<dim,dim>& Jacobian_inverse (const Vec<dim,albertCtype>& local);
+ Mat<dim,dim>& Jacobian_inverse (const FieldVector<albertCtype, dim>& local);
//***********************************************************************
//! Methods that not belong to the Interface, but have to be public
@@ -271,7 +271,7 @@ namespace Dune
//! built the jacobian inverse and store the volume
void buildJacobianInverse ();
- Vec<dim+1,albertCtype> tmpVec_;
+ FieldVector<albertCtype, dim+1> tmpVec_;
//! maps a global coordinate within the elements local barycentric
//! koordinates
//Vec<dim+1,albertCtype> localBary (const Vec<dimworld,albertCtype>& global);
@@ -287,10 +287,10 @@ namespace Dune
Mat<dimworld,dim+1,albertCtype> coord_;
//! storage for global coords
- Vec<dimworld,albertCtype> globalCoord_;
+ FieldVector<albertCtype, dimworld> globalCoord_;
//! storage for local coords
- Vec<dim,albertCtype> localCoord_;
+ FieldVector<albertCtype, dim> localCoord_;
// make empty EL_INFO
ALBERT EL_INFO * makeEmptyElInfo();
@@ -369,7 +369,7 @@ namespace Dune
AlbertGridLevelIterator<0,dim,dimworld,All_Partition> father ();
//! local coordinates within father
- Vec<dim,albertCtype>& local ();
+ FieldVector<albertCtype, dim>& local ();
//! no interface method
//! returns the global vertex number as default
@@ -399,7 +399,7 @@ namespace Dune
AlbertGridElement<dim-codim,dimworld> geo_;
bool builtgeometry_; //!< true if geometry has been constructed
- Vec<dim,albertCtype> localFatherCoords_;
+ FieldVector<albertCtype, dim> localFatherCoords_;
//! element number
int elNum_;
@@ -798,10 +798,10 @@ namespace Dune
//! return unit outer normal, this should be dependent on local
//! coordinates for higher order boundary
- Vec<dimworld,albertCtype>& unit_outer_normal (Vec<dim-1,albertCtype>& local);
+ FieldVector<albertCtype, dimworld>& unit_outer_normal (FieldVector<albertCtype, dim-1>& local);
//! return unit outer normal, if you know it is constant use this function instead
- Vec<dimworld,albertCtype>& unit_outer_normal ();
+ FieldVector<albertCtype, dimworld>& unit_outer_normal ();
//! intersection of codimension 1 of this neighbor with element where
//! iteration started.
@@ -829,10 +829,10 @@ namespace Dune
//! return outer normal, this should be dependent on local
//! coordinates for higher order boundary
- Vec<dimworld,albertCtype>& outer_normal (Vec<dim-1,albertCtype>& local);
+ FieldVector<albertCtype, dimworld>& outer_normal (FieldVector<albertCtype, dim-1>& local);
//! return unit outer normal, if you know it is constant use this function instead
- Vec<dimworld,albertCtype>& outer_normal ();
+ FieldVector<albertCtype, dimworld>& outer_normal ();
private:
//**********************************************************
@@ -894,11 +894,11 @@ namespace Dune
//! EL_INFO th store the information of the neighbor if needed
ALBERT EL_INFO * neighElInfo_;
- Vec<dimworld,albertCtype> outNormal_;
+ FieldVector<albertCtype, dimworld> outNormal_;
// tmp memory for normal calculation
- Vec<dimworld,albertCtype> tmpU_;
- Vec<dimworld,albertCtype> tmpV_;
+ FieldVector<albertCtype, dimworld> tmpU_;
+ FieldVector<albertCtype, dimworld> tmpV_;
};
@@ -1171,7 +1171,7 @@ namespace Dune
int myProcessor () const { return myProc_; };
//! transform grid N = scalar * x + trans
- void setNewCoords(const Vec<dimworld,albertCtype> & trans, const albertCtype scalar);
+ void setNewCoords(const FieldVector<albertCtype, dimworld> & trans, const albertCtype scalar);
//! calculate max edge length
double calcGridWidth ();
diff --git a/grid/albertgrid/albertgrid.cc b/grid/albertgrid/albertgrid.cc
index 64191a559..950e57dd0 100644
--- a/grid/albertgrid/albertgrid.cc
+++ b/grid/albertgrid/albertgrid.cc
@@ -359,7 +359,7 @@ namespace Dune
///////////////////////////////////////////////////////////////////////
template< int dim, int dimworld>
- inline Vec<dimworld,albertCtype>& AlbertGridElement<dim,dimworld>::
+ inline FieldVector<albertCtype, dimworld>& AlbertGridElement<dim,dimworld>::
operator [](int i)
{
return coord_(i);
@@ -401,8 +401,8 @@ namespace Dune
}
template< int dim, int dimworld>
- inline Vec<dimworld,albertCtype> AlbertGridElement<dim,dimworld>::
- global(const Vec<dim>& local)
+ inline FieldVector<albertCtype, dimworld> AlbertGridElement<dim,dimworld>::
+ global(const FieldVector<albertCtype, dim>& local)
{
// 1hecked, works
@@ -463,7 +463,7 @@ namespace Dune
enum { dimworld = 3 };
if( !builtElMat_)
{
- Vec<dimworld,albertCtype> & coord0 = coord_(0);
+ FieldVector<albertCtype, dimworld> & coord0 = coord_(0);
for(int i=0 ; i<dimworld; i++)
{
elMat_(i,0) = coord_(i,1) - coord0(i);
@@ -477,8 +477,8 @@ namespace Dune
// uses the element matrix, because faster
template<>
- inline Vec<2,albertCtype> AlbertGridElement<2,2>::
- global(const Vec<2>& local)
+ inline FieldVector<albertCtype, 2> AlbertGridElement<2,2>::
+ global(const FieldVector<albertCtype, 2>& local)
{
calcElMatrix();
globalCoord_ = elMat_ * local;
@@ -487,8 +487,8 @@ namespace Dune
}
template<>
- inline Vec<3,albertCtype> AlbertGridElement<3,3>::
- global(const Vec<3>& local)
+ inline FieldVector<albertCtype, 3> AlbertGridElement<3,3>::
+ global(const FieldVector<albertCtype, 3>& local)
{
calcElMatrix();
globalCoord_ = elMat_ * local;
@@ -498,17 +498,17 @@ namespace Dune
template< int dim, int dimworld>
- inline Vec<dim> AlbertGridElement<dim,dimworld>::
- local(const Vec<dimworld>& global)
+ inline FieldVector<albertCtype, dim> AlbertGridElement<dim,dimworld>::
+ local(const FieldVector<albertCtype, dimworld>& global)
{
- std::cerr << "local for dim != dimworld not implemented! \n";
+ std::cerr << "global for dim != dimworld not implemented! \n";
abort();
return localCoord_;
}
template <>
- inline Vec<2> AlbertGridElement<2,2>::
- local(const Vec<2>& global)
+ inline FieldVector<albertCtype, 2> AlbertGridElement<2,2>::
+ local(const FieldVector<albertCtype, 2>& global)
{
if(!builtinverse_)
buildJacobianInverse();
@@ -518,8 +518,8 @@ namespace Dune
}
template <>
- inline Vec<3> AlbertGridElement<3,3>::
- local(const Vec<3>& global)
+ inline FieldVector<albertCtype, 3> AlbertGridElement<3,3>::
+ local(const FieldVector<albertCtype, 3>& global)
{
if(!builtinverse_)
buildJacobianInverse();
@@ -580,7 +580,7 @@ namespace Dune
buildJacobianInverse()
{
// volume is length of edge
- Vec<2,albertCtype> vec = coord_(0) - coord_(1);
+ FieldVector<albertCtype, 2> vec = coord_(0) - coord_(1);
elDet_ = vec.norm2();
builtinverse_ = true;
@@ -612,7 +612,7 @@ namespace Dune
template< int dim, int dimworld>
inline albertCtype AlbertGridElement<dim,dimworld>::
- integration_element (const Vec<dim,albertCtype>& local)
+ integration_element (const FieldVector<albertCtype, dim>& local)
{
// if inverse was built, volume was calced already
if(builtinverse_)
@@ -624,7 +624,7 @@ namespace Dune
template <>
inline Mat<1,1>& AlbertGridElement<1,2>::
- Jacobian_inverse (const Vec<1,albertCtype>& global)
+ Jacobian_inverse (const FieldVector<albertCtype, 1>& global)
{
std::cout << "Jaconbian_inverse for dim=1,dimworld=2 not implemented yet! \n";
return Jinv_;
@@ -632,7 +632,7 @@ namespace Dune
template< int dim, int dimworld>
inline Mat<dim,dim>& AlbertGridElement<dim,dimworld>::
- Jacobian_inverse (const Vec<dim,albertCtype>& global)
+ Jacobian_inverse (const FieldVector<albertCtype, dim>& global)
{
if(builtinverse_)
return Jinv_;
@@ -659,11 +659,11 @@ namespace Dune
// checks if F^-1 (x_i) == xref_i
enum { dim =2 };
- Vec<dim> & coord = coord_(loc);
- Vec<dim> & refcoord = refelem()[loc];
+ FieldVector<albertCtype, dim> & coord = coord_(loc);
+ FieldVector<albertCtype, dim> & refcoord = refelem()[loc];
buildJacobianInverse();
- Vec<dim> tmp2 = coord - coord_(0);
+ FieldVector<albertCtype, dim> tmp2 = coord - coord_(0);
tmp2 = Jinv_ * tmp2;
for(int j=0; j<dim; j++)
@@ -681,11 +681,11 @@ namespace Dune
// checks if F^-1 (x_i) == xref_i
enum { dim = 3 };
- Vec<dim> & coord = coord_(loc);
- Vec<dim> & refcoord = refelem()[loc];
+ FieldVector<albertCtype, dim> & coord = coord_(loc);
+ FieldVector<albertCtype, dim> & refcoord = refelem()[loc];
buildJacobianInverse();
- Vec<dim> tmp2 = coord - coord_(0);
+ FieldVector<albertCtype, dim> tmp2 = coord - coord_(0);
tmp2 = Jinv_ * tmp2;
for(int j=0; j<dim; j++)
@@ -715,10 +715,10 @@ namespace Dune
calcElMatrix ();
- Vec<dim> & coord = coord_(loc);
- Vec<dim> & refcoord = refelem()[loc];
+ FieldVector<albertCtype, dim> & coord = coord_(loc);
+ FieldVector<albertCtype, dim> & refcoord = refelem()[loc];
- Vec<dim> tmp2 = elMat_ * refcoord;
+ FieldVector<albertCtype, dim> tmp2 = elMat_ * refcoord;
tmp2 += coord_(0);
for(int j=0; j<dim; j++)
@@ -741,10 +741,10 @@ namespace Dune
calcElMatrix ();
- Vec<dim> & coord = coord_(loc);
- Vec<dim> & refcoord = refelem()[loc];
+ FieldVector<albertCtype, dim> & coord = coord_(loc);
+ FieldVector<albertCtype, dim> & refcoord = refelem()[loc];
- Vec<dim> tmp2 = elMat_ * refcoord;
+ FieldVector<albertCtype, dim> tmp2 = elMat_ * refcoord;
tmp2 += coord_(0);
for(int j=0; j<dim; j++)
@@ -764,7 +764,7 @@ namespace Dune
template<int dim, int dimworld>
inline bool AlbertGridElement <dim ,dimworld >::
- checkInside(const Vec<dim,albertCtype> &local)
+ checkInside(const FieldVector<albertCtype, dim> &local)
{
albertCtype sum = 0.0;
@@ -945,7 +945,7 @@ namespace Dune
}
template<int codim, int dim, int dimworld>
- inline Vec<dim,albertCtype>&
+ inline FieldVector<albertCtype, dim>&
AlbertGridEntity < codim, dim ,dimworld >::local()
{
return localFatherCoords_;
@@ -979,7 +979,7 @@ namespace Dune
inline albertCtype AlbertGridEntity < 0, dim ,dimworld >::
diam()
{
- Vec<dim> tmp;
+ FieldVector<albertCtype, dim> tmp;
return sqrt(geometry().integration_element(tmp));
}
@@ -994,8 +994,8 @@ namespace Dune
InterIt it = ibegin();
InterIt endit = iend ();
- Vec<dim> tmp;
- Vec<dim-1> tmp1;
+ FieldVector<albertCtype, dim> tmp;
+ FieldVector<albertCtype, dim-1> tmp1;
albertCtype vol = 0.5*geometry().integration_element(tmp);
albertCtype fak = 2.0;
@@ -1804,11 +1804,11 @@ namespace Dune
}
template< int dim, int dimworld>
- inline Vec<dimworld,albertCtype>& AlbertGridIntersectionIterator<dim,dimworld>::
- unit_outer_normal(Vec<dim-1,albertCtype>& local)
+ inline FieldVector<albertCtype, dimworld>& AlbertGridIntersectionIterator<dim,dimworld>::
+ unit_outer_normal(FieldVector<albertCtype, dim-1>& local)
{
// calculates the outer_normal
- Vec<dimworld,albertCtype>& tmp = outer_normal(local);
+ FieldVector<albertCtype, dimworld>& tmp = outer_normal(local);
double norm_1 = (1.0/tmp.norm2());
assert(norm_1 > 0.0);
@@ -1818,11 +1818,11 @@ namespace Dune
}
template< int dim, int dimworld>
- inline Vec<dimworld,albertCtype>& AlbertGridIntersectionIterator<dim,dimworld>::
+ inline FieldVector<albertCtype, dimworld>& AlbertGridIntersectionIterator<dim,dimworld>::
unit_outer_normal()
{
// calculates the outer_normal
- Vec<dimworld,albertCtype>& tmp = outer_normal();
+ FieldVector<albertCtype, dimworld>& tmp = outer_normal();
double norm_1 = (1.0/tmp.norm2());
assert(norm_1 > 0.0);
@@ -1832,8 +1832,8 @@ namespace Dune
}
template< int dim, int dimworld>
- inline Vec<dimworld,albertCtype>& AlbertGridIntersectionIterator<dim,dimworld>::
- outer_normal(Vec<dim-1,albertCtype>& local)
+ inline FieldVector<albertCtype, dimworld>& AlbertGridIntersectionIterator<dim,dimworld>::
+ outer_normal(FieldVector<albertCtype, dim-1>& local)
{
// we dont have curved boundary
// therefore return outer_normal
@@ -1841,7 +1841,7 @@ namespace Dune
}
template< int dim, int dimworld>
- inline Vec<dimworld,albertCtype>& AlbertGridIntersectionIterator<dim,dimworld>::
+ inline FieldVector<albertCtype, dimworld>& AlbertGridIntersectionIterator<dim,dimworld>::
outer_normal()
{
std::cout << "outer_normal() not correctly implemented yet! \n";
@@ -1852,7 +1852,7 @@ namespace Dune
}
template <>
- inline Vec<2,albertCtype>& AlbertGridIntersectionIterator<2,2>::
+ inline FieldVector<albertCtype, 2>& AlbertGridIntersectionIterator<2,2>::
outer_normal()
{
// seems to work
@@ -1865,7 +1865,7 @@ namespace Dune
}
template <>
- inline Vec<3,albertCtype>& AlbertGridIntersectionIterator<3,3>::
+ inline FieldVector<albertCtype, 3>& AlbertGridIntersectionIterator<3,3>::
outer_normal()
{
enum { dim = 3 };
@@ -3679,7 +3679,7 @@ namespace Dune
nb[elNum][i] = nit->index();
//std::cout << nb[elNum][i] << " Neigh \n";
- Vec<dimworld>& vec = (it->geometry())[i];
+ FieldVector<albertCtype, dimworld>& vec = (it->geometry())[i];
for (int j = 0; j < dimworld; j++)
coord[k][j] = vec(j);
@@ -3766,7 +3766,7 @@ namespace Dune
vertex[elNum][i] = k;
- Vec<dimworld> vec = (it->geometry())[i];
+ FieldVector<albertCtype, dimworld> vec = (it->geometry())[i];
for (int j = 0; j < dimworld; j++)
coord[k][j] = vec(j);
}
@@ -4489,9 +4489,9 @@ namespace Dune
template < int dim, int dimworld >
inline void AlbertGrid < dim, dimworld >::setNewCoords
- (const Vec<dimworld,albertCtype> & trans, const albertCtype scalar)
+ (const FieldVector<albertCtype, dimworld> & trans, const albertCtype scalar)
{
- static Vec<dimworld,albertCtype> trans_(0.0);
+ static FieldVector<albertCtype, dimworld> trans_(0.0);
static albertCtype scalar_ (1.0);
for(int i=0; i<macroVertices_.size(); i++)
diff --git a/grid/bsgrid.hh b/grid/bsgrid.hh
index dab04f5fb..8c5dd038a 100644
--- a/grid/bsgrid.hh
+++ b/grid/bsgrid.hh
@@ -77,7 +77,7 @@ namespace Dune
int corners ();
//! access to coordinates of corners. Index is the number of the corner
- Vec<dimworld,bs_ctype>& operator[] (int i);
+ FieldVector<bs_ctype, dimworld>& operator[] (int i);
/*! return reference element corresponding to this element. If this is
a reference element then self is returned.
@@ -86,14 +86,14 @@ namespace Dune
//! maps a local coordinate within reference element to
//! global coordinate in element
- Vec<dimworld,bs_ctype> global (const Vec<dim,bs_ctype>& local);
+ FieldVector<bs_ctype, dimworld> global (const FieldVector<bs_ctype, dim>& local);
//! maps a global coordinate within the element to a
//! local coordinate in its reference element
- Vec<dim,bs_ctype> local (const Vec<dimworld,bs_ctype>& global);
+ FieldVector<bs_ctype, dim> local (const FieldVector<bs_ctype, dimworld>& global);
//! returns true if the point in local coordinates is inside reference element
- bool checkInside(const Vec<dim,bs_ctype>& local);
+ bool checkInside(const FieldVector<bs_ctype, dim>& local);
/*!
Copy from grid.hh:
@@ -119,10 +119,10 @@ namespace Dune
*/
//! A(l)
- bs_ctype integration_element (const Vec<dim,bs_ctype>& local);
+ bs_ctype integration_element (const FieldVector<bs_ctype, dim>& local);
//! can only be called for dim=dimworld!
- Mat<dim,dim>& Jacobian_inverse (const Vec<dim,bs_ctype>& local);
+ Mat<dim,dim>& Jacobian_inverse (const FieldVector<bs_ctype, dim>& local);
//***********************************************************************
//! Methods that not belong to the Interface, but have to be public
@@ -158,7 +158,7 @@ namespace Dune
bs_ctype detDF_; //!< storage of integration_element
Mat<dimworld,dim,bs_ctype> A_; //!< for map2world NORMALE Nr.(ZEILE/SPALTE)
- Vec<dimworld,bs_ctype> tmpVec_;
+ FieldVector<bs_ctype, dimworld> tmpVec_;
};
@@ -208,7 +208,7 @@ namespace Dune
BSGridLevelIterator<0,dim,dimworld,All_Partition> father ();
//! local coordinates within father
- Vec<dim,bs_ctype>& local ();
+ FieldVector<bs_ctype, dim>& local ();
private:
BSGrid<dim,dimworld> &grid_;
@@ -512,10 +512,10 @@ namespace Dune
//! return unit outer normal, this should be dependent on local
//! coordinates for higher order boundary
- Vec<dimworld,bs_ctype>& unit_outer_normal (Vec<dim-1,bs_ctype>& local);
+ FieldVector<bs_ctype, dimworld>& unit_outer_normal (FieldVector<bs_ctype, dim-1>& local);
//! return unit outer normal, if you know it is constant use this function instead
- Vec<dimworld,bs_ctype>& unit_outer_normal ();
+ FieldVector<bs_ctype, dimworld>& unit_outer_normal ();
//! intersection of codimension 1 of this neighbor with element where
//! iteration started.
@@ -543,10 +543,10 @@ namespace Dune
//! return outer normal, this should be dependent on local
//! coordinates for higher order boundary
- Vec<dimworld,bs_ctype>& outer_normal (Vec<dim-1,bs_ctype>& local);
+ FieldVector<bs_ctype, dimworld>& outer_normal (FieldVector<bs_ctype, dim-1>& local);
//! return unit outer normal, if you know it is constant use this function instead
- Vec<dimworld,bs_ctype>& outer_normal ();
+ FieldVector<bs_ctype, dimworld>& outer_normal ();
private:
// if neighbour exists , do setup of new neighbour
@@ -569,8 +569,8 @@ namespace Dune
bool theSituation_; //! true if the "situation" occurs :-)
- Vec<dimworld,bs_ctype> outerNormal_; //! outerNormal ro current intersection
- Vec<dimworld,bs_ctype> unitOuterNormal_;
+ FieldVector<bs_ctype, dimworld> outerNormal_; //! outerNormal ro current intersection
+ FieldVector<bs_ctype, dimworld> unitOuterNormal_;
bool needSetup_; //! true if setup is needed
bool needNormal_;
diff --git a/grid/bsgrid/bsgrid.cc b/grid/bsgrid/bsgrid.cc
index f2b0843e6..0268f5b21 100644
--- a/grid/bsgrid/bsgrid.cc
+++ b/grid/bsgrid/bsgrid.cc
@@ -774,15 +774,15 @@ namespace Dune {
}
template< int dim, int dimworld>
- inline Vec<dimworld,bs_ctype>&
+ inline FieldVector<bs_ctype, dimworld>&
BSGridIntersectionIterator<dim,dimworld>::
- outer_normal(Vec<dim-1,bs_ctype>& local)
+ outer_normal(FieldVector<bs_ctype, dim-1>& local)
{
return this->outer_normal();
}
template< int dim, int dimworld>
- inline Vec<dimworld,bs_ctype>&
+ inline FieldVector<bs_ctype, dimworld>&
BSGridIntersectionIterator<dim,dimworld>::outer_normal()
{
assert(item_ != 0);
@@ -795,15 +795,15 @@ namespace Dune {
}
template< int dim, int dimworld>
- inline Vec<dimworld,bs_ctype>&
+ inline FieldVector<bs_ctype, dimworld>&
BSGridIntersectionIterator<dim,dimworld>::
- unit_outer_normal(Vec<dim-1,bs_ctype>& local)
+ unit_outer_normal(FieldVector<bs_ctype, dim-1>& local)
{
return this->unit_outer_normal();
}
template< int dim, int dimworld>
- inline Vec<dimworld,bs_ctype>&
+ inline FieldVector<bs_ctype, dimworld>&
BSGridIntersectionIterator<dim,dimworld>::unit_outer_normal()
{
unitOuterNormal_ = this->outer_normal();
@@ -1140,7 +1140,7 @@ namespace Dune {
}
template<int dim, int dimworld>
- inline Vec<dimworld,bs_ctype>& BSGridElement<dim,dimworld> :: operator[] (int i)
+ inline FieldVector<bs_ctype, dimworld>& BSGridElement<dim,dimworld> :: operator[] (int i)
{
assert((i>=0) && (i < dim+1));
return coord_(i);
@@ -1150,23 +1150,23 @@ namespace Dune {
// G L O B A L - - -
template<int dim, int dimworld> // dim = 1,2,3 dimworld = 3
- inline Vec<dimworld,bs_ctype> BSGridElement<dim,dimworld>::
- global(const Vec<dim,bs_ctype>& local)
+ inline FieldVector<bs_ctype, dimworld> BSGridElement<dim,dimworld>::
+ global(const FieldVector<bs_ctype, dim>& local)
{
if(!builtA_) calcElMatrix();
return (A_ * local) + coord_(0);
}
template<> // dim = dimworld = 3
- inline Vec<3,bs_ctype> BSGridElement<3,3>::
- local(const Vec<3,bs_ctype>& global)
+ inline FieldVector<bs_ctype, 3> BSGridElement<3,3>::
+ local(const FieldVector<bs_ctype, 3>& global)
{
if (!builtinverse_) buildJacobianInverse();
return Jinv_ * ( global - coord_(0));
}
template<int dim, int dimworld>
- inline bool BSGridElement<dim,dimworld> :: checkInside(const Vec<dim,bs_ctype>& local)
+ inline bool BSGridElement<dim,dimworld> :: checkInside(const FieldVector<bs_ctype, dim>& local)
{
bs_ctype sum = 0.0;
@@ -1192,7 +1192,7 @@ namespace Dune {
}
template<int dim, int dimworld>
- inline bs_ctype BSGridElement<dim,dimworld> :: integration_element (const Vec<dim,bs_ctype>& local)
+ inline bs_ctype BSGridElement<dim,dimworld> :: integration_element (const FieldVector<bs_ctype, dim>& local)
{
if(builtDetDF_)
return detDF_;
@@ -1207,7 +1207,7 @@ namespace Dune {
// J A C O B I A N _ I N V E R S E - - -
template<> // dim = dimworld = 3
- inline Mat<3,3>& BSGridElement<3,3> :: Jacobian_inverse (const Vec<3,bs_ctype>& local)
+ inline Mat<3,3>& BSGridElement<3,3> :: Jacobian_inverse (const FieldVector<bs_ctype, 3>& local)
{
if (!builtinverse_) buildJacobianInverse();
return Jinv_;
diff --git a/grid/bsgrid/bsinclude.hh b/grid/bsgrid/bsinclude.hh
index ed63b701a..ceb4bbc80 100644
--- a/grid/bsgrid/bsinclude.hh
+++ b/grid/bsgrid/bsinclude.hh
@@ -14,7 +14,7 @@
#define _BSGRID_USE_INDEX_
// define parameter type for normal method
-#define BSGridVecType Dune::Vec<3,double>
+#define BSGridVecType Dune::FieldVector<double, 3>
#include <dune/grid/common/indexstack.hh>
diff --git a/grid/common/grid.cc b/grid/common/grid.cc
index 1fc9e52b6..c5c098707 100644
--- a/grid/common/grid.cc
+++ b/grid/common/grid.cc
@@ -21,13 +21,13 @@ namespace Dune {
}
template<int dim, class ct>
- inline Vec<dim,ct>& ReferenceTopology<dim,ct>::center_codim0_local (int elemtype)
+ inline FieldVector<ct, dim>& ReferenceTopology<dim,ct>::center_codim0_local (int elemtype)
{
DUNE_THROW(GridError, "dimension too large");
}
template<int dim, class ct>
- inline Vec<dim-1,ct>& ReferenceTopology<dim,ct>::center_codim1_local (int elemtype, int i)
+ inline FieldVector<ct, dim-1>& ReferenceTopology<dim,ct>::center_codim1_local (int elemtype, int i)
{
DUNE_THROW(GridError, "dimension too large");
}
@@ -36,17 +36,17 @@ namespace Dune {
template<class ct>
inline ReferenceTopology<1,ct>::ReferenceTopology ()
{
- center0_local[0] = Vec<1,ct>(0.5);
+ center0_local[0] = FieldVector<ct, 1>(0.5);
}
template<class ct>
- inline Vec<1,ct>& ReferenceTopology<1,ct>::center_codim0_local (int elemtype)
+ inline FieldVector<ct, 1>& ReferenceTopology<1,ct>::center_codim0_local (int elemtype)
{
return center0_local[0];
}
template<class ct>
- inline Vec<0,ct>& ReferenceTopology<1,ct>::center_codim1_local (int elemtype, int i)
+ inline FieldVector<ct, 0>& ReferenceTopology<1,ct>::center_codim1_local (int elemtype, int i)
{
return center1_local[0];
}
@@ -55,19 +55,19 @@ namespace Dune {
template<class ct>
inline ReferenceTopology<2,ct>::ReferenceTopology ()
{
- center0_local[0] = Vec<2,ct>(1.0/3.0);
- center0_local[1] = Vec<2,ct>(0.5);
- center1_local[0] = Vec<1,ct>(0.5);
+ center0_local[0] = FieldVector<ct, 2>(1.0/3.0);
+ center0_local[1] = FieldVector<ct, 2>(0.5);
+ center1_local[0] = FieldVector<ct, 1>(0.5);
}
template<class ct>
- inline Vec<2,ct>& ReferenceTopology<2,ct>::center_codim0_local (int elemtype)
+ inline FieldVector<ct, 2>& ReferenceTopology<2,ct>::center_codim0_local (int elemtype)
{
return center0_local[elemtype-2];
}
template<class ct>
- inline Vec<1,ct>& ReferenceTopology<2,ct>::center_codim1_local (int elemtype, int i)
+ inline FieldVector<ct, 1>& ReferenceTopology<2,ct>::center_codim1_local (int elemtype, int i)
{
return center1_local[0];
}
@@ -76,24 +76,24 @@ namespace Dune {
template<class ct>
inline ReferenceTopology<3,ct>::ReferenceTopology ()
{
- center0_local[0] = Vec<3,ct>(0.25);
- center0_local[1] = Vec<3,ct>(0.0); // pyramid is missing !
- center0_local[2] = Vec<3,ct>(0.0); // prism is missing
- center0_local[3] = Vec<3,ct>(0.5);
- for (int i=0; i<6; i++) center1_local[0][i] = Vec<2,ct>(1.0/3.0);
- for (int i=0; i<6; i++) center1_local[1][i] = Vec<2,ct>(0.0);
- for (int i=0; i<6; i++) center1_local[2][i] = Vec<2,ct>(0.0);
- for (int i=0; i<6; i++) center1_local[3][i] = Vec<2,ct>(0.5);
+ center0_local[0] = FieldVector<ct, 3>(0.25);
+ center0_local[1] = FieldVector<ct, 3>(0.0); // pyramid is missing !
+ center0_local[2] = FieldVector<ct, 3>(0.0); // prism is missing
+ center0_local[3] = FieldVector<ct, 3>(0.5);
+ for (int i=0; i<6; i++) center1_local[0][i] = FieldVector<ct, 2>(1.0/3.0);
+ for (int i=0; i<6; i++) center1_local[1][i] = FieldVector<ct, 2>(0.0);
+ for (int i=0; i<6; i++) center1_local[2][i] = FieldVector<ct, 2>(0.0);
+ for (int i=0; i<6; i++) center1_local[3][i] = FieldVector<ct, 2>(0.5);
}
template<class ct>
- inline Vec<3,ct>& ReferenceTopology<3,ct>::center_codim0_local (int elemtype)
+ inline FieldVector<ct, 3>& ReferenceTopology<3,ct>::center_codim0_local (int elemtype)
{
return center0_local[elemtype-4];
}
template<class ct>
- inline Vec<2,ct>& ReferenceTopology<3,ct>::center_codim1_local (int elemtype, int i)
+ inline FieldVector<ct, 2>& ReferenceTopology<3,ct>::center_codim1_local (int elemtype, int i)
{
return center1_local[elemtype-4][i];
}
@@ -118,7 +118,7 @@ namespace Dune {
}
template<int dim, int dimworld, class ct,template<int,int> class ElementImp>
- inline Vec<dimworld,ct>& Element<dim,dimworld,ct,ElementImp>::operator[] (int i)
+ inline FieldVector<ct, dimworld>& Element<dim,dimworld,ct,ElementImp>::operator[] (int i)
{
return asImp().operator[](i);
}
@@ -130,31 +130,31 @@ namespace Dune {
}
template<int dim, int dimworld, class ct,template<int,int> class ElementImp>
- inline Vec<dimworld,ct> Element<dim,dimworld,ct,ElementImp>::global (const Vec<dim,ct>& local)
+ inline FieldVector<ct, dimworld> Element<dim,dimworld,ct,ElementImp>::global (const FieldVector<ct, dim>& local)
{
return asImp().global(local);
}
template<int dim, int dimworld, class ct,template<int,int> class ElementImp>
- inline Vec<dim,ct> Element<dim,dimworld,ct,ElementImp>::local (const Vec<dimworld,ct>& global)
+ inline FieldVector<ct, dim> Element<dim,dimworld,ct,ElementImp>::local (const FieldVector<ct, dimworld>& global)
{
return asImp().local(global);
}
template<int dim, int dimworld, class ct,template<int,int> class ElementImp>
- inline bool Element<dim,dimworld,ct,ElementImp>::checkInside (const Vec<dim,ct>& local)
+ inline bool Element<dim,dimworld,ct,ElementImp>::checkInside (const FieldVector<ct, dim>& local)
{
return asImp().checkInside(local);
}
template<int dim, int dimworld, class ct,template<int,int> class ElementImp>
- inline ct Element<dim,dimworld,ct,ElementImp>::integration_element (const Vec<dim,ct>& local)
+ inline ct Element<dim,dimworld,ct,ElementImp>::integration_element (const FieldVector<ct, dim>& local)
{
return asImp().integration_element(local);
}
template<int dim, int dimworld, class ct,template<int,int> class ElementImp>
- inline Mat<dim,dim>& Element<dim,dimworld,ct,ElementImp>::Jacobian_inverse (const Vec<dim,ct>& local)
+ inline Mat<dim,dim>& Element<dim,dimworld,ct,ElementImp>::Jacobian_inverse (const FieldVector<ct, dim>& local)
{
return asImp().Jacobian_inverse(local);
}
@@ -165,13 +165,13 @@ namespace Dune {
// simply call all members, this forces compiler to compile them ...
type();
- Vec<dim,ct> l;
- Vec<dimworld,ct> g;
+ FieldVector<ct, dim> l;
+ FieldVector<ct, dimworld> g;
// check the methods local,global, refelem and checkInside
for(int i=0; i<corners(); i++)
{
- Vec<dimworld,ct> & coord = operator[] (i);
+ FieldVector<ct, dimworld> & coord = operator[] (i);
l = local ( coord );
if( !checkInside (l) )
{
@@ -219,7 +219,7 @@ namespace Dune {
}
template<int dimworld, class ct,template<int,int> class ElementImp>
- inline Vec<dimworld,ct>& Element<0,dimworld,ct,ElementImp>::operator[] (int i)
+ inline FieldVector<ct, dimworld>& Element<0,dimworld,ct,ElementImp>::operator[] (int i)
{
return asImp().operator[](i);
}
@@ -314,7 +314,7 @@ namespace Dune {
template<int,int> class ElementImp
, template<int,int> class BoundaryEntityImp
>
- inline Vec<dimworld,ct>& IntersectionIterator<dim,dimworld,ct,IntersectionIteratorImp,EntityImp,ElementImp,BoundaryEntityImp>::unit_outer_normal (Vec<dim-1,ct>& local)
+ inline FieldVector<ct, dimworld>& IntersectionIterator<dim,dimworld,ct,IntersectionIteratorImp,EntityImp,ElementImp,BoundaryEntityImp>::unit_outer_normal (FieldVector<ct, dim-1>& local)
{
return asImp().unit_outer_normal(local);
}
@@ -325,7 +325,7 @@ namespace Dune {
template<int,int> class ElementImp
, template<int,int> class BoundaryEntityImp
>
- inline Vec<dimworld,ct>& IntersectionIterator<dim,dimworld,ct,IntersectionIteratorImp,EntityImp,ElementImp,BoundaryEntityImp>::unit_outer_normal ()
+ inline FieldVector<ct, dimworld>& IntersectionIterator<dim,dimworld,ct,IntersectionIteratorImp,EntityImp,ElementImp,BoundaryEntityImp>::unit_outer_normal ()
{
return asImp().unit_outer_normal();
}
@@ -410,7 +410,7 @@ namespace Dune {
operator*();
operator->();
boundary();
- Vec<dim-1,ct> l;
+ FieldVector<ct, dim-1> l;
unit_outer_normal(l);
unit_outer_normal();
intersection_self_local();
@@ -430,7 +430,7 @@ namespace Dune {
template<int,int> class ElementImp
, template<int,int> class BoundaryEntityImp
>
- inline Vec<dimworld,ct>& IntersectionIteratorDefault<dim,dimworld,ct,IntersectionIteratorImp,EntityImp,ElementImp,BoundaryEntityImp>::outer_normal ()
+ inline FieldVector<ct, dimworld>& IntersectionIteratorDefault<dim,dimworld,ct,IntersectionIteratorImp,EntityImp,ElementImp,BoundaryEntityImp>::outer_normal ()
{
// make a copy, is nessasary
outerNormal_ = asImp().unit_outer_normal();
@@ -444,8 +444,8 @@ namespace Dune {
template<int,int> class ElementImp
, template<int,int> class BoundaryEntityImp
>
- inline Vec<dimworld,ct>& IntersectionIteratorDefault<dim,dimworld,ct,IntersectionIteratorImp,EntityImp,ElementImp,BoundaryEntityImp>::
- outer_normal (Vec<dim-1,ct>& local)
+ inline FieldVector<ct, dimworld>& IntersectionIteratorDefault<dim,dimworld,ct,IntersectionIteratorImp,EntityImp,ElementImp,BoundaryEntityImp>::
+ outer_normal (FieldVector<ct, dim-1>& local)
{
// make a copy, is nessasary
outerNormal_ = asImp().unit_outer_normal(local);
@@ -828,7 +828,7 @@ namespace Dune {
template<int,int> class IntersectionIteratorImp,
template<int,int> class HierarchicIteratorImp
>
- inline Vec<dim,ct>& Entity<dim,dim,dimworld,ct,EntityImp,ElementImp,LevelIteratorImp,IntersectionIteratorImp,HierarchicIteratorImp>::local ()
+ inline FieldVector<ct, dim>& Entity<dim,dim,dimworld,ct,EntityImp,ElementImp,LevelIteratorImp,IntersectionIteratorImp,HierarchicIteratorImp>::local ()
{
return asImp().local();
}
diff --git a/grid/common/grid.hh b/grid/common/grid.hh
index 7aa0c2b3b..505bbbffc 100644
--- a/grid/common/grid.hh
+++ b/grid/common/grid.hh
@@ -168,10 +168,10 @@ namespace Dune {
ReferenceTopology ();
//! return local coordinates of center in reference element
- Vec<dim,ct>& center_codim0_local (int elemtype);
+ FieldVector<ct, dim>& center_codim0_local (int elemtype);
//! return local coordinates of center of ith codim 1 subentity
- Vec<dim-1,ct>& center_codim1_local (int elemtype, int i);
+ FieldVector<ct, dim-1>& center_codim1_local (int elemtype, int i);
};
// Specialization dim=1
@@ -182,14 +182,14 @@ namespace Dune {
ReferenceTopology ();
//! return local coordinates of center in reference element
- Vec<1,ct>& center_codim0_local (int elemtype);
+ FieldVector<ct, 1>& center_codim0_local (int elemtype);
//! return local coordinates of center of ith codim 1 subentity
- Vec<0,ct>& center_codim1_local (int elemtype, int i);
+ FieldVector<ct, 0>& center_codim1_local (int elemtype, int i);
private:
- Vec<1,ct> center0_local[1]; // ElementType
- Vec<0,ct> center1_local[1]; // ElementType x faces
+ FieldVector<ct, 1> center0_local[1]; // ElementType
+ FieldVector<ct, 0> center1_local[1]; // ElementType x faces
};
// Specialization dim=2
@@ -200,14 +200,14 @@ namespace Dune {
ReferenceTopology ();
//! return local coordinates of center in reference element
- Vec<2,ct>& center_codim0_local (int elemtype);
+ FieldVector<ct, 2>& center_codim0_local (int elemtype);
//! return local coordinates of center of ith codim 1 subentity
- Vec<1,ct>& center_codim1_local (int elemtype, int i);
+ FieldVector<ct, 1>& center_codim1_local (int elemtype, int i);
private:
- Vec<2,ct> center0_local[2]; // ElementType
- Vec<1,ct> center1_local[2]; // ElementType x faces
+ FieldVector<ct, 2> center0_local[2]; // ElementType
+ FieldVector<ct, 1> center1_local[2]; // ElementType x faces
};
@@ -219,14 +219,14 @@ namespace Dune {
ReferenceTopology ();
//! return local coordinates of center in reference element
- Vec<3,ct>& center_codim0_local (int elemtype);
+ FieldVector<ct, 3>& center_codim0_local (int elemtype);
//! return local coordinates of center of ith codim 1 subentity
- Vec<2,ct>& center_codim1_local (int elemtype, int i);
+ FieldVector<ct, 2>& center_codim1_local (int elemtype, int i);
private:
- Vec<3,ct> center0_local[4]; // ElementType
- Vec<2,ct> center1_local[4][6]; // ElementType x faces
+ FieldVector<ct, 3> center0_local[4]; // ElementType
+ FieldVector<ct, 2> center1_local[4][6]; // ElementType x faces
};
@@ -292,7 +292,7 @@ namespace Dune {
int corners ();
//! access to coordinates of corners. Index is the number of the corner
- Vec<dimworld,ct>& operator[] (int i);
+ FieldVector<ct, dimworld>& operator[] (int i);
/*! return reference element corresponding to this element. If this is
a reference element then self is returned. A reference to a reference
@@ -302,13 +302,13 @@ namespace Dune {
ElementImp<dim,dim>& refelem ();
//! maps a local coordinate within reference element to global coordinate in element
- Vec<dimworld,ct> global (const Vec<dim,ct>& local);
+ FieldVector<ct, dimworld> global (const FieldVector<ct, dim>& local);
//! maps a global coordinate within the element to a local coordinate in its reference element
- Vec<dim,ct> local (const Vec<dimworld,ct>& global);
+ FieldVector<ct, dim> local (const FieldVector<ct, dimworld>& global);
//! return true if the point in local coordinates lies inside the reference element
- bool checkInside (const Vec<dim,ct>& local);
+ bool checkInside (const FieldVector<ct, dim>& local);
/*! Integration over a general element is done by integrating over the reference element
and using the transformation from the reference element to the global element as follows:
@@ -329,10 +329,10 @@ namespace Dune {
will directly translate in substantial savings in the computation of finite element
stiffness matrices.
*/
- ct integration_element (const Vec<dim,ct>& local);
+ ct integration_element (const FieldVector<ct, dim>& local);
//! can only be called for dim=dimworld!
- Mat<dim,dim>& Jacobian_inverse (const Vec<dim,ct>& local);
+ Mat<dim,dim>& Jacobian_inverse (const FieldVector<ct, dim>& local);
/*! \internal
Checking presence and format of all interface functions. With
@@ -391,7 +391,7 @@ namespace Dune {
int corners ();
//! access to coordinates of corners. Index is the number of the corner
- Vec<dimworld,ct>& operator[] (int i);
+ FieldVector<ct, dimworld>& operator[] (int i);
/*! \internal Checking presence and format of all interface functions. With
this method all derived classes can check their correct definition.
@@ -475,7 +475,7 @@ namespace Dune {
ElementImp<dim,dimworld> & geometry ();
//! return barycenter of ghostcell
- Vec<dimworld,ct>& outerPoint ();
+ FieldVector<ct, dimworld>& outerPoint ();
private:
//! Barton-Nackman trick
@@ -568,10 +568,10 @@ namespace Dune {
BoundaryEntityImp<dim,dimworld> & boundaryEntity ();
//! return unit outer normal, this should be dependent on local coordinates for higher order boundary
- Vec<dimworld,ct>& unit_outer_normal (Vec<dim-1,ct>& local);
+ FieldVector<ct, dimworld>& unit_outer_normal (FieldVector<ct, dim-1>& local);
//! return unit outer normal, if you know it is constant use this function instead
- Vec<dimworld,ct>& unit_outer_normal ();
+ FieldVector<ct, dimworld>& unit_outer_normal ();
/*! intersection of codimension 1 of this neighbor with element where iteration started.
Here returned element is in LOCAL coordinates of the element where iteration started.
@@ -631,19 +631,19 @@ namespace Dune {
public:
//! return outer normal, which is the unit_outer_normal() scaled with the
//! volume of the intersection_self_global ()
- Vec<dimworld,ct>& outer_normal ();
+ FieldVector<ct, dimworld>& outer_normal ();
//! return unit outer normal, this should be dependent on
//! local coordinates for higher order boundary
//! the normal is scaled with the integration element
- Vec<dimworld,ct>& outer_normal (Vec<dim-1,ct>& local);
+ FieldVector<ct, dimworld>& outer_normal (FieldVector<ct, dim-1>& local);
protected:
//! the outer normal
- Vec<dimworld,ct> outerNormal_;
+ FieldVector<ct, dimworld> outerNormal_;
//! tmp Vec for integration_element
- Vec<dim-1,ct> tmp_;
+ FieldVector<ct, dim-1> tmp_;
private:
//! Barton-Nackman trick
@@ -1104,7 +1104,7 @@ namespace Dune {
LevelIteratorImp<0,dim,dimworld,All_Partition> father ();
//! local coordinates within father
- Vec<dim,ct>& local ();
+ FieldVector<ct, dim>& local ();
/*! \internal Checking presence and format of all interface functions. With
this method all derived classes can check their correct definition.
diff --git a/grid/sgrid.hh b/grid/sgrid.hh
index f0458c86b..0556f5ad7 100644
--- a/grid/sgrid.hh
+++ b/grid/sgrid.hh
@@ -116,7 +116,7 @@ namespace Dune {
int corners ();
//! access to coordinates of corners. Index is the number of the corner
- Vec<dimworld,sgrid_ctype>& operator[] (int i);
+ FieldVector<sgrid_ctype, dimworld>& operator[] (int i);
/*! return reference element corresponding to this element. If this is
a reference element then self is returned. A reference to a reference
@@ -126,13 +126,13 @@ namespace Dune {
static SElement<dim,dim>& refelem ();
//! maps a local coordinate within reference element to global coordinate in element
- Vec<dimworld,sgrid_ctype> global (const Vec<dim,sgrid_ctype>& local);
+ FieldVector<sgrid_ctype, dimworld> global (const FieldVector<sgrid_ctype, dim>& local);
//! maps a global coordinate within the element to a local coordinate in its reference element
- Vec<dim,sgrid_ctype> local (const Vec<dimworld,sgrid_ctype>& global);
+ FieldVector<sgrid_ctype, dim> local (const FieldVector<sgrid_ctype, dimworld>& global);
//! returns true if the point in local coordinates is located within the refelem
- bool checkInside (const Vec<dim,sgrid_ctype>& local);
+ bool checkInside (const FieldVector<sgrid_ctype, dim>& local);
/*! Integration over a general element is done by integrating over the reference element
and using the transformation from the reference element to the global element as follows:
@@ -153,10 +153,10 @@ namespace Dune {
will directly translate in substantial savings in the computation of finite element
stiffness matrices.
*/
- sgrid_ctype integration_element (const Vec<dim,sgrid_ctype>& local);
+ sgrid_ctype integration_element (const FieldVector<sgrid_ctype, dim>& local);
//! can only be called for dim=dimworld!
- Mat<dim,dim>& Jacobian_inverse (const Vec<dim,sgrid_ctype>& local);
+ Mat<dim,dim>& Jacobian_inverse (const FieldVector<sgrid_ctype, dim>& local);
//! print internal data
void print (std::ostream& ss, int indent);
@@ -171,9 +171,9 @@ namespace Dune {
SElement (bool b);
private:
- Vec<dimworld,sgrid_ctype> s; //!< position of element
+ FieldVector<sgrid_ctype, dimworld> s; //!< position of element
Mat<dimworld,dim,sgrid_ctype> A; //!< direction vectors as matrix
- Vec<dimworld,sgrid_ctype> c[1<<dim]; //!< coordinate vectors of corners
+ FieldVector<sgrid_ctype, dimworld> c[1<<dim]; //!< coordinate vectors of corners
Mat<dim,dim,sgrid_ctype> Jinv; //!< storage for inverse of jacobian
bool builtinverse;
};
@@ -193,7 +193,7 @@ namespace Dune {
int corners ();
//! access to coordinates of corners. Index is the number of the corner
- Vec<dimworld,sgrid_ctype>& operator[] (int i);
+ FieldVector<sgrid_ctype, dimworld>& operator[] (int i);
//! print internal data
void print (std::ostream& ss, int indent);
@@ -205,7 +205,7 @@ namespace Dune {
SElement (bool b);
protected:
- Vec<dimworld,sgrid_ctype> s; //!< position of element
+ FieldVector<sgrid_ctype, dimworld> s; //!< position of element
};
template <int dim, int dimworld>
@@ -233,9 +233,9 @@ namespace Dune {
SElement<dim,dimworld> & geometry() { return elem_; };
//! return outer barycenter of ghost cell
- Vec<dimworld,sgrid_ctype> & outerPoint () { return outerPoint_; };
+ FieldVector<sgrid_ctype, dimworld> & outerPoint () { return outerPoint_; };
private:
- Vec<dimworld> outerPoint_;
+ FieldVector<sgrid_ctype, dimworld> outerPoint_;
SElement<dim,dimworld> elem_;
};
@@ -278,10 +278,10 @@ namespace Dune {
SEntity<0,dim,dimworld>* operator->();
//! return unit outer normal, this should be dependent on local coordinates for higher order boundary
- Vec<dimworld,sgrid_ctype>& unit_outer_normal (Vec<dim-1,sgrid_ctype>& local);
+ FieldVector<sgrid_ctype, dimworld>& unit_outer_normal (FieldVector<sgrid_ctype, dim-1>& local);
//! return unit outer normal, if you know it is constant use this function instead
- Vec<dimworld,sgrid_ctype>& unit_outer_normal ();
+ FieldVector<sgrid_ctype, dimworld>& unit_outer_normal ();
/*! intersection of codimension 1 of this neighbor with element where iteration started.
Here returned element is in LOCAL coordinates of the element where iteration started.
@@ -326,7 +326,7 @@ namespace Dune {
bool valid_count; //!< true if count is in range
bool is_on_boundary; //!< true if neighbor is otside the domain
SEntity<0,dim,dimworld> e; //!< virtual neighbor entity
- Vec<dimworld,sgrid_ctype> normal; //!< outer unit normal direction
+ FieldVector<sgrid_ctype, dimworld> normal; //!< outer unit normal direction
bool built_intersections; //!< true if all intersections have been built
SElement<dim-1,dim> is_self_local; //!< intersection in own local coordinates
SElement<dim-1,dimworld> is_global; //!< intersection in global coordinates, map consistent with is_self_local
@@ -631,7 +631,7 @@ namespace Dune {
SLevelIterator<0,dim,dimworld,All_Partition> father ();
//! local coordinates within father
- Vec<dim,sgrid_ctype>& local ();
+ FieldVector<sgrid_ctype, dim>& local ();
// members specific to SEntity
//! constructor
@@ -650,7 +650,7 @@ namespace Dune {
private:
bool built_father;
int father_id;
- Vec<dim,sgrid_ctype> in_father_local;
+ FieldVector<sgrid_ctype, dim> in_father_local;
void make_father();
};
@@ -815,8 +815,23 @@ namespace Dune {
*/
void globalRefine (int refCount);
+ /** \brief Get number of elements in each coordinate direction */
+ const FixedArray<int, dim>& dims(int level) const {
+ return N[level];
+ }
+
+ /** \brief Get lower left corner */
+ const FieldVector<sgrid_ctype, dimworld>& lowerLeft() const {
+ return low;
+ }
+
+ /** \brief Get upper right corner */
+ FieldVector<sgrid_ctype, dimworld> upperRight() const {
+ return low+H;
+ }
+
//! map expanded coordinates to position
- Vec<dim,sgrid_ctype> pos (int level, FixedArray<int,dim>& z);
+ FieldVector<sgrid_ctype, dim> pos (int level, FixedArray<int,dim>& z);
//! compute codim from coordinate
int codim (int level, FixedArray<int,dim>& z);
@@ -846,10 +861,10 @@ namespace Dune {
void makeSGrid (const int* N_, const sgrid_ctype* L_, const sgrid_ctype* H_);
int L; // number of levels in hierarchic mesh 0<=level<L
- FixedArray<sgrid_ctype,dim> low; // lower left corner of the grid
- FixedArray<sgrid_ctype,dim> H; // length of cube per direction
+ FieldVector<sgrid_ctype, dim> low; // lower left corner of the grid
+ FieldVector<sgrid_ctype, dim> H; // length of cube per direction
FixedArray<int,dim> N[MAXL]; // number of elements per direction
- Vec<dim,sgrid_ctype> h[MAXL]; // mesh size per direction
+ FieldVector<sgrid_ctype, dim> h[MAXL]; // mesh size per direction
CubeMapper<dim> mapper[MAXL]; // a mapper for each level
// faster implemantation od subIndex
diff --git a/grid/sgrid/numbering.hh b/grid/sgrid/numbering.hh
index 73abc0602..32f63cef4 100644
--- a/grid/sgrid/numbering.hh
+++ b/grid/sgrid/numbering.hh
@@ -1,7 +1,9 @@
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
-#ifndef __NUMBERING_HH__
-#define __NUMBERING_HH__
+#ifndef __DUNE_NUMBERING_HH__
+#define __DUNE_NUMBERING_HH__
+
+#include <dune/common/fixedarray.hh>
namespace Dune {
diff --git a/grid/sgrid/sgrid.cc b/grid/sgrid/sgrid.cc
index 542207ac6..0ff2f9a79 100644
--- a/grid/sgrid/sgrid.cc
+++ b/grid/sgrid/sgrid.cc
@@ -33,7 +33,11 @@ namespace Dune {
// copy arguments
s = 0.0;
- for (int j=0; j<dim; j++) A(j) = Vec<dimworld,sgrid_ctype>(j,1.0); // make unit vectors
+ for (int j=0; j<dim; j++) {
+ // make unit vectors
+ A(j) = FieldVector<sgrid_ctype, dimworld>(0.0);
+ A(j,j) = 1.0;
+ }
// make corners
for (int i=0; i<(1<<dim); i++) // there are 2^d corners
@@ -94,29 +98,29 @@ namespace Dune {
}
template<int dim, int dimworld>
- inline Vec<dimworld,sgrid_ctype>& SElement<dim,dimworld>::operator[] (int i)
+ inline FieldVector<sgrid_ctype, dimworld>& SElement<dim,dimworld>::operator[] (int i)
{
return c[i];
}
template<int dim, int dimworld>
- inline Vec<dimworld,sgrid_ctype> SElement<dim,dimworld>::global (const Vec<dim,sgrid_ctype>& local)
+ inline FieldVector<sgrid_ctype, dimworld> SElement<dim,dimworld>::global (const FieldVector<sgrid_ctype, dim>& local)
{
return s+(A*local);
}
template<int dim, int dimworld>
- inline Vec<dim,sgrid_ctype> SElement<dim,dimworld>::local (const Vec<dimworld,sgrid_ctype>& global)
+ inline FieldVector<sgrid_ctype, dim> SElement<dim,dimworld>::local (const FieldVector<sgrid_ctype, dimworld>& global)
{
- Vec<dim,sgrid_ctype> l; // result
- Vec<dimworld,sgrid_ctype> rhs = global-s;
+ FieldVector<sgrid_ctype, dim> l; // result
+ FieldVector<sgrid_ctype, dimworld> rhs = global-s;
for (int k=0; k<dim; k++)
- l(k) = (rhs*A(k)) / (A(k)*A(k));
+ l[k] = (rhs*A(k)) / (A(k)*A(k));
return l;
}
template<int dim, int dimworld>
- inline bool SElement<dim,dimworld>::checkInside (const Vec<dim,sgrid_ctype>& local)
+ inline bool SElement<dim,dimworld>::checkInside (const FieldVector<sgrid_ctype, dim>& local)
{
// check wether they are in the reference element
for(int i=0; i<dim; i++)
@@ -128,16 +132,16 @@ namespace Dune {
}
template<int dim, int dimworld>
- inline sgrid_ctype SElement<dim,dimworld>::integration_element (const Vec<dim,sgrid_ctype>& local)
+ inline sgrid_ctype SElement<dim,dimworld>::integration_element (const FieldVector<sgrid_ctype, dim>& local)
{
sgrid_ctype s = 1.0;
- for (int j=0; j<dim; j++) s *= A(j).norm1();
+ for (int j=0; j<dim; j++) s *= A(j).one_norm();
return s;
}
template<int dim, int dimworld>
- inline Mat<dim,dim>& SElement<dim,dimworld>::Jacobian_inverse (const Vec<dim,sgrid_ctype>& local)
+ inline Mat<dim,dim>& SElement<dim,dimworld>::Jacobian_inverse (const FieldVector<sgrid_ctype, dim>& local)
{
assert(dim==dimworld);
@@ -204,7 +208,7 @@ namespace Dune {
}
template<int dimworld>
- inline Vec<dimworld,sgrid_ctype>& SElement<0,dimworld>::operator[] (int i)
+ inline FieldVector<sgrid_ctype, dimworld>& SElement<0,dimworld>::operator[] (int i)
{
return s;
}
@@ -295,7 +299,7 @@ namespace Dune {
// count number of direction vectors found
int dir=0;
- Vec<dim,sgrid_ctype> p1,p2;
+ FieldVector<sgrid_ctype, dim> p1,p2;
FixedArray<int,dim> t=z;
// check all directions
@@ -425,7 +429,7 @@ namespace Dune {
FixedArray<int,dim> zz = this->grid->compress(this->l,this->z);
// look for odd coordinates
- Vec<dim,sgrid_ctype> delta;
+ FieldVector<sgrid_ctype, dim> delta;
for (int i=0; i<dim; i++)
if (zz[i]%2)
{
@@ -447,7 +451,7 @@ namespace Dune {
// now make a subcube of size 1/2 in each direction
Mat<dim,dim+1,sgrid_ctype> As;
- Vec<dim,sgrid_ctype> v;
+ FieldVector<sgrid_ctype, dim> v;
for (int i=0; i<dim; i++)
{
v = 0.0; v(i) = 0.5;
@@ -525,7 +529,7 @@ namespace Dune {
FixedArray<int,dim> zz = this->grid->compress(this->l,this->z);
// to find father, make all coordinates odd
- Vec<dim,sgrid_ctype> delta;
+ FieldVector<sgrid_ctype, dim> delta;
for (int i=0; i<dim; i++)
if (zz[i]%2)
{
@@ -568,7 +572,7 @@ namespace Dune {
}
template<int dim, int dimworld>
- inline Vec<dim,sgrid_ctype>& SEntity<dim,dim,dimworld>::local ()
+ inline FieldVector<sgrid_ctype, dim>& SEntity<dim,dim,dimworld>::local ()
{
if (!built_father) make_father();
return in_father_local;
@@ -702,9 +706,9 @@ namespace Dune {
// while we are at it, compute normal direction
normal = 0.0;
if (count%2)
- normal(count/2) = 1.0; // odd
+ normal[count/2] = 1.0; // odd
else
- normal(count/2) = -1.0; // even
+ normal[count/2] = -1.0; // even
// now check if neighbor exists
is_on_boundary = !grid->exists(self->level(),zrednb);
@@ -819,7 +823,7 @@ namespace Dune {
// z1 is even in direction dir, all others must be odd because it is codim 1
Mat<dim,dim,sgrid_ctype> As;
- Vec<dim,sgrid_ctype> p1,p2;
+ FieldVector<sgrid_ctype, dim> p1,p2;
int t;
// local coordinates in self
@@ -918,14 +922,14 @@ namespace Dune {
}
template<int dim, int dimworld>
- inline Vec<dimworld,sgrid_ctype>&
- SIntersectionIterator<dim,dimworld>::unit_outer_normal (Vec<dim-1,sgrid_ctype>& local)
+ inline FieldVector<sgrid_ctype, dimworld>&
+ SIntersectionIterator<dim,dimworld>::unit_outer_normal (FieldVector<sgrid_ctype, dim-1>& local)
{
return normal;
}
template<int dim, int dimworld>
- inline Vec<dimworld,sgrid_ctype>&
+ inline FieldVector<sgrid_ctype, dimworld>&
SIntersectionIterator<dim,dimworld>::unit_outer_normal ()
{
return normal;
@@ -994,7 +998,7 @@ namespace Dune {
// define coarse mesh
mapper[0].make(N[0]);
for (int i=0; i<dim; i++)
- h[0](i) = (H[i]-low[i])/((sgrid_ctype)N[0][i]);
+ h[0][i] = (H[i]-low[i])/((sgrid_ctype)N[0][i]);
std::cout << "level=" << L-1 << " size=(" << N[L-1][0];
for (int i=1; i<dim; i++) std::cout << "," << N[L-1][i];
@@ -1038,7 +1042,8 @@ namespace Dune {
mapper[L].make(N[L]);
// compute mesh size
- for (int i=0; i<dim; i++) h[L](i) = (H[i]-low[i])/((sgrid_ctype)N[L][i]);
+ for (int i=0; i<dim; i++)
+ h[L][i] = (H[i]-low[i])/((sgrid_ctype)N[L][i]);
L++;
std::cout << "level=" << L-1 << " size=(" << N[L-1][0];
@@ -1166,11 +1171,11 @@ namespace Dune {
}
template<int dim, int dimworld>
- inline Vec<dim,sgrid_ctype> SGrid<dim,dimworld>::pos (int level, FixedArray<int,dim>& z)
+ inline FieldVector<sgrid_ctype, dim> SGrid<dim,dimworld>::pos (int level, FixedArray<int,dim>& z)
{
- Vec<dim,sgrid_ctype> x;
+ FieldVector<sgrid_ctype, dim> x;
for (int k=0; k<dim; k++)
- x(k) = (z[k]*h[level](k))*0.5 + low[k];
+ x[k] = (z[k]*h[level][k])*0.5 + low[k];
return x;
}
diff --git a/grid/simplegrid.hh b/grid/simplegrid.hh
index 303e6ab61..0159980ef 100644
--- a/grid/simplegrid.hh
+++ b/grid/simplegrid.hh
@@ -144,7 +144,7 @@ namespace Dune {
}
//! access to coordinates of corners. Index is the number of the corner
- Vec<dimworld,simplegrid_ctype>& operator[] (int i)
+ FieldVector<simplegrid_ctype, dimworld>& operator[] (int i)
{
int lk=0;
@@ -172,9 +172,9 @@ namespace Dune {
}
//! maps a local coordinate within reference element to global coordinate in element
- Vec<dimworld,simplegrid_ctype> global (const Vec<dim,simplegrid_ctype>& local)
+ FieldVector<simplegrid_ctype, dimworld> global (const FieldVector<simplegrid_ctype, dim>& local)
{
- Vec<dimworld,simplegrid_ctype> g;
+ FieldVector<simplegrid_ctype, dimworld> g;
int lk=0;
for (int k=0; k<dimworld; k++)
{
@@ -190,10 +190,10 @@ namespace Dune {
}
//! maps a global coordinate within the element to a local coordinate in its reference element
- Vec<dim,simplegrid_ctype> local (const Vec<dimworld,simplegrid_ctype>& global)
+ FieldVector<simplegrid_ctype, dim> local (const FieldVector<simplegrid_ctype, dimworld>& global)
{
- Vec<dim,simplegrid_ctype> l; // result
- Vec<dimworld,simplegrid_ctype> rhs;
+ FieldVector<simplegrid_ctype, dim> l; // result
+ FieldVector<simplegrid_ctype, dimworld> rhs;
int k;
for (k=0; k<dimworld; k++)
if (k!=c) rhs(k) = global(k)-e.s[k];
@@ -209,7 +209,7 @@ namespace Dune {
return l;
}
- bool checkInside (const Vec<dim,simplegrid_ctype>& local)
+ bool checkInside (const FieldVector<simplegrid_ctype, dim>& local)
{
// check wether they are in the reference element
for(int i=0; i<dim; i++)
@@ -223,7 +223,7 @@ namespace Dune {
/*! Integration over a general element is done by integrating over the reference element
*/
- simplegrid_ctype integration_element (const Vec<dim,simplegrid_ctype>& local)
+ simplegrid_ctype integration_element (const FieldVector<simplegrid_ctype, dim>& local)
{
return e.li->facevol[c];
}
@@ -242,7 +242,7 @@ namespace Dune {
//! print function
void print (std::ostream& ss)
{
- Vec<dim,simplegrid_ctype> local = 0.5;
+ FieldVector<simplegrid_ctype, dim> local = 0.5;
ss << "SimpleElement<"<<dim<<","<<dimworld<< "> ";
ss << "position " << s << " mesh size " << h << " directions " << dir << " vol " << integration_element(local);
}
@@ -251,7 +251,7 @@ namespace Dune {
SimpleElement<dim+1,dimworld>& e; // my element
int c; // direction
int d; // face number in direction d=0,1
- Vec<dimworld,simplegrid_ctype> corn; //!< coordinate vector of corner
+ FieldVector<simplegrid_ctype, dimworld> corn; //!< coordinate vector of corner
};
@@ -286,7 +286,7 @@ namespace Dune {
}
//! access to coordinates of corners. Index is the number of the corner
- Vec<dim,simplegrid_ctype>& operator[] (int i)
+ FieldVector<simplegrid_ctype, dim>& operator[] (int i)
{
for (int k=0; k<dim; k++)
{
@@ -307,30 +307,30 @@ namespace Dune {
}
//! maps a local coordinate within reference element to global coordinate in element
- Vec<dim,simplegrid_ctype> global (const Vec<dim,simplegrid_ctype>& local)
+ FieldVector<simplegrid_ctype, dim> global (const FieldVector<simplegrid_ctype, dim>& local)
{
- Vec<dim,simplegrid_ctype> g;
+ FieldVector<simplegrid_ctype, dim> g;
for (int k=0; k<dim; k++) g(k) = s[k] + local(k)*li->h[k];
return g;
}
//! maps a global coordinate within the element to a local coordinate in its reference element
- Vec<dim,simplegrid_ctype> local (const Vec<dim,simplegrid_ctype>& global)
+ FieldVector<simplegrid_ctype, dim> local (const FieldVector<simplegrid_ctype, dim>& global)
{
- Vec<dim,simplegrid_ctype> l; // result
+ FieldVector<simplegrid_ctype, dim> l; // result
for (int k=0; k<dim; k++) l(k) = (global(k)-s[k])/li->h[k];
return l;
}
/*! Integration over a general element is done by integrating over the reference element
*/
- simplegrid_ctype integration_element (const Vec<dim,simplegrid_ctype>& local)
+ simplegrid_ctype integration_element (const FieldVector<simplegrid_ctype, dim>& local)
{
return li->volume;
}
//! can only be called for dim=dim!
- Mat<dim,dim>& Jacobian_inverse (const Vec<dim,simplegrid_ctype>& local)
+ Mat<dim,dim>& Jacobian_inverse (const FieldVector<simplegrid_ctype, dim>& local)
{
for (int i=0; i<dim; ++i)
{
@@ -380,7 +380,7 @@ namespace Dune {
//! print function
void print (std::ostream& ss)
{
- Vec<dim,simplegrid_ctype> local = 0.5;
+ FieldVector<simplegrid_ctype, dim> local = 0.5;
ss << "SimpleElement<"<<dim<<","<<dim<< "> ";
ss << "position ";
for (int i=0; i<dim; i++) ss << s[i] << " ";
@@ -457,7 +457,7 @@ namespace Dune {
levelinfo<dim>* li; // access to information common to all elements on a level
Mat<dim,dim,sgrid_ctype> Jinv; //!< storage for inverse of jacobian
- Vec<dim,simplegrid_ctype> c; // needed for returning references
+ FieldVector<simplegrid_ctype, dim> c; // needed for returning references
};
@@ -484,7 +484,7 @@ namespace Dune {
}
//! access to coordinates of corners. Index is the number of the corner
- Vec<dimworld,simplegrid_ctype>& operator[] (int i)
+ FieldVector<simplegrid_ctype, dimworld>& operator[] (int i)
{
return s;
}
@@ -556,7 +556,7 @@ namespace Dune {
private:
int id; //!< my id
int coord[dimworld]; //!< my integer position
- Vec<dimworld,simplegrid_ctype> s; //!< pos
+ FieldVector<simplegrid_ctype, dimworld> s; //!< pos
levelinfo<dimworld>& li; //!< common information for all vertices
};
@@ -592,7 +592,7 @@ namespace Dune {
{};
//! return outer barycenter of ghost cell
- Vec<dimworld,simplegrid_ctype> & outerPoint () {};
+ FieldVector<simplegrid_ctype, dimworld> & outerPoint () {};
private:
};
@@ -684,13 +684,13 @@ namespace Dune {
}
//! return unit outer normal, this should be dependent on local coordinates for higher order boundary
- Vec<dimworld,sgrid_ctype>& unit_outer_normal (Vec<dim-1,sgrid_ctype>& local)
+ FieldVector<simplegrid_ctype, dimworld>& unit_outer_normal (FieldVector<simplegrid_ctype, dim-1>& local)
{
return normal;
}
//! return unit outer normal, if you know it is constant use this function instead
- Vec<dimworld,sgrid_ctype>& unit_outer_normal ()
+ FieldVector<simplegrid_ctype, dimworld>& unit_outer_normal ()
{
return normal;
}
@@ -790,7 +790,7 @@ namespace Dune {
SimpleElement<dim-1,dim> is_self_local; //!< intersection in own local coordinates
SimpleElement<dim-1,dim> is_nb_local; //!< intersection in neighbors local coordinates
SimpleElement<dim-1,dimworld> is_global; //!< intersection in global coordinates
- Vec<dimworld,sgrid_ctype> normal;
+ FieldVector<simplegrid_ctype, dimworld> normal;
};
/*!
@@ -1050,7 +1050,7 @@ namespace Dune {
}
//! local coordinates within father
- Vec<dim,simplegrid_ctype>& local () {return loc;} // not implemented yet
+ FieldVector<simplegrid_ctype, dim>& local () {return loc;} // not implemented yet
// members specific to SimpleEntity
SimpleEntity (SimpleElement<0,dimworld>& _geo) : geo(_geo)
@@ -1058,7 +1058,7 @@ namespace Dune {
private:
SimpleElement<0,dimworld>& geo; // reference to geometry which is updated outside of here
- Vec<dim,simplegrid_ctype> loc; // needed for returning a reference
+ FieldVector<simplegrid_ctype, dim> loc; // needed for returning a reference
};
//************************************************************************
diff --git a/grid/uggrid/ugfunctions.hh b/grid/uggrid/ugfunctions.hh
index 3849e76c8..3034720b6 100644
--- a/grid/uggrid/ugfunctions.hh
+++ b/grid/uggrid/ugfunctions.hh
@@ -79,8 +79,8 @@ namespace Dune {
}
static void Local_To_Global(int n, DOUBLE** y,
- const Dune::Vec<dim, double>& local,
- Dune::Vec<dim, double>& global) {
+ const FieldVector<double, dim>& local,
+ FieldVector<double, dim>& global) {
LOCAL_TO_GLOBAL(n,y,local,global);
}
@@ -91,8 +91,7 @@ namespace Dune {
*
*/
static void Transformation(int n, double** x,
- const Vec<dim, double>& local, Mat<dim,dim,double>& mat) {
- //TRANSFORMATION(n, x, local, mat);
+ const FieldVector<double, dim>& local, Mat<dim,dim,double>& mat) {
typedef DOUBLE DOUBLE_VECTOR[dim];
double det;
INVERSE_TRANSFORMATION(n, x, local, mat, det);
@@ -145,6 +144,18 @@ namespace Dune {
return UG2d::GetMultigrid(name);
#endif
}
+
+ static void SetSubdomain(typename TargetType<0,dim>::T* theElement, int id) {
+#ifdef _2
+ using UG2d::control_entries;
+ using UG2d::SUBDOMAIN_CE;
+#else
+ using UG3d::control_entries;
+ using UG3d::SUBDOMAIN_CE;
+#endif
+ SETSUBDOMAIN(theElement, id);
+ }
+
};
diff --git a/grid/uggrid/uggridelement.cc b/grid/uggrid/uggridelement.cc
index 1c2181d78..4fec7779f 100644
--- a/grid/uggrid/uggridelement.cc
+++ b/grid/uggrid/uggridelement.cc
@@ -117,7 +117,7 @@ inline int UGGridElement<dim,dimworld>::corners()
///////////////////////////////////////////////////////////////////////
template<int dim, int dimworld>
-inline const Vec<dimworld,UGCtype>& UGGridElement<dim,dimworld>::
+inline const FieldVector<UGCtype, dimworld>& UGGridElement<dim,dimworld>::
operator [](int i)
{
std::cerr << "UGGridElement<" << dim << "," << dimworld << ">::operator[]:\n"
@@ -127,7 +127,7 @@ operator [](int i)
#ifdef _3
template<>
-inline const Vec<3,UGCtype>& UGGridElement<0,3>::
+inline const FieldVector<UGCtype, 3>& UGGridElement<0,3>::
operator [](int i)
{
const UG3d::VERTEX* vertex = target_->myvertex;
@@ -140,7 +140,7 @@ operator [](int i)
}
template<>
-inline const Vec<3,UGCtype>& UGGridElement<3,3>::
+inline const FieldVector<UGCtype, 3>& UGGridElement<3,3>::
operator [](int i)
{
assert(0<=i && i<corners());
@@ -160,7 +160,7 @@ operator [](int i)
#ifdef _2
template<>
-inline const Vec<2,UGCtype>& UGGridElement<0,2>::
+inline const FieldVector<UGCtype, 2>& UGGridElement<0,2>::
operator [](int i)
{
const UG2d::VERTEX* vertex = target_->myvertex;
@@ -172,7 +172,7 @@ operator [](int i)
}
template<>
-inline const Vec<2,UGCtype>& UGGridElement<2,2>::
+inline const FieldVector<UGCtype, 2>& UGGridElement<2,2>::
operator [](int i)
{
assert(0<=i && i<corners());
@@ -253,10 +253,10 @@ refelem()
template< int dim, int dimworld>
-inline Vec<dimworld,UGCtype> UGGridElement<dim,dimworld>::
-global(const Vec<dim>& local)
+inline FieldVector<UGCtype, dimworld> UGGridElement<dim,dimworld>::
+global(const FieldVector<UGCtype, dim>& local)
{
- Vec<dimworld, UGCtype> globalCoord;
+ FieldVector<UGCtype, dimworld> globalCoord;
// dimworld*dimworld is an upper bound for the number of vertices
UGCtype* cornerCoords[dimworld*dimworld];
@@ -271,10 +271,10 @@ global(const Vec<dim>& local)
// Maps a global coordinate within the element to a
// local coordinate in its reference element
template< int dim, int dimworld>
-inline Vec<dim, UGCtype> UGGridElement<dim,dimworld>::
-local (const Vec<dimworld, UGCtype>& global)
+inline FieldVector<UGCtype, dim> UGGridElement<dim,dimworld>::
+local (const FieldVector<UGCtype, dimworld>& global)
{
- Vec<dim, UGCtype> result;
+ FieldVector<UGCtype, dim> result;
UGCtype localCoords[dim];
// Copy input ADT into C-style array
@@ -304,7 +304,7 @@ local (const Vec<dimworld, UGCtype>& global)
template< int dim, int dimworld>
inline UGCtype UGGridElement<dim,dimworld>::
-integration_element (const Vec<dim,UGCtype>& local)
+integration_element (const FieldVector<UGCtype, dim>& local)
{
//std::cout << "integration element: " << ABS(Jacobian_inverse(local).determinant()) << std::endl;
return ABS(1/Jacobian_inverse(local).determinant());
@@ -312,7 +312,7 @@ integration_element (const Vec<dim,UGCtype>& local)
template< int dim, int dimworld>
inline const Mat<dim,dim>& UGGridElement<dim,dimworld>::
-Jacobian_inverse (const Vec<dim,UGCtype>& local)
+Jacobian_inverse (const FieldVector<UGCtype, dim>& local)
{
// dimworld*dimworld is an upper bound for the number of vertices
UGCtype* cornerCoords[dimworld*dimworld];
diff --git a/grid/uggrid/uggridelement.hh b/grid/uggrid/uggridelement.hh
index 12fbd7421..aa6cb4c20 100644
--- a/grid/uggrid/uggridelement.hh
+++ b/grid/uggrid/uggridelement.hh
@@ -29,9 +29,6 @@ namespace Dune {
//friend class UGGridBoundaryEntity<dim,dimworld>;
public:
- //! know dimension of world
- enum { dimbary=dim+1};
-
//! for makeRefElement == true a Element with the coordinates of the
//! reference element is made
UGGridElement(bool makeRefElement=false);
@@ -44,7 +41,7 @@ namespace Dune {
int corners ();
//! access to coordinates of corners. Index is the number of the corner
- const Vec<dimworld, UGCtype>& operator[] (int i);
+ const FieldVector<UGCtype, dimworld>& operator[] (int i);
/*! return reference element corresponding to this element. If this is
a reference element then self is returned.
@@ -53,14 +50,14 @@ namespace Dune {
//! maps a local coordinate within reference element to
//! global coordinate in element
- Vec<dimworld,UGCtype> global (const Vec<dim, UGCtype>& local);
+ FieldVector<UGCtype, dimworld> global (const FieldVector<UGCtype, dim>& local);
//! Maps a global coordinate within the element to a
//! local coordinate in its reference element
- Vec<dim, UGCtype> local (const Vec<dimworld, UGCtype>& global);
+ FieldVector<UGCtype, dim> local (const FieldVector<UGCtype, dimworld>& global);
//! Returns true if the point is in the current element
- bool checkInside(const Vec<dimworld, UGCtype> &global);
+ bool checkInside(const FieldVector<UGCtype, dimworld> &global);
/*!
Copy from sgrid.hh:
@@ -86,10 +83,10 @@ namespace Dune {
*/
// A(l)
- UGCtype integration_element (const Vec<dim, UGCtype>& local);
+ UGCtype integration_element (const FieldVector<UGCtype, dim>& local);
//! can only be called for dim=dimworld!
- const Mat<dim,dim>& Jacobian_inverse (const Vec<dim, UGCtype>& local);
+ const Mat<dim,dim>& Jacobian_inverse (const FieldVector<UGCtype, dim>& local);
//***********************************************************************
// Methods that not belong to the Interface, but have to be public
@@ -103,14 +100,6 @@ namespace Dune {
//! built the reference element
void makeRefElemCoords();
- //! maps a global coordinate within the elements local barycentric
- //! koordinates
- Vec<dim+1, UGCtype> localBary (const Vec<dimworld, UGCtype>& global);
-
- // template method for map the vertices of EL_INFO to the actual
- // coords with face_,edge_ and vertex_ , needes for operator []
- int mapVertices (int i) const;
-
//! the vertex coordinates
Mat<dimworld,dim+1, UGCtype> coord_;
@@ -118,18 +107,179 @@ namespace Dune {
Mat<dimworld,dimworld> jac_inverse_;
//! storage for global coords
- Vec<dimworld, UGCtype> globalCoord_;
+ FieldVector<UGCtype, dimworld> globalCoord_;
//! storage for local coords
- Vec<dim, UGCtype> localCoord_;
-
- //! storage for barycentric coords
- Vec<dimbary, UGCtype> localBary_;
+ FieldVector<UGCtype, dim> localCoord_;
typename TargetType<dimworld-dim,dimworld>::T* target_;
};
+
+ /****************************************************************/
+ /* Specialization for faces in 3d */
+ /****************************************************************/
+
+
+ template<>
+ class UGGridElement<2, 3> :
+ public ElementDefault <2, 3, UGCtype,UGGridElement>
+ {
+ //friend class UGGridBoundaryEntity<dim,dimworld>;
+ public:
+
+ //! for makeRefElement == true a Element with the coordinates of the
+ //! reference element is made
+ UGGridElement(bool makeRefElement=false){
+ std::cout << "UGGridElement<2,3> created!" << std::endl;
+ }
+
+ //! return the element type identifier (triangle or quadrilateral)
+ ElementType type ();
+
+ //! return the number of corners of this element. Corners are numbered 0...n-1
+ int corners ();
+
+ //! access to coordinates of corners. Index is the number of the corner
+ const FieldVector<UGCtype, 3>& operator[] (int i);
+
+ /*! return reference element corresponding to this element. If this is
+ a reference element then self is returned.
+ */
+ UGGridElement<2,2>& refelem ();
+
+ //! maps a local coordinate within reference element to
+ //! global coordinate in element
+ FieldVector<UGCtype, 3> global (const FieldVector<UGCtype, 2>& local);
+
+ //! Maps a global coordinate within the element to a
+ //! local coordinate in its reference element
+ FieldVector<UGCtype, 2> local (const FieldVector<UGCtype, 3>& global);
+
+ //! Returns true if the point is in the current element
+ bool checkInside(const FieldVector<UGCtype, 3> &global);
+
+ // A(l)
+ UGCtype integration_element (const FieldVector<UGCtype, 2>& local);
+
+ //! can only be called for dim=dimworld!
+ const Mat<2,2>& Jacobian_inverse (const FieldVector<UGCtype, 2>& local);
+
+ //***********************************************************************
+ // Methods that not belong to the Interface, but have to be public
+ //***********************************************************************
+
+ //void setToTarget(typename TargetType<dimworld-dim,dimworld>::T* target) {target_ = target;}
+ void setToTarget(TargetType<2,3>::T* target) {
+ DUNE_THROW(GridError, "UGGridElement<2,3>::setToTarget called!");
+ }
+
+ private:
+
+
+ //! built the reference element
+ void makeRefElemCoords();
+
+ //! the vertex coordinates
+ Mat<3,3, UGCtype> coord_;
+
+ //! The jacobian inverse
+ Mat<3,3> jac_inverse_;
+
+ //! storage for global coords
+ FieldVector<UGCtype, 3> globalCoord_;
+
+ //! storage for local coords
+ FieldVector<UGCtype, 2> localCoord_;
+
+ //typename TargetType<dimworld-dim,dimworld>::T* target_;
+
+ };
+
+
+ /****************************************************************/
+ /* Specialization for faces in 2d */
+ /****************************************************************/
+
+
+ template<>
+ class UGGridElement <1, 2> :
+ public ElementDefault <1, 2, UGCtype,UGGridElement>
+ {
+ //friend class UGGridBoundaryEntity<dim,dimworld>;
+ public:
+
+ //! for makeRefElement == true a Element with the coordinates of the
+ //! reference element is made
+ UGGridElement(bool makeRefElement=false) {
+ std::cout << "UGGridElement<1,2> created!" << std::endl;
+ }
+
+ //! return the element type identifier
+ //! line , triangle or tetrahedron, depends on dim
+ ElementType type () const {return line;}
+
+ //! return the number of corners of this element. Corners are numbered 0...n-1
+ int corners () const {return 2;}
+
+ //! access to coordinates of corners. Index is the number of the corner
+ const FieldVector<UGCtype, 2>& operator[] (int i) const;
+
+ /*! return reference element corresponding to this element. If this is
+ a reference element then self is returned.
+ */
+ UGGridElement<1,1>& refelem ();
+
+ //! maps a local coordinate within reference element to
+ //! global coordinate in element
+ FieldVector<UGCtype, 2> global (const FieldVector<UGCtype, 1>& local);
+
+ //! Maps a global coordinate within the element to a
+ //! local coordinate in its reference element
+ FieldVector<UGCtype, 1> local (const FieldVector<UGCtype, 2>& global);
+
+ //! Returns true if the point is in the current element
+ bool checkInside(const FieldVector<UGCtype, 2> &global);
+
+ // A(l)
+ UGCtype integration_element (const FieldVector<UGCtype, 1>& local);
+
+ //! can only be called for dim=dimworld!
+ const Mat<1,1>& Jacobian_inverse (const FieldVector<UGCtype, 1>& local);
+
+ //***********************************************************************
+ // Methods that not belong to the Interface, but have to be public
+ //***********************************************************************
+
+ //void setToTarget(typename TargetType<dimworld-dim,dimworld>::T* target) {target_ = target;}
+ void setToTarget(TargetType<1,2>::T* target) {
+ DUNE_THROW(GridError, "UGGridElement<1,2>::setToTarget called!");
+ }
+
+ private:
+
+
+ //! built the reference element
+ void makeRefElemCoords();
+
+ //! the vertex coordinates
+ Mat<2,2, UGCtype> coord_;
+
+ //! The jacobian inverse
+ Mat<2,2> jac_inverse_;
+
+ //! storage for global coords
+ FieldVector<UGCtype, 2> globalCoord_;
+
+ //! storage for local coords
+ FieldVector<UGCtype, 1> localCoord_;
+
+ //typename TargetType<dimworld-dim,dimworld>::T* target_;
+
+ };
+
+
} // namespace Dune
#endif
diff --git a/grid/uggrid/uggridentity.cc b/grid/uggrid/uggridentity.cc
index fdbb3061f..aba454083 100644
--- a/grid/uggrid/uggridentity.cc
+++ b/grid/uggrid/uggridentity.cc
@@ -96,15 +96,11 @@ geometry()
}
#ifdef _3
-//template<int codim, int dim, int dimworld>
template<>
-inline UGGridLevelIterator<0,3,3>
+inline UGGridLevelIterator<0,3,3, All_Partition>
UGGridEntity < 0, 3 ,3>::father()
{
- // This currently only works for elements
- //assert(codim==0);
-
- UGGridLevelIterator<0,3,3> it(level()-1);
+ UGGridLevelIterator<0,3,3, All_Partition> it(level()-1);
#define TAG(p) ReadCW(p, UG3d::TAG_CE)
#define EFATHER(p) ((UG3d::ELEMENT *) (p)->ge.refs[UG3d::father_offset[TAG(p)]])
UG3d::ELEMENT* fatherTarget = EFATHER(target_);
@@ -153,15 +149,9 @@ geometry()
//*****************************************************************8
// count
-/** \todo So far only works in 3d */
template <int codim, int dim, int dimworld> template <int cc>
inline int UGGridEntity<codim,dim,dimworld>::count ()
{
- // #define TAG(p) ReadCW(p, UG3d::TAG_CE)
- // #define SIDES_OF_ELEM(p) (UG3d::element_descriptors[TAG(p)]->sides_of_elem)
- // #define EDGES_OF_ELEM(p) (UG3d::element_descriptors[TAG(p)]->edges_of_elem)
- // #define CORNERS_OF_ELEM(p)(UG3d::element_descriptors[TAG(p)]->corners_of_elem)
-
if (dim==3) {
switch (cc) {
@@ -189,34 +179,14 @@ inline int UGGridEntity<codim,dim,dimworld>::count ()
}
return -1;
- // #undef SIDES_OF_ELEM
- // #undef EDGES_OF_ELEM
- // #undef CORNERS_OF_ELEM
- // #undef TAG
}
template <int codim, int dim, int dimworld>
template <int cc>
inline int UGGridEntity<codim, dim, dimworld>::subIndex(int i)
{
-#if 0
-
-#define TAG(p) ReadCW(p, UG3d::TAG_CE)
-#define CORNER(p,i) ((UG3d::node *) (p)->ge.refs[UG3d::n_offset[TAG(p)]+(i)])
- UG3d::node* node = CORNER(target_,i);
-#undef CORNER
-#undef TAG
- UG3d::vertex* vertex = node->myvertex;
- return vertex->iv.id;
-
-#else
-
typename TargetType<dim,dim>::T* node = CORNER(target_,i);
- //UG3d::vertex* vertex = node->myvertex;
return node->myvertex->iv.id;
-
-#endif
-
}
@@ -249,11 +219,19 @@ template <int dim, int dimworld>
template <int cc>
inline int UGGridEntity<0, dim, dimworld>::subIndex(int i)
{
-#ifdef _3
- UG3d::node* node = (UG3d::node*)UG_NS<dimworld>::Corner(target_,i);
-#else
- UG2d::node* node = (UG2d::node*)UG_NS<dimworld>::Corner(target_,i);
-#endif
+ //assert(i>=0 && i<count());
+
+ if (cc!=dim)
+ DUNE_THROW(GridError, "UGGrid::subIndex isn't implemented for cc != dim");
+
+ if (geometry().type()==hexahedron) {
+ // Dune numbers the vertices of a hexahedron differently than UG.
+ // The following two lines do the transformation
+ const int renumbering[8] = {0, 1, 3, 2, 4, 5, 7, 6};
+ i = renumbering[i];
+ }
+
+ typename TargetType<dim,dim>::T* node = UG_NS<dimworld>::Corner(target_,i);
return node->myvertex->iv.id;
}
@@ -315,9 +293,17 @@ UGGridEntity < 0, dim ,dimworld >::hbegin(int maxlevel)
// The 30 is the macro MAX_SONS from ug/gm/gm.h
UGElementType* sonList[30];
+#ifdef _2
UG2d::GetSons(target_,sonList);
+#else
+ UG3d::GetSons(target_,sonList);
+#endif
+#ifdef _2
#define NSONS(p) UG2d::ReadCW(p, UG2d::NSONS_CE)
+#else
+#define NSONS(p) UG3d::ReadCW(p, UG3d::NSONS_CE)
+#endif
// Load sons of entity into the iterator
for (unsigned int i=0; i<NSONS(target_); i++) {
typename UGGridHierarchicIterator<dim,dimworld>::StackEntry se;
diff --git a/grid/uggrid/uggridentity.hh b/grid/uggrid/uggridentity.hh
index 761313a8e..28e58d199 100644
--- a/grid/uggrid/uggridentity.hh
+++ b/grid/uggrid/uggridentity.hh
@@ -79,7 +79,7 @@ namespace Dune {
UGGridLevelIterator<0,dim,dimworld,All_Partition> father ();
//! local coordinates within father
- Vec<dim, UGCtype>& local ();
+ FieldVector<UGCtype, dim>& local ();
@@ -97,7 +97,7 @@ namespace Dune {
UGGridElement<dim-codim,dimworld> geo_;
bool builtgeometry_; //!< true if geometry has been constructed
- Vec<dim,UGCtype> localFatherCoords_;
+ FieldVector<UGCtype, dim> localFatherCoords_;
public:
//! element number
int elNum_;
diff --git a/grid/uggrid/uggridhieriterator.cc b/grid/uggrid/uggridhieriterator.cc
index cd2aa0ea1..ddf169841 100644
--- a/grid/uggrid/uggridhieriterator.cc
+++ b/grid/uggrid/uggridhieriterator.cc
@@ -32,9 +32,17 @@ UGGridHierarchicIterator< dim,dimworld >::operator ++()
// The 30 is the macro MAX_SONS from ug/gm/gm.h
UGElementType* sonList[30];
+#ifdef _2
UG2d::GetSons(old_target.element,sonList);
+#else
+ UG3d::GetSons(old_target.element,sonList);
+#endif
+#ifdef _2
#define NSONS(p) UG2d::ReadCW(p, UG2d::NSONS_CE)
+#else
+#define NSONS(p) UG3d::ReadCW(p, UG3d::NSONS_CE)
+#endif
// Load sons of old target onto the iterator stack
for (unsigned int i=0; i<NSONS(old_target.element); i++) {
StackEntry se;
diff --git a/grid/uggrid/uggridleveliterator.cc b/grid/uggrid/uggridleveliterator.cc
index d366c0ca3..2cf2a7969 100644
--- a/grid/uggrid/uggridleveliterator.cc
+++ b/grid/uggrid/uggridleveliterator.cc
@@ -31,8 +31,8 @@ UGGridLevelIterator(UGGrid<dim,dimworld> &grid, int travLevel) :
#ifdef _3
template<>
-inline UGGridLevelIterator < 3,3,3 >&
-UGGridLevelIterator < 3,3,3 >::operator++()
+inline UGGridLevelIterator < 3,3,3, All_Partition >&
+UGGridLevelIterator < 3,3,3,All_Partition >::operator++()
{
target_ = target_->succ;
@@ -44,8 +44,8 @@ UGGridLevelIterator < 3,3,3 >::operator++()
}
template<>
-inline UGGridLevelIterator < 0,3,3 >&
-UGGridLevelIterator < 0,3,3 >::operator++()
+inline UGGridLevelIterator < 0,3,3, All_Partition >&
+UGGridLevelIterator < 0,3,3, All_Partition >::operator++()
{
setToTarget(target_->ge.succ);
virtualEntity_.elNum_++;
diff --git a/grid/uggrid/ugintersectionit.cc b/grid/uggrid/ugintersectionit.cc
index 5ad83acc4..5c85d1455 100644
--- a/grid/uggrid/ugintersectionit.cc
+++ b/grid/uggrid/ugintersectionit.cc
@@ -112,7 +112,7 @@ namespace Dune {
#endif
template<>
- inline Vec<3,UGCtype>&
+ inline FieldVector<UGCtype, 3>&
UGGridIntersectionIterator < 3,3 >::unit_outer_normal ()
{
#ifdef _3
@@ -124,7 +124,7 @@ namespace Dune {
UG3d::VERTEX* c = UG_NS<3>::Corner(center_,CORNER_OF_SIDE(center_, neighborCount_, 2))->myvertex;
#undef CORNER_OF_SIDE
- Vec<3,UGCtype> aPos, bPos, cPos;
+ FieldVector<UGCtype, 3> aPos, bPos, cPos;
#define CVECT(p) ((p)->iv.x)
#define V3_COPY(A,C) {(C)[0] = (A)[0]; (C)[1] = (A)[1]; (C)[2] = (A)[2];\
@@ -135,8 +135,8 @@ namespace Dune {
#undef CVECT
#undef V3_COPY
- Vec<3> ba = bPos - aPos;
- Vec<3> ca = cPos - aPos;
+ FieldVector<UGCtype, 3> ba = bPos - aPos;
+ FieldVector<UGCtype, 3> ca = cPos - aPos;
// std::cout << "aPos: " << aPos << std::endl;
// std::cout << "bPos: " << bPos << std::endl;
@@ -157,7 +157,7 @@ namespace Dune {
}
template<>
- inline Vec<2,UGCtype>&
+ inline FieldVector<UGCtype, 2>&
UGGridIntersectionIterator < 2,2 >::unit_outer_normal ()
{
#ifdef _2
@@ -172,15 +172,15 @@ namespace Dune {
outerNormal_[1] = aPos[0] - bPos[0];
// normalize
- outerNormal_ *= (1/outerNormal_.norm2());
+ outerNormal_ *= (1/outerNormal_.two_norm());
#endif
return outerNormal_;
}
template<int dim, int dimworld>
- inline Vec<dimworld,UGCtype>&
+ inline FieldVector<UGCtype, dimworld>&
UGGridIntersectionIterator < dim,dimworld >::
- unit_outer_normal (Vec<dim-1,UGCtype>& local)
+ unit_outer_normal (const FieldVector<UGCtype, dim-1>& local)
{
return unit_outer_normal();
}
diff --git a/grid/uggrid/ugintersectionit.hh b/grid/uggrid/ugintersectionit.hh
index 73c1a74a1..ac48a093d 100644
--- a/grid/uggrid/ugintersectionit.hh
+++ b/grid/uggrid/ugintersectionit.hh
@@ -57,10 +57,10 @@ namespace Dune {
//! return unit outer normal, this should be dependent on local
//! coordinates for higher order boundary
- Vec<dimworld,UGCtype>& unit_outer_normal (Vec<dim-1,UGCtype>& local);
+ FieldVector<UGCtype, dimworld>& unit_outer_normal (const FieldVector<UGCtype, dim-1>& local);
//! return unit outer normal, if you know it is constant use this function instead
- Vec<dimworld,UGCtype>& unit_outer_normal ();
+ FieldVector<UGCtype, dimworld>& unit_outer_normal ();
//! intersection of codimension 1 of this neighbor with element where
//! iteration started.
@@ -88,10 +88,10 @@ namespace Dune {
//! return outer normal, this should be dependent on local
//! coordinates for higher order boundary
- Vec<dimworld,UGCtype>& outer_normal (Vec<dim-1,UGCtype>& local);
+ FieldVector<UGCtype, dimworld>& outer_normal (const FieldVector<UGCtype, dim-1>& local);
//! return unit outer normal, if you know it is constant use this function instead
- Vec<dimworld,UGCtype>& outer_normal ();
+ FieldVector<UGCtype, dimworld>& outer_normal ();
private:
//**********************************************************
@@ -114,7 +114,7 @@ namespace Dune {
UGGridEntity<0,dim,dimworld> virtualEntity_;
//! vector storing the outer normal
- Vec<dimworld,UGCtype> outerNormal_;
+ FieldVector<UGCtype, dimworld> outerNormal_;
//! pointer to element holding the self_local and self_global information.
//! This element is created on demand.
diff --git a/grid/yaspgrid.hh b/grid/yaspgrid.hh
index c519fbacc..e60aada9b 100644
--- a/grid/yaspgrid.hh
+++ b/grid/yaspgrid.hh
@@ -60,8 +60,8 @@ namespace Dune {
/** Singleton holding reference elements */
template<int dim>
struct YaspReferenceElement {
- static Vec<dim,yaspgrid_ctype> midpoint; // data neded for the refelem below
- static Vec<dim,yaspgrid_ctype> extension; // data needed for the refelem below
+ static FieldVector<yaspgrid_ctype, dim> midpoint; // data neded for the refelem below
+ static FieldVector<yaspgrid_ctype, dim> extension; // data needed for the refelem below
static YaspElement<dim,dim> refelem;
};
@@ -70,16 +70,16 @@ namespace Dune {
YaspElement<dim,dim> YaspReferenceElement<dim>::refelem(YaspReferenceElement<dim>::midpoint,
YaspReferenceElement<dim>::extension);
template<int dim>
- Vec<dim,yaspgrid_ctype> YaspReferenceElement<dim>::midpoint(0.5);
+ FieldVector<yaspgrid_ctype, dim> YaspReferenceElement<dim>::midpoint(0.5);
template<int dim>
- Vec<dim,yaspgrid_ctype> YaspReferenceElement<dim>::extension(1.0);
+ FieldVector<yaspgrid_ctype, dim> YaspReferenceElement<dim>::extension(1.0);
template<int dim>
class YaspFatherRelativeLocalElement {
public:
- static Vec<dim,yaspgrid_ctype> midpoint; // data neded for the refelem below
- static Vec<dim,yaspgrid_ctype> extension; // data needed for the refelem below
+ static FieldVector<yaspgrid_ctype, dim> midpoint; // data neded for the refelem below
+ static FieldVector<yaspgrid_ctype, dim> extension; // data needed for the refelem below
static YaspElement<dim,dim> element;
static YaspElement<dim,dim>& getson (int i)
{
@@ -97,10 +97,10 @@ namespace Dune {
YaspElement<dim,dim> YaspFatherRelativeLocalElement<dim>::element(YaspFatherRelativeLocalElement<dim>::midpoint,
YaspFatherRelativeLocalElement<dim>::extension);
template<int dim>
- Vec<dim,yaspgrid_ctype> YaspFatherRelativeLocalElement<dim>::midpoint(0.25);
+ FieldVector<yaspgrid_ctype, dim> YaspFatherRelativeLocalElement<dim>::midpoint(0.25);
template<int dim>
- Vec<dim,yaspgrid_ctype> YaspFatherRelativeLocalElement<dim>::extension(0.5);
+ FieldVector<yaspgrid_ctype, dim> YaspFatherRelativeLocalElement<dim>::extension(0.5);
//========================================================================
@@ -140,7 +140,7 @@ namespace Dune {
}
//! access to coordinates of corners. Index is the number of the corner
- Vec<dim,yaspgrid_ctype>& operator[] (int i)
+ FieldVector<yaspgrid_ctype, dim>& operator[] (int i)
{
int bit=0;
for (int k=0; k<dimworld; k++) // run over all directions in world
@@ -172,9 +172,9 @@ namespace Dune {
}
//! maps a local coordinate within reference element to global coordinate in element
- Vec<dimworld,yaspgrid_ctype> global (const Vec<dim,yaspgrid_ctype>& local)
+ FieldVector<yaspgrid_ctype,dimworld> global (const FieldVector<yaspgrid_ctype, dim>& local)
{
- Vec<dimworld,yaspgrid_ctype> g;
+ FieldVector<yaspgrid_ctype,dimworld> g;
int bit=0;
for (int k=0; k<dimworld; k++)
if (k==missing)
@@ -188,9 +188,9 @@ namespace Dune {
}
//! maps a global coordinate within the element to a local coordinate in its reference element
- Vec<dim,yaspgrid_ctype> local (const Vec<dimworld,yaspgrid_ctype>& global)
+ FieldVector<yaspgrid_ctype,dim> local (const FieldVector<yaspgrid_ctype,dimworld>& global)
{
- Vec<dim,yaspgrid_ctype> l; // result
+ FieldVector<yaspgrid_ctype,dim> l; // result
int bit=0;
for (int k=0; k<dimworld; k++)
if (k!=missing)
@@ -203,7 +203,7 @@ namespace Dune {
/*! determinant of the jacobian of the mapping
*/
- yaspgrid_ctype integration_element (const Vec<dim,yaspgrid_ctype>& local)
+ yaspgrid_ctype integration_element (const FieldVector<yaspgrid_ctype,dim>& local)
{
yaspgrid_ctype volume=1.0;
for (int k=0; k<dimworld; k++)
@@ -212,7 +212,7 @@ namespace Dune {
}
//! check whether local is inside reference element
- bool checkInside (const Vec<dim,yaspgrid_ctype>& local)
+ bool checkInside (const FieldVector<yaspgrid_ctype,dim>& local)
{
for (int i=0; i<dim; i++)
if (local(i)<-yasptolerance || local(i)>1+yasptolerance) return false;
@@ -220,7 +220,7 @@ namespace Dune {
}
//! constructor from (storage for) midpoint and extension and missing direction number
- YaspElement (Vec<dimworld,yaspgrid_ctype>& p, Vec<dimworld,yaspgrid_ctype>& h, int& m)
+ YaspElement (FieldVector<yaspgrid_ctype, dimworld>& p, FieldVector<yaspgrid_ctype, dimworld>& h, int& m)
: midpoint(p), extension(h), missing(m)
{
if (dimworld!=dim+1)
@@ -246,13 +246,13 @@ namespace Dune {
// References are used because this information
// is known outside the element in many cases.
// Note dimworld==dim+1
- Vec<dimworld,yaspgrid_ctype>& midpoint; // the midpoint
- Vec<dimworld,yaspgrid_ctype>& extension; // the extension
+ FieldVector<yaspgrid_ctype, dimworld>& midpoint; // the midpoint
+ FieldVector<yaspgrid_ctype, dimworld>& extension; // the extension
int& missing; // the missing, i.e. constant direction
// In addition we need memory in order to return references.
// Possibly we should change this in the interface ...
- Vec<dim,yaspgrid_ctype> c; // a point
+ FieldVector<yaspgrid_ctype, dim> c; // a point
};
@@ -284,7 +284,7 @@ namespace Dune {
}
//! access to coordinates of corners. Index is the number of the corner
- Vec<dim,yaspgrid_ctype>& operator[] (int i)
+ FieldVector<yaspgrid_ctype, dim>& operator[] (int i)
{
for (int k=0; k<dim; k++)
if (i&(1<<k))
@@ -305,18 +305,18 @@ namespace Dune {
}
//! maps a local coordinate within reference element to global coordinate in element
- Vec<dim,yaspgrid_ctype> global (const Vec<dim,yaspgrid_ctype>& local)
+ FieldVector<yaspgrid_ctype, dim> global (const FieldVector<yaspgrid_ctype, dim>& local)
{
- Vec<dim,yaspgrid_ctype> g;
+ FieldVector<yaspgrid_ctype,dim> g;
for (int k=0; k<dim; k++)
g(k) = midpoint(k) + (local(k)-0.5)*extension(k);
return g;
}
//! maps a global coordinate within the element to a local coordinate in its reference element
- Vec<dim,yaspgrid_ctype> local (const Vec<dim,yaspgrid_ctype>& global)
+ FieldVector<yaspgrid_ctype, dim> local (const FieldVector<yaspgrid_ctype,dim>& global)
{
- Vec<dim,yaspgrid_ctype> l; // result
+ FieldVector<yaspgrid_ctype, dim> l; // result
for (int k=0; k<dim; k++)
l(k) = (global(k)-midpoint(k))/extension(k) + 0.5;
return l;
@@ -324,7 +324,7 @@ namespace Dune {
/*! determinant of the jacobian of the mapping
*/
- yaspgrid_ctype integration_element (const Vec<dim,yaspgrid_ctype>& local)
+ yaspgrid_ctype integration_element (const FieldVector<yaspgrid_ctype, dim>& local)
{
yaspgrid_ctype volume=1.0;
for (int k=0; k<dim; k++) volume *= extension(k);
@@ -332,7 +332,7 @@ namespace Dune {
}
//! can only be called for dim=dim!
- Mat<dim,dim>& Jacobian_inverse (const Vec<dim,yaspgrid_ctype>& local)
+ Mat<dim,dim>& Jacobian_inverse (const FieldVector<yaspgrid_ctype, dim>& local)
{
for (int i=0; i<dim; ++i)
{
@@ -343,7 +343,7 @@ namespace Dune {
}
//! check whether local is inside reference element
- bool checkInside (const Vec<dim,yaspgrid_ctype>& local)
+ bool checkInside (const FieldVector<yaspgrid_ctype, dim>& local)
{
for (int i=0; i<dim; i++)
if (local(i)<-yasptolerance || local(i)>1+yasptolerance) return false;
@@ -351,7 +351,7 @@ namespace Dune {
}
//! constructor from (storage for) midpoint and extension
- YaspElement (Vec<dim,yaspgrid_ctype>& p, Vec<dim,yaspgrid_ctype>& h)
+ YaspElement (FieldVector<yaspgrid_ctype, dim>& p, FieldVector<yaspgrid_ctype, dim>& h)
: midpoint(p), extension(h)
{}
@@ -372,13 +372,13 @@ namespace Dune {
// in each direction. References are used because this information
// is known outside the element in many cases.
// Note dim==dimworld
- Vec<dim,yaspgrid_ctype>& midpoint; // the midpoint
- Vec<dim,yaspgrid_ctype>& extension; // the extension
+ FieldVector<yaspgrid_ctype, dim>& midpoint; // the midpoint
+ FieldVector<yaspgrid_ctype, dim>& extension; // the extension
// In addition we need memory in order to return references.
// Possibly we should change this in the interface ...
Mat<dim,dim,yaspgrid_ctype> Jinv; // the jacobian inverse
- Vec<dim,yaspgrid_ctype> c; // a point
+ FieldVector<yaspgrid_ctype, dim> c; // a point
};
@@ -405,13 +405,13 @@ namespace Dune {
}
//! access to coordinates of corners. Index is the number of the corner
- Vec<dimworld,yaspgrid_ctype>& operator[] (int i)
+ FieldVector<yaspgrid_ctype, dimworld>& operator[] (int i)
{
return position;
}
//! constructor
- YaspElement (Vec<dimworld,yaspgrid_ctype>& p) : position(p)
+ YaspElement (FieldVector<yaspgrid_ctype, dimworld>& p) : position(p)
{}
//! print function
@@ -422,7 +422,7 @@ namespace Dune {
}
private:
- Vec<dimworld,yaspgrid_ctype>& position; //!< where the vertex is
+ FieldVector<yaspgrid_ctype, dimworld>& position; //!< where the vertex is
};
// operator<< for all YaspElements
@@ -693,7 +693,7 @@ namespace Dune {
}
//! local coordinates within father
- Vec<dim,ctype>& local ()
+ FieldVector<ctype, dim>& local ()
{
// check if coarse level exists
if (_g.level<=0)
@@ -728,7 +728,7 @@ namespace Dune {
TSI& _it; // position in the grid level
YGLI& _g; // access to grid level
YaspElement<0,dimworld> _element; // the element geometry
- Vec<dim,ctype> loc; // always computed before being returned
+ FieldVector<ctype, dim> loc; // always computed before being returned
};
@@ -872,13 +872,13 @@ namespace Dune {
}
//! return unit outer normal, this should be dependent on local coordinates for higher order boundary
- Vec<dimworld,yaspgrid_ctype>& unit_outer_normal (Vec<dim-1,yaspgrid_ctype>& local)
+ FieldVector<yaspgrid_ctype, dimworld>& unit_outer_normal (FieldVector<yaspgrid_ctype, dim-1>& local)
{
return _normal;
}
//! return unit outer normal, if you know it is constant use this function instead
- Vec<dimworld,yaspgrid_ctype>& unit_outer_normal ()
+ FieldVector<yaspgrid_ctype, dimworld>& unit_outer_normal ()
{
return _normal;
}
@@ -973,14 +973,14 @@ namespace Dune {
TSI _itnb; //!< position of nb in the grid level
YaspEntity<0,dim,dimworld>& _myself; //!< reference to myself
YaspEntity<0,dim,dimworld> _nb; //!< virtual neighbor entity, built on the fly
- Vec<dim,yaspgrid_ctype> _pos_self_local; //!< center of face in own local coordinates
- Vec<dim,yaspgrid_ctype> _pos_nb_local; //!< center of face in neighbors local coordinates
- Vec<dim,yaspgrid_ctype> _pos_world; //!< center of face in world coordinates
- Vec<dim,yaspgrid_ctype> _ext_local; //!< extension of face in local coordinates
+ FieldVector<yaspgrid_ctype, dim> _pos_self_local; //!< center of face in own local coordinates
+ FieldVector<yaspgrid_ctype, dim> _pos_nb_local; //!< center of face in neighbors local coordinates
+ FieldVector<yaspgrid_ctype, dim> _pos_world; //!< center of face in world coordinates
+ FieldVector<yaspgrid_ctype, dim> _ext_local; //!< extension of face in local coordinates
YaspElement<dim-1,dim> _is_self_local; //!< intersection in own local coordinates
YaspElement<dim-1,dim> _is_nb_local; //!< intersection in neighbors local coordinates
YaspElement<dim-1,dimworld> _is_global; //!< intersection in global coordinates
- Vec<dimworld,yaspgrid_ctype> _normal; //!< for returning outer normal
+ FieldVector<yaspgrid_ctype, dimworld> _normal; //!< for returning outer normal
};
@@ -1223,7 +1223,9 @@ namespace Dune {
@param periodic tells if direction is periodic or not
@param size of overlap on coarsest grid (same in all directions)
*/
- YaspGrid (MPI_Comm comm, Dune::Vec<dim,ctype> L, Dune::Vec<dim,int> s, Dune::Vec<dim,bool> periodic, int overlap)
+ YaspGrid (MPI_Comm comm, Dune::FieldVector<ctype, dim> L,
+ Dune::FieldVector<int, dim> s,
+ Dune::FieldVector<bool, dim> periodic, int overlap)
: _mg(comm,L,s,periodic,overlap)
{ }
diff --git a/grid/yaspgrid/grids.hh b/grid/yaspgrid/grids.hh
index 77b54011f..4f020a86c 100644
--- a/grid/yaspgrid/grids.hh
+++ b/grid/yaspgrid/grids.hh
@@ -16,7 +16,7 @@
// local includes
#include "dune/common/array.hh"
-#include "dune/common/fixedvector.hh"
+#include "dune/common/fvector.hh"
/*! \file grids.hh
This is the basis for the yaspgrid implementation of the Dune grid interface.
@@ -71,9 +71,9 @@ namespace Dune {
class YGrid {
public:
//! define types used for arguments
- typedef Dune::Vec<d,int> iTupel;
- typedef Dune::Vec<d,ct> fTupel;
- typedef Dune::Vec<d,bool> bTupel;
+ typedef FieldVector<int, d> iTupel;
+ typedef FieldVector<ct, d> fTupel;
+ typedef FieldVector<bool, d> bTupel;
//! Make an empty YGrid with origin 0
YGrid ()
@@ -985,8 +985,8 @@ namespace Dune {
class Torus {
public:
//! type used to pass tupels in and out
- typedef Dune::Vec<d,int> iTupel;
- typedef Dune::Vec<d,bool> bTupel;
+ typedef FieldVector<int, d> iTupel;
+ typedef FieldVector<bool, d> bTupel;
private:
@@ -1335,7 +1335,7 @@ namespace Dune {
std::cout << "[" << rank() << "]: ERROR: local sends/receives do not match in exchange!" << std::endl;
return;
}
- for (int i=0; i<_localsendrequests.size(); i++)
+ for (unsigned int i=0; i<_localsendrequests.size(); i++)
{
if (_localsendrequests[i].size!=_localrecvrequests[i].size)
{
@@ -1352,7 +1352,7 @@ namespace Dune {
int recvs=0;
// issue sends to foreign processes
- for (int i=0; i<_sendrequests.size(); i++)
+ for (unsigned int i=0; i<_sendrequests.size(); i++)
if (_sendrequests[i].rank!=rank())
{
MPI_Isend(_sendrequests[i].buffer, _sendrequests[i].size, MPI_BYTE,
@@ -1362,7 +1362,7 @@ namespace Dune {
}
// issue receives from foreign processes
- for (int i=0; i<_recvrequests.size(); i++)
+ for (unsigned int i=0; i<_recvrequests.size(); i++)
if (_recvrequests[i].rank!=rank())
{
MPI_Irecv(_recvrequests[i].buffer, _recvrequests[i].size, MPI_BYTE,
@@ -1374,7 +1374,7 @@ namespace Dune {
// poll sends
while (sends>0)
{
- for (int i=0; i<_sendrequests.size(); i++)
+ for (unsigned int i=0; i<_sendrequests.size(); i++)
if (!_sendrequests[i].flag)
{
MPI_Status status;
@@ -1387,7 +1387,7 @@ namespace Dune {
// poll receives
while (recvs>0)
{
- for (int i=0; i<_recvrequests.size(); i++)
+ for (unsigned int i=0; i<_recvrequests.size(); i++)
if (!_recvrequests[i].flag)
{
MPI_Status status;
@@ -1613,9 +1613,9 @@ namespace Dune {
};
//! define types used for arguments
- typedef Dune::Vec<d,int> iTupel;
- typedef Dune::Vec<d,ct> fTupel;
- typedef Dune::Vec<d,bool> bTupel;
+ typedef FieldVector<int, d> iTupel;
+ typedef FieldVector<ct, d> fTupel;
+ typedef FieldVector<bool, d> bTupel;
// communication tag used by multigrid
enum { tag = 17 };
@@ -1645,7 +1645,7 @@ namespace Dune {
if (_torus.rank()==0) std::cout << "MultiYGrid<" << d
<< ">: coarse grid with size " << s
<< " imbalance=" << (imbal-1)*100 << "%" << std::endl;
- // print(std::cout);
+ // print(std::cout);
}
//! do a global mesh refinement; true: keep overlap in absolute size; false: keep overlap in mesh cells
--
GitLab