Skip to content
Snippets Groups Projects
Commit 94ac96f6 authored by Sander, Oliver's avatar Sander, Oliver
Browse files

Stop using LocalIndexSet if dune-functions is 2.7 or newer

parent bd02d57d
No related branches found
No related tags found
No related merge requests found
......@@ -174,18 +174,23 @@ public:
}
std::vector<int> connectivity(connectivitySize);
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
auto localIndexSet = basis.localIndexSet();
#endif
size_t i=0;
for (const auto& element : elements(gridView, Dune::Partitions::interior))
{
localView.bind(element);
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
localIndexSet.bind(localView);
#endif
if (element.type().isQuadrilateral())
{
if (vtkOrder==2)
{
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
connectivity[i++] = localIndexSet.index(0);
connectivity[i++] = localIndexSet.index(2);
connectivity[i++] = localIndexSet.index(8);
......@@ -195,31 +200,64 @@ public:
connectivity[i++] = localIndexSet.index(5);
connectivity[i++] = localIndexSet.index(7);
connectivity[i++] = localIndexSet.index(3);
#else
connectivity[i++] = localView.index(0);
connectivity[i++] = localView.index(2);
connectivity[i++] = localView.index(8);
connectivity[i++] = localView.index(6);
connectivity[i++] = localView.index(1);
connectivity[i++] = localView.index(5);
connectivity[i++] = localView.index(7);
connectivity[i++] = localView.index(3);
#endif
}
else // first order
{
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
connectivity[i++] = localIndexSet.index(0);
connectivity[i++] = localIndexSet.index(1);
connectivity[i++] = localIndexSet.index(3);
connectivity[i++] = localIndexSet.index(2);
#else
connectivity[i++] = localView.index(0);
connectivity[i++] = localView.index(1);
connectivity[i++] = localView.index(3);
connectivity[i++] = localView.index(2);
#endif
}
}
if (element.type().isTriangle())
{
if (vtkOrder==2)
{
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
connectivity[i++] = localIndexSet.index(0);
connectivity[i++] = localIndexSet.index(2);
connectivity[i++] = localIndexSet.index(5);
connectivity[i++] = localIndexSet.index(1);
connectivity[i++] = localIndexSet.index(4);
connectivity[i++] = localIndexSet.index(3);
#else
connectivity[i++] = localView.index(0);
connectivity[i++] = localView.index(2);
connectivity[i++] = localView.index(5);
connectivity[i++] = localView.index(1);
connectivity[i++] = localView.index(4);
connectivity[i++] = localView.index(3);
#endif
}
else // first order
{
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
connectivity[i++] = localIndexSet.index(0);
connectivity[i++] = localIndexSet.index(1);
connectivity[i++] = localIndexSet.index(2);
#else
connectivity[i++] = localView.index(0);
connectivity[i++] = localView.index(1);
connectivity[i++] = localView.index(2);
#endif
}
}
}
......@@ -375,15 +413,20 @@ public:
outFile << " <Cells>" << std::endl;
outFile << " <DataArray type=\"Int32\" Name=\"connectivity\" NumberOfComponents=\"1\" format=\"ascii\">" << std::endl;
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
auto localIndexSet = p2DeformationBasis.localIndexSet();
#endif
for (const auto& element : elements(gridView, Dune::Partitions::interior))
{
localView.bind(element);
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
localIndexSet.bind(localView);
#endif
outFile << " ";
if (element.type().isQuadrilateral())
{
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
outFile << localIndexSet.index(0) << " ";
outFile << localIndexSet.index(2) << " ";
outFile << localIndexSet.index(8) << " ";
......@@ -393,6 +436,17 @@ public:
outFile << localIndexSet.index(5) << " ";
outFile << localIndexSet.index(7) << " ";
outFile << localIndexSet.index(3) << " ";
#else
outFile << localView.index(0) << " ";
outFile << localView.index(2) << " ";
outFile << localView.index(8) << " ";
outFile << localView.index(6) << " ";
outFile << localView.index(1) << " ";
outFile << localView.index(5) << " ";
outFile << localView.index(7) << " ";
outFile << localView.index(3) << " ";
#endif
outFile << std::endl;
}
}
......
......@@ -68,17 +68,25 @@ public:
typename TargetSpace::CoordinateType operator()(const Element& element, const Dune::FieldVector<ctype,gridDim>& local) const
{
auto localView = basis_.localView();
auto localIndexSet = basis_.localIndexSet();
localView.bind(element);
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
auto localIndexSet = basis_.localIndexSet();
localIndexSet.bind(localView);
auto numOfBaseFct = localIndexSet.size();
#else
auto numOfBaseFct = localView.size();
#endif
// Extract local coefficients
std::vector<TargetSpace> localCoeff(numOfBaseFct);
for (size_t i=0; i<numOfBaseFct; i++)
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
localCoeff[i] = coefficients_[localIndexSet.index(i)];
#else
localCoeff[i] = coefficients_[localView.index(i)];
#endif
// create local gfe function
LocalInterpolationRule localInterpolationRule(localView.tree().finiteElement(),localCoeff);
......@@ -96,17 +104,26 @@ public:
Dune::FieldMatrix<ctype, embeddedDim, gridDim> derivative(const Element& element, const Dune::FieldVector<ctype,gridDim>& local) const
{
auto localView = basis_.localView();
auto localIndexSet = basis_.localIndexSet();
localView.bind(element);
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
auto localIndexSet = basis_.localIndexSet();
localIndexSet.bind(localView);
int numOfBaseFct = localIndexSet.size();
#else
auto numOfBaseFct = localView.size();
#endif
// Extract local coefficients
std::vector<TargetSpace> localCoeff(numOfBaseFct);
for (int i=0; i<numOfBaseFct; i++)
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
localCoeff[i] = coefficients_[localIndexSet.index(i)];
#else
localCoeff[i] = coefficients_[localView.index(i)];
#endif
// create local gfe function
LocalInterpolationRule localInterpolationRule(localView.tree().finiteElement(),localCoeff);
......
......@@ -76,7 +76,9 @@ getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const
// A view on the FE basis on a single element
auto localView = basis_.localView();
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
auto localIndexSet = basis_.localIndexSet();
#endif
ElementIterator it = basis_.gridView().template begin<0,Dune::Interior_Partition>();
ElementIterator endit = basis_.gridView().template end<0,Dune::Interior_Partition> ();
......@@ -85,7 +87,9 @@ getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const
// Bind the local FE basis view to the current element
localView.bind(*it);
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
localIndexSet.bind(localView);
#endif
const auto& lfe = localView.tree().finiteElement();
......@@ -93,8 +97,13 @@ getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const
for (size_t j=0; j<lfe.size(); j++) {
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
auto iIdx = localIndexSet.index(i)[0];
auto jIdx = localIndexSet.index(j)[0];
#else
auto iIdx = localView.index(i);
auto jIdx = localView.index(j);
#endif
nb.add(iIdx, jIdx);
......@@ -128,7 +137,9 @@ assembleGradientAndHessian(const std::vector<TargetSpace>& sol,
// A view on the FE basis on a single element
auto localView = basis_.localView();
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
auto localIndexSet = basis_.localIndexSet();
#endif
ElementIterator it = basis_.gridView().template begin<0,Dune::Interior_Partition>();
ElementIterator endit = basis_.gridView().template end<0,Dune::Interior_Partition> ();
......@@ -136,7 +147,9 @@ assembleGradientAndHessian(const std::vector<TargetSpace>& sol,
for( ; it != endit; ++it ) {
localView.bind(*it);
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
localIndexSet.bind(localView);
#endif
const int numOfBaseFct = localView.tree().size();
......@@ -144,7 +157,11 @@ assembleGradientAndHessian(const std::vector<TargetSpace>& sol,
std::vector<TargetSpace> localSolution(numOfBaseFct);
for (int i=0; i<numOfBaseFct; i++)
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
localSolution[i] = sol[localIndexSet.index(i)[0]];
#else
localSolution[i] = sol[localView.index(i)];
#endif
std::vector<Dune::FieldVector<double,blocksize> > localGradient(numOfBaseFct);
......@@ -154,11 +171,19 @@ assembleGradientAndHessian(const std::vector<TargetSpace>& sol,
// Add element matrix to global stiffness matrix
for(int i=0; i<numOfBaseFct; i++) {
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
auto row = localIndexSet.index(i)[0];
#else
auto row = localView.index(i);
#endif
for (int j=0; j<numOfBaseFct; j++ ) {
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
auto col = localIndexSet.index(j)[0];
#else
auto col = localView.index(j);
#endif
hessian[row][col] += localStiffness_->A_[i][j];
}
......@@ -166,7 +191,11 @@ assembleGradientAndHessian(const std::vector<TargetSpace>& sol,
// Add local gradient to global gradient
for (int i=0; i<numOfBaseFct; i++)
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
gradient[localIndexSet.index(i)[0]] += localGradient[i];
#else
gradient[localView.index(i)] += localGradient[i];
#endif
}
......@@ -185,7 +214,9 @@ assembleGradient(const std::vector<TargetSpace>& sol,
// A view on the FE basis on a single element
auto localView = basis_.localView();
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
auto localIndexSet = basis_.localIndexSet();
#endif
ElementIterator it = basis_.gridView().template begin<0,Dune::Interior_Partition>();
ElementIterator endIt = basis_.gridView().template end<0,Dune::Interior_Partition>();
......@@ -194,7 +225,9 @@ assembleGradient(const std::vector<TargetSpace>& sol,
for (; it!=endIt; ++it) {
localView.bind(*it);
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
localIndexSet.bind(localView);
#endif
// A 1d grid has two vertices
const auto nDofs = localView.tree().size();
......@@ -203,7 +236,11 @@ assembleGradient(const std::vector<TargetSpace>& sol,
std::vector<TargetSpace> localSolution(nDofs);
for (size_t i=0; i<nDofs; i++)
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
localSolution[i] = sol[localIndexSet.index(i)[0]];
#else
localSolution[i] = sol[localView.index(i)];
#endif
// Assemble local gradient
std::vector<Dune::FieldVector<double,blocksize> > localGradient(nDofs);
......@@ -212,8 +249,11 @@ assembleGradient(const std::vector<TargetSpace>& sol,
// Add to global gradient
for (size_t i=0; i<nDofs; i++)
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
grad[localIndexSet.index(i)[0]] += localGradient[i];
#else
grad[localView.index(i)[0]] += localGradient[i];
#endif
}
}
......@@ -230,7 +270,9 @@ computeEnergy(const std::vector<TargetSpace>& sol) const
// A view on the FE basis on a single element
auto localView = basis_.localView();
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
auto localIndexSet = basis_.localIndexSet();
#endif
ElementIterator it = basis_.gridView().template begin<0,Dune::Interior_Partition>();
ElementIterator endIt = basis_.gridView().template end<0,Dune::Interior_Partition>();
......@@ -239,7 +281,9 @@ computeEnergy(const std::vector<TargetSpace>& sol) const
for (; it!=endIt; ++it) {
localView.bind(*it);
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
localIndexSet.bind(localView);
#endif
// Number of degrees of freedom on this element
size_t nDofs = localView.tree().size();
......@@ -247,7 +291,11 @@ computeEnergy(const std::vector<TargetSpace>& sol) const
std::vector<TargetSpace> localSolution(nDofs);
for (size_t i=0; i<nDofs; i++)
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
localSolution[i] = sol[localIndexSet.index(i)[0]];
#else
localSolution[i] = sol[localView.index(i)[0]];
#endif
energy += localStiffness_->energy(localView, localSolution);
......
......@@ -95,13 +95,16 @@ getMatrixPattern(Dune::MatrixIndexSet& nb00,
// A view on the FE basis on a single element
auto localView = basis_.localView();
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
auto localIndexSet = basis_.localIndexSet();
#endif
// Loop over grid elements
for (const auto& element : elements(basis_.gridView(), Dune::Partitions::interior))
{
// Bind the local FE basis view to the current element
localView.bind(element);
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
localIndexSet.bind(localView);
// Add element stiffness matrix onto the global stiffness matrix
......@@ -114,6 +117,18 @@ getMatrixPattern(Dune::MatrixIndexSet& nb00,
{
// The global index of the j-th local degree of freedom of the element 'e'
auto col = localIndexSet.index(j);
#else
// Add element stiffness matrix onto the global stiffness matrix
for (size_t i=0; i<localView.size(); i++)
{
// The global index of the i-th local degree of freedom of the element 'e'
auto row = localView.index(i);
for (size_t j=0; j<localView.size(); j++ )
{
// The global index of the j-th local degree of freedom of the element 'e'
auto col = localView.index(j);
#endif
if (row[0]==0 and col[0]==0)
nb00.add(row[1],col[1]);
......@@ -167,13 +182,17 @@ assembleGradientAndHessian(const std::vector<TargetSpace0>& configuration0,
// A view on the FE basis on a single element
auto localView = basis_.localView();
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
auto localIndexSet = basis_.localIndexSet();
#endif
for (const auto& element : elements(basis_.gridView(), Dune::Partitions::interior))
{
// Bind the local FE basis view to the current element
localView.bind(element);
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
localIndexSet.bind(localView);
#endif
using namespace Dune::TypeTree::Indices;
const int nDofs0 = localView.tree().child(_0).finiteElement().size();
......@@ -185,10 +204,17 @@ assembleGradientAndHessian(const std::vector<TargetSpace0>& configuration0,
for (int i=0; i<nDofs0+nDofs1; i++)
{
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
if (localIndexSet.index(i)[0] == 0)
localConfiguration0[i] = configuration0[localIndexSet.index(i)[1]];
else
localConfiguration1[i-nDofs0] = configuration1[localIndexSet.index(i)[1]];
#else
if (localView.index(i)[0] == 0)
localConfiguration0[i] = configuration0[localView.index(i)[1]];
else
localConfiguration1[i-nDofs0] = configuration1[localView.index(i)[1]];
#endif
}
std::vector<Dune::FieldVector<double,blocksize0> > localGradient0(nDofs0);
......@@ -202,11 +228,19 @@ assembleGradientAndHessian(const std::vector<TargetSpace0>& configuration0,
// Add element matrix to global stiffness matrix
for (int i=0; i<nDofs0+nDofs1; i++)
{
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
auto row = localIndexSet.index(i);
#else
auto row = localView.index(i);
#endif
for (int j=0; j<nDofs0+nDofs1; j++ )
{
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
auto col = localIndexSet.index(j);
#else
auto col = localView.index(j);
#endif
if (row[0]==0 and col[0]==0)
hessian[_0][_0][row[1]][col[1]] += localStiffness_->A00_[i][j];
......@@ -222,10 +256,17 @@ assembleGradientAndHessian(const std::vector<TargetSpace0>& configuration0,
}
// Add local gradient to global gradient
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
if (localIndexSet.index(i)[0] == 0)
gradient0[localIndexSet.index(i)[1]] += localGradient0[i];
else
gradient1[localIndexSet.index(i)[1]] += localGradient1[i-nDofs0];
#else
if (localView.index(i)[0] == 0)
gradient0[localView.index(i)[1]] += localGradient0[i];
else
gradient1[localView.index(i)[1]] += localGradient1[i-nDofs0];
#endif
}
}
......@@ -287,14 +328,18 @@ computeEnergy(const std::vector<TargetSpace0>& configuration0,
// A view on the FE basis on a single element
auto localView = basis_.localView();
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
auto localIndexSet = basis_.localIndexSet();
#endif
// Loop over all elements
for (const auto& element : elements(basis_.gridView(), Dune::Partitions::interior))
{
// Bind the local FE basis view to the current element
localView.bind(element);
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
localIndexSet.bind(localView);
#endif
// Number of degrees of freedom on this element
using namespace Dune::TypeTree::Indices;
......@@ -306,10 +351,17 @@ computeEnergy(const std::vector<TargetSpace0>& configuration0,
for (int i=0; i<nDofs0+nDofs1; i++)
{
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
if (localIndexSet.index(i)[0] == 0)
localConfiguration0[i] = configuration0[localIndexSet.index(i)[1]];
else
localConfiguration1[i-nDofs0] = configuration1[localIndexSet.index(i)[1]];
#else
if (localView.index(i)[0] == 0)
localConfiguration0[i] = configuration0[localView.index(i)[1]];
else
localConfiguration1[i-nDofs0] = configuration1[localView.index(i)[1]];
#endif
}
energy += localStiffness_->energy(localView,
......
......@@ -46,21 +46,33 @@ namespace Dune {
size_ = globalVertexIndex.size(2) + globalEdgeIndex.size(1) + globalElementIndex.size(0);
auto localView = p2Mapper_.localView();
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
auto localIndexSet = p2Mapper_.localIndexSet();
#endif
// Determine
for (const auto& element : elements(gridView))
{
localView.bind(element);
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
localIndexSet.bind(localView);
#endif
// Loop over all local degrees of freedom
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
for (size_t i=0; i<localIndexSet.size(); i++)
#else
for (size_t i=0; i<localView.size(); i++)
#endif
{
int codim = localView.tree().finiteElement().localCoefficients().localKey(i).codim();
int entity = localView.tree().finiteElement().localCoefficients().localKey(i).subEntity();
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
int localIndex = localIndexSet.index(i);
#else
auto localIndex = localView.index(i);
#endif
int globalIndex;
switch (codim)
{
......@@ -105,19 +117,29 @@ namespace Dune {
bool contains(const Entity& entity, uint subEntity, uint codim, Index& result) const
{
auto localView = p2Mapper_.localView();
auto localIndexSet = p2Mapper_.localIndexSet();
localView.bind(entity);
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
auto localIndexSet = p2Mapper_.localIndexSet();
localIndexSet.bind(localView);
#endif
Index localIndex;
bool dofFound = false;
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
for (size_t i=0; i<localIndexSet.size(); i++)
#else
for (size_t i=0; i<localView.size(); i++)
#endif
{
if (localView.tree().finiteElement().localCoefficients().localKey(i).subEntity() == subEntity
and localView.tree().finiteElement().localCoefficients().localKey(i).codim() == codim)
{
dofFound = true;
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
localIndex = localIndexSet.index(i);
#else
localIndex = localView.index(i);
#endif
break;
}
}
......
......@@ -31,16 +31,26 @@ public:
bool contains(const Entity& entity, uint subEntity, uint codim, Index& result) const
{
auto localView = p2Basis_.localView();
auto localIndexSet = p2Basis_.localIndexSet();
localView.bind(entity);
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
auto localIndexSet = p2Basis_.localIndexSet();
localIndexSet.bind(localView);
#endif
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
for (size_t i=0; i<localIndexSet.size(); i++)
#else
for (size_t i=0; i<localView.size(); i++)
#endif
{
if (localView.tree().finiteElement().localCoefficients().localKey(i).subEntity() == subEntity
and localView.tree().finiteElement().localCoefficients().localKey(i).codim() == codim)
{
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
result = localIndexSet.index(i)[0];
#else
result = localView.index(i);
#endif
return true;
}
}
......
......@@ -26,7 +26,9 @@ assembleGradient(const std::vector<RigidBodyMotion<double,3> >& sol,
// A view on the FE basis on a single element
auto localView = this->basis_.localView();
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
auto localIndexSet = this->basis_.localIndexSet();
#endif
ElementIterator it = this->basis_.gridView().template begin<0>();
ElementIterator endIt = this->basis_.gridView().template end<0>();
......@@ -35,7 +37,9 @@ assembleGradient(const std::vector<RigidBodyMotion<double,3> >& sol,
for (; it!=endIt; ++it) {
localView.bind(*it);
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
localIndexSet.bind(localView);
#endif
// A 1d grid has two vertices
static const int nDofs = 2;
......@@ -44,7 +48,11 @@ assembleGradient(const std::vector<RigidBodyMotion<double,3> >& sol,
std::vector<RigidBodyMotion<double,3> > localSolution(nDofs);
for (int i=0; i<nDofs; i++)
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
localSolution[i] = sol[localIndexSet.index(i)[0]];
#else
localSolution[i] = sol[localView.index(i)];
#endif
// Assemble local gradient
std::vector<FieldVector<double,blocksize> > localGradient(nDofs);
......@@ -55,7 +63,11 @@ assembleGradient(const std::vector<RigidBodyMotion<double,3> >& sol,
// Add to global gradient
for (int i=0; i<nDofs; i++)
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
grad[localIndexSet.index(i)[0]] += localGradient[i];
#else
grad[localView.index(i)] += localGradient[i];
#endif
}
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment