Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision
  • bugfix/any_hh-use_std_remove_const
  • bugfix/matrix-size-in-dune-functions-assembler
  • bugfix/test-memory-leak
  • debug/functionspacebasistest
  • experimental/test-core-ci
  • face_hierarchy
  • feature/DG-Transferoperator-Assembler
  • feature/dune-functions-assembler-with-skeleton
  • feature/dune-functions-assemblers
  • feature/elastictensor-interface
  • feature/exploit-tensor-quadrature
  • feature/handout-wrapped-functionsbasis
  • feature/matrix-free-wip
  • feature/no-auto-ptr-via-boost
  • feature/powerBasis
  • feature/thread-parallel-assembler
  • feature/update-ci
  • improve_grid
  • introduce-periodic-basis
  • master
  • maxka/conformingbasis
  • missing-include
  • move-makefunction-to-namespace-dune-fufem
  • pipping/dune-fufem-feature/heterogeneous-second-order
  • releases/2.0-1
  • releases/2.1-1
  • releases/2.2-1
  • releases/2.3-1
  • releases/2.3-2
  • releases/2.4-1
  • releases/2.4-2
  • releases/2.5-1
  • releases/2.6-1
  • releases/2.7
  • releases/2.8
  • releases/2.9
  • subgridassembler-rework
  • temp/test-CI-with-subgrid-enabled
  • test/use-default-ci-setup
  • tobias-patches
  • tobias-patches-membrane
  • work-matrixindexset
  • subversion->git
43 results

Target

Select target project
  • lisa_julia.nebel_at_tu-dresden.de/dune-fufem
  • akbib/dune-fufem
  • patrick.jaap_at_tu-dresden.de/dune-fufem
  • burchardt_at_igpm.rwth-aachen.de/dune-fufem
4 results
Select Git revision
  • bugfix/matrix-size-in-dune-functions-assembler
  • debug/functionspacebasistest
  • experimental/test-core-ci
  • feature/DG-Transferoperator-Assembler
  • feature/assemble-transfer-operator-with-powerbasis
  • feature/dune-functions-assembler-with-skeleton
  • feature/dune-functions-assemblers
  • feature/elastictensor-interface
  • feature/matrix-free-wip
  • feature/no-auto-ptr-via-boost
  • feature/powerBasis
  • improve_grid
  • lisa-master
  • master
  • maxka/conformingbasis
  • missing-include
  • move-makefunction-to-namespace-dune-fufem
  • releases/2.0-1
  • releases/2.1-1
  • releases/2.2-1
  • releases/2.3-1
  • releases/2.3-2
  • releases/2.4-1
  • releases/2.4-2
  • releases/2.5-1
  • releases/2.6-1
  • releases/2.7
  • subgridassembler-rework
  • temp/test-CI-with-subgrid-enabled
  • subversion->git
