Add option for adaptive refinement, write the grid using an unstructured grid...

Add option for adaptive refinement, write the grid using an unstructured grid writer and create a file containing the final deformation
The grid and the file can then be read by dune-gfe

Make dune-vtk optional

dune.module

@@ -8,4 +8,4 @@ Version: 2.7-git
 Maintainer: oliver.sander@tu-dresden.de, youett@mi.fu-berlin.de
 #depending on
 Depends:dune-common dune-grid dune-istl dune-solvers dune-fufem dune-geometry dune-functions
-Suggests: dune-uggrid dune-parmg
+Suggests: dune-uggrid dune-parmg dune-vtk

problems/finite-strain-elasticity-bending.parset

@@ -3,14 +3,21 @@
 #############################################
 
 structuredGrid = true
+
+# bounding box
 lower = 0 0 0
-upper = 5 1 1
-elements = 20 5 5
+upper = 200 100 200
+
+elements = 10 5 5
+
+# Number of grid levels: refinement levels of the surface shell part, refined using hierarchic refinement
+numLevels = 2
+
 #Overlap indicator for the grid when doing parallel calculations; when set to false, all communication will be done assuming there is no overlap when assembling the matrix and the rhs
 overlap = false
 
-# Number of grid levels
-numLevels = 1
+# Adaptive Refinement
+adaptiveRefinement = true
 
 #############################################
 #  Solver parameters
@@ -23,13 +30,13 @@ numHomotopySteps = 1
 tolerance = 1e-6
 
 # Max number of steps of the trust region solver
-maxTrustRegionSteps = 500
+maxTrustRegionSteps = 200
 
 # Initial trust-region radius
-initialTrustRegionRadius = 0.1
+initialTrustRegionRadius = 10
 
 # Number of multigrid iterations per trust-region step
-numIt = 1000
+numIt = 100
 
 # Number of presmoothing steps
 nu1 = 3
@@ -56,31 +63,31 @@ baseTolerance = 1e-8
 #   Material parameters
 ############################
 
-energy = stvenantkirchhoff
+energy = mooneyrivlin
 
 ## For the Wriggers L-shape example
 [materialParameters]
 
-# Lame parameters
+# Lame parameters for stvenantkirchhoff
 # corresponds to E = 71240 N/mm^2, nu=0.31
 # However, we use units N/m^2
-mu = 2.7191e+6
-lambda = 4.4364e+6
+mu = 2.7191e+4
+lambda = 4.4364e+4
 
 mooneyrivlin_a = 2.7191e+6
 mooneyrivlin_b = 2.7191e+6
 mooneyrivlin_c = 2.7191e+6 
-mooneyrivlin_10 = -7.28e+5
-mooneyrivlin_01 = 9.17e+5
-mooneyrivlin_20 = 1.23e+5
-mooneyrivlin_02 = 8.15e+5
-mooneyrivlin_11 = -5.14e+5
+mooneyrivlin_10 = -7.28e+4
+mooneyrivlin_01 = 9.17e+4
+mooneyrivlin_20 = 1.23e+4
+mooneyrivlin_02 = 8.15e+4
+mooneyrivlin_11 = -5.14e+4
 mooneyrivlin_30 = 0
 mooneyrivlin_21 = 0
 mooneyrivlin_12 = 0
 mooneyrivlin_03 = 0
 
-mooneyrivlin_k = 1e+6	# volume-preserving parameter
+mooneyrivlin_k = 1e+5	# volume-preserving parameter
 mooneyrivlin_energy = log # log, square or ciarlet; different ways to compute the Mooney-Rivlin-Energy
 
 []
@@ -95,12 +102,16 @@ problem = wriggers-l-shape
 # x is the vertex coordinate
 dirichletVerticesPredicate = "x[0] < 0.01"
 
-###  Python predicate specifying all Dirichlet grid vertices
+###  Python predicate specifying all surface shell grid vertices
+# x is the vertex coordinate
+adaptiveRefinementPredicate = "x[2] > 199.99"
+
+###  Python predicate specifying all Neumann grid vertices
 # x is the vertex coordinate
-neumannVerticesPredicate = "x[0] > 4.99"
+neumannVerticesPredicate = "x[0] > 199.99"
 
 ###  Neumann values, if needed
-neumannValues = 0 5e4 0
+neumannValues = 15e3 0 0
 
 # Initial deformation
-initialDeformation = "[x[0], x[1], x[2]]"
+initialDeformation = "x"

src/finite-strain-elasticity.cc

@@ -42,6 +42,13 @@
 #include <dune/elasticity/materials/neumannenergy.hh>
 #include <dune/elasticity/materials/sumenergy.hh>
 
+#if HAVE_DUNE_VTK
+#include <dune/vtk/writers/vtkunstructuredgridwriter.hh>
+#endif
+
+#include <iostream>
+#include <fstream>
+
 // grid dimension
 const int dim = 3;
 
@@ -97,7 +104,7 @@ int main (int argc, char *argv[]) try
   ParameterTreeParser::readOptions(argc, argv, parameterSet);
 
   // read solver settings
-  const int numLevels                   = parameterSet.get<int>("numLevels");
+  int numLevels                         = parameterSet.get<int>("numLevels");
   int numHomotopySteps                  = parameterSet.get<int>("numHomotopySteps");
   const double tolerance                = parameterSet.get<double>("tolerance");
   const int maxTrustRegionSteps         = parameterSet.get<int>("maxTrustRegionSteps");
