Geometry.hpp

Go to the documentation of this file.

//===========================================================================
//
// File: Geometry.hpp
//
// Created: Fri May 29 23:29:24 2009
//
// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
//            B�rd Skaflestad     <bard.skaflestad@sintef.no>
//            Antonella Ritorto   <antonella.ritorto@opm-op.com>
//
// $Date$
//
// $Revision$
//
//===========================================================================
 
/*
  Copyright 2009, 2010, 2011 SINTEF ICT, Applied Mathematics.
  Copyright 2009, 2010, 2011, 2022 Equinor ASA.
 
  This file is part of the Open Porous Media project (OPM).
 
  OPM is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation, either version 3 of the License, or
  (at your option) any later version.
 
  OPM is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
 
  You should have received a copy of the GNU General Public License
  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
*/
 
#ifndef OPM_GEOMETRY_HEADER
#define OPM_GEOMETRY_HEADER
 
#include <cmath>
 
// Warning suppression for Dune includes.
#include <opm/grid/utility/platform_dependent/disable_warnings.h>
 
#include <dune/common/version.hh>
#include <dune/geometry/referenceelements.hh>
#include <dune/grid/common/geometry.hh>
 
#include <dune/geometry/type.hh>
 
#include <opm/grid/cpgrid/EntityRep.hpp>
#include <opm/grid/cpgrid/DefaultGeometryPolicy.hpp>
#include <opm/grid/cpgrid/OrientedEntityTable.hpp>
#include <opm/grid/cpgpreprocess/preprocess.h>
#include <opm/grid/common/Volumes.hpp>
#include <opm/grid/utility/platform_dependent/reenable_warnings.h>
#include <opm/common/utility/SparseTable.hpp>
 
#include <opm/grid/utility/ErrorMacros.hpp>
 
namespace Dune
{
    namespace cpgrid
    {
 
        template <int mydim, int cdim>
        class Geometry
        {
        };
 
 
 
 
        template <int cdim> // GridImp arg never used
        class Geometry<0, cdim>
        {
            static_assert(cdim == 3, "");
        public:
            enum { dimension = 3 };
            enum { mydimension = 0};
            enum { coorddimension = cdim };
            enum { dimensionworld = 3 };
 
            using ctype = double;
            using Volume = ctype;
 
            typedef FieldVector<ctype, mydimension> LocalCoordinate;
            typedef FieldVector<ctype, coorddimension> GlobalCoordinate;
 
            typedef FieldMatrix< ctype, coorddimension, mydimension >         Jacobian;
            typedef FieldMatrix< ctype, coorddimension, mydimension >         JacobianInverse;
            typedef FieldMatrix< ctype, mydimension, coorddimension >         JacobianTransposed;
            typedef FieldMatrix< ctype, coorddimension, mydimension >         JacobianInverseTransposed;
 
 
            Geometry(const GlobalCoordinate& pos)
                : pos_(pos)
            {
            }
 
            Geometry()
                : pos_(0.0)
            {
            }
 
            const GlobalCoordinate& global(const LocalCoordinate&) const
            {
                return pos_;
            }
 
            LocalCoordinate local(const GlobalCoordinate&) const
            {
                // return 0 to make the geometry check happy.
                return LocalCoordinate(0.0);
            }
 
            Volume integrationElement(const LocalCoordinate&) const
            {
                return volume();
            }
 
            GeometryType type() const
            {
                return Dune::GeometryTypes::cube(mydimension);
            }
 
            int corners() const
            {
                return 1;
            }
 
            GlobalCoordinate corner(int cor) const
            {
                static_cast<void>(cor);
                assert(cor == 0);
                return pos_;
            }
 
            Volume volume() const
            {
                return 1.0;
            }
 
            const GlobalCoordinate& center() const
            {
                return pos_;
            }
 
            JacobianTransposed
            jacobianTransposed(const LocalCoordinate& /* local */) const
            {
 
                // Meaningless to call jacobianTransposed() on singular geometries. But we need to make DUNE happy.
                return {};
            }
 
