diff --git a/src/assemblers.cc b/src/assemblers.cc
new file mode 100644
index 0000000000000000000000000000000000000000..f8a6ca19ea74ead1bc08509f5b1b303f4494b46e
--- /dev/null
+++ b/src/assemblers.cc
@@ -0,0 +1,92 @@
+#include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/functions/constantfunction.hh>
+#include <dune/fufem/assemblers/localassemblers/neumannboundaryassembler.hh>
+#include <dune/fufem/assemblers/boundaryfunctionalassembler.hh>
+
+#include <dune/tectonic/globallaursennonlinearity.hh>
+#include <dune/tectonic/globalruinanonlinearity.hh>
+
+// Assembles Neumann boundary term in f
+template <class GridType, class GridView, class LocalVectorType, class FEBasis>
+void assemble_neumann(GridView const &gridView, FEBasis const &feBasis,
+ Dune::BitSetVector<1> const &neumannNodes,
+ Dune::BlockVector<LocalVectorType> &f,
+ Dune::VirtualFunction<double, double> const &neumann,
+ double time) { // constant sample function on neumann
+ // boundary
+ BoundaryPatch<GridView> const neumannBoundary(gridView, neumannNodes);
+ LocalVectorType SampleVector(0);
+ neumann.evaluate(time, SampleVector[0]);
+ SampleVector[1] = 0;
+ ConstantFunction<LocalVectorType, LocalVectorType> const fNeumann(
+ SampleVector);
+ NeumannBoundaryAssembler<GridType, LocalVectorType> neumannBoundaryAssembler(
+ fNeumann);
+
+ BoundaryFunctionalAssembler<FEBasis>(feBasis, neumannBoundary).assemble(
+ neumannBoundaryAssembler, f, true); // resize and zero output vector
+}
+
+// Assembles constant 1-function on frictional boundary in nodalIntegrals
+template <class GridType, class GridView, class LocalVectorType, class FEBasis>
+Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>>
+assemble_frictional(GridView const &gridView, FEBasis const &feBasis,
+ Dune::BitSetVector<1> const &frictionalNodes) {
+ typedef Dune::FieldVector<double, 1> Singleton;
+ BoundaryPatch<GridView> const frictionalBoundary(gridView, frictionalNodes);
+ ConstantFunction<LocalVectorType, Singleton> const constantOneFunction(1);
+ NeumannBoundaryAssembler<GridType, Singleton> frictionalBoundaryAssembler(
+ constantOneFunction);
+
+ auto const nodalIntegrals = Dune::make_shared<Dune::BlockVector<Singleton>>();
+ BoundaryFunctionalAssembler<FEBasis>(feBasis, frictionalBoundary)
+ .assemble(frictionalBoundaryAssembler, *nodalIntegrals,
+ true); // resize and zero output vector
+ return nodalIntegrals;
+}
+
+template <class VectorType, class MatrixType>
+Dune::shared_ptr<Dune::GlobalNonlinearity<VectorType, MatrixType> const>
+assemble_nonlinearity(
+ int size, Dune::ParameterTree const &parset,
+ Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>>
+ nodalIntegrals,
+ Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>> state,
+ double h) {
+ typedef Dune::BlockVector<Dune::FieldVector<double, 1>> SingletonVectorType;
+ // {{{ Assemble terms for the nonlinearity
+ auto mu = Dune::make_shared<SingletonVectorType>(size);
+ *mu = parset.get<double>("boundary.friction.mu");
+
+ auto normalStress = Dune::make_shared<SingletonVectorType>(size);
+ *normalStress = parset.get<double>("boundary.friction.normalstress");
+
+ std::string const friction_model =
+ parset.get<std::string>("boundary.friction.model");
+ if (friction_model == std::string("Ruina")) {
+ auto a = Dune::make_shared<SingletonVectorType>(size);
+ *a = parset.get<double>("boundary.friction.ruina.a");
+
+ auto eta = Dune::make_shared<SingletonVectorType>(size);
+ *eta = parset.get<double>("boundary.friction.eta");
+
+ auto b = Dune::make_shared<SingletonVectorType>(size);
+ *b = parset.get<double>("boundary.friction.ruina.b");
+
+ auto L = Dune::make_shared<SingletonVectorType>(size);
+ *L = parset.get<double>("boundary.friction.ruina.L");
+
+ return Dune::make_shared<
+ Dune::GlobalRuinaNonlinearity<VectorType, MatrixType> const>(
+ nodalIntegrals, a, mu, eta, normalStress, b, state, L, h);
+ } else if (friction_model == std::string("Laursen")) {
+ return
+ // TODO: take state and h into account
+ Dune::make_shared<
+ Dune::GlobalLaursenNonlinearity<Dune::LinearFunction, VectorType,
+ MatrixType> const>(mu, normalStress,
+ nodalIntegrals);
+ } else {
+ assert(false);
+ }
+}
diff --git a/src/assemblers.hh b/src/assemblers.hh
new file mode 100644
index 0000000000000000000000000000000000000000..196c6b7b31c77974d47d4b6a222c431f83e328c4
--- /dev/null
+++ b/src/assemblers.hh
@@ -0,0 +1,35 @@
+#ifndef ASSEMBLERS_HH
+#define ASSEMBLERS_HH
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/function.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/parametertree.hh>
+#include <dune/common/shared_ptr.hh>
+#include <dune/istl/bvector.hh>
+
+#include <dune/tectonic/globalnonlinearity.hh>
+
+template <class GridType, class GridView, class LocalVectorType, class FEBasis>
+void assemble_neumann(GridView const &gridView, FEBasis const &feBasis,
+ Dune::BitSetVector<1> const &neumannNodes,
+ Dune::BlockVector<LocalVectorType> &f,
+ Dune::VirtualFunction<double, double> const &neumann,
+ double time);
+
+template <class GridType, class GridView, class LocalVectorType, class FEBasis>
+Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>>
+assemble_frictional(GridView const &gridView, FEBasis const &feBasis,
+ Dune::BitSetVector<1> const &frictionalNodes);
+
+template <class VectorType, class MatrixType>
+Dune::shared_ptr<Dune::GlobalNonlinearity<VectorType, MatrixType> const>
+assemble_nonlinearity(
+ int size, Dune::ParameterTree const &parset,
+ Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>>
+ nodalIntegrals,
+ Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>> state,
+ double h);
+
+#include "assemblers.cc"
+#endif
diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc
index 464e8e08360bbd73d3088aa669472e743bc9a525..2c8be4f7bc50007615e2e50c5fc63ded8672e578 100644
--- a/src/one-body-sample.cc
+++ b/src/one-body-sample.cc
@@ -41,17 +41,13 @@
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>
-#include <dune/fufem/assemblers/boundaryfunctionalassembler.hh>
#include <dune/fufem/assemblers/functionalassembler.hh>
-#include <dune/fufem/assemblers/localassemblers/neumannboundaryassembler.hh>
#include <dune/fufem/assemblers/localassemblers/stvenantkirchhoffassembler.hh>
#include <dune/fufem/assemblers/localassemblers/vonmisesstressassembler.hh>
#include <dune/fufem/assemblers/operatorassembler.hh>
#include <dune/fufem/assemblers/transferoperatorassembler.hh>
-#include <dune/fufem/boundarypatch.hh>
#include <dune/fufem/dunepython.hh>
#include <dune/fufem/functions/basisgridfunction.hh>
-#include <dune/fufem/functions/constantfunction.hh>
#include <dune/fufem/functions/vtkbasisgridfunction.hh>
#include <dune/fufem/functionspacebases/p0basis.hh>
#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
@@ -66,12 +62,12 @@
#include <dune/tnnmg/iterationsteps/genericnonlineargs.hh>
#include <dune/tnnmg/iterationsteps/tnnmgstep.hh>
-#include <dune/tectonic/globallaursennonlinearity.hh>
-#include <dune/tectonic/globalruinanonlinearity.hh>
#include <dune/tectonic/myblockproblem.hh>
#include <dune/tectonic/myconvexproblem.hh>
#include <dune/tectonic/mydirectionalconvexfunction.hh>
+#include "assemblers.hh"
+
int const dim = 2;
template <class GridView>
@@ -112,91 +108,6 @@ void setup_boundary(GridView const &gridView,
std::cout << "Number of frictional nodes: " << dirichlet_nodes << std::endl;
}
-// Assembles Neumann boundary term in f
-template <class GridType, class GridView, class LocalVectorType, class FEBasis>
-void assemble_neumann(GridView const &gridView, FEBasis const &feBasis,
- Dune::BitSetVector<1> const &neumannNodes,
- Dune::BlockVector<LocalVectorType> &f,
- Dune::VirtualFunction<double, double> const &neumann,
- double time) { // constant sample function on neumann
- // boundary
- BoundaryPatch<GridView> const neumannBoundary(gridView, neumannNodes);
- LocalVectorType SampleVector(0);
- neumann.evaluate(time, SampleVector[0]);
- SampleVector[1] = 0;
- ConstantFunction<LocalVectorType, LocalVectorType> const fNeumann(
- SampleVector);
- NeumannBoundaryAssembler<GridType, LocalVectorType> neumannBoundaryAssembler(
- fNeumann);
-
- BoundaryFunctionalAssembler<FEBasis>(feBasis, neumannBoundary).assemble(
- neumannBoundaryAssembler, f, true); // resize and zero output vector
-}
-
-// Assembles constant 1-function on frictional boundary in nodalIntegrals
-template <class GridType, class GridView, class LocalVectorType, class FEBasis>
-Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>>
-assemble_frictional(GridView const &gridView, FEBasis const &feBasis,
- Dune::BitSetVector<1> const &frictionalNodes) {
- typedef Dune::FieldVector<double, 1> Singleton;
- BoundaryPatch<GridView> const frictionalBoundary(gridView, frictionalNodes);
- ConstantFunction<LocalVectorType, Singleton> const constantOneFunction(1);
- NeumannBoundaryAssembler<GridType, Singleton> frictionalBoundaryAssembler(
- constantOneFunction);
-
- auto const nodalIntegrals = Dune::make_shared<Dune::BlockVector<Singleton>>();
- BoundaryFunctionalAssembler<FEBasis>(feBasis, frictionalBoundary)
- .assemble(frictionalBoundaryAssembler, *nodalIntegrals,
- true); // resize and zero output vector
- return nodalIntegrals;
-}
-
-template <class VectorType, class MatrixType>
-Dune::shared_ptr<Dune::GlobalNonlinearity<VectorType, MatrixType> const>
-assemble_nonlinearity(
- int size, Dune::ParameterTree const &parset,
- Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>>
- nodalIntegrals,
- Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>> state,
- double h) {
- typedef Dune::BlockVector<Dune::FieldVector<double, 1>> SingletonVectorType;
- // {{{ Assemble terms for the nonlinearity
- auto mu = Dune::make_shared<SingletonVectorType>(size);
- *mu = parset.get<double>("boundary.friction.mu");
-
- auto normalStress = Dune::make_shared<SingletonVectorType>(size);
- *normalStress = parset.get<double>("boundary.friction.normalstress");
-
- std::string const friction_model =
- parset.get<std::string>("boundary.friction.model");
- if (friction_model == std::string("Ruina")) {
- auto a = Dune::make_shared<SingletonVectorType>(size);
- *a = parset.get<double>("boundary.friction.ruina.a");
-
- auto eta = Dune::make_shared<SingletonVectorType>(size);
- *eta = parset.get<double>("boundary.friction.eta");
-
- auto b = Dune::make_shared<SingletonVectorType>(size);
- *b = parset.get<double>("boundary.friction.ruina.b");
-
- auto L = Dune::make_shared<SingletonVectorType>(size);
- *L = parset.get<double>("boundary.friction.ruina.L");
-
- return Dune::make_shared<
- Dune::GlobalRuinaNonlinearity<VectorType, MatrixType> const>(
- nodalIntegrals, a, mu, eta, normalStress, b, state, L, h);
- } else if (friction_model == std::string("Laursen")) {
- return
- // TODO: take state and h into account
- Dune::make_shared<
- Dune::GlobalLaursenNonlinearity<Dune::LinearFunction, VectorType,
- MatrixType> const>(mu, normalStress,
- nodalIntegrals);
- } else {
- assert(false);
- }
-}
-
// The nonlinearity exp(-x)
class DecayingExponential {
public: