diff --git a/configure.ac b/configure.ac index ab82a6594b7e599c78acea477f11b601d22fb7f2..48a891eff1ab4bd554ea5b1299fb2ff2e3c24c03 100644 --- a/configure.ac +++ b/configure.ac @@ -79,6 +79,8 @@ AC_CONFIG_FILES([dumux.pc dumux/material/fluidsystems/Makefile dumux/material/constraintsolvers/Makefile dumux/material/eos/Makefile + dumux/multidomain/Makefile + dumux/multidomain/common/Makefile dumux/nonlinear/Makefile dumux/parallel/Makefile m4/Makefile diff --git a/dumux/Makefile.am b/dumux/Makefile.am index bc835ce36fe34be06d1404ed62b5bbf2e60e0d39..640d0c72079152851d3f87e7662c01ed4ade1cc4 100644 --- a/dumux/Makefile.am +++ b/dumux/Makefile.am @@ -7,6 +7,7 @@ SUBDIRS = \ io \ linear \ material \ + multidomain \ nonlinear \ parallel diff --git a/dumux/multidomain/Makefile.am b/dumux/multidomain/Makefile.am new file mode 100644 index 0000000000000000000000000000000000000000..f1de35a999ec8752a2b014375be2c176d35212e7 --- /dev/null +++ b/dumux/multidomain/Makefile.am @@ -0,0 +1,4 @@ +SUBDIRS = common +multidomaindir = $(includedir)/dumux/multidomain + +include $(top_srcdir)/am/global-rules diff --git a/dumux/multidomain/common/Makefile.am b/dumux/multidomain/common/Makefile.am new file mode 100644 index 0000000000000000000000000000000000000000..1c56fe2ecadbe9cc1c8a185370bcdcfbb131c523 --- /dev/null +++ b/dumux/multidomain/common/Makefile.am @@ -0,0 +1,4 @@ +commondir = $(includedir)/dumux/multidomain/common +common_HEADERS = *.hh + +include $(top_srcdir)/am/global-rules diff --git a/dumux/multidomain/common/splitandmerge.hh b/dumux/multidomain/common/splitandmerge.hh new file mode 100644 index 0000000000000000000000000000000000000000..80529da88be142cd4415109356ffc884e6af3886 --- /dev/null +++ b/dumux/multidomain/common/splitandmerge.hh @@ -0,0 +1,254 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see <http://www.gnu.org/licenses/>. * + *****************************************************************************/ +/** + * \file + * \brief Some methods required by several classes of the coupling model. + */ +#ifndef DUMUX_SPLIT_AND_MERGE_HH +#define DUMUX_SPLIT_AND_MERGE_HH + +#include "coupledproperties.hh" +#include <dumux/common/valgrind.hh> + +namespace Dumux +{ +/*! + * \addtogroup ModelCoupling + */ +// \{ + +/*! + * \brief Some methods required by several classes of the coupling model. + */ +template<class TypeTag> +class SplitAndMerge +{ + typedef typename GET_PROP_TYPE(TypeTag, SubProblem1TypeTag) SubTypeTag1; + typedef typename GET_PROP_TYPE(TypeTag, SubProblem2TypeTag) SubTypeTag2; + + typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector; + typedef typename GET_PROP_TYPE(SubTypeTag1, SolutionVector) SolutionVector1; + typedef typename GET_PROP_TYPE(SubTypeTag2, SolutionVector) SolutionVector2; + + typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) JacobianMatrix; + typedef typename GET_PROP_TYPE(SubTypeTag1, JacobianMatrix) JacobianMatrix1; + typedef typename GET_PROP_TYPE(SubTypeTag2, JacobianMatrix) JacobianMatrix2; + + typedef typename GET_PROP_TYPE(TypeTag, Model) Model; + typedef typename GET_PROP_TYPE(SubTypeTag1, Model) Model1; + typedef typename GET_PROP_TYPE(SubTypeTag2, Model) Model2; + + enum { + numEq1 = GET_PROP_VALUE(SubTypeTag1, NumEq), + numEq2 = GET_PROP_VALUE(SubTypeTag2, NumEq) + }; +public: + /*! + * \brief Merge two solution vectors of the sub models into a + * global vector: nonoverlapping case. + * + * \param vec1 docme + * \param vec2 docme + * \param dest docme + * + */ + static void mergeSolVectors(const SolutionVector1 &vec1, + const SolutionVector2 &vec2, + SolutionVector &dest) + { +// printvector(std::cout, vec1, "vec1", "row", 200, 1, 3); +// printvector(std::cout, vec2, "vec2", "row", 200, 1, 3); + + int nDofs1 = vec1.size(); + int nDofs2 = vec2.size(); + + // merge the global vector into two individual ones + for (int i = 0; i < nDofs1; ++i) + for (int j = 0; j < numEq1; j++) + { + dest[i*numEq1 + j][0] = vec1[i][j]; + + Valgrind::CheckDefined(dest[i*numEq1 + j][0]); + } + + for (int i = 0; i < nDofs2; ++i) + for (int j = 0; j < numEq2; j++) + { + dest[nDofs1*numEq1 + i*numEq2 + j][0] = vec2[i][j]; + + Valgrind::CheckDefined(dest[nDofs1*numEq1 + i*numEq2 + j][0]); + } + +// printvector(std::cout, dest, "dest", "row", 200, 1, 3); + + } + + /*! + * \brief Split a global solution vector into two solution vectors + * of the sub models: nonoverlapping case. + * + * \param vec docme + * \param dest1 docme + * \param dest2 docme + * + */ + static void splitSolVector(const SolutionVector &vec, + SolutionVector1 &dest1, + SolutionVector2 &dest2) + { +// printvector(std::cout, vec, "vec", "row", 200, 1, 3); + + int nDofs1 = dest1.size(); + int nDofs2 = dest2.size(); + + for (int i = 0; i < nDofs1; ++i) + for (int j = 0; j < numEq1; j++) + dest1[i][j] = vec[i*numEq1 + j][0]; + + for (int i = 0; i < nDofs2; ++i) + for (int j = 0; j < numEq2; j++) + dest2[i][j] = vec[nDofs1*numEq1 + i*numEq2 + j][0]; + +// printvector(std::cout, dest1, "dest1", "row", 200, 1, 3); +// printvector(std::cout, dest2, "dest2", "row", 200, 1, 3); + + } + + /*! + * \brief Merge two solution vectors of the sub models into a + * global vector: more general case. + * + * \param vec1 docme + * \param vec2 docme + * \param dest docme + * \param subDOFToCoupledDOF docme + * + */ + static void mergeSolVectors(const SolutionVector1 &vec1, + const SolutionVector2 &vec2, + SolutionVector &dest, + const std::vector<int>& subDOFToCoupledDOF) + { + int nDofs1 = vec1.size(); + int nDofs2 = vec2.size(); +// printvector(std::cout, vec1, "vec1", "row", 200, 1, 3); +// printvector(std::cout, vec2, "vec2", "row", 200, 1, 3); + + for (int i = 0; i < nDofs1; ++i) + for (int j = 0; j < numEq1; j++) + dest[i*numEq1 + j] = vec1[i][j]; + + for (int i = 0; i < nDofs2; ++i) + { + for (int j = numEq1; j < numEq2; j++) + dest[nDofs1*numEq1 + i*(numEq2-numEq1) + (j - numEq1)] = vec2[i][j]; + } +// printvector(std::cout, dest, "dest", "row", 200, 1, 3); + } + + /*! + * \brief Split a global solution vector into two solution vectors + * of the sub models: more general case. + * + * \param vec docme + * \param dest1 docme + * \param dest2 docme + * \param subDOFToCoupledDOF docme + * + */ + static void splitSolVector(const SolutionVector &vec, + SolutionVector1 &dest1, + SolutionVector2 &dest2, + const std::vector<int>& subDOFToCoupledDOF) + { + int nDofs1 = dest1.size(); + int nDofs2 = dest2.size(); + +// printvector(std::cout, vec, "vec", "row", 200, 1, 3); + for (int i = 0; i < nDofs1; ++i) + for (int j = 0; j < numEq1; j++) + dest1[i][j] = vec[i*numEq1 + j]; + + for (int i = 0; i < nDofs2; ++i) + { + int blockIdxCoupled = subDOFToCoupledDOF[i]; + for (int j = 0; j < numEq1; j++) + { + dest2[i][j] = vec[blockIdxCoupled*numEq1 + j]; + } + + for (int j = numEq1; j < numEq2; j++) + dest2[i][j] = vec[nDofs1*numEq1 + i*(numEq2-numEq1) + (j - numEq1)]; + } +// printvector(std::cout, dest1, "dest1", "row", 200, 1, 3); +// printvector(std::cout, dest2, "dest2", "row", 200, 1, 3); + + } + + + /*! + * \brief Merge individual jacobian matrices of the sub models + * into a global jacobian matrix. + * + * \param M1 docme + * \param M2 docme + * \param M docme + * + */ + static void mergeMatrices(const JacobianMatrix1 &M1, + const JacobianMatrix2 &M2, + JacobianMatrix &M) + { + DUNE_THROW(Dune::NotImplemented, "mergeMatrices in coupled common"); + } + + /*! + * \brief Copy a sub matrix into into the main diagonal of the global matrix. + * + * \param Msub docme + * \param M docme + * \param offset docme? + * + */ + template <class SubMatrix> + static void copyMatrix(const SubMatrix &Msub, + JacobianMatrix &M, + size_t offset) + { + // loop over all rows of the submatrix + typedef typename SubMatrix::ConstRowIterator RowIterator; + typedef typename SubMatrix::ConstColIterator ColIterator; + RowIterator endRow = Msub.end(); + for (RowIterator row = Msub.begin(); row != endRow; ++row) { + // loop over columns of the current row + ColIterator endCol = row->end(); + for (ColIterator col = row->begin(); col != endCol; ++ col) { + // copy entry in the global matrix + M[row.index() + offset][col.index() + offset] + = *col; + } + } + } + +}; +// \} + +} // namespace Dumux + +#endif