Add exceptions when we cannot downcast a VectorSpaceVector object to Vector.

diff --git a/include/deal.II/lac/la_vector.h b/include/deal.II/lac/la_vector.h
index 20d972f2a5bcacd5efa927aa6a3bcdfba48f20e6..0e0ea0acb7dde5fcee18a33f8d3f4af5ba67bafd 100644 (file)
--- a/include/deal.II/lac/la_vector.h
+++ b/include/deal.II/lac/la_vector.h
@@ -51,6 +51,8 @@ namespace LinearAlgebra
   class Vector : public ReadWriteVector<Number>, public VectorSpaceVector<Number>
   {
   public:
+    typedef types::global_dof_index size_type;
+
     /**
      * Constructor. Create a vector of dimension zero.
      */
@@ -72,11 +74,11 @@ namespace LinearAlgebra
      * <tt>v=Vector@<Number@>(0);</tt>, i.e. the vector is replaced by one of
      * length zero.
      */
-    explicit Vector(const typename ReadWriteVector<Number>::size_type n);
+    explicit Vector(const size_type n);
 
     /**
      * Initialize the vector with a given range of values pointed to by the
-     * iterators. This function is there in analogy to the @p std::vector class.
+     * iterators. This function exists in analogy to the @p std::vector class.
      */
     template <typename InputIterator>
     Vector(const InputIterator first, const InputIterator last);
@@ -84,88 +86,88 @@ namespace LinearAlgebra
     /**
      * Destructor, deallocates memory.
      */
-    ~Vector();
+    virtual ~Vector();
 
     /**
      * Multiply the entire vector by a fixed factor.
      */
-    VectorSpaceVector<Number> &operator*= (const Number factor);
+    virtual VectorSpaceVector<Number> &operator*= (const Number factor);
 
     /**
      * Divide the entire vector by a fixed factor.
      */
-    VectorSpaceVector<Number> &operator/= (const Number factor);
+    virtual VectorSpaceVector<Number> &operator/= (const Number factor);
 
     /**
      * Add the vector @p V to the present one.
      */
-    VectorSpaceVector<Number> &operator+= (const VectorSpaceVector<Number> &V);
+    virtual VectorSpaceVector<Number> &operator+= (const VectorSpaceVector<Number> &V);
 
     /**
      * Substract the vector @p V from the present one.
      */
-    VectorSpaceVector<Number> &operator-= (const VectorSpaceVector<Number> &V);
+    virtual VectorSpaceVector<Number> &operator-= (const VectorSpaceVector<Number> &V);
 
     /**
      * Return the scalar product of two vectors.
      */
-    Number operator* (const VectorSpaceVector<Number> &V);
+    virtual Number operator* (const VectorSpaceVector<Number> &V);
 
     /**
      * Simple addition of a multiple of a vector, i.e. <tt>*this += a*V</tt>.
      */
-    void add(const Number a, const VectorSpaceVector<Number> &V);
+    virtual void add(const Number a, const VectorSpaceVector<Number> &V);
 
     /**
      * Multiple addition of a multiple of a vector, i.e. <tt>*this += a*V+b*W</tt>.
      */
-    void add(const Number a, const VectorSpaceVector<Number> &V,
-             const Number b, const VectorSpaceVector<Number> &W);
+    virtual void add(const Number a, const VectorSpaceVector<Number> &V,
+                     const Number b, const VectorSpaceVector<Number> &W);
 
     /**
      * Scaling and simple vector addition, i.e. <tt>*this = s*(*this)+V</tt>.
      */
-    void sadd(const Number s, const VectorSpaceVector<Number> &V);
+    virtual void sadd(const Number s, const VectorSpaceVector<Number> &V);
 
     /**
      * Scaling and simple addition of a multiple of a vector, i.e. <tt>*this =
      * s*(*this)+a*V<tt>.
      */
-    void sadd(const Number s, const Number a,
-              const VectorSpaceVector<Number> &V);
+    virtual void sadd(const Number s, const Number a,
+                      const VectorSpaceVector<Number> &V);
 
     /**
      * Scale each element of this vector by the corresponding element in the
      * argument. This function is mostly meant to simulate multiplication (and
      * immediate re-assignement) by a diagonal scaling matrix.
      */
-    void scale(const VectorSpaceVector<Number> &scaling_factors);
+    virtual void scale(const VectorSpaceVector<Number> &scaling_factors);
 
     /**
      * Assignement <tt>*this = a*V</tt>.
      */
-    void equ(const Number a, const VectorSpaceVector<Number> &V);
+    virtual void equ(const Number a, const VectorSpaceVector<Number> &V);
 
     /**
-     * Returns the l<sub>1</sub> norm of the vector (i.e., the sum of the
+     * Return the l<sub>1</sub> norm of the vector (i.e., the sum of the
      * absolute values of all entries).
      */
-    typename VectorSpaceVector<Number>::real_type l1_norm();
+    virtual typename VectorSpaceVector<Number>::real_type l1_norm();
 
     /**
-     * Returns the l<sub>2</sub> norm of the vector (i.e., the square root of
+     * Return the l<sub>2</sub> norm of the vector (i.e., the square root of
      * the sum of the square of all entries among all processors).
      */
-    typename VectorSpaceVector<Number>::real_type l2_norm();
+    virtual typename VectorSpaceVector<Number>::real_type l2_norm();
 
     /**
-     * Returns the maximum norm of the vector (i.e., the maximum absolute
+     * Return the maximum norm of the vector (i.e., the maximum absolute
      * value among all entries and among all processors).
      */
-    typename VectorSpaceVector<Number>::real_type linfty_norm();
+    virtual typename VectorSpaceVector<Number>::real_type linfty_norm();
 
     /**
-     * Performs a combined operation of a vector addition and a subsequent
+     * Perform a combined operation of a vector addition and a subsequent
      * inner product, returning the value of the inner product. In other
      * words, the result of this function is the same as if the user called
      * @code
@@ -173,18 +175,18 @@ namespace LinearAlgebra
      * return_value = *this * W;
      * @endcode
      */
-    Number add_and_dot(const Number a,
-                       const VectorSpaceVector<Number> &V,
-                       const VectorSpaceVector<Number> &W);
+    virtual Number add_and_dot(const Number a,
+                               const VectorSpaceVector<Number> &V,
+                               const VectorSpaceVector<Number> &W);
 
     /**
-     * Returns the global size of the vector, equal to the sum of the number
+     * Return the global size of the vector, equal to the sum of the number
      * of locally owned indices among all processors.
      */
-    typename ReadWriteVector<Number>::size_type size() const;
+    size_type size() const;
 
     /**
-     * Returns an index set that describes which elements of this vector are
+     * Return an index set that describes which elements of this vector are
      * owned by the current processor. As a consequence, the index sets returned
      * on different procesors if this is a distributed vector will form disjoint
      * sets that add up to the complete index set. Obviously, if a vector is
@@ -198,16 +200,23 @@ namespace LinearAlgebra
     /**
      * Prints the vector to the output stream @p out.
      */
-    void print(std::ostream &out,
-               const unsigned int precision=3,
-               const bool scientific=true,
-               const bool across=true) const;
+    virtual void print(std::ostream &out,
+                       const unsigned int precision=3,
+                       const bool scientific=true,
+                       const bool across=true) const;
 
     /**
      * Returns the memory consumption of this class in bytes.
      */
     std::size_t memory_consumption() const;
 
+    /**
+     * Attempt to perform an operation between two incompatible vector types.
+     *
+     * @ingroup Exceptions
+     */
+    DeclException0(ExcVectorTypeNotCompatible);
+
   private:
     /**
      * Compute the L1 norm in a recursive way by dividing the vector on smaller
@@ -217,11 +226,13 @@ namespace LinearAlgebra
         unsigned int j);
 
     /**
-     * Compute the L2 norm in a recursive way by dividing the vector on smaller
-     * and smaller intervals. This reduces the numerical error on large vector.
+     * Compute the squared L2 norm in a recursive way by dividing the vector on
+     * smaller and smaller intervals. This reduces the numerical error on large
+     * vector.
      */
-    typename VectorSpaceVector<Number>::real_type l2_norm_recursive(unsigned int i,
-        unsigned int j);
+    typename VectorSpaceVector<Number>::real_type l2_norm_squared_recursive(
+      unsigned int i,
+      unsigned int j);
   };
 
   /*@}*/
