From af1156f606b6ed5985434e0701549bb1dfecd0e7 Mon Sep 17 00:00:00 2001
From: "Dennis.Glaeser" <dennis.glaeser@iws.uni-stuttgart.de>
Date: Fri, 24 May 2019 16:02:31 +0200
Subject: [PATCH] [disc] add implementation of l2-projector
Functionality to create l2 projectors that project the solution
living in a function space basis on a domain grid into the function
space basis living on a target grid. Currently, the algorithm
assumes the grids to be non-conforming and expects a glue object.
An overload for different bases operating on the same grid can
be provided in the future, however, the current algorithm would
work as well (with the overhead of computing trivial intersections).
---
dumux/discretization/CMakeLists.txt | 1 +
.../discretization/projection/CMakeLists.txt | 3 +
dumux/discretization/projection/projector.hh | 337 ++++++++++++++++++
3 files changed, 341 insertions(+)
create mode 100644 dumux/discretization/projection/CMakeLists.txt
create mode 100644 dumux/discretization/projection/projector.hh
diff --git a/dumux/discretization/CMakeLists.txt b/dumux/discretization/CMakeLists.txt
index 3a2711c6a9..812b30de3f 100644
--- a/dumux/discretization/CMakeLists.txt
+++ b/dumux/discretization/CMakeLists.txt
@@ -2,6 +2,7 @@ add_subdirectory(box)
add_subdirectory(cellcentered)
add_subdirectory(fem)
add_subdirectory(staggered)
+add_subdirectory(projection)
install(FILES
basefvgridgeometry.hh
diff --git a/dumux/discretization/projection/CMakeLists.txt b/dumux/discretization/projection/CMakeLists.txt
new file mode 100644
index 0000000000..f3f288c353
--- /dev/null
+++ b/dumux/discretization/projection/CMakeLists.txt
@@ -0,0 +1,3 @@
+install(FILES
+projector.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/discretization/projection)
diff --git a/dumux/discretization/projection/projector.hh b/dumux/discretization/projection/projector.hh
new file mode 100644
index 0000000000..6124204b46
--- /dev/null
+++ b/dumux/discretization/projection/projector.hh
@@ -0,0 +1,337 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup Discretization
+ * \brief Contains functionality for L2-projections from one function
+ * space into another, which can live both on the same or different
+ * grids of potentially different dimensionality.
+ * \note In the case of the function space bases living on the same grid
+ * an optimized implementation could be implemented. Currently, we
+ * require a glue object containing the intersections between two
+ * grids to be provided to the free factory functions.
+ */
+#ifndef DUMUX_DISCRETIZATION_PROJECTOR_HH
+#define DUMUX_DISCRETIZATION_PROJECTOR_HH
+
+#include <string>
+#include <utility>
+#include <type_traits>
+
+#include <dune/common/timer.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/exceptions.hh>
+#include <dune/common/promotiontraits.hh>
+
+#include <dune/geometry/quadraturerules.hh>
+#include <dune/istl/matrixindexset.hh>
+
+#include <dumux/common/parameters.hh>
+#include <dumux/linear/seqsolverbackend.hh>
+#include <dumux/assembly/jacobianpattern.hh>
+#include <dumux/discretization/functionspacebasis.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup Discretization
+ * \brief Does an L2-projection from one discrete function space
+ * into another. The convenience functions makeProjectorPair
+ * or makeProjector can be used to create such a projection.
+ */
+template<class Scalar>
+class Projector
+{
+ using BlockType = Dune::FieldMatrix<Scalar, 1, 1>;
+
+public:
+ //! Export the type of the projection matrices
+ using Matrix = Dune::BCRSMatrix< BlockType >;
+
+ //! delete default constructor
+ Projector() = delete;
+
+ //! Constructor
+ Projector(Matrix&& M, Matrix&& P)
+ : M_(std::move(M)) , P_(std::move(P))
+ {}
+
+ /*!
+ * \brief Project a solution u into up
+ * \param u The solution living on the domain space
+ * \param up The projection of u into the target space
+ */
+ template< class DomainSolution, class TargetSolution>
+ void project(const DomainSolution& u, TargetSolution& up) const
+ {
+ // be picky about size of u
+ if ( u.size() != P_.M())
+ DUNE_THROW(Dune::InvalidStateException, "Vector size mismatch" );
+
+ up.resize(M_.N());
+
+ auto rhs = up;
+ P_.mv(u, rhs);
+
+ SSORCGBackend solver;
+ solver.solve(M_, up, rhs);
+ }
+
+private:
+ Matrix M_;
+ Matrix P_;
+};
+
+/*!
+ * \brief Traits class stating the type of projector between to bases
+ */
+template<class FEBasisDomain, class FEBasisTarget>
+class ProjectorTraits
+{
+ using FiniteElementDomain = typename FEBasisDomain::LocalView::Tree::FiniteElement;
+ using FiniteElementTarget = typename FEBasisTarget::LocalView::Tree::FiniteElement;
+ using ScalarDomain = typename FiniteElementDomain::Traits::LocalBasisType::Traits::RangeFieldType;
+ using ScalarTarget = typename FiniteElementTarget::Traits::LocalBasisType::Traits::RangeFieldType;
+ using Scalar = typename Dune::PromotionTraits<ScalarDomain, ScalarTarget>::PromotedType;
+public:
+ using Projector = Dumux::Projector< Scalar >;
+};
+
+
+// Projector construction details
+namespace Impl {
+
+/*!
+ * \brief Creates a projector class between two function space bases
+ * \tparam doBidirectional If false, the backward projection matrix is not assembled
+ */
+template<bool doBidirectional, class FEBasisDomain, class FEBasisTarget, class GlueType>
+std::pair< typename ProjectorTraits<FEBasisDomain, FEBasisTarget>::Projector,
+ typename ProjectorTraits<FEBasisTarget, FEBasisDomain>::Projector >
+makeProjectorPair(const FEBasisDomain& feBasisDomain,
+ const FEBasisTarget& feBasisTarget,
+ const GlueType& glue)
+{
+ // we assume that target dim <= domain dimension
+ static constexpr int domainDim = FEBasisDomain::GridView::dimension;
+ static constexpr int targetDim = FEBasisTarget::GridView::dimension;
+ static_assert(targetDim <= domainDim, "This expects target dim < domain dim, please swap arguments");
+
+ using ForwardProjector = typename ProjectorTraits<FEBasisDomain, FEBasisTarget>::Projector;
+ using BackwardProjector = typename ProjectorTraits<FEBasisTarget, FEBasisDomain>::Projector;
+
+ using ForwardProjectionMatrix = typename ForwardProjector::Matrix;
+ using BackwardProjectionMatrix = typename BackwardProjector::Matrix;
+
+ auto domainLocalView = feBasisDomain.localView();
+ auto targetLocalView = feBasisTarget.localView();
+
+ // determine mass matrix patterns (standard FE scheme pattern)
+ Dune::MatrixIndexSet backwardPatternM, forwardPatternM;
+ forwardPatternM = getFEJacobianPattern(feBasisTarget);
+ if (doBidirectional) backwardPatternM = getFEJacobianPattern(feBasisDomain);
+
+ // determine projection matrix patterns
+ Dune::MatrixIndexSet backwardPatternP, forwardPatternP;
+ forwardPatternP.resize(feBasisTarget.size(), feBasisDomain.size());
+ if (doBidirectional) backwardPatternP.resize(feBasisDomain.size(), feBasisTarget.size());
+
+ using std::max;
+ unsigned int maxBasisOrder = 0;
+ for (const auto& is : intersections(glue))
+ {
+ // "outside" element is of target domain
+ // since target dim <= domain dim there is maximum one!
+ targetLocalView.bind( is.outside(0) );
+ const auto& targetLocalBasis = targetLocalView.tree().finiteElement().localBasis();
+
+ for (unsigned int nIdx = 0; nIdx < is.neighbor(/*targetSide*/1); ++nIdx)
+ {
+ domainLocalView.bind( is.inside(nIdx) );
+ const auto& domainLocalBasis = domainLocalView.tree().finiteElement().localBasis();
+
+ // keep track of maximum basis order (used in integration)
+ maxBasisOrder = max(maxBasisOrder, max(domainLocalBasis.order(), targetLocalBasis.order()));
+
+ for (unsigned int i = 0; i < domainLocalBasis.size(); ++i)
+ for (unsigned int j = 0; j < targetLocalBasis.size(); ++j)
+ {
+ forwardPatternP.add(targetLocalView.index(j), domainLocalView.index(i));
+ if (doBidirectional) backwardPatternP.add(domainLocalView.index(i), targetLocalView.index(j));
+ }
+ }
+ }
+
+ // assemble matrices
+ ForwardProjectionMatrix forwardM, forwardP;
+ forwardPatternM.exportIdx(forwardM); forwardM = 0.0;
+ forwardPatternP.exportIdx(forwardP); forwardP = 0.0;
+
+ BackwardProjectionMatrix backwardM, backwardP;
+ if (doBidirectional)
+ {
+ backwardPatternM.exportIdx(backwardM); backwardM = 0.0;
+ backwardPatternP.exportIdx(backwardP); backwardP = 0.0;
+ }
+
+ for (const auto& is : intersections(glue))
+ {
+ const auto& targetElement = is.outside(0);
+ const auto& targetElementGeometry = targetElement.geometry();
+
+ targetLocalView.bind( targetElement );
+ const auto& targetLocalBasis = targetLocalView.tree().finiteElement().localBasis();
+
+ // perform integration over intersection geometry
+ using IsGeometry = typename std::decay_t<decltype(is.geometry())>;
+ using ctype = typename IsGeometry::ctype;
+
+ const auto& isGeometry = is.geometry();
+ const int intOrder = maxBasisOrder + 1;
+ const auto& quad = Dune::QuadratureRules<ctype, IsGeometry::mydimension>::rule(isGeometry.type(), intOrder);
+ for (auto&& qp : quad)
+ {
+ const auto weight = qp.weight();
+ const auto ie = isGeometry.integrationElement(qp.position());
+ const auto globalPos = isGeometry.global(qp.position());
+
+ std::vector< Dune::FieldVector<ctype, 1> > targetShapeVals;
+ targetLocalBasis.evaluateFunction(targetElementGeometry.local(globalPos), targetShapeVals);
+
+ // mass matrix entries target domain
+ for (unsigned int i = 0; i < targetLocalBasis.size(); ++i)
+ {
+ const auto dofIdx = targetLocalView.index(i);
+ forwardM[dofIdx][dofIdx] += ie*weight*targetShapeVals[i];
+ }
+
+ // If targetDim < domainDim, there can be several "neighbors" if
+ // targetElement is aligned with a facet of domainElement. In
+ // this case make sure the basis functions are not added
+ // multiple times! (division by numNeighbors)
+ const auto numNeighbors = is.neighbor(/*targetSide*/1);
+ for (unsigned int nIdx = 0; nIdx < numNeighbors; ++nIdx)
+ {
+ const auto& domainElement = is.inside(nIdx);
+ domainLocalView.bind( domainElement );
+ const auto& domainLocalBasis = domainLocalView.tree().finiteElement().localBasis();
+
+ std::vector< Dune::FieldVector<ctype, 1> > domainShapeVals;
+ domainLocalBasis.evaluateFunction(domainElement.geometry().local(globalPos), domainShapeVals);
+
+ // add entries in matrices
+ for (unsigned int i = 0; i < domainLocalBasis.size(); ++i)
+ {
+ const auto dofIdxDomain = domainLocalView.index(i);
+ const auto domainShapeVal = domainShapeVals[i];
+ backwardM[dofIdxDomain][dofIdxDomain] += ie*weight*domainShapeVal;
+
+ for (unsigned int j = 0; j < targetLocalBasis.size(); ++j)
+ {
+ const auto dofIdxTarget = targetLocalView.index(j);
+ const auto entry = ie*weight*domainShapeVal*targetShapeVals[j];
+
+ forwardP[dofIdxTarget][dofIdxDomain] += entry/numNeighbors;
+ if (doBidirectional)
+ backwardP[dofIdxDomain][dofIdxTarget] += entry;
+ }
+ }
+ }
+ }
+ }
+
+ // On those dofs that to not take part in any intersection we will have zeroes
+ // in the mass matrices. The right hand side will be zero because the projection
+ // matrix has no entries for those dofs as there was no intersection. Thus, we set
+ // 1.0 on the diagonal for those dofs, such that after projection, the projected
+ // solution is 0 on them.
+ for (std::size_t dofIdx = 0; dofIdx < forwardM.N(); ++dofIdx)
+ if (forwardM[dofIdx][dofIdx] == 0.0)
+ forwardM[dofIdx][dofIdx] = 1.0;
+
+ if (doBidirectional)
+ {
+ for (std::size_t dofIdx = 0; dofIdx < backwardM.N(); ++dofIdx)
+ if (backwardM[dofIdx][dofIdx] == 0.0)
+ backwardM[dofIdx][dofIdx] = 1.0;
+ }
+
+ ForwardProjector fw(std::move(forwardM), std::move(forwardP));
+ BackwardProjector bw(std::move(backwardM), std::move(backwardP));
+
+ return std::make_pair(std::move(fw), std::move(bw));
+}
+
+} // end namespace Implementation
+
+
+/*!
+ * \ingroup Discretization
+ * \brief Creates a pair of projectors between the space with
+ * basis feBasisDomain to the space with basis feBasisTarget.
+ * \param feBasisDomain The domain finite element space basis
+ * \param feBasisTarget The target finite element space basis
+ * \param glue The glue object containing the intersections between the grids.
+ * \return A pair of projectors where the first is the forward
+ * projector from the space with basis feBasisDomain to
+ * the space with basis feBasisTarget and the second
+ * does the backward projection.
+ */
+template< class FEBasisDomain, class FEBasisTarget, class GlueType >
+std::pair< typename ProjectorTraits<FEBasisDomain, FEBasisTarget>::Projector,
+ typename ProjectorTraits<FEBasisTarget, FEBasisDomain>::Projector >
+makeProjectorPair(const FEBasisDomain& feBasisDomain,
+ const FEBasisTarget& feBasisTarget,
+ GlueType glue)
+{
+ // we assume that target dim <= domain dimension
+ static constexpr int domainDim = FEBasisDomain::GridView::dimension;
+ static constexpr int targetDim = FEBasisTarget::GridView::dimension;
+ static_assert(targetDim <= domainDim, "makeProjectorPair() expects targetDim < domainDim, please swap arguments");
+
+ return Impl::makeProjectorPair<true>(feBasisDomain, feBasisTarget, glue);
+}
+
+/*!
+ * \ingroup Discretization
+ * \brief Creates a forward projector from the space feBasisDomain
+ * to the space with basis feBasisTarget.
+ * \param feBasisDomain The domain finite element space basis
+ * \param feBasisTarget The target finite element space basis
+ * \param glue The glue object containing the intersections between the grids.
+ * \return The forward projector from the space with basis feBasisDomain
+ * to the space with basis feBasisTarget.
+ */
+template< class FEBasisDomain, class FEBasisTarget, class GlueType >
+typename ProjectorTraits<FEBasisDomain, FEBasisTarget>::Projector
+makeProjector(const FEBasisDomain& feBasisDomain,
+ const FEBasisTarget& feBasisTarget,
+ GlueType glue)
+{
+ // we assume that target dim <= domain dimension
+ static constexpr int domainDim = FEBasisDomain::GridView::dimension;
+ static constexpr int targetDim = FEBasisTarget::GridView::dimension;
+ static_assert(targetDim <= domainDim, "makeProjectorPair() expects targetDim < domainDim, please swap arguments");
+
+ return Impl::makeProjectorPair<false>(feBasisDomain, feBasisTarget, glue).first;
+}
+
+} // end namespace Dumux
+
+#endif // DUMUX_PROJECTOR_HH
--
GitLab