contactmerge.hh

Go to the documentation of this file.

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-GPL-2.0-only-with-dune-grid-glue-exception
#ifndef DUNE_GRIDGLUE_MERGING_CONTACTMERGE_HH
#define DUNE_GRIDGLUE_MERGING_CONTACTMERGE_HH
 
 
#include <iostream>
#include <fstream>
#include <iomanip>
#include <vector>
#include <algorithm>
#include <limits>
#include <memory>
#include <functional>
 
#include <dune/common/fvector.hh>
#include <dune/common/exceptions.hh>
#include <dune/common/bitsetvector.hh>
#include <dune/common/deprecated.hh>
 
#include <dune/grid/common/grid.hh>
 
#include <dune/grid-glue/merging/standardmerge.hh>
#include <dune/grid-glue/gridglue.hh>
 
namespace Dune {
namespace GridGlue {
 
template<int dimworld, typename T = double>
class ContactMerge
: public StandardMerge<T,dimworld-1,dimworld-1,dimworld>
{
    static constexpr int dim = dimworld-1;
 
    static_assert( dim==1 || dim==2,
            "ContactMerge yet only handles the cases dim==1 and dim==2!");
 
    typedef StandardMerge<T,dim,dim,dimworld> Base;
public:
 
    /*   E X P O R T E D   T Y P E S   A N D   C O N S T A N T S   */
 
    typedef T ctype;
 
    typedef Dune::FieldVector<T, dimworld>  WorldCoords;
 
    typedef Dune::FieldVector<T, dim>  LocalCoords;
 
    enum ProjectionType {OUTER_NORMAL, CLOSEST_POINT};
    ContactMerge(const T allowedOverlap=T(0),
                 std::function<WorldCoords(WorldCoords)> domainDirections=nullptr,
                 std::function<WorldCoords(WorldCoords)> targetDirections=nullptr,
                 ProjectionType type = OUTER_NORMAL)
        : domainDirections_(domainDirections), targetDirections_(targetDirections),
          overlap_(allowedOverlap), type_(type)
    {}
 
    ContactMerge(const T allowedOverlap, ProjectionType type)
        : overlap_(allowedOverlap),
          type_(type)
    {}
 
    inline
    void setSurfaceDirections(std::function<WorldCoords(WorldCoords)> domainDirections,
                              std::function<WorldCoords(WorldCoords)> targetDirections)
    {
        domainDirections_ = domainDirections;
        targetDirections_ = targetDirections;
        this->valid = false;
    }
 
    void setOverlap(T overlap)
    {
        overlap_ = overlap;
    }
 
    T getOverlap() const
    {
        return overlap_;
    }
 
    void minNormalAngle(T angle)
    {
        using std::cos;
        maxNormalProduct_ = cos(angle);
    }
 
    T minNormalAngle() const
    {
        using std::acos;
        return acos(maxNormalProduct_);
    }
 
protected:
  typedef typename StandardMerge<T,dimworld-1,dimworld-1,dimworld>::SimplicialIntersection SimplicialIntersection;
 
private:
    std::function<WorldCoords(WorldCoords)> domainDirections_;
    std::vector<WorldCoords> nodalDomainDirections_;
 
    std::function<WorldCoords(WorldCoords)> targetDirections_;
    std::vector<WorldCoords> nodalTargetDirections_;
 
    T overlap_;
 
    ProjectionType type_;
 
    T maxNormalProduct_ = T(-0.1);
 
    void computeIntersections(const Dune::GeometryType& grid1ElementType,
            const std::vector<Dune::FieldVector<T,dimworld> >& grid1ElementCorners,
            std::bitset<(1<<dim)>& neighborIntersects1,
            unsigned int grid1Index,
            const Dune::GeometryType& grid2ElementType,
            const std::vector<Dune::FieldVector<T,dimworld> >& grid2ElementCorners,
            std::bitset<(1<<dim)>& neighborIntersects2,
            unsigned int grid2Index,
            std::vector<SimplicialIntersection>& intersections) override;
 
protected:
    void build(const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
        const std::vector<unsigned int>& grid1Elements,
        const std::vector<Dune::GeometryType>& grid1ElementTypes,
        const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
        const std::vector<unsigned int>& grid2Elements,
        const std::vector<Dune::GeometryType>& grid2ElementTypes) override
    {
        std::cout<<"ContactMerge building grid!\n";
        // setup the nodal direction vectors
        setupNodalDirections(grid1Coords, grid1Elements, grid1ElementTypes,
                grid2Coords, grid2Elements, grid2ElementTypes);
 
        Base::build(grid1Coords, grid1Elements, grid1ElementTypes,
                   grid2Coords, grid2Elements, grid2ElementTypes);
 
    }
 
private:
 
    static LocalCoords localCornerCoords(int i, const Dune::GeometryType& gt)
    {
        const auto& ref = Dune::ReferenceElements<T,dim>::general(gt);
        return ref.position(i,dim);
    }
 
protected:
 
    void computeCyclicOrder(const std::vector<std::array<LocalCoords,2> >& polytopeCorners,
            const LocalCoords& center, std::vector<int>& ordering) const;
 
    void setupNodalDirections(const std::vector<WorldCoords>& coords1,
            const std::vector<unsigned int>& elements1,
            const std::vector<Dune::GeometryType>& elementTypes1,
            const std::vector<WorldCoords>& coords2,
            const std::vector<unsigned int>& elements2,
            const std::vector<Dune::GeometryType>& elementTypes2);
 
    void computeOuterNormalField(const std::vector<WorldCoords>& coords,
            const std::vector<unsigned int>& elements,
            const std::vector<Dune::GeometryType>& elementTypes,
            std::vector<WorldCoords>& normals);
 
    void removeDoubles(std::vector<std::array<LocalCoords,2> >& polytopeCorners);
};
 
} /* namespace GridGlue */
} /* namespace Dune */
 
#include "contactmerge.cc"
 
#endif // DUNE_GRIDGLUE_MERGING_CONTACTMERGE_HH