DUNE-FEM (unstable)

brezzidouglasmarini1cube2dlocalbasis.hh
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
4// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
5#ifndef DUNE_LOCALFUNCTIONS_BREZZIDOUGLASMARINI1_CUBE2D_LOCALBASIS_HH
6#define DUNE_LOCALFUNCTIONS_BREZZIDOUGLASMARINI1_CUBE2D_LOCALBASIS_HH
7
8#include <array>
9#include <bitset>
10#include <numeric>
11#include <vector>
12
14
15#include "../../common/localbasis.hh"
16
17namespace Dune
18{
28 template<class D, class R>
30 {
31
32 public:
35
38 {
39 for (size_t i=0; i<4; i++)
40 sign_[i] = 1.0;
41 }
42
48 BDM1Cube2DLocalBasis (std::bitset<4> s)
49 {
50 for (size_t i=0; i<4; i++)
51 sign_[i] = s[i] ? -1.0 : 1.0;
52 }
53
55 unsigned int size () const
56 {
57 return 8;
58 }
59
66 inline void evaluateFunction (const typename Traits::DomainType& in,
67 std::vector<typename Traits::RangeType>& out) const
68 {
69 out.resize(8);
70
71 out[0][0] = sign_[0]*(in[0] - 1.0);
72 out[0][1] = 0.0;
73 out[1][0] = 6.0*in[0]*in[1] - 3.0*in[0]-6*in[1] + 3.0;
74 out[1][1] = -3.0*in[1]*in[1] + 3.0*in[1];
75 out[2][0] = sign_[1]*(in[0]);
76 out[2][1] = 0.0;
77 out[3][0] = -6.0*in[0]*in[1] + 3.0*in[0];
78 out[3][1] = 3.0*in[1]*in[1] - 3.0*in[1];
79 out[4][0] = 0.0;
80 out[4][1] = sign_[2]*(in[1] - 1.0);
81 out[5][0] = 3.0*in[0]*in[0] - 3.0*in[0];
82 out[5][1] = -6.0*in[0]*in[1] + 6.0*in[0] + 3.0*in[1] - 3.0;
83 out[6][0] = 0.0;
84 out[6][1] = sign_[3]*(in[1]);
85 out[7][0] = -3.0*in[0]*in[0] + 3.0*in[0];
86 out[7][1] = 6.0*in[0]*in[1] - 3.0*in[1];
87 }
88
95 inline void evaluateJacobian (const typename Traits::DomainType& in,
96 std::vector<typename Traits::JacobianType>& out) const
97 {
98 out.resize(8);
99
100 out[0][0][0] = sign_[0];
101 out[0][0][1] = 0.0;
102 out[0][1][0] = 0.0;
103 out[0][1][1] = 0.0;
104
105 out[1][0][0] = 6.0*in[1] - 3.0;
106 out[1][0][1] = 6.0*in[0] - 6.0;
107 out[1][1][0] = 0.0;
108 out[1][1][1] = -6.0*in[1] + 3.0;
109
110 out[2][0][0] = sign_[1];
111 out[2][0][1] = 0.0;
112 out[2][1][0] = 0.0;
113 out[2][1][1] = 0.0;
114
115 out[3][0][0] = -6.0*in[1] + 3.0;
116 out[3][0][1] = -6.0*in[0];
117 out[3][1][0] = 0.0;
118 out[3][1][1] = 6.0*in[1] - 3.0;
119
120 out[4][0][0] = 0.0;
121 out[4][0][1] = 0.0;
122 out[4][1][0] = 0.0;
123 out[4][1][1] = sign_[2];
124
125 out[5][0][0] = 6.0*in[0] - 3.0;
126 out[5][0][1] = 0.0;
127 out[5][1][0] = -6.0*in[1] + 6.0;
128 out[5][1][1] = -6.0*in[0] + 3.0;
129
130 out[6][0][0] = 0.0;
131 out[6][0][1] = 0.0;
132 out[6][1][0] = 0.0;
133 out[6][1][1] = sign_[3];
134
135 out[7][0][0] = -6.0*in[0] + 3.0;
136 out[7][0][1] = 0.0;
137 out[7][1][0] = 6.0*in[1];
138 out[7][1][1] = 6.0*in[0] - 3.0;
139 }
140
142 void partial (const std::array<unsigned int, 2>& order,
143 const typename Traits::DomainType& in, // position
144 std::vector<typename Traits::RangeType>& out) const // return value
145 {
146 auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
147 if (totalOrder == 0) {
148 evaluateFunction(in, out);
149 } else if (totalOrder == 1) {
150 out.resize(size());
151 auto const direction = std::distance(order.begin(), std::find(order.begin(), order.end(), 1));
152
153 switch (direction) {
154 case 0:
155 out[0][0] = sign_[0];
156 out[0][1] = 0.0;
157
158 out[1][0] = 6.0*in[1] - 3.0;
159 out[1][1] = 0.0;
160
161 out[2][0] = sign_[1];
162 out[2][1] = 0.0;
163
164 out[3][0] = -6.0*in[1] + 3.0;
165 out[3][1] = 0.0;
166
167 out[4][0] = 0.0;
168 out[4][1] = 0.0;
169
170 out[5][0] = 6.0*in[0] - 3.0;
171 out[5][1] = -6.0*in[1] + 6.0;
172
173 out[6][0] = 0.0;
174 out[6][1] = 0.0;
175
176 out[7][0] = -6.0*in[0] + 3.0;
177 out[7][1] = 6.0*in[1];
178 break;
179 case 1:
180 out[0][0] = 0.0;
181 out[0][1] = 0.0;
182
183 out[1][0] = 6.0*in[0] - 6.0;
184 out[1][1] = -6.0*in[1] + 3.0;
185
186 out[2][0] = 0.0;
187 out[2][1] = 0.0;
188
189 out[3][0] = -6.0*in[0];
190 out[3][1] = 6.0*in[1] - 3.0;
191
192 out[4][0] = 0.0;
193 out[4][1] = sign_[2];
194
195 out[5][0] = 0.0;
196 out[5][1] = -6.0*in[0] + 3.0;
197
198 out[6][0] = 0.0;
199 out[6][1] = sign_[3];
200
201 out[7][0] = 0.0;
202 out[7][1] = 6.0*in[0] - 3.0;
203 break;
204 default:
205 DUNE_THROW(RangeError, "Component out of range.");
206 }
207 } else {
208 DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
209 }
210 }
211
213 unsigned int order () const
214 {
215 return 2;
216 }
217
218 private:
219 std::array<R,4> sign_;
220 };
221}
222#endif // DUNE_LOCALFUNCTIONS_BREZZIDOUGLASMARINI1_CUBE2D_LOCALBASIS_HH
First order Brezzi-Douglas-Marini shape functions on the reference quadrilateral.
Definition: brezzidouglasmarini1cube2dlocalbasis.hh:30
void evaluateJacobian(const typename Traits::DomainType &in, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: brezzidouglasmarini1cube2dlocalbasis.hh:95
void partial(const std::array< unsigned int, 2 > &order, const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate partial derivatives of all shape functions.
Definition: brezzidouglasmarini1cube2dlocalbasis.hh:142
BDM1Cube2DLocalBasis()
Standard constructor.
Definition: brezzidouglasmarini1cube2dlocalbasis.hh:37
unsigned int size() const
number of shape functions
Definition: brezzidouglasmarini1cube2dlocalbasis.hh:55
BDM1Cube2DLocalBasis(std::bitset< 4 > s)
Make set number s, where 0 <= s < 16.
Definition: brezzidouglasmarini1cube2dlocalbasis.hh:48
unsigned int order() const
Polynomial order of the shape functions.
Definition: brezzidouglasmarini1cube2dlocalbasis.hh:213
void evaluateFunction(const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: brezzidouglasmarini1cube2dlocalbasis.hh:66
A dense n x m matrix.
Definition: fmatrix.hh:117
vector space out of a tensor product of fields.
Definition: fvector.hh:91
Default exception for dummy implementations.
Definition: exceptions.hh:263
Default exception class for range errors.
Definition: exceptions.hh:254
Implements a matrix constructed from a given type representing a field and compile-time given number ...
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
constexpr T accumulate(Range &&range, T value, F &&f)
Accumulate values.
Definition: hybridutilities.hh:279
Dune namespace.
Definition: alignedallocator.hh:13
Type traits for LocalBasisVirtualInterface.
Definition: localbasis.hh:35
D DomainType
domain type
Definition: localbasis.hh:43
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 24, 23:30, 2024)