From 9f9c95744c0717c4df253c0d1201d7ca7326092a Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Wed, 11 Jun 2008 14:23:02 +0000
Subject: [PATCH] use the BoundaryPatch iterator
[[Imported from SVN: r2421]]
---
src/averageinterface.hh | 61 +++++++++++++----------------------------
1 file changed, 19 insertions(+), 42 deletions(-)
diff --git a/src/averageinterface.hh b/src/averageinterface.hh
index 36e615d4..dd68ba54 100644
--- a/src/averageinterface.hh
+++ b/src/averageinterface.hh
@@ -449,23 +449,15 @@ void computeAveragePressure(const Dune::FieldVector<double,GridType::dimension>&
Dune::BlockVector<Dune::FieldVector<double,1> > nodalWeights(interface.numVertices());
nodalWeights = 0;
- typename GridType::template Codim<0>::LevelIterator eIt = indexSet.template begin<0,Dune::All_Partition>();
- typename GridType::template Codim<0>::LevelIterator eEndIt = indexSet.template end<0,Dune::All_Partition>();
+ typename BoundaryPatch<GridType>::iterator it = interface.begin();
+ typename BoundaryPatch<GridType>::iterator endIt = interface.end();
- for (; eIt!=eEndIt; ++eIt) {
-
- typename GridType::template Codim<0>::Entity::LevelIntersectionIterator nIt = eIt->ilevelbegin();
- typename GridType::template Codim<0>::Entity::LevelIntersectionIterator nEndIt = eIt->ilevelend();
-
- for (; nIt!=nEndIt; ++nIt) {
-
- if (!interface.contains(*eIt,nIt))
- continue;
+ for (; it!=endIt; ++it) {
const Dune::LagrangeShapeFunctionSet<ctype, field_type, dim-1>& baseSet
- = Dune::LagrangeShapeFunctions<ctype, field_type, dim-1>::general(nIt->intersectionGlobal().type(),1);
+ = Dune::LagrangeShapeFunctions<ctype, field_type, dim-1>::general(it->intersectionGlobal().type(),1);
- const Dune::ReferenceElement<double,dim>& refElement = Dune::ReferenceElements<double, dim>::general(eIt->type());
+ const Dune::ReferenceElement<double,dim>& refElement = Dune::ReferenceElements<double, dim>::general(it->inside()->type());
// four rows because a face may have no more than four vertices
Dune::FieldVector<double,4> mu(0);
@@ -475,23 +467,23 @@ void computeAveragePressure(const Dune::FieldVector<double,GridType::dimension>&
for (int j=0; j<3; j++)
mu_tilde[i][j] = 0;
- for (int i=0; i<nIt->intersectionGlobal().corners(); i++) {
+ for (int i=0; i<it->intersectionGlobal().corners(); i++) {
const Dune::QuadratureRule<double, dim-1>& quad
- = Dune::QuadratureRules<double, dim-1>::rule(nIt->intersectionGlobal().type(), dim-1);
+ = Dune::QuadratureRules<double, dim-1>::rule(it->intersectionGlobal().type(), dim-1);
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 double integrationElement = nIt->intersectionGlobal().integrationElement(quadPos);
+ const double integrationElement = it->intersectionGlobal().integrationElement(quadPos);
// \mu_i = \int_t \varphi_i \ds
mu[i] += quad[qp].weight() * integrationElement * baseSet[i].evaluateFunction(0,quadPos);
// \tilde{\mu}_i^j = \int_t \varphi_i \times (x - x_0) \ds
- Dune::FieldVector<double,dim> worldPos = nIt->intersectionGlobal().global(quadPos);
+ Dune::FieldVector<double,dim> worldPos = it->intersectionGlobal().global(quadPos);
for (int j=0; j<dim; j++) {
@@ -511,8 +503,8 @@ void computeAveragePressure(const Dune::FieldVector<double,GridType::dimension>&
// Set up matrix
for (int i=0; i<baseSet.size(); i++) {
- int faceIdxi = refElement.subEntity(nIt->numberInSelf(), 1, i, dim);
- int subIndex = globalToLocal[indexSet.template subIndex<dim>(*eIt, faceIdxi)];
+ int faceIdxi = refElement.subEntity(it->numberInSelf(), 1, i, dim);
+ int subIndex = globalToLocal[indexSet.template subIndex<dim>(*it->inside(), faceIdxi)];
nodalWeights[subIndex] += mu[i];
for (int j=0; j<dim; j++)
@@ -524,8 +516,6 @@ void computeAveragePressure(const Dune::FieldVector<double,GridType::dimension>&
}
- }
-
}
//printmatrix(std::cout, constraints, "jacobian", "--");
@@ -588,41 +578,30 @@ void computeAveragePressure(const Dune::FieldVector<double,GridType::dimension>&
#if 1
Dune::FieldVector<double,3> outputForce(0), outputTorque(0);
- eIt = indexSet.template begin<0,Dune::All_Partition>();
- eEndIt = indexSet.template end<0,Dune::All_Partition>();
-
- for (; eIt!=eEndIt; ++eIt) {
-
- typename GridType::template Codim<0>::Entity::LevelIntersectionIterator nIt = eIt->ilevelbegin();
- typename GridType::template Codim<0>::Entity::LevelIntersectionIterator nEndIt = eIt->ilevelend();
-
- for (; nIt!=nEndIt; ++nIt) {
-
- if (!interface.contains(*eIt,nIt))
- continue;
+ for (it=interface.begin(); it!=endIt; ++it) {
const Dune::LagrangeShapeFunctionSet<double, double, dim-1>& baseSet
- = Dune::LagrangeShapeFunctions<double, double, dim-1>::general(nIt->intersectionGlobal().type(),1);
+ = Dune::LagrangeShapeFunctions<double, double, dim-1>::general(it->intersectionGlobal().type(),1);
const Dune::QuadratureRule<double, dim-1>& quad
- = Dune::QuadratureRules<double, dim-1>::rule(nIt->intersectionGlobal().type(), dim-1);
+ = Dune::QuadratureRules<double, dim-1>::rule(it->intersectionGlobal().type(), dim-1);
- const Dune::ReferenceElement<double,dim>& refElement = Dune::ReferenceElements<double, dim>::general(eIt->type());
+ const Dune::ReferenceElement<double,dim>& refElement = Dune::ReferenceElements<double, dim>::general(it->inside()->type());
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 double integrationElement = nIt->intersectionGlobal().integrationElement(quadPos);
+ const double integrationElement = it->intersectionGlobal().integrationElement(quadPos);
// Evaluate function
Dune::FieldVector<double,dim> localPressure(0);
for (size_t i=0; i<baseSet.size(); i++) {
- int faceIdxi = refElement.subEntity(nIt->numberInSelf(), 1, i, dim);
- int subIndex = indexSet.template subIndex<dim>(*eIt, faceIdxi);
+ int faceIdxi = refElement.subEntity(it->numberInSelf(), 1, i, dim);
+ int subIndex = indexSet.template subIndex<dim>(*it->inside(), faceIdxi);
localPressure.axpy(baseSet[i].evaluateFunction(0,quadPos),
pressure[subIndex]);
@@ -633,14 +612,12 @@ void computeAveragePressure(const Dune::FieldVector<double,GridType::dimension>&
outputForce.axpy(quad[qp].weight()*integrationElement, localPressure);
// Sum up the total torque \int (x - x_0) \times f dx
- Dune::FieldVector<double,dim> worldPos = nIt->intersectionGlobal().global(quadPos);
+ Dune::FieldVector<double,dim> worldPos = it->intersectionGlobal().global(quadPos);
outputTorque.axpy(quad[qp].weight()*integrationElement,
crossProduct(worldPos - crossSection.r, localPressure));
}
- }
-
}
outputForce -= resultantForce;
--
GitLab