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