Pass and use global nonlinearity

dune/tectonic/myblockproblem.hh

@@ -19,6 +19,7 @@
 template <class MyConvexProblemTypeTEMPLATE> class MyBlockProblem {
 public:
   typedef MyConvexProblemTypeTEMPLATE MyConvexProblemType;
+  typedef typename MyConvexProblemType::FunctionType FunctionType;
   typedef typename MyConvexProblemType::NonlinearityType NonlinearityType;
   typedef typename MyConvexProblemType::VectorType VectorType;
   typedef typename MyConvexProblemType::MatrixType MatrixType;
@@ -98,9 +99,9 @@ class MyBlockProblem<MyConvexProblemTypeTEMPLATE>::IterateObject {
       }
       assert(localA != NULL);
 
-      // FIXME: Hardcoding a fixed function here for now
-      Dune::TrivialFunction func;
-      Dune::MyNonlinearity<block_size> phi(func);
+      FunctionType f;
+      problem.newphi.restriction(m, f);
+      Dune::MyNonlinearity<block_size> const phi(f);
       Dune::SampleFunctional<block_size> localJ(*localA, localb, phi);
 
       LocalVectorType correction;

dune/tectonic/myconvexproblem.hh

@@ -12,11 +12,11 @@
    part
 */
 template <class NonlinearityTypeTEMPLATE, class MatrixTypeTEMPLATE,
-          class VectorTypeTEMPLATE>
+          class VectorTypeTEMPLATE, class FunctionTypeTEMPLATE>
 class MyConvexProblem {
 public:
   typedef NonlinearityTypeTEMPLATE NonlinearityType;
-
+  typedef FunctionTypeTEMPLATE FunctionType;
   typedef VectorTypeTEMPLATE VectorType;
   typedef MatrixTypeTEMPLATE MatrixType;
   typedef typename VectorType::block_type LocalVectorType;
@@ -31,12 +31,14 @@ class MyConvexProblem {
       \param f The linear functional
       \param u The solution vector
   */
-  MyConvexProblem(MatrixType const &A, NonlinearityType &phi,
-                  VectorType const &f, VectorType &u)
-      : A(A), phi(phi), f(f), u(u) {};
+  MyConvexProblem(MatrixType const &A,
+                  Dune::MyGlobalNonlinearity<block_size, FunctionType> &newphi,
+                  NonlinearityType &phi, VectorType const &f, VectorType &u)
+      : A(A), newphi(newphi), phi(phi), f(f), u(u) {};
 
   MatrixType const &A;
   NonlinearityType &phi;
+  Dune::MyGlobalNonlinearity<block_size, FunctionType> const &newphi;
 
   VectorType const &f;
 

dune/tectonic/myglobalnonlinearity.hh

@@ -18,7 +18,7 @@ template <int dim, class OuterFunctionType> class MyGlobalNonlinearity {
         normalStress(normalStress),
         nodalIntegrals(nodalIntegrals) {}
 
-  void restriction(int i, OuterFunctionType &f) {
+  void restriction(int i, OuterFunctionType &f) const {
     double coefficient = normalStress[i] * nodalIntegrals[i];
     // FIXME: Assume Gamma = id and h_{n+1} = 1 for now;
     // We then only have to evaluate (1 + F_xi)(|x|)

src/one-body-sample.cc

@@ -36,6 +36,7 @@
 #include <dune/solvers/norms/energynorm.hh>
 #include <dune/solvers/solvers/loopsolver.hh>
 
+#include <dune/tectonic/myglobalnonlinearity.hh>
 #include <dune/tectonic/myconvexproblem.hh>
 #include <dune/tectonic/myblockproblem.hh>
 
@@ -184,6 +185,10 @@ int main() {
           neumannBoundaryAssembler, nodalIntegrals,
           true); // resize the output vector and zero all of its entries
 
+      // FIXME: We should be using S_F instead of S_N here
+      Dune::MyGlobalNonlinearity<dim, Dune::LinearFunction>
+      myGlobalNonlinearity(coefficientOfFriction, normalStress, nodalIntegrals);
+
       { // constant 2D function
         std::vector<SmallVector> b;
         SmallVector SampleVector;
@@ -205,10 +210,11 @@ int main() {
       typedef ZeroNonlinearity<SmallVector, SmallMatrix> NonlinearityType;
       NonlinearityType phi;
 
-      typedef MyConvexProblem<NonlinearityType, OperatorType, VectorType>
-      MyConvexProblemType;
+      typedef MyConvexProblem<NonlinearityType, OperatorType, VectorType,
+                              Dune::LinearFunction> MyConvexProblemType;
 
-      MyConvexProblemType myConvexProblem(stiffnessMatrix, phi, f, u1);
+      MyConvexProblemType myConvexProblem(stiffnessMatrix, myGlobalNonlinearity,
+                                          phi, f, u1);
 
       typedef MyBlockProblem<MyConvexProblemType> MyBlockProblemType;
       MyBlockProblemType myBlockProblem(myConvexProblem);