diff --git a/problems/gradient-flow-curve-shortening.parset b/problems/gradient-flow-curve-shortening.parset index 2412d9f1566b0a46d2680ffbd23d2643f20c7150..cf5f31f492d398fe59d167365f950dcb87647e8e 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 02dc27238f21f626c5ded16b20186c68a2378274..9861fafb3cdf7168590eac9c28069d9a231b0eb1 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); ////////////////////////////