            JacobianInverseTransposed
            jacobianInverseTransposed(const LocalCoordinate& /*local*/) const
            {
                // Meaningless to call jacobianInverseTransposed() on singular geometries. But we need to make DUNE happy.
                return {};
            }
 
            Jacobian
            jacobian(const LocalCoordinate& /*local*/) const
            {
                return {};
            }
 
            JacobianInverse
            jacobianInverse(const LocalCoordinate& /*local*/) const
            {
                return {};
            }
 
            bool affine() const
            {
                return true;
            }
 
        private:
            GlobalCoordinate pos_;
        };  // class Geometry<0,cdim>
 
 
 
 
        template <int cdim> // GridImp arg never used
        class Geometry<2, cdim>
        {
            static_assert(cdim == 3, "");
        public:
            enum { dimension = 3 };
            enum { mydimension = 2 };
            enum { coorddimension = cdim };
            enum { dimensionworld = 3 };
 
            using ctype = double;
            using Volume = ctype;
 
            typedef FieldVector<ctype, mydimension> LocalCoordinate;
            typedef FieldVector<ctype, coorddimension> GlobalCoordinate;
 
            typedef FieldMatrix< ctype, coorddimension, mydimension >         Jacobian;
            typedef FieldMatrix< ctype, coorddimension, mydimension >         JacobianInverse;
            typedef FieldMatrix< ctype, mydimension, coorddimension >         JacobianTransposed;
            typedef FieldMatrix< ctype, coorddimension, mydimension >         JacobianInverseTransposed;
 
            Geometry(const GlobalCoordinate& pos,
                     ctype vol)
                : pos_(pos), vol_(vol)
            {
            }
 
            Geometry()
                : pos_(0.0), vol_(0.0)
            {
            }
 
            const GlobalCoordinate& global(const LocalCoordinate&) const
            {
                OPM_THROW(std::runtime_error, "Geometry::global() meaningless on singular geometry.");
            }
 
            LocalCoordinate local(const GlobalCoordinate&) const
            {
                OPM_THROW(std::runtime_error, "Geometry::local() meaningless on singular geometry.");
            }
 
            Volume integrationElement(const LocalCoordinate&) const
            {
                return vol_;
            }
 
            GeometryType type() const
            {
                return Dune::GeometryTypes::none(mydimension);
            }
 
            int corners() const
            {
                return 0;
            }
 
            GlobalCoordinate corner(int /* cor */) const
            {
                // Meaningless call to cpgrid::Geometry::corner(int):
                //"singular geometry has no corners.
                // But the DUNE tests assume at least one corner.
                return GlobalCoordinate( 0.0 );
            }
 
            Volume volume() const
            {
                return vol_;
            }
 
            const GlobalCoordinate& center() const
            {
                return pos_;
            }
 
            const FieldMatrix<ctype, mydimension, coorddimension>&
            jacobianTransposed(const LocalCoordinate& /* local */) const
            {
                OPM_THROW(std::runtime_error, "Meaningless to call jacobianTransposed() on singular geometries.");
            }
 
            const FieldMatrix<ctype, coorddimension, mydimension>&
            jacobianInverseTransposed(const LocalCoordinate& /*local*/) const
            {
                OPM_THROW(std::runtime_error, "Meaningless to call jacobianInverseTransposed() on singular geometries.");
            }
 
            Jacobian
            jacobian(const LocalCoordinate& /*local*/) const
            {
                return jacobianTransposed({}).transposed();
            }
 
            JacobianInverse
            jacobianInverse(const LocalCoordinate& /*local*/) const
            {
                return jacobianInverseTransposed({}).transposed();
            }
 
            bool affine() const
            {
                return true;
            }
 
        private:
            GlobalCoordinate pos_;
            ctype vol_;
        };
 
 
 
 
 
        template <int cdim>
        class Geometry<3, cdim>
        {
            static_assert(cdim == 3, "");
        public:
            enum { dimension = 3 };
            enum { mydimension = 3 };
            enum { coorddimension = cdim };
            enum { dimensionworld = 3 };
 
