diff --git a/dune/faultnetworks/assemblers/CMakeLists.txt b/dune/faultnetworks/assemblers/CMakeLists.txt
index e3766c4a27eb5692101cee6972d81960f80155fc..84472309798e67593d4c6e6aec083739ad02063e 100644
--- a/dune/faultnetworks/assemblers/CMakeLists.txt
+++ b/dune/faultnetworks/assemblers/CMakeLists.txt
@@ -1,23 +1,14 @@
-add_subdirectory("test" EXCLUDE_FROM_ALL)
add_custom_target(faultnetworks_assemblers_src SOURCES
globalfaultassembler.hh
- oscassembler.hh
interfaceoperatorassembler.hh
intersectionjumpmassassembler.hh
- oscfunctionalassemblereff.hh
- oscl2functionalassembler.hh
osclocalassembler.hh
- oscoperatorassemblereff.hh
)
install(FILES
globalfaultassembler.hh
- oscassembler.hh
interfaceoperatorassembler.hh
intersectionjumpmassassembler.hh
- oscfunctionalassemblereff.hh
- oscl2functionalassembler.hh
osclocalassembler.hh
- oscoperatorassemblereff.hh
-DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/faultnetworks/assemblers)
\ No newline at end of file
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/faultnetworks/assemblers)
diff --git a/dune/faultnetworks/assemblers/oscassembler.hh b/dune/faultnetworks/assemblers/oscassembler.hh
deleted file mode 100644
index 599f24a1ceefab63bac6bd4fb27b52f833f77842..0000000000000000000000000000000000000000
--- a/dune/faultnetworks/assemblers/oscassembler.hh
+++ /dev/null
@@ -1,58 +0,0 @@
-#ifndef OSC_ASSEMBLER_HH
-#define OSC_ASSEMBLER_HH
-
-#include "dune/osccoeff/oscoperatorassemblereff.hh"
-#include "dune/osccoeff/oscfunctionalassemblereff.hh"
-
-//#include "dune/fufem/assemblers/boundaryfunctionalassembler.hh"
-//#include "dune/fufem/boundarypatch.hh"
-
-//! Generic global assembler for operators and functionals
-template <class TrialBasis, class AnsatzBasis>
-class OscAssembler
-{
-
- private:
- typedef typename TrialBasis::GridView::Grid GridType;
- typedef typename GridType::template Codim<0>::Entity EntityType;
-
- const TrialBasis& tBasis_;
- const AnsatzBasis& aBasis_;
-
- std::vector<EntityType>& regionElements_;
-
- public:
- //! create assembler
- OscAssembler(const TrialBasis& tBasis, const AnsatzBasis& aBasis, std::vector<EntityType>& regionElements) :
- tBasis_(tBasis),
- aBasis_(aBasis),
- regionElements_(regionElements) {}
-
- template <class LocalAssemblerType, class GlobalMatrixType>
- void assembleOperator(LocalAssemblerType& localAssembler, GlobalMatrixType& A, const bool lumping=false) const
- {
- const OscOperatorAssembler<TrialBasis, AnsatzBasis> operatorAssembler(tBasis_, aBasis_, regionElements_);
- operatorAssembler.assemble(localAssembler, A, lumping);
- }
-
- template <class LocalFunctionalAssemblerType, class GlobalVectorType>
- void assembleFunctional(LocalFunctionalAssemblerType& localAssembler, GlobalVectorType& b, bool initializeVector=true) const
- {
- const OscFunctionalAssembler<TrialBasis> functionalAssembler(tBasis_, regionElements_);
- functionalAssembler.assemble(localAssembler, b, initializeVector);
- }
-
- /*
- template <class LocalBoundaryFunctionalAssemblerType, class GlobalVectorType>
- void assembleBoundaryFunctional(LocalBoundaryFunctionalAssemblerType& localAssembler,
- GlobalVectorType& b,
- const BoundaryPatch<typename TrialBasis::GridView>& boundaryPatch,
- bool initializeVector=true) const
- {
- const BoundaryFunctionalAssembler<TrialBasis> boundaryFunctionalAssembler(tBasis_, boundaryPatch);
- boundaryFunctionalAssembler.assemble(localAssembler, b, initializeVector);
- }*/
-};
-
-#endif
-
diff --git a/dune/faultnetworks/assemblers/oscfunctionalassemblereff.hh b/dune/faultnetworks/assemblers/oscfunctionalassemblereff.hh
deleted file mode 100644
index 48b102dea19dd5b721f06447ea957741b0535c4f..0000000000000000000000000000000000000000
--- a/dune/faultnetworks/assemblers/oscfunctionalassemblereff.hh
+++ /dev/null
@@ -1,81 +0,0 @@
-#ifndef OSC_FUNCTIONAL_ASSEMBLER_HH
-#define OSC_FUNCTIONAL_ASSEMBLER_HH
-
-#include <dune/istl/bvector.hh>
-
-#include "dune/fufem/functionspacebases/functionspacebasis.hh"
-
-//! Generic global assembler for functionals on a gridview
-template <class TrialBasis>
-class OscFunctionalAssembler
-{
- private:
- typedef typename TrialBasis::GridView GridView;
- typedef typename GridView::Grid GridType;
- typedef typename GridType::template Codim<0>::Entity EntityType;
-
- public:
- //! create assembler for gridview
- OscFunctionalAssembler(const TrialBasis& tBasis, std::vector<EntityType>& regionElements) :
- tBasis_(tBasis),
- regionElements_(regionElements)
- {}
-
- /** \brief Assemble
- *
- * \param localAssembler local assembler
- * \param[out] b target vector
- * \param initializeVector If this is set the output vector is
- * set to the correct length and initialized to zero before
- * assembly. Otherwise the assembled values are just added to
- * the vector.
- */
- template <class LocalFunctionalAssemblerType, class GlobalVectorType>
- void assemble(LocalFunctionalAssemblerType& localAssembler, GlobalVectorType& b, bool initializeVector=true) const
- {
- typedef typename LocalFunctionalAssemblerType::LocalVector LocalVector;
- typedef typename TrialBasis::LinearCombination LinearCombination;
-
- int rows = tBasis_.size();
-
- if (initializeVector)
- {
- b.resize(rows);
- b=0.0;
- }
-
-
- for (size_t i=0; i<regionElements_.size(); ++i)
- {
- EntityType it = regionElements_[i];
-
- // get shape functions
- const typename TrialBasis::LocalFiniteElement& tFE = tBasis_.getLocalFiniteElement(it);
-
- LocalVector localB(tFE.localBasis().size());
- localAssembler.assemble(it, localB, tFE);
-
- for (size_t i=0; i<tFE.localBasis().size(); ++i)
- {
- int idx = tBasis_.index(it, i);
- const LinearCombination& constraints = tBasis_.constraints(idx);
- bool isConstrained = tBasis_.isConstrained(idx);
- if (isConstrained)
- {
- for (size_t w=0; w<constraints.size(); ++w)
- b[constraints[w].index].axpy(constraints[w].factor, localB[i]);
- }
- else
- b[idx] += localB[i];
- }
- }
- return;
- }
-
- private:
- const TrialBasis& tBasis_;
- std::vector<EntityType>& regionElements_;
-};
-
-#endif
-
diff --git a/dune/faultnetworks/assemblers/oscl2functionalassembler.hh b/dune/faultnetworks/assemblers/oscl2functionalassembler.hh
deleted file mode 100644
index 626e58068432ad3fcb66816c4c2a74a6aac3bb57..0000000000000000000000000000000000000000
--- a/dune/faultnetworks/assemblers/oscl2functionalassembler.hh
+++ /dev/null
@@ -1,147 +0,0 @@
-#ifndef OSC_L2_FUNCTIONAL_ASSEMBLER_HH
-#define OSC_L2_FUNCTIONAL_ASSEMBLER_HH
-
-
-#include <dune/common/fvector.hh>
-#include <dune/common/function.hh>
-#include <dune/istl/bvector.hh>
-
-#include <dune/fufem/quadraturerules/quadraturerulecache.hh>
-#include <dune/fufem/functions/virtualgridfunction.hh>
-
-#include "dune/fufem/assemblers/localfunctionalassembler.hh"
-
-#include <dune/faultnetworks/oscdata.hh>
-
-
-// TODO: oscData_ not implemented
-
-/** \brief Local finite element assembler for L^2 functionals. */
-template <class GridType, class TrialLocalFE, class T=Dune::FieldVector<double,1> >
-class OscL2FunctionalAssembler :
- public LocalFunctionalAssembler<GridType, TrialLocalFE, T>
-
-{
- private:
- typedef LocalFunctionalAssembler<GridType, TrialLocalFE, T> Base;
- static const int dim = GridType::dimension;
- static const int dimworld = GridType::dimensionworld;
-
- typedef typename GridType::ctype ctype;
- typedef typename GridType::template Codim<0>::Geometry::GlobalCoordinate GlobalCoordinate;
- typedef typename GridType::template Codim<0>::Geometry::LocalCoordinate LocalCoordinate;
- typedef VirtualGridFunction<GridType, T> GridFunction;
-
- public:
- typedef typename Base::Element Element;
- typedef typename Base::Element::Geometry Geometry;
- typedef typename Base::LocalVector LocalVector;
-
- typedef typename Dune::VirtualFunction<GlobalCoordinate, T> Function;
-
- /**
- * \brief Create L2FunctionalAssembler
- *
- * Creates a local functional assembler for an L2-functional.
- * The QuadratureRuleKey given here does only specify what is
- * needed to integrate f. Depending on this the quadrature rule
- * is selected automatically such that fv for test functions v
- * can be integrated.
- * If you, e.g., specify quadrature order 1 for f and the test
- * functions need order 2 the selected quadrature rule will
- * have order 3.
- *
- * \param f The L2 function
- * \param fQuadKey A QuadratureRuleKey that specifies how to integrate f
- */
- OscL2FunctionalAssembler(const Function& f, const QuadratureRuleKey& fQuadKey) :
- f_(f),
- quadOrder_(-1),
- functionQuadKey_(fQuadKey)
- {}
-
- /**
- * \brief Create L2FunctionalAssembler
- *
- * Using this constructor the quadrature order will be selected
- * such that the test functions can be integrated exactly.
- *
- * \param f The L2 function
- */
- OscL2FunctionalAssembler(const Function& f, OscData<> oscData) :
- f_(f),
- quadOrder_(-1),
- functionQuadKey_(dim, 0),
- region_(region)
- {}
-
- //! create assembler for grid
- DUNE_DEPRECATED_MSG("Please select quadrature for L2FunctionalAssembler using a QuadratureRuleKey.")
- OscL2FunctionalAssembler(const Function& f, int order) :
- f_(f),
- quadOrder_(order),
- functionQuadKey_(dim, order)
- {}
-
- virtual void assemble(const Element& element, LocalVector& localVector, const TrialLocalFE& tFE) const
- {
- typedef typename TrialLocalFE::Traits::LocalBasisType::Traits::RangeType RangeType;
-
- // get geometry and store it
- const Geometry geometry = element.geometry();
-
- localVector = 0.0;
-
- bool inBox;
- region_.evaluate(element, inBox);
- if (!inBox)
- return;
-
- // get quadrature rule
- QuadratureRuleKey tFEquad(tFE);
- QuadratureRuleKey quadKey = tFEquad.product(functionQuadKey_);
- if (quadOrder_>=0)
- quadKey.setOrder(quadOrder_);
- const Dune::template QuadratureRule<double, dim>& quad = QuadratureRuleCache<double, dim>::rule(quadKey);
-
- // store values of shape functions
- std::vector<RangeType> values(tFE.localBasis().size());
-
- // loop over quadrature points
- for (size_t pt=0; pt < quad.size(); ++pt)
- {
- // get quadrature point
- const LocalCoordinate& quadPos = quad[pt].position();
-
- // get integration factor
- const double integrationElement = geometry.integrationElement(quadPos);
-
- // evaluate basis functions
- tFE.localBasis().evaluateFunction(quadPos, values);
-
- // compute values of function
- T f_pos;
- const GridFunction* gf = dynamic_cast<const GridFunction*>(&f_);
- if (gf and gf->isDefinedOn(element))
- gf->evaluateLocal(element, quadPos, f_pos);
- else
- f_.evaluate(geometry.global(quadPos), f_pos);
-
- // and vector entries
- for(size_t i=0; i<values.size(); ++i)
- {
- localVector[i].axpy(values[i]*quad[pt].weight()*integrationElement, f_pos);
- }
- }
- return;
- }
-
- private:
- const Function& f_;
- const int quadOrder_;
- const QuadratureRuleKey functionQuadKey_;
-
- OscData<typename GridType::template Codim<0>::Entity, Dune::FieldVector<ctype, dimworld>>& region_;
-};
-#endif
-
diff --git a/dune/faultnetworks/assemblers/oscoperatorassemblereff.hh b/dune/faultnetworks/assemblers/oscoperatorassemblereff.hh
deleted file mode 100644
index 6fdc30b34b7ac1c451fdfd5ac2a076063ea42df1..0000000000000000000000000000000000000000
--- a/dune/faultnetworks/assemblers/oscoperatorassemblereff.hh
+++ /dev/null
@@ -1,214 +0,0 @@
-#ifndef OSC_OPERATOR_ASSEMBLER_HH
-#define OSC_OPERATOR_ASSEMBLER_HH
-
-#include <dune/istl/matrix.hh>
-#include <dune/istl/matrixindexset.hh>
-
-#include "dune/fufem/functionspacebases/functionspacebasis.hh"
-#include "dune/fufem/staticmatrixtools.hh"
-
-//! Generic global assembler for operators on a gridview
-template <class TrialBasis, class AnsatzBasis>
-class OscOperatorAssembler
-{
- private:
- typedef typename TrialBasis::GridView GridView;
- typedef typename GridView::Grid GridType;
- typedef typename GridType::template Codim<0>::Entity EntityType;
-
- public:
- //! create assembler for grid
- OscOperatorAssembler(const TrialBasis& tBasis, const AnsatzBasis& aBasis, std::vector<EntityType>& regionElements) :
- tBasis_(tBasis),
- aBasis_(aBasis),
- regionElements_(regionElements)
- {}
-
-
- template <class LocalAssemblerType>
- void addIndices(LocalAssemblerType& localAssembler, Dune::MatrixIndexSet& indices, const bool lumping=false) const
- {
- if (lumping)
- addIndicesStaticLumping<LocalAssemblerType,true>(localAssembler, indices);
- else
- addIndicesStaticLumping<LocalAssemblerType,false>(localAssembler, indices);
- }
-
-
- template <class LocalAssemblerType, class GlobalMatrixType>
- void addEntries(LocalAssemblerType& localAssembler, GlobalMatrixType& A, const bool lumping=false) const
- {
- if (lumping)
- addEntriesStaticLumping<LocalAssemblerType,GlobalMatrixType,true>(localAssembler, A);
- else
- addEntriesStaticLumping<LocalAssemblerType,GlobalMatrixType,false>(localAssembler, A);
- }
-
-
- template <class LocalAssemblerType, class GlobalMatrixType>
- void assemble(LocalAssemblerType& localAssembler, GlobalMatrixType& A, const bool lumping=false) const
- {
- int rows = tBasis_.size();
- int cols = aBasis_.size();
-
- Dune::MatrixIndexSet indices(rows, cols);
-
- addIndices(localAssembler, indices, lumping);
-
- indices.exportIdx(A);
- A=0.0;
-
- addEntries(localAssembler, A, lumping);
-
- return;
- }
-
-
- protected:
-
- template <class LocalAssemblerType, bool lumping>
- void addIndicesStaticLumping(LocalAssemblerType& localAssembler, Dune::MatrixIndexSet& indices) const
- {
- typedef typename LocalAssemblerType::BoolMatrix BoolMatrix;
- typedef typename TrialBasis::LinearCombination LinearCombination;
-
- for (size_t i=0; i<regionElements_.size(); ++i)
- {
- EntityType it = regionElements_[i];
-
- // get shape functions
- const typename TrialBasis::LocalFiniteElement& tFE = tBasis_.getLocalFiniteElement(it);
- const typename AnsatzBasis::LocalFiniteElement& aFE = aBasis_.getLocalFiniteElement(it);
-
- BoolMatrix localIndices(tFE.localBasis().size(), aFE.localBasis().size());
- localAssembler.indices(it, localIndices, tFE, aFE);
-
- for (size_t i=0; i<tFE.localBasis().size(); ++i)
- {
- int rowIndex = tBasis_.index(it, i);
- const LinearCombination& rowConstraints = tBasis_.constraints(rowIndex);
- bool rowIsConstrained = tBasis_.isConstrained(rowIndex);
- for (size_t j=0; j<aFE.localBasis().size(); ++j)
- {
- if (localIndices[i][j])
- {
- if (lumping)
- {
- if (rowIsConstrained)
- {
- for (size_t rw=0; rw<rowConstraints.size(); ++rw)
- indices.add(rowConstraints[rw].index, rowConstraints[rw].index);
- }
- else
- indices.add(rowIndex, rowIndex);
- }
- else
- {
- int colIndex = aBasis_.index(it, j);
- const LinearCombination& colConstraints = aBasis_.constraints(colIndex);
- bool colIsConstrained = aBasis_.isConstrained(colIndex);
- if (rowIsConstrained and colIsConstrained)
- {
- for (size_t rw=0; rw<rowConstraints.size(); ++rw)
- {
- for (size_t cw=0; cw<colConstraints.size(); ++cw)
- indices.add(rowConstraints[rw].index, colConstraints[cw].index);
- }
- }
- else if (rowIsConstrained)
- {
- for (size_t rw=0; rw<rowConstraints.size(); ++rw)
- indices.add(rowConstraints[rw].index, colIndex);
- }
- else if (colIsConstrained)
- {
- for (size_t cw=0; cw<colConstraints.size(); ++cw)
- indices.add(rowIndex, colConstraints[cw].index);
- }
- else
- indices.add(rowIndex, colIndex);
- }
- }
- }
- }
- }
- }
-
-
- template <class LocalAssemblerType, class GlobalMatrixType, bool lumping>
- void addEntriesStaticLumping(LocalAssemblerType& localAssembler, GlobalMatrixType& A) const
- {
- typedef typename LocalAssemblerType::LocalMatrix LocalMatrix;
- typedef typename TrialBasis::LinearCombination LinearCombination;
-
- for (size_t i=0; i<regionElements_.size(); ++i)
- {
- EntityType it = regionElements_[i];
-
- // get shape functions
- const typename TrialBasis::LocalFiniteElement& tFE = tBasis_.getLocalFiniteElement(it);
- const typename AnsatzBasis::LocalFiniteElement& aFE = aBasis_.getLocalFiniteElement(it);
-
- LocalMatrix localA(tFE.localBasis().size(), aFE.localBasis().size());
- localAssembler.assemble(it, localA, tFE, aFE);
-
- for (size_t i=0; i<tFE.localBasis().size(); ++i)
- {
- int rowIndex = tBasis_.index(it, i);
- const LinearCombination& rowConstraints = tBasis_.constraints(rowIndex);
- bool rowIsConstrained = tBasis_.isConstrained(rowIndex);
- for (size_t j=0; j<aFE.localBasis().size(); ++j)
- {
-
- if (localA[i][j].infinity_norm()!=0.0)
- {
- if (lumping)
- {
- if (rowIsConstrained)
- {
- for (size_t rw=0; rw<rowConstraints.size(); ++rw)
- Arithmetic::addProduct(A[rowConstraints[rw].index][rowConstraints[rw].index], rowConstraints[rw].factor, localA[i][j]);
- }
- else
- Arithmetic::addProduct(A[rowIndex][rowIndex], 1.0, localA[i][j]);
- }
- else
- {
- int colIndex = aBasis_.index(it, j);
- const LinearCombination& colConstraints = aBasis_.constraints(colIndex);
- bool colIsConstrained = aBasis_.isConstrained(colIndex);
- if (rowIsConstrained and colIsConstrained)
- {
- for (size_t rw=0; rw<rowConstraints.size(); ++rw)
- {
- for (size_t cw=0; cw<colConstraints.size(); ++cw)
- Arithmetic::addProduct(A[rowConstraints[rw].index][colConstraints[cw].index], rowConstraints[rw].factor * colConstraints[cw].factor, localA[i][j]);
- }
- }
- else if (rowIsConstrained)
- {
- for (size_t rw=0; rw<rowConstraints.size(); ++rw)
- Arithmetic::addProduct(A[rowConstraints[rw].index][colIndex], rowConstraints[rw].factor, localA[i][j]);
- }
- else if (colIsConstrained)
- {
- for (size_t cw=0; cw<colConstraints.size(); ++cw)
- Arithmetic::addProduct(A[rowIndex][colConstraints[cw].index], colConstraints[cw].factor, localA[i][j]);
- }
- else
- Arithmetic::addProduct(A[rowIndex][colIndex], 1.0, localA[i][j]);
- }
- }
- }
- }
- }
- }
-
-
- const TrialBasis& tBasis_;
- const AnsatzBasis& aBasis_;
- std::vector<EntityType>& regionElements_;
-};
-
-#endif
-
diff --git a/dune/faultnetworks/assemblers/test/CMakeLists.txt b/dune/faultnetworks/assemblers/test/CMakeLists.txt
deleted file mode 100644
index c2a45e21ef5681312921e943d7f6f51ff326cc99..0000000000000000000000000000000000000000
--- a/dune/faultnetworks/assemblers/test/CMakeLists.txt
+++ /dev/null
@@ -1,12 +0,0 @@
-# Put your test in here if it needs access to external grids
-set(FAULTNETWORKS_ASSEMBLERS_TESTS
- globalproblemassemblertest
- )
-
-set(TESTS ${FAULTNETWORKS_ASSEMBLERS_TESTS})
-
-# Setup targets, register tests, and add dune flags
-foreach(_test ${TESTS})
- #dune_add_test(SOURCES ${_test}.cc)
-endforeach()
-
diff --git a/dune/faultnetworks/assemblers/test/Makefile.am b/dune/faultnetworks/assemblers/test/Makefile.am
deleted file mode 100644
index a233502edfe296766368b94f0ae3913464aafffa..0000000000000000000000000000000000000000
--- a/dune/faultnetworks/assemblers/test/Makefile.am
+++ /dev/null
@@ -1,42 +0,0 @@
-# which tests where program to build and run are equal
-NORMALTESTS = globalproblemassemblertest
-
-# list of tests to run
-TESTS = $(NORMALTESTS)
-
-# programs just to build when "make check" is used
-check_PROGRAMS = $(NORMALTESTS)
-
-# define the programs
-globalproblemassemblertest_SOURCES = globalproblemassemblertest.cc
-
-
-globalproblemassemblertest_CPPFLAGS = $(AM_CPPFLAGS) \
- $(DUNEMPICPPFLAGS) \
- $(UG_CPPFLAGS)
-
-# The libraries have to be given in reverse order (most basic libraries
-# last). Also, due to some misunderstanding, a lot of libraries include the
-# -L option in LDFLAGS instead of LIBS -- so we have to include the LDFLAGS
-# here as well.
-globalproblemassemblertest_LDADD = \
- $(DUNE_LDFLAGS) $(DUNE_LIBS) \
- $(UG_LDFLAGS) $(UG_LIBS) \
- $(DUNEMPILIBS) \
- $(LDADD)
-globalproblemassemblertest_LDFLAGS = $(AM_LDFLAGS) \
- $(DUNEMPILDFLAGS) \
- $(UG_LDFLAGS) \
- $(DUNE_LDFLAGS)
-
-# don't follow the full GNU-standard
-# we need automake 1.9
-AUTOMAKE_OPTIONS = foreign 1.9
-
-# pass most important options when "make distcheck" is used
-DISTCHECK_CONFIGURE_FLAGS = --with-dune-common=$(DUNE_COMMON_ROOT) --with-dune-geometry=$(DUNE_GEOMETRY_ROOT) --with-dune-grid=$(DUNE_GRID_ROOT) --with-dune-istl=$(DUNE_ISTL_ROOT) --with-dune-localfunctions=$(DUNE_LOCALFUNCTIONS_ROOT) --with-dune-solvers=$(DUNE_SOLVERS_ROOT) --with-dune-fufem=$(DUNE_FUFEM_ROOT) CXX="$(CXX)" CC="$(CC)"
-
-EXTRA_DIST = CMakeLists.txt
-
-include $(top_srcdir)/am/global-rules
-
diff --git a/dune/faultnetworks/assemblers/test/globalproblemassemblertest.cc b/dune/faultnetworks/assemblers/test/globalproblemassemblertest.cc
deleted file mode 100644
index e9e28fb00b2b37fc4a587d0795eb138cd3257b89..0000000000000000000000000000000000000000
--- a/dune/faultnetworks/assemblers/test/globalproblemassemblertest.cc
+++ /dev/null
@@ -1,304 +0,0 @@
-#include <config.h>
-
-#include <dune/common/parallel/mpihelper.hh>
-
-#include <dune/grid/uggrid.hh>
-#include <dune/grid/common/mcmgmapper.hh>
-
-#include <dune/faultnetworks/assemblers/globalproblemassembler.hh>
-#include <dune/faultnetworks/faultp1nodalbasis.hh>
-
-#include <dune/faultnetworks/assemblers/intersectionjumpmassassembler.hh>
-#include <dune/fufem/assemblers/localassemblers/l2functionalassembler.hh>
-#include <dune/fufem/assemblers/localassemblers/laplaceassembler.hh>
-#include <dune/fufem/boundarypatch.hh>
-#include <dune/fufem/functiontools/boundarydofs.hh>
-
-#include <dune/faultnetworks/interfacenetwork.hh>
-#include <dune/faultnetworks/oscdata.hh>
-#include <dune/faultnetworks/oscdatahierarchy.hh>
-#include <dune/faultnetworks/oscrhs.hh>
-
-#include <dune/faultnetworks/faultfactories/multileveluniformfaultfactory.hh>
-#include <dune/faultnetworks/utils/debugutils.hh>
-
-const int dim = 2;
-
-const int coarseResolution = 2;
-const int fineResolution = 3;
-const int exactResolution = 4;
-
-const int faultCount = 3;
-
-typedef typename Dune::UGGrid<dim> GridType;
-typedef typename GridType::LevelGridView GridView;
-typedef typename GridType::LevelIntersection Intersection;
-typedef typename GridType::template Codim<0>::Entity Element;
-
-typedef FaultP1NodalBasis<GridView, double> DGBasis;
-typedef typename DGBasis::LocalFiniteElement LFE;
-
-typedef Dune::FieldVector<double, dim> WorldVectorType;
-typedef Dune::BlockVector<Dune::FieldVector<double,1>> VectorType;
-typedef Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>> MatrixType;
-typedef Dune::Matrix<Dune::FieldMatrix<double,1,1>, std::allocator<Dune::FieldMatrix<double,1,1>>> OscMatrixType;
-
-typedef IntersectionJumpMassAssembler<GridView, LFE, LFE> LocalInterfaceAssembler;
-
-static const int dimElement = Element::dimension;
-
-void gridLevelOutput(const GridType& grid, const int level) {
- std::cout << "Grid elements on lvl " << level << std::endl << std::endl;
-
- typedef typename GridView::template Codim<0>::Iterator ElementLevelIterator;
- GridView gridView(grid.levelGridView(level));
-
- //Dune::LevelMultipleCodimMultipleGeomTypeMapper <GridType, Dune::MCMGElementLayout > mapper(grid, level);
- Dune::LevelMultipleCodimMultipleGeomTypeMapper <GridType, Dune::MCMGVertexLayout > vMapper(grid, level);
-
- Dune::BitSetVector<1> vertexPrinted(vMapper.size(), false);
-
- ElementLevelIterator elemIt = gridView.template begin<0>();
- ElementLevelIterator elemEndIt = gridView.template end<0>();
- for (; elemIt!=elemEndIt; ++elemIt) {
- const Element& elem = *elemIt;
- const auto& geometry = elem.geometry();
-
- for(int i=0; i<geometry.corners(); ++i) {
- size_t globalIdx = gridView.indexSet().subIndex(elem, i, Element::dimension);
- if (!vertexPrinted[globalIdx][0]) {
- const auto& vertex = geometry.corner(i);
- print(vertex, "elem corner:" + std::to_string(globalIdx));
- vertexPrinted[globalIdx][0] = true;
- }
- }
- }
- std::cout << "----------------" << std::endl << std::endl;
-
- /*
- ElementLevelIterator elemIt = gridView.template begin<0>();
- ElementLevelIterator elemEndIt = gridView.template end<0>();
- for (; elemIt!=elemEndIt; ++elemIt) {
- const Element& elem = *elemIt;
- const auto& geometry = elem.geometry();
-
- std::cout << "Element idx: " << mapper.index(elem) << std::endl;
-
- for(int i=0; i<geometry.corners(); ++i) {
- size_t globalIdx = gridView.indexSet().subIndex(elem, i, Element::dimension);
- const auto& vertex = geometry.corner(i);
- print(vertex, "elem corner:" + std::to_string(globalIdx));
- }
- std::cout << "----------------" << std::endl << std::endl;
- }*/
-
-}
-
-
-void faultInterfaceOutput(const FaultInterface<GridView>& faultInterface) {
- std::cout << "Level: " << faultInterface.level() << " global faultID: " << faultInterface.id() << std::endl;
-
- const std::set<size_t>& interfaceDofs = faultInterface.getInterfaceDofs();
- print(interfaceDofs, "");
-}
-
-void interfaceNetworkLevelOutput(const InterfaceNetwork<GridType>& interfaceNetwork, const int level) {
- const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork = interfaceNetwork.levelInterfaceNetwork(level);
-
- std::cout << "LevelInterfaceNetwork level: " << level << std::endl;
- for (size_t i=0; i<levelInterfaceNetwork.size(); i++) {
- std::cout << "local faultID : " << i << " ";
- faultInterfaceOutput(*(levelInterfaceNetwork.getInterface(i)));
- }
-
- std::cout << "------------------------------" << std::endl << std::endl;
-}
-
-void computeIntersectionDofs(const GridView& gridView, const Intersection& intersection, std::set<size_t>& intersectionDofs) {
- intersectionDofs.clear();
-
- // loop over all vertices of the intersection
- const Element& insideElement = intersection.inside();
-
- const auto& refElement = Dune::ReferenceElements<double,dimElement>::general(insideElement.type());
- for (int i=0; i<refElement.size(intersection.indexInInside(), 1, dimElement); i++) {
- size_t idxInElement = refElement.subEntity(intersection.indexInInside(), 1, i, dimElement);
- size_t globalIdx = gridView.indexSet().subIndex(insideElement, idxInElement, dimElement);
- intersectionDofs.insert(globalIdx);
- }
-}
-
-void prolongFaultDofs(const GridType& grid, const FaultInterface<GridView>& faultInterface, const int toLevel, std::set<size_t>& prolongedDofs) {
- prolongedDofs.clear();
-
- const std::vector<Intersection>& faces = faultInterface.faces();
- const GridView gridView(grid.levelGridView(toLevel));
-
- for (size_t i=0; i<faces.size(); i++) {
- const Intersection& intersection = faces[i];
- Face<GridType> face(grid, intersection.inside(), intersection.indexInInside());
-
- typename Face<GridType>::HierarchicIterator exactIt = face.hbegin(toLevel);
- typename Face<GridType>::HierarchicIterator exactEnd = face.hend(toLevel);
- for(; exactIt!=exactEnd; ++exactIt) {
- if (exactIt->level() == toLevel) {
- std::set<size_t> faceDofs;
- computeIntersectionDofs(gridView, *(exactIt->intersection()), faceDofs);
-
- prolongedDofs.insert(faceDofs.begin(), faceDofs.end());
- }
- }
- }
-}
-
-void prolongInterfaceNetworkLevelDofs(const GridType& grid, const InterfaceNetwork<GridType>& interfaceNetwork, const int level, const int toLevel, std::vector<std::set<size_t>>& prolongedFaultDofs) {
- prolongedFaultDofs.resize(0);
-
- const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork = interfaceNetwork.levelInterfaceNetwork(level);
-
- std::cout << "Correct prolonged LevelInterfaceNetwork level " << level << " to level " << toLevel << std::endl;
- for (size_t i=0; i<levelInterfaceNetwork.size(); i++) {
- std::set<size_t> prolongedDofs;
- prolongFaultDofs(grid, *(levelInterfaceNetwork.getInterface(i)), toLevel, prolongedDofs);
-
- print(prolongedDofs, "local faultID: " + std::to_string(i));
- prolongedFaultDofs.push_back(prolongedDofs);
- }
-
- std::cout << "------------------------------" << std::endl << std::endl;
-}
-
-int main(int argc, char** argv) try
-{
- Dune::MPIHelper::instance(argc, argv);
-
- std::ofstream out("/home/mi/podlesny/software/dune/faultnetworks/dune/faultnetworks/assemblers/test/globalproblemassemblertest.log");
- std::streambuf *coutbuf = std::cout.rdbuf(); //save old buf
- std::cout.rdbuf(out.rdbuf()); //redirect std::cout to out.txt!
-
- const int fineLevelIdx = fineResolution - coarseResolution;
- const int exactLevelIdx = exactResolution - coarseResolution;
-
- std::vector<int> faultToLevel(faultCount);
- for (size_t i=0; i<faultToLevel.size(); ) {
- for (int j=fineLevelIdx; j>0 && i<faultToLevel.size(); j--) {
- faultToLevel[i] = j;
- i++;
- }
-
- for (int j=0; j<fineLevelIdx && i<faultToLevel.size(); j++) {
- faultToLevel[i] = j;
- i++;
- }
- }
-
- print(faultToLevel, "faultToLevel: ");
-
- MultilevelUniformFaultFactory<GridType> faultFactory(faultCount, faultToLevel, coarseResolution, exactLevelIdx);
- const GridType& grid = faultFactory.grid();
- const InterfaceNetwork<GridType>& interfaceNetwork = faultFactory.interfaceNetwork();
-
- // interface integral jump penalty, B:
- OscMatrixType B(2,2);
- B[0][0] = 1;
- B[1][1] = 1;
- B[0][1] = -1;
- B[1][0] = -1;
-
- // setting local assemblers for laplace operator
- OscRhs<WorldVectorType, Dune::FieldVector<double,1>> f;
- L2FunctionalAssembler<GridType, LFE> l2FunctionalAssembler(f);
- LocalInterfaceAssembler localInterfaceAssembler(B);
- LaplaceAssembler<GridType, LFE, LFE> localAssembler;
-
- GlobalProblemAssembler<DGBasis, MatrixType, VectorType> globalAssembler(interfaceNetwork, exactLevelIdx);
-
- // constructor test
- // ---------------------------------------------------------------
-
- // levelInterfaceNetwork test
-
- const LevelInterfaceNetwork<GridView>& leafInterfaceNetwork = globalAssembler.levelInterfaceNetwork();
-
- std::cout << "Computed LeafInterfaceNetwork level " << exactLevelIdx << std::endl;
- for (size_t i=0; i<leafInterfaceNetwork.size(); i++) {
- print(leafInterfaceNetwork.getInterface(i)->getInterfaceDofs(), "local faultID: " + std::to_string(i));
- }
-
- std::vector<std::set<size_t>> correctProlongedFaultDofs;
- prolongInterfaceNetworkLevelDofs(grid, interfaceNetwork, 0, exactLevelIdx, correctProlongedFaultDofs);
- prolongInterfaceNetworkLevelDofs(grid, interfaceNetwork, fineLevelIdx, exactLevelIdx, correctProlongedFaultDofs);
-
- size_t allProlongedFaultsCorrect = 0;
- for (size_t i=0; i<correctProlongedFaultDofs.size(); i++) {
- for (size_t j=0; j<leafInterfaceNetwork.size(); j++) {
- const std::set<size_t>& faultDofs = leafInterfaceNetwork.getInterface(j)->getInterfaceDofs();
-
- if (correctProlongedFaultDofs[i]==faultDofs) {
- allProlongedFaultsCorrect++;
- break;
- }
- }
- }
-
- assert(allProlongedFaultsCorrect == correctProlongedFaultDofs.size());
- print(globalAssembler.levelInterfaceNetwork().getInterfaceNetworkDofs(), "GlobalProblemAssembler::levelInterfaceNetwork::getInterfaceNetworkDofs :");
- // ---------------------------------------------------------------
-
-
- // assemble() test
- // ---------------------------------------------------------------
- globalAssembler.assemble(localAssembler, localInterfaceAssembler, l2FunctionalAssembler);
-
- // basis test
- const DGBasis& basis = globalAssembler.basis();
- assert(basis.size() == leafInterfaceNetwork.dofCount() + leafInterfaceNetwork.levelGridView().size(dim));
-
- GridView gridView(grid.levelGridView(exactLevelIdx));
- std::cout << "regular dofs: " << 0 << "..." << gridView.size(dim)-1 << " fault dofs : " << gridView.size(dim) << ".." << basis.size()-1 << std::endl;
-
- // boundary dofs test
- Dune::BitSetVector<1> boundaryDofs;
- BoundaryPatch<GridView> boundaryPatch(basis.getGridView(), true);
- constructBoundaryDofs(boundaryPatch, basis, boundaryDofs);
- print(boundaryDofs, "boundaryDofs: ");
-
- bool hasDiagonalZero = false;
- const MatrixType& matrix = globalAssembler.matrix();
- for (size_t i=0; i<matrix.N(); i++) {
- if (!matrix.exists(i,i)) {
- std::cout << i << " ";
- hasDiagonalZero = true;
- }
- }
- std::cout << std::endl;
-
- assert(!hasDiagonalZero);
-
-
-
- // ---------------------------------------------------------------
-
-
- // solve() test
- // ---------------------------------------------------------------
- // globalAssembler.solve();
-
- // const VectorType& sol = globalAssembler.getSol();
- // std::cout << "exactBasis.size() " << exactBasis.size() << std::endl;
-
- /*
- gridLevelOutput(grid, 0);
- gridLevelOutput(grid, fineLevelIdx);
- gridLevelOutput(grid, exactLevelIdx);
-
- interfaceNetworkLevelOutput(interfaceNetwork, 0);
- interfaceNetworkLevelOutput(interfaceNetwork, fineLevelIdx);
- interfaceNetworkLevelOutput(interfaceNetwork, exactLevelIdx);
-
-*/
-
- std::cout.rdbuf(coutbuf);
-} catch (Dune::Exception e){
- std::cout << e << std::endl;
-}
diff --git a/dune/faultnetworks/assemblers/test/globalproblemassemblertest.log b/dune/faultnetworks/assemblers/test/globalproblemassemblertest.log
deleted file mode 100644
index 4067b35ea0beab6ffc23a2e5001472234b31c550..0000000000000000000000000000000000000000
--- a/dune/faultnetworks/assemblers/test/globalproblemassemblertest.log
+++ /dev/null
@@ -1,45 +0,0 @@
-faultToLevel:
-1 0 1
-
-Coarse and fine grid were generated!
- Step 1!
- Step 3!
-fault 0 faces.size(): 8
-fault 1 faces.size(): 4
-fault 2 faces.size(): 8
-4
-GlobalProblemAssembler: levelInterfaceNetwork.size() 3
-GlobalProblemAssembler: levelInterfaceNetwork.getInterface(0).size() 16
-Computed LeafInterfaceNetwork level 2
-local faultID: 0
-11 15 17 20 23 35 39 41 43 55 59 61 63 75 79 81 83
-
-local faultID: 1
-91 95 97 100 103 111 115 117 119 127 131 133 135 143 147 149 151
-
-local faultID: 2
-159 163 165 168 171 179 183 185 187 195 199 201 203 211 215 217 219
-
-Correct prolonged LevelInterfaceNetwork level 0 to level 2
-local faultID: 0
-91 95 97 100 103 111 115 117 119 127 131 133 135 143 147 149 151
-
-------------------------------
-
-Correct prolonged LevelInterfaceNetwork level 1 to level 2
-local faultID: 0
-11 15 17 20 23 35 39 41 43 55 59 61 63 75 79 81 83
-
-local faultID: 1
-159 163 165 168 171 179 183 185 187 195 199 201 203 211 215 217 219
-
-------------------------------
-
-GlobalProblemAssembler::levelInterfaceNetwork::getInterfaceNetworkDofs :
-11 15 17 20 23 35 39 41 43 55 59 61 63 75 79 81 83 91 95 97 100 103 111 115 117 119 127 131 133 135 143 147 149 151 159 163 165 168 171 179 183 185 187 195 199 201 203 211 215 217 219
-
-regular dofs: 0...288 fault dofs : 289..339
-boundaryDofs:
-1 1 0 1 0 0 1 0 1 0 0 0 0 0 0 1 1 0 1 0 0 1 0 0 0 1 0 1 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 1 1 1 1 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 1 1 1 0 1 1 0 1 0 0 0 0 0 0 0 1 0 0 0 1 0 1 0 1 0 0 0 0 0 0 0 1 0 0 0 1 0 1 0 1 0 0 0 0 1 1 0 1 1 0 0 1 0 1 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0
-
-
diff --git a/dune/faultnetworks/faultfactories/CMakeLists.txt b/dune/faultnetworks/faultfactories/CMakeLists.txt
index 2b49f6f00465244fb90842c3c2ce2d313c6f136c..2b3253225373ea179fd3463b3603c2fbe84a1bed 100644
--- a/dune/faultnetworks/faultfactories/CMakeLists.txt
+++ b/dune/faultnetworks/faultfactories/CMakeLists.txt
@@ -1,13 +1,9 @@
-add_subdirectory("test" EXCLUDE_FROM_ALL)
add_custom_target(faultnetworks_faultfactories_src SOURCES
oscunitcube.hh
geofaultfactory.hh
cantorconvergencefaultfactory.hh
cantorfaultfactory.hh
- multileveluniformfaultfactory.hh
- twofaultfactory.hh
- uniformfaultfactory.hh
)
install(FILES
@@ -15,7 +11,4 @@ install(FILES
geofaultfactory.hh
cantorconvergencefaultfactory.hh
cantorfaultfactory.hh
- multileveluniformfaultfactory.hh
- twofaultfactory.hh
- uniformfaultfactory.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/faultnetworks/faultfactories)
diff --git a/dune/faultnetworks/faultfactories/multileveluniformfaultfactory.hh b/dune/faultnetworks/faultfactories/multileveluniformfaultfactory.hh
deleted file mode 100644
index a196429741aa29223b83c1a46bd496fef4bc8bbc..0000000000000000000000000000000000000000
--- a/dune/faultnetworks/faultfactories/multileveluniformfaultfactory.hh
+++ /dev/null
@@ -1,242 +0,0 @@
-#ifndef UNIFORM_FAULT_FACTORY_HH
-#define UNIFORM_FAULT_FACTORY_HH
-
-#include <math.h>
-
-#include <dune/faultnetworks/facehierarchy.hh>
-#include <dune/faultnetworks/utils/debugutils.hh>
-#include <dune/faultnetworks/faultfactories/oscunitcube.hh>
-#include <dune/faultnetworks/levelinterfacenetwork.hh>
-#include <dune/faultnetworks/interfacenetwork.hh>
-#include <dune/faultnetworks/faultinterface.hh>
-#include <dune/faultnetworks/hierarchicleveliterator.hh>
-
-#include <dune/grid/common/mcmgmapper.hh>
-#include <dune/grid/io/file/vtk/vtkwriter.hh>
-
-
-template <class GridType>
-class MultilevelUniformFaultFactory
-{
- //! Parameter for mapper class
- template<int dim>
- struct FaceMapperLayout
- {
- bool contains (Dune::GeometryType gt)
- {
- return gt.dim() == dim-1;
- }
- };
-
- protected:
- typedef OscUnitCube<GridType, 2> GridOb;
- static const int dimworld = GridType::dimensionworld ;
- static const int dim = GridType::dimension;
- typedef typename GridType::ctype ctype;
-
- typedef typename GridType::LevelGridView GV;
- typedef typename GridType::LevelIntersection Intersection;
- typedef typename Intersection::Geometry IntersectionGeometry;
-
- typedef typename GV::IntersectionIterator NeighborIterator;
- typedef HierarchicLevelIterator<GridType> HierarchicLevelIteratorType;
-
- typedef typename GridType::template Codim<0>::Entity Element;
-
- typedef typename Dune::MultipleCodimMultipleGeomTypeMapper<GV, FaceMapperLayout > FaceMapper;
-
- private:
- const int faultCount_;
- const std::vector<int>& faultToLevel_;
- const int coarseResolution_;
- const int maxLevel_;
- const int coarseGridN_;
-
- std::shared_ptr<GridType> grid_;
-
- std::shared_ptr<InterfaceNetwork<GridType>> interfaceNetwork_;
-
- std::vector<std::shared_ptr<FaultInterface<GV>>> faults_;
-
-private:
-
-typename std::enable_if<!std::numeric_limits<ctype>::is_integer, bool>::type
- almost_equal(ctype x, ctype y, int ulp)
- {
- return std::abs(x-y) < std::numeric_limits<ctype>::epsilon() * std::abs(x+y) * ulp
- || std::abs(x-y) < std::numeric_limits<ctype>::min();
- }
-
-// returns fault id
-int isFaultIntersection(const Intersection& intersection, const std::set<int>& interfaceYCoords) {
- Dune::FieldVector<ctype, dimworld> faceCenter = intersection.geometry().center();
-
- ctype yID = faceCenter[1]*coarseGridN_;
- ctype intpart;
-
- if (almost_equal(std::modf(yID, &intpart), 0.0, 2)) {
- std::set<int>::iterator it = interfaceYCoords.find((int) intpart);
- if (it!=interfaceYCoords.end())
- return std::distance(interfaceYCoords.begin(), it);
- else
- return -1;
- } else
- return -1;
-}
-
-
-public:
- //setup
- MultilevelUniformFaultFactory(const int faultCount, const std::vector<int>& faultToLevel, const int coarseResolution, const int maxLevel) :
- faultCount_(faultCount),
- faultToLevel_(faultToLevel),
- coarseResolution_(coarseResolution),
- maxLevel_(maxLevel),
- coarseGridN_(std::pow(2, coarseResolution_)),
- interfaceNetwork_(nullptr)
- {
- if (faultCount_>coarseGridN_-1)
- DUNE_THROW(Dune::Exception, "Coarse grid resolution too low to support " + std::to_string(faultCount_) + " fault interfaces!");
-
- assert((size_t) faultCount == faultToLevel_.size());
-
- Dune::UGGrid<dim>::setDefaultHeapSize(4000);
- GridOb unitCube(coarseGridN_);
- grid_ = unitCube.grid();
-
- grid_->globalRefine(maxLevel_);
-
- // init interface network
- interfaceNetwork_ = std::make_shared<InterfaceNetwork<GridType>>(*grid_);
-
- faults_.resize(faultCount);
- for (size_t i=0; i<faults_.size(); i++) {
- const int level = faultToLevel_[i];
- faults_[i] = std::make_shared<FaultInterface<GV>>(interfaceNetwork_->levelGridView(level), level);
- }
-
- std::set<int> interfaceYCoords;
- int yStep = coarseGridN_/(faultCount+1);
- int offSet = (coarseGridN_ - (faultCount-1)*yStep)/2;
- for (int i=0; i<faultCount; i++) {
- interfaceYCoords.insert(i*yStep + offSet);
- }
-
- std::cout << "yStep: " << yStep << " offSet: " << offSet << std::endl;
- print(interfaceYCoords, "interfaceYCoords: ");
-
- std::cout << "Coarse and fine grid were generated!" << std::endl;
-
-
- typedef typename GV::template Codim <0>::Iterator ElementLevelIterator;
-
- const GV& coarseGridView = grid_->levelGridView(0);
-
- // TODO: remove after debugging
- /*Dune::LevelMultipleCodimMultipleGeomTypeMapper <GridType, Dune::MCMGElementLayout > fineMapper(*grid_, fineLevelIdx);
- ElementLevelIterator fineElemIt = fineGridView.template begin<0>();
- ElementLevelIterator fineElemEndIt = fineGridView.template end<0>();
- for (; fineElemIt!=fineElemEndIt; ++fineElemIt) {
- const Element& fineElem = *fineElemIt;
- const auto& geometry = fineElem.geometry();
- std::cout << "Element idx: " << fineMapper.index(fineElem) << std::endl;
-
- for(int i=0; i<geometry.corners(); ++i) {
- size_t globalIdx = fineGridView.indexSet().subIndex(fineElem, i, Element::dimension);
- const auto& vertex = geometry.corner(i);
- print(vertex, "elem corner:" + std::to_string(globalIdx));
- }
- std::cout << "----------------" << std::endl << std::endl;
- }*/
- // end remove
-
- FaceMapper coarseFaceMapper(coarseGridView);
- std::vector<bool> coarseHandled(coarseFaceMapper.size(),false);
-
- ElementLevelIterator coarseElemIt = coarseGridView.template begin<0>();
- ElementLevelIterator coarseElemEndIt = coarseGridView.template end<0>();
- for (; coarseElemIt!=coarseElemEndIt; ++coarseElemIt) {
- const Element& coarseElem = *coarseElemIt;
-
- bool hasFaultIntersection = false;
-
- int faceID, faultID = -1;
-
- NeighborIterator coarseNIt = coarseGridView.ibegin(coarseElem);
- NeighborIterator coarseNEnd = coarseGridView.iend(coarseElem);
- for (; coarseNIt != coarseNEnd; ++coarseNIt){
-
- if (coarseHandled[coarseFaceMapper.subIndex(coarseElem,coarseNIt->indexInInside(),1)])
- continue;
-
- coarseHandled[coarseFaceMapper.subIndex(coarseElem,coarseNIt->indexInInside(),1)]=true;
-
- faultID = isFaultIntersection(*coarseNIt, interfaceYCoords);
-
- if (faultID >= 0) {
- hasFaultIntersection = true;
- faceID = coarseNIt->indexInInside();
- break;
- }
- }
-
- if (hasFaultIntersection) {
-
- Face<GridType> face(*grid_, coarseElem, faceID);
- const int faultLevel = faults_[faultID]->level();
-
- if (faultLevel == 0) {
- faults_[faultID]->addFace(*(face.intersection()));
- } else {
- typename Face<GridType>::HierarchicIterator exactIt = face.hbegin(faultLevel);
- typename Face<GridType>::HierarchicIterator exactEnd = face.hend(faultLevel);
-
- for (; exactIt!=exactEnd; ++exactIt) {
- if (exactIt->level() == faultLevel) {
- faults_[faultID]->addFace(*(exactIt->intersection()));
- }
- }
- }
- }
- }
-
- for (size_t i=0; i<faults_.size(); i++) {
- interfaceNetwork_->addInterface(faults_[i]);
- }
- }
-
- void prolongToAll() {
- // prolong all faults to all subsequent levels
- for (int i=maxLevel_-1; i>=0; i--) {
- if (interfaceNetwork_->size(i)>0) {
- std::set<int> toLevels;
- for (int j=i+1; j<=maxLevel_; j++) {
- toLevels.insert(j);
- }
-
- interfaceNetwork_->prolongLevels(i, toLevels);
- }
- }
- }
-
- /*void prolongToAll() {
- // prolong all faults to all subsequent levels
- for (int i=interfaceNetwork_->size()-1; i>=0; i--) {
- interfaceNetwork_->prolongLevelInterfaces(i, maxLevel_);
- }
- interfaceNetwork_->build();
- }*/
-
- const GridType& grid() const {
- return *grid_;
- }
-
- /*const InterfaceNetwork<GridType>& interfaceNetwork() {
- return *interfaceNetwork_;
- }*/
-
- InterfaceNetwork<GridType>& interfaceNetwork() {
- return *interfaceNetwork_;
- }
-};
-#endif
diff --git a/dune/faultnetworks/faultfactories/test/CMakeLists.txt b/dune/faultnetworks/faultfactories/test/CMakeLists.txt
deleted file mode 100644
index 012342bfba217823125c0b44af0d61ab07a0bbd2..0000000000000000000000000000000000000000
--- a/dune/faultnetworks/faultfactories/test/CMakeLists.txt
+++ /dev/null
@@ -1,12 +0,0 @@
-# Put your test in here if it needs access to external grids
-set(FAULTNETWORKS_FAULTFACTORY_TESTS
- multileveluniformfaultfactorytest
- )
-
-set(TESTS ${FAULTNETWORKS_FAULTFACTORY_TESTS})
-
-# Setup targets, register tests, and add dune flags
-foreach(_test ${TESTS})
- #dune_add_test(SOURCES ${_test}.cc)
-endforeach()
-
diff --git a/dune/faultnetworks/faultfactories/test/Makefile.am b/dune/faultnetworks/faultfactories/test/Makefile.am
deleted file mode 100644
index fad6cc38ace2886f0e29b0100b1516f65b754e83..0000000000000000000000000000000000000000
--- a/dune/faultnetworks/faultfactories/test/Makefile.am
+++ /dev/null
@@ -1,42 +0,0 @@
-# which tests where program to build and run are equal
-NORMALTESTS = multileveluniformfaultfactorytest
-
-# list of tests to run
-TESTS = $(NORMALTESTS)
-
-# programs just to build when "make check" is used
-check_PROGRAMS = $(NORMALTESTS)
-
-# define the programs
-multileveluniformfaultfactorytest_SOURCES = multileveluniformfaultfactorytest.cc
-
-
-multileveluniformfaultfactorytest_CPPFLAGS = $(AM_CPPFLAGS) \
- $(DUNEMPICPPFLAGS) \
- $(UG_CPPFLAGS)
-
-# The libraries have to be given in reverse order (most basic libraries
-# last). Also, due to some misunderstanding, a lot of libraries include the
-# -L option in LDFLAGS instead of LIBS -- so we have to include the LDFLAGS
-# here as well.
-multileveluniformfaultfactorytest_LDADD = \
- $(DUNE_LDFLAGS) $(DUNE_LIBS) \
- $(UG_LDFLAGS) $(UG_LIBS) \
- $(DUNEMPILIBS) \
- $(LDADD)
-multileveluniformfaultfactorytest_LDFLAGS = $(AM_LDFLAGS) \
- $(DUNEMPILDFLAGS) \
- $(UG_LDFLAGS) \
- $(DUNE_LDFLAGS)
-
-# don't follow the full GNU-standard
-# we need automake 1.9
-AUTOMAKE_OPTIONS = foreign 1.9
-
-# pass most important options when "make distcheck" is used
-DISTCHECK_CONFIGURE_FLAGS = --with-dune-common=$(DUNE_COMMON_ROOT) --with-dune-geometry=$(DUNE_GEOMETRY_ROOT) --with-dune-grid=$(DUNE_GRID_ROOT) --with-dune-istl=$(DUNE_ISTL_ROOT) --with-dune-localfunctions=$(DUNE_LOCALFUNCTIONS_ROOT) --with-dune-solvers=$(DUNE_SOLVERS_ROOT) --with-dune-fufem=$(DUNE_FUFEM_ROOT) CXX="$(CXX)" CC="$(CC)"
-
-EXTRA_DIST = CMakeLists.txt
-
-include $(top_srcdir)/am/global-rules
-
diff --git a/dune/faultnetworks/faultfactories/test/multileveluniformfaultfactorytest.cc b/dune/faultnetworks/faultfactories/test/multileveluniformfaultfactorytest.cc
deleted file mode 100644
index d9a6582d2e461773c829d3042af642439b866cbb..0000000000000000000000000000000000000000
--- a/dune/faultnetworks/faultfactories/test/multileveluniformfaultfactorytest.cc
+++ /dev/null
@@ -1,224 +0,0 @@
-#include <config.h>
-
-#include <dune/common/parallel/mpihelper.hh>
-
-#include <dune/grid/uggrid.hh>
-#include <dune/grid/common/mcmgmapper.hh>
-
-#include <dune/faultnetworks/interfacenetwork.hh>
-#include <dune/faultnetworks/faultfactories/multileveluniformfaultfactory.hh>
-#include <dune/faultnetworks/utils/debugutils.hh>
-
-const int dim = 2;
-
-const int coarseResolution = 2;
-const int fineResolution = 3;
-const int exactResolution = 4;
-
-const int faultCount = 3;
-
-typedef double ctype;
-typedef typename Dune::UGGrid<dim> GridType;
-typedef typename GridType::LevelGridView GridView;
-typedef typename GridType::LevelIntersection Intersection;
-typedef typename GridType::template Codim<0>::Entity Element;
-
-static const int dimElement = Element::dimension;
-
-void gridLevelOutput(const GridType& grid, const int level) {
- std::cout << "Grid elements on lvl " << level << std::endl << std::endl;
-
- typedef typename GridView::template Codim<0>::Iterator ElementLevelIterator;
- GridView gridView(grid.levelGridView(level));
-
- //Dune::LevelMultipleCodimMultipleGeomTypeMapper <GridType, Dune::MCMGElementLayout > mapper(grid, level);
- Dune::LevelMultipleCodimMultipleGeomTypeMapper <GridType, Dune::MCMGVertexLayout > vMapper(grid, level);
-
- Dune::BitSetVector<1> vertexPrinted(vMapper.size(), false);
-
- ElementLevelIterator elemIt = gridView.template begin<0>();
- ElementLevelIterator elemEndIt = gridView.template end<0>();
- for (; elemIt!=elemEndIt; ++elemIt) {
- const Element& elem = *elemIt;
- const auto& geometry = elem.geometry();
-
- for(int i=0; i<geometry.corners(); ++i) {
- size_t globalIdx = gridView.indexSet().subIndex(elem, i, Element::dimension);
- if (!vertexPrinted[globalIdx][0]) {
- const auto& vertex = geometry.corner(i);
- print(vertex, "elem corner:" + std::to_string(globalIdx));
- vertexPrinted[globalIdx][0] = true;
- }
- }
- }
- std::cout << "----------------" << std::endl << std::endl;
-
- /*
- ElementLevelIterator elemIt = gridView.template begin<0>();
- ElementLevelIterator elemEndIt = gridView.template end<0>();
- for (; elemIt!=elemEndIt; ++elemIt) {
- const Element& elem = *elemIt;
- const auto& geometry = elem.geometry();
-
- std::cout << "Element idx: " << mapper.index(elem) << std::endl;
-
- for(int i=0; i<geometry.corners(); ++i) {
- size_t globalIdx = gridView.indexSet().subIndex(elem, i, Element::dimension);
- const auto& vertex = geometry.corner(i);
- print(vertex, "elem corner:" + std::to_string(globalIdx));
- }
- std::cout << "----------------" << std::endl << std::endl;
- }*/
-
-}
-
-
-void faultInterfaceOutput(const FaultInterface<GridView>& faultInterface) {
- std::cout << "Level: " << faultInterface.level() << " global faultID: " << faultInterface.id() << std::endl;
-
- const std::set<size_t>& interfaceDofs = faultInterface.getInterfaceDofs();
- print(interfaceDofs, "");
-}
-
-void interfaceNetworkLevelOutput(const InterfaceNetwork<GridType>& interfaceNetwork, const int level) {
- const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork = interfaceNetwork.levelInterfaceNetwork(level);
-
- std::cout << "LevelInterfaceNetwork level: " << level << std::endl;
- for (size_t i=0; i<levelInterfaceNetwork.size(); i++) {
- std::cout << "local faultID : " << i << " ";
- faultInterfaceOutput(*(levelInterfaceNetwork.getInterface(i)));
- }
-
- std::cout << "------------------------------" << std::endl << std::endl;
-}
-
-void computeIntersectionDofs(const GridView& gridView, const Intersection& intersection, std::set<size_t>& intersectionDofs) {
- intersectionDofs.clear();
-
- // loop over all vertices of the intersection
- const Element& insideElement = intersection.inside();
-
- const auto& refElement = Dune::ReferenceElements<double,dimElement>::general(insideElement.type());
- for (int i=0; i<refElement.size(intersection.indexInInside(), 1, dimElement); i++) {
- size_t idxInElement = refElement.subEntity(intersection.indexInInside(), 1, i, dimElement);
- size_t globalIdx = gridView.indexSet().subIndex(insideElement, idxInElement, dimElement);
- intersectionDofs.insert(globalIdx);
- }
-}
-
-void prolongFaultDofs(const GridType& grid, const FaultInterface<GridView>& faultInterface, const int toLevel, std::set<size_t>& prolongedDofs) {
- prolongedDofs.clear();
-
- const std::vector<Intersection>& faces = faultInterface.faces();
- const GridView gridView(grid.levelGridView(toLevel));
-
- for (size_t i=0; i<faces.size(); i++) {
- const Intersection& intersection = faces[i];
- Face<GridType> face(grid, intersection.inside(), intersection.indexInInside());
-
- typename Face<GridType>::HierarchicIterator exactIt = face.hbegin(toLevel);
- typename Face<GridType>::HierarchicIterator exactEnd = face.hend(toLevel);
- for(; exactIt!=exactEnd; ++exactIt) {
- if (exactIt->level() == toLevel) {
- std::set<size_t> faceDofs;
- computeIntersectionDofs(gridView, *(exactIt->intersection()), faceDofs);
-
- prolongedDofs.insert(faceDofs.begin(), faceDofs.end());
- }
- }
- }
-}
-
-void prolongInterfaceNetworkLevelDofs(const GridType& grid, const InterfaceNetwork<GridType>& interfaceNetwork, const int level, const int toLevel, std::vector<std::set<size_t>>& prolongedFaultDofs) {
- prolongedFaultDofs.resize(0);
-
- const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork = interfaceNetwork.levelInterfaceNetwork(level);
-
- std::cout << "Correct prolonged LevelInterfaceNetwork level " << level << " to level " << toLevel << std::endl;
- for (size_t i=0; i<levelInterfaceNetwork.size(); i++) {
- std::set<size_t> prolongedDofs;
- prolongFaultDofs(grid, *(levelInterfaceNetwork.getInterface(i)), toLevel, prolongedDofs);
-
- print(prolongedDofs, "local faultID: " + std::to_string(i));
- prolongedFaultDofs.push_back(prolongedDofs);
- }
-
- std::cout << "------------------------------" << std::endl << std::endl;
-}
-
-int main(int argc, char** argv) try
-{
-
- Dune::MPIHelper::instance(argc, argv);
-
- std::ofstream out("/home/mi/podlesny/software/dune/faultnetworks/dune/faultnetworks/faultfactories/test/multileveluniformfaultfactorytest.log");
- std::streambuf *coutbuf = std::cout.rdbuf(); //save old buf
- std::cout.rdbuf(out.rdbuf()); //redirect std::cout to out.txt!
-
- const int fineLevelIdx = fineResolution - coarseResolution;
- const int exactLevelIdx = exactResolution - coarseResolution;
-
- std::vector<int> faultToLevel(faultCount);
- for (size_t i=0; i<faultToLevel.size(); ) {
- for (int j=fineLevelIdx; j>0 && i<faultToLevel.size(); j--) {
- faultToLevel[i] = j;
- i++;
- }
-
- for (int j=0; j<fineLevelIdx && i<faultToLevel.size(); j++) {
- faultToLevel[i] = j;
- i++;
- }
- }
-
- print(faultToLevel, "faultToLevel: ");
-
- MultilevelUniformFaultFactory<GridType> faultFactory(faultCount, faultToLevel, coarseResolution, exactLevelIdx);
- const GridType& grid = faultFactory.grid();
- const InterfaceNetwork<GridType>& interfaceNetwork = faultFactory.interfaceNetwork();
-
- gridLevelOutput(grid, 0);
- gridLevelOutput(grid, fineLevelIdx);
- gridLevelOutput(grid, exactLevelIdx);
-
- interfaceNetworkLevelOutput(interfaceNetwork, 0);
- interfaceNetworkLevelOutput(interfaceNetwork, fineLevelIdx);
- interfaceNetworkLevelOutput(interfaceNetwork, exactLevelIdx);
-
- std::vector<std::set<size_t>> correctProlongedFaultDofs;
- prolongInterfaceNetworkLevelDofs(grid, interfaceNetwork, 0, exactLevelIdx, correctProlongedFaultDofs);
- prolongInterfaceNetworkLevelDofs(grid, interfaceNetwork, fineLevelIdx, exactLevelIdx, correctProlongedFaultDofs);
-
-
- const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork = interfaceNetwork.levelInterfaceNetwork(exactLevelIdx);
-
- // build leaf view of all faults
- LevelInterfaceNetwork<GridView> leafInterfaceNetwork(levelInterfaceNetwork.grid(), levelInterfaceNetwork.level(), levelInterfaceNetwork.referencePoint());
-
- for (int i=0; i<=exactLevelIdx; i++) {
- interfaceNetwork.prolongLevelInterfaces(i, leafInterfaceNetwork);
- }
-
- std::cout << "Computed LeafInterfaceNetwork level " << exactLevelIdx << std::endl;
- for (size_t i=0; i<leafInterfaceNetwork.size(); i++) {
- print(leafInterfaceNetwork.getInterface(i)->getInterfaceDofs(), "local faultID: " + std::to_string(i));
- }
-
- int allProlongedFaultsCorrect = 0;
- for (size_t i=0; i<correctProlongedFaultDofs.size(); i++) {
- for (size_t j=0; j<leafInterfaceNetwork.size(); j++) {
- const std::set<size_t>& faultDofs = leafInterfaceNetwork.getInterface(j)->getInterfaceDofs();
-
- if (correctProlongedFaultDofs[i]==faultDofs) {
- allProlongedFaultsCorrect++;
- break;
- }
- }
- }
-
- assert(allProlongedFaultsCorrect == (int) correctProlongedFaultDofs.size());
-
- std::cout.rdbuf(coutbuf);
-} catch (Dune::Exception e){
- std::cout << e << std::endl;
-}
diff --git a/dune/faultnetworks/faultfactories/test/multileveluniformfaultfactorytest.log b/dune/faultnetworks/faultfactories/test/multileveluniformfaultfactorytest.log
deleted file mode 100644
index 8a6552ccd3a95cf4f0689321742ba44e7454233f..0000000000000000000000000000000000000000
--- a/dune/faultnetworks/faultfactories/test/multileveluniformfaultfactorytest.log
+++ /dev/null
@@ -1,1247 +0,0 @@
-faultToLevel:
-1 0 1
-
-Coarse and fine grid were generated!
-Grid elements on lvl 0
-
-elem corner:0
-0 0
-
-elem corner:1
-0.25 0
-
-elem corner:6
-0.25 0.25
-
-elem corner:5
-0 0.25
-
-elem corner:2
-0.5 0
-
-elem corner:7
-0.5 0.25
-
-elem corner:3
-0.75 0
-
-elem corner:8
-0.75 0.25
-
-elem corner:4
-1 0
-
-elem corner:9
-1 0.25
-
-elem corner:11
-0.25 0.5
-
-elem corner:10
-0 0.5
-
-elem corner:12
-0.5 0.5
-
-elem corner:13
-0.75 0.5
-
-elem corner:14
-1 0.5
-
-elem corner:16
-0.25 0.75
-
-elem corner:15
-0 0.75
-
-elem corner:17
-0.5 0.75
-
-elem corner:18
-0.75 0.75
-
-elem corner:19
-1 0.75
-
-elem corner:21
-0.25 1
-
-elem corner:20
-0 1
-
-elem corner:22
-0.5 1
-
-elem corner:23
-0.75 1
-
-elem corner:24
-1 1
-
-----------------
-
-Grid elements on lvl 1
-
-elem corner:0
-0 0
-
-elem corner:3
-0.125 0
-
-elem corner:5
-0.125 0.125
-
-elem corner:1
-0.25 0
-
-elem corner:4
-0.25 0.125
-
-elem corner:2
-0.25 0.25
-
-elem corner:6
-0 0.25
-
-elem corner:7
-0 0.125
-
-elem corner:8
-0.125 0.25
-
-elem corner:11
-0.375 0
-
-elem corner:13
-0.375 0.125
-
-elem corner:9
-0.5 0
-
-elem corner:12
-0.5 0.125
-
-elem corner:10
-0.5 0.25
-
-elem corner:14
-0.375 0.25
-
-elem corner:17
-0.625 0
-
-elem corner:19
-0.625 0.125
-
-elem corner:15
-0.75 0
-
-elem corner:18
-0.75 0.125
-
-elem corner:16
-0.75 0.25
-
-elem corner:20
-0.625 0.25
-
-elem corner:23
-0.875 0
-
-elem corner:25
-0.875 0.125
-
-elem corner:21
-1 0
-
-elem corner:24
-1 0.125
-
-elem corner:22
-1 0.25
-
-elem corner:26
-0.875 0.25
-
-elem corner:29
-0.125 0.375
-
-elem corner:28
-0.25 0.375
-
-elem corner:27
-0.25 0.5
-
-elem corner:30
-0 0.5
-
-elem corner:31
-0 0.375
-
-elem corner:32
-0.125 0.5
-
-elem corner:35
-0.375 0.375
-
-elem corner:34
-0.5 0.375
-
-elem corner:33
-0.5 0.5
-
-elem corner:36
-0.375 0.5
-
-elem corner:39
-0.625 0.375
-
-elem corner:38
-0.75 0.375
-
-elem corner:37
-0.75 0.5
-
-elem corner:40
-0.625 0.5
-
-elem corner:43
-0.875 0.375
-
-elem corner:42
-1 0.375
-
-elem corner:41
-1 0.5
-
-elem corner:44
-0.875 0.5
-
-elem corner:47
-0.125 0.625
-
-elem corner:46
-0.25 0.625
-
-elem corner:45
-0.25 0.75
-
-elem corner:48
-0 0.75
-
-elem corner:49
-0 0.625
-
-elem corner:50
-0.125 0.75
-
-elem corner:53
-0.375 0.625
-
-elem corner:52
-0.5 0.625
-
-elem corner:51
-0.5 0.75
-
-elem corner:54
-0.375 0.75
-
-elem corner:57
-0.625 0.625
-
-elem corner:56
-0.75 0.625
-
-elem corner:55
-0.75 0.75
-
-elem corner:58
-0.625 0.75
-
-elem corner:61
-0.875 0.625
-
-elem corner:60
-1 0.625
-
-elem corner:59
-1 0.75
-
-elem corner:62
-0.875 0.75
-
-elem corner:65
-0.125 0.875
-
-elem corner:64
-0.25 0.875
-
-elem corner:63
-0.25 1
-
-elem corner:66
-0 1
-
-elem corner:67
-0 0.875
-
-elem corner:68
-0.125 1
-
-elem corner:71
-0.375 0.875
-
-elem corner:70
-0.5 0.875
-
-elem corner:69
-0.5 1
-
-elem corner:72
-0.375 1
-
-elem corner:75
-0.625 0.875
-
-elem corner:74
-0.75 0.875
-
-elem corner:73
-0.75 1
-
-elem corner:76
-0.625 1
-
-elem corner:79
-0.875 0.875
-
-elem corner:78
-1 0.875
-
-elem corner:77
-1 1
-
-elem corner:80
-0.875 1
-
-----------------
-
-Grid elements on lvl 2
-
-elem corner:0
-0 0
-
-elem corner:3
-0.0625 0
-
-elem corner:5
-0.0625 0.0625
-
-elem corner:1
-0.125 0
-
-elem corner:4
-0.125 0.0625
-
-elem corner:2
-0.125 0.125
-
-elem corner:8
-0.1875 0
-
-elem corner:10
-0.1875 0.0625
-
-elem corner:6
-0.25 0
-
-elem corner:9
-0.25 0.0625
-
-elem corner:7
-0.25 0.125
-
-elem corner:12
-0.25 0.1875
-
-elem corner:14
-0.1875 0.125
-
-elem corner:11
-0.25 0.25
-
-elem corner:13
-0.1875 0.1875
-
-elem corner:15
-0 0.25
-
-elem corner:18
-0 0.1875
-
-elem corner:20
-0.0625 0.25
-
-elem corner:16
-0 0.125
-
-elem corner:19
-0.0625 0.1875
-
-elem corner:17
-0.125 0.25
-
-elem corner:21
-0 0.0625
-
-elem corner:22
-0.0625 0.125
-
-elem corner:24
-0.125 0.1875
-
-elem corner:23
-0.1875 0.25
-
-elem corner:27
-0.3125 0
-
-elem corner:29
-0.3125 0.0625
-
-elem corner:25
-0.375 0
-
-elem corner:28
-0.375 0.0625
-
-elem corner:26
-0.375 0.125
-
-elem corner:32
-0.4375 0
-
-elem corner:34
-0.4375 0.0625
-
-elem corner:30
-0.5 0
-
-elem corner:33
-0.5 0.0625
-
-elem corner:31
-0.5 0.125
-
-elem corner:36
-0.5 0.1875
-
-elem corner:38
-0.4375 0.125
-
-elem corner:35
-0.5 0.25
-
-elem corner:37
-0.4375 0.1875
-
-elem corner:41
-0.3125 0.25
-
-elem corner:40
-0.3125 0.1875
-
-elem corner:39
-0.375 0.25
-
-elem corner:42
-0.3125 0.125
-
-elem corner:44
-0.375 0.1875
-
-elem corner:43
-0.4375 0.25
-
-elem corner:47
-0.5625 0
-
-elem corner:49
-0.5625 0.0625
-
-elem corner:45
-0.625 0
-
-elem corner:48
-0.625 0.0625
-
-elem corner:46
-0.625 0.125
-
-elem corner:52
-0.6875 0
-
-elem corner:54
-0.6875 0.0625
-
-elem corner:50
-0.75 0
-
-elem corner:53
-0.75 0.0625
-
-elem corner:51
-0.75 0.125
-
-elem corner:56
-0.75 0.1875
-
-elem corner:58
-0.6875 0.125
-
-elem corner:55
-0.75 0.25
-
-elem corner:57
-0.6875 0.1875
-
-elem corner:61
-0.5625 0.25
-
-elem corner:60
-0.5625 0.1875
-
-elem corner:59
-0.625 0.25
-
-elem corner:62
-0.5625 0.125
-
-elem corner:64
-0.625 0.1875
-
-elem corner:63
-0.6875 0.25
-
-elem corner:67
-0.8125 0
-
-elem corner:69
-0.8125 0.0625
-
-elem corner:65
-0.875 0
-
-elem corner:68
-0.875 0.0625
-
-elem corner:66
-0.875 0.125
-
-elem corner:72
-0.9375 0
-
-elem corner:74
-0.9375 0.0625
-
-elem corner:70
-1 0
-
-elem corner:73
-1 0.0625
-
-elem corner:71
-1 0.125
-
-elem corner:76
-1 0.1875
-
-elem corner:78
-0.9375 0.125
-
-elem corner:75
-1 0.25
-
-elem corner:77
-0.9375 0.1875
-
-elem corner:81
-0.8125 0.25
-
-elem corner:80
-0.8125 0.1875
-
-elem corner:79
-0.875 0.25
-
-elem corner:82
-0.8125 0.125
-
-elem corner:84
-0.875 0.1875
-
-elem corner:83
-0.9375 0.25
-
-elem corner:87
-0.0625 0.3125
-
-elem corner:86
-0.125 0.3125
-
-elem corner:85
-0.125 0.375
-
-elem corner:90
-0.1875 0.3125
-
-elem corner:89
-0.25 0.3125
-
-elem corner:88
-0.25 0.375
-
-elem corner:92
-0.25 0.4375
-
-elem corner:94
-0.1875 0.375
-
-elem corner:91
-0.25 0.5
-
-elem corner:93
-0.1875 0.4375
-
-elem corner:95
-0 0.5
-
-elem corner:98
-0 0.4375
-
-elem corner:100
-0.0625 0.5
-
-elem corner:96
-0 0.375
-
-elem corner:99
-0.0625 0.4375
-
-elem corner:97
-0.125 0.5
-
-elem corner:101
-0 0.3125
-
-elem corner:102
-0.0625 0.375
-
-elem corner:104
-0.125 0.4375
-
-elem corner:103
-0.1875 0.5
-
-elem corner:107
-0.3125 0.3125
-
-elem corner:106
-0.375 0.3125
-
-elem corner:105
-0.375 0.375
-
-elem corner:110
-0.4375 0.3125
-
-elem corner:109
-0.5 0.3125
-
-elem corner:108
-0.5 0.375
-
-elem corner:112
-0.5 0.4375
-
-elem corner:114
-0.4375 0.375
-
-elem corner:111
-0.5 0.5
-
-elem corner:113
-0.4375 0.4375
-
-elem corner:117
-0.3125 0.5
-
-elem corner:116
-0.3125 0.4375
-
-elem corner:115
-0.375 0.5
-
-elem corner:118
-0.3125 0.375
-
-elem corner:120
-0.375 0.4375
-
-elem corner:119
-0.4375 0.5
-
-elem corner:123
-0.5625 0.3125
-
-elem corner:122
-0.625 0.3125
-
-elem corner:121
-0.625 0.375
-
-elem corner:126
-0.6875 0.3125
-
-elem corner:125
-0.75 0.3125
-
-elem corner:124
-0.75 0.375
-
-elem corner:128
-0.75 0.4375
-
-elem corner:130
-0.6875 0.375
-
-elem corner:127
-0.75 0.5
-
-elem corner:129
-0.6875 0.4375
-
-elem corner:133
-0.5625 0.5
-
-elem corner:132
-0.5625 0.4375
-
-elem corner:131
-0.625 0.5
-
-elem corner:134
-0.5625 0.375
-
-elem corner:136
-0.625 0.4375
-
-elem corner:135
-0.6875 0.5
-
-elem corner:139
-0.8125 0.3125
-
-elem corner:138
-0.875 0.3125
-
-elem corner:137
-0.875 0.375
-
-elem corner:142
-0.9375 0.3125
-
-elem corner:141
-1 0.3125
-
-elem corner:140
-1 0.375
-
-elem corner:144
-1 0.4375
-
-elem corner:146
-0.9375 0.375
-
-elem corner:143
-1 0.5
-
-elem corner:145
-0.9375 0.4375
-
-elem corner:149
-0.8125 0.5
-
-elem corner:148
-0.8125 0.4375
-
-elem corner:147
-0.875 0.5
-
-elem corner:150
-0.8125 0.375
-
-elem corner:152
-0.875 0.4375
-
-elem corner:151
-0.9375 0.5
-
-elem corner:155
-0.0625 0.5625
-
-elem corner:154
-0.125 0.5625
-
-elem corner:153
-0.125 0.625
-
-elem corner:158
-0.1875 0.5625
-
-elem corner:157
-0.25 0.5625
-
-elem corner:156
-0.25 0.625
-
-elem corner:160
-0.25 0.6875
-
-elem corner:162
-0.1875 0.625
-
-elem corner:159
-0.25 0.75
-
-elem corner:161
-0.1875 0.6875
-
-elem corner:163
-0 0.75
-
-elem corner:166
-0 0.6875
-
-elem corner:168
-0.0625 0.75
-
-elem corner:164
-0 0.625
-
-elem corner:167
-0.0625 0.6875
-
-elem corner:165
-0.125 0.75
-
-elem corner:169
-0 0.5625
-
-elem corner:170
-0.0625 0.625
-
-elem corner:172
-0.125 0.6875
-
-elem corner:171
-0.1875 0.75
-
-elem corner:175
-0.3125 0.5625
-
-elem corner:174
-0.375 0.5625
-
-elem corner:173
-0.375 0.625
-
-elem corner:178
-0.4375 0.5625
-
-elem corner:177
-0.5 0.5625
-
-elem corner:176
-0.5 0.625
-
-elem corner:180
-0.5 0.6875
-
-elem corner:182
-0.4375 0.625
-
-elem corner:179
-0.5 0.75
-
-elem corner:181
-0.4375 0.6875
-
-elem corner:185
-0.3125 0.75
-
-elem corner:184
-0.3125 0.6875
-
-elem corner:183
-0.375 0.75
-
-elem corner:186
-0.3125 0.625
-
-elem corner:188
-0.375 0.6875
-
-elem corner:187
-0.4375 0.75
-
-elem corner:191
-0.5625 0.5625
-
-elem corner:190
-0.625 0.5625
-
-elem corner:189
-0.625 0.625
-
-elem corner:194
-0.6875 0.5625
-
-elem corner:193
-0.75 0.5625
-
-elem corner:192
-0.75 0.625
-
-elem corner:196
-0.75 0.6875
-
-elem corner:198
-0.6875 0.625
-
-elem corner:195
-0.75 0.75
-
-elem corner:197
-0.6875 0.6875
-
-elem corner:201
-0.5625 0.75
-
-elem corner:200
-0.5625 0.6875
-
-elem corner:199
-0.625 0.75
-
-elem corner:202
-0.5625 0.625
-
-elem corner:204
-0.625 0.6875
-
-elem corner:203
-0.6875 0.75
-
-elem corner:207
-0.8125 0.5625
-
-elem corner:206
-0.875 0.5625
-
-elem corner:205
-0.875 0.625
-
-elem corner:210
-0.9375 0.5625
-
-elem corner:209
-1 0.5625
-
-elem corner:208
-1 0.625
-
-elem corner:212
-1 0.6875
-
-elem corner:214
-0.9375 0.625
-
-elem corner:211
-1 0.75
-
-elem corner:213
-0.9375 0.6875
-
-elem corner:217
-0.8125 0.75
-
-elem corner:216
-0.8125 0.6875
-
-elem corner:215
-0.875 0.75
-
-elem corner:218
-0.8125 0.625
-
-elem corner:220
-0.875 0.6875
-
-elem corner:219
-0.9375 0.75
-
-elem corner:223
-0.0625 0.8125
-
-elem corner:222
-0.125 0.8125
-
-elem corner:221
-0.125 0.875
-
-elem corner:226
-0.1875 0.8125
-
-elem corner:225
-0.25 0.8125
-
-elem corner:224
-0.25 0.875
-
-elem corner:228
-0.25 0.9375
-
-elem corner:230
-0.1875 0.875
-
-elem corner:227
-0.25 1
-
-elem corner:229
-0.1875 0.9375
-
-elem corner:231
-0 1
-
-elem corner:234
-0 0.9375
-
-elem corner:236
-0.0625 1
-
-elem corner:232
-0 0.875
-
-elem corner:235
-0.0625 0.9375
-
-elem corner:233
-0.125 1
-
-elem corner:237
-0 0.8125
-
-elem corner:238
-0.0625 0.875
-
-elem corner:240
-0.125 0.9375
-
-elem corner:239
-0.1875 1
-
-elem corner:243
-0.3125 0.8125
-
-elem corner:242
-0.375 0.8125
-
-elem corner:241
-0.375 0.875
-
-elem corner:246
-0.4375 0.8125
-
-elem corner:245
-0.5 0.8125
-
-elem corner:244
-0.5 0.875
-
-elem corner:248
-0.5 0.9375
-
-elem corner:250
-0.4375 0.875
-
-elem corner:247
-0.5 1
-
-elem corner:249
-0.4375 0.9375
-
-elem corner:253
-0.3125 1
-
-elem corner:252
-0.3125 0.9375
-
-elem corner:251
-0.375 1
-
-elem corner:254
-0.3125 0.875
-
-elem corner:256
-0.375 0.9375
-
-elem corner:255
-0.4375 1
-
-elem corner:259
-0.5625 0.8125
-
-elem corner:258
-0.625 0.8125
-
-elem corner:257
-0.625 0.875
-
-elem corner:262
-0.6875 0.8125
-
-elem corner:261
-0.75 0.8125
-
-elem corner:260
-0.75 0.875
-
-elem corner:264
-0.75 0.9375
-
-elem corner:266
-0.6875 0.875
-
-elem corner:263
-0.75 1
-
-elem corner:265
-0.6875 0.9375
-
-elem corner:269
-0.5625 1
-
-elem corner:268
-0.5625 0.9375
-
-elem corner:267
-0.625 1
-
-elem corner:270
-0.5625 0.875
-
-elem corner:272
-0.625 0.9375
-
-elem corner:271
-0.6875 1
-
-elem corner:275
-0.8125 0.8125
-
-elem corner:274
-0.875 0.8125
-
-elem corner:273
-0.875 0.875
-
-elem corner:278
-0.9375 0.8125
-
-elem corner:277
-1 0.8125
-
-elem corner:276
-1 0.875
-
-elem corner:280
-1 0.9375
-
-elem corner:282
-0.9375 0.875
-
-elem corner:279
-1 1
-
-elem corner:281
-0.9375 0.9375
-
-elem corner:285
-0.8125 1
-
-elem corner:284
-0.8125 0.9375
-
-elem corner:283
-0.875 1
-
-elem corner:286
-0.8125 0.875
-
-elem corner:288
-0.875 0.9375
-
-elem corner:287
-0.9375 1
-
-----------------
-
-LevelInterfaceNetwork level: 0
-local faultID : 0 Level: 0 global faultID: 1
-
-10 11 12 13 14
-
-------------------------------
-
-LevelInterfaceNetwork level: 1
-local faultID : 0 Level: 1 global faultID: 0
-
-2 6 8 10 14 16 20 22 26
-
-local faultID : 1 Level: 1 global faultID: 2
-
-45 48 50 51 54 55 58 59 62
-
-------------------------------
-
-LevelInterfaceNetwork level: 2
-------------------------------
-
-Correct prolonged LevelInterfaceNetwork level 0 to level 2
-local faultID: 0
-91 95 97 100 103 111 115 117 119 127 131 133 135 143 147 149 151
-
-------------------------------
-
-Correct prolonged LevelInterfaceNetwork level 1 to level 2
-local faultID: 0
-11 15 17 20 23 35 39 41 43 55 59 61 63 75 79 81 83
-
-local faultID: 1
-159 163 165 168 171 179 183 185 187 195 199 201 203 211 215 217 219
-
-------------------------------
-
-Computed LeafInterfaceNetwork level 2
-local faultID: 0
-11 15 17 20 23 35 39 41 43 55 59 61 63 75 79 81 83
-
-local faultID: 1
-91 95 97 100 103 111 115 117 119 127 131 133 135 143 147 149 151
-
-local faultID: 2
-159 163 165 168 171 179 183 185 187 195 199 201 203 211 215 217 219
-
diff --git a/dune/faultnetworks/faultfactories/twofaultfactory.hh b/dune/faultnetworks/faultfactories/twofaultfactory.hh
deleted file mode 100644
index c201968cabbcb12384f216e8f949d8b2fa953ee2..0000000000000000000000000000000000000000
--- a/dune/faultnetworks/faultfactories/twofaultfactory.hh
+++ /dev/null
@@ -1,261 +0,0 @@
-#ifndef TWO_FAULT_FACTORY_HH
-#define TWO_FAULT_FACTORY_HH
-
-#include <math.h>
-
-#include <dune/faultnetworks/facehierarchy.hh>
-#include <dune/faultnetworks/utils/debugutils.hh>
-#include <dune/faultnetworks/faultfactories/oscunitcube.hh>
-#include <dune/faultnetworks/levelinterfacenetwork.hh>
-#include <dune/faultnetworks/faultinterface.hh>
-#include <dune/faultnetworks/hierarchicleveliterator.hh>
-
-#include <dune/grid/common/mcmgmapper.hh>
-#include <dune/grid/io/file/vtk/vtkwriter.hh>
-
-
-template <class GridType>
-class TwoFaultFactory
-{
- //! Parameter for mapper class
- template<int dim>
- struct FaceMapperLayout
- {
- bool contains (Dune::GeometryType gt)
- {
- return gt.dim() == dim-1;
- }
- };
-
- protected:
- typedef OscUnitCube<GridType, 2> GridOb;
- static const int dimworld = GridType::dimensionworld ;
- static const int dim = GridType::dimension;
- static const int faultCount = 2;
- typedef typename GridType::ctype ctype;
-
- typedef typename GridType::LevelGridView GV;
- typedef typename GridType::LevelIntersection Intersection;
- typedef typename Intersection::Geometry IntersectionGeometry;
-
- typedef typename GV::IntersectionIterator NeighborIterator;
- typedef HierarchicLevelIterator<GridType> HierarchicLevelIteratorType;
-
- typedef typename GridType::template Codim<0>::Entity Element;
-
- typedef typename Dune::MultipleCodimMultipleGeomTypeMapper<GV, FaceMapperLayout > FaceMapper;
-
- private:
- const int faultDistance_;
- const int coarseResolution_;
- const int fineResolution_;
- const int exactResolution_;
-
- int coarseGridN_;
-
- std::shared_ptr<GridType> grid_;
-
- std::shared_ptr<LevelInterfaceNetwork<GV>> coarseInterfaceNetwork_;
- std::shared_ptr<LevelInterfaceNetwork<GV>> fineInterfaceNetwork_;
- std::shared_ptr<LevelInterfaceNetwork<GV>> exactInterfaceNetwork_;
-
-private:
-
-typename std::enable_if<!std::numeric_limits<ctype>::is_integer, bool>::type
- almost_equal(ctype x, ctype y, int ulp)
- {
- return std::abs(x-y) < std::numeric_limits<ctype>::epsilon() * std::abs(x+y) * ulp
- || std::abs(x-y) < std::numeric_limits<ctype>::min();
- }
-
-// returns fault id
-int isFaultIntersection(const Intersection& intersection, const std::set<int>& interfaceYCoords) {
- Dune::FieldVector<ctype, dimworld> faceCenter = intersection.geometry().center();
-
- ctype yID = faceCenter[1]*coarseGridN_;
- ctype intpart;
-
- if (almost_equal(std::modf(yID, &intpart), 0.0, 2)) {
- std::set<int>::iterator it = interfaceYCoords.find((int) intpart);
- if (it!=interfaceYCoords.end())
- return std::distance(interfaceYCoords.begin(), it);
- else
- return -1;
- } else
- return -1;
-}
-
-
- public:
- //setup
- TwoFaultFactory(const int faultDistance, const int coarseResolution, const int fineResolution, const int exactResolution) :
- faultDistance_(faultDistance),
- coarseResolution_(coarseResolution),
- fineResolution_(fineResolution),
- exactResolution_(exactResolution),
- coarseGridN_(std::pow(2, coarseResolution_)),
- coarseInterfaceNetwork_(nullptr),
- fineInterfaceNetwork_(nullptr),
- exactInterfaceNetwork_(nullptr)
- {
- if (faultDistance_<1)
- DUNE_THROW(Dune::Exception, "faultDistance has to be >0");
- if (faultDistance_>coarseGridN_-2)
- DUNE_THROW(Dune::Exception, "Coarse grid resolution too low to support fault distance of " + std::to_string(faultDistance_));
-
- Dune::UGGrid<dim>::setDefaultHeapSize(4000);
- GridOb unitCube(coarseGridN_);
- grid_ = unitCube.grid();
-
- const int fineLevelIdx = fineResolution_ - coarseResolution_;
- const int exactLevelIdx = exactResolution_ - coarseResolution_;
- grid_->globalRefine(exactLevelIdx);
-
- std::set<int> interfaceYCoords;
- int yStep = (coarseGridN_-faultDistance_)/2;
- interfaceYCoords.insert(yStep);
- interfaceYCoords.insert(yStep+faultDistance_);
-
- std::cout << "Coarse and fine grid were generated!" << std::endl;
-
- /*
- GV gridView(grid.levelGridView(0));
-
-
- VTKWriter<GV> vtkwriter(gridView);
- vtkwriter.pwrite("coarseGrid", "/home/mi/podlesjo/data/", "/home/mi/podlesjo/data/");//("~/home/mi/podlesjo/data/coarseGrid");
-
- GV gridViewFine(grid.levelGridView(fineLevelIdx));
-
-
- VTKWriter<GV> vtkwriterFine(gridViewFine);
- vtkwriterFine.pwrite("fineGrid", "/home/mi/podlesjo/data/", "/home/mi/podlesjo/data/");//("~/home/mi/podlesjo/data/coarseGrid");
- */
-
- typedef typename GV::template Codim <0>::Iterator ElementLevelIterator;
-
- GV coarseGridView(grid_->levelGridView(0));
- GV fineGridView(grid_->levelGridView(fineLevelIdx));
- GV exactGridView(grid_->levelGridView(exactLevelIdx));
-
- coarseInterfaceNetwork_ = std::make_shared<LevelInterfaceNetwork<GV>>(*grid_, 0);
- fineInterfaceNetwork_ = std::make_shared<LevelInterfaceNetwork<GV>>(*grid_, fineLevelIdx);
- exactInterfaceNetwork_ = std::make_shared<LevelInterfaceNetwork<GV>>(*grid_, exactLevelIdx);
-
- // init interfaces
- std::vector<std::shared_ptr<FaultInterface<GV>>> coarseInterfaces(faultCount);
- std::vector<std::shared_ptr<FaultInterface<GV>>> fineInterfaces(faultCount);
- std::vector<std::shared_ptr<FaultInterface<GV>>> exactInterfaces(faultCount);
- for (int i=0; i<faultCount; i++) {
- coarseInterfaces[i] = std::make_shared<FaultInterface<GV>>(FaultInterface<GV>(coarseGridView, 0));
- fineInterfaces[i] = std::make_shared<FaultInterface<GV>>(FaultInterface<GV>(fineGridView, fineLevelIdx));
- exactInterfaces[i] = std::make_shared<FaultInterface<GV>>(FaultInterface<GV>(exactGridView, exactLevelIdx));
- }
-
-
- /*
- VTKWriter<GV> vtkwriter(fineGridView);
- vtkwriter.pwrite("fineGrid", "/home/mi/podlesny/data/", "/home/mi/podlesny/data/");//("~/home/mi/podlesny/data/fineGrid");
-
- VTKWriter<GV> vtkwriterFine(exactGridView);
- vtkwriterFine.pwrite("exactGrid", "/home/mi/podlesny/data/", "/home/mi/podlesny/data/");//("~/home/mi/podlesny/data/exactGrid");
- */
-
- // TODO: remove after debugging
- /*Dune::LevelMultipleCodimMultipleGeomTypeMapper <GridType, Dune::MCMGElementLayout > fineMapper(*grid_, fineLevelIdx);
- ElementLevelIterator fineElemIt = fineGridView.template begin<0>();
- ElementLevelIterator fineElemEndIt = fineGridView.template end<0>();
- for (; fineElemIt!=fineElemEndIt; ++fineElemIt) {
- const Element& fineElem = *fineElemIt;
- const auto& geometry = fineElem.geometry();
- std::cout << "Element idx: " << fineMapper.index(fineElem) << std::endl;
-
- for(int i=0; i<geometry.corners(); ++i) {
- size_t globalIdx = fineGridView.indexSet().subIndex(fineElem, i, Element::dimension);
- const auto& vertex = geometry.corner(i);
- print(vertex, "elem corner:" + std::to_string(globalIdx));
- }
- std::cout << "----------------" << std::endl << std::endl;
- }*/
- // end remove
-
-
- FaceMapper coarseFaceMapper(coarseGridView);
- std::vector<bool> coarseHandled(coarseFaceMapper.size(),false);
-
- ElementLevelIterator coarseElemIt = coarseGridView.template begin<0>();
- ElementLevelIterator coarseElemEndIt = coarseGridView.template end<0>();
- for (; coarseElemIt!=coarseElemEndIt; ++coarseElemIt) {
- const Element& coarseElem = *coarseElemIt;
-
- bool hasFaultIntersection = false;
-
- int faceID, faultID = -1;
-
- NeighborIterator coarseNIt = coarseGridView.ibegin(coarseElem);
- NeighborIterator coarseNEnd = coarseGridView.iend(coarseElem);
- for (; coarseNIt != coarseNEnd; ++coarseNIt){
-
- if (coarseHandled[coarseFaceMapper.subIndex(coarseElem,coarseNIt->indexInInside(),1)])
- continue;
-
- coarseHandled[coarseFaceMapper.subIndex(coarseElem,coarseNIt->indexInInside(),1)]=true;
-
- faultID = isFaultIntersection(*coarseNIt, interfaceYCoords);
-
- if (faultID >= 0) {
- hasFaultIntersection = true;
- faceID = coarseNIt->indexInInside();
- break;
- }
- }
-
- if (hasFaultIntersection) {
- Face<GridType> face(*grid_, coarseElem, faceID);
-
- typename Face<GridType>::HierarchicIterator exactIt = face.hbegin(exactLevelIdx);
- typename Face<GridType>::HierarchicIterator exactEnd = face.hend(exactLevelIdx);
- for(; exactIt!=exactEnd; ++exactIt) {
- if (exactIt->level()==exactLevelIdx)
- exactInterfaces[faultID]->addFace(*(exactIt->intersection()));
- else if (exactIt->level()==fineLevelIdx)
- fineInterfaces[faultID]->addFace(*(exactIt->intersection()));
- else if (exactIt->level()==0)
- coarseInterfaces[faultID]->addFace(*(exactIt->intersection()));
- }
- }
- }
-
- for (int i=0; i<faultCount; i++) {
- coarseInterfaces[i]->build();
- coarseInterfaceNetwork_->addInterface(coarseInterfaces[i]);
-
- fineInterfaces[i]->build();
- fineInterfaceNetwork_->addInterface(fineInterfaces[i]);
- //print(fineInterfaceNetwork_->getInterfaceNetworkDofs(), "fineInterface dofs:");
-
- exactInterfaces[i]->build();
- exactInterfaceNetwork_->addInterface(exactInterfaces[i]);
- //print(exactInterfaceNetwork_->getInterfaceNetworkDofs(), "exactInterface dofs:");
- }
- }
-
- const GridType& grid() const {
- return *grid_;
- }
-
- const LevelInterfaceNetwork<GV>& coarseInterfaceNetwork() {
- return *coarseInterfaceNetwork_;
- }
-
-
- const LevelInterfaceNetwork<GV>& fineInterfaceNetwork() {
- return *fineInterfaceNetwork_;
- }
-
- const LevelInterfaceNetwork<GV>& exactInterfaceNetwork() {
- return *exactInterfaceNetwork_;
- }
-
-};
-#endif
diff --git a/dune/faultnetworks/faultfactories/uniformfaultfactory.hh b/dune/faultnetworks/faultfactories/uniformfaultfactory.hh
deleted file mode 100644
index c66d72ddca2bffe5372203f99b5419a425d7ccfb..0000000000000000000000000000000000000000
--- a/dune/faultnetworks/faultfactories/uniformfaultfactory.hh
+++ /dev/null
@@ -1,259 +0,0 @@
-#ifndef UNIFORM_FAULT_FACTORY_HH
-#define UNIFORM_FAULT_FACTORY_HH
-
-#include <math.h>
-
-#include <dune/faultnetworks/facehierarchy.hh>
-#include <dune/faultnetworks/utils/debugutils.hh>
-#include <dune/faultnetworks/faultfactories/oscunitcube.hh>
-#include <dune/faultnetworks/levelinterfacenetwork.hh>
-#include <dune/faultnetworks/faultinterface.hh>
-#include <dune/faultnetworks/hierarchicleveliterator.hh>
-
-#include <dune/grid/common/mcmgmapper.hh>
-#include <dune/grid/io/file/vtk/vtkwriter.hh>
-
-
-template <class GridType>
-class UniformFaultFactory
-{
- //! Parameter for mapper class
- template<int dim>
- struct FaceMapperLayout
- {
- bool contains (Dune::GeometryType gt)
- {
- return gt.dim() == dim-1;
- }
- };
-
- protected:
- typedef OscUnitCube<GridType, 2> GridOb;
- static const int dimworld = GridType::dimensionworld ;
- static const int dim = GridType::dimension;
- typedef typename GridType::ctype ctype;
-
- typedef typename GridType::LevelGridView GV;
- typedef typename GridType::LevelIntersection Intersection;
- typedef typename Intersection::Geometry IntersectionGeometry;
-
- typedef typename GV::IntersectionIterator NeighborIterator;
- typedef HierarchicLevelIterator<GridType> HierarchicLevelIteratorType;
-
- typedef typename GridType::template Codim<0>::Entity Element;
-
- typedef typename Dune::MultipleCodimMultipleGeomTypeMapper<GV, FaceMapperLayout > FaceMapper;
-
- private:
- const int faultCount_;
- const int coarseResolution_;
- const int fineResolution_;
- const int exactResolution_;
-
- int coarseGridN_;
-
- std::shared_ptr<GridType> grid_;
-
- std::shared_ptr<LevelInterfaceNetwork<GV>> coarseInterfaceNetwork_;
- std::shared_ptr<LevelInterfaceNetwork<GV>> fineInterfaceNetwork_;
- std::shared_ptr<LevelInterfaceNetwork<GV>> exactInterfaceNetwork_;
-
-private:
-
-typename std::enable_if<!std::numeric_limits<ctype>::is_integer, bool>::type
- almost_equal(ctype x, ctype y, int ulp)
- {
- return std::abs(x-y) < std::numeric_limits<ctype>::epsilon() * std::abs(x+y) * ulp
- || std::abs(x-y) < std::numeric_limits<ctype>::min();
- }
-
-// returns fault id
-int isFaultIntersection(const Intersection& intersection, const std::set<int>& interfaceYCoords) {
- Dune::FieldVector<ctype, dimworld> faceCenter = intersection.geometry().center();
-
- ctype yID = faceCenter[1]*coarseGridN_;
- ctype intpart;
-
- if (almost_equal(std::modf(yID, &intpart), 0.0, 2)) {
- std::set<int>::iterator it = interfaceYCoords.find((int) intpart);
- if (it!=interfaceYCoords.end())
- return std::distance(interfaceYCoords.begin(), it);
- else
- return -1;
- } else
- return -1;
-}
-
-
- public:
- //setup
- UniformFaultFactory(const int faultCount, const int coarseResolution, const int fineResolution, const int exactResolution) :
- faultCount_(faultCount),
- coarseResolution_(coarseResolution),
- fineResolution_(fineResolution),
- exactResolution_(exactResolution),
- coarseGridN_(std::pow(2, coarseResolution_)),
- coarseInterfaceNetwork_(nullptr),
- fineInterfaceNetwork_(nullptr),
- exactInterfaceNetwork_(nullptr)
- {
- if (faultCount_>coarseGridN_-1)
- DUNE_THROW(Dune::Exception, "Coarse grid resolution too low to support " + std::to_string(faultCount_) + " fault interfaces!");
-
- Dune::UGGrid<dim>::setDefaultHeapSize(4000);
- GridOb unitCube(coarseGridN_);
- grid_ = unitCube.grid();
-
- const int fineLevelIdx = fineResolution_ - coarseResolution_;
- const int exactLevelIdx = exactResolution_ - coarseResolution_;
- grid_->globalRefine(exactLevelIdx);
-
- std::set<int> interfaceYCoords;
- int yStep = coarseGridN_/(faultCount+1);
- for (int i=1; i<=faultCount; i++) {
- interfaceYCoords.insert(i*yStep);
- }
-
- std::cout << "Coarse and fine grid were generated!" << std::endl;
-
- /*
- GV gridView(grid.levelGridView(0));
-
-
- VTKWriter<GV> vtkwriter(gridView);
- vtkwriter.pwrite("coarseGrid", "/home/mi/podlesjo/data/", "/home/mi/podlesjo/data/");//("~/home/mi/podlesjo/data/coarseGrid");
-
- GV gridViewFine(grid.levelGridView(fineLevelIdx));
-
-
- VTKWriter<GV> vtkwriterFine(gridViewFine);
- vtkwriterFine.pwrite("fineGrid", "/home/mi/podlesjo/data/", "/home/mi/podlesjo/data/");//("~/home/mi/podlesjo/data/coarseGrid");
- */
-
- typedef typename GV::template Codim <0>::Iterator ElementLevelIterator;
-
- GV coarseGridView(grid_->levelGridView(0));
- GV fineGridView(grid_->levelGridView(fineLevelIdx));
- GV exactGridView(grid_->levelGridView(exactLevelIdx));
-
- coarseInterfaceNetwork_ = std::make_shared<LevelInterfaceNetwork<GV>>(*grid_, 0);
- fineInterfaceNetwork_ = std::make_shared<LevelInterfaceNetwork<GV>>(*grid_, fineLevelIdx);
- exactInterfaceNetwork_ = std::make_shared<LevelInterfaceNetwork<GV>>(*grid_, exactLevelIdx);
-
- // init interfaces
- std::vector<std::shared_ptr<FaultInterface<GV>>> coarseInterfaces(faultCount);
- std::vector<std::shared_ptr<FaultInterface<GV>>> fineInterfaces(faultCount);
- std::vector<std::shared_ptr<FaultInterface<GV>>> exactInterfaces(faultCount);
- for (int i=0; i<faultCount; i++) {
- coarseInterfaces[i] = std::make_shared<FaultInterface<GV>>(FaultInterface<GV>(coarseGridView, 0));
- fineInterfaces[i] = std::make_shared<FaultInterface<GV>>(FaultInterface<GV>(fineGridView, fineLevelIdx));
- exactInterfaces[i] = std::make_shared<FaultInterface<GV>>(FaultInterface<GV>(exactGridView, exactLevelIdx));
- }
-
-
- /*
- VTKWriter<GV> vtkwriter(fineGridView);
- vtkwriter.pwrite("fineGrid", "/home/mi/podlesny/data/", "/home/mi/podlesny/data/");//("~/home/mi/podlesny/data/fineGrid");
-
- VTKWriter<GV> vtkwriterFine(exactGridView);
- vtkwriterFine.pwrite("exactGrid", "/home/mi/podlesny/data/", "/home/mi/podlesny/data/");//("~/home/mi/podlesny/data/exactGrid");
- */
-
- // TODO: remove after debugging
- /*Dune::LevelMultipleCodimMultipleGeomTypeMapper <GridType, Dune::MCMGElementLayout > fineMapper(*grid_, fineLevelIdx);
- ElementLevelIterator fineElemIt = fineGridView.template begin<0>();
- ElementLevelIterator fineElemEndIt = fineGridView.template end<0>();
- for (; fineElemIt!=fineElemEndIt; ++fineElemIt) {
- const Element& fineElem = *fineElemIt;
- const auto& geometry = fineElem.geometry();
- std::cout << "Element idx: " << fineMapper.index(fineElem) << std::endl;
-
- for(int i=0; i<geometry.corners(); ++i) {
- size_t globalIdx = fineGridView.indexSet().subIndex(fineElem, i, Element::dimension);
- const auto& vertex = geometry.corner(i);
- print(vertex, "elem corner:" + std::to_string(globalIdx));
- }
- std::cout << "----------------" << std::endl << std::endl;
- }*/
- // end remove
-
-
- FaceMapper coarseFaceMapper(coarseGridView);
- std::vector<bool> coarseHandled(coarseFaceMapper.size(),false);
-
- ElementLevelIterator coarseElemIt = coarseGridView.template begin<0>();
- ElementLevelIterator coarseElemEndIt = coarseGridView.template end<0>();
- for (; coarseElemIt!=coarseElemEndIt; ++coarseElemIt) {
- const Element& coarseElem = *coarseElemIt;
-
- bool hasFaultIntersection = false;
-
- int faceID, faultID = -1;
-
- NeighborIterator coarseNIt = coarseGridView.ibegin(coarseElem);
- NeighborIterator coarseNEnd = coarseGridView.iend(coarseElem);
- for (; coarseNIt != coarseNEnd; ++coarseNIt){
-
- if (coarseHandled[coarseFaceMapper.subIndex(coarseElem,coarseNIt->indexInInside(),1)])
- continue;
-
- coarseHandled[coarseFaceMapper.subIndex(coarseElem,coarseNIt->indexInInside(),1)]=true;
-
- faultID = isFaultIntersection(*coarseNIt, interfaceYCoords);
-
- if (faultID >= 0) {
- hasFaultIntersection = true;
- faceID = coarseNIt->indexInInside();
- break;
- }
- }
-
- if (hasFaultIntersection) {
- Face<GridType> face(*grid_, coarseElem, faceID);
-
- typename Face<GridType>::HierarchicIterator exactIt = face.hbegin(exactLevelIdx);
- typename Face<GridType>::HierarchicIterator exactEnd = face.hend(exactLevelIdx);
- for(; exactIt!=exactEnd; ++exactIt) {
- if (exactIt->level()==exactLevelIdx)
- exactInterfaces[faultID]->addFace(*(exactIt->intersection()));
- else if (exactIt->level()==fineLevelIdx)
- fineInterfaces[faultID]->addFace(*(exactIt->intersection()));
- else if (exactIt->level()==0)
- coarseInterfaces[faultID]->addFace(*(exactIt->intersection()));
- }
- }
- }
-
- for (int i=0; i<faultCount; i++) {
- coarseInterfaces[i]->build();
- coarseInterfaceNetwork_->addInterface(coarseInterfaces[i]);
-
- fineInterfaces[i]->build();
- fineInterfaceNetwork_->addInterface(fineInterfaces[i]);
- //print(fineInterfaceNetwork_->getInterfaceNetworkDofs(), "fineInterface dofs:");
-
- exactInterfaces[i]->build();
- exactInterfaceNetwork_->addInterface(exactInterfaces[i]);
- //print(exactInterfaceNetwork_->getInterfaceNetworkDofs(), "exactInterface dofs:");
- }
- }
-
- const GridType& grid() const {
- return *grid_;
- }
-
- const LevelInterfaceNetwork<GV>& coarseInterfaceNetwork() {
- return *coarseInterfaceNetwork_;
- }
-
-
- const LevelInterfaceNetwork<GV>& fineInterfaceNetwork() {
- return *fineInterfaceNetwork_;
- }
-
- const LevelInterfaceNetwork<GV>& exactInterfaceNetwork() {
- return *exactInterfaceNetwork_;
- }
-
-};
-#endif
diff --git a/dune/faultnetworks/faultinterface_old.hh b/dune/faultnetworks/faultinterface_old.hh
deleted file mode 100644
index 2120b2417a48675486aafcc55f67fe8edd04a2e0..0000000000000000000000000000000000000000
--- a/dune/faultnetworks/faultinterface_old.hh
+++ /dev/null
@@ -1,309 +0,0 @@
-#ifndef FAULT_INTERFACE_HH
-#define FAULT_INTERFACE_HH
-
-#include <queue>
-#include <map>
-
-#include <dune/common/bitsetvector.hh>
-#include <dune/common/fvector.hh>
-
-#include <dune/grid/common/mcmgmapper.hh>
-#include <dune/grid/common/intersectioniterator.hh>
-
-//template <typename GridView> class FaultInterface;
-
-//------------------------------------------------------
-//
-// old version, where inside and outside region seeds
-// are constructed for fault networks that connect
-// boundary to boundary and do not intersect
-//
-//------------------------------------------------------
-template <typename GridView>
-class FaultInterface : public std::enable_shared_from_this<FaultInterface<GridView>>
-{
-protected:
- static int idCounter;
-
- typedef typename GridView::Grid GridType;
- typedef typename GridType::LevelIntersection Intersection;
- typedef typename GridType::template Codim<0>::Entity Element;
- typedef typename Intersection::Geometry::GlobalCoordinate GlobalCoordinate;
- static const int dimElement = Element::dimension;
-
- typedef typename GridType::LevelIntersectionIterator LevelIntersectionIterator;
-
- typedef typename Dune::LevelMultipleCodimMultipleGeomTypeMapper <GridType, Dune::MCMGElementLayout > ElementMapper;
-
- const GridView& gridView_;
- const GridType& grid_;
- const int level_;
- const int id_;
- const std::shared_ptr<FaultInterface> father_;
-
- ElementMapper elementMapper_;
-
- std::vector<Intersection> faces_;
- std::map<std::pair<size_t, size_t>, bool> isInterfaceIntersectionMapper_;
-
- std::set<size_t> interfaceDofs_;
-
- bool requireBuild_;
-
- bool intersectionsIntersect(std::set<size_t>& intersection1Dofs, std::set<size_t>& intersection2Dofs) {
- std::vector<int> result(std::min<size_t>(intersection1Dofs.size(), intersection2Dofs.size()));
- std::vector<int>::iterator it;
-
- it = std::set_intersection (intersection1Dofs.begin(), intersection1Dofs.end(), intersection2Dofs.begin(), intersection2Dofs.end(), result.begin());
- result.resize(it-result.begin());
-
- return result.size()>0;
- }
-
- void computeIntersectionDofs(const Intersection& intersection, std::set<size_t>& intersectionDofs) {
- intersectionDofs.clear();
-
- // loop over all vertices of the intersection
- const Element& insideElement = intersection.inside();
-
- const auto& refElement = Dune::ReferenceElements<double,dimElement>::general(insideElement.type());
- for (int i=0; i<refElement.size(intersection.indexInInside(), 1, dimElement); i++) {
- size_t idxInElement = refElement.subEntity(intersection.indexInInside(), 1, i, dimElement);
- size_t globalIdx = gridView_.indexSet().subIndex(insideElement, idxInElement, dimElement);
- intersectionDofs.insert(globalIdx);
- }
- }
-
-public:
- FaultInterface(const GridView& gridView, const int level, std::shared_ptr<FaultInterface> father = nullptr) :
- gridView_(gridView), grid_(gridView_.grid()), level_(level), id_(idCounter++), father_(father), elementMapper_(grid_, level_), requireBuild_(false)
- {}
-
- bool isInterfaceIntersection(size_t elementIdx, size_t intersectionIdx) {
- return isInterfaceIntersectionMapper_.count(std::make_pair(elementIdx, intersectionIdx));
- }
-
- size_t getElementIndex(const Element& elem) {
- return elementMapper_.index(elem);
- }
-
- void addFace(const Intersection& face) {
- faces_.push_back(face);
-
- // add to interfaceDofs_
- std::set<size_t> faceDofs;
- computeIntersectionDofs(face, faceDofs);
-
- interfaceDofs_.insert(faceDofs.begin(), faceDofs.end());
-
- // add to isInterfaceIntersectionMapper_
- size_t insideElementIdx = elementMapper_.index(face.inside());
- isInterfaceIntersectionMapper_[std::make_pair(insideElementIdx, face.indexInInside())] = true;
-
- if (face.neighbor()) {
- size_t outsideElementIdx = elementMapper_.index(face.outside());
- isInterfaceIntersectionMapper_[std::make_pair(outsideElementIdx, face.indexInOutside())] = true;
- }
-
- requireBuild_ = true;
- }
-
- void build() {
- if (requireBuild_) {
- // define inside and outside regions by the respective tag of first intersection
- Intersection& intersection = faces_[0];
-
-
- const Element& insideElement = intersection.inside();
- size_t insideElementIdx = elementMapper_.index(insideElement);
- insideInterfaceRegion_.insert(insideElementIdx);
-
- std::set<size_t> visited;
- std::queue<Element> elementQueue;
-
- visited.insert(insideElementIdx);
- elementQueue.push(insideElement);
-
- while (!elementQueue.empty()) {
- const Element& elem = elementQueue.front();
- elementQueue.pop();
-
-
- LevelIntersectionIterator neighborIt = gridView_.ibegin(elem);
- LevelIntersectionIterator neighborItEnd = gridView_.iend(elem);
- for (; neighborIt!=neighborItEnd; ++neighborIt) {
- const Intersection& intersection = *neighborIt;
-
- if (!intersection.neighbor())
- continue;
-
- const Element& neighbor = intersection.outside();
- const size_t neighborIdx = elementMapper_.index(neighbor);
-
- if (visited.count(neighborIdx))
- continue;
-
- if (isInterfaceIntersection(elementMapper_.index(elem), intersection.indexInInside()))
- continue;
-
- std::set<size_t> intersectionDofs;
- computeIntersectionDofs(intersection, intersectionDofs);
-
- if (intersectionsIntersect(intersectionDofs, interfaceDofs_)) {
- elementQueue.push(neighbor);
- visited.insert(neighborIdx);
- insideInterfaceRegion_.insert(neighborIdx);
- }
-
- }
- }
-
- if (intersection.neighbor()) {
- const Element& outsideElement = intersection.outside();
- size_t outsideElementIdx = elementMapper_.index(outsideElement);
- outsideInterfaceRegion_.insert(outsideElementIdx);
- //std::cout << "outside region elem: " << outsideElementIdx << std::endl;
-
- visited.clear();
- visited.insert(outsideElementIdx);
- elementQueue.push(outsideElement);
-
- while (!elementQueue.empty()) {
- const Element& elem = elementQueue.front();
- elementQueue.pop();
-
- //std::cout << "intersection iteration" << std::endl;
-
- LevelIntersectionIterator neighborIt = gridView_.ibegin(elem);
- LevelIntersectionIterator neighborItEnd = gridView_.iend(elem);
- for (; neighborIt!=neighborItEnd; ++neighborIt) {
- const Intersection& intersection = *neighborIt;
-
- if (!intersection.neighbor())
- continue;
-
- const Element& neighbor = intersection.outside();
- const size_t neighborIdx = elementMapper_.index(neighbor);
-
- if (visited.count(neighborIdx))
- continue;
-
- if (isInterfaceIntersection(elementMapper_.index(elem), intersection.indexInInside()))
- continue;
-
- std::set<size_t> intersectionDofs;
- computeIntersectionDofs(intersection, intersectionDofs);
-
- // std::cout << "neighborIdx: " << neighborIdx << std::endl;
- // print(intersectionDofs, "intersectionDofs:");
-
- if (intersectionsIntersect(intersectionDofs, interfaceDofs_)) {
- // std::cout << "inserted" << std::endl;
- elementQueue.push(neighbor);
- visited.insert(neighborIdx);
- outsideInterfaceRegion_.insert(neighborIdx);
- }
-
- }
- }
- }
- }
-
- requireBuild_ = false;
- }
-
- double minDistance(const GlobalCoordinate refPoint) const {
- if (requireBuild_) {
- DUNE_THROW(Dune::Exception, "FaultInterface is not up to date, call build() before minDistance()!");
- } else {
- double result = std::numeric_limits<double>::infinity();
-
- for (size_t i=0; i<faces_.size(); i++) {
- // loop over all vertices of the intersection
- const Intersection& intersection = faces_[i];
- const typename Intersection::Geometry& geometry = intersection.geometry();
-
- for (int i=0; i<geometry.corners(); i++) {
- GlobalCoordinate candidate = refPoint;
- candidate -= geometry.corner(i);
-
- result = std::min(result, candidate.two_norm());
- }
- }
-
- return result;
- }
- }
-
- void averageUnitNormal(GlobalCoordinate normal) const {
- if (requireBuild_) {
- DUNE_THROW(Dune::Exception, "FaultInterface is not up to date, call build() before averageUnitNormal()!");
- } else {
- normal = 0;
- for (int i=0; i<faces_.size(); i++) {
- normal += faces_[i].centerUnitOuterNormal();
- }
- normal *= 1.0/faces_.size();
- }
- }
-
- bool inInside(const Element& elem) const {
- size_t elemID = elementMapper_.index(elem);
- return insideInterfaceRegion_.count(elemID);
- }
-
- size_t dofCount() const {
- return interfaceDofs_.size();
- }
-
- // local index of interface dof
- size_t localIndex(const size_t i) const {
- std::set<size_t>::iterator it;
-
- it = interfaceDofs_.find(i);
-
- if (it == interfaceDofs_.end())
- DUNE_THROW(Dune::Exception, "Global index of dof not part of the interface.");
-
- return std::distance(interfaceDofs_.begin(), it);
- }
-
- std::shared_ptr<const FaultInterface> coarsestFather() const {
- std::shared_ptr<const FaultInterface> res = this->shared_from_this();
-
- while (res->father() != nullptr)
- res = res->father();
-
- return res;
- }
-
- const std::set<size_t>& getInterfaceDofs() const {
- return interfaceDofs_;
- }
-
- const ElementMapper& getElementMapper() const {
- return elementMapper_;
- }
-
- const std::vector<Intersection>& faces() const {
- return faces_;
- }
-
- int level() const {
- return level_;
- }
-
- int id() const {
- return id_;
- }
-
- std::shared_ptr<FaultInterface> father() const {
- return father_;
- }
-
- const GridType& grid() const {
- return grid_;
- }
-};
-template <typename GridView> int FaultInterface<GridView>::idCounter = 0;
-#endif
diff --git a/dune/faultnetworks/faultp1nodalbasis_old.hh b/dune/faultnetworks/faultp1nodalbasis_old.hh
deleted file mode 100644
index fdb7225e7c11b455caf5b74ed835ed9a9fae9bc1..0000000000000000000000000000000000000000
--- a/dune/faultnetworks/faultp1nodalbasis_old.hh
+++ /dev/null
@@ -1,151 +0,0 @@
-#ifndef FAULT_P1_NODALBASIS_HH
-#define FAULT_P1_NODALBASIS_HH
-
-/**
- @file
- @brief
-
- @author
- */
-
-#include <dune/localfunctions/lagrange/pqkfactory.hh>
-#include <dune/fufem/functionspacebases/functionspacebasis.hh>
-#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
-#include <dune/istl/matrixindexset.hh>
-
-#include <dune/faultnetworks/levelinterfacenetwork.hh>
-
-template <class GV, class RT=double>
-class FaultP1NodalBasis :
- public P1NodalBasis<
- GV,
- RT>
-{
- protected:
- typedef typename GV::Grid::ctype ctype;
-
- typedef typename Dune::PQkLocalFiniteElementCache<typename GV::Grid::ctype, RT, GV::dimension, 1> FiniteElementCache;
- typedef typename FiniteElementCache::FiniteElementType LFE;
-
- typedef P1NodalBasis<GV, RT> Base;
- typedef typename Base::Element Element;
-
- using Base::dim;
- using Base::gridview_;
- using Base::cache_;
-
- typedef typename Dune::BCRSMatrix<Dune::FieldMatrix<RT, 1, 1>> MatrixType;
-
- MatrixType prolongationMatrix_;
- MatrixType restrictionMatrix_;
-
-
- const LevelInterfaceNetwork<GV>& faultNetwork_;
-
- public:
- typedef typename Base::GridView GridView;
- typedef typename Base::ReturnType ReturnType;
- typedef typename Base::LocalFiniteElement LocalFiniteElement;
- typedef typename Base::LinearCombination LinearCombination;
-
- FaultP1NodalBasis(const LevelInterfaceNetwork<GV>& faultNetwork) :
- Base(faultNetwork.levelGridView()),
- faultNetwork_(faultNetwork)
- {
- // compute prolongation and restriction matrices
- const std::set<size_t>& interfaceNetworkDofs = faultNetwork_.getInterfaceNetworkDofs();
-
- const size_t conDim = gridview_.indexSet().size(dim);
- const size_t disconDim = conDim + faultNetwork_.dofCount();
-
- Dune::MatrixIndexSet prolongIdxSet(disconDim, conDim);
- Dune::MatrixIndexSet restrictIdxSet(conDim, disconDim);
-
- // compute respective index sets
- for (size_t i=0; i<conDim; i++) {
- prolongIdxSet.add(i, i);
- restrictIdxSet.add(i, i);
- }
-
- std::set<size_t>::iterator beginIt = interfaceNetworkDofs.begin();
- for (size_t i=conDim; i<disconDim; i++) {
- std::set<size_t>::iterator it = beginIt;
- std::advance(it, i-conDim);
- const size_t idx = *it;
-
- prolongIdxSet.add(i, idx);
- restrictIdxSet.add(idx, i);
- }
-
- // set entries of prolongation
- prolongIdxSet.exportIdx(prolongationMatrix_);
- prolongationMatrix_ = 1;
-
- // set entries of restriction
- restrictIdxSet.exportIdx(restrictionMatrix_);
-
- typedef typename MatrixType::row_type RowType;
- typedef typename RowType::ConstIterator ColumnIterator;
-
- for(size_t rowIdx = 0; rowIdx<conDim; rowIdx++) {
- RowType& row = restrictionMatrix_[rowIdx];
-
- RT entry = 1.0/row.size();
-
- ColumnIterator colIt = row.begin();
- ColumnIterator colEndIt = row.end();
- for(; colIt!=colEndIt; ++colIt) {
- row[colIt.index()] = entry;
- }
- }
- }
-
- const LevelInterfaceNetwork<GV>& faultNetwork() const {
- return faultNetwork_;
- }
-
- size_t size() const
- {
- return gridview_.indexSet().size(dim) + faultNetwork_.dofCount();
- }
-
- const LocalFiniteElement& getLocalFiniteElement(const Element& e) const
- {
- return cache_.get(e.type());
- }
-
- int index(const Element& e, const int i) const
- {
- size_t globalIdx = gridview_.indexSet().subIndex(e, getLocalFiniteElement(e).localCoefficients().localKey(i).subEntity(), dim);
-
- if (faultNetwork_.inInside(e, globalIdx))
- return gridview_.indexSet().size(dim) + faultNetwork_.localIndex(globalIdx);
- else
- return globalIdx;
- }
-
- int indexInGridView(const Element& e, const int i) const
- {
- return gridview_.indexSet().subIndex(e, getLocalFiniteElement(e).localCoefficients().localKey(i).subEntity(), dim);
- }
-
- // prolong continuous p1 vector to discontinuous p1 vector
- template<class VectorType>
- void prolong(const VectorType& x, VectorType& res) {
- res.resize(prolongationMatrix_.N());
-
- prolongationMatrix_.mv(x, res);
- }
-
- // restrict discontinuous p1 vector to continuous p1 vector
- template<class VectorType>
- void restrict(const VectorType& x, VectorType& res) {
- res.resize(restrictionMatrix_.N());
-
- restrictionMatrix_.mv(x, res);
- }
-};
-
-#endif
-
-
diff --git a/dune/faultnetworks/interfacenetwork_old.hh b/dune/faultnetworks/interfacenetwork_old.hh
deleted file mode 100644
index f634453733ac2763ac26fc52e8661081b05b67d8..0000000000000000000000000000000000000000
--- a/dune/faultnetworks/interfacenetwork_old.hh
+++ /dev/null
@@ -1,264 +0,0 @@
-#ifndef INTERFACE_NETWORK_HH
-#define INTERFACE_NETWORK_HH
-
-
-#include <dune/common/bitsetvector.hh>
-#include <dune/common/fvector.hh>
-
-#include <dune/faultnetworks/facehierarchy.hh>
-#include <dune/faultnetworks/faultinterface.hh>
-#include <dune/faultnetworks/levelinterfacenetwork.hh>
-
-//------------------------------------------------------
-//
-// old version, where inside and outside region seeds
-// are constructed for fault networks that connect
-// boundary to boundary and do not intersect
-//
-//------------------------------------------------------
-template <class GridType>
-class InterfaceNetwork
-{
- protected:
- typedef typename GridType::LevelGridView GridView;
- typedef typename GridType::LevelIntersection Intersection;
- typedef typename GridType::ctype ctype;
- static const int dim = GridType::dimension;
- typedef typename GridType::template Codim<0>::Entity Element;
- typedef typename GridType::template Codim<dim>::Entity Vertex;
-
- private:
- const GridType& grid_;
- std::vector<std::shared_ptr<LevelInterfaceNetwork<GridView>>> levelInterfaceNetworks_;
-
- bool requireBuild_;
-
- public:
- InterfaceNetwork(const GridType& grid) : grid_(grid) {
-
- // init levelInterfaceNetworks
- levelInterfaceNetworks_.resize(grid.maxLevel()+1);
-
- for (size_t i=0; i<levelInterfaceNetworks_.size(); i++) {
- levelInterfaceNetworks_[i] = std::make_shared<LevelInterfaceNetwork<GridView>>(grid_, i);
- }
-
- requireBuild_ = false;
- }
-
- size_t size() const {
- return levelInterfaceNetworks_.size();
- }
-
- size_t size(int level) const {
- return levelInterfaceNetworks_[level]->size();
- }
-
- const std::shared_ptr<FaultInterface<GridView>> getInterface(const size_t interfaceID, const int level) const {
- return levelInterfaceNetworks_[level]->getInterface(interfaceID);
- }
-
- void build() {
- if (requireBuild_) {
- for (size_t i=0; i<levelInterfaceNetworks_.size(); i++) {
- levelInterfaceNetworks_[i]->build();
- }
-
- requireBuild_ = false;
- }
- }
-
- void prolongInterface(const std::shared_ptr<FaultInterface<GridView>> interfacePtr, LevelInterfaceNetwork<GridView>& toLevelInterfaceNetwork) {
- interfacePtr->build();
-
- const GridType& grid = interfacePtr->grid();
- const int toLevel = toLevelInterfaceNetwork.level();
-
- // init prolonged interfaces
- std::shared_ptr<FaultInterface<GridView>> newInterface = std::make_shared<FaultInterface<GridView>>(toLevelInterfaceNetwork.levelGridView(), toLevel, interfacePtr);
-
- // iterate over interface intersections
- const std::vector<Intersection>& faces = interfacePtr->faces();
- for (size_t i=0; i<faces.size(); i++) {
- const Intersection& intersection = faces[i];
- Face<GridType> face(grid, intersection.inside(), intersection.indexInInside());
-
- if (toLevel == interfacePtr->level()) {
- newInterface->addFace(intersection);
- } else {
- typename Face<GridType>::HierarchicIterator exactIt = face.hbegin(toLevel);
- typename Face<GridType>::HierarchicIterator exactEnd = face.hend(toLevel);
-
- for(; exactIt!=exactEnd; ++exactIt) {
- if (exactIt->level() == toLevel) {
- newInterface->addFace(*(exactIt->intersection()));
- }
- }
- }
- }
-
- // build new interface and add to levelInterfaceNetwork
- newInterface->build();
- toLevelInterfaceNetwork.addInterface(newInterface);
- toLevelInterfaceNetwork.build();
- }
-
- // CANNOT be used to copy interface to same level (prolong to level of interfacePtr->level())
- void prolongInterface(const std::shared_ptr<FaultInterface<GridView>> interfacePtr, const std::set<int> toLevels) {
- interfacePtr->build();
-
- const GridType& grid = interfacePtr->grid();
- const int maxLevel = *(toLevels.rbegin());
-
- // init prolonged interfaces
- std::vector<std::shared_ptr<FaultInterface<GridView>>> newInterfaces(toLevels.size());
- std::set<int>::iterator it = toLevels.begin();
- for (size_t i=0; i<newInterfaces.size(); i++) {
- const int newLevel = *it;
-
- // set father interface appropriately
- if (i==0)
- newInterfaces[i] = std::make_shared<FaultInterface<GridView>>(levelInterfaceNetworks_[newLevel]->levelGridView(), newLevel, interfacePtr);
- else
- newInterfaces[i] = std::make_shared<FaultInterface<GridView>>(levelInterfaceNetworks_[newLevel]->levelGridView(), newLevel, newInterfaces[i-1]);
-
- it++;
- }
-
- // iterate over interface intersections
- const std::vector<Intersection>& faces = interfacePtr->faces();
- for (size_t i=0; i<faces.size(); i++) {
- const Intersection& intersection = faces[i];
- Face<GridType> face(grid, intersection.inside(), intersection.indexInInside());
-
- typename Face<GridType>::HierarchicIterator exactIt = face.hbegin(maxLevel);
- typename Face<GridType>::HierarchicIterator exactEnd = face.hend(maxLevel);
- for(; exactIt!=exactEnd; ++exactIt) {
- std::set<int>::iterator it = toLevels.find(exactIt->level());
- if (it != toLevels.end()) {
- newInterfaces[std::distance(toLevels.begin(), it)]->addFace(*(exactIt->intersection()));
- }
- }
- }
-
- // build new interfaces and add to levelInterfaceNetworks
- for (size_t i=0; i<newInterfaces.size(); i++) {
- std::shared_ptr<FaultInterface<GridView>> newInterface = newInterfaces[i];
- newInterface->build();
-
- if (levelInterfaceNetworks_[newInterface->level()]->size()==0) {
- levelInterfaceNetworks_[newInterface->level()]->setReferencePoint(levelInterfaceNetworks_[interfacePtr->level()]->referencePoint());
- }
-
- levelInterfaceNetworks_[newInterface->level()]->addInterface(newInterface);
- levelInterfaceNetworks_[newInterface->level()]->build();
- }
- }
-
- void prolongLevelInterfaces(const int level, const int toLevel) {
- std::shared_ptr<LevelInterfaceNetwork<GridView>>& levelInterfaceNetworkPtr = levelInterfaceNetworks_[level];
- levelInterfaceNetworkPtr->build();
-
- for (size_t i=0; i<levelInterfaceNetworkPtr->size(); i++) {
- prolongInterface(levelInterfaceNetworkPtr->getInterface(i), *levelInterfaceNetworks_[toLevel]);
- }
- }
-
- void prolongLevelInterfaces(const int level, const std::set<int> toLevels) {
- std::shared_ptr<LevelInterfaceNetwork<GridView>>& levelInterfaceNetworkPtr = levelInterfaceNetworks_[level];
- levelInterfaceNetworkPtr->build();
-
- for (size_t i=0; i<levelInterfaceNetworkPtr->size(); i++) {
- prolongInterface(levelInterfaceNetworkPtr->getInterface(i), toLevels);
- }
- }
-
- void prolongInterface(const std::shared_ptr<FaultInterface<GridView>> interfacePtr, LevelInterfaceNetwork<GridView>& toLevelInterfaceNetwork) const {
- interfacePtr->build();
-
- const GridType& grid = interfacePtr->grid();
- const int toLevel = toLevelInterfaceNetwork.level();
-
- // init prolonged interfaces
- std::shared_ptr<FaultInterface<GridView>> newInterface = std::make_shared<FaultInterface<GridView>>(toLevelInterfaceNetwork.levelGridView(), toLevel, interfacePtr);
-
- // iterate over interface intersections
- const std::vector<Intersection>& faces = interfacePtr->faces();
- for (size_t i=0; i<faces.size(); i++) {
- const Intersection& intersection = faces[i];
- Face<GridType> face(grid, intersection.inside(), intersection.indexInInside());
-
- if (toLevel == interfacePtr->level()) {
- newInterface->addFace(intersection);
- } else {
- typename Face<GridType>::HierarchicIterator exactIt = face.hbegin(toLevel);
- typename Face<GridType>::HierarchicIterator exactEnd = face.hend(toLevel);
-
- for(; exactIt!=exactEnd; ++exactIt) {
- if (exactIt->level() == toLevel) {
- newInterface->addFace(*(exactIt->intersection()));
- }
- }
- }
- }
-
- // build new interface and add to levelInterfaceNetwork
- newInterface->build();
- toLevelInterfaceNetwork.addInterface(newInterface);
- toLevelInterfaceNetwork.build();
- }
-
- void prolongLevelInterfaces(const int level, LevelInterfaceNetwork<GridView>& toLevelInterfaceNetwork) const {
- const std::shared_ptr<LevelInterfaceNetwork<GridView>>& levelInterfaceNetworkPtr = levelInterfaceNetworks_[level];
- levelInterfaceNetworkPtr->build();
-
- for (size_t i=0; i<levelInterfaceNetworkPtr->size(); i++) {
- prolongInterface(levelInterfaceNetworkPtr->getInterface(i), toLevelInterfaceNetwork);
- }
- }
-
- void addInterface(std::shared_ptr<FaultInterface<GridView>>& interface){
- const int level = interface->level();
- levelInterfaceNetworks_[level]->addInterface(interface);
-
- requireBuild_ = true;
- }
-
- void addLevelInterfaceNetwork(std::shared_ptr<LevelInterfaceNetwork<GridView>> interfaceNetwork){
- const int level = interfaceNetwork->level();
- levelInterfaceNetworks_[level] = interfaceNetwork;
-
- requireBuild_ = true;
- }
-
-
- size_t localIndex(const size_t i, const int level) const {
- return levelInterfaceNetworks_[level]->localIndex(i);
- }
-
- size_t dofCount(const int level) const {
- return levelInterfaceNetworks_[level]->dofCount();
- }
-
- const std::set<size_t>& getInterfaceNetworkDofs(const int level) const{
- return levelInterfaceNetworks_[level]->getInterfaceNetworkDofs();
- }
-
- int maxLevel() const {
- return levelInterfaceNetworks_.size();
- }
-
- const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork(const int i) const {
- return *(levelInterfaceNetworks_[i]);
- }
-
- const GridType& grid() {
- return grid_;
- }
-
- const GridView& levelGridView(const int i) const {
- return levelInterfaceNetworks_[i]->levelGridView();
- }
-};
-
-#endif
diff --git a/dune/faultnetworks/levelinterfacenetwork_old.hh b/dune/faultnetworks/levelinterfacenetwork_old.hh
deleted file mode 100644
index 838d43a14fb310af30ebe5f5ae66017dc4a54a89..0000000000000000000000000000000000000000
--- a/dune/faultnetworks/levelinterfacenetwork_old.hh
+++ /dev/null
@@ -1,239 +0,0 @@
-#ifndef LEVEL_INTERFACE_NETWORK_HH
-#define LEVEL_INTERFACE_NETWORK_HH
-
-
-#include <dune/common/bitsetvector.hh>
-#include <dune/common/fvector.hh>
-
-#include <dune/faultnetworks/faultinterface.hh>
-
-//------------------------------------------------------
-//
-// old version, where inside and outside region seeds
-// are constructed for fault networks that connect
-// boundary to boundary and do not intersect
-//
-//------------------------------------------------------
-template <typename GridView>
-class LevelInterfaceNetwork
-{
- protected:
- // typedef typename BasisType::GridView GV;
- // typedef typename BasisType::LocalFiniteElement LFE;
-
- typedef typename GridView::Grid GridType;
- typedef typename GridType::ctype ctype;
- static const int dim = GridType::dimension;
- typedef typename GridType::template Codim<0>::Entity Element;
- typedef typename GridType::template Codim<dim>::Entity Vertex;
- typedef typename GridType::LevelIntersection::Geometry::GlobalCoordinate GlobalCoordinate;
-
- private:
- // const GridType& grid_;
- GridView gridView_;
- const int level_;
-
- const bool orderInterfaces_;
- GlobalCoordinate referencePoint_;
- std::map<int, double> interfaceDistances_;
-
- std::vector<std::shared_ptr<FaultInterface<GridView>>> interfaces_;
-
- std::map<size_t, std::shared_ptr<FaultInterface<GridView>>> dofToInterface_;
-
- std::set<size_t> interfaceNetworkDofs_;
-
- bool requireBuild_;
-
- class InterfaceDistanceCompare {
- private:
- std::map<int, double>& interfaceDistances;
-
- public:
- InterfaceDistanceCompare(std::map<int, double>& interfaceDistancesIn) : interfaceDistances(interfaceDistancesIn) {}
-
- bool operator()(const std::shared_ptr<FaultInterface<GridView>>& a, const std::shared_ptr<FaultInterface<GridView>>& b) const {
- return interfaceDistances[a->id()] < interfaceDistances[b->id()];
- }
- };
-
- public:
- LevelInterfaceNetwork(const GridType& grid, const int level, const GlobalCoordinate& refPoint) :
- gridView_(grid.levelGridView(level)), level_(level), orderInterfaces_(true), referencePoint_(refPoint) {
-
- interfaces_.resize(0);
- dofToInterface_.clear();
-
- requireBuild_ = false;
- }
-
- LevelInterfaceNetwork(const GridType& grid, const int level) :
- gridView_(grid.levelGridView(level)), level_(level), orderInterfaces_(true), referencePoint_(GlobalCoordinate()) {
-
- interfaces_.resize(0);
- dofToInterface_.clear();
-
- requireBuild_ = false;
- }
-
- size_t size() const {
- return interfaces_.size();
- }
-
- bool inInside(const Element& elem, const size_t globalIdx) const {
- typename std::map<size_t, std::shared_ptr<FaultInterface<GridView>>>::const_iterator it = dofToInterface_.find(globalIdx);
-
- if (it == dofToInterface_.end())
- return false;
- else
- return (it->second)->inInside(elem);
- }
-
- int isInterfaceIntersection(const Element& elem, const size_t intersectionIdx) const {
- int res = -1;
-
- for (size_t i=0; i<interfaces_.size(); ++i) {
- std::shared_ptr<FaultInterface<GridView>> interfacePtr = interfaces_[i];
- const size_t elemIdx = interfacePtr->getElementIndex(elem);
-
- if (interfacePtr->isInterfaceIntersection(elemIdx, intersectionIdx)) {
- res = i;
- break;
- }
- }
-
- return res;
- }
-
-
- int isInterfaceIntersection(const Element& elem, const size_t intersectionIdx, const int level) {
- int res = -1;
-
- for (size_t i=0; i<interfaces_.size(); ++i) {
- std::shared_ptr<FaultInterface<GridView>> interfacePtr = interfaces_[i];
-
- if (interfacePtr->getLevel()==level) {
- const size_t elemIdx = interfacePtr->getElementIndex(elem);
-
- if (interfacePtr->isInterfaceIntersection(elemIdx, intersectionIdx)) {
- res = i;
- break;
- }
- }
- }
-
- return res;
- }
-
- const std::shared_ptr<FaultInterface<GridView>> getInterface(const size_t interfaceID) const {
- return interfaces_[interfaceID];
- }
-
- void build() {
- if (requireBuild_) {
- if (orderInterfaces_) {
- std::sort(interfaces_.begin(), interfaces_.end(), InterfaceDistanceCompare(interfaceDistances_));
- }
-
- for (size_t i=0; i<interfaces_.size(); i++) {
- interfaces_[i]->build();
- }
-
- requireBuild_ = false;
- }
- }
-
- void addInterface(std::shared_ptr<FaultInterface<GridView>> interface){
-
- // for the moment: no (self) intersections allowed
- bool admissible = true;
- int errorCode = 0;
-
- // TODO: check for self intersection
- //std::set<size_t> interfaceDofs;
- //interface->getInterfaceDofs(interfaceDofs);
- const std::set<size_t>& interfaceDofs = interface->getInterfaceDofs();
-
- std::set<size_t>::iterator it = interfaceDofs.begin();
- std::set<size_t>::iterator endIt = interfaceDofs.end();
- for (; it != endIt; ++it) {
- size_t newDof = *it;
-
- // checking for intersection with other (already added) faults
- if (dofToInterface_.find(newDof) != dofToInterface_.end()) {
- admissible = false;
- errorCode = 1;
- break;
- }
- }
-
- // checking for self intersection
- if (false) {
- admissible = false;
- errorCode = 2;
- }
-
- size_t interfaceID = interfaces_.size();
- //if (admissible) {
- if (true) {
- for (it = interfaceDofs.begin(); it != endIt; ++it) {
- // build dofToInterfaceID_
- dofToInterface_[*it] = interface;
- interfaceNetworkDofs_.insert(*it);
- }
-
- interfaces_.push_back(interface);
-
- if (orderInterfaces_) {
- interfaceDistances_[interface->id()] = interface->minDistance(referencePoint_);
- }
-
- requireBuild_ = true;
- } else {
- DUNE_THROW(Dune::Exception, "Adding fault with ID " << interfaceID << " failed! " << ((errorCode==1) ? "Intersection with other fault!" : "Self intersection!"));
- }
- }
-
-
- // running local index of interface dof
- size_t localIndex(const size_t i) const {
- std::set<size_t>::iterator it;
-
- it = interfaceNetworkDofs_.find(i);
-
- if (it == interfaceNetworkDofs_.end())
- DUNE_THROW(Dune::Exception, "Global index of dof not part of the interface network.");
-
- return std::distance(interfaceNetworkDofs_.begin(), it);
- }
-
- size_t dofCount() const {
- return interfaceNetworkDofs_.size();
- }
-
- const std::set<size_t>& getInterfaceNetworkDofs() const{
- return interfaceNetworkDofs_;
- }
-
- int level() const {
- return level_;
- }
-
- const GridView& levelGridView() const {
- return gridView_;
- }
-
- const GridType& grid() const {
- return gridView_.grid();
- }
-
- const GlobalCoordinate& referencePoint() const {
- return referencePoint_;
- }
-
- void setReferencePoint(const GlobalCoordinate& newRefPoint) {
- referencePoint_ = newRefPoint;
- }
-};
-
-#endif
diff --git a/dune/faultnetworks/patchfactories/CMakeLists.txt b/dune/faultnetworks/patchfactories/CMakeLists.txt
index 8dd976ff4ba76769cc95fb3e31be29f22f07c1d3..146bebeca7469f31bbb2375df1e22822ba8fd2e7 100644
--- a/dune/faultnetworks/patchfactories/CMakeLists.txt
+++ b/dune/faultnetworks/patchfactories/CMakeLists.txt
@@ -1,11 +1,9 @@
#add_subdirectory("test" EXCLUDE_FROM_ALL)
add_custom_target(faultnetworks_patchfactories_src SOURCES
- localpatchfactory.hh
supportpatchfactory.hh
)
install(FILES
- localpatchfactory.hh
supportpatchfactory.hh
-DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/faultnetworks/patchfactories)
\ No newline at end of file
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/faultnetworks/patchfactories)
diff --git a/dune/faultnetworks/patchfactories/localpatchfactory.hh b/dune/faultnetworks/patchfactories/localpatchfactory.hh
deleted file mode 100644
index 380ea5840fe24a5f312f9fac654898e5b0ae1390..0000000000000000000000000000000000000000
--- a/dune/faultnetworks/patchfactories/localpatchfactory.hh
+++ /dev/null
@@ -1,299 +0,0 @@
-#ifndef LOCAL_PATCH_FACTORY_HH
-#define LOCAL_PATCH_FACTORY_HH
-
-#include<queue>
-
-#include <dune/common/bitsetvector.hh>
-#include <dune/common/fvector.hh>
-
-//#include <dune/fufem/functionspacebases/dofconstraints.hh>
-#include <dune/fufem/referenceelementhelper.hh>
-
-#include <dune/grid/common/mcmgmapper.hh>
-#include <dune/faultnetworks/hierarchicleveliterator.hh>
-#include <dune/faultnetworks/levelinterfacenetwork.hh>
-
-
-/* --------------------------------------------------------------------------------
- *
- * Deprecated: LevelInterfaceNetwork does not support individual interfaces anymore
- *
- * --------------------------------------------------------------------------------
- */
-
-template <class BasisType>
-class LocalPatchFactory
-{
- protected:
- typedef typename BasisType::GridView GV;
- typedef typename BasisType::LocalFiniteElement LFE;
-
- typedef typename GV::Grid GridType;
- typedef typename GridType::LevelIntersection Intersection;
- typedef typename GridType::LevelIntersectionIterator LevelIntersectionIterator;
-
- typedef typename GridType::ctype ctype;
- static const int dim = GridType::dimension;
- typedef typename GridType::template Codim<0>::Entity Element;
- typedef typename GridType::template Codim<dim>::Entity Vertex;
- static const int dimElement = Element::dimension;
-
- typedef typename Dune::LevelMultipleCodimMultipleGeomTypeMapper<GridType, Dune::MCMGElementLayout > ElementMapper;
-
- private:
-
- const BasisType& basis_;
-
- const int level_;
-
- const LevelInterfaceNetwork<GV>& interfaceNetwork_;
- const Dune::BitSetVector<1> activeFaults_;
- const int interfaceID_;
-
- const ElementMapper& elementMapper_;
-
- const GV& gridView_;
- const GridType& grid_;
-
-
- std::vector<int> localToGlobal_;
- Dune::BitSetVector<1> boundaryDofs_;
- std::vector<Element> patchElements_;
-
- public:
- //setup
- LocalPatchFactory(const BasisType& basis,
- const LevelInterfaceNetwork<GV>& interfaceNetwork,
- const int interfaceID,
- const Dune::BitSetVector<1> activeFaults) :
- basis_(basis),
- level_(interfaceNetwork.getInterface(interfaceID)->level()),
- interfaceNetwork_(interfaceNetwork),
- activeFaults_(activeFaults),
- interfaceID_(interfaceID),
- elementMapper_(interfaceNetwork_.getInterface(interfaceID_)->getElementMapper()),
- gridView_(basis_.getGridView()),
- grid_(gridView_.grid())
- {
-
- assert(interfaceNetwork_.size() == activeFaults_.size());
-
- patchElements_.resize(0);
-
- std::set<int> localDofs;
- std::set<int> localBoundaryDofs;
-
- std::shared_ptr<FaultInterface<GV>> interface = interfaceNetwork_.getInterface(interfaceID_);
- const std::vector<Intersection>& faces = interface->faces();
- const Element& seedElement = faces[0].inside();
-
- std::set<size_t> visited;
- std::queue<Element> elementQueue;
-
- visited.insert(elementMapper_.index(seedElement));
- elementQueue.push(seedElement);
- patchElements_.push_back(seedElement);
-
- while (!elementQueue.empty()) {
- const Element& elem = elementQueue.front();
- elementQueue.pop();
-
- // insert element dofs into localDofs
- const LFE& lfe = basis_.getLocalFiniteElement(elem);
-
- for (size_t i=0; i<lfe.localBasis().size(); ++i) {
- int dofIndex = basis_.index(elem, i);
- localDofs.insert(dofIndex);
- }
-
- LevelIntersectionIterator neighborIt = gridView_.ibegin(elem);
- LevelIntersectionIterator neighborItEnd = gridView_.iend(elem);
- for (; neighborIt!=neighborItEnd; ++neighborIt) {
- const Intersection& intersection = *neighborIt;
-
- bool isBoundaryIntersection = false;
-
-
- if (intersection.boundary()) {
- isBoundaryIntersection = true;
- } else {
-
- int isInterfaceBoundary = interfaceNetwork_.isInterfaceIntersection(elem, intersection.indexInInside());
-
- if (isInterfaceBoundary != -1 && isInterfaceBoundary != interfaceID_ && activeFaults_[isInterfaceBoundary][0]) {
- isBoundaryIntersection = true;
- } else {
- const Element& neighbor = intersection.outside();
- const size_t neighborIdx = elementMapper_.index(neighbor);
-
- if (!visited.count(neighborIdx)) {
- elementQueue.push(neighbor);
- visited.insert(neighborIdx);
- patchElements_.push_back(elem);
- }
- }
- }
-
- if (isBoundaryIntersection) {
- const auto& refElement = Dune::ReferenceElements<double,dimElement>::general(elem.type());
- for (int i=0; i<refElement.size(intersection.indexInInside(), 1, dimElement); i++) {
- size_t idxInElement = refElement.subEntity(intersection.indexInInside(), 1, i, dimElement);
- localBoundaryDofs.insert(basis_.index(elem, idxInElement));
- }
-
- /* typedef typename LFE::Traits::LocalCoefficientsType LocalCoefficients;
- const LocalCoefficients* localCoefficients = &basis_.getLocalFiniteElement(elem).localCoefficients();
-
- for (size_t i=0; i<localCoefficients->size(); i++) {
- unsigned int entity = localCoefficients->localKey(i).subEntity();
- unsigned int codim = localCoefficients->localKey(i).codim();
-
- if (containsInsideSubentity(elem, intersection, entity, codim))
- localBoundaryDofs.insert(basis_.index(elem, i));
- }*/
- }
- }
- }
-
- localToGlobal_.resize(localDofs.size());
- boundaryDofs_.resize(localDofs.size());
- boundaryDofs_.unsetAll();
-
- std::set<int>::iterator dofIt = localDofs.begin();
- std::set<int>::iterator dofEndIt = localDofs.end();
- size_t i=0;
- for (; dofIt != dofEndIt; ++dofIt) {
- int localDof = *dofIt;
- localToGlobal_[i] = localDof;
-
- if (localBoundaryDofs.count(localDof)) {
- boundaryDofs_[i][0] = true;
- }
- i++;
- }
- }
-
- LocalPatchFactory(const BasisType& basis,
- const LevelInterfaceNetwork<GV>& interfaceNetwork,
- const int interfaceID) : LocalPatchFactory(basis, interfaceNetwork, interfaceID, Dune::BitSetVector<1>(interfaceNetwork.size(), true)){}
-
- LocalPatchFactory(const BasisType& basis,
- const LevelInterfaceNetwork<GV>& interfaceNetwork,
- const Element& seedElement,
- const Dune::BitSetVector<1> activeFaults) :
- basis_(basis),
- level_(interfaceNetwork.level()),
- interfaceNetwork_(interfaceNetwork),
- activeFaults_(activeFaults),
- interfaceID_(-1),
- elementMapper_(interfaceNetwork_.getInterface(0)->getElementMapper()),
- gridView_(basis_.getGridView()),
- grid_(gridView_.grid())
- {
-
- patchElements_.resize(0);
-
- std::set<int> localDofs;
- std::set<int> localBoundaryDofs;
-
- std::set<size_t> visited;
- std::queue<Element> elementQueue;
-
- visited.insert(elementMapper_.index(seedElement));
- elementQueue.push(seedElement);
- patchElements_.push_back(seedElement);
-
- while (!elementQueue.empty()) {
- const Element& elem = elementQueue.front();
- elementQueue.pop();
-
- // insert element dofs into localDofs
- const LFE& lfe = basis_.getLocalFiniteElement(elem);
-
- for (size_t i=0; i<lfe.localBasis().size(); ++i) {
- int dofIndex = basis_.index(elem, i);
- localDofs.insert(dofIndex);
- }
-
- LevelIntersectionIterator neighborIt = gridView_.ibegin(elem);
- LevelIntersectionIterator neighborItEnd = gridView_.iend(elem);
- for (; neighborIt!=neighborItEnd; ++neighborIt) {
- const Intersection& intersection = *neighborIt;
-
- bool isBoundaryIntersection = false;
-
- if (intersection.boundary()) {
- isBoundaryIntersection = true;
- } else {
- int isInterfaceBoundary = interfaceNetwork_.isInterfaceIntersection(elem, intersection.indexInInside());
-
- if (isInterfaceBoundary != -1 && activeFaults_[isInterfaceBoundary][0]) {
- isBoundaryIntersection = true;
- } else {
- const Element& neighbor = intersection.outside();
- const size_t neighborIdx = elementMapper_.index(neighbor);
-
- if (!visited.count(neighborIdx)) {
- elementQueue.push(neighbor);
- visited.insert(neighborIdx);
- patchElements_.push_back(elem);
- }
- }
- }
-
- if (isBoundaryIntersection) {
- const auto& refElement = Dune::ReferenceElements<double,dimElement>::general(elem.type());
- for (int i=0; i<refElement.size(intersection.indexInInside(), 1, dimElement); i++) {
- size_t idxInElement = refElement.subEntity(intersection.indexInInside(), 1, i, dimElement);
- localBoundaryDofs.insert(basis_.index(elem, idxInElement));
- }
-
- /* typedef typename LFE::Traits::LocalCoefficientsType LocalCoefficients;
- const LocalCoefficients* localCoefficients = &basis_.getLocalFiniteElement(elem).localCoefficients();
-
- for (size_t i=0; i<localCoefficients->size(); i++) {
- unsigned int entity = localCoefficients->localKey(i).subEntity();
- unsigned int codim = localCoefficients->localKey(i).codim();
-
- if (containsInsideSubentity(elem, intersection, entity, codim))
- localBoundaryDofs.insert(basis_.index(elem, i));
- }*/
- }
- }
- }
-
- localToGlobal_.resize(localDofs.size());
- boundaryDofs_.resize(localDofs.size());
- boundaryDofs_.unsetAll();
-
- std::set<int>::iterator dofIt = localDofs.begin();
- std::set<int>::iterator dofEndIt = localDofs.end();
- size_t i=0;
- for (; dofIt != dofEndIt; ++dofIt) {
- int localDof = *dofIt;
- localToGlobal_[i] = localDof;
-
- if (localBoundaryDofs.count(localDof)) {
- boundaryDofs_[i][0] = true;
- }
- i++;
- }
- }
-
- LocalPatchFactory(const BasisType& basis,
- const LevelInterfaceNetwork<GV>& interfaceNetwork,
- const Element& seedElement) : LocalPatchFactory(basis, interfaceNetwork, seedElement, Dune::BitSetVector<1>(interfaceNetwork.size(), true)) {}
-
- std::vector<int>& getLocalToGlobal() {
- return localToGlobal_;
- }
-
- std::vector<Element>& getPatchElements() {
- return patchElements_;
- }
-
- Dune::BitSetVector<1>& getBoundaryDofs() {
- return boundaryDofs_;
- }
-};
-#endif
diff --git a/dune/faultnetworks/preconditioners/CMakeLists.txt b/dune/faultnetworks/preconditioners/CMakeLists.txt
index 497bb345dbd166fee2658d54211fbbac9ccb6ef0..f4274985a25e585c4acf41c0a6506fd6ce0256a5 100644
--- a/dune/faultnetworks/preconditioners/CMakeLists.txt
+++ b/dune/faultnetworks/preconditioners/CMakeLists.txt
@@ -1,21 +1,13 @@
-add_subdirectory("test" EXCLUDE_FROM_ALL)
+
add_custom_target(faultnetworks_preconditioners_src SOURCES
- levelfaultpreconditioner.hh
levelglobalpreconditioner.hh
levelpatchpreconditioner.hh
- multilevelfaultpreconditioner.hh
multilevelpatchpreconditioner.hh
- faultpreconditioner.hh
- faultpreconditioner2.hh
)
install(FILES
- levelfaultpreconditioner.hh
levelglobalpreconditioner.hh
levelpatchpreconditioner.hh
- multilevelfaultpreconditioner.hh
multilevelpatchpreconditioner.hh
- faultpreconditioner.hh
- faultpreconditioner2.hh
-DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/faultnetworks/preconditioners)
\ No newline at end of file
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/faultnetworks/preconditioners)
diff --git a/dune/faultnetworks/preconditioners/faultpreconditioner.hh b/dune/faultnetworks/preconditioners/faultpreconditioner.hh
deleted file mode 100644
index 9b045b7918aa8418ae99de0d60b28270aeadd529..0000000000000000000000000000000000000000
--- a/dune/faultnetworks/preconditioners/faultpreconditioner.hh
+++ /dev/null
@@ -1,331 +0,0 @@
-#ifndef FAULT_PRECONDITIONER_HH
-#define FAULT_PRECONDITIONER_HH
-
-#include <queue>
-#include <string>
-
-#include <dune/common/timer.hh>
-#include <dune/common/fvector.hh>
-
-#include <dune/solvers/iterationsteps/lineariterationstep.hh>
-#include <dune/faultnetworks/compressedmultigridtransfer.hh>
-
-#include <dune/faultnetworks/patchfactories/localpatchfactory.hh>
-//#include <dune/osccoeff/osclocalproblemeff_latestWorking.hh>
-#include <dune/faultnetworks/localproblem.hh>
-#include <dune/faultnetworks/assemblers/osclocalassembler.hh>
-#include <dune/faultnetworks/assemblers/oscl2functionalassembler.hh>
-#include <dune/faultnetworks/hierarchicleveliterator.hh>
-#include <dune/faultnetworks/oscrhs.hh>
-#include <dune/faultnetworks/utils/debugutils.hh>
-
-#include <dune/fufem/assemblers/assembler.hh>
-#include <dune/fufem/boundarypatch.hh>
-#include <dune/fufem/functiontools/boundarydofs.hh>
-
-#include <dune/grid/common/mcmgmapper.hh>
-#include <dune/grid/common/intersectioniterator.hh>
-
-#include <dune/istl/superlu.hh>
-
-#include <dune/faultnetworks/levelinterfacenetwork.hh>
-
-/* --------------------------------------------------------------------------------
- *
- * Deprecated: LevelInterfaceNetwork does not support individual interfaces anymore
- *
- * --------------------------------------------------------------------------------
- */
-
-template <class GridType, class CoarseBasisType, class DGBasisType, class MatrixType, class VectorType>
-class FaultPreconditioner : public LinearIterationStep<MatrixType, VectorType> {
-
-public:
- typedef GridType grid_type;
- typedef MatrixType matrix_type;
-
-private:
- static const int dim= GridType::dimension;
- static const int dimworld = GridType::dimensionworld;
- typedef typename GridType::ctype ctype;
- typedef typename Dune::FieldVector<ctype, dimworld> WorldVectorType;
- typedef typename Dune::LevelMultipleCodimMultipleGeomTypeMapper <GridType, Dune::MCMGElementLayout > MapperType;
- typedef typename GridType::template Codim<0>::Entity EntityType;
- typedef HierarchicLevelIterator<GridType> HierarchicLevelIteratorType;
-
-
- const GridType& grid_;
- std::vector<ctype>& oscData_;
-
- const int coarseLevel_;
- const int fineLevel_;
-
- MatrixType coarseStiffMat_;
- VectorType coarseRhs_;
-
- MapperType& mapper_;
-
- MatrixType& fineMat_;
- VectorType& fineRhs_;
-
- const LevelInterfaceNetwork<typename GridType::LevelGridView>& interfaceNetwork_;
-
- const bool debugMode_;
-
- std::vector<OscLocalProblem<MatrixType, VectorType>* > localProblems_;
-
- DGBasisType fineBasis_;
- Dune::BitSetVector<1> coarseBoundaryDofs_;
- CompressedMultigridTransfer<VectorType> mgTransfer_;
-
-
-public:
-
- FaultPreconditioner(const GridType& grid, int coarseLevel, int fineLevel, std::vector<ctype>& oscData, MapperType& mapper, MatrixType& fineMat, VectorType& fineRhs, const LevelInterfaceNetwork<typename GridType::LevelGridView>& interfaceNetwork, const bool debugMode)
- : grid_(grid), oscData_(oscData), coarseLevel_(coarseLevel), fineLevel_(fineLevel), mapper_(mapper), fineMat_(fineMat), fineRhs_(fineRhs), interfaceNetwork_(interfaceNetwork), debugMode_(debugMode), fineBasis_(grid_.levelGridView(fineLevel_), interfaceNetwork) {
-
- typedef typename GridType::LevelGridView GV;
- typedef typename CoarseBasisType::LocalFiniteElement CoarseFE;
-
- GV coarseGridView(grid_.levelGridView(coarseLevel_));
- GV fineGridView(grid_.levelGridView(fineLevel_));
-
- CoarseBasisType coarseBasis(coarseGridView);
-
- std::cout << "Setting coarse grid oscData:" << std::endl;
- Dune::LevelMultipleCodimMultipleGeomTypeMapper <GridType, Dune::MCMGElementLayout > coarseMapper(grid_, coarseLevel_);
-
- // allocate a vector for the data
- std::vector<ctype> coarseOscData(coarseMapper.size());
-
-
-
- typedef typename GV::template Codim <0>::Iterator ElementLevelIterator;
- ElementLevelIterator endElemIt = coarseGridView.template end <0>();
- for (ElementLevelIterator it = coarseGridView. template begin <0>(); it!=endElemIt; ++it) {
-
- // set coarseOscData
- int counter = 0;
- ctype sum = 0;
-
- HierarchicLevelIteratorType endHierIt(grid_, *it, HierarchicLevelIteratorType::PositionFlag::end, fineLevel_);
- for (HierarchicLevelIteratorType hierIt(grid_, *it, HierarchicLevelIteratorType::PositionFlag::begin, fineLevel_); hierIt!=endHierIt; ++hierIt) {
-
- sum += oscData_[mapper_.index(*hierIt)];
- counter++;
- }
- coarseOscData[coarseMapper.index(*it)] = sum/counter;
- }
-
-
- std::cout << "Assembling coarse grid stiffness matrix!" << std::endl;
-
- // assemble coarse stiffness matrix
- Assembler<CoarseBasisType, CoarseBasisType> coarseAssembler(coarseBasis, coarseBasis);
- OscLocalAssembler<GridType, CoarseFE, CoarseFE> coarseLocalAssembler(coarseOscData, coarseMapper);
- coarseAssembler.assembleOperator(coarseLocalAssembler, coarseStiffMat_);
-
-
- std::cout << "Setting coarse grid boundary conditions!" << std::endl;
-
- typedef typename MatrixType::row_type RowType;
- typedef typename RowType::ConstIterator ColumnIterator;
-
- // set coarse boundary conditions
- BoundaryPatch<typename GridType::LevelGridView> coarseBoundaryPatch(coarseGridView, true);
- constructBoundaryDofs(coarseBoundaryPatch, coarseBasis, coarseBoundaryDofs_);
-
- for(size_t i=0; i<coarseBoundaryDofs_.size(); i++) {
- if(!coarseBoundaryDofs_[i][0])
- continue;
-
- const RowType& row = coarseStiffMat_[i];
-
- ColumnIterator cIt = row.begin();
- ColumnIterator cEndIt = row.end();
-
- for(; cIt!=cEndIt; ++cIt) {
- coarseStiffMat_[i][cIt.index()] = 0;
- }
- coarseStiffMat_[i][i]=1;
- }
-
- mgTransfer_.setup(coarseGridView, fineGridView);
-
- localProblems_.resize(interfaceNetwork_.size()+2);
-
- std::cout << "Initializing local fine grid corrections! " << std::endl;
-
- typedef typename GridType::template Codim<0>::Entity Element;
- typedef typename GridType::LevelIntersection Intersection;
- std::vector<Element> faultZoneSeeds;
- std::vector<int> faultZoneSeedsDofs;
- std::vector<int> faultZoneSeedsCounter;
- for (size_t i=0; i<interfaceNetwork_.size(); i++) {
- const std::vector<Intersection>& faces = interfaceNetwork_.getIntersections();
-
- const Element& insideElem = faces[0].inside();
- faultZoneSeeds.push_back(insideElem);
- faultZoneSeedsDofs.push_back(fineBasis_.index(insideElem, 0));
-
- const Element& outsideElem = faces[0].outside();
- faultZoneSeeds.push_back(outsideElem);
- faultZoneSeedsDofs.push_back(fineBasis_.index(outsideElem, 0));
- }
- faultZoneSeedsCounter.resize(faultZoneSeeds.size(), 0);
-
- // init local fine level corrections
- Timer timer;
- timer.reset();
- timer.start();
-
-
- for (size_t i=0; i<interfaceNetwork_.size(); i++) {
- // std::cout << i << std::endl;
- // std::cout << "---------------" << std::endl;
- LocalPatchFactory<DGBasisType> patchFactory(fineBasis_, interfaceNetwork_, i);
-
- std::vector<int>& localToGlobal = patchFactory.getLocalToGlobal();
- Dune::BitSetVector<1>& boundaryDofs = patchFactory.getBoundaryDofs();
-
- for (size_t j=0; j<faultZoneSeeds.size(); j++) {
- if (std::find(localToGlobal.begin(), localToGlobal.end(), faultZoneSeedsDofs[j]) != localToGlobal.end()) {
- faultZoneSeedsCounter[j]++;
- }
- }
-
-
- // print(localToGlobal, "localToGlobal: ");
- // print(boundaryDofs, "boundaryDofs: ");
-
-/* std::vector<EntityType>& regionElements = patchFactory.getRegionElements();
-
-
- Dune::BitSetVector<1> fineBoundaryDofs(fineRhs_.size());
- fineBoundaryDofs.unsetAll();
- for (size_t j=0; j<localToGlobal.size(); j++) {
- if (boundaryDofs[j][0]) {
- fineBoundaryDofs[localToGlobal[j]][0] = true;
- }
- }
-
- OscRhs<WorldVectorType, Dune::FieldVector<ctype,1>> f;
- OscLocalAssembler<GridType, FE, FE> oscLocalAssembler(oscData_, mapper_, coarseLocalProblemRegion);
- OscL2FunctionalAssembler<GridType, FE> oscLocalFunctionalAssembler(f, coarseLocalProblemRegion);
- OscLocalProblem<BasisType, MatrixType, VectorType, VectorType>* localProblem = new OscLocalProblem<BasisType, MatrixType, VectorType, VectorType>(regionElements, fineBasis, fineBoundaryDofs, oscLocalAssembler, fineRhs_, oscLocalFunctionalAssembler);
-*/
-
-
- VectorType rhs;
- rhs.resize(fineMat_.N());
- rhs = 0;
-
- localProblems_[i] = new OscLocalProblem<MatrixType, VectorType>(fineMat_, rhs, localToGlobal, boundaryDofs, false);
-
- //print(localToGlobal, "localToGlobal: ");
- //print(boundaryDofs, "boundaryDofs: ");
-
- if ((i+1) % 10 == 0) {
- std::cout << '\r' << std::floor((i+1.0)/coarseBasis.size()*100) << " % done. Elapsed time: " << timer.elapsed() << " seconds. Predicted total setup time: " << timer.elapsed()/(i+1.0)*coarseBasis.size() << " seconds." << std::flush;
- }
- }
-
- std::vector<Element> seedElements;
- for (size_t j=0; j<faultZoneSeedsCounter.size(); j++) {
- if (faultZoneSeedsCounter[j]==1)
- seedElements.push_back(faultZoneSeeds[j]);
- }
-
- assert(seedElements.size()==2);
-
- VectorType rhs1, rhs2;
- rhs1.resize(fineMat_.N());
- rhs2.resize(fineMat_.N());
- rhs1 = 0;
- rhs2 = 0;
-
- LocalPatchFactory<DGBasisType> patchFactory1(fineBasis_, interfaceNetwork_, seedElements[0]);
- std::vector<int>& localToGlobal1 = patchFactory1.getLocalToGlobal();
- Dune::BitSetVector<1>& boundaryDofs1 = patchFactory1.getBoundaryDofs();
- localProblems_[localProblems_.size()-2] = new OscLocalProblem<MatrixType, VectorType>(fineMat_, rhs1, localToGlobal1, boundaryDofs1, false);
-
- LocalPatchFactory<DGBasisType> patchFactory2(fineBasis_, interfaceNetwork_, seedElements[1]);
- std::vector<int>& localToGlobal2 = patchFactory2.getLocalToGlobal();
- Dune::BitSetVector<1>& boundaryDofs2 = patchFactory2.getBoundaryDofs();
- localProblems_[localProblems_.size()-1] = new OscLocalProblem<MatrixType, VectorType>(fineMat_, rhs2, localToGlobal2, boundaryDofs2, false);
-
- std::cout << std::endl;
- std::cout << "Total setup time: " << timer.elapsed() << " seconds." << std::endl;
- timer.stop();
- }
-
-
- virtual void setProblem(const MatrixType& mat, VectorType& x, const VectorType& rhs) {
- this->x_ = &x;
- this->rhs_ = &rhs;
- this->mat_ = Dune::stackobject_to_shared_ptr(mat);
-
- VectorType contFineRhs;
- fineBasis_.restrict(*(this->rhs_), contFineRhs);
- mgTransfer_.restrict(contFineRhs, coarseRhs_);
-
- for(size_t i=0; i<coarseBoundaryDofs_.size(); i++) {
- if(!coarseBoundaryDofs_[i][0])
- continue;
-
- coarseRhs_[i] = 0;
- }
-
- for (size_t i=0; i<localProblems_.size(); i++) {
- localProblems_[i]->updateRhs(*(this->rhs_));
- }
- }
-
- virtual void iterate() {
-
- // compute exact coarse solution
- *(this->x_) = 0;
-
- VectorType coarseSol(coarseStiffMat_.M());
- #if HAVE_SUPERLU
- Dune::InverseOperatorResult res;
-
- Dune::SuperLU<MatrixType> coarseSolver(coarseStiffMat_);
- coarseSolver.apply(coarseSol, coarseRhs_, res);
- #else
- #error No SuperLU!
- #endif
-
- // print(coarseRhs_, "coarseRhs_: ");
- // print(coarseSol, "coarseSol just on coarse lvl: ");
-
- VectorType contFineSol;
- mgTransfer_.prolong(coarseSol, contFineSol);
-
- //prolong to discontinous fineLevel basis
- fineBasis_.prolong(contFineSol, *(this->x_));
-
-
- // print(*(this->x_), "coarseSol: ");
-
- VectorType it, x;
- for (size_t i=0; i<localProblems_.size(); i++) {
-
- localProblems_[i]->solve(it);
- localProblems_[i]->prolong(it, x);
-
- *(this->x_) += x;
- }
-
- //print(*(this->x_), "globalX final: ");
-
- }
-
- DGBasisType& fineBasis() const {
- return fineBasis_;
- }
-
-};
-
-#endif
-
diff --git a/dune/faultnetworks/preconditioners/faultpreconditioner2.hh b/dune/faultnetworks/preconditioners/faultpreconditioner2.hh
deleted file mode 100644
index 91ad5ed59bb5fecb9efbf1e12d17ddd5a8fe6a55..0000000000000000000000000000000000000000
--- a/dune/faultnetworks/preconditioners/faultpreconditioner2.hh
+++ /dev/null
@@ -1,284 +0,0 @@
-#ifndef FAULT_PRECONDITIONER_HH
-#define FAULT_PRECONDITIONER_HH
-
-#include <queue>
-#include <string>
-
-#include <dune/common/timer.hh>
-#include <dune/common/fvector.hh>
-
-#include <dune/solvers/iterationsteps/lineariterationstep.hh>
-#include <dune/faultnetworks/compressedmultigridtransfer.hh>
-
-#include <dune/faultnetworks/patchfactories/localpatchfactory.hh>
-//#include <dune/osccoeff/osclocalproblemeff_latestWorking.hh>
-#include <dune/faultnetworks/localproblem.hh>
-#include <dune/faultnetworks/assemblers/osclocalassembler.hh>
-#include <dune/faultnetworks/assemblers/oscl2functionalassembler.hh>
-#include <dune/faultnetworks/hierarchicleveliterator.hh>
-#include <dune/faultnetworks/oscrhs.hh>
-#include <dune/faultnetworks/utils/debugutils.hh>
-#include <dune/faultnetworks/dgmgtransfer.hh>
-
-#include <dune/fufem/assemblers/assembler.hh>
-#include <dune/fufem/boundarypatch.hh>
-#include <dune/fufem/functiontools/boundarydofs.hh>
-
-#include <dune/grid/common/intersectioniterator.hh>
-
-#include <dune/istl/superlu.hh>
-
-#include <dune/faultnetworks/levelinterfacenetwork.hh>
-
-/* --------------------------------------------------------------------------------
- *
- * Deprecated: LevelInterfaceNetwork does not support individual interfaces anymore
- *
- * --------------------------------------------------------------------------------
- */
-
-template <class GridType, class CoarseBasisType, class DGBasisType, class MatrixType, class VectorType>
-class FaultPreconditioner : public LinearIterationStep<MatrixType, VectorType> {
-
-public:
- typedef GridType grid_type;
- typedef MatrixType matrix_type;
-
-private:
- static const int dim= GridType::dimension;
- static const int dimworld = GridType::dimensionworld;
- typedef typename GridType::ctype ctype;
- typedef typename Dune::FieldVector<ctype, dimworld> WorldVectorType;
- typedef typename GridType::template Codim<0>::Entity EntityType;
- typedef HierarchicLevelIterator<GridType> HierarchicLevelIteratorType;
-
-
- const GridType& grid_;
-
- const int coarseLevel_;
- const int fineLevel_;
-
- MatrixType& fineMat_;
- MatrixType& coarseStiffMat_;
-
- const LevelInterfaceNetwork<typename GridType::LevelGridView>& interfaceNetwork_;
-
- const DGMGTransfer<DGBasisType>& mgTransfer_;
-
- const bool debugMode_;
-
- //std::vector<OscLocalProblem<BasisType, MatrixType, VectorType, VectorType>* > localProblems_;
- std::vector<OscLocalProblem<MatrixType, VectorType>* > localProblems_;
-
- DGBasisType fineBasis_;
-
- VectorType coarseRhs_;
- Dune::BitSetVector<1> coarseBoundaryDofs_;
-
- void print(const VectorType& x, std::string message){
- std::cout << message << std::endl;
- for (size_t i=0; i<x.size(); i++) {
- std::cout << x[i] << " ";
- }
- std::cout << std::endl << std::endl;
- }
-
- void print(const std::vector<int>& x, std::string message){
- std::cout << message << std::endl;
- for (size_t i=0; i<x.size(); i++) {
- std::cout << x[i] << " ";
- }
- std::cout << std::endl << std::endl;
- }
-
-
- void print(const Dune::BitSetVector<1>& x, std::string message){
- std::cout << message << std::endl;
- for (size_t i=0; i<x.size(); i++) {
- std::cout << x[i][0] << " ";
- }
- std::cout << std::endl << std::endl;
- }
-
-public:
-
- FaultPreconditioner(const GridType& grid, int coarseLevel, int fineLevel, MatrixType& fineMat,
- const LevelInterfaceNetwork<typename GridType::LevelGridView>& interfaceNetwork, MatrixType& coarseMat,
- const DGMGTransfer<DGBasisType>& mgTransfer, const bool debugMode)
- : grid_(grid), coarseLevel_(coarseLevel), fineLevel_(fineLevel), fineMat_(fineMat),
- coarseStiffMat_(coarseMat), interfaceNetwork_(interfaceNetwork),
- mgTransfer_(mgTransfer), debugMode_(debugMode), fineBasis_(grid_.levelGridView(fineLevel_), interfaceNetwork) {
-
- // GV coarseGridView(grid_.levelGridView(coarseLevel_));
- // GV fineGridView(grid_.levelGridView(fineLevel_));
-
- localProblems_.resize(interfaceNetwork_.size()+2);
-
- std::cout << "Initializing local fine grid corrections! " << std::endl;
-
- typedef typename GridType::template Codim<0>::Entity Element;
- typedef typename GridType::LevelIntersection Intersection;
- std::vector<Element> faultZoneSeeds;
- std::vector<int> faultZoneSeedsDofs;
- std::vector<int> faultZoneSeedsCounter;
- for (size_t i=0; i<interfaceNetwork_.size(); i++) {
- const std::vector<Intersection>& faces = interfaceNetwork_.getInterface(i)->faces();
-
- const Element& insideElem = faces[0].inside();
- faultZoneSeeds.push_back(insideElem);
- faultZoneSeedsDofs.push_back(fineBasis_.index(insideElem, 0));
-
- const Element& outsideElem = faces[0].outside();
- faultZoneSeeds.push_back(outsideElem);
- faultZoneSeedsDofs.push_back(fineBasis_.index(outsideElem, 0));
- }
- faultZoneSeedsCounter.resize(faultZoneSeeds.size(), 0);
-
- // init local fine level corrections
- Timer timer;
- timer.reset();
- timer.start();
-
-
- for (size_t i=0; i<interfaceNetwork_.size(); i++) {
- // std::cout << i << std::endl;
- // std::cout << "---------------" << std::endl;
- LocalPatchFactory<DGBasisType> patchFactory(fineBasis_, interfaceNetwork_, i);
-
- std::vector<int>& localToGlobal = patchFactory.getLocalToGlobal();
- Dune::BitSetVector<1>& boundaryDofs = patchFactory.getBoundaryDofs();
-
- for (size_t j=0; j<faultZoneSeeds.size(); j++) {
- if (std::find(localToGlobal.begin(), localToGlobal.end(), faultZoneSeedsDofs[j]) != localToGlobal.end()) {
- faultZoneSeedsCounter[j]++;
- }
- }
-
-
- // print(localToGlobal, "localToGlobal: ");
- // print(boundaryDofs, "boundaryDofs: ");
-
-/* std::vector<EntityType>& regionElements = patchFactory.getRegionElements();
-
-
- Dune::BitSetVector<1> fineBoundaryDofs(fineRhs_.size());
- fineBoundaryDofs.unsetAll();
- for (size_t j=0; j<localToGlobal.size(); j++) {
- if (boundaryDofs[j][0]) {
- fineBoundaryDofs[localToGlobal[j]][0] = true;
- }
- }
-
- OscRhs<WorldVectorType, Dune::FieldVector<ctype,1>> f;
- OscLocalAssembler<GridType, FE, FE> oscLocalAssembler(oscData_, mapper_, coarseLocalProblemRegion);
- OscL2FunctionalAssembler<GridType, FE> oscLocalFunctionalAssembler(f, coarseLocalProblemRegion);
- OscLocalProblem<BasisType, MatrixType, VectorType, VectorType>* localProblem = new OscLocalProblem<BasisType, MatrixType, VectorType, VectorType>(regionElements, fineBasis, fineBoundaryDofs, oscLocalAssembler, fineRhs_, oscLocalFunctionalAssembler);
-*/
-
-
- VectorType rhs;
- rhs.resize(fineMat_.N());
- rhs = 0;
-
- localProblems_[i] = new OscLocalProblem<MatrixType, VectorType>(fineMat_, rhs, localToGlobal, boundaryDofs, false);
-
- //print(localToGlobal, "localToGlobal: ");
- //print(boundaryDofs, "boundaryDofs: ");
-
- if ((i+1) % 10 == 0) {
- std::cout << '\r' << std::floor((i+1.0)/interfaceNetwork_.size()*100) << " % done. Elapsed time: " << timer.elapsed() << " seconds. Predicted total setup time: " << timer.elapsed()/(i+1.0)*interfaceNetwork_.size() << " seconds." << std::flush;
- }
- }
-
- std::vector<Element> seedElements;
- for (size_t j=0; j<faultZoneSeedsCounter.size(); j++) {
- if (faultZoneSeedsCounter[j]==1)
- seedElements.push_back(faultZoneSeeds[j]);
- }
-
- assert(seedElements.size()==2);
-
- VectorType rhs1, rhs2;
- rhs1.resize(fineMat_.N());
- rhs2.resize(fineMat_.N());
- rhs1 = 0;
- rhs2 = 0;
-
- LocalPatchFactory<DGBasisType> patchFactory1(fineBasis_, interfaceNetwork_, seedElements[0]);
- std::vector<int>& localToGlobal1 = patchFactory1.getLocalToGlobal();
- Dune::BitSetVector<1>& boundaryDofs1 = patchFactory1.getBoundaryDofs();
- localProblems_[localProblems_.size()-2] = new OscLocalProblem<MatrixType, VectorType>(fineMat_, rhs1, localToGlobal1, boundaryDofs1, false);
-
- LocalPatchFactory<DGBasisType> patchFactory2(fineBasis_, interfaceNetwork_, seedElements[1]);
- std::vector<int>& localToGlobal2 = patchFactory2.getLocalToGlobal();
- Dune::BitSetVector<1>& boundaryDofs2 = patchFactory2.getBoundaryDofs();
- localProblems_[localProblems_.size()-1] = new OscLocalProblem<MatrixType, VectorType>(fineMat_, rhs2, localToGlobal2, boundaryDofs2, false);
-
- std::cout << std::endl;
- std::cout << "Total setup time: " << timer.elapsed() << " seconds." << std::endl;
- timer.stop();
- }
-
-
- virtual void setProblem(const MatrixType& mat, VectorType& x, const VectorType& rhs) {
- this->x_ = &x;
- this->rhs_ = &rhs;
- this->mat_ = Dune::stackobject_to_shared_ptr(mat);
-
- mgTransfer_.restrict(*(this->rhs_), coarseRhs_);
-
- for(size_t i=0; i<coarseBoundaryDofs_.size(); i++) {
- if(!coarseBoundaryDofs_[i][0])
- continue;
-
- coarseRhs_[i] = 0;
- }
-
- for (size_t i=0; i<localProblems_.size(); i++) {
- localProblems_[i]->updateRhs(*(this->rhs_));
- }
- }
-
- virtual void iterate() {
-
- // compute exact coarse solution
- *(this->x_) = 0;
-
- VectorType coarseSol(coarseStiffMat_.M());
- #if HAVE_SUPERLU
- Dune::InverseOperatorResult res;
-
- Dune::SuperLU<MatrixType> coarseSolver(coarseStiffMat_);
- coarseSolver.apply(coarseSol, coarseRhs_, res);
- #else
- #error No SuperLU!
- #endif
-
- // print(coarseRhs_, "coarseRhs_: ");
- // print(coarseSol, "coarseSol just on coarse lvl: ");
-
- mgTransfer_.prolong(coarseSol, *(this->x_));
-
- // print(*(this->x_), "coarseSol: ");
-
- VectorType it, x;
- for (size_t i=0; i<localProblems_.size(); i++) {
-
- localProblems_[i]->solve(it);
- localProblems_[i]->prolong(it, x);
-
- *(this->x_) += x;
- }
-
- //print(*(this->x_), "globalX final: ");
-
- }
-
- DGBasisType& fineBasis() const {
- return fineBasis_;
- }
-
-};
-
-#endif
-
diff --git a/dune/faultnetworks/preconditioners/levelfaultpreconditioner.hh b/dune/faultnetworks/preconditioners/levelfaultpreconditioner.hh
deleted file mode 100644
index 85827f593f68e519c5fc91da2f1e37fc2e754279..0000000000000000000000000000000000000000
--- a/dune/faultnetworks/preconditioners/levelfaultpreconditioner.hh
+++ /dev/null
@@ -1,223 +0,0 @@
-#ifndef LEVEL_FAULT_PRECONDITIONER_HH
-#define LEVEL_FAULT_PRECONDITIONER_HH
-
-#include <string>
-
-#include <dune/common/timer.hh>
-#include <dune/common/fvector.hh>
-#include <dune/common/bitsetvector.hh>
-
-#include <dune/solvers/iterationsteps/lineariterationstep.hh>
-
-#include <dune/faultnetworks/assemblers/globalfaultassembler.hh>
-#include <dune/faultnetworks/patchfactories/localpatchfactory.hh>
-#include <dune/faultnetworks/localproblem.hh>
-#include <dune/faultnetworks/levelinterfacenetwork.hh>
-#include <dune/faultnetworks/utils/debugutils.hh>
-
-#include <dune/fufem/boundarypatch.hh>
-#include <dune/fufem/functiontools/boundarydofs.hh>
-
-
-/* --------------------------------------------------------------------------------
- *
- * Deprecated: LevelInterfaceNetwork does not support individual interfaces anymore
- *
- * --------------------------------------------------------------------------------
- */
-
-template <class BasisType, class LocalAssembler, class LocalInterfaceAssembler, class MatrixType, class VectorType>
-class LevelFaultPreconditioner : public LinearIterationStep<MatrixType, VectorType> {
-
-public:
- enum BoundaryMode {homogeneous, fromIterate};
-
-private:
- typedef typename BasisType::GridView GridView;
- typedef typename GridView::Grid GridType;
-
- const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork_;
- const LocalAssembler& localAssembler_;
- const std::vector<std::shared_ptr<LocalInterfaceAssembler>>& localInterfaceAssemblers_;
-
- const GridType& grid_;
- const int level_;
- const BasisType basis_;
-
- Dune::BitSetVector<1> activeFaults_;
- size_t patchDepth_;
- BoundaryMode boundaryMode_;
-
- MatrixType matrix_;
- std::vector<std::shared_ptr<OscLocalProblem<MatrixType, VectorType>> > localProblems_;
-
-public:
-
- // for each active fault in levelInterfaceNetwork: set up local problem, compute local correction
- LevelFaultPreconditioner(const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork,
- const LocalAssembler& localAssembler,
- const std::vector<std::shared_ptr<LocalInterfaceAssembler>>& localInterfaceAssemblers) :
- levelInterfaceNetwork_(levelInterfaceNetwork),
- localAssembler_(localAssembler),
- localInterfaceAssemblers_(localInterfaceAssemblers),
- grid_(levelInterfaceNetwork_.grid()),
- level_(levelInterfaceNetwork_.level()),
- basis_(levelInterfaceNetwork_)
- {
- setActiveFaults();
- setPatchDepth();
- setBoundaryMode();
-
- assert(localInterfaceAssemblers_.size() == levelInterfaceNetwork_.size());
- }
-
- void build() {
- // assemble stiffness matrix for entire level including all faults
- GlobalFaultAssembler<BasisType, BasisType> globalFaultAssembler(basis_, basis_, levelInterfaceNetwork_);
- globalFaultAssembler.assembleOperator(localAssembler_, localInterfaceAssemblers_, matrix_);
-
- // set boundary conditions
- Dune::BitSetVector<1> globalBoundaryDofs;
- BoundaryPatch<GridView> boundaryPatch(levelInterfaceNetwork_.levelGridView(), true);
- constructBoundaryDofs(boundaryPatch, basis_, globalBoundaryDofs);
-
- typedef typename MatrixType::row_type RowType;
- typedef typename RowType::ConstIterator ColumnIterator;
-
- for(size_t i=0; i<globalBoundaryDofs.size(); i++) {
- if(!globalBoundaryDofs[i][0])
- continue;
-
- RowType& row = matrix_[i];
-
- ColumnIterator cIt = row.begin();
- ColumnIterator cEndIt = row.end();
-
- for(; cIt!=cEndIt; ++cIt) {
- row[cIt.index()] = 0;
- }
- row[i] = 1;
- }
-
- localProblems_.resize(levelInterfaceNetwork_.size());
-
- std::cout << std::endl;
- std::cout << "---------------------------------------------" << std::endl;
- std::cout << "Initializing local fine grid corrections! " << std::endl;
-
- // init local fine level corrections
- Dune::Timer timer;
- timer.reset();
- timer.start();
-
- for (size_t i=0; i<levelInterfaceNetwork_.size(); i++) {
- if (activeFaults_[i][0]) {
- Dune::BitSetVector<1> locallyActiveFaults(activeFaults_);
-
- // turn off neighboring two active faults for local problem
- int depthCount = 0;
- for (size_t j= (size_t) i +1; j<locallyActiveFaults.size() && depthCount < (int) patchDepth_; ++j) {
- if (locallyActiveFaults[j][0]) {
- locallyActiveFaults[j][0] = false;
- depthCount++;
- }
- }
-
- depthCount = 0;
- for (int j=i-1; j>=0 && depthCount < (int) patchDepth_; --j) {
- if (locallyActiveFaults[j][0]) {
- locallyActiveFaults[j][0] = false;
- depthCount++;
- }
- }
-
- LocalPatchFactory<BasisType> patchFactory(basis_, levelInterfaceNetwork_, i, locallyActiveFaults);
-
- std::vector<int>& localToGlobal = patchFactory.getLocalToGlobal();
- Dune::BitSetVector<1>& boundaryDofs = patchFactory.getBoundaryDofs();
-
- VectorType rhs;
- rhs.resize(matrix_.N());
- rhs = 0;
-
- localProblems_[i] = std::make_shared<OscLocalProblem<MatrixType, VectorType>>(matrix_, rhs, localToGlobal, boundaryDofs);
-
- if ((i+1) % 10 == 0) {
- std::cout << '\r' << std::floor((i+1.0)/levelInterfaceNetwork_.size()*100) << " % done. Elapsed time: " << timer.elapsed() << " seconds. Predicted total setup time: " << timer.elapsed()/(i+1.0)*levelInterfaceNetwork_.size() << " seconds." << std::flush;
- }
- }
-
- }
-
- std::cout << std::endl;
- std::cout << "Total setup time: " << timer.elapsed() << " seconds." << std::endl;
- std::cout << "---------------------------------------------" << std::endl;
- timer.stop();
- }
-
- void setActiveFaults() {
- Dune::BitSetVector<1> activeFaults(levelInterfaceNetwork_.size(), true);
- setActiveFaults(activeFaults);
- }
-
- void setActiveFaults(const Dune::BitSetVector<1> activeFaults) {
- assert(levelInterfaceNetwork_.size() == activeFaults.size());
- activeFaults_ = activeFaults;
- }
-
- void setPatchDepth(const size_t patchDepth = 0) {
- patchDepth_ = patchDepth;
- }
-
- void setBoundaryMode(const BoundaryMode boundaryMode = LevelFaultPreconditioner::BoundaryMode::homogeneous) {
- boundaryMode_ = boundaryMode;
- }
-
- virtual void setProblem(const MatrixType& mat, VectorType& x, const VectorType& rhs) {
- this->x_ = &x;
- this->rhs_ = &rhs;
- this->mat_ = Dune::stackobject_to_shared_ptr(mat);
-
- for (size_t i=0; i<localProblems_.size(); i++) {
- if (boundaryMode_ == BoundaryMode::homogeneous)
- localProblems_[i]->updateRhs(rhs);
- else
- localProblems_[i]->updateRhsAndBoundary(rhs, x);
- }
- }
-
- virtual void iterate() {
- *(this->x_) = 0;
-
- VectorType it, x;
- for (size_t i=0; i<localProblems_.size(); i++) {
- localProblems_[i]->solve(it);
- localProblems_[i]->prolong(it, x);
-
- if (i==5) {
- writeToVTK(basis_, x, "/storage/mi/podlesny/data/faultnetworks/iterates/", "exactvertexdata_patchDepth_"+std::to_string(patchDepth_));
- }
-
- *(this->x_) += x;
- }
- }
-
- const BasisType& basis() const {
- return basis_;
- }
-
- const GridView& gridView() const {
- return basis_.getGridView();
- }
-
- const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork() const {
- return levelInterfaceNetwork_;
- }
-
- size_t size() const {
- return localProblems_.size();
- }
-};
-
-#endif
-
diff --git a/dune/faultnetworks/preconditioners/levelfaultpreconditioner_latestworking.hh b/dune/faultnetworks/preconditioners/levelfaultpreconditioner_latestworking.hh
deleted file mode 100644
index e789eb57bd8e2388b60aea20037c68fc4ddda186..0000000000000000000000000000000000000000
--- a/dune/faultnetworks/preconditioners/levelfaultpreconditioner_latestworking.hh
+++ /dev/null
@@ -1,322 +0,0 @@
-#ifndef LEVEL_FAULT_PRECONDITIONER_HH
-#define LEVEL_FAULT_PRECONDITIONER_HH
-
-#include <string>
-
-#include <dune/common/timer.hh>
-#include <dune/common/fvector.hh>
-
-#include <dune/solvers/iterationsteps/lineariterationstep.hh>
-#include <dune/faultnetworks/compressedmultigridtransfer.hh>
-
-#include <dune/faultnetworks/patchfactories/localpatchfactory.hh>
-#include <dune/faultnetworks/localproblem.hh>
-
-#include <dune/faultnetworks/assemblers/osclocalassembler.hh>
-//#include <dune/faultnetworks/assemblers/oscl2functionalassembler.hh>
-#include <dune/faultnetworks/assemblers/interfaceoperatorassembler.hh>
-
-#include <dune/faultnetworks/hierarchicleveliterator.hh>
-#include <dune/faultnetworks/oscrhs.hh>
-#include <dune/faultnetworks/utils/debugutils.hh>
-
-#include <dune/fufem/assemblers/assembler.hh>
-#include <dune/fufem/boundarypatch.hh>
-#include <dune/fufem/functiontools/boundarydofs.hh>
-
-#include <dune/grid/common/mcmgmapper.hh>
-#include <dune/grid/common/intersectioniterator.hh>
-
-#include <dune/istl/superlu.hh>
-#include <dune/istl/matrixindexset.hh>
-
-#include <dune/faultnetworks/levelinterfacenetwork.hh>
-#include <dune/faultnetworks/oscdata.hh>
-
-
-template <class BasisType, class LocalAssembler, class LocalInterfaceAssembler, class MatrixType, class VectorType>
-class LevelFaultPreconditioner : public LinearIterationStep<MatrixType, VectorType> {
-
-private:
- typedef typename BasisType::GridView GridView;
- typedef typename GridView::Grid GridType;
-
- static const int dim= GridType::dimension;
- static const int dimworld = GridType::dimensionworld;
-
- typedef typename GridType::ctype ctype;
- typedef typename Dune::FieldVector<ctype, dimworld> WorldVectorType;
- typedef typename Dune::LevelMultipleCodimMultipleGeomTypeMapper <GridType, Dune::MCMGElementLayout > MapperType;
- typedef typename GridType::template Codim<0>::Entity EntityType;
- typedef HierarchicLevelIterator<GridType> HierarchicLevelIteratorType;
-
-
- const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork_;
- const LocalAssembler& localAssembler_;
- const LocalInterfaceAssembler& localInterfaceAssembler_;
- const GridType& grid_;
- const int level_;
-
-
- MatrixType matrix_;
-
- const bool debugMode_;
-
- const GridView gridView_;
- const BasisType basis_;
-
- std::vector<std::shared_ptr<OscLocalProblem<MatrixType, VectorType>> > localProblems_;
-
- Dune::BitSetVector<1> boundaryDofs_;
-
-
-protected:
-
- // res = A + B
- void add(MatrixType& A, MatrixType& B, MatrixType& res) {
- if (A.N()!=B.N() || A.M()!=B.M())
- DUNE_THROW(Dune::RangeError, "levelFaultPreconditioner::add() Matrices must have same dimensions!");
-
- typedef typename std::map<std::pair<size_t, size_t>, ctype> MatEntriesType;
- typedef typename MatEntriesType::iterator MatEntriesIterator;
- typedef typename MatrixType::row_type RowType;
- typedef typename RowType::ConstIterator ColumnIterator;
-
- MatEntriesType matEntries;
- MatEntriesIterator entryIt;
-
- Dune::MatrixIndexSet joinedIdxSet(A.N(), A.M());
- joinedIdxSet.import(A);
-
- for(size_t rowIdx = 0; rowIdx<B.N(); rowIdx++) {
- const RowType& BRow = B[rowIdx];
- const RowType& ARow = A[rowIdx];
-
- ColumnIterator cIt = BRow.begin();
- ColumnIterator cEndIt = BRow.end();
-
- for(; cIt!=cEndIt; ++cIt) {
- const size_t colIdx = cIt.index();
-
- if (A.exists(rowIdx, colIdx)) {
- matEntries[std::make_pair(rowIdx, colIdx)] = ARow[colIdx] + BRow[colIdx];
- } else {
- matEntries[std::make_pair(rowIdx, colIdx)] = BRow[colIdx];
- joinedIdxSet.add(rowIdx, colIdx);
- }
- }
- }
-
- joinedIdxSet.exportIdx(res);
-
- for(size_t rowIdx = 0; rowIdx<res.N(); rowIdx++) {
- RowType& row = res[rowIdx];
- const RowType& ARow = A[rowIdx];
-
- ColumnIterator cIt = row.begin();
- ColumnIterator cEndIt = row.end();
-
- for(; cIt!=cEndIt; ++cIt) {
- const size_t colIdx = cIt.index();
- entryIt = matEntries.find(std::make_pair(rowIdx, colIdx));
-
- if (entryIt==matEntries.end()) {
- row[colIdx] = ARow[colIdx];
- } else {
- row[colIdx] = entryIt->second;
- }
- }
- }
- }
-
-public:
-
- // all faults in levelInterfaceNetwork are active and seeds for a local problem, employs 2 additional smooth problems at boundary without faults for improved convergence
- LevelFaultPreconditioner(const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork, const LocalAssembler& localAssembler, const LocalInterfaceAssembler& localInterfaceAssembler, const bool debugMode=false)
- : levelInterfaceNetwork_(levelInterfaceNetwork),
- localAssembler_(localAssembler),
- localInterfaceAssembler_(localInterfaceAssembler),
- grid_(levelInterfaceNetwork_.grid()),
- level_(levelInterfaceNetwork_.level()),
- debugMode_(debugMode),
- gridView_(grid_.levelGridView(level_)),
- basis_(gridView_, levelInterfaceNetwork_) {
-
-
- // assemble stiffness matrix
- MatrixType bilinStiffMat, interfaceMats;
- Assembler<BasisType, BasisType> assembler(basis_, basis_);
- assembler.assembleOperator(localAssembler_, bilinStiffMat);
-
- // assemble interface integrals
- Dune::MatrixIndexSet idxSet(bilinStiffMat.N(), bilinStiffMat.M());
- idxSet.exportIdx(interfaceMats);
-
- for (size_t i=0; i<levelInterfaceNetwork_.size(); i++) {
- MatrixType interfaceMat, res;
- const FaultInterface<GridView>& fault = *(levelInterfaceNetwork_.getInterface(i));
-
- InterfaceOperatorAssembler<BasisType, BasisType> interfaceAssembler(basis_, basis_, fault);
- interfaceAssembler.assemble(localInterfaceAssembler_, interfaceMat);
-
- //print(exactInterfaceMat, "exactInterfaceMat:");
-
- add(interfaceMats, interfaceMat, res);
- interfaceMats = res;
- }
-
- add(bilinStiffMat, interfaceMats, matrix_);
-
- // set boundary conditions
- BoundaryPatch<GridView> boundaryPatch(gridView_, true);
- constructBoundaryDofs(boundaryPatch, basis_, boundaryDofs_);
-
- typedef typename MatrixType::row_type RowType;
- typedef typename RowType::ConstIterator ColumnIterator;
-
- for(size_t i=0; i<boundaryDofs_.size(); i++) {
- if(!boundaryDofs_[i][0])
- continue;
-
- RowType& row = matrix_[i];
-
- ColumnIterator cIt = row.begin();
- ColumnIterator cEndIt = row.end();
-
- for(; cIt!=cEndIt; ++cIt) {
- row[cIt.index()] = 0;
- }
-
- row[i] = 1;
- }
-
- localProblems_.resize(levelInterfaceNetwork_.size()+2);
-
- std::cout << "Initializing local fine grid corrections! " << std::endl;
- typedef typename GridType::template Codim<0>::Entity Element;
- typedef typename GridType::LevelIntersection Intersection;
- std::vector<Element> faultZoneSeeds;
- std::vector<int> faultZoneSeedsDofs;
- std::vector<int> faultZoneSeedsCounter;
-
- for (size_t i=0; i<levelInterfaceNetwork_.size(); i++) {
- const std::vector<Intersection>& faces = levelInterfaceNetwork_.getInterface(i)->faces();
-
- const Element& insideElem = faces[0].inside();
- faultZoneSeeds.push_back(insideElem);
- faultZoneSeedsDofs.push_back(basis_.index(insideElem, 0));
-
- const Element& outsideElem = faces[0].outside();
- faultZoneSeeds.push_back(outsideElem);
- faultZoneSeedsDofs.push_back(basis_.index(outsideElem, 0));
- }
- faultZoneSeedsCounter.resize(faultZoneSeeds.size(), 0);
-
- // init local fine level corrections
- Dune::Timer timer;
- timer.reset();
- timer.start();
-
- for (size_t i=0; i<levelInterfaceNetwork_.size(); i++) {
- // std::cout << i << std::endl;
- // std::cout << "---------------" << std::endl;
- LocalPatchFactory<BasisType> patchFactory(basis_, levelInterfaceNetwork_, i);
-
- std::vector<int>& localToGlobal = patchFactory.getLocalToGlobal();
- Dune::BitSetVector<1>& boundaryDofs = patchFactory.getBoundaryDofs();
-
- for (size_t j=0; j<faultZoneSeeds.size(); j++) {
- if (std::find(localToGlobal.begin(), localToGlobal.end(), faultZoneSeedsDofs[j]) != localToGlobal.end()) {
- faultZoneSeedsCounter[j]++;
- }
- }
-
- // print(localToGlobal, "localToGlobal: ");
- // print(boundaryDofs, "boundaryDofs: ");
-
- VectorType rhs;
- rhs.resize(matrix_.N());
- rhs = 0;
-
- localProblems_[i] = std::make_shared<OscLocalProblem<MatrixType, VectorType>>(matrix_, rhs, localToGlobal, boundaryDofs, false);
-
- //print(localToGlobal, "localToGlobal: ");
- //print(boundaryDofs, "boundaryDofs: ");
-
- if ((i+1) % 10 == 0) {
- std::cout << '\r' << std::floor((i+1.0)/levelInterfaceNetwork_.size()*100) << " % done. Elapsed time: " << timer.elapsed() << " seconds. Predicted total setup time: " << timer.elapsed()/(i+1.0)*levelInterfaceNetwork_.size() << " seconds." << std::flush;
- }
- }
-
- std::vector<Element> seedElements;
- for (size_t j=0; j<faultZoneSeedsCounter.size(); j++) {
- if (faultZoneSeedsCounter[j]==1)
- seedElements.push_back(faultZoneSeeds[j]);
- }
-
- assert(seedElements.size()==2);
-
- VectorType rhs1, rhs2;
- rhs1.resize(matrix_.N());
- rhs2.resize(matrix_.N());
- rhs1 = 0;
- rhs2 = 0;
-
- LocalPatchFactory<BasisType> patchFactory1(basis_, levelInterfaceNetwork_, seedElements[0]);
- std::vector<int>& localToGlobal1 = patchFactory1.getLocalToGlobal();
- Dune::BitSetVector<1>& boundaryDofs1 = patchFactory1.getBoundaryDofs();
- localProblems_[localProblems_.size()-2] = std::make_shared<OscLocalProblem<MatrixType, VectorType>>(matrix_, rhs1, localToGlobal1, boundaryDofs1, false);
-
- LocalPatchFactory<BasisType> patchFactory2(basis_, levelInterfaceNetwork_, seedElements[1]);
- std::vector<int>& localToGlobal2 = patchFactory2.getLocalToGlobal();
- Dune::BitSetVector<1>& boundaryDofs2 = patchFactory2.getBoundaryDofs();
- localProblems_[localProblems_.size()-1] = std::make_shared<OscLocalProblem<MatrixType, VectorType>>(matrix_, rhs2, localToGlobal2, boundaryDofs2, false);
-
- std::cout << std::endl;
- std::cout << "Total setup time: " << timer.elapsed() << " seconds." << std::endl;
- timer.stop();
- }
-
-
- virtual void setProblem(const MatrixType& mat, VectorType& x, const VectorType& rhs) {
- this->x_ = &x;
- this->rhs_ = &rhs;
- //this->mat_ = Dune::stackobject_to_shared_ptr(mat);
-
- for (size_t i=0; i<localProblems_.size(); i++) {
- localProblems_[i]->updateRhs(rhs);
- }
- }
-
- virtual void iterate() {
- *(this->x_) = 0;
-
- VectorType it, x;
- for (size_t i=0; i<localProblems_.size(); i++) {
- localProblems_[i]->solve(it);
- localProblems_[i]->prolong(it, x);
-
- //print(it, "it: ");
- //print(x, "x");
-
- *(this->x_) += x;
- }
- //print(*(this->x_), "globalX final: ");
- //print(getSol(), "getSol() output: ");
- }
-
- const BasisType& basis() const {
- return basis_;
- }
-
- const GridView& gridView() const {
- return gridView_;
- }
-
- /*VectorType getSol(){
- return *(this->x_);
- }*/
-};
-
-#endif
-
diff --git a/dune/faultnetworks/preconditioners/multilevelfaultpreconditioner.hh b/dune/faultnetworks/preconditioners/multilevelfaultpreconditioner.hh
deleted file mode 100644
index 730f5056edc3a99a9853f0f7366ed36942c794a3..0000000000000000000000000000000000000000
--- a/dune/faultnetworks/preconditioners/multilevelfaultpreconditioner.hh
+++ /dev/null
@@ -1,179 +0,0 @@
-#ifndef MULTILEVEL_FAULT_PRECONDITIONER_HH
-#define MULTILEVEL_FAULT_PRECONDITIONER_HH
-
-#include <string>
-
-#include <dune/common/timer.hh>
-#include <dune/common/fvector.hh>
-#include <dune/common/bitsetvector.hh>
-
-#include <dune/solvers/iterationsteps/lineariterationstep.hh>
-
-#include <dune/faultnetworks/preconditioners/levelfaultpreconditioner.hh>
-#include <dune/faultnetworks/levelinterfacenetwork.hh>
-#include <dune/faultnetworks/interfacenetwork.hh>
-#include <dune/faultnetworks/dgmgtransfer.hh>
-
-#include <dune/faultnetworks/utils/debugutils.hh>
-
-/* --------------------------------------------------------------------------------
- *
- * Deprecated: LevelInterfaceNetwork does not support individual interfaces anymore
- *
- * --------------------------------------------------------------------------------
- */
-
-template <class BasisType, class LocalOperatorAssembler, class LocalInterfaceAssembler, class MatrixType, class VectorType>
-class MultilevelFaultPreconditioner : public LinearIterationStep<MatrixType, VectorType> {
-
-public:
- enum Mode {additive, multiplicative};
- typedef LevelFaultPreconditioner<BasisType, LocalOperatorAssembler, LocalInterfaceAssembler, MatrixType, VectorType> LevelFaultPreconditionerType;
-
-private:
- typedef typename BasisType::GridView GridView;
- typedef typename GridView::Grid GridType;
-
-
- typedef typename Dune::BitSetVector<1> BitVector;
-
- const InterfaceNetwork<GridType>& interfaceNetwork_;
- const BitVector& activeLevels_;
- const std::vector<std::shared_ptr<LocalOperatorAssembler>>& localAssemblers_;
- const std::vector<std::vector<std::shared_ptr<LocalInterfaceAssembler>>>& localInterfaceAssemblers_;
- const Mode mode_;
-
- std::vector<std::shared_ptr<LevelFaultPreconditionerType>> levelFaultPreconditioners_;
- std::vector<VectorType> levelX_;
- std::vector<VectorType> levelRhs_;
-
- std::shared_ptr<LevelInterfaceNetwork<GridView>> allFaultLevelInterfaceNetwork_;
- std::shared_ptr<BasisType> allFaultBasis_;
-
- std::vector<std::shared_ptr<DGMGTransfer<BasisType>>> mgTransfer_;
-
-public:
- MultilevelFaultPreconditioner(const InterfaceNetwork<GridType>& interfaceNetwork,
- const BitVector& activeLevels,
- const std::vector<std::shared_ptr<LocalOperatorAssembler>>& localAssemblers,
- const std::vector<std::vector<std::shared_ptr<LocalInterfaceAssembler>>>& localInterfaceAssemblers,
- const Mode mode = MultilevelFaultPreconditioner::Mode::additive) :
- interfaceNetwork_(interfaceNetwork),
- activeLevels_(activeLevels),
- localAssemblers_(localAssemblers),
- localInterfaceAssemblers_(localInterfaceAssemblers),
- mode_(mode)
- {
-
- if (activeLevels_.size() > (size_t) interfaceNetwork.maxLevel() +1)
- DUNE_THROW(Dune::Exception, "MultilevelFaultPreconditioner: too many active levels; preconditioner supports at most (grid.maxLevel + 1) levels!");
-
- assert(activeLevels.size() == localAssemblers_.size() && activeLevels.size() == localInterfaceAssemblers_.size());
-
- // init level fault preconditioners and multigrid transfer
- levelFaultPreconditioners_.resize(0);
- mgTransfer_.resize(0);
-
- int maxLevel = 0;
- for (size_t i=0; i<activeLevels_.size(); i++) {
- if (activeLevels_[i][0]) {
- const LevelInterfaceNetwork<GridView>& levelNetwork = interfaceNetwork.levelInterfaceNetwork(i);
-
- levelFaultPreconditioners_.push_back(std::make_shared<LevelFaultPreconditionerType>(levelNetwork, *localAssemblers_[i], localInterfaceAssemblers_[i]));
- maxLevel = i;
- }
- }
-
- levelX_.resize(levelFaultPreconditioners_.size());
- levelRhs_.resize(levelFaultPreconditioners_.size());
-
- allFaultBasis_ = std::make_shared<BasisType>(interfaceNetwork_.levelInterfaceNetwork(maxLevel));
-
- // init multigrid transfer
- for (size_t i=0; i<levelFaultPreconditioners_.size(); i++) {
- mgTransfer_.push_back(std::make_shared<DGMGTransfer<BasisType>>(levelFaultPreconditioners_[i]->basis(), *allFaultBasis_));
- }
- }
-
-
- virtual void setProblem(const MatrixType& mat, VectorType& x, const VectorType& rhs) {
- this->x_ = &x;
- this->rhs_ = &rhs;
- this->mat_ = Dune::stackobject_to_shared_ptr(mat);
-
- for (size_t i=0; i<levelFaultPreconditioners_.size(); i++) {
- mgTransfer_[i]->restrict(x, levelX_[i]);
- mgTransfer_[i]->restrict(rhs, levelRhs_[i]);
-
- levelFaultPreconditioners_[i]->setProblem(mat, levelX_[i], levelRhs_[i]);
- }
- }
-
-
- virtual void iterate() {
- if (mode_ == additive)
- iterateAdd();
- else
- iterateMult();
- }
-
- void iterateAdd() {
- *(this->x_) = 0;
-
- VectorType x;
- for (size_t i=0; i<levelFaultPreconditioners_.size(); i++) {
- levelFaultPreconditioners_[i]->iterate();
- const VectorType& it = levelFaultPreconditioners_[i]->getSol();
-
- mgTransfer_[i]->prolong(it, x);
-
- //writeToVTK(*allFaultBasis_, x, "/home/mi/podlesny/data/faultnetworks/sol/", "preconditioner_step_"+std::to_string(i));
-
- *(this->x_) += x;
- }
- }
-
-
- void iterateMult() {
- *(this->x_) = 0;
-
- VectorType x;
- for (size_t i=0; i<levelFaultPreconditioners_.size(); i++) {
- VectorType updatedRhs(*(this->rhs_));
- this->mat_->mmv(*(this->x_), updatedRhs);
-
- mgTransfer_[i]->restrict(*(this->x_), levelX_[i]);
- mgTransfer_[i]->restrict(updatedRhs, levelRhs_[i]);
-
- levelFaultPreconditioners_[i]->setProblem(*(this->mat_), levelX_[i], levelRhs_[i]);
-
- levelFaultPreconditioners_[i]->iterate();
- const VectorType& it = levelFaultPreconditioners_[i]->getSol();
-
- mgTransfer_[i]->prolong(it, x);
-
- *(this->x_) += x;
- }
- }
-
- void build() {
- for (size_t i=0; i<levelFaultPreconditioners_.size(); i++) {
- levelFaultPreconditioners_[i]->build();
- }
- }
-
- std::shared_ptr<LevelFaultPreconditionerType> levelFaultPreconditioner(const int level) {
- return levelFaultPreconditioners_[level];
- }
-
- const BasisType& basis() const {
- return *allFaultBasis_;
- }
-
- size_t size() const {
- return levelFaultPreconditioners_.size();
- }
-};
-
-#endif
-
diff --git a/dune/faultnetworks/preconditioners/old/oscpreconditionereff_latestWorking.hh b/dune/faultnetworks/preconditioners/old/oscpreconditionereff_latestWorking.hh
deleted file mode 100644
index 4108d3addbd4c71d5f4d82717454c1f21b647010..0000000000000000000000000000000000000000
--- a/dune/faultnetworks/preconditioners/old/oscpreconditionereff_latestWorking.hh
+++ /dev/null
@@ -1,361 +0,0 @@
-#ifndef OSC_PRECONDITIONER_HH
-#define OSC_PRECONDITIONER_HH
-
-#include <queue>
-#include <string>
-
-#include <dune/common/timer.hh>
-#include <dune/common/fvector.hh>
-
-#include <dune/solvers/iterationsteps/lineariterationstep.hh>
-#include <dune/osccoeff/compressedmultigridtransfer.hh>
-
-#include <dune/osccoeff/osclocalproblemeff_latestWorking.hh>
-#include <dune/osccoeff/osclocalproblemregion.hh>
-#include <dune/osccoeff/oscl2functionalassembler.hh>
-#include <dune/osccoeff/hierarchicleveliterator.hh>
-
-#include <dune/fufem/assemblers/assembler.hh>
-#include <dune/fufem/boundarypatch.hh>
-
-#include <dune/grid/common/mcmgmapper.hh>
-#include <dune/grid/common/intersectioniterator.hh>
-
-#include <dune/istl/superlu.hh>
-
-
-
-
-template <class GridType, class BasisType, class MatrixType, class VectorType>
-class OscPreconditioner : public LinearIterationStep<MatrixType, VectorType> {
-
-public:
- typedef GridType grid_type;
- typedef MatrixType matrix_type;
-
-private:
- static const int dim= GridType::dimension;
- static const int dimworld = GridType::dimensionworld;
- typedef typename GridType::ctype ctype;
- typedef typename Dune::FieldVector<ctype, dimworld> WorldVectorType;
- typedef typename Dune::LevelMultipleCodimMultipleGeomTypeMapper <GridType, Dune::MCMGElementLayout > MapperType;
- typedef typename GridType::template Codim<0>::Entity EntityType;
-
- typedef HierarchicLevelIterator<GridType> HierarchicLevelIteratorType;
-
- GridType& grid_;
- std::vector<ctype>& oscData_;
-
- const int coarseLevel_;
- const int fineLevel_;
-
- MapperType& mapper_;
-
- const ctype boxLength_;
-
- VectorType& fineRhs_;
-
- std::vector<OscLocalProblem<BasisType, MatrixType, VectorType, VectorType>* > localProblems_;
-
- MatrixType coarseStiffMat_;
- VectorType coarseRhs_;
- Dune::BitSetVector<1> coarseBoundaryDofs_;
- CompressedMultigridTransfer<VectorType> mgTransfer_;
-
- struct classcomp {
- bool operator() (const WorldVectorType& lhs, const WorldVectorType rhs) const {
-
- for (size_t i=0; i<lhs.size(); ++i) {
- if (lhs[i]>rhs[i])
- return false;
- else if (lhs[i]<rhs[i])
- return true;
- }
- return false;
- }
- };
-
- void print(const VectorType& x, std::string message){
- std::cout << message << std::endl;
- for (size_t i=0; i<x.size(); i++) {
- std::cout << x[i] << " ";
- }
- std::cout << std::endl << std::endl;
- }
-
- void addRegionElements(const EntityType seed, std::vector<EntityType>& regionElements, int& elemCounter){
-
- HierarchicLevelIteratorType endHierIt(grid_, seed, HierarchicLevelIteratorType::PositionFlag::end, fineLevel_);
- for (HierarchicLevelIteratorType hierIt(grid_, seed, HierarchicLevelIteratorType::PositionFlag::begin, fineLevel_); hierIt!=endHierIt; ++hierIt) {
-
- EntityType entity(*hierIt);
- regionElements[elemCounter] = entity;
- elemCounter++;
- }
- }
-
-
-
-public:
-
- OscPreconditioner(GridType& grid, int coarseLevel, int fineLevel, std::vector<ctype>& oscData, MapperType& mapper, ctype boxLength, VectorType& fineRhs)
- : grid_(grid), oscData_(oscData), coarseLevel_(coarseLevel), fineLevel_(fineLevel), mapper_(mapper), boxLength_(boxLength), fineRhs_(fineRhs) {
-
- typedef typename GridType::LevelGridView GV;
- typedef typename BasisType::LocalFiniteElement FE;
-
- GV coarseGridView(grid_.levelGridView(coarseLevel_));
- GV fineGridView(grid_.levelGridView(fineLevel_));
-
- BasisType coarseBasis(coarseGridView);
- BasisType fineBasis(fineGridView);
-
-
-
- std::cout << "Setting coarse grid oscData:" << std::endl;
- Dune::LevelMultipleCodimMultipleGeomTypeMapper <GridType, Dune::MCMGElementLayout > coarseMapper(grid_, coarseLevel_);
-
- // allocate a vector for the data
- std::vector<ctype> coarseOscData(coarseMapper.size());
-
- typedef typename GV::template Codim <0>::Iterator ElementLevelIterator;
-
- std::map<WorldVectorType, EntityType, classcomp> regionSeeds;
-
- ElementLevelIterator endElemIt = coarseGridView.template end <0>();
- for (ElementLevelIterator it = coarseGridView. template begin <0>(); it!=endElemIt; ++it) {
- const typename EntityType::Geometry geometry = it->geometry();
-
- // set region seed elements
- for (int i=0; i<geometry.corners(); i++) {
- WorldVectorType corner = geometry.corner(i);
- typename std::map<WorldVectorType, EntityType>::iterator hasCorner = regionSeeds.find(corner);
-
- if (hasCorner==regionSeeds.end()) {
- EntityType entity(*it);
- regionSeeds.insert(std::pair<WorldVectorType, EntityType>(corner, entity));
- }
- }
-
-
- // set coarseOscData
- int counter = 0;
- ctype sum = 0;
-
- HierarchicLevelIteratorType endHierIt(grid_, *it, HierarchicLevelIteratorType::PositionFlag::end, fineLevel_);
- for (HierarchicLevelIteratorType hierIt(grid_, *it, HierarchicLevelIteratorType::PositionFlag::begin, fineLevel_); hierIt!=endHierIt; ++hierIt) {
-
- sum += oscData_[mapper_.index(*hierIt)];
- counter++;
- }
- coarseOscData[coarseMapper.index(*it)] = sum/counter;
- // std::cout << "elem: " << (int) mapper.map(*it) << "| "<< coarseOscData[coarseMapper.map(*it)] << std::endl;
- }
-
-
- WorldVectorType zero(2);
- zero = 0;
- typename std::map<WorldVectorType, EntityType>::iterator has = regionSeeds.find(zero);
- std::cout << "region size: "<<regionSeeds.size() << std::endl;
- //std::cout << (has==regionSeeds.end())<< " " << std::endl;
-
- std::cout << "Assembling coarse grid stiffness matrix!" << std::endl;
-
- // assemble coarse stiffness matrix
- WorldVectorType midPoint(1/2);
- ctype box = 1;
-
- Assembler<BasisType, BasisType> coarseAssembler(coarseBasis, coarseBasis);
- OscLocalProblemBox<typename GridType::template Codim<0>::Entity, WorldVectorType> coarseLocalProblemRegion(midPoint, box);
- OscLocalAssembler<GridType, FE, FE> coarseLocalAssembler(coarseOscData, coarseMapper, coarseLocalProblemRegion);
- coarseAssembler.assembleOperator(coarseLocalAssembler, coarseStiffMat_);
-
-
- std::cout << "Setting coarse grid boundary conditions!" << std::endl;
-
- typedef typename MatrixType::row_type RowType;
- typedef typename RowType::ConstIterator ColumnIterator;
-
- // set coarse boundary conditions
- BoundaryPatch<typename GridType::LevelGridView> coarseBoundaryPatch(coarseGridView, true);
- constructBoundaryDofs(coarseBoundaryPatch, coarseBasis, coarseBoundaryDofs_);
-
- for(size_t i=0; i<coarseBoundaryDofs_.size(); i++) {
- if(!coarseBoundaryDofs_[i][0])
- continue;
-
- const RowType& row = coarseStiffMat_[i];
-
- ColumnIterator cIt = row.begin();
- ColumnIterator cEndIt = row.end();
-
- for(; cIt!=cEndIt; ++cIt) {
- coarseStiffMat_[i][cIt.index()] = 0;
- }
- coarseStiffMat_[i][i]=1;
- }
-
- mgTransfer_.setup(coarseGridView, fineGridView);
-
- localProblems_.resize(grid_.size(coarseLevel_, dim));
-
- std::cout << "Initializing local fine grid corrections! " << std::endl;
-
-
-
- // init local fine level corrections
- int count = 0;
- OscRhs<WorldVectorType, Dune::FieldVector<ctype,1>> f;
- Dune::BitSetVector<1> fineBoundaryDofs;
- BoundaryPatch<typename GridType::LevelGridView> fineBoundaryPatch(fineGridView, true);
- constructBoundaryDofs(fineBoundaryPatch, fineBasis, fineBoundaryDofs);
-
- typedef typename GV::template Codim <dim>::Iterator VertexLevelIterator;
-
- Timer timer;
- timer.reset();
- timer.start();
- VertexLevelIterator endIt = coarseGridView.template end <dim>();
- for (VertexLevelIterator it = coarseGridView.template begin <dim>(); it!=endIt; ++it) {
-
- const WorldVectorType vertex = it->geometry().corner(0);
- EntityType regionSeed(regionSeeds.at(vertex));
-
- /* std::cout << vertex[0] << " " << vertex[1] << " " << std::endl;
- std::cout << regionSeed->geometry().center().two_norm() << std::endl;*/
-
- OscLocalProblemBox<EntityType, WorldVectorType> fineLocalProblemRegion(vertex, boxLength_);
-
- std::vector<EntityType> regionElements(mapper_.size(), regionSeed);
- int elemCounter = 0;
- addRegionElements(regionSeed, regionElements, elemCounter);
-
- std::queue<EntityType> elementQueue;
- elementQueue.push(regionSeed);
-
- Dune::BitSetVector<1> visited(coarseMapper.size());
- visited.unsetAll();
-
- visited[coarseMapper.index(regionSeed)][0] = true;
-
- while (!elementQueue.empty()) {
-
- const EntityType elem = elementQueue.front();
- elementQueue.pop();
-
- typename GridType::LevelIntersectionIterator neighborIt = coarseGridView.ibegin(elem);
- typename GridType::LevelIntersectionIterator neighborItEnd = coarseGridView.iend(elem);
- for (; neighborIt!=neighborItEnd; ++neighborIt) {
- const typename GridType::LevelIntersection& intersection = *neighborIt;
-
- if (intersection.neighbor()) {
- EntityType neighbor = intersection.outside();
-
- bool inRegion;
- fineLocalProblemRegion.evaluate(neighbor, inRegion);
- int neighborId = coarseMapper.index(neighbor);
-
-
- if (inRegion and !(visited[neighborId][0])) {
- elementQueue.push(neighbor);
- visited[neighborId][0] = true;
- addRegionElements(neighbor, regionElements, elemCounter);
- }
- }
- }
- }
-
- regionElements.erase(regionElements.begin()+elemCounter, regionElements.end());
- //std::cout << regionElements.size() << " " << elemCounter << std::endl;
-
- OscLocalAssembler<GridType, FE, FE> oscLocalAssembler(oscData_, mapper_, fineLocalProblemRegion);
- OscL2FunctionalAssembler<GridType, FE> oscLocalFunctionalAssembler(f, fineLocalProblemRegion);
- OscLocalProblem<BasisType, MatrixType, VectorType, VectorType>* localProblem = new OscLocalProblem<BasisType, MatrixType, VectorType, VectorType>(regionElements, fineBasis, fineBoundaryDofs, oscLocalAssembler, fineRhs_, oscLocalFunctionalAssembler);
-
- localProblems_[count] = localProblem;
-
- if ((count+1) % 10 == 0) {
- std::cout << '\r' << std::floor((count+1.0)/coarseBasis.size()*100) << " % done. Elapsed time: " << timer.elapsed() << " seconds. Predicted total setup time: " << timer.elapsed()/(count+1.0)*coarseBasis.size() << " seconds." << std::flush;
- }
- count++;
- }
-
- std::cout << std::endl;
- std::cout << "Total setup time: " << timer.elapsed() << " seconds." << std::endl;
- timer.stop();
- }
-
-
- virtual void setProblem(const MatrixType& mat, VectorType& x, const VectorType& rhs) {
- this->x_ = &x;
- this->rhs_ = &rhs;
- this->mat_ = Dune::stackobject_to_shared_ptr(mat);
-
- mgTransfer_.restrict(*(this->rhs_), coarseRhs_);
- for(size_t i=0; i<coarseBoundaryDofs_.size(); i++) {
- if(!coarseBoundaryDofs_[i][0])
- continue;
-
- coarseRhs_[i] = 0;
- }
-
- for (size_t i=0; i<localProblems_.size(); i++) {
- localProblems_[i]->updateRhs(*(this->rhs_));
- }
- }
-
- virtual void iterate() {
-
- // compute exact coarse solution
- *(this->x_) = 0;
-
- VectorType coarseSol(coarseStiffMat_.M());
- #if HAVE_SUPERLU
- Dune::InverseOperatorResult res;
-
- Dune::SuperLU<MatrixType> coarseSolver(coarseStiffMat_);
- coarseSolver.apply(coarseSol, coarseRhs_, res);
- #else
- #error No SuperLU!
- #endif
-
- /*std::cout << "coarseStiffMat_: " << std::endl;
- for(size_t rowIdx = 0; rowIdx<localProblems_[6]->stiffnessMat_.N(); rowIdx++) {
- for(size_t colIdx = 0; colIdx<localProblems_[6]->stiffnessMat_.M(); colIdx++) {
- if (localProblems_[6]->stiffnessMat_.exists(rowIdx, colIdx))
- std::cout << localProblems_[6]->stiffnessMat_[rowIdx][colIdx] << " ";
- else
- std::cout << "0 ";
- }
- std::cout << std::endl;
- }
- std::cout << std::endl;*/
-/*
- std::cout << "coarseRhs_: " << std::endl;
- for (size_t i = 0; i<coarseRhs_.size(); i++) {
- std::cout << coarseRhs_[i] << " ";
- }
- std::cout << std::endl;
-
- std::cout << "coarseSol: " << std::endl;
- for (size_t i = 0; i<coarseSol.size(); i++) {
- std::cout << coarseSol[i] << " ";
- }
- std::cout << std::endl;
-*/
- mgTransfer_.prolong(coarseSol, *(this->x_));
-
- // print(*(this->x_), "coarseSol: ");
-
- VectorType x;
- for (size_t i=0; i<localProblems_.size(); i++) {
-
- localProblems_[i]->solve(x);
- *(this->x_) += x;
- }
- }
-
-};
-
-#endif
-
diff --git a/dune/faultnetworks/preconditioners/old/oscpreconditionereffparallel.hh b/dune/faultnetworks/preconditioners/old/oscpreconditionereffparallel.hh
deleted file mode 100644
index 13ac9e830629d3bc09080c6b8fa66c20236024da..0000000000000000000000000000000000000000
--- a/dune/faultnetworks/preconditioners/old/oscpreconditionereffparallel.hh
+++ /dev/null
@@ -1,359 +0,0 @@
-#ifndef OSC_PRECONDITIONER_HH
-#define OSC_PRECONDITIONER_HH
-
-#include <queue>
-#include <thread>
-
-#include <dune/common/timer.hh>
-#include <dune/common/fvector.hh>
-
-#include <dune/solvers/iterationsteps/lineariterationstep.hh>
-#include <osccoeffImpl/compressedmultigridtransfer.hh>
-
-#include <osccoeffImpl/osclocalproblemeff.hh>
-#include <osccoeffImpl/osclocalproblemregion.hh>
-#include <osccoeffImpl/oscl2functionalassembler.hh>
-#include <osccoeffImpl/hierarchicleveliterator.hh>
-
-#include <dune/fufem/assemblers/assembler.hh>
-#include <dune/fufem/assemblers/localassemblers/laplaceassembler.hh>
-#include <dune/fufem/boundarypatch.hh>
-
-#include <dune/grid/common/mcmgmapper.hh>
-#include <dune/grid/common/intersectioniterator.hh>
-
-#include <dune/istl/umfpack.hh>
-
-
-
-
-template <class GridType, class BasisType, class MatrixType, class VectorType>
-class OscPreconditioner : public LinearIterationStep<MatrixType, VectorType> {
-
-public:
- typedef GridType grid_type;
- typedef MatrixType matrix_type;
-
-private:
- static const int dim= GridType::dimension;
- static const int dimworld = GridType::dimensionworld;
- typedef typename GridType::ctype ctype;
- typedef typename Dune::FieldVector<ctype, dimworld> WorldVectorType;
- typedef typename Dune::LevelMultipleCodimMultipleGeomTypeMapper <GridType, Dune::MCMGElementLayout > MapperType;
- typedef typename GridType::template Codim<0>::Entity EntityType;
- typedef typename GridType::template Codim<0>::EntityPointer EntityPointerType;
- typedef typename GridType::LevelGridView GV;
- typedef typename BasisType::LocalFiniteElement FE;
-
- typedef typename GV::template Codim <dim>::Iterator VertexLevelIterator;
- typedef HierarchicLevelIterator<GridType> HierarchicLevelIteratorType;
-
- GridType& grid_;
- std::vector<ctype>& oscData_;
- MapperType& mapper_;
- std::vector<OscLocalProblem<BasisType, MatrixType, VectorType, VectorType>* > localProblems_;
- const int coarseLevel_;
- const int fineLevel_;
- const ctype boxLength_;
-
- VectorType& fineRhs_;
-
- MatrixType coarseStiffMat_;
- VectorType coarseRhs_;
- Dune::BitSetVector<1> coarseBoundaryDofs_;
- CompressedMultigridTransfer<VectorType> mgTransfer_;
-
- GV coarseGridView;
- GV fineGridView;
-
- BasisType coarseBasis;
- BasisType fineBasis;
-
- Dune::LevelMultipleCodimMultipleGeomTypeMapper <GridType, Dune::MCMGElementLayout > coarseMapper;
-
- Dune::BitSetVector<1> fineBoundaryDofs;
-
- struct classcomp {
- bool operator() (const WorldVectorType& lhs, const WorldVectorType rhs) const {
-
- for (size_t i=0; i<lhs.size(); ++i) {
- if (lhs[i]>rhs[i])
- return false;
- else if (lhs[i]<rhs[i])
- return true;
- }
- return false;
- }
- };
-
- std::map<WorldVectorType, EntityPointerType, classcomp> regionSeeds;
-
- void addRegionElements(const EntityPointerType seed, std::vector<EntityPointerType>& regionElements, int& elemCounter){
-
- HierarchicLevelIteratorType endHierIt(grid_, *seed, HierarchicLevelIteratorType::PositionFlag::end, fineLevel_);
- for (HierarchicLevelIteratorType hierIt(grid_, *seed, HierarchicLevelIteratorType::PositionFlag::begin, fineLevel_); hierIt!=endHierIt; ++hierIt) {
-
- EntityPointerType entityPointer(*hierIt);
- regionElements[elemCounter] = entityPointer;
- elemCounter++;
- }
- }
-
- void addLocalProblem(VertexLevelIterator& it /*, int& count*/) {
- OscRhs<WorldVectorType, Dune::FieldVector<ctype,1>> f;
-
- int count =0;
-
- const WorldVectorType vertex = it->geometry().corner(0);
- EntityPointerType regionSeed(regionSeeds.at(vertex));
-
- /* std::cout << vertex[0] << " " << vertex[1] << " " << std::endl;
- std::cout << regionSeed->geometry().center().two_norm() << std::endl;*/
-
- OscLocalProblemBox<EntityType, WorldVectorType> fineLocalProblemRegion(vertex, boxLength_);
-
- std::vector<EntityPointerType> regionElements(mapper_.size(), regionSeed);
- int elemCounter = 0;
- addRegionElements(regionSeed, regionElements, elemCounter);
-
- std::queue<EntityPointerType> elementQueue;
- elementQueue.push(regionSeed);
-
- Dune::BitSetVector<1> visited(coarseMapper.size());
- visited.unsetAll();
-
- visited[coarseMapper.map(*regionSeed)][0] = true;
-
- while (!elementQueue.empty()) {
-
- const EntityPointerType elem = elementQueue.front();
- elementQueue.pop();
-
- typename GridType::LevelIntersectionIterator neighborIt = coarseGridView.ibegin(*elem);
- typename GridType::LevelIntersectionIterator neighborItEnd = coarseGridView.iend(*elem);
- for (; neighborIt!=neighborItEnd; ++neighborIt) {
- const typename GridType::LevelIntersection& intersection = *neighborIt;
-
- if (intersection.neighbor()) {
- EntityType& neighbor = *intersection.outside();
-
- bool inRegion;
- fineLocalProblemRegion.evaluate(neighbor, inRegion);
- int neighborId = coarseMapper.map(neighbor);
-
-
- if (inRegion and !(visited[neighborId][0])) {
- EntityPointerType neighborPointer(neighbor);
- elementQueue.push(neighborPointer);
- visited[neighborId][0] = true;
- addRegionElements(neighborPointer, regionElements, elemCounter);
- }
- }
- }
- }
-
- regionElements.erase(regionElements.begin()+elemCounter, regionElements.end());
- //std::cout << regionElements.size() << " " << elemCounter << std::endl;
-
- OscLocalAssembler<GridType, FE, FE> oscLocalAssembler(oscData_, mapper_, fineLocalProblemRegion);
- OscL2FunctionalAssembler<GridType> oscLocalFunctionalAssembler(f, fineLocalProblemRegion);
- OscLocalProblem<BasisType, MatrixType, VectorType, VectorType>* localProblem = new OscLocalProblem<BasisType, MatrixType, VectorType, VectorType>(regionElements, fineBasis, fineBoundaryDofs, oscLocalAssembler, fineRhs_, oscLocalFunctionalAssembler);
-
- localProblems_[count] = localProblem;
-
- count++;
- }
-
-public:
-
- OscPreconditioner(GridType& grid, int coarseLevel, int fineLevel, std::vector<ctype>& oscData, MapperType& mapper, ctype boxLength, VectorType& fineRhs)
- : grid_(grid), oscData_(oscData), coarseLevel_(coarseLevel), fineLevel_(fineLevel), mapper_(mapper), boxLength_(boxLength), fineRhs_(fineRhs),
- coarseGridView(grid_.levelView(coarseLevel_)), fineGridView(grid_.levelView(fineLevel_)), coarseBasis(coarseGridView), fineBasis(fineGridView),
- coarseMapper(grid_, coarseLevel_)
- {
-
-
- std::cout << "Setting coarse grid oscData:" << std::endl;
-
-
- // allocate a vector for the data
- std::vector<ctype> coarseOscData(coarseMapper.size());
-
- typedef typename GV::template Codim <0>::Iterator ElementLevelIterator;
-
- ElementLevelIterator endElemIt = coarseGridView.template end <0>();
- for (ElementLevelIterator it = coarseGridView. template begin <0>(); it!=endElemIt; ++it) {
- const typename EntityType::Geometry geometry = it->geometry();
-
- // set region seed elements
- for (size_t i=0; i<geometry.corners(); i++) {
- WorldVectorType corner = geometry.corner(i);
- typename std::map<WorldVectorType, EntityPointerType>::iterator hasCorner = regionSeeds.find(corner);
-
- if (hasCorner==regionSeeds.end()) {
- EntityPointerType entityPointer(*it);
- regionSeeds.insert(std::pair<WorldVectorType, EntityPointerType>(corner, entityPointer));
- }
- }
-
-
- // set coarseOscData
- int counter = 0;
- ctype sum = 0;
-
- HierarchicLevelIteratorType endHierIt(grid_, *it, HierarchicLevelIteratorType::PositionFlag::end, fineLevel_);
- for (HierarchicLevelIteratorType hierIt(grid_, *it, HierarchicLevelIteratorType::PositionFlag::begin, fineLevel_); hierIt!=endHierIt; ++hierIt) {
-
- sum += oscData_[mapper_.map(*hierIt)];
- counter++;
- }
- coarseOscData[coarseMapper.map(*it)] = sum/counter;
- // std::cout << "elem: " << (int) mapper.map(*it) << "| "<< coarseOscData[coarseMapper.map(*it)] << std::endl;
- }
-
-
- WorldVectorType zero(2);
- zero = 0;
- typename std::map<WorldVectorType, EntityPointerType>::iterator has = regionSeeds.find(zero);
- std::cout << "region size: "<<regionSeeds.size() << std::endl;
- //std::cout << (has==regionSeeds.end())<< " " << std::endl;
-
- std::cout << "Assembling coarse grid stiffness matrix!" << std::endl;
-
- // assemble coarse stiffness matrix
- WorldVectorType midPoint(1/2);
- ctype box = 1;
- Assembler<BasisType, BasisType> coarseAssembler(coarseBasis, coarseBasis);
- OscLocalProblemBox<typename GridType::template Codim<0>::Entity, WorldVectorType> coarseLocalProblemRegion(midPoint, box);
- OscLocalAssembler<GridType, FE, FE> coarseLocalAssembler(coarseOscData, coarseMapper, coarseLocalProblemRegion);
- coarseAssembler.assembleOperator(coarseLocalAssembler, coarseStiffMat_);
-
-
- std::cout << "Setting coarse grid boundary conditions!" << std::endl;
-
- // set coarse boundary conditions
- BoundaryPatch<typename GridType::LevelGridView> coarseBoundaryPatch(coarseGridView, true);
- constructBoundaryDofs(coarseBoundaryPatch, coarseBasis, coarseBoundaryDofs_);
-
- for(size_t i=0; i<coarseBoundaryDofs_.size(); i++) {
- if(!coarseBoundaryDofs_[i][0])
- continue;
-
- for(size_t j=0; j<coarseStiffMat_.M(); j++) {
- if(coarseStiffMat_.exists(i,j)) {
- coarseStiffMat_[i][j] = 0;
- }
- }
- coarseStiffMat_[i][i]=1;
- }
-
- mgTransfer_.setup(coarseGridView, fineGridView);
-
- localProblems_.resize(grid_.size(coarseLevel_, dim));
-
- std::cout << "Initializing local fine grid corrections! " << std::endl;
-
-
-
- // init local fine level corrections
- int count = 0;
- BoundaryPatch<typename GridType::LevelGridView> fineBoundaryPatch(fineGridView, true);
- constructBoundaryDofs(fineBoundaryPatch, fineBasis, fineBoundaryDofs);
-
-
-
- Timer timer;
- timer.reset();
- timer.start();
- std::list<std::thread*> threads;
- VertexLevelIterator endIt = coarseGridView.template end <dim>();
- for (VertexLevelIterator it = coarseGridView.template begin <dim>(); it!=endIt; ++it) {
- if (count==10) {
- std::cout << "Predicted total setup time: " << timer.elapsed()/10*coarseBasis.size() << " seconds." << std::endl;
- }
- std::thread* thread = new std::thread(&OscPreconditioner::addLocalProblem, this /*, count*/);
- // threads.push_back(thread);
- }
-
- // std::list<std::thread*>::iterator threadEndIt = threads.end();
- // for (std::list<std::thread*>::iterator threadIt = threads.begin(); threadIt!=threadEndIt; ++threadIt){
- // (*threadIt)->join();
- // }
-
-
- std::cout << "Total setup time: " << timer.elapsed() << " seconds." << std::endl;
- timer.stop();
- }
-
-
- virtual void setProblem(const MatrixType& mat, VectorType& x, const VectorType& rhs) {
- this->x_ = &x;
- this->rhs_ = &rhs;
- this->mat_ = Dune::stackobject_to_shared_ptr(mat);
-
- mgTransfer_.restrict(*(this->rhs_), coarseRhs_);
- for(size_t i=0; i<coarseBoundaryDofs_.size(); i++) {
- if(!coarseBoundaryDofs_[i][0])
- continue;
-
- coarseRhs_[i] = 0;
- }
-
- for (size_t i=0; i<localProblems_.size(); i++) {
- localProblems_[i]->updateRhs(*(this->rhs_));
- }
- }
-
- virtual void iterate() {
-
- // compute exact coarse solution
- *(this->x_) = 0;
-
- VectorType coarseSol(coarseStiffMat_.M());
- #if HAVE_UMFPACK
- Dune::InverseOperatorResult res;
-
- Dune::UMFPack<MatrixType> coarseSolver(coarseStiffMat_);
- coarseSolver.apply(coarseSol, coarseRhs_, res);
- #else
- #error No UMFPack!
- #endif
-
- /*std::cout << "coarseStiffMat_: " << std::endl;
- for(size_t rowIdx = 0; rowIdx<localProblems_[6]->stiffnessMat_.N(); rowIdx++) {
- for(size_t colIdx = 0; colIdx<localProblems_[6]->stiffnessMat_.M(); colIdx++) {
- if (localProblems_[6]->stiffnessMat_.exists(rowIdx, colIdx))
- std::cout << localProblems_[6]->stiffnessMat_[rowIdx][colIdx] << " ";
- else
- std::cout << "0 ";
- }
- std::cout << std::endl;
- }
- std::cout << std::endl;*/
-/*
- std::cout << "coarseRhs_: " << std::endl;
- for (size_t i = 0; i<coarseRhs_.size(); i++) {
- std::cout << coarseRhs_[i] << " ";
- }
- std::cout << std::endl;
-
- std::cout << "coarseSol: " << std::endl;
- for (size_t i = 0; i<coarseSol.size(); i++) {
- std::cout << coarseSol[i] << " ";
- }
- std::cout << std::endl;
-*/
- mgTransfer_.prolong(coarseSol, *(this->x_));
-
- VectorType x;
- for (size_t i=0; i<localProblems_.size(); i++) {
-
- localProblems_[i]->solve(x);
- *(this->x_) += x;
- }
- }
-
-};
-
-#endif
-
diff --git a/dune/faultnetworks/preconditioners/test/CMakeLists.txt b/dune/faultnetworks/preconditioners/test/CMakeLists.txt
deleted file mode 100644
index 3a5964adeb30007351d7b9a5c14236191cbe0222..0000000000000000000000000000000000000000
--- a/dune/faultnetworks/preconditioners/test/CMakeLists.txt
+++ /dev/null
@@ -1,12 +0,0 @@
-# Put your test in here if it needs access to external grids
-set(FAULTNETWORKS_PRECONDITIONERS_TESTS
- levelfaultpreconditionertest
- )
-
-set(TESTS ${FAULTNETWORKS_PRECONDITIONERS_TESTS})
-
-# Setup targets, register tests, and add dune flags
-foreach(_test ${TESTS})
- #dune_add_test(SOURCES ${_test}.cc)
-endforeach()
-
diff --git a/dune/faultnetworks/preconditioners/test/Makefile.am b/dune/faultnetworks/preconditioners/test/Makefile.am
deleted file mode 100644
index a4c19c2c84cac95071500f9fa9c79aa2677057d1..0000000000000000000000000000000000000000
--- a/dune/faultnetworks/preconditioners/test/Makefile.am
+++ /dev/null
@@ -1,42 +0,0 @@
-# which tests where program to build and run are equal
-NORMALTESTS = levelfaultpreconditionertest
-
-# list of tests to run
-TESTS = $(NORMALTESTS)
-
-# programs just to build when "make check" is used
-check_PROGRAMS = $(NORMALTESTS)
-
-# define the programs
-levelfaultpreconditionertest_SOURCES = levelfaultpreconditionertest.cc
-
-
-levelfaultpreconditionertest_CPPFLAGS = $(AM_CPPFLAGS) \
- $(DUNEMPICPPFLAGS) \
- $(UG_CPPFLAGS)
-
-# The libraries have to be given in reverse order (most basic libraries
-# last). Also, due to some misunderstanding, a lot of libraries include the
-# -L option in LDFLAGS instead of LIBS -- so we have to include the LDFLAGS
-# here as well.
-levelfaultpreconditionertest_LDADD = \
- $(DUNE_LDFLAGS) $(DUNE_LIBS) \
- $(UG_LDFLAGS) $(UG_LIBS) \
- $(DUNEMPILIBS) \
- $(LDADD)
-levelfaultpreconditionertest_LDFLAGS = $(AM_LDFLAGS) \
- $(DUNEMPILDFLAGS) \
- $(UG_LDFLAGS) \
- $(DUNE_LDFLAGS)
-
-# don't follow the full GNU-standard
-# we need automake 1.9
-AUTOMAKE_OPTIONS = foreign 1.9
-
-# pass most important options when "make distcheck" is used
-DISTCHECK_CONFIGURE_FLAGS = --with-dune-common=$(DUNE_COMMON_ROOT) --with-dune-geometry=$(DUNE_GEOMETRY_ROOT) --with-dune-grid=$(DUNE_GRID_ROOT) --with-dune-istl=$(DUNE_ISTL_ROOT) --with-dune-localfunctions=$(DUNE_LOCALFUNCTIONS_ROOT) --with-dune-solvers=$(DUNE_SOLVERS_ROOT) --with-dune-fufem=$(DUNE_FUFEM_ROOT) CXX="$(CXX)" CC="$(CC)"
-
-EXTRA_DIST = CMakeLists.txt
-
-include $(top_srcdir)/am/global-rules
-
diff --git a/dune/faultnetworks/preconditioners/test/levelfaultpreconditionertest.cc b/dune/faultnetworks/preconditioners/test/levelfaultpreconditionertest.cc
deleted file mode 100644
index 69b2283a34652d187aa57bf7bb0623e0168cb519..0000000000000000000000000000000000000000
--- a/dune/faultnetworks/preconditioners/test/levelfaultpreconditionertest.cc
+++ /dev/null
@@ -1,274 +0,0 @@
-#include <config.h>
-
-#include <dune/common/parallel/mpihelper.hh>
-
-#include <dune/grid/uggrid.hh>
-#include <dune/grid/common/mcmgmapper.hh>
-
-#include <dune/fufem/assemblers/localassemblers/laplaceassembler.hh>
-
-#include <dune/faultnetworks/interfacenetwork.hh>
-#include <dune/faultnetworks/faultfactories/multileveluniformfaultfactory.hh>
-#include <dune/faultnetworks/preconditioners/levelfaultpreconditioner.hh>
-#include <dune/faultnetworks/preconditioners/multilevelfaultpreconditioner.hh>
-#include <dune/faultnetworks/assemblers/intersectionjumpmassassembler.hh>
-#include <dune/faultnetworks/faultp1nodalbasis.hh>
-#include <dune/faultnetworks/utils/debugutils.hh>
-
-const int dim = 2;
-
-const int coarseResolution = 2;
-const int fineResolution = 3;
-const int exactResolution = 4;
-
-const int faultCount = 3;
-
-
-typedef typename Dune::UGGrid<dim> GridType;
-typedef typename GridType::LevelGridView GridView;
-typedef typename GridType::LevelIntersection Intersection;
-typedef typename GridType::template Codim<0>::Entity Element;
-
-static const int dimElement = Element::dimension;
-
-void gridLevelOutput(const GridType& grid, const int level) {
- std::cout << "Grid elements on lvl " << level << std::endl << std::endl;
-
- typedef typename GridView::template Codim<0>::Iterator ElementLevelIterator;
- GridView gridView(grid.levelGridView(level));
-
- //Dune::LevelMultipleCodimMultipleGeomTypeMapper <GridType, Dune::MCMGElementLayout > mapper(grid, level);
- Dune::LevelMultipleCodimMultipleGeomTypeMapper <GridType, Dune::MCMGVertexLayout > vMapper(grid, level);
-
- Dune::BitSetVector<1> vertexPrinted(vMapper.size(), false);
-
- ElementLevelIterator elemIt = gridView.template begin<0>();
- ElementLevelIterator elemEndIt = gridView.template end<0>();
- for (; elemIt!=elemEndIt; ++elemIt) {
- const Element& elem = *elemIt;
- const auto& geometry = elem.geometry();
-
- for(int i=0; i<geometry.corners(); ++i) {
- size_t globalIdx = gridView.indexSet().subIndex(elem, i, Element::dimension);
- if (!vertexPrinted[globalIdx][0]) {
- const auto& vertex = geometry.corner(i);
- print(vertex, "elem corner:" + std::to_string(globalIdx));
- vertexPrinted[globalIdx][0] = true;
- }
- }
- }
- std::cout << "----------------" << std::endl << std::endl;
-
- /*
- ElementLevelIterator elemIt = gridView.template begin<0>();
- ElementLevelIterator elemEndIt = gridView.template end<0>();
- for (; elemIt!=elemEndIt; ++elemIt) {
- const Element& elem = *elemIt;
- const auto& geometry = elem.geometry();
-
- std::cout << "Element idx: " << mapper.index(elem) << std::endl;
-
- for(int i=0; i<geometry.corners(); ++i) {
- size_t globalIdx = gridView.indexSet().subIndex(elem, i, Element::dimension);
- const auto& vertex = geometry.corner(i);
- print(vertex, "elem corner:" + std::to_string(globalIdx));
- }
- std::cout << "----------------" << std::endl << std::endl;
- }*/
-
-}
-
-
-void faultInterfaceOutput(const FaultInterface<GridView>& faultInterface) {
- std::cout << "Level: " << faultInterface.level() << " global faultID: " << faultInterface.id() << std::endl;
-
- const std::set<size_t>& interfaceDofs = faultInterface.getInterfaceDofs();
- print(interfaceDofs, "");
-}
-
-void interfaceNetworkLevelOutput(const InterfaceNetwork<GridType>& interfaceNetwork, const int level) {
- const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork = interfaceNetwork.levelInterfaceNetwork(level);
-
- std::cout << "LevelInterfaceNetwork level: " << level << std::endl;
- for (size_t i=0; i<levelInterfaceNetwork.size(); i++) {
- std::cout << "local faultID : " << i << " ";
- faultInterfaceOutput(*(levelInterfaceNetwork.getInterface(i)));
- }
-
- std::cout << "------------------------------" << std::endl << std::endl;
-}
-
-void computeIntersectionDofs(const GridView& gridView, const Intersection& intersection, std::set<size_t>& intersectionDofs) {
- intersectionDofs.clear();
-
- // loop over all vertices of the intersection
- const Element& insideElement = intersection.inside();
-
- const auto& refElement = Dune::ReferenceElements<double,dimElement>::general(insideElement.type());
- for (int i=0; i<refElement.size(intersection.indexInInside(), 1, dimElement); i++) {
- size_t idxInElement = refElement.subEntity(intersection.indexInInside(), 1, i, dimElement);
- size_t globalIdx = gridView.indexSet().subIndex(insideElement, idxInElement, dimElement);
- intersectionDofs.insert(globalIdx);
- }
-}
-
-void prolongFaultDofs(const GridType& grid, const FaultInterface<GridView>& faultInterface, const int toLevel, std::set<size_t>& prolongedDofs) {
- prolongedDofs.clear();
-
- const std::vector<Intersection>& faces = faultInterface.faces();
- const GridView gridView(grid.levelGridView(toLevel));
-
- for (size_t i=0; i<faces.size(); i++) {
- const Intersection& intersection = faces[i];
- Face<GridType> face(grid, intersection.inside(), intersection.indexInInside());
-
- typename Face<GridType>::HierarchicIterator exactIt = face.hbegin(toLevel);
- typename Face<GridType>::HierarchicIterator exactEnd = face.hend(toLevel);
- for(; exactIt!=exactEnd; ++exactIt) {
- if (exactIt->level() == toLevel) {
- std::set<size_t> faceDofs;
- computeIntersectionDofs(gridView, *(exactIt->intersection()), faceDofs);
-
- prolongedDofs.insert(faceDofs.begin(), faceDofs.end());
- }
- }
- }
-}
-
-void prolongInterfaceNetworkLevelDofs(const GridType& grid, const InterfaceNetwork<GridType>& interfaceNetwork, const int level, const int toLevel, std::vector<std::set<size_t>>& prolongedFaultDofs) {
- prolongedFaultDofs.resize(0);
-
- const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork = interfaceNetwork.levelInterfaceNetwork(level);
-
- std::cout << "Correct prolonged LevelInterfaceNetwork level " << level << " to level " << toLevel << std::endl;
- for (size_t i=0; i<levelInterfaceNetwork.size(); i++) {
- std::set<size_t> prolongedDofs;
- prolongFaultDofs(grid, *(levelInterfaceNetwork.getInterface(i)), toLevel, prolongedDofs);
-
- print(prolongedDofs, "local faultID: " + std::to_string(i));
- prolongedFaultDofs.push_back(prolongedDofs);
- }
-
- std::cout << "------------------------------" << std::endl << std::endl;
-}
-
-int main(int argc, char** argv) try
-{
- typedef double ctype;
-
- Dune::MPIHelper::instance(argc, argv);
-
- std::ofstream out("/home/mi/podlesny/software/dune/faultnetworks/dune/faultnetworks/preconditioners/test/levelfaultpreconditionertest.log");
- std::streambuf *coutbuf = std::cout.rdbuf(); //save old buf
- std::cout.rdbuf(out.rdbuf()); //redirect std::cout to out.txt!
-
- const int fineLevelIdx = fineResolution - coarseResolution;
- const int exactLevelIdx = exactResolution - coarseResolution;
-
- std::vector<int> faultToLevel(faultCount);
- for (size_t i=0; i<faultToLevel.size(); ) {
- for (int j=fineLevelIdx; j>0 && i<faultToLevel.size(); j--) {
- faultToLevel[i] = j;
- i++;
- }
-
- for (int j=0; j<fineLevelIdx && i<faultToLevel.size(); j++) {
- faultToLevel[i] = j;
- i++;
- }
- }
-
- print(faultToLevel, "faultToLevel: ");
-
- MultilevelUniformFaultFactory<GridType> faultFactory(faultCount, faultToLevel, coarseResolution, exactLevelIdx);
- const InterfaceNetwork<GridType>& interfaceNetwork = faultFactory.interfaceNetwork();
-
- /*
- const GridType& grid = faultFactory.grid();
-
- gridLevelOutput(grid, 0);
- gridLevelOutput(grid, fineLevelIdx);
- gridLevelOutput(grid, exactLevelIdx);
-
- interfaceNetworkLevelOutput(interfaceNetwork, 0);
- interfaceNetworkLevelOutput(interfaceNetwork, fineLevelIdx);
- interfaceNetworkLevelOutput(interfaceNetwork, exactLevelIdx); */
-
- typedef typename Dune::BlockVector<Dune::FieldVector<ctype,1>> VectorType;
- typedef typename Dune::BCRSMatrix<Dune::FieldMatrix<ctype, 1, 1>> MatrixType;
-
- typedef FaultP1NodalBasis<GridView, ctype> DGBasis;
- typedef typename DGBasis::LocalFiniteElement DGFE;
-
- typedef IntersectionJumpMassAssembler<GridView, DGFE, DGFE> LocalInterfaceAssembler;
- typedef Dune::Matrix<Dune::FieldMatrix<double,1,1>, std::allocator<Dune::FieldMatrix<double,1,1>>> OscMatrixType;
- OscMatrixType B(2,2);
- B[0][0] = 1;
- B[1][1] = 1;
- B[0][1] = -1;
- B[1][0] = -1;
-
- //--------------------------------------------------------------
- // Multilevelfaultpreconditioner
- //--------------------------------------------------------------
-
- // setting local assemblers for laplace operator
- LocalInterfaceAssembler localInterfaceAssembler(B);
- typedef LaplaceAssembler<GridType, DGFE, DGFE> LocalLaplaceAssembler;
- LocalLaplaceAssembler localAssembler;
-
- Dune::BitSetVector<1> activeLevels(fineLevelIdx+1, true);
- print(activeLevels, "activeLevels: ");
-
- std::vector<std::shared_ptr<LocalLaplaceAssembler>> localAssemblers;
- std::vector<std::shared_ptr<LocalInterfaceAssembler>> localIntersectionAssemblers;
-
- localAssemblers.resize(activeLevels.size());
- localIntersectionAssemblers.resize(activeLevels.size());
-
- for (size_t i=0; i<activeLevels.size(); i++) {
- localAssemblers[i] = std::make_shared<LocalLaplaceAssembler>(localAssembler);
- localIntersectionAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(localInterfaceAssembler);
- }
-
- MultilevelFaultPreconditioner<DGBasis, LocalLaplaceAssembler, LocalInterfaceAssembler, MatrixType, VectorType> preconditioner(interfaceNetwork, activeLevels, localAssemblers, localIntersectionAssemblers);
-
- MatrixType matrix;
- VectorType mgX(preconditioner.basis().size());
- VectorType mgRhs(preconditioner.basis().size());
- mgRhs = 0;
-
- preconditioner.setProblem(matrix, mgX, mgRhs);
- preconditioner.iterate();
-
- const VectorType& mgSol = preconditioner.getSol();
- print(mgSol, "mgSol: ");
-
- step(0);
-
- const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork = interfaceNetwork.levelInterfaceNetwork(1);
- LevelFaultPreconditioner<DGBasis, LocalLaplaceAssembler, LocalInterfaceAssembler, MatrixType, VectorType> levelFaultPreconditioner(levelInterfaceNetwork, localAssembler, localInterfaceAssembler);
-
- step(1);
-
- VectorType x(levelFaultPreconditioner.basis().size());
- VectorType rhs(levelFaultPreconditioner.basis().size());
- rhs = 0;
-
- step(2);
-
-
- levelFaultPreconditioner.setProblem(matrix, x, rhs);
-
- step(3);
-
- levelFaultPreconditioner.iterate();
-
- const VectorType& sol = levelFaultPreconditioner.getSol();
-
- print(sol, "sol: ");
-
- std::cout.rdbuf(coutbuf);
-} catch (Dune::Exception e){
- std::cout << e << std::endl;
-}
diff --git a/dune/faultnetworks/preconditioners/test/levelfaultpreconditionertest.log b/dune/faultnetworks/preconditioners/test/levelfaultpreconditionertest.log
deleted file mode 100644
index a66d14bfdf5e4aed77911c4071d9eda95e757075..0000000000000000000000000000000000000000
--- a/dune/faultnetworks/preconditioners/test/levelfaultpreconditionertest.log
+++ /dev/null
@@ -1,26 +0,0 @@
-faultToLevel:
-1 0 1
-
-Coarse and fine grid were generated!
-activeLevels:
-1 1
-
-Initializing local fine grid corrections!
-
-Total setup time: 0.000599407 seconds.
-Initializing local fine grid corrections!
-
-Total setup time: 0.00250133 seconds.
-mgSol:
-0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
-
- Step 0!
-Initializing local fine grid corrections!
-
-Total setup time: 0.00260553 seconds.
- Step 1!
- Step 2!
- Step 3!
-sol:
-0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
-
diff --git a/dune/faultnetworks/utils/exactvertexdata.vtu b/dune/faultnetworks/utils/exactvertexdata.vtu
deleted file mode 100644
index 55ffbafe9382f0ae7d3e8e18909004dec6eb0a23..0000000000000000000000000000000000000000
Binary files a/dune/faultnetworks/utils/exactvertexdata.vtu and /dev/null differ
diff --git a/dune/faultnetworks/utils/finevertexdata.vtu b/dune/faultnetworks/utils/finevertexdata.vtu
deleted file mode 100644
index d76c79a48e53887f500cb29231c76c3a8ba7eebe..0000000000000000000000000000000000000000
Binary files a/dune/faultnetworks/utils/finevertexdata.vtu and /dev/null differ