From 55cc41d88b7636765ac2721a7959f90777eb2502 Mon Sep 17 00:00:00 2001 From: Oliver Sander <sander@igpm.rwth-aachen.de> Date: Mon, 6 Apr 2009 20:39:35 +0000 Subject: [PATCH] moved to new module dune-solvers [[Imported from SVN: r2360]] --- dune-solvers/norms/blocknorm.hh | 47 +++++++++++++++++ dune-solvers/norms/diagnorm.hh | 44 ++++++++++++++++ dune-solvers/norms/energynorm.hh | 79 ++++++++++++++++++++++++++++ dune-solvers/norms/fullnorm.hh | 45 ++++++++++++++++ dune-solvers/norms/h1seminorm.hh | 60 +++++++++++++++++++++ dune-solvers/norms/norm.hh | 21 ++++++++ dune-solvers/norms/pnorm.hh | 90 ++++++++++++++++++++++++++++++++ dune-solvers/norms/sumnorm.hh | 45 ++++++++++++++++ 8 files changed, 431 insertions(+) create mode 100644 dune-solvers/norms/blocknorm.hh create mode 100644 dune-solvers/norms/diagnorm.hh create mode 100644 dune-solvers/norms/energynorm.hh create mode 100644 dune-solvers/norms/fullnorm.hh create mode 100644 dune-solvers/norms/h1seminorm.hh create mode 100644 dune-solvers/norms/norm.hh create mode 100644 dune-solvers/norms/pnorm.hh create mode 100644 dune-solvers/norms/sumnorm.hh diff --git a/dune-solvers/norms/blocknorm.hh b/dune-solvers/norms/blocknorm.hh new file mode 100644 index 00000000..b6f88cc1 --- /dev/null +++ b/dune-solvers/norms/blocknorm.hh @@ -0,0 +1,47 @@ +#ifndef BLOCK_NORM_HH +#define BLOCK_NORM_HH + +#include <vector> + +#include "norm.hh" + +//! A norm for blocks of vectors +template <class VectorType> +class BlockNorm: public Norm<VectorType> +{ + public: + typedef typename std::vector<const Norm<typename VectorType::block_type>* > NormPointerVector; + + BlockNorm(const NormPointerVector& norms) : + norms_(norms) + {} + + //! Compute the norm of the given vector + double operator()(const VectorType &v) const + { + double r = 0.0; + for (int i=0; i<norms_.size(); ++i) + { + double ri = (*norms_[i])(v[i]); + r += ri*ri; + } + return sqrt(r); + } + + //! Compute the norm of the difference of two vectors + double diff(const VectorType &v1, const VectorType &v2) const + { + double r = 0.0; + for (int i=0; i<norms_.size(); ++i) + { + double ri = (*norms_[i]).diff(v1[i], v2[i]); + r += ri*ri; + } + return sqrt(r); + } + + private: + const NormPointerVector& norms_; +}; + +#endif diff --git a/dune-solvers/norms/diagnorm.hh b/dune-solvers/norms/diagnorm.hh new file mode 100644 index 00000000..1b4f5b8c --- /dev/null +++ b/dune-solvers/norms/diagnorm.hh @@ -0,0 +1,44 @@ +#ifndef __DIAGNORM__HH__ +#define __DIAGNORM__HH__ + +#include "norm.hh" + +typedef Dune::BlockVector<Dune::FieldVector <double,1> > VectorType; + +class DiagNorm: public Norm<VectorType> +{ + public: + DiagNorm(double alpha, const VectorType &d) : + d(d), + alpha(alpha) + {} + + //! Compute the norm of the given vector + double operator()(const VectorType &v) const { + double r = 0.0; + + for(int row = 0; row < v.size(); ++row) + r += d[row] * v[row] * v[row]; + + return sqrt(fabs(alpha * r)); + } + + //! Compute the norm of the difference of two vectors + double diff(const VectorType &v1, const VectorType &v2) const + { + double r = 0.0; + + for(int row = 0; row < v1.size(); ++row) + r += (double)d[row] * (v1[row]-v2[row]) * (v1[row] - v2[row]); + + return sqrt(fabs(alpha * r)); + } + + private: + const VectorType &d; + + const double alpha; + +}; + +#endif diff --git a/dune-solvers/norms/energynorm.hh b/dune-solvers/norms/energynorm.hh new file mode 100644 index 00000000..b24d2f7b --- /dev/null +++ b/dune-solvers/norms/energynorm.hh @@ -0,0 +1,79 @@ +#ifndef ENERGY_NORM_HH +#define ENERGY_NORM_HH + +#include "norm.hh" +#include <dune/ag-common/lineariterationstep.hh> + + template<class OperatorType, class DiscFuncType> + class EnergyNorm : public Norm<DiscFuncType> + { + public: + + /** \brief The type used for the result */ + typedef typename DiscFuncType::field_type field_type; + + EnergyNorm() : iterationStep_(NULL), matrix_(NULL) {} + + EnergyNorm(LinearIterationStep<OperatorType, DiscFuncType>& it) + : iterationStep_(&it), matrix_(NULL) + {} + + EnergyNorm(const OperatorType& matrix) + : iterationStep_(NULL), matrix_(&matrix) + {} + + void setMatrix(const OperatorType* matrix) { + matrix_ = matrix; + } + + void setIterationStep(LinearIterationStep<OperatorType, DiscFuncType>* iterationStep) { + iterationStep_ = iterationStep; + } + + //! Compute the norm of the difference of two vectors + double diff(const DiscFuncType& f1, const DiscFuncType& f2) const { + if (iterationStep_ == NULL && matrix_ == NULL) + DUNE_THROW(Dune::Exception, "You have not supplied neither a matrix nor an IterationStep to the EnergyNorm routine!"); + + DiscFuncType tmp_f = f1; + tmp_f -= f2; + return (*this)(tmp_f); + } + + //! Compute the norm of the given vector + double operator()(const DiscFuncType& f) const + { + if (iterationStep_ == NULL && matrix_ == NULL) + DUNE_THROW(Dune::Exception, "You have not supplied neither a matrix nor an IterationStep to the EnergyNorm routine!"); + + const OperatorType& A = (iterationStep_) + ? *(iterationStep_->getMatrix()) + : *matrix_; + + DiscFuncType tmp(f.size()); + tmp = 0; + A.umv(f, tmp); + + return sqrt(f * tmp); + } + + /** \brief Compute the squared norm for a given vector and matrix + \todo This could be implemented without the temporary. */ + static field_type normSquared(const DiscFuncType& u, + const OperatorType& A) + { + DiscFuncType tmp(u.size()); + tmp = 0; + A.umv(u, tmp); + return u*tmp; + } + + protected: + + LinearIterationStep<OperatorType, DiscFuncType>* iterationStep_; + + const OperatorType* matrix_; + + }; + +#endif diff --git a/dune-solvers/norms/fullnorm.hh b/dune-solvers/norms/fullnorm.hh new file mode 100644 index 00000000..a3559372 --- /dev/null +++ b/dune-solvers/norms/fullnorm.hh @@ -0,0 +1,45 @@ +#ifndef __FULLNORM__HH__ +#define __FULLNORM__HH__ + +#include "norm.hh" + +typedef Dune::BlockVector<Dune::FieldVector<double,1> > Vector; + +class FullNorm: public Norm<Vector> +{ + public: + FullNorm(double alpha, const Vector &m) : + m(m), + alpha(alpha) + {} + + //! Compute the norm of the given vector + double operator()(const Vector &v) const + { + double r = 0.0; + + for(int row = 0; row < v.size(); ++row) + r += m[row] * v[row]; + + return sqrt(fabs(alpha*r*r)); + } + + //! Compute the norm of the difference of two vectors + double diff(const Vector &v1, const Vector &v2) const + { + double r = 0.0; + + for(int row = 0; row < v1.size(); ++row) + r += (double)m[row] * (v1[row] - v2[row]); + + return sqrt(fabs(alpha*r*r)); + } + + private: + const Vector &m; + + const double alpha; + +}; + +#endif diff --git a/dune-solvers/norms/h1seminorm.hh b/dune-solvers/norms/h1seminorm.hh new file mode 100644 index 00000000..12e37562 --- /dev/null +++ b/dune-solvers/norms/h1seminorm.hh @@ -0,0 +1,60 @@ +#ifndef DUNE_H1_SEMINORM_HH +#define DUNE_H1_SEMINORM_HH + +#include "norm.hh" + +template<class VectorType> +class H1SemiNorm : public Norm<VectorType> +{ +public: + H1SemiNorm() : matrix_(NULL) {} + + H1SemiNorm(const Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> >& matrix) + : matrix_(&matrix) + {} + + //! Compute the norm of the difference of two vectors + double diff(const VectorType& u1, const VectorType& u2) const + { + double sum = 0; + + for (int i=0; i<matrix_->N(); i++) + { + typename Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> >::row_type::const_iterator cIt = (*matrix_)[i].begin(); + typename Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> >::row_type::const_iterator cEndIt = (*matrix_)[i].end(); + + typename VectorType::block_type differ_i = u1[i] - u2[i]; + for (; cIt!=cEndIt; ++cIt) + sum += differ_i * (u1[cIt.index()]-u2[cIt.index()]) * (*cIt); + } + + return std::sqrt(sum); + } + + //! Compute the norm of the given vector + double operator()(const VectorType& u) const + { + if (matrix_ == NULL) + DUNE_THROW(Dune::Exception, "You have not supplied neither a matrix nor an IterationStep to the EnergyNorm routine!"); + // we compute sqrt{uAu^t}. We have implemented this by hand because the matrix is + // always scalar but the vectors may not be + double sum = 0; + + for (int i=0; i<matrix_->N(); i++) { + + typename Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> >::row_type::const_iterator cIt = (*matrix_)[i].begin(); + typename Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> >::row_type::const_iterator cEndIt = (*matrix_)[i].end(); + + for (; cIt!=cEndIt; ++cIt) + sum += u[i]*u[cIt.index()] * (*cIt); + + } + + return std::sqrt(sum); + } + + const Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> >* matrix_; + +}; + +#endif diff --git a/dune-solvers/norms/norm.hh b/dune-solvers/norms/norm.hh new file mode 100644 index 00000000..324d22f0 --- /dev/null +++ b/dune-solvers/norms/norm.hh @@ -0,0 +1,21 @@ +#ifndef NORM_HH +#define NORM_HH + +//! Abstract base for classes computing norms of discrete functions +template <class VectorType> +class Norm { + + public: + + /** \brief Destructor, doing nothing */ + virtual ~Norm() {}; + + //! Compute the norm of the given vector + virtual double operator()(const VectorType& f) const = 0; + + //! Compute the norm of the difference of two vectors + virtual double diff(const VectorType& f1, const VectorType& f2) const = 0; + +}; + +#endif diff --git a/dune-solvers/norms/pnorm.hh b/dune-solvers/norms/pnorm.hh new file mode 100644 index 00000000..e728b8f0 --- /dev/null +++ b/dune-solvers/norms/pnorm.hh @@ -0,0 +1,90 @@ +#ifndef __PNORM__HH__ +#define __PNORM__HH__ + +#include "norm.hh" + +typedef BlockVector<FieldVector <double,1> > Vector; + +class PNorm: public Norm<Vector> +{ + public: + PNorm(int p=2, double alpha=1.0): + p(p), + alpha(alpha) + {} + + double operator()(const Vector &v) const + { + double r = 0.0; + + if (p<1) + { + for(int row = 0; row < v.size(); ++row) + { + double z = fabs(v[row]); + if (r<z) + r = z; + } + } + else if (p==1) + { + for(int row = 0; row < v.size(); ++row) + r += fabs(v[row]); + } + else if (p==2) + { + for(int row = 0; row < v.size(); ++row) + r += v[row] * v[row]; + r = sqrt(r); + } + else + { + for(int row = 0; row < v.size(); ++row) + r += pow(fabs(v[row]), p); + r = pow(r, 1.0/p); + } + + return alpha*r; + } + + double diff(const Vector &v1, const Vector &v2) const + { + double r = 0.0; + + if (p<1) + { + for(int row = 0; row < v1.size(); ++row) + { + double z = fabs(v1[row]-v2[row]); + if (z>r) + r = z; + } + } + else if (p==1) + { + for(int row = 0; row < v1.size(); ++row) + r += fabs(v1[row]-v2[row]); + } + else if (p==2) + { + for(int row = 0; row < v1.size(); ++row) + r += (v1[row]-v2[row]) * (v1[row]-v2[row]); + r = sqrt(r); + } + else + { + for(int row = 0; row < v1.size(); ++row) + r += pow(fabs(v1[row]-v2[row]), p); + r = pow(r, 1.0/p); + } + + return alpha*r; + } + + private: + double alpha; + double p; + +}; + +#endif diff --git a/dune-solvers/norms/sumnorm.hh b/dune-solvers/norms/sumnorm.hh new file mode 100644 index 00000000..95e4fb0e --- /dev/null +++ b/dune-solvers/norms/sumnorm.hh @@ -0,0 +1,45 @@ +#ifndef __SUMNORM__HH__ +#define __SUMNORM__HH__ + +#include "norm.hh" + +template <class VectorType> +class SumNorm: public Norm<VectorType> +{ + public: + SumNorm(double _alpha1, const Norm<VectorType> &_norm1, double _alpha2, const Norm<VectorType> &_norm2) : + alpha1(_alpha1), + norm1(_norm1), + alpha2(_alpha2), + norm2(_norm2) + {} + + //! Compute the norm of the given vector + double operator()(const VectorType &v) const + { + double r1 = norm1(v); + double r2 = norm2(v); + + return sqrt(alpha1 * r1 * r1 + alpha2 * r2 *r2); + } + + //! Compute the norm of the difference of two vectors + double diff(const VectorType &v1, const VectorType &v2) const + { + double r1 = norm1.diff(v1,v2); + double r2 = norm2.diff(v1,v2); + + return sqrt(alpha1 * r1 * r1 + alpha2 * r2 *r2); + } + + private: + const double alpha1; + const Norm<VectorType> &norm1; + + const double alpha2; + const Norm<VectorType> &norm2; + + +}; + +#endif -- GitLab