From 063194fa4b94869e190ef72fe5a91afc21d5a7d6 Mon Sep 17 00:00:00 2001
From: Markus Blatt <mblatt@dune-project.org>
Date: Tue, 15 Nov 2005 10:40:49 +0000
Subject: [PATCH] Laplacian finite difference 5 point stencil.
[[Imported from SVN: r402]]
---
istl/test/laplacian.hh | 65 ++++++++++++++++++++++++++++++++++++++++++
1 file changed, 65 insertions(+)
create mode 100644 istl/test/laplacian.hh
diff --git a/istl/test/laplacian.hh b/istl/test/laplacian.hh
new file mode 100644
index 000000000..6e3020f57
--- /dev/null
+++ b/istl/test/laplacian.hh
@@ -0,0 +1,65 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef LAPLACIAN_HH
+#define LAPLACIAN_HH
+
+template<int N, class B>
+void setupSparsityPattern(Dune::BCRSMatrix<B>& A)
+{
+ for (typename Dune::BCRSMatrix<B>::CreateIterator i = A.createbegin(); i != A.createend(); ++i) {
+ int x = i.index()%N; // x coordinate in the 2d field
+ int y = i.index()/N; // y coordinate in the 2d field
+
+ if(y>0)
+ // insert lower neighbour
+ i.insert(i.index()-N);
+ if(x>0)
+ // insert left neighbour
+ i.insert(i.index()-1);
+
+ // insert diagonal value
+ i.insert(i.index());
+
+ if(x<N-1)
+ //insert right neighbour
+ i.insert(i.index()+1);
+ if(y<N-1)
+ // insert upper neighbour
+ i.insert(i.index()+N);
+ }
+}
+
+template<int N, class B>
+void setupLaplacian(Dune::BCRSMatrix<B>& A)
+{
+ setupSparsityPattern<N>(A);
+
+ B diagonal = 0, bone=0, beps=0;
+ for(typename B::RowIterator b = diagonal.begin(); b != diagonal.end(); ++b)
+ b->operator[](b.index())=4;
+
+
+ for(typename B::RowIterator b=bone.begin(); b != bone.end(); ++b)
+ b->operator[](b.index())=-1.0;
+
+
+ for (typename Dune::BCRSMatrix<B>::RowIterator i = A.begin(); i != A.end(); ++i) {
+ int x = i.index()%N; // x coordinate in the 2d field
+ int y = i.index()/N; // y coordinate in the 2d field
+
+ i->operator[](i.index())=diagonal;
+
+ if(y>0)
+ i->operator[](i.index()-N)=bone;
+
+ if(y<N-1)
+ i->operator[](i.index()+N)=bone;
+
+ if(x>0)
+ i->operator[](i.index()-1)=bone;
+
+ if(x < N-1)
+ i->operator[](i.index()+1)=bone;
+ }
+}
+#endif
--
GitLab