localrestrictprolong.hh 8.77 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13
#ifndef DUNE_FEM_LFESPACE_LOCALRESTRICTPROLONG_HH
#define DUNE_FEM_LFESPACE_LOCALRESTRICTPROLONG_HH

#include <dune/common/dynmatrix.hh>
#include <dune/fem/function/localfunction/localfunction.hh>
#include <dune/fem/space/common/localrestrictprolong.hh>

namespace Dune
{

  namespace Fem
  {

14 15
    namespace Impl
    {
16 17
      template< class LocalGeometry, class LF>
      struct FatherWrapper
18
      {
19
        typedef std::remove_reference_t< LF > LocalFunctionType;
20 21 22
        typedef std::remove_reference_t< LocalGeometry > LocalGeometryType;
        struct Traits
        {
23 24
          typedef typename LocalFunctionType::DomainType DomainType;
          typedef typename LocalFunctionType::RangeType RangeType;
25
        };
26 27 28 29 30 31 32 33 34 35 36
        typedef typename LocalFunctionType::EntityType EntityType;
        typedef typename LocalFunctionType::FunctionSpaceType FunctionSpaceType;
        typedef typename LocalFunctionType::DomainType DomainType;
        typedef typename LocalFunctionType::RangeType RangeType;
        typedef typename LocalFunctionType::JacobianRangeType JacobianRangeType;
        typedef typename LocalFunctionType::HessianRangeType HessianRangeType;

        FatherWrapper ( const LocalGeometryType &localGeo, const LocalFunctionType &lfFather)
          : localGeo_(localGeo), lfFather_(lfFather) {}
        template <class Point>
        void evaluate ( const Point &x, RangeType &y ) const
37
        {
38
          lfFather_.evaluate( localGeo_.global(x), y);
39 40 41 42
        }

      private:
        const LocalGeometryType &localGeo_;
43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94
        const LocalFunctionType &lfFather_;
      };
      template< class BasisFunctionSet, class LF>
      struct SonsWrapper
      {
        typedef std::remove_reference_t< LF > LocalFunctionType;
        typedef std::remove_reference_t< BasisFunctionSet > BasisFunctionSetType;
        struct Traits
        {
          typedef typename LocalFunctionType::DomainType DomainType;
          typedef typename LocalFunctionType::RangeType RangeType;
        };
        typedef typename LocalFunctionType::FunctionSpaceType FunctionSpaceType;
        typedef typename LocalFunctionType::DomainType DomainType;
        typedef typename LocalFunctionType::RangeType RangeType;
        typedef typename LocalFunctionType::JacobianRangeType JacobianRangeType;
        typedef typename LocalFunctionType::HessianRangeType HessianRangeType;
        typedef typename LocalFunctionType::EntityType EntityType;
        typedef typename EntityType::LocalGeometry LocalGeometryType;

        SonsWrapper (
          const std::vector< EntityType >& childEntities,
          const std::vector< BasisFunctionSetType >& childBasisSets,
          const std::vector< std::vector<double> >& childDofs )
          : childEntities_(childEntities), childBasisSets_(childBasisSets), childDofs_(childDofs)
        {}
        template <class Point>
        void evaluate ( const Point &x, RangeType &val ) const
        {
          val = RangeType(0);
          RangeType tmp;
          double weight = 0;
          for (unsigned int i=0; i<childEntities_.size();++i)
          {
            const auto &refSon = Dune::ReferenceElements< typename LocalGeometryType::ctype, LocalGeometryType::mydimension >
              ::general( childEntities_[i].type() );
            auto y = childEntities_[i].geometryInFather().local(x);
            if( refSon.checkInside( y ) )
            {
              childBasisSets_[i].evaluateAll( y, childDofs_[i], tmp );
              val += tmp;
              weight += 1.;
            }
          }
          assert( weight > 0); // weight==0 would mean that point was found in none of the children
          val /= weight;
        }

