From 873fa08b3c943c10d2d713a1e46b06836a88b175 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Wed, 22 Aug 2007 16:11:10 +0000
Subject: [PATCH] make this work for interface which consist of more than a
single segment
[[Imported from SVN: r1563]]
---
src/averageinterface.hh | 18 +++++++++++++-----
1 file changed, 13 insertions(+), 5 deletions(-)
diff --git a/src/averageinterface.hh b/src/averageinterface.hh
index b5db8571..36d2fe20 100644
--- a/src/averageinterface.hh
+++ b/src/averageinterface.hh
@@ -23,9 +23,13 @@ void computeAveragePressure(const Dune::FieldVector<double,GridType::dimension>&
const typename GridType::Traits::LevelIndexSet& indexSet = grid.levelIndexSet(level);
const int dim = GridType::dimension;
typedef typename GridType::ctype ctype;
+ typedef double field_type;
typedef typename GridType::template Codim<dim>::LevelIterator VertexIterator;
+ // Get total interface area
+ ctype area = interface.area();
+
// set up output array
DGIndexSet<GridType> dgIndexSet(grid,level);
dgIndexSet.setup(grid,level);
@@ -45,8 +49,8 @@ void computeAveragePressure(const Dune::FieldVector<double,GridType::dimension>&
if (!interface.contains(*eIt,nIt))
continue;
- const Dune::LagrangeShapeFunctionSet<double, double, dim-1>& baseSet
- = Dune::LagrangeShapeFunctions<double, double, dim-1>::general(nIt.intersectionGlobal().type(),1);
+ const Dune::LagrangeShapeFunctionSet<ctype, field_type, dim-1>& baseSet
+ = Dune::LagrangeShapeFunctions<ctype, field_type, dim-1>::general(nIt.intersectionGlobal().type(),1);
if (baseSet.size() != 4)
DUNE_THROW(Dune::NotImplemented, "average interface only for quad faces");
@@ -113,9 +117,13 @@ void computeAveragePressure(const Dune::FieldVector<double,GridType::dimension>&
Dune::FieldVector<double,12> u;
Dune::FieldVector<double,6> b;
+ // Scale the resultant force and torque with this segments area percentage.
+ // That way the resulting pressure gets distributed fairly uniformly.
+ ctype segmentArea = nIt.intersectionGlobal().volume() / area;
+
for (int i=0; i<3; i++) {
- b[i] = resultantForce[i];
- b[i+3] = resultantTorque[i];
+ b[i] = resultantForce[i] * segmentArea;
+ b[i+3] = resultantTorque[i] * segmentArea;
}
// std::cout << b << std::endl;
@@ -333,7 +341,7 @@ void computeAverageInterface(const BoundaryPatch<GridType>& interface,
// average deformation is the integral over the deformation gradient
// divided by its area
deformationGradient /= interfaceArea;
- std::cout << "deformationGradient: " << std::endl << deformationGradient << std::endl;
+ //std::cout << "deformationGradient: " << std::endl << deformationGradient << std::endl;
// Get the rotational part of the deformation gradient by performing a
// polar composition.
--
GitLab