Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
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
#include <config.h>
#include <iostream>
#include <dune/common/fvector.hh>
#include <dune/grid/common/quadraturerules.hh>
#include <dune/src/rotation.hh>
#include <dune/src/localgeodesicfefunction.hh>
using namespace Dune;
int main()
{
typedef Rotation<3,double> TargetSpace;
const int dim = 2;
FieldVector<double,3> xAxis(0);
xAxis[0] = 1;
FieldVector<double,3> yAxis(0);
yAxis[1] = 1;
FieldVector<double,3> zAxis(0);
zAxis[2] = 1;
TargetSpace v0 = Rotation<3,double>(xAxis,1);
TargetSpace v1 = Rotation<3,double>(yAxis,2);
TargetSpace v2 = Rotation<3,double>(zAxis,3);
std::vector<TargetSpace> coef0(3), coef1(3), coef2(3);
coef0[0] = v0; coef0[1] = v1; coef0[2] = v2;
coef1[0] = v2; coef1[1] = v0; coef1[2] = v1;
coef2[0] = v1; coef2[1] = v2; coef2[2] = v0;
GeometryType triangle;
triangle.makeTriangle();
LocalGeodesicFEFunction<2,double,TargetSpace> f0(triangle, coef0);
LocalGeodesicFEFunction<2,double,TargetSpace> f1(triangle, coef1);
LocalGeodesicFEFunction<2,double,TargetSpace> f2(triangle, coef2);
int quadOrder = 1;
const Dune::QuadratureRule<double, dim>& quad
= Dune::QuadratureRules<double, dim>::rule(triangle, quadOrder);
for (size_t pt=0; pt<quad.size(); pt++) {
const Dune::FieldVector<double,dim>& quadPos = quad[pt].position();
Dune::FieldVector<double,dim> l0 = quadPos;
Dune::FieldVector<double,dim> l1, l2;
l1[0] = 1-quadPos[0]-quadPos[1];
l1[1] = quadPos[0];
l2[0] = quadPos[1];
l2[1] = 1-quadPos[0]-quadPos[1];
std::cout << f0.evaluate(l0) << std::endl;
std::cout << f1.evaluate(l1) << std::endl;
std::cout << f2.evaluate(l2) << std::endl;
}
}