DUNE-FEM (unstable)

equidistantpoints.hh
1// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
2// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
3#ifndef DUNE_LOCALFUNCTIONS_LAGRANGE_EQUIDISTANTPOINTS_HH
4#define DUNE_LOCALFUNCTIONS_LAGRANGE_EQUIDISTANTPOINTS_HH
5
6#include <cstddef>
7
8#include <algorithm>
9#include <vector>
10
11#include <dune/geometry/referenceelements.hh>
12#include <dune/geometry/type.hh>
13
14#include <dune/localfunctions/lagrange/emptypoints.hh>
15#include <dune/localfunctions/utility/field.hh>
16
17namespace Dune
18{
19
20 // numLagrangePoints
21 // -----------------
22
23 inline std::size_t numLagrangePoints ( const GeometryType& gt, std::size_t order )
24 {
25 const int dim = gt.dim();
26 if( dim > 0 )
27 {
28 const GeometryType baseGeometryType = Impl::getBase( gt );
29 if( gt.isConical() )
30 {
31 std::size_t size = 0;
32 for( unsigned int o = 0; o <= order; ++o )
33 size += numLagrangePoints( baseGeometryType, o );
34 return size;
35 }
36 else
37 return numLagrangePoints( baseGeometryType, order ) * (order+1);
38 }
39 else
40 return 1;
41 }
42
43
44
45 // equidistantLagrangePoints
46 // -------------------------
47
48 template< class ct, unsigned int cdim >
49 inline static unsigned int equidistantLagrangePoints ( const GeometryType& gt, unsigned int codim, std::size_t order, unsigned int *count, LagrangePoint< ct, cdim > *points )
50 {
51 const unsigned int dim = gt.dim();
52 assert( (0 <= codim) && (codim <= dim) && (dim <= cdim) );
53
54 if( dim > 0 )
55 {
56 const GeometryType baseGeometryType = Impl::getBase( gt );
57 const unsigned int numBaseN = (codim < dim ? Geo::Impl::size( baseGeometryType.id(), baseGeometryType.dim(), codim ) : 0);
58 const unsigned int numBaseM = (codim > 0 ? Geo::Impl::size( baseGeometryType.id(), baseGeometryType.dim(), codim-1 ) : 0);
59
60 if( gt.isPrismatic() )
61 {
62 unsigned int size = 0;
63 if( codim < dim )
64 {
65 for( unsigned int i = 1; i < order; ++i )
66 {
67 const unsigned int n = equidistantLagrangePoints( baseGeometryType, codim, order, count, points );
68 for( unsigned int j = 0; j < n; ++j )
69 {
70 LocalKey &key = points->localKey_;
71 key = LocalKey( key.subEntity(), codim, key.index() );
72 points->point_[ dim-1 ] = ct( i ) / ct( order );
73 ++points;
74 }
75 size += n;
76 }
77 }
78
79 if( codim > 0 )
80 {
81 const unsigned int n = equidistantLagrangePoints( baseGeometryType, codim-1, order, count+numBaseN, points );
82 for( unsigned int j = 0; j < n; ++j )
83 {
84 LocalKey &key = points[ j ].localKey_;
85 key = LocalKey( key.subEntity() + numBaseN, codim, key.index() );
86
87 points[ j + n ].point_ = points[ j ].point_;
88 points[ j + n ].point_[ dim-1 ] = ct( 1 );
89 points[ j + n ].localKey_ = LocalKey( key.subEntity() + numBaseM, codim, key.index() );
90 ++count[ key.subEntity() + numBaseM ];
91 }
92 size += 2*n;
93 }
94
95 return size;
96 }
97 else
98 {
99 unsigned int size = (codim > 0 ? equidistantLagrangePoints( baseGeometryType, codim-1, order, count, points ) : 0);
100 LagrangePoint< ct, cdim > *const end = points + size;
101 for( ; points != end; ++points )
102 points->localKey_ = LocalKey( points->localKey_.subEntity(), codim, points->localKey_.index() );
103
104 if( codim < dim )
105 {
106 for( unsigned int i = order-1; i > 0; --i )
107 {
108 const unsigned int n = equidistantLagrangePoints( baseGeometryType, codim, i, count+numBaseM, points );
109 LagrangePoint< ct, cdim > *const end = points + n;
110 for( ; points != end; ++points )
111 {
112 points->localKey_ = LocalKey( points->localKey_.subEntity()+numBaseM, codim, points->localKey_.index() );
113 for( unsigned int j = 0; j < dim-1; ++j )
114 points->point_[ j ] *= ct( i ) / ct( order );
115 points->point_[ dim-1 ] = ct( order - i ) / ct( order );
116 }
117 size += n;
118 }
119 }
120 else
121 {
122 points->localKey_ = LocalKey( numBaseM, dim, count[ numBaseM ]++ );
123 points->point_ = 0;
124 points->point_[ dim-1 ] = ct( 1 );
125 ++size;
126 }
127
128 return size;
129 }
130 }
131 else
132 {
133 points->localKey_ = LocalKey( 0, 0, count[ 0 ]++ );
134 points->point_ = 0;
135 return 1;
136 }
137 }
138
139
140
141 // EquidistantPointSet
142 // -------------------
143
144 template< class F, unsigned int dim >
145 class EquidistantPointSet
146 : public EmptyPointSet< F, dim >
147 {
148 typedef EmptyPointSet< F, dim > Base;
149
150 public:
151 static const unsigned int dimension = dim;
152
153 using Base::order;
154
155 EquidistantPointSet ( std::size_t order ) : Base( order ) {}
156
157 void build ( GeometryType gt )
158 {
159 assert( gt.dim() == dimension );
160 points_.resize( numLagrangePoints( gt, order() ) );
161
162 typename Base::LagrangePoint *p = points_.data();
163 std::vector< unsigned int > count;
164 for( unsigned int mydim = 0; mydim <= dimension; ++mydim )
165 {
166 count.resize( Geo::Impl::size( gt.id(), dimension, dimension-mydim ) );
167 std::fill( count.begin(), count.end(), 0u );
168 p += equidistantLagrangePoints( gt, dimension-mydim, order(), count.data(), p );
169 }
170 const auto &refElement = referenceElement<F,dimension>(gt);
171 F weight = refElement.volume()/F(double(points_.size()));
172 for (auto &p : points_)
173 p.weight_ = weight;
174 }
175
176 template< GeometryType::Id geometryId >
177 bool build ()
178 {
179 build( GeometryType( geometryId ) );
180 return true;
181 }
182
183 bool buildCube ()
184 {
185 return build< GeometryTypes::cube(dim) > ();
186 }
187
188 static bool supports ( GeometryType, std::size_t /*order*/ ) { return true; }
189 template< GeometryType::Id geometryId>
190 static bool supports ( std::size_t order ) {
191 return supports( GeometryType( geometryId ), order );
192 }
193
194 private:
195 using Base::points_;
196 };
197
198} // namespace Dune
199
200#endif // #ifndef DUNE_LOCALFUNCTIONS_LAGRANGE_EQUIDISTANTPOINTS_HH
GeometryType
Type representing VTK's entity geometry types.
Definition: common.hh:132
bool gt(const T &first, const T &second, typename EpsilonType< T >::Type epsilon)
test if first greater than second
Definition: float_cmp.cc:158
Dune namespace.
Definition: alignedallocator.hh:13
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
A unique label for each type of element that can occur in a grid.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 21, 23:30, 2024)