From 208f032e99b047eec07214e8c2c095989d40dd9f Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Thu, 5 Jun 2008 15:50:46 +0000
Subject: [PATCH] remove all code that uses lapackcpp
[[Imported from SVN: r2362]]
---
src/averageinterface.hh | 246 ----------------------------------------
1 file changed, 246 deletions(-)
diff --git a/src/averageinterface.hh b/src/averageinterface.hh
index f8beb972..c8a80609 100644
--- a/src/averageinterface.hh
+++ b/src/averageinterface.hh
@@ -9,11 +9,6 @@
#include <dune/ag-common/surfmassmatrix.hh>
#include "svd.hh"
-#ifdef HAVE_LAPACKPP
-#include "lapackpp.h"
-#undef max
-#endif
-
template <class GridType>
class PressureAverager : public Ipopt::TNLP
{
@@ -658,247 +653,6 @@ void computeAveragePressureIPOpt(const Dune::FieldVector<double,GridType::dimens
}
-
-#ifdef HAVE_LAPACKPP
-// Given a resultant force and torque (from a rod problem), this method computes the corresponding
-// Neumann data for a 3d elasticity problem.
-template <class GridType>
-void computeAveragePressure(const Dune::FieldVector<double,GridType::dimension>& resultantForce,
- const Dune::FieldVector<double,GridType::dimension>& resultantTorque,
- const BoundaryPatch<GridType>& interface,
- const Configuration& crossSection,
- Dune::BlockVector<Dune::FieldVector<double, GridType::dimension> >& pressure)
-{
- const GridType& grid = interface.getGrid();
- const int level = interface.level();
- 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);
- pressure.resize(dgIndexSet.size());
- pressure = 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>();
-
- 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;
-
- const Dune::LagrangeShapeFunctionSet<ctype, field_type, dim-1>& baseSet
- = Dune::LagrangeShapeFunctions<ctype, field_type, dim-1>::general(nIt->intersectionGlobal().type(),1);
-
- // four rows because a face may have no more than four vertices
- Dune::FieldVector<double,4> mu(0);
- Dune::FieldVector<double,3> mu_tilde[4][3];
-
- for (int i=0; i<4; i++)
- for (int j=0; j<3; j++)
- mu_tilde[i][j] = 0;
-
- for (int i=0; i<nIt->intersectionGlobal().corners(); i++) {
-
- const Dune::QuadratureRule<double, dim-1>& quad
- = Dune::QuadratureRules<double, dim-1>::rule(nIt->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);
-
- // \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);
-
- for (int j=0; j<dim; j++) {
-
- // Vector-valued basis function
- Dune::FieldVector<double,dim> phi_i(0);
- phi_i[j] = baseSet[i].evaluateFunction(0,quadPos);
-
- mu_tilde[i][j].axpy(quad[qp].weight() * integrationElement,
- crossProduct(worldPos-crossSection.r, phi_i));
-
- }
-
- }
-
- }
-
- // Set up matrix
- // LaPack++ style
- LaGenMatDouble matrix(6, 3*baseSet.size());
- matrix = 0;
- for (int i=0; i<baseSet.size(); i++)
- for (int j=0; j<3; j++)
- matrix(j, i*3+j) = mu[i];
-
- for (int i=0; i<baseSet.size(); i++)
- for (int j=0; j<3; j++)
- for (int k=0; k<3; k++)
- matrix(3+k, 3*i+j) = mu_tilde[i][j][k];
-
- LaVectorDouble u(3*baseSet.size());
- LaVectorDouble b(6);
-
- // 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] * segmentArea;
- b(i+3) = resultantTorque[i] * segmentArea;
- }
-
- LaLinearSolve(matrix, u, b);
-
- for (int i=0; i<baseSet.size(); i++)
- for (int j=0; j<3; j++)
- pressure[dgIndexSet(*eIt, nIt->numberInSelf())+i][j] = u(i*3+j);
-
- }
-
- }
-
-
- // /////////////////////////////////////////////////////////////////////////////////////
- // Compute the overall force and torque to see whether the preceding code is correct
- // /////////////////////////////////////////////////////////////////////////////////////
-#ifdef NDEBUG
- 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;
-
- const Dune::LagrangeShapeFunctionSet<double, double, dim-1>& baseSet
- = Dune::LagrangeShapeFunctions<double, double, dim-1>::general(nIt->intersectionGlobal().type(),1);
-
- const Dune::QuadratureRule<double, dim-1>& quad
- = Dune::QuadratureRules<double, dim-1>::rule(nIt->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);
-
- // Evaluate function
- Dune::FieldVector<double,dim> localPressure(0);
-
- for (size_t i=0; i<baseSet.size(); i++)
- localPressure.axpy(baseSet[i].evaluateFunction(0,quadPos),
- pressure[dgIndexSet(*eIt,nIt->numberInSelf())+i]);
-
-
- // Sum up the total force
- 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);
- outputTorque.axpy(quad[qp].weight()*integrationElement,
- crossProduct(worldPos - crossSection.r, localPressure));
-
- }
-
- }
-
- }
-
- outputForce -= resultantForce;
- outputTorque -= resultantTorque;
- assert( outputForce.infinity_norm() < 1e-6 );
- assert( outputTorque.infinity_norm() < 1e-6 );
-// std::cout << "Output force: " << outputForce << std::endl;
-// std::cout << "Output torque: " << outputTorque << " " << resultantTorque[0]/outputTorque[0] << std::endl;
-#endif
-
-}
-#endif
-
-
-template <class GridType>
-void averageSurfaceDGFunction(const GridType& grid,
- const Dune::BlockVector<Dune::FieldVector<double,GridType::dimension> >& dgFunction,
- Dune::BlockVector<Dune::FieldVector<double,GridType::dimension> >& p1Function,
- const DGIndexSet<GridType>& dgIndexSet)
-{
- const int dim = GridType::dimension;
-
- const typename GridType::Traits::LeafIndexSet& indexSet = grid.leafIndexSet();
-
- p1Function.resize(indexSet.size(dim));
- p1Function = 0;
-
- std::vector<int> counter(indexSet.size(dim), 0);
-
- typename GridType::template Codim<0>::LeafIterator eIt = grid.template leafbegin<0>();
- typename GridType::template Codim<0>::LeafIterator eEndIt = grid.template leafend<0>();
-
- for (; eIt!=eEndIt; ++eIt) {
-
- typename GridType::template Codim<0>::Entity::LeafIntersectionIterator nIt = eIt->ileafbegin();
- typename GridType::template Codim<0>::Entity::LeafIntersectionIterator nEndIt = eIt->ileafend();
-
- for (; nIt!=nEndIt; ++nIt) {
-
- if (!nIt->boundary())
- continue;
-
- const Dune::ReferenceElement<double,dim>& refElement
- = Dune::ReferenceElements<double, dim>::general(eIt->type());
-
- for (int i=0; i<refElement.size(nIt->numberInSelf(),1,dim); i++) {
-
- int idxInElement = refElement.subEntity(nIt->numberInSelf(),1, i, dim);
-
- p1Function[indexSet.template subIndex<dim>(*eIt,idxInElement)]
- += dgFunction[dgIndexSet(*eIt,nIt->numberInSelf())+i];
-
- counter[indexSet.template subIndex<dim>(*eIt,idxInElement)]++;
-
- }
-
- }
-
- }
-
- for (int i=0; i<p1Function.size(); i++)
- if (counter[i]!=0)
- p1Function[i] /= counter[i];
-
-}
-
-
template <class GridType>
void computeAverageInterface(const BoundaryPatch<GridType>& interface,
const Dune::BlockVector<Dune::FieldVector<double,GridType::dimension> > deformation,
--
GitLab