diff --git a/dune/fufem/arithmetic.hh b/dune/fufem/arithmetic.hh
index 9272244aa5f192ee2d3ac1fae4554dd7c2db094c..86c5f20fc083230f5e9dacda473b601c99345b7c 100644
--- a/dune/fufem/arithmetic.hh
+++ b/dune/fufem/arithmetic.hh
@@ -11,6 +11,8 @@
#include <dune/istl/matrixindexset.hh>
#include <dune/istl/scaledidmatrix.hh>
+#include <dune/matrix-vector/axpy.hh>
+
/** \brief Namespace containing helper classes and functions for arithmetic operations
*
* Everything in this namespace is experimental and might change
@@ -138,386 +140,6 @@ namespace Arithmetic
template <bool Condition, typename Type=void>
using enable_if_t = typename std::enable_if<Condition, Type>::type;
- /** \brief Internal helper class for product operations
- *
- */
- template<class A, class B, class C, bool AisScalar, bool BisScalar, bool CisScalar>
- struct ProductHelper
- {
- template <class ADummy=A, enable_if_t<!MatrixTraits<ADummy>::isMatrix,int> SFINAE_Dummy=0>
- static void addProduct(A& a, const B& b, const C& c)
- {
- b.umv(c, a);
- }
-
- template <class ADummy=A, enable_if_t<MatrixTraits<ADummy>::isMatrix,int> SFINAE_Dummy=0>
- static void addProduct(A& a, const B& b, const C& c)
- {
- typename B::ConstRowIterator bi = b.begin();
- typename B::ConstRowIterator bEnd = b.end();
- for(; bi!=bEnd; ++bi)
- {
- typename B::ConstColIterator bik = bi->begin();
- typename B::ConstColIterator biEnd= bi->end();
- for(; bik!=biEnd; ++bik)
- {
- typename C::ConstColIterator ckj =c[bik.index()].begin();
- typename C::ConstColIterator ckEnd=c[bik.index()].end();
- for(; ckj!=ckEnd; ++ckj)
- a[bi.index()][ckj.index()] += (*bik) * (*ckj);
- }
- }
- }
- };
-
- // Internal helper class for scaled product operations (i.e., b is always a scalar)
- template<class A, class Scalar, class B, class C, bool AisScalar, bool BisScalar, bool CisScalar>
- struct ScaledProductHelper
- {
- template <class ADummy=A, enable_if_t<!MatrixTraits<ADummy>::isMatrix,int> SFINAE_Dummy=0>
- static void addProduct(A& a, const Scalar& scalar, const B& b, const C& c)
- {
- b.usmv(scalar, c, a);
- }
-
- template <class ADummy=A, enable_if_t<MatrixTraits<ADummy>::isMatrix,int> SFINAE_Dummy=0>
- static void addProduct(A& a, const Scalar& scalar, const B& b, const C& c)
- {
- typename B::ConstRowIterator bi = b.begin();
- typename B::ConstRowIterator bEnd = b.end();
- for(; bi!=bEnd; ++bi)
- {
- typename B::ConstColIterator bik =b[bi.index()].begin();
- typename B::ConstColIterator biEnd=b[bi.index()].end();
- for(; bik!=biEnd; ++bik)
- {
- typename C::ConstColIterator ckj =c[bik.index()].begin();
- typename C::ConstColIterator ckEnd=c[bik.index()].end();
- for(; ckj!=ckEnd; ++ckj)
- a[bi.index()][ckj.index()] += scalar * (*bik) * (*ckj);
- }
- }
- }
- };
-
- template<class T, int n, int m, int k>
- struct ProductHelper<Dune::FieldMatrix<T,n,k>, Dune::FieldMatrix<T,n,m>, Dune::FieldMatrix<T,m,k>, false, false, false>
- {
- typedef Dune::FieldMatrix<T,n,k> A;
- typedef Dune::FieldMatrix<T,n,m> B;
- typedef Dune::FieldMatrix<T,m,k> C;
- static void addProduct(A& a, const B& b, const C& c)
- {
- for (size_t row = 0; row < n; ++row)
- {
- for (size_t col = 0 ; col < k; ++col)
- {
- for (size_t i = 0; i < m; ++i)
- a[row][col] += b[row][i]*c[i][col];
- }
- }
- }
- };
-
- template<class T, int n, int m, int k>
- struct ScaledProductHelper<Dune::FieldMatrix<T,n,k>, T, Dune::FieldMatrix<T,n,m>, Dune::FieldMatrix<T,m,k>, false, false, false>
- {
- typedef Dune::FieldMatrix<T,n,k> A;
- typedef Dune::FieldMatrix<T,n,m> C;
- typedef Dune::FieldMatrix<T,m,k> D;
- static void addProduct(A& a, const T& b, const C& c, const D& d)
- {
- for (size_t row = 0; row < n; ++row) {
- for (size_t col = 0 ; col < k; ++col) {
- for (size_t i = 0; i < m; ++i)
- a[row][col] += b * c[row][i]*d[i][col];
- }
- }
- }
- };
-
- template<class T, int n>
- struct ProductHelper<Dune::DiagonalMatrix<T,n>, Dune::DiagonalMatrix<T,n>, Dune::DiagonalMatrix<T,n>, false, false, false>
- {
- typedef Dune::DiagonalMatrix<T,n> A;
- typedef A B;
- typedef B C;
- static void addProduct(A& a, const B& b, const C& c)
- {
- for (size_t i=0; i<n; ++i)
- a.diagonal(i) += b.diagonal(i)*c.diagonal(i);
- }
- };
-
- template<class T, int n>
- struct ScaledProductHelper<Dune::DiagonalMatrix<T,n>, T, Dune::DiagonalMatrix<T,n>, Dune::DiagonalMatrix<T,n>, false, false, false>
- {
- typedef Dune::DiagonalMatrix<T,n> A;
- typedef A C;
- typedef C D;
- static void addProduct(A& a, const T& b, const C& c, const D& d)
- {
- for (size_t i=0; i<n; ++i)
- a.diagonal(i) += b * c.diagonal(i)*d.diagonal(i);
- }
- };
-
- template<class T, int n>
- struct ProductHelper<Dune::FieldMatrix<T,n,n>, Dune::DiagonalMatrix<T,n>, Dune::DiagonalMatrix<T,n>, false, false, false>
- {
- typedef Dune::FieldMatrix<T,n,n> A;
- typedef Dune::DiagonalMatrix<T,n> B;
- typedef B C;
- static void addProduct(A& a, const B& b, const C& c)
- {
- for (size_t i=0; i<n; ++i)
- a[i][i] += b.diagonal(i)*c.diagonal(i);
- }
- };
-
- template<class T, int n>
- struct ScaledProductHelper<Dune::FieldMatrix<T,n,n>, T, Dune::DiagonalMatrix<T,n>, Dune::DiagonalMatrix<T,n>, false, false, false>
- {
- typedef Dune::FieldMatrix<T,n,n> A;
- typedef Dune::DiagonalMatrix<T,n> C;
- typedef C D;
- static void addProduct(A& a, const T& b, const C& c, const D& d)
- {
- for (size_t i=0; i<n; ++i)
- a[i][i] += b * c.diagonal(i)*d.diagonal(i);
- }
- };
-
- template<class T, int n, int m>
- struct ProductHelper<Dune::FieldMatrix<T,n,m>, Dune::DiagonalMatrix<T,n>, Dune::FieldMatrix<T,n,m>, false, false, false>
- {
- typedef Dune::FieldMatrix<T,n,m> A;
- typedef Dune::DiagonalMatrix<T,n> B;
- typedef A C;
- static void addProduct(A& a, const B& b, const C& c)
- {
- for (size_t i=0; i<n; ++i)
- a[i].axpy(b.diagonal(i),c[i]);
- }
- };
-
- template<class T, int n, int m>
- struct ScaledProductHelper<Dune::FieldMatrix<T,n,m>, T, Dune::DiagonalMatrix<T,n>, Dune::FieldMatrix<T,n,m>, false, false, false>
- {
- typedef Dune::FieldMatrix<T,n,m> A;
- typedef Dune::DiagonalMatrix<T,n> C;
- typedef A D;
- static void addProduct(A& a, const T& b, const C& c, const D& d)
- {
- for (size_t i=0; i<n; ++i)
- a[i].axpy(b*c.diagonal(i),d[i]);
- }
- };
-
- template<class T, int n, int m>
- struct ProductHelper<Dune::FieldMatrix<T,n,m>, Dune::FieldMatrix<T,n,m>, Dune::DiagonalMatrix<T,m>, false, false, false>
- {
- typedef Dune::FieldMatrix<T,n,m> A;
- typedef Dune::DiagonalMatrix<T,m> C;
- typedef A B;
-
- static void addProduct(A& a, const B& b, const C& c)
- {
- for (size_t i=0; i<n; ++i)
- for (size_t j=0; j<m; j++)
- a[i][j] += c.diagonal(j)*b[i][j];
- }
- };
-
- template<class T, int n, int m>
- struct ScaledProductHelper<Dune::FieldMatrix<T,n,m>, T, Dune::FieldMatrix<T,n,m>, Dune::DiagonalMatrix<T,m>, false, false, false>
- {
- typedef Dune::FieldMatrix<T,n,m> A;
- typedef Dune::DiagonalMatrix<T,m> D;
- typedef A C;
-
- static void addProduct(A& a, const T& b, const C& c, const D& d)
- {
- for (size_t i=0; i<n; ++i)
- for (size_t j=0; j<m; j++)
- a[i][j] += b*d.diagonal(j)*c[i][j];
- }
- };
-
- template<class T, int n>
- struct ProductHelper<Dune::ScaledIdentityMatrix<T,n>, Dune::ScaledIdentityMatrix<T,n>, Dune::ScaledIdentityMatrix<T,n>, false, false, false>
- {
- typedef Dune::ScaledIdentityMatrix<T,n> A;
- typedef A B;
- typedef B C;
-
- static void addProduct(A& a, const B& b, const C& c)
- {
- a.scalar() += b.scalar()*c.scalar();
- }
- };
-
- template<class T, class Scalar, int n>
- struct ScaledProductHelper<Dune::ScaledIdentityMatrix<T,n>, Scalar, Dune::ScaledIdentityMatrix<T,n>, Dune::ScaledIdentityMatrix<T,n>, false, false, false>
- {
- typedef Dune::ScaledIdentityMatrix<T,n> A;
- typedef A C;
- typedef C D;
-
- static void addProduct(A& a, const Scalar& b, const C& c, const D& d)
- {
- a.scalar() += b * c.scalar()*d.scalar();
- }
- };
-
- /** \brief Specialization for b being a scalar type
- */
- template<class A, class ScalarB, class C, bool AisScalar, bool CisScalar>
- struct ProductHelper<A, ScalarB, C, AisScalar, true, CisScalar>
- {
- typedef ScalarB B;
-
- template <class ADummy=A, enable_if_t<!MatrixTraits<ADummy>::isMatrix,int> SFINAE_Dummy=0>
- static void addProduct(A& a, const B& b, const C& c)
- {
- a.axpy(b, c);
- }
-
- template <class ADummy=A, enable_if_t<MatrixTraits<ADummy>::isMatrix,int> SFINAE_Dummy=0>
- static void addProduct(A& a, const B& b, const C& c)
- {
- typename C::ConstRowIterator ci = c.begin();
- typename C::ConstRowIterator cEnd = c.end();
- for(; ci!=cEnd; ++ci)
- {
- typename C::ConstColIterator cik =c[ci.index()].begin();
- typename C::ConstColIterator ciEnd=c[ci.index()].end();
- for(; cik!=ciEnd; ++cik)
- {
- a[ci.index()][cik.index()] += b*(*cik);
- }
- }
- }
- };
-
- template<class A, class Scalar, class ScalarB, class C, bool AisScalar, bool CisScalar>
- struct ScaledProductHelper<A, Scalar, ScalarB, C, AisScalar, true, CisScalar>
- {
- typedef ScalarB B;
-
- template <class ADummy=A, enable_if_t<!MatrixTraits<ADummy>::isMatrix,int> SFINAE_Dummy=0>
- static void addProduct(A& a, const Scalar& scalar, const B& b, const C& c)
- {
- a.axpy(scalar*b, c);
- }
-
- template <class ADummy=A, enable_if_t<MatrixTraits<ADummy>::isMatrix,int> SFINAE_Dummy=0>
- static void addProduct(A& a, const Scalar& scalar, const B& b, const C& c)
- {
- typename C::ConstRowIterator ci = c.begin();
- typename C::ConstRowIterator cEnd = c.end();
- for(; ci!=cEnd; ++ci)
- {
- typename C::ConstColIterator cik =c[ci.index()].begin();
- typename C::ConstColIterator ciEnd=c[ci.index()].end();
- for(; cik!=ciEnd; ++cik)
- {
- a[ci.index()][cik.index()] += scalar*b*(*cik);
- }
- }
- }
- };
-
- template<class ABlock, class ScalarB, class CBlock>
- struct ProductHelper<Dune::BCRSMatrix<ABlock>, ScalarB, Dune::BCRSMatrix<CBlock>, false, true, false>
- {
- typedef Dune::BCRSMatrix<ABlock> A;
- typedef ScalarB B;
- typedef Dune::BCRSMatrix<CBlock> C;
-
- static void addProduct(A& a, const B& b, const C& c)
- {
- // Workaround for the fact that Dune::BCRSMatrix assumes
- // its blocks to have an axpy() implementation, yet
- // Dune::DiagonalMatrix does not.
-
- //a.axpy(b, c);
-
- for (typename C::ConstIterator i=c.begin(); i!=c.end(); ++i) {
- const size_t iindex = i.index();
- for (typename C::row_type::ConstIterator j=i->begin(); j!=i->end(); ++j)
- ProductHelper<typename A::block_type, B, typename C::block_type, false, true, false>::addProduct(a[iindex][j.index()],b,*j);
- }
- }
- };
-
- template<class ABlock, class Scalar, class ScalarB, class CBlock>
- struct ScaledProductHelper<Dune::BCRSMatrix<ABlock>, Scalar, ScalarB, Dune::BCRSMatrix<CBlock>, false, true, false>
- {
- typedef Dune::BCRSMatrix<ABlock> A;
- typedef ScalarB B;
- typedef Dune::BCRSMatrix<CBlock> C;
-
- static void addProduct(A& a, const Scalar& scalar, const B& b, const C& c)
- {
- // Workaround for the fact that Dune::BCRSMatrix assumes
- // its blocks to have an axpy() implementation, yet
- // Dune::DiagonalMatrix does not.
-
- //a.axpy(scalar*b, c);
-
- for (typename C::ConstIterator i=c.begin(); i!=c.end(); ++i) {
- const size_t iindex = i.index();
- for (typename C::row_type::ConstIterator j=i->begin(); j!=i->end(); ++j)
- ProductHelper<typename A::block_type, B, typename C::block_type, false, true, false>::addProduct(a[iindex][j.index()],scalar*b,*j);
- }
- }
- };
-
- template<class T, int n, class ScalarB>
- struct ProductHelper<Dune::ScaledIdentityMatrix<T,n>, ScalarB, Dune::ScaledIdentityMatrix<T,n>, false, true, false>
- {
- typedef Dune::ScaledIdentityMatrix<T,n> A;
- typedef ScalarB B;
- typedef Dune::ScaledIdentityMatrix<T,n> C;
-
- static void addProduct(A& a, const B& b, const C& c)
- {
- a.scalar() += b*c.scalar();
- }
- };
-
- template<class T, int n, class Scalar, class ScalarB>
- struct ScaledProductHelper<Dune::ScaledIdentityMatrix<T,n>, Scalar, ScalarB, Dune::ScaledIdentityMatrix<T,n>, false, true, false>
- {
- typedef Dune::ScaledIdentityMatrix<T,n> A;
- typedef ScalarB B;
- typedef Dune::ScaledIdentityMatrix<T,n> C;
-
- static void addProduct(A& a, const Scalar& scalar, const B& b, const C& c)
- {
- a.scalar() += scalar*b*c.scalar();
- }
- };
-
- template<class A, class B, class C>
- struct ProductHelper<A, B, C, true, true, true>
- {
- static void addProduct(A& a, const B& b, const C& c)
- {
- a += b*c;
- }
- };
-
- template<class A, class B, class C, class D>
- struct ScaledProductHelper<A, B, C, D, true, true, true>
- {
- static void addProduct(A& a, const B& b, const C& c, const D& d)
- {
- a += b*c*d;
- }
- };
-
//! Helper class for computing the cross product
template <class T, int n>
struct CrossProductHelper
@@ -651,7 +273,7 @@ namespace Arithmetic
template<class A, class B, class C>
void addProduct(A& a, const B& b, const C& c)
{
- ProductHelper<A,B,C,ScalarTraits<A>::isScalar, ScalarTraits<B>::isScalar, ScalarTraits<C>::isScalar>::addProduct(a,b,c);
+ Dune::MatrixVector::addProduct(a,b,c);
}
/** \brief Subtract a product from some matrix or vector
@@ -668,7 +290,7 @@ namespace Arithmetic
template<class A, class B, class C>
void subtractProduct(A& a, const B& b, const C& c)
{
- ScaledProductHelper<A,int,B,C,ScalarTraits<A>::isScalar, ScalarTraits<B>::isScalar, ScalarTraits<C>::isScalar>::addProduct(a,-1,b,c);
+ Dune::MatrixVector::subtractProduct(a,b,c);
}
//! Compute the cross product of two vectors. Only works for n==3
@@ -700,7 +322,7 @@ namespace Arithmetic
typename std::enable_if<ScalarTraits<B>::isScalar, void>::type
addProduct(A& a, const B& b, const C& c, const D& d)
{
- ScaledProductHelper<A,B,C,D,ScalarTraits<A>::isScalar, ScalarTraits<C>::isScalar, ScalarTraits<D>::isScalar>::addProduct(a,b,c,d);
+ Dune::MatrixVector::addProduct(a,b,c,d);
}
/** \brief Subtract a scaled product from some matrix or vector
@@ -718,7 +340,7 @@ namespace Arithmetic
typename std::enable_if<ScalarTraits<B>::isScalar, void>::type
subtractProduct(A& a, const B& b, const C& c, const D& d)
{
- ScaledProductHelper<A,B,C,D,ScalarTraits<A>::isScalar, ScalarTraits<C>::isScalar, ScalarTraits<D>::isScalar>::addProduct(a,-b,c,d);
+ Dune::MatrixVector::subtractProduct(a,b,c,d);
}
/** \brief Internal helper class for Matrix operations
diff --git a/dune/fufem/assemblers/boundaryoperatorassembler.hh b/dune/fufem/assemblers/boundaryoperatorassembler.hh
index ab9766f732819280811ad6c46bbcc24943bb1b97..25a9e11bd54ee44d0ecdbd1a0bfd6d1e66e59289 100644
--- a/dune/fufem/assemblers/boundaryoperatorassembler.hh
+++ b/dune/fufem/assemblers/boundaryoperatorassembler.hh
@@ -4,8 +4,9 @@
#include <dune/istl/matrix.hh>
#include <dune/istl/matrixindexset.hh>
+#include <dune/matrix-vector/axpy.hh>
+
#include "dune/fufem/functionspacebases/functionspacebasis.hh"
-#include "dune/fufem/arithmetic.hh"
#include "dune/fufem/boundarypatch.hh"
//! Generic global assembler for operators on a gridview
@@ -174,11 +175,11 @@ class BoundaryOperatorAssembler
{
for (size_t rw=0; rw<rowConstraints.size(); ++rw)
A[rowConstraints[rw].index][rowConstraints[rw].index].axpy(rowConstraints[rw].factor, localA[i][j]);
-// Arithmetic::addProduct(A[rowConstraints[rw].index][rowConstraints[rw].index], rowConstraints[rw].factor, localA[i][j]);
+// Dune::MatrixVector::addProduct(A[rowConstraints[rw].index][rowConstraints[rw].index], rowConstraints[rw].factor, localA[i][j]);
}
else
A[rowIndex][rowIndex].axpy(1.0, localA[i][j]);
-// Arithmetic::addProduct(A[rowIndex][rowIndex], 1.0, localA[i][j]);
+// Dune::MatrixVector::addProduct(A[rowIndex][rowIndex], 1.0, localA[i][j]);
}
else
{
@@ -191,24 +192,24 @@ class BoundaryOperatorAssembler
{
for (size_t cw=0; cw<colConstraints.size(); ++cw)
A[rowConstraints[rw].index][colConstraints[cw].index].axpy(rowConstraints[rw].factor * colConstraints[cw].factor, localA[i][j]);
-// Arithmetic::addProduct(A[rowConstraints[rw].index][colConstraints[cw].index], rowConstraints[rw].factor * colConstraints[cw].factor, localA[i][j]);
+// Dune::MatrixVector::addProduct(A[rowConstraints[rw].index][colConstraints[cw].index], rowConstraints[rw].factor * colConstraints[cw].factor, localA[i][j]);
}
}
else if (rowIsConstrained)
{
for (size_t rw=0; rw<rowConstraints.size(); ++rw)
A[rowConstraints[rw].index][colIndex].axpy(rowConstraints[rw].factor, localA[i][j]);
-// Arithmetic::addProduct(A[rowConstraints[rw].index][colIndex], rowConstraints[rw].factor, localA[i][j]);
+// Dune::MatrixVector::addProduct(A[rowConstraints[rw].index][colIndex], rowConstraints[rw].factor, localA[i][j]);
}
else if (colIsConstrained)
{
for (size_t cw=0; cw<colConstraints.size(); ++cw)
A[rowIndex][colConstraints[cw].index].axpy(colConstraints[cw].factor, localA[i][j]);
-// Arithmetic::addProduct(A[rowIndex][colConstraints[cw].index], colConstraints[cw].factor, localA[i][j]);
+// Dune::MatrixVector::addProduct(A[rowIndex][colConstraints[cw].index], colConstraints[cw].factor, localA[i][j]);
}
else
A[rowIndex][colIndex].axpy(1.0, localA[i][j]);
-// Arithmetic::addProduct(A[rowIndex][colIndex], 1.0, localA[i][j]);
+// Dune::MatrixVector::addProduct(A[rowIndex][colIndex], 1.0, localA[i][j]);
}
}
}
diff --git a/dune/fufem/assemblers/integraloperatorassembler.hh b/dune/fufem/assemblers/integraloperatorassembler.hh
index 0198ff55f678f04921d9fabb636fbc45f98186c1..d85cc464e8aacd8201f757a2b8d9b29a9d760ce9 100644
--- a/dune/fufem/assemblers/integraloperatorassembler.hh
+++ b/dune/fufem/assemblers/integraloperatorassembler.hh
@@ -4,7 +4,8 @@
#include <dune/istl/matrix.hh>
#include <dune/istl/matrixindexset.hh>
-#include "dune/fufem/arithmetic.hh"
+#include <dune/matrix-vector/axpy.hh>
+
#include "dune/fufem/functionspacebases/functionspacebasis.hh"
//! Generic global assembler for operators on a gridview
@@ -89,21 +90,21 @@ class IntegralOperatorAssembler
for (size_t rw=0; rw<rowConstraints.size(); ++rw)
{
for (size_t cw=0; cw<colConstraints.size(); ++cw)
- Arithmetic::addProduct(A[rowConstraints[rw].index][colConstraints[cw].index], rowConstraints[rw].factor * colConstraints[cw].factor, localA[i][j]);
+ Dune::MatrixVector::addProduct(A[rowConstraints[rw].index][colConstraints[cw].index], rowConstraints[rw].factor * colConstraints[cw].factor, localA[i][j]);
}
}
else if (rowIsConstrained)
{
for (size_t rw=0; rw<rowConstraints.size(); ++rw)
- Arithmetic::addProduct(A[rowConstraints[rw].index][colIndex], rowConstraints[rw].factor, localA[i][j]);
+ Dune::MatrixVector::addProduct(A[rowConstraints[rw].index][colIndex], rowConstraints[rw].factor, localA[i][j]);
}
else if (colIsConstrained)
{
for (size_t cw=0; cw<colConstraints.size(); ++cw)
- Arithmetic::addProduct(A[rowIndex][colConstraints[cw].index], colConstraints[cw].factor, localA[i][j]);
+ Dune::MatrixVector::addProduct(A[rowIndex][colConstraints[cw].index], colConstraints[cw].factor, localA[i][j]);
}
else
- Arithmetic::addProduct(A[rowIndex][colIndex], 1.0, localA[i][j]);
+ Dune::MatrixVector::addProduct(A[rowIndex][colIndex], 1.0, localA[i][j]);
}
}
}
diff --git a/dune/fufem/assemblers/intersectionoperatorassembler.hh b/dune/fufem/assemblers/intersectionoperatorassembler.hh
index 469b7429dc40f2b4d36b80acac8fcb5213dd4002..60eb394ef08726d97dc26676348a916f1b01e468 100644
--- a/dune/fufem/assemblers/intersectionoperatorassembler.hh
+++ b/dune/fufem/assemblers/intersectionoperatorassembler.hh
@@ -8,7 +8,7 @@
#include <dune/grid/common/mcmgmapper.hh>
-#include "dune/fufem/arithmetic.hh"
+#include <dune/matrix-vector/axpy.hh>
// TODO Modify such that the global assembler also works for nonconforming grids!
@@ -279,10 +279,10 @@ class IntersectionOperatorAssembler
if (rowIsConstrained)
{
for (size_t rw=0; rw<rowConstraints.size(); ++rw)
- Arithmetic::addProduct(A[rowConstraints[rw].index][rowConstraints[rw].index], rowConstraints[rw].factor, localA[i][j]);
+ Dune::MatrixVector::addProduct(A[rowConstraints[rw].index][rowConstraints[rw].index], rowConstraints[rw].factor, localA[i][j]);
}
else
- Arithmetic::addProduct(A[rowIndex][rowIndex], 1.0, localA[i][j]);
+ Dune::MatrixVector::addProduct(A[rowIndex][rowIndex], 1.0, localA[i][j]);
}
else
{
@@ -294,21 +294,21 @@ class IntersectionOperatorAssembler
for (size_t rw=0; rw<rowConstraints.size(); ++rw)
{
for (size_t cw=0; cw<colConstraints.size(); ++cw)
- Arithmetic::addProduct(A[rowConstraints[rw].index][colConstraints[cw].index], rowConstraints[rw].factor * colConstraints[cw].factor, localA[i][j]);
+ Dune::MatrixVector::addProduct(A[rowConstraints[rw].index][colConstraints[cw].index], rowConstraints[rw].factor * colConstraints[cw].factor, localA[i][j]);
}
}
else if (rowIsConstrained)
{
for (size_t rw=0; rw<rowConstraints.size(); ++rw)
- Arithmetic::addProduct(A[rowConstraints[rw].index][colIndex], rowConstraints[rw].factor, localA[i][j]);
+ Dune::MatrixVector::addProduct(A[rowConstraints[rw].index][colIndex], rowConstraints[rw].factor, localA[i][j]);
}
else if (colIsConstrained)
{
for (size_t cw=0; cw<colConstraints.size(); ++cw)
- Arithmetic::addProduct(A[rowIndex][colConstraints[cw].index], colConstraints[cw].factor, localA[i][j]);
+ Dune::MatrixVector::addProduct(A[rowIndex][colConstraints[cw].index], colConstraints[cw].factor, localA[i][j]);
}
else
- Arithmetic::addProduct(A[rowIndex][colIndex], 1.0, localA[i][j]);
+ Dune::MatrixVector::addProduct(A[rowIndex][colIndex], 1.0, localA[i][j]);
}
}
}
@@ -365,10 +365,10 @@ class IntersectionOperatorAssembler
if (rowIsConstrained)
{
for (size_t rw=0; rw<rowConstraints.size(); ++rw)
- Arithmetic::addProduct(A[rowConstraints[rw].index][rowConstraints[rw].index], rowConstraints[rw].factor, localA[i][j]);
+ Dune::MatrixVector::addProduct(A[rowConstraints[rw].index][rowConstraints[rw].index], rowConstraints[rw].factor, localA[i][j]);
}
else
- Arithmetic::addProduct(A[rowIndex][rowIndex], 1.0, localA[i][j]);
+ Dune::MatrixVector::addProduct(A[rowIndex][rowIndex], 1.0, localA[i][j]);
}
else
{
@@ -380,21 +380,21 @@ class IntersectionOperatorAssembler
for (size_t rw=0; rw<rowConstraints.size(); ++rw)
{
for (size_t cw=0; cw<colConstraints.size(); ++cw)
- Arithmetic::addProduct(A[rowConstraints[rw].index][colConstraints[cw].index], rowConstraints[rw].factor * colConstraints[cw].factor, localA[i][j]);
+ Dune::MatrixVector::addProduct(A[rowConstraints[rw].index][colConstraints[cw].index], rowConstraints[rw].factor * colConstraints[cw].factor, localA[i][j]);
}
}
else if (rowIsConstrained)
{
for (size_t rw=0; rw<rowConstraints.size(); ++rw)
- Arithmetic::addProduct(A[rowConstraints[rw].index][colIndex], rowConstraints[rw].factor, localA[i][j]);
+ Dune::MatrixVector::addProduct(A[rowConstraints[rw].index][colIndex], rowConstraints[rw].factor, localA[i][j]);
}
else if (colIsConstrained)
{
for (size_t cw=0; cw<colConstraints.size(); ++cw)
- Arithmetic::addProduct(A[rowIndex][colConstraints[cw].index], colConstraints[cw].factor, localA[i][j]);
+ Dune::MatrixVector::addProduct(A[rowIndex][colConstraints[cw].index], colConstraints[cw].factor, localA[i][j]);
}
else
- Arithmetic::addProduct(A[rowIndex][colIndex], 1.0, localA[i][j]);
+ Dune::MatrixVector::addProduct(A[rowIndex][colIndex], 1.0, localA[i][j]);
}
}
}
diff --git a/dune/fufem/assemblers/localassemblers/generalizedlaplaceassembler.hh b/dune/fufem/assemblers/localassemblers/generalizedlaplaceassembler.hh
index 27dc1a305fad8cd2e21ae36590998e5de7c3470f..8aabb2f235fe3121efa655e976ce28077c78d72e 100644
--- a/dune/fufem/assemblers/localassemblers/generalizedlaplaceassembler.hh
+++ b/dune/fufem/assemblers/localassemblers/generalizedlaplaceassembler.hh
@@ -9,7 +9,9 @@
#include <dune/geometry/quadraturerules.hh>
-#include <dune/fufem/arithmetic.hh>
+#include <dune/matrix-vector/axpy.hh>
+#include <dune/matrix-vector/matrixtraits.hh>
+#include <dune/matrix-vector/scalartraits.hh>
#include <dune/fufem/quadraturerules/quadraturerulecache.hh>
@@ -35,15 +37,15 @@ class GeneralizedLaplaceAssembler : public LocalOperatorAssembler < GridType, Tr
static void checkCoefficientType()
{
typedef typename FunctionType::RangeType RangeType;
- if (not(Arithmetic::ScalarTraits<RangeType>::isScalar))
+ if (not(Dune::MatrixVector::ScalarTraits<RangeType>::isScalar))
{
// check if non-scalar coefficient type is a matrix type
- if (not(Arithmetic::MatrixTraits<RangeType>::isMatrix))
+ if (not(Dune::MatrixVector::MatrixTraits<RangeType>::isMatrix))
DUNE_THROW(Dune::Exception,"Coefficient function for generalized Laplace s neither scalar nor matrix valued.");
// check if matrix coefficient type has proper dimensions
- int crows = Arithmetic::MatrixTraits<RangeType>::rows;
- int ccols = Arithmetic::MatrixTraits<RangeType>::cols;
+ int crows = Dune::MatrixVector::MatrixTraits<RangeType>::rows;
+ int ccols = Dune::MatrixVector::MatrixTraits<RangeType>::cols;
if (not(crows==GridType::dimension and ccols==GridType::dimension))
DUNE_THROW(Dune::Exception,"Coefficient function for generalized Laplace has wrong dimension of range.\n Should be 1x1 or " << GridType::dimension << "x" << GridType::dimension << ", is " << crows << "x" << ccols << ".");
}
@@ -149,7 +151,7 @@ class GeneralizedLaplaceAssembler : public LocalOperatorAssembler < GridType, Tr
{
GlobalCoordinate dummy(0.0);
- Arithmetic::addProduct(dummy, coeffs, gradients[i]);
+ Dune::MatrixVector::addProduct(dummy, coeffs, gradients[i]);
for (int j=i+1; j<cols; ++j)
{
diff --git a/dune/fufem/assemblers/localassemblers/scaledsumfunctionalassembler.hh b/dune/fufem/assemblers/localassemblers/scaledsumfunctionalassembler.hh
index 6fa0859ca76bbd36051dc9fd972b5b4458feb83f..14f1cca78640cdd196f4a229c2fd8d8866243e68 100644
--- a/dune/fufem/assemblers/localassemblers/scaledsumfunctionalassembler.hh
+++ b/dune/fufem/assemblers/localassemblers/scaledsumfunctionalassembler.hh
@@ -1,8 +1,9 @@
#ifndef DUNE_FUFEM_ASSEMBLERS_LOCALASSEMBLERS_SCALEDSUMFUNCTIONALASSEMBLER_HH
#define DUNE_FUFEM_ASSEMBLERS_LOCALASSEMBLERS_SCALEDSUMFUNCTIONALASSEMBLER_HH
+#include <dune/matrix-vector/axpy.hh>
+
#include <dune/fufem/assemblers/localfunctionalassembler.hh>
-#include <dune/fufem/arithmetic.hh>
template<typename Grid, typename FE, typename T>
class LocalSumFunctionalAssembler
@@ -27,7 +28,7 @@ public:
LocalVector v(localVector.size());
scaledAssembler.assembler().assemble(element, v = 0.0, tFE);
for(uint i=0; i<v.size(); ++i)
- Arithmetic::addProduct(localVector[i], scaledAssembler.scale(), v[i]);
+ Dune::MatrixVector::addProduct(localVector[i], scaledAssembler.scale(), v[i]);
}
}
diff --git a/dune/fufem/assemblers/localassemblers/secondorderoperatorassembler.hh b/dune/fufem/assemblers/localassemblers/secondorderoperatorassembler.hh
index 4aa112bc1e13d391e4155addbd94d416d89f60b1..2629dd718d5092fde06a7124f9fd235eea3b3794 100644
--- a/dune/fufem/assemblers/localassemblers/secondorderoperatorassembler.hh
+++ b/dune/fufem/assemblers/localassemblers/secondorderoperatorassembler.hh
@@ -16,6 +16,8 @@
#include <dune/functions/common/utility.hh>
#include <dune/functions/gridfunctions/gridfunction.hh>
+#include <dune/matrix-vector/axpy.hh>
+
#include "dune/fufem/staticmatrixtools.hh"
#include "dune/fufem/arithmetic.hh"
#include "dune/fufem/quadraturerules/quadraturerulecache.hh"
@@ -54,7 +56,7 @@ namespace Concept {
template<class F>
auto require(F&& f) -> decltype(
- Arithmetic::addProduct(r, s, Dune::Std::apply(f, args))
+ Dune::MatrixVector::addProduct(r, s, Dune::Std::apply(f, args))
);
};
@@ -164,7 +166,7 @@ public:
typename std::result_of<FunctionTypes(WorldCoordinate)>::type...>;
static_assert(
Dune::Fufem::Concept::models<ContractionConcept, Contraction>(),
- "Contraction type does not support: Arithmetic::addProduct(MatrixEntry, double, contraction(WorldCoordinate, WorldCoordinate, ...))");
+ "Contraction type does not support: Dune::MatrixVector::addProduct(MatrixEntry, double, contraction(WorldCoordinate, WorldCoordinate, ...))");
}
SecondOrderOperatorAssembler(const Contraction& contraction,
bool isSymmetric,
@@ -253,7 +255,7 @@ public:
const auto boundContraction = [&](auto... args) {
return contraction_(gradients_[i], gradients_[j], args...);
};
- Arithmetic::addProduct(
+ Dune::MatrixVector::addProduct(
localMatrix[i][j], z,
Dune::Std::apply(boundContraction, returnValues));
}
diff --git a/dune/fufem/assemblers/localassemblers/subgridgeneralizedlaplaceassembler.hh b/dune/fufem/assemblers/localassemblers/subgridgeneralizedlaplaceassembler.hh
index 1a9dea4d16e3ae43a93c95ac8b49d7370940c690..7809b4d04af2692201d3a7088872eb8e72d44884 100644
--- a/dune/fufem/assemblers/localassemblers/subgridgeneralizedlaplaceassembler.hh
+++ b/dune/fufem/assemblers/localassemblers/subgridgeneralizedlaplaceassembler.hh
@@ -9,7 +9,9 @@
#include <dune/geometry/quadraturerules.hh>
-#include <dune/fufem/arithmetic.hh>
+#include <dune/matrix-vector/arithmetic.hh>
+#include <dune/matrix-vector/matrixtraits.hh>
+#include <dune/matrix-vector/scalartraits.hh>
#include <dune/fufem/quadraturerules/quadraturerulecache.hh>
@@ -52,15 +54,15 @@ class SubgridGeneralizedLaplaceAssembler : public LocalOperatorAssembler < GridT
static void checkCoefficientType()
{
typedef typename FunctionType::RangeType RangeType;
- if (not(Arithmetic::ScalarTraits<RangeType>::isScalar))
+ if (not(Dune::MatrixVector::ScalarTraits<RangeType>::isScalar))
{
// check if non-scalar coefficient type is a matrix type
- if (not(Arithmetic::MatrixTraits<RangeType>::isMatrix))
+ if (not(Dune::MatrixVector::MatrixTraits<RangeType>::isMatrix))
DUNE_THROW(Dune::Exception,"Coefficient function for generalized Laplace is neither scalar nor matrix valued.");
// check if matrix coefficient type has proper dimensions
- int crows = Arithmetic::MatrixTraits<RangeType>::rows;
- int ccols = Arithmetic::MatrixTraits<RangeType>::cols;
+ int crows = Dune::MatrixVector::MatrixTraits<RangeType>::rows;
+ int ccols = Dune::MatrixVector::MatrixTraits<RangeType>::cols;
if (not(crows==GridType::dimension and ccols==GridType::dimension))
DUNE_THROW(Dune::Exception,"Coefficient function for generalized Laplace has wrong dimension of range.\n Should be 1x1 or " << GridType::dimension << "x" << GridType::dimension << ", is " << crows << "x" << ccols << ".");
}
@@ -168,7 +170,7 @@ class SubgridGeneralizedLaplaceAssembler : public LocalOperatorAssembler < GridT
{
GlobalCoordinate dummy(0.0);
- Arithmetic::addProduct(dummy, coeffs, gradients[i]);
+ Dune::MatrixVector::addProduct(dummy, coeffs, gradients[i]);
for (int j=i+1; j<cols; ++j)
{
@@ -227,7 +229,7 @@ class SubgridGeneralizedLaplaceAssembler : public LocalOperatorAssembler < GridT
{
GlobalCoordinate dummy(0.0);
- Arithmetic::addProduct(dummy, coeffs, gradients[i]);
+ Dune::MatrixVector::addProduct(dummy, coeffs, gradients[i]);
for (int j=i+1; j<cols; ++j)
{
diff --git a/dune/fufem/assemblers/localassemblers/vvlaplaceassembler.hh b/dune/fufem/assemblers/localassemblers/vvlaplaceassembler.hh
index 99e283dcd2deb2e656cc7c460545ac290928f5db..0fc310ccb627ec1276f777d48ff4469810f6c331 100644
--- a/dune/fufem/assemblers/localassemblers/vvlaplaceassembler.hh
+++ b/dune/fufem/assemblers/localassemblers/vvlaplaceassembler.hh
@@ -6,7 +6,8 @@
#include <dune/istl/matrix.hh>
-#include <dune/fufem/arithmetic.hh>
+#include <dune/matrix-vector/axpy.hh>
+
#include <dune/fufem/quadraturerules/quadraturerulecache.hh>
#include <dune/fufem/assemblers/localoperatorassembler.hh>
@@ -105,10 +106,10 @@ class VVLaplaceAssembler : public LocalOperatorAssembler <GridType, TrialLocalFE
for (int j=i+1; j<cols; ++j)
{
double zij = (gradients[i] * gradients[j]) * z;
- Arithmetic::addProduct(localMatrix[i][j], zij, one_);
- Arithmetic::addProduct(localMatrix[j][i], zij, one_);
+ Dune::MatrixVector::addProduct(localMatrix[i][j], zij, one_);
+ Dune::MatrixVector::addProduct(localMatrix[j][i], zij, one_);
}
- Arithmetic::addProduct(localMatrix[i][i], gradients[i] * gradients[i] * z, one_);
+ Dune::MatrixVector::addProduct(localMatrix[i][i], gradients[i] * gradients[i] * z, one_);
}
}
return;
@@ -193,10 +194,10 @@ class VVLaplaceAssembler : public LocalOperatorAssembler <GridType, TrialLocalFE
for (int j=i+1; j<cols; ++j)
{
double zij = (gradients[i] * gradients[j]) * z;
- Arithmetic::addProduct(localMatrix[i][j], zij, one_);
- Arithmetic::addProduct(localMatrix[j][i], zij, one_);
+ Dune::MatrixVector::addProduct(localMatrix[i][j], zij, one_);
+ Dune::MatrixVector::addProduct(localMatrix[j][i], zij, one_);
}
- Arithmetic::addProduct(localMatrix[i][i], gradients[i] * gradients[i] * z, one_);
+ Dune::MatrixVector::addProduct(localMatrix[i][i], gradients[i] * gradients[i] * z, one_);
}
}
diff --git a/dune/fufem/assemblers/localassemblers/vvmassassembler.hh b/dune/fufem/assemblers/localassemblers/vvmassassembler.hh
index 0f1d57f5da8c7898968b433b6a2ef303cdf4a3dd..2f8e14eac80baa693f9194159168705308929d64 100644
--- a/dune/fufem/assemblers/localassemblers/vvmassassembler.hh
+++ b/dune/fufem/assemblers/localassemblers/vvmassassembler.hh
@@ -6,7 +6,8 @@
#include <dune/istl/matrix.hh>
-#include <dune/fufem/arithmetic.hh>
+#include <dune/matrix-vector/axpy.hh>
+
#include <dune/fufem/quadraturerules/quadraturerulecache.hh>
#include <dune/fufem/assemblers/localoperatorassembler.hh>
@@ -91,10 +92,10 @@ class VVMassAssembler : public LocalOperatorAssembler < GridType, TrialLocalFE,
for (int j=i+1; j<cols; ++j)
{
double zij = values[j] * zi;
- Arithmetic::addProduct(localMatrix[i][j], zij, one_);
- Arithmetic::addProduct(localMatrix[j][i], zij, one_);
+ Dune::MatrixVector::addProduct(localMatrix[i][j], zij, one_);
+ Dune::MatrixVector::addProduct(localMatrix[j][i], zij, one_);
}
- Arithmetic::addProduct(localMatrix[i][i], values[i] * zi, one_);
+ Dune::MatrixVector::addProduct(localMatrix[i][i], values[i] * zi, one_);
}
}
}
diff --git a/dune/fufem/assemblers/operatorassembler.hh b/dune/fufem/assemblers/operatorassembler.hh
index 6e31ce8dc8c8860e5284fc6732847935286e3adc..3a4cdaafa5bc28df6b0a6f007c67dbb761605abd 100644
--- a/dune/fufem/assemblers/operatorassembler.hh
+++ b/dune/fufem/assemblers/operatorassembler.hh
@@ -4,8 +4,9 @@
#include <dune/istl/matrix.hh>
#include <dune/istl/matrixindexset.hh>
+#include <dune/matrix-vector/axpy.hh>
+
#include "dune/fufem/functionspacebases/functionspacebasis.hh"
-#include "dune/fufem/arithmetic.hh"
//! Generic global assembler for operators on a gridview
template <class TrialBasis, class AnsatzBasis>
@@ -166,10 +167,10 @@ class OperatorAssembler
if (rowIsConstrained)
{
for (size_t rw=0; rw<rowConstraints.size(); ++rw)
- Arithmetic::addProduct(A[rowConstraints[rw].index][rowConstraints[rw].index], rowConstraints[rw].factor, localA[i][j]);
+ Dune::MatrixVector::addProduct(A[rowConstraints[rw].index][rowConstraints[rw].index], rowConstraints[rw].factor, localA[i][j]);
}
else
- Arithmetic::addProduct(A[rowIndex][rowIndex], 1.0, localA[i][j]);
+ Dune::MatrixVector::addProduct(A[rowIndex][rowIndex], 1.0, localA[i][j]);
}
else
{
@@ -181,21 +182,21 @@ class OperatorAssembler
for (size_t rw=0; rw<rowConstraints.size(); ++rw)
{
for (size_t cw=0; cw<colConstraints.size(); ++cw)
- Arithmetic::addProduct(A[rowConstraints[rw].index][colConstraints[cw].index], rowConstraints[rw].factor * colConstraints[cw].factor, localA[i][j]);
+ Dune::MatrixVector::addProduct(A[rowConstraints[rw].index][colConstraints[cw].index], rowConstraints[rw].factor * colConstraints[cw].factor, localA[i][j]);
}
}
else if (rowIsConstrained)
{
for (size_t rw=0; rw<rowConstraints.size(); ++rw)
- Arithmetic::addProduct(A[rowConstraints[rw].index][colIndex], rowConstraints[rw].factor, localA[i][j]);
+ Dune::MatrixVector::addProduct(A[rowConstraints[rw].index][colIndex], rowConstraints[rw].factor, localA[i][j]);
}
else if (colIsConstrained)
{
for (size_t cw=0; cw<colConstraints.size(); ++cw)
- Arithmetic::addProduct(A[rowIndex][colConstraints[cw].index], colConstraints[cw].factor, localA[i][j]);
+ Dune::MatrixVector::addProduct(A[rowIndex][colConstraints[cw].index], colConstraints[cw].factor, localA[i][j]);
}
else
- Arithmetic::addProduct(A[rowIndex][colIndex], 1.0, localA[i][j]);
+ Dune::MatrixVector::addProduct(A[rowIndex][colIndex], 1.0, localA[i][j]);
}
}
}
diff --git a/dune/fufem/assemblers/preconditioneddefectassembler.hh b/dune/fufem/assemblers/preconditioneddefectassembler.hh
index 2bc6f8b0ea718b0a88a1ef45a22814019678b0f7..69c53853e1ffa018f641a5361cb438d1fe0bf0c0 100644
--- a/dune/fufem/assemblers/preconditioneddefectassembler.hh
+++ b/dune/fufem/assemblers/preconditioneddefectassembler.hh
@@ -5,8 +5,9 @@
#include <dune/istl/matrix.hh>
+#include <dune/matrix-vector/axpy.hh>
+
#include "dune/fufem/functionspacebases/functionspacebasis.hh"
-#include <dune/fufem/arithmetic.hh>
//! \todo Please doc me!
template <class ExtensionBasisType, class ExtendedBasisType>
@@ -124,7 +125,7 @@ class PreconditionedDefectAssembler
for (size_t k=0; k<alpha.size(); ++k)
localA[i][j].usmv(-rowConstraints[rw].factor*alpha[k], (*x[k])[rowConstraints[rw].index], (*b[k])[rowConstraints[rw].index]);
// A[rowConstraints[rw].index][rowConstraints[rw].index].axpy(rowConstraints[rw].factor, localA[i][j]);
- Arithmetic::addProduct(A[rowConstraints[rw].index][rowConstraints[rw].index], rowConstraints[rw].factor, localA[i][j]);
+ Dune::MatrixVector::addProduct(A[rowConstraints[rw].index][rowConstraints[rw].index], rowConstraints[rw].factor, localA[i][j]);
}
}
else
@@ -132,7 +133,7 @@ class PreconditionedDefectAssembler
for (size_t k=0; k<alpha.size(); ++k)
localA[i][j].usmv(-alpha[k], (*x[k])[rowIndex], (*b[k])[rowIndex]);
// A[rowIndex][rowIndex].axpy(1.0, localA[i][j]);
- Arithmetic::addProduct(A[rowIndex][rowIndex], 1.0, localA[i][j]);
+ Dune::MatrixVector::addProduct(A[rowIndex][rowIndex], 1.0, localA[i][j]);
}
}
else
@@ -153,7 +154,7 @@ class PreconditionedDefectAssembler
(*b[k])[rowConstraints[rw].index]);
if (rowConstraints[rw].index == colConstraints[cw].index)
// A[rowConstraints[rw].index][colConstraints[cw].index].axpy(rowConstraints[rw].factor * colConstraints[cw].factor, localA[i][j]);
- Arithmetic::addProduct(A[rowConstraints[rw].index][colConstraints[cw].index], rowConstraints[rw].factor * colConstraints[cw].factor, localA[i][j]);
+ Dune::MatrixVector::addProduct(A[rowConstraints[rw].index][colConstraints[cw].index], rowConstraints[rw].factor * colConstraints[cw].factor, localA[i][j]);
}
}
}
@@ -165,7 +166,7 @@ class PreconditionedDefectAssembler
localA[i][j].usmv(-rowConstraints[rw].factor*alpha[k], (*x[k])[colIndex], (*b[k])[rowConstraints[rw].index]);
if (rowConstraints[rw].index == colIndex)
// A[rowConstraints[rw].index][colIndex].axpy(rowConstraints[rw].factor, localA[i][j]);
- Arithmetic::addProduct(A[rowConstraints[rw].index][colIndex], rowConstraints[rw].factor, localA[i][j]);
+ Dune::MatrixVector::addProduct(A[rowConstraints[rw].index][colIndex], rowConstraints[rw].factor, localA[i][j]);
}
}
else if (colIsConstrained)
@@ -176,7 +177,7 @@ class PreconditionedDefectAssembler
localA[i][j].usmv(-colConstraints[cw].factor*alpha[k], (*x[k])[colConstraints[cw].index], (*b[k])[rowIndex]);
if (rowIndex == colConstraints[cw].index)
// A[rowIndex][colConstraints[cw].index].axpy(colConstraints[cw].factor, localA[i][j]);
- Arithmetic::addProduct(A[rowIndex][colConstraints[cw].index], colConstraints[cw].factor, localA[i][j]);
+ Dune::MatrixVector::addProduct(A[rowIndex][colConstraints[cw].index], colConstraints[cw].factor, localA[i][j]);
}
}
else
@@ -185,7 +186,7 @@ class PreconditionedDefectAssembler
localA[i][j].usmv(-alpha[k], (*x[k])[colIndex], (*b[k])[rowIndex]);
if (rowIndex == colIndex)
// A[rowIndex][colIndex].axpy(1.0, localA[i][j]);
- Arithmetic::addProduct(A[rowIndex][colIndex], 1.0, localA[i][j]);
+ Dune::MatrixVector::addProduct(A[rowIndex][colIndex], 1.0, localA[i][j]);
}
}
diff --git a/dune/fufem/functions/basisgridfunction.hh b/dune/fufem/functions/basisgridfunction.hh
index 4caa556124de9ba7e0abc9e12c9fcd04dfbbc210..122646782aabcd7d39a09cabc804b1a8d45c806a 100644
--- a/dune/fufem/functions/basisgridfunction.hh
+++ b/dune/fufem/functions/basisgridfunction.hh
@@ -11,7 +11,8 @@
#include <dune/istl/bvector.hh>
-#include <dune/fufem/arithmetic.hh>
+#include <dune/matrix-vector/axpy.hh>
+
#include <dune/fufem/quadraturerules/quadraturerulecache.hh>
#include <dune/fufem/functions/virtualgridfunction.hh>
@@ -51,10 +52,10 @@ struct LocalEvaluator
{
const LinearCombination& constraints = basis.constraints(index);
for (size_t w=0; w<constraints.size(); ++w)
- Arithmetic::addProduct(y, values[i] * constraints[w].factor, coefficients[constraints[w].index]);
+ Dune::MatrixVector::addProduct(y, values[i] * constraints[w].factor, coefficients[constraints[w].index]);
}
else
- Arithmetic::addProduct(y, values[i], coefficients[index]);
+ Dune::MatrixVector::addProduct(y, values[i], coefficients[index]);
}
}
};
diff --git a/dune/fufem/functiontools/basisidmapper.hh b/dune/fufem/functiontools/basisidmapper.hh
index 590229b9d49cef007175fef7927cc3a6eb8b5fef..3007d9b5455a2c9011e790d3199a701f2a1d7885 100644
--- a/dune/fufem/functiontools/basisidmapper.hh
+++ b/dune/fufem/functiontools/basisidmapper.hh
@@ -10,7 +10,8 @@
#include <dune/localfunctions/common/localkey.hh>
-#include <dune/fufem/staticmatrixtools.hh>
+#include <dune/matrix-vector/axpy.hh>
+
#include <dune/fufem/functions/virtualgridfunction.hh>
#include <dune/fufem/functionspacebases/functionspacebasis.hh>
@@ -283,10 +284,10 @@ class BasisIdGridFunction :
{
const LinearCombination& constraints = mapper_.constraints(index);
for (size_t w=0; w<constraints.size(); ++w)
- Arithmetic::addProduct(y, values[i] * constraints[w].factor, coefficients_[constraints[w].index]);
+ Dune::MatrixVector::addProduct(y, values[i] * constraints[w].factor, coefficients_[constraints[w].index]);
}
else
- Arithmetic::addProduct(y, values[i], coefficients_[index]);
+ Dune::MatrixVector::addProduct(y, values[i], coefficients_[index]);
}
}
diff --git a/dune/fufem/functiontools/basisinterpolator.hh b/dune/fufem/functiontools/basisinterpolator.hh
index 498c272bdb94724721aaf5c4797e5255a8133a29..17144e76be8446163b2d926e4ca6d5aed7ce2c3f 100644
--- a/dune/fufem/functiontools/basisinterpolator.hh
+++ b/dune/fufem/functiontools/basisinterpolator.hh
@@ -9,7 +9,8 @@
#include <dune/localfunctions/common/virtualinterface.hh>
-#include <dune/fufem/arithmetic.hh>
+#include <dune/matrix-vector/axpy.hh>
+
#include <dune/fufem/functions/virtualgridfunction.hh>
#include <dune/fufem/functionspacebases/conformingbasis.hh>
@@ -215,7 +216,7 @@ struct BasisInterpolator<B, C, BasisGridFunction<ConformingBasis<B>, C> >
{
coeff[i] = 0.0;
for(size_t k=0; k<fBasis.constraints(i).size(); ++k)
- Arithmetic::addProduct(coeff[i], fBasis.constraints(i)[k].factor, coeff[fBasis.constraints(i)[k].index]);
+ Dune::MatrixVector::addProduct(coeff[i], fBasis.constraints(i)[k].factor, coeff[fBasis.constraints(i)[k].index]);
}
}
}
diff --git a/dune/fufem/geometry/convexpolyhedron.hh b/dune/fufem/geometry/convexpolyhedron.hh
index 4dbdfffe16ba95343adfd3e38adad5467909e22f..7b6efead5314b22f64499b6f9ccdd01d3b2f08e6 100644
--- a/dune/fufem/geometry/convexpolyhedron.hh
+++ b/dune/fufem/geometry/convexpolyhedron.hh
@@ -1,7 +1,7 @@
#ifndef DUNE_FUFEM_GEOMETRY_CONVEXPOLYHEDRON_HH
#define DUNE_FUFEM_GEOMETRY_CONVEXPOLYHEDRON_HH
-#include <dune/fufem/arithmetic.hh>
+#include <dune/matrix-vector/axpy.hh>
template <class Coordinate> struct ConvexPolyhedron {
std::vector<Coordinate> vertices;
@@ -11,7 +11,7 @@ template <class Coordinate> struct ConvexPolyhedron {
assert(b.size() == vertices.size());
Coordinate r(0);
for (size_t i = 0; i < b.size(); i++)
- Arithmetic::addProduct(r, b[i], vertices[i]);
+ Dune::MatrixVector::addProduct(r, b[i], vertices[i]);
return r;
}
diff --git a/dune/fufem/staticmatrixtools.hh b/dune/fufem/staticmatrixtools.hh
index d9ef27ec83abf68aa2e6934e5588c207b79e7128..ec9d5b78fa72bbd3968c69f3fe16c0b6035cec89 100644
--- a/dune/fufem/staticmatrixtools.hh
+++ b/dune/fufem/staticmatrixtools.hh
@@ -7,6 +7,9 @@
#include "dune/istl/bcrsmatrix.hh"
#include "dune/istl/matrixindexset.hh"
+#include <dune/matrix-vector/axpy.hh>
+#include <dune/matrix-vector/scalartraits.hh>
+
#include "arithmetic.hh"
namespace StaticMatrix
@@ -167,7 +170,7 @@ namespace StaticMatrix
{
static void addTransformedMatrix(MatrixA& A, const ScalarTransform& T1, const MatrixB& B, const ScalarTransform& T2)
{
- Arithmetic::addProduct(A, T1*T2, B);
+ Dune::MatrixVector::addProduct(A, T1*T2, B);
}
};
@@ -281,7 +284,7 @@ namespace StaticMatrix
void addTransformedMatrix(MatrixA& A, const TransformationMatrix& T1, const MatrixB& B, const TransformationMatrix& T2)
{
TransformMatrixHelper<MatrixA,MatrixB,TransformationMatrix,
- Arithmetic::ScalarTraits<MatrixA>::isScalar, Arithmetic::ScalarTraits<MatrixB>::isScalar, Arithmetic::ScalarTraits<TransformationMatrix>::isScalar>
+ Dune::MatrixVector::ScalarTraits<MatrixA>::isScalar, Dune::MatrixVector::ScalarTraits<MatrixB>::isScalar, Dune::MatrixVector::ScalarTraits<TransformationMatrix>::isScalar>
::addTransformedMatrix(A,T1,B,T2);
}
@@ -290,7 +293,7 @@ namespace StaticMatrix
{
A = 0;
TransformMatrixHelper<MatrixA,MatrixB,TransformationMatrix,
- Arithmetic::ScalarTraits<MatrixA>::isScalar, Arithmetic::ScalarTraits<MatrixB>::isScalar, Arithmetic::ScalarTraits<TransformationMatrix>::isScalar>
+ Dune::MatrixVector::ScalarTraits<MatrixA>::isScalar, Dune::MatrixVector::ScalarTraits<MatrixB>::isScalar, Dune::MatrixVector::ScalarTraits<TransformationMatrix>::isScalar>
::addTransformedMatrix(A,T1,B,T2);
}
diff --git a/dune/fufem/test/secondorderassemblertest.cc b/dune/fufem/test/secondorderassemblertest.cc
index 7b3a258fbdea38631ac7db18297a24b712744ea5..02a24f64331cfd47fff5de5baf9e22f286d2b495 100644
--- a/dune/fufem/test/secondorderassemblertest.cc
+++ b/dune/fufem/test/secondorderassemblertest.cc
@@ -12,6 +12,8 @@
#include <dune/istl/io.hh>
#include <dune/istl/scaledidmatrix.hh>
+#include <dune/matrix-vector/axpy.hh>
+
#include <dune/fufem/assemblers/operatorassembler.hh>
#include <dune/fufem/assemblers/localassemblers/secondorderoperatorassembler.hh>
@@ -130,7 +132,7 @@ int main (int argc, char *argv[])
{
LocalMatrix L(1);
- Arithmetic::addProduct(L, 1, LocalIdentityMatrix(1));
+ Dune::MatrixVector::addProduct(L, 1, LocalIdentityMatrix(1));
auto const leafView = grid.leafGridView();
auto funcL = [L](GlobalCoordinate const &) {return L;};