            using ctype = double;
            using Volume = ctype;
 
            typedef FieldVector<ctype, mydimension> LocalCoordinate;
            typedef FieldVector<ctype, coorddimension> GlobalCoordinate;
 
            typedef FieldMatrix< ctype, coorddimension, mydimension >         Jacobian;
            typedef FieldMatrix< ctype, coorddimension, mydimension >         JacobianInverse;
            typedef FieldMatrix< ctype, mydimension, coorddimension >         JacobianTransposed;
            typedef FieldMatrix< ctype, coorddimension, mydimension >         JacobianInverseTransposed;
 
            typedef Dune::Impl::FieldMatrixHelper< double >  MatrixHelperType;
 
            Geometry(const GlobalCoordinate& pos,
                     ctype vol,
                     std::shared_ptr<const EntityVariable<cpgrid::Geometry<0, 3>, 3>> allcorners_ptr,
                     const int* corner_indices)
                : pos_(pos), vol_(vol),
                  allcorners_(allcorners_ptr), cor_idx_(corner_indices)
            {
                assert(allcorners_ && corner_indices);
            }
 
            Geometry()
                : pos_(0.0), vol_(0.0), allcorners_(0), cor_idx_(0)
            {
            }
 
            GlobalCoordinate global(const LocalCoordinate& local_coord) const
            {
                static_assert(mydimension == 3, "");
                static_assert(coorddimension == 3, "");
                // uvw = { (1-u, 1-v, 1-w), (u, v, w) }
                LocalCoordinate uvw[2] = { LocalCoordinate(1.0), local_coord };
                uvw[0] -= local_coord;
                // Access pattern for uvw matching ordering of corners.
                const int pat[8][3] = { { 0, 0, 0 },
                                        { 1, 0, 0 },
                                        { 0, 1, 0 },
                                        { 1, 1, 0 },
                                        { 0, 0, 1 },
                                        { 1, 0, 1 },
                                        { 0, 1, 1 },
                                        { 1, 1, 1 } };
                GlobalCoordinate xyz(0.0);
                for (int i = 0; i < 8; ++i) {
                    GlobalCoordinate corner_contrib = corner(i);
                    double factor = 1.0;
                    for (int j = 0; j < 3; ++j) {
                        factor *= uvw[pat[i][j]][j];
                    }
                    corner_contrib *= factor;
                    xyz += corner_contrib;
                }
                return xyz;
            }
 
            LocalCoordinate local(const GlobalCoordinate& y) const
            {
                static_assert(mydimension == 3, "");
                static_assert(coorddimension == 3, "");
                // This code is modified from dune/grid/genericgeometry/mapping.hh
                // \todo: Implement direct computation.
                const ctype epsilon = 1e-12;
                auto refElement = Dune::ReferenceElements<ctype, 3>::cube();
                LocalCoordinate x = refElement.position(0,0);
                LocalCoordinate dx;
                do {
                    // DF^n dx^n = F^n, x^{n+1} -= dx^n
                    JacobianTransposed JT = jacobianTransposed(x);
                    GlobalCoordinate z = global(x);
                    z -= y;
                    MatrixHelperType::template xTRightInvA<3, 3>(JT, z, dx );
                    x -= dx;
                } while (dx.two_norm2() > epsilon*epsilon);
                return x;
            }
 
            Volume integrationElement(const LocalCoordinate& local_coord) const
            {
                JacobianTransposed Jt = jacobianTransposed(local_coord);
                return MatrixHelperType::template sqrtDetAAT<3, 3>(Jt);
            }
 
            GeometryType type() const
            {
                return Dune::GeometryTypes::cube(mydimension);
            }
 
            int corners() const
            {
                return 8;
            }
 
            GlobalCoordinate corner(int cor) const
            {
                assert(allcorners_ && cor_idx_);
                return (allcorners_->data())[cor_idx_[cor]].center();
            }
 
            Volume volume() const
            {
                return vol_;
            }
 
            void set_volume(ctype volume) {
                vol_ = volume;
            }
 
            const GlobalCoordinate& center() const
            {
                return pos_;
            }
 
