computeintersection.hh

// 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_COMPUTEINTERSECTION_HH
#define DUNE_GRIDGLUE_MERGING_COMPUTEINTERSECTION_HH
 
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
 
namespace Dune {
namespace GridGlue {
 
template<int dimWorld, int dim1, int dim2, typename T = double>
class ComputationMethod {
public:
    typedef FieldVector<T, dimWorld> Vector;
    static const int grid1Dimension = dim1;
    static const int grid2Dimension = dim2;
    static const int intersectionDimension = (dim1 < dim2)?(dim1):(dim2);
 
    static bool computeIntersectionPoints(const std::vector<Vector> X,
                                                  const std::vector<Vector> Y,
                                                  std::vector<std::vector<int> >& SX,
                                                  std::vector<std::vector<int> >& SY,
                                                  std::vector<Vector>& P);
    static void grid1_subdivisions(const std::vector<Vector> elementCorners,
                                     std::vector<std::vector<unsigned int> >& subElements,
                                     std::vector<std::vector<int> >& faceIds);
    static void grid2_subdivisions(const std::vector<Vector> elementCorners,
                                     std::vector<std::vector<unsigned int> >& subElements,
                                     std::vector<std::vector<int> >& faceIds);
};
 
template<class CM>
class IntersectionComputation {
private:
    typedef typename CM::Vector V;
    const int dimWorld = V::dimension;
    const int dim1 = CM::grid1Dimension;
    const int dim2 = CM::grid2Dimension;
public:
    static bool computeIntersection(const std::vector<V>& X,
                                    const std::vector<V>& Y,
                                    std::vector<std::vector<int> >& SX,
                                    std::vector<std::vector<int> >& SY,
                                    std::vector<V>& P);
 
    template<int isDim, int dW>
    static void orderPoints(const V& centroid,
                            const std::vector<std::vector<int> >& SX,
                            const std::vector<std::vector<int> >& SY,
                            const std::vector<V>& P,
                            std::vector<std::vector<int> >& H)
    {
        if (isDim > 1)
            orderPoints_(std::integral_constant<int,isDim>(),std::integral_constant<int,dW>(),
                    centroid, SX, SY, P,H);
    }
 
private:
    static void orderPoints_(std::integral_constant<int,1>,
                             std::integral_constant<int,1>,
                             const V& centroid,
                             const std::vector<std::vector<int> >& SX,
                             const std::vector<std::vector<int> >& SY,
                             const std::vector<V>& P,
                             std::vector<std::vector<int> >& H) {}
    static void orderPoints_(std::integral_constant<int,1>,
                             std::integral_constant<int,2>,
                             const V& centroid,
                             const std::vector<std::vector<int> >& SX,
                             const std::vector<std::vector<int> >& SY,
                             const std::vector<V>& P,
                             std::vector<std::vector<int> >& H) {}
    static void orderPoints_(std::integral_constant<int,1>,
                             std::integral_constant<int,3>,
                             const V& centroid,
                             const std::vector<std::vector<int> >& SX,
                             const std::vector<std::vector<int> >& SY,
                             const std::vector<V>& P,
                             std::vector<std::vector<int> >& H) {}
    static void orderPoints_(std::integral_constant<int,2>,
                             std::integral_constant<int,2>,
                             const V& centroid,
                             const std::vector<std::vector<int> >& SX,
                             const std::vector<std::vector<int> >& SY,
                             const std::vector<V>& P,
                             std::vector<std::vector<int> >& H);
    static void orderPoints_(std::integral_constant<int,2>,
                             std::integral_constant<int,3>,
                             const V& centroid,
                             const std::vector<std::vector<int> >& SX,
                             const std::vector<std::vector<int> >& SY,
                             const std::vector<V>& P,
                             std::vector<std::vector<int> >& H);
    static void orderPoints_(std::integral_constant<int,3>,
                             std::integral_constant<int,3>,
                             const V& centroid,
                             const std::vector<std::vector<int> >& SX,
                             const std::vector<std::vector<int> >& SY,
                             const std::vector<V>& P,
                             std::vector<std::vector<int> > & H);
 
    static void orderPointsCC(std::integral_constant<int,2>,
                              const V& centroid,
                              std::vector<int> &id,
                              const std::vector<V>& P);
    static void orderPointsCC(std::integral_constant<int,3>,
                              const V& centroid,
                              std::vector<int> &id,
                              const std::vector<V>& P);
 
    static void removeDuplicates( std::vector<int> & p);
 
    static bool newFace3D(const std::vector<int>& id,
                          const std::vector<std::vector<int> >& H);
};
 
template<class V>
inline int insertPoint(const V p, std::vector<V>& P)
{
    double eps= 1e-8;     // tolerance for identical nodes
    std::size_t k=0 ;
 
    if (P.size()>0) {
 
        while ((k<P.size())&&
               ((p - P[k]).infinity_norm()>eps*(P[k].infinity_norm()) &&
                (p - P[k]).infinity_norm()>eps*(p.infinity_norm())) &&
                !(p.infinity_norm() < eps && P[k].infinity_norm() <eps &&
                   (p - P[k]).infinity_norm() <  eps))
            k++ ;
 
        if (k>=P.size())
            P.push_back(p) ;        //  new node is not contained in P
 
    }
    else
        P.push_back(p);
 
    return k ;
}
 
 
} /* namespace Dune::GridGlue */
} /* namespace Dune */
 
#include "simplexintersection.cc"
#include "computeintersection.cc"
 
#endif