Fix neumann boundary in y z

src/film-on-substrate.cc

@@ -253,7 +253,9 @@ int main (int argc, char *argv[]) try
   //                      BOUNDARY DATA
   /////////////////////////////////////////////////////////////
 
-  BitSetVector<dim> 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);
 
@@ -262,8 +264,9 @@ int main (int argc, char *argv[]) try
   for (auto&& v : vertices(gridView))
   {
     FieldVector<bool,dim> isDirichlet = pythonDirichletVertices(v.geometry().corner(0));
-    for (int i = 0; i )
-    dirichletVertices[indexSet.index(v)] = isDirichlet;
+    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;
@@ -272,19 +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<dim> dirichletNodes(compositeBasis.size({0}),false);
-  BitSetVector<1> neumannNodes(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,deformationPowerBasis,dirichletNodes);
-  constructBoundaryDofs(*neumannBoundary,deformationFEBasis,neumannNodes);
+  constructBoundaryDofs(dirichletBoundaryX,deformationFEBasis,dirichletNodesX);
+  constructBoundaryDofs(dirichletBoundaryY,deformationFEBasis,dirichletNodesY);
+  constructBoundaryDofs(dirichletBoundaryZ,deformationFEBasis,dirichletNodesZ);
   constructBoundaryDofs(surfaceShellBoundary,orientationFEBasis,surfaceShellNodes);
 
   //Create BitVector matching the tangential space
@@ -295,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][j]; //dirichlet nodes and also fix the area where we can pull in y- and z-direction, as we only want movement in x-direction!
-    }
+  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][j]);
+      dirichletDofs[_1][i][j] = not surfaceShellNodes[i][0];
     }
   }