From 53629bf4424e41297fc3c6b33dca42834b136224 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Fri, 21 Mar 2014 21:49:25 +0000
Subject: [PATCH] [cleanup] Use a BoundaryPatch for the Dirichlet boundary
[[Imported from SVN: r9655]]
---
cosserat-continuum.cc | 51 +++++++++++++++++++------------------------
1 file changed, 22 insertions(+), 29 deletions(-)
diff --git a/cosserat-continuum.cc b/cosserat-continuum.cc
index 000dead2..db6f18c1 100644
--- a/cosserat-continuum.cc
+++ b/cosserat-continuum.cc
@@ -170,13 +170,11 @@ int main (int argc, char *argv[]) try
typedef P1NodalBasis<GridView,double> FEBasis;
FEBasis feBasis(gridView);
- SolutionType x(feBasis.size());
-
// /////////////////////////////////////////
// Read Dirichlet values
// /////////////////////////////////////////
- BitSetVector<blocksize> dirichletNodes(feBasis.size(), false);
+ BitSetVector<1> dirichletVertices(feBasis.size(), false);
BitSetVector<1> neumannNodes(feBasis.size(), false);
GridType::Codim<dim>::LeafIterator vIt = gridView.begin<dim>();
@@ -187,51 +185,46 @@ int main (int argc, char *argv[]) try
for (; vIt!=vEndIt; ++vIt) {
#if 0 // Boundary conditions for the cantilever example
- if (vIt->geometry().corner(0)[0] < 1.0+1e-3 /* or vIt->geometry().corner(0)[0] > upper[0]-1e-3*/ ) {
- // Only translation dofs are Dirichlet
- for (int j=0; j<3; j++)
- dirichletNodes[indexSet.index(*vIt)][j] = true;
- }
+ if (vIt->geometry().corner(0)[0] < 1.0+1e-3 /* or vIt->geometry().corner(0)[0] > upper[0]-1e-3*/ )
+ dirichletVertices[indexSet.index(*vIt)] = true;
if (vIt->geometry().corner(0)[0] > upper[0]-1e-3 )
- neumannNodes[indexSet.index(*vIt)][0] = true;
+ neumannNodes[indexSet.index(*vIt)] = true;
#endif
#if 1 // Boundary conditions for the shearing/wrinkling example
- if (vIt->geometry().corner(0)[1] < 1e-4 or vIt->geometry().corner(0)[1] > upper[1]-1e-4 ) {
- // Only translation dofs are Dirichlet
- for (int j=0; j<5; j++)
- dirichletNodes[indexSet.index(*vIt)][j] = true;
- }
+ if (vIt->geometry().corner(0)[1] < 1e-4 or vIt->geometry().corner(0)[1] > upper[1]-1e-4 )
+ dirichletVertices[indexSet.index(*vIt)] = true;
#endif
#if 0 // Boundary conditions for the twisted-strip example
- if (vIt->geometry().corner(0)[0] < lower[0]+1e-3 or vIt->geometry().corner(0)[0] > upper[0]-1e-3 ) {
- // Only translation dofs are Dirichlet
- for (int j=0; j<3; j++)
- dirichletNodes[indexSet.index(*vIt)][j] = true;
- }
+ if (vIt->geometry().corner(0)[0] < lower[0]+1e-3 or vIt->geometry().corner(0)[0] > upper[0]-1e-3 )
+ dirichletVertices[indexSet.index(*vIt)] = true;
#endif
#if 0 // Boundary conditions for the L-shape example
- if (vIt->geometry().corner(0)[0] < 1.0) {
- // Only translation dofs are Dirichlet
- for (int j=0; j<3; j++)
- dirichletNodes[indexSet.index(*vIt)][j] = true;
- }
+ if (vIt->geometry().corner(0)[0] < 1.0)
+ dirichletVertices[indexSet.index(*vIt)] = true;
if (vIt->geometry().corner(0)[1] < -239 )
neumannNodes[indexSet.index(*vIt)][0] = true;
#endif
}
- //////////////////////////////////////////////////////////////////////////////
- // Assemble Neumann term
- //////////////////////////////////////////////////////////////////////////////
-
+ BoundaryPatch<GridView> dirichletBoundary(gridView, dirichletVertices);
BoundaryPatch<GridView> neumannBoundary(gridView, neumannNodes);
std::cout << "Neumann boundary has " << neumannBoundary.numFaces() << " faces\n";
+
+
+ BitSetVector<blocksize> dirichletNodes(feBasis.size(), false);
+ for (size_t i=0; i<feBasis.size(); i++)
+ if (dirichletVertices[i][0])
+ for (int j=0; j<5; j++)
+ dirichletNodes[i][j] = true;
+
// //////////////////////////
- // Initial solution
+ // Initial iterate
// //////////////////////////
+ SolutionType x(feBasis.size());
+
vIt = gridView.begin<dim>();
for (; vIt!=vEndIt; ++vIt) {
--
GitLab