From add56cb21b61245df065893a40adf8568f9b7b25 Mon Sep 17 00:00:00 2001
From: Lisa Julia Nebel <lisa_julia.nebel@tu-dresden.de>
Date: Thu, 16 Jul 2020 20:38:45 +0200
Subject: [PATCH] Fix Neumann boundary in y-and-z diretion
---
problems/film-on-substrate.parset | 2 +-
src/film-on-substrate.cc | 38 ++++++++++++++++++++-----------
2 files changed, 26 insertions(+), 14 deletions(-)
diff --git a/problems/film-on-substrate.parset b/problems/film-on-substrate.parset
index fb745548..dd3d24ac 100644
--- a/problems/film-on-substrate.parset
+++ b/problems/film-on-substrate.parset
@@ -30,7 +30,7 @@ dirichletValues = identity-dirichlet-values
### Python predicate specifying all Dirichlet grid vertices
# x is the vertex coordinate
-dirichletVerticesPredicate = "x[0] < 0.01"
+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
# x is the vertex coordinate
diff --git a/src/film-on-substrate.cc b/src/film-on-substrate.cc
index 94c64915..235be15b 100644
--- a/src/film-on-substrate.cc
+++ b/src/film-on-substrate.cc
@@ -177,7 +177,7 @@ int main (int argc, char *argv[]) try
// 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));
+ PythonFunction<FieldVector<double,dim>, FieldVector<bool,dim>> pythonDirichletVertices(Python::evaluate(lambda));
// Same for the Neumann boundary
lambda = std::string("lambda x: (") + parameterSet.get<std::string>("neumannVerticesPredicate", "0") + std::string(")");
@@ -253,7 +253,9 @@ int main (int argc, char *argv[]) try
// BOUNDARY DATA
/////////////////////////////////////////////////////////////
- BitSetVector<1> dirichletVertices(gridView.size(dim), false);
+ BitSetVector<1> dirichletVerticesX(gridView.size(dim), false);
+ BitSetVector<1> dirichletVerticesY(gridView.size(dim), false);
+ BitSetVector<1> dirichletVerticesZ(gridView.size(dim), false);
BitSetVector<1> neumannVertices(gridView.size(dim), false);
BitSetVector<1> surfaceShellVertices(gridView.size(dim), false);
@@ -261,8 +263,10 @@ int main (int argc, char *argv[]) try
for (auto&& v : vertices(gridView))
{
- bool isDirichlet = pythonDirichletVertices(v.geometry().corner(0));
- dirichletVertices[indexSet.index(v)] = isDirichlet;
+ FieldVector<bool,dim> isDirichlet = pythonDirichletVertices(v.geometry().corner(0));
+ dirichletVerticesX[indexSet.index(v)] = isDirichlet[0];
+ dirichletVerticesY[indexSet.index(v)] = isDirichlet[1];
+ dirichletVerticesZ[indexSet.index(v)] = isDirichlet[2];
bool isNeumann = pythonNeumannVertices(v.geometry().corner(0));
neumannVertices[indexSet.index(v)] = isNeumann;
@@ -271,17 +275,25 @@ int main (int argc, char *argv[]) try
surfaceShellVertices[indexSet.index(v)] = isSurfaceShell;
}
- BoundaryPatch<GridView> dirichletBoundary(gridView, dirichletVertices);
+ BoundaryPatch<GridView> dirichletBoundaryX(gridView, dirichletVerticesX);
+ BoundaryPatch<GridView> dirichletBoundaryY(gridView, dirichletVerticesY);
+ BoundaryPatch<GridView> dirichletBoundaryZ(gridView, dirichletVerticesZ);
auto neumannBoundary = std::make_shared<BoundaryPatch<GridView>>(gridView, neumannVertices);
BoundaryPatch<GridView> surfaceShellBoundary(gridView, surfaceShellVertices);
- std::cout << "On rank " << mpiHelper.rank() << ": Dirichlet boundary has " << dirichletBoundary.numFaces() << " faces\n";
+ std::cout << "On rank " << mpiHelper.rank() << ": Dirichlet boundary has [" << dirichletBoundaryX.numFaces() <<
+ ", " << dirichletBoundaryY.numFaces() <<
+ ", " << dirichletBoundaryZ.numFaces() <<"] faces\n";
std::cout << "On rank " << mpiHelper.rank() << ": Neumann boundary has " << neumannBoundary->numFaces() << " faces\n";
std::cout << "On rank " << mpiHelper.rank() << ": Shell boundary has " << surfaceShellBoundary.numFaces() << " faces\n";
- BitSetVector<1> dirichletNodes(compositeBasis.size({0}),false);
+ BitSetVector<1> dirichletNodesX(compositeBasis.size({0}),false);
+ BitSetVector<1> dirichletNodesY(compositeBasis.size({0}),false);
+ BitSetVector<1> dirichletNodesZ(compositeBasis.size({0}),false);
BitSetVector<1> surfaceShellNodes(compositeBasis.size({1}),false);
- constructBoundaryDofs(dirichletBoundary,deformationFEBasis,dirichletNodes);
+ constructBoundaryDofs(dirichletBoundaryX,deformationFEBasis,dirichletNodesX);
+ constructBoundaryDofs(dirichletBoundaryY,deformationFEBasis,dirichletNodesY);
+ constructBoundaryDofs(dirichletBoundaryZ,deformationFEBasis,dirichletNodesZ);
constructBoundaryDofs(surfaceShellBoundary,orientationFEBasis,surfaceShellNodes);
//Create BitVector matching the tangential space
@@ -292,14 +304,14 @@ int main (int argc, char *argv[]) try
BitVector dirichletDofs;
dirichletDofs[_0].resize(compositeBasis.size({0}));
dirichletDofs[_1].resize(compositeBasis.size({1}));
- for (size_t i = 0; i < compositeBasis.size({0}); i++){
- for (int j = 0; j < dim; j++){
- dirichletDofs[_0][i][j] = dirichletNodes[i][0];
- }
+ for (size_t i = 0; i < compositeBasis.size({0}); i++) {
+ dirichletDofs[_0][i][0] = dirichletNodesX[i][0];
+ dirichletDofs[_0][i][1] = dirichletNodesY[i][0];
+ dirichletDofs[_0][i][2] = dirichletNodesZ[i][0];
}
for (size_t i = 0; i < compositeBasis.size({1}); i++) {
for (int j = 0; j < dimRotationTangent; j++){
- dirichletDofs[_1][i][j] = (not surfaceShellNodes[i][0] or dirichletNodes[i][0]);
+ dirichletDofs[_1][i][j] = not surfaceShellNodes[i][0];
}
}
--
GitLab