~johan-hake/dolfin/general-rk-solver
» Revision
7386
: demo/undocumented/nonmatching-interpolation/cpp/P3.h
« back to all changes in this revision
Viewing changes to demo/undocumented/nonmatching-interpolation/cpp/P3.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
demo/undocumented/nonmatching-interpolation/cpp/P3.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 __P3_H
31
#define __P3_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 p3_finite_element_0: public ufc::finite_element
41
{
42
public:
43
44
/// Constructor
45
p3_finite_element_0() : ufc::finite_element()
46
{
47
// Do nothing
48
}
49
50
/// Destructor
51
virtual ~p3_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('triangle', 2), 'triangle_multiverse', 2, 2), 3, None)";
60
}
61
62
/// Return the cell shape
63
virtual ufc::shape cell_shape() const
64
{
65
return ufc::triangle;
66
}
67
68
/// Return the topological dimension of the cell shape
69
virtual std::size_t topological_dimension() const
70
{
71
return 2;
72
}
73
74
/// Return the geometric dimension of the cell shape
75
virtual std::size_t geometric_dimension() const
76
{
77
return 2;
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[4];
107
compute_jacobian_triangle_2d(J, vertex_coordinates);
108
109
// Compute Jacobian inverse and determinant
110
double K[4];
111
double detJ;
112
compute_jacobian_inverse_triangle_2d(K, detJ, J);
113
114
115
// Compute constants
116
const double C0 = vertex_coordinates[2] + vertex_coordinates[4];
117
const double C1 = vertex_coordinates[3] + vertex_coordinates[5];
118
119
// Get coordinates and map to the reference (FIAT) element
120
double X = (J[1]*(C1 - 2.0*x[1]) + J[3]*(2.0*x[0] - C0)) / detJ;
121
double Y = (J[0]*(2.0*x[1] - C1) + J[2]*(C0 - 2.0*x[0])) / detJ;
122
123
// Reset values
124
*values = 0.0;
125
switch (i)
126
{
127
case 0:
128
{
129
130
// Array of basisvalues
131
double basisvalues[10] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
132
133
// Declare helper variables
134
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
135
double tmp1 = (1.0 - Y)/2.0;
136
double tmp2 = tmp1*tmp1;
137
138
// Compute basisvalues
139
basisvalues[0] = 1.0;
140
basisvalues[1] = tmp0;
141
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
142
basisvalues[6] = basisvalues[3]*1.66666666666667*tmp0 - basisvalues[1]*0.666666666666667*tmp2;
143
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
144
basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
145
basisvalues[7] = basisvalues[3]*(2.5 + 3.5*Y);
146
basisvalues[5] = basisvalues[2]*(0.111111111111111 + Y*1.66666666666667) - basisvalues[0]*0.555555555555556;
147
basisvalues[9] = basisvalues[5]*(0.05 + Y*1.75) - basisvalues[2]*0.7;
148
basisvalues[8] = basisvalues[4]*(0.54 + Y*2.1) - basisvalues[1]*0.56;
149
basisvalues[0] *= std::sqrt(0.5);
150
basisvalues[2] *= std::sqrt(1.0);
151
basisvalues[5] *= std::sqrt(1.5);
152
basisvalues[9] *= std::sqrt(2.0);
153
basisvalues[1] *= std::sqrt(3.0);
154
basisvalues[4] *= std::sqrt(4.5);
155
basisvalues[8] *= std::sqrt(6.0);
156
basisvalues[3] *= std::sqrt(7.5);
157
basisvalues[7] *= std::sqrt(10.0);
158
basisvalues[6] *= std::sqrt(14.0);
159
160
// Table(s) of coefficients
161
static const double coefficients0[10] = \
162
{0.0471404520791032, -0.0288675134594813, -0.0166666666666667, 0.0782460796435951, 0.0606091526731326, 0.0349927106111883, -0.0601337794302955, -0.0508223195384204, -0.0393667994375868, -0.0227284322524248};
163
164
// Compute value(s)
165
for (unsigned int r = 0; r < 10; r++)
166
{
167
*values += coefficients0[r]*basisvalues[r];
168
}// end loop over 'r'
169
break;
170
}
171
case 1:
172
{
173
174
// Array of basisvalues
175
double basisvalues[10] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
176
177
// Declare helper variables
178
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
179
double tmp1 = (1.0 - Y)/2.0;
180
double tmp2 = tmp1*tmp1;
181
182
// Compute basisvalues
183
basisvalues[0] = 1.0;
184
basisvalues[1] = tmp0;
185
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
186
basisvalues[6] = basisvalues[3]*1.66666666666667*tmp0 - basisvalues[1]*0.666666666666667*tmp2;
187
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
188
basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
189
basisvalues[7] = basisvalues[3]*(2.5 + 3.5*Y);
190
basisvalues[5] = basisvalues[2]*(0.111111111111111 + Y*1.66666666666667) - basisvalues[0]*0.555555555555556;
191
basisvalues[9] = basisvalues[5]*(0.05 + Y*1.75) - basisvalues[2]*0.7;
192
basisvalues[8] = basisvalues[4]*(0.54 + Y*2.1) - basisvalues[1]*0.56;
193
basisvalues[0] *= std::sqrt(0.5);
194
basisvalues[2] *= std::sqrt(1.0);
195
basisvalues[5] *= std::sqrt(1.5);
196
basisvalues[9] *= std::sqrt(2.0);
197
basisvalues[1] *= std::sqrt(3.0);
198
basisvalues[4] *= std::sqrt(4.5);
199
basisvalues[8] *= std::sqrt(6.0);
200
basisvalues[3] *= std::sqrt(7.5);
201
basisvalues[7] *= std::sqrt(10.0);
202
basisvalues[6] *= std::sqrt(14.0);
203
204
// Table(s) of coefficients
205
static const double coefficients0[10] = \
206
{0.0471404520791032, 0.0288675134594813, -0.0166666666666666, 0.0782460796435952, -0.0606091526731327, 0.0349927106111883, 0.0601337794302955, -0.0508223195384204, 0.0393667994375868, -0.0227284322524248};
207
208
// Compute value(s)
209
for (unsigned int r = 0; r < 10; r++)
210
{
211
*values += coefficients0[r]*basisvalues[r];
212
}// end loop over 'r'
213
break;
214
}
215
case 2:
216
{
217
218
// Array of basisvalues
219
double basisvalues[10] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
220
221
// Declare helper variables
222
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
223
double tmp1 = (1.0 - Y)/2.0;
224
double tmp2 = tmp1*tmp1;
225
226
// Compute basisvalues
227
basisvalues[0] = 1.0;
228
basisvalues[1] = tmp0;
229
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
230
basisvalues[6] = basisvalues[3]*1.66666666666667*tmp0 - basisvalues[1]*0.666666666666667*tmp2;
231
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
232
basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
233
basisvalues[7] = basisvalues[3]*(2.5 + 3.5*Y);
234
basisvalues[5] = basisvalues[2]*(0.111111111111111 + Y*1.66666666666667) - basisvalues[0]*0.555555555555556;
235
basisvalues[9] = basisvalues[5]*(0.05 + Y*1.75) - basisvalues[2]*0.7;
236
basisvalues[8] = basisvalues[4]*(0.54 + Y*2.1) - basisvalues[1]*0.56;
237
basisvalues[0] *= std::sqrt(0.5);
238
basisvalues[2] *= std::sqrt(1.0);
239
basisvalues[5] *= std::sqrt(1.5);
240
basisvalues[9] *= std::sqrt(2.0);
241
basisvalues[1] *= std::sqrt(3.0);
242
basisvalues[4] *= std::sqrt(4.5);
243
basisvalues[8] *= std::sqrt(6.0);
244
basisvalues[3] *= std::sqrt(7.5);
245
basisvalues[7] *= std::sqrt(10.0);
246
basisvalues[6] *= std::sqrt(14.0);
247
248
// Table(s) of coefficients
249
static const double coefficients0[10] = \
250
{0.0471404520791032, 0.0, 0.0333333333333334, 0.0, 0.0, 0.104978131833565, 0.0, 0.0, 0.0, 0.0909137290096989};
251
252
// Compute value(s)
253
for (unsigned int r = 0; r < 10; r++)
254
{
255
*values += coefficients0[r]*basisvalues[r];
256
}// end loop over 'r'
257
break;
258
}
259
case 3:
260
{
261
262
// Array of basisvalues
263
double basisvalues[10] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
264
265
// Declare helper variables
266
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
267
double tmp1 = (1.0 - Y)/2.0;
268
double tmp2 = tmp1*tmp1;
269
270
// Compute basisvalues
271
basisvalues[0] = 1.0;
272
basisvalues[1] = tmp0;
273
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
274
basisvalues[6] = basisvalues[3]*1.66666666666667*tmp0 - basisvalues[1]*0.666666666666667*tmp2;
275
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
276
basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
277
basisvalues[7] = basisvalues[3]*(2.5 + 3.5*Y);
278
basisvalues[5] = basisvalues[2]*(0.111111111111111 + Y*1.66666666666667) - basisvalues[0]*0.555555555555556;
279
basisvalues[9] = basisvalues[5]*(0.05 + Y*1.75) - basisvalues[2]*0.7;
280
basisvalues[8] = basisvalues[4]*(0.54 + Y*2.1) - basisvalues[1]*0.56;
281
basisvalues[0] *= std::sqrt(0.5);
282
basisvalues[2] *= std::sqrt(1.0);
283
basisvalues[5] *= std::sqrt(1.5);
284
basisvalues[9] *= std::sqrt(2.0);
285
basisvalues[1] *= std::sqrt(3.0);
286
basisvalues[4] *= std::sqrt(4.5);
287
basisvalues[8] *= std::sqrt(6.0);
288
basisvalues[3] *= std::sqrt(7.5);
289
basisvalues[7] *= std::sqrt(10.0);
290
basisvalues[6] *= std::sqrt(14.0);
291
292
// Table(s) of coefficients
293
static const double coefficients0[10] = \
294
{0.106066017177982, 0.259807621135332, -0.15, 0.117369119465393, 0.0606091526731327, -0.0787335988751736, 0.0, 0.101644639076841, -0.131222664791956, 0.090913729009699};
295
296
// Compute value(s)
297
for (unsigned int r = 0; r < 10; r++)
298
{
299
*values += coefficients0[r]*basisvalues[r];
300
}// end loop over 'r'
301
break;
302
}
303
case 4:
304
{
305
306
// Array of basisvalues
307
double basisvalues[10] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
308
309
// Declare helper variables
310
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
311
double tmp1 = (1.0 - Y)/2.0;
312
double tmp2 = tmp1*tmp1;
313
314
// Compute basisvalues
315
basisvalues[0] = 1.0;
316
basisvalues[1] = tmp0;
317
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
318
basisvalues[6] = basisvalues[3]*1.66666666666667*tmp0 - basisvalues[1]*0.666666666666667*tmp2;
319
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
320
basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
321
basisvalues[7] = basisvalues[3]*(2.5 + 3.5*Y);
322
basisvalues[5] = basisvalues[2]*(0.111111111111111 + Y*1.66666666666667) - basisvalues[0]*0.555555555555556;
323
basisvalues[9] = basisvalues[5]*(0.05 + Y*1.75) - basisvalues[2]*0.7;
324
basisvalues[8] = basisvalues[4]*(0.54 + Y*2.1) - basisvalues[1]*0.56;
325
basisvalues[0] *= std::sqrt(0.5);
326
basisvalues[2] *= std::sqrt(1.0);
327
basisvalues[5] *= std::sqrt(1.5);
328
basisvalues[9] *= std::sqrt(2.0);
329
basisvalues[1] *= std::sqrt(3.0);
330
basisvalues[4] *= std::sqrt(4.5);
331
basisvalues[8] *= std::sqrt(6.0);
332
basisvalues[3] *= std::sqrt(7.5);
333
basisvalues[7] *= std::sqrt(10.0);
334
basisvalues[6] *= std::sqrt(14.0);
335
336
// Table(s) of coefficients
337
static const double coefficients0[10] = \
338
{0.106066017177982, 0.0, 0.3, 0.0, 0.151522881682832, 0.0262445329583912, 0.0, 0.0, 0.131222664791956, -0.136370593514548};
339
340
// Compute value(s)
341
for (unsigned int r = 0; r < 10; r++)
342
{
343
*values += coefficients0[r]*basisvalues[r];
344
}// end loop over 'r'
345
break;
346
}
347
case 5:
348
{
349
350
// Array of basisvalues
351
double basisvalues[10] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
352
353
// Declare helper variables
354
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
355
double tmp1 = (1.0 - Y)/2.0;
356
double tmp2 = tmp1*tmp1;
357
358
// Compute basisvalues
359
basisvalues[0] = 1.0;
360
basisvalues[1] = tmp0;
361
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
362
basisvalues[6] = basisvalues[3]*1.66666666666667*tmp0 - basisvalues[1]*0.666666666666667*tmp2;
363
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
364
basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
365
basisvalues[7] = basisvalues[3]*(2.5 + 3.5*Y);
366
basisvalues[5] = basisvalues[2]*(0.111111111111111 + Y*1.66666666666667) - basisvalues[0]*0.555555555555556;
367
basisvalues[9] = basisvalues[5]*(0.05 + Y*1.75) - basisvalues[2]*0.7;
368
basisvalues[8] = basisvalues[4]*(0.54 + Y*2.1) - basisvalues[1]*0.56;
369
basisvalues[0] *= std::sqrt(0.5);
370
basisvalues[2] *= std::sqrt(1.0);
371
basisvalues[5] *= std::sqrt(1.5);
372
basisvalues[9] *= std::sqrt(2.0);
373
basisvalues[1] *= std::sqrt(3.0);
374
basisvalues[4] *= std::sqrt(4.5);
375
basisvalues[8] *= std::sqrt(6.0);
376
basisvalues[3] *= std::sqrt(7.5);
377
basisvalues[7] *= std::sqrt(10.0);
378
basisvalues[6] *= std::sqrt(14.0);
379
380
// Table(s) of coefficients
381
static const double coefficients0[10] = \
382
{0.106066017177982, -0.259807621135332, -0.15, 0.117369119465393, -0.0606091526731326, -0.0787335988751736, 0.0, 0.101644639076841, 0.131222664791956, 0.090913729009699};
383
384
// Compute value(s)
385
for (unsigned int r = 0; r < 10; r++)
386
{
387
*values += coefficients0[r]*basisvalues[r];
388
}// end loop over 'r'
389
break;
390
}
391
case 6:
392
{
393
394
// Array of basisvalues
395
double basisvalues[10] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
396
397
// Declare helper variables
398
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
399
double tmp1 = (1.0 - Y)/2.0;
400
double tmp2 = tmp1*tmp1;
401
402
// Compute basisvalues
403
basisvalues[0] = 1.0;
404
basisvalues[1] = tmp0;
405
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
406
basisvalues[6] = basisvalues[3]*1.66666666666667*tmp0 - basisvalues[1]*0.666666666666667*tmp2;
407
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
408
basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
409
basisvalues[7] = basisvalues[3]*(2.5 + 3.5*Y);
410
basisvalues[5] = basisvalues[2]*(0.111111111111111 + Y*1.66666666666667) - basisvalues[0]*0.555555555555556;
411
basisvalues[9] = basisvalues[5]*(0.05 + Y*1.75) - basisvalues[2]*0.7;
412
basisvalues[8] = basisvalues[4]*(0.54 + Y*2.1) - basisvalues[1]*0.56;
413
basisvalues[0] *= std::sqrt(0.5);
414
basisvalues[2] *= std::sqrt(1.0);
415
basisvalues[5] *= std::sqrt(1.5);
416
basisvalues[9] *= std::sqrt(2.0);
417
basisvalues[1] *= std::sqrt(3.0);
418
basisvalues[4] *= std::sqrt(4.5);
419
basisvalues[8] *= std::sqrt(6.0);
420
basisvalues[3] *= std::sqrt(7.5);
421
basisvalues[7] *= std::sqrt(10.0);
422
basisvalues[6] *= std::sqrt(14.0);
423
424
// Table(s) of coefficients
425
static const double coefficients0[10] = \
426
{0.106066017177982, 0.0, 0.3, 0.0, -0.151522881682832, 0.0262445329583912, 0.0, 0.0, -0.131222664791956, -0.136370593514548};
427
428
// Compute value(s)
429
for (unsigned int r = 0; r < 10; r++)
430
{
431
*values += coefficients0[r]*basisvalues[r];
432
}// end loop over 'r'
433
break;
434
}
435
case 7:
436
{
437
438
// Array of basisvalues
439
double basisvalues[10] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
440
441
// Declare helper variables
442
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
443
double tmp1 = (1.0 - Y)/2.0;
444
double tmp2 = tmp1*tmp1;
445
446
// Compute basisvalues
447
basisvalues[0] = 1.0;
448
basisvalues[1] = tmp0;
449
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
450
basisvalues[6] = basisvalues[3]*1.66666666666667*tmp0 - basisvalues[1]*0.666666666666667*tmp2;
451
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
452
basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
453
basisvalues[7] = basisvalues[3]*(2.5 + 3.5*Y);
454
basisvalues[5] = basisvalues[2]*(0.111111111111111 + Y*1.66666666666667) - basisvalues[0]*0.555555555555556;
455
basisvalues[9] = basisvalues[5]*(0.05 + Y*1.75) - basisvalues[2]*0.7;
456
basisvalues[8] = basisvalues[4]*(0.54 + Y*2.1) - basisvalues[1]*0.56;
457
basisvalues[0] *= std::sqrt(0.5);
458
basisvalues[2] *= std::sqrt(1.0);
459
basisvalues[5] *= std::sqrt(1.5);
460
basisvalues[9] *= std::sqrt(2.0);
461
basisvalues[1] *= std::sqrt(3.0);
462
basisvalues[4] *= std::sqrt(4.5);
463
basisvalues[8] *= std::sqrt(6.0);
464
basisvalues[3] *= std::sqrt(7.5);
465
basisvalues[7] *= std::sqrt(10.0);
466
basisvalues[6] *= std::sqrt(14.0);
467
468
// Table(s) of coefficients
469
static const double coefficients0[10] = \
470
{0.106066017177982, -0.259807621135332, -0.15, -0.0782460796435952, 0.090913729009699, 0.0962299541807677, 0.180401338290886, 0.0508223195384204, -0.0131222664791956, -0.0227284322524247};
471
472
// Compute value(s)
473
for (unsigned int r = 0; r < 10; r++)
474
{
475
*values += coefficients0[r]*basisvalues[r];
476
}// end loop over 'r'
477
break;
478
}
479
case 8:
480
{
481
482
// Array of basisvalues
483
double basisvalues[10] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
484
485
// Declare helper variables
486
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
487
double tmp1 = (1.0 - Y)/2.0;
488
double tmp2 = tmp1*tmp1;
489
490
// Compute basisvalues
491
basisvalues[0] = 1.0;
492
basisvalues[1] = tmp0;
493
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
494
basisvalues[6] = basisvalues[3]*1.66666666666667*tmp0 - basisvalues[1]*0.666666666666667*tmp2;
495
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
496
basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
497
basisvalues[7] = basisvalues[3]*(2.5 + 3.5*Y);
498
basisvalues[5] = basisvalues[2]*(0.111111111111111 + Y*1.66666666666667) - basisvalues[0]*0.555555555555556;
499
basisvalues[9] = basisvalues[5]*(0.05 + Y*1.75) - basisvalues[2]*0.7;
500
basisvalues[8] = basisvalues[4]*(0.54 + Y*2.1) - basisvalues[1]*0.56;
501
basisvalues[0] *= std::sqrt(0.5);
502
basisvalues[2] *= std::sqrt(1.0);
503
basisvalues[5] *= std::sqrt(1.5);
504
basisvalues[9] *= std::sqrt(2.0);
505
basisvalues[1] *= std::sqrt(3.0);
506
basisvalues[4] *= std::sqrt(4.5);
507
basisvalues[8] *= std::sqrt(6.0);
508
basisvalues[3] *= std::sqrt(7.5);
509
basisvalues[7] *= std::sqrt(10.0);
510
basisvalues[6] *= std::sqrt(14.0);
511
512
// Table(s) of coefficients
513
static const double coefficients0[10] = \
514
{0.106066017177982, 0.259807621135332, -0.15, -0.0782460796435952, -0.090913729009699, 0.0962299541807678, -0.180401338290886, 0.0508223195384204, 0.0131222664791956, -0.0227284322524248};
515
516
// Compute value(s)
517
for (unsigned int r = 0; r < 10; r++)
518
{
519
*values += coefficients0[r]*basisvalues[r];
520
}// end loop over 'r'
521
break;
522
}
523
case 9:
524
{
525
526
// Array of basisvalues
527
double basisvalues[10] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
528
529
// Declare helper variables
530
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
531
double tmp1 = (1.0 - Y)/2.0;
532
double tmp2 = tmp1*tmp1;
533
534
// Compute basisvalues
535
basisvalues[0] = 1.0;
536
basisvalues[1] = tmp0;
537
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
538
basisvalues[6] = basisvalues[3]*1.66666666666667*tmp0 - basisvalues[1]*0.666666666666667*tmp2;
539
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
540
basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
541
basisvalues[7] = basisvalues[3]*(2.5 + 3.5*Y);
542
basisvalues[5] = basisvalues[2]*(0.111111111111111 + Y*1.66666666666667) - basisvalues[0]*0.555555555555556;
543
basisvalues[9] = basisvalues[5]*(0.05 + Y*1.75) - basisvalues[2]*0.7;
544
basisvalues[8] = basisvalues[4]*(0.54 + Y*2.1) - basisvalues[1]*0.56;
545
basisvalues[0] *= std::sqrt(0.5);
546
basisvalues[2] *= std::sqrt(1.0);
547
basisvalues[5] *= std::sqrt(1.5);
548
basisvalues[9] *= std::sqrt(2.0);
549
basisvalues[1] *= std::sqrt(3.0);
550
basisvalues[4] *= std::sqrt(4.5);
551
basisvalues[8] *= std::sqrt(6.0);
552
basisvalues[3] *= std::sqrt(7.5);
553
basisvalues[7] *= std::sqrt(10.0);
554
basisvalues[6] *= std::sqrt(14.0);
555
556
// Table(s) of coefficients
557
static const double coefficients0[10] = \
558
{0.636396103067893, 0.0, 0.0, -0.234738238930785, 0.0, -0.262445329583912, 0.0, -0.203289278153682, 0.0, 0.090913729009699};
559
560
// Compute value(s)
561
for (unsigned int r = 0; r < 10; r++)
562
{
563
*values += coefficients0[r]*basisvalues[r];
564
}// end loop over 'r'
565
break;
566
}
567
}
568
569
}
570
571
/// Evaluate all basis functions at given point x in cell
572
virtual void evaluate_basis_all(double* values,
573
const double* x,
574
const double* vertex_coordinates,
575
int cell_orientation) const
576
{
577
// Helper variable to hold values of a single dof.
578
double dof_values = 0.0;
579
580
// Loop dofs and call evaluate_basis
581
for (unsigned int r = 0; r < 10; r++)
582
{
583
evaluate_basis(r, &dof_values, x, vertex_coordinates, cell_orientation);
584
values[r] = dof_values;
585
}// end loop over 'r'
586
}
587
588
/// Evaluate order n derivatives of basis function i at given point x in cell
589
virtual void evaluate_basis_derivatives(std::size_t i,
590
std::size_t n,
591
double* values,
592
const double* x,
593
const double* vertex_coordinates,
594
int cell_orientation) const
595
{
596
// Compute Jacobian
597
double J[4];
598
compute_jacobian_triangle_2d(J, vertex_coordinates);
599
600
// Compute Jacobian inverse and determinant
601
double K[4];
602
double detJ;
603
compute_jacobian_inverse_triangle_2d(K, detJ, J);
604
605
606
// Compute constants
607
const double C0 = vertex_coordinates[2] + vertex_coordinates[4];
608
const double C1 = vertex_coordinates[3] + vertex_coordinates[5];
609
610
// Get coordinates and map to the reference (FIAT) element
611
double X = (J[1]*(C1 - 2.0*x[1]) + J[3]*(2.0*x[0] - C0)) / detJ;
612
double Y = (J[0]*(2.0*x[1] - C1) + J[2]*(C0 - 2.0*x[0])) / detJ;
613
614
// Compute number of derivatives.
615
unsigned int num_derivatives = 1;
616
for (unsigned int r = 0; r < n; r++)
617
{
618
num_derivatives *= 2;
619
}// end loop over 'r'
620
621
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
622
unsigned int **combinations = new unsigned int *[num_derivatives];
623
for (unsigned int row = 0; row < num_derivatives; row++)
624
{
625
combinations[row] = new unsigned int [n];
626
for (unsigned int col = 0; col < n; col++)
627
combinations[row][col] = 0;
628
}
629
630
// Generate combinations of derivatives
631
for (unsigned int row = 1; row < num_derivatives; row++)
632
{
633
for (unsigned int num = 0; num < row; num++)
634
{
635
for (unsigned int col = n-1; col+1 > 0; col--)
636
{
637
if (combinations[row][col] + 1 > 1)
638
combinations[row][col] = 0;
639
else
640
{
641
combinations[row][col] += 1;
642
break;
643
}
644
}
645
}
646
}
647
648
// Compute inverse of Jacobian
649
const double Jinv[2][2] = {{K[0], K[1]}, {K[2], K[3]}};
650
651
// Declare transformation matrix
652
// Declare pointer to two dimensional array and initialise
653
double **transform = new double *[num_derivatives];
654
655
for (unsigned int j = 0; j < num_derivatives; j++)
656
{
657
transform[j] = new double [num_derivatives];
658
for (unsigned int k = 0; k < num_derivatives; k++)
659
transform[j][k] = 1;
660
}
661
662
// Construct transformation matrix
663
for (unsigned int row = 0; row < num_derivatives; row++)
664
{
665
for (unsigned int col = 0; col < num_derivatives; col++)
666
{
667
for (unsigned int k = 0; k < n; k++)
668
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
669
}
670
}
671
672
// Reset values. Assuming that values is always an array.
673
for (unsigned int r = 0; r < num_derivatives; r++)
674
{
675
values[r] = 0.0;
676
}// end loop over 'r'
677
678
switch (i)
679
{
680
case 0:
681
{
682
683
// Array of basisvalues
684
double basisvalues[10] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
685
686
// Declare helper variables
687
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
688
double tmp1 = (1.0 - Y)/2.0;
689
double tmp2 = tmp1*tmp1;
690
691
// Compute basisvalues
692
basisvalues[0] = 1.0;
693
basisvalues[1] = tmp0;
694
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
695
basisvalues[6] = basisvalues[3]*1.66666666666667*tmp0 - basisvalues[1]*0.666666666666667*tmp2;
696
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
697
basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
698
basisvalues[7] = basisvalues[3]*(2.5 + 3.5*Y);
699
basisvalues[5] = basisvalues[2]*(0.111111111111111 + Y*1.66666666666667) - basisvalues[0]*0.555555555555556;
700
basisvalues[9] = basisvalues[5]*(0.05 + Y*1.75) - basisvalues[2]*0.7;
701
basisvalues[8] = basisvalues[4]*(0.54 + Y*2.1) - basisvalues[1]*0.56;
702
basisvalues[0] *= std::sqrt(0.5);
703
basisvalues[2] *= std::sqrt(1.0);
704
basisvalues[5] *= std::sqrt(1.5);
705
basisvalues[9] *= std::sqrt(2.0);
706
basisvalues[1] *= std::sqrt(3.0);
707
basisvalues[4] *= std::sqrt(4.5);
708
basisvalues[8] *= std::sqrt(6.0);
709
basisvalues[3] *= std::sqrt(7.5);
710
basisvalues[7] *= std::sqrt(10.0);
711
basisvalues[6] *= std::sqrt(14.0);
712
713
// Table(s) of coefficients
714
static const double coefficients0[10] = \
715
{0.0471404520791032, -0.0288675134594813, -0.0166666666666667, 0.0782460796435951, 0.0606091526731326, 0.0349927106111883, -0.0601337794302955, -0.0508223195384204, -0.0393667994375868, -0.0227284322524248};
716
717
// Tables of derivatives of the polynomial base (transpose).
718
static const double dmats0[10][10] = \
719
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
720
{4.89897948556636, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
721
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
722
{0.0, 9.48683298050514, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
723
{4, 0.0, 7.07106781186548, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
724
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
725
{5.29150262212918, 0.0, -2.99332590941915, 13.6626010212795, 0.0, 0.611010092660778, 0.0, 0.0, 0.0, 0.0},
726
{0.0, 4.38178046004133, 0.0, 0.0, 12.5219806739988, 0.0, 0.0, 0.0, 0.0, 0.0},
727
{3.46410161513775, 0.0, 7.83836717690617, 0.0, 0.0, 8.4, 0.0, 0.0, 0.0, 0.0},
728
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
729
730
static const double dmats1[10][10] = \
731
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
732
{2.44948974278318, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
733
{4.24264068711928, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
734
{2.58198889747161, 4.74341649025257, -0.912870929175277, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
735
{2.0, 6.12372435695795, 3.53553390593274, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
736
{-2.30940107675851, 0.0, 8.16496580927726, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
737
{2.64575131106459, 5.18459255872629, -1.49666295470958, 6.83130051063973, -1.05830052442584, 0.30550504633039, 0.0, 0.0, 0.0, 0.0},
738
{2.23606797749979, 2.19089023002067, 2.52982212813471, 8.08290376865476, 6.26099033699941, -1.80739222823013, 0.0, 0.0, 0.0, 0.0},
739
{1.73205080756887, -5.09116882454314, 3.91918358845309, 0.0, 9.69948452238572, 4.2, 0.0, 0.0, 0.0, 0.0},
740
{5.0, 0.0, -2.82842712474619, 0.0, 0.0, 12.1243556529821, 0.0, 0.0, 0.0, 0.0}};
741
742
// Compute reference derivatives.
743
// Declare pointer to array of derivatives on FIAT element.
744
double *derivatives = new double[num_derivatives];
745
for (unsigned int r = 0; r < num_derivatives; r++)
746
{
747
derivatives[r] = 0.0;
748
}// end loop over 'r'
749
750
// Declare derivative matrix (of polynomial basis).
751
double dmats[10][10] = \
752
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
753
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
754
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
755
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
756
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
757
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
758
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
759
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
760
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
761
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
762
763
// Declare (auxiliary) derivative matrix (of polynomial basis).
764
double dmats_old[10][10] = \
765
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
766
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
767
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
768
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
769
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
770
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
771
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
772
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
773
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
774
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
775
776
// Loop possible derivatives.
777
for (unsigned int r = 0; r < num_derivatives; r++)
778
{
779
// Resetting dmats values to compute next derivative.
780
for (unsigned int t = 0; t < 10; t++)
781
{
782
for (unsigned int u = 0; u < 10; u++)
783
{
784
dmats[t][u] = 0.0;
785
if (t == u)
786
{
787
dmats[t][u] = 1.0;
788
}
789
790
}// end loop over 'u'
791
}// end loop over 't'
792
793
// Looping derivative order to generate dmats.
794
for (unsigned int s = 0; s < n; s++)
795
{
796
// Updating dmats_old with new values and resetting dmats.
797
for (unsigned int t = 0; t < 10; t++)
798
{
799
for (unsigned int u = 0; u < 10; u++)
800
{
801
dmats_old[t][u] = dmats[t][u];
802
dmats[t][u] = 0.0;
803
}// end loop over 'u'
804
}// end loop over 't'
805
806
// Update dmats using an inner product.
807
if (combinations[r][s] == 0)
808
{
809
for (unsigned int t = 0; t < 10; t++)
810
{
811
for (unsigned int u = 0; u < 10; u++)
812
{
813
for (unsigned int tu = 0; tu < 10; tu++)
814
{
815
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
816
}// end loop over 'tu'
817
}// end loop over 'u'
818
}// end loop over 't'
819
}
820
821
if (combinations[r][s] == 1)
822
{
823
for (unsigned int t = 0; t < 10; t++)
824
{
825
for (unsigned int u = 0; u < 10; u++)
826
{
827
for (unsigned int tu = 0; tu < 10; tu++)
828
{
829
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
830
}// end loop over 'tu'
831
}// end loop over 'u'
832
}// end loop over 't'
833
}
834
835
}// end loop over 's'
836
for (unsigned int s = 0; s < 10; s++)
837
{
838
for (unsigned int t = 0; t < 10; t++)
839
{
840
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
841
}// end loop over 't'
842
}// end loop over 's'
843
}// end loop over 'r'
844
845
// Transform derivatives back to physical element
846
for (unsigned int r = 0; r < num_derivatives; r++)
847
{
848
for (unsigned int s = 0; s < num_derivatives; s++)
849
{
850
values[r] += transform[r][s]*derivatives[s];
851
}// end loop over 's'
852
}// end loop over 'r'
853
854
// Delete pointer to array of derivatives on FIAT element
855
delete [] derivatives;
856
857
// Delete pointer to array of combinations of derivatives and transform
858
for (unsigned int r = 0; r < num_derivatives; r++)
859
{
860
delete [] combinations[r];
861
}// end loop over 'r'
862
delete [] combinations;
863
for (unsigned int r = 0; r < num_derivatives; r++)
864
{
865
delete [] transform[r];
866
}// end loop over 'r'
867
delete [] transform;
868
break;
869
}
870
case 1:
871
{
872
873
// Array of basisvalues
874
double basisvalues[10] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
875
876
// Declare helper variables
877
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
878
double tmp1 = (1.0 - Y)/2.0;
879
double tmp2 = tmp1*tmp1;
880
881
// Compute basisvalues
882
basisvalues[0] = 1.0;
883
basisvalues[1] = tmp0;
884
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
885
basisvalues[6] = basisvalues[3]*1.66666666666667*tmp0 - basisvalues[1]*0.666666666666667*tmp2;
886
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
887
basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
888
basisvalues[7] = basisvalues[3]*(2.5 + 3.5*Y);
889
basisvalues[5] = basisvalues[2]*(0.111111111111111 + Y*1.66666666666667) - basisvalues[0]*0.555555555555556;
890
basisvalues[9] = basisvalues[5]*(0.05 + Y*1.75) - basisvalues[2]*0.7;
891
basisvalues[8] = basisvalues[4]*(0.54 + Y*2.1) - basisvalues[1]*0.56;
892
basisvalues[0] *= std::sqrt(0.5);
893
basisvalues[2] *= std::sqrt(1.0);
894
basisvalues[5] *= std::sqrt(1.5);
895
basisvalues[9] *= std::sqrt(2.0);
896
basisvalues[1] *= std::sqrt(3.0);
897
basisvalues[4] *= std::sqrt(4.5);
898
basisvalues[8] *= std::sqrt(6.0);
899
basisvalues[3] *= std::sqrt(7.5);
900
basisvalues[7] *= std::sqrt(10.0);
901
basisvalues[6] *= std::sqrt(14.0);
902
903
// Table(s) of coefficients
904
static const double coefficients0[10] = \
905
{0.0471404520791032, 0.0288675134594813, -0.0166666666666666, 0.0782460796435952, -0.0606091526731327, 0.0349927106111883, 0.0601337794302955, -0.0508223195384204, 0.0393667994375868, -0.0227284322524248};
906
907
// Tables of derivatives of the polynomial base (transpose).
908
static const double dmats0[10][10] = \
909
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
910
{4.89897948556636, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
911
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
912
{0.0, 9.48683298050514, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
913
{4, 0.0, 7.07106781186548, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
914
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
915
{5.29150262212918, 0.0, -2.99332590941915, 13.6626010212795, 0.0, 0.611010092660778, 0.0, 0.0, 0.0, 0.0},
916
{0.0, 4.38178046004133, 0.0, 0.0, 12.5219806739988, 0.0, 0.0, 0.0, 0.0, 0.0},
917
{3.46410161513775, 0.0, 7.83836717690617, 0.0, 0.0, 8.4, 0.0, 0.0, 0.0, 0.0},
918
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
919
920
static const double dmats1[10][10] = \
921
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
922
{2.44948974278318, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
923
{4.24264068711928, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
924
{2.58198889747161, 4.74341649025257, -0.912870929175277, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
925
{2.0, 6.12372435695795, 3.53553390593274, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
926
{-2.30940107675851, 0.0, 8.16496580927726, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
927
{2.64575131106459, 5.18459255872629, -1.49666295470958, 6.83130051063973, -1.05830052442584, 0.30550504633039, 0.0, 0.0, 0.0, 0.0},
928
{2.23606797749979, 2.19089023002067, 2.52982212813471, 8.08290376865476, 6.26099033699941, -1.80739222823013, 0.0, 0.0, 0.0, 0.0},
929
{1.73205080756887, -5.09116882454314, 3.91918358845309, 0.0, 9.69948452238572, 4.2, 0.0, 0.0, 0.0, 0.0},
930
{5.0, 0.0, -2.82842712474619, 0.0, 0.0, 12.1243556529821, 0.0, 0.0, 0.0, 0.0}};
931
932
// Compute reference derivatives.
933
// Declare pointer to array of derivatives on FIAT element.
934
double *derivatives = new double[num_derivatives];
935
for (unsigned int r = 0; r < num_derivatives; r++)
936
{
937
derivatives[r] = 0.0;
938
}// end loop over 'r'
939
940
// Declare derivative matrix (of polynomial basis).
941
double dmats[10][10] = \
942
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
943
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
944
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
945
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
946
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
947
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
948
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
949
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
950
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
951
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
952
953
// Declare (auxiliary) derivative matrix (of polynomial basis).
954
double dmats_old[10][10] = \
955
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
956
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
957
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
958
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
959
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
960
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
961
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
962
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
963
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
964
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
965
966
// Loop possible derivatives.
967
for (unsigned int r = 0; r < num_derivatives; r++)
968
{
969
// Resetting dmats values to compute next derivative.
970
for (unsigned int t = 0; t < 10; t++)
971
{
972
for (unsigned int u = 0; u < 10; u++)
973
{
974
dmats[t][u] = 0.0;
975
if (t == u)
976
{
977
dmats[t][u] = 1.0;
978
}
979
980
}// end loop over 'u'
981
}// end loop over 't'
982
983
// Looping derivative order to generate dmats.
984
for (unsigned int s = 0; s < n; s++)
985
{
986
// Updating dmats_old with new values and resetting dmats.
987
for (unsigned int t = 0; t < 10; t++)
988
{
989
for (unsigned int u = 0; u < 10; u++)
990
{
991
dmats_old[t][u] = dmats[t][u];
992
dmats[t][u] = 0.0;
993
}// end loop over 'u'
994
}// end loop over 't'
995
996
// Update dmats using an inner product.
997
if (combinations[r][s] == 0)
998
{
999
for (unsigned int t = 0; t < 10; t++)
1000
{
1001
for (unsigned int u = 0; u < 10; u++)
1002
{
1003
for (unsigned int tu = 0; tu < 10; tu++)
1004
{
1005
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1006
}// end loop over 'tu'
1007
}// end loop over 'u'
1008
}// end loop over 't'
1009
}
1010
1011
if (combinations[r][s] == 1)
1012
{
1013
for (unsigned int t = 0; t < 10; t++)
1014
{
1015
for (unsigned int u = 0; u < 10; u++)
1016
{
1017
for (unsigned int tu = 0; tu < 10; tu++)
1018
{
1019
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1020
}// end loop over 'tu'
1021
}// end loop over 'u'
1022
}// end loop over 't'
1023
}
1024
1025
}// end loop over 's'
1026
for (unsigned int s = 0; s < 10; s++)
1027
{
1028
for (unsigned int t = 0; t < 10; t++)
1029
{
1030
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1031
}// end loop over 't'
1032
}// end loop over 's'
1033
}// end loop over 'r'
1034
1035
// Transform derivatives back to physical element
1036
for (unsigned int r = 0; r < num_derivatives; r++)
1037
{
1038
for (unsigned int s = 0; s < num_derivatives; s++)
1039
{
1040
values[r] += transform[r][s]*derivatives[s];
1041
}// end loop over 's'
1042
}// end loop over 'r'
1043
1044
// Delete pointer to array of derivatives on FIAT element
1045
delete [] derivatives;
1046
1047
// Delete pointer to array of combinations of derivatives and transform
1048
for (unsigned int r = 0; r < num_derivatives; r++)
1049
{
1050
delete [] combinations[r];
1051
}// end loop over 'r'
1052
delete [] combinations;
1053
for (unsigned int r = 0; r < num_derivatives; r++)
1054
{
1055
delete [] transform[r];
1056
}// end loop over 'r'
1057
delete [] transform;
1058
break;
1059
}
1060
case 2:
1061
{
1062
1063
// Array of basisvalues
1064
double basisvalues[10] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
1065
1066
// Declare helper variables
1067
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1068
double tmp1 = (1.0 - Y)/2.0;
1069
double tmp2 = tmp1*tmp1;
1070
1071
// Compute basisvalues
1072
basisvalues[0] = 1.0;
1073
basisvalues[1] = tmp0;
1074
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
1075
basisvalues[6] = basisvalues[3]*1.66666666666667*tmp0 - basisvalues[1]*0.666666666666667*tmp2;
1076
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1077
basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
1078
basisvalues[7] = basisvalues[3]*(2.5 + 3.5*Y);
1079
basisvalues[5] = basisvalues[2]*(0.111111111111111 + Y*1.66666666666667) - basisvalues[0]*0.555555555555556;
1080
basisvalues[9] = basisvalues[5]*(0.05 + Y*1.75) - basisvalues[2]*0.7;
1081
basisvalues[8] = basisvalues[4]*(0.54 + Y*2.1) - basisvalues[1]*0.56;
1082
basisvalues[0] *= std::sqrt(0.5);
1083
basisvalues[2] *= std::sqrt(1.0);
1084
basisvalues[5] *= std::sqrt(1.5);
1085
basisvalues[9] *= std::sqrt(2.0);
1086
basisvalues[1] *= std::sqrt(3.0);
1087
basisvalues[4] *= std::sqrt(4.5);
1088
basisvalues[8] *= std::sqrt(6.0);
1089
basisvalues[3] *= std::sqrt(7.5);
1090
basisvalues[7] *= std::sqrt(10.0);
1091
basisvalues[6] *= std::sqrt(14.0);
1092
1093
// Table(s) of coefficients
1094
static const double coefficients0[10] = \
1095
{0.0471404520791032, 0.0, 0.0333333333333334, 0.0, 0.0, 0.104978131833565, 0.0, 0.0, 0.0, 0.0909137290096989};
1096
1097
// Tables of derivatives of the polynomial base (transpose).
1098
static const double dmats0[10][10] = \
1099
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1100
{4.89897948556636, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1101
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1102
{0.0, 9.48683298050514, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1103
{4, 0.0, 7.07106781186548, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1104
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1105
{5.29150262212918, 0.0, -2.99332590941915, 13.6626010212795, 0.0, 0.611010092660778, 0.0, 0.0, 0.0, 0.0},
1106
{0.0, 4.38178046004133, 0.0, 0.0, 12.5219806739988, 0.0, 0.0, 0.0, 0.0, 0.0},
1107
{3.46410161513775, 0.0, 7.83836717690617, 0.0, 0.0, 8.4, 0.0, 0.0, 0.0, 0.0},
1108
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
1109
1110
static const double dmats1[10][10] = \
1111
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1112
{2.44948974278318, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1113
{4.24264068711928, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1114
{2.58198889747161, 4.74341649025257, -0.912870929175277, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1115
{2.0, 6.12372435695795, 3.53553390593274, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1116
{-2.30940107675851, 0.0, 8.16496580927726, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1117
{2.64575131106459, 5.18459255872629, -1.49666295470958, 6.83130051063973, -1.05830052442584, 0.30550504633039, 0.0, 0.0, 0.0, 0.0},
1118
{2.23606797749979, 2.19089023002067, 2.52982212813471, 8.08290376865476, 6.26099033699941, -1.80739222823013, 0.0, 0.0, 0.0, 0.0},
1119
{1.73205080756887, -5.09116882454314, 3.91918358845309, 0.0, 9.69948452238572, 4.2, 0.0, 0.0, 0.0, 0.0},
1120
{5.0, 0.0, -2.82842712474619, 0.0, 0.0, 12.1243556529821, 0.0, 0.0, 0.0, 0.0}};
1121
1122
// Compute reference derivatives.
1123
// Declare pointer to array of derivatives on FIAT element.
1124
double *derivatives = new double[num_derivatives];
1125
for (unsigned int r = 0; r < num_derivatives; r++)
1126
{
1127
derivatives[r] = 0.0;
1128
}// end loop over 'r'
1129
1130
// Declare derivative matrix (of polynomial basis).
1131
double dmats[10][10] = \
1132
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1133
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1134
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1135
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1136
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1137
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
1138
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
1139
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
1140
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
1141
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
1142
1143
// Declare (auxiliary) derivative matrix (of polynomial basis).
1144
double dmats_old[10][10] = \
1145
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1146
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1147
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1148
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1149
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1150
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
1151
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
1152
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
1153
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
1154
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
1155
1156
// Loop possible derivatives.
1157
for (unsigned int r = 0; r < num_derivatives; r++)
1158
{
1159
// Resetting dmats values to compute next derivative.
1160
for (unsigned int t = 0; t < 10; t++)
1161
{
1162
for (unsigned int u = 0; u < 10; u++)
1163
{
1164
dmats[t][u] = 0.0;
1165
if (t == u)
1166
{
1167
dmats[t][u] = 1.0;
1168
}
1169
1170
}// end loop over 'u'
1171
}// end loop over 't'
1172
1173
// Looping derivative order to generate dmats.
1174
for (unsigned int s = 0; s < n; s++)
1175
{
1176
// Updating dmats_old with new values and resetting dmats.
1177
for (unsigned int t = 0; t < 10; t++)
1178
{
1179
for (unsigned int u = 0; u < 10; u++)
1180
{
1181
dmats_old[t][u] = dmats[t][u];
1182
dmats[t][u] = 0.0;
1183
}// end loop over 'u'
1184
}// end loop over 't'
1185
1186
// Update dmats using an inner product.
1187
if (combinations[r][s] == 0)
1188
{
1189
for (unsigned int t = 0; t < 10; t++)
1190
{
1191
for (unsigned int u = 0; u < 10; u++)
1192
{
1193
for (unsigned int tu = 0; tu < 10; tu++)
1194
{
1195
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1196
}// end loop over 'tu'
1197
}// end loop over 'u'
1198
}// end loop over 't'
1199
}
1200
1201
if (combinations[r][s] == 1)
1202
{
1203
for (unsigned int t = 0; t < 10; t++)
1204
{
1205
for (unsigned int u = 0; u < 10; u++)
1206
{
1207
for (unsigned int tu = 0; tu < 10; tu++)
1208
{
1209
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1210
}// end loop over 'tu'
1211
}// end loop over 'u'
1212
}// end loop over 't'
1213
}
1214
1215
}// end loop over 's'
1216
for (unsigned int s = 0; s < 10; s++)
1217
{
1218
for (unsigned int t = 0; t < 10; t++)
1219
{
1220
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1221
}// end loop over 't'
1222
}// end loop over 's'
1223
}// end loop over 'r'
1224
1225
// Transform derivatives back to physical element
1226
for (unsigned int r = 0; r < num_derivatives; r++)
1227
{
1228
for (unsigned int s = 0; s < num_derivatives; s++)
1229
{
1230
values[r] += transform[r][s]*derivatives[s];
1231
}// end loop over 's'
1232
}// end loop over 'r'
1233
1234
// Delete pointer to array of derivatives on FIAT element
1235
delete [] derivatives;
1236
1237
// Delete pointer to array of combinations of derivatives and transform
1238
for (unsigned int r = 0; r < num_derivatives; r++)
1239
{
1240
delete [] combinations[r];
1241
}// end loop over 'r'
1242
delete [] combinations;
1243
for (unsigned int r = 0; r < num_derivatives; r++)
1244
{
1245
delete [] transform[r];
1246
}// end loop over 'r'
1247
delete [] transform;
1248
break;
1249
}
1250
case 3:
1251
{
1252
1253
// Array of basisvalues
1254
double basisvalues[10] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
1255
1256
// Declare helper variables
1257
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1258
double tmp1 = (1.0 - Y)/2.0;
1259
double tmp2 = tmp1*tmp1;
1260
1261
// Compute basisvalues
1262
basisvalues[0] = 1.0;
1263
basisvalues[1] = tmp0;
1264
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
1265
basisvalues[6] = basisvalues[3]*1.66666666666667*tmp0 - basisvalues[1]*0.666666666666667*tmp2;
1266
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1267
basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
1268
basisvalues[7] = basisvalues[3]*(2.5 + 3.5*Y);
1269
basisvalues[5] = basisvalues[2]*(0.111111111111111 + Y*1.66666666666667) - basisvalues[0]*0.555555555555556;
1270
basisvalues[9] = basisvalues[5]*(0.05 + Y*1.75) - basisvalues[2]*0.7;
1271
basisvalues[8] = basisvalues[4]*(0.54 + Y*2.1) - basisvalues[1]*0.56;
1272
basisvalues[0] *= std::sqrt(0.5);
1273
basisvalues[2] *= std::sqrt(1.0);
1274
basisvalues[5] *= std::sqrt(1.5);
1275
basisvalues[9] *= std::sqrt(2.0);
1276
basisvalues[1] *= std::sqrt(3.0);
1277
basisvalues[4] *= std::sqrt(4.5);
1278
basisvalues[8] *= std::sqrt(6.0);
1279
basisvalues[3] *= std::sqrt(7.5);
1280
basisvalues[7] *= std::sqrt(10.0);
1281
basisvalues[6] *= std::sqrt(14.0);
1282
1283
// Table(s) of coefficients
1284
static const double coefficients0[10] = \
1285
{0.106066017177982, 0.259807621135332, -0.15, 0.117369119465393, 0.0606091526731327, -0.0787335988751736, 0.0, 0.101644639076841, -0.131222664791956, 0.090913729009699};
1286
1287
// Tables of derivatives of the polynomial base (transpose).
1288
static const double dmats0[10][10] = \
1289
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1290
{4.89897948556636, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1291
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1292
{0.0, 9.48683298050514, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1293
{4, 0.0, 7.07106781186548, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1294
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1295
{5.29150262212918, 0.0, -2.99332590941915, 13.6626010212795, 0.0, 0.611010092660778, 0.0, 0.0, 0.0, 0.0},
1296
{0.0, 4.38178046004133, 0.0, 0.0, 12.5219806739988, 0.0, 0.0, 0.0, 0.0, 0.0},
1297
{3.46410161513775, 0.0, 7.83836717690617, 0.0, 0.0, 8.4, 0.0, 0.0, 0.0, 0.0},
1298
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
1299
1300
static const double dmats1[10][10] = \
1301
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1302
{2.44948974278318, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1303
{4.24264068711928, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1304
{2.58198889747161, 4.74341649025257, -0.912870929175277, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1305
{2.0, 6.12372435695795, 3.53553390593274, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1306
{-2.30940107675851, 0.0, 8.16496580927726, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1307
{2.64575131106459, 5.18459255872629, -1.49666295470958, 6.83130051063973, -1.05830052442584, 0.30550504633039, 0.0, 0.0, 0.0, 0.0},
1308
{2.23606797749979, 2.19089023002067, 2.52982212813471, 8.08290376865476, 6.26099033699941, -1.80739222823013, 0.0, 0.0, 0.0, 0.0},
1309
{1.73205080756887, -5.09116882454314, 3.91918358845309, 0.0, 9.69948452238572, 4.2, 0.0, 0.0, 0.0, 0.0},
1310
{5.0, 0.0, -2.82842712474619, 0.0, 0.0, 12.1243556529821, 0.0, 0.0, 0.0, 0.0}};
1311
1312
// Compute reference derivatives.
1313
// Declare pointer to array of derivatives on FIAT element.
1314
double *derivatives = new double[num_derivatives];
1315
for (unsigned int r = 0; r < num_derivatives; r++)
1316
{
1317
derivatives[r] = 0.0;
1318
}// end loop over 'r'
1319
1320
// Declare derivative matrix (of polynomial basis).
1321
double dmats[10][10] = \
1322
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1323
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1324
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1325
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1326
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1327
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
1328
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
1329
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
1330
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
1331
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
1332
1333
// Declare (auxiliary) derivative matrix (of polynomial basis).
1334
double dmats_old[10][10] = \
1335
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1336
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1337
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1338
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1339
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1340
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
1341
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
1342
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
1343
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
1344
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
1345
1346
// Loop possible derivatives.
1347
for (unsigned int r = 0; r < num_derivatives; r++)
1348
{
1349
// Resetting dmats values to compute next derivative.
1350
for (unsigned int t = 0; t < 10; t++)
1351
{
1352
for (unsigned int u = 0; u < 10; u++)
1353
{
1354
dmats[t][u] = 0.0;
1355
if (t == u)
1356
{
1357
dmats[t][u] = 1.0;
1358
}
1359
1360
}// end loop over 'u'
1361
}// end loop over 't'
1362
1363
// Looping derivative order to generate dmats.
1364
for (unsigned int s = 0; s < n; s++)
1365
{
1366
// Updating dmats_old with new values and resetting dmats.
1367
for (unsigned int t = 0; t < 10; t++)
1368
{
1369
for (unsigned int u = 0; u < 10; u++)
1370
{
1371
dmats_old[t][u] = dmats[t][u];
1372
dmats[t][u] = 0.0;
1373
}// end loop over 'u'
1374
}// end loop over 't'
1375
1376
// Update dmats using an inner product.
1377
if (combinations[r][s] == 0)
1378
{
1379
for (unsigned int t = 0; t < 10; t++)
1380
{
1381
for (unsigned int u = 0; u < 10; u++)
1382
{
1383
for (unsigned int tu = 0; tu < 10; tu++)
1384
{
1385
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1386
}// end loop over 'tu'
1387
}// end loop over 'u'
1388
}// end loop over 't'
1389
}
1390
1391
if (combinations[r][s] == 1)
1392
{
1393
for (unsigned int t = 0; t < 10; t++)
1394
{
1395
for (unsigned int u = 0; u < 10; u++)
1396
{
1397
for (unsigned int tu = 0; tu < 10; tu++)
1398
{
1399
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1400
}// end loop over 'tu'
1401
}// end loop over 'u'
1402
}// end loop over 't'
1403
}
1404
1405
}// end loop over 's'
1406
for (unsigned int s = 0; s < 10; s++)
1407
{
1408
for (unsigned int t = 0; t < 10; t++)
1409
{
1410
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1411
}// end loop over 't'
1412
}// end loop over 's'
1413
}// end loop over 'r'
1414
1415
// Transform derivatives back to physical element
1416
for (unsigned int r = 0; r < num_derivatives; r++)
1417
{
1418
for (unsigned int s = 0; s < num_derivatives; s++)
1419
{
1420
values[r] += transform[r][s]*derivatives[s];
1421
}// end loop over 's'
1422
}// end loop over 'r'
1423
1424
// Delete pointer to array of derivatives on FIAT element
1425
delete [] derivatives;
1426
1427
// Delete pointer to array of combinations of derivatives and transform
1428
for (unsigned int r = 0; r < num_derivatives; r++)
1429
{
1430
delete [] combinations[r];
1431
}// end loop over 'r'
1432
delete [] combinations;
1433
for (unsigned int r = 0; r < num_derivatives; r++)
1434
{
1435
delete [] transform[r];
1436
}// end loop over 'r'
1437
delete [] transform;
1438
break;
1439
}
1440
case 4:
1441
{
1442
1443
// Array of basisvalues
1444
double basisvalues[10] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
1445
1446
// Declare helper variables
1447
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1448
double tmp1 = (1.0 - Y)/2.0;
1449
double tmp2 = tmp1*tmp1;
1450
1451
// Compute basisvalues
1452
basisvalues[0] = 1.0;
1453
basisvalues[1] = tmp0;
1454
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
1455
basisvalues[6] = basisvalues[3]*1.66666666666667*tmp0 - basisvalues[1]*0.666666666666667*tmp2;
1456
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1457
basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
1458
basisvalues[7] = basisvalues[3]*(2.5 + 3.5*Y);
1459
basisvalues[5] = basisvalues[2]*(0.111111111111111 + Y*1.66666666666667) - basisvalues[0]*0.555555555555556;
1460
basisvalues[9] = basisvalues[5]*(0.05 + Y*1.75) - basisvalues[2]*0.7;
1461
basisvalues[8] = basisvalues[4]*(0.54 + Y*2.1) - basisvalues[1]*0.56;
1462
basisvalues[0] *= std::sqrt(0.5);
1463
basisvalues[2] *= std::sqrt(1.0);
1464
basisvalues[5] *= std::sqrt(1.5);
1465
basisvalues[9] *= std::sqrt(2.0);
1466
basisvalues[1] *= std::sqrt(3.0);
1467
basisvalues[4] *= std::sqrt(4.5);
1468
basisvalues[8] *= std::sqrt(6.0);
1469
basisvalues[3] *= std::sqrt(7.5);
1470
basisvalues[7] *= std::sqrt(10.0);
1471
basisvalues[6] *= std::sqrt(14.0);
1472
1473
// Table(s) of coefficients
1474
static const double coefficients0[10] = \
1475
{0.106066017177982, 0.0, 0.3, 0.0, 0.151522881682832, 0.0262445329583912, 0.0, 0.0, 0.131222664791956, -0.136370593514548};
1476
1477
// Tables of derivatives of the polynomial base (transpose).
1478
static const double dmats0[10][10] = \
1479
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1480
{4.89897948556636, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1481
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1482
{0.0, 9.48683298050514, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1483
{4, 0.0, 7.07106781186548, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1484
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1485
{5.29150262212918, 0.0, -2.99332590941915, 13.6626010212795, 0.0, 0.611010092660778, 0.0, 0.0, 0.0, 0.0},
1486
{0.0, 4.38178046004133, 0.0, 0.0, 12.5219806739988, 0.0, 0.0, 0.0, 0.0, 0.0},
1487
{3.46410161513775, 0.0, 7.83836717690617, 0.0, 0.0, 8.4, 0.0, 0.0, 0.0, 0.0},
1488
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
1489
1490
static const double dmats1[10][10] = \
1491
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1492
{2.44948974278318, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1493
{4.24264068711928, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1494
{2.58198889747161, 4.74341649025257, -0.912870929175277, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1495
{2.0, 6.12372435695795, 3.53553390593274, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1496
{-2.30940107675851, 0.0, 8.16496580927726, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1497
{2.64575131106459, 5.18459255872629, -1.49666295470958, 6.83130051063973, -1.05830052442584, 0.30550504633039, 0.0, 0.0, 0.0, 0.0},
1498
{2.23606797749979, 2.19089023002067, 2.52982212813471, 8.08290376865476, 6.26099033699941, -1.80739222823013, 0.0, 0.0, 0.0, 0.0},
1499
{1.73205080756887, -5.09116882454314, 3.91918358845309, 0.0, 9.69948452238572, 4.2, 0.0, 0.0, 0.0, 0.0},
1500
{5.0, 0.0, -2.82842712474619, 0.0, 0.0, 12.1243556529821, 0.0, 0.0, 0.0, 0.0}};
1501
1502
// Compute reference derivatives.
1503
// Declare pointer to array of derivatives on FIAT element.
1504
double *derivatives = new double[num_derivatives];
1505
for (unsigned int r = 0; r < num_derivatives; r++)
1506
{
1507
derivatives[r] = 0.0;
1508
}// end loop over 'r'
1509
1510
// Declare derivative matrix (of polynomial basis).
1511
double dmats[10][10] = \
1512
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1513
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1514
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1515
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1516
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1517
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
1518
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
1519
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
1520
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
1521
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
1522
1523
// Declare (auxiliary) derivative matrix (of polynomial basis).
1524
double dmats_old[10][10] = \
1525
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1526
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1527
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1528
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1529
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1530
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
1531
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
1532
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
1533
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
1534
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
1535
1536
// Loop possible derivatives.
1537
for (unsigned int r = 0; r < num_derivatives; r++)
1538
{
1539
// Resetting dmats values to compute next derivative.
1540
for (unsigned int t = 0; t < 10; t++)
1541
{
1542
for (unsigned int u = 0; u < 10; u++)
1543
{
1544
dmats[t][u] = 0.0;
1545
if (t == u)
1546
{
1547
dmats[t][u] = 1.0;
1548
}
1549
1550
}// end loop over 'u'
1551
}// end loop over 't'
1552
1553
// Looping derivative order to generate dmats.
1554
for (unsigned int s = 0; s < n; s++)
1555
{
1556
// Updating dmats_old with new values and resetting dmats.
1557
for (unsigned int t = 0; t < 10; t++)
1558
{
1559
for (unsigned int u = 0; u < 10; u++)
1560
{
1561
dmats_old[t][u] = dmats[t][u];
1562
dmats[t][u] = 0.0;
1563
}// end loop over 'u'
1564
}// end loop over 't'
1565
1566
// Update dmats using an inner product.
1567
if (combinations[r][s] == 0)
1568
{
1569
for (unsigned int t = 0; t < 10; t++)
1570
{
1571
for (unsigned int u = 0; u < 10; u++)
1572
{
1573
for (unsigned int tu = 0; tu < 10; tu++)
1574
{
1575
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1576
}// end loop over 'tu'
1577
}// end loop over 'u'
1578
}// end loop over 't'
1579
}
1580
1581
if (combinations[r][s] == 1)
1582
{
1583
for (unsigned int t = 0; t < 10; t++)
1584
{
1585
for (unsigned int u = 0; u < 10; u++)
1586
{
1587
for (unsigned int tu = 0; tu < 10; tu++)
1588
{
1589
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1590
}// end loop over 'tu'
1591
}// end loop over 'u'
1592
}// end loop over 't'
1593
}
1594
1595
}// end loop over 's'
1596
for (unsigned int s = 0; s < 10; s++)
1597
{
1598
for (unsigned int t = 0; t < 10; t++)
1599
{
1600
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1601
}// end loop over 't'
1602
}// end loop over 's'
1603
}// end loop over 'r'
1604
1605
// Transform derivatives back to physical element
1606
for (unsigned int r = 0; r < num_derivatives; r++)
1607
{
1608
for (unsigned int s = 0; s < num_derivatives; s++)
1609
{
1610
values[r] += transform[r][s]*derivatives[s];
1611
}// end loop over 's'
1612
}// end loop over 'r'
1613
1614
// Delete pointer to array of derivatives on FIAT element
1615
delete [] derivatives;
1616
1617
// Delete pointer to array of combinations of derivatives and transform
1618
for (unsigned int r = 0; r < num_derivatives; r++)
1619
{
1620
delete [] combinations[r];
1621
}// end loop over 'r'
1622
delete [] combinations;
1623
for (unsigned int r = 0; r < num_derivatives; r++)
1624
{
1625
delete [] transform[r];
1626
}// end loop over 'r'
1627
delete [] transform;
1628
break;
1629
}
1630
case 5:
1631
{
1632
1633
// Array of basisvalues
1634
double basisvalues[10] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
1635
1636
// Declare helper variables
1637
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1638
double tmp1 = (1.0 - Y)/2.0;
1639
double tmp2 = tmp1*tmp1;
1640
1641
// Compute basisvalues
1642
basisvalues[0] = 1.0;
1643
basisvalues[1] = tmp0;
1644
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
1645
basisvalues[6] = basisvalues[3]*1.66666666666667*tmp0 - basisvalues[1]*0.666666666666667*tmp2;
1646
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1647
basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
1648
basisvalues[7] = basisvalues[3]*(2.5 + 3.5*Y);
1649
basisvalues[5] = basisvalues[2]*(0.111111111111111 + Y*1.66666666666667) - basisvalues[0]*0.555555555555556;
1650
basisvalues[9] = basisvalues[5]*(0.05 + Y*1.75) - basisvalues[2]*0.7;
1651
basisvalues[8] = basisvalues[4]*(0.54 + Y*2.1) - basisvalues[1]*0.56;
1652
basisvalues[0] *= std::sqrt(0.5);
1653
basisvalues[2] *= std::sqrt(1.0);
1654
basisvalues[5] *= std::sqrt(1.5);
1655
basisvalues[9] *= std::sqrt(2.0);
1656
basisvalues[1] *= std::sqrt(3.0);
1657
basisvalues[4] *= std::sqrt(4.5);
1658
basisvalues[8] *= std::sqrt(6.0);
1659
basisvalues[3] *= std::sqrt(7.5);
1660
basisvalues[7] *= std::sqrt(10.0);
1661
basisvalues[6] *= std::sqrt(14.0);
1662
1663
// Table(s) of coefficients
1664
static const double coefficients0[10] = \
1665
{0.106066017177982, -0.259807621135332, -0.15, 0.117369119465393, -0.0606091526731326, -0.0787335988751736, 0.0, 0.101644639076841, 0.131222664791956, 0.090913729009699};
1666
1667
// Tables of derivatives of the polynomial base (transpose).
1668
static const double dmats0[10][10] = \
1669
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1670
{4.89897948556636, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1671
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1672
{0.0, 9.48683298050514, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1673
{4, 0.0, 7.07106781186548, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1674
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1675
{5.29150262212918, 0.0, -2.99332590941915, 13.6626010212795, 0.0, 0.611010092660778, 0.0, 0.0, 0.0, 0.0},
1676
{0.0, 4.38178046004133, 0.0, 0.0, 12.5219806739988, 0.0, 0.0, 0.0, 0.0, 0.0},
1677
{3.46410161513775, 0.0, 7.83836717690617, 0.0, 0.0, 8.4, 0.0, 0.0, 0.0, 0.0},
1678
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
1679
1680
static const double dmats1[10][10] = \
1681
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1682
{2.44948974278318, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1683
{4.24264068711928, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1684
{2.58198889747161, 4.74341649025257, -0.912870929175277, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1685
{2.0, 6.12372435695795, 3.53553390593274, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1686
{-2.30940107675851, 0.0, 8.16496580927726, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1687
{2.64575131106459, 5.18459255872629, -1.49666295470958, 6.83130051063973, -1.05830052442584, 0.30550504633039, 0.0, 0.0, 0.0, 0.0},
1688
{2.23606797749979, 2.19089023002067, 2.52982212813471, 8.08290376865476, 6.26099033699941, -1.80739222823013, 0.0, 0.0, 0.0, 0.0},
1689
{1.73205080756887, -5.09116882454314, 3.91918358845309, 0.0, 9.69948452238572, 4.2, 0.0, 0.0, 0.0, 0.0},
1690
{5.0, 0.0, -2.82842712474619, 0.0, 0.0, 12.1243556529821, 0.0, 0.0, 0.0, 0.0}};
1691
1692
// Compute reference derivatives.
1693
// Declare pointer to array of derivatives on FIAT element.
1694
double *derivatives = new double[num_derivatives];
1695
for (unsigned int r = 0; r < num_derivatives; r++)
1696
{
1697
derivatives[r] = 0.0;
1698
}// end loop over 'r'
1699
1700
// Declare derivative matrix (of polynomial basis).
1701
double dmats[10][10] = \
1702
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1703
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1704
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1705
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1706
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1707
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
1708
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
1709
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
1710
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
1711
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
1712
1713
// Declare (auxiliary) derivative matrix (of polynomial basis).
1714
double dmats_old[10][10] = \
1715
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1716
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1717
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1718
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1719
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1720
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
1721
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
1722
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
1723
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
1724
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
1725
1726
// Loop possible derivatives.
1727
for (unsigned int r = 0; r < num_derivatives; r++)
1728
{
1729
// Resetting dmats values to compute next derivative.
1730
for (unsigned int t = 0; t < 10; t++)
1731
{
1732
for (unsigned int u = 0; u < 10; u++)
1733
{
1734
dmats[t][u] = 0.0;
1735
if (t == u)
1736
{
1737
dmats[t][u] = 1.0;
1738
}
1739
1740
}// end loop over 'u'
1741
}// end loop over 't'
1742
1743
// Looping derivative order to generate dmats.
1744
for (unsigned int s = 0; s < n; s++)
1745
{
1746
// Updating dmats_old with new values and resetting dmats.
1747
for (unsigned int t = 0; t < 10; t++)
1748
{
1749
for (unsigned int u = 0; u < 10; u++)
1750
{
1751
dmats_old[t][u] = dmats[t][u];
1752
dmats[t][u] = 0.0;
1753
}// end loop over 'u'
1754
}// end loop over 't'
1755
1756
// Update dmats using an inner product.
1757
if (combinations[r][s] == 0)
1758
{
1759
for (unsigned int t = 0; t < 10; t++)
1760
{
1761
for (unsigned int u = 0; u < 10; u++)
1762
{
1763
for (unsigned int tu = 0; tu < 10; tu++)
1764
{
1765
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1766
}// end loop over 'tu'
1767
}// end loop over 'u'
1768
}// end loop over 't'
1769
}
1770
1771
if (combinations[r][s] == 1)
1772
{
1773
for (unsigned int t = 0; t < 10; t++)
1774
{
1775
for (unsigned int u = 0; u < 10; u++)
1776
{
1777
for (unsigned int tu = 0; tu < 10; tu++)
1778
{
1779
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1780
}// end loop over 'tu'
1781
}// end loop over 'u'
1782
}// end loop over 't'
1783
}
1784
1785
}// end loop over 's'
1786
for (unsigned int s = 0; s < 10; s++)
1787
{
1788
for (unsigned int t = 0; t < 10; t++)
1789
{
1790
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1791
}// end loop over 't'
1792
}// end loop over 's'
1793
}// end loop over 'r'
1794
1795
// Transform derivatives back to physical element
1796
for (unsigned int r = 0; r < num_derivatives; r++)
1797
{
1798
for (unsigned int s = 0; s < num_derivatives; s++)
1799
{
1800
values[r] += transform[r][s]*derivatives[s];
1801
}// end loop over 's'
1802
}// end loop over 'r'
1803
1804
// Delete pointer to array of derivatives on FIAT element
1805
delete [] derivatives;
1806
1807
// Delete pointer to array of combinations of derivatives and transform
1808
for (unsigned int r = 0; r < num_derivatives; r++)
1809
{
1810
delete [] combinations[r];
1811
}// end loop over 'r'
1812
delete [] combinations;
1813
for (unsigned int r = 0; r < num_derivatives; r++)
1814
{
1815
delete [] transform[r];
1816
}// end loop over 'r'
1817
delete [] transform;
1818
break;
1819
}
1820
case 6:
1821
{
1822
1823
// Array of basisvalues
1824
double basisvalues[10] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
1825
1826
// Declare helper variables
1827
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1828
double tmp1 = (1.0 - Y)/2.0;
1829
double tmp2 = tmp1*tmp1;
1830
1831
// Compute basisvalues
1832
basisvalues[0] = 1.0;
1833
basisvalues[1] = tmp0;
1834
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
1835
basisvalues[6] = basisvalues[3]*1.66666666666667*tmp0 - basisvalues[1]*0.666666666666667*tmp2;
1836
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1837
basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
1838
basisvalues[7] = basisvalues[3]*(2.5 + 3.5*Y);
1839
basisvalues[5] = basisvalues[2]*(0.111111111111111 + Y*1.66666666666667) - basisvalues[0]*0.555555555555556;
1840
basisvalues[9] = basisvalues[5]*(0.05 + Y*1.75) - basisvalues[2]*0.7;
1841
basisvalues[8] = basisvalues[4]*(0.54 + Y*2.1) - basisvalues[1]*0.56;
1842
basisvalues[0] *= std::sqrt(0.5);
1843
basisvalues[2] *= std::sqrt(1.0);
1844
basisvalues[5] *= std::sqrt(1.5);
1845
basisvalues[9] *= std::sqrt(2.0);
1846
basisvalues[1] *= std::sqrt(3.0);
1847
basisvalues[4] *= std::sqrt(4.5);
1848
basisvalues[8] *= std::sqrt(6.0);
1849
basisvalues[3] *= std::sqrt(7.5);
1850
basisvalues[7] *= std::sqrt(10.0);
1851
basisvalues[6] *= std::sqrt(14.0);
1852
1853
// Table(s) of coefficients
1854
static const double coefficients0[10] = \
1855
{0.106066017177982, 0.0, 0.3, 0.0, -0.151522881682832, 0.0262445329583912, 0.0, 0.0, -0.131222664791956, -0.136370593514548};
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
{4.89897948556636, 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, 9.48683298050514, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1863
{4, 0.0, 7.07106781186548, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1864
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1865
{5.29150262212918, 0.0, -2.99332590941915, 13.6626010212795, 0.0, 0.611010092660778, 0.0, 0.0, 0.0, 0.0},
1866
{0.0, 4.38178046004133, 0.0, 0.0, 12.5219806739988, 0.0, 0.0, 0.0, 0.0, 0.0},
1867
{3.46410161513775, 0.0, 7.83836717690617, 0.0, 0.0, 8.4, 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
{2.44948974278318, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1873
{4.24264068711928, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1874
{2.58198889747161, 4.74341649025257, -0.912870929175277, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1875
{2.0, 6.12372435695795, 3.53553390593274, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1876
{-2.30940107675851, 0.0, 8.16496580927726, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1877
{2.64575131106459, 5.18459255872629, -1.49666295470958, 6.83130051063973, -1.05830052442584, 0.30550504633039, 0.0, 0.0, 0.0, 0.0},
1878
{2.23606797749979, 2.19089023002067, 2.52982212813471, 8.08290376865476, 6.26099033699941, -1.80739222823013, 0.0, 0.0, 0.0, 0.0},
1879
{1.73205080756887, -5.09116882454314, 3.91918358845309, 0.0, 9.69948452238572, 4.2, 0.0, 0.0, 0.0, 0.0},
1880
{5.0, 0.0, -2.82842712474619, 0.0, 0.0, 12.1243556529821, 0.0, 0.0, 0.0, 0.0}};
1881
1882
// Compute reference derivatives.
1883
// Declare pointer to array of derivatives on FIAT element.
1884
double *derivatives = new double[num_derivatives];
1885
for (unsigned int r = 0; r < num_derivatives; r++)
1886
{
1887
derivatives[r] = 0.0;
1888
}// end loop over 'r'
1889
1890
// Declare derivative matrix (of polynomial basis).
1891
double dmats[10][10] = \
1892
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1893
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1894
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1895
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1896
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1897
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
1898
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
1899
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
1900
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
1901
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
1902
1903
// Declare (auxiliary) derivative matrix (of polynomial basis).
1904
double dmats_old[10][10] = \
1905
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1906
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1907
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1908
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1909
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
1910
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
1911
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
1912
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
1913
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
1914
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
1915
1916
// Loop possible derivatives.
1917
for (unsigned int r = 0; r < num_derivatives; r++)
1918
{
1919
// Resetting dmats values to compute next derivative.
1920
for (unsigned int t = 0; t < 10; t++)
1921
{
1922
for (unsigned int u = 0; u < 10; u++)
1923
{
1924
dmats[t][u] = 0.0;
1925
if (t == u)
1926
{
1927
dmats[t][u] = 1.0;
1928
}
1929
1930
}// end loop over 'u'
1931
}// end loop over 't'
1932
1933
// Looping derivative order to generate dmats.
1934
for (unsigned int s = 0; s < n; s++)
1935
{
1936
// Updating dmats_old with new values and resetting dmats.
1937
for (unsigned int t = 0; t < 10; t++)
1938
{
1939
for (unsigned int u = 0; u < 10; u++)
1940
{
1941
dmats_old[t][u] = dmats[t][u];
1942
dmats[t][u] = 0.0;
1943
}// end loop over 'u'
1944
}// end loop over 't'
1945
1946
// Update dmats using an inner product.
1947
if (combinations[r][s] == 0)
1948
{
1949
for (unsigned int t = 0; t < 10; t++)
1950
{
1951
for (unsigned int u = 0; u < 10; u++)
1952
{
1953
for (unsigned int tu = 0; tu < 10; tu++)
1954
{
1955
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1956
}// end loop over 'tu'
1957
}// end loop over 'u'
1958
}// end loop over 't'
1959
}
1960
1961
if (combinations[r][s] == 1)
1962
{
1963
for (unsigned int t = 0; t < 10; t++)
1964
{
1965
for (unsigned int u = 0; u < 10; u++)
1966
{
1967
for (unsigned int tu = 0; tu < 10; tu++)
1968
{
1969
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1970
}// end loop over 'tu'
1971
}// end loop over 'u'
1972
}// end loop over 't'
1973
}
1974
1975
}// end loop over 's'
1976
for (unsigned int s = 0; s < 10; s++)
1977
{
1978
for (unsigned int t = 0; t < 10; t++)
1979
{
1980
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1981
}// end loop over 't'
1982
}// end loop over 's'
1983
}// end loop over 'r'
1984
1985
// Transform derivatives back to physical element
1986
for (unsigned int r = 0; r < num_derivatives; r++)
1987
{
1988
for (unsigned int s = 0; s < num_derivatives; s++)
1989
{
1990
values[r] += transform[r][s]*derivatives[s];
1991
}// end loop over 's'
1992
}// end loop over 'r'
1993
1994
// Delete pointer to array of derivatives on FIAT element
1995
delete [] derivatives;
1996
1997
// Delete pointer to array of combinations of derivatives and transform
1998
for (unsigned int r = 0; r < num_derivatives; r++)
1999
{
2000
delete [] combinations[r];
2001
}// end loop over 'r'
2002
delete [] combinations;
2003
for (unsigned int r = 0; r < num_derivatives; r++)
2004
{
2005
delete [] transform[r];
2006
}// end loop over 'r'
2007
delete [] transform;
2008
break;
2009
}
2010
case 7:
2011
{
2012
2013
// Array of basisvalues
2014
double basisvalues[10] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
2015
2016
// Declare helper variables
2017
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2018
double tmp1 = (1.0 - Y)/2.0;
2019
double tmp2 = tmp1*tmp1;
2020
2021
// Compute basisvalues
2022
basisvalues[0] = 1.0;
2023
basisvalues[1] = tmp0;
2024
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
2025
basisvalues[6] = basisvalues[3]*1.66666666666667*tmp0 - basisvalues[1]*0.666666666666667*tmp2;
2026
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
2027
basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
2028
basisvalues[7] = basisvalues[3]*(2.5 + 3.5*Y);
2029
basisvalues[5] = basisvalues[2]*(0.111111111111111 + Y*1.66666666666667) - basisvalues[0]*0.555555555555556;
2030
basisvalues[9] = basisvalues[5]*(0.05 + Y*1.75) - basisvalues[2]*0.7;
2031
basisvalues[8] = basisvalues[4]*(0.54 + Y*2.1) - basisvalues[1]*0.56;
2032
basisvalues[0] *= std::sqrt(0.5);
2033
basisvalues[2] *= std::sqrt(1.0);
2034
basisvalues[5] *= std::sqrt(1.5);
2035
basisvalues[9] *= std::sqrt(2.0);
2036
basisvalues[1] *= std::sqrt(3.0);
2037
basisvalues[4] *= std::sqrt(4.5);
2038
basisvalues[8] *= std::sqrt(6.0);
2039
basisvalues[3] *= std::sqrt(7.5);
2040
basisvalues[7] *= std::sqrt(10.0);
2041
basisvalues[6] *= std::sqrt(14.0);
2042
2043
// Table(s) of coefficients
2044
static const double coefficients0[10] = \
2045
{0.106066017177982, -0.259807621135332, -0.15, -0.0782460796435952, 0.090913729009699, 0.0962299541807677, 0.180401338290886, 0.0508223195384204, -0.0131222664791956, -0.0227284322524247};
2046
2047
// Tables of derivatives of the polynomial base (transpose).
2048
static const double dmats0[10][10] = \
2049
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2050
{4.89897948556636, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2051
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2052
{0.0, 9.48683298050514, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2053
{4, 0.0, 7.07106781186548, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2054
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2055
{5.29150262212918, 0.0, -2.99332590941915, 13.6626010212795, 0.0, 0.611010092660778, 0.0, 0.0, 0.0, 0.0},
2056
{0.0, 4.38178046004133, 0.0, 0.0, 12.5219806739988, 0.0, 0.0, 0.0, 0.0, 0.0},
2057
{3.46410161513775, 0.0, 7.83836717690617, 0.0, 0.0, 8.4, 0.0, 0.0, 0.0, 0.0},
2058
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
2059
2060
static const double dmats1[10][10] = \
2061
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2062
{2.44948974278318, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2063
{4.24264068711928, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2064
{2.58198889747161, 4.74341649025257, -0.912870929175277, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2065
{2.0, 6.12372435695795, 3.53553390593274, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2066
{-2.30940107675851, 0.0, 8.16496580927726, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2067
{2.64575131106459, 5.18459255872629, -1.49666295470958, 6.83130051063973, -1.05830052442584, 0.30550504633039, 0.0, 0.0, 0.0, 0.0},
2068
{2.23606797749979, 2.19089023002067, 2.52982212813471, 8.08290376865476, 6.26099033699941, -1.80739222823013, 0.0, 0.0, 0.0, 0.0},
2069
{1.73205080756887, -5.09116882454314, 3.91918358845309, 0.0, 9.69948452238572, 4.2, 0.0, 0.0, 0.0, 0.0},
2070
{5.0, 0.0, -2.82842712474619, 0.0, 0.0, 12.1243556529821, 0.0, 0.0, 0.0, 0.0}};
2071
2072
// Compute reference derivatives.
2073
// Declare pointer to array of derivatives on FIAT element.
2074
double *derivatives = new double[num_derivatives];
2075
for (unsigned int r = 0; r < num_derivatives; r++)
2076
{
2077
derivatives[r] = 0.0;
2078
}// end loop over 'r'
2079
2080
// Declare derivative matrix (of polynomial basis).
2081
double dmats[10][10] = \
2082
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2083
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2084
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2085
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2086
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2087
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
2088
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
2089
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
2090
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
2091
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
2092
2093
// Declare (auxiliary) derivative matrix (of polynomial basis).
2094
double dmats_old[10][10] = \
2095
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2096
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2097
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2098
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2099
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2100
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
2101
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
2102
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
2103
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
2104
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
2105
2106
// Loop possible derivatives.
2107
for (unsigned int r = 0; r < num_derivatives; r++)
2108
{
2109
// Resetting dmats values to compute next derivative.
2110
for (unsigned int t = 0; t < 10; t++)
2111
{
2112
for (unsigned int u = 0; u < 10; u++)
2113
{
2114
dmats[t][u] = 0.0;
2115
if (t == u)
2116
{
2117
dmats[t][u] = 1.0;
2118
}
2119
2120
}// end loop over 'u'
2121
}// end loop over 't'
2122
2123
// Looping derivative order to generate dmats.
2124
for (unsigned int s = 0; s < n; s++)
2125
{
2126
// Updating dmats_old with new values and resetting dmats.
2127
for (unsigned int t = 0; t < 10; t++)
2128
{
2129
for (unsigned int u = 0; u < 10; u++)
2130
{
2131
dmats_old[t][u] = dmats[t][u];
2132
dmats[t][u] = 0.0;
2133
}// end loop over 'u'
2134
}// end loop over 't'
2135
2136
// Update dmats using an inner product.
2137
if (combinations[r][s] == 0)
2138
{
2139
for (unsigned int t = 0; t < 10; t++)
2140
{
2141
for (unsigned int u = 0; u < 10; u++)
2142
{
2143
for (unsigned int tu = 0; tu < 10; tu++)
2144
{
2145
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2146
}// end loop over 'tu'
2147
}// end loop over 'u'
2148
}// end loop over 't'
2149
}
2150
2151
if (combinations[r][s] == 1)
2152
{
2153
for (unsigned int t = 0; t < 10; t++)
2154
{
2155
for (unsigned int u = 0; u < 10; u++)
2156
{
2157
for (unsigned int tu = 0; tu < 10; tu++)
2158
{
2159
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2160
}// end loop over 'tu'
2161
}// end loop over 'u'
2162
}// end loop over 't'
2163
}
2164
2165
}// end loop over 's'
2166
for (unsigned int s = 0; s < 10; s++)
2167
{
2168
for (unsigned int t = 0; t < 10; t++)
2169
{
2170
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2171
}// end loop over 't'
2172
}// end loop over 's'
2173
}// end loop over 'r'
2174
2175
// Transform derivatives back to physical element
2176
for (unsigned int r = 0; r < num_derivatives; r++)
2177
{
2178
for (unsigned int s = 0; s < num_derivatives; s++)
2179
{
2180
values[r] += transform[r][s]*derivatives[s];
2181
}// end loop over 's'
2182
}// end loop over 'r'
2183
2184
// Delete pointer to array of derivatives on FIAT element
2185
delete [] derivatives;
2186
2187
// Delete pointer to array of combinations of derivatives and transform
2188
for (unsigned int r = 0; r < num_derivatives; r++)
2189
{
2190
delete [] combinations[r];
2191
}// end loop over 'r'
2192
delete [] combinations;
2193
for (unsigned int r = 0; r < num_derivatives; r++)
2194
{
2195
delete [] transform[r];
2196
}// end loop over 'r'
2197
delete [] transform;
2198
break;
2199
}
2200
case 8:
2201
{
2202
2203
// Array of basisvalues
2204
double basisvalues[10] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
2205
2206
// Declare helper variables
2207
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2208
double tmp1 = (1.0 - Y)/2.0;
2209
double tmp2 = tmp1*tmp1;
2210
2211
// Compute basisvalues
2212
basisvalues[0] = 1.0;
2213
basisvalues[1] = tmp0;
2214
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
2215
basisvalues[6] = basisvalues[3]*1.66666666666667*tmp0 - basisvalues[1]*0.666666666666667*tmp2;
2216
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
2217
basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
2218
basisvalues[7] = basisvalues[3]*(2.5 + 3.5*Y);
2219
basisvalues[5] = basisvalues[2]*(0.111111111111111 + Y*1.66666666666667) - basisvalues[0]*0.555555555555556;
2220
basisvalues[9] = basisvalues[5]*(0.05 + Y*1.75) - basisvalues[2]*0.7;
2221
basisvalues[8] = basisvalues[4]*(0.54 + Y*2.1) - basisvalues[1]*0.56;
2222
basisvalues[0] *= std::sqrt(0.5);
2223
basisvalues[2] *= std::sqrt(1.0);
2224
basisvalues[5] *= std::sqrt(1.5);
2225
basisvalues[9] *= std::sqrt(2.0);
2226
basisvalues[1] *= std::sqrt(3.0);
2227
basisvalues[4] *= std::sqrt(4.5);
2228
basisvalues[8] *= std::sqrt(6.0);
2229
basisvalues[3] *= std::sqrt(7.5);
2230
basisvalues[7] *= std::sqrt(10.0);
2231
basisvalues[6] *= std::sqrt(14.0);
2232
2233
// Table(s) of coefficients
2234
static const double coefficients0[10] = \
2235
{0.106066017177982, 0.259807621135332, -0.15, -0.0782460796435952, -0.090913729009699, 0.0962299541807678, -0.180401338290886, 0.0508223195384204, 0.0131222664791956, -0.0227284322524248};
2236
2237
// Tables of derivatives of the polynomial base (transpose).
2238
static const double dmats0[10][10] = \
2239
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2240
{4.89897948556636, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2241
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2242
{0.0, 9.48683298050514, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2243
{4, 0.0, 7.07106781186548, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2244
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2245
{5.29150262212918, 0.0, -2.99332590941915, 13.6626010212795, 0.0, 0.611010092660778, 0.0, 0.0, 0.0, 0.0},
2246
{0.0, 4.38178046004133, 0.0, 0.0, 12.5219806739988, 0.0, 0.0, 0.0, 0.0, 0.0},
2247
{3.46410161513775, 0.0, 7.83836717690617, 0.0, 0.0, 8.4, 0.0, 0.0, 0.0, 0.0},
2248
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
2249
2250
static const double dmats1[10][10] = \
2251
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2252
{2.44948974278318, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2253
{4.24264068711928, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2254
{2.58198889747161, 4.74341649025257, -0.912870929175277, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2255
{2.0, 6.12372435695795, 3.53553390593274, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2256
{-2.30940107675851, 0.0, 8.16496580927726, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2257
{2.64575131106459, 5.18459255872629, -1.49666295470958, 6.83130051063973, -1.05830052442584, 0.30550504633039, 0.0, 0.0, 0.0, 0.0},
2258
{2.23606797749979, 2.19089023002067, 2.52982212813471, 8.08290376865476, 6.26099033699941, -1.80739222823013, 0.0, 0.0, 0.0, 0.0},
2259
{1.73205080756887, -5.09116882454314, 3.91918358845309, 0.0, 9.69948452238572, 4.2, 0.0, 0.0, 0.0, 0.0},
2260
{5.0, 0.0, -2.82842712474619, 0.0, 0.0, 12.1243556529821, 0.0, 0.0, 0.0, 0.0}};
2261
2262
// Compute reference derivatives.
2263
// Declare pointer to array of derivatives on FIAT element.
2264
double *derivatives = new double[num_derivatives];
2265
for (unsigned int r = 0; r < num_derivatives; r++)
2266
{
2267
derivatives[r] = 0.0;
2268
}// end loop over 'r'
2269
2270
// Declare derivative matrix (of polynomial basis).
2271
double dmats[10][10] = \
2272
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2273
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2274
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2275
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2276
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2277
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
2278
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
2279
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
2280
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
2281
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
2282
2283
// Declare (auxiliary) derivative matrix (of polynomial basis).
2284
double dmats_old[10][10] = \
2285
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2286
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2287
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2288
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2289
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2290
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
2291
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
2292
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
2293
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
2294
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
2295
2296
// Loop possible derivatives.
2297
for (unsigned int r = 0; r < num_derivatives; r++)
2298
{
2299
// Resetting dmats values to compute next derivative.
2300
for (unsigned int t = 0; t < 10; t++)
2301
{
2302
for (unsigned int u = 0; u < 10; u++)
2303
{
2304
dmats[t][u] = 0.0;
2305
if (t == u)
2306
{
2307
dmats[t][u] = 1.0;
2308
}
2309
2310
}// end loop over 'u'
2311
}// end loop over 't'
2312
2313
// Looping derivative order to generate dmats.
2314
for (unsigned int s = 0; s < n; s++)
2315
{
2316
// Updating dmats_old with new values and resetting dmats.
2317
for (unsigned int t = 0; t < 10; t++)
2318
{
2319
for (unsigned int u = 0; u < 10; u++)
2320
{
2321
dmats_old[t][u] = dmats[t][u];
2322
dmats[t][u] = 0.0;
2323
}// end loop over 'u'
2324
}// end loop over 't'
2325
2326
// Update dmats using an inner product.
2327
if (combinations[r][s] == 0)
2328
{
2329
for (unsigned int t = 0; t < 10; t++)
2330
{
2331
for (unsigned int u = 0; u < 10; u++)
2332
{
2333
for (unsigned int tu = 0; tu < 10; tu++)
2334
{
2335
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2336
}// end loop over 'tu'
2337
}// end loop over 'u'
2338
}// end loop over 't'
2339
}
2340
2341
if (combinations[r][s] == 1)
2342
{
2343
for (unsigned int t = 0; t < 10; t++)
2344
{
2345
for (unsigned int u = 0; u < 10; u++)
2346
{
2347
for (unsigned int tu = 0; tu < 10; tu++)
2348
{
2349
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2350
}// end loop over 'tu'
2351
}// end loop over 'u'
2352
}// end loop over 't'
2353
}
2354
2355
}// end loop over 's'
2356
for (unsigned int s = 0; s < 10; s++)
2357
{
2358
for (unsigned int t = 0; t < 10; t++)
2359
{
2360
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2361
}// end loop over 't'
2362
}// end loop over 's'
2363
}// end loop over 'r'
2364
2365
// Transform derivatives back to physical element
2366
for (unsigned int r = 0; r < num_derivatives; r++)
2367
{
2368
for (unsigned int s = 0; s < num_derivatives; s++)
2369
{
2370
values[r] += transform[r][s]*derivatives[s];
2371
}// end loop over 's'
2372
}// end loop over 'r'
2373
2374
// Delete pointer to array of derivatives on FIAT element
2375
delete [] derivatives;
2376
2377
// Delete pointer to array of combinations of derivatives and transform
2378
for (unsigned int r = 0; r < num_derivatives; r++)
2379
{
2380
delete [] combinations[r];
2381
}// end loop over 'r'
2382
delete [] combinations;
2383
for (unsigned int r = 0; r < num_derivatives; r++)
2384
{
2385
delete [] transform[r];
2386
}// end loop over 'r'
2387
delete [] transform;
2388
break;
2389
}
2390
case 9:
2391
{
2392
2393
// Array of basisvalues
2394
double basisvalues[10] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
2395
2396
// Declare helper variables
2397
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2398
double tmp1 = (1.0 - Y)/2.0;
2399
double tmp2 = tmp1*tmp1;
2400
2401
// Compute basisvalues
2402
basisvalues[0] = 1.0;
2403
basisvalues[1] = tmp0;
2404
basisvalues[3] = basisvalues[1]*1.5*tmp0 - basisvalues[0]*0.5*tmp2;
2405
basisvalues[6] = basisvalues[3]*1.66666666666667*tmp0 - basisvalues[1]*0.666666666666667*tmp2;
2406
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
2407
basisvalues[4] = basisvalues[1]*(1.5 + 2.5*Y);
2408
basisvalues[7] = basisvalues[3]*(2.5 + 3.5*Y);
2409
basisvalues[5] = basisvalues[2]*(0.111111111111111 + Y*1.66666666666667) - basisvalues[0]*0.555555555555556;
2410
basisvalues[9] = basisvalues[5]*(0.05 + Y*1.75) - basisvalues[2]*0.7;
2411
basisvalues[8] = basisvalues[4]*(0.54 + Y*2.1) - basisvalues[1]*0.56;
2412
basisvalues[0] *= std::sqrt(0.5);
2413
basisvalues[2] *= std::sqrt(1.0);
2414
basisvalues[5] *= std::sqrt(1.5);
2415
basisvalues[9] *= std::sqrt(2.0);
2416
basisvalues[1] *= std::sqrt(3.0);
2417
basisvalues[4] *= std::sqrt(4.5);
2418
basisvalues[8] *= std::sqrt(6.0);
2419
basisvalues[3] *= std::sqrt(7.5);
2420
basisvalues[7] *= std::sqrt(10.0);
2421
basisvalues[6] *= std::sqrt(14.0);
2422
2423
// Table(s) of coefficients
2424
static const double coefficients0[10] = \
2425
{0.636396103067893, 0.0, 0.0, -0.234738238930785, 0.0, -0.262445329583912, 0.0, -0.203289278153682, 0.0, 0.090913729009699};
2426
2427
// Tables of derivatives of the polynomial base (transpose).
2428
static const double dmats0[10][10] = \
2429
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2430
{4.89897948556636, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2431
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2432
{0.0, 9.48683298050514, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2433
{4, 0.0, 7.07106781186548, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2434
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2435
{5.29150262212918, 0.0, -2.99332590941915, 13.6626010212795, 0.0, 0.611010092660778, 0.0, 0.0, 0.0, 0.0},
2436
{0.0, 4.38178046004133, 0.0, 0.0, 12.5219806739988, 0.0, 0.0, 0.0, 0.0, 0.0},
2437
{3.46410161513775, 0.0, 7.83836717690617, 0.0, 0.0, 8.4, 0.0, 0.0, 0.0, 0.0},
2438
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
2439
2440
static const double dmats1[10][10] = \
2441
{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2442
{2.44948974278318, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2443
{4.24264068711928, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2444
{2.58198889747161, 4.74341649025257, -0.912870929175277, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2445
{2.0, 6.12372435695795, 3.53553390593274, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2446
{-2.30940107675851, 0.0, 8.16496580927726, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2447
{2.64575131106459, 5.18459255872629, -1.49666295470958, 6.83130051063973, -1.05830052442584, 0.30550504633039, 0.0, 0.0, 0.0, 0.0},
2448
{2.23606797749979, 2.19089023002067, 2.52982212813471, 8.08290376865476, 6.26099033699941, -1.80739222823013, 0.0, 0.0, 0.0, 0.0},
2449
{1.73205080756887, -5.09116882454314, 3.91918358845309, 0.0, 9.69948452238572, 4.2, 0.0, 0.0, 0.0, 0.0},
2450
{5.0, 0.0, -2.82842712474619, 0.0, 0.0, 12.1243556529821, 0.0, 0.0, 0.0, 0.0}};
2451
2452
// Compute reference derivatives.
2453
// Declare pointer to array of derivatives on FIAT element.
2454
double *derivatives = new double[num_derivatives];
2455
for (unsigned int r = 0; r < num_derivatives; r++)
2456
{
2457
derivatives[r] = 0.0;
2458
}// end loop over 'r'
2459
2460
// Declare derivative matrix (of polynomial basis).
2461
double dmats[10][10] = \
2462
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2463
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2464
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2465
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2466
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2467
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
2468
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
2469
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
2470
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
2471
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
2472
2473
// Declare (auxiliary) derivative matrix (of polynomial basis).
2474
double dmats_old[10][10] = \
2475
{{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2476
{0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2477
{0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2478
{0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2479
{0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2480
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0},
2481
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
2482
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
2483
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
2484
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}};
2485
2486
// Loop possible derivatives.
2487
for (unsigned int r = 0; r < num_derivatives; r++)
2488
{
2489
// Resetting dmats values to compute next derivative.
2490
for (unsigned int t = 0; t < 10; t++)
2491
{
2492
for (unsigned int u = 0; u < 10; u++)
2493
{
2494
dmats[t][u] = 0.0;
2495
if (t == u)
2496
{
2497
dmats[t][u] = 1.0;
2498
}
2499
2500
}// end loop over 'u'
2501
}// end loop over 't'
2502
2503
// Looping derivative order to generate dmats.
2504
for (unsigned int s = 0; s < n; s++)
2505
{
2506
// Updating dmats_old with new values and resetting dmats.
2507
for (unsigned int t = 0; t < 10; t++)
2508
{
2509
for (unsigned int u = 0; u < 10; u++)
2510
{
2511
dmats_old[t][u] = dmats[t][u];
2512
dmats[t][u] = 0.0;
2513
}// end loop over 'u'
2514
}// end loop over 't'
2515
2516
// Update dmats using an inner product.
2517
if (combinations[r][s] == 0)
2518
{
2519
for (unsigned int t = 0; t < 10; t++)
2520
{
2521
for (unsigned int u = 0; u < 10; u++)
2522
{
2523
for (unsigned int tu = 0; tu < 10; tu++)
2524
{
2525
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2526
}// end loop over 'tu'
2527
}// end loop over 'u'
2528
}// end loop over 't'
2529
}
2530
2531
if (combinations[r][s] == 1)
2532
{
2533
for (unsigned int t = 0; t < 10; t++)
2534
{
2535
for (unsigned int u = 0; u < 10; u++)
2536
{
2537
for (unsigned int tu = 0; tu < 10; tu++)
2538
{
2539
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2540
}// end loop over 'tu'
2541
}// end loop over 'u'
2542
}// end loop over 't'
2543
}
2544
2545
}// end loop over 's'
2546
for (unsigned int s = 0; s < 10; s++)
2547
{
2548
for (unsigned int t = 0; t < 10; t++)
2549
{
2550
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2551
}// end loop over 't'
2552
}// end loop over 's'
2553
}// end loop over 'r'
2554
2555
// Transform derivatives back to physical element
2556
for (unsigned int r = 0; r < num_derivatives; r++)
2557
{
2558
for (unsigned int s = 0; s < num_derivatives; s++)
2559
{
2560
values[r] += transform[r][s]*derivatives[s];
2561
}// end loop over 's'
2562
}// end loop over 'r'
2563
2564
// Delete pointer to array of derivatives on FIAT element
2565
delete [] derivatives;
2566
2567
// Delete pointer to array of combinations of derivatives and transform
2568
for (unsigned int r = 0; r < num_derivatives; r++)
2569
{
2570
delete [] combinations[r];
2571
}// end loop over 'r'
2572
delete [] combinations;
2573
for (unsigned int r = 0; r < num_derivatives; r++)
2574
{
2575
delete [] transform[r];
2576
}// end loop over 'r'
2577
delete [] transform;
2578
break;
2579
}
2580
}
2581
2582
}
2583
2584
/// Evaluate order n derivatives of all basis functions at given point x in cell
2585
virtual void evaluate_basis_derivatives_all(std::size_t n,
2586
double* values,
2587
const double* x,
2588
const double* vertex_coordinates,
2589
int cell_orientation) const
2590
{
2591
// Compute number of derivatives.
2592
unsigned int num_derivatives = 1;
2593
for (unsigned int r = 0; r < n; r++)
2594
{
2595
num_derivatives *= 2;
2596
}// end loop over 'r'
2597
2598
// Helper variable to hold values of a single dof.
2599
double *dof_values = new double[num_derivatives];
2600
for (unsigned int r = 0; r < num_derivatives; r++)
2601
{
2602
dof_values[r] = 0.0;
2603
}// end loop over 'r'
2604
2605
// Loop dofs and call evaluate_basis_derivatives.
2606
for (unsigned int r = 0; r < 10; r++)
2607
{
2608
evaluate_basis_derivatives(r, n, dof_values, x, vertex_coordinates, cell_orientation);
2609
for (unsigned int s = 0; s < num_derivatives; s++)
2610
{
2611
values[r*num_derivatives + s] = dof_values[s];
2612
}// end loop over 's'
2613
}// end loop over 'r'
2614
2615
// Delete pointer.
2616
delete [] dof_values;
2617
}
2618
2619
/// Evaluate linear functional for dof i on the function f
2620
virtual double evaluate_dof(std::size_t i,
2621
const ufc::function& f,
2622
const double* vertex_coordinates,
2623
int cell_orientation,
2624
const ufc::cell& c) const
2625
{
2626
// Declare variables for result of evaluation
2627
double vals[1];
2628
2629
// Declare variable for physical coordinates
2630
double y[2];
2631
switch (i)
2632
{
2633
case 0:
2634
{
2635
y[0] = vertex_coordinates[0];
2636
y[1] = vertex_coordinates[1];
2637
f.evaluate(vals, y, c);
2638
return vals[0];
2639
break;
2640
}
2641
case 1:
2642
{
2643
y[0] = vertex_coordinates[2];
2644
y[1] = vertex_coordinates[3];
2645
f.evaluate(vals, y, c);
2646
return vals[0];
2647
break;
2648
}
2649
case 2:
2650
{
2651
y[0] = vertex_coordinates[4];
2652
y[1] = vertex_coordinates[5];
2653
f.evaluate(vals, y, c);
2654
return vals[0];
2655
break;
2656
}
2657
case 3:
2658
{
2659
y[0] = 0.666666666666667*vertex_coordinates[2] + 0.333333333333333*vertex_coordinates[4];
2660
y[1] = 0.666666666666667*vertex_coordinates[3] + 0.333333333333333*vertex_coordinates[5];
2661
f.evaluate(vals, y, c);
2662
return vals[0];
2663
break;
2664
}
2665
case 4:
2666
{
2667
y[0] = 0.333333333333333*vertex_coordinates[2] + 0.666666666666667*vertex_coordinates[4];
2668
y[1] = 0.333333333333333*vertex_coordinates[3] + 0.666666666666667*vertex_coordinates[5];
2669
f.evaluate(vals, y, c);
2670
return vals[0];
2671
break;
2672
}
2673
case 5:
2674
{
2675
y[0] = 0.666666666666667*vertex_coordinates[0] + 0.333333333333333*vertex_coordinates[4];
2676
y[1] = 0.666666666666667*vertex_coordinates[1] + 0.333333333333333*vertex_coordinates[5];
2677
f.evaluate(vals, y, c);
2678
return vals[0];
2679
break;
2680
}
2681
case 6:
2682
{
2683
y[0] = 0.333333333333333*vertex_coordinates[0] + 0.666666666666667*vertex_coordinates[4];
2684
y[1] = 0.333333333333333*vertex_coordinates[1] + 0.666666666666667*vertex_coordinates[5];
2685
f.evaluate(vals, y, c);
2686
return vals[0];
2687
break;
2688
}
2689
case 7:
2690
{
2691
y[0] = 0.666666666666667*vertex_coordinates[0] + 0.333333333333333*vertex_coordinates[2];
2692
y[1] = 0.666666666666667*vertex_coordinates[1] + 0.333333333333333*vertex_coordinates[3];
2693
f.evaluate(vals, y, c);
2694
return vals[0];
2695
break;
2696
}
2697
case 8:
2698
{
2699
y[0] = 0.333333333333333*vertex_coordinates[0] + 0.666666666666667*vertex_coordinates[2];
2700
y[1] = 0.333333333333333*vertex_coordinates[1] + 0.666666666666667*vertex_coordinates[3];
2701
f.evaluate(vals, y, c);
2702
return vals[0];
2703
break;
2704
}
2705
case 9:
2706
{
2707
y[0] = 0.333333333333333*vertex_coordinates[0] + 0.333333333333333*vertex_coordinates[2] + 0.333333333333333*vertex_coordinates[4];
2708
y[1] = 0.333333333333333*vertex_coordinates[1] + 0.333333333333333*vertex_coordinates[3] + 0.333333333333333*vertex_coordinates[5];
2709
f.evaluate(vals, y, c);
2710
return vals[0];
2711
break;
2712
}
2713
}
2714
2715
return 0.0;
2716
}
2717
2718
/// Evaluate linear functionals for all dofs on the function f
2719
virtual void evaluate_dofs(double* values,
2720
const ufc::function& f,
2721
const double* vertex_coordinates,
2722
int cell_orientation,
2723
const ufc::cell& c) const
2724
{
2725
// Declare variables for result of evaluation
2726
double vals[1];
2727
2728
// Declare variable for physical coordinates
2729
double y[2];
2730
y[0] = vertex_coordinates[0];
2731
y[1] = vertex_coordinates[1];
2732
f.evaluate(vals, y, c);
2733
values[0] = vals[0];
2734
y[0] = vertex_coordinates[2];
2735
y[1] = vertex_coordinates[3];
2736
f.evaluate(vals, y, c);
2737
values[1] = vals[0];
2738
y[0] = vertex_coordinates[4];
2739
y[1] = vertex_coordinates[5];
2740
f.evaluate(vals, y, c);
2741
values[2] = vals[0];
2742
y[0] = 0.666666666666667*vertex_coordinates[2] + 0.333333333333333*vertex_coordinates[4];
2743
y[1] = 0.666666666666667*vertex_coordinates[3] + 0.333333333333333*vertex_coordinates[5];
2744
f.evaluate(vals, y, c);
2745
values[3] = vals[0];
2746
y[0] = 0.333333333333333*vertex_coordinates[2] + 0.666666666666667*vertex_coordinates[4];
2747
y[1] = 0.333333333333333*vertex_coordinates[3] + 0.666666666666667*vertex_coordinates[5];
2748
f.evaluate(vals, y, c);
2749
values[4] = vals[0];
2750
y[0] = 0.666666666666667*vertex_coordinates[0] + 0.333333333333333*vertex_coordinates[4];
2751
y[1] = 0.666666666666667*vertex_coordinates[1] + 0.333333333333333*vertex_coordinates[5];
2752
f.evaluate(vals, y, c);
2753
values[5] = vals[0];
2754
y[0] = 0.333333333333333*vertex_coordinates[0] + 0.666666666666667*vertex_coordinates[4];
2755
y[1] = 0.333333333333333*vertex_coordinates[1] + 0.666666666666667*vertex_coordinates[5];
2756
f.evaluate(vals, y, c);
2757
values[6] = vals[0];
2758
y[0] = 0.666666666666667*vertex_coordinates[0] + 0.333333333333333*vertex_coordinates[2];
2759
y[1] = 0.666666666666667*vertex_coordinates[1] + 0.333333333333333*vertex_coordinates[3];
2760
f.evaluate(vals, y, c);
2761
values[7] = vals[0];
2762
y[0] = 0.333333333333333*vertex_coordinates[0] + 0.666666666666667*vertex_coordinates[2];
2763
y[1] = 0.333333333333333*vertex_coordinates[1] + 0.666666666666667*vertex_coordinates[3];
2764
f.evaluate(vals, y, c);
2765
values[8] = vals[0];
2766
y[0] = 0.333333333333333*vertex_coordinates[0] + 0.333333333333333*vertex_coordinates[2] + 0.333333333333333*vertex_coordinates[4];
2767
y[1] = 0.333333333333333*vertex_coordinates[1] + 0.333333333333333*vertex_coordinates[3] + 0.333333333333333*vertex_coordinates[5];
2768
f.evaluate(vals, y, c);
2769
values[9] = vals[0];
2770
}
2771
2772
/// Interpolate vertex values from dof values
2773
virtual void interpolate_vertex_values(double* vertex_values,
2774
const double* dof_values,
2775
const double* vertex_coordinates,
2776
int cell_orientation,
2777
const ufc::cell& c) const
2778
{
2779
// Evaluate function and change variables
2780
vertex_values[0] = dof_values[0];
2781
vertex_values[1] = dof_values[1];
2782
vertex_values[2] = dof_values[2];
2783
}
2784
2785
/// Map coordinate xhat from reference cell to coordinate x in cell
2786
virtual void map_from_reference_cell(double* x,
2787
const double* xhat,
2788
const ufc::cell& c) const
2789
{
2790
throw std::runtime_error("map_from_reference_cell not yet implemented.");
2791
}
2792
2793
/// Map from coordinate x in cell to coordinate xhat in reference cell
2794
virtual void map_to_reference_cell(double* xhat,
2795
const double* x,
2796
const ufc::cell& c) const
2797
{
2798
throw std::runtime_error("map_to_reference_cell not yet implemented.");
2799
}
2800
2801
/// Return the number of sub elements (for a mixed element)
2802
virtual std::size_t num_sub_elements() const
2803
{
2804
return 0;
2805
}
2806
2807
/// Create a new finite element for sub element i (for a mixed element)
2808
virtual ufc::finite_element* create_sub_element(std::size_t i) const
2809
{
2810
return 0;
2811
}
2812
2813
/// Create a new class instance
2814
virtual ufc::finite_element* create() const
2815
{
2816
return new p3_finite_element_0();
2817
}
2818
2819
};
2820
2821
/// This class defines the interface for a local-to-global mapping of
2822
/// degrees of freedom (dofs).
2823
2824
class p3_dofmap_0: public ufc::dofmap
2825
{
2826
public:
2827
2828
/// Constructor
2829
p3_dofmap_0() : ufc::dofmap()
2830
{
2831
// Do nothing
2832
}
2833
2834
/// Destructor
2835
virtual ~p3_dofmap_0()
2836
{
2837
// Do nothing
2838
}
2839
2840
/// Return a string identifying the dofmap
2841
virtual const char* signature() const
2842
{
2843
return "FFC dofmap for FiniteElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 3, None)";
2844
}
2845
2846
/// Return true iff mesh entities of topological dimension d are needed
2847
virtual bool needs_mesh_entities(std::size_t d) const
2848
{
2849
switch (d)
2850
{
2851
case 0:
2852
{
2853
return true;
2854
break;
2855
}
2856
case 1:
2857
{
2858
return true;
2859
break;
2860
}
2861
case 2:
2862
{
2863
return true;
2864
break;
2865
}
2866
}
2867
2868
return false;
2869
}
2870
2871
/// Return the topological dimension of the associated cell shape
2872
virtual std::size_t topological_dimension() const
2873
{
2874
return 2;
2875
}
2876
2877
/// Return the geometric dimension of the associated cell shape
2878
virtual std::size_t geometric_dimension() const
2879
{
2880
return 2;
2881
}
2882
2883
/// Return the dimension of the global finite element function space
2884
virtual std::size_t global_dimension(const std::vector<std::size_t>&
2885
num_global_entities) const
2886
{
2887
return num_global_entities[0] + 2*num_global_entities[1] + num_global_entities[2];
2888
}
2889
2890
/// Return the dimension of the local finite element function space for a cell
2891
virtual std::size_t local_dimension() const
2892
{
2893
return 10;
2894
}
2895
2896
/// Return the number of dofs on each cell facet
2897
virtual std::size_t num_facet_dofs() const
2898
{
2899
return 4;
2900
}
2901
2902
/// Return the number of dofs associated with each cell entity of dimension d
2903
virtual std::size_t num_entity_dofs(std::size_t d) const
2904
{
2905
switch (d)
2906
{
2907
case 0:
2908
{
2909
return 1;
2910
break;
2911
}
2912
case 1:
2913
{
2914
return 2;
2915
break;
2916
}
2917
case 2:
2918
{
2919
return 1;
2920
break;
2921
}
2922
}
2923
2924
return 0;
2925
}
2926
2927
/// Tabulate the local-to-global mapping of dofs on a cell
2928
virtual void tabulate_dofs(std::size_t* dofs,
2929
const std::vector<std::size_t>& num_global_entities,
2930
const ufc::cell& c) const
2931
{
2932
unsigned int offset = 0;
2933
dofs[0] = offset + c.entity_indices[0][0];
2934
dofs[1] = offset + c.entity_indices[0][1];
2935
dofs[2] = offset + c.entity_indices[0][2];
2936
offset += num_global_entities[0];
2937
dofs[3] = offset + 2*c.entity_indices[1][0];
2938
dofs[4] = offset + 2*c.entity_indices[1][0] + 1;
2939
dofs[5] = offset + 2*c.entity_indices[1][1];
2940
dofs[6] = offset + 2*c.entity_indices[1][1] + 1;
2941
dofs[7] = offset + 2*c.entity_indices[1][2];
2942
dofs[8] = offset + 2*c.entity_indices[1][2] + 1;
2943
offset += 2*num_global_entities[1];
2944
dofs[9] = offset + c.entity_indices[2][0];
2945
offset += num_global_entities[2];
2946
}
2947
2948
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
2949
virtual void tabulate_facet_dofs(std::size_t* dofs,
2950
std::size_t facet) const
2951
{
2952
switch (facet)
2953
{
2954
case 0:
2955
{
2956
dofs[0] = 1;
2957
dofs[1] = 2;
2958
dofs[2] = 3;
2959
dofs[3] = 4;
2960
break;
2961
}
2962
case 1:
2963
{
2964
dofs[0] = 0;
2965
dofs[1] = 2;
2966
dofs[2] = 5;
2967
dofs[3] = 6;
2968
break;
2969
}
2970
case 2:
2971
{
2972
dofs[0] = 0;
2973
dofs[1] = 1;
2974
dofs[2] = 7;
2975
dofs[3] = 8;
2976
break;
2977
}
2978
}
2979
2980
}
2981
2982
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
2983
virtual void tabulate_entity_dofs(std::size_t* dofs,
2984
std::size_t d, std::size_t i) const
2985
{
2986
if (d > 2)
2987
{
2988
throw std::runtime_error("d is larger than dimension (2)");
2989
}
2990
2991
switch (d)
2992
{
2993
case 0:
2994
{
2995
if (i > 2)
2996
{
2997
throw std::runtime_error("i is larger than number of entities (2)");
2998
}
2999
3000
switch (i)
3001
{
3002
case 0:
3003
{
3004
dofs[0] = 0;
3005
break;
3006
}
3007
case 1:
3008
{
3009
dofs[0] = 1;
3010
break;
3011
}
3012
case 2:
3013
{
3014
dofs[0] = 2;
3015
break;
3016
}
3017
}
3018
3019
break;
3020
}
3021
case 1:
3022
{
3023
if (i > 2)
3024
{
3025
throw std::runtime_error("i is larger than number of entities (2)");
3026
}
3027
3028
switch (i)
3029
{
3030
case 0:
3031
{
3032
dofs[0] = 3;
3033
dofs[1] = 4;
3034
break;
3035
}
3036
case 1:
3037
{
3038
dofs[0] = 5;
3039
dofs[1] = 6;
3040
break;
3041
}
3042
case 2:
3043
{
3044
dofs[0] = 7;
3045
dofs[1] = 8;
3046
break;
3047
}
3048
}
3049
3050
break;
3051
}
3052
case 2:
3053
{
3054
if (i > 0)
3055
{
3056
throw std::runtime_error("i is larger than number of entities (0)");
3057
}
3058
3059
dofs[0] = 9;
3060
break;
3061
}
3062
}
3063
3064
}
3065
3066
/// Tabulate the coordinates of all dofs on a cell
3067
virtual void tabulate_coordinates(double** dof_coordinates,
3068
const double* vertex_coordinates) const
3069
{
3070
dof_coordinates[0][0] = vertex_coordinates[0];
3071
dof_coordinates[0][1] = vertex_coordinates[1];
3072
dof_coordinates[1][0] = vertex_coordinates[2];
3073
dof_coordinates[1][1] = vertex_coordinates[3];
3074
dof_coordinates[2][0] = vertex_coordinates[4];
3075
dof_coordinates[2][1] = vertex_coordinates[5];
3076
dof_coordinates[3][0] = 0.666666666666667*vertex_coordinates[2] + 0.333333333333333*vertex_coordinates[4];
3077
dof_coordinates[3][1] = 0.666666666666667*vertex_coordinates[3] + 0.333333333333333*vertex_coordinates[5];
3078
dof_coordinates[4][0] = 0.333333333333333*vertex_coordinates[2] + 0.666666666666667*vertex_coordinates[4];
3079
dof_coordinates[4][1] = 0.333333333333333*vertex_coordinates[3] + 0.666666666666667*vertex_coordinates[5];
3080
dof_coordinates[5][0] = 0.666666666666667*vertex_coordinates[0] + 0.333333333333333*vertex_coordinates[4];
3081
dof_coordinates[5][1] = 0.666666666666667*vertex_coordinates[1] + 0.333333333333333*vertex_coordinates[5];
3082
dof_coordinates[6][0] = 0.333333333333333*vertex_coordinates[0] + 0.666666666666667*vertex_coordinates[4];
3083
dof_coordinates[6][1] = 0.333333333333333*vertex_coordinates[1] + 0.666666666666667*vertex_coordinates[5];
3084
dof_coordinates[7][0] = 0.666666666666667*vertex_coordinates[0] + 0.333333333333333*vertex_coordinates[2];
3085
dof_coordinates[7][1] = 0.666666666666667*vertex_coordinates[1] + 0.333333333333333*vertex_coordinates[3];
3086
dof_coordinates[8][0] = 0.333333333333333*vertex_coordinates[0] + 0.666666666666667*vertex_coordinates[2];
3087
dof_coordinates[8][1] = 0.333333333333333*vertex_coordinates[1] + 0.666666666666667*vertex_coordinates[3];
3088
dof_coordinates[9][0] = 0.333333333333333*vertex_coordinates[0] + 0.333333333333333*vertex_coordinates[2] + 0.333333333333333*vertex_coordinates[4];
3089
dof_coordinates[9][1] = 0.333333333333333*vertex_coordinates[1] + 0.333333333333333*vertex_coordinates[3] + 0.333333333333333*vertex_coordinates[5];
3090
}
3091
3092
/// Return the number of sub dofmaps (for a mixed element)
3093
virtual std::size_t num_sub_dofmaps() const
3094
{
3095
return 0;
3096
}
3097
3098
/// Create a new dofmap for sub dofmap i (for a mixed element)
3099
virtual ufc::dofmap* create_sub_dofmap(std::size_t i) const
3100
{
3101
return 0;
3102
}
3103
3104
/// Create a new class instance
3105
virtual ufc::dofmap* create() const
3106
{
3107
return new p3_dofmap_0();
3108
}
3109
3110
};
3111
3112
// DOLFIN wrappers
3113
3114
// Standard library includes
3115
#include <string>
3116
3117
// DOLFIN includes
3118
#include <dolfin/common/NoDeleter.h>
3119
#include <dolfin/mesh/Restriction.h>
3120
#include <dolfin/fem/FiniteElement.h>
3121
#include <dolfin/fem/DofMap.h>
3122
#include <dolfin/fem/Form.h>
3123
#include <dolfin/function/FunctionSpace.h>
3124
#include <dolfin/function/GenericFunction.h>
3125
#include <dolfin/function/CoefficientAssigner.h>
3126
#include <dolfin/adaptivity/ErrorControl.h>
3127
#include <dolfin/adaptivity/GoalFunctional.h>
3128
3129
namespace P3
3130
{
3131
3132
class FunctionSpace: public dolfin::FunctionSpace
3133
{
3134
public:
3135
3136
//--- Constructors for standard function space, 2 different versions ---
3137
3138
// Create standard function space (reference version)
3139
FunctionSpace(const dolfin::Mesh& mesh):
3140
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
3141
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new p3_finite_element_0()))),
3142
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new p3_dofmap_0()), mesh)))
3143
{
3144
// Do nothing
3145
}
3146
3147
// Create standard function space (shared pointer version)
3148
FunctionSpace(boost::shared_ptr<const dolfin::Mesh> mesh):
3149
dolfin::FunctionSpace(mesh,
3150
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new p3_finite_element_0()))),
3151
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new p3_dofmap_0()), *mesh)))
3152
{
3153
// Do nothing
3154
}
3155
3156
//--- Constructors for constrained function space, 2 different versions ---
3157
3158
// Create standard function space (reference version)
3159
FunctionSpace(const dolfin::Mesh& mesh, const dolfin::SubDomain& constrained_domain):
3160
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
3161
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new p3_finite_element_0()))),
3162
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new p3_dofmap_0()), mesh,
3163
dolfin::reference_to_no_delete_pointer(constrained_domain))))
3164
{
3165
// Do nothing
3166
}
3167
3168
// Create standard function space (shared pointer version)
3169
FunctionSpace(boost::shared_ptr<const dolfin::Mesh> mesh, boost::shared_ptr<const dolfin::SubDomain> constrained_domain):
3170
dolfin::FunctionSpace(mesh,
3171
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new p3_finite_element_0()))),
3172
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new p3_dofmap_0()), *mesh, constrained_domain)))
3173
{
3174
// Do nothing
3175
}
3176
3177
//--- Constructors for restricted function space, 2 different versions ---
3178
3179
// Create restricted function space (reference version)
3180
FunctionSpace(const dolfin::Restriction& restriction):
3181
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(restriction.mesh()),
3182
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new p3_finite_element_0()))),
3183
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new p3_dofmap_0()),
3184
reference_to_no_delete_pointer(restriction))))
3185
{
3186
// Do nothing
3187
}
3188
3189
// Create restricted function space (shared pointer version)
3190
FunctionSpace(boost::shared_ptr<const dolfin::Restriction> restriction):
3191
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(restriction->mesh()),
3192
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new p3_finite_element_0()))),
3193
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new p3_dofmap_0()),
3194
restriction)))
3195
{
3196
// Do nothing
3197
}
3198
3199
// Copy constructor
3200
~FunctionSpace()
3201
{
3202
}
3203
3204
};
3205
3206
}
3207
3208
#endif
Loggerhead is a web-based interface for Breezy
Version: 2.0.1