« back to all changes in this revision
Viewing changes to test/unit/function/cpp/Projection.h
- browse files at revision 7386
- compare with another revision
- download diff
- download tarball
- view history from revision 7386
- Committer: Johan Hake
- Date: 2013-03-27 15:18:08 UTC
- mfrom: (7352.1.227 working)
- Revision ID: hake.dev@gmail.com-20130327151808-b73d4pueq1n432hg
Merge with trunk
- files added:
- cmake/scripts/download-demo-data
- files removed:
- bench/fem/assembly/cpp/forms/Elasticity3D.h
- data/graphs
- data/meshes
- site-packages/dolfin_utils/wrappers
- files renamed:
- scripts/dolfin-order/mesh0.xml => scripts/dolfin-order/mesh0.xml.gz
- scripts/dolfin-order/mesh1.xml => scripts/dolfin-order/mesh1.xml.gz
- scripts/dolfin-plot/mesh.xml => scripts/dolfin-plot/mesh.xml.gz
- files modified:
- CMakeLists.txt
added
removed
test/unit/function/cpp/Projection.h
1
// This code conforms with the UFC specification version 2.1.0+
2
// and was automatically generated by FFC version 1.1.0+.
3
//
4
// This code was generated with the option '-l dolfin' and
5
// contains DOLFIN-specific wrappers that depend on DOLFIN.
6
//
7
// This code was generated with the following parameters:
8
//
9
// cache_dir: ''
10
// convert_exceptions_to_warnings: False
11
// cpp_optimize: False
12
// cpp_optimize_flags: '-O2'
13
// epsilon: 1e-14
14
// error_control: False
15
// form_postfix: True
16
// format: 'dolfin'
17
// log_level: 10
18
// log_prefix: ''
19
// no_ferari: True
20
// optimize: True
21
// output_dir: '.'
22
// precision: 15
23
// quadrature_degree: 'auto'
24
// quadrature_rule: 'auto'
25
// representation: 'auto'
26
// split: False
27
// swig_binary: 'swig'
28
// swig_path: ''
29
30
#ifndef __PROJECTION_H
31
#define __PROJECTION_H
32
33
#include <cmath>
34
#include <stdexcept>
35
#include <fstream>
36
#include <ufc.h>
37
38
/// This class defines the interface for a finite element.
39
40
class projection_finite_element_0: public ufc::finite_element
41
{
42
public:
43
44
/// Constructor
45
projection_finite_element_0() : ufc::finite_element()
46
{
47
// Do nothing
48
}
49
50
/// Destructor
51
virtual ~projection_finite_element_0()
52
{
53
// Do nothing
54
}
55
56
/// Return a string identifying the finite element
57
virtual const char* signature() const
58
{
59
return "FiniteElement('Lagrange', Domain(Cell('tetrahedron', 3), 'tetrahedron_multiverse', 3, 3), 2, None)";
60
}
61
62
/// Return the cell shape
63
virtual ufc::shape cell_shape() const
64
{
65
return ufc::tetrahedron;
66
}
67
68
/// Return the topological dimension of the cell shape
69
virtual std::size_t topological_dimension() const
70
{
71
return 3;
72
}
73
74
/// Return the geometric dimension of the cell shape
75
virtual std::size_t geometric_dimension() const
76
{
77
return 3;
78
}
79
80
/// Return the dimension of the finite element function space
81
virtual std::size_t space_dimension() const
82
{
83
return 10;
84
}
85
86
/// Return the rank of the value space
87
virtual std::size_t value_rank() const
88
{
89
return 0;
90
}
91
92
/// Return the dimension of the value space for axis i
93
virtual std::size_t value_dimension(std::size_t i) const
94
{
95
return 1;
96
}
97
98
/// Evaluate basis function i at given point x in cell
99
virtual void evaluate_basis(std::size_t i,
100
double* values,
101
const double* x,
102
const double* vertex_coordinates,
103
int cell_orientation) const
104
{
105
// Compute Jacobian
106
double J[9];
107
compute_jacobian_tetrahedron_3d(J, vertex_coordinates);
108
109
// Compute Jacobian inverse and determinant
110
double K[9];
111
double detJ;
112
compute_jacobian_inverse_tetrahedron_3d(K, detJ, J);
113
114
115
// Compute constants
116
const double C0 = vertex_coordinates[9] + vertex_coordinates[6] + vertex_coordinates[3] - vertex_coordinates[0];
117
const double C1 = vertex_coordinates[10] + vertex_coordinates[7] + vertex_coordinates[4] - vertex_coordinates[1];
118
const double C2 = vertex_coordinates[11] + vertex_coordinates[8] + vertex_coordinates[5] - vertex_coordinates[2];
119
120
// Compute subdeterminants
121
const double d_00 = J[4]*J[8] - J[5]*J[7];
122
const double d_01 = J[5]*J[6] - J[3]*J[8];
123
const double d_02 = J[3]*J[7] - J[4]*J[6];
124
const double d_10 = J[2]*J[7] - J[1]*J[8];
125
const double d_11 = J[0]*J[8] - J[2]*J[6];
126
const double d_12 = J[1]*J[6] - J[0]*J[7];
127
const double d_20 = J[1]*J[5] - J[2]*J[4];
128
const double d_21 = J[2]*J[3] - J[0]*J[5];
129
const double d_22 = J[0]*J[4] - J[1]*J[3];
130
131
// Get coordinates and map to the reference (FIAT) element
132
double X = (d_00*(2.0*x[0] - C0) + d_10*(2.0*x[1] - C1) + d_20*(2.0*x[2] - C2)) / detJ;
133
double Y = (d_01*(2.0*x[0] - C0) + d_11*(2.0*x[1] - C1) + d_21*(2.0*x[2] - C2)) / detJ;
134
double Z = (d_02*(2.0*x[0] - C0) + d_12*(2.0*x[1] - C1) + d_22*(2.0*x[2] - C2)) / detJ;
135
136
137
// Reset values
138
*values = 0.0;
139
switch (i)
140
{
141
case 0:
142
{
143
144
// Array of basisvalues
145
double basisvalues[10] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
146
147
// Declare helper variables
148
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
149
double tmp1 = 0.25*(Y + Z)*(Y + Z);
150
double tmp2 = 0.5*(1.0 + Z + 2.0*Y);
151
double tmp3 = 0.5*(1.0 - Z);
152
double tmp4 = tmp3*tmp3;
153
154
// Compute basisvalues
155
basisvalues[0] = 1.0;
156
basisvalues[1] = tmp0;
157
basisvalues[4] = 1.5*tmp0*basisvalues[1] - 0.5*tmp1*basisvalues[0];
158
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
159
basisvalues[5] = (0.5*(2.0 + 3.0*Y + Z) + 1.0*(1.0 + Y))*basisvalues[1];
160
basisvalues[7] = (1.66666666666667*tmp2 + 0.111111111111111*tmp3)*basisvalues[2] - 0.555555555555556*tmp4*basisvalues[0];
161
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
162
basisvalues[8] = (3.0*Z + 2.0)*basisvalues[2];
163
basisvalues[6] = (3.0*Z + 2.0)*basisvalues[1];
164
basisvalues[9] = basisvalues[3]*(0.3125 + 1.875*Z) - 0.5625*basisvalues[0];
165
basisvalues[0] *= std::sqrt(0.75);
166
basisvalues[3] *= std::sqrt(1.25);
167
basisvalues[9] *= std::sqrt(1.75);
168
basisvalues[2] *= std::sqrt(2.5);
169
basisvalues[8] *= std::sqrt(3.5);
170
basisvalues[7] *= std::sqrt(5.25);
171
basisvalues[1] *= std::sqrt(7.5);
172
basisvalues[6] *= std::sqrt(10.5);
173
basisvalues[5] *= std::sqrt(15.75);
174
basisvalues[4] *= std::sqrt(26.25);
175
176
// Table(s) of coefficients
177
static const double coefficients0[10] = \
178
{-0.0577350269189625, -0.0608580619450185, -0.0351364184463153, -0.0248451997499977, 0.0650600048632355, 0.050395263067897, 0.0411475599898912, 0.0290957186981323, 0.0237565548366599, 0.0167984210226323};
179
180
// Compute value(s)
181
for (unsigned int r = 0; r < 10; r++)
182
{
183
*values += coefficients0[r]*basisvalues[r];
184
}// end loop over 'r'
185
break;
186
}
187
case 1:
188
{
189
190
// Array of basisvalues
191
double basisvalues[10] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
192
193
// Declare helper variables
194
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
195
double tmp1 = 0.25*(Y + Z)*(Y + Z);
196
double tmp2 = 0.5*(1.0 + Z + 2.0*Y);
197
double tmp3 = 0.5*(1.0 - Z);
198
double tmp4 = tmp3*tmp3;
199
200
// Compute basisvalues
201
basisvalues[0] = 1.0;
202
basisvalues[1] = tmp0;
203
basisvalues[4] = 1.5*tmp0*basisvalues[1] - 0.5*tmp1*basisvalues[0];
204
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
205
basisvalues[5] = (0.5*(2.0 + 3.0*Y + Z) + 1.0*(1.0 + Y))*basisvalues[1];
206
basisvalues[7] = (1.66666666666667*tmp2 + 0.111111111111111*tmp3)*basisvalues[2] - 0.555555555555556*tmp4*basisvalues[0];
207
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
208
basisvalues[8] = (3.0*Z + 2.0)*basisvalues[2];
209
basisvalues[6] = (3.0*Z + 2.0)*basisvalues[1];
210
basisvalues[9] = basisvalues[3]*(0.3125 + 1.875*Z) - 0.5625*basisvalues[0];
211
basisvalues[0] *= std::sqrt(0.75);
212
basisvalues[3] *= std::sqrt(1.25);
213
basisvalues[9] *= std::sqrt(1.75);
214
basisvalues[2] *= std::sqrt(2.5);
215
basisvalues[8] *= std::sqrt(3.5);
216
basisvalues[7] *= std::sqrt(5.25);
217
basisvalues[1] *= std::sqrt(7.5);
218
basisvalues[6] *= std::sqrt(10.5);
219
basisvalues[5] *= std::sqrt(15.75);
220
basisvalues[4] *= std::sqrt(26.25);
221
222
// Table(s) of coefficients
223
static const double coefficients0[10] = \
224
{-0.0577350269189625, 0.0608580619450185, -0.0351364184463153, -0.0248451997499977, 0.0650600048632355, -0.050395263067897, -0.0411475599898912, 0.0290957186981323, 0.0237565548366599, 0.0167984210226323};
225
226
// Compute value(s)
227
for (unsigned int r = 0; r < 10; r++)
228
{
229
*values += coefficients0[r]*basisvalues[r];
230
}// end loop over 'r'
231
break;
232
}
233
case 2:
234
{
235
236
// Array of basisvalues
237
double basisvalues[10] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
238
239
// Declare helper variables
240
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
241
double tmp1 = 0.25*(Y + Z)*(Y + Z);
242
double tmp2 = 0.5*(1.0 + Z + 2.0*Y);
243
double tmp3 = 0.5*(1.0 - Z);
244
double tmp4 = tmp3*tmp3;
245
246
// Compute basisvalues
247
basisvalues[0] = 1.0;
248
basisvalues[1] = tmp0;
249
basisvalues[4] = 1.5*tmp0*basisvalues[1] - 0.5*tmp1*basisvalues[0];
250
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
251
basisvalues[5] = (0.5*(2.0 + 3.0*Y + Z) + 1.0*(1.0 + Y))*basisvalues[1];
252
basisvalues[7] = (1.66666666666667*tmp2 + 0.111111111111111*tmp3)*basisvalues[2] - 0.555555555555556*tmp4*basisvalues[0];
253
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
254
basisvalues[8] = (3.0*Z + 2.0)*basisvalues[2];
255
basisvalues[6] = (3.0*Z + 2.0)*basisvalues[1];
256
basisvalues[9] = basisvalues[3]*(0.3125 + 1.875*Z) - 0.5625*basisvalues[0];
257
basisvalues[0] *= std::sqrt(0.75);
258
basisvalues[3] *= std::sqrt(1.25);
259
basisvalues[9] *= std::sqrt(1.75);
260
basisvalues[2] *= std::sqrt(2.5);
261
basisvalues[8] *= std::sqrt(3.5);
262
basisvalues[7] *= std::sqrt(5.25);
263
basisvalues[1] *= std::sqrt(7.5);
264
basisvalues[6] *= std::sqrt(10.5);
265
basisvalues[5] *= std::sqrt(15.75);
266
basisvalues[4] *= std::sqrt(26.25);
267
268
// Table(s) of coefficients
269
static const double coefficients0[10] = \
270
{-0.0577350269189626, 0.0, 0.0702728368926307, -0.0248451997499977, 0.0, 0.0, 0.0, 0.0872871560943969, -0.0475131096733199, 0.0167984210226323};
271
272
// Compute value(s)
273
for (unsigned int r = 0; r < 10; r++)
274
{
275
*values += coefficients0[r]*basisvalues[r];
276
}// end loop over 'r'
277
break;
278
}
279
case 3:
280
{
281
282
// Array of basisvalues
283
double basisvalues[10] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
284
285
// Declare helper variables
286
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
287
double tmp1 = 0.25*(Y + Z)*(Y + Z);
288
double tmp2 = 0.5*(1.0 + Z + 2.0*Y);
289
double tmp3 = 0.5*(1.0 - Z);
290
double tmp4 = tmp3*tmp3;
291
292
// Compute basisvalues
293
basisvalues[0] = 1.0;
294
basisvalues[1] = tmp0;
295
basisvalues[4] = 1.5*tmp0*basisvalues[1] - 0.5*tmp1*basisvalues[0];
296
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
297
basisvalues[5] = (0.5*(2.0 + 3.0*Y + Z) + 1.0*(1.0 + Y))*basisvalues[1];
298
basisvalues[7] = (1.66666666666667*tmp2 + 0.111111111111111*tmp3)*basisvalues[2] - 0.555555555555556*tmp4*basisvalues[0];
299
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
300
basisvalues[8] = (3.0*Z + 2.0)*basisvalues[2];
301
basisvalues[6] = (3.0*Z + 2.0)*basisvalues[1];
302
basisvalues[9] = basisvalues[3]*(0.3125 + 1.875*Z) - 0.5625*basisvalues[0];
303
basisvalues[0] *= std::sqrt(0.75);
304
basisvalues[3] *= std::sqrt(1.25);
305
basisvalues[9] *= std::sqrt(1.75);
306
basisvalues[2] *= std::sqrt(2.5);
307
basisvalues[8] *= std::sqrt(3.5);
308
basisvalues[7] *= std::sqrt(5.25);
309
basisvalues[1] *= std::sqrt(7.5);
310
basisvalues[6] *= std::sqrt(10.5);
311
basisvalues[5] *= std::sqrt(15.75);
312
basisvalues[4] *= std::sqrt(26.25);
313
314
// Table(s) of coefficients
315
static const double coefficients0[10] = \
316
{-0.0577350269189626, 0.0, 0.0, 0.074535599249993, 0.0, 0.0, 0.0, 0.0, 0.0, 0.100790526135794};
317
318
// Compute value(s)
319
for (unsigned int r = 0; r < 10; r++)
320
{
321
*values += coefficients0[r]*basisvalues[r];
322
}// end loop over 'r'
323
break;
324
}
325
case 4:
326
{
327
328
// Array of basisvalues
329
double basisvalues[10] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
330
331
// Declare helper variables
332
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
333
double tmp1 = 0.25*(Y + Z)*(Y + Z);
334
double tmp2 = 0.5*(1.0 + Z + 2.0*Y);
335
double tmp3 = 0.5*(1.0 - Z);
336
double tmp4 = tmp3*tmp3;
337
338
// Compute basisvalues
339
basisvalues[0] = 1.0;
340
basisvalues[1] = tmp0;
341
basisvalues[4] = 1.5*tmp0*basisvalues[1] - 0.5*tmp1*basisvalues[0];
342
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
343
basisvalues[5] = (0.5*(2.0 + 3.0*Y + Z) + 1.0*(1.0 + Y))*basisvalues[1];
344
basisvalues[7] = (1.66666666666667*tmp2 + 0.111111111111111*tmp3)*basisvalues[2] - 0.555555555555556*tmp4*basisvalues[0];
345
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
346
basisvalues[8] = (3.0*Z + 2.0)*basisvalues[2];
347
basisvalues[6] = (3.0*Z + 2.0)*basisvalues[1];
348
basisvalues[9] = basisvalues[3]*(0.3125 + 1.875*Z) - 0.5625*basisvalues[0];
349
basisvalues[0] *= std::sqrt(0.75);
350
basisvalues[3] *= std::sqrt(1.25);
351
basisvalues[9] *= std::sqrt(1.75);
352
basisvalues[2] *= std::sqrt(2.5);
353
basisvalues[8] *= std::sqrt(3.5);
354
basisvalues[7] *= std::sqrt(5.25);
355
basisvalues[1] *= std::sqrt(7.5);
356
basisvalues[6] *= std::sqrt(10.5);
357
basisvalues[5] *= std::sqrt(15.75);
358
basisvalues[4] *= std::sqrt(26.25);
359
360
// Table(s) of coefficients
361
static const double coefficients0[10] = \
362
{0.23094010767585, 0.0, 0.140545673785261, 0.0993807989999907, 0.0, 0.0, 0.0, 0.0, 0.1187827741833, -0.0671936840905293};
363
364
// Compute value(s)
365
for (unsigned int r = 0; r < 10; r++)
366
{
367
*values += coefficients0[r]*basisvalues[r];
368
}// end loop over 'r'
369
break;
370
}
371
case 5:
372
{
373
374
// Array of basisvalues
375
double basisvalues[10] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
376
377
// Declare helper variables
378
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
379
double tmp1 = 0.25*(Y + Z)*(Y + Z);
380
double tmp2 = 0.5*(1.0 + Z + 2.0*Y);
381
double tmp3 = 0.5*(1.0 - Z);
382
double tmp4 = tmp3*tmp3;
383
384
// Compute basisvalues
385
basisvalues[0] = 1.0;
386
basisvalues[1] = tmp0;
387
basisvalues[4] = 1.5*tmp0*basisvalues[1] - 0.5*tmp1*basisvalues[0];
388
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
389
basisvalues[5] = (0.5*(2.0 + 3.0*Y + Z) + 1.0*(1.0 + Y))*basisvalues[1];
390
basisvalues[7] = (1.66666666666667*tmp2 + 0.111111111111111*tmp3)*basisvalues[2] - 0.555555555555556*tmp4*basisvalues[0];
391
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
392
basisvalues[8] = (3.0*Z + 2.0)*basisvalues[2];
393
basisvalues[6] = (3.0*Z + 2.0)*basisvalues[1];
394
basisvalues[9] = basisvalues[3]*(0.3125 + 1.875*Z) - 0.5625*basisvalues[0];
395
basisvalues[0] *= std::sqrt(0.75);
396
basisvalues[3] *= std::sqrt(1.25);
397
basisvalues[9] *= std::sqrt(1.75);
398
basisvalues[2] *= std::sqrt(2.5);
399
basisvalues[8] *= std::sqrt(3.5);
400
basisvalues[7] *= std::sqrt(5.25);
401
basisvalues[1] *= std::sqrt(7.5);
402
basisvalues[6] *= std::sqrt(10.5);
403
basisvalues[5] *= std::sqrt(15.75);
404
basisvalues[4] *= std::sqrt(26.25);
405
406
// Table(s) of coefficients
407
static const double coefficients0[10] = \
408
{0.23094010767585, 0.121716123890037, -0.0702728368926307, 0.0993807989999907, 0.0, 0.0, 0.102868899974728, 0.0, -0.0593913870916499, -0.0671936840905293};
409
410
// Compute value(s)
411
for (unsigned int r = 0; r < 10; r++)
412
{
413
*values += coefficients0[r]*basisvalues[r];
414
}// end loop over 'r'
415
break;
416
}
417
case 6:
418
{
419
420
// Array of basisvalues
421
double basisvalues[10] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
422
423
// Declare helper variables
424
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
425
double tmp1 = 0.25*(Y + Z)*(Y + Z);
426
double tmp2 = 0.5*(1.0 + Z + 2.0*Y);
427
double tmp3 = 0.5*(1.0 - Z);
428
double tmp4 = tmp3*tmp3;
429
430
// Compute basisvalues
431
basisvalues[0] = 1.0;
432
basisvalues[1] = tmp0;
433
basisvalues[4] = 1.5*tmp0*basisvalues[1] - 0.5*tmp1*basisvalues[0];
434
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
435
basisvalues[5] = (0.5*(2.0 + 3.0*Y + Z) + 1.0*(1.0 + Y))*basisvalues[1];
436
basisvalues[7] = (1.66666666666667*tmp2 + 0.111111111111111*tmp3)*basisvalues[2] - 0.555555555555556*tmp4*basisvalues[0];
437
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
438
basisvalues[8] = (3.0*Z + 2.0)*basisvalues[2];
439
basisvalues[6] = (3.0*Z + 2.0)*basisvalues[1];
440
basisvalues[9] = basisvalues[3]*(0.3125 + 1.875*Z) - 0.5625*basisvalues[0];
441
basisvalues[0] *= std::sqrt(0.75);
442
basisvalues[3] *= std::sqrt(1.25);
443
basisvalues[9] *= std::sqrt(1.75);
444
basisvalues[2] *= std::sqrt(2.5);
445
basisvalues[8] *= std::sqrt(3.5);
446
basisvalues[7] *= std::sqrt(5.25);
447
basisvalues[1] *= std::sqrt(7.5);
448
basisvalues[6] *= std::sqrt(10.5);
449
basisvalues[5] *= std::sqrt(15.75);
450
basisvalues[4] *= std::sqrt(26.25);
451
452
// Table(s) of coefficients
453
static const double coefficients0[10] = \
454
{0.23094010767585, 0.121716123890037, 0.0702728368926307, -0.0993807989999906, 0.0, 0.100790526135794, -0.0205737799949456, -0.087287156094397, -0.01187827741833, 0.0167984210226323};
455
456
// Compute value(s)
457
for (unsigned int r = 0; r < 10; r++)
458
{
459
*values += coefficients0[r]*basisvalues[r];
460
}// end loop over 'r'
461
break;
462
}
463
case 7:
464
{
465
466
// Array of basisvalues
467
double basisvalues[10] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
468
469
// Declare helper variables
470
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
471
double tmp1 = 0.25*(Y + Z)*(Y + Z);
472
double tmp2 = 0.5*(1.0 + Z + 2.0*Y);
473
double tmp3 = 0.5*(1.0 - Z);
474
double tmp4 = tmp3*tmp3;
475
476
// Compute basisvalues
477
basisvalues[0] = 1.0;
478
basisvalues[1] = tmp0;
479
basisvalues[4] = 1.5*tmp0*basisvalues[1] - 0.5*tmp1*basisvalues[0];
480
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
481
basisvalues[5] = (0.5*(2.0 + 3.0*Y + Z) + 1.0*(1.0 + Y))*basisvalues[1];
482
basisvalues[7] = (1.66666666666667*tmp2 + 0.111111111111111*tmp3)*basisvalues[2] - 0.555555555555556*tmp4*basisvalues[0];
483
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
484
basisvalues[8] = (3.0*Z + 2.0)*basisvalues[2];
485
basisvalues[6] = (3.0*Z + 2.0)*basisvalues[1];
486
basisvalues[9] = basisvalues[3]*(0.3125 + 1.875*Z) - 0.5625*basisvalues[0];
487
basisvalues[0] *= std::sqrt(0.75);
488
basisvalues[3] *= std::sqrt(1.25);
489
basisvalues[9] *= std::sqrt(1.75);
490
basisvalues[2] *= std::sqrt(2.5);
491
basisvalues[8] *= std::sqrt(3.5);
492
basisvalues[7] *= std::sqrt(5.25);
493
basisvalues[1] *= std::sqrt(7.5);
494
basisvalues[6] *= std::sqrt(10.5);
495
basisvalues[5] *= std::sqrt(15.75);
496
basisvalues[4] *= std::sqrt(26.25);
497
498
// Table(s) of coefficients
499
static const double coefficients0[10] = \
500
{0.23094010767585, -0.121716123890037, -0.0702728368926306, 0.0993807989999906, 0.0, 0.0, -0.102868899974728, 0.0, -0.0593913870916499, -0.0671936840905293};
501
502
// Compute value(s)
503
for (unsigned int r = 0; r < 10; r++)
504
{
505
*values += coefficients0[r]*basisvalues[r];
506
}// end loop over 'r'
507
break;
508
}
509
case 8:
510
{
511
512
// Array of basisvalues
513
double basisvalues[10] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
514
515
// Declare helper variables
516
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
517
double tmp1 = 0.25*(Y + Z)*(Y + Z);
518
double tmp2 = 0.5*(1.0 + Z + 2.0*Y);
519
double tmp3 = 0.5*(1.0 - Z);
520
double tmp4 = tmp3*tmp3;
521
522
// Compute basisvalues
523
basisvalues[0] = 1.0;
524
basisvalues[1] = tmp0;
525
basisvalues[4] = 1.5*tmp0*basisvalues[1] - 0.5*tmp1*basisvalues[0];
526
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
527
basisvalues[5] = (0.5*(2.0 + 3.0*Y + Z) + 1.0*(1.0 + Y))*basisvalues[1];
528
basisvalues[7] = (1.66666666666667*tmp2 + 0.111111111111111*tmp3)*basisvalues[2] - 0.555555555555556*tmp4*basisvalues[0];
529
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
530
basisvalues[8] = (3.0*Z + 2.0)*basisvalues[2];
531
basisvalues[6] = (3.0*Z + 2.0)*basisvalues[1];
532
basisvalues[9] = basisvalues[3]*(0.3125 + 1.875*Z) - 0.5625*basisvalues[0];
533
basisvalues[0] *= std::sqrt(0.75);
534
basisvalues[3] *= std::sqrt(1.25);
535
basisvalues[9] *= std::sqrt(1.75);
536
basisvalues[2] *= std::sqrt(2.5);
537
basisvalues[8] *= std::sqrt(3.5);
538
basisvalues[7] *= std::sqrt(5.25);
539
basisvalues[1] *= std::sqrt(7.5);
540
basisvalues[6] *= std::sqrt(10.5);
541
basisvalues[5] *= std::sqrt(15.75);
542
basisvalues[4] *= std::sqrt(26.25);
543
544
// Table(s) of coefficients
545
static const double coefficients0[10] = \
546
{0.23094010767585, -0.121716123890037, 0.0702728368926306, -0.0993807989999906, 0.0, -0.100790526135794, 0.0205737799949456, -0.0872871560943969, -0.01187827741833, 0.0167984210226323};
547
548
// Compute value(s)
549
for (unsigned int r = 0; r < 10; r++)
550
{
551
*values += coefficients0[r]*basisvalues[r];
552
}// end loop over 'r'
553
break;
554
}
555
case 9:
556
{
557
558
// Array of basisvalues
559
double basisvalues[10] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
560
561
// Declare helper variables
562
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
563
double tmp1 = 0.25*(Y + Z)*(Y + Z);
564
double tmp2 = 0.5*(1.0 + Z + 2.0*Y);
565
double tmp3 = 0.5*(1.0 - Z);
566
double tmp4 = tmp3*tmp3;
567
568
// Compute basisvalues
569
basisvalues[0] = 1.0;
570
basisvalues[1] = tmp0;
571
basisvalues[4] = 1.5*tmp0*basisvalues[1] - 0.5*tmp1*basisvalues[0];
572
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
573
basisvalues[5] = (0.5*(2.0 + 3.0*Y + Z) + 1.0*(1.0 + Y))*basisvalues[1];
574
basisvalues[7] = (1.66666666666667*tmp2 + 0.111111111111111*tmp3)*basisvalues[2] - 0.555555555555556*tmp4*basisvalues[0];
575
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
576
basisvalues[8] = (3.0*Z + 2.0)*basisvalues[2];
577
basisvalues[6] = (3.0*Z + 2.0)*basisvalues[1];
578
basisvalues[9] = basisvalues[3]*(0.3125 + 1.875*Z) - 0.5625*basisvalues[0];
579
basisvalues[0] *= std::sqrt(0.75);
580
basisvalues[3] *= std::sqrt(1.25);
581
basisvalues[9] *= std::sqrt(1.75);
582
basisvalues[2] *= std::sqrt(2.5);
583
basisvalues[8] *= std::sqrt(3.5);
584
basisvalues[7] *= std::sqrt(5.25);
585
basisvalues[1] *= std::sqrt(7.5);
586
basisvalues[6] *= std::sqrt(10.5);
587
basisvalues[5] *= std::sqrt(15.75);
588
basisvalues[4] *= std::sqrt(26.25);
589
590
// Table(s) of coefficients
591
static const double coefficients0[10] = \
592
{0.23094010767585, 0.0, -0.140545673785261, -0.0993807989999906, -0.130120009726471, 0.0, 0.0, 0.0290957186981323, 0.02375655483666, 0.0167984210226323};
593
594
// Compute value(s)
595
for (unsigned int r = 0; r < 10; r++)
596
{
597
*values += coefficients0[r]*basisvalues[r];
598
}// end loop over 'r'
599
break;
600
}
601
}
602
603
}
604
605
/// Evaluate all basis functions at given point x in cell
606
virtual void evaluate_basis_all(double* values,
607
const double* x,
608
const double* vertex_coordinates,
609
int cell_orientation) const
610
{
611
// Helper variable to hold values of a single dof.
612
double dof_values = 0.0;
613
614
// Loop dofs and call evaluate_basis
615
for (unsigned int r = 0; r < 10; r++)
616
{
617
evaluate_basis(r, &dof_values, x, vertex_coordinates, cell_orientation);
618
values[r] = dof_values;
619
}// end loop over 'r'
620
}
621
622
/// Evaluate order n derivatives of basis function i at given point x in cell
623
virtual void evaluate_basis_derivatives(std::size_t i,
624
std::size_t n,
625
double* values,
626
const double* x,
627
const double* vertex_coordinates,
628
int cell_orientation) const
629
{
630
// Compute Jacobian
631
double J[9];
632
compute_jacobian_tetrahedron_3d(J, vertex_coordinates);
633
634
// Compute Jacobian inverse and determinant
635
double K[9];
636
double detJ;
637
compute_jacobian_inverse_tetrahedron_3d(K, detJ, J);
638
639
640
// Compute constants
641
const double C0 = vertex_coordinates[9] + vertex_coordinates[6] + vertex_coordinates[3] - vertex_coordinates[0];
642
const double C1 = vertex_coordinates[10] + vertex_coordinates[7] + vertex_coordinates[4] - vertex_coordinates[1];
643
const double C2 = vertex_coordinates[11] + vertex_coordinates[8] + vertex_coordinates[5] - vertex_coordinates[2];
644
645
// Compute subdeterminants
646
const double d_00 = J[4]*J[8] - J[5]*J[7];
647
const double d_01 = J[5]*J[6] - J[3]*J[8];
648
const double d_02 = J[3]*J[7] - J[4]*J[6];
649
const double d_10 = J[2]*J[7] - J[1]*J[8];
650
const double d_11 = J[0]*J[8] - J[2]*J[6];
651
const double d_12 = J[1]*J[6] - J[0]*J[7];
652
const double d_20 = J[1]*J[5] - J[2]*J[4];
653
const double d_21 = J[2]*J[3] - J[0]*J[5];
654
const double d_22 = J[0]*J[4] - J[1]*J[3];
655
656
// Get coordinates and map to the reference (FIAT) element
657
double X = (d_00*(2.0*x[0] - C0) + d_10*(2.0*x[1] - C1) + d_20*(2.0*x[2] - C2)) / detJ;
658
double Y = (d_01*(2.0*x[0] - C0) + d_11*(2.0*x[1] - C1) + d_21*(2.0*x[2] - C2)) / detJ;
659
double Z = (d_02*(2.0*x[0] - C0) + d_12*(2.0*x[1] - C1) + d_22*(2.0*x[2] - C2)) / detJ;
660
661
662
// Compute number of derivatives.
663
unsigned int num_derivatives = 1;
664
for (unsigned int r = 0; r < n; r++)
665
{
666
num_derivatives *= 3;
667
}// end loop over 'r'
668
669
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
670
unsigned int **combinations = new unsigned int *[num_derivatives];
671
for (unsigned int row = 0; row < num_derivatives; row++)
672
{
673
combinations[row] = new unsigned int [n];
674
for (unsigned int col = 0; col < n; col++)
675
combinations[row][col] = 0;
676
}
677
678
// Generate combinations of derivatives
679
for (unsigned int row = 1; row < num_derivatives; row++)
680
{
681
for (unsigned int num = 0; num < row; num++)
682
{
683
for (unsigned int col = n-1; col+1 > 0; col--)
684
{
685
if (combinations[row][col] + 1 > 2)
686
combinations[row][col] = 0;
687
else
688
{
689
combinations[row][col] += 1;
690
break;
691
}
692
}
693
}
694
}
695
696
// Compute inverse of Jacobian
697
const double Jinv[3][3] = {{K[0], K[1], K[2]}, {K[3], K[4], K[5]}, {K[6], K[7], K[8]}};
698
699
// Declare transformation matrix
700
// Declare pointer to two dimensional array and initialise
701
double **transform = new double *[num_derivatives];
702
703
for (unsigned int j = 0; j < num_derivatives; j++)
704
{
705
transform[j] = new double [num_derivatives];
706
for (unsigned int k = 0; k < num_derivatives; k++)
707
transform[j][k] = 1;
708
}
709
710
// Construct transformation matrix
711
for (unsigned int row = 0; row < num_derivatives; row++)
712
{
713
for (unsigned int col = 0; col < num_derivatives; col++)
714
{
715
for (unsigned int k = 0; k < n; k++)
716
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
717
}
718
}
719
720
// Reset values. Assuming that values is always an array.
721
for (unsigned int r = 0; r < num_derivatives; r++)
722
{
723
values[r] = 0.0;
724
}// end loop over 'r'
725
726
switch (i)
727
{
728
case 0:
729
{
730
731
// Array of basisvalues
732
double basisvalues[10] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
733
734
// Declare helper variables
735
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
736
double tmp1 = 0.25*(Y + Z)*(Y + Z);
737
double tmp2 = 0.5*(1.0 + Z + 2.0*Y);
738
double tmp3 = 0.5*(1.0 - Z);
739
double tmp4 = tmp3*tmp3;
740
741
// Compute basisvalues
742
basisvalues[0] = 1.0;
743
basisvalues[1] = tmp0;
744
basisvalues[4] = 1.5*tmp0*basisvalues[1] - 0.5*tmp1*basisvalues[0];
745
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
746
basisvalues[5] = (0.5*(2.0 + 3.0*Y + Z) + 1.0*(1.0 + Y))*basisvalues[1];
747
basisvalues[7] = (1.66666666666667*tmp2 + 0.111111111111111*tmp3)*basisvalues[2] - 0.555555555555556*tmp4*basisvalues[0];
748
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
749
basisvalues[8] = (3.0*Z + 2.0)*basisvalues[2];
750
basisvalues[6] = (3.0*Z + 2.0)*basisvalues[1];
751
basisvalues[9] = basisvalues[3]*(0.3125 + 1.875*Z) - 0.5625*basisvalues[0];
752
basisvalues[0] *= std::sqrt(0.75);
753
basisvalues[3] *= std::sqrt(1.25);
754
basisvalues[9] *= std::sqrt(1.75);
755
basisvalues[2] *= std::sqrt(2.5);
756
basisvalues[8] *= std::sqrt(3.5);
757
basisvalues[7] *= std::sqrt(5.25);
758
basisvalues[1] *= std::sqrt(7.5);
759
basisvalues[6] *= std::sqrt(10.5);
760
basisvalues[5] *= std::sqrt(15.75);
761
basisvalues[4] *= std::sqrt(26.25);
762
763
// Table(s) of coefficients
764
static const double coefficients0[10] = \
765
{-0.0577350269189625, -0.0608580619450185, -0.0351364184463153, -0.0248451997499977, 0.0650600048632355, 0.050395263067897, 0.0411475599898912, 0.0290957186981323, 0.0237565548366599, 0.0167984210226323};
766
767
// Tables of derivatives of the polynomial base (transpose).
768
static const double dmats0[10][10] = \
769
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
770
{6.32455532033676, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
771
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
772
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
773
{0.0, 11.2249721603218, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
774
{4.58257569495584, 0.0, 8.36660026534075, -1.18321595661992, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
775
{3.74165738677394, 0.0, 0.0, 8.69482604771366, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
776
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
777
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
778
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
779
780
static const double dmats1[10][10] = \
781
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
782
{3.16227766016838, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
783
{5.47722557505166, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
784
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
785
{2.95803989154981, 5.61248608016091, -1.08012344973464, -0.763762615825973, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
786
{2.29128784747792, 7.24568837309472, 4.18330013267038, -0.591607978309961, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
787
{1.87082869338697, 0.0, 0.0, 4.34741302385683, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
788
{-2.64575131106459, 0.0, 9.66091783079296, 0.683130051063973, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
789
{3.24037034920393, 0.0, 0.0, 7.52994023880668, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
790
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
791
792
static const double dmats2[10][10] = \
793
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
794
{3.16227766016838, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
795
{1.82574185835055, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
796
{5.16397779494322, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
797
{2.95803989154981, 5.61248608016091, -1.08012344973464, -0.763762615825973, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
798
{2.29128784747792, 1.44913767461894, 4.18330013267038, -0.591607978309962, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
799
{1.87082869338697, 7.09929573971954, 0.0, 4.34741302385683, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
800
{1.3228756555323, 0.0, 3.86436713231718, -0.341565025531987, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
801
{1.08012344973464, 0.0, 7.09929573971954, 2.50998007960223, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
802
{-3.81881307912987, 0.0, 0.0, 8.87411967464942, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
803
804
// Compute reference derivatives.
805
// Declare pointer to array of derivatives on FIAT element.
806
double *derivatives = new double[num_derivatives];
807
for (unsigned int r = 0; r < num_derivatives; r++)
808
{
809
derivatives[r] = 0.0;
810
}// end loop over 'r'
811
812
// Declare derivative matrix (of polynomial basis).
813
double dmats[10][10] = \
814
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
815
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
816
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
817
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
818
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
819
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
820
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
821
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
822
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
823
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
824
825
// Declare (auxiliary) derivative matrix (of polynomial basis).
826
double dmats_old[10][10] = \
827
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
828
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
829
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
830
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
831
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
832
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
833
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
834
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
835
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
836
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
837
838
// Loop possible derivatives.
839
for (unsigned int r = 0; r < num_derivatives; r++)
840
{
841
// Resetting dmats values to compute next derivative.
842
for (unsigned int t = 0; t < 10; t++)
843
{
844
for (unsigned int u = 0; u < 10; u++)
845
{
846
dmats[t][u] = 0.0;
847
if (t == u)
848
{
849
dmats[t][u] = 1.0;
850
}
851
852
}// end loop over 'u'
853
}// end loop over 't'
854
855
// Looping derivative order to generate dmats.
856
for (unsigned int s = 0; s < n; s++)
857
{
858
// Updating dmats_old with new values and resetting dmats.
859
for (unsigned int t = 0; t < 10; t++)
860
{
861
for (unsigned int u = 0; u < 10; u++)
862
{
863
dmats_old[t][u] = dmats[t][u];
864
dmats[t][u] = 0.0;
865
}// end loop over 'u'
866
}// end loop over 't'
867
868
// Update dmats using an inner product.
869
if (combinations[r][s] == 0)
870
{
871
for (unsigned int t = 0; t < 10; t++)
872
{
873
for (unsigned int u = 0; u < 10; u++)
874
{
875
for (unsigned int tu = 0; tu < 10; tu++)
876
{
877
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
878
}// end loop over 'tu'
879
}// end loop over 'u'
880
}// end loop over 't'
881
}
882
883
if (combinations[r][s] == 1)
884
{
885
for (unsigned int t = 0; t < 10; t++)
886
{
887
for (unsigned int u = 0; u < 10; u++)
888
{
889
for (unsigned int tu = 0; tu < 10; tu++)
890
{
891
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
892
}// end loop over 'tu'
893
}// end loop over 'u'
894
}// end loop over 't'
895
}
896
897
if (combinations[r][s] == 2)
898
{
899
for (unsigned int t = 0; t < 10; t++)
900
{
901
for (unsigned int u = 0; u < 10; u++)
902
{
903
for (unsigned int tu = 0; tu < 10; tu++)
904
{
905
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
906
}// end loop over 'tu'
907
}// end loop over 'u'
908
}// end loop over 't'
909
}
910
911
}// end loop over 's'
912
for (unsigned int s = 0; s < 10; s++)
913
{
914
for (unsigned int t = 0; t < 10; t++)
915
{
916
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
917
}// end loop over 't'
918
}// end loop over 's'
919
}// end loop over 'r'
920
921
// Transform derivatives back to physical element
922
for (unsigned int r = 0; r < num_derivatives; r++)
923
{
924
for (unsigned int s = 0; s < num_derivatives; s++)
925
{
926
values[r] += transform[r][s]*derivatives[s];
927
}// end loop over 's'
928
}// end loop over 'r'
929
930
// Delete pointer to array of derivatives on FIAT element
931
delete [] derivatives;
932
933
// Delete pointer to array of combinations of derivatives and transform
934
for (unsigned int r = 0; r < num_derivatives; r++)
935
{
936
delete [] combinations[r];
937
}// end loop over 'r'
938
delete [] combinations;
939
for (unsigned int r = 0; r < num_derivatives; r++)
940
{
941
delete [] transform[r];
942
}// end loop over 'r'
943
delete [] transform;
944
break;
945
}
946
case 1:
947
{
948
949
// Array of basisvalues
950
double basisvalues[10] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
951
952
// Declare helper variables
953
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
954
double tmp1 = 0.25*(Y + Z)*(Y + Z);
955
double tmp2 = 0.5*(1.0 + Z + 2.0*Y);
956
double tmp3 = 0.5*(1.0 - Z);
957
double tmp4 = tmp3*tmp3;
958
959
// Compute basisvalues
960
basisvalues[0] = 1.0;
961
basisvalues[1] = tmp0;
962
basisvalues[4] = 1.5*tmp0*basisvalues[1] - 0.5*tmp1*basisvalues[0];
963
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
964
basisvalues[5] = (0.5*(2.0 + 3.0*Y + Z) + 1.0*(1.0 + Y))*basisvalues[1];
965
basisvalues[7] = (1.66666666666667*tmp2 + 0.111111111111111*tmp3)*basisvalues[2] - 0.555555555555556*tmp4*basisvalues[0];
966
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
967
basisvalues[8] = (3.0*Z + 2.0)*basisvalues[2];
968
basisvalues[6] = (3.0*Z + 2.0)*basisvalues[1];
969
basisvalues[9] = basisvalues[3]*(0.3125 + 1.875*Z) - 0.5625*basisvalues[0];
970
basisvalues[0] *= std::sqrt(0.75);
971
basisvalues[3] *= std::sqrt(1.25);
972
basisvalues[9] *= std::sqrt(1.75);
973
basisvalues[2] *= std::sqrt(2.5);
974
basisvalues[8] *= std::sqrt(3.5);
975
basisvalues[7] *= std::sqrt(5.25);
976
basisvalues[1] *= std::sqrt(7.5);
977
basisvalues[6] *= std::sqrt(10.5);
978
basisvalues[5] *= std::sqrt(15.75);
979
basisvalues[4] *= std::sqrt(26.25);
980
981
// Table(s) of coefficients
982
static const double coefficients0[10] = \
983
{-0.0577350269189625, 0.0608580619450185, -0.0351364184463153, -0.0248451997499977, 0.0650600048632355, -0.050395263067897, -0.0411475599898912, 0.0290957186981323, 0.0237565548366599, 0.0167984210226323};
984
985
// Tables of derivatives of the polynomial base (transpose).
986
static const double dmats0[10][10] = \
987
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
988
{6.32455532033676, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
989
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
990
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
991
{0.0, 11.2249721603218, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
992
{4.58257569495584, 0.0, 8.36660026534075, -1.18321595661992, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
993
{3.74165738677394, 0.0, 0.0, 8.69482604771366, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
994
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
995
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
996
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
997
998
static const double dmats1[10][10] = \
999
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1000
{3.16227766016838, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1001
{5.47722557505166, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1002
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1003
{2.95803989154981, 5.61248608016091, -1.08012344973464, -0.763762615825973, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1004
{2.29128784747792, 7.24568837309472, 4.18330013267038, -0.591607978309961, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1005
{1.87082869338697, 0.0, 0.0, 4.34741302385683, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1006
{-2.64575131106459, 0.0, 9.66091783079296, 0.683130051063973, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1007
{3.24037034920393, 0.0, 0.0, 7.52994023880668, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1008
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
1009
1010
static const double dmats2[10][10] = \
1011
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1012
{3.16227766016838, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1013
{1.82574185835055, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1014
{5.16397779494322, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1015
{2.95803989154981, 5.61248608016091, -1.08012344973464, -0.763762615825973, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1016
{2.29128784747792, 1.44913767461894, 4.18330013267038, -0.591607978309962, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1017
{1.87082869338697, 7.09929573971954, 0.0, 4.34741302385683, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1018
{1.3228756555323, 0.0, 3.86436713231718, -0.341565025531987, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1019
{1.08012344973464, 0.0, 7.09929573971954, 2.50998007960223, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1020
{-3.81881307912987, 0.0, 0.0, 8.87411967464942, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
1021
1022
// Compute reference derivatives.
1023
// Declare pointer to array of derivatives on FIAT element.
1024
double *derivatives = new double[num_derivatives];
1025
for (unsigned int r = 0; r < num_derivatives; r++)
1026
{
1027
derivatives[r] = 0.0;
1028
}// end loop over 'r'
1029
1030
// Declare derivative matrix (of polynomial basis).
1031
double dmats[10][10] = \
1032
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1033
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1034
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1035
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1036
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1037
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
1038
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
1039
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
1040
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
1041
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
1042
1043
// Declare (auxiliary) derivative matrix (of polynomial basis).
1044
double dmats_old[10][10] = \
1045
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1046
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1047
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1048
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1049
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1050
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
1051
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
1052
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
1053
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
1054
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
1055
1056
// Loop possible derivatives.
1057
for (unsigned int r = 0; r < num_derivatives; r++)
1058
{
1059
// Resetting dmats values to compute next derivative.
1060
for (unsigned int t = 0; t < 10; t++)
1061
{
1062
for (unsigned int u = 0; u < 10; u++)
1063
{
1064
dmats[t][u] = 0.0;
1065
if (t == u)
1066
{
1067
dmats[t][u] = 1.0;
1068
}
1069
1070
}// end loop over 'u'
1071
}// end loop over 't'
1072
1073
// Looping derivative order to generate dmats.
1074
for (unsigned int s = 0; s < n; s++)
1075
{
1076
// Updating dmats_old with new values and resetting dmats.
1077
for (unsigned int t = 0; t < 10; t++)
1078
{
1079
for (unsigned int u = 0; u < 10; u++)
1080
{
1081
dmats_old[t][u] = dmats[t][u];
1082
dmats[t][u] = 0.0;
1083
}// end loop over 'u'
1084
}// end loop over 't'
1085
1086
// Update dmats using an inner product.
1087
if (combinations[r][s] == 0)
1088
{
1089
for (unsigned int t = 0; t < 10; t++)
1090
{
1091
for (unsigned int u = 0; u < 10; u++)
1092
{
1093
for (unsigned int tu = 0; tu < 10; tu++)
1094
{
1095
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1096
}// end loop over 'tu'
1097
}// end loop over 'u'
1098
}// end loop over 't'
1099
}
1100
1101
if (combinations[r][s] == 1)
1102
{
1103
for (unsigned int t = 0; t < 10; t++)
1104
{
1105
for (unsigned int u = 0; u < 10; u++)
1106
{
1107
for (unsigned int tu = 0; tu < 10; tu++)
1108
{
1109
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1110
}// end loop over 'tu'
1111
}// end loop over 'u'
1112
}// end loop over 't'
1113
}
1114
1115
if (combinations[r][s] == 2)
1116
{
1117
for (unsigned int t = 0; t < 10; t++)
1118
{
1119
for (unsigned int u = 0; u < 10; u++)
1120
{
1121
for (unsigned int tu = 0; tu < 10; tu++)
1122
{
1123
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
1124
}// end loop over 'tu'
1125
}// end loop over 'u'
1126
}// end loop over 't'
1127
}
1128
1129
}// end loop over 's'
1130
for (unsigned int s = 0; s < 10; s++)
1131
{
1132
for (unsigned int t = 0; t < 10; t++)
1133
{
1134
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1135
}// end loop over 't'
1136
}// end loop over 's'
1137
}// end loop over 'r'
1138
1139
// Transform derivatives back to physical element
1140
for (unsigned int r = 0; r < num_derivatives; r++)
1141
{
1142
for (unsigned int s = 0; s < num_derivatives; s++)
1143
{
1144
values[r] += transform[r][s]*derivatives[s];
1145
}// end loop over 's'
1146
}// end loop over 'r'
1147
1148
// Delete pointer to array of derivatives on FIAT element
1149
delete [] derivatives;
1150
1151
// Delete pointer to array of combinations of derivatives and transform
1152
for (unsigned int r = 0; r < num_derivatives; r++)
1153
{
1154
delete [] combinations[r];
1155
}// end loop over 'r'
1156
delete [] combinations;
1157
for (unsigned int r = 0; r < num_derivatives; r++)
1158
{
1159
delete [] transform[r];
1160
}// end loop over 'r'
1161
delete [] transform;
1162
break;
1163
}
1164
case 2:
1165
{
1166
1167
// Array of basisvalues
1168
double basisvalues[10] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
1169
1170
// Declare helper variables
1171
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
1172
double tmp1 = 0.25*(Y + Z)*(Y + Z);
1173
double tmp2 = 0.5*(1.0 + Z + 2.0*Y);
1174
double tmp3 = 0.5*(1.0 - Z);
1175
double tmp4 = tmp3*tmp3;
1176
1177
// Compute basisvalues
1178
basisvalues[0] = 1.0;
1179
basisvalues[1] = tmp0;
1180
basisvalues[4] = 1.5*tmp0*basisvalues[1] - 0.5*tmp1*basisvalues[0];
1181
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
1182
basisvalues[5] = (0.5*(2.0 + 3.0*Y + Z) + 1.0*(1.0 + Y))*basisvalues[1];
1183
basisvalues[7] = (1.66666666666667*tmp2 + 0.111111111111111*tmp3)*basisvalues[2] - 0.555555555555556*tmp4*basisvalues[0];
1184
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
1185
basisvalues[8] = (3.0*Z + 2.0)*basisvalues[2];
1186
basisvalues[6] = (3.0*Z + 2.0)*basisvalues[1];
1187
basisvalues[9] = basisvalues[3]*(0.3125 + 1.875*Z) - 0.5625*basisvalues[0];
1188
basisvalues[0] *= std::sqrt(0.75);
1189
basisvalues[3] *= std::sqrt(1.25);
1190
basisvalues[9] *= std::sqrt(1.75);
1191
basisvalues[2] *= std::sqrt(2.5);
1192
basisvalues[8] *= std::sqrt(3.5);
1193
basisvalues[7] *= std::sqrt(5.25);
1194
basisvalues[1] *= std::sqrt(7.5);
1195
basisvalues[6] *= std::sqrt(10.5);
1196
basisvalues[5] *= std::sqrt(15.75);
1197
basisvalues[4] *= std::sqrt(26.25);
1198
1199
// Table(s) of coefficients
1200
static const double coefficients0[10] = \
1201
{-0.0577350269189626, 0.0, 0.0702728368926307, -0.0248451997499977, 0.0, 0.0, 0.0, 0.0872871560943969, -0.0475131096733199, 0.0167984210226323};
1202
1203
// Tables of derivatives of the polynomial base (transpose).
1204
static const double dmats0[10][10] = \
1205
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1206
{6.32455532033676, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1207
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1208
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1209
{0.0, 11.2249721603218, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1210
{4.58257569495584, 0.0, 8.36660026534075, -1.18321595661992, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1211
{3.74165738677394, 0.0, 0.0, 8.69482604771366, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1212
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1213
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1214
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
1215
1216
static const double dmats1[10][10] = \
1217
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1218
{3.16227766016838, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1219
{5.47722557505166, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1220
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1221
{2.95803989154981, 5.61248608016091, -1.08012344973464, -0.763762615825973, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1222
{2.29128784747792, 7.24568837309472, 4.18330013267038, -0.591607978309961, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1223
{1.87082869338697, 0.0, 0.0, 4.34741302385683, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1224
{-2.64575131106459, 0.0, 9.66091783079296, 0.683130051063973, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1225
{3.24037034920393, 0.0, 0.0, 7.52994023880668, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1226
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
1227
1228
static const double dmats2[10][10] = \
1229
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1230
{3.16227766016838, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1231
{1.82574185835055, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1232
{5.16397779494322, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1233
{2.95803989154981, 5.61248608016091, -1.08012344973464, -0.763762615825973, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1234
{2.29128784747792, 1.44913767461894, 4.18330013267038, -0.591607978309962, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1235
{1.87082869338697, 7.09929573971954, 0.0, 4.34741302385683, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1236
{1.3228756555323, 0.0, 3.86436713231718, -0.341565025531987, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1237
{1.08012344973464, 0.0, 7.09929573971954, 2.50998007960223, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1238
{-3.81881307912987, 0.0, 0.0, 8.87411967464942, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
1239
1240
// Compute reference derivatives.
1241
// Declare pointer to array of derivatives on FIAT element.
1242
double *derivatives = new double[num_derivatives];
1243
for (unsigned int r = 0; r < num_derivatives; r++)
1244
{
1245
derivatives[r] = 0.0;
1246
}// end loop over 'r'
1247
1248
// Declare derivative matrix (of polynomial basis).
1249
double dmats[10][10] = \
1250
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1251
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1252
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1253
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1254
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1255
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
1256
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
1257
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
1258
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
1259
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
1260
1261
// Declare (auxiliary) derivative matrix (of polynomial basis).
1262
double dmats_old[10][10] = \
1263
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1264
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1265
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1266
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1267
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1268
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
1269
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
1270
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
1271
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
1272
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
1273
1274
// Loop possible derivatives.
1275
for (unsigned int r = 0; r < num_derivatives; r++)
1276
{
1277
// Resetting dmats values to compute next derivative.
1278
for (unsigned int t = 0; t < 10; t++)
1279
{
1280
for (unsigned int u = 0; u < 10; u++)
1281
{
1282
dmats[t][u] = 0.0;
1283
if (t == u)
1284
{
1285
dmats[t][u] = 1.0;
1286
}
1287
1288
}// end loop over 'u'
1289
}// end loop over 't'
1290
1291
// Looping derivative order to generate dmats.
1292
for (unsigned int s = 0; s < n; s++)
1293
{
1294
// Updating dmats_old with new values and resetting dmats.
1295
for (unsigned int t = 0; t < 10; t++)
1296
{
1297
for (unsigned int u = 0; u < 10; u++)
1298
{
1299
dmats_old[t][u] = dmats[t][u];
1300
dmats[t][u] = 0.0;
1301
}// end loop over 'u'
1302
}// end loop over 't'
1303
1304
// Update dmats using an inner product.
1305
if (combinations[r][s] == 0)
1306
{
1307
for (unsigned int t = 0; t < 10; t++)
1308
{
1309
for (unsigned int u = 0; u < 10; u++)
1310
{
1311
for (unsigned int tu = 0; tu < 10; tu++)
1312
{
1313
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1314
}// end loop over 'tu'
1315
}// end loop over 'u'
1316
}// end loop over 't'
1317
}
1318
1319
if (combinations[r][s] == 1)
1320
{
1321
for (unsigned int t = 0; t < 10; t++)
1322
{
1323
for (unsigned int u = 0; u < 10; u++)
1324
{
1325
for (unsigned int tu = 0; tu < 10; tu++)
1326
{
1327
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1328
}// end loop over 'tu'
1329
}// end loop over 'u'
1330
}// end loop over 't'
1331
}
1332
1333
if (combinations[r][s] == 2)
1334
{
1335
for (unsigned int t = 0; t < 10; t++)
1336
{
1337
for (unsigned int u = 0; u < 10; u++)
1338
{
1339
for (unsigned int tu = 0; tu < 10; tu++)
1340
{
1341
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
1342
}// end loop over 'tu'
1343
}// end loop over 'u'
1344
}// end loop over 't'
1345
}
1346
1347
}// end loop over 's'
1348
for (unsigned int s = 0; s < 10; s++)
1349
{
1350
for (unsigned int t = 0; t < 10; t++)
1351
{
1352
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1353
}// end loop over 't'
1354
}// end loop over 's'
1355
}// end loop over 'r'
1356
1357
// Transform derivatives back to physical element
1358
for (unsigned int r = 0; r < num_derivatives; r++)
1359
{
1360
for (unsigned int s = 0; s < num_derivatives; s++)
1361
{
1362
values[r] += transform[r][s]*derivatives[s];
1363
}// end loop over 's'
1364
}// end loop over 'r'
1365
1366
// Delete pointer to array of derivatives on FIAT element
1367
delete [] derivatives;
1368
1369
// Delete pointer to array of combinations of derivatives and transform
1370
for (unsigned int r = 0; r < num_derivatives; r++)
1371
{
1372
delete [] combinations[r];
1373
}// end loop over 'r'
1374
delete [] combinations;
1375
for (unsigned int r = 0; r < num_derivatives; r++)
1376
{
1377
delete [] transform[r];
1378
}// end loop over 'r'
1379
delete [] transform;
1380
break;
1381
}
1382
case 3:
1383
{
1384
1385
// Array of basisvalues
1386
double basisvalues[10] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
1387
1388
// Declare helper variables
1389
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
1390
double tmp1 = 0.25*(Y + Z)*(Y + Z);
1391
double tmp2 = 0.5*(1.0 + Z + 2.0*Y);
1392
double tmp3 = 0.5*(1.0 - Z);
1393
double tmp4 = tmp3*tmp3;
1394
1395
// Compute basisvalues
1396
basisvalues[0] = 1.0;
1397
basisvalues[1] = tmp0;
1398
basisvalues[4] = 1.5*tmp0*basisvalues[1] - 0.5*tmp1*basisvalues[0];
1399
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
1400
basisvalues[5] = (0.5*(2.0 + 3.0*Y + Z) + 1.0*(1.0 + Y))*basisvalues[1];
1401
basisvalues[7] = (1.66666666666667*tmp2 + 0.111111111111111*tmp3)*basisvalues[2] - 0.555555555555556*tmp4*basisvalues[0];
1402
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
1403
basisvalues[8] = (3.0*Z + 2.0)*basisvalues[2];
1404
basisvalues[6] = (3.0*Z + 2.0)*basisvalues[1];
1405
basisvalues[9] = basisvalues[3]*(0.3125 + 1.875*Z) - 0.5625*basisvalues[0];
1406
basisvalues[0] *= std::sqrt(0.75);
1407
basisvalues[3] *= std::sqrt(1.25);
1408
basisvalues[9] *= std::sqrt(1.75);
1409
basisvalues[2] *= std::sqrt(2.5);
1410
basisvalues[8] *= std::sqrt(3.5);
1411
basisvalues[7] *= std::sqrt(5.25);
1412
basisvalues[1] *= std::sqrt(7.5);
1413
basisvalues[6] *= std::sqrt(10.5);
1414
basisvalues[5] *= std::sqrt(15.75);
1415
basisvalues[4] *= std::sqrt(26.25);
1416
1417
// Table(s) of coefficients
1418
static const double coefficients0[10] = \
1419
{-0.0577350269189626, 0.0, 0.0, 0.074535599249993, 0.0, 0.0, 0.0, 0.0, 0.0, 0.100790526135794};
1420
1421
// Tables of derivatives of the polynomial base (transpose).
1422
static const double dmats0[10][10] = \
1423
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1424
{6.32455532033676, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1425
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1426
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1427
{0.0, 11.2249721603218, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1428
{4.58257569495584, 0.0, 8.36660026534075, -1.18321595661992, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1429
{3.74165738677394, 0.0, 0.0, 8.69482604771366, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1430
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1431
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1432
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
1433
1434
static const double dmats1[10][10] = \
1435
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1436
{3.16227766016838, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1437
{5.47722557505166, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1438
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1439
{2.95803989154981, 5.61248608016091, -1.08012344973464, -0.763762615825973, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1440
{2.29128784747792, 7.24568837309472, 4.18330013267038, -0.591607978309961, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1441
{1.87082869338697, 0.0, 0.0, 4.34741302385683, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1442
{-2.64575131106459, 0.0, 9.66091783079296, 0.683130051063973, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1443
{3.24037034920393, 0.0, 0.0, 7.52994023880668, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1444
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
1445
1446
static const double dmats2[10][10] = \
1447
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1448
{3.16227766016838, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1449
{1.82574185835055, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1450
{5.16397779494322, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1451
{2.95803989154981, 5.61248608016091, -1.08012344973464, -0.763762615825973, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1452
{2.29128784747792, 1.44913767461894, 4.18330013267038, -0.591607978309962, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1453
{1.87082869338697, 7.09929573971954, 0.0, 4.34741302385683, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1454
{1.3228756555323, 0.0, 3.86436713231718, -0.341565025531987, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1455
{1.08012344973464, 0.0, 7.09929573971954, 2.50998007960223, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1456
{-3.81881307912987, 0.0, 0.0, 8.87411967464942, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
1457
1458
// Compute reference derivatives.
1459
// Declare pointer to array of derivatives on FIAT element.
1460
double *derivatives = new double[num_derivatives];
1461
for (unsigned int r = 0; r < num_derivatives; r++)
1462
{
1463
derivatives[r] = 0.0;
1464
}// end loop over 'r'
1465
1466
// Declare derivative matrix (of polynomial basis).
1467
double dmats[10][10] = \
1468
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1469
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1470
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1471
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1472
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1473
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
1474
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
1475
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
1476
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
1477
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
1478
1479
// Declare (auxiliary) derivative matrix (of polynomial basis).
1480
double dmats_old[10][10] = \
1481
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1482
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1483
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1484
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1485
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1486
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
1487
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
1488
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
1489
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
1490
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
1491
1492
// Loop possible derivatives.
1493
for (unsigned int r = 0; r < num_derivatives; r++)
1494
{
1495
// Resetting dmats values to compute next derivative.
1496
for (unsigned int t = 0; t < 10; t++)
1497
{
1498
for (unsigned int u = 0; u < 10; u++)
1499
{
1500
dmats[t][u] = 0.0;
1501
if (t == u)
1502
{
1503
dmats[t][u] = 1.0;
1504
}
1505
1506
}// end loop over 'u'
1507
}// end loop over 't'
1508
1509
// Looping derivative order to generate dmats.
1510
for (unsigned int s = 0; s < n; s++)
1511
{
1512
// Updating dmats_old with new values and resetting dmats.
1513
for (unsigned int t = 0; t < 10; t++)
1514
{
1515
for (unsigned int u = 0; u < 10; u++)
1516
{
1517
dmats_old[t][u] = dmats[t][u];
1518
dmats[t][u] = 0.0;
1519
}// end loop over 'u'
1520
}// end loop over 't'
1521
1522
// Update dmats using an inner product.
1523
if (combinations[r][s] == 0)
1524
{
1525
for (unsigned int t = 0; t < 10; t++)
1526
{
1527
for (unsigned int u = 0; u < 10; u++)
1528
{
1529
for (unsigned int tu = 0; tu < 10; tu++)
1530
{
1531
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1532
}// end loop over 'tu'
1533
}// end loop over 'u'
1534
}// end loop over 't'
1535
}
1536
1537
if (combinations[r][s] == 1)
1538
{
1539
for (unsigned int t = 0; t < 10; t++)
1540
{
1541
for (unsigned int u = 0; u < 10; u++)
1542
{
1543
for (unsigned int tu = 0; tu < 10; tu++)
1544
{
1545
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1546
}// end loop over 'tu'
1547
}// end loop over 'u'
1548
}// end loop over 't'
1549
}
1550
1551
if (combinations[r][s] == 2)
1552
{
1553
for (unsigned int t = 0; t < 10; t++)
1554
{
1555
for (unsigned int u = 0; u < 10; u++)
1556
{
1557
for (unsigned int tu = 0; tu < 10; tu++)
1558
{
1559
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
1560
}// end loop over 'tu'
1561
}// end loop over 'u'
1562
}// end loop over 't'
1563
}
1564
1565
}// end loop over 's'
1566
for (unsigned int s = 0; s < 10; s++)
1567
{
1568
for (unsigned int t = 0; t < 10; t++)
1569
{
1570
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1571
}// end loop over 't'
1572
}// end loop over 's'
1573
}// end loop over 'r'
1574
1575
// Transform derivatives back to physical element
1576
for (unsigned int r = 0; r < num_derivatives; r++)
1577
{
1578
for (unsigned int s = 0; s < num_derivatives; s++)
1579
{
1580
values[r] += transform[r][s]*derivatives[s];
1581
}// end loop over 's'
1582
}// end loop over 'r'
1583
1584
// Delete pointer to array of derivatives on FIAT element
1585
delete [] derivatives;
1586
1587
// Delete pointer to array of combinations of derivatives and transform
1588
for (unsigned int r = 0; r < num_derivatives; r++)
1589
{
1590
delete [] combinations[r];
1591
}// end loop over 'r'
1592
delete [] combinations;
1593
for (unsigned int r = 0; r < num_derivatives; r++)
1594
{
1595
delete [] transform[r];
1596
}// end loop over 'r'
1597
delete [] transform;
1598
break;
1599
}
1600
case 4:
1601
{
1602
1603
// Array of basisvalues
1604
double basisvalues[10] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
1605
1606
// Declare helper variables
1607
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
1608
double tmp1 = 0.25*(Y + Z)*(Y + Z);
1609
double tmp2 = 0.5*(1.0 + Z + 2.0*Y);
1610
double tmp3 = 0.5*(1.0 - Z);
1611
double tmp4 = tmp3*tmp3;
1612
1613
// Compute basisvalues
1614
basisvalues[0] = 1.0;
1615
basisvalues[1] = tmp0;
1616
basisvalues[4] = 1.5*tmp0*basisvalues[1] - 0.5*tmp1*basisvalues[0];
1617
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
1618
basisvalues[5] = (0.5*(2.0 + 3.0*Y + Z) + 1.0*(1.0 + Y))*basisvalues[1];
1619
basisvalues[7] = (1.66666666666667*tmp2 + 0.111111111111111*tmp3)*basisvalues[2] - 0.555555555555556*tmp4*basisvalues[0];
1620
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
1621
basisvalues[8] = (3.0*Z + 2.0)*basisvalues[2];
1622
basisvalues[6] = (3.0*Z + 2.0)*basisvalues[1];
1623
basisvalues[9] = basisvalues[3]*(0.3125 + 1.875*Z) - 0.5625*basisvalues[0];
1624
basisvalues[0] *= std::sqrt(0.75);
1625
basisvalues[3] *= std::sqrt(1.25);
1626
basisvalues[9] *= std::sqrt(1.75);
1627
basisvalues[2] *= std::sqrt(2.5);
1628
basisvalues[8] *= std::sqrt(3.5);
1629
basisvalues[7] *= std::sqrt(5.25);
1630
basisvalues[1] *= std::sqrt(7.5);
1631
basisvalues[6] *= std::sqrt(10.5);
1632
basisvalues[5] *= std::sqrt(15.75);
1633
basisvalues[4] *= std::sqrt(26.25);
1634
1635
// Table(s) of coefficients
1636
static const double coefficients0[10] = \
1637
{0.23094010767585, 0.0, 0.140545673785261, 0.0993807989999907, 0.0, 0.0, 0.0, 0.0, 0.1187827741833, -0.0671936840905293};
1638
1639
// Tables of derivatives of the polynomial base (transpose).
1640
static const double dmats0[10][10] = \
1641
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1642
{6.32455532033676, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1643
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1644
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1645
{0.0, 11.2249721603218, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1646
{4.58257569495584, 0.0, 8.36660026534075, -1.18321595661992, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1647
{3.74165738677394, 0.0, 0.0, 8.69482604771366, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1648
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1649
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1650
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
1651
1652
static const double dmats1[10][10] = \
1653
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1654
{3.16227766016838, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1655
{5.47722557505166, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1656
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1657
{2.95803989154981, 5.61248608016091, -1.08012344973464, -0.763762615825973, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1658
{2.29128784747792, 7.24568837309472, 4.18330013267038, -0.591607978309961, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1659
{1.87082869338697, 0.0, 0.0, 4.34741302385683, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1660
{-2.64575131106459, 0.0, 9.66091783079296, 0.683130051063973, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1661
{3.24037034920393, 0.0, 0.0, 7.52994023880668, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1662
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
1663
1664
static const double dmats2[10][10] = \
1665
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1666
{3.16227766016838, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1667
{1.82574185835055, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1668
{5.16397779494322, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1669
{2.95803989154981, 5.61248608016091, -1.08012344973464, -0.763762615825973, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1670
{2.29128784747792, 1.44913767461894, 4.18330013267038, -0.591607978309962, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1671
{1.87082869338697, 7.09929573971954, 0.0, 4.34741302385683, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1672
{1.3228756555323, 0.0, 3.86436713231718, -0.341565025531987, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1673
{1.08012344973464, 0.0, 7.09929573971954, 2.50998007960223, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1674
{-3.81881307912987, 0.0, 0.0, 8.87411967464942, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
1675
1676
// Compute reference derivatives.
1677
// Declare pointer to array of derivatives on FIAT element.
1678
double *derivatives = new double[num_derivatives];
1679
for (unsigned int r = 0; r < num_derivatives; r++)
1680
{
1681
derivatives[r] = 0.0;
1682
}// end loop over 'r'
1683
1684
// Declare derivative matrix (of polynomial basis).
1685
double dmats[10][10] = \
1686
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1687
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1688
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1689
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1690
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1691
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
1692
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
1693
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
1694
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
1695
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
1696
1697
// Declare (auxiliary) derivative matrix (of polynomial basis).
1698
double dmats_old[10][10] = \
1699
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1700
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1701
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1702
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1703
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1704
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
1705
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
1706
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
1707
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
1708
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
1709
1710
// Loop possible derivatives.
1711
for (unsigned int r = 0; r < num_derivatives; r++)
1712
{
1713
// Resetting dmats values to compute next derivative.
1714
for (unsigned int t = 0; t < 10; t++)
1715
{
1716
for (unsigned int u = 0; u < 10; u++)
1717
{
1718
dmats[t][u] = 0.0;
1719
if (t == u)
1720
{
1721
dmats[t][u] = 1.0;
1722
}
1723
1724
}// end loop over 'u'
1725
}// end loop over 't'
1726
1727
// Looping derivative order to generate dmats.
1728
for (unsigned int s = 0; s < n; s++)
1729
{
1730
// Updating dmats_old with new values and resetting dmats.
1731
for (unsigned int t = 0; t < 10; t++)
1732
{
1733
for (unsigned int u = 0; u < 10; u++)
1734
{
1735
dmats_old[t][u] = dmats[t][u];
1736
dmats[t][u] = 0.0;
1737
}// end loop over 'u'
1738
}// end loop over 't'
1739
1740
// Update dmats using an inner product.
1741
if (combinations[r][s] == 0)
1742
{
1743
for (unsigned int t = 0; t < 10; t++)
1744
{
1745
for (unsigned int u = 0; u < 10; u++)
1746
{
1747
for (unsigned int tu = 0; tu < 10; tu++)
1748
{
1749
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1750
}// end loop over 'tu'
1751
}// end loop over 'u'
1752
}// end loop over 't'
1753
}
1754
1755
if (combinations[r][s] == 1)
1756
{
1757
for (unsigned int t = 0; t < 10; t++)
1758
{
1759
for (unsigned int u = 0; u < 10; u++)
1760
{
1761
for (unsigned int tu = 0; tu < 10; tu++)
1762
{
1763
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1764
}// end loop over 'tu'
1765
}// end loop over 'u'
1766
}// end loop over 't'
1767
}
1768
1769
if (combinations[r][s] == 2)
1770
{
1771
for (unsigned int t = 0; t < 10; t++)
1772
{
1773
for (unsigned int u = 0; u < 10; u++)
1774
{
1775
for (unsigned int tu = 0; tu < 10; tu++)
1776
{
1777
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
1778
}// end loop over 'tu'
1779
}// end loop over 'u'
1780
}// end loop over 't'
1781
}
1782
1783
}// end loop over 's'
1784
for (unsigned int s = 0; s < 10; s++)
1785
{
1786
for (unsigned int t = 0; t < 10; t++)
1787
{
1788
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1789
}// end loop over 't'
1790
}// end loop over 's'
1791
}// end loop over 'r'
1792
1793
// Transform derivatives back to physical element
1794
for (unsigned int r = 0; r < num_derivatives; r++)
1795
{
1796
for (unsigned int s = 0; s < num_derivatives; s++)
1797
{
1798
values[r] += transform[r][s]*derivatives[s];
1799
}// end loop over 's'
1800
}// end loop over 'r'
1801
1802
// Delete pointer to array of derivatives on FIAT element
1803
delete [] derivatives;
1804
1805
// Delete pointer to array of combinations of derivatives and transform
1806
for (unsigned int r = 0; r < num_derivatives; r++)
1807
{
1808
delete [] combinations[r];
1809
}// end loop over 'r'
1810
delete [] combinations;
1811
for (unsigned int r = 0; r < num_derivatives; r++)
1812
{
1813
delete [] transform[r];
1814
}// end loop over 'r'
1815
delete [] transform;
1816
break;
1817
}
1818
case 5:
1819
{
1820
1821
// Array of basisvalues
1822
double basisvalues[10] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
1823
1824
// Declare helper variables
1825
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
1826
double tmp1 = 0.25*(Y + Z)*(Y + Z);
1827
double tmp2 = 0.5*(1.0 + Z + 2.0*Y);
1828
double tmp3 = 0.5*(1.0 - Z);
1829
double tmp4 = tmp3*tmp3;
1830
1831
// Compute basisvalues
1832
basisvalues[0] = 1.0;
1833
basisvalues[1] = tmp0;
1834
basisvalues[4] = 1.5*tmp0*basisvalues[1] - 0.5*tmp1*basisvalues[0];
1835
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
1836
basisvalues[5] = (0.5*(2.0 + 3.0*Y + Z) + 1.0*(1.0 + Y))*basisvalues[1];
1837
basisvalues[7] = (1.66666666666667*tmp2 + 0.111111111111111*tmp3)*basisvalues[2] - 0.555555555555556*tmp4*basisvalues[0];
1838
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
1839
basisvalues[8] = (3.0*Z + 2.0)*basisvalues[2];
1840
basisvalues[6] = (3.0*Z + 2.0)*basisvalues[1];
1841
basisvalues[9] = basisvalues[3]*(0.3125 + 1.875*Z) - 0.5625*basisvalues[0];
1842
basisvalues[0] *= std::sqrt(0.75);
1843
basisvalues[3] *= std::sqrt(1.25);
1844
basisvalues[9] *= std::sqrt(1.75);
1845
basisvalues[2] *= std::sqrt(2.5);
1846
basisvalues[8] *= std::sqrt(3.5);
1847
basisvalues[7] *= std::sqrt(5.25);
1848
basisvalues[1] *= std::sqrt(7.5);
1849
basisvalues[6] *= std::sqrt(10.5);
1850
basisvalues[5] *= std::sqrt(15.75);
1851
basisvalues[4] *= std::sqrt(26.25);
1852
1853
// Table(s) of coefficients
1854
static const double coefficients0[10] = \
1855
{0.23094010767585, 0.121716123890037, -0.0702728368926307, 0.0993807989999907, 0.0, 0.0, 0.102868899974728, 0.0, -0.0593913870916499, -0.0671936840905293};
1856
1857
// Tables of derivatives of the polynomial base (transpose).
1858
static const double dmats0[10][10] = \
1859
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1860
{6.32455532033676, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1861
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1862
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1863
{0.0, 11.2249721603218, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1864
{4.58257569495584, 0.0, 8.36660026534075, -1.18321595661992, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1865
{3.74165738677394, 0.0, 0.0, 8.69482604771366, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1866
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1867
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1868
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
1869
1870
static const double dmats1[10][10] = \
1871
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1872
{3.16227766016838, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1873
{5.47722557505166, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1874
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1875
{2.95803989154981, 5.61248608016091, -1.08012344973464, -0.763762615825973, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1876
{2.29128784747792, 7.24568837309472, 4.18330013267038, -0.591607978309961, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1877
{1.87082869338697, 0.0, 0.0, 4.34741302385683, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1878
{-2.64575131106459, 0.0, 9.66091783079296, 0.683130051063973, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1879
{3.24037034920393, 0.0, 0.0, 7.52994023880668, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1880
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
1881
1882
static const double dmats2[10][10] = \
1883
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1884
{3.16227766016838, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1885
{1.82574185835055, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1886
{5.16397779494322, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1887
{2.95803989154981, 5.61248608016091, -1.08012344973464, -0.763762615825973, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1888
{2.29128784747792, 1.44913767461894, 4.18330013267038, -0.591607978309962, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1889
{1.87082869338697, 7.09929573971954, 0.0, 4.34741302385683, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1890
{1.3228756555323, 0.0, 3.86436713231718, -0.341565025531987, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1891
{1.08012344973464, 0.0, 7.09929573971954, 2.50998007960223, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1892
{-3.81881307912987, 0.0, 0.0, 8.87411967464942, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
1893
1894
// Compute reference derivatives.
1895
// Declare pointer to array of derivatives on FIAT element.
1896
double *derivatives = new double[num_derivatives];
1897
for (unsigned int r = 0; r < num_derivatives; r++)
1898
{
1899
derivatives[r] = 0.0;
1900
}// end loop over 'r'
1901
1902
// Declare derivative matrix (of polynomial basis).
1903
double dmats[10][10] = \
1904
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1905
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1906
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1907
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1908
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1909
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
1910
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
1911
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
1912
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
1913
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
1914
1915
// Declare (auxiliary) derivative matrix (of polynomial basis).
1916
double dmats_old[10][10] = \
1917
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1918
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1919
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1920
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1921
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1922
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
1923
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
1924
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
1925
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
1926
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
1927
1928
// Loop possible derivatives.
1929
for (unsigned int r = 0; r < num_derivatives; r++)
1930
{
1931
// Resetting dmats values to compute next derivative.
1932
for (unsigned int t = 0; t < 10; t++)
1933
{
1934
for (unsigned int u = 0; u < 10; u++)
1935
{
1936
dmats[t][u] = 0.0;
1937
if (t == u)
1938
{
1939
dmats[t][u] = 1.0;
1940
}
1941
1942
}// end loop over 'u'
1943
}// end loop over 't'
1944
1945
// Looping derivative order to generate dmats.
1946
for (unsigned int s = 0; s < n; s++)
1947
{
1948
// Updating dmats_old with new values and resetting dmats.
1949
for (unsigned int t = 0; t < 10; t++)
1950
{
1951
for (unsigned int u = 0; u < 10; u++)
1952
{
1953
dmats_old[t][u] = dmats[t][u];
1954
dmats[t][u] = 0.0;
1955
}// end loop over 'u'
1956
}// end loop over 't'
1957
1958
// Update dmats using an inner product.
1959
if (combinations[r][s] == 0)
1960
{
1961
for (unsigned int t = 0; t < 10; t++)
1962
{
1963
for (unsigned int u = 0; u < 10; u++)
1964
{
1965
for (unsigned int tu = 0; tu < 10; tu++)
1966
{
1967
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1968
}// end loop over 'tu'
1969
}// end loop over 'u'
1970
}// end loop over 't'
1971
}
1972
1973
if (combinations[r][s] == 1)
1974
{
1975
for (unsigned int t = 0; t < 10; t++)
1976
{
1977
for (unsigned int u = 0; u < 10; u++)
1978
{
1979
for (unsigned int tu = 0; tu < 10; tu++)
1980
{
1981
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1982
}// end loop over 'tu'
1983
}// end loop over 'u'
1984
}// end loop over 't'
1985
}
1986
1987
if (combinations[r][s] == 2)
1988
{
1989
for (unsigned int t = 0; t < 10; t++)
1990
{
1991
for (unsigned int u = 0; u < 10; u++)
1992
{
1993
for (unsigned int tu = 0; tu < 10; tu++)
1994
{
1995
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
1996
}// end loop over 'tu'
1997
}// end loop over 'u'
1998
}// end loop over 't'
1999
}
2000
2001
}// end loop over 's'
2002
for (unsigned int s = 0; s < 10; s++)
2003
{
2004
for (unsigned int t = 0; t < 10; t++)
2005
{
2006
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2007
}// end loop over 't'
2008
}// end loop over 's'
2009
}// end loop over 'r'
2010
2011
// Transform derivatives back to physical element
2012
for (unsigned int r = 0; r < num_derivatives; r++)
2013
{
2014
for (unsigned int s = 0; s < num_derivatives; s++)
2015
{
2016
values[r] += transform[r][s]*derivatives[s];
2017
}// end loop over 's'
2018
}// end loop over 'r'
2019
2020
// Delete pointer to array of derivatives on FIAT element
2021
delete [] derivatives;
2022
2023
// Delete pointer to array of combinations of derivatives and transform
2024
for (unsigned int r = 0; r < num_derivatives; r++)
2025
{
2026
delete [] combinations[r];
2027
}// end loop over 'r'
2028
delete [] combinations;
2029
for (unsigned int r = 0; r < num_derivatives; r++)
2030
{
2031
delete [] transform[r];
2032
}// end loop over 'r'
2033
delete [] transform;
2034
break;
2035
}
2036
case 6:
2037
{
2038
2039
// Array of basisvalues
2040
double basisvalues[10] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
2041
2042
// Declare helper variables
2043
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
2044
double tmp1 = 0.25*(Y + Z)*(Y + Z);
2045
double tmp2 = 0.5*(1.0 + Z + 2.0*Y);
2046
double tmp3 = 0.5*(1.0 - Z);
2047
double tmp4 = tmp3*tmp3;
2048
2049
// Compute basisvalues
2050
basisvalues[0] = 1.0;
2051
basisvalues[1] = tmp0;
2052
basisvalues[4] = 1.5*tmp0*basisvalues[1] - 0.5*tmp1*basisvalues[0];
2053
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
2054
basisvalues[5] = (0.5*(2.0 + 3.0*Y + Z) + 1.0*(1.0 + Y))*basisvalues[1];
2055
basisvalues[7] = (1.66666666666667*tmp2 + 0.111111111111111*tmp3)*basisvalues[2] - 0.555555555555556*tmp4*basisvalues[0];
2056
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
2057
basisvalues[8] = (3.0*Z + 2.0)*basisvalues[2];
2058
basisvalues[6] = (3.0*Z + 2.0)*basisvalues[1];
2059
basisvalues[9] = basisvalues[3]*(0.3125 + 1.875*Z) - 0.5625*basisvalues[0];
2060
basisvalues[0] *= std::sqrt(0.75);
2061
basisvalues[3] *= std::sqrt(1.25);
2062
basisvalues[9] *= std::sqrt(1.75);
2063
basisvalues[2] *= std::sqrt(2.5);
2064
basisvalues[8] *= std::sqrt(3.5);
2065
basisvalues[7] *= std::sqrt(5.25);
2066
basisvalues[1] *= std::sqrt(7.5);
2067
basisvalues[6] *= std::sqrt(10.5);
2068
basisvalues[5] *= std::sqrt(15.75);
2069
basisvalues[4] *= std::sqrt(26.25);
2070
2071
// Table(s) of coefficients
2072
static const double coefficients0[10] = \
2073
{0.23094010767585, 0.121716123890037, 0.0702728368926307, -0.0993807989999906, 0.0, 0.100790526135794, -0.0205737799949456, -0.087287156094397, -0.01187827741833, 0.0167984210226323};
2074
2075
// Tables of derivatives of the polynomial base (transpose).
2076
static const double dmats0[10][10] = \
2077
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2078
{6.32455532033676, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2079
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2080
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2081
{0.0, 11.2249721603218, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2082
{4.58257569495584, 0.0, 8.36660026534075, -1.18321595661992, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2083
{3.74165738677394, 0.0, 0.0, 8.69482604771366, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2084
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2085
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2086
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
2087
2088
static const double dmats1[10][10] = \
2089
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2090
{3.16227766016838, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2091
{5.47722557505166, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2092
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2093
{2.95803989154981, 5.61248608016091, -1.08012344973464, -0.763762615825973, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2094
{2.29128784747792, 7.24568837309472, 4.18330013267038, -0.591607978309961, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2095
{1.87082869338697, 0.0, 0.0, 4.34741302385683, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2096
{-2.64575131106459, 0.0, 9.66091783079296, 0.683130051063973, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2097
{3.24037034920393, 0.0, 0.0, 7.52994023880668, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2098
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
2099
2100
static const double dmats2[10][10] = \
2101
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2102
{3.16227766016838, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2103
{1.82574185835055, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2104
{5.16397779494322, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2105
{2.95803989154981, 5.61248608016091, -1.08012344973464, -0.763762615825973, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2106
{2.29128784747792, 1.44913767461894, 4.18330013267038, -0.591607978309962, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2107
{1.87082869338697, 7.09929573971954, 0.0, 4.34741302385683, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2108
{1.3228756555323, 0.0, 3.86436713231718, -0.341565025531987, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2109
{1.08012344973464, 0.0, 7.09929573971954, 2.50998007960223, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2110
{-3.81881307912987, 0.0, 0.0, 8.87411967464942, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
2111
2112
// Compute reference derivatives.
2113
// Declare pointer to array of derivatives on FIAT element.
2114
double *derivatives = new double[num_derivatives];
2115
for (unsigned int r = 0; r < num_derivatives; r++)
2116
{
2117
derivatives[r] = 0.0;
2118
}// end loop over 'r'
2119
2120
// Declare derivative matrix (of polynomial basis).
2121
double dmats[10][10] = \
2122
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2123
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2124
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2125
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2126
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2127
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
2128
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
2129
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
2130
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
2131
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
2132
2133
// Declare (auxiliary) derivative matrix (of polynomial basis).
2134
double dmats_old[10][10] = \
2135
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2136
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2137
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2138
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2139
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2140
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
2141
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
2142
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
2143
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
2144
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
2145
2146
// Loop possible derivatives.
2147
for (unsigned int r = 0; r < num_derivatives; r++)
2148
{
2149
// Resetting dmats values to compute next derivative.
2150
for (unsigned int t = 0; t < 10; t++)
2151
{
2152
for (unsigned int u = 0; u < 10; u++)
2153
{
2154
dmats[t][u] = 0.0;
2155
if (t == u)
2156
{
2157
dmats[t][u] = 1.0;
2158
}
2159
2160
}// end loop over 'u'
2161
}// end loop over 't'
2162
2163
// Looping derivative order to generate dmats.
2164
for (unsigned int s = 0; s < n; s++)
2165
{
2166
// Updating dmats_old with new values and resetting dmats.
2167
for (unsigned int t = 0; t < 10; t++)
2168
{
2169
for (unsigned int u = 0; u < 10; u++)
2170
{
2171
dmats_old[t][u] = dmats[t][u];
2172
dmats[t][u] = 0.0;
2173
}// end loop over 'u'
2174
}// end loop over 't'
2175
2176
// Update dmats using an inner product.
2177
if (combinations[r][s] == 0)
2178
{
2179
for (unsigned int t = 0; t < 10; t++)
2180
{
2181
for (unsigned int u = 0; u < 10; u++)
2182
{
2183
for (unsigned int tu = 0; tu < 10; tu++)
2184
{
2185
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2186
}// end loop over 'tu'
2187
}// end loop over 'u'
2188
}// end loop over 't'
2189
}
2190
2191
if (combinations[r][s] == 1)
2192
{
2193
for (unsigned int t = 0; t < 10; t++)
2194
{
2195
for (unsigned int u = 0; u < 10; u++)
2196
{
2197
for (unsigned int tu = 0; tu < 10; tu++)
2198
{
2199
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2200
}// end loop over 'tu'
2201
}// end loop over 'u'
2202
}// end loop over 't'
2203
}
2204
2205
if (combinations[r][s] == 2)
2206
{
2207
for (unsigned int t = 0; t < 10; t++)
2208
{
2209
for (unsigned int u = 0; u < 10; u++)
2210
{
2211
for (unsigned int tu = 0; tu < 10; tu++)
2212
{
2213
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
2214
}// end loop over 'tu'
2215
}// end loop over 'u'
2216
}// end loop over 't'
2217
}
2218
2219
}// end loop over 's'
2220
for (unsigned int s = 0; s < 10; s++)
2221
{
2222
for (unsigned int t = 0; t < 10; t++)
2223
{
2224
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2225
}// end loop over 't'
2226
}// end loop over 's'
2227
}// end loop over 'r'
2228
2229
// Transform derivatives back to physical element
2230
for (unsigned int r = 0; r < num_derivatives; r++)
2231
{
2232
for (unsigned int s = 0; s < num_derivatives; s++)
2233
{
2234
values[r] += transform[r][s]*derivatives[s];
2235
}// end loop over 's'
2236
}// end loop over 'r'
2237
2238
// Delete pointer to array of derivatives on FIAT element
2239
delete [] derivatives;
2240
2241
// Delete pointer to array of combinations of derivatives and transform
2242
for (unsigned int r = 0; r < num_derivatives; r++)
2243
{
2244
delete [] combinations[r];
2245
}// end loop over 'r'
2246
delete [] combinations;
2247
for (unsigned int r = 0; r < num_derivatives; r++)
2248
{
2249
delete [] transform[r];
2250
}// end loop over 'r'
2251
delete [] transform;
2252
break;
2253
}
2254
case 7:
2255
{
2256
2257
// Array of basisvalues
2258
double basisvalues[10] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
2259
2260
// Declare helper variables
2261
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
2262
double tmp1 = 0.25*(Y + Z)*(Y + Z);
2263
double tmp2 = 0.5*(1.0 + Z + 2.0*Y);
2264
double tmp3 = 0.5*(1.0 - Z);
2265
double tmp4 = tmp3*tmp3;
2266
2267
// Compute basisvalues
2268
basisvalues[0] = 1.0;
2269
basisvalues[1] = tmp0;
2270
basisvalues[4] = 1.5*tmp0*basisvalues[1] - 0.5*tmp1*basisvalues[0];
2271
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
2272
basisvalues[5] = (0.5*(2.0 + 3.0*Y + Z) + 1.0*(1.0 + Y))*basisvalues[1];
2273
basisvalues[7] = (1.66666666666667*tmp2 + 0.111111111111111*tmp3)*basisvalues[2] - 0.555555555555556*tmp4*basisvalues[0];
2274
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
2275
basisvalues[8] = (3.0*Z + 2.0)*basisvalues[2];
2276
basisvalues[6] = (3.0*Z + 2.0)*basisvalues[1];
2277
basisvalues[9] = basisvalues[3]*(0.3125 + 1.875*Z) - 0.5625*basisvalues[0];
2278
basisvalues[0] *= std::sqrt(0.75);
2279
basisvalues[3] *= std::sqrt(1.25);
2280
basisvalues[9] *= std::sqrt(1.75);
2281
basisvalues[2] *= std::sqrt(2.5);
2282
basisvalues[8] *= std::sqrt(3.5);
2283
basisvalues[7] *= std::sqrt(5.25);
2284
basisvalues[1] *= std::sqrt(7.5);
2285
basisvalues[6] *= std::sqrt(10.5);
2286
basisvalues[5] *= std::sqrt(15.75);
2287
basisvalues[4] *= std::sqrt(26.25);
2288
2289
// Table(s) of coefficients
2290
static const double coefficients0[10] = \
2291
{0.23094010767585, -0.121716123890037, -0.0702728368926306, 0.0993807989999906, 0.0, 0.0, -0.102868899974728, 0.0, -0.0593913870916499, -0.0671936840905293};
2292
2293
// Tables of derivatives of the polynomial base (transpose).
2294
static const double dmats0[10][10] = \
2295
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2296
{6.32455532033676, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2297
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2298
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2299
{0.0, 11.2249721603218, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2300
{4.58257569495584, 0.0, 8.36660026534075, -1.18321595661992, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2301
{3.74165738677394, 0.0, 0.0, 8.69482604771366, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2302
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2303
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2304
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
2305
2306
static const double dmats1[10][10] = \
2307
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2308
{3.16227766016838, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2309
{5.47722557505166, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2310
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2311
{2.95803989154981, 5.61248608016091, -1.08012344973464, -0.763762615825973, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2312
{2.29128784747792, 7.24568837309472, 4.18330013267038, -0.591607978309961, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2313
{1.87082869338697, 0.0, 0.0, 4.34741302385683, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2314
{-2.64575131106459, 0.0, 9.66091783079296, 0.683130051063973, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2315
{3.24037034920393, 0.0, 0.0, 7.52994023880668, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2316
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
2317
2318
static const double dmats2[10][10] = \
2319
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2320
{3.16227766016838, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2321
{1.82574185835055, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2322
{5.16397779494322, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2323
{2.95803989154981, 5.61248608016091, -1.08012344973464, -0.763762615825973, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2324
{2.29128784747792, 1.44913767461894, 4.18330013267038, -0.591607978309962, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2325
{1.87082869338697, 7.09929573971954, 0.0, 4.34741302385683, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2326
{1.3228756555323, 0.0, 3.86436713231718, -0.341565025531987, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2327
{1.08012344973464, 0.0, 7.09929573971954, 2.50998007960223, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2328
{-3.81881307912987, 0.0, 0.0, 8.87411967464942, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
2329
2330
// Compute reference derivatives.
2331
// Declare pointer to array of derivatives on FIAT element.
2332
double *derivatives = new double[num_derivatives];
2333
for (unsigned int r = 0; r < num_derivatives; r++)
2334
{
2335
derivatives[r] = 0.0;
2336
}// end loop over 'r'
2337
2338
// Declare derivative matrix (of polynomial basis).
2339
double dmats[10][10] = \
2340
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2341
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2342
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2343
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2344
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2345
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
2346
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
2347
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
2348
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
2349
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
2350
2351
// Declare (auxiliary) derivative matrix (of polynomial basis).
2352
double dmats_old[10][10] = \
2353
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2354
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2355
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2356
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2357
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2358
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
2359
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
2360
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
2361
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
2362
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
2363
2364
// Loop possible derivatives.
2365
for (unsigned int r = 0; r < num_derivatives; r++)
2366
{
2367
// Resetting dmats values to compute next derivative.
2368
for (unsigned int t = 0; t < 10; t++)
2369
{
2370
for (unsigned int u = 0; u < 10; u++)
2371
{
2372
dmats[t][u] = 0.0;
2373
if (t == u)
2374
{
2375
dmats[t][u] = 1.0;
2376
}
2377
2378
}// end loop over 'u'
2379
}// end loop over 't'
2380
2381
// Looping derivative order to generate dmats.
2382
for (unsigned int s = 0; s < n; s++)
2383
{
2384
// Updating dmats_old with new values and resetting dmats.
2385
for (unsigned int t = 0; t < 10; t++)
2386
{
2387
for (unsigned int u = 0; u < 10; u++)
2388
{
2389
dmats_old[t][u] = dmats[t][u];
2390
dmats[t][u] = 0.0;
2391
}// end loop over 'u'
2392
}// end loop over 't'
2393
2394
// Update dmats using an inner product.
2395
if (combinations[r][s] == 0)
2396
{
2397
for (unsigned int t = 0; t < 10; t++)
2398
{
2399
for (unsigned int u = 0; u < 10; u++)
2400
{
2401
for (unsigned int tu = 0; tu < 10; tu++)
2402
{
2403
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2404
}// end loop over 'tu'
2405
}// end loop over 'u'
2406
}// end loop over 't'
2407
}
2408
2409
if (combinations[r][s] == 1)
2410
{
2411
for (unsigned int t = 0; t < 10; t++)
2412
{
2413
for (unsigned int u = 0; u < 10; u++)
2414
{
2415
for (unsigned int tu = 0; tu < 10; tu++)
2416
{
2417
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2418
}// end loop over 'tu'
2419
}// end loop over 'u'
2420
}// end loop over 't'
2421
}
2422
2423
if (combinations[r][s] == 2)
2424
{
2425
for (unsigned int t = 0; t < 10; t++)
2426
{
2427
for (unsigned int u = 0; u < 10; u++)
2428
{
2429
for (unsigned int tu = 0; tu < 10; tu++)
2430
{
2431
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
2432
}// end loop over 'tu'
2433
}// end loop over 'u'
2434
}// end loop over 't'
2435
}
2436
2437
}// end loop over 's'
2438
for (unsigned int s = 0; s < 10; s++)
2439
{
2440
for (unsigned int t = 0; t < 10; t++)
2441
{
2442
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2443
}// end loop over 't'
2444
}// end loop over 's'
2445
}// end loop over 'r'
2446
2447
// Transform derivatives back to physical element
2448
for (unsigned int r = 0; r < num_derivatives; r++)
2449
{
2450
for (unsigned int s = 0; s < num_derivatives; s++)
2451
{
2452
values[r] += transform[r][s]*derivatives[s];
2453
}// end loop over 's'
2454
}// end loop over 'r'
2455
2456
// Delete pointer to array of derivatives on FIAT element
2457
delete [] derivatives;
2458
2459
// Delete pointer to array of combinations of derivatives and transform
2460
for (unsigned int r = 0; r < num_derivatives; r++)
2461
{
2462
delete [] combinations[r];
2463
}// end loop over 'r'
2464
delete [] combinations;
2465
for (unsigned int r = 0; r < num_derivatives; r++)
2466
{
2467
delete [] transform[r];
2468
}// end loop over 'r'
2469
delete [] transform;
2470
break;
2471
}
2472
case 8:
2473
{
2474
2475
// Array of basisvalues
2476
double basisvalues[10] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
2477
2478
// Declare helper variables
2479
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
2480
double tmp1 = 0.25*(Y + Z)*(Y + Z);
2481
double tmp2 = 0.5*(1.0 + Z + 2.0*Y);
2482
double tmp3 = 0.5*(1.0 - Z);
2483
double tmp4 = tmp3*tmp3;
2484
2485
// Compute basisvalues
2486
basisvalues[0] = 1.0;
2487
basisvalues[1] = tmp0;
2488
basisvalues[4] = 1.5*tmp0*basisvalues[1] - 0.5*tmp1*basisvalues[0];
2489
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
2490
basisvalues[5] = (0.5*(2.0 + 3.0*Y + Z) + 1.0*(1.0 + Y))*basisvalues[1];
2491
basisvalues[7] = (1.66666666666667*tmp2 + 0.111111111111111*tmp3)*basisvalues[2] - 0.555555555555556*tmp4*basisvalues[0];
2492
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
2493
basisvalues[8] = (3.0*Z + 2.0)*basisvalues[2];
2494
basisvalues[6] = (3.0*Z + 2.0)*basisvalues[1];
2495
basisvalues[9] = basisvalues[3]*(0.3125 + 1.875*Z) - 0.5625*basisvalues[0];
2496
basisvalues[0] *= std::sqrt(0.75);
2497
basisvalues[3] *= std::sqrt(1.25);
2498
basisvalues[9] *= std::sqrt(1.75);
2499
basisvalues[2] *= std::sqrt(2.5);
2500
basisvalues[8] *= std::sqrt(3.5);
2501
basisvalues[7] *= std::sqrt(5.25);
2502
basisvalues[1] *= std::sqrt(7.5);
2503
basisvalues[6] *= std::sqrt(10.5);
2504
basisvalues[5] *= std::sqrt(15.75);
2505
basisvalues[4] *= std::sqrt(26.25);
2506
2507
// Table(s) of coefficients
2508
static const double coefficients0[10] = \
2509
{0.23094010767585, -0.121716123890037, 0.0702728368926306, -0.0993807989999906, 0.0, -0.100790526135794, 0.0205737799949456, -0.0872871560943969, -0.01187827741833, 0.0167984210226323};
2510
2511
// Tables of derivatives of the polynomial base (transpose).
2512
static const double dmats0[10][10] = \
2513
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2514
{6.32455532033676, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2515
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2516
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2517
{0.0, 11.2249721603218, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2518
{4.58257569495584, 0.0, 8.36660026534075, -1.18321595661992, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2519
{3.74165738677394, 0.0, 0.0, 8.69482604771366, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2520
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2521
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2522
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
2523
2524
static const double dmats1[10][10] = \
2525
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2526
{3.16227766016838, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2527
{5.47722557505166, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2528
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2529
{2.95803989154981, 5.61248608016091, -1.08012344973464, -0.763762615825973, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2530
{2.29128784747792, 7.24568837309472, 4.18330013267038, -0.591607978309961, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2531
{1.87082869338697, 0.0, 0.0, 4.34741302385683, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2532
{-2.64575131106459, 0.0, 9.66091783079296, 0.683130051063973, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2533
{3.24037034920393, 0.0, 0.0, 7.52994023880668, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2534
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
2535
2536
static const double dmats2[10][10] = \
2537
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2538
{3.16227766016838, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2539
{1.82574185835055, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2540
{5.16397779494322, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2541
{2.95803989154981, 5.61248608016091, -1.08012344973464, -0.763762615825973, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2542
{2.29128784747792, 1.44913767461894, 4.18330013267038, -0.591607978309962, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2543
{1.87082869338697, 7.09929573971954, 0.0, 4.34741302385683, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2544
{1.3228756555323, 0.0, 3.86436713231718, -0.341565025531987, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2545
{1.08012344973464, 0.0, 7.09929573971954, 2.50998007960223, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2546
{-3.81881307912987, 0.0, 0.0, 8.87411967464942, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
2547
2548
// Compute reference derivatives.
2549
// Declare pointer to array of derivatives on FIAT element.
2550
double *derivatives = new double[num_derivatives];
2551
for (unsigned int r = 0; r < num_derivatives; r++)
2552
{
2553
derivatives[r] = 0.0;
2554
}// end loop over 'r'
2555
2556
// Declare derivative matrix (of polynomial basis).
2557
double dmats[10][10] = \
2558
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2559
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2560
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2561
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2562
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2563
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
2564
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
2565
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
2566
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
2567
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
2568
2569
// Declare (auxiliary) derivative matrix (of polynomial basis).
2570
double dmats_old[10][10] = \
2571
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2572
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2573
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2574
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2575
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2576
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
2577
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
2578
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
2579
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
2580
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
2581
2582
// Loop possible derivatives.
2583
for (unsigned int r = 0; r < num_derivatives; r++)
2584
{
2585
// Resetting dmats values to compute next derivative.
2586
for (unsigned int t = 0; t < 10; t++)
2587
{
2588
for (unsigned int u = 0; u < 10; u++)
2589
{
2590
dmats[t][u] = 0.0;
2591
if (t == u)
2592
{
2593
dmats[t][u] = 1.0;
2594
}
2595
2596
}// end loop over 'u'
2597
}// end loop over 't'
2598
2599
// Looping derivative order to generate dmats.
2600
for (unsigned int s = 0; s < n; s++)
2601
{
2602
// Updating dmats_old with new values and resetting dmats.
2603
for (unsigned int t = 0; t < 10; t++)
2604
{
2605
for (unsigned int u = 0; u < 10; u++)
2606
{
2607
dmats_old[t][u] = dmats[t][u];
2608
dmats[t][u] = 0.0;
2609
}// end loop over 'u'
2610
}// end loop over 't'
2611
2612
// Update dmats using an inner product.
2613
if (combinations[r][s] == 0)
2614
{
2615
for (unsigned int t = 0; t < 10; t++)
2616
{
2617
for (unsigned int u = 0; u < 10; u++)
2618
{
2619
for (unsigned int tu = 0; tu < 10; tu++)
2620
{
2621
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2622
}// end loop over 'tu'
2623
}// end loop over 'u'
2624
}// end loop over 't'
2625
}
2626
2627
if (combinations[r][s] == 1)
2628
{
2629
for (unsigned int t = 0; t < 10; t++)
2630
{
2631
for (unsigned int u = 0; u < 10; u++)
2632
{
2633
for (unsigned int tu = 0; tu < 10; tu++)
2634
{
2635
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2636
}// end loop over 'tu'
2637
}// end loop over 'u'
2638
}// end loop over 't'
2639
}
2640
2641
if (combinations[r][s] == 2)
2642
{
2643
for (unsigned int t = 0; t < 10; t++)
2644
{
2645
for (unsigned int u = 0; u < 10; u++)
2646
{
2647
for (unsigned int tu = 0; tu < 10; tu++)
2648
{
2649
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
2650
}// end loop over 'tu'
2651
}// end loop over 'u'
2652
}// end loop over 't'
2653
}
2654
2655
}// end loop over 's'
2656
for (unsigned int s = 0; s < 10; s++)
2657
{
2658
for (unsigned int t = 0; t < 10; t++)
2659
{
2660
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2661
}// end loop over 't'
2662
}// end loop over 's'
2663
}// end loop over 'r'
2664
2665
// Transform derivatives back to physical element
2666
for (unsigned int r = 0; r < num_derivatives; r++)
2667
{
2668
for (unsigned int s = 0; s < num_derivatives; s++)
2669
{
2670
values[r] += transform[r][s]*derivatives[s];
2671
}// end loop over 's'
2672
}// end loop over 'r'
2673
2674
// Delete pointer to array of derivatives on FIAT element
2675
delete [] derivatives;
2676
2677
// Delete pointer to array of combinations of derivatives and transform
2678
for (unsigned int r = 0; r < num_derivatives; r++)
2679
{
2680
delete [] combinations[r];
2681
}// end loop over 'r'
2682
delete [] combinations;
2683
for (unsigned int r = 0; r < num_derivatives; r++)
2684
{
2685
delete [] transform[r];
2686
}// end loop over 'r'
2687
delete [] transform;
2688
break;
2689
}
2690
case 9:
2691
{
2692
2693
// Array of basisvalues
2694
double basisvalues[10] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
2695
2696
// Declare helper variables
2697
double tmp0 = 0.5*(2.0 + Y + Z + 2.0*X);
2698
double tmp1 = 0.25*(Y + Z)*(Y + Z);
2699
double tmp2 = 0.5*(1.0 + Z + 2.0*Y);
2700
double tmp3 = 0.5*(1.0 - Z);
2701
double tmp4 = tmp3*tmp3;
2702
2703
// Compute basisvalues
2704
basisvalues[0] = 1.0;
2705
basisvalues[1] = tmp0;
2706
basisvalues[4] = 1.5*tmp0*basisvalues[1] - 0.5*tmp1*basisvalues[0];
2707
basisvalues[2] = 0.5*(2.0 + 3.0*Y + Z)*basisvalues[0];
2708
basisvalues[5] = (0.5*(2.0 + 3.0*Y + Z) + 1.0*(1.0 + Y))*basisvalues[1];
2709
basisvalues[7] = (1.66666666666667*tmp2 + 0.111111111111111*tmp3)*basisvalues[2] - 0.555555555555556*tmp4*basisvalues[0];
2710
basisvalues[3] = (2.0*Z + 1.0)*basisvalues[0];
2711
basisvalues[8] = (3.0*Z + 2.0)*basisvalues[2];
2712
basisvalues[6] = (3.0*Z + 2.0)*basisvalues[1];
2713
basisvalues[9] = basisvalues[3]*(0.3125 + 1.875*Z) - 0.5625*basisvalues[0];
2714
basisvalues[0] *= std::sqrt(0.75);
2715
basisvalues[3] *= std::sqrt(1.25);
2716
basisvalues[9] *= std::sqrt(1.75);
2717
basisvalues[2] *= std::sqrt(2.5);
2718
basisvalues[8] *= std::sqrt(3.5);
2719
basisvalues[7] *= std::sqrt(5.25);
2720
basisvalues[1] *= std::sqrt(7.5);
2721
basisvalues[6] *= std::sqrt(10.5);
2722
basisvalues[5] *= std::sqrt(15.75);
2723
basisvalues[4] *= std::sqrt(26.25);
2724
2725
// Table(s) of coefficients
2726
static const double coefficients0[10] = \
2727
{0.23094010767585, 0.0, -0.140545673785261, -0.0993807989999906, -0.130120009726471, 0.0, 0.0, 0.0290957186981323, 0.02375655483666, 0.0167984210226323};
2728
2729
// Tables of derivatives of the polynomial base (transpose).
2730
static const double dmats0[10][10] = \
2731
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2732
{6.32455532033676, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2733
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2734
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2735
{0.0, 11.2249721603218, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2736
{4.58257569495584, 0.0, 8.36660026534075, -1.18321595661992, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2737
{3.74165738677394, 0.0, 0.0, 8.69482604771366, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2738
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2739
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2740
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
2741
2742
static const double dmats1[10][10] = \
2743
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2744
{3.16227766016838, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2745
{5.47722557505166, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2746
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2747
{2.95803989154981, 5.61248608016091, -1.08012344973464, -0.763762615825973, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2748
{2.29128784747792, 7.24568837309472, 4.18330013267038, -0.591607978309961, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2749
{1.87082869338697, 0.0, 0.0, 4.34741302385683, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2750
{-2.64575131106459, 0.0, 9.66091783079296, 0.683130051063973, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2751
{3.24037034920393, 0.0, 0.0, 7.52994023880668, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2752
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
2753
2754
static const double dmats2[10][10] = \
2755
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2756
{3.16227766016838, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2757
{1.82574185835055, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2758
{5.16397779494322, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2759
{2.95803989154981, 5.61248608016091, -1.08012344973464, -0.763762615825973, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2760
{2.29128784747792, 1.44913767461894, 4.18330013267038, -0.591607978309962, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2761
{1.87082869338697, 7.09929573971954, 0.0, 4.34741302385683, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2762
{1.3228756555323, 0.0, 3.86436713231718, -0.341565025531987, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2763
{1.08012344973464, 0.0, 7.09929573971954, 2.50998007960223, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2764
{-3.81881307912987, 0.0, 0.0, 8.87411967464942, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
2765
2766
// Compute reference derivatives.
2767
// Declare pointer to array of derivatives on FIAT element.
2768
double *derivatives = new double[num_derivatives];
2769
for (unsigned int r = 0; r < num_derivatives; r++)
2770
{
2771
derivatives[r] = 0.0;
2772
}// end loop over 'r'
2773
2774
// Declare derivative matrix (of polynomial basis).
2775
double dmats[10][10] = \
2776
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2777
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2778
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2779
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2780
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2781
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
2782
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
2783
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
2784
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
2785
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
2786
2787
// Declare (auxiliary) derivative matrix (of polynomial basis).
2788
double dmats_old[10][10] = \
2789
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2790
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2791
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2792
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2793
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2794
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
2795
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
2796
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
2797
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
2798
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
2799
2800
// Loop possible derivatives.
2801
for (unsigned int r = 0; r < num_derivatives; r++)
2802
{
2803
// Resetting dmats values to compute next derivative.
2804
for (unsigned int t = 0; t < 10; t++)
2805
{
2806
for (unsigned int u = 0; u < 10; u++)
2807
{
2808
dmats[t][u] = 0.0;
2809
if (t == u)
2810
{
2811
dmats[t][u] = 1.0;
2812
}
2813
2814
}// end loop over 'u'
2815
}// end loop over 't'
2816
2817
// Looping derivative order to generate dmats.
2818
for (unsigned int s = 0; s < n; s++)
2819
{
2820
// Updating dmats_old with new values and resetting dmats.
2821
for (unsigned int t = 0; t < 10; t++)
2822
{
2823
for (unsigned int u = 0; u < 10; u++)
2824
{
2825
dmats_old[t][u] = dmats[t][u];
2826
dmats[t][u] = 0.0;
2827
}// end loop over 'u'
2828
}// end loop over 't'
2829
2830
// Update dmats using an inner product.
2831
if (combinations[r][s] == 0)
2832
{
2833
for (unsigned int t = 0; t < 10; t++)
2834
{
2835
for (unsigned int u = 0; u < 10; u++)
2836
{
2837
for (unsigned int tu = 0; tu < 10; tu++)
2838
{
2839
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2840
}// end loop over 'tu'
2841
}// end loop over 'u'
2842
}// end loop over 't'
2843
}
2844
2845
if (combinations[r][s] == 1)
2846
{
2847
for (unsigned int t = 0; t < 10; t++)
2848
{
2849
for (unsigned int u = 0; u < 10; u++)
2850
{
2851
for (unsigned int tu = 0; tu < 10; tu++)
2852
{
2853
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2854
}// end loop over 'tu'
2855
}// end loop over 'u'
2856
}// end loop over 't'
2857
}
2858
2859
if (combinations[r][s] == 2)
2860
{
2861
for (unsigned int t = 0; t < 10; t++)
2862
{
2863
for (unsigned int u = 0; u < 10; u++)
2864
{
2865
for (unsigned int tu = 0; tu < 10; tu++)
2866
{
2867
dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
2868
}// end loop over 'tu'
2869
}// end loop over 'u'
2870
}// end loop over 't'
2871
}
2872
2873
}// end loop over 's'
2874
for (unsigned int s = 0; s < 10; s++)
2875
{
2876
for (unsigned int t = 0; t < 10; t++)
2877
{
2878
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2879
}// end loop over 't'
2880
}// end loop over 's'
2881
}// end loop over 'r'
2882
2883
// Transform derivatives back to physical element
2884
for (unsigned int r = 0; r < num_derivatives; r++)
2885
{
2886
for (unsigned int s = 0; s < num_derivatives; s++)
2887
{
2888
values[r] += transform[r][s]*derivatives[s];
2889
}// end loop over 's'
2890
}// end loop over 'r'
2891
2892
// Delete pointer to array of derivatives on FIAT element
2893
delete [] derivatives;
2894
2895
// Delete pointer to array of combinations of derivatives and transform
2896
for (unsigned int r = 0; r < num_derivatives; r++)
2897
{
2898
delete [] combinations[r];
2899
}// end loop over 'r'
2900
delete [] combinations;
2901
for (unsigned int r = 0; r < num_derivatives; r++)
2902
{
2903
delete [] transform[r];
2904
}// end loop over 'r'
2905
delete [] transform;
2906
break;
2907
}
2908
}
2909
2910
}
2911
2912
/// Evaluate order n derivatives of all basis functions at given point x in cell
2913
virtual void evaluate_basis_derivatives_all(std::size_t n,
2914
double* values,
2915
const double* x,
2916
const double* vertex_coordinates,
2917
int cell_orientation) const
2918
{
2919
// Compute number of derivatives.
2920
unsigned int num_derivatives = 1;
2921
for (unsigned int r = 0; r < n; r++)
2922
{
2923
num_derivatives *= 3;
2924
}// end loop over 'r'
2925
2926
// Helper variable to hold values of a single dof.
2927
double *dof_values = new double[num_derivatives];
2928
for (unsigned int r = 0; r < num_derivatives; r++)
2929
{
2930
dof_values[r] = 0.0;
2931
}// end loop over 'r'
2932
2933
// Loop dofs and call evaluate_basis_derivatives.
2934
for (unsigned int r = 0; r < 10; r++)
2935
{
2936
evaluate_basis_derivatives(r, n, dof_values, x, vertex_coordinates, cell_orientation);
2937
for (unsigned int s = 0; s < num_derivatives; s++)
2938
{
2939
values[r*num_derivatives + s] = dof_values[s];
2940
}// end loop over 's'
2941
}// end loop over 'r'
2942
2943
// Delete pointer.
2944
delete [] dof_values;
2945
}
2946
2947
/// Evaluate linear functional for dof i on the function f
2948
virtual double evaluate_dof(std::size_t i,
2949
const ufc::function& f,
2950
const double* vertex_coordinates,
2951
int cell_orientation,
2952
const ufc::cell& c) const
2953
{
2954
// Declare variables for result of evaluation
2955
double vals[1];
2956
2957
// Declare variable for physical coordinates
2958
double y[3];
2959
switch (i)
2960
{
2961
case 0:
2962
{
2963
y[0] = vertex_coordinates[0];
2964
y[1] = vertex_coordinates[1];
2965
y[2] = vertex_coordinates[2];
2966
f.evaluate(vals, y, c);
2967
return vals[0];
2968
break;
2969
}
2970
case 1:
2971
{
2972
y[0] = vertex_coordinates[3];
2973
y[1] = vertex_coordinates[4];
2974
y[2] = vertex_coordinates[5];
2975
f.evaluate(vals, y, c);
2976
return vals[0];
2977
break;
2978
}
2979
case 2:
2980
{
2981
y[0] = vertex_coordinates[6];
2982
y[1] = vertex_coordinates[7];
2983
y[2] = vertex_coordinates[8];
2984
f.evaluate(vals, y, c);
2985
return vals[0];
2986
break;
2987
}
2988
case 3:
2989
{
2990
y[0] = vertex_coordinates[9];
2991
y[1] = vertex_coordinates[10];
2992
y[2] = vertex_coordinates[11];
2993
f.evaluate(vals, y, c);
2994
return vals[0];
2995
break;
2996
}
2997
case 4:
2998
{
2999
y[0] = 0.5*vertex_coordinates[6] + 0.5*vertex_coordinates[9];
3000
y[1] = 0.5*vertex_coordinates[7] + 0.5*vertex_coordinates[10];
3001
y[2] = 0.5*vertex_coordinates[8] + 0.5*vertex_coordinates[11];
3002
f.evaluate(vals, y, c);
3003
return vals[0];
3004
break;
3005
}
3006
case 5:
3007
{
3008
y[0] = 0.5*vertex_coordinates[3] + 0.5*vertex_coordinates[9];
3009
y[1] = 0.5*vertex_coordinates[4] + 0.5*vertex_coordinates[10];
3010
y[2] = 0.5*vertex_coordinates[5] + 0.5*vertex_coordinates[11];
3011
f.evaluate(vals, y, c);
3012
return vals[0];
3013
break;
3014
}
3015
case 6:
3016
{
3017
y[0] = 0.5*vertex_coordinates[3] + 0.5*vertex_coordinates[6];
3018
y[1] = 0.5*vertex_coordinates[4] + 0.5*vertex_coordinates[7];
3019
y[2] = 0.5*vertex_coordinates[5] + 0.5*vertex_coordinates[8];
3020
f.evaluate(vals, y, c);
3021
return vals[0];
3022
break;
3023
}
3024
case 7:
3025
{
3026
y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[9];
3027
y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[10];
3028
y[2] = 0.5*vertex_coordinates[2] + 0.5*vertex_coordinates[11];
3029
f.evaluate(vals, y, c);
3030
return vals[0];
3031
break;
3032
}
3033
case 8:
3034
{
3035
y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[6];
3036
y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[7];
3037
y[2] = 0.5*vertex_coordinates[2] + 0.5*vertex_coordinates[8];
3038
f.evaluate(vals, y, c);
3039
return vals[0];
3040
break;
3041
}
3042
case 9:
3043
{
3044
y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[3];
3045
y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[4];
3046
y[2] = 0.5*vertex_coordinates[2] + 0.5*vertex_coordinates[5];
3047
f.evaluate(vals, y, c);
3048
return vals[0];
3049
break;
3050
}
3051
}
3052
3053
return 0.0;
3054
}
3055
3056
/// Evaluate linear functionals for all dofs on the function f
3057
virtual void evaluate_dofs(double* values,
3058
const ufc::function& f,
3059
const double* vertex_coordinates,
3060
int cell_orientation,
3061
const ufc::cell& c) const
3062
{
3063
// Declare variables for result of evaluation
3064
double vals[1];
3065
3066
// Declare variable for physical coordinates
3067
double y[3];
3068
y[0] = vertex_coordinates[0];
3069
y[1] = vertex_coordinates[1];
3070
y[2] = vertex_coordinates[2];
3071
f.evaluate(vals, y, c);
3072
values[0] = vals[0];
3073
y[0] = vertex_coordinates[3];
3074
y[1] = vertex_coordinates[4];
3075
y[2] = vertex_coordinates[5];
3076
f.evaluate(vals, y, c);
3077
values[1] = vals[0];
3078
y[0] = vertex_coordinates[6];
3079
y[1] = vertex_coordinates[7];
3080
y[2] = vertex_coordinates[8];
3081
f.evaluate(vals, y, c);
3082
values[2] = vals[0];
3083
y[0] = vertex_coordinates[9];
3084
y[1] = vertex_coordinates[10];
3085
y[2] = vertex_coordinates[11];
3086
f.evaluate(vals, y, c);
3087
values[3] = vals[0];
3088
y[0] = 0.5*vertex_coordinates[6] + 0.5*vertex_coordinates[9];
3089
y[1] = 0.5*vertex_coordinates[7] + 0.5*vertex_coordinates[10];
3090
y[2] = 0.5*vertex_coordinates[8] + 0.5*vertex_coordinates[11];
3091
f.evaluate(vals, y, c);
3092
values[4] = vals[0];
3093
y[0] = 0.5*vertex_coordinates[3] + 0.5*vertex_coordinates[9];
3094
y[1] = 0.5*vertex_coordinates[4] + 0.5*vertex_coordinates[10];
3095
y[2] = 0.5*vertex_coordinates[5] + 0.5*vertex_coordinates[11];
3096
f.evaluate(vals, y, c);
3097
values[5] = vals[0];
3098
y[0] = 0.5*vertex_coordinates[3] + 0.5*vertex_coordinates[6];
3099
y[1] = 0.5*vertex_coordinates[4] + 0.5*vertex_coordinates[7];
3100
y[2] = 0.5*vertex_coordinates[5] + 0.5*vertex_coordinates[8];
3101
f.evaluate(vals, y, c);
3102
values[6] = vals[0];
3103
y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[9];
3104
y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[10];
3105
y[2] = 0.5*vertex_coordinates[2] + 0.5*vertex_coordinates[11];
3106
f.evaluate(vals, y, c);
3107
values[7] = vals[0];
3108
y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[6];
3109
y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[7];
3110
y[2] = 0.5*vertex_coordinates[2] + 0.5*vertex_coordinates[8];
3111
f.evaluate(vals, y, c);
3112
values[8] = vals[0];
3113
y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[3];
3114
y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[4];
3115
y[2] = 0.5*vertex_coordinates[2] + 0.5*vertex_coordinates[5];
3116
f.evaluate(vals, y, c);
3117
values[9] = vals[0];
3118
}
3119
3120
/// Interpolate vertex values from dof values
3121
virtual void interpolate_vertex_values(double* vertex_values,
3122
const double* dof_values,
3123
const double* vertex_coordinates,
3124
int cell_orientation,
3125
const ufc::cell& c) const
3126
{
3127
// Evaluate function and change variables
3128
vertex_values[0] = dof_values[0];
3129
vertex_values[1] = dof_values[1];
3130
vertex_values[2] = dof_values[2];
3131
vertex_values[3] = dof_values[3];
3132
}
3133
3134
/// Map coordinate xhat from reference cell to coordinate x in cell
3135
virtual void map_from_reference_cell(double* x,
3136
const double* xhat,
3137
const ufc::cell& c) const
3138
{
3139
throw std::runtime_error("map_from_reference_cell not yet implemented.");
3140
}
3141
3142
/// Map from coordinate x in cell to coordinate xhat in reference cell
3143
virtual void map_to_reference_cell(double* xhat,
3144
const double* x,
3145
const ufc::cell& c) const
3146
{
3147
throw std::runtime_error("map_to_reference_cell not yet implemented.");
3148
}
3149
3150
/// Return the number of sub elements (for a mixed element)
3151
virtual std::size_t num_sub_elements() const
3152
{
3153
return 0;
3154
}
3155
3156
/// Create a new finite element for sub element i (for a mixed element)
3157
virtual ufc::finite_element* create_sub_element(std::size_t i) const
3158
{
3159
return 0;
3160
}
3161
3162
/// Create a new class instance
3163
virtual ufc::finite_element* create() const
3164
{
3165
return new projection_finite_element_0();
3166
}
3167
3168
};
3169
3170
/// This class defines the interface for a local-to-global mapping of
3171
/// degrees of freedom (dofs).
3172
3173
class projection_dofmap_0: public ufc::dofmap
3174
{
3175
public:
3176
3177
/// Constructor
3178
projection_dofmap_0() : ufc::dofmap()
3179
{
3180
// Do nothing
3181
}
3182
3183
/// Destructor
3184
virtual ~projection_dofmap_0()
3185
{
3186
// Do nothing
3187
}
3188
3189
/// Return a string identifying the dofmap
3190
virtual const char* signature() const
3191
{
3192
return "FFC dofmap for FiniteElement('Lagrange', Domain(Cell('tetrahedron', 3), 'tetrahedron_multiverse', 3, 3), 2, None)";
3193
}
3194
3195
/// Return true iff mesh entities of topological dimension d are needed
3196
virtual bool needs_mesh_entities(std::size_t d) const
3197
{
3198
switch (d)
3199
{
3200
case 0:
3201
{
3202
return true;
3203
break;
3204
}
3205
case 1:
3206
{
3207
return true;
3208
break;
3209
}
3210
case 2:
3211
{
3212
return false;
3213
break;
3214
}
3215
case 3:
3216
{
3217
return false;
3218
break;
3219
}
3220
}
3221
3222
return false;
3223
}
3224
3225
/// Return the topological dimension of the associated cell shape
3226
virtual std::size_t topological_dimension() const
3227
{
3228
return 3;
3229
}
3230
3231
/// Return the geometric dimension of the associated cell shape
3232
virtual std::size_t geometric_dimension() const
3233
{
3234
return 3;
3235
}
3236
3237
/// Return the dimension of the global finite element function space
3238
virtual std::size_t global_dimension(const std::vector<std::size_t>&
3239
num_global_entities) const
3240
{
3241
return num_global_entities[0] + num_global_entities[1];
3242
}
3243
3244
/// Return the dimension of the local finite element function space for a cell
3245
virtual std::size_t local_dimension() const
3246
{
3247
return 10;
3248
}
3249
3250
/// Return the number of dofs on each cell facet
3251
virtual std::size_t num_facet_dofs() const
3252
{
3253
return 6;
3254
}
3255
3256
/// Return the number of dofs associated with each cell entity of dimension d
3257
virtual std::size_t num_entity_dofs(std::size_t d) const
3258
{
3259
switch (d)
3260
{
3261
case 0:
3262
{
3263
return 1;
3264
break;
3265
}
3266
case 1:
3267
{
3268
return 1;
3269
break;
3270
}
3271
case 2:
3272
{
3273
return 0;
3274
break;
3275
}
3276
case 3:
3277
{
3278
return 0;
3279
break;
3280
}
3281
}
3282
3283
return 0;
3284
}
3285
3286
/// Tabulate the local-to-global mapping of dofs on a cell
3287
virtual void tabulate_dofs(std::size_t* dofs,
3288
const std::vector<std::size_t>& num_global_entities,
3289
const ufc::cell& c) const
3290
{
3291
unsigned int offset = 0;
3292
dofs[0] = offset + c.entity_indices[0][0];
3293
dofs[1] = offset + c.entity_indices[0][1];
3294
dofs[2] = offset + c.entity_indices[0][2];
3295
dofs[3] = offset + c.entity_indices[0][3];
3296
offset += num_global_entities[0];
3297
dofs[4] = offset + c.entity_indices[1][0];
3298
dofs[5] = offset + c.entity_indices[1][1];
3299
dofs[6] = offset + c.entity_indices[1][2];
3300
dofs[7] = offset + c.entity_indices[1][3];
3301
dofs[8] = offset + c.entity_indices[1][4];
3302
dofs[9] = offset + c.entity_indices[1][5];
3303
offset += num_global_entities[1];
3304
}
3305
3306
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
3307
virtual void tabulate_facet_dofs(std::size_t* dofs,
3308
std::size_t facet) const
3309
{
3310
switch (facet)
3311
{
3312
case 0:
3313
{
3314
dofs[0] = 1;
3315
dofs[1] = 2;
3316
dofs[2] = 3;
3317
dofs[3] = 4;
3318
dofs[4] = 5;
3319
dofs[5] = 6;
3320
break;
3321
}
3322
case 1:
3323
{
3324
dofs[0] = 0;
3325
dofs[1] = 2;
3326
dofs[2] = 3;
3327
dofs[3] = 4;
3328
dofs[4] = 7;
3329
dofs[5] = 8;
3330
break;
3331
}
3332
case 2:
3333
{
3334
dofs[0] = 0;
3335
dofs[1] = 1;
3336
dofs[2] = 3;
3337
dofs[3] = 5;
3338
dofs[4] = 7;
3339
dofs[5] = 9;
3340
break;
3341
}
3342
case 3:
3343
{
3344
dofs[0] = 0;
3345
dofs[1] = 1;
3346
dofs[2] = 2;
3347
dofs[3] = 6;
3348
dofs[4] = 8;
3349
dofs[5] = 9;
3350
break;
3351
}
3352
}
3353
3354
}
3355
3356
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
3357
virtual void tabulate_entity_dofs(std::size_t* dofs,
3358
std::size_t d, std::size_t i) const
3359
{
3360
if (d > 3)
3361
{
3362
throw std::runtime_error("d is larger than dimension (3)");
3363
}
3364
3365
switch (d)
3366
{
3367
case 0:
3368
{
3369
if (i > 3)
3370
{
3371
throw std::runtime_error("i is larger than number of entities (3)");
3372
}
3373
3374
switch (i)
3375
{
3376
case 0:
3377
{
3378
dofs[0] = 0;
3379
break;
3380
}
3381
case 1:
3382
{
3383
dofs[0] = 1;
3384
break;
3385
}
3386
case 2:
3387
{
3388
dofs[0] = 2;
3389
break;
3390
}
3391
case 3:
3392
{
3393
dofs[0] = 3;
3394
break;
3395
}
3396
}
3397
3398
break;
3399
}
3400
case 1:
3401
{
3402
if (i > 5)
3403
{
3404
throw std::runtime_error("i is larger than number of entities (5)");
3405
}
3406
3407
switch (i)
3408
{
3409
case 0:
3410
{
3411
dofs[0] = 4;
3412
break;
3413
}
3414
case 1:
3415
{
3416
dofs[0] = 5;
3417
break;
3418
}
3419
case 2:
3420
{
3421
dofs[0] = 6;
3422
break;
3423
}
3424
case 3:
3425
{
3426
dofs[0] = 7;
3427
break;
3428
}
3429
case 4:
3430
{
3431
dofs[0] = 8;
3432
break;
3433
}
3434
case 5:
3435
{
3436
dofs[0] = 9;
3437
break;
3438
}
3439
}
3440
3441
break;
3442
}
3443
case 2:
3444
{
3445
3446
break;
3447
}
3448
case 3:
3449
{
3450
3451
break;
3452
}
3453
}
3454
3455
}
3456
3457
/// Tabulate the coordinates of all dofs on a cell
3458
virtual void tabulate_coordinates(double** dof_coordinates,
3459
const double* vertex_coordinates) const
3460
{
3461
dof_coordinates[0][0] = vertex_coordinates[0];
3462
dof_coordinates[0][1] = vertex_coordinates[1];
3463
dof_coordinates[0][2] = vertex_coordinates[2];
3464
dof_coordinates[1][0] = vertex_coordinates[3];
3465
dof_coordinates[1][1] = vertex_coordinates[4];
3466
dof_coordinates[1][2] = vertex_coordinates[5];
3467
dof_coordinates[2][0] = vertex_coordinates[6];
3468
dof_coordinates[2][1] = vertex_coordinates[7];
3469
dof_coordinates[2][2] = vertex_coordinates[8];
3470
dof_coordinates[3][0] = vertex_coordinates[9];
3471
dof_coordinates[3][1] = vertex_coordinates[10];
3472
dof_coordinates[3][2] = vertex_coordinates[11];
3473
dof_coordinates[4][0] = 0.5*vertex_coordinates[6] + 0.5*vertex_coordinates[9];
3474
dof_coordinates[4][1] = 0.5*vertex_coordinates[7] + 0.5*vertex_coordinates[10];
3475
dof_coordinates[4][2] = 0.5*vertex_coordinates[8] + 0.5*vertex_coordinates[11];
3476
dof_coordinates[5][0] = 0.5*vertex_coordinates[3] + 0.5*vertex_coordinates[9];
3477
dof_coordinates[5][1] = 0.5*vertex_coordinates[4] + 0.5*vertex_coordinates[10];
3478
dof_coordinates[5][2] = 0.5*vertex_coordinates[5] + 0.5*vertex_coordinates[11];
3479
dof_coordinates[6][0] = 0.5*vertex_coordinates[3] + 0.5*vertex_coordinates[6];
3480
dof_coordinates[6][1] = 0.5*vertex_coordinates[4] + 0.5*vertex_coordinates[7];
3481
dof_coordinates[6][2] = 0.5*vertex_coordinates[5] + 0.5*vertex_coordinates[8];
3482
dof_coordinates[7][0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[9];
3483
dof_coordinates[7][1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[10];
3484
dof_coordinates[7][2] = 0.5*vertex_coordinates[2] + 0.5*vertex_coordinates[11];
3485
dof_coordinates[8][0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[6];
3486
dof_coordinates[8][1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[7];
3487
dof_coordinates[8][2] = 0.5*vertex_coordinates[2] + 0.5*vertex_coordinates[8];
3488
dof_coordinates[9][0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[3];
3489
dof_coordinates[9][1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[4];
3490
dof_coordinates[9][2] = 0.5*vertex_coordinates[2] + 0.5*vertex_coordinates[5];
3491
}
3492
3493
/// Return the number of sub dofmaps (for a mixed element)
3494
virtual std::size_t num_sub_dofmaps() const
3495
{
3496
return 0;
3497
}
3498
3499
/// Create a new dofmap for sub dofmap i (for a mixed element)
3500
virtual ufc::dofmap* create_sub_dofmap(std::size_t i) const
3501
{
3502
return 0;
3503
}
3504
3505
/// Create a new class instance
3506
virtual ufc::dofmap* create() const
3507
{
3508
return new projection_dofmap_0();
3509
}
3510
3511
};
3512
3513
/// This class defines the interface for the tabulation of the cell
3514
/// tensor corresponding to the local contribution to a form from
3515
/// the integral over a cell.
3516
3517
class projection_cell_integral_0_otherwise: public ufc::cell_integral
3518
{
3519
public:
3520
3521
/// Constructor
3522
projection_cell_integral_0_otherwise() : ufc::cell_integral()
3523
{
3524
// Do nothing
3525
}
3526
3527
/// Destructor
3528
virtual ~projection_cell_integral_0_otherwise()
3529
{
3530
// Do nothing
3531
}
3532
3533
/// Tabulate the tensor for the contribution from a local cell
3534
virtual void tabulate_tensor(double* A,
3535
const double * const * w,
3536
const double* vertex_coordinates,
3537
int cell_orientation) const
3538
{
3539
// Number of operations (multiply-add pairs) for Jacobian data: 3
3540
// Number of operations (multiply-add pairs) for geometry tensor: 0
3541
// Number of operations (multiply-add pairs) for tensor contraction: 50
3542
// Total number of operations (multiply-add pairs): 53
3543
3544
// Compute Jacobian
3545
double J[9];
3546
compute_jacobian_tetrahedron_3d(J, vertex_coordinates);
3547
3548
// Compute Jacobian inverse and determinant
3549
double K[9];
3550
double detJ;
3551
compute_jacobian_inverse_tetrahedron_3d(K, detJ, J);
3552
3553
// Set scale factor
3554
const double det = std::abs(detJ);
3555
3556
// Compute geometry tensor
3557
const double G0_ = det;
3558
3559
// Compute element tensor
3560
A[0] = 0.00238095238095238*G0_;
3561
A[1] = 0.000396825396825395*G0_;
3562
A[2] = 0.000396825396825396*G0_;
3563
A[3] = 0.000396825396825396*G0_;
3564
A[4] = -0.00238095238095238*G0_;
3565
A[5] = -0.00238095238095238*G0_;
3566
A[6] = -0.00238095238095238*G0_;
3567
A[7] = -0.00158730158730158*G0_;
3568
A[8] = -0.00158730158730158*G0_;
3569
A[9] = -0.00158730158730158*G0_;
3570
A[10] = 0.000396825396825395*G0_;
3571
A[11] = 0.00238095238095238*G0_;
3572
A[12] = 0.000396825396825396*G0_;
3573
A[13] = 0.000396825396825396*G0_;
3574
A[14] = -0.00238095238095238*G0_;
3575
A[15] = -0.00158730158730158*G0_;
3576
A[16] = -0.00158730158730158*G0_;
3577
A[17] = -0.00238095238095238*G0_;
3578
A[18] = -0.00238095238095238*G0_;
3579
A[19] = -0.00158730158730158*G0_;
3580
A[20] = 0.000396825396825396*G0_;
3581
A[21] = 0.000396825396825396*G0_;
3582
A[22] = 0.00238095238095238*G0_;
3583
A[23] = 0.000396825396825396*G0_;
3584
A[24] = -0.00158730158730159*G0_;
3585
A[25] = -0.00238095238095238*G0_;
3586
A[26] = -0.00158730158730158*G0_;
3587
A[27] = -0.00238095238095238*G0_;
3588
A[28] = -0.00158730158730159*G0_;
3589
A[29] = -0.00238095238095238*G0_;
3590
A[30] = 0.000396825396825396*G0_;
3591
A[31] = 0.000396825396825396*G0_;
3592
A[32] = 0.000396825396825396*G0_;
3593
A[33] = 0.00238095238095238*G0_;
3594
A[34] = -0.00158730158730159*G0_;
3595
A[35] = -0.00158730158730159*G0_;
3596
A[36] = -0.00238095238095238*G0_;
3597
A[37] = -0.00158730158730159*G0_;
3598
A[38] = -0.00238095238095238*G0_;
3599
A[39] = -0.00238095238095238*G0_;
3600
A[40] = -0.00238095238095238*G0_;
3601
A[41] = -0.00238095238095238*G0_;
3602
A[42] = -0.00158730158730159*G0_;
3603
A[43] = -0.00158730158730159*G0_;
3604
A[44] = 0.0126984126984127*G0_;
3605
A[45] = 0.00634920634920635*G0_;
3606
A[46] = 0.00634920634920635*G0_;
3607
A[47] = 0.00634920634920635*G0_;
3608
A[48] = 0.00634920634920634*G0_;
3609
A[49] = 0.00317460317460317*G0_;
3610
A[50] = -0.00238095238095238*G0_;
3611
A[51] = -0.00158730158730158*G0_;
3612
A[52] = -0.00238095238095238*G0_;
3613
A[53] = -0.00158730158730159*G0_;
3614
A[54] = 0.00634920634920635*G0_;
3615
A[55] = 0.0126984126984127*G0_;
3616
A[56] = 0.00634920634920635*G0_;
3617
A[57] = 0.00634920634920635*G0_;
3618
A[58] = 0.00317460317460317*G0_;
3619
A[59] = 0.00634920634920635*G0_;
3620
A[60] = -0.00238095238095238*G0_;
3621
A[61] = -0.00158730158730158*G0_;
3622
A[62] = -0.00158730158730159*G0_;
3623
A[63] = -0.00238095238095238*G0_;
3624
A[64] = 0.00634920634920635*G0_;
3625
A[65] = 0.00634920634920635*G0_;
3626
A[66] = 0.0126984126984127*G0_;
3627
A[67] = 0.00317460317460317*G0_;
3628
A[68] = 0.00634920634920634*G0_;
3629
A[69] = 0.00634920634920635*G0_;
3630
A[70] = -0.00158730158730158*G0_;
3631
A[71] = -0.00238095238095238*G0_;
3632
A[72] = -0.00238095238095238*G0_;
3633
A[73] = -0.00158730158730159*G0_;
3634
A[74] = 0.00634920634920635*G0_;
3635
A[75] = 0.00634920634920635*G0_;
3636
A[76] = 0.00317460317460317*G0_;
3637
A[77] = 0.0126984126984127*G0_;
3638
A[78] = 0.00634920634920635*G0_;
3639
A[79] = 0.00634920634920635*G0_;
3640
A[80] = -0.00158730158730158*G0_;
3641
A[81] = -0.00238095238095238*G0_;
3642
A[82] = -0.00158730158730159*G0_;
3643
A[83] = -0.00238095238095238*G0_;
3644
A[84] = 0.00634920634920634*G0_;
3645
A[85] = 0.00317460317460317*G0_;
3646
A[86] = 0.00634920634920634*G0_;
3647
A[87] = 0.00634920634920635*G0_;
3648
A[88] = 0.0126984126984127*G0_;
3649
A[89] = 0.00634920634920634*G0_;
3650
A[90] = -0.00158730158730158*G0_;
3651
A[91] = -0.00158730158730158*G0_;
3652
A[92] = -0.00238095238095238*G0_;
3653
A[93] = -0.00238095238095238*G0_;
3654
A[94] = 0.00317460317460317*G0_;
3655
A[95] = 0.00634920634920635*G0_;
3656
A[96] = 0.00634920634920635*G0_;
3657
A[97] = 0.00634920634920635*G0_;
3658
A[98] = 0.00634920634920634*G0_;
3659
A[99] = 0.0126984126984127*G0_;
3660
}
3661
3662
};
3663
3664
/// This class defines the interface for the tabulation of the cell
3665
/// tensor corresponding to the local contribution to a form from
3666
/// the integral over a cell.
3667
3668
class projection_cell_integral_1_otherwise: public ufc::cell_integral
3669
{
3670
public:
3671
3672
/// Constructor
3673
projection_cell_integral_1_otherwise() : ufc::cell_integral()
3674
{
3675
// Do nothing
3676
}
3677
3678
/// Destructor
3679
virtual ~projection_cell_integral_1_otherwise()
3680
{
3681
// Do nothing
3682
}
3683
3684
/// Tabulate the tensor for the contribution from a local cell
3685
virtual void tabulate_tensor(double* A,
3686
const double * const * w,
3687
const double* vertex_coordinates,
3688
int cell_orientation) const
3689
{
3690
// Number of operations (multiply-add pairs) for Jacobian data: 3
3691
// Number of operations (multiply-add pairs) for geometry tensor: 10
3692
// Number of operations (multiply-add pairs) for tensor contraction: 95
3693
// Total number of operations (multiply-add pairs): 108
3694
3695
// Compute Jacobian
3696
double J[9];
3697
compute_jacobian_tetrahedron_3d(J, vertex_coordinates);
3698
3699
// Compute Jacobian inverse and determinant
3700
double K[9];
3701
double detJ;
3702
compute_jacobian_inverse_tetrahedron_3d(K, detJ, J);
3703
3704
// Set scale factor
3705
const double det = std::abs(detJ);
3706
3707
// Compute geometry tensor
3708
const double G0_0 = det*w[0][0]*(1.0);
3709
const double G0_1 = det*w[0][1]*(1.0);
3710
const double G0_2 = det*w[0][2]*(1.0);
3711
const double G0_3 = det*w[0][3]*(1.0);
3712
const double G0_4 = det*w[0][4]*(1.0);
3713
const double G0_5 = det*w[0][5]*(1.0);
3714
const double G0_6 = det*w[0][6]*(1.0);
3715
const double G0_7 = det*w[0][7]*(1.0);
3716
const double G0_8 = det*w[0][8]*(1.0);
3717
const double G0_9 = det*w[0][9]*(1.0);
3718
3719
// Compute element tensor
3720
A[0] = 0.00238095238095238*G0_0 + 0.000396825396825395*G0_1 + 0.000396825396825396*G0_2 + 0.000396825396825396*G0_3 - 0.00238095238095238*G0_4 - 0.00238095238095238*G0_5 - 0.00238095238095238*G0_6 - 0.00158730158730158*G0_7 - 0.00158730158730158*G0_8 - 0.00158730158730158*G0_9;
3721
A[1] = 0.000396825396825395*G0_0 + 0.00238095238095238*G0_1 + 0.000396825396825396*G0_2 + 0.000396825396825396*G0_3 - 0.00238095238095238*G0_4 - 0.00158730158730158*G0_5 - 0.00158730158730158*G0_6 - 0.00238095238095238*G0_7 - 0.00238095238095238*G0_8 - 0.00158730158730158*G0_9;
3722
A[2] = 0.000396825396825396*G0_0 + 0.000396825396825396*G0_1 + 0.00238095238095238*G0_2 + 0.000396825396825396*G0_3 - 0.00158730158730159*G0_4 - 0.00238095238095238*G0_5 - 0.00158730158730158*G0_6 - 0.00238095238095238*G0_7 - 0.00158730158730159*G0_8 - 0.00238095238095238*G0_9;
3723
A[3] = 0.000396825396825396*G0_0 + 0.000396825396825396*G0_1 + 0.000396825396825396*G0_2 + 0.00238095238095238*G0_3 - 0.00158730158730159*G0_4 - 0.00158730158730159*G0_5 - 0.00238095238095238*G0_6 - 0.00158730158730159*G0_7 - 0.00238095238095238*G0_8 - 0.00238095238095238*G0_9;
3724
A[4] = -0.00238095238095238*G0_0 - 0.00238095238095238*G0_1 - 0.00158730158730159*G0_2 - 0.00158730158730159*G0_3 + 0.0126984126984127*G0_4 + 0.00634920634920635*G0_5 + 0.00634920634920635*G0_6 + 0.00634920634920635*G0_7 + 0.00634920634920634*G0_8 + 0.00317460317460317*G0_9;
3725
A[5] = -0.00238095238095238*G0_0 - 0.00158730158730158*G0_1 - 0.00238095238095238*G0_2 - 0.00158730158730159*G0_3 + 0.00634920634920635*G0_4 + 0.0126984126984127*G0_5 + 0.00634920634920635*G0_6 + 0.00634920634920635*G0_7 + 0.00317460317460317*G0_8 + 0.00634920634920635*G0_9;
3726
A[6] = -0.00238095238095238*G0_0 - 0.00158730158730158*G0_1 - 0.00158730158730159*G0_2 - 0.00238095238095238*G0_3 + 0.00634920634920635*G0_4 + 0.00634920634920635*G0_5 + 0.0126984126984127*G0_6 + 0.00317460317460317*G0_7 + 0.00634920634920634*G0_8 + 0.00634920634920635*G0_9;
3727
A[7] = -0.00158730158730158*G0_0 - 0.00238095238095238*G0_1 - 0.00238095238095238*G0_2 - 0.00158730158730159*G0_3 + 0.00634920634920635*G0_4 + 0.00634920634920635*G0_5 + 0.00317460317460317*G0_6 + 0.0126984126984127*G0_7 + 0.00634920634920635*G0_8 + 0.00634920634920635*G0_9;
3728
A[8] = -0.00158730158730158*G0_0 - 0.00238095238095238*G0_1 - 0.00158730158730159*G0_2 - 0.00238095238095238*G0_3 + 0.00634920634920634*G0_4 + 0.00317460317460317*G0_5 + 0.00634920634920634*G0_6 + 0.00634920634920635*G0_7 + 0.0126984126984127*G0_8 + 0.00634920634920634*G0_9;
3729
A[9] = -0.00158730158730158*G0_0 - 0.00158730158730158*G0_1 - 0.00238095238095238*G0_2 - 0.00238095238095238*G0_3 + 0.00317460317460317*G0_4 + 0.00634920634920635*G0_5 + 0.00634920634920635*G0_6 + 0.00634920634920635*G0_7 + 0.00634920634920634*G0_8 + 0.0126984126984127*G0_9;
3730
}
3731
3732
};
3733
3734
/// This class defines the interface for the assembly of the global
3735
/// tensor corresponding to a form with r + n arguments, that is, a
3736
/// mapping
3737
///
3738
/// a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
3739
///
3740
/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
3741
/// global tensor A is defined by
3742
///
3743
/// A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
3744
///
3745
/// where each argument Vj represents the application to the
3746
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
3747
/// fixed functions (coefficients).
3748
3749
class projection_form_0: public ufc::form
3750
{
3751
public:
3752
3753
/// Constructor
3754
projection_form_0() : ufc::form()
3755
{
3756
// Do nothing
3757
}
3758
3759
/// Destructor
3760
virtual ~projection_form_0()
3761
{
3762
// Do nothing
3763
}
3764
3765
/// Return a string identifying the form
3766
virtual const char* signature() const
3767
{
3768
return "d784c0e2a973a24797958f6509b56c50f2578240181e8ff00098270cf2fa2972b67951ef90f3a43d309d2ab5f41d7fefc9ce77a68c5d4fca8664b99cf5789190";
3769
}
3770
3771
/// Return the rank of the global tensor (r)
3772
virtual std::size_t rank() const
3773
{
3774
return 2;
3775
}
3776
3777
/// Return the number of coefficients (n)
3778
virtual std::size_t num_coefficients() const
3779
{
3780
return 0;
3781
}
3782
3783
/// Return the number of cell domains
3784
virtual std::size_t num_cell_domains() const
3785
{
3786
return 0;
3787
}
3788
3789
/// Return the number of exterior facet domains
3790
virtual std::size_t num_exterior_facet_domains() const
3791
{
3792
return 0;
3793
}
3794
3795
/// Return the number of interior facet domains
3796
virtual std::size_t num_interior_facet_domains() const
3797
{
3798
return 0;
3799
}
3800
3801
/// Return the number of point domains
3802
virtual std::size_t num_point_domains() const
3803
{
3804
return 0;
3805
}
3806
3807
/// Return whether the form has any cell integrals
3808
virtual bool has_cell_integrals() const
3809
{
3810
return true;
3811
}
3812
3813
/// Return whether the form has any exterior facet integrals
3814
virtual bool has_exterior_facet_integrals() const
3815
{
3816
return false;
3817
}
3818
3819
/// Return whether the form has any interior facet integrals
3820
virtual bool has_interior_facet_integrals() const
3821
{
3822
return false;
3823
}
3824
3825
/// Return whether the form has any point integrals
3826
virtual bool has_point_integrals() const
3827
{
3828
return false;
3829
}
3830
3831
/// Create a new finite element for argument function i
3832
virtual ufc::finite_element* create_finite_element(std::size_t i) const
3833
{
3834
switch (i)
3835
{
3836
case 0:
3837
{
3838
return new projection_finite_element_0();
3839
break;
3840
}
3841
case 1:
3842
{
3843
return new projection_finite_element_0();
3844
break;
3845
}
3846
}
3847
3848
return 0;
3849
}
3850
3851
/// Create a new dofmap for argument function i
3852
virtual ufc::dofmap* create_dofmap(std::size_t i) const
3853
{
3854
switch (i)
3855
{
3856
case 0:
3857
{
3858
return new projection_dofmap_0();
3859
break;
3860
}
3861
case 1:
3862
{
3863
return new projection_dofmap_0();
3864
break;
3865
}
3866
}
3867
3868
return 0;
3869
}
3870
3871
/// Create a new cell integral on sub domain i
3872
virtual ufc::cell_integral* create_cell_integral(std::size_t i) const
3873
{
3874
return 0;
3875
}
3876
3877
/// Create a new exterior facet integral on sub domain i
3878
virtual ufc::exterior_facet_integral* create_exterior_facet_integral(std::size_t i) const
3879
{
3880
return 0;
3881
}
3882
3883
/// Create a new interior facet integral on sub domain i
3884
virtual ufc::interior_facet_integral* create_interior_facet_integral(std::size_t i) const
3885
{
3886
return 0;
3887
}
3888
3889
/// Create a new point integral on sub domain i
3890
virtual ufc::point_integral* create_point_integral(std::size_t i) const
3891
{
3892
return 0;
3893
}
3894
3895
/// Create a new cell integral on everywhere else
3896
virtual ufc::cell_integral* create_default_cell_integral() const
3897
{
3898
return new projection_cell_integral_0_otherwise();
3899
}
3900
3901
/// Create a new exterior facet integral on everywhere else
3902
virtual ufc::exterior_facet_integral* create_default_exterior_facet_integral() const
3903
{
3904
return 0;
3905
}
3906
3907
/// Create a new interior facet integral on everywhere else
3908
virtual ufc::interior_facet_integral* create_default_interior_facet_integral() const
3909
{
3910
return 0;
3911
}
3912
3913
/// Create a new point integral on everywhere else
3914
virtual ufc::point_integral* create_default_point_integral() const
3915
{
3916
return 0;
3917
}
3918
3919
};
3920
3921
/// This class defines the interface for the assembly of the global
3922
/// tensor corresponding to a form with r + n arguments, that is, a
3923
/// mapping
3924
///
3925
/// a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
3926
///
3927
/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
3928
/// global tensor A is defined by
3929
///
3930
/// A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
3931
///
3932
/// where each argument Vj represents the application to the
3933
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
3934
/// fixed functions (coefficients).
3935
3936
class projection_form_1: public ufc::form
3937
{
3938
public:
3939
3940
/// Constructor
3941
projection_form_1() : ufc::form()
3942
{
3943
// Do nothing
3944
}
3945
3946
/// Destructor
3947
virtual ~projection_form_1()
3948
{
3949
// Do nothing
3950
}
3951
3952
/// Return a string identifying the form
3953
virtual const char* signature() const
3954
{
3955
return "5e05fa9312df6ed7e318c7a07aff232acbdcbfa3737e1a69ccbada9505317bbff5ea736d79a7854db7a559f13450590edc4c8477ff3cec5390f13f6a1f5804bd";
3956
}
3957
3958
/// Return the rank of the global tensor (r)
3959
virtual std::size_t rank() const
3960
{
3961
return 1;
3962
}
3963
3964
/// Return the number of coefficients (n)
3965
virtual std::size_t num_coefficients() const
3966
{
3967
return 1;
3968
}
3969
3970
/// Return the number of cell domains
3971
virtual std::size_t num_cell_domains() const
3972
{
3973
return 0;
3974
}
3975
3976
/// Return the number of exterior facet domains
3977
virtual std::size_t num_exterior_facet_domains() const
3978
{
3979
return 0;
3980
}
3981
3982
/// Return the number of interior facet domains
3983
virtual std::size_t num_interior_facet_domains() const
3984
{
3985
return 0;
3986
}
3987
3988
/// Return the number of point domains
3989
virtual std::size_t num_point_domains() const
3990
{
3991
return 0;
3992
}
3993
3994
/// Return whether the form has any cell integrals
3995
virtual bool has_cell_integrals() const
3996
{
3997
return true;
3998
}
3999
4000
/// Return whether the form has any exterior facet integrals
4001
virtual bool has_exterior_facet_integrals() const
4002
{
4003
return false;
4004
}
4005
4006
/// Return whether the form has any interior facet integrals
4007
virtual bool has_interior_facet_integrals() const
4008
{
4009
return false;
4010
}
4011
4012
/// Return whether the form has any point integrals
4013
virtual bool has_point_integrals() const
4014
{
4015
return false;
4016
}
4017
4018
/// Create a new finite element for argument function i
4019
virtual ufc::finite_element* create_finite_element(std::size_t i) const
4020
{
4021
switch (i)
4022
{
4023
case 0:
4024
{
4025
return new projection_finite_element_0();
4026
break;
4027
}
4028
case 1:
4029
{
4030
return new projection_finite_element_0();
4031
break;
4032
}
4033
}
4034
4035
return 0;
4036
}
4037
4038
/// Create a new dofmap for argument function i
4039
virtual ufc::dofmap* create_dofmap(std::size_t i) const
4040
{
4041
switch (i)
4042
{
4043
case 0:
4044
{
4045
return new projection_dofmap_0();
4046
break;
4047
}
4048
case 1:
4049
{
4050
return new projection_dofmap_0();
4051
break;
4052
}
4053
}
4054
4055
return 0;
4056
}
4057
4058
/// Create a new cell integral on sub domain i
4059
virtual ufc::cell_integral* create_cell_integral(std::size_t i) const
4060
{
4061
return 0;
4062
}
4063
4064
/// Create a new exterior facet integral on sub domain i
4065
virtual ufc::exterior_facet_integral* create_exterior_facet_integral(std::size_t i) const
4066
{
4067
return 0;
4068
}
4069
4070
/// Create a new interior facet integral on sub domain i
4071
virtual ufc::interior_facet_integral* create_interior_facet_integral(std::size_t i) const
4072
{
4073
return 0;
4074
}
4075
4076
/// Create a new point integral on sub domain i
4077
virtual ufc::point_integral* create_point_integral(std::size_t i) const
4078
{
4079
return 0;
4080
}
4081
4082
/// Create a new cell integral on everywhere else
4083
virtual ufc::cell_integral* create_default_cell_integral() const
4084
{
4085
return new projection_cell_integral_1_otherwise();
4086
}
4087
4088
/// Create a new exterior facet integral on everywhere else
4089
virtual ufc::exterior_facet_integral* create_default_exterior_facet_integral() const
4090
{
4091
return 0;
4092
}
4093
4094
/// Create a new interior facet integral on everywhere else
4095
virtual ufc::interior_facet_integral* create_default_interior_facet_integral() const
4096
{
4097
return 0;
4098
}
4099
4100
/// Create a new point integral on everywhere else
4101
virtual ufc::point_integral* create_default_point_integral() const
4102
{
4103
return 0;
4104
}
4105
4106
};
4107
4108
// DOLFIN wrappers
4109
4110
// Standard library includes
4111
#include <string>
4112
4113
// DOLFIN includes
4114
#include <dolfin/common/NoDeleter.h>
4115
#include <dolfin/mesh/Restriction.h>
4116
#include <dolfin/fem/FiniteElement.h>
4117
#include <dolfin/fem/DofMap.h>
4118
#include <dolfin/fem/Form.h>
4119
#include <dolfin/function/FunctionSpace.h>
4120
#include <dolfin/function/GenericFunction.h>
4121
#include <dolfin/function/CoefficientAssigner.h>
4122
#include <dolfin/adaptivity/ErrorControl.h>
4123
#include <dolfin/adaptivity/GoalFunctional.h>
4124
4125
namespace Projection
4126
{
4127
4128
class CoefficientSpace_f: public dolfin::FunctionSpace
4129
{
4130
public:
4131
4132
//--- Constructors for standard function space, 2 different versions ---
4133
4134
// Create standard function space (reference version)
4135
CoefficientSpace_f(const dolfin::Mesh& mesh):
4136
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
4137
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new projection_finite_element_0()))),
4138
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new projection_dofmap_0()), mesh)))
4139
{
4140
// Do nothing
4141
}
4142
4143
// Create standard function space (shared pointer version)
4144
CoefficientSpace_f(boost::shared_ptr<const dolfin::Mesh> mesh):
4145
dolfin::FunctionSpace(mesh,
4146
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new projection_finite_element_0()))),
4147
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new projection_dofmap_0()), *mesh)))
4148
{
4149
// Do nothing
4150
}
4151
4152
//--- Constructors for constrained function space, 2 different versions ---
4153
4154
// Create standard function space (reference version)
4155
CoefficientSpace_f(const dolfin::Mesh& mesh, const dolfin::SubDomain& constrained_domain):
4156
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
4157
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new projection_finite_element_0()))),
4158
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new projection_dofmap_0()), mesh,
4159
dolfin::reference_to_no_delete_pointer(constrained_domain))))
4160
{
4161
// Do nothing
4162
}
4163
4164
// Create standard function space (shared pointer version)
4165
CoefficientSpace_f(boost::shared_ptr<const dolfin::Mesh> mesh, boost::shared_ptr<const dolfin::SubDomain> constrained_domain):
4166
dolfin::FunctionSpace(mesh,
4167
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new projection_finite_element_0()))),
4168
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new projection_dofmap_0()), *mesh, constrained_domain)))
4169
{
4170
// Do nothing
4171
}
4172
4173
//--- Constructors for restricted function space, 2 different versions ---
4174
4175
// Create restricted function space (reference version)
4176
CoefficientSpace_f(const dolfin::Restriction& restriction):
4177
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(restriction.mesh()),
4178
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new projection_finite_element_0()))),
4179
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new projection_dofmap_0()),
4180
reference_to_no_delete_pointer(restriction))))
4181
{
4182
// Do nothing
4183
}
4184
4185
// Create restricted function space (shared pointer version)
4186
CoefficientSpace_f(boost::shared_ptr<const dolfin::Restriction> restriction):
4187
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(restriction->mesh()),
4188
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new projection_finite_element_0()))),
4189
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new projection_dofmap_0()),
4190
restriction)))
4191
{
4192
// Do nothing
4193
}
4194
4195
// Copy constructor
4196
~CoefficientSpace_f()
4197
{
4198
}
4199
4200
};
4201
4202
class Form_a_FunctionSpace_0: public dolfin::FunctionSpace
4203
{
4204
public:
4205
4206
//--- Constructors for standard function space, 2 different versions ---
4207
4208
// Create standard function space (reference version)
4209
Form_a_FunctionSpace_0(const dolfin::Mesh& mesh):
4210
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
4211
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new projection_finite_element_0()))),
4212
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new projection_dofmap_0()), mesh)))
4213
{
4214
// Do nothing
4215
}
4216
4217
// Create standard function space (shared pointer version)
4218
Form_a_FunctionSpace_0(boost::shared_ptr<const dolfin::Mesh> mesh):
4219
dolfin::FunctionSpace(mesh,
4220
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new projection_finite_element_0()))),
4221
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new projection_dofmap_0()), *mesh)))
4222
{
4223
// Do nothing
4224
}
4225
4226
//--- Constructors for constrained function space, 2 different versions ---
4227
4228
// Create standard function space (reference version)
4229
Form_a_FunctionSpace_0(const dolfin::Mesh& mesh, const dolfin::SubDomain& constrained_domain):
4230
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
4231
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new projection_finite_element_0()))),
4232
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new projection_dofmap_0()), mesh,
4233
dolfin::reference_to_no_delete_pointer(constrained_domain))))
4234
{
4235
// Do nothing
4236
}
4237
4238
// Create standard function space (shared pointer version)
4239
Form_a_FunctionSpace_0(boost::shared_ptr<const dolfin::Mesh> mesh, boost::shared_ptr<const dolfin::SubDomain> constrained_domain):
4240
dolfin::FunctionSpace(mesh,
4241
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new projection_finite_element_0()))),
4242
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new projection_dofmap_0()), *mesh, constrained_domain)))
4243
{
4244
// Do nothing
4245
}
4246
4247
//--- Constructors for restricted function space, 2 different versions ---
4248
4249
// Create restricted function space (reference version)
4250
Form_a_FunctionSpace_0(const dolfin::Restriction& restriction):
4251
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(restriction.mesh()),
4252
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new projection_finite_element_0()))),
4253
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new projection_dofmap_0()),
4254
reference_to_no_delete_pointer(restriction))))
4255
{
4256
// Do nothing
4257
}
4258
4259
// Create restricted function space (shared pointer version)
4260
Form_a_FunctionSpace_0(boost::shared_ptr<const dolfin::Restriction> restriction):
4261
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(restriction->mesh()),
4262
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new projection_finite_element_0()))),
4263
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new projection_dofmap_0()),
4264
restriction)))
4265
{
4266
// Do nothing
4267
}
4268
4269
// Copy constructor
4270
~Form_a_FunctionSpace_0()
4271
{
4272
}
4273
4274
};
4275
4276
class Form_a_FunctionSpace_1: public dolfin::FunctionSpace
4277
{
4278
public:
4279
4280
//--- Constructors for standard function space, 2 different versions ---
4281
4282
// Create standard function space (reference version)
4283
Form_a_FunctionSpace_1(const dolfin::Mesh& mesh):
4284
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
4285
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new projection_finite_element_0()))),
4286
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new projection_dofmap_0()), mesh)))
4287
{
4288
// Do nothing
4289
}
4290
4291
// Create standard function space (shared pointer version)
4292
Form_a_FunctionSpace_1(boost::shared_ptr<const dolfin::Mesh> mesh):
4293
dolfin::FunctionSpace(mesh,
4294
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new projection_finite_element_0()))),
4295
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new projection_dofmap_0()), *mesh)))
4296
{
4297
// Do nothing
4298
}
4299
4300
//--- Constructors for constrained function space, 2 different versions ---
4301
4302
// Create standard function space (reference version)
4303
Form_a_FunctionSpace_1(const dolfin::Mesh& mesh, const dolfin::SubDomain& constrained_domain):
4304
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
4305
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new projection_finite_element_0()))),
4306
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new projection_dofmap_0()), mesh,
4307
dolfin::reference_to_no_delete_pointer(constrained_domain))))
4308
{
4309
// Do nothing
4310
}
4311
4312
// Create standard function space (shared pointer version)
4313
Form_a_FunctionSpace_1(boost::shared_ptr<const dolfin::Mesh> mesh, boost::shared_ptr<const dolfin::SubDomain> constrained_domain):
4314
dolfin::FunctionSpace(mesh,
4315
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new projection_finite_element_0()))),
4316
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new projection_dofmap_0()), *mesh, constrained_domain)))
4317
{
4318
// Do nothing
4319
}
4320
4321
//--- Constructors for restricted function space, 2 different versions ---
4322
4323
// Create restricted function space (reference version)
4324
Form_a_FunctionSpace_1(const dolfin::Restriction& restriction):
4325
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(restriction.mesh()),
4326
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new projection_finite_element_0()))),
4327
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new projection_dofmap_0()),
4328
reference_to_no_delete_pointer(restriction))))
4329
{
4330
// Do nothing
4331
}
4332
4333
// Create restricted function space (shared pointer version)
4334
Form_a_FunctionSpace_1(boost::shared_ptr<const dolfin::Restriction> restriction):
4335
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(restriction->mesh()),
4336
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new projection_finite_element_0()))),
4337
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new projection_dofmap_0()),
4338
restriction)))
4339
{
4340
// Do nothing
4341
}
4342
4343
// Copy constructor
4344
~Form_a_FunctionSpace_1()
4345
{
4346
}
4347
4348
};
4349
4350
class Form_a: public dolfin::Form
4351
{
4352
public:
4353
4354
// Constructor
4355
Form_a(const dolfin::FunctionSpace& V1, const dolfin::FunctionSpace& V0):
4356
dolfin::Form(2, 0)
4357
{
4358
_function_spaces[0] = reference_to_no_delete_pointer(V0);
4359
_function_spaces[1] = reference_to_no_delete_pointer(V1);
4360
4361
_ufc_form = boost::shared_ptr<const ufc::form>(new projection_form_0());
4362
}
4363
4364
// Constructor
4365
Form_a(boost::shared_ptr<const dolfin::FunctionSpace> V1, boost::shared_ptr<const dolfin::FunctionSpace> V0):
4366
dolfin::Form(2, 0)
4367
{
4368
_function_spaces[0] = V0;
4369
_function_spaces[1] = V1;
4370
4371
_ufc_form = boost::shared_ptr<const ufc::form>(new projection_form_0());
4372
}
4373
4374
// Destructor
4375
~Form_a()
4376
{}
4377
4378
/// Return the number of the coefficient with this name
4379
virtual std::size_t coefficient_number(const std::string& name) const
4380
{
4381
4382
dolfin::dolfin_error("generated code for class Form",
4383
"access coefficient data",
4384
"There are no coefficients");
4385
return 0;
4386
}
4387
4388
/// Return the name of the coefficient with this number
4389
virtual std::string coefficient_name(std::size_t i) const
4390
{
4391
4392
dolfin::dolfin_error("generated code for class Form",
4393
"access coefficient data",
4394
"There are no coefficients");
4395
return "unnamed";
4396
}
4397
4398
// Typedefs
4399
typedef Form_a_FunctionSpace_0 TestSpace;
4400
typedef Form_a_FunctionSpace_1 TrialSpace;
4401
4402
// Coefficients
4403
};
4404
4405
class Form_L_FunctionSpace_0: public dolfin::FunctionSpace
4406
{
4407
public:
4408
4409
//--- Constructors for standard function space, 2 different versions ---
4410
4411
// Create standard function space (reference version)
4412
Form_L_FunctionSpace_0(const dolfin::Mesh& mesh):
4413
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
4414
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new projection_finite_element_0()))),
4415
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new projection_dofmap_0()), mesh)))
4416
{
4417
// Do nothing
4418
}
4419
4420
// Create standard function space (shared pointer version)
4421
Form_L_FunctionSpace_0(boost::shared_ptr<const dolfin::Mesh> mesh):
4422
dolfin::FunctionSpace(mesh,
4423
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new projection_finite_element_0()))),
4424
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new projection_dofmap_0()), *mesh)))
4425
{
4426
// Do nothing
4427
}
4428
4429
//--- Constructors for constrained function space, 2 different versions ---
4430
4431
// Create standard function space (reference version)
4432
Form_L_FunctionSpace_0(const dolfin::Mesh& mesh, const dolfin::SubDomain& constrained_domain):
4433
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
4434
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new projection_finite_element_0()))),
4435
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new projection_dofmap_0()), mesh,
4436
dolfin::reference_to_no_delete_pointer(constrained_domain))))
4437
{
4438
// Do nothing
4439
}
4440
4441
// Create standard function space (shared pointer version)
4442
Form_L_FunctionSpace_0(boost::shared_ptr<const dolfin::Mesh> mesh, boost::shared_ptr<const dolfin::SubDomain> constrained_domain):
4443
dolfin::FunctionSpace(mesh,
4444
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new projection_finite_element_0()))),
4445
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new projection_dofmap_0()), *mesh, constrained_domain)))
4446
{
4447
// Do nothing
4448
}
4449
4450
//--- Constructors for restricted function space, 2 different versions ---
4451
4452
// Create restricted function space (reference version)
4453
Form_L_FunctionSpace_0(const dolfin::Restriction& restriction):
4454
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(restriction.mesh()),
4455
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new projection_finite_element_0()))),
4456
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new projection_dofmap_0()),
4457
reference_to_no_delete_pointer(restriction))))
4458
{
4459
// Do nothing
4460
}
4461
4462
// Create restricted function space (shared pointer version)
4463
Form_L_FunctionSpace_0(boost::shared_ptr<const dolfin::Restriction> restriction):
4464
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(restriction->mesh()),
4465
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new projection_finite_element_0()))),
4466
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new projection_dofmap_0()),
4467
restriction)))
4468
{
4469
// Do nothing
4470
}
4471
4472
// Copy constructor
4473
~Form_L_FunctionSpace_0()
4474
{
4475
}
4476
4477
};
4478
4479
typedef CoefficientSpace_f Form_L_FunctionSpace_1;
4480
4481
class Form_L: public dolfin::Form
4482
{
4483
public:
4484
4485
// Constructor
4486
Form_L(const dolfin::FunctionSpace& V0):
4487
dolfin::Form(1, 1), f(*this, 0)
4488
{
4489
_function_spaces[0] = reference_to_no_delete_pointer(V0);
4490
4491
_ufc_form = boost::shared_ptr<const ufc::form>(new projection_form_1());
4492
}
4493
4494
// Constructor
4495
Form_L(const dolfin::FunctionSpace& V0, const dolfin::GenericFunction& f):
4496
dolfin::Form(1, 1), f(*this, 0)
4497
{
4498
_function_spaces[0] = reference_to_no_delete_pointer(V0);
4499
4500
this->f = f;
4501
4502
_ufc_form = boost::shared_ptr<const ufc::form>(new projection_form_1());
4503
}
4504
4505
// Constructor
4506
Form_L(const dolfin::FunctionSpace& V0, boost::shared_ptr<const dolfin::GenericFunction> f):
4507
dolfin::Form(1, 1), f(*this, 0)
4508
{
4509
_function_spaces[0] = reference_to_no_delete_pointer(V0);
4510
4511
this->f = *f;
4512
4513
_ufc_form = boost::shared_ptr<const ufc::form>(new projection_form_1());
4514
}
4515
4516
// Constructor
4517
Form_L(boost::shared_ptr<const dolfin::FunctionSpace> V0):
4518
dolfin::Form(1, 1), f(*this, 0)
4519
{
4520
_function_spaces[0] = V0;
4521
4522
_ufc_form = boost::shared_ptr<const ufc::form>(new projection_form_1());
4523
}
4524
4525
// Constructor
4526
Form_L(boost::shared_ptr<const dolfin::FunctionSpace> V0, const dolfin::GenericFunction& f):
4527
dolfin::Form(1, 1), f(*this, 0)
4528
{
4529
_function_spaces[0] = V0;
4530
4531
this->f = f;
4532
4533
_ufc_form = boost::shared_ptr<const ufc::form>(new projection_form_1());
4534
}
4535
4536
// Constructor
4537
Form_L(boost::shared_ptr<const dolfin::FunctionSpace> V0, boost::shared_ptr<const dolfin::GenericFunction> f):
4538
dolfin::Form(1, 1), f(*this, 0)
4539
{
4540
_function_spaces[0] = V0;
4541
4542
this->f = *f;
4543
4544
_ufc_form = boost::shared_ptr<const ufc::form>(new projection_form_1());
4545
}
4546
4547
// Destructor
4548
~Form_L()
4549
{}
4550
4551
/// Return the number of the coefficient with this name
4552
virtual std::size_t coefficient_number(const std::string& name) const
4553
{
4554
if (name == "f")
4555
return 0;
4556
4557
dolfin::dolfin_error("generated code for class Form",
4558
"access coefficient data",
4559
"Invalid coefficient");
4560
return 0;
4561
}
4562
4563
/// Return the name of the coefficient with this number
4564
virtual std::string coefficient_name(std::size_t i) const
4565
{
4566
switch (i)
4567
{
4568
case 0:
4569
return "f";
4570
}
4571
4572
dolfin::dolfin_error("generated code for class Form",
4573
"access coefficient data",
4574
"Invalid coefficient");
4575
return "unnamed";
4576
}
4577
4578
// Typedefs
4579
typedef Form_L_FunctionSpace_0 TestSpace;
4580
typedef Form_L_FunctionSpace_1 CoefficientSpace_f;
4581
4582
// Coefficients
4583
dolfin::CoefficientAssigner f;
4584
};
4585
4586
// Class typedefs
4587
typedef Form_a BilinearForm;
4588
typedef Form_a JacobianForm;
4589
typedef Form_L LinearForm;
4590
typedef Form_L ResidualForm;
4591
typedef Form_a::TestSpace FunctionSpace;
4592
4593
}
4594
4595
#endif
Loggerhead is a web-based interface for Breezy
Version: 2.0.1