@@ -244,7 +255,7 @@ namespace LinearAlgebra
 
   template <typename Number>
   inline
-  Vector<Number>::Vector(const typename ReadWriteVector<Number>::size_type n)
+  Vector<Number>::Vector(const size_type n)
     :
     ReadWriteVector<Number>(n)
   {}
@@ -270,7 +281,7 @@ namespace LinearAlgebra
 
   template <typename Number>
   inline
-  typename ReadWriteVector<Number>::size_type Vector<Number>::size() const
+  typename Vector<Number>::size_type Vector<Number>::size() const
   {
     return ReadWriteVector<Number>::size();
   }

diff --git a/include/deal.II/lac/la_vector.templates.h b/include/deal.II/lac/la_vector.templates.h
index 9057e5424989bcadb3a3cf68f6bd5ef83fe5fd5e..f436504d48dabd603906271c1fd4dac225d569b4 100644 (file)
--- a/include/deal.II/lac/la_vector.templates.h
+++ b/include/deal.II/lac/la_vector.templates.h
@@ -49,6 +49,9 @@ namespace LinearAlgebra
   template <typename Number>
   VectorSpaceVector<Number> &Vector<Number>::operator+= (const VectorSpaceVector<Number> &V)
   {
+    // Check that casting will work.
+    Assert(dynamic_cast<const Vector<Number>*>(&V)!=NULL, ExcVectorTypeNotCompatible());
+
     // Downcast V. If fails, throws an exception.
     const Vector<Number> &down_V = dynamic_cast<const Vector<Number>&>(V);
     Assert(down_V.size()==this->size(),
@@ -67,6 +70,9 @@ namespace LinearAlgebra
   template <typename Number>
   VectorSpaceVector<Number> &Vector<Number>::operator-= (const VectorSpaceVector<Number> &V)
   {
+    // Check that casting will work.
+    Assert(dynamic_cast<const Vector<Number>*>(&V)!=NULL, ExcVectorTypeNotCompatible());
+
     // Downcast V. If fails, throws an exception.
     const Vector<Number> &down_V = dynamic_cast<const Vector<Number>&>(V);
     Assert(down_V.size()==this->size(),
@@ -85,6 +91,10 @@ namespace LinearAlgebra
   template <typename Number>
   Number Vector<Number>::operator* (const VectorSpaceVector<Number> &V)
   {
+    // Check that casting will work.
+    Assert(dynamic_cast<const Vector<Number>*>(&V)!=NULL,
+           ExcVectorTypeNotCompatible());
+
     // Downcast V. If fails, throws an exception.
     const Vector<Number> &down_V = dynamic_cast<const Vector<Number>&>(V);
     Assert(down_V.size()==this->size(),
@@ -105,6 +115,10 @@ namespace LinearAlgebra
   template <typename Number>
   void Vector<Number>::add(const Number a, const VectorSpaceVector<Number> &V)
   {
+    // Check that casting will work.
+    Assert(dynamic_cast<const Vector<Number>*>(&V)!=NULL,
+           ExcVectorTypeNotCompatible());
+
     // Downcast V. If fails, throws an exception.
     const Vector<Number> &down_V = dynamic_cast<const Vector<Number>&>(V);
     AssertIsFinite(a);
@@ -123,6 +137,13 @@ namespace LinearAlgebra
   void Vector<Number>::add(const Number a, const VectorSpaceVector<Number> &V,
                            const Number b, const VectorSpaceVector<Number> &W)
   {
+    // Check that casting will work.
+    Assert(dynamic_cast<const Vector<Number>*>(&V)!=NULL,
+           ExcVectorTypeNotCompatible());
+    // Check that casting will work.
+    Assert(dynamic_cast<const Vector<Number>*>(&W)!=NULL,
+           ExcVectorTypeNotCompatible());
+
     // Downcast V. If fails, throws an exception.
     const Vector<Number> &down_V = dynamic_cast<const Vector<Number>&>(V);
     // Downcast W. If fails, throws an exception.
@@ -156,6 +177,10 @@ namespace LinearAlgebra
   void Vector<Number>::sadd(const Number s, const Number a,
                             const VectorSpaceVector<Number> &V)
   {
+    // Check that casting will work.
+    Assert(dynamic_cast<const Vector<Number>*>(&V)!=NULL,
+           ExcVectorTypeNotCompatible());
+
     *this *= s;
     // Downcast V. It fails, throws an exception.
     const Vector<Number> &down_V = dynamic_cast<const Vector<Number>&>(V);
@@ -169,6 +194,10 @@ namespace LinearAlgebra
   template <typename Number>
   void Vector<Number>::scale(const VectorSpaceVector<Number> &scaling_factors)
   {
+    // Check that casting will work.
+    Assert(dynamic_cast<const Vector<Number>*>(&scaling_factors)!=NULL,
+           ExcVectorTypeNotCompatible());
+
     // Downcast scaling_factors. If fails, throws an exception.
     const Vector<Number> &down_scaling_factors =
       dynamic_cast<const Vector<Number>&>(scaling_factors);
@@ -186,6 +215,10 @@ namespace LinearAlgebra
   template <typename Number>
   void Vector<Number>::equ(const Number a, const VectorSpaceVector<Number> &V)
   {
+    // Check that casting will work.
+    Assert(dynamic_cast<const Vector<Number>*>(&V)!=NULL,
+           ExcVectorTypeNotCompatible());
+
     // Downcast V. If fails, throws an exception.
     const Vector<Number> &down_V = dynamic_cast<const Vector<Number>&>(V);
     *this = down_V;
@@ -220,7 +253,7 @@ namespace LinearAlgebra
     // precision) using the BLAS approach with a weight, see e.g. dnrm2.f.
     typedef typename ReadWriteVector<Number>::real_type real_type;
     real_type norm_square = 0.;
-    norm_square = l2_norm_recursive(0,this->size()-1);
+    norm_square = l2_norm_squared_recursive(0,this->size()-1);
     if (numbers::is_finite(norm_square) &&
         norm_square>=std::numeric_limits<real_type>::min())
       return std::sqrt(norm_square);
@@ -255,8 +288,7 @@ namespace LinearAlgebra
   {
     typename ReadWriteVector<Number>::real_type norm = 0.;
     for (unsigned int i=0; i<this->size(); ++i)
-      if (std::abs(this->val[i])>norm)
-        norm = std::abs(this->val[i]);
+      norm = std::max(std::abs(this->val[i]),norm);
 
     return norm;
   }
@@ -276,8 +308,9 @@ namespace LinearAlgebra
 
 
   template <typename Number>
-  typename VectorSpaceVector<Number>::real_type Vector<Number>::l1_norm_recursive(unsigned int i,
-      unsigned int j)
+  typename VectorSpaceVector<Number>::real_type Vector<Number>::l1_norm_recursive(
+    unsigned int i,
+    unsigned int j)
   {
     Assert(j>=i, ExcInternalError());
     typename ReadWriteVector<Number>::real_type norm = 0.;
@@ -296,16 +329,17 @@ namespace LinearAlgebra
 
 
   template <typename Number>
-  typename VectorSpaceVector<Number>::real_type Vector<Number>::l2_norm_recursive(unsigned int i,
-      unsigned int j)
+  typename VectorSpaceVector<Number>::real_type Vector<Number>::l2_norm_squared_recursive(
+    unsigned int i,
+    unsigned int j)
   {
     Assert(j>=i, ExcInternalError());
 
     typename ReadWriteVector<Number>::real_type norm = 0.;
     if ((j-i)!=0)
       {
-        norm += l2_norm_recursive(i,(i+j)/2);
-        norm += l2_norm_recursive((i+j)/2+1,j);
+        norm += l2_norm_squared_recursive(i,(i+j)/2);
+        norm += l2_norm_squared_recursive((i+j)/2+1,j);
       }
     else
       norm += std::pow(std::abs(this->val[i]),2);

diff --git a/include/deal.II/lac/vector_space_vector.h b/include/deal.II/lac/vector_space_vector.h
index 60e47ee77164a60fd45645841b4108b65c998694..31ba1b1c7263f1b18793b79aff9f2983612a5d28 100644 (file)
--- a/include/deal.II/lac/vector_space_vector.h
+++ b/include/deal.II/lac/vector_space_vector.h
@@ -113,25 +113,25 @@ namespace LinearAlgebra
     virtual void equ(const Number a, const VectorSpaceVector<Number> &V) = 0;
 
     /**
-     * Returns the l<sub>1</sub> norm of the vector (i.e., the sum of the
+     * Return the l<sub>1</sub> norm of the vector (i.e., the sum of the
      * absolute values of all entries among all processors).
      */
     virtual real_type l1_norm() = 0;
 
     /**
-     * Returns the l<sub>2</sub> norm of the vector (i.e., the square root of
+     * Return the l<sub>2</sub> norm of the vector (i.e., the square root of
      * the sum of the square of all entries among all processors).
      */
     virtual real_type l2_norm() = 0;
 
     /**
-     * Returns the maximum norm of the vector (i.e., the maximum absolute value
+     * Return the maximum norm of the vector (i.e., the maximum absolute value
      * among all entries and among all processors).
      */
     virtual real_type linfty_norm() = 0;
 
     /**
-     * Performs a combined operation of a vector addition and a subsequent
+     * Perform a combined operation of a vector addition and a subsequent
      * inner product, returning the value of the inner product. In other
      * words, the result of this function is the same as if the user called
      * @code
@@ -151,13 +151,13 @@ namespace LinearAlgebra
                                const VectorSpaceVector<Number> &W) = 0;
 
     /**
-     * Returns the global size of the vector, equal to the sum of the number of
+     * Return the global size of the vector, equal to the sum of the number of
      * locally owned indices among all processors.
      */
     virtual size_type size() const = 0;
 
     /**
-     * Returns an index set that describes which elements of this vector are
+     * Return an index set that describes which elements of this vector are
      * owned by the current processor. As a consequence, the index sets returned
      * on different procesors if this is a distributed vector will form disjoint
      * sets that add up to the complete index set. Obviously, if a vector is
@@ -169,7 +169,7 @@ namespace LinearAlgebra
     virtual const dealii::IndexSet &locally_owned_elements() const = 0;
 
     /**
-     * Prints the vector to the output stream @p out.
+     * Print the vector to the output stream @p out.
      */
     virtual void print(std::ostream &out,
                        const unsigned int precision=3,
@@ -177,7 +177,7 @@ namespace LinearAlgebra
                        const bool across=true) const = 0;
 
     /**
-     * Returns the memory consumption of this class in bytes.
+     * Return the memory consumption of this class in bytes.
      */
     virtual std::size_t memory_consumption() const = 0;
   };