From 499f5092abb89afcd93d0429c7d12c9367b9f3e1 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Tue, 20 Nov 2012 10:17:48 +0100 Subject: [PATCH] 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];