From 122b19650249e0c8e8da9c1fb8cacdf24e8525da Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Wed, 6 Apr 2011 15:14:13 +0000
Subject: [PATCH] bugfix in computeAveragePressure: use the proper embedding of
intersections into elements when integrating over the boundary
[[Imported from SVN: r7105]]
---
dune/gfe/averageinterface.hh | 19 ++++++++++---------
1 file changed, 10 insertions(+), 9 deletions(-)
diff --git a/dune/gfe/averageinterface.hh b/dune/gfe/averageinterface.hh
index 0cef863b..bc05223c 100644
--- a/dune/gfe/averageinterface.hh
+++ b/dune/gfe/averageinterface.hh
@@ -519,7 +519,7 @@ void computeAveragePressure(const typename RigidBodyMotion<GridView::dimension>:
typedef typename GridView::ctype ctype;
typedef double field_type;
- Dune::PQkLocalFiniteElementCache<double,double, dim-1, 1> finiteElementCache;
+ Dune::PQkLocalFiniteElementCache<double,double, dim, 1> finiteElementCache;
typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
@@ -564,8 +564,8 @@ void computeAveragePressure(const typename RigidBodyMotion<GridView::dimension>:
for (; it!=endIt; ++it) {
// Get shape functions
- const typename Dune::PQkLocalFiniteElementCache<double,double, dim-1, 1>::FiniteElementType&
- localFiniteElement = finiteElementCache.get(it->type());
+ const typename Dune::PQkLocalFiniteElementCache<double,double, dim, 1>::FiniteElementType&
+ localFiniteElement = finiteElementCache.get(it->inside()->type());
const Dune::GenericReferenceElement<double,dim>& refElement = Dune::GenericReferenceElements<double, dim>::general(it->inside()->type());
@@ -585,24 +585,25 @@ void computeAveragePressure(const typename RigidBodyMotion<GridView::dimension>:
for (size_t qp=0; qp<quad.size(); qp++) {
// Local position of the quadrature point
- const Dune::FieldVector<double,dim-1>& quadPos = quad[qp].position();
+ const Dune::FieldVector<double,dim>& quadPos = it->geometryInInside().global(quad[qp].position());
- const double integrationElement = it->geometry().integrationElement(quadPos);
+ const double integrationElement = it->geometry().integrationElement(quad[qp].position());
std::vector<Dune::FieldVector<double,1> > shapeFunctionValues;
localFiniteElement.localBasis().evaluateFunction(quadPos, shapeFunctionValues);
// \mu_i = \int_t \varphi_i \ds
- mu[i] += quad[qp].weight() * integrationElement * shapeFunctionValues[i];
+ int indexInFace = refElement.subEntity(it->indexInInside(), 1, i, dim);
+ mu[i] += quad[qp].weight() * integrationElement * shapeFunctionValues[indexInFace];
// \tilde{\mu}_i^j = \int_t \varphi_i \times (x - x_0) \ds
- Dune::FieldVector<double,dim> worldPos = it->geometry().global(quadPos);
+ Dune::FieldVector<double,dim> worldPos = it->geometry().global(quad[qp].position());
for (int j=0; j<dim; j++) {
// Vector-valued basis function
Dune::FieldVector<double,dim> phi_i(0);
- phi_i[j] = shapeFunctionValues[i];
+ phi_i[j] = shapeFunctionValues[indexInFace];
mu_tilde[i][j].axpy(quad[qp].weight() * integrationElement,
crossProduct(Dune::FieldVector<double,dim>(worldPos-centerOfTorque), phi_i));
@@ -614,7 +615,7 @@ void computeAveragePressure(const typename RigidBodyMotion<GridView::dimension>:
}
// Set up matrix
- for (int i=0; i<localFiniteElement.localCoefficients().size(); i++) {
+ for (int i=0; i<refElement.size(it->indexInInside(),1,dim); i++) {
int faceIdxi = refElement.subEntity(it->indexInInside(), 1, i, dim);
int subIndex = globalToLocal[indexSet.subIndex(*it->inside(), faceIdxi, dim)];
--
GitLab