From 6bd24612173df3ced4ec05d5836f5a975cee4bd5 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Tue, 20 Nov 2012 09:48:16 +0100 Subject: [PATCH 01/17] added: utilities to find complex eigenvalues of dyn and field matrices this is missing in dune. process has been started to get this integrated upstream http://www.dune-project.org/flyspray/index.php?do=details&task_id=1186&project=1 --- dune/elasticity/Makefile.am | 1 + dune/elasticity/dynmatrixev.cpp | 152 +++++++++++++++++++++++++++++++ dune/elasticity/dynmatrixev.hh | 67 ++++++++++++++ dune/elasticity/fmatrixev_ext.cc | 152 +++++++++++++++++++++++++++++++ dune/elasticity/fmatrixev_ext.hh | 70 ++++++++++++++ 5 files changed, 442 insertions(+) create mode 100644 dune/elasticity/dynmatrixev.cpp create mode 100644 dune/elasticity/dynmatrixev.hh create mode 100644 dune/elasticity/fmatrixev_ext.cc create mode 100644 dune/elasticity/fmatrixev_ext.hh diff --git a/dune/elasticity/Makefile.am b/dune/elasticity/Makefile.am index 19010a7..70e2d68 100644 --- a/dune/elasticity/Makefile.am +++ b/dune/elasticity/Makefile.am @@ -2,6 +2,7 @@ lib_LTLIBRARIES = libelasticityupscale.la libelasticityupscale_la_SOURCES = \ boundarygrid.cpp \ + dynmatrixev.cpp \ material.cpp \ materials.cpp \ mpc.cpp diff --git a/dune/elasticity/dynmatrixev.cpp b/dune/elasticity/dynmatrixev.cpp new file mode 100644 index 0000000..d5880a9 --- /dev/null +++ b/dune/elasticity/dynmatrixev.cpp @@ -0,0 +1,152 @@ +#ifndef DUNE_FMATRIXEIGENVALUES_EXT_CC +#define DUNE_FMATRIXEIGENVALUES_EXT_CC + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include +#include +#include + +#include + +#if HAVE_LAPACK + +// nonsymetric matrices +#define DGEEV_FORTRAN FC_FUNC (dgeev, DGEEV) + +// dsyev declaration (in liblapack) +extern "C" { + +/* + * + ** purpose + ** ======= + ** + ** xgeev computes for an N-by-N BASE DATA TYPE nonsymmetric matrix A, the + ** eigenvalues and, optionally, the left and/or right eigenvectors. + ** + ** The right eigenvector v(j) of A satisfies + ** A * v(j) = lambda(j) * v(j) + ** where lambda(j) is its eigenvalue. + ** The left eigenvector u(j) of A satisfies + ** u(j)**T * A = lambda(j) * u(j)**T + ** where u(j)**T denotes the transpose of u(j). + ** + ** The computed eigenvectors are normalized to have Euclidean norm + ** equal to 1 and largest component real. + ** + ** arguments + ** ========= + ** + ** jobvl (input) char + ** = 'n': left eigenvectors of a are not computed; + ** = 'v': left eigenvectors of a are computed. + ** + ** jobvr (input) char + ** = 'n': right eigenvectors of a are not computed; + ** = 'v': right eigenvectors of a are computed. + ** + ** n (input) long int + ** the order of the matrix v. v >= 0. + ** + ** a (input/output) BASE DATA TYPE array, dimension (lda,n) + ** on entry, the n-by-n matrix a. + ** on exit, a has been overwritten. + ** + ** lda (input) long int + ** the leading dimension of the array a. lda >= max(1,n). + ** + ** wr (output) BASE DATA TYPE array, dimension (n) + ** wi (output) BASE DATA TYPE array, dimension (n) + ** wr and wi contain the real and imaginary parts, + ** respectively, of the computed eigenvalues. complex + ** conjugate pairs of eigenvalues appear consecutively + ** with the eigenvalue having the positive imaginary part + ** first. + ** + ** vl (output) COMPLEX DATA TYPE array, dimension (ldvl,n) + ** if jobvl = 'v', the left eigenvectors u(j) are stored one + ** after another in the columns of vl, in the same order + ** as their eigenvalues. + ** if jobvl = 'n', vl is not referenced. + ** if the j-th eigenvalue is real, then u(j) = vl(:,j), + ** the j-th column of vl. + ** if the j-th and (j+1)-st eigenvalues form a complex + ** conjugate pair, then u(j) = vl(:,j) + i*vl(:,j+1) and + ** u(j+1) = vl(:,j) - i*vl(:,j+1). + ** + ** ldvl (input) long int + ** the leading dimension of the array vl. ldvl >= 1; if + ** jobvl = 'v', ldvl >= n. + ** + ** vr (output) COMPLEX DATA TYPE array, dimension (ldvr,n) + ** if jobvr = 'v', the right eigenvectors v(j) are stored one + ** after another in the columns of vr, in the same order + ** as their eigenvalues. + ** if jobvr = 'n', vr is not referenced. + ** if the j-th eigenvalue is real, then v(j) = vr(:,j), + ** the j-th column of vr. + ** if the j-th and (j+1)-st eigenvalues form a complex + ** conjugate pair, then v(j) = vr(:,j) + i*vr(:,j+1) and + ** v(j+1) = vr(:,j) - i*vr(:,j+1). + ** + ** ldvr (input) long int + ** the leading dimension of the array vr. ldvr >= 1; if + ** jobvr = 'v', ldvr >= n. + ** + ** work (workspace/output) BASE DATA TYPE array, dimension (max(1,lwork)) + ** on exit, if info = 0, work(1) returns the optimal lwork. + ** + ** lwork (input) long int + ** the dimension of the array work. lwork >= max(1,3*n), and + ** if jobvl = 'v' or jobvr = 'v', lwork >= 4*n. for good + ** performance, lwork must generally be larger. + ** + ** if lwork = -1, then a workspace query is assumed; the routine + ** only calculates the optimal size of the work array, returns + ** this value as the first entry of the work array, and no error + ** message related to lwork is issued by xerbla. + ** + ** info (output) long int + ** = 0: successful exit + ** < 0: if info = -i, the i-th argument had an illegal value. + ** > 0: if info = i, the qr algorithm failed to compute all the + ** eigenvalues, and no eigenvectors have been computed; + ** elements i+1:n of wr and wi contain eigenvalues which + ** have converged. + ** +**/ + +extern void DGEEV_FORTRAN(const char* jobvl, const char* jobvr, const long + int* n, double* a, const long int* lda, double* wr, double* wi, double* vl, + const long int* ldvl, double* vr, const long int* ldvr, double* work, + const long int* lwork, const long int* info); + +} // end extern C +#endif + +namespace Dune { + +namespace DynamicMatrixHelp { + +void eigenValuesNonsymLapackCall( + const char* jobvl, const char* jobvr, const long + int* n, double* a, const long int* lda, double* wr, double* wi, double* vl, + const long int* ldvl, double* vr, const long int* ldvr, double* work, + const long int* lwork, const long int* info) +{ +#if HAVE_LAPACK + // call LAPACK dgeev + DGEEV_FORTRAN(jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr, ldvr, + work, lwork, info); +#else + DUNE_THROW(NotImplemented,"eigenValuesNonsymLapackCall: LAPACK not found!"); +#endif +} + +} // end namespace FMatrixHelp + +} // end namespace Dune +#endif diff --git a/dune/elasticity/dynmatrixev.hh b/dune/elasticity/dynmatrixev.hh new file mode 100644 index 0000000..570c5cc --- /dev/null +++ b/dune/elasticity/dynmatrixev.hh @@ -0,0 +1,67 @@ +#ifndef DUNE_DYNMATRIXEIGENVALUES_HH +#define DUNE_DYNMATRIXEIGENVALUES_HH + +#include + +namespace Dune { + +namespace DynamicMatrixHelp { + +// defined in fmatrixev_ext.cpp +extern void eigenValuesNonsymLapackCall( + const char* jobvl, const char* jobvr, const long + int* n, double* a, const long int* lda, double* wr, double* wi, double* vl, + const long int* ldvl, double* vr, const long int* ldvr, double* work, + const long int* lwork, const long int* info); + +template +static void eigenValuesNonSym(const DynamicMatrix& matrix, + DynamicVector& eigenValues) +{ + { + const long int N = matrix.rows(); + const char jobvl = 'n'; + const char jobvr = 'n'; + + // matrix to put into dgeev + double matrixVector[N * N]; + + // copy matrix + int row = 0; + for(int i=0; i(eigenR[i], eigenI[i]); + } +} + +} + +} + +#endif diff --git a/dune/elasticity/fmatrixev_ext.cc b/dune/elasticity/fmatrixev_ext.cc new file mode 100644 index 0000000..a80c060 --- /dev/null +++ b/dune/elasticity/fmatrixev_ext.cc @@ -0,0 +1,152 @@ +#ifndef DUNE_FMATRIXEIGENVALUES_EXT_CC +#define DUNE_FMATRIXEIGENVALUES_EXT_CC + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include +#include +#include + +#include + +#if HAVE_LAPACK + +// nonsymetric matrices +#define DGEEV_FORTRAN FC_FUNC (dgeev, DGEEV) + +// dsyev declaration (in liblapack) +extern "C" { + +/* + * + ** purpose + ** ======= + ** + ** xgeev computes for an N-by-N BASE DATA TYPE nonsymmetric matrix A, the + ** eigenvalues and, optionally, the left and/or right eigenvectors. + ** + ** The right eigenvector v(j) of A satisfies + ** A * v(j) = lambda(j) * v(j) + ** where lambda(j) is its eigenvalue. + ** The left eigenvector u(j) of A satisfies + ** u(j)**T * A = lambda(j) * u(j)**T + ** where u(j)**T denotes the transpose of u(j). + ** + ** The computed eigenvectors are normalized to have Euclidean norm + ** equal to 1 and largest component real. + ** + ** arguments + ** ========= + ** + ** jobvl (input) char + ** = 'n': left eigenvectors of a are not computed; + ** = 'v': left eigenvectors of a are computed. + ** + ** jobvr (input) char + ** = 'n': right eigenvectors of a are not computed; + ** = 'v': right eigenvectors of a are computed. + ** + ** n (input) long int + ** the order of the matrix v. v >= 0. + ** + ** a (input/output) BASE DATA TYPE array, dimension (lda,n) + ** on entry, the n-by-n matrix a. + ** on exit, a has been overwritten. + ** + ** lda (input) long int + ** the leading dimension of the array a. lda >= max(1,n). + ** + ** wr (output) BASE DATA TYPE array, dimension (n) + ** wi (output) BASE DATA TYPE array, dimension (n) + ** wr and wi contain the real and imaginary parts, + ** respectively, of the computed eigenvalues. complex + ** conjugate pairs of eigenvalues appear consecutively + ** with the eigenvalue having the positive imaginary part + ** first. + ** + ** vl (output) COMPLEX DATA TYPE array, dimension (ldvl,n) + ** if jobvl = 'v', the left eigenvectors u(j) are stored one + ** after another in the columns of vl, in the same order + ** as their eigenvalues. + ** if jobvl = 'n', vl is not referenced. + ** if the j-th eigenvalue is real, then u(j) = vl(:,j), + ** the j-th column of vl. + ** if the j-th and (j+1)-st eigenvalues form a complex + ** conjugate pair, then u(j) = vl(:,j) + i*vl(:,j+1) and + ** u(j+1) = vl(:,j) - i*vl(:,j+1). + ** + ** ldvl (input) long int + ** the leading dimension of the array vl. ldvl >= 1; if + ** jobvl = 'v', ldvl >= n. + ** + ** vr (output) COMPLEX DATA TYPE array, dimension (ldvr,n) + ** if jobvr = 'v', the right eigenvectors v(j) are stored one + ** after another in the columns of vr, in the same order + ** as their eigenvalues. + ** if jobvr = 'n', vr is not referenced. + ** if the j-th eigenvalue is real, then v(j) = vr(:,j), + ** the j-th column of vr. + ** if the j-th and (j+1)-st eigenvalues form a complex + ** conjugate pair, then v(j) = vr(:,j) + i*vr(:,j+1) and + ** v(j+1) = vr(:,j) - i*vr(:,j+1). + ** + ** ldvr (input) long int + ** the leading dimension of the array vr. ldvr >= 1; if + ** jobvr = 'v', ldvr >= n. + ** + ** work (workspace/output) BASE DATA TYPE array, dimension (max(1,lwork)) + ** on exit, if info = 0, work(1) returns the optimal lwork. + ** + ** lwork (input) long int + ** the dimension of the array work. lwork >= max(1,3*n), and + ** if jobvl = 'v' or jobvr = 'v', lwork >= 4*n. for good + ** performance, lwork must generally be larger. + ** + ** if lwork = -1, then a workspace query is assumed; the routine + ** only calculates the optimal size of the work array, returns + ** this value as the first entry of the work array, and no error + ** message related to lwork is issued by xerbla. + ** + ** info (output) long int + ** = 0: successful exit + ** < 0: if info = -i, the i-th argument had an illegal value. + ** > 0: if info = i, the qr algorithm failed to compute all the + ** eigenvalues, and no eigenvectors have been computed; + ** elements i+1:n of wr and wi contain eigenvalues which + ** have converged. + ** +**/ + +extern void DGEEV_FORTRAN(const char* jobvl, const char* jobvr, const long + int* n, double* a, const long int* lda, double* wr, double* wi, double* vl, + const long int* ldvl, double* vr, const long int* ldvr, double* work, + const long int* lwork, const long int* info); + +} // end extern C +#endif + +namespace Dune { + +namespace FMatrixHelp { + +void eigenValuesNonsymLapackCall( + const char* jobvl, const char* jobvr, const long + int* n, double* a, const long int* lda, double* wr, double* wi, double* vl, + const long int* ldvl, double* vr, const long int* ldvr, double* work, + const long int* lwork, const long int* info) +{ +#if HAVE_LAPACK + // call LAPACK dgeev + DGEEV_FORTRAN(jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr, ldvr, + work, lwork, info); +#else + DUNE_THROW(NotImplemented,"eigenValuesNonsymLapackCall: LAPACK not found!"); +#endif +} + +} // end namespace FMatrixHelp + +} // end namespace Dune +#endif diff --git a/dune/elasticity/fmatrixev_ext.hh b/dune/elasticity/fmatrixev_ext.hh new file mode 100644 index 0000000..bfdebcc --- /dev/null +++ b/dune/elasticity/fmatrixev_ext.hh @@ -0,0 +1,70 @@ +#ifndef DUNE_FMATRIXEIGENVALUES_EXT_HH +#define DUNE_FMATRIXEIGENVALUES_EXT_HH + +namespace Dune { + +namespace FMatrixHelp { + +// defined in fmatrixev_ext.cpp +extern void eigenValuesNonsymLapackCall( + const char* jobvl, const char* jobvr, const long + int* n, double* a, const long int* lda, double* wr, double* wi, double* vl, + const long int* ldvl, double* vr, const long int* ldvr, double* work, + const long int* lwork, const long int* info); + +template +static void eigenValuesNonSym(const FieldMatrix& matrix, + FieldVector& eigenValues) +{ + { + const long int N = dim ; + const char jobvl = 'n'; + const char jobvr = 'n'; + + // length of matrix vector + const long int w = N * N ; + + // matrix to put into dgeev + double matrixVector[dim * dim]; + + // copy matrix + int row = 0; + for(int i=0; i Date: Tue, 20 Nov 2012 09:51:07 +0100 Subject: [PATCH 02/17] fixed: move Vector type under the correct namespace --- dune/elasticity/asmhandler.hpp | 3 --- dune/elasticity/matrixops.hpp | 3 +++ examples/upscale_elasticity.cpp | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/dune/elasticity/asmhandler.hpp b/dune/elasticity/asmhandler.hpp index 7259eb1..a190263 100644 --- a/dune/elasticity/asmhandler.hpp +++ b/dune/elasticity/asmhandler.hpp @@ -11,9 +11,6 @@ //============================================================================== #pragma once -//! \brief A vector holding our RHS -typedef Dune::BlockVector > Vector; - #include #include #include diff --git a/dune/elasticity/matrixops.hpp b/dune/elasticity/matrixops.hpp index e7ef807..a255947 100644 --- a/dune/elasticity/matrixops.hpp +++ b/dune/elasticity/matrixops.hpp @@ -22,6 +22,9 @@ typedef Dune::BCRSMatrix > Matrix; //! \brief For storing matrix adjacency/sparsity patterns typedef std::vector< std::set > AdjacencyPattern; +//! \brief A vector holding our RHS +typedef Dune::BlockVector > Vector; + //! \brief Helper class with some matrix operations class MatrixOps { public: diff --git a/examples/upscale_elasticity.cpp b/examples/upscale_elasticity.cpp index e26a12a..36eaa73 100644 --- a/examples/upscale_elasticity.cpp +++ b/examples/upscale_elasticity.cpp @@ -258,7 +258,7 @@ int main(int argc, char** argv) } Dune::FieldMatrix C; Dune::VTKWriter vtkwriter(grid.leafView()); - Vector field[6]; + Opm::Elasticity::Vector field[6]; std::cout << "assembling..." << "\n"; upscale.assemble(-1,true); upscale.setupSolvers(p.solver); From e4ad0d62df529ac3a5cdcde183a28e5b53ed24fc Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Tue, 20 Nov 2012 09:56:02 +0100 Subject: [PATCH 03/17] added: MatrixOps::extractDiagonal extracts the diagonal of a matrix as a new matrix --- dune/elasticity/matrixops.hpp | 5 +++++ dune/elasticity/matrixops_impl.hpp | 14 ++++++++++++++ 2 files changed, 19 insertions(+) diff --git a/dune/elasticity/matrixops.hpp b/dune/elasticity/matrixops.hpp index a255947..03ed906 100644 --- a/dune/elasticity/matrixops.hpp +++ b/dune/elasticity/matrixops.hpp @@ -55,6 +55,11 @@ class MatrixOps { //! \param[in] symmetric If true, augment symmetrically static Matrix augment(const Matrix& A, const Matrix& B, size_t r0, size_t c0, bool symmetric); + + //! \brief Extract the diagonal of a matrix into a new matrix + //! \param[in] The matrix to extract the diagonal from + //! \returns M = diag(A) + static Matrix extractDiagonal(const Matrix& A); }; #include "matrixops_impl.hpp" diff --git a/dune/elasticity/matrixops_impl.hpp b/dune/elasticity/matrixops_impl.hpp index 7f6253c..7da51d7 100644 --- a/dune/elasticity/matrixops_impl.hpp +++ b/dune/elasticity/matrixops_impl.hpp @@ -144,3 +144,17 @@ Matrix MatrixOps::augment(const Matrix& A, const Matrix& B, return result; } + +Matrix MatrixOps::extractDiagonal(const Matrix& A) +{ + AdjacencyPattern adj; + adj.resize(A.M()); + for (size_t i=0;i Date: Tue, 20 Nov 2012 09:57:17 +0100 Subject: [PATCH 04/17] added: MatrixOps::saveAsc save a matrix as a dense asc file for easy import into e.g. octave note this is only useful for debugging small systems, as the asc files grow very big and unmanageable --- dune/elasticity/matrixops.hpp | 6 ++++++ dune/elasticity/matrixops_impl.hpp | 30 ++++++++++++++++++++++++++++++ 2 files changed, 36 insertions(+) diff --git a/dune/elasticity/matrixops.hpp b/dune/elasticity/matrixops.hpp index 03ed906..e2aefc1 100644 --- a/dune/elasticity/matrixops.hpp +++ b/dune/elasticity/matrixops.hpp @@ -60,6 +60,12 @@ class MatrixOps { //! \param[in] The matrix to extract the diagonal from //! \returns M = diag(A) static Matrix extractDiagonal(const Matrix& A); + + //! \brief Save a matrix as a dense asc file + //! \param[in] A The matrix to save + //! \param[in] file File name + //! \details This is only useful for debugging as the files grow very big + static void saveAsc(const Matrix& A, const std::string& file); }; #include "matrixops_impl.hpp" diff --git a/dune/elasticity/matrixops_impl.hpp b/dune/elasticity/matrixops_impl.hpp index 7da51d7..a4c6de5 100644 --- a/dune/elasticity/matrixops_impl.hpp +++ b/dune/elasticity/matrixops_impl.hpp @@ -158,3 +158,33 @@ Matrix MatrixOps::extractDiagonal(const Matrix& A) return result; } + +void MatrixOps::saveAsc(const Matrix& A, const std::string& file) +{ + std::ofstream f; + f.open(file.c_str()); + f << "% " << A.N() << " " << A.M() << std::endl; + int prevrow=-1; + for (Matrix::ConstIterator it = A.begin(); + it != A.end(); ++it) { + for (int i=0;ibegin(); + it2 != it->end();++it2) { + for (int j=0;j Date: Tue, 20 Nov 2012 09:58:15 +0100 Subject: [PATCH 05/17] fixed: make sure we do not dereference null ptrs --- dune/elasticity/asmhandler_impl.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dune/elasticity/asmhandler_impl.hpp b/dune/elasticity/asmhandler_impl.hpp index 8a4dbbe..73b5e16 100644 --- a/dune/elasticity/asmhandler_impl.hpp +++ b/dune/elasticity/asmhandler_impl.hpp @@ -70,7 +70,7 @@ void ASMHandler::addDOF(int row, int erow, } } } - if (S) + if (S && b) (*b)[row] += scale*(*S)[erow]; } From 499f5092abb89afcd93d0429c7d12c9367b9f3e1 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Tue, 20 Nov 2012 10:17:48 +0100 Subject: [PATCH 06/17] added: support for higher order lagrangian multipliers - add PN shape function support - make the way we group boundary elements into pillars more robust - use the new shape functions in mortar matrix assembly - add options to specify multiplier order (defaults to linears for now) --- dune/elasticity/boundarygrid.hh | 18 ++ dune/elasticity/elasticity_upscale.hpp | 10 +- dune/elasticity/elasticity_upscale_impl.hpp | 125 +++++++--- dune/elasticity/shapefunctions.hpp | 245 +++++++++++++++++++- examples/upscale_elasticity.cpp | 8 +- 5 files changed, 360 insertions(+), 46 deletions(-) diff --git a/dune/elasticity/boundarygrid.hh b/dune/elasticity/boundarygrid.hh index b4a274e..6978f65 100644 --- a/dune/elasticity/boundarygrid.hh +++ b/dune/elasticity/boundarygrid.hh @@ -113,6 +113,13 @@ class BoundaryGrid { //! \param[in] quad The quad to add void add(const Quad& quad); + void addToColumn(size_t col, const Quad& quad) + { + if (col >= colGrids.size()) + colGrids.resize(col+1); + colGrids[col].push_back(quad); + } + //! \brief Obtain a reference to a quad //! \param[in] index The index of the requested quad //! \returns A reference to the requested quad @@ -129,12 +136,22 @@ class BoundaryGrid { return grid[index]; } + const Quad& getQuad(int col, int index) const + { + return colGrids[col][index]; + } + //! \brief Obtain the number of quads in the grid size_t size() const { return grid.size(); } + size_t colSize(int i) const + { + return colGrids[i].size(); + } + //! \brief Return the total number of nodes on the grid when known //! \sa uniform size_t totalNodes() const @@ -209,6 +226,7 @@ class BoundaryGrid { protected: //! \brief Our quadrilateral elements std::vector grid; + std::vector > colGrids; //! \brief Whether or not a given node is marked as fixed std::vector fixNodes; diff --git a/dune/elasticity/elasticity_upscale.hpp b/dune/elasticity/elasticity_upscale.hpp index d2a4c4f..20545cf 100644 --- a/dune/elasticity/elasticity_upscale.hpp +++ b/dune/elasticity/elasticity_upscale.hpp @@ -151,8 +151,11 @@ class ElasticityUpscale //! \param[in] max The maximum coordinates of the grid //! \param[in] n1 The number of elements on the lambda grid in the X direction //! \param[in] n2 The number of elements on the lambda grid in the Y direction + //! \param[in] p1 The order of multipliers in the X direction + //! \param[in] p2 The order of multipliers in the Y direction void periodicBCsMortar(const double* min, - const double* max, int n1, int n2); + const double* max, int n1, int n2, + int p1, int p2); //! \brief Assemble (optionally) stiffness matrix A and load vector //! \param[in] loadcase The strain load case. Set to -1 to skip @@ -283,8 +286,11 @@ class ElasticityUpscale //! \param[in] b1 The primal boundary to match //! \param[in] interface The interface/multiplier grid //! \param[in] dir The normal direction on the boundary (0/1) + //! \param[in] n1 The multipler order in the first direction + //! \param[in] n2 The multipler order in the second direction Matrix findLMatrixMortar(const BoundaryGrid& b1, - const BoundaryGrid& interface, int dir); + const BoundaryGrid& interface, int dir, + int n1, int n2); //! \brief This function loads and maps materials to active grid cells //! \param[in] file The eclipse grid to read materials from diff --git a/dune/elasticity/elasticity_upscale_impl.hpp b/dune/elasticity/elasticity_upscale_impl.hpp index 5b8f52d..b619687 100644 --- a/dune/elasticity/elasticity_upscale_impl.hpp +++ b/dune/elasticity/elasticity_upscale_impl.hpp @@ -56,6 +56,9 @@ BoundaryGrid ElasticityUpscale::extractMasterFace(Direction dir, int c = 0; int i = log2(dir); BoundaryGrid result; + // we first group nodes into this map through the coordinate of lower left + // vertex. we then split this up into pillars for easy processing later + std::map > nodeMap; for (LeafIterator cell = gv.leafView().template begin<0>(); cell != cellend; ++cell, ++c) { std::vector verts; @@ -99,10 +102,27 @@ BoundaryGrid ElasticityUpscale::extractMasterFace(Direction dir, q.v[2] = maxXmaxY(verts); q.v[3] = minXmaxY(verts); } + std::map >::iterator it; + for (it = nodeMap.begin(); it != nodeMap.end(); ++it) { + if (fabs(it->first-q.v[0].c[0]) < 1.e-7) { + it->second.push_back(q); + break; + } + } + if (it == nodeMap.end()) + nodeMap[q.v[0].c[0]].push_back(q); + result.add(q); } } + int p=0; + std::map >::const_iterator it; + for (it = nodeMap.begin(); it != nodeMap.end(); ++it, ++p) { + for (size_t i=0;isecond.size();++i) + result.addToColumn(p,it->second[i]); + } + return result; } @@ -301,34 +321,65 @@ Matrix ElasticityUpscale::findLMatrixLLM(const SlaveGrid& slave, return L; } +static std::vector< std::vector > renumber(const BoundaryGrid& b, + int n1, int n2) +{ + std::vector > nodes; + nodes.resize(b.size()); + // loop over elements + int ofs = 0; + for (size_t e=0; e < b.size(); ++e) { + // first direction major ordered nodes within each element + for (int i2=0; i2 < n2; ++i2) { + if (e != 0) + nodes[e].push_back(nodes[e-1][i2*n1+n1-1]); + for (int i1=(e==0?0:1); i1 < n1; ++i1) + nodes[e].push_back(ofs++); + } + } + + return nodes; +} + template Matrix ElasticityUpscale::findLMatrixMortar(const BoundaryGrid& b1, - const BoundaryGrid& interface, - int dir) + const BoundaryGrid& interface, + int dir, int n1, int n2) { std::vector< std::set > adj; adj.resize(A.getEqns()); +#define QINDEX(p,q,dir) (q) + + // get a set of P1 shape functions for the displacements + P1ShapeFunctionSet ubasis = + P1ShapeFunctionSet::instance(); + + // get a set of PN shape functions for the multipliers + PNShapeFunctionSet<2> lbasis(n1+1, n2+1); + + std::vector > lnodes = renumber(interface,n1,n2); + int totalNodes = lnodes.back().back()+1; + // process pillar by pillar - size_t per_pillar = b1.size()/interface.size(); for (size_t p=0;pgetNoMaster();++n) { + for (size_t n=0;ngetNoMaster();++n) { int dof = A.getEquationForDof(mpc->getMaster(n).node,d); if (dof > -1) { - for (int j=0;j<4;++j) - adj[dof].insert(3*interface[p].v[j].i+d); + for (size_t j=0;j -1) { - for (int j=0;j<4;++j) - adj[dof].insert(3*interface[p].v[j].i+d); + for (size_t j=0;j::findLMatrixMortar(const BoundaryGrid& b1, } Matrix B; - MatrixOps::fromAdjacency(B,adj,A.getEqns(),3*interface.totalNodes()); + MatrixOps::fromAdjacency(B,adj,A.getEqns(),3*totalNodes); - // get a set of P1 shape functions for the velocity - P1ShapeFunctionSet ubasis = P1ShapeFunctionSet::instance(); - // get a set of P1 shape functions for multipliers - P1ShapeFunctionSet lbasis = P1ShapeFunctionSet::instance(); // get a reference element Dune::GeometryType gt; gt.makeCube(2); - const Dune::template GenericReferenceElement &ref = - Dune::GenericReferenceElements::general(gt); // get a quadrature rule + int quadorder = std::max((1.0+n1+0.5)/2.0,(1.0+n2+0.5)/2.0); + quadorder = std::max(quadorder, 2); const Dune::QuadratureRule& rule = - Dune::QuadratureRules::rule(gt,2); + Dune::QuadratureRules::rule(gt,quadorder); // do the assembly loop typename Dune::QuadratureRule::const_iterator r; + Dune::DynamicMatrix E(ubasis.n,(n1+1)*(n2+1),0.0); for (size_t p=0;p lg(qi); - for (size_t q=0;q hex(qu,gv,dir); - Dune::FieldMatrix E; E = 0; for (r = rule.begin(); r != rule.end();++r) { ctype detJ = hex.integrationElement(r->position()); - if (detJ < 0) - continue; + if (detJ < 1.e-4 && r == rule.begin()) + std::cerr << "warning: interface cell (close to) degenerated, |J|=" << detJ << std::endl; + typename HexGeometry<2,2,GridType>::LocalCoordinate loc = lg.local(hex.global(r->position())); - for (int i=0;i<4;++i) { - for (int j=0;j<4;++j) + assert(loc[0] <= 1.0+1.e-4 && loc[0] >= 0.0 && loc[1] <= 1.0+1.e-4 && loc[1] >= 0.0); + for (int i=0;iposition())* lbasis[j].evaluateFunction(loc)*detJ*r->weight(); + } } } // and assemble element contributions @@ -379,11 +429,11 @@ Matrix ElasticityUpscale::findLMatrixMortar(const BoundaryGrid& b1, for (int i=0;i<4;++i) { MPC* mpc = A.getMPC(qu.v[i].i,d); if (mpc) { - for (int n=0;ngetNoMaster();++n) { + for (size_t n=0;ngetNoMaster();++n) { int indexi = A.getEquationForDof(mpc->getMaster(n).node,d); if (indexi > -1) { - for (int j=0;j<4;++j) { - int indexj = qi.v[j].i*3+d; + for (size_t j=0;j::findLMatrixMortar(const BoundaryGrid& b1, } else { int indexi = A.getEquationForDof(qu.v[i].i,d); if (indexi > -1) { - for (int j=0;j<4;++j) { - int indexj = qi.v[j].i*3+d; + for (size_t j=0;j::periodicBCs(const double* min, template void ElasticityUpscale::periodicBCsMortar(const double* min, const double* max, - int n1, int n2) + int n1, int n2, + int p1, int p2) { // this method // 1. fixes the primal corner dofs @@ -844,13 +895,13 @@ void ElasticityUpscale::periodicBCsMortar(const double* min, BoundaryGrid lambday = BoundaryGrid::uniform(fmin,fmax,n1,1,true); // step 5 - Matrix L1 = findLMatrixMortar(master[0],lambdax,0); - Matrix L2 = findLMatrixMortar(master[1],lambdax,0); + Matrix L1 = findLMatrixMortar(master[0],lambdax,0, p2, 1); + Matrix L2 = findLMatrixMortar(master[1],lambdax,0, p2, 1); L.push_back(MatrixOps::Axpy(L1,L2,-1)); // step 6 - Matrix L3 = findLMatrixMortar(master[2],lambday,1); - Matrix L4 = findLMatrixMortar(master[3],lambday,1); + Matrix L3 = findLMatrixMortar(master[2],lambday,1, p1, 1); + Matrix L4 = findLMatrixMortar(master[3],lambday,1, p1, 1); L.push_back(MatrixOps::Axpy(L3,L4,-1)); } diff --git a/dune/elasticity/shapefunctions.hpp b/dune/elasticity/shapefunctions.hpp index b757c09..3047c9d 100644 --- a/dune/elasticity/shapefunctions.hpp +++ b/dune/elasticity/shapefunctions.hpp @@ -6,12 +6,15 @@ //! //! \author Arne Morten Kvarving / SINTEF //! -//! \brief Class for linear shape functions. Loosely based on code in dune-grid-howto +//! \brief Classes for shape functions. Loosely based on code in dune-grid-howto //! //============================================================================== #pragma once #include +#include "dynmatrixev.hh" + +#include namespace Opm { namespace Elasticity { @@ -80,6 +83,96 @@ class LinearShapeFunction Dune::FieldVector coeff1; }; +//! \brief Represents a cardinal function on a line + template +class LagrangeCardinalFunction +{ + public: + //! \brief Empty default constructor + LagrangeCardinalFunction() {} + + //! \brief Construct a cardinal function with the given nodes + //! \param[in] nodes_ The nodes + //! \param[in] i The node this function is associated with + LagrangeCardinalFunction(const std::vector& nodes_, + size_t i) + : nodes(nodes_), node(i) {} + + //! \brief Evaluate the shape function + //! \param[in] local The local coordinates + rtype evaluateFunction(const ctype& local) const + { + rtype result = 1; + for (size_t i=0; i < nodes.size(); ++i) { + if (i != node) + result *= (local-nodes[i])/(nodes[node]-nodes[i]); + } + + return result; + } + + //! \brief Evaluate the derivative of the cardinal function + //! \param[in] local The local coordinates + rtype evaluateGradient(const ctype& local) const + { + rtype result = 0; + for (size_t i=0; i < nodes.size(); ++i) { + rtype f = 1; + for (int j=0; j < nodes.size(); ++j) { + if (i != j && j != node) + f *= (local-nodes[j])/(nodes[node]-nodes[j]); + } + result += f/(nodes[node]-nodes[i]); + } + + return result; + } + + private: + //! \brief The nodes + std::vector nodes; + + size_t node; +}; + +//! \brief Represents a tensor-product of 1D functions + template +class TensorProductFunction +{ + public: + //! \brief The dimension of the function + enum { dimension = dim }; + + //! \brief Empty default constructor + TensorProductFunction() {} + + //! \brief Construct a tensor-product function + //! \param[in] funcs_ The functions + TensorProductFunction(const Dune::FieldVector& funcs_) + : funcs(funcs_) {} + + //! \brief Evaluate the function + //! \param[in] local The local coordinates + rtype evaluateFunction(const Dune::FieldVector& local) const + { + rtype result = 1; + for (int i=0; i < dim; ++i) + result *= funcs[i].evaluateFunction(local[i]); + + return result; + } + + Dune::FieldVector + evaluateGradient(const Dune::FieldVector& local) const + { + Dune::FieldVector result; + for (int i=0; i < dim; ++i) + result[i] = funcs[i].evaluateGradient(local[i]); + } + private: + Dune::FieldVector funcs; +}; + //! \brief Singleton handler for the set of LinearShapeFunction template class P1ShapeFunctionSet @@ -121,8 +214,8 @@ private: 0, 1, 1, 0, 0, 0}; - static rtype coeffs22[] = {-1,-1, - 1,-1, + static rtype coeffs22[] = {-1,-1, + 1,-1, -1, 1, 1, 1}; @@ -134,8 +227,8 @@ private: 0, 1, 0, 1, 0, 0, 0, 0, 0}; - static rtype coeffs32[] = {-1,-1,-1, - 1,-1,-1, + static rtype coeffs32[] = {-1,-1,-1, + 1,-1,-1, -1, 1,-1, 1, 1,-1, -1,-1, 1, @@ -173,5 +266,147 @@ private: ShapeFunction f[n]; }; + template +class PNShapeFunctionSet +{ +public: + typedef LagrangeCardinalFunction CardinalFunction; + + typedef TensorProductFunction + ShapeFunction; + + PNShapeFunctionSet(int n1, int n2, int n3=0) + { + int dims[3] = {n1, n2, n3}; + cfuncs.resize(dim); + for (int i=0; i < dim; ++i) { + std::vector grid; + grid = gaussLobattoLegendreGrid(dims[i]); + for (int j=0;j fs; + if (dim == 3) { + f.resize(n1*n2*n3); + for (int k=0; k < n3; ++k) { + for (int j=0; j < n2; ++j) + for (int i=0; i< n1; ++i) { + fs[0] = cfuncs[0][i]; + fs[1] = cfuncs[1][j]; + fs[2] = cfuncs[2][k]; + f[l++] = ShapeFunction(fs); + } + } + } else { + f.resize(n1*n2); + for (int j=0; j < n2; ++j) { + for (int i=0; i< n1; ++i) { + fs[0] = cfuncs[0][i]; + fs[1] = cfuncs[1][j]; + f[l++] = ShapeFunction(fs); + } + } + } + } + + //! \brief Obtain a given shape function + //! \param[in] i The requested shape function + const ShapeFunction& operator[](int i) const + { + return f[i]; + } + + int size() + { + return f.size(); + } +protected: + std::vector< std::vector > cfuncs; + std::vector f; + + double legendre(double x, int n) + { + std::vector Ln; + Ln.resize(n+1); + Ln[0] = 1.f; + Ln[1] = x; + if( n > 1 ) { + for( int i=1;i Ln; + Ln.resize(n+1); + + Ln[0] = 1.0; Ln[1] = x; + + if( (x == 1.0) || (x == -1.0) ) + return( pow(x,n-1)*n*(n+1.0)/2.0 ); + else { + for( int i=1;i gaussLegendreGrid(int n) + { + Dune::DynamicMatrix A(n,n,0.0); + + A[0][1] = 1.f; + for (int i=1;i > eigenValues(n); + Dune::DynamicMatrixHelp::eigenValuesNonSym(A, eigenValues); + + std::vector result(n); + for (int i=0; i < n; ++i) + result[i] = std::real(eigenValues[i]); + std::sort(result.begin(),result.begin()+n); + + return result; + } + + std::vector gaussLobattoLegendreGrid(int n) + { + assert(n > 1); + const double tolerance = 1.e-15; + + std::vector result(n); + result[0] = 0.0; + result[n-1] = 1.0; + if (n == 3) + result[1] = 0.5; + + if (n < 4) + return result; + + std::vector glgrid = gaussLegendreGrid(n-1); + for (int i=1;i tolerance) { + old = result[i]; + double L = legendre(old, n-1); + double Ld = legendreDerivative(old, n-1); + result[i] += (1.0-old*old)*Ld/((n-1.0)*n*L); + } + result[i] = (result[i]+1.0)/2.0; + } + + return result; + } +}; + } } diff --git a/examples/upscale_elasticity.cpp b/examples/upscale_elasticity.cpp index 36eaa73..0b632d9 100644 --- a/examples/upscale_elasticity.cpp +++ b/examples/upscale_elasticity.cpp @@ -91,6 +91,8 @@ struct Params { int cellsy; //! \brief cell in z int cellsz; + //! \brief Polynomial order of multipliers in first / second direction + int lambda[2]; }; //! \brief Parse the command line arguments @@ -104,6 +106,8 @@ void parseCommandLine(int argc, char** argv, Params& p) p.min[1] = param.getDefault("ymin",-1); p.min[2] = param.getDefault("zmin",-1); //std::string method = param.getDefault("method","mpc"); + p.lambda[0] = param.getDefault("lambdax", 1); + p.lambda[1] = param.getDefault("lambday", 1); std::string method = param.getDefault("method","mortar"); p.LLM = (strcasecmp(method.c_str(),"llm") == 0); //p.mortar = (strcasecmp(method.c_str(),"mortar") == 0); @@ -254,8 +258,8 @@ int main(int argc, char** argv) upscale.A.initForAssembly(); } else { std::cout << "using Mortar couplings.." << std::endl; - upscale.periodicBCsMortar(p.min,p.max,p.n1,p.n2); - } + upscale.periodicBCsMortar(p.min,p.max,p.n1,p.n2,p.lambda[0], p.lambda[1]); + } Dune::FieldMatrix C; Dune::VTKWriter vtkwriter(grid.leafView()); Opm::Elasticity::Vector field[6]; From da7729c54f28edb20713bb0b5a37796260ac37c2 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Tue, 20 Nov 2012 11:06:09 +0100 Subject: [PATCH 07/17] changed: refactor method setting into an enum instead of a set of bools --- examples/upscale_elasticity.cpp | 37 +++++++++++++++++++-------------- 1 file changed, 21 insertions(+), 16 deletions(-) diff --git a/examples/upscale_elasticity.cpp b/examples/upscale_elasticity.cpp index 0b632d9..711a74c 100644 --- a/examples/upscale_elasticity.cpp +++ b/examples/upscale_elasticity.cpp @@ -51,6 +51,13 @@ void syntax(char** argv) << "\t linsolver_type=slu - use the SuperLU sparse direct solver" << std::endl; } + +enum UpscaleMethod { + UPSCALE_MPC = 1, + UPSCALE_LLM = 2, + UPSCALE_MORTAR = 3 +}; + //! \brief Structure holding parameters configurable from command line struct Params { //! \brief The eclipse grid file @@ -61,12 +68,8 @@ struct Params { std::string vtufile; //! \brief Text output file std::string output; - //! \brief Whether or not to use LLM couplings - bool LLM; - //! \brief Whether or not to use Mortar couplings - bool mortar; - //! \brief Whether or not to use MPC couplings - bool mpc; + //! \brief Method + UpscaleMethod method; //! \brief A scaling parameter for the E-modulus to avoid numerical issues double Emin; //! \brief The tolerance for collapsing nodes in the Z direction @@ -105,14 +108,16 @@ void parseCommandLine(int argc, char** argv, Params& p) p.min[0] = param.getDefault("xmin",-1); p.min[1] = param.getDefault("ymin",-1); p.min[2] = param.getDefault("zmin",-1); - //std::string method = param.getDefault("method","mpc"); p.lambda[0] = param.getDefault("lambdax", 1); p.lambda[1] = param.getDefault("lambday", 1); std::string method = param.getDefault("method","mortar"); - p.LLM = (strcasecmp(method.c_str(),"llm") == 0); - //p.mortar = (strcasecmp(method.c_str(),"mortar") == 0); - p.mpc = (strcasecmp(method.c_str(),"mpc") == 0); - p.Emin = param.getDefault("Emin",0.f); + if (!strcasecmp(method.c_str(),"mpc")) + p.method = UPSCALE_MPC; + if (!strcasecmp(method.c_str(),"llm")) + p.method = UPSCALE_LLM; + if (!strcasecmp(method.c_str(),"mortar")) + p.method = UPSCALE_MORTAR; + p.Emin = param.getDefault("Emin",0.0); p.ctol = param.getDefault("ctol",1.e-8); p.ltol = param.getDefault("ltol",1.e-10); p.file = param.get("gridfilename"); @@ -153,9 +158,9 @@ void writeOutput(const Params& p, Opm::time::StopWatch& watch, int cells, gethostname(hostname,1024); std::string method = "mortar"; - if (p.LLM) + if (p.method == UPSCALE_LLM) method = "llm"; - if (p.mpc) + if (p.method == UPSCALE_MPC) method = "mpc"; // write log @@ -248,15 +253,15 @@ int main(int argc, char** argv) p.n2 = grid.logicalCartesianSize()[1]; } - if (p.LLM) { + if (p.method == UPSCALE_LLM) { std::cout << "using LLM couplings..." << std::endl; upscale.periodicBCsLLM(p.min,p.max,p.n1,p.n2); - } else if (p.mpc) { + } else if (p.method == UPSCALE_MPC) { std::cout << "using MPC couplings in all directions..." << std::endl; upscale.periodicBCs(p.min,p.max); std::cout << "preprocessing grid..." << std::endl; upscale.A.initForAssembly(); - } else { + } else if (p.method == UPSCALE_MORTAR) { std::cout << "using Mortar couplings.." << std::endl; upscale.periodicBCsMortar(p.min,p.max,p.n1,p.n2,p.lambda[0], p.lambda[1]); } From 6f1e511531c5486bfe2dbbeac0148a6081e67660 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Tue, 20 Nov 2012 11:14:09 +0100 Subject: [PATCH 08/17] added: method=none use this to apply no periodicity approach. useful for debugging --- dune/elasticity/elasticity_upscale.hpp | 10 +++++----- examples/upscale_elasticity.cpp | 21 +++++++++++++++------ 2 files changed, 20 insertions(+), 11 deletions(-) diff --git a/dune/elasticity/elasticity_upscale.hpp b/dune/elasticity/elasticity_upscale.hpp index 20545cf..79f3124 100644 --- a/dune/elasticity/elasticity_upscale.hpp +++ b/dune/elasticity/elasticity_upscale.hpp @@ -157,6 +157,11 @@ class ElasticityUpscale const double* max, int n1, int n2, int p1, int p2); + //! \brief Fix corner nodes + //! \param[in] min The minimum coordinates on the grid + //! \param[in] max The maximum coordinates on the grid + void fixCorners(const double* min, const double* max); + //! \brief Assemble (optionally) stiffness matrix A and load vector //! \param[in] loadcase The strain load case. Set to -1 to skip //! \param[in] matrix Whether or not to assemble the matrix @@ -213,11 +218,6 @@ class ElasticityUpscale //! \brief Vector holding material parameters for each active grid cell std::vector materials; - //! \brief Fix corner nodes - //! \param[in] min The minimum coordinates on the grid - //! \param[in] max The maximum coordinates on the grid - void fixCorners(const double* min, const double* max); - //! \brief Extract the vertices on a given face //! \param[in] dir The direction of the face normal //! \param[in] coord The coordinate of the face plane diff --git a/examples/upscale_elasticity.cpp b/examples/upscale_elasticity.cpp index 711a74c..6fcdca1 100644 --- a/examples/upscale_elasticity.cpp +++ b/examples/upscale_elasticity.cpp @@ -53,6 +53,7 @@ void syntax(char** argv) enum UpscaleMethod { + UPSCALE_NONE = 0, UPSCALE_MPC = 1, UPSCALE_LLM = 2, UPSCALE_MORTAR = 3 @@ -102,12 +103,12 @@ struct Params { void parseCommandLine(int argc, char** argv, Params& p) { Opm::parameter::ParameterGroup param(argc, argv); - p.max[0] = param.getDefault("xmax",-1); - p.max[1] = param.getDefault("ymax",-1); - p.max[2] = param.getDefault("zmax",-1); - p.min[0] = param.getDefault("xmin",-1); - p.min[1] = param.getDefault("ymin",-1); - p.min[2] = param.getDefault("zmin",-1); + p.max[0] = param.getDefault("xmax",-1); + p.max[1] = param.getDefault("ymax",-1); + p.max[2] = param.getDefault("zmax",-1); + p.min[0] = param.getDefault("xmin",-1); + p.min[1] = param.getDefault("ymin",-1); + p.min[2] = param.getDefault("zmin",-1); p.lambda[0] = param.getDefault("lambdax", 1); p.lambda[1] = param.getDefault("lambday", 1); std::string method = param.getDefault("method","mortar"); @@ -117,6 +118,8 @@ void parseCommandLine(int argc, char** argv, Params& p) p.method = UPSCALE_LLM; if (!strcasecmp(method.c_str(),"mortar")) p.method = UPSCALE_MORTAR; + if (!strcasecmp(method.c_str(),"none")) + p.method = UPSCALE_NONE; p.Emin = param.getDefault("Emin",0.0); p.ctol = param.getDefault("ctol",1.e-8); p.ltol = param.getDefault("ltol",1.e-10); @@ -162,6 +165,8 @@ void writeOutput(const Params& p, Opm::time::StopWatch& watch, int cells, method = "llm"; if (p.method == UPSCALE_MPC) method = "mpc"; + if (p.method == UPSCALE_NONE) + method = "none"; // write log std::ofstream f; @@ -264,6 +269,10 @@ int main(int argc, char** argv) } else if (p.method == UPSCALE_MORTAR) { std::cout << "using Mortar couplings.." << std::endl; upscale.periodicBCsMortar(p.min,p.max,p.n1,p.n2,p.lambda[0], p.lambda[1]); + } else if (p.method == UPSCALE_NONE) { + std::cout << "no periodicity approach applied.." << std::endl; + upscale.fixCorners(p.min, p.max); + upscale.A.initForAssembly(); } Dune::FieldMatrix C; Dune::VTKWriter vtkwriter(grid.leafView()); From 343411cc9b7ad3b6b903f7dfb4235604777fd4ef Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Tue, 20 Nov 2012 11:22:04 +0100 Subject: [PATCH 09/17] quell compiler warnings (sign/unsigned mismatch, unused variables) --- dune/elasticity/asmhandler.hpp | 2 +- dune/elasticity/asmhandler_impl.hpp | 20 ++++----- dune/elasticity/boundarygrid.cpp | 8 ++-- dune/elasticity/boundarygrid.hh | 5 +-- dune/elasticity/elasticity_impl.hpp | 6 +-- dune/elasticity/elasticity_upscale_impl.hpp | 49 +++++++-------------- dune/elasticity/matrixops_impl.hpp | 6 +-- examples/upscale_elasticity.cpp | 5 +-- 8 files changed, 42 insertions(+), 59 deletions(-) diff --git a/dune/elasticity/asmhandler.hpp b/dune/elasticity/asmhandler.hpp index a190263..4f29be0 100644 --- a/dune/elasticity/asmhandler.hpp +++ b/dune/elasticity/asmhandler.hpp @@ -55,7 +55,7 @@ class ASMHandler { //! \brief Get the number of equations in the system //! \returns The number of equations in the system - int getEqns() const + size_t getEqns() const { return maxeqn; } diff --git a/dune/elasticity/asmhandler_impl.hpp b/dune/elasticity/asmhandler_impl.hpp index 73b5e16..7ca9068 100644 --- a/dune/elasticity/asmhandler_impl.hpp +++ b/dune/elasticity/asmhandler_impl.hpp @@ -59,7 +59,7 @@ void ASMHandler::addDOF(int row, int erow, for (int l=0;lgetNoMaster();++n) { + for (size_t n=0;ngetNoMaster();++n) { int idx = meqn[mpc->getMaster(n).node*dim+mpc->getMaster(n).dof-1]; if (idx != -1) A[row][idx] += scale*mpc->getMaster(n).coeff*(*K)[erow][j*dim+l]; @@ -93,7 +93,7 @@ void ASMHandler::addElement( for (int k=0;kgetNoMaster();++n) { + for (size_t n=0;ngetNoMaster();++n) { int idx = meqn[mpc->getMaster(n).node*dim+mpc->getMaster(n).dof-1]; addDOF(idx,i*dim+k,K,S,set,cell,b2,mpc->getMaster(n).coeff); } @@ -124,7 +124,7 @@ void ASMHandler::extractValues(Dune::FieldVector& v, if (it2 != fixedNodes.end() && it2->second.first & (1 << n)) v[l++] = it2->second.second[n]; else if (mpc) { - for (int m=0;mgetNoMaster();++m) { + for (size_t m=0;mgetNoMaster();++m) { int idx = meqn[mpc->getMaster(m).node*dim+mpc->getMaster(m).dof-1]; if (idx != -1) v[l] += u[idx]*mpc->getMaster(m).coeff; @@ -143,7 +143,7 @@ void ASMHandler::expandSolution(Vector& result, const Vector& u) result.resize(nodes*dim); result = 0; int l=0; - for (size_t i=0;i::expandSolution(Vector& result, const Vector& u) } // second loop - handle MPC couplings l = 0; - for (size_t i=0;igetNoMaster();++n) { + for (size_t n=0;ngetNoMaster();++n) { int idx = mpc->getMaster(n).node*dim+mpc->getMaster(n).dof-1; if (meqn[idx] != -1) result[l] += u[meqn[idx]]*mpc->getMaster(n).coeff; @@ -188,7 +188,7 @@ void ASMHandler::addMPC(MPC* mpc) fixIt it = fixedNodes.find(mpc->getSlave().node); int flag = 1 << (mpc->getSlave().dof-1); if (it == fixedNodes.end() || - !(it->second.first & (1 << (mpc->getSlave().dof-1)))) { + !(it->second.first & flag)) { mpcs.insert(std::make_pair(slaveNode,mpc)); return; } @@ -257,7 +257,7 @@ void ASMHandler::resolveMPCChain(MPC* mpc) coeff[master.dof-1] = master.coeff; int removeOld = 0; - for (size_t d = 1; d <= dim; d++) + for (int d = 1; d <= dim; d++) if (fabs(coeff[d-1]) > 1.0e-8) { MPC* mpc2 = getMPC(mpc->getMaster(i).node,d-1); @@ -354,7 +354,7 @@ void ASMHandler::nodeAdjacency(const LeafIterator& it, for (int l=0;lgetNoMaster();++i) { + for (size_t i=0;igetNoMaster();++i) { int idx = meqn[mpc->getMaster(i).node*dim+ mpc->getMaster(i).dof-1]; if (idx != -1) @@ -386,7 +386,7 @@ void ASMHandler::determineAdjacencyPattern() for (int k=0;kgetNoMaster();++l) { + for (size_t l=0;lgetNoMaster();++l) { nodeAdjacency(it,vertexsize, meqn[mpc->getMaster(l).node*dim+ mpc->getMaster(l).dof-1]); diff --git a/dune/elasticity/boundarygrid.cpp b/dune/elasticity/boundarygrid.cpp index d167a08..04bb408 100644 --- a/dune/elasticity/boundarygrid.cpp +++ b/dune/elasticity/boundarygrid.cpp @@ -175,7 +175,7 @@ int BoundaryGrid::Q4inv(FaceCoord& res, const Quad& q, // check that obtained solutions are inside element double tol = 1+epsOut; int nInside=0; - for (int i=0;i 1) { std::cout << "multiple solutions" << std::endl; @@ -232,7 +232,7 @@ bool BoundaryGrid::bilinearSolve(double epsilon, double order, double Q2 = A[0]*B[1]-B[0]*A[1]; std::vector Z; cubicSolve(epsilon,0,Q2,Q1,Q0,Z); - for (int i=0;i tol) { X.push_back(Z[i]); @@ -244,10 +244,10 @@ bool BoundaryGrid::bilinearSolve(double epsilon, double order, Q2 = A[0]*B[2]-B[0]*A[2]; Z.clear(); cubicSolve(epsilon,0,Q2,Q1,Q0,Z); - for (int i=0;i tol) { - int j=0; + size_t j=0; for (j=0;j const typename GridImp::LeafGridView::template Codim<3>::Iterator itend = gv.leafView().template end<3>(); for (int i=0;i<4;++i) { typename GridImp::LeafGridView::template Codim<3>::Iterator it=start; - for (it ; it != itend; ++it) { + for (; it != itend; ++it) { if (mapper.map(*it) == q.v[i].i) break; } @@ -414,7 +414,6 @@ class HexGeometry<2, cdim, GridImp> Dune::GenericReferenceElements< ctype, 2 >::general(type()); LocalCoordinate x = refElement.position(0,0); LocalCoordinate dx; - int i=0; do { using namespace Dune::GenericGeometry; // DF^n dx^n = F^n, x^{n+1} -= dx^n diff --git a/dune/elasticity/elasticity_impl.hpp b/dune/elasticity/elasticity_impl.hpp index 3aa3ff8..34743b8 100644 --- a/dune/elasticity/elasticity_impl.hpp +++ b/dune/elasticity/elasticity_impl.hpp @@ -30,7 +30,7 @@ void Elasticity::getBmatrix(Dune::FieldMatrix::getBmatrix(Dune::FieldMatrix::getBmatrix(Dune::FieldMatrix ElasticityUpscale::extractFace(Direction dir, ctype coord) { std::vector result; - const LeafIndexSet& set = gv.leafView().indexSet(); const LeafVertexIterator itend = gv.leafView().template end(); // make a mapper for codim dim entities in the leaf grid @@ -68,7 +67,7 @@ BoundaryGrid ElasticityUpscale::extractMasterFace(Direction dir, else if (side == RIGHT) idx = set.subIndex(*cell,V2[i][0],dim); LeafVertexIterator it=start; - for (it ; it != itend; ++it) { + for (; it != itend; ++it) { if (mapper.map(*it) == idx) break; } @@ -79,7 +78,7 @@ BoundaryGrid ElasticityUpscale::extractMasterFace(Direction dir, if (side == RIGHT) idx = set.subIndex(*cell,V2[i][j],dim); LeafVertexIterator it=start; - for (it ; it != itend; ++it) { + for (; it != itend; ++it) { if (mapper.map(*it) == idx) break; } @@ -145,7 +144,6 @@ void ElasticityUpscale::findBoundaries(double* min, { max[0] = max[1] = max[2] = -1e5; min[0] = min[1] = min[2] = 1e5; - const LeafIndexSet& set = gv.leafView().indexSet(); const LeafVertexIterator itend = gv.leafView().template end(); // iterate over vertices and find slaves @@ -179,7 +177,7 @@ void ElasticityUpscale::periodicPlane(Direction plane, Direction dir, const std::vector& slave, const BoundaryGrid& master) { - for (int i=0;i::findBMatrixLLM(const SlaveGrid& slave) continue; MPC* mpc = A.getMPC(it->first,k); if (mpc) { - for (int n=0;ngetNoMaster();++n) { + for (size_t n=0;ngetNoMaster();++n) { int idx = A.getEquationForDof(mpc->getMaster(n).node, mpc->getMaster(n).dof-1); if (idx != -1) @@ -227,7 +225,7 @@ Matrix ElasticityUpscale::findBMatrixLLM(const SlaveGrid& slave) continue; MPC* mpc = A.getMPC(it->first,k); if (mpc) { - for (int n=0;ngetNoMaster();++n) { + for (size_t n=0;ngetNoMaster();++n) { int idx = A.getEquationForDof(mpc->getMaster(n).node, mpc->getMaster(n).dof-1); if (idx != -1) @@ -251,7 +249,7 @@ Matrix ElasticityUpscale::findLMatrixLLM(const SlaveGrid& slave, int nbeqn=0; std::vector dofmap(master.totalNodes()*dim,-1); int col=0; - for (int i=0;i::assemble(int loadcase, bool matrix) { const int comp = 3+(dim-2)*3; static const int bfunc = 4+(dim-2)*4; - const int N = gv.size(dim); - const LeafIndexSet& set = gv.leafView().indexSet(); const LeafIterator itend = gv.leafView().template end<0>(); Dune::FieldMatrix C; @@ -567,8 +563,6 @@ void ElasticityUpscale::assemble(int loadcase, bool matrix) materials[m++]->getConstitutiveMatrix(C); // determine geometry type of the current element and get the matching reference element Dune::GeometryType gt = it->type(); - const Dune::template GenericReferenceElement &ref = - Dune::GenericReferenceElements::general(gt); Dune::FieldMatrix Aq; K = 0; @@ -616,14 +610,9 @@ void ElasticityUpscale::averageStress(Dune::FieldVector& s return; static const int bfunc = 4+(dim-2)*4; - const int N = gv.size(dim); - const LeafIndexSet& set = gv.leafView().indexSet(); const LeafIterator itend = gv.leafView().template end<0>(); - // get a set of P1 shape functions - P1ShapeFunctionSet basis = P1ShapeFunctionSet::instance(); - Dune::FieldMatrix C; Dune::FieldVector eps0; eps0 = 0; @@ -635,8 +624,6 @@ void ElasticityUpscale::averageStress(Dune::FieldVector& s materials[m++]->getConstitutiveMatrix(C); // determine geometry type of the current element and get the matching reference element Dune::GeometryType gt = it->type(); - const Dune::template GenericReferenceElement &ref = - Dune::GenericReferenceElements::general(gt); Dune::FieldVector v; A.extractValues(v,u,it); @@ -710,7 +697,7 @@ void ElasticityUpscale::loadMaterialsFromGrid(const std::string& file) // scale E modulus of materials if (Escale > 0) { Emin = *std::min_element(Emod.begin(),Emod.end()); - for (int i=0;i::loadMaterialsFromRocklist(const std::string& f // scale E modulus of materials if (Escale > 0) { Emin=1e10; - for (int i=0;igetE()); - for (int i=0;igetE(); ((Isotropic*)cache[i])->setE(E*Escale/Emin); } @@ -780,7 +767,7 @@ void ElasticityUpscale::loadMaterialsFromRocklist(const std::string& f std::vector volume; volume.resize(cache.size()); if (file == "uniform") { - for (size_t i=0;i::periodicBCsLLM(const double* min, Dune::LeafMultipleCodimMultipleGeomTypeMapper mapper(gv); - LeafVertexIterator start=gv.leafView().template begin(); - LeafVertexIterator itend = gv.leafView().template end(); BoundaryGrid::FaceCoord fmin,fmax; // YZ plane @@ -969,10 +954,10 @@ void ElasticityUpscale::periodicBCsLLM(const double* min, // step 5 std::map m; // find matching coefficients - for (int i=0;i::setupPreconditioner() if (cell / gv.logicalCartesianSize()[2] > 0 && cell % gv.logicalCartesianSize()[2] == 0) currdomain++; - for (int i=0;i<8;++i) { + for (size_t i=0;i<8;++i) { int idx=set.subIndex(*it,i,dim); for (int d=0;d<3;++d) { MPC* mpc = A.getMPC(idx,d); if (mpc) { - for (int n=0;ngetNoMaster();++n) { + for (size_t n=0;ngetNoMaster();++n) { int row = A.getEquationForDof(mpc->getMaster(n).node,d); if (row > -1) rows[currdomain].insert(row); @@ -1034,11 +1019,11 @@ void ElasticityUpscale::setupSolvers(Solver solver) // [ 0 0 0 L3' L4' 0 0] [ub_2] [0] int r = A.getOperator().N(); int c = A.getOperator().M(); - for (int i=0;i::setupSolvers(Solver solver) // [L1' 0 0] [l_1] = [0] // [L2' 0 0] [l_2] [0] int c = A.getOperator().M(); - for (int i=0;i >& adj A.setSize(rows, cols, sum); A.setBuildMode(Matrix::random); - for (int i = 0; i < adj.size(); i++) + for (size_t i = 0; i < adj.size(); i++) A.setrowsize(i,adj[i].size()); A.endrowsizes(); @@ -53,7 +53,7 @@ Matrix MatrixOps::Axpy(const Matrix& A, const Matrix& B, double alpha) assert(A.M() == B.M() && A.N() == B.N()); // establish union adjacency pattern - std::vector > adj; + AdjacencyPattern adj; adj.resize(A.N()); for (Matrix::ConstIterator it = A.begin(); it != A.end(); ++it) { @@ -120,7 +120,7 @@ Matrix MatrixOps::augment(const Matrix& A, const Matrix& B, } if (symmetric) { // always establish diagonal elements or superLU crashes - for (int i=0;i Date: Tue, 20 Nov 2012 11:23:05 +0100 Subject: [PATCH 10/17] fixed: make sure we load vectors --- dune/elasticity/elasticity_upscale_impl.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/dune/elasticity/elasticity_upscale_impl.hpp b/dune/elasticity/elasticity_upscale_impl.hpp index 33574e8..19bf41d 100644 --- a/dune/elasticity/elasticity_upscale_impl.hpp +++ b/dune/elasticity/elasticity_upscale_impl.hpp @@ -552,6 +552,7 @@ void ElasticityUpscale::assemble(int loadcase, bool matrix) EP = &ES; eps0[loadcase] = 1; A.getLoadVector() = 0; + b[loadcase] = 0; } int m=0; Dune::FieldMatrix* KP=0; From 29ec9824e9ed5d63d263b85ab9368f3c64cdd3cd Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Tue, 20 Nov 2012 11:23:52 +0100 Subject: [PATCH 11/17] added: warn if volume cell is degenerated also gives the collapse tolerance that will get rid of the cell --- dune/elasticity/elasticity_upscale_impl.hpp | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/dune/elasticity/elasticity_upscale_impl.hpp b/dune/elasticity/elasticity_upscale_impl.hpp index 19bf41d..93e0ccc 100644 --- a/dune/elasticity/elasticity_upscale_impl.hpp +++ b/dune/elasticity/elasticity_upscale_impl.hpp @@ -578,8 +578,13 @@ void ElasticityUpscale::assemble(int loadcase, bool matrix) it->geometry().jacobianInverseTransposed(r->position()); ctype detJ = it->geometry().integrationElement(r->position()); - if (detJ <= 0) - continue; + if (detJ <= 1.e-4) { + std::cout << "cell " << m << " is (close to) degenerated, detJ " << detJ << std::endl; + double zdiff=0.0; + for (int i=0;i<4;++i) + zdiff = std::max(zdiff, it->geometry().corner(i+4)[2]-it->geometry().corner(i)[2]); + std::cout << " - Consider setting ctol larger than " << zdiff << std::endl; + } Dune::FieldMatrix B; E.getBmatrix(B,r->position(),jacInvTra); From cd6c93a28b3d6d9e7c3bc0ed66998bda8a95d150 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Tue, 20 Nov 2012 11:24:51 +0100 Subject: [PATCH 12/17] fixed: warn and exit instead of segfault if model uses more materials than is in rock list --- dune/elasticity/elasticity_upscale_impl.hpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/dune/elasticity/elasticity_upscale_impl.hpp b/dune/elasticity/elasticity_upscale_impl.hpp index 93e0ccc..2d96ebc 100644 --- a/dune/elasticity/elasticity_upscale_impl.hpp +++ b/dune/elasticity/elasticity_upscale_impl.hpp @@ -782,6 +782,10 @@ void ElasticityUpscale::loadMaterialsFromRocklist(const std::string& f std::vector cells = gv.globalCell(); for (size_t i=0;i= cache.size()) { + std::cerr << "Material " << satnum[k] << " referenced but not available. Check your rocklist." << std::endl; + exit(1); + } materials.push_back(cache[satnum[k]-1]); volume[satnum[k]-1] += gv.cellVolume(i); } From a7dc718f2caf9f9be293b72f3da5b15437192b9d Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Tue, 20 Nov 2012 11:27:21 +0100 Subject: [PATCH 13/17] changed: use human friendly 1 based loadcase numbering --- examples/upscale_elasticity.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/upscale_elasticity.cpp b/examples/upscale_elasticity.cpp index e70d706..fb6bc3d 100644 --- a/examples/upscale_elasticity.cpp +++ b/examples/upscale_elasticity.cpp @@ -282,7 +282,7 @@ int main(int argc, char** argv) #pragma omp parallel for schedule(static) for (int i=0;i<6;++i) { upscale.assemble(i,false); - std::cout << "solving case " << i << "..." << "\n"; + std::cout << "solving case " << i+1 << "..." << "\n"; upscale.solve(p.solver,p.ltol,i); upscale.A.expandSolution(field[i],upscale.u[i]); #define CLAMP(x) (fabs(x)<1.e-5?0.f:x) From eaca1c97ad7d610af41fb62fb250e489777eae29 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Wed, 2 Jan 2013 14:02:09 +0100 Subject: [PATCH 14/17] added: improve the linear algebra - change linsolver_type=cg to linsolver_type=iterative - rename linsolver_type=slu to linsolver_type=direct - adds a schur based preconditioner for the mortar system - can use a symmetric iterative scheme (minres) or an asymmetric (gmres) - allows employing a uzawa approach (mostly for debugging) the elasticity block and (through schur-decomposition) multiplier blocks are both preconditioned using an AMG. a diagonal preconditioner for the multipliers are also available. --- dune/elasticity/applier.hpp | 47 ++ dune/elasticity/asmhandler.hpp | 8 + dune/elasticity/asmhandler_impl.hpp | 22 +- dune/elasticity/elasticity_upscale.hpp | 250 +++++--- dune/elasticity/elasticity_upscale_impl.hpp | 655 ++++++++------------ dune/elasticity/logutils.hpp | 31 + dune/elasticity/matrixops.hpp | 15 + dune/elasticity/matrixops_impl.hpp | 75 ++- dune/elasticity/mortar_evaluator.hpp | 79 +++ dune/elasticity/mortar_schur.hpp | 73 +++ dune/elasticity/mortar_schur_precond.hpp | 123 ++++ dune/elasticity/mortar_utils.hpp | 31 + dune/elasticity/seqlu.hpp | 77 +++ dune/elasticity/uzawa_solver.hpp | 56 ++ examples/upscale_elasticity.cpp | 140 +++-- 15 files changed, 1146 insertions(+), 536 deletions(-) create mode 100644 dune/elasticity/applier.hpp create mode 100644 dune/elasticity/logutils.hpp create mode 100644 dune/elasticity/mortar_evaluator.hpp create mode 100644 dune/elasticity/mortar_schur.hpp create mode 100644 dune/elasticity/mortar_schur_precond.hpp create mode 100644 dune/elasticity/mortar_utils.hpp create mode 100644 dune/elasticity/seqlu.hpp create mode 100644 dune/elasticity/uzawa_solver.hpp diff --git a/dune/elasticity/applier.hpp b/dune/elasticity/applier.hpp new file mode 100644 index 0000000..7db2c4c --- /dev/null +++ b/dune/elasticity/applier.hpp @@ -0,0 +1,47 @@ +//============================================================================== +//! +//! \file applier.hpp +//! +//! \date Dec 22 2012 +//! +//! \author Arne Morten Kvarving / SINTEF +//! +//! \brief Helper class for abstracting differences between inverse operators +//! and preconditioners in DUNE +//! +//============================================================================== +#pragma once + +namespace Opm { + namespace Elasticity { + + template +struct OperatorApplier +{ + OperatorApplier(T& t) : A(t) + { + } + + void apply(Vector& v, Vector& d); + + T& A; +}; + +typedef OperatorApplier > InverseApplier; +typedef OperatorApplier > PreApplier; + + template<> +void InverseApplier::apply(Vector& v, Vector& d) +{ + Dune::InverseOperatorResult r; + A.apply(v, d, r); +} + + template<> +void PreApplier::apply(Vector& v, Vector& d) +{ + A.apply(v, d); +} + +} +} diff --git a/dune/elasticity/asmhandler.hpp b/dune/elasticity/asmhandler.hpp index 4f29be0..b37842e 100644 --- a/dune/elasticity/asmhandler.hpp +++ b/dune/elasticity/asmhandler.hpp @@ -15,6 +15,7 @@ #include #include #include +#include #include #include @@ -139,6 +140,13 @@ class ASMHandler { //! \brief Print the current load vector void printLoadVector() const; + + //! \brief Access current adjacency pattern + //! \details Can be used to add extra entries, such as other blocks + AdjacencyPattern& getAdjacencyPattern() + { + return adjacencyPattern; + } protected: //! \brief Resolve chained MPCs void resolveMPCChains() diff --git a/dune/elasticity/asmhandler_impl.hpp b/dune/elasticity/asmhandler_impl.hpp index 7ca9068..664ded4 100644 --- a/dune/elasticity/asmhandler_impl.hpp +++ b/dune/elasticity/asmhandler_impl.hpp @@ -21,14 +21,16 @@ void ASMHandler::initForAssembly() A.setBuildMode(Matrix::random); A.endrowsizes(); - MatrixOps::fromAdjacency(A,adjacencyPattern,maxeqn,maxeqn); - b.resize(maxeqn); + MatrixOps::fromAdjacency(A,adjacencyPattern, + adjacencyPattern.size(),adjacencyPattern.size()); + b.resize(adjacencyPattern.size()); b = 0; + adjacencyPattern.clear(); // print some information - std::cout << "Number of nodes: " << gv.size(dim) << std::endl; - std::cout << "Number of elements: " << gv.size(0) << std::endl; - std::cout << "Number of constraints: " << mpcs.size() << std::endl; + std::cout << "\tNumber of nodes: " << gv.size(dim) << std::endl; + std::cout << "\tNumber of elements: " << gv.size(0) << std::endl; + std::cout << "\tNumber of constraints: " << mpcs.size() << std::endl; int fixedDofs=0; for (fixIt it = fixedNodes.begin(); it != fixedNodes.end(); ++it) { if (it->second.first & X) @@ -38,7 +40,7 @@ void ASMHandler::initForAssembly() if (it->second.first & Z) fixedDofs++; } - std::cout << "Number of fixed dofs: " << fixedDofs << std::endl; + std::cout << "\tNumber of fixed dofs: " << fixedDofs << std::endl; } template @@ -339,7 +341,7 @@ void ASMHandler::preprocess() } } } - std::cout << "number of equations: " << maxeqn << std::endl; + std::cout << "\tnumber of equations: " << maxeqn << std::endl; } template @@ -370,12 +372,15 @@ void ASMHandler::nodeAdjacency(const LeafIterator& it, void ASMHandler::determineAdjacencyPattern() { adjacencyPattern.resize(maxeqn); + std::cout << "\tsetting up sparsity pattern..." << std::endl; + LoggerHelper help(gv.size(0), 5, 50000); const LeafIndexSet& set = gv.leafView().indexSet(); LeafIterator itend = gv.leafView().template end<0>(); // iterate over cells - for (LeafIterator it = gv.leafView().template begin<0>(); it != itend; ++it) { + int cell=0; + for (LeafIterator it = gv.leafView().template begin<0>(); it != itend; ++it, ++cell) { Dune::GeometryType gt = it->type(); const Dune::template GenericReferenceElement& ref = Dune::GenericReferenceElements::general(gt); @@ -395,5 +400,6 @@ void ASMHandler::determineAdjacencyPattern() nodeAdjacency(it,vertexsize,meqn[indexi*dim+k]); } } + help.log(cell, "\t\t... still processing ... cell "); } } diff --git a/dune/elasticity/elasticity_upscale.hpp b/dune/elasticity/elasticity_upscale.hpp index 79f3124..7320e5d 100644 --- a/dune/elasticity/elasticity_upscale.hpp +++ b/dune/elasticity/elasticity_upscale.hpp @@ -17,30 +17,101 @@ #include #include #include -#if HAVE_SUPERLU -#include -#endif #include -#include +#include #include #include #include #include #include +#include #include #include +#include +#include +#include +#include +#include + +#include namespace Opm { namespace Elasticity { -//! \brief An enumeration of available linear solvers + +//! \brief An enumeration of available linear solver classes enum Solver { - SLU, - CG + DIRECT, + ITERATIVE }; -typedef Dune::SeqOverlappingSchwarz OverlappingSchwarz; +//! \brief An enumeration of the available preconditioners +enum Preconditioner { + SCHURAMG, + SCHURDIAG +}; + +struct LinSolParams { + //! \brief The linear solver to employ + Opm::Elasticity::Solver type; + + //! \brief Number of iterations in GMRES before restart + int restart; + + //! \brief Max number of iterations + int maxit; + + //! \brief The tolerance for the iterative linear solver + double tol; + + //! \brief Use MINRES instead of GMRES (and thus symmetric preconditioning) + bool symmetric; + + //! \brief Use a Uzawa approach + bool uzawa; + + //! \brief The number of pre/post steps in the AMG + int steps[2]; + + //! \brief Coarsening target in the AMG + int coarsen_target; + + //! \brief Number of cells in z to collapse in each cell + int zcells; + + //! \brief Preconditioner for mortar block + Opm::Elasticity::Preconditioner mortarpre; + + //! \brief Parse command line parameters + //! \param[in] param The parameter group to parse + void parse(Opm::parameter::ParameterGroup& param) + { + std::string solver = param.getDefault("linsolver_type","iterative"); + if (solver == "iterative") + type = Opm::Elasticity::ITERATIVE; + else + type = Opm::Elasticity::DIRECT; + restart = param.getDefault("linsolver_restart", 1000); + tol = param.getDefault("ltol",1.e-8); + maxit = param.getDefault("linsolver_maxit", 10000); + steps[0] = param.getDefault("linsolver_presteps", 2); + steps[1] = param.getDefault("linsolver_poststeps", 2); + coarsen_target = param.getDefault("linsolver_coarsen", 5000); + symmetric = param.getDefault("linsolver_symmetric", true); + solver = param.getDefault("linsolver_mortarpre","schuramg"); + if (solver == "schuramg") + mortarpre = Opm::Elasticity::SCHURAMG; + else + mortarpre = Opm::Elasticity::SCHURDIAG; + uzawa = param.getDefault("linsolver_uzawa", false); + + zcells = param.getDefault("linsolver_zcells", 2); + + if (symmetric) + steps[1] = steps[0]; + } +}; //! \brief The main driver class template @@ -81,22 +152,24 @@ class ElasticityUpscale //! \param[in] tol_ The tolerance to use when deciding whether or not a coordinate falls on a plane/line/point. \sa tol //! \param[in] Escale_ A scale value for E-moduluses to avoid numerical issues //! \param[in] file The eclipse grid file - //! \param[in] rocklist If true, file is a rocklist + //! \param[in] rocklist If not blank, file is a rocklist + //! \param[in] verbose If true, give verbose output ElasticityUpscale(const GridType& gv_, ctype tol_, ctype Escale_, - const std::string& file, const std::string& rocklist) - : gv(gv_), tol(tol_), A(gv_), E(gv_), Escale(Escale_) + const std::string& file, const std::string& rocklist, + bool verbose_) + : A(gv_), gv(gv_), tol(tol_), Escale(Escale_), E(gv_), verbose(verbose_) { if (rocklist.empty()) loadMaterialsFromGrid(file); else loadMaterialsFromRocklist(file,rocklist); -#if HAVE_SUPERLU - slu = 0; -#endif - cgsolver = 0; + solver = 0; op = 0; - ilu = 0; - ovl = 0; + mpre = 0; + upre = 0; + op2 = 0; + meval = 0; + lpre = 0; } //! \brief The destructor @@ -111,13 +184,13 @@ class ElasticityUpscale it != itend; ++it) delete *it; -#if HAVE_SUPERLU - delete slu; -#endif - delete cgsolver; - delete ilu; + delete solver; delete op; - delete ovl; + delete op2; + delete meval; + delete upre; + delete lpre; + delete mpre; } //! \brief Find boundary coordinates @@ -137,15 +210,6 @@ class ElasticityUpscale //! \param[in] max The maximum coordinates of the grid void periodicBCs(const double* min, const double* max); - //! \brief Establish periodic boundaries using the LLM approach - //! \param[in] min The minimum coordinates of the grid - //! \param[in] max The maximum coordinates of the grid - //! \param[in] n1 The number of elements on the interface grid in the X direction - //! \param[in] n2 The number of elements on the interface grid in the Y direction - void periodicBCsLLM(const double* min, - const double* max, - int n1, int n2); - //! \brief Establish periodic boundaries using the mortar approach //! \param[in] min The minimum coordinates of the grid //! \param[in] max The maximum coordinates of the grid @@ -176,15 +240,16 @@ class ElasticityUpscale const Vector& u, int loadcase); //! \brief Solve Au = b for u - //! \param[in] solver The linear equation solver to employ - //! \param[in] tol The tolerance for iterative solvers - void solve(Solver solver, double tol, int loadcase); + //! \param[in] loadcase The load case to solve + void solve(int loadcase); + + //! \param[in] params The linear solver parameters + void setupSolvers(const LinSolParams& params); - void setupSolvers(Solver solver); private: //! \brief An iterator over grid vertices typedef typename GridType::LeafGridView::template Codim::Iterator LeafVertexIterator; - + //! \brief A reference to our grid const GridType& gv; @@ -271,26 +336,29 @@ class ElasticityUpscale //! \brief A point grid typedef std::map SlaveGrid; - //! \brief Find the B matrix associated with a particular slave grid (LLM) - //! \param[in] slave The slave grid to extract - //! \returns The B matrix associated with the given slave grid - Matrix findBMatrixLLM(const SlaveGrid& slave); - - //! \brief Find the L matrix associated with a particular slave grid (LLM) - //! \param[in] slave The slave grid we want to couple - //! \param[in] master The interface grid we want to couple to - //! \returns The L matrix describing the coupling between slave and master - Matrix findLMatrixLLM(const SlaveGrid& slave, const BoundaryGrid& master); - - //! \brief Find the L matrix associated with mortar couplings + //! \brief Add a block to the B matrix associated with mortar couplings //! \param[in] b1 The primal boundary to match //! \param[in] interface The interface/multiplier grid //! \param[in] dir The normal direction on the boundary (0/1) //! \param[in] n1 The multipler order in the first direction //! \param[in] n2 The multipler order in the second direction - Matrix findLMatrixMortar(const BoundaryGrid& b1, - const BoundaryGrid& interface, int dir, - int n1, int n2); + //! \param[in] colofs The column offset (multiplier number) + //! \returns Number of multipliers DOFs added + int addBBlockMortar(const BoundaryGrid& b1, + const BoundaryGrid& interface, int dir, + int n1, int n2, int colofs); + + //! \brief Assemble part of the B block associated with mortar couplings + //! \param[in] b1 The primal boundary to match + //! \param[in] interface The interface/multiplier grid + //! \param[in] dir The normal direction on the boundary (0/1) + //! \param[in] n1 The multipler order in the first direction + //! \param[in] n2 The multipler order in the second direction + //! \param[in] colofs The column offset (first multiplier unknown) + //! \param[in] alpha Scaling for matrix elements + void assembleBBlockMortar(const BoundaryGrid& b1, + const BoundaryGrid& interface, int dir, + int n1, int n2, int colofs, double alpha=1.0); //! \brief This function loads and maps materials to active grid cells //! \param[in] file The eclipse grid to read materials from @@ -302,9 +370,12 @@ class ElasticityUpscale void loadMaterialsFromRocklist(const std::string& file, const std::string& rocklist); - //! \brief Setup an overlapping schwarz preconditioner with one - //! subdomain per pillar. - void setupPreconditioner(); + //! \brief Setup AMG preconditioner + //! \param[in] pre The number of pre-smoothing steps + //! \param[in] post The number of post-smoothing steps + //! \param[in] target The coarsening target + //! \param[in] zcells The wanted number of cells to collapse in z per level + void setupAMG(int pre, int post, int target, int zcells); //! \brief Master grids std::vector master; @@ -312,30 +383,63 @@ class ElasticityUpscale //! \brief Slave point grids std::vector< std::vector > slave; - //! \brief Vector of matrices holding boolean coupling matrices (LLM) - std::vector B; - - //! \brief Vector of matrices holding Lagrangian multipliers - std::vector L; + //! \brief Lagrangian multiplier block + AdjacencyPattern Bpatt; + Matrix B; -#if HAVE_SUPERLU - //! \brief SuperLU solver - Dune::SuperLU* slu; -#endif - //! \brief CG solver - Dune::CGSolver* cgsolver; - //! \brief The linear operator - Dune::MatrixAdapter* op; + Matrix P; //!< Preconditioner for multiplier block - //! \brief ILU preconditioner - Dune::SeqILU0* ilu; - //! \brief Overlapping Schwarz preconditioner - OverlappingSchwarz* ovl; + //! \brief Linear solver + Dune::InverseOperator* solver; + + //! \brief The smoother used in the AMG + typedef Dune::SeqSSOR Smoother; + + //! \brief The coupling metric used in the AMG + typedef Dune::Amg::RowSum CouplingMetric; + + //! \brief The coupling criterion used in the AMG + typedef Dune::Amg::SymmetricCriterion CritBase; + + //! \brief The coarsening criterion used in the AMG + typedef Dune::Amg::CoarsenCriterion Criterion; + + //! \brief A linear operator + typedef Dune::MatrixAdapter Operator; + + //! \brief A preconditioner for an elasticity operator + typedef Dune::Amg::AMG ElasticityAMG; + + //! \brief Matrix adaptor for the elasticity block + Operator* op; + + //! \brief Preconditioner for multiplier block + typedef MortarBlockEvaluator > SchurPreconditioner; + + //! \brief Evaluator for multiplier block + typedef MortarBlockEvaluator > SchurEvaluator; + + //! \brief Outer evaluator, used with uzawa + SchurEvaluator* op2; + + //! \brief The preconditioner for the elasticity operator + ElasticityAMG* upre; + + //! \brief The preconditioner for the multiplier block (used with uzawa) + SeqLU* lpre; + + //! \brief Preconditioner for the Mortar system + MortarSchurPre* mpre; + + //! \brief Evaluator for the Mortar system + MortarEvaluator* meval; //! \brief Elasticity helper class Elasticity E; -}; + //! \brief Verbose output + bool verbose; +}; #include "elasticity_upscale_impl.hpp" } diff --git a/dune/elasticity/elasticity_upscale_impl.hpp b/dune/elasticity/elasticity_upscale_impl.hpp index 2d96ebc..4338965 100644 --- a/dune/elasticity/elasticity_upscale_impl.hpp +++ b/dune/elasticity/elasticity_upscale_impl.hpp @@ -45,13 +45,7 @@ BoundaryGrid ElasticityUpscale::extractMasterFace(Direction dir, {2,3,6,7}, {4,5,6,7}}; const LeafIndexSet& set = gv.leafView().indexSet(); - const LeafVertexIterator itend = gv.leafView().template end(); - // make a mapper for codim dim entities in the leaf grid - Dune::LeafMultipleCodimMultipleGeomTypeMapper mapper(gv); - LeafVertexIterator start=gv.leafView().template begin(); - LeafIterator cellend = gv.leafView().template end<0>(); int c = 0; int i = log2(dir); BoundaryGrid result; @@ -59,33 +53,25 @@ BoundaryGrid ElasticityUpscale::extractMasterFace(Direction dir, // vertex. we then split this up into pillars for easy processing later std::map > nodeMap; for (LeafIterator cell = gv.leafView().template begin<0>(); - cell != cellend; ++cell, ++c) { + cell != gv.leafView().template end<0>(); ++cell, ++c) { std::vector verts; - int idx; + int idx=0; if (side == LEFT) idx = set.subIndex(*cell,V1[i][0],dim); else if (side == RIGHT) idx = set.subIndex(*cell,V2[i][0],dim); - LeafVertexIterator it=start; - for (; it != itend; ++it) { - if (mapper.map(*it) == idx) - break; - } - if (isOnPlane(dir,it->geometry().corner(0),coord)) { + Dune::FieldVector pos = gv.vertexPosition(idx); + if (isOnPlane(dir,pos,coord)) { for (int j=0;j<4;++j) { if (side == LEFT) idx = set.subIndex(*cell,V1[i][j],dim); if (side == RIGHT) idx = set.subIndex(*cell,V2[i][j],dim); - LeafVertexIterator it=start; - for (; it != itend; ++it) { - if (mapper.map(*it) == idx) - break; - } - if (!isOnPlane(dir,it->geometry().corner(0),coord)) + pos = gv.vertexPosition(idx); + if (!isOnPlane(dir,pos,coord)) continue; BoundaryGrid::Vertex v; - BoundaryGrid::extract(v,it->geometry().corner(0),i); + BoundaryGrid::extract(v,pos,i); v.i = idx; verts.push_back(v); } @@ -173,9 +159,10 @@ void ElasticityUpscale::addMPC(Direction dir, int slave, } template -void ElasticityUpscale::periodicPlane(Direction plane, Direction dir, - const std::vector& slave, - const BoundaryGrid& master) +void ElasticityUpscale::periodicPlane(Direction plane, + Direction dir, + const std::vector& slave, + const BoundaryGrid& master) { for (size_t i=0;i::periodicPlane(Direction plane, Direction dir, } } - template -Matrix ElasticityUpscale::findBMatrixLLM(const SlaveGrid& slave) -{ - std::vector< std::set > adj; - adj.resize(A.getEqns()); - int cols=0; - std::vector colmap; - cols=0; - for (std::map::const_iterator it = slave.begin(); - it != slave.end();++it) { - for (int k=0;kfirst)) - continue; - MPC* mpc = A.getMPC(it->first,k); - if (mpc) { - for (size_t n=0;ngetNoMaster();++n) { - int idx = A.getEquationForDof(mpc->getMaster(n).node, - mpc->getMaster(n).dof-1); - if (idx != -1) - adj[idx].insert(cols++); - } - } else { - int idx = A.getEquationForDof(it->first,k); - if (idx != -1) - adj[idx].insert(cols++); - } - } - } - Matrix B; - MatrixOps::fromAdjacency(B,adj,A.getEqns(),cols); - int col=0; - for (std::map::const_iterator it = slave.begin(); - it != slave.end();++it) { - for (int k=0;kfirst)) - continue; - MPC* mpc = A.getMPC(it->first,k); - if (mpc) { - for (size_t n=0;ngetNoMaster();++n) { - int idx = A.getEquationForDof(mpc->getMaster(n).node, - mpc->getMaster(n).dof-1); - if (idx != -1) - B[idx][col++] += 1; - } - } else { - int idx = A.getEquationForDof(it->first,k); - if (idx != -1) - B[idx][col++] += 1; - } - } - } - - return B; -} - - template -Matrix ElasticityUpscale::findLMatrixLLM(const SlaveGrid& slave, - const BoundaryGrid& master) -{ - int nbeqn=0; - std::vector dofmap(master.totalNodes()*dim,-1); - int col=0; - for (size_t i=0;ifirst)) - nbeqn += 3; - } - std::vector< std::set > adj; - adj.resize(nbeqn); - - int row=0; - for (SlaveGrid::const_iterator it = slave.begin(); - it != slave.end();++it) { - if (A.isFixed(it->first)) - continue; - for (int k=0;ksecond.i == -1) { - for (int i=0;i<4;++i) { - int idx = dofmap[it->second.q->v[i].i*dim+k]; - if (idx > -1) - adj[row].insert(idx); - } - } - else { - int idx = dofmap[it->second.i*dim+k]; - if (idx > -1) - adj[row].insert(idx); - } - row++; - } - } - Matrix L; - MatrixOps::fromAdjacency(L,adj,nbeqn,col); - row = 0; - for (SlaveGrid::const_iterator it = slave.begin(); - it != slave.end();++it) { - if (A.isFixed(it->first)) - continue; - for (int k=0;ksecond.i == -1) { - std::vector N = it->second.q->evalBasis(it->second.c[0], - it->second.c[1]); - for (int i=0;i<4;++i) { - int idx = dofmap[it->second.q->v[i].i*dim+k]; - if (idx > -1) - L[row][idx] -= N[i]; - } - } - else { - int idx = dofmap[it->second.i*dim+k]; - if (idx > -1) - L[row][idx] -= 1; - } - row++; - } - } - - return L; -} - +//! \brief Static helper to renumber a Q4 mesh to a P_1/P_N lag mesh. +//! \param[in] b The boundary grid describing the Q4 mesh +//! \param[in] n1 Number of DOFS in the first direction for each element +//! \param[in] n2 Number of DOFS in the first direction for each element +//! \param[out] totalDOFs The total number of free DOFs static std::vector< std::vector > renumber(const BoundaryGrid& b, - int n1, int n2) + int n1, int n2, + int& totalDOFs) { std::vector > nodes; nodes.resize(b.size()); // loop over elements - int ofs = 0; + totalDOFs = 0; + + // fix lower left multiplicator. + // will be "transfered" to all corners through periodic conditions + nodes[0].push_back(-1); for (size_t e=0; e < b.size(); ++e) { // first direction major ordered nodes within each element for (int i2=0; i2 < n2; ++i2) { if (e != 0) nodes[e].push_back(nodes[e-1][i2*n1+n1-1]); - for (int i1=(e==0?0:1); i1 < n1; ++i1) - nodes[e].push_back(ofs++); + + int start = (e==0 && i2 != 0)?0:1; + + // slave the buggers + if (i2 == n2-1 && n2 > 2) { + for (int i1=(e==0?0:1); i1 < n1; ++i1) { + nodes[e].push_back(nodes[e][i1]); + } + } else { + for (int i1=start; i1 < n1; ++i1) { + if (e == b.size()-1) + nodes[e].push_back(nodes[0][i2*n1]); + else + nodes[e].push_back(totalDOFs++); + } + } } } @@ -340,44 +219,43 @@ static std::vector< std::vector > renumber(const BoundaryGrid& b, } template -Matrix ElasticityUpscale::findLMatrixMortar(const BoundaryGrid& b1, - const BoundaryGrid& interface, - int dir, int n1, int n2) +int ElasticityUpscale::addBBlockMortar(const BoundaryGrid& b1, + const BoundaryGrid& interface, + int dir, int n1, int n2, + int colofs) { - std::vector< std::set > adj; - adj.resize(A.getEqns()); - -#define QINDEX(p,q,dir) (q) - - // get a set of P1 shape functions for the displacements - P1ShapeFunctionSet ubasis = - P1ShapeFunctionSet::instance(); - - // get a set of PN shape functions for the multipliers - PNShapeFunctionSet<2> lbasis(n1+1, n2+1); - - std::vector > lnodes = renumber(interface,n1,n2); - int totalNodes = lnodes.back().back()+1; + // renumber the linear grid to the real multiplier grid + int totalEqns; + std::vector > lnodes = renumber(interface, n1+1, + n2+1, totalEqns); + if (Bpatt.empty()) + Bpatt.resize(A.getEqns()); // process pillar by pillar for (size_t p=0;pgetNoMaster();++n) { int dof = A.getEquationForDof(mpc->getMaster(n).node,d); if (dof > -1) { - for (size_t j=0;j -1) + Bpatt[dof].insert(indexj+colofs); + } } } } else { - int dof = A.getEquationForDof(b1.getQuad(p,QINDEX(p,q,dir)).v[i].i,d); + int dof = A.getEquationForDof(b1.getQuad(p,q).v[i].i,d); if (dof > -1) { - for (size_t j=0;j -1) + Bpatt[dof].insert(indexj+colofs); + } } } } @@ -385,9 +263,27 @@ Matrix ElasticityUpscale::findLMatrixMortar(const BoundaryGrid& b1, } } - Matrix B; - MatrixOps::fromAdjacency(B,adj,A.getEqns(),3*totalNodes); + return 3*totalEqns; +} + template +void ElasticityUpscale::assembleBBlockMortar(const BoundaryGrid& b1, + const BoundaryGrid& interface, + int dir, int n1, + int n2, int colofs, + double alpha) +{ + // get a set of P1 shape functions for the displacements + P1ShapeFunctionSet ubasis = + P1ShapeFunctionSet::instance(); + + // get a set of PN shape functions for the multipliers + PNShapeFunctionSet<2> lbasis(n1+1, n2+1); + + // renumber the linear grid to the real multiplier grid + int totalEqns; + std::vector > lnodes = renumber(interface, n1+1, + n2+1, totalEqns); // get a reference element Dune::GeometryType gt; gt.makeCube(2); @@ -400,6 +296,7 @@ Matrix ElasticityUpscale::findLMatrixMortar(const BoundaryGrid& b1, // do the assembly loop typename Dune::QuadratureRule::const_iterator r; Dune::DynamicMatrix E(ubasis.n,(n1+1)*(n2+1),0.0); + LoggerHelper help(interface.size(), 5, 1000); for (size_t p=0;p lg(qi); @@ -409,8 +306,8 @@ Matrix ElasticityUpscale::findLMatrixMortar(const BoundaryGrid& b1, E = 0; for (r = rule.begin(); r != rule.end();++r) { ctype detJ = hex.integrationElement(r->position()); - if (detJ < 1.e-4 && r == rule.begin()) - std::cerr << "warning: interface cell (close to) degenerated, |J|=" << detJ << std::endl; + if (detJ < 0) + assert(0); typename HexGeometry<2,2,GridType>::LocalCoordinate loc = lg.local(hex.global(r->position())); @@ -422,6 +319,7 @@ Matrix ElasticityUpscale::findLMatrixMortar(const BoundaryGrid& b1, } } } + // and assemble element contributions for (int d=0;d<3;++d) { for (int i=0;i<4;++i) { @@ -432,7 +330,8 @@ Matrix ElasticityUpscale::findLMatrixMortar(const BoundaryGrid& b1, if (indexi > -1) { for (size_t j=0;j -1) + B[indexi][indexj+colofs] += alpha*E[i][j]; } } } @@ -441,16 +340,16 @@ Matrix ElasticityUpscale::findLMatrixMortar(const BoundaryGrid& b1, if (indexi > -1) { for (size_t j=0;j -1) + B[indexi][indexj+colofs] += alpha*E[i][j]; } } } } } } + help.log(p, "\t\t\t... still processing ... pillar "); } - - return B; } template @@ -557,9 +456,11 @@ void ElasticityUpscale::assemble(int loadcase, bool matrix) int m=0; Dune::FieldMatrix* KP=0; if (matrix) { - A.getOperator() = 0; KP = &K; + A.getOperator() = 0; } + + LoggerHelper help(gv.size(0), 5, 50000); for (LeafIterator it = gv.leafView().template begin<0>(); it != itend; ++it) { materials[m++]->getConstitutiveMatrix(C); // determine geometry type of the current element and get the matching reference element @@ -578,7 +479,7 @@ void ElasticityUpscale::assemble(int loadcase, bool matrix) it->geometry().jacobianInverseTransposed(r->position()); ctype detJ = it->geometry().integrationElement(r->position()); - if (detJ <= 1.e-4) { + if (detJ <= 1.e-5 && verbose) { std::cout << "cell " << m << " is (close to) degenerated, detJ " << detJ << std::endl; double zdiff=0.0; for (int i=0;i<4;++i) @@ -604,6 +505,7 @@ void ElasticityUpscale::assemble(int loadcase, bool matrix) } A.addElement(KP,EP,it,(loadcase > -1)?&b[loadcase]:NULL); + help.log(m, "\t\t... still processing ... cell "); } } @@ -643,8 +545,6 @@ void ElasticityUpscale::averageStress(Dune::FieldVector& s it->geometry().jacobianInverseTransposed(r->position()); ctype detJ = it->geometry().integrationElement(r->position()); - if (detJ <= 0) // wtf ? - continue; volume += detJ*r->weight(); @@ -720,12 +620,14 @@ void ElasticityUpscale::loadMaterialsFromGrid(const std::string& file) MaterialMap::iterator it; if ((it = cache.find(std::make_pair(Emod[k],Poiss[k]))) != cache.end()) { + assert(gv.cellVolume(i) > 0); volume[it->second] += gv.cellVolume(i); materials.push_back(it->second); } else { Material* mat = new Isotropic(j++,Emod[k],Poiss[k]); cache.insert(std::make_pair(std::make_pair(Emod[k],Poiss[k]),mat)); + assert(gv.cellVolume(i) > 0); volume.insert(std::make_pair(mat,gv.cellVolume(i))); materials.push_back(mat); } @@ -782,7 +684,7 @@ void ElasticityUpscale::loadMaterialsFromRocklist(const std::string& f std::vector cells = gv.globalCell(); for (size_t i=0;i= cache.size()) { + if (satnum[k]-1 >= (int)cache.size()) { std::cerr << "Material " << satnum[k] << " referenced but not available. Check your rocklist." << std::endl; exit(1); } @@ -857,254 +759,215 @@ void ElasticityUpscale::periodicBCsMortar(const double* min, // 2. establishes strong couplings (MPC) on top and bottom // 3. extracts and establishes a quad grid for the left/right/front/back sides // 4. establishes grids for the dual dofs - // 5. calculates the coupling matrix L1 between the left/right sides - // 6. calculates the coupling matrix L2 between the front/back sides + // 5. calculates the coupling matrix between the left/right sides + // 6. calculates the coupling matrix between the front/back sides + // + // The Mortar linear system is of the form + // [A B] [u] [b] + // [B' 0] [l] = [0] // step 1 fixCorners(min,max); // step 2 + std::cout << "\textracting nodes on top face..." << std::endl; slave.push_back(extractFace(Z,max[2])); + std::cout << "\treconstructing bottom face..." << std::endl; BoundaryGrid bottom = extractMasterFace(Z,min[2]); + std::cout << "\testablishing couplings on top/bottom..." << std::endl; periodicPlane(Z,XYZ,slave[0],bottom); + std::cout << "\tinitializing matrix..." << std::endl; A.initForAssembly(); // step 3 - master.push_back(extractMasterFace(X,min[0],LEFT,true)); - master.push_back(extractMasterFace(X,max[0],RIGHT,true)); - master.push_back(extractMasterFace(Y,min[1],LEFT,true)); - master.push_back(extractMasterFace(Y,max[1],RIGHT,true)); + std::cout << "\treconstructing left face..." << std::endl; + master.push_back(extractMasterFace(X, min[0], LEFT, true)); + std::cout << "\treconstructing right face..." << std::endl; + master.push_back(extractMasterFace(X, max[0], RIGHT, true)); + std::cout << "\treconstructing front face..." << std::endl; + master.push_back(extractMasterFace(Y, min[1], LEFT, true)); + std::cout << "\treconstructing back face..." << std::endl; + master.push_back(extractMasterFace(Y, max[1], RIGHT, true)); - std::cout << "Xsides " << master[0].size() << " " << master[1].size() << std::endl - << "Ysides " << master[2].size() << " " << master[3].size() << std::endl - << "Zmaster " << bottom.size() << std::endl; - std::cout << "Establish YZ multiplier grid with " << n2 << "x1" << " elements" << std::endl; + std::cout << "\testablished YZ multiplier grid with " << n2 << "x1" << " elements" << std::endl; // step 4 BoundaryGrid::FaceCoord fmin,fmax; fmin[0] = min[1]; fmin[1] = min[2]; fmax[0] = max[1]; fmax[1] = max[2]; - BoundaryGrid lambdax = BoundaryGrid::uniform(fmin,fmax,n2,1,true); + BoundaryGrid lambdax = BoundaryGrid::uniform(fmin, fmax, n2, 1, true); fmin[0] = min[0]; fmin[1] = min[2]; fmax[0] = max[0]; fmax[1] = max[2]; - std::cout << "Establish XZ multiplier grid with " << n1 << "x1" << " elements" << std::endl; - BoundaryGrid lambday = BoundaryGrid::uniform(fmin,fmax,n1,1,true); + std::cout << "\testablished XZ multiplier grid with " << n1 << "x1" << " elements" << std::endl; + BoundaryGrid lambday = BoundaryGrid::uniform(fmin, fmax, n1, 1, true); + + addBBlockMortar(master[0], lambdax, 0, 1, p2, 0); + int eqns = addBBlockMortar(master[1], lambdax, 0, 1, p2, 0); + addBBlockMortar(master[2], lambday, 1, 1, p1, eqns); + int eqns2 = addBBlockMortar(master[3], lambday, 1, 1, p1, eqns); + + MatrixOps::fromAdjacency(B, Bpatt, A.getEqns(), eqns+eqns2); + Bpatt.clear(); // step 5 - Matrix L1 = findLMatrixMortar(master[0],lambdax,0, p2, 1); - Matrix L2 = findLMatrixMortar(master[1],lambdax,0, p2, 1); - L.push_back(MatrixOps::Axpy(L1,L2,-1)); + std::cout << "\tassembling YZ mortar matrix..." << std::endl; + assembleBBlockMortar(master[0], lambdax, 0, 1, p2, 0); + assembleBBlockMortar(master[1], lambdax, 0, 1, p2, 0, -1.0); // step 6 - Matrix L3 = findLMatrixMortar(master[2],lambday,1, p1, 1); - Matrix L4 = findLMatrixMortar(master[3],lambday,1, p1, 1); - L.push_back(MatrixOps::Axpy(L3,L4,-1)); + std::cout << "\tassembling XZ mortar matrix..." << std::endl; + assembleBBlockMortar(master[2], lambday, 1, 1, p1, eqns); + assembleBBlockMortar(master[3], lambday, 1, 1, p1, eqns, -1.0); + + master.clear(); + slave.clear(); } template -void ElasticityUpscale::periodicBCsLLM(const double* min, - const double* max, - int n1, int n2) +void ElasticityUpscale::setupAMG(int pre, int post, + int target, int zcells) { - // this method - // 1. fixes the primal corner dofs - // 2. fixes the primal dofs on the skeleton - // 3. establishes strong couplings (MPC) on top and bottom - // 4. extracts a point grid for the left/right/front/back sides, - // and establishes a uniform interface grid in each direction - // 5. calculates the coupling matrices B1-4 and L1-4 + Criterion crit; + ElasticityAMG::SmootherArgs args; + args.relaxationFactor = 1.0; + crit.setCoarsenTarget(target); + crit.setGamma(1); + crit.setNoPreSmoothSteps(pre); + crit.setNoPostSmoothSteps(post); + crit.setDefaultValuesIsotropic(3, zcells); - // step 1 - fixCorners(min,max); + std::cout << "\t collapsing 2x2x" << zcells << " cells per level" << std::endl; + op = new Operator(A.getOperator()); + upre = new ElasticityAMG(*op, crit, args); - // step 2 - fixLine(X,min[1],min[2]); - fixLine(X,max[1],min[2]); - fixLine(X,min[1],max[2]); - fixLine(X,max[1],max[2]); - - fixLine(Y,min[0],min[2]); - fixLine(Y,max[0],min[2]); - fixLine(Y,min[0],max[2]); - fixLine(Y,max[0],max[2]); - - fixLine(Z,min[0],min[1]); - fixLine(Z,max[0],min[1]); - fixLine(Z,min[0],max[1]); - fixLine(Z,max[0],max[1]); - - // step 3 - std::vector Zs = extractFace(Z,max[2]); - BoundaryGrid Zm = extractMasterFace(Z,min[2]); - periodicPlane(Z,XYZ,Zs,Zm); - A.initForAssembly(); - - // step 4 - slave.push_back(extractFace(X,min[0])); - slave.push_back(extractFace(X,max[0])); - slave.push_back(extractFace(Y,min[1])); - slave.push_back(extractFace(Y,max[1])); - - Dune::LeafMultipleCodimMultipleGeomTypeMapper mapper(gv); - - BoundaryGrid::FaceCoord fmin,fmax; - // YZ plane - fmin[0] = min[1]; fmin[1] = min[2]; - fmax[0] = max[1]; fmax[1] = max[2]; - std::cout << "Establish YZ interface grid with " << n2 << "x1" << " elements" << std::endl; - master.push_back(BoundaryGrid::uniform(fmin,fmax,n2,1)); - // XZ plane - fmin[0] = min[0]; fmin[1] = min[2]; - fmax[0] = max[0]; fmax[1] = max[2]; - std::cout << "Establish XZ interface grid with " << n1 << "x1" << " elements" << std::endl; - master.push_back(BoundaryGrid::uniform(fmin,fmax,n1,1)); - - // step 5 - std::map m; - // find matching coefficients - for (size_t i=0;iaddContext("Apre"); + amg->setContext("Apre"); + Vector x,y; + // this is done here to make sure we are in a single-threaded section + // will have to be redone when AMG is refactored upstream + amg->pre(x,y); + */ } - template -void ElasticityUpscale::setupPreconditioner() + template +void ElasticityUpscale::setupSolvers(const LinSolParams& params) { - // setup subdomain maps - typename OverlappingSchwarz::subdomain_vector rows; - const LeafIterator itend = gv.leafView().template end<0>(); - const LeafIndexSet& set = gv.leafView().indexSet(); - - rows.resize(gv.logicalCartesianSize()[0]*gv.logicalCartesianSize()[1]); - int cell=0, currdomain=0; - for (LeafIterator it = gv.leafView().template begin<0>(); it != itend; ++it, ++cell) { - if (cell / gv.logicalCartesianSize()[2] > 0 - && cell % gv.logicalCartesianSize()[2] == 0) - currdomain++; - for (size_t i=0;i<8;++i) { - int idx=set.subIndex(*it,i,dim); - for (int d=0;d<3;++d) { - MPC* mpc = A.getMPC(idx,d); - if (mpc) { - for (size_t n=0;ngetNoMaster();++n) { - int row = A.getEquationForDof(mpc->getMaster(n).node,d); - if (row > -1) - rows[currdomain].insert(row); + int siz = A.getOperator().N(); // system size + if (params.type == ITERATIVE) { + setupAMG(params.steps[0], params.steps[1], params.coarsen_target, + params.zcells); + + // Mortar in use + if (B.N()) { + siz += B.M(); + + // schur system: B'*diag(A)^-1*B + if (params.mortarpre == SCHURAMG) { + Vector v, v2, v3; + v.resize(B.N()); + v2.resize(B.N()); + v = 0; + v2 = 0; + Dune::DynamicMatrix T(B.M(), B.M()); + upre->pre(v, v); + std::cout << "\tBuilding preconditioner for multipliers..." << std::endl; + MortarBlockEvaluator > pre(*upre, B); + LoggerHelper help(B.M(), 10, 100); + for (size_t i=0; i < B.M(); ++i) { + v[i] = 1; + pre.apply(v, v2); + for (size_t j=0; j < B.M(); ++j) + T[j][i] = v2[j]; + + v[i] = 0; + help.log(i, "\t\t... still processing ... multiplier "); + } + upre->post(v); + P = MatrixOps::fromDense(T); + } else if (params.mortarpre == SCHURDIAG) { + Matrix D = MatrixOps::diagonal(A.getEqns()); + + // scale by row sums + size_t row=0; + for (Matrix::ConstRowIterator it = A.getOperator().begin(); + it != A.getOperator().end(); ++it, ++row) { + double alpha=0; + for (Matrix::ConstColIterator it2 = it->begin(); + it2 != it->end(); ++it2) { + if (it2.index() != row) + alpha += fabs(*it2); } + D[row][row] = 1.0/(A.getOperator()[row][row]/alpha); + } + + Matrix t1; + // t1 = Ad*B + Dune::matMultMat(t1, D, B); + // P = B'*t1 = B'*Ad*B + Dune::transposeMatMultMat(P, B, t1); + } + + if (params.uzawa) { + Dune::CGSolver* innersolver = + new Dune::CGSolver(*op, *upre, params.tol, + params.maxit, verbose?2:0); + op2 = new SchurEvaluator(*innersolver, B); + lpre = new SeqLU(P); + Dune::CGSolver* outersolver = + new Dune::CGSolver(*op2, *lpre, params.tol*10, + params.maxit, verbose?2:0); + solver = new UzawaSolver(innersolver, outersolver, B); + } else { + mpre = new MortarSchurPre(P, B, *upre, params.symmetric); + meval = new MortarEvaluator(A.getOperator(), B); + if (params.symmetric) { + solver = new Dune::MINRESSolver(*meval, *mpre, + params.tol, + params.maxit, + verbose?2:0); } else { - int row = A.getEquationForDof(idx,d); - if (row > -1) - rows[currdomain].insert(row); + solver = new Dune::RestartedGMResSolver(*meval, *mpre, + params.tol, + params.restart, + params.maxit, + verbose?2:0, true); } } + } else { + solver = new Dune::CGSolver(*op, *upre, params.tol, + params.maxit, verbose?2:0); } - } - ovl = new OverlappingSchwarz(A.getOperator(),rows,1,false); -} - - template -void ElasticityUpscale::setupSolvers(Solver solver) -{ - if (!B.empty() && A.getOperator().N() == A.getEqns()) { // LLM in use - // The LLM linear system is of the form - // [A B1 B2 B3 B4 0 0] [ u] [b] - // [B1' 0 0 0 0 L1 0] [ l_1] [0] - // [B2' 0 0 0 0 L2 0] [ l_2] [0] - // [B3' 0 0 0 0 0 L3] [ l_3] = [0] - // [B4' 0 0 0 0 0 L4] [ l_4] [0] - // [ 0 L1' L2' 0 0 0 0] [ub_1] [0] - // [ 0 0 0 L3' L4' 0 0] [ub_2] [0] - int r = A.getOperator().N(); - int c = A.getOperator().M(); - for (size_t i=0;i(A.getOperator(),false); + solver = new Dune::SuperLU(A.getOperator(),verbose); #else std::cerr << "SuperLU solver not enabled" << std::endl; exit(1); #endif - } else if (solver == CG) { - if (!op) - op = new Dune::MatrixAdapter(A.getOperator()); - if (!ovl && !ilu) - ilu = new Dune::SeqILU0(A.getOperator(),0.92); - if (!cgsolver) { - if (ovl) - cgsolver = new Dune::CGSolver(*op, *ovl, tol, 5000, 1); - else - cgsolver = new Dune::CGSolver(*op, *ilu, tol, 5000, 1); - } + siz = A.getOperator().N(); } + + for (int i=0;i<6;++i) + b[i].resize(siz); } template -void ElasticityUpscale::solve(Solver solver, double tol, int loadcase) +void ElasticityUpscale::solve(int loadcase) { - Dune::InverseOperatorResult r; + try { + Dune::InverseOperatorResult r; + u[loadcase].resize(b[loadcase].size(), false); + u[loadcase] = 0; - // initialize u to some arbitrary value - u[loadcase].resize(A.getOperator().N(), false); - u[loadcase] = 2.0; + solver->apply(u[loadcase], b[loadcase], r); -#if HAVE_SUPERLU - if (solver == SLU) { - slu->apply(u[loadcase], b[loadcase], r); + std::cout << "\tsolution norm: " << u[loadcase].two_norm() << std::endl; + } catch (Dune::ISTLError& e) { + std::cerr << "exception thrown " << e << std::endl; } - else -#endif - if (solver == CG) { - cgsolver->apply(u[loadcase], b[loadcase], r); - } - std::cout << "\t solution norm: " << u[loadcase].two_norm() << std::endl; } diff --git a/dune/elasticity/logutils.hpp b/dune/elasticity/logutils.hpp new file mode 100644 index 0000000..2dd2415 --- /dev/null +++ b/dune/elasticity/logutils.hpp @@ -0,0 +1,31 @@ +//============================================================================== +//! +//! \file logutils.hpp +//! +//! \date Nov 9 2011 +//! +//! \author Arne Morten Kvarving / SINTEF +//! +//! \brief Logging helper utilities +//! +//============================================================================== +#pragma once + +class LoggerHelper { + public: + LoggerHelper(int max_, int chunks, int minsize) + : per_chunk(max_/chunks), max(max_) + { + if (max < minsize) + per_chunk=-1; + } + + void log(int it, const std::string& prefix) + { + if (per_chunk > 0 && it > 0 && (it % per_chunk) == 0 && max-it > per_chunk) + std::cout << prefix << it << '/' << max << std::endl; + } + protected: + int per_chunk; //!< Will log for each per_chunk processed + int max; //!< Total number of its +}; diff --git a/dune/elasticity/matrixops.hpp b/dune/elasticity/matrixops.hpp index e2aefc1..73f7a82 100644 --- a/dune/elasticity/matrixops.hpp +++ b/dune/elasticity/matrixops.hpp @@ -36,6 +36,11 @@ class MatrixOps { static void fromAdjacency(Matrix& A, const AdjacencyPattern& adj, int rows, int cols); + //! \brief Create a sparse matrix from a dense matrix + //! \param[in] T The dense matrix + //! \returns The sparse matrix + static Matrix fromDense(const Dune::DynamicMatrix& T); + //! \brief Print a matrix to stdout //! \param[in] A The matrix to print static void print(const Matrix& A); @@ -61,6 +66,16 @@ class MatrixOps { //! \returns M = diag(A) static Matrix extractDiagonal(const Matrix& A); + //! \brief Returns a diagonal matrix + //! \param[in] N The dimension of the matrix + static Matrix diagonal(size_t N); + + //! \brief Extract a subblock of a matrix into a new matrix + //! \param[in] The matrix to extract from + //! \returns The subblock + static Matrix extractBlock(const Matrix& A, + size_t r0, size_t N, size_t c0, size_t M); + //! \brief Save a matrix as a dense asc file //! \param[in] A The matrix to save //! \param[in] file File name diff --git a/dune/elasticity/matrixops_impl.hpp b/dune/elasticity/matrixops_impl.hpp index 57533c2..83e3001 100644 --- a/dune/elasticity/matrixops_impl.hpp +++ b/dune/elasticity/matrixops_impl.hpp @@ -12,7 +12,7 @@ void MatrixOps::fromAdjacency(Matrix& A, const std::vector< std::set >& adj, int rows, int cols) -{ +{ size_t sum=0; for (size_t i=0;i >& adj A = 0; } +Matrix MatrixOps::fromDense(const Dune::DynamicMatrix& T) +{ + AdjacencyPattern a; + a.resize(T.N()); + for (size_t i=0; i < T.N(); ++i) { + for (size_t j=0; j < T.M(); ++j) { + if (fabs(T[i][j]) > 1.e-14) + a[i].insert(j); + } + } + + Matrix result; + fromAdjacency(result, a, T.N(), T.M()); + + for (Matrix::ConstIterator it = result.begin(); + it != result.end(); ++it) { + for (Matrix::ConstColIterator it2 = it->begin(); + it2 != it->end();++it2) + result[it.index()][it2.index()] = T[it.index()][it2.index()]; + } + + return result; +} + void MatrixOps::print(const Matrix& A) { for (Matrix::ConstIterator it = A.begin(); @@ -159,6 +183,18 @@ Matrix MatrixOps::extractDiagonal(const Matrix& A) return result; } +Matrix MatrixOps::diagonal(size_t N) +{ + AdjacencyPattern adj; + adj.resize(N); + for (size_t i=0;i= N+r0) + continue; + for (Matrix::ConstColIterator it2 = it->begin(); + it2 != it->end(); ++it2) { + if (it2.index() < c0 || it2.index() >= M+c0) + continue; + adj[it.index()-r0].insert(it2.index()-c0); + } + } + + Matrix result; + fromAdjacency(result, adj, N, M); + + // now insert elements from A + for (Matrix::ConstIterator it = A.begin(); + it != A.end(); ++it) { + if (it.index() < r0 || it.index() >= N+r0) + continue; + for (Matrix::ConstColIterator it2 = it->begin(); + it2 != it->end(); ++it2) { + if (it2.index() < c0 || it2.index() >= M+c0) + continue; + result[it.index()-r0][it2.index()-c0] = *it2; + } + } + + return result; +} diff --git a/dune/elasticity/mortar_evaluator.hpp b/dune/elasticity/mortar_evaluator.hpp new file mode 100644 index 0000000..c5713de --- /dev/null +++ b/dune/elasticity/mortar_evaluator.hpp @@ -0,0 +1,79 @@ +//============================================================================== +//! +//! \file mortar_evaluator.hpp +//! +//! \date Nov 15 2012 +//! +//! \author Arne Morten Kvarving / SINTEF +//! +//! \brief Linear operator for a Mortar block +//! +//============================================================================== +#pragma once + +#include +#include +#include + +namespace Opm { +namespace Elasticity { + +/*! This implements a operator evaluation for the + * schur mortar-block S = B^T*A^-1*B +!*/ + +class MortarEvaluator : public Dune::LinearOperator { + public: + // define the category + enum { + //! \brief The category the preconditioner is part of. + category=Dune::SolverCategory::sequential + }; + + //! \brief Constructor + //! \param[in] Ai Evaluator for A^-1 + //! \param[in] B The mortar coupling matrix + MortarEvaluator(const Matrix& A_, + const Matrix& B_) : + A(A_), B(B_) + { + } + + //! \brief Apply the multiplier block + //! \param[in] x The vector to apply the operator to + //! \param[out] y The result of the operator evaluation + void apply(const Vector& x, Vector& y) const + { + Vector lambda, l2; + + MortarUtils::extractBlock(lambda, x, B.M(), A.N()); + A.mv(x, y); + B.umv(lambda, y); + l2.resize(lambda.size()); + B.mtv(x, l2); + MortarUtils::injectBlock(y, l2, B.M(), A.N()); + } + + //! \brief Apply the multiplier block with an embedded axpy + //! \param[in] alpha The scalar to scale with + //! \param[in] x The vector to apply the operator to + //! \param[out] y The result of the operator evaluation + void applyscaleadd(field_type alpha, const Vector& x, Vector& y) const + { + Vector lambda, l2; + + A.usmv(alpha, x, y); + MortarUtils::extractBlock(lambda, x, B.M(), A.N()); + B.umv(lambda, y); + l2.resize(lambda.size()); + B.mtv(x, l2); + for (size_t i=A.N(); i < y.size(); ++i) + y[i] += alpha*l2[i-A.N()]; + } + protected: + const Matrix& A; //!< Reference to A matrix + const Matrix& B; //!< Reference to the mortar coupling matrix +}; + +} +} diff --git a/dune/elasticity/mortar_schur.hpp b/dune/elasticity/mortar_schur.hpp new file mode 100644 index 0000000..25b87f7 --- /dev/null +++ b/dune/elasticity/mortar_schur.hpp @@ -0,0 +1,73 @@ +//============================================================================== +//! +//! \file mortar_schur.hpp +//! +//! \date Nov 15 2012 +//! +//! \author Arne Morten Kvarving / SINTEF +//! +//! \brief Schur based linear operator for a Mortar block +//! +//============================================================================== +#pragma once + +#include +#include +#include + +namespace Opm { +namespace Elasticity { + +/*! This implements a operator evaluation for the + * schur mortar-block S = B^T*A^-1*B +!*/ + + template +class MortarBlockEvaluator : public Dune::LinearOperator { + public: + // define the category + enum { + //! \brief The category the preconditioner is part of. + category=Dune::SolverCategory::sequential + }; + + //! \brief Constructor + //! \param[in] Ai Solver or preconditioner for A^-1 + //! \param[in] B The mortar coupling matrix + MortarBlockEvaluator(T& Ai_, + const Matrix& B_) : + Ai(Ai_), B(B_), op(Ai) + { + } + + //! \brief Apply the multiplier block + //! \param[in] x The vector to apply the operator to + //! \param[out] y The result of the operator evaluation + void apply(const Vector& x, Vector& y) const + { + y = 0; + applyscaleadd(1.0, x, y); + } + + //! \brief Apply the multiplier block with an embedded axpy + //! \param[in] alpha The scalar to scale with + //! \param[in] x The vector to apply the operator to + //! \param[out] y The result of the operator evaluation + void applyscaleadd(field_type alpha, const Vector& x, Vector& y) const + { + Vector temp1(B.N()); + B.mv(x, temp1); + Vector temp(temp1.size()); + Dune::InverseOperatorResult res; + temp = 0; + op.apply(temp, temp1); + B.usmtv(alpha, temp, y); + } + protected: + T& Ai; //!< Reference to solver or evaluator for inverse operator + const Matrix& B; //!< Reference to the mortar coupling matrix + mutable OperatorApplier op; //!< Applier for the preconditioner / inverse solver +}; + +} +} diff --git a/dune/elasticity/mortar_schur_precond.hpp b/dune/elasticity/mortar_schur_precond.hpp new file mode 100644 index 0000000..c780422 --- /dev/null +++ b/dune/elasticity/mortar_schur_precond.hpp @@ -0,0 +1,123 @@ +//============================================================================== +//! +//! \file schur_precond.hpp +//! +//! \date Nov 15 2012 +//! +//! \author Arne Morten Kvarving / SINTEF +//! +//! \brief Schur based preconditioner for saddle point systems +//! +//============================================================================== +#pragma once + +#include +#include +#include +#include + +namespace Opm { +namespace Elasticity { + +/*! This implements a Schur-decomposition based preconditioner for the + * mortar-elasticity system + * [A B] + * [B' ] + * + * The preconditioner is + * [Apre B] + * [ P] + * Here Apre is some preconditioner for A and P some preconditioner for + * S = B^TA^-1B +!*/ + + template +class MortarSchurPre : public Dune::Preconditioner { + public: + // define the category + enum { + //! \brief The category the preconditioner is part of. + category=Dune::SolverCategory::sequential + }; + + //! \brief Constructor + //! \param[in] P The multiplier block with diagonal A approximation + //! \param[in] B The mortar coupling matrix + //! \param[in] Apre_ A preconfigured preconditioner for A + //! \param[in] symmetric If true, use symmetric preconditioning + MortarSchurPre(const Matrix& P_, const Matrix& B_, + PrecondElasticityBlock& Apre_, bool symmetric_=false) : + Apre(Apre_), B(B_), N(B.N()), M(B.M()), + Lpre(P_, false, false), symmetric(symmetric_) + { + } + + //! \brief Destructor + virtual ~MortarSchurPre() + { + } + + //! \brief Preprocess preconditioner + virtual void pre(Vector& x, Vector& b) + { + // extract displacement DOFs + Vector tempx, tempb; + MortarUtils::extractBlock(tempx, x, N); + MortarUtils::extractBlock(tempb, b, N); + Apre.pre(tempx, tempb); + MortarUtils::injectBlock(x, tempx, N); + MortarUtils::injectBlock(b, tempb, N); + } + + //! \brief Applies the preconditioner + //! \param[out] v The resulting vector + //! \param[in] d The vector to apply the preconditioner to + virtual void apply(Vector& v, const Vector& d) + { + // multiplier block (second row) + Vector temp; + MortarUtils::extractBlock(temp, d, M, N); + Dune::InverseOperatorResult r; + Vector temp2; + temp2.resize(temp.size(), false); + Lpre.apply(temp2, temp, r); + MortarUtils::injectBlock(v, temp2, M, N); + + // elasticity block (first row) + MortarUtils::extractBlock(temp, d, N); + if (!symmetric) + B.mmv(temp2, temp); + + temp2.resize(N, false); + temp2 = 0; + Apre.apply(temp2, temp); + MortarUtils::injectBlock(v, temp2, N); + } + + //! \brief Dummy post-process function + virtual void post(Vector& x) + { + Apre.post(x); + } + protected: + //! \brief The preconditioner for the elasticity operator + PrecondElasticityBlock& Apre; + + //! \brief The mortar coupling matrix + const Matrix& B; + + //! \brief Number of displacement DOFs + int N; + + //! \brief Number of multiplier DOFs + int M; + + //! \brief Linear solver for the multiplier block + Dune::SuperLU Lpre; + + //! \brief Whether or not to use a symmetric preconditioner + bool symmetric; +}; + +} +} diff --git a/dune/elasticity/mortar_utils.hpp b/dune/elasticity/mortar_utils.hpp new file mode 100644 index 0000000..1e1c552 --- /dev/null +++ b/dune/elasticity/mortar_utils.hpp @@ -0,0 +1,31 @@ +#pragma once + +namespace Opm { + namespace Elasticity { + +class MortarUtils { + public: + //! \brief Extract a range of indices from a vector + //! \param[out] x The vector with the extracted data + //! \param[in] y The vector + //! \param[in] len The number of indices to extract + //! \param[in] start The first index in the range + static void extractBlock(Vector& x, const Vector& y, int len, int start=0) + { + x.resize(len,false); + std::copy(y.begin()+start,y.begin()+len+start,x.begin()); + } + + //! \brief Inject a range of indices into a vector + //! \param[in/out] x The vector to inject into + //! \param[in] y The vector with the data to inject + //! \param[in] len The number of indices to inject + //! \param[in] start The first index in the range + static void injectBlock(Vector& x, const Vector& y, int len, int start=0) + { + std::copy(y.begin(),y.begin()+len,x.begin()+start); + } +}; + +} +} diff --git a/dune/elasticity/seqlu.hpp b/dune/elasticity/seqlu.hpp new file mode 100644 index 0000000..e85c5bf --- /dev/null +++ b/dune/elasticity/seqlu.hpp @@ -0,0 +1,77 @@ +//============================================================================== +//! +//! \file seqlu.hpp +//! +//! \date Dec 22 2012 +//! +//! \author Arne Morten Kvarving / SINTEF +//! +//! \brief SuperLU preconditioner interface +//! +//============================================================================== +#pragma once + +#if HAVE_SUPERLU +#include +#endif + +#if HAVE_SUPERLU + template +class SeqLU : public Dune::Preconditioner { + public: + //! \brief The matrix type the preconditioner is for. + typedef typename Dune::remove_const::type matrix_type; + //! \brief The domain type of the preconditioner. + typedef X domain_type; + //! \brief The range type of the preconditioner. + typedef Y range_type; + //! \brief The field type of the preconditioner. + typedef typename X::field_type field_type; + + // define the category + enum { + //! \brief The category the preconditioner is part of. + category=Dune::SolverCategory::sequential + }; + + /*! \brief Constructor. + + Constructor gets all parameters to operate the prec. + \param A The matrix to operate on. + */ + SeqLU (const M& A) : + slu(A, false) + { + } + + /*! + \brief Prepare the preconditioner. + + \copydoc Preconditioner::pre(X&,Y&) + */ + virtual void pre (X& x, Y& b) {} + + /*! + \brief Apply the precondioner. + + \copydoc Preconditioner::apply(X&,const Y&) + */ + virtual void apply (X& v, const Y& d) + { + Dune::InverseOperatorResult res; + Y t(d); + slu.apply(v, t, res); + } + + /*! + \brief Clean up. + + \copydoc Preconditioner::post(X&) + */ + virtual void post (X& x) {} + + private: + Dune::SuperLU slu; + }; +#endif + diff --git a/dune/elasticity/uzawa_solver.hpp b/dune/elasticity/uzawa_solver.hpp new file mode 100644 index 0000000..312787b --- /dev/null +++ b/dune/elasticity/uzawa_solver.hpp @@ -0,0 +1,56 @@ +#pragma once + +namespace Opm { + namespace Elasticity { + + template +class UzawaSolver : public Dune::InverseOperator +{ + public: + UzawaSolver(Dune::InverseOperator* innersolver_, + Dune::InverseOperator* outersolver_, + const Matrix& B_) : + innersolver(innersolver_), outersolver(outersolver_), B(B_) + { + } + + virtual ~UzawaSolver() + { + delete innersolver; + delete outersolver; + } + + void apply(X& x, Y& b, double reduction, + Dune::InverseOperatorResult& res) + { + apply(x, b, res); + } + + void apply(X& x, Y& b, Dune::InverseOperatorResult& res) + { + Vector lambda, lambda2, u, u2; + MortarUtils::extractBlock(u, b, B.N()); + Dune::InverseOperatorResult res2; + u2.resize(u.size()); + u2 = 0; + innersolver->apply(u2, u, res2); + lambda.resize(B.M()); + B.mtv(u2, lambda); + lambda2.resize(lambda.size()); + lambda2 = 0; + outersolver->apply(lambda2, lambda, res2); + MortarUtils::extractBlock(u, b, B.N()); + B.usmv(-1.0, lambda2, u); + u2 = 0; + innersolver->apply(u2, u, res); + MortarUtils::injectBlock(x, u2, B.N()); + MortarUtils::injectBlock(x, lambda2, B.M(), B.N()); + } + protected: + Dune::InverseOperator* innersolver; + Dune::InverseOperator* outersolver; + const Matrix& B; +}; + +} +} diff --git a/examples/upscale_elasticity.cpp b/examples/upscale_elasticity.cpp index fb6bc3d..d579503 100644 --- a/examples/upscale_elasticity.cpp +++ b/examples/upscale_elasticity.cpp @@ -27,6 +27,7 @@ #endif #include +#include //! \brief Display the available command line parameters @@ -34,29 +35,36 @@ void syntax(char** argv) { std::cerr << "Usage: " << argv[0] << " gridfilename=filename.grdecl [method=]" << std::endl << "\t[xmax=] [ymax=] [zmax=] [xmin=] [ymin=] [zmin=] [linsolver_type=]" << std::endl - <<" \t[Emin=] [ctol=] [ltol=] [rock_list=] [vtufilename=] [output=]" << std::endl << std::endl - << "\t gridfilename - the grid file. can be 'uniform'" << std::endl - << "\t vtufilename - save results to vtu file" << std::endl - << "\t rock_list - use material information from given rocklist" << std::endl - << "\t output - output results in text report format" << std::endl - << "\t method - can be MPC, LLM or Mortar" << std::endl + <<" \t[Emin=] [ctol=] [ltol=] [rock_list=] [vtufilename=] [output=] [linsolver_verbose=]" << std::endl << std::endl + << "\t gridfilename - the grid file. can be 'uniform'" << std::endl + << "\t vtufilename - save results to vtu file" << std::endl + << "\t rock_list - use material information from given rocklist" << std::endl + << "\t output - output results in text report format" << std::endl + << "\t method - can be MPC or Mortar" << std::endl << "\t\t if not specified, Mortar couplings are used" << std::endl << "\t (x|y|z)max/(x|y|z)min - maximum/minimum grid coordinates." << std::endl << "\t\t if not specified, coordinates are found by grid traversal" << std::endl - << "\t cells(x|y|z) - number of cells if gridfilename=uniform." << std::endl - << "\t Emin - minimum E modulus (to avoid numerical issues)" << std::endl - << "\t ctol - collapse tolerance in grid parsing" << std::endl - << "\t ltol - tolerance in iterative linear solvers" << std::endl - << "\t linsolver_type=cg - use the conjugate gradient method" << std::endl - << "\t linsolver_type=slu - use the SuperLU sparse direct solver" << std::endl; + << "\t cells(x|y|z) - number of cells if gridfilename=uniform." << std::endl + << "\t lambda(x|y) - order of lagrangian multipliers in first/second direction" << std::endl + << "\t Emin - minimum E modulus (to avoid numerical issues)" << std::endl + << "\t ctol - collapse tolerance in grid parsing" << std::endl + << "\t ltol - tolerance in iterative linear solvers" << std::endl + << "\t linsolver_type=iterative - use a suitable iterative method (cg or gmres)" << std::endl + << "\t linsolver_type=direct - use the SuperLU sparse direct solver" << std::endl + << "\t verbose - set to true to get verbose output" << std::endl + << "\t linsolver_restart - number of iterations before gmres is restarted" << std::endl + << "\t linsolver_presteps - number of pre-smooth steps in the AMG" << std::endl + << "\t linsolver_poststeps - number of post-smooth steps in the AMG" << std::endl + << "\t\t affects memory usage" << std::endl + << "\t linsolver_symmetric - use symmetric linear solver. Defaults to true" << std::endl + << "\t mortar_precond - preconditioner for mortar block. Defaults to schur-amg" << std::endl; } enum UpscaleMethod { UPSCALE_NONE = 0, UPSCALE_MPC = 1, - UPSCALE_LLM = 2, - UPSCALE_MORTAR = 3 + UPSCALE_MORTAR = 2 }; //! \brief Structure holding parameters configurable from command line @@ -75,28 +83,28 @@ struct Params { double Emin; //! \brief The tolerance for collapsing nodes in the Z direction double ctol; - //! \brief The tolerance for the iterative linear solver - double ltol; //! \brief Minimum grid coordinates double min[3]; //! \brief Maximum grid coordinates double max[3]; - //! \brief The linear solver to employ - Opm::Elasticity::Solver solver; + //! \brief Polynomial order of multipliers in first / second direction + int lambda[2]; + //! \brief Number of elements on interface grid in the x direction + int n1; + //! \brief Number of elements on interface grid in the y direction + int n2; + Opm::Elasticity::LinSolParams linsolver; // Debugging options - //! \brief Number of elements on interface grid in the x direction (LLM) - int n1; - //! \brief Number of elements on interface grid in the y direction (LLM) - int n2; + //! \brief cell in x int cellsx; //! \brief cell in y int cellsy; //! \brief cell in z int cellsz; - //! \brief Polynomial order of multipliers in first / second direction - int lambda[2]; + //! \brief verbose output + bool verbose; }; //! \brief Parse the command line arguments @@ -109,38 +117,32 @@ void parseCommandLine(int argc, char** argv, Params& p) p.min[0] = param.getDefault("xmin",-1); p.min[1] = param.getDefault("ymin",-1); p.min[2] = param.getDefault("zmin",-1); - p.lambda[0] = param.getDefault("lambdax", 1); - p.lambda[1] = param.getDefault("lambday", 1); + p.lambda[0] = param.getDefault("lambdax", 2); + p.lambda[1] = param.getDefault("lambday", 2); std::string method = param.getDefault("method","mortar"); if (!strcasecmp(method.c_str(),"mpc")) p.method = UPSCALE_MPC; - if (!strcasecmp(method.c_str(),"llm")) - p.method = UPSCALE_LLM; if (!strcasecmp(method.c_str(),"mortar")) p.method = UPSCALE_MORTAR; if (!strcasecmp(method.c_str(),"none")) p.method = UPSCALE_NONE; p.Emin = param.getDefault("Emin",0.0); - p.ctol = param.getDefault("ctol",1.e-8); - p.ltol = param.getDefault("ltol",1.e-10); + p.ctol = param.getDefault("ctol",1.e-6); p.file = param.get("gridfilename"); p.rocklist = param.getDefault("rock_list",""); p.vtufile = param.getDefault("vtufilename",""); p.output = param.getDefault("output",""); + p.verbose = param.getDefault("verbose",false); size_t i; if ((i=p.vtufile.find(".vtu")) != std::string::npos) p.vtufile = p.vtufile.substr(0,i); - //std::string solver = param.getDefault("linsolver_type","slu"); - std::string solver = param.getDefault("linsolver_type","cg"); - if (solver == "cg") - p.solver = Opm::Elasticity::CG; - else - p.solver = Opm::Elasticity::SLU; + if (p.file == "uniform") { p.cellsx = param.getDefault("cellsx",3); p.cellsy = param.getDefault("cellsy",3); p.cellsz = param.getDefault("cellsz",3); } + p.linsolver.parse(param); p.n1 = -1; p.n2 = -1; } @@ -161,8 +163,6 @@ void writeOutput(const Params& p, Opm::time::StopWatch& watch, int cells, gethostname(hostname,1024); std::string method = "mortar"; - if (p.method == UPSCALE_LLM) - method = "llm"; if (p.method == UPSCALE_MPC) method = "mpc"; if (p.method == UPSCALE_NONE) @@ -194,10 +194,10 @@ void writeOutput(const Params& p, Opm::time::StopWatch& watch, int cells, } f << "# Options used:" << std::endl << "#\t method: " << method << std::endl - << "#\t linsolver_type: " << (p.solver==Opm::Elasticity::SLU?"slu":"cg") + << "#\t linsolver_type: " << (p.linsolver.type==Opm::Elasticity::DIRECT?"direct":"iterative") << std::endl; - if (p.solver == Opm::Elasticity::CG) - f << "#\t ltol: " << p.ltol << std::endl; + if (p.linsolver.type == Opm::Elasticity::ITERATIVE) + f << "#\t ltol: " << p.linsolver.tol << std::endl; if (p.file == "uniform") { f << "#\t cellsx: " << p.cellsx << std::endl << "#\t cellsy: " << p.cellsy << std::endl @@ -245,7 +245,10 @@ int main(int argc, char** argv) grid.readEclipseFormat(p.file,p.ctol,false); typedef GridType::ctype ctype; - Opm::Elasticity::ElasticityUpscale upscale(grid, p.ctol, p.Emin, p.file, p.rocklist); + Opm::Elasticity::ElasticityUpscale upscale(grid, p.ctol, + p.Emin, p.file, + p.rocklist, + p.verbose); if (p.max[0] < 0 || p.min[0] < 0) { std::cout << "determine side coordinates..." << std::endl; upscale.findBoundaries(p.min,p.max); @@ -256,36 +259,54 @@ int main(int argc, char** argv) p.n1 = grid.logicalCartesianSize()[0]; p.n2 = grid.logicalCartesianSize()[1]; } + if (p.linsolver.zcells == -1) { + double lz = p.max[2]-p.min[2]; + int nz = grid.logicalCartesianSize()[2]; + double hz = lz/nz; + double lp = sqrt((double)(p.max[0]-p.min[0])*(p.max[1]-p.min[1])); + int np = std::max(grid.logicalCartesianSize()[0], + grid.logicalCartesianSize()[1]); + double hp = lp/np; + p.linsolver.zcells = (int)(2*hp/hz+0.5); + } + std::cout << "logical dimension: " << grid.logicalCartesianSize()[0] + << "x" << grid.logicalCartesianSize()[1] + << "x" << grid.logicalCartesianSize()[2] + << std::endl; - if (p.method == UPSCALE_LLM) { - std::cout << "using LLM couplings..." << std::endl; - upscale.periodicBCsLLM(p.min,p.max,p.n1,p.n2); - } else if (p.method == UPSCALE_MPC) { + if (p.method == UPSCALE_MPC) { std::cout << "using MPC couplings in all directions..." << std::endl; - upscale.periodicBCs(p.min,p.max); + upscale.periodicBCs(p.min, p.max); std::cout << "preprocessing grid..." << std::endl; upscale.A.initForAssembly(); } else if (p.method == UPSCALE_MORTAR) { std::cout << "using Mortar couplings.." << std::endl; - upscale.periodicBCsMortar(p.min,p.max,p.n1,p.n2,p.lambda[0], p.lambda[1]); + upscale.periodicBCsMortar(p.min, p.max, p.n1, p.n2, + p.lambda[0], p.lambda[1]); } else if (p.method == UPSCALE_NONE) { std::cout << "no periodicity approach applied.." << std::endl; upscale.fixCorners(p.min, p.max); upscale.A.initForAssembly(); } Dune::FieldMatrix C; - Dune::VTKWriter vtkwriter(grid.leafView()); + Dune::VTKWriter* vtkwriter=0; + if (!p.vtufile.empty()) + vtkwriter = new Dune::VTKWriter(grid.leafView()); Opm::Elasticity::Vector field[6]; - std::cout << "assembling..." << "\n"; + std::cout << "assembling elasticity operator..." << "\n"; upscale.assemble(-1,true); - upscale.setupSolvers(p.solver); -#pragma omp parallel for schedule(static) + std::cout << "setting up linear solver..." << std::endl; + upscale.setupSolvers(p.linsolver); + +//#pragma omp parallel for schedule(static) for (int i=0;i<6;++i) { + std::cout << "processing case " << i+1 << "..." << std::endl; + std::cout << "\tassembling load vector..." << std::endl; upscale.assemble(i,false); - std::cout << "solving case " << i+1 << "..." << "\n"; - upscale.solve(p.solver,p.ltol,i); + std::cout << "\tsolving..." << std::endl; + upscale.solve(i); upscale.A.expandSolution(field[i],upscale.u[i]); -#define CLAMP(x) (fabs(x)<1.e-5?0.f:x) +#define CLAMP(x) (fabs(x)<1.e-4?0.0:x) for (size_t j=0;jaddVertexData(field[i], str.str().c_str(), dim); + } + if (vtkwriter) { + vtkwriter->write(p.vtufile); + delete vtkwriter; } - if (!p.vtufile.empty()) - vtkwriter.write(p.vtufile); // voigt notation for (int j=0;j<6;++j) std::swap(C[3][j],C[5][j]); From cd8dfb11b703dd261086324867eab538112c32ef Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Mon, 7 Jan 2013 09:43:35 +0100 Subject: [PATCH 15/17] fixed: don't rely on a patched istl --- dune/elasticity/mortar_schur_precond.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dune/elasticity/mortar_schur_precond.hpp b/dune/elasticity/mortar_schur_precond.hpp index c780422..c3875ce 100644 --- a/dune/elasticity/mortar_schur_precond.hpp +++ b/dune/elasticity/mortar_schur_precond.hpp @@ -48,7 +48,7 @@ class MortarSchurPre : public Dune::Preconditioner { MortarSchurPre(const Matrix& P_, const Matrix& B_, PrecondElasticityBlock& Apre_, bool symmetric_=false) : Apre(Apre_), B(B_), N(B.N()), M(B.M()), - Lpre(P_, false, false), symmetric(symmetric_) + Lpre(P_, false), symmetric(symmetric_) { } From b3b547a4f10c1d506fccc833694fdc4f25250dbe Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Mon, 7 Jan 2013 09:54:34 +0100 Subject: [PATCH 16/17] changed: get rid of case-insensitive comparisons in command line parameter parsing the paramgroup parser is case sensitive, the data may as well be as well --- examples/upscale_elasticity.cpp | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/examples/upscale_elasticity.cpp b/examples/upscale_elasticity.cpp index d579503..b1584e2 100644 --- a/examples/upscale_elasticity.cpp +++ b/examples/upscale_elasticity.cpp @@ -17,6 +17,7 @@ #include #include +#include #include // We use exceptions #include #include @@ -120,11 +121,11 @@ void parseCommandLine(int argc, char** argv, Params& p) p.lambda[0] = param.getDefault("lambdax", 2); p.lambda[1] = param.getDefault("lambday", 2); std::string method = param.getDefault("method","mortar"); - if (!strcasecmp(method.c_str(),"mpc")) + if (method == "mpc") p.method = UPSCALE_MPC; - if (!strcasecmp(method.c_str(),"mortar")) + if (method == "mortar") p.method = UPSCALE_MORTAR; - if (!strcasecmp(method.c_str(),"none")) + if (method == "none") p.method = UPSCALE_NONE; p.Emin = param.getDefault("Emin",0.0); p.ctol = param.getDefault("ctol",1.e-6); @@ -220,9 +221,9 @@ int main(int argc, char** argv) typedef Dune::CpGrid GridType; - if (argc < 2 || strcasecmp(argv[1],"-h") == 0 - || strcasecmp(argv[1],"--help") == 0 - || strcasecmp(argv[1],"-?") == 0) { + if (argc < 2 || strcmp(argv[1],"-h") == 0 + || strcmp(argv[1],"--help") == 0 + || strcmp(argv[1],"-?") == 0) { syntax(argv); exit(1); } From 42b271cfd17b4a94747a3b574ef286ebcfbd80de Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Mon, 7 Jan 2013 10:01:20 +0100 Subject: [PATCH 17/17] changed: certain elements do not like pragma once --- dune/elasticity/applier.hpp | 5 ++++- dune/elasticity/asmhandler.hpp | 5 ++++- dune/elasticity/boundarygrid.hh | 5 ++++- dune/elasticity/elasticity.hpp | 5 ++++- dune/elasticity/elasticity_upscale.hpp | 5 ++++- dune/elasticity/logutils.hpp | 5 ++++- dune/elasticity/material.hh | 5 ++++- dune/elasticity/materials.hh | 5 ++++- dune/elasticity/matrixops.hpp | 5 ++++- dune/elasticity/mortar_evaluator.hpp | 5 ++++- dune/elasticity/mortar_schur.hpp | 5 ++++- dune/elasticity/mortar_schur_precond.hpp | 5 ++++- dune/elasticity/mortar_utils.hpp | 16 +++++++++++++++- dune/elasticity/mpc.hh | 5 ++++- dune/elasticity/seqlu.hpp | 4 +++- dune/elasticity/shapefunctions.hpp | 5 ++++- dune/elasticity/uzawa_solver.hpp | 16 +++++++++++++++- 17 files changed, 89 insertions(+), 17 deletions(-) diff --git a/dune/elasticity/applier.hpp b/dune/elasticity/applier.hpp index 7db2c4c..4442bfe 100644 --- a/dune/elasticity/applier.hpp +++ b/dune/elasticity/applier.hpp @@ -10,7 +10,8 @@ //! and preconditioners in DUNE //! //============================================================================== -#pragma once +#ifndef APPLIER_H_ +#define APPLIER_H_ namespace Opm { namespace Elasticity { @@ -45,3 +46,5 @@ void PreApplier::apply(Vector& v, Vector& d) } } + +#endif diff --git a/dune/elasticity/asmhandler.hpp b/dune/elasticity/asmhandler.hpp index b37842e..f34b095 100644 --- a/dune/elasticity/asmhandler.hpp +++ b/dune/elasticity/asmhandler.hpp @@ -9,7 +9,8 @@ //! \brief Class handling finite element assembly //! //============================================================================== -#pragma once +#ifndef ASMHANDLER_HPP_ +#define ASMHANDLER_HPP_ #include #include @@ -233,3 +234,5 @@ class ASMHandler { } } + +#endif diff --git a/dune/elasticity/boundarygrid.hh b/dune/elasticity/boundarygrid.hh index 2fb054d..75b8fc2 100644 --- a/dune/elasticity/boundarygrid.hh +++ b/dune/elasticity/boundarygrid.hh @@ -9,7 +9,8 @@ //! \brief Class describing 2D quadrilateral grids //! //============================================================================== -#pragma once +#ifndef BOUNDARYGRID_HH_ +#define BOUNDARYGRID_HH_ #include #include @@ -500,3 +501,5 @@ BoundaryGrid::Vertex minXmaxY(std::vector& in); } } + +#endif diff --git a/dune/elasticity/elasticity.hpp b/dune/elasticity/elasticity.hpp index ecc75f7..cd2bd99 100644 --- a/dune/elasticity/elasticity.hpp +++ b/dune/elasticity/elasticity.hpp @@ -9,7 +9,8 @@ //! \brief Elasticity helper class //! //============================================================================== -#pragma once +#ifndef ELASTICITY_HPP_ +#define ELASTICITY_HPP_ namespace Opm { namespace Elasticity { @@ -69,3 +70,5 @@ class Elasticity { } } + +#endif diff --git a/dune/elasticity/elasticity_upscale.hpp b/dune/elasticity/elasticity_upscale.hpp index 7320e5d..98605a1 100644 --- a/dune/elasticity/elasticity_upscale.hpp +++ b/dune/elasticity/elasticity_upscale.hpp @@ -9,7 +9,8 @@ //! \brief Elasticity upscale class //! //============================================================================== -#pragma once +#ifndef ELASTICITY_UPSCALE_HPP_ +#define ELASTICITY_UPSCALE_HPP_ #include #include @@ -444,3 +445,5 @@ class ElasticityUpscale } } + +#endif diff --git a/dune/elasticity/logutils.hpp b/dune/elasticity/logutils.hpp index 2dd2415..2bd92ea 100644 --- a/dune/elasticity/logutils.hpp +++ b/dune/elasticity/logutils.hpp @@ -9,7 +9,8 @@ //! \brief Logging helper utilities //! //============================================================================== -#pragma once +#ifndef LOGUTILS_HPP_ +#define LOGUTILS_HPP_ class LoggerHelper { public: @@ -29,3 +30,5 @@ class LoggerHelper { int per_chunk; //!< Will log for each per_chunk processed int max; //!< Total number of its }; + +#endif diff --git a/dune/elasticity/material.hh b/dune/elasticity/material.hh index 8fe0634..92ee786 100644 --- a/dune/elasticity/material.hh +++ b/dune/elasticity/material.hh @@ -9,7 +9,8 @@ //! \brief Material interface. //! //============================================================================== -#pragma once +#ifndef MATERIAL_HH_ +#define MATERIAL_HH_ #include #include @@ -97,3 +98,5 @@ private: } } + +#endif diff --git a/dune/elasticity/materials.hh b/dune/elasticity/materials.hh index 9fd3d83..2eaca60 100644 --- a/dune/elasticity/materials.hh +++ b/dune/elasticity/materials.hh @@ -9,7 +9,8 @@ //! \brief Material properties. //! //============================================================================== -#pragma once +#ifndef MATERIALS_HH_ +#define MATERIALS_HH_ #include "material.hh" @@ -173,3 +174,5 @@ private: } } + +#endif diff --git a/dune/elasticity/matrixops.hpp b/dune/elasticity/matrixops.hpp index 73f7a82..373d289 100644 --- a/dune/elasticity/matrixops.hpp +++ b/dune/elasticity/matrixops.hpp @@ -9,7 +9,8 @@ //! \brief Helper class with some matrix operations //! //============================================================================== -#pragma once +#ifndef MATRIXOPS_HPP_ +#define MATRIXOPS_HPP_ #include @@ -87,3 +88,5 @@ class MatrixOps { } } + +#endif diff --git a/dune/elasticity/mortar_evaluator.hpp b/dune/elasticity/mortar_evaluator.hpp index c5713de..68c7dd8 100644 --- a/dune/elasticity/mortar_evaluator.hpp +++ b/dune/elasticity/mortar_evaluator.hpp @@ -9,7 +9,8 @@ //! \brief Linear operator for a Mortar block //! //============================================================================== -#pragma once +#ifndef MORTAR_EVALUATOR_HPP_ +#define MORTAR_EVALUATOR_HPP_ #include #include @@ -77,3 +78,5 @@ class MortarEvaluator : public Dune::LinearOperator { } } + +#endif diff --git a/dune/elasticity/mortar_schur.hpp b/dune/elasticity/mortar_schur.hpp index 25b87f7..7093873 100644 --- a/dune/elasticity/mortar_schur.hpp +++ b/dune/elasticity/mortar_schur.hpp @@ -9,7 +9,8 @@ //! \brief Schur based linear operator for a Mortar block //! //============================================================================== -#pragma once +#ifndef MORTAR_SCHUR_HPP_ +#define MORTAR_SCHUR_HPP_ #include #include @@ -71,3 +72,5 @@ class MortarBlockEvaluator : public Dune::LinearOperator { } } + +#endif diff --git a/dune/elasticity/mortar_schur_precond.hpp b/dune/elasticity/mortar_schur_precond.hpp index c3875ce..9dfe7a7 100644 --- a/dune/elasticity/mortar_schur_precond.hpp +++ b/dune/elasticity/mortar_schur_precond.hpp @@ -9,7 +9,8 @@ //! \brief Schur based preconditioner for saddle point systems //! //============================================================================== -#pragma once +#ifndef SCHUR_PRECOND_HPP_ +#define SCHUR_PRECOND_HPP_ #include #include @@ -121,3 +122,5 @@ class MortarSchurPre : public Dune::Preconditioner { } } + +#endif diff --git a/dune/elasticity/mortar_utils.hpp b/dune/elasticity/mortar_utils.hpp index 1e1c552..23c3c78 100644 --- a/dune/elasticity/mortar_utils.hpp +++ b/dune/elasticity/mortar_utils.hpp @@ -1,4 +1,16 @@ -#pragma once +//============================================================================== +//! +//! \file mortar_utils.hpp +//! +//! \date Nov 9 2012 +//! +//! \author Arne Morten Kvarving / SINTEF +//! +//! \brief Mortar helper class +//! +//============================================================================== +#ifndef MORTAR_UTILS_HPP_ +#define MORTAR_UTILS_HPP_ namespace Opm { namespace Elasticity { @@ -29,3 +41,5 @@ class MortarUtils { } } + +#endif diff --git a/dune/elasticity/mpc.hh b/dune/elasticity/mpc.hh index 0108ef6..3986ecf 100644 --- a/dune/elasticity/mpc.hh +++ b/dune/elasticity/mpc.hh @@ -9,7 +9,8 @@ //! \brief Representation of multi-point constraint (MPC) equations. //! //============================================================================== -#pragma once +#ifndef MPC_HH_ +#define MPC_HH_ #include #include @@ -159,3 +160,5 @@ typedef std::map MPCMap; } } + +#endif diff --git a/dune/elasticity/seqlu.hpp b/dune/elasticity/seqlu.hpp index e85c5bf..9c6020a 100644 --- a/dune/elasticity/seqlu.hpp +++ b/dune/elasticity/seqlu.hpp @@ -9,7 +9,8 @@ //! \brief SuperLU preconditioner interface //! //============================================================================== -#pragma once +#ifndef SEQLU_HPP_ +#define SEQLU_HPP_ #if HAVE_SUPERLU #include @@ -75,3 +76,4 @@ class SeqLU : public Dune::Preconditioner { }; #endif +#endif diff --git a/dune/elasticity/shapefunctions.hpp b/dune/elasticity/shapefunctions.hpp index 3047c9d..eef6723 100644 --- a/dune/elasticity/shapefunctions.hpp +++ b/dune/elasticity/shapefunctions.hpp @@ -9,7 +9,8 @@ //! \brief Classes for shape functions. Loosely based on code in dune-grid-howto //! //============================================================================== -#pragma once +#ifndef SHAPEFUNCTIONS_HPP_ +#define SHAPEFUNCTIONS_HPP_ #include #include "dynmatrixev.hh" @@ -410,3 +411,5 @@ protected: } } + +#endif diff --git a/dune/elasticity/uzawa_solver.hpp b/dune/elasticity/uzawa_solver.hpp index 312787b..890c8de 100644 --- a/dune/elasticity/uzawa_solver.hpp +++ b/dune/elasticity/uzawa_solver.hpp @@ -1,4 +1,16 @@ -#pragma once +#ifndef UZAWA_SOLVER_HPP_ +#define UZAWA_SOLVER_HPP_ +//============================================================================== +//! +//! \file uzawa_solver.hpp +//! +//! \date Nov 9 2012 +//! +//! \author Arne Morten Kvarving / SINTEF +//! +//! \brief Uzawa scheme helper class +//! +//============================================================================== namespace Opm { namespace Elasticity { @@ -54,3 +66,5 @@ class UzawaSolver : public Dune::InverseOperator } } + +#endif