Play around with the StVKA

src/Makefile.am

-
 SUBDIRS =
 
 check_PROGRAMS = \
@@ -6,7 +5,11 @@ check_PROGRAMS = \
 	test-gradient-method
 
 bin_PROGRAMS = \
-	gauss-seidel-sample
+	gauss-seidel-sample \
+	one-body-sample
+
+one_body_sample_SOURCES = \
+	one-body-sample.cc
 
 gauss_seidel_sample_SOURCES = \
 	gauss-seidel-sample.cc gaussseidel.hh

src/one-body-sample.cc

+/* -*- mode:c++; mode:semantic -*- */
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/common/exceptions.hh>
+
+#include <dune/common/fmatrix.hh>
+#include <dune/istl/bcrsmatrix.hh>
+
+#include <dune/grid/yaspgrid.hh>
+
+#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
+#include <dune/fufem/assemblers/operatorassembler.hh>
+#include <dune/fufem/assemblers/localassemblers/stvenantkirchhoffassembler.hh>
+
+#include <exception>
+
+int const dim = 2;
+
+int main() {
+  try {
+    typedef Dune::BCRSMatrix<Dune::FieldMatrix<double, dim, dim>> OperatorType;
+
+    // FIXME: Random values
+    double const E = 1e8;
+    double const nu = 0.3;
+    int const refinements = 5;
+
+    typedef Dune::YaspGrid<dim> GridType;
+
+    Dune::FieldVector<double, dim> const end_points(
+        1); // nth dimension (zero-indexed) goes from 0 to end_points[n]
+    Dune::FieldVector<int, dim> const elements(
+        2); // number of elements in each direction
+    Dune::FieldVector<bool, dim> const periodic(false);
+
+    GridType grid(end_points, elements, periodic, 0);
+
+    grid.globalRefine(refinements);
+
+    typedef P1NodalBasis<GridType::LeafGridView, double> P1Basis;
+    P1Basis const p1Basis(grid.leafView());
+
+    OperatorAssembler<P1Basis, P1Basis> const globalAssembler(p1Basis, p1Basis);
+
+    StVenantKirchhoffAssembler<GridType, P1Basis::LocalFiniteElement,
+                               P1Basis::LocalFiniteElement> const
+    localStiffness(E, nu);
+
+    OperatorType stiffnessMatrix;
+    globalAssembler.assemble(localStiffness, stiffnessMatrix);
+
+    return 0;
+  }
+  catch (Dune::Exception& e) {
+    Dune::derr << "Dune reported error: " << e << std::endl;
+  }
+  catch (std::exception& e) {
+    std::cout << "Standard exception: " << e.what() << std::endl;
+  }
+}