Merge branch 'feature/add-option-to-run-film-on-substrate-without-any-shell' into 'master'

Introduce python function for adaptive refinement

See merge request !86

problems/film-on-substrate.parset

@@ -40,7 +40,10 @@ dirichletValues = identity-dirichlet-values
 # x is the vertex coordinate
 dirichletVerticesPredicate = "[x[0] < 0.01, x[0] < 0.01 or x[0] > 199.99, x[0] < 0.01 or x[0] > 199.99]"
 
-###  Python predicate specifying all surfaceshell grid vertices, elements conataining these vertices will get adaptively refined
+###  Python predicate specifying the vertices whose elements will get adaptively refined 
+adaptiveRefinementVerticesPredicate = "x[2] > 199.99 and x[0] > 49.99 and x[0] < 150.01"
+
+###  Python predicate specifying all surfaceshell grid vertices, elements conataining these vertices will get ALSO get adaptively refined
 # x is the vertex coordinate
 surfaceShellVerticesPredicate = "x[2] > 199.99 and x[0] > 49.99 and x[0] < 150.01"
 
@@ -136,24 +139,53 @@ b1 = 1
 b2 = 1
 b3 = 1
 
-# 
-mooneyrivlin_10 = -1.67e+6 #184 2:1
+
+mooneyrivlin_energy = log
+# log, square or ciarlet; different ways to compute the Mooney-Rivlin-Energy
+
+# ciarlet: Formula from "Ciarlet: Three-Dimensional Elasticity": not tested thoroughly yet
+# log: Generalized Rivlin model or polynomial hyperelastic model, using  0.5*mooneyrivlin_k*log(det(∇φ)) as the volume-preserving penalty term
+# square: Generalized Rivlin model or polynomial hyperelastic model, using mooneyrivlin_k*(det(∇φ)-1)² as the volume-preserving penalty term
+
+# volume-preserving parameter:
+# 184 2:1, mooneyrivlin_k = 57e+6 and mooneyrivlin_energy = log, the neumannValues = 27e4 0 0 result in a stretch of 30% in x-direction
+# 184 10:1, mooneyrivlin_k = 90e+6 and mooneyrivlin_energy = log, the neumannValues = 27e4 0 0 result in a stretch of 30% in x-direction
+
+#182 2:1
+#mooneyrivlin_10 = -3.01e+6 
+#mooneyrivlin_01 = 3.36e+6
+#mooneyrivlin_20 = 5.03e+6
+#mooneyrivlin_02 = 13.1e+6
+#mooneyrivlin_11 = -15.2e+6
+#mooneyrivlin_k = ???
+
+#182 20:1
+#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_k = ???
+
+#184 2:1 ANDRE
+mooneyrivlin_10 = -1.67e+6 
 mooneyrivlin_01 = 1.94e+6
 mooneyrivlin_20 = 2.42e+6
 mooneyrivlin_02 = 6.52e+6
 mooneyrivlin_11 = -7.34e+6
+mooneyrivlin_k = 57e+6
+
+#184 10:1 ANIK
+#mooneyrivlin_10 = -1.44e+6 
+#mooneyrivlin_01 = 1.95e+6
+#mooneyrivlin_20 = 2.65e+6
+#mooneyrivlin_02 = 5.06e+6
+#mooneyrivlin_11 = -6.56e+6
+#mooneyrivlin_k = 90e+6
 
 mooneyrivlin_30 = 0
 mooneyrivlin_21 = 0
 mooneyrivlin_12 = 0
 mooneyrivlin_03 = 0
 
-# volume-preserving parameter 
-mooneyrivlin_k = 57e+6	# 184 2:1, mooneyrivlin_k = 57e+6 and mooneyrivlin_energy = log, the neumannValues = 27e4 0 0 result in a stretch of 30% of 45e4 10e4 2e4 in x-direction, so a stretch of 45e4*0.3 = 13.5e4
-mooneyrivlin_energy = log # log, square or ciarlet; different ways to compute the Mooney-Rivlin-Energy
-
-# ciarlet: Fomula from "Ciarlet: Three-Dimensional Elasticity", here no penalty term is 
-# log: Generalized Rivlin model or polynomial hyperelastic model, using  0.5*mooneyrivlin_k*log(det(∇φ)) as the volume-preserving penalty term
-# square: Generalized Rivlin model or polynomial hyperelastic model, using mooneyrivlin_k*(det(∇φ)-1)² as the volume-preserving penalty term
-
 []

src/film-on-substrate.cc

@@ -176,13 +176,17 @@ int main (int argc, char *argv[]) try
   lambda = std::string("lambda x: (") + parameterSet.get<std::string>("surfaceShellVerticesPredicate", "0") + std::string(")");
   auto pythonSurfaceShellVertices = Python::make_function<bool>(Python::evaluate(lambda));
 
+  // Same for adaptive refinement vertices
+  lambda = std::string("lambda x: (") + parameterSet.get<std::string>("adaptiveRefinementVerticesPredicate", "0") + std::string(")");
+  auto pythonAdaptiveRefinementVertices = Python::make_function<bool>(Python::evaluate(lambda));
+
   while (numLevels > 0) {
     for (auto&& e : elements(grid->leafGridView())){
-      bool isSurfaceShell = false;
+      bool refineHere = false;
       for (int i = 0; i < e.geometry().corners(); i++) {
-          isSurfaceShell = isSurfaceShell || pythonSurfaceShellVertices(e.geometry().corner(i));
+          refineHere = refineHere || pythonSurfaceShellVertices(e.geometry().corner(i)) || pythonAdaptiveRefinementVertices(e.geometry().corner(i));
       }
-      grid->mark(isSurfaceShell ? 1 : 0,e);
+      grid->mark(refineHere ? 1 : 0, e);
     }
 
     grid->adapt();
@@ -329,9 +333,12 @@ int main (int argc, char *argv[]) try
 
   BlockVector<FieldVector<double,dim> > identity(compositeBasis.size({0}));
   Dune::Functions::interpolate(deformationPowerBasis, identity, [](FieldVector<double,dim> x){ return x; });
-
+  
+  double max_x = 0;
+  double initial_max_x = 0;
   for (int i = 0; i < displacement.size(); i++) {
     x[_0][i] = displacement[i]; //Copy over the initial deformation to the deformation part of x
+    initial_max_x = std::max(x[_0][i][0], initial_max_x);
     displacement[i] -= identity[i]; //Subtract identity to get the initial displacement as a function
   }
   auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(deformationPowerBasis, displacement);
@@ -650,6 +657,14 @@ int main (int argc, char *argv[]) try
     vtkWriter.addVertexData(displacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::scalar, dim));
     vtkWriter.write(resultPath + "finite-strain_homotopy_" + parameterSet.get<std::string>("energy") + "_" + std::to_string(neumannValues[0]) + "_" + std::to_string(i+1));
   }
+  for (int i = 0; i < x[_0].size(); i++) {
+    max_x = std::max(x[_0][i][0], max_x);
+  }
+
+  if (mpiHelper.rank()==0) {
+    std::cout << "Maximal value in x-direction: " << max_x << ", this is a stretch in of " << 100*(max_x - initial_max_x)/initial_max_x << " %." << std::endl;
+  }
+
   std::string ending = grid->leafGridView().comm().size() > 1 ? std::to_string(mpiHelper.rank()) : "";
   std::ofstream file;
   std::string pathToOutput = parameterSet.hasKey("pathToOutput") ?  parameterSet.get<std::string>("pathToOutput") : "./";