conbase_fork.hh

#ifndef CONBASE_FORK_HH
#define CONBASE_FORK_HH
 
#include<dune/pdelab/boilerplate/pdelab.hh>
 
namespace Dune {
    namespace PDELab {
 
    template<typename P, typename GFS, typename CG, bool isFunction>
    struct ConstraintsAssemblerHelperExterior
    {
 
      static void
      assemble(const P& p, const GFS& gfs, CG& cg, const bool verbose)
      {
        // get some types
        using ES = typename GFS::Traits::EntitySet;
        //using Element = typename ES::Traits::Element;
        using Intersection = typename ES::Traits::Intersection;
 
        ES es = gfs.entitySet();
 
        // make local function space
        using LFS = LocalFunctionSpace<GFS>;
        LFS lfs_e(gfs);
        LFSIndexCache<LFS> lfs_cache_e(lfs_e);
        LFS lfs_f(gfs);
        LFSIndexCache<LFS> lfs_cache_f(lfs_f);
 
        // get index set
        //auto& is = es.indexSet();
 
        // loop once over the grid
        for (const auto& element : elements(es))
        {
 
          //auto id = is.uniqueIndex(element);
 
          // bind local function space to element
          lfs_e.bind(element);
 
          using CL = typename CG::LocalTransformation;
 
          CL cl_self;
 
          //using ElementWrapper = ElementGeometry<Element>;
          using IntersectionWrapper = IntersectionGeometry<Intersection>;
 
          //TypeTree::applyToTreePair(p,lfs_e,VolumeConstraints<ElementWrapper,CL>(ElementWrapper(element),cl_self));
 
          // iterate over intersections and call metaprogram
          unsigned int intersection_index = 0;
          for (const auto& intersection : intersections(es,element))
          {
 
            auto intersection_data = classifyIntersection(es,intersection);
            auto intersection_type = std::get<0>(intersection_data);
            //auto& outside_element = std::get<1>(intersection_data);
 
            switch (intersection_type) {
 
            case IntersectionType::skeleton:
            case IntersectionType::periodic:
              {
                /*auto idn = is.uniqueIndex(outside_element);
 
                if(id > idn){
                  // bind local function space to element in neighbor
                  lfs_f.bind(outside_element);
 
                  CL cl_neighbor;
 
                  TypeTree::applyToTreePair(lfs_e,lfs_f,SkeletonConstraints<IntersectionWrapper,CL>(IntersectionWrapper(intersection,intersection_index),cl_self,cl_neighbor));
 
                  if (!cl_neighbor.empty())
                    {
                      lfs_cache_f.update();
                      cg.import_local_transformation(cl_neighbor,lfs_cache_f);
                    }
 
                }*/
                break;
              }
 
            case IntersectionType::boundary:
              //if (intersection.boundary())
                TypeTree::applyToTreePair(p,lfs_e,BoundaryConstraints<IntersectionWrapper,CL>(IntersectionWrapper(intersection,intersection_index),cl_self));
              break;
 
            case IntersectionType::processor:
              //TypeTree::applyToTree(lfs_e,ProcessorConstraints<IntersectionWrapper,CL>(IntersectionWrapper(intersection,intersection_index),cl_self));
              break;
 
            }
            ++intersection_index;
          }
 
          if (!cl_self.empty())
            {
              lfs_cache_e.update();
              cg.import_local_transformation(cl_self,lfs_cache_e);
            }
 
        }
 
        // print result
        if(verbose){
          std::cout << "constraints:" << std::endl;
 
          std::cout << cg.size() << " constrained degrees of freedom" << std::endl;
 
          for (const auto& col : cg)
          {
            std::cout << col.first << ": "; // col index
            for (const auto& row : col.second)
              std::cout << "(" << row.first << "," << row.second << ") "; // row index , value
            std::cout << std::endl;
          }
        }
      }
    }; // end ConstraintsAssemblerHelper
 
    /*template<typename GFS, typename CG>
    void constraints_exterior(const GFS& gfs, CG& cg,
                     const bool verbose = false)
    {
      NoConstraintsParameters p;
      ConstraintsAssemblerHelperExterior<NoConstraintsParameters, GFS, CG, false>::assemble(p,gfs,cg,verbose);
    }*/
    template<typename P, typename GFS, typename CG>
    void constraints_exterior(const P& p, const GFS& gfs, CG& cg,
                     const bool verbose = false)
    {
      // clear global constraints
      cg.clear();
      ConstraintsAssemblerHelperExterior<P, GFS, CG, IsGridFunction<P>::value>::assemble(p,gfs,cg,verbose);
    }
 
        template<typename Grid, unsigned int degree, Dune::GeometryType::BasicType gt,typename BCType, typename GV = typename Grid::LeafGridView>
        class CGCONBaseExterior//<Grid,degree,gt,MeshType::conforming,SolverCategory::overlapping,BCType,GV>
        {
        public:
            typedef ConformingDirichletConstraints CON;
 
            CGCONBaseExterior (Grid& grid, const BCType& bctype, const GV& gv)
            {
                conp = std::shared_ptr<CON>(new CON());
            }
 
            CGCONBaseExterior (Grid& grid, const BCType& bctype)
            {
                conp = std::shared_ptr<CON>(new CON());
            }
 
            template<typename GFS>
            void postGFSHook (const GFS& gfs) {}
            CON& getCON() {return *conp;}
            const CON& getCON() const {return *conp;}
            template<typename GFS, typename DOF>
            void make_consistent (const GFS& gfs, DOF& x) const
            {
                // make vector consistent; this is needed for all overlapping solvers
                istl::ParallelHelper<GFS> helper(gfs);
                helper.maskForeignDOFs(Backend::native(x));
                Dune::PDELab::AddDataHandle<GFS,DOF> adddh(gfs,x);
                if (gfs.gridView().comm().size()>1)
                    gfs.gridView().communicate(adddh,Dune::InteriorBorder_All_Interface,Dune::ForwardCommunication);
            }
        private:
            std::shared_ptr<CON> conp;
        };
    }
}
 
#endif