30 results
Show changes
Commits on Source (201)
101 additional commits have been omitted to prevent performance issues.
Hide whitespace changes
Inline Side-by-side
Showing
with 592 additions and 878 deletions
.gitignore
View file @ b215e298
# ignore all build folders
/build*/
# ignore backup files
*~
# ignore Python files
*.pyc
# ignore files generated during python setup.py sdist
MANIFEST
_skbuild/
dist
# ignore macOS filesystem
.DS_Store
# /
/Makefile
/Makefile.in
/config.h
/configure
/aclocal.m4
/dependencies.m4
/autom4te.cache
/depcomp
/install-sh
/missing
/mkinstalldirs
/libtool
/dune-fufem.pc
/semantic.cache
/configure.lineno
/stamp-h1
/dune-fufem-*.tar.gz
/dune-fufem-?.?
/ltmain.sh
/am
/.libs
/compile
/test-driver
# /bin/
/bin/*.pyc
# /cmake/
/cmake/Makefile
/cmake/Makefile.in
/cmake/modules/Makefile
/cmake/modules/Makefile.in
# Standard cmake build directory
/build-cmake
# /doc/
/doc/.deps
/doc/Makefile
/doc/Makefile.in
/doc/semantic.cache
/doc/*.html
/doc/*.out
/doc/*.pdf
/doc/*.ps
/doc/*.toc
/doc/*.aux
/doc/*.bbl
/doc/*.blg
/doc/*.log
/doc/*.dvi
/doc/Makefile.dist.in
/doc/Makefile.dist
# /doc/doxygen/
/doc/doxygen/html
/doc/doxygen/html-dist
/doc/doxygen/Makefile
/doc/doxygen/Makefile.in
/doc/doxygen/semantic.cache
/doc/doxygen/.deps
/doc/doxygen/Doxyfile.in
/doc/doxygen/Doxyfile
/doc/doxygen/doxygen-tag
/doc/doxygen/doxygen.log
/doc/doxygen/doxyerr.log
# /dune/
/dune/Makefile.in
/dune/Makefile
# /dune/fufem/
/dune/fufem/Makefile
/dune/fufem/Makefile.in
/dune/fufem/semantic.cache
/dune/fufem/.*.swp
# /dune/fufem/assemblers/
/dune/fufem/assemblers/.*.swp
/dune/fufem/assemblers/Makefile.in
/dune/fufem/assemblers/Makefile
# /dune/fufem/assemblers/localassemblers/
/dune/fufem/assemblers/localassemblers/Makefile.in
/dune/fufem/assemblers/localassemblers/Makefile
# /dune/fufem/estimators/
/dune/fufem/estimators/Makefile.in
/dune/fufem/estimators/Makefile
# /dune/fufem/functions/
/dune/fufem/functions/Makefile.in
/dune/fufem/functions/Makefile
# /dune/fufem/functionspacebases/
/dune/fufem/functionspacebases/Makefile.in
/dune/fufem/functionspacebases/Makefile
# /dune/fufem/functiontools/
/dune/fufem/functiontools/Makefile.in
/dune/fufem/functiontools/Makefile
# /dune/fufem/mechanics/
/dune/fufem/mechanics/Makefile.in
/dune/fufem/mechanics/Makefile
# /dune/fufem/python/
/dune/fufem/python/Makefile.in
/dune/fufem/python/Makefile
# /dune/fufem/quadraturerules/
/dune/fufem/quadraturerules/Makefile.in
/dune/fufem/quadraturerules/Makefile
# /dune/fufem/test/
/dune/fufem/test/# Test executables
/dune/fufem/test/arithmetictest
/dune/fufem/test/basisgridfunctiontest
/dune/fufem/test/basisinterpolationmatrixtest
/dune/fufem/test/basisinterpolatortest
/dune/fufem/test/bsplinebasistest
/dune/fufem/test/boundarypatchtest
/dune/fufem/test/coarsegridfunctionwrappertest
/dune/fufem/test/composedfunctiontest
/dune/fufem/test/constantfunctiontest
/dune/fufem/test/dunepythontest
/dune/fufem/test/functionintegratortest
/dune/fufem/test/functionspacebasistest
/dune/fufem/test/generalizedlaplaceassemblertest
/dune/fufem/test/gradientassemblertest
/dune/fufem/test/gridfunctiontest
/dune/fufem/test/gridfunctionadaptortest
/dune/fufem/test/h1functionalassemblertest
/dune/fufem/test/integraloperatorassemblertest
/dune/fufem/test/laplaceassemblertest
/dune/fufem/test/pgmtest
/dune/fufem/test/polynomialtest
/dune/fufem/test/ppmtest
/dune/fufem/test/serializationtest
/dune/fufem/test/staticmatrixtoolstest
/dune/fufem/test/subgridxyfunctionalassemblertest
/dune/fufem/test/sumfunctiontest
/dune/fufem/test/tensortest
/dune/fufem/test/# Build system stuff
/dune/fufem/test/Makefile
/dune/fufem/test/Makefile.in
/dune/fufem/test/.deps
/dune/fufem/test/*.log
/dune/fufem/test/*.trs
/dune/fufem/test/*.o
# /dune/fufem/utilities/
/dune/fufem/utilities/Makefile.in
/dune/fufem/utilities/Makefile
# /m4/
/m4/Makefile.in
/m4/Makefile
.gitlab-ci.yml
View file @ b215e298
......@@ -2,13 +2,31 @@
# Install external dependencies
before_script:
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-alugrid.git
- duneci-install-module https://git.imp.fu-berlin.de/agnumpde/dune-subgrid.git
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-subgrid.git
- duneci-install-module https://git.imp.fu-berlin.de/agnumpde/dune-matrix-vector.git
dune:git clang C++17:
image: registry.dune-project.org/docker/ci/dune:git-debian-10-clang-7-libcpp-17
dune:2.8 debian-11 gcc-10 C++20:
variables:
DUNECI_BRANCH: releases/2.8
image: registry.dune-project.org/docker/ci/dune:2.8-debian-11-gcc-10-20
before_script:
script: duneci-standard-test
dune:git gcc-8 C++17:
image: registry.dune-project.org/docker/ci/dune:git-debian-10-gcc-8-17
dune:2.9 debian-11 gcc-10 C++20:
variables:
DUNECI_BRANCH: releases/2.9
image: registry.dune-project.org/docker/ci/dune:2.9-debian-11-gcc-10-20
before_script:
script: duneci-standard-test
# dune:git debian-10 clang-7 C++17:
# image: registry.dune-project.org/docker/ci/dune:git-debian-10-clang-7-libcpp-17
# script: duneci-standard-test
dune:git ubuntu-20-04 clang-10 C++20:
image: registry.dune-project.org/docker/ci/dune:git-ubuntu-20.04-clang-10-20
script: duneci-standard-test
dune:git debian-11 gcc-10 C++20:
image: registry.dune-project.org/docker/ci/dune:git-debian-11-gcc-10-20
script: duneci-standard-test
CHANGELOG.md
View file @ b215e298
# Master (will become release 2.8)
# Master (will become release 2.10)
- ...
## Deprecations
- Objects of type `VirtualFunction` and `DifferentiableVirtualFunction` cannot be created
from embedded Python anymore. Both have been deprecated since before the release of
`dune-common` 2.7.
- The following implementations of the deprecated `Dune::VirtualFunction` interface
have been removed: `ConstantFunction`, `SumFunction`, `SumGridFunction`,
`ComposedFunction`, `ComposedGridFunction`, `Polynomial`.
- The header `dune/fufem/concept.hh` has been deprecated. Please use the concept utilities
from `dune/common/concept.hh` instead.
## Removals
- The method `assembleBasisInterpolationMatrix` has been removed. It only worked
for old `dune-fufem`-style function space bases.
- The class `AmiraMeshBasisWriter` has been removed. It relied on AmiraMesh support
in `dune-grid`, which has been removed before the 2.9 release.
- The class `VTKBasisWriter` has been removed. It only worked
for old `dune-fufem`-style function space bases.
- The deprecated file `ulis_tools.hh` has been removed. It contained a short list of
basic linear algebra methods that are now covered by the `dune-matrix-vector` module.
- The deprecated method `Dune::Fufem::istlVectorBackend` has been removed.
Please use `Dune::Functions::istlVectorBackend` instead!
- The deprecated class `MassAssembler` (from the global namespace) has been removed.
It still relied on the old `dune-fufem` function space basis interface.
Please use `Dune::Fufem::MassAssembler` instead!
- The class `NamedFunctionMap` based on `Dune::VirtualFunction` has been removed.
Use `std::map<:std::string, std::function<Range(Domain)>>` as a replacement.
# 2.9 Release
- Various improvements to the `MappedMatrix` class
- `MappedMatrix` objects can now be printed using `printmatrix` (from `dune-istl`).
- Iteration over rows is now implemented. This implies that range-based `for`-loops
over rows also work.
# 2.8 Release
- constructBoundaryDofs:
- Small interface change in the template parameters: drop `blocksize` and replace it by `BitSetVector`
- The method can now handle generic `dune-functions` basis types, as long as we have consistency in the data types
## Deprecations and removals
- assembleGlobalBasisTransferMatrix:
- Support for all bases in `dune-functions` compatible form added
- ...
- `istlMatrixBackend` can now be used with `MultiTypeBlockMatrix`
- The class `DuneFunctionsLocalMassAssembler` has been renamed to `MassAssembler`,
moved into the namespace `Dune::Fufem`, and to the file `massassembler.hh`.
The old class is still there, but it is deprecated now.
CMakeLists.txt
View file @ b215e298
......@@ -11,6 +11,7 @@ include(DuneMacros)
dune_project()
add_subdirectory("cmake")
add_subdirectory("doc")
add_subdirectory("dune")
......
cmake/CMakeLists.txt 0 → 100644
View file @ b215e298
add_subdirectory(modules)
cmake/modules/AddAdolcFlags.cmake
View file @ b215e298
......@@ -11,7 +11,7 @@
function(add_dune_adolc_flags _targets)
if(ADOLC_FOUND)
foreach(_target ${_targets})
target_link_libraries(${_target} ${ADOLC_LIBRARIES})
target_link_libraries(${_target} PUBLIC ${ADOLC_LIBRARIES})
get_target_property(_props ${_target} COMPILE_FLAGS)
string(REPLACE "_props-NOTFOUND" "" _props "${_props}")
set_target_properties(${_target} PROPERTIES COMPILE_FLAGS
......
cmake/modules/AddHDF5Flags.cmake
View file @ b215e298
......@@ -11,7 +11,7 @@
function(add_dune_hdf5_flags _targets)
if(HDF5_FOUND)
foreach(_target ${_targets})
target_link_libraries(${_target} ${HDF5_LIBRARIES} ${HDF5_HL_LIBRARIES})
target_link_libraries(${_target} PUBLIC ${HDF5_LIBRARIES} ${HDF5_HL_LIBRARIES})
set_property(TARGET ${_target} APPEND PROPERTY INCLUDE_DIRECTORIES ${HDF5_INCLUDE_DIRS})
endforeach(_target ${_targets})
endif(HDF5_FOUND)
......
cmake/modules/AddPythonLibsFlags.cmake
View file @ b215e298
......@@ -9,15 +9,16 @@
#
function(add_dune_pythonlibs_flags _targets)
if(PYTHONLIBS_FOUND)
foreach(_target ${_targets})
target_link_libraries(${_target} ${PYTHON_LIBRARIES})
# target_include_directories(${_target} PRIVATE ${PYTHON_INCLUDE_DIRS})
# target_compile_definitions(${_target} PRIVATE "-DHAVE_PYTHON")
# target_compile_options(${_target} PRIVATE "-fno-strict-aliasing")
foreach(_target ${_targets})
if(Python3_FOUND)
target_link_libraries(${_target} PUBLIC ${Python3_LIBRARIES})
set_property(TARGET ${_target} APPEND PROPERTY INCLUDE_DIRECTORIES ${Python3_INCLUDE_DIRS})
# backward compatibility with 2.7 core tests
elseif(PYTHONLIBS_FOUND)
target_link_libraries(${_target} PUBLIC ${PYTHON_LIBRARIES})
set_property(TARGET ${_target} APPEND PROPERTY INCLUDE_DIRECTORIES ${PYTHON_INCLUDE_DIRS})
endif()
set_property(TARGET ${_target} APPEND PROPERTY COMPILE_DEFINITIONS "HAVE_PYTHON")
set_property(TARGET ${_target} APPEND PROPERTY COMPILE_OPTIONS "-fno-strict-aliasing")
endforeach(_target ${_targets})
endif(PYTHONLIBS_FOUND)
endforeach(_target ${_targets})
endfunction(add_dune_pythonlibs_flags)
cmake/modules/CMakeLists.txt
View file @ b215e298
set(modules "DuneFufemMacros.cmake")
set(modules
AddAdolcFlags.cmake
AddHDF5Flags.cmake
AddPythonLibsFlags.cmake
DuneFufemMacros.cmake
FindAdolc.cmake
)
install(FILES ${modules} DESTINATION ${DUNE_INSTALL_MODULEDIR})
cmake/modules/DuneFufemMacros.cmake
View file @ b215e298
......@@ -36,11 +36,12 @@ macro(dune_fufem_add_copy_dependency target file_name)
add_dependencies(${target} "${target}_copy_${file_name}")
endmacro(dune_fufem_add_copy_dependency)
# Backward compatibility: Test here for PYTHONLIBS if dune-common is < 2.8
# Can be removed once tests against 2.7 core are dropped
if(${dune-common_VERSION} VERSION_LESS_EQUAL "2.7.9")
find_package(PythonLibs)
endif()
set(Python_ADDITIONAL_VERSIONS 2.7 3 3.2)
find_package(PythonLibs)
find_package(Adolc)
find_package(HDF5 COMPONENTS C HL)
......
dune.module
View file @ b215e298
Module: dune-fufem
Depends: dune-common (>= 2.7) dune-geometry (>= 2.7) dune-grid (>= 2.7) dune-istl (>= 2.7) dune-localfunctions (>= 2.7) dune-functions (>= 2.7) dune-matrix-vector (>= 2.7)
Version: 2.8-git
Maintainer: Carsten Gräser (graeser@math.fu-berlin.de)
Depends: dune-common (>= 2.8) dune-geometry (>= 2.8) dune-grid (>= 2.8) dune-istl (>= 2.8) dune-localfunctions (>= 2.8) dune-functions (>= 2.8) dune-matrix-vector (>= 2.8)
Version: 2.10-git
Maintainer: Carsten Gräser (graeser@math.fau.de)
Suggests: dune-alugrid dune-subgrid
dune/fufem/CMakeLists.txt
View file @ b215e298
......@@ -13,7 +13,6 @@ add_subdirectory(utilities)
add_subdirectory(test)
install(FILES
any.hh
arcofcircle.hh
arithmetic.hh
basistraceindexset.hh
......@@ -52,5 +51,4 @@ install(FILES
staticmatrixtools.hh
surfmassmatrix.hh
symmetrictensor.hh
ulis_tools.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/fufem)
dune/fufem/any.hh deleted 100644 → 0
View file @ 7a3c7508
#ifndef ANY_HH
#define ANY_HH
#include <list>
#include <memory>
#include <dune/common/typetraits.hh>
#include <dune/common/classname.hh>
/** \brief Class to encapsulate objects of any type
*
* This class is inspired by boost::any
*/
class any
{
private:
class AnyObjectInterface
{
public:
virtual ~AnyObjectInterface() {}
virtual AnyObjectInterface* clone() const
{
return new AnyObjectInterface(*this);
}
virtual std::string className() const
{
return "";
}
};
template <class T>
class AnyObjectWrapper :
public AnyObjectInterface
{
public:
AnyObjectWrapper(const T& t) : t_(t) {}
virtual ~AnyObjectWrapper() {}
virtual AnyObjectWrapper* clone() const
{
return new AnyObjectWrapper(*this);
}
const T& get() const
{
return t_;
}
T& get()
{
return t_;
}
std::string className() const
{
return Dune::className<T>();
}
private:
T t_;
};
public:
/**
* \brief Create an empty any object.
*/
any() :
p_(0)
{}
/**
* \brief Create an any object storing a copy of some value.
*
* \tparam T Type of the stored value
*
* \param t Store a copy of this value
*/
template <class T>
any(const T& t)
{
p_ = new AnyObjectWrapper<T>(t);
}
/**
* \brief Copy content of another any object.
*
* \param other Copy content of this any object.
*
* The current content of this any object is destroyed.
*/
any(const any& other)
{
if (other.p_)
p_ = other.p_->clone();
else
p_ = 0;
}
/**
* \brief Copy content of another any object.
*
* \param other Copy content of this any object.
*
* The current content of this any object is destroyed.
*/
any& operator= (const any & other)
{
AnyObjectInterface* pOld = p_;
if (other.p_)
p_ = other.p_->clone();
else
p_ = 0;
delete pOld;
return *this;
}
/**
* \brief Destroy any object and its content.
*/
~any()
{
delete p_;
}
template<class T>
T* get()
{
AnyObjectWrapper<T>* wrappedT = dynamic_cast<AnyObjectWrapper<T>*>(p_);
if (wrappedT==0)
return 0;
return &(wrappedT->get());
}
template<class T>
const T* get() const
{
const AnyObjectWrapper<T>* wrappedT = dynamic_cast<const AnyObjectWrapper<T>*>(p_);
if (wrappedT==0)
return 0;
return &(wrappedT->get());
}
std::string className() const
{
if (p_)
return p_->className();
return "";
}
private:
AnyObjectInterface* p_;
};
/**
* \brief Exception for failed any_cast calls
*/
struct bad_any_cast {};
/**
* \brief Try to cast reference to an any object to a const value or reference
*
* If the any object does not contain an object of this type a bad_any_cast exception is thrown.
*/
template<typename T>
T any_cast(const any& operand)
{
using RawT = std::decay_t<T>;
const RawT* t = operand.template get<RawT>();
if (t==0)
throw bad_any_cast();
return *t;
}
/**
* \brief Try to cast reference to an any object to a value or reference
*
* If the any object does not contain an object of this type a bad_any_cast exception is thrown.
*/
template<typename T>
T any_cast(any& operand)
{
using RawT = std::decay_t<T>;
RawT* t = operand.template get<RawT>();
if (t==0)
throw bad_any_cast();
return *t;
}
/**
* \brief Try to cast a pointer to a const any object to a pointer to a const value
*
* Returns 0 if the cast fails.
*/
template<typename T>
const T* any_cast(const any* operand)
{
return operand->template get<typename Dune::remove_const<T>::type>();
}
/**
* \brief Try to cast a pointer to a any object to a pointer to a value
*
* Returns 0 if the cast fails.
*/
template<typename T>
T* any_cast(any* operand)
{
return operand->template get<typename Dune::remove_const<T>::type>();
}
/**
* \brief Container for any type of shared_ptr objects
*
* This class stores shared_ptr<T> objects for any type T in a list.
* Its purpose is to control life time of objects allocated on the heap
* without storing a separate shared_ptr for each object.
*/
class any_shared_ptr_list :
public std::list<any>
{
public:
/**
* \brief Store a copy of the shared_ptr in the container.
*
* \returns Reference to the stored shared_ptr
*/
template<class T>
std::shared_ptr<T>& push_back_ptr(std::shared_ptr<T>& t)
{
this->push_back(any(t));
return any_cast<typename std::shared_ptr<T>&>(this->back());
}
/**
* \brief Create a shared_ptr and store it in the container.
*
* This creates a shared_ptr encapsulating the raw pointer and stores
* it in the container. This e.g. allows to do
* \code
* any_shared_ptr_list c;
* A* p1 = c.push_back_ptr(new SomeClassDerivedFromA).get();
* B* p2 = c.push_back_ptr(new SomeClassDerivedFromB).get();
* B* p3 = c.push_back_ptr(new AnotherClassDerivedFromB).get();
* \endcode
* Here the container takes ownership for the created objects and thus
* controls their life time. If you want to share ownership you can just
* copy the returned shared_ptr.
*
* \returns Reference to the stored shared_ptr
*/
template<class T>
std::shared_ptr<T>& push_back_ptr(T* t)
{
std::shared_ptr<T> p(t);
return push_back_ptr<T>(p);
}
};
#endif
dune/fufem/assemblers/basisinterpolationmatrixassembler.hh
View file @ b215e298
......@@ -5,6 +5,7 @@
#include <type_traits>
#include <typeinfo>
#include <unordered_set>
#include <dune/common/bitsetvector.hh>
......@@ -14,7 +15,13 @@
#include <dune/matrix-vector/addtodiagonal.hh>
#include <dune/fufem/functions/cachedcomponentwrapper.hh>
#include <dune/localfunctions/common/virtualinterface.hh>
#include <dune/fufem/assemblers/istlbackend.hh>
#include <dune/functions/functionspacebases/basistags.hh>
#include <dune/typetree/visitor.hh>
#include <dune/typetree/pairtraversal.hh>
/** \brief Wrapper that extracts a single local basis function. */
template<class LocalBasis, class Base>
......@@ -48,179 +55,6 @@ class LocalBasisComponentWrapper :
const LocalBasis& localBasis_;
};
/** \brief Method that computes the interpolation matrix mapping coefficients from one basis to
* the coefficients of another basis.
* Note that the interpolation only works exactly if one basis can be represented by the other.
* (Here denoted by coarse and fine basis)
* In this case multiplication of the matrix from the left represents the interpolation.
*/
template<class MatrixType, class BasisType0, class BasisType1>
static void assembleBasisInterpolationMatrix(MatrixType& matrix, const BasisType0& coarseBasis, const BasisType1& fineBasis)
{
typedef typename BasisType0::GridView GridView;
static_assert(std::is_same<GridView, typename BasisType1::GridView>::value, "GridView's don't match!");
const GridView& gridView = coarseBasis.getGridView();
int rows = fineBasis.size();
int cols = coarseBasis.size();
typedef typename BasisType0::LocalFiniteElement FEType0;
typedef typename BasisType1::LocalFiniteElement FEType1;
typedef typename Dune::LocalFiniteElementFunctionBase<FEType0>::type FunctionBaseClass;
typedef LocalBasisComponentWrapper<typename FEType0::Traits::LocalBasisType, FunctionBaseClass> LocalBasisWrapper;
matrix.setSize(rows,cols);
matrix.setBuildMode(MatrixType::random);
auto eIt = gridView.template begin<0>();
auto eEndIt = gridView.template end<0>();
// ///////////////////////////////////////////
// Determine the indices present in the matrix
// /////////////////////////////////////////////////
// Only handle every dof once
Dune::BitSetVector<1> processed(fineBasis.size());
for (size_t i=0; i<processed.size();i++)
if (fineBasis.isConstrained(i))
processed[i] = true;
Dune::MatrixIndexSet indices(rows, cols);
for (; eIt != eEndIt; ++eIt) {
// Get local finite elements
const FEType0& lfe0 = coarseBasis.getLocalFiniteElement(*eIt);
const FEType1& lfe1 = fineBasis.getLocalFiniteElement(*eIt);
const size_t numBaseFct0 = lfe0.localBasis().size();
const size_t numBaseFct1 = lfe1.localBasis().size();
// check if all components have already been processed
bool allProcessed = true;
for(size_t i=0; i<numBaseFct1; ++i)
allProcessed = allProcessed and processed[fineBasis.index(*eIt, i)][0];
// preallocate vector for function evaluations
std::vector<typename FEType0::Traits::LocalBasisType::Traits::RangeType> values(numBaseFct0);
// Extract single basis functions into a format that can be used within local interpolation
LocalBasisWrapper basisFctj(lfe0.localBasis(),0);
for (size_t j=0; j<numBaseFct0; j++)
{
// wrap each local basis function as a local function.
basisFctj.setIndex(j);
// Interpolate j^th base function by the fine basis
lfe1.localInterpolation().interpolate(basisFctj, values);
int globalCoarse = coarseBasis.index(*eIt,j);
for (size_t i=0; i<numBaseFct1; i++) {
if (std::abs(values[i])<1e-12)
continue;
int globalFine = fineBasis.index(*eIt,i);
if (processed[globalFine][0])
continue;
if (coarseBasis.isConstrained(globalCoarse)) {
const auto& lin = coarseBasis.constraints(globalCoarse);
for (size_t k=0; k<lin.size(); k++)
indices.add(globalFine, lin[k].index);
} else
indices.add(globalFine, globalCoarse);
}
}
for (size_t i=0; i<numBaseFct1; i++)
processed[fineBasis.index(*eIt,i)]=true;
}
indices.exportIdx(matrix);
matrix = 0;
// /////////////////////////////////////////////
// Compute the matrix
// /////////////////////////////////////////////
// Only handle every dof once
processed.unsetAll();
for (size_t i=0; i<processed.size();i++)
if (fineBasis.isConstrained(i))
processed[i] = true;
eIt = gridView.template begin<0>();
for (; eIt != eEndIt; ++eIt) {
// Get local finite element
const FEType0& lfe0 = coarseBasis.getLocalFiniteElement(*eIt);
const FEType1& lfe1 = fineBasis.getLocalFiniteElement(*eIt);
const size_t numBaseFct0 = lfe0.localBasis().size();
const size_t numBaseFct1 = lfe1.localBasis().size();
// check if all components have already been processed
bool allProcessed = true;
for(size_t i=0; i<numBaseFct1; ++i)
allProcessed = allProcessed and processed[fineBasis.index(*eIt, i)][0];
if (not allProcessed) {
// preallocate vector for function evaluations
std::vector<typename FEType0::Traits::LocalBasisType::Traits::RangeType> values(numBaseFct0);
// Extract single basis functions into a format that can be used within local interpolation
LocalBasisWrapper basisFctj(lfe0.localBasis(),0);
for (size_t j=0; j<numBaseFct0; j++)
{
// wrap each local basis function as a local function.
basisFctj.setIndex(j);
int globalCoarse = coarseBasis.index(*eIt,j);
// Interpolate j^th base function by the fine basis
lfe1.localInterpolation().interpolate(basisFctj, values);
for (size_t i=0; i<numBaseFct1; i++) {
if (std::abs(values[i])<1e-12)
continue;
int globalFine = fineBasis.index(*eIt,i);
if (processed[globalFine][0])
continue;
if (coarseBasis.isConstrained(globalCoarse)) {
const auto& lin = coarseBasis.constraints(globalCoarse);
for (size_t k=0; k<lin.size(); k++)
Dune::MatrixVector::addToDiagonal(matrix[globalFine][lin[k].index],lin[k].factor*values[i]);
} else
Dune::MatrixVector::addToDiagonal(matrix[globalFine][globalCoarse],values[i]);
}
}
for (size_t i=0; i<numBaseFct1; i++)
processed[fineBasis.index(*eIt,i)] = true;
}
}
}
/** \brief Wrapper for evaluating local functions in another geometry */
template<class Function, class Geometry>
......@@ -229,20 +63,20 @@ class OtherGeometryWrapper
public:
using Traits = typename Function::Traits;
OtherGeometryWrapper(const Function& function, const Geometry& geometry)
: function_(function)
, geometry_(geometry)
{}
// forward the evaluation into the other geometry
template<class X, class Y>
void evaluate(const X& x, Y& y) const
template<class X>
auto operator()(const X& x) const
{
function_.evaluate(geometry_.global(x), y);
return function_(geometry_.global(x));
}
private:
const Function function_;
......@@ -280,16 +114,144 @@ private:
namespace Impl{
/** \brief Visitor for interpolating the local basis of a leaf node by a fine local basis
*
* The LeafNodeInterpolationVisitor traverses the local tree of the coarse basis at an
* element bound to the coarse localView. At every leaf node the local basis is extracted and
* interpolated by the corresponding local basis of the fine element (bound to the fine localView).
*
* This works for all dune-functions compatible local basis trees.
*
* \todo In case of powerNodes we could recycle the interpolation for all child nodes to save
* some assembly time.
*/
template <class CLV, class FLV, class G, class B, class M, class I>
class LeafNodeInterpolationVisitor
: public Dune::TypeTree::TreePairVisitor
, public Dune::TypeTree::DynamicTraversal
{
public:
using CoarseLocalView = CLV;
using FineLocalView = FLV;
using Geometry = G;
using BitSet = B;
using Matrix = M;
using IndexSet = I;
/** \brief Constructor with all necessary context information
*
* \param [in] clv (coarseLocalView) coarse element information and localBasis
* \param [in] flv (fineLocalView) fine element information and localBasis
* \param [in] g (geometry) mapping fine -> coarse
* \param [in] p (processed) set of already processed dof's
* \param [in] t (tolerance) threshold for considering an interpolated value as zero
* \param [out] m (matrix) resulting tranferoperator matrix (coarse basis) -> (fine basis) wrapped in ISTLBackend
* \param [out] i (indices) index set of the matrix
* \param [in] e (exportIndices) boolean switch: true = add indices to the index set
* false = increase the matrix entries
*/
LeafNodeInterpolationVisitor(const CLV& clv, const FLV& flv, const G& g, const B& p, double t, M& m, I& i, bool e)
: coarseLocalView_(clv)
, fineLocalView_(flv)
, geometry_(g)
, processed_(p)
, tolerance_(t)
, matrix_(m)
, indices_(i)
, exportIndices_(e)
{}
template<class CoarseNode, class FineNode, class TreePath>
void leaf(CoarseNode& coarseNode, FineNode& fineNode, TreePath treePath)
{
const auto& coarseFiniteElement = coarseNode.finiteElement();
const auto& fineFiniteElement = fineNode.finiteElement();
const auto& coarseLocalBasis = coarseFiniteElement.localBasis();
const auto& fineLocalBasis = fineFiniteElement.localBasis();
// Hack: The RangeFieldType is the best guess of a suitable type for coefficients we have here
using CoefficientType = typename std::decay_t<decltype(coarseLocalBasis)>::Traits::RangeFieldType;
// interpolated coefficients for the fine local basis
std::vector<CoefficientType> coefficients(fineLocalBasis.size());
using CoarseFEType = std::decay_t<decltype(coarseFiniteElement)>;
using FunctionTraits = typename CoarseFEType::Traits::LocalBasisType::Traits;
using LocalBasisWrapper = LocalBasisComponentWrapper<typename CoarseFEType::Traits::LocalBasisType, FunctionTraits>;
LocalBasisWrapper basisFctj(coarseLocalBasis,0);
for (size_t j=0; j<coarseLocalBasis.size(); j++)
{
// wrap each local basis function as a local function.
basisFctj.setIndex(j);
// transform the local fine function to a local coarse function
const auto baseFctj = OtherGeometryWrapper(basisFctj, geometry_);
// Interpolate j^th base function by the fine basis
fineFiniteElement.localInterpolation().interpolate(baseFctj, coefficients);
// get the matrix col index
const auto& globalCoarseIndex = coarseLocalView_.index( coarseNode.localIndex( j ) );
// set the matrix indices
for (size_t i=0; i<fineLocalBasis.size(); i++)
{
// some values may be practically zero -- no need to store those
if ( std::abs(coefficients[i]) < tolerance_ )
continue;
// get the matrix row index
const auto& globalFineIndex = fineLocalView_.index( fineNode.localIndex( i ) );
if (processed_.find(globalFineIndex) != processed_.end() )
continue;
if ( exportIndices_ )
{
// add only the index of the top level
indices_.insertEntry(globalFineIndex, globalCoarseIndex);
}
else
{
// use ISTLBackend wrapper with ()-operator
matrix_(globalFineIndex,globalCoarseIndex) += coefficients[i];
}
}
}
}
private:
const CoarseLocalView& coarseLocalView_;
const FineLocalView& fineLocalView_;
const Geometry& geometry_;
const BitSet& processed_;
double tolerance_;
Matrix& matrix_;
IndexSet& indices_;
const bool exportIndices_;
};
} // namespace Impl
/** \brief Assemble the transfer matrix for multigrid methods
*
* The coarse and fine basis has to be related: in a way that the fine grid is descendant of the
* coarse grid (not necessarily direct descandant).
* The matrix is assumed to be in block structure with correct block size.
* The two basis are assumed to be scalar dune-functions basis (i.e. leaf nodes in the basis tree).
* The two bases are assumed to implement the dune-functions basis interface.
*
* \param [out] matrix the block matrix corresponding to the basis transfer operator
* \param [in] coarseBasis scalar global basis on the coarse grid
* \param [in] fineasis scalar global basis on the fine grid
* \param [in] coarseBasis global basis on the coarse grid
* \param [in] fineBasis global basis on the fine grid
* \param [in] tolerance (optional) interpolation threshold. Default is 1e-12.
*/
template<class MatrixType, class CoarseBasis, class FineBasis>
......@@ -304,28 +266,21 @@ static void assembleGlobalBasisTransferMatrix(MatrixType& matrix,
auto coarseLocalView = coarseBasis.localView();
auto fineLocalView = fineBasis.localView();
const auto rows = fineBasis.size();
const auto cols = coarseBasis.size();
if ( not fineLocalView.tree().isLeaf or not coarseLocalView.tree().isLeaf )
DUNE_THROW(Dune::Exception, "currently, we can handle power nodes only");
using FineLocalView = std::decay_t<decltype(fineLocalView)>;
using CoarseFEType = typename decltype(coarseLocalView)::Tree::FiniteElement;
using FunctionBaseClass = typename Dune::LocalFiniteElementFunctionBase<CoarseFEType>::type;
using LocalBasisWrapper = LocalBasisComponentWrapper<typename CoarseFEType::Traits::LocalBasisType, FunctionBaseClass>;
matrix.setSize(rows,cols);
matrix.setBuildMode(MatrixType::random);
// for multi index access wrap the matrix to the ISTLBackend
auto matrixBackend = Dune::Fufem::ISTLMatrixBackend(matrix);
// prepare index set with multi index access
auto indices = matrixBackend.patternBuilder();
indices.resize(fineBasis,coarseBasis);
// ///////////////////////////////////////////
// Determine the indices present in the matrix
// ///////////////////////////////////////////
// Only handle every dof once
Dune::BitSetVector<1> processed(fineBasis.size());
Dune::MatrixIndexSet indices(rows, cols);
// Only handle every fine dof once
std::unordered_set<typename FineLocalView::MultiIndex> processed;
// loop over all fine grid elements
for ( const auto& fineElement : elements(fineGridView) )
......@@ -349,69 +304,34 @@ static void assembleGlobalBasisTransferMatrix(MatrixType& matrix,
const auto& coarseElement = fE;
const auto geometryInFathers = GeometryInFathers<decltype(geometryInFathersVector)>(geometryInFathersVector);
// geometric mapping fine -> coarse
const auto geometryInFathers = GeometryInFathers(geometryInFathersVector);
coarseLocalView.bind(coarseElement);
fineLocalView.bind(fineElement);
// get the root as reference
const auto& node = fineLocalView.tree();
// visit all children of the coarse node and interpolate with the fine basis functions
::Impl::LeafNodeInterpolationVisitor visitor( coarseLocalView, fineLocalView, geometryInFathers,
processed, tolerance, matrixBackend, indices,
true /* = set indexSet */ );
// get the local basis
const auto& localBasis = node.finiteElement().localBasis();
// Hack: The RangeFieldType is the best guess of a suitable type for coefficients we have here
using CoefficientType = typename std::decay_t<decltype(localBasis)>::Traits::RangeFieldType;
std::vector<CoefficientType> coefficients(localBasis.size());
LocalBasisWrapper basisFctj(node.finiteElement().localBasis(),0);
Dune::TypeTree::applyToTreePair(coarseLocalView.tree(),fineLocalView.tree(),visitor);
const auto& coarseLocalBasis = coarseLocalView.tree().finiteElement().localBasis();
for (size_t j=0; j<coarseLocalBasis.size(); j++)
{
// wrap each local basis function as a local function.
basisFctj.setIndex(j);
// transform the local fine function to a local coarse function
const auto coarseBaseFctj = OtherGeometryWrapper<decltype(basisFctj),decltype(geometryInFathers)>(basisFctj, geometryInFathers);
// Interpolate j^th base function by the fine basis
node.finiteElement().localInterpolation().interpolate(coarseBaseFctj, coefficients);
// get the matrix col index
const auto& globalCoarseIndex = coarseLocalView.index(j);
// set the matrix indices
for (size_t i=0; i<localBasis.size(); i++) {
// some values may be practically zero -- no need to store those
if ( std::abs(coefficients[i]) < tolerance )
continue;
// get the matrix row index
const auto& globalFineIndex = fineLocalView.index(i);
if (processed[globalFineIndex][0])
continue;
indices.add(globalFineIndex, globalCoarseIndex);
}
}
// now the all dof's of this fine element to 'processed'
for (size_t i=0; i<localBasis.size(); i++)
for (size_t i=0; i<fineLocalView.size(); i++)
{
const auto& globalFineIndex = fineLocalView.index(i);
processed[globalFineIndex] = true;
const auto& globalFineIndex = fineLocalView.index( i );
processed.insert(globalFineIndex);
}
}
indices.exportIdx(matrix);
// apply the index set to the matrix
indices.setupMatrix();
matrix = 0;
// reset all dof's
processed.unsetAll();
processed.clear();
//////////////////////
// Now set the values
......@@ -432,66 +352,26 @@ static void assembleGlobalBasisTransferMatrix(MatrixType& matrix,
fE = fE.father();
}
// did we really end up in coarseElement?
if ( not coarseGridView.contains(fE) )
DUNE_THROW( Dune::Exception, "There is a fine element without a parent in the coarse grid!");
const auto& coarseElement = fE;
const auto geometryInFathers = GeometryInFathers<decltype(geometryInFathersVector)>(geometryInFathersVector);
// geometric mapping fine -> coarse
const auto geometryInFathers = GeometryInFathers(geometryInFathersVector);
coarseLocalView.bind(coarseElement);
fineLocalView.bind(fineElement);
// get the root as reference
const auto& node = fineLocalView.tree();
// get the local basis
const auto& localBasis = node.finiteElement().localBasis();
using Coefficient = typename std::decay_t<decltype(localBasis)>::Traits::RangeFieldType;
std::vector<Coefficient> coefficients(localBasis.size());
LocalBasisWrapper basisFctj(node.finiteElement().localBasis(),0);
const auto& coarseLocalBasis = coarseLocalView.tree().finiteElement().localBasis();
// visit all children of the coarse node and interpolate with the fine basis functions
::Impl::LeafNodeInterpolationVisitor visitor( coarseLocalView, fineLocalView, geometryInFathers,
processed, tolerance, matrixBackend, indices,
false /* = increase matrix entries */ );
for (size_t j=0; j<coarseLocalBasis.size(); j++)
{
// wrap each local basis function as a local function.
basisFctj.setIndex(j);
// transform the local fine function to a local coarse function
const auto coarseBaseFctj = OtherGeometryWrapper<decltype(basisFctj),decltype(geometryInFathers)>(basisFctj, geometryInFathers);
// Interpolate j^th base function by the fine basis
node.finiteElement().localInterpolation().interpolate(coarseBaseFctj, coefficients);
// get the matrix col index
const auto& globalCoarseIndex = coarseLocalView.index(j);
// set the matrix indices
for (size_t i=0; i<localBasis.size(); i++) {
// some values may be practically zero -- no need to store those
if ( std::abs(coefficients[i]) < tolerance )
continue;
Dune::TypeTree::applyToTreePair(coarseLocalView.tree(),fineLocalView.tree(),visitor);
// get the matrix row index
const auto& globalFineIndex = fineLocalView.index(i);
if (processed[globalFineIndex][0])
continue;
Dune::MatrixVector::addToDiagonal(matrix[globalFineIndex][globalCoarseIndex],coefficients[i]);
}
}
// now the all dof's of this fine element to 'processed'
for (size_t i=0; i<localBasis.size(); i++)
// now the all dof's of this fine element to 'processed'
for (size_t i=0; i<fineLocalView.size(); i++)
{
const auto& globalFineIndex = fineLocalView.index(i);
processed[globalFineIndex] = true;
const auto& globalFineIndex = fineLocalView.index( i );
processed.insert(globalFineIndex);
}
}
}
......
dune/fufem/assemblers/dunefunctionsboundaryfunctionalassembler.hh
View file @ b215e298
......@@ -36,7 +36,7 @@ public:
{
auto trialLocalView = trialBasis_.localView();
using Field = std::decay_t<decltype(vectorBackend(trialLocalView.index(0)))>;
using Field = std::decay_t<decltype(vectorBackend[trialLocalView.index(0)])>;
using LocalVector = Dune::BlockVector<Dune::FieldVector<Field,1>>;
auto localVector = LocalVector(trialLocalView.maxSize());
......@@ -47,12 +47,19 @@ public:
const auto& element = it->inside();
trialLocalView.bind(element);
for (size_t i=0; i<trialLocalView.tree().size(); ++i)
{
auto localRow = trialLocalView.tree().localIndex(i);
localVector[localRow] = 0;
}
localAssembler(it, localVector, trialLocalView);
for (size_t localRow=0; localRow<trialLocalView.size(); ++localRow)
for (size_t i=0; i<trialLocalView.tree().size(); ++i)
{
auto localRow = trialLocalView.tree().localIndex(i);
auto row = trialLocalView.index(localRow);
vectorBackend(row) += localVector[localRow];
vectorBackend[row] += localVector[localRow];
}
}
}
......
dune/fufem/assemblers/dunefunctionsfunctionalassembler.hh
View file @ b215e298
......@@ -43,11 +43,18 @@ public:
{
trialLocalView.bind(element);
for (size_t i=0; i<trialLocalView.tree().size(); ++i)
{
auto localRow = trialLocalView.tree().localIndex(i);
localVector[localRow] = 0;
}
localAssembler(element, localVector, trialLocalView);
// Add element stiffness matrix onto the global stiffness matrix
for (size_t localRow=0; localRow<trialLocalView.size(); ++localRow)
for (size_t i=0; i<trialLocalView.tree().size(); ++i)
{
auto localRow = trialLocalView.tree().localIndex(i);
auto row = trialLocalView.index(localRow);
vectorBackend[row] += localVector[localRow];
}
......
dune/fufem/assemblers/dunefunctionsoperatorassembler.hh
View file @ b215e298
......@@ -3,6 +3,8 @@
#ifndef DUNE_FUFEM_ASSEMBLERS_DUNEFUNCTIONSOPERATORASSEMBLER_HH
#define DUNE_FUFEM_ASSEMBLERS_DUNEFUNCTIONSOPERATORASSEMBLER_HH
#include <type_traits>
#include <dune/istl/matrix.hh>
#include <dune/istl/matrixindexset.hh>
......@@ -23,6 +25,61 @@ class DuneFunctionsOperatorAssembler
{
using GridView = typename TrialBasis::GridView;
// Zero-initialize local entries contained in the tree.
// In case of a subspace basis, this does not touch
// entries not contained in the subspace.
template<class TrialLocalView, class AnsatzLocalView, class LocalMatrix>
static void zeroInitializeLocal(const TrialLocalView& trialLocalView, const AnsatzLocalView& ansatzLocalView, LocalMatrix& localMatrix)
{
for (size_t i=0; i<trialLocalView.tree().size(); ++i)
{
auto localRow = trialLocalView.tree().localIndex(i);
for (size_t j=0; j<ansatzLocalView.tree().size(); ++j)
{
auto localCol = ansatzLocalView.tree().localIndex(j);
localMatrix[localRow][localCol] = 0;
}
}
}
// Distribute local matrix pattern to global pattern builder
// In case of a subspace basis, this does not touch
// entries not contained in the subspace.
template<class TrialLocalView, class AnsatzLocalView, class PatternBuilder>
static void distributeLocalPattern(const TrialLocalView& trialLocalView, const AnsatzLocalView& ansatzLocalView, PatternBuilder& patternBuilder)
{
for (size_t i=0; i<trialLocalView.tree().size(); ++i)
{
auto localRow = trialLocalView.tree().localIndex(i);
auto row = trialLocalView.index(localRow);
for (size_t j=0; j<ansatzLocalView.tree().size(); ++j)
{
auto localCol = ansatzLocalView.tree().localIndex(j);
auto col = ansatzLocalView.index(localCol);
patternBuilder.insertEntry(row,col);
}
}
}
// Distribute local matrix entries to global matrix
// In case of a subspace basis, this does not touch
// entries not contained in the subspace.
template<class TrialLocalView, class AnsatzLocalView, class LocalMatrix, class MatrixBackend>
static void distributeLocalEntries(const TrialLocalView& trialLocalView, const AnsatzLocalView& ansatzLocalView, const LocalMatrix& localMatrix, MatrixBackend& matrixBackend)
{
for (size_t i=0; i<trialLocalView.tree().size(); ++i)
{
auto localRow = trialLocalView.tree().localIndex(i);
auto row = trialLocalView.index(localRow);
for (size_t j=0; j<ansatzLocalView.tree().size(); ++j)
{
auto localCol = ansatzLocalView.tree().localIndex(j);
auto col = ansatzLocalView.index(localCol);
matrixBackend(row,col) += localMatrix[localRow][localCol];
}
}
}
public:
//! create assembler for grid
DuneFunctionsOperatorAssembler(const TrialBasis& tBasis, const AnsatzBasis& aBasis) :
......@@ -37,26 +94,17 @@ public:
{
patternBuilder.resize(trialBasis_, ansatzBasis_);
auto trialLocalView = trialBasis_.localView();
auto ansatzLocalView = ansatzBasis_.localView();
// create two localViews but use only one if bases are the same
auto ansatzLocalView = ansatzBasis_.localView();
auto separateTrialLocalView = trialBasis_.localView();
auto& trialLocalView = selectTrialLocalView(ansatzLocalView, separateTrialLocalView);
for (const auto& element : elements(trialBasis_.gridView()))
{
trialLocalView.bind(element);
ansatzLocalView.bind(element);
// bind the localViews to the element
bind(ansatzLocalView, trialLocalView, element);
// Add element stiffness matrix onto the global stiffness matrix
for (size_t i=0; i<trialLocalView.size(); ++i)
{
auto row = trialLocalView.index(i);
for (size_t j=0; j<ansatzLocalView.size(); ++j)
{
auto col = ansatzLocalView.index(j);
patternBuilder.insertEntry(row,col);
}
}
distributeLocalPattern(trialLocalView, ansatzLocalView, patternBuilder);
}
}
......@@ -80,13 +128,7 @@ public:
insideAnsatzLocalView.bind(element);
/* Coupling on the same element */
for (size_t i = 0; i < insideTrialLocalView.size(); i++) {
auto row = insideTrialLocalView.index(i);
for (size_t j = 0; j < insideAnsatzLocalView.size(); j++) {
auto col = insideAnsatzLocalView.index(j);
patternBuilder.insertEntry(row,col);
}
}
distributeLocalPattern(insideTrialLocalView, insideAnsatzLocalView, patternBuilder);
for (const auto& is : intersections(trialBasis_.gridView(), element))
{
......@@ -100,15 +142,7 @@ public:
outsideAnsatzLocalView.bind(outsideElement);
// We assume that all basis functions of the inner element couple with all basis functions from the outer one
for (size_t i=0; i<insideTrialLocalView.size(); ++i)
{
auto row = insideTrialLocalView.index(i);
for (size_t j=0; j<outsideAnsatzLocalView.size(); ++j)
{
auto col = outsideAnsatzLocalView.index(j);
patternBuilder.insertEntry(row,col);
}
}
distributeLocalPattern(insideTrialLocalView, outsideAnsatzLocalView, patternBuilder);
}
}
}
......@@ -119,35 +153,29 @@ public:
template <class MatrixBackend, class LocalAssembler>
void assembleBulkEntries(MatrixBackend&& matrixBackend, LocalAssembler&& localAssembler) const
{
auto trialLocalView = trialBasis_.localView();
auto ansatzLocalView = ansatzBasis_.localView();
// create two localViews but use only one if bases are the same
auto ansatzLocalView = ansatzBasis_.localView();
auto separateTrialLocalView = trialBasis_.localView();
auto& trialLocalView = selectTrialLocalView(ansatzLocalView, separateTrialLocalView);
using Field = std::decay_t<decltype(matrixBackend(trialLocalView.index(0), ansatzLocalView.index(0)))>;
using LocalMatrix = Dune::Matrix<Dune::FieldMatrix<Field,1,1>>;
auto localMatrix = LocalMatrix();
localMatrix.setSize(trialLocalView.maxSize(), ansatzLocalView.maxSize());
for (const auto& element : elements(trialBasis_.gridView()))
{
trialLocalView.bind(element);
ansatzLocalView.bind(element);
// bind the localViews to the element
bind(ansatzLocalView, trialLocalView, element);
localMatrix.setSize(trialLocalView.size(), ansatzLocalView.size());
zeroInitializeLocal(trialLocalView, ansatzLocalView, localMatrix);
localAssembler(element, localMatrix, trialLocalView, ansatzLocalView);
// Add element stiffness matrix onto the global stiffness matrix
for (size_t localRow=0; localRow<trialLocalView.size(); ++localRow)
{
auto row = trialLocalView.index(localRow);
for (size_t localCol=0; localCol<ansatzLocalView.size(); ++localCol)
{
auto col = ansatzLocalView.index(localCol);
matrixBackend(row,col) += localMatrix[localRow][localCol];
}
}
distributeLocalEntries(trialLocalView, ansatzLocalView, localMatrix, matrixBackend);
}
}
......@@ -178,15 +206,17 @@ public:
using MatrixContainer = Dune::Matrix<LocalMatrix>;
auto matrixContrainer = MatrixContainer(2,2);
matrixContrainer[0][0].setSize(insideTrialLocalView.maxSize(), insideAnsatzLocalView.maxSize());
matrixContrainer[0][1].setSize(insideTrialLocalView.maxSize(), outsideAnsatzLocalView.maxSize());
matrixContrainer[1][0].setSize(outsideTrialLocalView.maxSize(), insideAnsatzLocalView.maxSize());
matrixContrainer[1][1].setSize(outsideTrialLocalView.maxSize(), outsideAnsatzLocalView.maxSize());
for (const auto& element : elements(trialBasis_.gridView()))
{
insideTrialLocalView.bind(element);
insideAnsatzLocalView.bind(element);
// Resize
matrixContrainer[0][0].setSize(insideTrialLocalView.size(), insideAnsatzLocalView.size());
for (const auto& is : intersections(trialBasis_.gridView(), element))
{
......@@ -194,12 +224,13 @@ public:
if (is.boundary())
{
auto& localMatrix = matrixContrainer[0][0];
zeroInitializeLocal(insideTrialLocalView, insideAnsatzLocalView, localMatrix);
localBoundaryAssembler(is, localMatrix, insideTrialLocalView, insideAnsatzLocalView);
for (size_t i=0; i< insideTrialLocalView.size(); i++) {
for (size_t j=0; j< insideAnsatzLocalView.size(); j++)
matrixBackend(insideTrialLocalView.index(i), insideAnsatzLocalView.index(j))+=localMatrix[i][j];
}
distributeLocalEntries(insideTrialLocalView, insideAnsatzLocalView, localMatrix, matrixBackend);
}
else if (is.neighbor())
{
......@@ -216,45 +247,17 @@ public:
outsideAnsatzLocalView.bind(outsideElement);
matrixContrainer[0][1].setSize(insideTrialLocalView.size(), outsideAnsatzLocalView.size());
matrixContrainer[1][0].setSize(outsideTrialLocalView.size(), insideAnsatzLocalView.size());
matrixContrainer[1][1].setSize(outsideTrialLocalView.size(), outsideAnsatzLocalView.size());
zeroInitializeLocal(insideTrialLocalView, insideAnsatzLocalView, matrixContrainer[0][0]);
zeroInitializeLocal(insideTrialLocalView, outsideAnsatzLocalView, matrixContrainer[0][1]);
zeroInitializeLocal(outsideTrialLocalView, insideAnsatzLocalView, matrixContrainer[1][0]);
zeroInitializeLocal(outsideTrialLocalView, outsideAnsatzLocalView, matrixContrainer[1][1]);
localAssembler(is, matrixContrainer, insideTrialLocalView, insideAnsatzLocalView, outsideTrialLocalView, outsideAnsatzLocalView);
for (size_t i=0; i < insideTrialLocalView.size(); i++)
{
auto row = insideTrialLocalView.index(i);
// inside x inside
for (size_t j=0; j < insideTrialLocalView.size(); j++)
{
auto col = insideTrialLocalView.index(j);
matrixBackend(row, col) += matrixContrainer[0][0][i][j];
}
// inside x outside
for (size_t j=0; j < outsideAnsatzLocalView.size(); j++)
{
auto col = outsideAnsatzLocalView.index(j);
matrixBackend(row, col) += matrixContrainer[0][1][i][j];
}
}
for (size_t i=0; i < outsideTrialLocalView.size(); i++)
{
auto row = outsideTrialLocalView.index(i);
// outside x inside
for (size_t j=0; j < insideAnsatzLocalView.size(); j++)
{
auto col = insideAnsatzLocalView.index(j);
matrixBackend(row, col) += matrixContrainer[1][0][i][j];
}
// outside x outside
for (size_t j=0; j < outsideAnsatzLocalView.size(); j++)
{
auto col = outsideAnsatzLocalView.index(j);
matrixBackend(row, col) += matrixContrainer[1][1][i][j];
}
}
distributeLocalEntries(insideTrialLocalView, insideAnsatzLocalView, matrixContrainer[0][0], matrixBackend);
distributeLocalEntries(insideTrialLocalView, outsideAnsatzLocalView, matrixContrainer[0][1], matrixBackend);
distributeLocalEntries(outsideTrialLocalView, insideAnsatzLocalView, matrixContrainer[1][0], matrixBackend);
distributeLocalEntries(outsideTrialLocalView, outsideAnsatzLocalView, matrixContrainer[1][1], matrixBackend);
}
}
}
......@@ -267,9 +270,34 @@ public:
auto patternBuilder = matrixBackend.patternBuilder();
assembleBulkPattern(patternBuilder);
patternBuilder.setupMatrix();
matrixBackend.assign(0);
assembleBulkEntries(matrixBackend, std::forward<LocalAssembler>(localAssembler));
}
private:
//! helper function to select a trialLocalView (possibly a reference to ansatzLocalView if bases are the same)
template<class AnsatzLocalView, class TrialLocalView>
TrialLocalView& selectTrialLocalView(AnsatzLocalView& ansatzLocalView, TrialLocalView& trialLocalView) const
{
if constexpr (std::is_same<TrialBasis,AnsatzBasis>::value)
if (&trialBasis_ == &ansatzBasis_)
return ansatzLocalView;
return trialLocalView;
}
//! small helper that checks whether two localViews are the same and binds one or both to an element
template<class AnsatzLocalView, class TrialLocalView, class E>
void bind(AnsatzLocalView& ansatzLocalView, TrialLocalView& trialLocalView, const E& e) const
{
ansatzLocalView.bind(e);
if constexpr (std::is_same<TrialBasis,AnsatzBasis>::value)
if (&trialLocalView == &ansatzLocalView)
return;
// localViews differ: bind trialLocalView too
trialLocalView.bind(e);
}
protected:
const TrialBasis& trialBasis_;
const AnsatzBasis& ansatzBasis_;
......
dune/fufem/assemblers/istlbackend.hh
View file @ b215e298
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_FUNCTIONS_COMMON_ISTL_BACKEND_HH
#define DUNE_FUNCTIONS_COMMON_ISTL_BACKEND_HH
#ifndef DUNE_FUFEM_ASSEMBLERS_ISTLBACKEND_HH
#define DUNE_FUFEM_ASSEMBLERS_ISTLBACKEND_HH
#include <dune/common/indices.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/hybridutilities.hh>
#include <dune/istl/matrix.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/matrixindexset.hh>
#include <dune/istl/multitypeblockvector.hh>
#include <dune/istl/multitypeblockmatrix.hh>
#include <dune/functions/common/indexaccess.hh>
#include <dune/functions/functionspacebases/hierarchicvectorwrapper.hh>
......@@ -134,11 +137,132 @@ private:
};
/**
* \brief Helper class for building matrix pattern
*
* Version for MultiTypeBlockMatrix with BCRSMatrix-like blocks.
* It uses a 2D std::array of MatrixIndexSets
*/
template<class Row0, class... Rows>
class MatrixBuilder<MultiTypeBlockMatrix<Row0,Rows...>>
{
public:
using Matrix = MultiTypeBlockMatrix<Row0,Rows...>;
MatrixBuilder(Matrix& matrix) :
matrix_(matrix)
{}
template<class RowSizeInfo, class ColSizeInfo>
void resize(const RowSizeInfo& rowSizeInfo, const ColSizeInfo& colSizeInfo)
{
for(size_t i=0; i<rowBlocks_; i++)
for(size_t j=0; j<colBlocks_; j++)
indices_[i][j].resize(rowSizeInfo.size({i}), colSizeInfo.size({j}));
}
template<class RowIndex, class ColIndex>
void insertEntry(const RowIndex& rowIndex, const ColIndex& colIndex)
{
indices_[rowIndex[0]][colIndex[0]].add(rowIndex[1], colIndex[1]);
}
void setupMatrix()
{
Hybrid::forEach(Hybrid::integralRange(std::integral_constant<std::size_t,rowBlocks_>()), [&](auto&& i) {
Hybrid::forEach(Hybrid::integralRange(std::integral_constant<std::size_t,colBlocks_>()), [&](auto&& j) {
indices_[i][j].exportIdx(matrix_[i][j]);
});
});
}
private:
//! number of multi type rows
static constexpr size_t rowBlocks_ = Matrix::N();
//! number of multi type cols (we assume the matrix is well-formed)
static constexpr size_t colBlocks_ = Matrix::M();
//! 2D array of IndexSets
std::array<std::array<Dune::MatrixIndexSet,colBlocks_>,rowBlocks_> indices_;
Matrix& matrix_;
};
template<class T, class A>
class MatrixBuilder<Dune::Matrix<Dune::BCRSMatrix<T, A>>>
{
public:
using Matrix = Dune::Matrix<Dune::BCRSMatrix<T, A>>;
MatrixBuilder(Matrix& matrix) :
matrix_(matrix)
{}
template<class RowSizeInfo, class ColSizeInfo>
void resize(const RowSizeInfo& rowSizeInfo, const ColSizeInfo& colSizeInfo)
{
auto rowBlocks = rowSizeInfo.size();
auto colBlocks = colSizeInfo.size();
indices_.resize(rowBlocks);
for(size_t i=0; i<rowBlocks; i++)
{
indices_[i].resize(colBlocks);
for(size_t j=0; j<colBlocks; j++)
indices_[i][j].resize(rowSizeInfo.size({i}), colSizeInfo.size({j}));
}
}
template<class RowIndex, class ColIndex>
void insertEntry(const RowIndex& rowIndex, const ColIndex& colIndex)
{
indices_[rowIndex[0]][colIndex[0]].add(rowIndex[1], colIndex[1]);
}
void setupMatrix()
{
auto rowBlocks = indices_.size();
auto colBlocks = indices_[0].size();
matrix_.setSize(rowBlocks, colBlocks);
for(size_t i=0; i<rowBlocks; i++)
for(size_t j=0; j<colBlocks; j++)
indices_[i][j].exportIdx(matrix_[i][j]);
}
private:
std::vector<std::vector<Dune::MatrixIndexSet>> indices_;
Matrix& matrix_;
};
template<class M, class E=typename M::field_type>
class ISTLMatrixBackend
{
// The following is essentially a copy of Dune::Functions::hybridIndexAccess
// But since matrices do not provide size(), we have to use N() instead.
template<class C, class I, class F,
typename std::enable_if< Dune::models<Dune::Functions::Imp::Concept::HasDynamicIndexAccess<I>, C>(), int>::type = 0>
static auto hybridRowIndexAccess(C&& c, const I& i, F&& f)
-> decltype(f(c[i]))
{
return f(c[i]);
}
template<class C, class I, class F,
typename std::enable_if< not Dune::models<Dune::Functions::Imp::Concept::HasDynamicIndexAccess<I>, C>(), int>::type = 0>
static decltype(auto) hybridRowIndexAccess(C&& c, const I& i, F&& f)
{
return Hybrid::switchCases(std::make_index_sequence<std::decay_t<C>::N()>(), i,
[&](const auto& ii) -> decltype(auto){
return f(c[ii]);
}, [&]() -> decltype(auto){
return f(c[Dune::Indices::_0]);
});
}
template<class Result, class RowIndex, class ColIndex>
struct MatrixMultiIndexResolver
{
......@@ -178,7 +302,7 @@ class ISTLMatrixBackend
auto&& subRowIndex = rowIndex_;
auto&& subColIndex = Dune::Functions::Imp::shiftedDynamicMultiIndex<1>(colIndex_);
auto&& subIndexResolver = MatrixMultiIndexResolver<Result, decltype(subRowIndex), decltype(subColIndex)>(subRowIndex, subColIndex);
return static_cast<Result>((hybridIndexAccess(c, _0, [&](auto&& ci) -> decltype(auto) {
return static_cast<Result>((hybridRowIndexAccess(c, _0, [&](auto&& ci) -> decltype(auto) {
return static_cast<Result>((hybridIndexAccess(ci, colIndex_[_0], subIndexResolver)));
})));
}
......@@ -193,7 +317,7 @@ class ISTLMatrixBackend
auto&& subRowIndex = Dune::Functions::Imp::shiftedDynamicMultiIndex<1>(rowIndex_);
auto&& subColIndex = colIndex_;
auto&& subIndexResolver = MatrixMultiIndexResolver<Result, decltype(subRowIndex), decltype(subColIndex)>(subRowIndex, subColIndex);
return static_cast<Result>((hybridIndexAccess(c, rowIndex_[_0], [&](auto&& ci) -> decltype(auto) {
return static_cast<Result>((hybridRowIndexAccess(c, rowIndex_[_0], [&](auto&& ci) -> decltype(auto) {
return static_cast<Result>((hybridIndexAccess(ci, _0, subIndexResolver)));
})));
}
......@@ -209,7 +333,7 @@ class ISTLMatrixBackend
auto&& subRowIndex = Dune::Functions::Imp::shiftedDynamicMultiIndex<1>(rowIndex_);
auto&& subColIndex = Dune::Functions::Imp::shiftedDynamicMultiIndex<1>(colIndex_);
auto&& subIndexResolver = MatrixMultiIndexResolver<Result, decltype(subRowIndex), decltype(subColIndex)>(subRowIndex, subColIndex);
return static_cast<Result>((hybridIndexAccess(c, rowIndex_[_0], [&](auto&& ci) -> decltype(auto) {
return static_cast<Result>((hybridRowIndexAccess(c, rowIndex_[_0], [&](auto&& ci) -> decltype(auto) {
return static_cast<Result>((hybridIndexAccess(ci, colIndex_[_0], subIndexResolver)));
})));
}
......@@ -256,6 +380,31 @@ public:
return i(*matrix_);
}
/**
* \brief Const access to wrapped matrix
*/
const Matrix& matrix() const
{
return *matrix_;
}
/**
* \brief Mutable access to wrapped matrix
*/
Matrix& matrix()
{
return *matrix_;
}
/**
* \brief Assign value to wrapped matrix
*/
template<class Value>
void assign(const Value& value)
{
matrix() = value;
}
protected:
Matrix* matrix_;
......@@ -277,23 +426,9 @@ auto istlMatrixBackend(Matrix& matrix)
template<class Entry, class Vector>
auto istlVectorBackend(Vector& vector)
{
return Dune::Functions::HierarchicVectorWrapper<Vector, Entry>(vector);
}
template<class Vector>
auto istlVectorBackend(Vector& vector)
{
return Dune::Functions::HierarchicVectorWrapper<Vector, typename Vector::field_type>(vector);
}
} // namespace Dune::Fufem
} // namespace Dune
#endif // DUNE_FUNCTIONS_COMMON_ISTL_BACKEND_HH
#endif // DUNE_FUFEM_ASSEMBLERS_ISTLBACKEND_HH
dune/fufem/assemblers/localassembler.hh
View file @ b215e298
......@@ -22,7 +22,7 @@ class LocalAssembler
public:
typedef typename GridType::template Codim<0>::Entity::Entity Element;
typedef typename GridType::template Codim<0>::Entity Element;
typedef Dune::Matrix< Dune::FieldMatrix<int,1,1> > BoolMatrix;
typedef typename Dune::Matrix<T> LocalMatrix;
typedef typename T::field_type field_type;
......
dune/fufem/assemblers/localassemblers/CMakeLists.txt
View file @ b215e298
......@@ -2,7 +2,6 @@ install(FILES
adolchessianassembler.hh
adolclinearizationassembler.hh
adolclocalenergy.hh
boundarymassassembler.hh
convolutionassembler.hh
dunefunctionsl2functionalassembler.hh
generalizedboundarymassassembler.hh
......