pqkfactory.hh 5.26 KB
Newer Older
1
2
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
3
4
#ifndef DUNE_LOCALFUNCTIONS_LAGRANGE_PQKFACTORY_HH
#define DUNE_LOCALFUNCTIONS_LAGRANGE_PQKFACTORY_HH
5
6
7

#include <map>

8
#include <dune/geometry/type.hh>
9
10
11
12
13
14
15

#include <dune/localfunctions/common/virtualinterface.hh>
#include <dune/localfunctions/common/virtualwrappers.hh>

#include <dune/localfunctions/lagrange/p0.hh>
#include <dune/localfunctions/lagrange/pk.hh>
#include <dune/localfunctions/lagrange/q1.hh>
16
#include <dune/localfunctions/lagrange/qk.hh>
17
18
#include <dune/localfunctions/lagrange/prismp1.hh>
#include <dune/localfunctions/lagrange/prismp2.hh>
Oliver Sander's avatar
Oliver Sander committed
19
#include <dune/localfunctions/lagrange/pyramidp1.hh>
20
#include <dune/localfunctions/lagrange/pyramidp2.hh>
21
22
23
24
25
26
27
28
29
30
31

namespace Dune
{

  /** \brief Factory that only creates dimension specific local finite elements
   *
   * Empty default implementation
   */
  template<class D, class R, int d, int k>
  struct DimSpecificPQkLocalFiniteElementFactory
  {
32
    typedef typename P0LocalFiniteElement<D,R,d>::Traits::LocalBasisType::Traits T;
33
34

    //! create finite element for given GeometryType
35
    static LocalFiniteElementVirtualInterface<T>* create(const GeometryType&)
36
    {
Oliver Sander's avatar
Oliver Sander committed
37
      return nullptr;
38
39
40
41
42
    }
  };

  /** \brief Factory that only creates dimension specific local finite elements
   *
Christian Engwer's avatar
Christian Engwer committed
43
   * Specialization for dim=3
44
45
46
47
   */
  template<class D, class R, int k>
  struct DimSpecificPQkLocalFiniteElementFactory<D,R,3,k>
  {
48
    typedef typename P0LocalFiniteElement<D,R,3>::Traits::LocalBasisType::Traits T;
49
50
    typedef PrismP1LocalFiniteElement<D,R> PrismP1;
    typedef PrismP2LocalFiniteElement<D,R> PrismP2;
51
    typedef PyramidP1LocalFiniteElement<D,R> PyramidP1;
52
    typedef PyramidP2LocalFiniteElement<D,R> PyramidP2;
53
54
55
56

    //! create finite element for given GeometryType
    static LocalFiniteElementVirtualInterface<T>* create(const GeometryType& gt)
    {
57
      if ((gt.isPrism())and (k==1))
58
        return new LocalFiniteElementVirtualImp<PrismP1>(PrismP1());
59
      if ((gt.isPrism())and (k==2))
60
        return new LocalFiniteElementVirtualImp<PrismP2>(PrismP2());
61
      if ((gt.isPyramid())and (k==1))
62
        return new LocalFiniteElementVirtualImp<PyramidP1>(PyramidP1());
63
64
      if ((gt.isPyramid())and (k==2))
        return new LocalFiniteElementVirtualImp<PyramidP2>(PyramidP2());
Oliver Sander's avatar
Oliver Sander committed
65
      return nullptr;
66
67
68
69
70
71
72
73
74
75
    }
  };


  /** \brief Factory to create any kind of Pk/Qk like element wrapped for the virtual interface
   *
   */
  template<class D, class R, int dim, int k>
  struct PQkLocalFiniteElementFactory
  {
76
    typedef typename P0LocalFiniteElement<D,R,dim>::Traits::LocalBasisType::Traits T;
77
78
79
    typedef LocalFiniteElementVirtualInterface<T> FiniteElementType;
    typedef P0LocalFiniteElement<D,R,dim> P0;
    typedef PkLocalFiniteElement<D,R,dim,k> Pk;
80
    typedef QkLocalFiniteElement<D,R,dim,k> Qk;
81
82
83
84
85
86


    //! create finite element for given GeometryType
    static FiniteElementType* create(const GeometryType& gt)
    {
      if (k==0)
87
        return new LocalFiniteElementVirtualImp<P0>(P0(gt));
88

89
      if (gt.isSimplex())
90
91
        return new LocalFiniteElementVirtualImp<Pk>(Pk());

92
93
      if (gt.isCube())
        return new LocalFiniteElementVirtualImp<Qk>(Qk());
94
95
96
97
98
99
100

      return DimSpecificPQkLocalFiniteElementFactory<D,R,dim,k>::create(gt);
    }
  };



101
  /** \brief A cache that stores all available Pk/Qk like local finite elements for the given dimension and order
102
103
104
   *
   * An interface for dealing with different vertex orders is currently missing.
   * So you can in general only use this for order=1,2 or with global DG spaces
105
106
107
108
109
   *
   * \tparam D Type used for domain coordinates
   * \tparam R Type used for shape function values
   * \tparam dim Element dimension
   * \tparam k Element order
110
111
112
113
114
   */
  template<class D, class R, int dim, int k>
  class PQkLocalFiniteElementCache
  {
  protected:
115
    typedef typename P0LocalFiniteElement<D,R,dim>::Traits::LocalBasisType::Traits T;
116
    typedef LocalFiniteElementVirtualInterface<T> FE;
Carsten Gräser's avatar
Carsten Gräser committed
117
    typedef typename std::map<GeometryType,FE*> FEMap;
118
119

  public:
120
    /** \brief Type of the finite elements stored in this cache */
121
122
    typedef FE FiniteElementType;

123
    /** \brief Default constructor */
124
125
    PQkLocalFiniteElementCache() {}

126
    /** \brief Copy constructor */
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
    PQkLocalFiniteElementCache(const PQkLocalFiniteElementCache& other)
    {
      typename FEMap::iterator it = other.cache_.begin();
      typename FEMap::iterator end = other.cache_.end();
      for(; it!=end; ++it)
        cache_[it->first] = (it->second)->clone();
    }

    ~PQkLocalFiniteElementCache()
    {
      typename FEMap::iterator it = cache_.begin();
      typename FEMap::iterator end = cache_.end();
      for(; it!=end; ++it)
        delete it->second;
    }

    //! Get local finite element for given GeometryType
    const FiniteElementType& get(const GeometryType& gt) const
    {
      typename FEMap::const_iterator it = cache_.find(gt);
      if (it==cache_.end())
      {
        FiniteElementType* fe = PQkLocalFiniteElementFactory<D,R,dim,k>::create(gt);
        if (fe==0)
          DUNE_THROW(Dune::NotImplemented,"No Pk/Qk like local finite element available for geometry type " << gt << " and order " << k);

        cache_[gt] = fe;
        return *fe;
      }
      return *(it->second);
    }

  protected:
    mutable FEMap cache_;

  };

}

#endif