moved to new module dune-solvers

[[Imported from SVN: r2360]]

dune-solvers/norms/blocknorm.hh

+#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

dune-solvers/norms/diagnorm.hh

+#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 

dune-solvers/norms/energynorm.hh

+#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

dune-solvers/norms/fullnorm.hh

+#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 

dune-solvers/norms/h1seminorm.hh

+#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

dune-solvers/norms/norm.hh

+#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

dune-solvers/norms/pnorm.hh

+#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 

dune-solvers/norms/sumnorm.hh

+#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