diff --git a/problems/film-on-substrate.parset b/problems/film-on-substrate.parset
index 8ccce0d1c12ba5b5145e813910dade3ba1281cea..e6cf6ce9164fe7216fc7b8392d764965e02aa2c0 100644
--- a/problems/film-on-substrate.parset
+++ b/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
-
[]
diff --git a/src/film-on-substrate.cc b/src/film-on-substrate.cc
index dfef58c8317ab30f3b84d97e3726488f65bb7e51..8978172ac9fe81b8ce9c8a0b68e5cddc08a9961f 100644
--- a/src/film-on-substrate.cc
+++ b/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") : "./";