DUNE-FEM (unstable)

test-fe.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
6// This header is not part of the official Dune API and might be subject
7// to change. You can use this header to test external finite element
8// implementations, but be warned that your tests might break with future
9// Dune versions.
10
11#ifndef DUNE_LOCALFUNCTIONS_TEST_TEST_FE_HH
12#define DUNE_LOCALFUNCTIONS_TEST_TEST_FE_HH
13
14#include <algorithm>
15#include <cmath>
16#include <cstddef>
17#include <cstdlib>
18#include <iomanip>
19#include <iostream>
20#include <ostream>
21#include <vector>
22
25
27
28// This class defines a local finite element function.
29// It is determined by a local finite element and
30// representing the local basis and a coefficient vector.
31// This provides the evaluate method needed by the interpolate()
32// method.
33template<class FE>
34class FEFunction
35{
36 const FE& fe;
37
38public:
39 using RangeType = typename FE::Traits::Basis::Traits::Range;
40 using DomainType = typename FE::Traits::Basis::Traits::DomainLocal;
41
42 typedef typename FE::Traits::Basis::Traits::RangeField CT;
43
44 std::vector<CT> coeff;
45
46 FEFunction(const FE& fe_) : fe(fe_) { resetCoefficients(); }
47
48 void resetCoefficients() {
49 coeff.resize(fe.basis().size());
50 for(std::size_t i=0; i<coeff.size(); ++i)
51 coeff[i] = 0;
52 }
53
54 void setRandom(double max) {
55 coeff.resize(fe.basis().size());
56 for(std::size_t i=0; i<coeff.size(); ++i)
57 coeff[i] = ((1.0*std::rand()) / RAND_MAX - 0.5)*2.0*max;
58 }
59
60 RangeType operator() (const DomainType& x) const {
61 RangeType y;
62 std::vector<RangeType> yy;
63 fe.basis().evaluateFunction(x, yy);
64 y = 0.0;
65 for (std::size_t i=0; i<yy.size(); ++i)
66 y.axpy(coeff[i], yy[i]);
67 return y;
68 }
69
70};
71
72
73// Check if interpolation is consistent with basis evaluation.
85template<class FE>
86bool testInterpolation(const FE& fe, double eps, int n=5)
87{
88 bool success = true;
89 FEFunction<FE> f(fe);
90
91 std::vector<typename FEFunction<FE>::CT> coeff;
92 for(int i=0; i<n && success; ++i) {
93 // Set random coefficient vector
94 f.setRandom(100);
95
96 // Compute interpolation weights
97
99 // Feed the shape functions to the 'interpolate' method in form of a callable.
101 fe.interpolation().interpolate(f, coeff);
102
103 // Check size of weight vector
104 if (coeff.size() != fe.basis().size()) {
105 std::cout << "Bug in LocalInterpolation for finite element type "
106 << Dune::className<FE>() << ":" << std::endl;
107 std::cout << " Interpolation vector has size " << coeff.size()
108 << std::endl;
109 std::cout << " Basis has size " << fe.basis().size() << std::endl;
110 std::cout << std::endl;
111 success = false;
112
113 // skip rest of loop since that depends on matching sizes
114 continue;
115 }
116
117 // Check if interpolation weights are equal to coefficients
118 for(std::size_t j=0; j<coeff.size() && success; ++j) {
119 if ( std::abs(coeff[j]-f.coeff[j]) >
120 (2*coeff.size()*eps)*(std::max(std::abs(f.coeff[j]), 1.0)) )
121 {
122 std::cout << std::setprecision(16);
123 std::cout << "Bug in LocalInterpolation for finite element type "
124 << Dune::className<FE>() << ":" << std::endl;
125 std::cout << " Interpolation weight " << j << " differs by "
126 << std::abs(coeff[j]-f.coeff[j]) << " from coefficient of "
127 << "linear combination." << std::endl;
128 std::cout << std::endl;
129 success = false;
130 }
131 }
132 }
133 return success;
134}
135
136// check whether Jacobian agrees with FD approximation
148template<class Geo, class FE>
149bool testJacobian(const Geo &geo, const FE& fe, double eps, double delta,
150 std::size_t order = 2)
151{
152 typedef typename FE::Traits::Basis Basis;
153
154 typedef typename Basis::Traits::DomainField DF;
155 static const std::size_t dimDLocal = Basis::Traits::dimDomainLocal;
156 typedef typename Basis::Traits::DomainLocal DomainLocal;
157 static const std::size_t dimDGlobal = Basis::Traits::dimDomainGlobal;
158
159 static const std::size_t dimR = Basis::Traits::dimRange;
160 typedef typename Basis::Traits::Range Range;
161
162 typedef typename Basis::Traits::Jacobian Jacobian;
163
164 bool success = true;
165
166 // ////////////////////////////////////////////////////////////
167 // Check the partial derivatives by comparing them
168 // to finite difference approximations
169 // ////////////////////////////////////////////////////////////
170
171 // A set of test points
174
175 // Loop over all quadrature points
176 for (std::size_t i=0; i < quad.size(); i++) {
177
178 // Get a test point
179 const DomainLocal& testPoint = quad[i].position();
180
181 // Get the shape function derivatives there
182 std::vector<Jacobian> jacobians;
183 fe.basis().evaluateJacobian(testPoint, jacobians);
184 if(jacobians.size() != fe.basis().size()) {
185 std::cout << "Bug in evaluateJacobianGlobal() for finite element type "
186 << Dune::className<FE>() << ":" << std::endl;
187 std::cout << " Jacobian vector has size " << jacobians.size()
188 << std::endl;
189 std::cout << " Basis has size " << fe.basis().size() << std::endl;
190 std::cout << std::endl;
191 return false;
192 }
193
195 geo.jacobianTransposed(testPoint);
196
197 // Loop over all shape functions in this set
198 for (std::size_t j=0; j<fe.basis().size(); ++j) {
199
200 // basis.evaluateJacobian returns global derivatives, however we can
201 // only do local derivatives, so transform the derivatives back into
202 // local coordinates
204 for(std::size_t k = 0; k < dimR; ++k)
205 for(std::size_t l = 0; l < dimDGlobal; ++l)
206 for(std::size_t m = 0; m < dimDLocal; ++m)
207 localJacobian[k][m] += jacobians[j][k][l] * geoJT[m][l];
208
209 // Loop over all local directions
210 for (std::size_t m = 0; m < dimDLocal; ++m) {
211
212 // Compute an approximation to the derivative by finite differences
213 DomainLocal upPos = testPoint;
214 DomainLocal downPos = testPoint;
215
216 upPos[m] += delta;
217 downPos[m] -= delta;
218
219 std::vector<Range> upValues, downValues;
220
221 fe.basis().evaluateFunction(upPos, upValues);
222 fe.basis().evaluateFunction(downPos, downValues);
223
224 //Loop over all components
225 for(std::size_t k = 0; k < dimR; ++k) {
226
227 // The current partial derivative, just for ease of notation
228 double derivative = localJacobian[k][m];
229
230 double finiteDiff = (upValues[j][k] - downValues[j][k]) / (2*delta);
231
232 // Check
233 if ( std::abs(derivative-finiteDiff) >
234 eps/delta*(std::max(std::abs(finiteDiff), 1.0)) )
235 {
236 std::cout << std::setprecision(16);
237 std::cout << "Bug in evaluateJacobian() for finite element type "
238 << Dune::className<FE>() << ":" << std::endl;
239 std::cout << " Shape function derivative does not agree with "
240 << "FD approximation" << std::endl;
241 std::cout << " Shape function " << j << " component " << k
242 << " at position " << testPoint << ": derivative in "
243 << "local direction " << m << " is "
244 << derivative << ", but " << finiteDiff << " is "
245 << "expected." << std::endl;
246 std::cout << std::endl;
247 success = false;
248 }
249 } //Loop over all components
250 } // Loop over all local directions
251 } // Loop over all shape functions in this set
252 } // Loop over all quadrature points
253
254 return success;
255}
256
257// call tests for given finite element
258template<class Geo, class FE>
259bool testFE(const Geo &geo, const FE& fe, double eps, double delta,
260 unsigned order = 2)
261{
262 bool success = true;
263
264 success = testInterpolation(fe, eps) and success;
265 success = testJacobian(geo, fe, eps, delta, order) and success;
266
267 return success;
268}
269
270#endif // DUNE_LOCALFUNCTIONS_TEST_TEST_FE_HH
A dense n x m matrix.
Definition: fmatrix.hh:117
Abstract base class for quadrature rules.
Definition: quadraturerules.hh:214
static const QuadratureRule & rule(const GeometryType &t, int p, QuadratureType::Enum qt=QuadratureType::GaussLegendre)
select the appropriate QuadratureRule for GeometryType t and order p
Definition: quadraturerules.hh:326
A free function to provide the demangled class name of a given object or type as a string.
Implements a matrix constructed from a given type representing a field and compile-time given number ...
constexpr auto max
Function object that returns the greater of the given values.
Definition: hybridutilities.hh:484
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 12, 23:30, 2024)