@@ -121,6 +128,8 @@ int main (int argc, char *argv[]) try
 
   std::shared_ptr<GridType> grid;
 
+  const bool adaptiveRefinement = parameterSet.get<bool>("adaptiveRefinement");
+
   FieldVector<double,dim> lower(0), upper(1);
 
   if (parameterSet.get<bool>("structuredGrid")) {
@@ -129,18 +138,50 @@ int main (int argc, char *argv[]) try
     upper = parameterSet.get<FieldVector<double,dim> >("upper");
 
     std::array<unsigned int,dim> elements = parameterSet.get<std::array<unsigned int,dim> >("elements");
-    grid = StructuredGridFactory<GridType>::createCubeGrid(lower, upper, elements);
 
+    grid = StructuredGridFactory<GridType>::createCubeGrid(lower, upper, elements);
   } else {
     std::string path                = parameterSet.get<std::string>("path");
     std::string gridFile            = parameterSet.get<std::string>("gridFile");
     grid = std::shared_ptr<GridType>(GmshReader<GridType>::read(path + "/" + gridFile));
   }
 
-  grid->globalRefine(numLevels-1);
+  // Make Python function that computes which vertices are on the Dirichlet boundary,
+  // based on the vertex positions.
+  std::string lambda = std::string("lambda x: (") + parameterSet.get<std::string>("dirichletVerticesPredicate") + std::string(")");
+  PythonFunction<FieldVector<double,dim>, bool> pythonDirichletVertices(Python::evaluate(lambda));
+
+  // Same for the Neumann boundary
+  lambda = std::string("lambda x: (") + parameterSet.get<std::string>("neumannVerticesPredicate", "0") + std::string(")");
+  PythonFunction<FieldVector<double,dim>, bool> pythonNeumannVertices(Python::evaluate(lambda));
+
+  if (adaptiveRefinement) {
+    // Same for the boundary that will get adaptively refined
+    lambda = std::string("lambda x: (") + parameterSet.get<std::string>("adaptiveRefinementPredicate", "0") + std::string(")");
+    PythonFunction<FieldVector<double,dim>, bool> pythonSurfaceShellVertices(Python::evaluate(lambda));
+
+    grid->setRefinementType(GridType::RefinementType::COPY);
+
+    while (numLevels > 0) {
+      for (auto&& e : elements(grid->leafGridView())){
+        bool isSurfaceShell = false;
+        for (int i = 0; i < e.geometry().corners(); i++) {
+            isSurfaceShell = isSurfaceShell || pythonSurfaceShellVertices(e.geometry().corner(i));
+        }
+        grid->mark(isSurfaceShell ? 1 : 0,e);
+      }
+
+      grid->adapt();
 
       grid->loadBalance();
 
+      numLevels--;
+    }
+  } else {
+    grid->globalRefine(numLevels-1);
+    grid->loadBalance();
+  }
+
   if (mpiHelper.rank()==0)
     std::cout << "There are " << grid->leafGridView().comm().size() << " processes" << std::endl;
   
@@ -165,15 +206,6 @@ int main (int argc, char *argv[]) try
 
   const GridView::IndexSet& indexSet = gridView.indexSet();
 
-  // Make Python function that computes which vertices are on the Dirichlet boundary,
-  // based on the vertex positions.
-  std::string lambda = std::string("lambda x: (") + parameterSet.get<std::string>("dirichletVerticesPredicate") + std::string(")");
-  PythonFunction<FieldVector<double,dim>, bool> pythonDirichletVertices(Python::evaluate(lambda));
-
-  // Same for the Neumann boundary
-  lambda = std::string("lambda x: (") + parameterSet.get<std::string>("neumannVerticesPredicate", "0") + std::string(")");
-  PythonFunction<FieldVector<double,dim>, bool> pythonNeumannVertices(Python::evaluate(lambda));
-
   for (auto&& v : vertices(gridView))
   {
     bool isDirichlet;
@@ -384,7 +416,25 @@ int main (int argc, char *argv[]) try
     vtkWriter.addVertexData(displacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::vector, dim));
     vtkWriter.write(resultPath + "finite-strain_homotopy_" + std::to_string(i+1));
   }
-
+#if HAVE_DUNE_VTK
+  if (grid->leafGridView().comm().size() == 1) {
+    displacement = x;
+    displacement -= identity;
+    auto finalDisplacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(feBasis, displacement);
+    //Write out the grid using an unstructured GridWriter
+    VtkUnstructuredGridWriter<GridView> vtkUnstructuredGridWriter(gridView, Vtk::ASCII);
+    vtkUnstructuredGridWriter.addPointData(finalDisplacementFunction, "displacement");
+    vtkUnstructuredGridWriter.write(resultPath + "finite-strain_homotopy_original");
+    std::ofstream file;
+    file.open("deformation");
+    for (auto& d : displacement){
+      file << d << "\n";
+    }
+    file.close();
+  } else if (mpiHelper.rank() == 0) {
+    std::cout << "Not writing out the deformation to a file, this currently only works for sequential calculations" << std::endl;
+  }
+#endif
   if (mpiHelper.rank()==0)
     std::cout << "Complete duration: " << homotopyTimer.elapsed() << " sec." << std::endl;