diff --git a/dune/fufem/assemblers/CMakeLists.txt b/dune/fufem/assemblers/CMakeLists.txt index 24a2f37497efed6a0fae42dfcb6da966f6170323..f588564a794a1c1050ed71f00d9a5c60e5d88490 100644 --- a/dune/fufem/assemblers/CMakeLists.txt +++ b/dune/fufem/assemblers/CMakeLists.txt @@ -5,6 +5,7 @@ install(FILES boundaryfunctionalassembler.hh boundaryoperatorassembler.hh dunefunctionsoperatorassembler.hh + dunefunctionsboundaryfunctionalassembler.hh dunefunctionsfunctionalassembler.hh functionalassembler.hh integraloperatorassembler.hh diff --git a/dune/fufem/assemblers/dunefunctionsboundaryfunctionalassembler.hh b/dune/fufem/assemblers/dunefunctionsboundaryfunctionalassembler.hh new file mode 100644 index 0000000000000000000000000000000000000000..af6cc24bec1db36a188e9440f65f878ecd3c9584 --- /dev/null +++ b/dune/fufem/assemblers/dunefunctionsboundaryfunctionalassembler.hh @@ -0,0 +1,93 @@ +// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*- +// vi: set et ts=4 sw=2 sts=2: +#ifndef DUNE_FUFEM_ASSEMBERS_DUNEFUNCTIONSBOUNDARYFUNCTIONALASSEMBLER_HH +#define DUNE_FUFEM_ASSEMBERS_DUNEFUNCTIONSBOUNDARYFUNCTIONALASSEMBLER_HH + +#include <dune/istl/matrix.hh> +#include <dune/istl/matrixindexset.hh> + +#include <dune/matrix-vector/axpy.hh> + +#include "dune/fufem/functionspacebases/functionspacebasis.hh" +#include "dune/fufem/boundarypatch.hh" + + +namespace Dune { +namespace Fufem { + + +//! Generic global assembler for functionals on a boundary +template <class TrialBasis> +class DuneFunctionsBoundaryFunctionalAssembler +{ + using GridView = typename TrialBasis::GridView; + +public: + //! create assembler for grid + DuneFunctionsBoundaryFunctionalAssembler(const TrialBasis& tBasis, const BoundaryPatch<GridView>& boundaryPatch) : + trialBasis_(tBasis), + boundaryPatch_(boundaryPatch) + {} + + + + template <class VectorBackend, class LocalAssembler> + void assembleBulkEntries(VectorBackend&& vectorBackend, LocalAssembler&& localAssembler) const + { + auto trialLocalView = trialBasis_.localView(); + auto trialLocalIndexSet = trialBasis_.localIndexSet(); + + using Field = std::decay_t<decltype(vectorBackend(trialLocalIndexSet.index(0)))>; + using LocalVector = Dune::BlockVector<Dune::FieldVector<Field,1>>; + + auto localVector = LocalVector(trialLocalView.maxSize()); + + // iterate over boundary + for (auto it =boundaryPatch_.begin(); it!= boundaryPatch_.end(); ++it) + { + const auto& element = it->inside(); + trialLocalView.bind(element); + trialLocalIndexSet.bind(trialLocalView); + + localAssembler(it, localVector, trialLocalView); + + for (size_t localRow=0; localRow<trialLocalIndexSet.size(); ++localRow) + { + auto row = trialLocalIndexSet.index(localRow); + vectorBackend(row) += localVector[localRow]; + } + } + } + + template <class VectorBackend, class LocalAssembler> + void assembleBulk(VectorBackend&& vectorBackend, LocalAssembler&& localAssembler) const + { + vectorBackend.resize(trialBasis_); + vectorBackend.vector() = 0.0; + assembleBulkEntries(vectorBackend, std::forward<LocalAssembler>(localAssembler)); + } + + +protected: + const TrialBasis& trialBasis_; + const BoundaryPatch<GridView>& boundaryPatch_; +}; + + +/** + * \brief Create DuneFunctionsBoundaryFunctionalAssembler + */ +template <class TrialBasis, class BP> +auto duneFunctionsBoundaryFunctionalAssembler(const TrialBasis& trialBasis, const BP& bp) +{ + return DuneFunctionsBoundaryFunctionalAssembler<TrialBasis>(trialBasis, bp); +} + + + +} // namespace Fufem +} // namespace Dune + + +#endif +