From 88209c04c0f6caf9b489381430f9ac13cd930dd7 Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Mon, 8 Feb 2016 13:51:07 +0100
Subject: [PATCH] Allow to set the Dirichlet vertices using a Python predicate
---
.../gradient-flow-curve-shortening.parset | 3 +++
src/gradient-flow.cc | 22 +++++++++++++++----
2 files changed, 21 insertions(+), 4 deletions(-)
diff --git a/problems/gradient-flow-curve-shortening.parset b/problems/gradient-flow-curve-shortening.parset
index 2412d9f1..cf5f31f4 100644
--- a/problems/gradient-flow-curve-shortening.parset
+++ b/problems/gradient-flow-curve-shortening.parset
@@ -54,6 +54,9 @@ energy = harmonic
# Inverse stereographic projection
initialIterate = initial-curve
+# Function of position returning '1' for each Dirichlet vertex
+dirichletVerticesPredicate = "x[0] < 0.0001 or x[0]>0.999"
+
# Number of time steps
numTimeSteps = 50
diff --git a/src/gradient-flow.cc b/src/gradient-flow.cc
index 02dc2723..9861fafb 100644
--- a/src/gradient-flow.cc
+++ b/src/gradient-flow.cc
@@ -151,11 +151,25 @@ int main (int argc, char *argv[]) try
typedef DuneFunctionsBasis<FEBasis> FufemFEBasis;
FufemFEBasis fufemFeBasis(feBasis);
- BitSetVector<1> allNodes(grid->size(dim));
- allNodes.setAll();
- BoundaryPatch<typename GridType::LeafGridView> dirichletBoundary(grid->leafGridView(), allNodes);
+ BitSetVector<1> dirichletVertices(grid->leafGridView().size(dim), false);
- BitSetVector<blocksize> dirichletNodes;
+ const auto& indexSet = grid->leafGridView().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,dimworld>, bool> pythonDirichletVertices(Python::evaluate(lambda));
+
+ for (auto&& vertex : vertices(grid->leafGridView()))
+ {
+ bool isDirichlet;
+ pythonDirichletVertices.evaluate(vertex.geometry().corner(0), isDirichlet);
+ dirichletVertices[indexSet.index(vertex)] = isDirichlet;
+ }
+
+ BoundaryPatch<GridType::LeafGridView> dirichletBoundary(grid->leafGridView(), dirichletVertices);
+
+ BitSetVector<blocksize> dirichletNodes(feBasis.indexSet().size(), false);
constructBoundaryDofs(dirichletBoundary,fufemFeBasis,dirichletNodes);
////////////////////////////
--
GitLab