gradientOf(DV<WorldVector>) corrected

AMDiS/src/DOFVector.cc

@@ -596,6 +596,7 @@ namespace AMDiS {
 
     // define result vector
     static WorldVector<DOFVector<double>*> *result = NULL; // TODO: REMOVE STATIC
+    DOFVector<double>* null_ptr = NULL;
 
     if (grad) {
       result = grad;
@@ -603,7 +604,7 @@ namespace AMDiS {
       if (!result) {
 	result = new WorldVector<DOFVector<double>*>;
 
-	result->set(NULL);
+	result->set(null_ptr);
       }
       for (int i = 0; i < dow; i++) {
 	if ((*result)[i] && (*result)[i]->getFeSpace() != feSpace) {

AMDiS/src/DOFVector.hh

@@ -1402,7 +1402,7 @@ namespace AMDiS {
       vecAtQPs.change_dim(nPoints);
 
       for (int iq = 0; iq < nPoints; iq++) {
-       	vecAtQPs[iq] = 0.0;
+       	nullify(vecAtQPs[iq]);
        	for (int j = 0; j < nBasFcts; j++)
        	  vecAtQPs[iq] += 
        	    localVec[j] * (*(basFcts->getPhi(j)))(quad->getLambda(iq));
@@ -1440,7 +1440,7 @@ namespace AMDiS {
     mtl::dense2D<double> &m = smallElInfo->getSubElemCoordsMat(basFcts->getDegree());
 
     for (int iq = 0; iq < nPoints; iq++) {
-      vecAtQPs[iq] = 0.0;
+      nullify(vecAtQPs[iq]);
       for (int j = 0; j < nBasFcts; j++) {
 	double val = 0.0;
 	for (int k = 0; k < nBasFcts; k++)

AMDiS/src/FixVec.h

@@ -412,7 +412,9 @@ namespace AMDiS {
     {}
 
     /// Sets all entries to d
-    inline const WorldVector<T>& operator=(const T& d)
+    template <typename S>
+    inline const WorldVector<T>&
+    operator=(const S& d)
     {
       this->set(d);
       return (*this);
@@ -643,7 +645,25 @@ namespace AMDiS {
     M1 -= M2;
     return M1;
   }
+
+  template<typename T>
+  WorldVector<WorldVector<T> > 
+  operator*(const WorldVector<WorldVector<T> >& A, const WorldVector<WorldVector<T> >& B)
+  {
+    WorldVector<WorldVector<T> > result;
+    nullify(result);
+    for (size_t r = 0; r < num_rows(A); r++)
+      for (size_t c = 0; c < num_cols(A); c++)
+	for (size_t i = 0; i < num_cols(A); i++)
+	  result[r][c] += A[r][i] * B[i][c];
+    return result;
+  }
   
+  template<typename T>
+  void set_to_zero(WorldVector<WorldVector<T> >& mat)
+  {
+    nullify(mat);
+  }
 
   template<typename T>
   WorldVector<T> operator-(WorldVector<T> v)
@@ -715,6 +735,45 @@ namespace AMDiS {
 }
 
 
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/utility/ashape.hpp>
+#include <boost/numeric/mtl/utility/category.hpp>
+#include <boost/numeric/mtl/utility/is_row_major.hpp>
+
+namespace mtl 
+{
+  template <typename T>
+  struct Collection<AMDiS::WorldVector<AMDiS::WorldVector<T> > >
+  {
+    typedef T      value_type;
+    typedef int    size_type;
+  };
+
+  namespace ashape 
+  {
+    template <typename T> 
+    struct ashape_aux<AMDiS::WorldVector<AMDiS::WorldVector<T> > > 
+    { 
+      typedef mat<typename ashape<T>::type> type;  // nonscal
+    };
+  } // end namespace ashape  
+  
+  namespace traits
+  {
+    template <typename T>
+    struct category<AMDiS::WorldVector<AMDiS::WorldVector<T> > > 
+    {
+	typedef tag::dense2D type;
+    };
+    
+    template <typename T>
+    struct is_row_major<AMDiS::WorldVector<AMDiS::WorldVector<T> > >
+      : boost::mpl::true_ {};
+  } // end namespace traits
+} // end namespace mtl
+
+
+
 #include "FixVec.hh"
 
 #endif // AMDIS_FIXVEC_H

AMDiS/src/MatrixVector.h

@@ -101,7 +101,9 @@ namespace AMDiS {
     }
 
     /// Assignement operator
-    inline const Vector<T>& operator=(const T& scal) 
+    template <typename S>
+    inline const Vector<T>&
+    operator=(const S& scal) 
     {
       for (T *thisIt = this->begin(); thisIt != this->end(); ++thisIt)
 	*thisIt = scal;
@@ -123,10 +125,18 @@ namespace AMDiS {
     }
 
     /// Sets all entries to scal.
+    template <typename S>
+    inline typename boost::disable_if<boost::mpl::and_<boost::is_integral<S>, boost::is_pointer<T> >, const Vector<T>& >::type
+    set(const S& scal) 
+    {
+      return *this = scal;
+    }
+    
     inline const Vector<T>& set(const T& scal) 
     {
       return *this = scal;
     }
+    
 
     /// Sets all entries.
     inline const Vector<T>& setValues(const T* values) 

AMDiS/src/Traits.h

@@ -134,14 +134,14 @@ namespace AMDiS
     
     /// Scalar*Vector => Vector
     template<template<class> class Container, typename T1, typename T2>
-    struct mult_type_dispatch<T1, Container<T2>, tag::scalar, tag::matrix>
+    struct mult_type_dispatch<T1, Container<T2>, tag::scalar, tag::vector>
     {
       typedef typename mult_type<T1, T2>::type value_type;
       typedef Container<value_type> type;
     };
     
     
-    /// Vec*Scalar => Vector
+    /// Matrix*Scalar => Matrix
     template<template<class> class Container, typename T1, typename T2>
     struct mult_type_dispatch<Container<T1>, T2, tag::matrix, tag::scalar>
     {
@@ -149,9 +149,9 @@ namespace AMDiS
       typedef Container<value_type> type;
     };
     
-    /// Scalar*Vector => Vector
+    /// Scalar*Matrix => Matrix
     template<template<class> class Container, typename T1, typename T2>
-    struct mult_type_dispatch<T1, Container<T2>, tag::scalar, tag::vector>
+    struct mult_type_dispatch<T1, Container<T2>, tag::scalar, tag::matrix>
     {
       typedef typename mult_type<T1, T2>::type value_type;
       typedef Container<value_type> type;

AMDiS/src/expressions/functor_expr.hpp

@@ -245,7 +245,7 @@ namespace AMDiS
   // call a function with 3 arguments
   // _____________________________________________________________________________
   template<typename F, typename Term1, typename Term2, typename Term3>
-  inline typename result_of::Function3<F, Term1, Term2, Term2>::type
+  inline typename result_of::Function3<F, Term1, Term2, Term3>::type
   function_(const F& f, const Term1& t1, const Term2& t2, const Term3& t3) 
   { 
     typedef typename traits::to_expr<Term1>::to Expr1;
@@ -257,7 +257,7 @@ namespace AMDiS
 
 
   template<typename F, typename Term1, typename Term2, typename Term3>
-  inline typename result_of::Function3<F, Term1, Term2, Term2>::type
+  inline typename result_of::Function3<F, Term1, Term2, Term3>::type
   function_(const Term1& t1, const Term2& t2, const Term3& t3) 
   { 
     typedef typename traits::to_expr<Term1>::to Expr1;

AMDiS/src/expressions/gradientOf.hpp

@@ -380,4 +380,24 @@ namespace AMDiS
 } // end namespace AMDiS
 
 
+// ------- something special needed for gradientOf(DOFVector<WorldVector>)
+#include <boost/numeric/mtl/operation/mult_result.hpp>
+namespace mtl {
+  namespace traits {
+  
+    typedef AMDiS::WorldVector<AMDiS::WorldVector<double> > WWMatrix;
+  
+    template <typename Op1, typename IsMatrix>
+    struct mult_result_WWMatrix {
+	typedef mat_cvec_times_expr<Op1, mtl::dense_vector<WWMatrix> > type;
+    };
+  
+    /// Multiply matrix with column vector
+    template <typename Op1>
+    struct mult_result<Op1, mtl::dense_vector<WWMatrix> > 
+      : public mult_result_WWMatrix<Op1, typename boost::enable_if<is_matrix<Op1> >::type >
+    {};
+  }
+}
+
 #endif // AMDIS_GRADIENT_OF_HPP

AMDiS/src/traits/scalar_types.hpp

+/******************************************************************************
+ *
+ * AMDiS - Adaptive multidimensional simulations
+ *
+ * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
+ * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
+ *
+ * Authors: 
+ * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
+ *
+ * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+ * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+ *
+ *
+ * This file is part of AMDiS
+ *
+ * See also license.opensource.txt in the distribution.
+ * 
+ ******************************************************************************/
+
+
+
+/** \file scalar_types.hpp */
+
+#ifndef AMDIS_TYPE_TRAITS_SCALAR_TYPES_HPP
+#define AMDIS_TYPE_TRAITS_SCALAR_TYPES_HPP
+
+namespace AMDiS 
+{
+  namespace traits 
+  {
+	
+    // type-traits
+    // _________________________________________________________________________
+    template<typename T>
+    struct is_integer : boost::mpl::or_
+      <
+	typename boost::is_signed<T>::type,
+	typename boost::is_unsigned<T>::type
+      >::type {};
+    
+    template<typename T>
+    struct is_numeric : boost::mpl::or_
+      <
+	typename boost::is_floating_point<T>::type,
+	typename is_integer<T>::type
+      >::type {};
+      
+    
+  } // end namespace traits
+  
+} // end namespace AMDiS
+
+#endif // AMDIS_TYPE_TRAITS_TYPES_HPP

AMDiS/src/traits/types.hpp

@@ -37,6 +37,8 @@
 
 #include "boost/numeric/ublas/detail/returntype_deduction.hpp"
 
+#include "scalar_types.hpp"
+
 namespace AMDiS 
 {
   namespace traits 
@@ -44,23 +46,6 @@ namespace AMDiS
   
     // dummy type
     typedef boost::numeric::ublas::error_cant_deduce_type no_valid_type;
-	
-    // type-traits
-    // _________________________________________________________________________
-    template<typename T>
-    struct is_integer : boost::mpl::or_
-      <
-	typename boost::is_signed<T>::type,
-	typename boost::is_unsigned<T>::type
-      >::type {};
-    
-    template<typename T>
-    struct is_numeric : boost::mpl::or_
-      <
-	typename boost::is_floating_point<T>::type,
-	typename is_integer<T>::type
-      >::type {};
-      
       
     // test for mtl4 types
     // _________________________________________________________________________