      private:
        const std::vector< EntityType >& childEntities_;
        const std::vector< BasisFunctionSetType >& childBasisSets_;
        const std::vector< std::vector<double> >& childDofs_;
95 96
      };

97 98 99 100
      // a detailed description is given in MR308
      // https://gitlab.dune-project.org/dune-fem/dune-fem/merge_requests/308
      template< class LFESpace >
      struct DefaultLocalRestrictProlongLFE
101
      {
102 103
        typedef LFESpace DiscreteFunctionSpaceType;
        typedef DefaultLocalRestrictProlongLFE< DiscreteFunctionSpaceType > ThisType;
104

105 106 107 108 109 110 111 112 113
        typedef typename DiscreteFunctionSpaceType::DomainFieldType DomainFieldType;
        typedef typename DiscreteFunctionSpaceType::BasisFunctionSetType BasisFunctionSetType;
        typedef typename DiscreteFunctionSpaceType::EntityType EntityType;
        typedef typename EntityType::LocalGeometry LocalGeometryType;
        typedef typename EntityType::EntitySeed EntitySeedType;

        DefaultLocalRestrictProlongLFE (const DiscreteFunctionSpaceType &space)
        : space_( space ), childSeeds_(0), childDofs_(0)
        {}
114

115 116 117 118 119 120 121
        /** \brief explicit set volume ratio of son and father
         *
         *  \param[in]  weight  volume of son / volume of father
         *
         *  \note If this ratio is set, it is assume to be constant.
         */
        void setFatherChildWeight ( const DomainFieldType &weight )
122 123
        {
        }
124

125 126 127 128 129 130 131 132 133 134 135 136 137
        //! restrict data to father
        template< class LFFather, class LFSon, class LocalGeometry >
        void restrictLocal ( LFFather &lfFather, const LFSon &lfSon,
                             const LocalGeometry &geometryInFather, bool initialize ) const
        {
          const int numDofs = lfFather.numDofs();
          assert( lfFather.numDofs() == lfSon.numDofs() );

          if (initialize)
          {
            childSeeds_.resize(0);
            childDofs_.resize(0);
          }
138

139 140 141 142 143 144 145 146
          childSeeds_.push_back(lfSon.entity().seed());
          childDofs_.push_back(std::vector<double>(lfSon.size()));
          for (unsigned int i=0;i<lfSon.size();++i)
            childDofs_.back()[i] = lfSon[i];
        }

        template <class LFFather>
        void restrictFinalize( LFFather &lfFather ) const
147
        {
148 149 150 151 152 153 154 155 156 157 158
          const int numDofs = lfFather.numDofs();
          std::vector< EntityType > childEntities(childSeeds_.size());
          std::vector< BasisFunctionSetType > childBasisSets(childSeeds_.size());
          for (unsigned int i=0; i<childSeeds_.size();++i)
          {
            childEntities[i] = space_.gridPart().entity( childSeeds_[i] );
            childBasisSets[i] = space_.basisFunctionSet( childEntities[i] );
          }
          space_.interpolation(lfFather.entity())
            ( Impl::SonsWrapper<BasisFunctionSetType, LFFather>( childEntities, childBasisSets, childDofs_ ),
              lfFather.localDofVector() );
159
        }
160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198

        //! prolong data to children
        template< class LFFather, class LFSon, class LocalGeometry >
        void prolongLocal ( const LFFather &lfFather, LFSon &lfSon,
                            const LocalGeometry &geometryInFather, bool initialize ) const
        {
          const int numDofs = lfFather.numDofs();
          assert( lfFather.numDofs() == lfSon.numDofs() );
          DynamicVector<double> sonDofs( numDofs );
          space_.interpolation(lfSon.entity())
            ( Impl::FatherWrapper<LocalGeometry,LFFather>(geometryInFather,lfFather),
              sonDofs );
          for (int i=0;i<numDofs;++i)
            lfSon[i] = sonDofs[i];
        }

        //! do discrete functions need a communication after restriction / prolongation?
        bool needCommunication () const { return true; }

      protected:
        template< class Entity >
        static DomainFieldType calcWeight ( const Entity &father, const Entity &son )
        {
          return son.geometry().volume() / father.geometry().volume();
        }
        const DiscreteFunctionSpaceType &space_;
        mutable std::vector< EntitySeedType > childSeeds_;
        mutable std::vector< std::vector<double> > childDofs_;
      };
    } // namespce Impl

    template< class LFEMap, class FunctionSpace, template< class > class Storage >
    struct DefaultLocalRestrictProlong< LocalFiniteElementSpace< LFEMap, FunctionSpace, Storage > >
    : public Impl::DefaultLocalRestrictProlongLFE
               < LocalFiniteElementSpace< LFEMap, FunctionSpace, Storage > >
    {
      typedef Impl::DefaultLocalRestrictProlongLFE
               < LocalFiniteElementSpace< LFEMap, FunctionSpace, Storage > > BaseType;
      using BaseType::BaseType;
199
    };
200 201 202 203
    template< class LFEMap, class FunctionSpace, template< class > class Storage >
    struct DefaultLocalRestrictProlong< DiscontinuousLocalFiniteElementSpace< LFEMap, FunctionSpace, Storage > >
    : public Impl::DefaultLocalRestrictProlongLFE
               < DiscontinuousLocalFiniteElementSpace< LFEMap, FunctionSpace, Storage > >
204

205 206 207 208 209
    {
      typedef Impl::DefaultLocalRestrictProlongLFE
               < DiscontinuousLocalFiniteElementSpace< LFEMap, FunctionSpace, Storage > > BaseType;
      using BaseType::BaseType;
    };
210 211 212 213 214
  } // namespace Fem

} // namespace Dune

#endif // #ifndef DUNE_FEM_LOCALRESTRICTPROLONG_HH