            const JacobianTransposed
            jacobianTransposed(const LocalCoordinate& local_coord) const
            {
                static_assert(mydimension == 3, "");
                static_assert(coorddimension == 3, "");
 
                // uvw = { (1-u, 1-v, 1-w), (u, v, w) }
                LocalCoordinate uvw[2] = { LocalCoordinate(1.0), local_coord };
                uvw[0] -= local_coord;
                // Access pattern for uvw matching ordering of corners.
                const int pat[8][3] = { { 0, 0, 0 },
                                        { 1, 0, 0 },
                                        { 0, 1, 0 },
                                        { 1, 1, 0 },
                                        { 0, 0, 1 },
                                        { 1, 0, 1 },
                                        { 0, 1, 1 },
                                        { 1, 1, 1 } };
                JacobianTransposed  Jt(0.0);
                for (int i = 0; i < 8; ++i) {
                    for (int deriv = 0; deriv < 3; ++deriv) {
                        // This part contributing to dg/du_{deriv}
                        double factor = 1.0;
                        for (int j = 0; j < 3; ++j) {
                            factor *= (j != deriv) ? uvw[pat[i][j]][j]
                                : (pat[i][j] == 0 ? -1.0 : 1.0);
                        }
                        GlobalCoordinate corner_contrib = corner(i);
                        corner_contrib *= factor;
                        Jt[deriv] += corner_contrib; // using FieldMatrix row access.
                    }
                }
                return Jt;
            }
 
            const JacobianInverseTransposed
            jacobianInverseTransposed(const LocalCoordinate& local_coord) const
            {
                JacobianInverseTransposed Jti = jacobianTransposed(local_coord);
                Jti.invert();
                return Jti;
            }
 
            Jacobian
            jacobian(const LocalCoordinate& local_coord) const
            {
                return jacobianTransposed(local_coord).transposed();
            }
 
            JacobianInverse
            jacobianInverse(const LocalCoordinate& local_coord) const
            {
                return jacobianInverseTransposed(local_coord).transposed();
            }
 
            bool affine() const
            {
                return false;
            }
 
            typedef Dune::FieldVector<double,3> PointType;
            void refineCellifiedPatch(const std::array<int,3>& cells_per_dim,
                                      DefaultGeometryPolicy& all_geom,
                                      std::vector<std::array<int,8>>&  refined_cell_to_point,
                                      cpgrid::OrientedEntityTable<0,1>& refined_cell_to_face,
                                      Opm::SparseTable<int>& refined_face_to_point,
                                      cpgrid::OrientedEntityTable<1,0>& refined_face_to_cell,
                                      cpgrid::EntityVariable<enum face_tag, 1>& refined_face_tags,
                                      cpgrid::SignedEntityVariable<PointType, 1>& refined_face_normals,
                                      const std::array<int,3>& patch_dim,
                                      const std::vector<double>& widthsX,
                                      const std::vector<double>& lengthsY,
                                      const std::vector<double>& heightsZ) const;
 
        private:
            GlobalCoordinate pos_;
            double vol_;
            std::shared_ptr<const EntityVariable<Geometry<0, 3>,3>> allcorners_; // For dimension 3 only
            const int* cor_idx_; // For dimension 3 only
 
            const std::tuple< enum face_tag, int,
                              std::array<int, 4>, std::vector<cpgrid::EntityRep<0>>,
                              LocalCoordinate>
            getIndicesFace(int l, int m, int n, int constant_direction, const std::array<int, 3>& cells_per_dim) const;
        };
 
        // Reference element for a given geometry
        template< int mydim, int cdim >
        auto referenceElement(const cpgrid::Geometry<mydim,cdim>& geo) -> decltype(referenceElement<double,mydim>(geo.type()))
        {
            return referenceElement<double,mydim>(geo.type());
        }
 
    } // namespace cpgrid
} // namespace Dune
 
#include <opm/grid/cpgrid/refinement/trilinear/GeometryRefinement.hpp>
 
#endif // OPM_GEOMETRY_HEADER