From 27006393436564e766f3670d2a1db7c735024c61 Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Tue, 2 Jan 2018 11:51:39 +0100
Subject: [PATCH] Determine Dirichlet boundary from a Python predicate

---
 src/harmonicmaps.cc | 21 +++++++++++++++++----
 1 file changed, 17 insertions(+), 4 deletions(-)

diff --git a/src/harmonicmaps.cc b/src/harmonicmaps.cc
index b6260d0d..09b50c94 100644
--- a/src/harmonicmaps.cc
+++ b/src/harmonicmaps.cc
@@ -170,12 +170,25 @@ int main (int argc, char *argv[])
     // /////////////////////////////////////////
     //   Read Dirichlet values
     // /////////////////////////////////////////
+    BitSetVector<1> dirichletVertices(gridView.size(dim), false);
+    const GridView::IndexSet& indexSet = gridView.indexSet();
 
-    BitSetVector<1> allNodes(grid->size(dim));
-    allNodes.setAll();
-    BoundaryPatch<GridView> dirichletBoundary(gridView, allNodes);
+    // 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));
 
-    BitSetVector<blocksize> dirichletNodes;
+    for (auto&& vertex : vertices(gridView))
+    {
+      //bool isDirichlet;
+      bool isDirichlet = pythonDirichletVertices(vertex.geometry().corner(0));
+      pythonDirichletVertices.evaluate(vertex.geometry().corner(0), isDirichlet);
+      dirichletVertices[indexSet.index(vertex)] = isDirichlet;
+    }
+
+    BoundaryPatch<GridView> dirichletBoundary(gridView, dirichletVertices);
+
+    BitSetVector<blocksize> dirichletNodes(feBasis.size(), false);
 
     typedef DuneFunctionsBasis<FEBasis> FufemFEBasis;
     FufemFEBasis fufemFeBasis(feBasis);
-- 
GitLab