8
8
#include <dolfin/Mesh.h>
9
9
#include <dolfin/Cell.h>
10
10
#include <dolfin/Point.h>
11
#include <dolfin/Vector.h>
12
11
#include <dolfin/AffineMap.h>
13
12
#include <dolfin/FiniteElement.h>
14
13
#include <dolfin/FiniteElementSpec.h>
19
class P4tet : public dolfin::FiniteElement
23
P4tet() : dolfin::FiniteElement(), tensordims(0), subelements(0)
25
// Element is scalar, don't need to initialize tensordims
27
// Element is simple, don't need to initialize subelements
32
if ( tensordims ) delete [] tensordims;
35
for (unsigned int i = 0; i < elementdim(); i++)
36
delete subelements[i];
37
delete [] subelements;
41
inline unsigned int spacedim() const
46
inline unsigned int shapedim() const
51
inline unsigned int tensordim(unsigned int i) const
53
dolfin_error("Element is scalar.");
57
inline unsigned int elementdim() const
62
inline unsigned int rank() const
67
void nodemap(int nodes[], const Cell& cell, const Mesh& mesh) const
69
static unsigned int edge_reordering_0[2][3] = {{0, 1, 2}, {2, 1, 0}};
70
static unsigned int face_reordering_0[6][3] = {{0, 1, 2}, {0, 2, 1}, {2, 0, 1}, {1, 0, 2}, {1, 2, 0}, {2, 1, 0}};
71
nodes[0] = cell.vertexID(0);
72
nodes[1] = cell.vertexID(1);
73
nodes[2] = cell.vertexID(2);
74
nodes[3] = cell.vertexID(3);
75
int alignment = cell.edgeAlignment(0);
76
int offset = mesh.numVertices();
77
nodes[4] = offset + 3*cell.edgeID(0) + edge_reordering_0[alignment][0];
78
nodes[5] = offset + 3*cell.edgeID(0) + edge_reordering_0[alignment][1];
79
nodes[6] = offset + 3*cell.edgeID(0) + edge_reordering_0[alignment][2];
80
alignment = cell.edgeAlignment(1);
81
nodes[7] = offset + 3*cell.edgeID(1) + edge_reordering_0[alignment][0];
82
nodes[8] = offset + 3*cell.edgeID(1) + edge_reordering_0[alignment][1];
83
nodes[9] = offset + 3*cell.edgeID(1) + edge_reordering_0[alignment][2];
84
alignment = cell.edgeAlignment(2);
85
nodes[10] = offset + 3*cell.edgeID(2) + edge_reordering_0[alignment][0];
86
nodes[11] = offset + 3*cell.edgeID(2) + edge_reordering_0[alignment][1];
87
nodes[12] = offset + 3*cell.edgeID(2) + edge_reordering_0[alignment][2];
88
alignment = cell.edgeAlignment(3);
89
nodes[13] = offset + 3*cell.edgeID(3) + edge_reordering_0[alignment][0];
90
nodes[14] = offset + 3*cell.edgeID(3) + edge_reordering_0[alignment][1];
91
nodes[15] = offset + 3*cell.edgeID(3) + edge_reordering_0[alignment][2];
92
alignment = cell.edgeAlignment(4);
93
nodes[16] = offset + 3*cell.edgeID(4) + edge_reordering_0[alignment][0];
94
nodes[17] = offset + 3*cell.edgeID(4) + edge_reordering_0[alignment][1];
95
nodes[18] = offset + 3*cell.edgeID(4) + edge_reordering_0[alignment][2];
96
alignment = cell.edgeAlignment(5);
97
nodes[19] = offset + 3*cell.edgeID(5) + edge_reordering_0[alignment][0];
98
nodes[20] = offset + 3*cell.edgeID(5) + edge_reordering_0[alignment][1];
99
nodes[21] = offset + 3*cell.edgeID(5) + edge_reordering_0[alignment][2];
100
alignment = cell.faceAlignment(0);
101
offset = offset + 3*mesh.numEdges();
102
nodes[22] = offset + 3*cell.faceID(0) + face_reordering_0[alignment][0];
103
nodes[23] = offset + 3*cell.faceID(0) + face_reordering_0[alignment][1];
104
nodes[24] = offset + 3*cell.faceID(0) + face_reordering_0[alignment][2];
105
alignment = cell.faceAlignment(1);
106
nodes[25] = offset + 3*cell.faceID(1) + face_reordering_0[alignment][0];
107
nodes[26] = offset + 3*cell.faceID(1) + face_reordering_0[alignment][1];
108
nodes[27] = offset + 3*cell.faceID(1) + face_reordering_0[alignment][2];
109
alignment = cell.faceAlignment(2);
110
nodes[28] = offset + 3*cell.faceID(2) + face_reordering_0[alignment][0];
111
nodes[29] = offset + 3*cell.faceID(2) + face_reordering_0[alignment][1];
112
nodes[30] = offset + 3*cell.faceID(2) + face_reordering_0[alignment][2];
113
alignment = cell.faceAlignment(3);
114
nodes[31] = offset + 3*cell.faceID(3) + face_reordering_0[alignment][0];
115
nodes[32] = offset + 3*cell.faceID(3) + face_reordering_0[alignment][1];
116
nodes[33] = offset + 3*cell.faceID(3) + face_reordering_0[alignment][2];
117
offset = offset + 3*mesh.numFaces();
118
nodes[34] = offset + cell.id();
121
void pointmap(Point points[], unsigned int components[], const AffineMap& map) const
123
points[0] = map(0.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00);
124
points[1] = map(1.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00);
125
points[2] = map(0.000000000000000e+00, 1.000000000000000e+00, 0.000000000000000e+00);
126
points[3] = map(0.000000000000000e+00, 0.000000000000000e+00, 1.000000000000000e+00);
127
points[4] = map(7.500000000000000e-01, 2.500000000000000e-01, 0.000000000000000e+00);
128
points[5] = map(5.000000000000000e-01, 5.000000000000000e-01, 0.000000000000000e+00);
129
points[6] = map(2.500000000000000e-01, 7.500000000000000e-01, 0.000000000000000e+00);
130
points[7] = map(0.000000000000000e+00, 7.500000000000000e-01, 0.000000000000000e+00);
131
points[8] = map(0.000000000000000e+00, 5.000000000000000e-01, 0.000000000000000e+00);
132
points[9] = map(0.000000000000000e+00, 2.500000000000000e-01, 0.000000000000000e+00);
133
points[10] = map(2.500000000000000e-01, 0.000000000000000e+00, 0.000000000000000e+00);
134
points[11] = map(5.000000000000000e-01, 0.000000000000000e+00, 0.000000000000000e+00);
135
points[12] = map(7.500000000000000e-01, 0.000000000000000e+00, 0.000000000000000e+00);
136
points[13] = map(0.000000000000000e+00, 0.000000000000000e+00, 2.500000000000000e-01);
137
points[14] = map(0.000000000000000e+00, 0.000000000000000e+00, 5.000000000000000e-01);
138
points[15] = map(0.000000000000000e+00, 0.000000000000000e+00, 7.500000000000000e-01);
139
points[16] = map(7.500000000000000e-01, 0.000000000000000e+00, 2.500000000000000e-01);
140
points[17] = map(5.000000000000000e-01, 0.000000000000000e+00, 5.000000000000000e-01);
141
points[18] = map(2.500000000000000e-01, 0.000000000000000e+00, 7.500000000000000e-01);
142
points[19] = map(0.000000000000000e+00, 7.500000000000000e-01, 2.500000000000000e-01);
143
points[20] = map(0.000000000000000e+00, 5.000000000000000e-01, 5.000000000000000e-01);
144
points[21] = map(0.000000000000000e+00, 2.500000000000000e-01, 7.500000000000000e-01);
145
points[22] = map(5.000000000000000e-01, 2.500000000000000e-01, 2.500000000000000e-01);
146
points[23] = map(2.500000000000000e-01, 2.500000000000000e-01, 5.000000000000000e-01);
147
points[24] = map(2.500000000000000e-01, 5.000000000000000e-01, 2.500000000000000e-01);
148
points[25] = map(0.000000000000000e+00, 5.000000000000000e-01, 2.500000000000000e-01);
149
points[26] = map(0.000000000000000e+00, 2.500000000000000e-01, 5.000000000000000e-01);
150
points[27] = map(0.000000000000000e+00, 2.500000000000000e-01, 2.500000000000000e-01);
151
points[28] = map(2.500000000000000e-01, 0.000000000000000e+00, 5.000000000000000e-01);
152
points[29] = map(5.000000000000000e-01, 0.000000000000000e+00, 2.500000000000000e-01);
153
points[30] = map(2.500000000000000e-01, 0.000000000000000e+00, 2.500000000000000e-01);
154
points[31] = map(2.500000000000000e-01, 2.500000000000000e-01, 0.000000000000000e+00);
155
points[32] = map(5.000000000000000e-01, 2.500000000000000e-01, 0.000000000000000e+00);
156
points[33] = map(2.500000000000000e-01, 5.000000000000000e-01, 0.000000000000000e+00);
157
points[34] = map(2.500000000000000e-01, 2.500000000000000e-01, 2.500000000000000e-01);
195
void vertexeval(real values[], unsigned int vertex, const real x[], const Mesh& mesh) const
197
// FIXME: Temporary fix for Lagrange elements
198
values[0] = x[vertex];
201
const FiniteElement& operator[] (unsigned int i) const
206
FiniteElement& operator[] (unsigned int i)
211
FiniteElementSpec spec() const
213
FiniteElementSpec s("Lagrange", "tetrahedron", 4);
219
unsigned int* tensordims;
220
FiniteElement** subelements;
18
class P4tet : public dolfin::FiniteElement
22
P4tet() : dolfin::FiniteElement(), tensordims(0), subelements(0)
24
// Element is scalar, don't need to initialize tensordims
26
// Element is simple, don't need to initialize subelements
31
if ( tensordims ) delete [] tensordims;
34
for (unsigned int i = 0; i < elementdim(); i++)
35
delete subelements[i];
36
delete [] subelements;
40
inline unsigned int spacedim() const
45
inline unsigned int shapedim() const
50
inline unsigned int tensordim(unsigned int i) const
52
dolfin_error("Element is scalar.");
56
inline unsigned int elementdim() const
61
inline unsigned int rank() const
66
void nodemap(int nodes[], const Cell& cell, const Mesh& mesh) const
68
static unsigned int edge_reordering_0[2][3] = {{0, 1, 2}, {2, 1, 0}};
69
static unsigned int face_reordering_0[6][3] = {{0, 1, 2}, {0, 2, 1}, {2, 0, 1}, {1, 0, 2}, {1, 2, 0}, {2, 1, 0}};
70
nodes[0] = cell.entities(0)[0];
71
nodes[1] = cell.entities(0)[1];
72
nodes[2] = cell.entities(0)[2];
73
nodes[3] = cell.entities(0)[3];
74
int alignment = cell.alignment(1, 0);
75
int offset = mesh.topology().size(0);
76
nodes[4] = offset + 3*cell.entities(1)[0] + edge_reordering_0[alignment][0];
77
nodes[5] = offset + 3*cell.entities(1)[0] + edge_reordering_0[alignment][1];
78
nodes[6] = offset + 3*cell.entities(1)[0] + edge_reordering_0[alignment][2];
79
alignment = cell.alignment(1, 1);
80
nodes[7] = offset + 3*cell.entities(1)[1] + edge_reordering_0[alignment][0];
81
nodes[8] = offset + 3*cell.entities(1)[1] + edge_reordering_0[alignment][1];
82
nodes[9] = offset + 3*cell.entities(1)[1] + edge_reordering_0[alignment][2];
83
alignment = cell.alignment(1, 2);
84
nodes[10] = offset + 3*cell.entities(1)[2] + edge_reordering_0[alignment][0];
85
nodes[11] = offset + 3*cell.entities(1)[2] + edge_reordering_0[alignment][1];
86
nodes[12] = offset + 3*cell.entities(1)[2] + edge_reordering_0[alignment][2];
87
alignment = cell.alignment(1, 3);
88
nodes[13] = offset + 3*cell.entities(1)[3] + edge_reordering_0[alignment][0];
89
nodes[14] = offset + 3*cell.entities(1)[3] + edge_reordering_0[alignment][1];
90
nodes[15] = offset + 3*cell.entities(1)[3] + edge_reordering_0[alignment][2];
91
alignment = cell.alignment(1, 4);
92
nodes[16] = offset + 3*cell.entities(1)[4] + edge_reordering_0[alignment][0];
93
nodes[17] = offset + 3*cell.entities(1)[4] + edge_reordering_0[alignment][1];
94
nodes[18] = offset + 3*cell.entities(1)[4] + edge_reordering_0[alignment][2];
95
alignment = cell.alignment(1, 5);
96
nodes[19] = offset + 3*cell.entities(1)[5] + edge_reordering_0[alignment][0];
97
nodes[20] = offset + 3*cell.entities(1)[5] + edge_reordering_0[alignment][1];
98
nodes[21] = offset + 3*cell.entities(1)[5] + edge_reordering_0[alignment][2];
99
alignment = cell.alignment(2, 0);
100
offset = offset + 3*mesh.topology().size(1);
101
nodes[22] = offset + 3*cell.entities(2)[0] + face_reordering_0[alignment][0];
102
nodes[23] = offset + 3*cell.entities(2)[0] + face_reordering_0[alignment][1];
103
nodes[24] = offset + 3*cell.entities(2)[0] + face_reordering_0[alignment][2];
104
alignment = cell.alignment(2, 1);
105
nodes[25] = offset + 3*cell.entities(2)[1] + face_reordering_0[alignment][0];
106
nodes[26] = offset + 3*cell.entities(2)[1] + face_reordering_0[alignment][1];
107
nodes[27] = offset + 3*cell.entities(2)[1] + face_reordering_0[alignment][2];
108
alignment = cell.alignment(2, 2);
109
nodes[28] = offset + 3*cell.entities(2)[2] + face_reordering_0[alignment][0];
110
nodes[29] = offset + 3*cell.entities(2)[2] + face_reordering_0[alignment][1];
111
nodes[30] = offset + 3*cell.entities(2)[2] + face_reordering_0[alignment][2];
112
alignment = cell.alignment(2, 3);
113
nodes[31] = offset + 3*cell.entities(2)[3] + face_reordering_0[alignment][0];
114
nodes[32] = offset + 3*cell.entities(2)[3] + face_reordering_0[alignment][1];
115
nodes[33] = offset + 3*cell.entities(2)[3] + face_reordering_0[alignment][2];
116
offset = offset + 3*mesh.topology().size(2);
117
nodes[34] = offset + cell.index();
120
void pointmap(Point points[], unsigned int components[], const AffineMap& map) const
122
points[0] = map(0.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00);
123
points[1] = map(1.000000000000000e+00, 0.000000000000000e+00, 0.000000000000000e+00);
124
points[2] = map(0.000000000000000e+00, 1.000000000000000e+00, 0.000000000000000e+00);
125
points[3] = map(0.000000000000000e+00, 0.000000000000000e+00, 1.000000000000000e+00);
126
points[4] = map(7.500000000000000e-01, 2.500000000000000e-01, 0.000000000000000e+00);
127
points[5] = map(5.000000000000000e-01, 5.000000000000000e-01, 0.000000000000000e+00);
128
points[6] = map(2.500000000000000e-01, 7.500000000000000e-01, 0.000000000000000e+00);
129
points[7] = map(0.000000000000000e+00, 7.500000000000000e-01, 0.000000000000000e+00);
130
points[8] = map(0.000000000000000e+00, 5.000000000000000e-01, 0.000000000000000e+00);
131
points[9] = map(0.000000000000000e+00, 2.500000000000000e-01, 0.000000000000000e+00);
132
points[10] = map(2.500000000000000e-01, 0.000000000000000e+00, 0.000000000000000e+00);
133
points[11] = map(5.000000000000000e-01, 0.000000000000000e+00, 0.000000000000000e+00);
134
points[12] = map(7.500000000000000e-01, 0.000000000000000e+00, 0.000000000000000e+00);
135
points[13] = map(0.000000000000000e+00, 0.000000000000000e+00, 2.500000000000000e-01);
136
points[14] = map(0.000000000000000e+00, 0.000000000000000e+00, 5.000000000000000e-01);
137
points[15] = map(0.000000000000000e+00, 0.000000000000000e+00, 7.500000000000000e-01);
138
points[16] = map(7.500000000000000e-01, 0.000000000000000e+00, 2.500000000000000e-01);
139
points[17] = map(5.000000000000000e-01, 0.000000000000000e+00, 5.000000000000000e-01);
140
points[18] = map(2.500000000000000e-01, 0.000000000000000e+00, 7.500000000000000e-01);
141
points[19] = map(0.000000000000000e+00, 7.500000000000000e-01, 2.500000000000000e-01);
142
points[20] = map(0.000000000000000e+00, 5.000000000000000e-01, 5.000000000000000e-01);
143
points[21] = map(0.000000000000000e+00, 2.500000000000000e-01, 7.500000000000000e-01);
144
points[22] = map(5.000000000000000e-01, 2.500000000000000e-01, 2.500000000000000e-01);
145
points[23] = map(2.500000000000000e-01, 2.500000000000000e-01, 5.000000000000000e-01);
146
points[24] = map(2.500000000000000e-01, 5.000000000000000e-01, 2.500000000000000e-01);
147
points[25] = map(0.000000000000000e+00, 5.000000000000000e-01, 2.500000000000000e-01);
148
points[26] = map(0.000000000000000e+00, 2.500000000000000e-01, 5.000000000000000e-01);
149
points[27] = map(0.000000000000000e+00, 2.500000000000000e-01, 2.500000000000000e-01);
150
points[28] = map(2.500000000000000e-01, 0.000000000000000e+00, 5.000000000000000e-01);
151
points[29] = map(5.000000000000000e-01, 0.000000000000000e+00, 2.500000000000000e-01);
152
points[30] = map(2.500000000000000e-01, 0.000000000000000e+00, 2.500000000000000e-01);
153
points[31] = map(2.500000000000000e-01, 2.500000000000000e-01, 0.000000000000000e+00);
154
points[32] = map(5.000000000000000e-01, 2.500000000000000e-01, 0.000000000000000e+00);
155
points[33] = map(2.500000000000000e-01, 5.000000000000000e-01, 0.000000000000000e+00);
156
points[34] = map(2.500000000000000e-01, 2.500000000000000e-01, 2.500000000000000e-01);
194
void vertexeval(uint vertex_nodes[], unsigned int vertex, const Mesh& mesh) const
196
// FIXME: Temporary fix for Lagrange elements
197
vertex_nodes[0] = vertex;
200
const FiniteElement& operator[] (unsigned int i) const
205
FiniteElement& operator[] (unsigned int i)
210
FiniteElementSpec spec() const
212
FiniteElementSpec s("Lagrange", "tetrahedron", 4);
218
unsigned int* tensordims;
219
FiniteElement** subelements;
added
removed
src/kernel/fem/dolfin/P4tet.h
