Skip to content
Snippets Groups Projects
Commit 208f032e authored by Oliver Sander's avatar Oliver Sander Committed by sander@PCPOOL.MI.FU-BERLIN.DE
Browse files

remove all code that uses lapackcpp 

[[Imported from SVN: r2362]]
parent 370770ae
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 0 additions and 246 deletions
src/averageinterface.hh
+ 0
− 246
View file @ 208f032e
...@@ -9,11 +9,6 @@ ...@@ -9,11 +9,6 @@
#include <dune/ag-common/surfmassmatrix.hh> #include <dune/ag-common/surfmassmatrix.hh>
#include "svd.hh" #include "svd.hh"
#ifdef HAVE_LAPACKPP
#include "lapackpp.h"
#undef max
#endif
template <class GridType> template <class GridType>
class PressureAverager : public Ipopt::TNLP class PressureAverager : public Ipopt::TNLP
{ {
...@@ -658,247 +653,6 @@ void computeAveragePressureIPOpt(const Dune::FieldVector<double,GridType::dimens ...@@ -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> template <class GridType>
void computeAverageInterface(const BoundaryPatch<GridType>& interface, void computeAverageInterface(const BoundaryPatch<GridType>& interface,
const Dune::BlockVector<Dune::FieldVector<double,GridType::dimension> > deformation, const Dune::BlockVector<Dune::FieldVector<double,GridType::dimension> > deformation,
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment