Skip to content
Snippets Groups Projects
Commit de9b228e authored by Oliver Sander's avatar Oliver Sander Committed by sander
Browse files

GeodesicFEAssembler and RiemannianTRSolver now directly accept dune-functions bases 

At least the harmonicmaps.cc example still works as before.
[caveat: only tested without HAVE_MPI]
Internally, at various places in the code, a dune-fufem basis is still constructed
from the dune-functions one.  These parts can be ported one by one, as dune-fufem
slowly replaces its own bases by the ones from dune-functions.

[[Imported from SVN: r10080]]
parent ef777be9
Branches
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dune/gfe/geodesicfeassembler.hh
+ 63
− 32
View file @ de9b228e
......@@ -28,19 +28,21 @@ class GeodesicFEAssembler {
public:
const Basis basis_;
const typename Basis::IndexSet basisIndexSet_;
protected:
LocalGeodesicFEStiffness<GridView,
typename Basis::LocalFiniteElement,
typename Basis::LocalView::Tree::FiniteElement,
TargetSpace>* localStiffness_;
public:
/** \brief Constructor for a given grid */
GeodesicFEAssembler(const Basis& basis,
LocalGeodesicFEStiffness<GridView,typename Basis::LocalFiniteElement, TargetSpace>* localStiffness)
LocalGeodesicFEStiffness<GridView,typename Basis::LocalView::Tree::FiniteElement, TargetSpace>* localStiffness)
: basis_(basis),
basisIndexSet_(basis_.indexSet()),
localStiffness_(localStiffness)
{}
......@@ -72,23 +74,31 @@ template <class Basis, class TargetSpace>
void GeodesicFEAssembler<Basis,TargetSpace>::
getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const
{
int n = basis_.size();
auto n = basisIndexSet_.size();
nb.resize(n, n);
ElementIterator it = basis_.getGridView().template begin<0,Dune::Interior_Partition>();
ElementIterator endit = basis_.getGridView().template end<0,Dune::Interior_Partition> ();
// A view on the FE basis on a single element
typename Basis::LocalView localView(&basis_);
auto localIndexSet = basisIndexSet_.localIndexSet();
ElementIterator it = basis_.gridView().template begin<0,Dune::Interior_Partition>();
ElementIterator endit = basis_.gridView().template end<0,Dune::Interior_Partition> ();
for (; it!=endit; ++it) {
const typename Basis::LocalFiniteElement& lfe = basis_.getLocalFiniteElement(*it);
// Bind the local FE basis view to the current element
localView.bind(*it);
localIndexSet.bind(localView);
const auto& lfe = localView.tree().finiteElement();
for (size_t i=0; i<lfe.localBasis().size(); i++) {
for (size_t i=0; i<lfe.size(); i++) {
for (size_t j=0; j<lfe.localBasis().size(); j++) {
for (size_t j=0; j<lfe.size(); j++) {
int iIdx = basis_.index(*it,i);
int jIdx = basis_.index(*it,j);
auto iIdx = localIndexSet.index(i)[0];
auto jIdx = localIndexSet.index(j)[0];
nb.add(iIdx, jIdx);
......@@ -120,32 +130,39 @@ assembleGradientAndHessian(const std::vector<TargetSpace>& sol,
gradient.resize(sol.size());
gradient = 0;
ElementIterator it = basis_.getGridView().template begin<0,Dune::Interior_Partition>();
ElementIterator endit = basis_.getGridView().template end<0,Dune::Interior_Partition> ();
// A view on the FE basis on a single element
typename Basis::LocalView localView(&basis_);
auto localIndexSet = basisIndexSet_.localIndexSet();
ElementIterator it = basis_.gridView().template begin<0,Dune::Interior_Partition>();
ElementIterator endit = basis_.gridView().template end<0,Dune::Interior_Partition> ();
for( ; it != endit; ++it ) {
const int numOfBaseFct = basis_.getLocalFiniteElement(*it).localBasis().size();
localView.bind(*it);
localIndexSet.bind(localView);
const int numOfBaseFct = localView.tree().size();
// Extract local solution
std::vector<TargetSpace> localSolution(numOfBaseFct);
for (int i=0; i<numOfBaseFct; i++)
localSolution[i] = sol[basis_.index(*it,i)];
localSolution[i] = sol[localIndexSet.index(i)[0]];
std::vector<Dune::FieldVector<double,blocksize> > localGradient(numOfBaseFct);
// setup local matrix and gradient
localStiffness_->assembleGradientAndHessian(*it, basis_.getLocalFiniteElement(*it), localSolution, localGradient);
localStiffness_->assembleGradientAndHessian(*it, localView.tree().finiteElement(), localSolution, localGradient);
// Add element matrix to global stiffness matrix
for(int i=0; i<numOfBaseFct; i++) {
int row = basis_.index(*it,i);
auto row = localIndexSet.index(i)[0];
for (int j=0; j<numOfBaseFct; j++ ) {
int col = basis_.index(*it,j);
auto col = localIndexSet.index(j)[0];
hessian[row][col] += localStiffness_->A_[i][j];
}
......@@ -153,7 +170,7 @@ assembleGradientAndHessian(const std::vector<TargetSpace>& sol,
// Add local gradient to global gradient
for (int i=0; i<numOfBaseFct; i++)
gradient[basis_.index(*it,i)] += localGradient[i];
gradient[localIndexSet.index(i)[0]] += localGradient[i];
}
......@@ -164,35 +181,42 @@ void GeodesicFEAssembler<Basis,TargetSpace>::
assembleGradient(const std::vector<TargetSpace>& sol,
Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const
{
if (sol.size()!=basis_.size())
if (sol.size()!=basisIndexSet_.size())
DUNE_THROW(Dune::Exception, "Solution vector doesn't match the grid!");
grad.resize(sol.size());
grad = 0;
ElementIterator it = basis_.getGridView().template begin<0,Dune::Interior_Partition>();
ElementIterator endIt = basis_.getGridView().template end<0,Dune::Interior_Partition>();
// A view on the FE basis on a single element
typename Basis::LocalView localView(&basis_);
auto localIndexSet = basisIndexSet_.localIndexSet();
ElementIterator it = basis_.gridView().template begin<0,Dune::Interior_Partition>();
ElementIterator endIt = basis_.gridView().template end<0,Dune::Interior_Partition>();
// Loop over all elements
for (; it!=endIt; ++it) {
localView.bind(*it);
localIndexSet.bind(localView);
// A 1d grid has two vertices
const int nDofs = basis_.getLocalFiniteElement(*it).localBasis().size();
const auto nDofs = localView.tree().size();
// Extract local solution
std::vector<TargetSpace> localSolution(nDofs);
for (int i=0; i<nDofs; i++)
localSolution[i] = sol[basis_.index(*it,i)];
localSolution[i] = sol[localIndexSet.index(i)[0]];
// Assemble local gradient
std::vector<Dune::FieldVector<double,blocksize> > localGradient(nDofs);
localStiffness_->assembleGradient(*it, basis_.getLocalFiniteElement(*it), localSolution, localGradient);
localStiffness_->assembleGradient(*it, localView.tree().finiteElement(), localSolution, localGradient);
// Add to global gradient
for (int i=0; i<nDofs; i++)
grad[basis_.index(*it,i)] += localGradient[i];
grad[localIndexSet.index(i)[0]] += localGradient[i];
}
......@@ -205,24 +229,31 @@ computeEnergy(const std::vector<TargetSpace>& sol) const
{
double energy = 0;
if (sol.size()!=basis_.size())
DUNE_THROW(Dune::Exception, "Solution vector doesn't match the grid!");
if (sol.size() != basisIndexSet_.size())
DUNE_THROW(Dune::Exception, "Coefficient vector doesn't match the function space basis!");
ElementIterator it = basis_.getGridView().template begin<0,Dune::Interior_Partition>();
ElementIterator endIt = basis_.getGridView().template end<0,Dune::Interior_Partition>();
// A view on the FE basis on a single element
typename Basis::LocalView localView(&basis_);
auto localIndexSet = basisIndexSet_.localIndexSet();
ElementIterator it = basis_.gridView().template begin<0,Dune::Interior_Partition>();
ElementIterator endIt = basis_.gridView().template end<0,Dune::Interior_Partition>();
// Loop over all elements
for (; it!=endIt; ++it) {
localView.bind(*it);
localIndexSet.bind(localView);
// Number of degrees of freedom on this element
size_t nDofs = basis_.getLocalFiniteElement(*it).localBasis().size();
size_t nDofs = localView.tree().size();
std::vector<TargetSpace> localSolution(nDofs);
for (size_t i=0; i<nDofs; i++)
localSolution[i] = sol[basis_.index(*it,i)];
localSolution[i] = sol[localIndexSet.index(i)[0]];
energy += localStiffness_->energy(*it, basis_.getLocalFiniteElement(*it), localSolution);
energy += localStiffness_->energy(*it, localView.tree().finiteElement(), localSolution);
}
......
dune/gfe/riemanniantrsolver.cc
+ 7
− 6
View file @ de9b228e
......@@ -120,10 +120,11 @@ setup(const GridType& grid,
// Assemble a Laplace matrix to create a norm that's equivalent to the H1-norm
// //////////////////////////////////////////////////////////////////////////////////////
Basis basis(grid.leafGridView());
OperatorAssembler<Basis,Basis> operatorAssembler(basis, basis);
typedef DuneFunctionsBasis<Basis> FufemBasis;
FufemBasis basis(grid.leafGridView());
OperatorAssembler<FufemBasis,FufemBasis> operatorAssembler(basis, basis);
LaplaceAssembler<GridType, typename Basis::LocalFiniteElement, typename Basis::LocalFiniteElement> laplaceStiffness;
LaplaceAssembler<GridType, typename FufemBasis::LocalFiniteElement, typename FufemBasis::LocalFiniteElement> laplaceStiffness;
typedef Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> > ScalarMatrixType;
ScalarMatrixType localA;
......@@ -158,7 +159,7 @@ setup(const GridType& grid,
// This will be used to monitor the gradient
// //////////////////////////////////////////////////////////////////////////////////////
MassAssembler<GridType, typename Basis::LocalFiniteElement, typename Basis::LocalFiniteElement> massStiffness;
MassAssembler<GridType, typename Basis::LocalView::Tree::FiniteElement, typename Basis::LocalView::Tree::FiniteElement> massStiffness;
ScalarMatrixType localMassMatrix;
operatorAssembler.assemble(massStiffness, localMassMatrix);
......@@ -216,7 +217,7 @@ setup(const GridType& grid,
TransferOperatorType pkToP1TransferMatrix;
assembleBasisInterpolationMatrix<TransferOperatorType,
P1NodalBasis<typename GridType::LeafGridView,double>,
Basis>(pkToP1TransferMatrix,p1Basis,basis);
FufemBasis>(pkToP1TransferMatrix,p1Basis,basis);
#if HAVE_MPI
// If we are on more than 1 processors, join all local transfer matrices on rank 0,
// and construct a single global transfer operator there.
......@@ -316,7 +317,7 @@ void RiemannianTrustRegionSolver<Basis,TargetSpace>::solve()
MaxNormTrustRegion<blocksize> trustRegion(globalMapper_->size(), initialTrustRegionRadius_);
#else
Basis basis(grid_->leafGridView());
MaxNormTrustRegion<blocksize> trustRegion(basis.size(), initialTrustRegionRadius_);
MaxNormTrustRegion<blocksize> trustRegion(basis.indexSet().size(), initialTrustRegionRadius_);
#endif
trustRegion.set(initialTrustRegionRadius_, scaling_);
......
src/harmonicmaps.cc
+ 16
− 12
View file @ de9b228e
......@@ -132,10 +132,14 @@ int main (int argc, char *argv[]) try
// Construct the scalar function space basis corresponding to the GFE space
//////////////////////////////////////////////////////////////////////////////////
typedef DuneFunctionsBasis<Dune::Functions::PQKNodalBasis<typename GridType::LeafGridView, 3> > FEBasis;
typedef Dune::Functions::PQKNodalBasis<typename GridType::LeafGridView, 3> FEBasis;
FEBasis feBasis(grid->leafGridView());
SolutionType x(feBasis.size());
typedef DuneFunctionsBasis<FEBasis> FufemFEBasis;
FufemFEBasis fufemFeBasis(feBasis);
SolutionType x(fufemFeBasis.size());
// /////////////////////////////////////////
// Read Dirichlet values
......@@ -146,7 +150,7 @@ int main (int argc, char *argv[]) try
BoundaryPatch<typename GridType::LeafGridView> dirichletBoundary(grid->leafGridView(), allNodes);
BitSetVector<blocksize> dirichletNodes;
constructBoundaryDofs(dirichletBoundary,feBasis,dirichletNodes);
constructBoundaryDofs(dirichletBoundary,fufemFeBasis,dirichletNodes);
// //////////////////////////
// Initial iterate
......@@ -159,7 +163,7 @@ int main (int argc, char *argv[]) try
auto pythonInitialIterate = module.get("fdf").toC<std::shared_ptr<FBase>>();
std::vector<TargetSpace::CoordinateType> v;
::Functions::interpolate(feBasis, v, *pythonInitialIterate);
::Functions::interpolate(fufemFeBasis, v, *pythonInitialIterate);
for (size_t i=0; i<x.size(); i++)
x[i] = v[i];
......@@ -173,24 +177,24 @@ int main (int argc, char *argv[]) try
// Assembler using ADOL-C
typedef TargetSpace::rebind<adouble>::other ATargetSpace;
std::shared_ptr<LocalGeodesicFEStiffness<GridType::LeafGridView,FEBasis::LocalFiniteElement,ATargetSpace> > localEnergy;
std::shared_ptr<LocalGeodesicFEStiffness<GridType::LeafGridView,FEBasis::LocalView::Tree::FiniteElement,ATargetSpace> > localEnergy;
std::string energy = parameterSet.get<std::string>("energy");
if (energy == "harmonic")
{
localEnergy.reset(new HarmonicEnergyLocalStiffness<GridType::LeafGridView, FEBasis::LocalFiniteElement, ATargetSpace>);
localEnergy.reset(new HarmonicEnergyLocalStiffness<GridType::LeafGridView, FEBasis::LocalView::Tree::FiniteElement, ATargetSpace>);
} else if (energy == "chiral_skyrmion")
{
localEnergy.reset(new GFE::ChiralSkyrmionEnergy<GridType::LeafGridView, FEBasis::LocalFiniteElement, adouble>(parameterSet.sub("energyParameters")));
localEnergy.reset(new GFE::ChiralSkyrmionEnergy<GridType::LeafGridView, FEBasis::LocalView::Tree::FiniteElement, adouble>(parameterSet.sub("energyParameters")));
} else
DUNE_THROW(Exception, "Unknown energy type '" << energy << "'");
LocalGeodesicFEADOLCStiffness<GridType::LeafGridView,
FEBasis::LocalFiniteElement,
FEBasis::LocalView::Tree::FiniteElement,
TargetSpace> localGFEADOLCStiffness(localEnergy.get());
GeodesicFEAssembler<FEBasis,TargetSpace> assembler(feBasis, &localGFEADOLCStiffness);
......@@ -239,9 +243,9 @@ int main (int argc, char *argv[]) try
}
VTKWriter<GridType::LeafGridView> vtkWriter(grid->leafGridView());
Dune::shared_ptr<VTKBasisGridFunction<FEBasis,EmbeddedVectorType> > vtkVectorField
= Dune::shared_ptr<VTKBasisGridFunction<FEBasis,EmbeddedVectorType> >
(new VTKBasisGridFunction<FEBasis,EmbeddedVectorType>(feBasis, xEmbedded, "orientation"));
Dune::shared_ptr<VTKBasisGridFunction<FufemFEBasis,EmbeddedVectorType> > vtkVectorField
= Dune::shared_ptr<VTKBasisGridFunction<FufemFEBasis,EmbeddedVectorType> >
(new VTKBasisGridFunction<FufemFEBasis,EmbeddedVectorType>(fufemFeBasis, xEmbedded, "orientation"));
vtkWriter.addVertexData(vtkVectorField);
vtkWriter.write(resultPath + "_" + energy + "_result");
......@@ -259,7 +263,7 @@ int main (int argc, char *argv[]) try
auto referenceSolution = module.get("fdf").toC<std::shared_ptr<FBase>>();
// The numerical solution, as a grid function
GFE::EmbeddedGlobalGFEFunction<FEBasis, TargetSpace> numericalSolution(feBasis, x);
GFE::EmbeddedGlobalGFEFunction<FufemFEBasis, TargetSpace> numericalSolution(feBasis, x);
// QuadratureRule for the integral of the L^2 error
QuadratureRuleKey quadKey(dim,6);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment