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

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
Show changes
Commits on Source (362)
262 additional commits have been omitted to prevent performance issues.
Hide whitespace changes
Inline Side-by-side
Showing
with 583 additions and 722 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 @@ ...@@ -2,13 +2,31 @@
# Install external dependencies # Install external dependencies
before_script: before_script:
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-alugrid.git - 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 - duneci-install-module https://git.imp.fu-berlin.de/agnumpde/dune-matrix-vector.git
dune:git--clang: dune:2.8 debian-11 gcc-10 C++20:
image: duneci/dune:git variables:
script: duneci-standard-test --opts=/duneci/opts.clang 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: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--gcc: dune:git debian-11 gcc-10 C++20:
image: duneci/dune:git image: registry.dune-project.org/docker/ci/dune:git-debian-11-gcc-10-20
script: duneci-standard-test script: duneci-standard-test
CHANGELOG.md 0 → 100644
View file @ b215e298
# 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
- 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) ...@@ -11,6 +11,7 @@ include(DuneMacros)
dune_project() dune_project()
add_subdirectory("cmake")
add_subdirectory("doc") add_subdirectory("doc")
add_subdirectory("dune") 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 @@ ...@@ -11,7 +11,7 @@
function(add_dune_adolc_flags _targets) function(add_dune_adolc_flags _targets)
if(ADOLC_FOUND) if(ADOLC_FOUND)
foreach(_target ${_targets}) foreach(_target ${_targets})
target_link_libraries(${_target} ${ADOLC_LIBRARIES}) target_link_libraries(${_target} PUBLIC ${ADOLC_LIBRARIES})
get_target_property(_props ${_target} COMPILE_FLAGS) get_target_property(_props ${_target} COMPILE_FLAGS)
string(REPLACE "_props-NOTFOUND" "" _props "${_props}") string(REPLACE "_props-NOTFOUND" "" _props "${_props}")
set_target_properties(${_target} PROPERTIES COMPILE_FLAGS set_target_properties(${_target} PROPERTIES COMPILE_FLAGS
......
cmake/modules/AddHDF5Flags.cmake
View file @ b215e298
...@@ -11,7 +11,7 @@ ...@@ -11,7 +11,7 @@
function(add_dune_hdf5_flags _targets) function(add_dune_hdf5_flags _targets)
if(HDF5_FOUND) if(HDF5_FOUND)
foreach(_target ${_targets}) 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}) set_property(TARGET ${_target} APPEND PROPERTY INCLUDE_DIRECTORIES ${HDF5_INCLUDE_DIRS})
endforeach(_target ${_targets}) endforeach(_target ${_targets})
endif(HDF5_FOUND) endif(HDF5_FOUND)
......
cmake/modules/AddPythonLibsFlags.cmake
View file @ b215e298
...@@ -9,15 +9,16 @@ ...@@ -9,15 +9,16 @@
# #
function(add_dune_pythonlibs_flags _targets) function(add_dune_pythonlibs_flags _targets)
if(PYTHONLIBS_FOUND) foreach(_target ${_targets})
foreach(_target ${_targets}) if(Python3_FOUND)
target_link_libraries(${_target} ${PYTHON_LIBRARIES}) target_link_libraries(${_target} PUBLIC ${Python3_LIBRARIES})
# target_include_directories(${_target} PRIVATE ${PYTHON_INCLUDE_DIRS}) set_property(TARGET ${_target} APPEND PROPERTY INCLUDE_DIRECTORIES ${Python3_INCLUDE_DIRS})
# target_compile_definitions(${_target} PRIVATE "-DHAVE_PYTHON") # backward compatibility with 2.7 core tests
# target_compile_options(${_target} PRIVATE "-fno-strict-aliasing") elseif(PYTHONLIBS_FOUND)
target_link_libraries(${_target} PUBLIC ${PYTHON_LIBRARIES})
set_property(TARGET ${_target} APPEND PROPERTY INCLUDE_DIRECTORIES ${PYTHON_INCLUDE_DIRS}) 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_DEFINITIONS "HAVE_PYTHON")
set_property(TARGET ${_target} APPEND PROPERTY COMPILE_OPTIONS "-fno-strict-aliasing") set_property(TARGET ${_target} APPEND PROPERTY COMPILE_OPTIONS "-fno-strict-aliasing")
endforeach(_target ${_targets}) endforeach(_target ${_targets})
endif(PYTHONLIBS_FOUND)
endfunction(add_dune_pythonlibs_flags) 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}) 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) ...@@ -36,11 +36,12 @@ macro(dune_fufem_add_copy_dependency target file_name)
add_dependencies(${target} "${target}_copy_${file_name}") add_dependencies(${target} "${target}_copy_${file_name}")
endmacro(dune_fufem_add_copy_dependency) 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(Adolc)
find_package(HDF5 COMPONENTS C HL) find_package(HDF5 COMPONENTS C HL)
......
cmake/modules/FindAdolc.cmake
View file @ b215e298
...@@ -4,6 +4,10 @@ ...@@ -4,6 +4,10 @@
# Add the path to your locally installed library to the variable # Add the path to your locally installed library to the variable
# CMAKE_PREFIX_PATH # CMAKE_PREFIX_PATH
# #
# Uses the following variables:
#
# ADOLC_ROOT Path to ADOL-C installation
#
# Sets the following variables: # Sets the following variables:
# #
# ADOLC_FOUND True if ADOL-C is available and usable. # ADOLC_FOUND True if ADOL-C is available and usable.
...@@ -13,12 +17,22 @@ ...@@ -13,12 +17,22 @@
find_path(ADOLC_INCLUDE_DIR find_path(ADOLC_INCLUDE_DIR
NAMES "adolc/adouble.h" NAMES "adolc/adouble.h"
PATHS ${ADOLC_ROOT}
PATH_SUFFIXES "include" PATH_SUFFIXES "include"
NO_DEFAULT_PATH
)
find_path(ADOLC_INCLUDE_DIR
NAMES "adolc/adouble.h"
) )
find_library(ADOLC_LIBRARY find_library(ADOLC_LIBRARY
NAMES adolc NAMES adolc
PATHS ${ADOLC_ROOT}
PATH_SUFFIXES "lib" "lib64" PATH_SUFFIXES "lib" "lib64"
NO_DEFAULT_PATH
)
find_library(ADOLC_LIBRARY
NAMES adolc
) )
include(FindPackageHandleStandardArgs) include(FindPackageHandleStandardArgs)
......
config.h.cmake
View file @ b215e298
...@@ -35,7 +35,7 @@ ...@@ -35,7 +35,7 @@
#define DUNE_FUFEM_VERSION_MAJOR @DUNE_FUFEM_VERSION_MAJOR@ #define DUNE_FUFEM_VERSION_MAJOR @DUNE_FUFEM_VERSION_MAJOR@
/* Define to the minor version of dune-fufem */ /* Define to the minor version of dune-fufem */
#define DUNE_FUFEM_VERSION_MINOR @DUNE_FUFE;_VERSION_MINOR@ #define DUNE_FUFEM_VERSION_MINOR @DUNE_FUFEM_VERSION_MINOR@
/* Define to the revision of dune-fufem */ /* Define to the revision of dune-fufem */
#define DUNE_FUFEM_VERSION_REVISION @DUNE_FUFEM_VERSION_REVISION@ #define DUNE_FUFEM_VERSION_REVISION @DUNE_FUFEM_VERSION_REVISION@
......
dune.module
View file @ b215e298
Module: dune-fufem Module: dune-fufem
Depends: dune-common (>= 2.5) dune-geometry (>= 2.5) dune-grid (>= 2.5) dune-istl (>= 2.5) dune-localfunctions (>= 2.5) dune-functions dune-matrix-vector 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.6-git Version: 2.10-git
Maintainer: Carsten Gräser (graeser@math.fu-berlin.de) Maintainer: Carsten Gräser (graeser@math.fau.de)
Suggests: dune-alugrid dune-subgrid Suggests: dune-alugrid dune-subgrid
dune/fufem/CMakeLists.txt
View file @ b215e298
...@@ -13,7 +13,6 @@ add_subdirectory(utilities) ...@@ -13,7 +13,6 @@ add_subdirectory(utilities)
add_subdirectory(test) add_subdirectory(test)
install(FILES install(FILES
any.hh
arcofcircle.hh arcofcircle.hh
arithmetic.hh arithmetic.hh
basistraceindexset.hh basistraceindexset.hh
...@@ -52,5 +51,4 @@ install(FILES ...@@ -52,5 +51,4 @@ install(FILES
staticmatrixtools.hh staticmatrixtools.hh
surfmassmatrix.hh surfmassmatrix.hh
symmetrictensor.hh symmetrictensor.hh
ulis_tools.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/fufem) DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/fufem)
dune/fufem/any.hh deleted 100644 → 0
View file @ 85bcbd04
#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/assembler.hh
View file @ b215e298
...@@ -19,7 +19,7 @@ class Assembler ...@@ -19,7 +19,7 @@ class Assembler
{} {}
template <class LocalAssemblerType, class GlobalMatrixType> template <class LocalAssemblerType, class GlobalMatrixType>
void assembleOperator(LocalAssemblerType& localAssembler, GlobalMatrixType& A, const bool lumping=false) const void assembleOperator(LocalAssemblerType&& localAssembler, GlobalMatrixType& A, const bool lumping=false) const
{ {
const OperatorAssembler<TrialBasis, AnsatzBasis> operatorAssembler(tBasis_, aBasis_); const OperatorAssembler<TrialBasis, AnsatzBasis> operatorAssembler(tBasis_, aBasis_);
operatorAssembler.assemble(localAssembler, A, lumping); operatorAssembler.assemble(localAssembler, A, lumping);
...@@ -47,5 +47,10 @@ class Assembler ...@@ -47,5 +47,10 @@ class Assembler
const AnsatzBasis& aBasis_; const AnsatzBasis& aBasis_;
}; };
template <class TrialBasis, class AnsatzBasis>
auto assembler(const TrialBasis& tBasis, const AnsatzBasis& aBasis) {
return Assembler<TrialBasis, AnsatzBasis>(tBasis, aBasis);
}
#endif #endif
dune/fufem/assemblers/basisinterpolationmatrixassembler.hh
View file @ b215e298
...@@ -4,6 +4,9 @@ ...@@ -4,6 +4,9 @@
#define DUNE_FUFEM_ASSEMBLERS_BASIS_INTERPOLATION_MATRIX_ASSEMBLER_HH #define DUNE_FUFEM_ASSEMBLERS_BASIS_INTERPOLATION_MATRIX_ASSEMBLER_HH
#include <type_traits> #include <type_traits>
#include <typeinfo>
#include <unordered_set>
#include <dune/common/bitsetvector.hh> #include <dune/common/bitsetvector.hh>
#include <dune/istl/matrixindexset.hh> #include <dune/istl/matrixindexset.hh>
...@@ -12,7 +15,13 @@ ...@@ -12,7 +15,13 @@
#include <dune/matrix-vector/addtodiagonal.hh> #include <dune/matrix-vector/addtodiagonal.hh>
#include <dune/fufem/functions/cachedcomponentwrapper.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. */ /** \brief Wrapper that extracts a single local basis function. */
template<class LocalBasis, class Base> template<class LocalBasis, class Base>
...@@ -46,173 +55,326 @@ class LocalBasisComponentWrapper : ...@@ -46,173 +55,326 @@ class LocalBasisComponentWrapper :
const LocalBasis& localBasis_; 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; /** \brief Wrapper for evaluating local functions in another geometry */
typedef typename BasisType1::LocalFiniteElement FEType1; template<class Function, class Geometry>
typedef typename Dune::LocalFiniteElementFunctionBase<FEType0>::type FunctionBaseClass; class OtherGeometryWrapper
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>();
// /////////////////////////////////////////// public:
// Determine the indices present in the matrix
// /////////////////////////////////////////////////
// Only handle every dof once OtherGeometryWrapper(const Function& function, const Geometry& geometry)
Dune::BitSetVector<1> processed(fineBasis.size()); : function_(function)
for (size_t i=0; i<processed.size();i++) , geometry_(geometry)
if (fineBasis.isConstrained(i)) {}
processed[i] = true;
Dune::MatrixIndexSet indices(rows, cols); // forward the evaluation into the other geometry
template<class X>
auto operator()(const X& x) const
{
return function_(geometry_.global(x));
}
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(); private:
const size_t numBaseFct1 = lfe1.localBasis().size();
// check if all components have already been processed const Function function_;
bool allProcessed = true; const Geometry geometry_;
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 /** \brief Wrapper chaining sequencial father geometries */
LocalBasisWrapper basisFctj(lfe0.localBasis(),0); template<class GeometryContainer>
class GeometryInFathers
{
for (size_t j=0; j<numBaseFct0; j++) public:
{
// wrap each local basis function as a local function.
basisFctj.setIndex(j);
// Interpolate j^th base function by the fine basis GeometryInFathers(const GeometryContainer& geometryContainer)
lfe1.localInterpolation().interpolate(basisFctj, values); : geometryContainer_(geometryContainer)
{}
int globalCoarse = coarseBasis.index(*eIt,j); // forward the global() method to all single geometries
template<class X>
auto global(const X& x) const
{
auto y = x;
for (const auto& g : geometryContainer_)
{
y = g.global(y);
}
return y;
}
for (size_t i=0; i<numBaseFct1; i++) { private:
if (std::abs(values[i])<1e-12) const GeometryContainer geometryContainer_;
continue; };
int globalFine = fineBasis.index(*eIt,i);
if (processed[globalFine][0])
continue;
if (coarseBasis.isConstrained(globalCoarse)) { namespace Impl{
const auto& lin = coarseBasis.constraints(globalCoarse);
for (size_t k=0; k<lin.size(); k++) /** \brief Visitor for interpolating the local basis of a leaf node by a fine local basis
indices.add(globalFine, lin[k].index); *
} else * The LeafNodeInterpolationVisitor traverses the local tree of the coarse basis at an
indices.add(globalFine, globalCoarse); * 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
for (size_t i=0; i<numBaseFct1; i++) {
processed[fineBasis.index(*eIt,i)]=true; // 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_;
};
indices.exportIdx(matrix); } // namespace Impl
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(); /** \brief Assemble the transfer matrix for multigrid methods
const size_t numBaseFct1 = lfe1.localBasis().size(); *
* 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 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 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>
static void assembleGlobalBasisTransferMatrix(MatrixType& matrix,
const CoarseBasis& coarseBasis,
const FineBasis& fineBasis,
const double tolerance = 1e-12)
{
const auto& coarseGridView = coarseBasis.gridView();
const auto& fineGridView = fineBasis.gridView();
auto coarseLocalView = coarseBasis.localView();
auto fineLocalView = fineBasis.localView();
using FineLocalView = std::decay_t<decltype(fineLocalView)>;
// 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 fine dof once
std::unordered_set<typename FineLocalView::MultiIndex> processed;
// loop over all fine grid elements
for ( const auto& fineElement : elements(fineGridView) )
{
// find a coarse element that is the parent of the fine element
// start in the fine grid and track the geometry mapping
auto fE = fineElement;
// collect all geometryInFather's on the way
std::vector<decltype(fineElement.geometryInFather())> geometryInFathersVector;
while ( not coarseGridView.contains(fE) and fE.hasFather() )
{
// add the geometry to the container
geometryInFathersVector.push_back(fE.geometryInFather());
// step up one level
fE = fE.father();
}
// check if all components have already been processed // did we really end up in coarseElement?
bool allProcessed = true; if ( not coarseGridView.contains(fE) )
for(size_t i=0; i<numBaseFct1; ++i) DUNE_THROW( Dune::Exception, "There is a fine element without a parent in the coarse grid!");
allProcessed = allProcessed and processed[fineBasis.index(*eIt, i)][0];
if (not allProcessed) { const auto& coarseElement = fE;
// preallocate vector for function evaluations // geometric mapping fine -> coarse
std::vector<typename FEType0::Traits::LocalBasisType::Traits::RangeType> values(numBaseFct0); const auto geometryInFathers = GeometryInFathers(geometryInFathersVector);
// Extract single basis functions into a format that can be used within local interpolation coarseLocalView.bind(coarseElement);
LocalBasisWrapper basisFctj(lfe0.localBasis(),0); fineLocalView.bind(fineElement);
for (size_t j=0; j<numBaseFct0; j++) // visit all children of the coarse node and interpolate with the fine basis functions
{ ::Impl::LeafNodeInterpolationVisitor visitor( coarseLocalView, fineLocalView, geometryInFathers,
// wrap each local basis function as a local function. processed, tolerance, matrixBackend, indices,
basisFctj.setIndex(j); true /* = set indexSet */ );
int globalCoarse = coarseBasis.index(*eIt,j); Dune::TypeTree::applyToTreePair(coarseLocalView.tree(),fineLocalView.tree(),visitor);
// Interpolate j^th base function by the fine basis
lfe1.localInterpolation().interpolate(basisFctj, values);
for (size_t i=0; i<numBaseFct1; 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.insert(globalFineIndex);
}
}
// apply the index set to the matrix
indices.setupMatrix();
matrix = 0;
// reset all dof's
processed.clear();
//////////////////////
// Now set the values
//////////////////////
for ( const auto& fineElement : elements(fineGridView) )
{
// find a coarse element that is the parent of the fine element
// start in the fine grid and track the geometry mapping
auto fE = fineElement;
// collect all geometryInFather's on the way
std::vector<decltype(fineElement.geometryInFather())> geometryInFathersVector;
while ( not coarseGridView.contains(fE) and fE.hasFather() )
{
// add the geometry to the container
geometryInFathersVector.push_back(fE.geometryInFather());
// step up one level
fE = fE.father();
}
if (std::abs(values[i])<1e-12) const auto& coarseElement = fE;
continue;
int globalFine = fineBasis.index(*eIt,i); // geometric mapping fine -> coarse
if (processed[globalFine][0]) const auto geometryInFathers = GeometryInFathers(geometryInFathersVector);
continue;
if (coarseBasis.isConstrained(globalCoarse)) { coarseLocalView.bind(coarseElement);
fineLocalView.bind(fineElement);
const auto& lin = coarseBasis.constraints(globalCoarse); // 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 k=0; k<lin.size(); k++) Dune::TypeTree::applyToTreePair(coarseLocalView.tree(),fineLocalView.tree(),visitor);
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++) // now the all dof's of this fine element to 'processed'
processed[fineBasis.index(*eIt,i)] = true; for (size_t i=0; i<fineLocalView.size(); i++)
} {
const auto& globalFineIndex = fineLocalView.index( i );
processed.insert(globalFineIndex);
} }
}
} }
#endif #endif
dune/fufem/assemblers/dunefunctionsboundaryfunctionalassembler.hh
View file @ b215e298
...@@ -35,9 +35,8 @@ public: ...@@ -35,9 +35,8 @@ public:
void assembleBulkEntries(VectorBackend&& vectorBackend, LocalAssembler&& localAssembler) const void assembleBulkEntries(VectorBackend&& vectorBackend, LocalAssembler&& localAssembler) const
{ {
auto trialLocalView = trialBasis_.localView(); auto trialLocalView = trialBasis_.localView();
auto trialLocalIndexSet = trialBasis_.localIndexSet();
using Field = std::decay_t<decltype(vectorBackend(trialLocalIndexSet.index(0)))>; using Field = std::decay_t<decltype(vectorBackend[trialLocalView.index(0)])>;
using LocalVector = Dune::BlockVector<Dune::FieldVector<Field,1>>; using LocalVector = Dune::BlockVector<Dune::FieldVector<Field,1>>;
auto localVector = LocalVector(trialLocalView.maxSize()); auto localVector = LocalVector(trialLocalView.maxSize());
...@@ -47,14 +46,20 @@ public: ...@@ -47,14 +46,20 @@ public:
{ {
const auto& element = it->inside(); const auto& element = it->inside();
trialLocalView.bind(element); trialLocalView.bind(element);
trialLocalIndexSet.bind(trialLocalView);
for (size_t i=0; i<trialLocalView.tree().size(); ++i)
{
auto localRow = trialLocalView.tree().localIndex(i);
localVector[localRow] = 0;
}
localAssembler(it, localVector, trialLocalView); localAssembler(it, localVector, trialLocalView);
for (size_t localRow=0; localRow<trialLocalIndexSet.size(); ++localRow) for (size_t i=0; i<trialLocalView.tree().size(); ++i)
{ {
auto row = trialLocalIndexSet.index(localRow); auto localRow = trialLocalView.tree().localIndex(i);
vectorBackend(row) += localVector[localRow]; auto row = trialLocalView.index(localRow);
vectorBackend[row] += localVector[localRow];
} }
} }
} }
......
dune/fufem/assemblers/dunefunctionsfunctionalassembler.hh
View file @ b215e298
...@@ -33,9 +33,8 @@ public: ...@@ -33,9 +33,8 @@ public:
void assembleBulkEntries(VectorBackend&& vectorBackend, LocalAssembler&& localAssembler) const void assembleBulkEntries(VectorBackend&& vectorBackend, LocalAssembler&& localAssembler) const
{ {
auto trialLocalView = trialBasis_.localView(); auto trialLocalView = trialBasis_.localView();
auto trialLocalIndexSet = trialBasis_.localIndexSet();
using Field = std::decay_t<decltype(vectorBackend(trialLocalIndexSet.index(0)))>; using Field = std::decay_t<decltype(vectorBackend[trialLocalView.index(0)])>;
using LocalVector = Dune::BlockVector<Dune::FieldVector<Field,1>>; using LocalVector = Dune::BlockVector<Dune::FieldVector<Field,1>>;
auto localVector = LocalVector(trialLocalView.maxSize()); auto localVector = LocalVector(trialLocalView.maxSize());
...@@ -43,15 +42,21 @@ public: ...@@ -43,15 +42,21 @@ public:
for (const auto& element : elements(trialBasis_.gridView())) for (const auto& element : elements(trialBasis_.gridView()))
{ {
trialLocalView.bind(element); trialLocalView.bind(element);
trialLocalIndexSet.bind(trialLocalView);
for (size_t i=0; i<trialLocalView.tree().size(); ++i)
{
auto localRow = trialLocalView.tree().localIndex(i);
localVector[localRow] = 0;
}
localAssembler(element, localVector, trialLocalView); localAssembler(element, localVector, trialLocalView);
// Add element stiffness matrix onto the global stiffness matrix // Add element stiffness matrix onto the global stiffness matrix
for (size_t localRow=0; localRow<trialLocalIndexSet.size(); ++localRow) for (size_t i=0; i<trialLocalView.tree().size(); ++i)
{ {
auto row = trialLocalIndexSet.index(localRow); auto localRow = trialLocalView.tree().localIndex(i);
vectorBackend(row) += localVector[localRow]; auto row = trialLocalView.index(localRow);
vectorBackend[row] += localVector[localRow];
} }
} }
} }
......
dune/fufem/assemblers/dunefunctionsoperatorassembler.hh
View file @ b215e298
...@@ -3,6 +3,8 @@ ...@@ -3,6 +3,8 @@
#ifndef DUNE_FUFEM_ASSEMBLERS_DUNEFUNCTIONSOPERATORASSEMBLER_HH #ifndef DUNE_FUFEM_ASSEMBLERS_DUNEFUNCTIONSOPERATORASSEMBLER_HH
#define DUNE_FUFEM_ASSEMBLERS_DUNEFUNCTIONSOPERATORASSEMBLER_HH #define DUNE_FUFEM_ASSEMBLERS_DUNEFUNCTIONSOPERATORASSEMBLER_HH
#include <type_traits>
#include <dune/istl/matrix.hh> #include <dune/istl/matrix.hh>
#include <dune/istl/matrixindexset.hh> #include <dune/istl/matrixindexset.hh>
...@@ -23,6 +25,61 @@ class DuneFunctionsOperatorAssembler ...@@ -23,6 +25,61 @@ class DuneFunctionsOperatorAssembler
{ {
using GridView = typename TrialBasis::GridView; 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: public:
//! create assembler for grid //! create assembler for grid
DuneFunctionsOperatorAssembler(const TrialBasis& tBasis, const AnsatzBasis& aBasis) : DuneFunctionsOperatorAssembler(const TrialBasis& tBasis, const AnsatzBasis& aBasis) :
...@@ -37,30 +94,17 @@ public: ...@@ -37,30 +94,17 @@ public:
{ {
patternBuilder.resize(trialBasis_, ansatzBasis_); patternBuilder.resize(trialBasis_, ansatzBasis_);
auto trialLocalView = trialBasis_.localView(); // create two localViews but use only one if bases are the same
auto trialLocalIndexSet = trialBasis_.localIndexSet(); auto ansatzLocalView = ansatzBasis_.localView();
auto separateTrialLocalView = trialBasis_.localView();
auto ansatzLocalView = ansatzBasis_.localView(); auto& trialLocalView = selectTrialLocalView(ansatzLocalView, separateTrialLocalView);
auto ansatzLocalIndexSet = ansatzBasis_.localIndexSet();
for (const auto& element : elements(trialBasis_.gridView())) for (const auto& element : elements(trialBasis_.gridView()))
{ {
trialLocalView.bind(element); // bind the localViews to the element
trialLocalIndexSet.bind(trialLocalView); bind(ansatzLocalView, trialLocalView, element);
ansatzLocalView.bind(element); distributeLocalPattern(trialLocalView, ansatzLocalView, patternBuilder);
ansatzLocalIndexSet.bind(ansatzLocalView);
// Add element stiffness matrix onto the global stiffness matrix
for (size_t i=0; i<trialLocalIndexSet.size(); ++i)
{
auto row = trialLocalIndexSet.index(i);
for (size_t j=0; j<ansatzLocalIndexSet.size(); ++j)
{
auto col = ansatzLocalIndexSet.index(j);
patternBuilder.insertEntry(row,col);
}
}
} }
} }
...@@ -72,59 +116,33 @@ public: ...@@ -72,59 +116,33 @@ public:
patternBuilder.resize(trialBasis_, ansatzBasis_); patternBuilder.resize(trialBasis_, ansatzBasis_);
auto insideTrialLocalView = trialBasis_.localView(); auto insideTrialLocalView = trialBasis_.localView();
auto insideTrialLocalIndexSet = trialBasis_.localIndexSet();
auto insideAnsatzLocalView = ansatzBasis_.localView(); auto insideAnsatzLocalView = ansatzBasis_.localView();
auto insideAnsatzLocalIndexSet = ansatzBasis_.localIndexSet();
auto outsideTrialLocalView = trialBasis_.localView();
auto outsideTrialLocalIndexSet = trialBasis_.localIndexSet();
auto outsideAnsatzLocalView = ansatzBasis_.localView(); auto outsideAnsatzLocalView = ansatzBasis_.localView();
auto outsideAnsatzLocalIndexSet = ansatzBasis_.localIndexSet();
for (const auto& element : elements(trialBasis_.gridView())) for (const auto& element : elements(trialBasis_.gridView()))
{ {
insideTrialLocalView.bind(element); insideTrialLocalView.bind(element);
insideTrialLocalIndexSet.bind(insideTrialLocalView);
insideAnsatzLocalView.bind(element); insideAnsatzLocalView.bind(element);
insideAnsatzLocalIndexSet.bind(insideAnsatzLocalView);
/* Coupling on the same element */ /* Coupling on the same element */
for (size_t i = 0; i < insideTrialLocalIndexSet.size(); i++) { distributeLocalPattern(insideTrialLocalView, insideAnsatzLocalView, patternBuilder);
auto row = insideTrialLocalIndexSet.index(i);
for (size_t j = 0; j < insideAnsatzLocalIndexSet.size(); j++) {
auto col = insideAnsatzLocalIndexSet.index(j);
patternBuilder.insertEntry(row,col);
}
}
for (const auto& is : intersections(trialBasis_.gridView(), element)) for (const auto& is : intersections(trialBasis_.gridView(), element))
{ {
/* Elements on the boundary have been treated above, iterate along the inner edges */ /* Elements on the boundary have been treated above, iterate along the inner edges */
if (!is.boundary()) if (is.neighbor())
{ {
// Current outside element // Current outside element
const auto& outsideElement = is.outside(); const auto& outsideElement = is.outside();
// Bind the outer parts to the outer element // Bind the outer parts to the outer element
outsideTrialLocalView.bind(outsideElement);
outsideTrialLocalIndexSet.bind(outsideTrialLocalView);
outsideAnsatzLocalView.bind(outsideElement); outsideAnsatzLocalView.bind(outsideElement);
outsideAnsatzLocalIndexSet.bind(outsideAnsatzLocalView);
// We assume that all basis functions of the inner element couple with all basis functions from the outer one // 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<insideTrialLocalIndexSet.size(); ++i) distributeLocalPattern(insideTrialLocalView, outsideAnsatzLocalView, patternBuilder);
{
auto row = insideTrialLocalIndexSet.index(i);
for (size_t j=0; j<outsideAnsatzLocalIndexSet.size(); ++j)
{
auto col = outsideAnsatzLocalIndexSet.index(j);
patternBuilder.insertEntry(row,col);
}
}
} }
} }
} }
...@@ -135,37 +153,29 @@ public: ...@@ -135,37 +153,29 @@ public:
template <class MatrixBackend, class LocalAssembler> template <class MatrixBackend, class LocalAssembler>
void assembleBulkEntries(MatrixBackend&& matrixBackend, LocalAssembler&& localAssembler) const void assembleBulkEntries(MatrixBackend&& matrixBackend, LocalAssembler&& localAssembler) const
{ {
auto trialLocalView = trialBasis_.localView(); // create two localViews but use only one if bases are the same
auto trialLocalIndexSet = trialBasis_.localIndexSet(); auto ansatzLocalView = ansatzBasis_.localView();
auto separateTrialLocalView = trialBasis_.localView();
auto ansatzLocalView = ansatzBasis_.localView(); auto& trialLocalView = selectTrialLocalView(ansatzLocalView, separateTrialLocalView);
auto ansatzLocalIndexSet = ansatzBasis_.localIndexSet();
using Field = std::decay_t<decltype(matrixBackend(trialLocalIndexSet.index(0), ansatzLocalIndexSet.index(0)))>; using Field = std::decay_t<decltype(matrixBackend(trialLocalView.index(0), ansatzLocalView.index(0)))>;
using LocalMatrix = Dune::Matrix<Dune::FieldMatrix<Field,1,1>>; using LocalMatrix = Dune::Matrix<Dune::FieldMatrix<Field,1,1>>;
auto localMatrix = LocalMatrix(trialLocalView.maxSize(), ansatzLocalView.maxSize()); auto localMatrix = LocalMatrix();
localMatrix.setSize(trialLocalView.maxSize(), ansatzLocalView.maxSize());
for (const auto& element : elements(trialBasis_.gridView())) for (const auto& element : elements(trialBasis_.gridView()))
{ {
trialLocalView.bind(element); // bind the localViews to the element
trialLocalIndexSet.bind(trialLocalView); bind(ansatzLocalView, trialLocalView, element);
ansatzLocalView.bind(element); zeroInitializeLocal(trialLocalView, ansatzLocalView, localMatrix);
ansatzLocalIndexSet.bind(ansatzLocalView);
localAssembler(element, localMatrix, trialLocalView, ansatzLocalView); localAssembler(element, localMatrix, trialLocalView, ansatzLocalView);
// Add element stiffness matrix onto the global stiffness matrix // Add element stiffness matrix onto the global stiffness matrix
for (size_t localRow=0; localRow<trialLocalIndexSet.size(); ++localRow) distributeLocalEntries(trialLocalView, ansatzLocalView, localMatrix, matrixBackend);
{
auto row = trialLocalIndexSet.index(localRow);
for (size_t localCol=0; localCol<ansatzLocalIndexSet.size(); ++localCol)
{
auto col = ansatzLocalIndexSet.index(localCol);
matrixBackend(row,col) += localMatrix[localRow][localCol];
}
}
} }
} }
...@@ -184,35 +194,28 @@ public: ...@@ -184,35 +194,28 @@ public:
Dune::MultipleCodimMultipleGeomTypeMapper<GridView> faceMapper(trialBasis_.gridView(), mcmgElementLayout()); Dune::MultipleCodimMultipleGeomTypeMapper<GridView> faceMapper(trialBasis_.gridView(), mcmgElementLayout());
auto insideTrialLocalView = trialBasis_.localView(); auto insideTrialLocalView = trialBasis_.localView();
auto insideTrialLocalIndexSet = trialBasis_.localIndexSet();
auto insideAnsatzLocalView = ansatzBasis_.localView(); auto insideAnsatzLocalView = ansatzBasis_.localView();
auto insideAnsatzLocalIndexSet = ansatzBasis_.localIndexSet();
auto outsideTrialLocalView = trialBasis_.localView(); auto outsideTrialLocalView = trialBasis_.localView();
auto outsideTrialLocalIndexSet = trialBasis_.localIndexSet();
auto outsideAnsatzLocalView = ansatzBasis_.localView(); auto outsideAnsatzLocalView = ansatzBasis_.localView();
auto outsideAnsatzLocalIndexSet = ansatzBasis_.localIndexSet();
using Field = std::decay_t<decltype(matrixBackend(insideTrialLocalIndexSet.index(0), insideAnsatzLocalIndexSet.index(0)))>; using Field = std::decay_t<decltype(matrixBackend(insideTrialLocalView.index(0), insideAnsatzLocalView.index(0)))>;
using LocalMatrix = Dune::Matrix<Dune::FieldMatrix<Field,1,1>>; using LocalMatrix = Dune::Matrix<Dune::FieldMatrix<Field,1,1>>;
using MatrixContainer = Dune::Matrix<LocalMatrix>; using MatrixContainer = Dune::Matrix<LocalMatrix>;
auto matrixContrainer = MatrixContainer(2,2); auto matrixContrainer = MatrixContainer(2,2);
// Resize matrixContrainer[0][0].setSize(insideTrialLocalView.maxSize(), insideAnsatzLocalView.maxSize());
matrixContrainer[0][0].setSize(insideTrialLocalView.maxSize(), insideAnsatzLocalView.maxSize()); matrixContrainer[0][1].setSize(insideTrialLocalView.maxSize(), outsideAnsatzLocalView.maxSize());
matrixContrainer[0][1].setSize(insideTrialLocalView.maxSize(), outsideAnsatzLocalView.maxSize());
matrixContrainer[1][0].setSize(outsideTrialLocalView.maxSize(), insideAnsatzLocalView.maxSize()); matrixContrainer[1][0].setSize(outsideTrialLocalView.maxSize(), insideAnsatzLocalView.maxSize());
matrixContrainer[1][1].setSize(outsideTrialLocalView.maxSize(), outsideAnsatzLocalView.maxSize()); matrixContrainer[1][1].setSize(outsideTrialLocalView.maxSize(), outsideAnsatzLocalView.maxSize());
for (const auto& element : elements(trialBasis_.gridView())) for (const auto& element : elements(trialBasis_.gridView()))
{ {
insideTrialLocalView.bind(element); insideTrialLocalView.bind(element);
insideTrialLocalIndexSet.bind(insideTrialLocalView);
insideAnsatzLocalView.bind(element); insideAnsatzLocalView.bind(element);
insideAnsatzLocalIndexSet.bind(insideAnsatzLocalView);
for (const auto& is : intersections(trialBasis_.gridView(), element)) for (const auto& is : intersections(trialBasis_.gridView(), element))
{ {
...@@ -221,14 +224,15 @@ public: ...@@ -221,14 +224,15 @@ public:
if (is.boundary()) if (is.boundary())
{ {
auto& localMatrix = matrixContrainer[0][0]; auto& localMatrix = matrixContrainer[0][0];
zeroInitializeLocal(insideTrialLocalView, insideAnsatzLocalView, localMatrix);
localBoundaryAssembler(is, localMatrix, insideTrialLocalView, insideAnsatzLocalView); localBoundaryAssembler(is, localMatrix, insideTrialLocalView, insideAnsatzLocalView);
for (size_t i=0; i< insideTrialLocalIndexSet.size(); i++) { distributeLocalEntries(insideTrialLocalView, insideAnsatzLocalView, localMatrix, matrixBackend);
for (size_t j=0; j< insideAnsatzLocalIndexSet.size(); j++)
matrixBackend(insideTrialLocalIndexSet.index(i), insideAnsatzLocalIndexSet.index(j))+=localMatrix[i][j];
}
} }
else else if (is.neighbor())
{ {
/* Only handle every intersection once; we choose the case where the inner element has the lower index /* Only handle every intersection once; we choose the case where the inner element has the lower index
* *
...@@ -240,46 +244,20 @@ public: ...@@ -240,46 +244,20 @@ public:
// Bind the outer parts to the outer element // Bind the outer parts to the outer element
outsideTrialLocalView.bind(outsideElement); outsideTrialLocalView.bind(outsideElement);
outsideTrialLocalIndexSet.bind(outsideTrialLocalView);
outsideAnsatzLocalView.bind(outsideElement); outsideAnsatzLocalView.bind(outsideElement);
outsideAnsatzLocalIndexSet.bind(outsideAnsatzLocalView);
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); localAssembler(is, matrixContrainer, insideTrialLocalView, insideAnsatzLocalView, outsideTrialLocalView, outsideAnsatzLocalView);
for (size_t i=0; i < insideTrialLocalIndexSet.size(); i++) distributeLocalEntries(insideTrialLocalView, insideAnsatzLocalView, matrixContrainer[0][0], matrixBackend);
{ distributeLocalEntries(insideTrialLocalView, outsideAnsatzLocalView, matrixContrainer[0][1], matrixBackend);
auto row = insideTrialLocalIndexSet.index(i); distributeLocalEntries(outsideTrialLocalView, insideAnsatzLocalView, matrixContrainer[1][0], matrixBackend);
// inside x inside distributeLocalEntries(outsideTrialLocalView, outsideAnsatzLocalView, matrixContrainer[1][1], matrixBackend);
for (size_t j=0; j < insideTrialLocalIndexSet.size(); j++)
{
auto col = insideTrialLocalIndexSet.index(j);
matrixBackend(row, col) += matrixContrainer[0][0][i][j];
}
// inside x outside
for (size_t j=0; j < outsideAnsatzLocalIndexSet.size(); j++)
{
auto col = outsideAnsatzLocalIndexSet.index(j);
matrixBackend(row, col) += matrixContrainer[0][1][i][j];
}
}
for (size_t i=0; i < outsideTrialLocalIndexSet.size(); i++)
{
auto row = outsideTrialLocalIndexSet.index(i);
// outside x inside
for (size_t j=0; j < insideAnsatzLocalIndexSet.size(); j++)
{
auto col = insideAnsatzLocalIndexSet.index(j);
matrixBackend(row, col) += matrixContrainer[1][0][i][j];
}
// outside x outside
for (size_t j=0; j < outsideAnsatzLocalIndexSet.size(); j++)
{
auto col = outsideAnsatzLocalIndexSet.index(j);
matrixBackend(row, col) += matrixContrainer[1][1][i][j];
}
}
} }
} }
} }
...@@ -292,9 +270,34 @@ public: ...@@ -292,9 +270,34 @@ public:
auto patternBuilder = matrixBackend.patternBuilder(); auto patternBuilder = matrixBackend.patternBuilder();
assembleBulkPattern(patternBuilder); assembleBulkPattern(patternBuilder);
patternBuilder.setupMatrix(); patternBuilder.setupMatrix();
matrixBackend.assign(0);
assembleBulkEntries(matrixBackend, std::forward<LocalAssembler>(localAssembler)); 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: protected:
const TrialBasis& trialBasis_; const TrialBasis& trialBasis_;
const AnsatzBasis& ansatzBasis_; const AnsatzBasis& ansatzBasis_;
......