From 4021b000c59b304a1fae3d1ef5c55f4439a9f847 Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Wed, 17 Oct 2012 14:53:19 +0200
Subject: [PATCH] Only write out data for a central node
---
src/one-body-sample.cc | 54 +++++++++++++++++++++++++-----------------
1 file changed, 32 insertions(+), 22 deletions(-)
diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc
index 3bb04f71..a10414ec 100644
--- a/src/one-body-sample.cc
+++ b/src/one-body-sample.cc
@@ -86,7 +86,8 @@ void setup_boundary(GridView const &gridView,
Dune::BitSetVector<dim> &ignoreNodes,
Dune::BitSetVector<1> &neumannNodes,
Dune::BitSetVector<1> &frictionalNodes,
- GridCorner const &lowerLeft, GridCorner const &upperRight) {
+ GridCorner const &lowerLeft, GridCorner const &upperRight,
+ size_t &specialNode) {
typedef typename GridView::template Codim<dim>::Iterator VertexLeafIterator;
Dune::MultipleCodimMultipleGeomTypeMapper<
@@ -97,6 +98,22 @@ void setup_boundary(GridView const &gridView,
assert(it->geometry().corners() == 1);
Dune::FieldVector<double, dim> const coordinates = it->geometry().corner(0);
size_t const id = myVertexMapper.map(*it);
+
+ // Find the center of the lower face
+ switch (dim) {
+ case 3:
+ if (coordinates[2] != lowerLeft[2])
+ break;
+ // fallthrough
+ case 2:
+ if (coordinates[1] == lowerLeft[1] &&
+ std::abs(coordinates[0] - lowerLeft[0]) < 1e-8)
+ specialNode = id;
+ break;
+ default:
+ assert(false);
+ }
+
if (coordinates[1] == upperRight[1]) // upper face
ignoreNodes[id] = true;
else if (coordinates[1] == lowerLeft[1]) { // lower face
@@ -234,11 +251,14 @@ int main(int argc, char *argv[]) {
EnergyNorm<MatrixType, VectorType> energyNorm(stiffnessMatrix);
// Set up the boundary
+ size_t specialNode = finestSize;
Dune::BitSetVector<dim> ignoreNodes(finestSize, false);
Dune::BitSetVector<1> neumannNodes(finestSize, false);
Dune::BitSetVector<1> frictionalNodes(finestSize, false);
setup_boundary(leafView, ignoreNodes, neumannNodes, frictionalNodes,
- lowerLeft, upperRight);
+ lowerLeft, upperRight, specialNode);
+ assert(specialNode != finestSize);
+ assert(frictionalNodes[specialNode][0]);
auto const nodalIntegrals =
assemble_frictional<GridType, GridView, SmallVector, P1Basis>(
@@ -425,26 +445,16 @@ int main(int argc, char *argv[]) {
alpha_old = alpha;
{ // Write all kinds of data
- for (size_t i = 0; i < frictionalNodes.size(); ++i)
- if (frictionalNodes[i][0])
- state_writer << alpha[i][0] << " ";
- state_writer << std::endl;
-
- for (size_t i = 0; i < frictionalNodes.size(); ++i)
- if (frictionalNodes[i][0])
- displacement_writer << u[i][0] << " ";
- displacement_writer << std::endl;
-
- for (size_t i = 0; i < frictionalNodes.size(); ++i)
- if (frictionalNodes[i][0])
- velocity_writer << ud[i][0] << " ";
- velocity_writer << std::endl;
-
- for (size_t i = 0; i < frictionalNodes.size(); ++i)
- if (frictionalNodes[i][0])
- coefficient_writer << mu + a *std::log(ud[i].two_norm() / V0) +
- b * (alpha[i] + std::log(V0 / L)) << " ";
- coefficient_writer << std::endl;
+ state_writer << alpha[specialNode][0] << " " << std::endl;
+
+ displacement_writer << u[specialNode][0] << " " << std::endl;
+
+ velocity_writer << ud[specialNode][0] << " " << std::endl;
+
+ coefficient_writer << mu +
+ a *std::log(ud[specialNode].two_norm() / V0) +
+ b * (alpha[specialNode] + std::log(V0 / L))
+ << " " << std::endl;
}
// Compute von Mises stress and write everything to a file
--
GitLab