Skip to content
Snippets Groups Projects
Commit a5e8e71c authored by Dr. Felix Tobias Schindler's avatar Dr. Felix Tobias Schindler
Browse files

[space.continuouslagrange.fem] updated to new SpaceInterface

parent ca380a96
No related branches found
No related tags found
No related merge requests found
...@@ -19,7 +19,7 @@ ...@@ -19,7 +19,7 @@
#include "../../mapper/fem.hh" #include "../../mapper/fem.hh"
#include "../../basefunctionset/fem.hh" #include "../../basefunctionset/fem.hh"
#include "../interface.hh" #include "../continuouslagrange.hh"
#include "../constraints.hh" #include "../constraints.hh"
namespace Dune { namespace Dune {
...@@ -46,12 +46,13 @@ class FemWrapperTraits ...@@ -46,12 +46,13 @@ class FemWrapperTraits
public: public:
typedef FemWrapper<GridPartImp, polynomialOrder, RangeFieldImp, rangeDim, rangeDimCols> derived_type; typedef FemWrapper<GridPartImp, polynomialOrder, RangeFieldImp, rangeDim, rangeDimCols> derived_type;
typedef GridPartImp GridPartType; typedef GridPartImp GridPartType;
typedef typename GridPartType::GridViewType GridViewType;
static const int polOrder = polynomialOrder; static const int polOrder = polynomialOrder;
static_assert(polOrder >= 1, "Wrong polOrder given!"); static_assert(polOrder >= 1, "Wrong polOrder given!");
static const unsigned int dimDomain = GridPartType::dimension;
private: private:
typedef typename GridPartType::ctype DomainFieldType; typedef typename GridPartType::ctype DomainFieldType;
static const unsigned int dimDomain = GridPartType::dimension;
public: public:
typedef RangeFieldImp RangeFieldType; typedef RangeFieldImp RangeFieldType;
...@@ -70,17 +71,21 @@ public: ...@@ -70,17 +71,21 @@ public:
}; // class SpaceWrappedFemContinuousLagrangeTraits }; // class SpaceWrappedFemContinuousLagrangeTraits
template <class GridPartImp, int polynomialOrder, class RangeFieldImp, int r> // untested for the vector-valued case, especially ContinuousLagrangeSpaceBase
class FemWrapper<GridPartImp, polynomialOrder, RangeFieldImp, r, 1> template <class GridPartImp, int polynomialOrder, class RangeFieldImp>
: public SpaceInterface<FemWrapperTraits<GridPartImp, polynomialOrder, RangeFieldImp, r, 1>> class FemWrapper<GridPartImp, polynomialOrder, RangeFieldImp, 1, 1>
: public ContinuousLagrangeSpaceBase<FemWrapperTraits<GridPartImp, polynomialOrder, RangeFieldImp, 1, 1>,
GridPartImp::dimension, RangeFieldImp, 1, 1>
{ {
typedef SpaceInterface<FemWrapperTraits<GridPartImp, polynomialOrder, RangeFieldImp, r, 1>> BaseType; typedef ContinuousLagrangeSpaceBase<FemWrapperTraits<GridPartImp, polynomialOrder, RangeFieldImp, 1, 1>,
typedef FemWrapper<GridPartImp, polynomialOrder, RangeFieldImp, r, 1> ThisType; GridPartImp::dimension, RangeFieldImp, 1, 1> BaseType;
typedef FemWrapper<GridPartImp, polynomialOrder, RangeFieldImp, 1, 1> ThisType;
public: public:
typedef FemWrapperTraits<GridPartImp, polynomialOrder, RangeFieldImp, r, 1> Traits; typedef FemWrapperTraits<GridPartImp, polynomialOrder, RangeFieldImp, 1, 1> Traits;
typedef typename Traits::GridPartType GridPartType; typedef typename Traits::GridPartType GridPartType;
typedef typename Traits::GridViewType GridViewType;
static const int polOrder = Traits::polOrder; static const int polOrder = Traits::polOrder;
typedef typename GridPartType::ctype DomainFieldType; typedef typename GridPartType::ctype DomainFieldType;
static const unsigned int dimDomain = GridPartType::dimension; static const unsigned int dimDomain = GridPartType::dimension;
...@@ -97,6 +102,7 @@ public: ...@@ -97,6 +102,7 @@ public:
FemWrapper(const std::shared_ptr<const GridPartType>& gridP) FemWrapper(const std::shared_ptr<const GridPartType>& gridP)
: gridPart_(gridP) : gridPart_(gridP)
, gridView_(std::make_shared<GridViewType>(gridPart_->gridView()))
, backend_(std::make_shared<BackendType>(const_cast<GridPartType&>(*(gridPart_)))) , backend_(std::make_shared<BackendType>(const_cast<GridPartType&>(*(gridPart_))))
, mapper_(std::make_shared<MapperType>(backend_->mapper())) , mapper_(std::make_shared<MapperType>(backend_->mapper()))
, tmpMappedRows_(mapper_->maxNumDofs()) , tmpMappedRows_(mapper_->maxNumDofs())
...@@ -106,6 +112,7 @@ public: ...@@ -106,6 +112,7 @@ public:
FemWrapper(const ThisType& other) FemWrapper(const ThisType& other)
: gridPart_(other.gridPart_) : gridPart_(other.gridPart_)
, gridView_(other.gridView_)
, backend_(other.backend_) , backend_(other.backend_)
, mapper_(other.mapper_) , mapper_(other.mapper_)
, tmpMappedRows_(mapper_->maxNumDofs()) , tmpMappedRows_(mapper_->maxNumDofs())
...@@ -117,6 +124,7 @@ public: ...@@ -117,6 +124,7 @@ public:
{ {
if (this != &other) { if (this != &other) {
gridPart_ = other.gridPart_; gridPart_ = other.gridPart_;
gridView_ = other.gridView_;
backend_ = other.backend_; backend_ = other.backend_;
mapper_ = other.mapper_; mapper_ = other.mapper_;
tmpMappedRows_.resize(mapper_->maxNumDofs()); tmpMappedRows_.resize(mapper_->maxNumDofs());
...@@ -129,11 +137,16 @@ public: ...@@ -129,11 +137,16 @@ public:
{ {
} }
const std::shared_ptr<const GridPartType>& gridPart() const const std::shared_ptr<const GridPartType>& grid_part() const
{ {
return gridPart_; return gridPart_;
} }
const std::shared_ptr<const GridViewType>& grid_view() const
{
return gridView_;
}
const BackendType& backend() const const BackendType& backend() const
{ {
return *backend_; return *backend_;
...@@ -149,135 +162,14 @@ public: ...@@ -149,135 +162,14 @@ public:
return *mapper_; return *mapper_;
} }
BaseFunctionSetType baseFunctionSet(const EntityType& entity) const BaseFunctionSetType base_function_set(const EntityType& entity) const
{ {
return BaseFunctionSetType(*backend_, entity); return BaseFunctionSetType(*backend_, entity);
} }
template <class R>
void localConstraints(const EntityType& /*entity*/, Constraints::LocalDefault<R>& /*ret*/) const
{
static_assert((Dune::AlwaysFalse<R>::value), "Not implemented for arbitrary constraints!");
}
void
localConstraints(const EntityType& entity,
Constraints::Dirichlet<typename GridPartType::IntersectionType, RangeFieldType, true>& ret) const
{
std::set<size_t> localDirichletDofs;
const auto& gridBoundary = ret.gridBoundary();
if (polOrder == 1) {
localDirichletDofs = this->findLocalDirichletDoFs(entity, gridBoundary);
} else {
const auto& lagrangePointSet = backend_->lagrangePointSet(entity);
// loop over all intersections
const auto intersectionEndIt = gridPart_->iend(entity);
for (auto intersectionIt = gridPart_->ibegin(entity); intersectionIt != intersectionEndIt; ++intersectionIt) {
const auto& intersection = *intersectionIt;
// only work on dirichlet intersections
if (gridBoundary.dirichlet(intersection)) {
// get local face number of boundary intersection
const int intersectionIndex = intersection.indexInInside();
// iterate over face dofs and set unit row
const auto faceDofEndIt = lagrangePointSet.template endSubEntity<1>(intersectionIndex);
for (auto faceDofIt = lagrangePointSet.template beginSubEntity<1>(intersectionIndex);
faceDofIt != faceDofEndIt;
++faceDofIt) {
const size_t localDofIndex = *faceDofIt;
localDirichletDofs.insert(localDofIndex);
} // iterate over face dofs and set unit row
} // only work on dirichlet intersections
} // loop over all intersections
}
const size_t numRows = localDirichletDofs.size();
if (numRows > 0) {
const size_t numCols = mapper_->numDofs(entity);
ret.setSize(numRows, numCols);
mapper_->globalIndices(entity, tmpMappedRows_);
mapper_->globalIndices(entity, tmpMappedCols_);
size_t localRow = 0;
const RangeFieldType zero(0);
const RangeFieldType one(1);
for (auto localDirichletDofIt = localDirichletDofs.begin(); localDirichletDofIt != localDirichletDofs.end();
++localDirichletDofIt) {
const size_t& localDirichletDofIndex = *localDirichletDofIt;
ret.globalRow(localRow) = tmpMappedRows_[localDirichletDofIndex];
for (size_t jj = 0; jj < ret.cols(); ++jj) {
ret.globalCol(jj) = tmpMappedCols_[jj];
if (tmpMappedCols_[jj] == tmpMappedRows_[localDirichletDofIndex])
ret.value(localRow, jj) = one;
else
ret.value(localRow, jj) = zero;
}
++localRow;
}
} else {
ret.setSize(0, 0);
}
} // ... localConstraints(..., Dirichlet< ..., true >)
void
localConstraints(const EntityType& entity,
Constraints::Dirichlet<typename GridPartType::IntersectionType, RangeFieldType, false>& ret) const
{
std::set<size_t> localDirichletDofs;
const auto& gridBoundary = ret.gridBoundary();
if (polOrder == 1) {
localDirichletDofs = this->findLocalDirichletDoFs(entity, gridBoundary);
} else {
const auto& lagrangePointSet = backend_->lagrangePointSet(entity);
// loop over all intersections
const auto intersectionEndIt = gridPart_->iend(entity);
for (auto intersectionIt = gridPart_->ibegin(entity); intersectionIt != intersectionEndIt; ++intersectionIt) {
const auto& intersection = *intersectionIt;
// only work on dirichlet intersections
if (gridBoundary.dirichlet(intersection)) {
// get local face number of boundary intersection
const int intersectionIndex = intersection.indexInInside();
// iterate over face dofs and set unit row
const auto faceDofEndIt = lagrangePointSet.template endSubEntity<1>(intersectionIndex);
for (auto faceDofIt = lagrangePointSet.template beginSubEntity<1>(intersectionIndex);
faceDofIt != faceDofEndIt;
++faceDofIt) {
const size_t localDofIndex = *faceDofIt;
localDirichletDofs.insert(localDofIndex);
} // iterate over face dofs and set unit row
} // only work on dirichlet intersections
} // loop over all intersections
}
const size_t numRows = localDirichletDofs.size();
if (numRows > 0) {
const size_t numCols = mapper_->numDofs(entity);
ret.setSize(numRows, numCols);
mapper_->globalIndices(entity, tmpMappedRows_);
mapper_->globalIndices(entity, tmpMappedCols_);
size_t localRow = 0;
const RangeFieldType zero(0);
for (auto localDirichletDofIt = localDirichletDofs.begin(); localDirichletDofIt != localDirichletDofs.end();
++localDirichletDofIt) {
const size_t& localDirichletDofIndex = *localDirichletDofIt;
ret.globalRow(localRow) = tmpMappedRows_[localDirichletDofIndex];
for (size_t jj = 0; jj < ret.cols(); ++jj) {
ret.globalCol(jj) = tmpMappedCols_[jj];
ret.value(localRow, jj) = zero;
}
++localRow;
}
} else {
ret.setSize(0, 0);
}
} // ... localConstraints(..., Dirichlet< ..., false >)
using BaseType::computePattern;
template <class LocalGridPartType, class OtherSpaceType>
PatternType* computePattern(const LocalGridPartType& localGridPart, const OtherSpaceType& otherSpace) const
{
return BaseType::computeCodim0Pattern(localGridPart, otherSpace);
}
private: private:
std::shared_ptr<const GridPartType> gridPart_; std::shared_ptr<const GridPartType> gridPart_;
std::shared_ptr<const GridViewType> gridView_;
std::shared_ptr<const BackendType> backend_; std::shared_ptr<const BackendType> backend_;
std::shared_ptr<const MapperType> mapper_; std::shared_ptr<const MapperType> mapper_;
mutable Dune::DynamicVector<size_t> tmpMappedRows_; mutable Dune::DynamicVector<size_t> tmpMappedRows_;
......
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