1
// This code conforms with the UFC specification version 2.2.0
2
// and was automatically generated by FFC version 1.2.0.
4
// This code was generated with the following parameters:
7
// convert_exceptions_to_warnings: True
9
// cpp_optimize_flags: '-O2'
11
// error_control: False
19
// quadrature_degree: 'auto'
20
// quadrature_rule: 'auto'
21
// representation: 'auto'
23
// swig_binary: 'swig'
26
#ifndef __PROJECTIONMANIFOLD_H
27
#define __PROJECTIONMANIFOLD_H
34
/// This class defines the interface for a finite element.
36
class projectionmanifold_finite_element_0: public ufc::finite_element
41
projectionmanifold_finite_element_0() : ufc::finite_element()
47
virtual ~projectionmanifold_finite_element_0()
52
/// Return a string identifying the finite element
53
virtual const char* signature() const
55
return "FiniteElement('Raviart-Thomas', Domain(Cell('triangle', 3), 'triangle_multiverse', 3, 2), 1, None)";
58
/// Return the cell shape
59
virtual ufc::shape cell_shape() const
64
/// Return the topological dimension of the cell shape
65
virtual std::size_t topological_dimension() const
70
/// Return the geometric dimension of the cell shape
71
virtual std::size_t geometric_dimension() const
76
/// Return the dimension of the finite element function space
77
virtual std::size_t space_dimension() const
82
/// Return the rank of the value space
83
virtual std::size_t value_rank() const
88
/// Return the dimension of the value space for axis i
89
virtual std::size_t value_dimension(std::size_t i) const
103
/// Evaluate basis function i at given point x in cell
104
virtual void evaluate_basis(std::size_t i,
107
const double* vertex_coordinates,
108
int cell_orientation) const
112
compute_jacobian_triangle_3d(J, vertex_coordinates);
114
// Compute Jacobian inverse and determinant
117
compute_jacobian_inverse_triangle_3d(K, detJ, J);
121
if (cell_orientation == -1)
122
throw std::runtime_error("cell orientation must be defined (not -1)");
123
// (If cell_orientation == 1 = down, multiply det(J) by -1)
124
else if (cell_orientation == 1)
128
const double b0 = vertex_coordinates[0];
129
const double b1 = vertex_coordinates[1];
130
const double b2 = vertex_coordinates[2];
132
// P_FFC = J^dag (p - b), P_FIAT = 2*P_FFC - (1, 1)
133
double X = 2*(K[0]*(x[0] - b0) + K[1]*(x[1] - b1) + K[2]*(x[2] - b2)) - 1.0;
134
double Y = 2*(K[3]*(x[0] - b0) + K[4]*(x[1] - b1) + K[5]*(x[2] - b2)) - 1.0;
145
// Array of basisvalues
146
double basisvalues[3] = {0.0, 0.0, 0.0};
148
// Declare helper variables
149
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
151
// Compute basisvalues
152
basisvalues[0] = 1.0;
153
basisvalues[1] = tmp0;
154
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
155
basisvalues[0] *= std::sqrt(0.5);
156
basisvalues[2] *= std::sqrt(1.0);
157
basisvalues[1] *= std::sqrt(3.0);
159
// Table(s) of coefficients
160
static const double coefficients0[3] = \
161
{0.47140452, 0.28867513, -0.16666667};
163
static const double coefficients1[3] = \
164
{0.47140452, 0.0, 0.33333333};
167
for (unsigned int r = 0; r < 3; r++)
169
values[0] += coefficients0[r]*basisvalues[r];
170
values[1] += coefficients1[r]*basisvalues[r];
171
}// end loop over 'r'
173
// Using contravariant Piola transform to map values back to the physical element
174
const double tmp_ref0 = values[0];
175
const double tmp_ref1 = values[1];
176
values[0] = (1.0/detJ)*(J[0]*tmp_ref0 + J[1]*tmp_ref1);
177
values[1] = (1.0/detJ)*(J[2]*tmp_ref0 + J[3]*tmp_ref1);
178
values[2] = (1.0/detJ)*(J[4]*tmp_ref0 + J[5]*tmp_ref1);
184
// Array of basisvalues
185
double basisvalues[3] = {0.0, 0.0, 0.0};
187
// Declare helper variables
188
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
190
// Compute basisvalues
191
basisvalues[0] = 1.0;
192
basisvalues[1] = tmp0;
193
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
194
basisvalues[0] *= std::sqrt(0.5);
195
basisvalues[2] *= std::sqrt(1.0);
196
basisvalues[1] *= std::sqrt(3.0);
198
// Table(s) of coefficients
199
static const double coefficients0[3] = \
200
{0.94280904, -0.28867513, 0.16666667};
202
static const double coefficients1[3] = \
203
{-0.47140452, 0.0, -0.33333333};
206
for (unsigned int r = 0; r < 3; r++)
208
values[0] += coefficients0[r]*basisvalues[r];
209
values[1] += coefficients1[r]*basisvalues[r];
210
}// end loop over 'r'
212
// Using contravariant Piola transform to map values back to the physical element
213
const double tmp_ref0 = values[0];
214
const double tmp_ref1 = values[1];
215
values[0] = (1.0/detJ)*(J[0]*tmp_ref0 + J[1]*tmp_ref1);
216
values[1] = (1.0/detJ)*(J[2]*tmp_ref0 + J[3]*tmp_ref1);
217
values[2] = (1.0/detJ)*(J[4]*tmp_ref0 + J[5]*tmp_ref1);
223
// Array of basisvalues
224
double basisvalues[3] = {0.0, 0.0, 0.0};
226
// Declare helper variables
227
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
229
// Compute basisvalues
230
basisvalues[0] = 1.0;
231
basisvalues[1] = tmp0;
232
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
233
basisvalues[0] *= std::sqrt(0.5);
234
basisvalues[2] *= std::sqrt(1.0);
235
basisvalues[1] *= std::sqrt(3.0);
237
// Table(s) of coefficients
238
static const double coefficients0[3] = \
239
{0.47140452, 0.28867513, -0.16666667};
241
static const double coefficients1[3] = \
242
{-0.94280904, 0.0, 0.33333333};
245
for (unsigned int r = 0; r < 3; r++)
247
values[0] += coefficients0[r]*basisvalues[r];
248
values[1] += coefficients1[r]*basisvalues[r];
249
}// end loop over 'r'
251
// Using contravariant Piola transform to map values back to the physical element
252
const double tmp_ref0 = values[0];
253
const double tmp_ref1 = values[1];
254
values[0] = (1.0/detJ)*(J[0]*tmp_ref0 + J[1]*tmp_ref1);
255
values[1] = (1.0/detJ)*(J[2]*tmp_ref0 + J[3]*tmp_ref1);
256
values[2] = (1.0/detJ)*(J[4]*tmp_ref0 + J[5]*tmp_ref1);
263
/// Evaluate all basis functions at given point x in cell
264
virtual void evaluate_basis_all(double* values,
266
const double* vertex_coordinates,
267
int cell_orientation) const
269
// Helper variable to hold values of a single dof.
270
double dof_values[3] = {0.0, 0.0, 0.0};
272
// Loop dofs and call evaluate_basis
273
for (unsigned int r = 0; r < 3; r++)
275
evaluate_basis(r, dof_values, x, vertex_coordinates, cell_orientation);
276
for (unsigned int s = 0; s < 3; s++)
278
values[r*3 + s] = dof_values[s];
279
}// end loop over 's'
280
}// end loop over 'r'
283
/// Evaluate order n derivatives of basis function i at given point x in cell
284
virtual void evaluate_basis_derivatives(std::size_t i,
288
const double* vertex_coordinates,
289
int cell_orientation) const
293
compute_jacobian_triangle_3d(J, vertex_coordinates);
295
// Compute Jacobian inverse and determinant
298
compute_jacobian_inverse_triangle_3d(K, detJ, J);
302
if (cell_orientation == -1)
303
throw std::runtime_error("cell orientation must be defined (not -1)");
304
// (If cell_orientation == 1 = down, multiply det(J) by -1)
305
else if (cell_orientation == 1)
309
const double b0 = vertex_coordinates[0];
310
const double b1 = vertex_coordinates[1];
311
const double b2 = vertex_coordinates[2];
313
// P_FFC = J^dag (p - b), P_FIAT = 2*P_FFC - (1, 1)
314
double X = 2*(K[0]*(x[0] - b0) + K[1]*(x[1] - b1) + K[2]*(x[2] - b2)) - 1.0;
315
double Y = 2*(K[3]*(x[0] - b0) + K[4]*(x[1] - b1) + K[5]*(x[2] - b2)) - 1.0;
318
// Compute number of derivatives.
319
unsigned int num_derivatives_t = 1;
320
for (unsigned int r = 0; r < n; r++)
322
num_derivatives_t *= 2;
323
}// end loop over 'r'
325
// Compute number of derivatives.
326
unsigned int num_derivatives_g = 1;
327
for (unsigned int r = 0; r < n; r++)
329
num_derivatives_g *= 3;
330
}// end loop over 'r'
332
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
333
unsigned int **combinations_t = new unsigned int *[num_derivatives_t];
334
for (unsigned int row = 0; row < num_derivatives_t; row++)
336
combinations_t[row] = new unsigned int [n];
337
for (unsigned int col = 0; col < n; col++)
338
combinations_t[row][col] = 0;
341
// Generate combinations of derivatives
342
for (unsigned int row = 1; row < num_derivatives_t; row++)
344
for (unsigned int num = 0; num < row; num++)
346
for (unsigned int col = n-1; col+1 > 0; col--)
348
if (combinations_t[row][col] + 1 > 1)
349
combinations_t[row][col] = 0;
352
combinations_t[row][col] += 1;
359
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
360
unsigned int **combinations_g = new unsigned int *[num_derivatives_g];
361
for (unsigned int row = 0; row < num_derivatives_g; row++)
363
combinations_g[row] = new unsigned int [n];
364
for (unsigned int col = 0; col < n; col++)
365
combinations_g[row][col] = 0;
368
// Generate combinations of derivatives
369
for (unsigned int row = 1; row < num_derivatives_g; row++)
371
for (unsigned int num = 0; num < row; num++)
373
for (unsigned int col = n-1; col+1 > 0; col--)
375
if (combinations_g[row][col] + 1 > 2)
376
combinations_g[row][col] = 0;
379
combinations_g[row][col] += 1;
386
// Compute inverse of Jacobian
387
const double Jinv[2][3] = {{K[0], K[1], K[2]}, {K[3], K[4], K[5]}};
389
// Declare transformation matrix
390
// Declare pointer to two dimensional array and initialise
391
double **transform = new double *[num_derivatives_g];
393
for (unsigned int j = 0; j < num_derivatives_g; j++)
395
transform[j] = new double [num_derivatives_t];
396
for (unsigned int k = 0; k < num_derivatives_t; k++)
400
// Construct transformation matrix
401
for (unsigned int row = 0; row < num_derivatives_g; row++)
403
for (unsigned int col = 0; col < num_derivatives_t; col++)
405
for (unsigned int k = 0; k < n; k++)
406
transform[row][col] *= Jinv[combinations_t[col][k]][combinations_g[row][k]];
410
// Reset values. Assuming that values is always an array.
411
for (unsigned int r = 0; r < 3*num_derivatives_g; r++)
414
}// end loop over 'r'
421
// Array of basisvalues
422
double basisvalues[3] = {0.0, 0.0, 0.0};
424
// Declare helper variables
425
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
427
// Compute basisvalues
428
basisvalues[0] = 1.0;
429
basisvalues[1] = tmp0;
430
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
431
basisvalues[0] *= std::sqrt(0.5);
432
basisvalues[2] *= std::sqrt(1.0);
433
basisvalues[1] *= std::sqrt(3.0);
435
// Table(s) of coefficients
436
static const double coefficients0[3] = \
437
{0.47140452, 0.28867513, -0.16666667};
439
static const double coefficients1[3] = \
440
{0.47140452, 0.0, 0.33333333};
442
// Tables of derivatives of the polynomial base (transpose).
443
static const double dmats0[3][3] = \
445
{4.8989795, 0.0, 0.0},
448
static const double dmats1[3][3] = \
450
{2.4494897, 0.0, 0.0},
451
{4.2426407, 0.0, 0.0}};
453
// Compute reference derivatives.
454
// Declare pointer to array of derivatives on FIAT element.
455
double *derivatives = new double[2*num_derivatives_t];
456
for (unsigned int r = 0; r < 2*num_derivatives_t; r++)
458
derivatives[r] = 0.0;
459
}// end loop over 'r'
461
// Declare pointer to array of reference derivatives on physical element.
462
double *derivatives_p = new double[3*num_derivatives_t];
463
for (unsigned int r = 0; r < 3*num_derivatives_t; r++)
465
derivatives_p[r] = 0.0;
466
}// end loop over 'r'
468
// Declare derivative matrix (of polynomial basis).
469
double dmats[3][3] = \
474
// Declare (auxiliary) derivative matrix (of polynomial basis).
475
double dmats_old[3][3] = \
480
// Loop possible derivatives.
481
for (unsigned int r = 0; r < num_derivatives_t; r++)
483
// Resetting dmats values to compute next derivative.
484
for (unsigned int t = 0; t < 3; t++)
486
for (unsigned int u = 0; u < 3; u++)
494
}// end loop over 'u'
495
}// end loop over 't'
497
// Looping derivative order to generate dmats.
498
for (unsigned int s = 0; s < n; s++)
500
// Updating dmats_old with new values and resetting dmats.
501
for (unsigned int t = 0; t < 3; t++)
503
for (unsigned int u = 0; u < 3; u++)
505
dmats_old[t][u] = dmats[t][u];
507
}// end loop over 'u'
508
}// end loop over 't'
510
// Update dmats using an inner product.
511
if (combinations_t[r][s] == 0)
513
for (unsigned int t = 0; t < 3; t++)
515
for (unsigned int u = 0; u < 3; u++)
517
for (unsigned int tu = 0; tu < 3; tu++)
519
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
520
}// end loop over 'tu'
521
}// end loop over 'u'
522
}// end loop over 't'
525
if (combinations_t[r][s] == 1)
527
for (unsigned int t = 0; t < 3; t++)
529
for (unsigned int u = 0; u < 3; u++)
531
for (unsigned int tu = 0; tu < 3; tu++)
533
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
534
}// end loop over 'tu'
535
}// end loop over 'u'
536
}// end loop over 't'
539
}// end loop over 's'
540
for (unsigned int s = 0; s < 3; s++)
542
for (unsigned int t = 0; t < 3; t++)
544
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
545
derivatives[num_derivatives_t + r] += coefficients1[s]*dmats[s][t]*basisvalues[t];
546
}// end loop over 't'
547
}// end loop over 's'
549
// Using contravariant Piola transform to map values back to the physical element.
550
const double tmp_ref0 = derivatives[r];
551
const double tmp_ref1 = derivatives[num_derivatives_t + r];
552
derivatives_p[r] = (1.0/detJ)*(J[0]*tmp_ref0 + J[1]*tmp_ref1);
553
derivatives_p[num_derivatives_t + r] = (1.0/detJ)*(J[2]*tmp_ref0 + J[3]*tmp_ref1);
554
derivatives_p[2*num_derivatives_t + r] = (1.0/detJ)*(J[4]*tmp_ref0 + J[5]*tmp_ref1);
555
}// end loop over 'r'
557
// Transform derivatives back to physical element
558
for (unsigned int r = 0; r < num_derivatives_g; r++)
560
for (unsigned int s = 0; s < num_derivatives_t; s++)
562
values[r] += transform[r][s]*derivatives_p[s];
563
values[num_derivatives_g + r] += transform[r][s]*derivatives_p[num_derivatives_t + s];
564
values[2*num_derivatives_g + r] += transform[r][s]*derivatives_p[2*num_derivatives_t + s];
565
}// end loop over 's'
566
}// end loop over 'r'
568
// Delete pointer to array of derivatives on FIAT element
569
delete [] derivatives;
571
// Delete pointer to array of reference derivatives on physical element.
572
delete [] derivatives_p;
574
// Delete pointer to array of combinations of derivatives and transform
575
for (unsigned int r = 0; r < num_derivatives_t; r++)
577
delete [] combinations_t[r];
578
}// end loop over 'r'
579
delete [] combinations_t;
580
for (unsigned int r = 0; r < num_derivatives_g; r++)
582
delete [] combinations_g[r];
583
}// end loop over 'r'
584
delete [] combinations_g;
585
for (unsigned int r = 0; r < num_derivatives_g; r++)
587
delete [] transform[r];
588
}// end loop over 'r'
595
// Array of basisvalues
596
double basisvalues[3] = {0.0, 0.0, 0.0};
598
// Declare helper variables
599
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
601
// Compute basisvalues
602
basisvalues[0] = 1.0;
603
basisvalues[1] = tmp0;
604
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
605
basisvalues[0] *= std::sqrt(0.5);
606
basisvalues[2] *= std::sqrt(1.0);
607
basisvalues[1] *= std::sqrt(3.0);
609
// Table(s) of coefficients
610
static const double coefficients0[3] = \
611
{0.94280904, -0.28867513, 0.16666667};
613
static const double coefficients1[3] = \
614
{-0.47140452, 0.0, -0.33333333};
616
// Tables of derivatives of the polynomial base (transpose).
617
static const double dmats0[3][3] = \
619
{4.8989795, 0.0, 0.0},
622
static const double dmats1[3][3] = \
624
{2.4494897, 0.0, 0.0},
625
{4.2426407, 0.0, 0.0}};
627
// Compute reference derivatives.
628
// Declare pointer to array of derivatives on FIAT element.
629
double *derivatives = new double[2*num_derivatives_t];
630
for (unsigned int r = 0; r < 2*num_derivatives_t; r++)
632
derivatives[r] = 0.0;
633
}// end loop over 'r'
635
// Declare pointer to array of reference derivatives on physical element.
636
double *derivatives_p = new double[3*num_derivatives_t];
637
for (unsigned int r = 0; r < 3*num_derivatives_t; r++)
639
derivatives_p[r] = 0.0;
640
}// end loop over 'r'
642
// Declare derivative matrix (of polynomial basis).
643
double dmats[3][3] = \
648
// Declare (auxiliary) derivative matrix (of polynomial basis).
649
double dmats_old[3][3] = \
654
// Loop possible derivatives.
655
for (unsigned int r = 0; r < num_derivatives_t; r++)
657
// Resetting dmats values to compute next derivative.
658
for (unsigned int t = 0; t < 3; t++)
660
for (unsigned int u = 0; u < 3; u++)
668
}// end loop over 'u'
669
}// end loop over 't'
671
// Looping derivative order to generate dmats.
672
for (unsigned int s = 0; s < n; s++)
674
// Updating dmats_old with new values and resetting dmats.
675
for (unsigned int t = 0; t < 3; t++)
677
for (unsigned int u = 0; u < 3; u++)
679
dmats_old[t][u] = dmats[t][u];
681
}// end loop over 'u'
682
}// end loop over 't'
684
// Update dmats using an inner product.
685
if (combinations_t[r][s] == 0)
687
for (unsigned int t = 0; t < 3; t++)
689
for (unsigned int u = 0; u < 3; u++)
691
for (unsigned int tu = 0; tu < 3; tu++)
693
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
694
}// end loop over 'tu'
695
}// end loop over 'u'
696
}// end loop over 't'
699
if (combinations_t[r][s] == 1)
701
for (unsigned int t = 0; t < 3; t++)
703
for (unsigned int u = 0; u < 3; u++)
705
for (unsigned int tu = 0; tu < 3; tu++)
707
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
708
}// end loop over 'tu'
709
}// end loop over 'u'
710
}// end loop over 't'
713
}// end loop over 's'
714
for (unsigned int s = 0; s < 3; s++)
716
for (unsigned int t = 0; t < 3; t++)
718
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
719
derivatives[num_derivatives_t + r] += coefficients1[s]*dmats[s][t]*basisvalues[t];
720
}// end loop over 't'
721
}// end loop over 's'
723
// Using contravariant Piola transform to map values back to the physical element.
724
const double tmp_ref0 = derivatives[r];
725
const double tmp_ref1 = derivatives[num_derivatives_t + r];
726
derivatives_p[r] = (1.0/detJ)*(J[0]*tmp_ref0 + J[1]*tmp_ref1);
727
derivatives_p[num_derivatives_t + r] = (1.0/detJ)*(J[2]*tmp_ref0 + J[3]*tmp_ref1);
728
derivatives_p[2*num_derivatives_t + r] = (1.0/detJ)*(J[4]*tmp_ref0 + J[5]*tmp_ref1);
729
}// end loop over 'r'
731
// Transform derivatives back to physical element
732
for (unsigned int r = 0; r < num_derivatives_g; r++)
734
for (unsigned int s = 0; s < num_derivatives_t; s++)
736
values[r] += transform[r][s]*derivatives_p[s];
737
values[num_derivatives_g + r] += transform[r][s]*derivatives_p[num_derivatives_t + s];
738
values[2*num_derivatives_g + r] += transform[r][s]*derivatives_p[2*num_derivatives_t + s];
739
}// end loop over 's'
740
}// end loop over 'r'
742
// Delete pointer to array of derivatives on FIAT element
743
delete [] derivatives;
745
// Delete pointer to array of reference derivatives on physical element.
746
delete [] derivatives_p;
748
// Delete pointer to array of combinations of derivatives and transform
749
for (unsigned int r = 0; r < num_derivatives_t; r++)
751
delete [] combinations_t[r];
752
}// end loop over 'r'
753
delete [] combinations_t;
754
for (unsigned int r = 0; r < num_derivatives_g; r++)
756
delete [] combinations_g[r];
757
}// end loop over 'r'
758
delete [] combinations_g;
759
for (unsigned int r = 0; r < num_derivatives_g; r++)
761
delete [] transform[r];
762
}// end loop over 'r'
769
// Array of basisvalues
770
double basisvalues[3] = {0.0, 0.0, 0.0};
772
// Declare helper variables
773
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
775
// Compute basisvalues
776
basisvalues[0] = 1.0;
777
basisvalues[1] = tmp0;
778
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
779
basisvalues[0] *= std::sqrt(0.5);
780
basisvalues[2] *= std::sqrt(1.0);
781
basisvalues[1] *= std::sqrt(3.0);
783
// Table(s) of coefficients
784
static const double coefficients0[3] = \
785
{0.47140452, 0.28867513, -0.16666667};
787
static const double coefficients1[3] = \
788
{-0.94280904, 0.0, 0.33333333};
790
// Tables of derivatives of the polynomial base (transpose).
791
static const double dmats0[3][3] = \
793
{4.8989795, 0.0, 0.0},
796
static const double dmats1[3][3] = \
798
{2.4494897, 0.0, 0.0},
799
{4.2426407, 0.0, 0.0}};
801
// Compute reference derivatives.
802
// Declare pointer to array of derivatives on FIAT element.
803
double *derivatives = new double[2*num_derivatives_t];
804
for (unsigned int r = 0; r < 2*num_derivatives_t; r++)
806
derivatives[r] = 0.0;
807
}// end loop over 'r'
809
// Declare pointer to array of reference derivatives on physical element.
810
double *derivatives_p = new double[3*num_derivatives_t];
811
for (unsigned int r = 0; r < 3*num_derivatives_t; r++)
813
derivatives_p[r] = 0.0;
814
}// end loop over 'r'
816
// Declare derivative matrix (of polynomial basis).
817
double dmats[3][3] = \
822
// Declare (auxiliary) derivative matrix (of polynomial basis).
823
double dmats_old[3][3] = \
828
// Loop possible derivatives.
829
for (unsigned int r = 0; r < num_derivatives_t; r++)
831
// Resetting dmats values to compute next derivative.
832
for (unsigned int t = 0; t < 3; t++)
834
for (unsigned int u = 0; u < 3; u++)
842
}// end loop over 'u'
843
}// end loop over 't'
845
// Looping derivative order to generate dmats.
846
for (unsigned int s = 0; s < n; s++)
848
// Updating dmats_old with new values and resetting dmats.
849
for (unsigned int t = 0; t < 3; t++)
851
for (unsigned int u = 0; u < 3; u++)
853
dmats_old[t][u] = dmats[t][u];
855
}// end loop over 'u'
856
}// end loop over 't'
858
// Update dmats using an inner product.
859
if (combinations_t[r][s] == 0)
861
for (unsigned int t = 0; t < 3; t++)
863
for (unsigned int u = 0; u < 3; u++)
865
for (unsigned int tu = 0; tu < 3; tu++)
867
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
868
}// end loop over 'tu'
869
}// end loop over 'u'
870
}// end loop over 't'
873
if (combinations_t[r][s] == 1)
875
for (unsigned int t = 0; t < 3; t++)
877
for (unsigned int u = 0; u < 3; u++)
879
for (unsigned int tu = 0; tu < 3; tu++)
881
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
882
}// end loop over 'tu'
883
}// end loop over 'u'
884
}// end loop over 't'
887
}// end loop over 's'
888
for (unsigned int s = 0; s < 3; s++)
890
for (unsigned int t = 0; t < 3; t++)
892
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
893
derivatives[num_derivatives_t + r] += coefficients1[s]*dmats[s][t]*basisvalues[t];
894
}// end loop over 't'
895
}// end loop over 's'
897
// Using contravariant Piola transform to map values back to the physical element.
898
const double tmp_ref0 = derivatives[r];
899
const double tmp_ref1 = derivatives[num_derivatives_t + r];
900
derivatives_p[r] = (1.0/detJ)*(J[0]*tmp_ref0 + J[1]*tmp_ref1);
901
derivatives_p[num_derivatives_t + r] = (1.0/detJ)*(J[2]*tmp_ref0 + J[3]*tmp_ref1);
902
derivatives_p[2*num_derivatives_t + r] = (1.0/detJ)*(J[4]*tmp_ref0 + J[5]*tmp_ref1);
903
}// end loop over 'r'
905
// Transform derivatives back to physical element
906
for (unsigned int r = 0; r < num_derivatives_g; r++)
908
for (unsigned int s = 0; s < num_derivatives_t; s++)
910
values[r] += transform[r][s]*derivatives_p[s];
911
values[num_derivatives_g + r] += transform[r][s]*derivatives_p[num_derivatives_t + s];
912
values[2*num_derivatives_g + r] += transform[r][s]*derivatives_p[2*num_derivatives_t + s];
913
}// end loop over 's'
914
}// end loop over 'r'
916
// Delete pointer to array of derivatives on FIAT element
917
delete [] derivatives;
919
// Delete pointer to array of reference derivatives on physical element.
920
delete [] derivatives_p;
922
// Delete pointer to array of combinations of derivatives and transform
923
for (unsigned int r = 0; r < num_derivatives_t; r++)
925
delete [] combinations_t[r];
926
}// end loop over 'r'
927
delete [] combinations_t;
928
for (unsigned int r = 0; r < num_derivatives_g; r++)
930
delete [] combinations_g[r];
931
}// end loop over 'r'
932
delete [] combinations_g;
933
for (unsigned int r = 0; r < num_derivatives_g; r++)
935
delete [] transform[r];
936
}// end loop over 'r'
944
/// Evaluate order n derivatives of all basis functions at given point x in cell
945
virtual void evaluate_basis_derivatives_all(std::size_t n,
948
const double* vertex_coordinates,
949
int cell_orientation) const
951
// Compute number of derivatives.
952
unsigned int num_derivatives_g = 1;
953
for (unsigned int r = 0; r < n; r++)
955
num_derivatives_g *= 3;
956
}// end loop over 'r'
958
// Helper variable to hold values of a single dof.
959
double *dof_values = new double[3*num_derivatives_g];
960
for (unsigned int r = 0; r < 3*num_derivatives_g; r++)
963
}// end loop over 'r'
965
// Loop dofs and call evaluate_basis_derivatives.
966
for (unsigned int r = 0; r < 3; r++)
968
evaluate_basis_derivatives(r, n, dof_values, x, vertex_coordinates, cell_orientation);
969
for (unsigned int s = 0; s < 3*num_derivatives_g; s++)
971
values[r*3*num_derivatives_g + s] = dof_values[s];
972
}// end loop over 's'
973
}// end loop over 'r'
976
delete [] dof_values;
979
/// Evaluate linear functional for dof i on the function f
980
virtual double evaluate_dof(std::size_t i,
981
const ufc::function& f,
982
const double* vertex_coordinates,
983
int cell_orientation,
984
const ufc::cell& c) const
986
// Declare variables for result of evaluation
989
// Declare variable for physical coordinates
995
compute_jacobian_triangle_3d(J, vertex_coordinates);
998
// Compute Jacobian inverse and determinant
1001
compute_jacobian_inverse_triangle_3d(K, detJ, J);
1005
// Check orientation
1006
if (cell_orientation == -1)
1007
throw std::runtime_error("cell orientation must be defined (not -1)");
1008
// (If cell_orientation == 1 = down, multiply det(J) by -1)
1009
else if (cell_orientation == 1)
1015
y[0] = 0.5*vertex_coordinates[3] + 0.5*vertex_coordinates[6];
1016
y[1] = 0.5*vertex_coordinates[4] + 0.5*vertex_coordinates[7];
1017
y[2] = 0.5*vertex_coordinates[5] + 0.5*vertex_coordinates[8];
1018
f.evaluate(vals, y, c);
1019
result = (detJ*(K[0]*vals[0] + K[1]*vals[1] + K[2]*vals[2])) + (detJ*(K[3]*vals[0] + K[4]*vals[1] + K[5]*vals[2]));
1025
y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[6];
1026
y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[7];
1027
y[2] = 0.5*vertex_coordinates[2] + 0.5*vertex_coordinates[8];
1028
f.evaluate(vals, y, c);
1029
result = (detJ*(K[0]*vals[0] + K[1]*vals[1] + K[2]*vals[2]));
1035
y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[3];
1036
y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[4];
1037
y[2] = 0.5*vertex_coordinates[2] + 0.5*vertex_coordinates[5];
1038
f.evaluate(vals, y, c);
1039
result = (-1.0)*(detJ*(K[3]*vals[0] + K[4]*vals[1] + K[5]*vals[2]));
1048
/// Evaluate linear functionals for all dofs on the function f
1049
virtual void evaluate_dofs(double* values,
1050
const ufc::function& f,
1051
const double* vertex_coordinates,
1052
int cell_orientation,
1053
const ufc::cell& c) const
1055
// Declare variables for result of evaluation
1058
// Declare variable for physical coordinates
1064
compute_jacobian_triangle_3d(J, vertex_coordinates);
1067
// Compute Jacobian inverse and determinant
1070
compute_jacobian_inverse_triangle_3d(K, detJ, J);
1074
// Check orientation
1075
if (cell_orientation == -1)
1076
throw std::runtime_error("cell orientation must be defined (not -1)");
1077
// (If cell_orientation == 1 = down, multiply det(J) by -1)
1078
else if (cell_orientation == 1)
1080
y[0] = 0.5*vertex_coordinates[3] + 0.5*vertex_coordinates[6];
1081
y[1] = 0.5*vertex_coordinates[4] + 0.5*vertex_coordinates[7];
1082
y[2] = 0.5*vertex_coordinates[5] + 0.5*vertex_coordinates[8];
1083
f.evaluate(vals, y, c);
1084
result = (detJ*(K[0]*vals[0] + K[1]*vals[1] + K[2]*vals[2])) + (detJ*(K[3]*vals[0] + K[4]*vals[1] + K[5]*vals[2]));
1086
y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[6];
1087
y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[7];
1088
y[2] = 0.5*vertex_coordinates[2] + 0.5*vertex_coordinates[8];
1089
f.evaluate(vals, y, c);
1090
result = (detJ*(K[0]*vals[0] + K[1]*vals[1] + K[2]*vals[2]));
1092
y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[3];
1093
y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[4];
1094
y[2] = 0.5*vertex_coordinates[2] + 0.5*vertex_coordinates[5];
1095
f.evaluate(vals, y, c);
1096
result = (-1.0)*(detJ*(K[3]*vals[0] + K[4]*vals[1] + K[5]*vals[2]));
1100
/// Interpolate vertex values from dof values
1101
virtual void interpolate_vertex_values(double* vertex_values,
1102
const double* dof_values,
1103
const double* vertex_coordinates,
1104
int cell_orientation,
1105
const ufc::cell& c) const
1109
compute_jacobian_triangle_3d(J, vertex_coordinates);
1112
// Compute Jacobian inverse and determinant
1115
compute_jacobian_inverse_triangle_3d(K, detJ, J);
1119
// Check orientation
1120
if (cell_orientation == -1)
1121
throw std::runtime_error("cell orientation must be defined (not -1)");
1122
// (If cell_orientation == 1 = down, multiply det(J) by -1)
1123
else if (cell_orientation == 1)
1126
// Evaluate function and change variables
1127
vertex_values[0] = dof_values[1]*(1.0/detJ)*J[0] + dof_values[2]*((1.0/detJ)*(J[1]*(-1.0)));
1128
vertex_values[2] = dof_values[0]*(1.0/detJ)*J[0] + dof_values[2]*((1.0/detJ)*(J[0] + J[1]*(-1.0)));
1129
vertex_values[4] = dof_values[0]*(1.0/detJ)*J[1] + dof_values[1]*((1.0/detJ)*(J[0] + J[1]*(-1.0)));
1130
vertex_values[1] = dof_values[1]*(1.0/detJ)*J[2] + dof_values[2]*((1.0/detJ)*(J[3]*(-1.0)));
1131
vertex_values[3] = dof_values[0]*(1.0/detJ)*J[2] + dof_values[2]*((1.0/detJ)*(J[2] + J[3]*(-1.0)));
1132
vertex_values[5] = dof_values[0]*(1.0/detJ)*J[3] + dof_values[1]*((1.0/detJ)*(J[2] + J[3]*(-1.0)));
1135
/// Map coordinate xhat from reference cell to coordinate x in cell
1136
virtual void map_from_reference_cell(double* x,
1138
const ufc::cell& c) const
1140
std::cerr << "*** FFC warning: " << "map_from_reference_cell not yet implemented." << std::endl;
1143
/// Map from coordinate x in cell to coordinate xhat in reference cell
1144
virtual void map_to_reference_cell(double* xhat,
1146
const ufc::cell& c) const
1148
std::cerr << "*** FFC warning: " << "map_to_reference_cell not yet implemented." << std::endl;
1151
/// Return the number of sub elements (for a mixed element)
1152
virtual std::size_t num_sub_elements() const
1157
/// Create a new finite element for sub element i (for a mixed element)
1158
virtual ufc::finite_element* create_sub_element(std::size_t i) const
1163
/// Create a new class instance
1164
virtual ufc::finite_element* create() const
1166
return new projectionmanifold_finite_element_0();
1171
/// This class defines the interface for a finite element.
1173
class projectionmanifold_finite_element_1: public ufc::finite_element
1178
projectionmanifold_finite_element_1() : ufc::finite_element()
1184
virtual ~projectionmanifold_finite_element_1()
1189
/// Return a string identifying the finite element
1190
virtual const char* signature() const
1192
return "FiniteElement('Discontinuous Lagrange', Domain(Cell('triangle', 3), 'triangle_multiverse', 3, 2), 0, None)";
1195
/// Return the cell shape
1196
virtual ufc::shape cell_shape() const
1198
return ufc::triangle;
1201
/// Return the topological dimension of the cell shape
1202
virtual std::size_t topological_dimension() const
1207
/// Return the geometric dimension of the cell shape
1208
virtual std::size_t geometric_dimension() const
1213
/// Return the dimension of the finite element function space
1214
virtual std::size_t space_dimension() const
1219
/// Return the rank of the value space
1220
virtual std::size_t value_rank() const
1225
/// Return the dimension of the value space for axis i
1226
virtual std::size_t value_dimension(std::size_t i) const
1231
/// Evaluate basis function i at given point x in cell
1232
virtual void evaluate_basis(std::size_t i,
1235
const double* vertex_coordinates,
1236
int cell_orientation) const
1240
compute_jacobian_triangle_3d(J, vertex_coordinates);
1242
// Compute Jacobian inverse and determinant
1245
compute_jacobian_inverse_triangle_3d(K, detJ, J);
1249
// P_FFC = J^dag (p - b), P_FIAT = 2*P_FFC - (1, 1)
1255
// Array of basisvalues
1256
double basisvalues[1] = {0.0};
1258
// Declare helper variables
1260
// Compute basisvalues
1261
basisvalues[0] = 1.0;
1263
// Table(s) of coefficients
1264
static const double coefficients0[1] = \
1268
for (unsigned int r = 0; r < 1; r++)
1270
*values += coefficients0[r]*basisvalues[r];
1271
}// end loop over 'r'
1274
/// Evaluate all basis functions at given point x in cell
1275
virtual void evaluate_basis_all(double* values,
1277
const double* vertex_coordinates,
1278
int cell_orientation) const
1280
// Element is constant, calling evaluate_basis.
1281
evaluate_basis(0, values, x, vertex_coordinates, cell_orientation);
1284
/// Evaluate order n derivatives of basis function i at given point x in cell
1285
virtual void evaluate_basis_derivatives(std::size_t i,
1289
const double* vertex_coordinates,
1290
int cell_orientation) const
1294
compute_jacobian_triangle_3d(J, vertex_coordinates);
1296
// Compute Jacobian inverse and determinant
1299
compute_jacobian_inverse_triangle_3d(K, detJ, J);
1303
// P_FFC = J^dag (p - b), P_FIAT = 2*P_FFC - (1, 1)
1306
// Compute number of derivatives.
1307
unsigned int num_derivatives_t = 1;
1308
for (unsigned int r = 0; r < n; r++)
1310
num_derivatives_t *= 2;
1311
}// end loop over 'r'
1313
// Compute number of derivatives.
1314
unsigned int num_derivatives_g = 1;
1315
for (unsigned int r = 0; r < n; r++)
1317
num_derivatives_g *= 3;
1318
}// end loop over 'r'
1320
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
1321
unsigned int **combinations_t = new unsigned int *[num_derivatives_t];
1322
for (unsigned int row = 0; row < num_derivatives_t; row++)
1324
combinations_t[row] = new unsigned int [n];
1325
for (unsigned int col = 0; col < n; col++)
1326
combinations_t[row][col] = 0;
1329
// Generate combinations of derivatives
1330
for (unsigned int row = 1; row < num_derivatives_t; row++)
1332
for (unsigned int num = 0; num < row; num++)
1334
for (unsigned int col = n-1; col+1 > 0; col--)
1336
if (combinations_t[row][col] + 1 > 1)
1337
combinations_t[row][col] = 0;
1340
combinations_t[row][col] += 1;
1347
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
1348
unsigned int **combinations_g = new unsigned int *[num_derivatives_g];
1349
for (unsigned int row = 0; row < num_derivatives_g; row++)
1351
combinations_g[row] = new unsigned int [n];
1352
for (unsigned int col = 0; col < n; col++)
1353
combinations_g[row][col] = 0;
1356
// Generate combinations of derivatives
1357
for (unsigned int row = 1; row < num_derivatives_g; row++)
1359
for (unsigned int num = 0; num < row; num++)
1361
for (unsigned int col = n-1; col+1 > 0; col--)
1363
if (combinations_g[row][col] + 1 > 2)
1364
combinations_g[row][col] = 0;
1367
combinations_g[row][col] += 1;
1374
// Compute inverse of Jacobian
1375
const double Jinv[2][3] = {{K[0], K[1], K[2]}, {K[3], K[4], K[5]}};
1377
// Declare transformation matrix
1378
// Declare pointer to two dimensional array and initialise
1379
double **transform = new double *[num_derivatives_g];
1381
for (unsigned int j = 0; j < num_derivatives_g; j++)
1383
transform[j] = new double [num_derivatives_t];
1384
for (unsigned int k = 0; k < num_derivatives_t; k++)
1385
transform[j][k] = 1;
1388
// Construct transformation matrix
1389
for (unsigned int row = 0; row < num_derivatives_g; row++)
1391
for (unsigned int col = 0; col < num_derivatives_t; col++)
1393
for (unsigned int k = 0; k < n; k++)
1394
transform[row][col] *= Jinv[combinations_t[col][k]][combinations_g[row][k]];
1398
// Reset values. Assuming that values is always an array.
1399
for (unsigned int r = 0; r < num_derivatives_g; r++)
1402
}// end loop over 'r'
1405
// Array of basisvalues
1406
double basisvalues[1] = {0.0};
1408
// Declare helper variables
1410
// Compute basisvalues
1411
basisvalues[0] = 1.0;
1413
// Table(s) of coefficients
1414
static const double coefficients0[1] = \
1417
// Tables of derivatives of the polynomial base (transpose).
1418
static const double dmats0[1][1] = \
1421
static const double dmats1[1][1] = \
1424
// Compute reference derivatives.
1425
// Declare pointer to array of derivatives on FIAT element.
1426
double *derivatives = new double[num_derivatives_t];
1427
for (unsigned int r = 0; r < num_derivatives_t; r++)
1429
derivatives[r] = 0.0;
1430
}// end loop over 'r'
1432
// Declare derivative matrix (of polynomial basis).
1433
double dmats[1][1] = \
1436
// Declare (auxiliary) derivative matrix (of polynomial basis).
1437
double dmats_old[1][1] = \
1440
// Loop possible derivatives.
1441
for (unsigned int r = 0; r < num_derivatives_t; r++)
1443
// Resetting dmats values to compute next derivative.
1444
for (unsigned int t = 0; t < 1; t++)
1446
for (unsigned int u = 0; u < 1; u++)
1454
}// end loop over 'u'
1455
}// end loop over 't'
1457
// Looping derivative order to generate dmats.
1458
for (unsigned int s = 0; s < n; s++)
1460
// Updating dmats_old with new values and resetting dmats.
1461
for (unsigned int t = 0; t < 1; t++)
1463
for (unsigned int u = 0; u < 1; u++)
1465
dmats_old[t][u] = dmats[t][u];
1467
}// end loop over 'u'
1468
}// end loop over 't'
1470
// Update dmats using an inner product.
1471
if (combinations_t[r][s] == 0)
1473
for (unsigned int t = 0; t < 1; t++)
1475
for (unsigned int u = 0; u < 1; u++)
1477
for (unsigned int tu = 0; tu < 1; tu++)
1479
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1480
}// end loop over 'tu'
1481
}// end loop over 'u'
1482
}// end loop over 't'
1485
if (combinations_t[r][s] == 1)
1487
for (unsigned int t = 0; t < 1; t++)
1489
for (unsigned int u = 0; u < 1; u++)
1491
for (unsigned int tu = 0; tu < 1; tu++)
1493
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1494
}// end loop over 'tu'
1495
}// end loop over 'u'
1496
}// end loop over 't'
1499
}// end loop over 's'
1500
for (unsigned int s = 0; s < 1; s++)
1502
for (unsigned int t = 0; t < 1; t++)
1504
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1505
}// end loop over 't'
1506
}// end loop over 's'
1507
}// end loop over 'r'
1509
// Transform derivatives back to physical element
1510
for (unsigned int r = 0; r < num_derivatives_g; r++)
1512
for (unsigned int s = 0; s < num_derivatives_t; s++)
1514
values[r] += transform[r][s]*derivatives[s];
1515
}// end loop over 's'
1516
}// end loop over 'r'
1518
// Delete pointer to array of derivatives on FIAT element
1519
delete [] derivatives;
1521
// Delete pointer to array of combinations of derivatives and transform
1522
for (unsigned int r = 0; r < num_derivatives_t; r++)
1524
delete [] combinations_t[r];
1525
}// end loop over 'r'
1526
delete [] combinations_t;
1527
for (unsigned int r = 0; r < num_derivatives_g; r++)
1529
delete [] combinations_g[r];
1530
}// end loop over 'r'
1531
delete [] combinations_g;
1532
for (unsigned int r = 0; r < num_derivatives_g; r++)
1534
delete [] transform[r];
1535
}// end loop over 'r'
1536
delete [] transform;
1539
/// Evaluate order n derivatives of all basis functions at given point x in cell
1540
virtual void evaluate_basis_derivatives_all(std::size_t n,
1543
const double* vertex_coordinates,
1544
int cell_orientation) const
1546
// Element is constant, calling evaluate_basis_derivatives.
1547
evaluate_basis_derivatives(0, n, values, x, vertex_coordinates, cell_orientation);
1550
/// Evaluate linear functional for dof i on the function f
1551
virtual double evaluate_dof(std::size_t i,
1552
const ufc::function& f,
1553
const double* vertex_coordinates,
1554
int cell_orientation,
1555
const ufc::cell& c) const
1557
// Declare variables for result of evaluation
1560
// Declare variable for physical coordinates
1566
y[0] = 0.33333333*vertex_coordinates[0] + 0.33333333*vertex_coordinates[3] + 0.33333333*vertex_coordinates[6];
1567
y[1] = 0.33333333*vertex_coordinates[1] + 0.33333333*vertex_coordinates[4] + 0.33333333*vertex_coordinates[7];
1568
y[2] = 0.33333333*vertex_coordinates[2] + 0.33333333*vertex_coordinates[5] + 0.33333333*vertex_coordinates[8];
1569
f.evaluate(vals, y, c);
1578
/// Evaluate linear functionals for all dofs on the function f
1579
virtual void evaluate_dofs(double* values,
1580
const ufc::function& f,
1581
const double* vertex_coordinates,
1582
int cell_orientation,
1583
const ufc::cell& c) const
1585
// Declare variables for result of evaluation
1588
// Declare variable for physical coordinates
1590
y[0] = 0.33333333*vertex_coordinates[0] + 0.33333333*vertex_coordinates[3] + 0.33333333*vertex_coordinates[6];
1591
y[1] = 0.33333333*vertex_coordinates[1] + 0.33333333*vertex_coordinates[4] + 0.33333333*vertex_coordinates[7];
1592
y[2] = 0.33333333*vertex_coordinates[2] + 0.33333333*vertex_coordinates[5] + 0.33333333*vertex_coordinates[8];
1593
f.evaluate(vals, y, c);
1594
values[0] = vals[0];
1597
/// Interpolate vertex values from dof values
1598
virtual void interpolate_vertex_values(double* vertex_values,
1599
const double* dof_values,
1600
const double* vertex_coordinates,
1601
int cell_orientation,
1602
const ufc::cell& c) const
1604
// Evaluate function and change variables
1605
vertex_values[0] = dof_values[0];
1606
vertex_values[1] = dof_values[0];
1607
vertex_values[2] = dof_values[0];
1610
/// Map coordinate xhat from reference cell to coordinate x in cell
1611
virtual void map_from_reference_cell(double* x,
1613
const ufc::cell& c) const
1615
std::cerr << "*** FFC warning: " << "map_from_reference_cell not yet implemented." << std::endl;
1618
/// Map from coordinate x in cell to coordinate xhat in reference cell
1619
virtual void map_to_reference_cell(double* xhat,
1621
const ufc::cell& c) const
1623
std::cerr << "*** FFC warning: " << "map_to_reference_cell not yet implemented." << std::endl;
1626
/// Return the number of sub elements (for a mixed element)
1627
virtual std::size_t num_sub_elements() const
1632
/// Create a new finite element for sub element i (for a mixed element)
1633
virtual ufc::finite_element* create_sub_element(std::size_t i) const
1638
/// Create a new class instance
1639
virtual ufc::finite_element* create() const
1641
return new projectionmanifold_finite_element_1();
1646
/// This class defines the interface for a finite element.
1648
class projectionmanifold_finite_element_2: public ufc::finite_element
1653
projectionmanifold_finite_element_2() : ufc::finite_element()
1659
virtual ~projectionmanifold_finite_element_2()
1664
/// Return a string identifying the finite element
1665
virtual const char* signature() const
1667
return "MixedElement(*[FiniteElement('Raviart-Thomas', Domain(Cell('triangle', 3), 'triangle_multiverse', 3, 2), 1, None), FiniteElement('Discontinuous Lagrange', Domain(Cell('triangle', 3), 'triangle_multiverse', 3, 2), 0, None)], **{'value_shape': (4,) })";
1670
/// Return the cell shape
1671
virtual ufc::shape cell_shape() const
1673
return ufc::triangle;
1676
/// Return the topological dimension of the cell shape
1677
virtual std::size_t topological_dimension() const
1682
/// Return the geometric dimension of the cell shape
1683
virtual std::size_t geometric_dimension() const
1688
/// Return the dimension of the finite element function space
1689
virtual std::size_t space_dimension() const
1694
/// Return the rank of the value space
1695
virtual std::size_t value_rank() const
1700
/// Return the dimension of the value space for axis i
1701
virtual std::size_t value_dimension(std::size_t i) const
1715
/// Evaluate basis function i at given point x in cell
1716
virtual void evaluate_basis(std::size_t i,
1719
const double* vertex_coordinates,
1720
int cell_orientation) const
1724
compute_jacobian_triangle_3d(J, vertex_coordinates);
1726
// Compute Jacobian inverse and determinant
1729
compute_jacobian_inverse_triangle_3d(K, detJ, J);
1732
// Check orientation
1733
if (cell_orientation == -1)
1734
throw std::runtime_error("cell orientation must be defined (not -1)");
1735
// (If cell_orientation == 1 = down, multiply det(J) by -1)
1736
else if (cell_orientation == 1)
1740
const double b0 = vertex_coordinates[0];
1741
const double b1 = vertex_coordinates[1];
1742
const double b2 = vertex_coordinates[2];
1744
// P_FFC = J^dag (p - b), P_FIAT = 2*P_FFC - (1, 1)
1745
double X = 2*(K[0]*(x[0] - b0) + K[1]*(x[1] - b1) + K[2]*(x[2] - b2)) - 1.0;
1746
double Y = 2*(K[3]*(x[0] - b0) + K[4]*(x[1] - b1) + K[5]*(x[2] - b2)) - 1.0;
1758
// Array of basisvalues
1759
double basisvalues[3] = {0.0, 0.0, 0.0};
1761
// Declare helper variables
1762
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1764
// Compute basisvalues
1765
basisvalues[0] = 1.0;
1766
basisvalues[1] = tmp0;
1767
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1768
basisvalues[0] *= std::sqrt(0.5);
1769
basisvalues[2] *= std::sqrt(1.0);
1770
basisvalues[1] *= std::sqrt(3.0);
1772
// Table(s) of coefficients
1773
static const double coefficients0[3] = \
1774
{0.47140452, 0.28867513, -0.16666667};
1776
static const double coefficients1[3] = \
1777
{0.47140452, 0.0, 0.33333333};
1780
for (unsigned int r = 0; r < 3; r++)
1782
values[0] += coefficients0[r]*basisvalues[r];
1783
values[1] += coefficients1[r]*basisvalues[r];
1784
}// end loop over 'r'
1786
// Using contravariant Piola transform to map values back to the physical element
1787
const double tmp_ref0 = values[0];
1788
const double tmp_ref1 = values[1];
1789
values[0] = (1.0/detJ)*(J[0]*tmp_ref0 + J[1]*tmp_ref1);
1790
values[1] = (1.0/detJ)*(J[2]*tmp_ref0 + J[3]*tmp_ref1);
1791
values[2] = (1.0/detJ)*(J[4]*tmp_ref0 + J[5]*tmp_ref1);
1797
// Array of basisvalues
1798
double basisvalues[3] = {0.0, 0.0, 0.0};
1800
// Declare helper variables
1801
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1803
// Compute basisvalues
1804
basisvalues[0] = 1.0;
1805
basisvalues[1] = tmp0;
1806
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1807
basisvalues[0] *= std::sqrt(0.5);
1808
basisvalues[2] *= std::sqrt(1.0);
1809
basisvalues[1] *= std::sqrt(3.0);
1811
// Table(s) of coefficients
1812
static const double coefficients0[3] = \
1813
{0.94280904, -0.28867513, 0.16666667};
1815
static const double coefficients1[3] = \
1816
{-0.47140452, 0.0, -0.33333333};
1819
for (unsigned int r = 0; r < 3; r++)
1821
values[0] += coefficients0[r]*basisvalues[r];
1822
values[1] += coefficients1[r]*basisvalues[r];
1823
}// end loop over 'r'
1825
// Using contravariant Piola transform to map values back to the physical element
1826
const double tmp_ref0 = values[0];
1827
const double tmp_ref1 = values[1];
1828
values[0] = (1.0/detJ)*(J[0]*tmp_ref0 + J[1]*tmp_ref1);
1829
values[1] = (1.0/detJ)*(J[2]*tmp_ref0 + J[3]*tmp_ref1);
1830
values[2] = (1.0/detJ)*(J[4]*tmp_ref0 + J[5]*tmp_ref1);
1836
// Array of basisvalues
1837
double basisvalues[3] = {0.0, 0.0, 0.0};
1839
// Declare helper variables
1840
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1842
// Compute basisvalues
1843
basisvalues[0] = 1.0;
1844
basisvalues[1] = tmp0;
1845
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1846
basisvalues[0] *= std::sqrt(0.5);
1847
basisvalues[2] *= std::sqrt(1.0);
1848
basisvalues[1] *= std::sqrt(3.0);
1850
// Table(s) of coefficients
1851
static const double coefficients0[3] = \
1852
{0.47140452, 0.28867513, -0.16666667};
1854
static const double coefficients1[3] = \
1855
{-0.94280904, 0.0, 0.33333333};
1858
for (unsigned int r = 0; r < 3; r++)
1860
values[0] += coefficients0[r]*basisvalues[r];
1861
values[1] += coefficients1[r]*basisvalues[r];
1862
}// end loop over 'r'
1864
// Using contravariant Piola transform to map values back to the physical element
1865
const double tmp_ref0 = values[0];
1866
const double tmp_ref1 = values[1];
1867
values[0] = (1.0/detJ)*(J[0]*tmp_ref0 + J[1]*tmp_ref1);
1868
values[1] = (1.0/detJ)*(J[2]*tmp_ref0 + J[3]*tmp_ref1);
1869
values[2] = (1.0/detJ)*(J[4]*tmp_ref0 + J[5]*tmp_ref1);
1875
// Array of basisvalues
1876
double basisvalues[1] = {0.0};
1878
// Declare helper variables
1880
// Compute basisvalues
1881
basisvalues[0] = 1.0;
1883
// Table(s) of coefficients
1884
static const double coefficients0[1] = \
1888
for (unsigned int r = 0; r < 1; r++)
1890
values[2] += coefficients0[r]*basisvalues[r];
1891
}// end loop over 'r'
1898
/// Evaluate all basis functions at given point x in cell
1899
virtual void evaluate_basis_all(double* values,
1901
const double* vertex_coordinates,
1902
int cell_orientation) const
1904
// Helper variable to hold values of a single dof.
1905
double dof_values[4] = {0.0, 0.0, 0.0, 0.0};
1907
// Loop dofs and call evaluate_basis
1908
for (unsigned int r = 0; r < 4; r++)
1910
evaluate_basis(r, dof_values, x, vertex_coordinates, cell_orientation);
1911
for (unsigned int s = 0; s < 4; s++)
1913
values[r*4 + s] = dof_values[s];
1914
}// end loop over 's'
1915
}// end loop over 'r'
1918
/// Evaluate order n derivatives of basis function i at given point x in cell
1919
virtual void evaluate_basis_derivatives(std::size_t i,
1923
const double* vertex_coordinates,
1924
int cell_orientation) const
1928
compute_jacobian_triangle_3d(J, vertex_coordinates);
1930
// Compute Jacobian inverse and determinant
1933
compute_jacobian_inverse_triangle_3d(K, detJ, J);
1936
// Check orientation
1937
if (cell_orientation == -1)
1938
throw std::runtime_error("cell orientation must be defined (not -1)");
1939
// (If cell_orientation == 1 = down, multiply det(J) by -1)
1940
else if (cell_orientation == 1)
1944
const double b0 = vertex_coordinates[0];
1945
const double b1 = vertex_coordinates[1];
1946
const double b2 = vertex_coordinates[2];
1948
// P_FFC = J^dag (p - b), P_FIAT = 2*P_FFC - (1, 1)
1949
double X = 2*(K[0]*(x[0] - b0) + K[1]*(x[1] - b1) + K[2]*(x[2] - b2)) - 1.0;
1950
double Y = 2*(K[3]*(x[0] - b0) + K[4]*(x[1] - b1) + K[5]*(x[2] - b2)) - 1.0;
1953
// Compute number of derivatives.
1954
unsigned int num_derivatives_t = 1;
1955
for (unsigned int r = 0; r < n; r++)
1957
num_derivatives_t *= 2;
1958
}// end loop over 'r'
1960
// Compute number of derivatives.
1961
unsigned int num_derivatives_g = 1;
1962
for (unsigned int r = 0; r < n; r++)
1964
num_derivatives_g *= 3;
1965
}// end loop over 'r'
1967
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
1968
unsigned int **combinations_t = new unsigned int *[num_derivatives_t];
1969
for (unsigned int row = 0; row < num_derivatives_t; row++)
1971
combinations_t[row] = new unsigned int [n];
1972
for (unsigned int col = 0; col < n; col++)
1973
combinations_t[row][col] = 0;
1976
// Generate combinations of derivatives
1977
for (unsigned int row = 1; row < num_derivatives_t; row++)
1979
for (unsigned int num = 0; num < row; num++)
1981
for (unsigned int col = n-1; col+1 > 0; col--)
1983
if (combinations_t[row][col] + 1 > 1)
1984
combinations_t[row][col] = 0;
1987
combinations_t[row][col] += 1;
1994
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
1995
unsigned int **combinations_g = new unsigned int *[num_derivatives_g];
1996
for (unsigned int row = 0; row < num_derivatives_g; row++)
1998
combinations_g[row] = new unsigned int [n];
1999
for (unsigned int col = 0; col < n; col++)
2000
combinations_g[row][col] = 0;
2003
// Generate combinations of derivatives
2004
for (unsigned int row = 1; row < num_derivatives_g; row++)
2006
for (unsigned int num = 0; num < row; num++)
2008
for (unsigned int col = n-1; col+1 > 0; col--)
2010
if (combinations_g[row][col] + 1 > 2)
2011
combinations_g[row][col] = 0;
2014
combinations_g[row][col] += 1;
2021
// Compute inverse of Jacobian
2022
const double Jinv[2][3] = {{K[0], K[1], K[2]}, {K[3], K[4], K[5]}};
2024
// Declare transformation matrix
2025
// Declare pointer to two dimensional array and initialise
2026
double **transform = new double *[num_derivatives_g];
2028
for (unsigned int j = 0; j < num_derivatives_g; j++)
2030
transform[j] = new double [num_derivatives_t];
2031
for (unsigned int k = 0; k < num_derivatives_t; k++)
2032
transform[j][k] = 1;
2035
// Construct transformation matrix
2036
for (unsigned int row = 0; row < num_derivatives_g; row++)
2038
for (unsigned int col = 0; col < num_derivatives_t; col++)
2040
for (unsigned int k = 0; k < n; k++)
2041
transform[row][col] *= Jinv[combinations_t[col][k]][combinations_g[row][k]];
2045
// Reset values. Assuming that values is always an array.
2046
for (unsigned int r = 0; r < 4*num_derivatives_g; r++)
2049
}// end loop over 'r'
2056
// Array of basisvalues
2057
double basisvalues[3] = {0.0, 0.0, 0.0};
2059
// Declare helper variables
2060
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2062
// Compute basisvalues
2063
basisvalues[0] = 1.0;
2064
basisvalues[1] = tmp0;
2065
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
2066
basisvalues[0] *= std::sqrt(0.5);
2067
basisvalues[2] *= std::sqrt(1.0);
2068
basisvalues[1] *= std::sqrt(3.0);
2070
// Table(s) of coefficients
2071
static const double coefficients0[3] = \
2072
{0.47140452, 0.28867513, -0.16666667};
2074
static const double coefficients1[3] = \
2075
{0.47140452, 0.0, 0.33333333};
2077
// Tables of derivatives of the polynomial base (transpose).
2078
static const double dmats0[3][3] = \
2080
{4.8989795, 0.0, 0.0},
2083
static const double dmats1[3][3] = \
2085
{2.4494897, 0.0, 0.0},
2086
{4.2426407, 0.0, 0.0}};
2088
// Compute reference derivatives.
2089
// Declare pointer to array of derivatives on FIAT element.
2090
double *derivatives = new double[2*num_derivatives_t];
2091
for (unsigned int r = 0; r < 2*num_derivatives_t; r++)
2093
derivatives[r] = 0.0;
2094
}// end loop over 'r'
2096
// Declare pointer to array of reference derivatives on physical element.
2097
double *derivatives_p = new double[3*num_derivatives_t];
2098
for (unsigned int r = 0; r < 3*num_derivatives_t; r++)
2100
derivatives_p[r] = 0.0;
2101
}// end loop over 'r'
2103
// Declare derivative matrix (of polynomial basis).
2104
double dmats[3][3] = \
2109
// Declare (auxiliary) derivative matrix (of polynomial basis).
2110
double dmats_old[3][3] = \
2115
// Loop possible derivatives.
2116
for (unsigned int r = 0; r < num_derivatives_t; r++)
2118
// Resetting dmats values to compute next derivative.
2119
for (unsigned int t = 0; t < 3; t++)
2121
for (unsigned int u = 0; u < 3; u++)
2129
}// end loop over 'u'
2130
}// end loop over 't'
2132
// Looping derivative order to generate dmats.
2133
for (unsigned int s = 0; s < n; s++)
2135
// Updating dmats_old with new values and resetting dmats.
2136
for (unsigned int t = 0; t < 3; t++)
2138
for (unsigned int u = 0; u < 3; u++)
2140
dmats_old[t][u] = dmats[t][u];
2142
}// end loop over 'u'
2143
}// end loop over 't'
2145
// Update dmats using an inner product.
2146
if (combinations_t[r][s] == 0)
2148
for (unsigned int t = 0; t < 3; t++)
2150
for (unsigned int u = 0; u < 3; u++)
2152
for (unsigned int tu = 0; tu < 3; tu++)
2154
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2155
}// end loop over 'tu'
2156
}// end loop over 'u'
2157
}// end loop over 't'
2160
if (combinations_t[r][s] == 1)
2162
for (unsigned int t = 0; t < 3; t++)
2164
for (unsigned int u = 0; u < 3; u++)
2166
for (unsigned int tu = 0; tu < 3; tu++)
2168
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2169
}// end loop over 'tu'
2170
}// end loop over 'u'
2171
}// end loop over 't'
2174
}// end loop over 's'
2175
for (unsigned int s = 0; s < 3; s++)
2177
for (unsigned int t = 0; t < 3; t++)
2179
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2180
derivatives[num_derivatives_t + r] += coefficients1[s]*dmats[s][t]*basisvalues[t];
2181
}// end loop over 't'
2182
}// end loop over 's'
2184
// Using contravariant Piola transform to map values back to the physical element.
2185
const double tmp_ref0 = derivatives[r];
2186
const double tmp_ref1 = derivatives[num_derivatives_t + r];
2187
derivatives_p[r] = (1.0/detJ)*(J[0]*tmp_ref0 + J[1]*tmp_ref1);
2188
derivatives_p[num_derivatives_t + r] = (1.0/detJ)*(J[2]*tmp_ref0 + J[3]*tmp_ref1);
2189
derivatives_p[2*num_derivatives_t + r] = (1.0/detJ)*(J[4]*tmp_ref0 + J[5]*tmp_ref1);
2190
}// end loop over 'r'
2192
// Transform derivatives back to physical element
2193
for (unsigned int r = 0; r < num_derivatives_g; r++)
2195
for (unsigned int s = 0; s < num_derivatives_t; s++)
2197
values[r] += transform[r][s]*derivatives_p[s];
2198
values[num_derivatives_g + r] += transform[r][s]*derivatives_p[num_derivatives_t + s];
2199
values[2*num_derivatives_g + r] += transform[r][s]*derivatives_p[2*num_derivatives_t + s];
2200
}// end loop over 's'
2201
}// end loop over 'r'
2203
// Delete pointer to array of derivatives on FIAT element
2204
delete [] derivatives;
2206
// Delete pointer to array of reference derivatives on physical element.
2207
delete [] derivatives_p;
2209
// Delete pointer to array of combinations of derivatives and transform
2210
for (unsigned int r = 0; r < num_derivatives_t; r++)
2212
delete [] combinations_t[r];
2213
}// end loop over 'r'
2214
delete [] combinations_t;
2215
for (unsigned int r = 0; r < num_derivatives_g; r++)
2217
delete [] combinations_g[r];
2218
}// end loop over 'r'
2219
delete [] combinations_g;
2220
for (unsigned int r = 0; r < num_derivatives_g; r++)
2222
delete [] transform[r];
2223
}// end loop over 'r'
2224
delete [] transform;
2230
// Array of basisvalues
2231
double basisvalues[3] = {0.0, 0.0, 0.0};
2233
// Declare helper variables
2234
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2236
// Compute basisvalues
2237
basisvalues[0] = 1.0;
2238
basisvalues[1] = tmp0;
2239
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
2240
basisvalues[0] *= std::sqrt(0.5);
2241
basisvalues[2] *= std::sqrt(1.0);
2242
basisvalues[1] *= std::sqrt(3.0);
2244
// Table(s) of coefficients
2245
static const double coefficients0[3] = \
2246
{0.94280904, -0.28867513, 0.16666667};
2248
static const double coefficients1[3] = \
2249
{-0.47140452, 0.0, -0.33333333};
2251
// Tables of derivatives of the polynomial base (transpose).
2252
static const double dmats0[3][3] = \
2254
{4.8989795, 0.0, 0.0},
2257
static const double dmats1[3][3] = \
2259
{2.4494897, 0.0, 0.0},
2260
{4.2426407, 0.0, 0.0}};
2262
// Compute reference derivatives.
2263
// Declare pointer to array of derivatives on FIAT element.
2264
double *derivatives = new double[2*num_derivatives_t];
2265
for (unsigned int r = 0; r < 2*num_derivatives_t; r++)
2267
derivatives[r] = 0.0;
2268
}// end loop over 'r'
2270
// Declare pointer to array of reference derivatives on physical element.
2271
double *derivatives_p = new double[3*num_derivatives_t];
2272
for (unsigned int r = 0; r < 3*num_derivatives_t; r++)
2274
derivatives_p[r] = 0.0;
2275
}// end loop over 'r'
2277
// Declare derivative matrix (of polynomial basis).
2278
double dmats[3][3] = \
2283
// Declare (auxiliary) derivative matrix (of polynomial basis).
2284
double dmats_old[3][3] = \
2289
// Loop possible derivatives.
2290
for (unsigned int r = 0; r < num_derivatives_t; r++)
2292
// Resetting dmats values to compute next derivative.
2293
for (unsigned int t = 0; t < 3; t++)
2295
for (unsigned int u = 0; u < 3; u++)
2303
}// end loop over 'u'
2304
}// end loop over 't'
2306
// Looping derivative order to generate dmats.
2307
for (unsigned int s = 0; s < n; s++)
2309
// Updating dmats_old with new values and resetting dmats.
2310
for (unsigned int t = 0; t < 3; t++)
2312
for (unsigned int u = 0; u < 3; u++)
2314
dmats_old[t][u] = dmats[t][u];
2316
}// end loop over 'u'
2317
}// end loop over 't'
2319
// Update dmats using an inner product.
2320
if (combinations_t[r][s] == 0)
2322
for (unsigned int t = 0; t < 3; t++)
2324
for (unsigned int u = 0; u < 3; u++)
2326
for (unsigned int tu = 0; tu < 3; tu++)
2328
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2329
}// end loop over 'tu'
2330
}// end loop over 'u'
2331
}// end loop over 't'
2334
if (combinations_t[r][s] == 1)
2336
for (unsigned int t = 0; t < 3; t++)
2338
for (unsigned int u = 0; u < 3; u++)
2340
for (unsigned int tu = 0; tu < 3; tu++)
2342
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2343
}// end loop over 'tu'
2344
}// end loop over 'u'
2345
}// end loop over 't'
2348
}// end loop over 's'
2349
for (unsigned int s = 0; s < 3; s++)
2351
for (unsigned int t = 0; t < 3; t++)
2353
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2354
derivatives[num_derivatives_t + r] += coefficients1[s]*dmats[s][t]*basisvalues[t];
2355
}// end loop over 't'
2356
}// end loop over 's'
2358
// Using contravariant Piola transform to map values back to the physical element.
2359
const double tmp_ref0 = derivatives[r];
2360
const double tmp_ref1 = derivatives[num_derivatives_t + r];
2361
derivatives_p[r] = (1.0/detJ)*(J[0]*tmp_ref0 + J[1]*tmp_ref1);
2362
derivatives_p[num_derivatives_t + r] = (1.0/detJ)*(J[2]*tmp_ref0 + J[3]*tmp_ref1);
2363
derivatives_p[2*num_derivatives_t + r] = (1.0/detJ)*(J[4]*tmp_ref0 + J[5]*tmp_ref1);
2364
}// end loop over 'r'
2366
// Transform derivatives back to physical element
2367
for (unsigned int r = 0; r < num_derivatives_g; r++)
2369
for (unsigned int s = 0; s < num_derivatives_t; s++)
2371
values[r] += transform[r][s]*derivatives_p[s];
2372
values[num_derivatives_g + r] += transform[r][s]*derivatives_p[num_derivatives_t + s];
2373
values[2*num_derivatives_g + r] += transform[r][s]*derivatives_p[2*num_derivatives_t + s];
2374
}// end loop over 's'
2375
}// end loop over 'r'
2377
// Delete pointer to array of derivatives on FIAT element
2378
delete [] derivatives;
2380
// Delete pointer to array of reference derivatives on physical element.
2381
delete [] derivatives_p;
2383
// Delete pointer to array of combinations of derivatives and transform
2384
for (unsigned int r = 0; r < num_derivatives_t; r++)
2386
delete [] combinations_t[r];
2387
}// end loop over 'r'
2388
delete [] combinations_t;
2389
for (unsigned int r = 0; r < num_derivatives_g; r++)
2391
delete [] combinations_g[r];
2392
}// end loop over 'r'
2393
delete [] combinations_g;
2394
for (unsigned int r = 0; r < num_derivatives_g; r++)
2396
delete [] transform[r];
2397
}// end loop over 'r'
2398
delete [] transform;
2404
// Array of basisvalues
2405
double basisvalues[3] = {0.0, 0.0, 0.0};
2407
// Declare helper variables
2408
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2410
// Compute basisvalues
2411
basisvalues[0] = 1.0;
2412
basisvalues[1] = tmp0;
2413
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
2414
basisvalues[0] *= std::sqrt(0.5);
2415
basisvalues[2] *= std::sqrt(1.0);
2416
basisvalues[1] *= std::sqrt(3.0);
2418
// Table(s) of coefficients
2419
static const double coefficients0[3] = \
2420
{0.47140452, 0.28867513, -0.16666667};
2422
static const double coefficients1[3] = \
2423
{-0.94280904, 0.0, 0.33333333};
2425
// Tables of derivatives of the polynomial base (transpose).
2426
static const double dmats0[3][3] = \
2428
{4.8989795, 0.0, 0.0},
2431
static const double dmats1[3][3] = \
2433
{2.4494897, 0.0, 0.0},
2434
{4.2426407, 0.0, 0.0}};
2436
// Compute reference derivatives.
2437
// Declare pointer to array of derivatives on FIAT element.
2438
double *derivatives = new double[2*num_derivatives_t];
2439
for (unsigned int r = 0; r < 2*num_derivatives_t; r++)
2441
derivatives[r] = 0.0;
2442
}// end loop over 'r'
2444
// Declare pointer to array of reference derivatives on physical element.
2445
double *derivatives_p = new double[3*num_derivatives_t];
2446
for (unsigned int r = 0; r < 3*num_derivatives_t; r++)
2448
derivatives_p[r] = 0.0;
2449
}// end loop over 'r'
2451
// Declare derivative matrix (of polynomial basis).
2452
double dmats[3][3] = \
2457
// Declare (auxiliary) derivative matrix (of polynomial basis).
2458
double dmats_old[3][3] = \
2463
// Loop possible derivatives.
2464
for (unsigned int r = 0; r < num_derivatives_t; r++)
2466
// Resetting dmats values to compute next derivative.
2467
for (unsigned int t = 0; t < 3; t++)
2469
for (unsigned int u = 0; u < 3; u++)
2477
}// end loop over 'u'
2478
}// end loop over 't'
2480
// Looping derivative order to generate dmats.
2481
for (unsigned int s = 0; s < n; s++)
2483
// Updating dmats_old with new values and resetting dmats.
2484
for (unsigned int t = 0; t < 3; t++)
2486
for (unsigned int u = 0; u < 3; u++)
2488
dmats_old[t][u] = dmats[t][u];
2490
}// end loop over 'u'
2491
}// end loop over 't'
2493
// Update dmats using an inner product.
2494
if (combinations_t[r][s] == 0)
2496
for (unsigned int t = 0; t < 3; t++)
2498
for (unsigned int u = 0; u < 3; u++)
2500
for (unsigned int tu = 0; tu < 3; tu++)
2502
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2503
}// end loop over 'tu'
2504
}// end loop over 'u'
2505
}// end loop over 't'
2508
if (combinations_t[r][s] == 1)
2510
for (unsigned int t = 0; t < 3; t++)
2512
for (unsigned int u = 0; u < 3; u++)
2514
for (unsigned int tu = 0; tu < 3; tu++)
2516
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2517
}// end loop over 'tu'
2518
}// end loop over 'u'
2519
}// end loop over 't'
2522
}// end loop over 's'
2523
for (unsigned int s = 0; s < 3; s++)
2525
for (unsigned int t = 0; t < 3; t++)
2527
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2528
derivatives[num_derivatives_t + r] += coefficients1[s]*dmats[s][t]*basisvalues[t];
2529
}// end loop over 't'
2530
}// end loop over 's'
2532
// Using contravariant Piola transform to map values back to the physical element.
2533
const double tmp_ref0 = derivatives[r];
2534
const double tmp_ref1 = derivatives[num_derivatives_t + r];
2535
derivatives_p[r] = (1.0/detJ)*(J[0]*tmp_ref0 + J[1]*tmp_ref1);
2536
derivatives_p[num_derivatives_t + r] = (1.0/detJ)*(J[2]*tmp_ref0 + J[3]*tmp_ref1);
2537
derivatives_p[2*num_derivatives_t + r] = (1.0/detJ)*(J[4]*tmp_ref0 + J[5]*tmp_ref1);
2538
}// end loop over 'r'
2540
// Transform derivatives back to physical element
2541
for (unsigned int r = 0; r < num_derivatives_g; r++)
2543
for (unsigned int s = 0; s < num_derivatives_t; s++)
2545
values[r] += transform[r][s]*derivatives_p[s];
2546
values[num_derivatives_g + r] += transform[r][s]*derivatives_p[num_derivatives_t + s];
2547
values[2*num_derivatives_g + r] += transform[r][s]*derivatives_p[2*num_derivatives_t + s];
2548
}// end loop over 's'
2549
}// end loop over 'r'
2551
// Delete pointer to array of derivatives on FIAT element
2552
delete [] derivatives;
2554
// Delete pointer to array of reference derivatives on physical element.
2555
delete [] derivatives_p;
2557
// Delete pointer to array of combinations of derivatives and transform
2558
for (unsigned int r = 0; r < num_derivatives_t; r++)
2560
delete [] combinations_t[r];
2561
}// end loop over 'r'
2562
delete [] combinations_t;
2563
for (unsigned int r = 0; r < num_derivatives_g; r++)
2565
delete [] combinations_g[r];
2566
}// end loop over 'r'
2567
delete [] combinations_g;
2568
for (unsigned int r = 0; r < num_derivatives_g; r++)
2570
delete [] transform[r];
2571
}// end loop over 'r'
2572
delete [] transform;
2578
// Array of basisvalues
2579
double basisvalues[1] = {0.0};
2581
// Declare helper variables
2583
// Compute basisvalues
2584
basisvalues[0] = 1.0;
2586
// Table(s) of coefficients
2587
static const double coefficients0[1] = \
2590
// Tables of derivatives of the polynomial base (transpose).
2591
static const double dmats0[1][1] = \
2594
static const double dmats1[1][1] = \
2597
// Compute reference derivatives.
2598
// Declare pointer to array of derivatives on FIAT element.
2599
double *derivatives = new double[num_derivatives_t];
2600
for (unsigned int r = 0; r < num_derivatives_t; r++)
2602
derivatives[r] = 0.0;
2603
}// end loop over 'r'
2605
// Declare derivative matrix (of polynomial basis).
2606
double dmats[1][1] = \
2609
// Declare (auxiliary) derivative matrix (of polynomial basis).
2610
double dmats_old[1][1] = \
2613
// Loop possible derivatives.
2614
for (unsigned int r = 0; r < num_derivatives_t; r++)
2616
// Resetting dmats values to compute next derivative.
2617
for (unsigned int t = 0; t < 1; t++)
2619
for (unsigned int u = 0; u < 1; u++)
2627
}// end loop over 'u'
2628
}// end loop over 't'
2630
// Looping derivative order to generate dmats.
2631
for (unsigned int s = 0; s < n; s++)
2633
// Updating dmats_old with new values and resetting dmats.
2634
for (unsigned int t = 0; t < 1; t++)
2636
for (unsigned int u = 0; u < 1; u++)
2638
dmats_old[t][u] = dmats[t][u];
2640
}// end loop over 'u'
2641
}// end loop over 't'
2643
// Update dmats using an inner product.
2644
if (combinations_t[r][s] == 0)
2646
for (unsigned int t = 0; t < 1; t++)
2648
for (unsigned int u = 0; u < 1; u++)
2650
for (unsigned int tu = 0; tu < 1; tu++)
2652
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2653
}// end loop over 'tu'
2654
}// end loop over 'u'
2655
}// end loop over 't'
2658
if (combinations_t[r][s] == 1)
2660
for (unsigned int t = 0; t < 1; t++)
2662
for (unsigned int u = 0; u < 1; u++)
2664
for (unsigned int tu = 0; tu < 1; tu++)
2666
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2667
}// end loop over 'tu'
2668
}// end loop over 'u'
2669
}// end loop over 't'
2672
}// end loop over 's'
2673
for (unsigned int s = 0; s < 1; s++)
2675
for (unsigned int t = 0; t < 1; t++)
2677
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2678
}// end loop over 't'
2679
}// end loop over 's'
2680
}// end loop over 'r'
2682
// Transform derivatives back to physical element
2683
for (unsigned int r = 0; r < num_derivatives_g; r++)
2685
for (unsigned int s = 0; s < num_derivatives_t; s++)
2687
values[2*num_derivatives_g + r] += transform[r][s]*derivatives[s];
2688
}// end loop over 's'
2689
}// end loop over 'r'
2691
// Delete pointer to array of derivatives on FIAT element
2692
delete [] derivatives;
2694
// Delete pointer to array of combinations of derivatives and transform
2695
for (unsigned int r = 0; r < num_derivatives_t; r++)
2697
delete [] combinations_t[r];
2698
}// end loop over 'r'
2699
delete [] combinations_t;
2700
for (unsigned int r = 0; r < num_derivatives_g; r++)
2702
delete [] combinations_g[r];
2703
}// end loop over 'r'
2704
delete [] combinations_g;
2705
for (unsigned int r = 0; r < num_derivatives_g; r++)
2707
delete [] transform[r];
2708
}// end loop over 'r'
2709
delete [] transform;
2716
/// Evaluate order n derivatives of all basis functions at given point x in cell
2717
virtual void evaluate_basis_derivatives_all(std::size_t n,
2720
const double* vertex_coordinates,
2721
int cell_orientation) const
2723
// Compute number of derivatives.
2724
unsigned int num_derivatives_g = 1;
2725
for (unsigned int r = 0; r < n; r++)
2727
num_derivatives_g *= 3;
2728
}// end loop over 'r'
2730
// Helper variable to hold values of a single dof.
2731
double *dof_values = new double[4*num_derivatives_g];
2732
for (unsigned int r = 0; r < 4*num_derivatives_g; r++)
2734
dof_values[r] = 0.0;
2735
}// end loop over 'r'
2737
// Loop dofs and call evaluate_basis_derivatives.
2738
for (unsigned int r = 0; r < 4; r++)
2740
evaluate_basis_derivatives(r, n, dof_values, x, vertex_coordinates, cell_orientation);
2741
for (unsigned int s = 0; s < 4*num_derivatives_g; s++)
2743
values[r*4*num_derivatives_g + s] = dof_values[s];
2744
}// end loop over 's'
2745
}// end loop over 'r'
2748
delete [] dof_values;
2751
/// Evaluate linear functional for dof i on the function f
2752
virtual double evaluate_dof(std::size_t i,
2753
const ufc::function& f,
2754
const double* vertex_coordinates,
2755
int cell_orientation,
2756
const ufc::cell& c) const
2758
// Declare variables for result of evaluation
2761
// Declare variable for physical coordinates
2767
compute_jacobian_triangle_3d(J, vertex_coordinates);
2770
// Compute Jacobian inverse and determinant
2773
compute_jacobian_inverse_triangle_3d(K, detJ, J);
2777
// Check orientation
2778
if (cell_orientation == -1)
2779
throw std::runtime_error("cell orientation must be defined (not -1)");
2780
// (If cell_orientation == 1 = down, multiply det(J) by -1)
2781
else if (cell_orientation == 1)
2787
y[0] = 0.5*vertex_coordinates[3] + 0.5*vertex_coordinates[6];
2788
y[1] = 0.5*vertex_coordinates[4] + 0.5*vertex_coordinates[7];
2789
y[2] = 0.5*vertex_coordinates[5] + 0.5*vertex_coordinates[8];
2790
f.evaluate(vals, y, c);
2791
result = (detJ*(K[0]*vals[0] + K[1]*vals[1] + K[2]*vals[2])) + (detJ*(K[3]*vals[0] + K[4]*vals[1] + K[5]*vals[2]));
2797
y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[6];
2798
y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[7];
2799
y[2] = 0.5*vertex_coordinates[2] + 0.5*vertex_coordinates[8];
2800
f.evaluate(vals, y, c);
2801
result = (detJ*(K[0]*vals[0] + K[1]*vals[1] + K[2]*vals[2]));
2807
y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[3];
2808
y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[4];
2809
y[2] = 0.5*vertex_coordinates[2] + 0.5*vertex_coordinates[5];
2810
f.evaluate(vals, y, c);
2811
result = (-1.0)*(detJ*(K[3]*vals[0] + K[4]*vals[1] + K[5]*vals[2]));
2817
y[0] = 0.33333333*vertex_coordinates[0] + 0.33333333*vertex_coordinates[3] + 0.33333333*vertex_coordinates[6];
2818
y[1] = 0.33333333*vertex_coordinates[1] + 0.33333333*vertex_coordinates[4] + 0.33333333*vertex_coordinates[7];
2819
y[2] = 0.33333333*vertex_coordinates[2] + 0.33333333*vertex_coordinates[5] + 0.33333333*vertex_coordinates[8];
2820
f.evaluate(vals, y, c);
2829
/// Evaluate linear functionals for all dofs on the function f
2830
virtual void evaluate_dofs(double* values,
2831
const ufc::function& f,
2832
const double* vertex_coordinates,
2833
int cell_orientation,
2834
const ufc::cell& c) const
2836
// Declare variables for result of evaluation
2839
// Declare variable for physical coordinates
2845
compute_jacobian_triangle_3d(J, vertex_coordinates);
2848
// Compute Jacobian inverse and determinant
2851
compute_jacobian_inverse_triangle_3d(K, detJ, J);
2855
// Check orientation
2856
if (cell_orientation == -1)
2857
throw std::runtime_error("cell orientation must be defined (not -1)");
2858
// (If cell_orientation == 1 = down, multiply det(J) by -1)
2859
else if (cell_orientation == 1)
2861
y[0] = 0.5*vertex_coordinates[3] + 0.5*vertex_coordinates[6];
2862
y[1] = 0.5*vertex_coordinates[4] + 0.5*vertex_coordinates[7];
2863
y[2] = 0.5*vertex_coordinates[5] + 0.5*vertex_coordinates[8];
2864
f.evaluate(vals, y, c);
2865
result = (detJ*(K[0]*vals[0] + K[1]*vals[1] + K[2]*vals[2])) + (detJ*(K[3]*vals[0] + K[4]*vals[1] + K[5]*vals[2]));
2867
y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[6];
2868
y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[7];
2869
y[2] = 0.5*vertex_coordinates[2] + 0.5*vertex_coordinates[8];
2870
f.evaluate(vals, y, c);
2871
result = (detJ*(K[0]*vals[0] + K[1]*vals[1] + K[2]*vals[2]));
2873
y[0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[3];
2874
y[1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[4];
2875
y[2] = 0.5*vertex_coordinates[2] + 0.5*vertex_coordinates[5];
2876
f.evaluate(vals, y, c);
2877
result = (-1.0)*(detJ*(K[3]*vals[0] + K[4]*vals[1] + K[5]*vals[2]));
2879
y[0] = 0.33333333*vertex_coordinates[0] + 0.33333333*vertex_coordinates[3] + 0.33333333*vertex_coordinates[6];
2880
y[1] = 0.33333333*vertex_coordinates[1] + 0.33333333*vertex_coordinates[4] + 0.33333333*vertex_coordinates[7];
2881
y[2] = 0.33333333*vertex_coordinates[2] + 0.33333333*vertex_coordinates[5] + 0.33333333*vertex_coordinates[8];
2882
f.evaluate(vals, y, c);
2883
values[3] = vals[3];
2886
/// Interpolate vertex values from dof values
2887
virtual void interpolate_vertex_values(double* vertex_values,
2888
const double* dof_values,
2889
const double* vertex_coordinates,
2890
int cell_orientation,
2891
const ufc::cell& c) const
2895
compute_jacobian_triangle_3d(J, vertex_coordinates);
2898
// Compute Jacobian inverse and determinant
2901
compute_jacobian_inverse_triangle_3d(K, detJ, J);
2905
// Check orientation
2906
if (cell_orientation == -1)
2907
throw std::runtime_error("cell orientation must be defined (not -1)");
2908
// (If cell_orientation == 1 = down, multiply det(J) by -1)
2909
else if (cell_orientation == 1)
2912
// Evaluate function and change variables
2913
vertex_values[0] = dof_values[1]*(1.0/detJ)*J[0] + dof_values[2]*((1.0/detJ)*(J[1]*(-1.0)));
2914
vertex_values[3] = dof_values[0]*(1.0/detJ)*J[0] + dof_values[2]*((1.0/detJ)*(J[0] + J[1]*(-1.0)));
2915
vertex_values[6] = dof_values[0]*(1.0/detJ)*J[1] + dof_values[1]*((1.0/detJ)*(J[0] + J[1]*(-1.0)));
2916
vertex_values[1] = dof_values[1]*(1.0/detJ)*J[2] + dof_values[2]*((1.0/detJ)*(J[3]*(-1.0)));
2917
vertex_values[4] = dof_values[0]*(1.0/detJ)*J[2] + dof_values[2]*((1.0/detJ)*(J[2] + J[3]*(-1.0)));
2918
vertex_values[7] = dof_values[0]*(1.0/detJ)*J[3] + dof_values[1]*((1.0/detJ)*(J[2] + J[3]*(-1.0)));
2919
// Evaluate function and change variables
2920
vertex_values[2] = dof_values[3];
2921
vertex_values[5] = dof_values[3];
2922
vertex_values[8] = dof_values[3];
2925
/// Map coordinate xhat from reference cell to coordinate x in cell
2926
virtual void map_from_reference_cell(double* x,
2928
const ufc::cell& c) const
2930
std::cerr << "*** FFC warning: " << "map_from_reference_cell not yet implemented." << std::endl;
2933
/// Map from coordinate x in cell to coordinate xhat in reference cell
2934
virtual void map_to_reference_cell(double* xhat,
2936
const ufc::cell& c) const
2938
std::cerr << "*** FFC warning: " << "map_to_reference_cell not yet implemented." << std::endl;
2941
/// Return the number of sub elements (for a mixed element)
2942
virtual std::size_t num_sub_elements() const
2947
/// Create a new finite element for sub element i (for a mixed element)
2948
virtual ufc::finite_element* create_sub_element(std::size_t i) const
2954
return new projectionmanifold_finite_element_0();
2959
return new projectionmanifold_finite_element_1();
2967
/// Create a new class instance
2968
virtual ufc::finite_element* create() const
2970
return new projectionmanifold_finite_element_2();
2975
/// This class defines the interface for a local-to-global mapping of
2976
/// degrees of freedom (dofs).
2978
class projectionmanifold_dofmap_0: public ufc::dofmap
2983
projectionmanifold_dofmap_0() : ufc::dofmap()
2989
virtual ~projectionmanifold_dofmap_0()
2994
/// Return a string identifying the dofmap
2995
virtual const char* signature() const
2997
return "FFC dofmap for FiniteElement('Raviart-Thomas', Domain(Cell('triangle', 3), 'triangle_multiverse', 3, 2), 1, None)";
3000
/// Return true iff mesh entities of topological dimension d are needed
3001
virtual bool needs_mesh_entities(std::size_t d) const
3025
/// Return the topological dimension of the associated cell shape
3026
virtual std::size_t topological_dimension() const
3031
/// Return the geometric dimension of the associated cell shape
3032
virtual std::size_t geometric_dimension() const
3037
/// Return the dimension of the global finite element function space
3038
virtual std::size_t global_dimension(const std::vector<std::size_t>&
3039
num_global_entities) const
3041
return num_global_entities[1];
3044
/// Return the dimension of the local finite element function space for a cell
3045
virtual std::size_t local_dimension(const ufc::cell& c) const
3050
/// Return the maximum dimension of the local finite element function space
3051
virtual std::size_t max_local_dimension() const
3056
/// Return the number of dofs on each cell facet
3057
virtual std::size_t num_facet_dofs() const
3062
/// Return the number of dofs associated with each cell entity of dimension d
3063
virtual std::size_t num_entity_dofs(std::size_t d) const
3087
/// Tabulate the local-to-global mapping of dofs on a cell
3088
virtual void tabulate_dofs(std::size_t* dofs,
3089
const std::vector<std::size_t>& num_global_entities,
3090
const ufc::cell& c) const
3092
dofs[0] = c.entity_indices[1][0];
3093
dofs[1] = c.entity_indices[1][1];
3094
dofs[2] = c.entity_indices[1][2];
3097
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
3098
virtual void tabulate_facet_dofs(std::size_t* dofs,
3099
std::size_t facet) const
3122
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
3123
virtual void tabulate_entity_dofs(std::size_t* dofs,
3124
std::size_t d, std::size_t i) const
3128
std::cerr << "*** FFC warning: " << "d is larger than dimension (2)" << std::endl;
3142
std::cerr << "*** FFC warning: " << "i is larger than number of entities (2)" << std::endl;
3175
/// Tabulate the coordinates of all dofs on a cell
3176
virtual void tabulate_coordinates(double** dof_coordinates,
3177
const double* vertex_coordinates) const
3179
dof_coordinates[0][0] = 0.5*vertex_coordinates[3] + 0.5*vertex_coordinates[6];
3180
dof_coordinates[0][1] = 0.5*vertex_coordinates[4] + 0.5*vertex_coordinates[7];
3181
dof_coordinates[0][2] = 0.5*vertex_coordinates[5] + 0.5*vertex_coordinates[8];
3182
dof_coordinates[1][0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[6];
3183
dof_coordinates[1][1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[7];
3184
dof_coordinates[1][2] = 0.5*vertex_coordinates[2] + 0.5*vertex_coordinates[8];
3185
dof_coordinates[2][0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[3];
3186
dof_coordinates[2][1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[4];
3187
dof_coordinates[2][2] = 0.5*vertex_coordinates[2] + 0.5*vertex_coordinates[5];
3190
/// Return the number of sub dofmaps (for a mixed element)
3191
virtual std::size_t num_sub_dofmaps() const
3196
/// Create a new dofmap for sub dofmap i (for a mixed element)
3197
virtual ufc::dofmap* create_sub_dofmap(std::size_t i) const
3202
/// Create a new class instance
3203
virtual ufc::dofmap* create() const
3205
return new projectionmanifold_dofmap_0();
3210
/// This class defines the interface for a local-to-global mapping of
3211
/// degrees of freedom (dofs).
3213
class projectionmanifold_dofmap_1: public ufc::dofmap
3218
projectionmanifold_dofmap_1() : ufc::dofmap()
3224
virtual ~projectionmanifold_dofmap_1()
3229
/// Return a string identifying the dofmap
3230
virtual const char* signature() const
3232
return "FFC dofmap for FiniteElement('Discontinuous Lagrange', Domain(Cell('triangle', 3), 'triangle_multiverse', 3, 2), 0, None)";
3235
/// Return true iff mesh entities of topological dimension d are needed
3236
virtual bool needs_mesh_entities(std::size_t d) const
3260
/// Return the topological dimension of the associated cell shape
3261
virtual std::size_t topological_dimension() const
3266
/// Return the geometric dimension of the associated cell shape
3267
virtual std::size_t geometric_dimension() const
3272
/// Return the dimension of the global finite element function space
3273
virtual std::size_t global_dimension(const std::vector<std::size_t>&
3274
num_global_entities) const
3276
return num_global_entities[2];
3279
/// Return the dimension of the local finite element function space for a cell
3280
virtual std::size_t local_dimension(const ufc::cell& c) const
3285
/// Return the maximum dimension of the local finite element function space
3286
virtual std::size_t max_local_dimension() const
3291
/// Return the number of dofs on each cell facet
3292
virtual std::size_t num_facet_dofs() const
3297
/// Return the number of dofs associated with each cell entity of dimension d
3298
virtual std::size_t num_entity_dofs(std::size_t d) const
3322
/// Tabulate the local-to-global mapping of dofs on a cell
3323
virtual void tabulate_dofs(std::size_t* dofs,
3324
const std::vector<std::size_t>& num_global_entities,
3325
const ufc::cell& c) const
3327
dofs[0] = c.entity_indices[2][0];
3330
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
3331
virtual void tabulate_facet_dofs(std::size_t* dofs,
3332
std::size_t facet) const
3355
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
3356
virtual void tabulate_entity_dofs(std::size_t* dofs,
3357
std::size_t d, std::size_t i) const
3361
std::cerr << "*** FFC warning: " << "d is larger than dimension (2)" << std::endl;
3380
std::cerr << "*** FFC warning: " << "i is larger than number of entities (0)" << std::endl;
3390
/// Tabulate the coordinates of all dofs on a cell
3391
virtual void tabulate_coordinates(double** dof_coordinates,
3392
const double* vertex_coordinates) const
3394
dof_coordinates[0][0] = 0.33333333*vertex_coordinates[0] + 0.33333333*vertex_coordinates[3] + 0.33333333*vertex_coordinates[6];
3395
dof_coordinates[0][1] = 0.33333333*vertex_coordinates[1] + 0.33333333*vertex_coordinates[4] + 0.33333333*vertex_coordinates[7];
3396
dof_coordinates[0][2] = 0.33333333*vertex_coordinates[2] + 0.33333333*vertex_coordinates[5] + 0.33333333*vertex_coordinates[8];
3399
/// Return the number of sub dofmaps (for a mixed element)
3400
virtual std::size_t num_sub_dofmaps() const
3405
/// Create a new dofmap for sub dofmap i (for a mixed element)
3406
virtual ufc::dofmap* create_sub_dofmap(std::size_t i) const
3411
/// Create a new class instance
3412
virtual ufc::dofmap* create() const
3414
return new projectionmanifold_dofmap_1();
3419
/// This class defines the interface for a local-to-global mapping of
3420
/// degrees of freedom (dofs).
3422
class projectionmanifold_dofmap_2: public ufc::dofmap
3427
projectionmanifold_dofmap_2() : ufc::dofmap()
3433
virtual ~projectionmanifold_dofmap_2()
3438
/// Return a string identifying the dofmap
3439
virtual const char* signature() const
3441
return "FFC dofmap for MixedElement(*[FiniteElement('Raviart-Thomas', Domain(Cell('triangle', 3), 'triangle_multiverse', 3, 2), 1, None), FiniteElement('Discontinuous Lagrange', Domain(Cell('triangle', 3), 'triangle_multiverse', 3, 2), 0, None)], **{'value_shape': (4,) })";
3444
/// Return true iff mesh entities of topological dimension d are needed
3445
virtual bool needs_mesh_entities(std::size_t d) const
3469
/// Return the topological dimension of the associated cell shape
3470
virtual std::size_t topological_dimension() const
3475
/// Return the geometric dimension of the associated cell shape
3476
virtual std::size_t geometric_dimension() const
3481
/// Return the dimension of the global finite element function space
3482
virtual std::size_t global_dimension(const std::vector<std::size_t>&
3483
num_global_entities) const
3485
return num_global_entities[1] + num_global_entities[2];
3488
/// Return the dimension of the local finite element function space for a cell
3489
virtual std::size_t local_dimension(const ufc::cell& c) const
3494
/// Return the maximum dimension of the local finite element function space
3495
virtual std::size_t max_local_dimension() const
3500
/// Return the number of dofs on each cell facet
3501
virtual std::size_t num_facet_dofs() const
3506
/// Return the number of dofs associated with each cell entity of dimension d
3507
virtual std::size_t num_entity_dofs(std::size_t d) const
3531
/// Tabulate the local-to-global mapping of dofs on a cell
3532
virtual void tabulate_dofs(std::size_t* dofs,
3533
const std::vector<std::size_t>& num_global_entities,
3534
const ufc::cell& c) const
3536
unsigned int offset = 0;
3537
dofs[0] = offset + c.entity_indices[1][0];
3538
dofs[1] = offset + c.entity_indices[1][1];
3539
dofs[2] = offset + c.entity_indices[1][2];
3540
offset += num_global_entities[1];
3541
dofs[3] = offset + c.entity_indices[2][0];
3542
offset += num_global_entities[2];
3545
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
3546
virtual void tabulate_facet_dofs(std::size_t* dofs,
3547
std::size_t facet) const
3570
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
3571
virtual void tabulate_entity_dofs(std::size_t* dofs,
3572
std::size_t d, std::size_t i) const
3576
std::cerr << "*** FFC warning: " << "d is larger than dimension (2)" << std::endl;
3590
std::cerr << "*** FFC warning: " << "i is larger than number of entities (2)" << std::endl;
3618
std::cerr << "*** FFC warning: " << "i is larger than number of entities (0)" << std::endl;
3628
/// Tabulate the coordinates of all dofs on a cell
3629
virtual void tabulate_coordinates(double** dof_coordinates,
3630
const double* vertex_coordinates) const
3632
dof_coordinates[0][0] = 0.5*vertex_coordinates[3] + 0.5*vertex_coordinates[6];
3633
dof_coordinates[0][1] = 0.5*vertex_coordinates[4] + 0.5*vertex_coordinates[7];
3634
dof_coordinates[0][2] = 0.5*vertex_coordinates[5] + 0.5*vertex_coordinates[8];
3635
dof_coordinates[1][0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[6];
3636
dof_coordinates[1][1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[7];
3637
dof_coordinates[1][2] = 0.5*vertex_coordinates[2] + 0.5*vertex_coordinates[8];
3638
dof_coordinates[2][0] = 0.5*vertex_coordinates[0] + 0.5*vertex_coordinates[3];
3639
dof_coordinates[2][1] = 0.5*vertex_coordinates[1] + 0.5*vertex_coordinates[4];
3640
dof_coordinates[2][2] = 0.5*vertex_coordinates[2] + 0.5*vertex_coordinates[5];
3641
dof_coordinates[3][0] = 0.33333333*vertex_coordinates[0] + 0.33333333*vertex_coordinates[3] + 0.33333333*vertex_coordinates[6];
3642
dof_coordinates[3][1] = 0.33333333*vertex_coordinates[1] + 0.33333333*vertex_coordinates[4] + 0.33333333*vertex_coordinates[7];
3643
dof_coordinates[3][2] = 0.33333333*vertex_coordinates[2] + 0.33333333*vertex_coordinates[5] + 0.33333333*vertex_coordinates[8];
3646
/// Return the number of sub dofmaps (for a mixed element)
3647
virtual std::size_t num_sub_dofmaps() const
3652
/// Create a new dofmap for sub dofmap i (for a mixed element)
3653
virtual ufc::dofmap* create_sub_dofmap(std::size_t i) const
3659
return new projectionmanifold_dofmap_0();
3664
return new projectionmanifold_dofmap_1();
3672
/// Create a new class instance
3673
virtual ufc::dofmap* create() const
3675
return new projectionmanifold_dofmap_2();
3680
/// This class defines the interface for the tabulation of the cell
3681
/// tensor corresponding to the local contribution to a form from
3682
/// the integral over a cell.
3684
class projectionmanifold_cell_integral_0_otherwise: public ufc::cell_integral
3689
projectionmanifold_cell_integral_0_otherwise() : ufc::cell_integral()
3695
virtual ~projectionmanifold_cell_integral_0_otherwise()
3700
/// Tabulate the tensor for the contribution from a local cell
3701
virtual void tabulate_tensor(double* A,
3702
const double * const * w,
3703
const double* vertex_coordinates,
3704
int cell_orientation) const
3706
// Number of operations (multiply-add pairs) for Jacobian data: 5
3707
// Number of operations (multiply-add pairs) for geometry tensor: 66
3708
// Number of operations (multiply-add pairs) for tensor contraction: 136
3709
// Total number of operations (multiply-add pairs): 207
3713
compute_jacobian_triangle_3d(J, vertex_coordinates);
3715
// Compute Jacobian inverse and determinant
3718
compute_jacobian_inverse_triangle_3d(K, detJ, J);
3720
// Check orientation
3721
if (cell_orientation == -1)
3722
throw std::runtime_error("cell orientation must be defined (not -1)");
3723
// (If cell_orientation == 1 = down, multiply det(J) by -1)
3724
else if (cell_orientation == 1)
3728
const double det = std::abs(detJ);
3730
// Compute geometry tensor
3731
const double G0_0_0 = 1.0 / (detJ)*det*J[4]*K[2]*(1.0);
3732
const double G0_1_1 = 1.0 / (detJ)*det*J[5]*K[5]*(1.0);
3733
const double G1_0_0 = 1.0 / (detJ)*det*J[0]*K[0]*(1.0);
3734
const double G1_1_1 = 1.0 / (detJ)*det*J[1]*K[3]*(1.0);
3735
const double G2_0_0 = 1.0 / (detJ)*det*J[2]*K[1]*(1.0);
3736
const double G2_1_1 = 1.0 / (detJ)*det*J[3]*K[4]*(1.0);
3737
const double G3_0_0 = 1.0 / (detJ)*det*J[4]*K[2]*(1.0);
3738
const double G3_1_1 = 1.0 / (detJ)*det*J[5]*K[5]*(1.0);
3739
const double G4_0_0 = 1.0 / (detJ)*det*J[0]*K[0]*(1.0);
3740
const double G4_1_1 = 1.0 / (detJ)*det*J[1]*K[3]*(1.0);
3741
const double G5_0_0 = 1.0 / (detJ)*det*J[2]*K[1]*(1.0);
3742
const double G5_1_1 = 1.0 / (detJ)*det*J[3]*K[4]*(1.0);
3743
const double G6_0_0 = 1.0 / (detJ*detJ)*det*J[4]*J[4]*(1.0);
3744
const double G6_0_1 = 1.0 / (detJ*detJ)*det*J[4]*J[5]*(1.0);
3745
const double G6_1_0 = 1.0 / (detJ*detJ)*det*J[5]*J[4]*(1.0);
3746
const double G6_1_1 = 1.0 / (detJ*detJ)*det*J[5]*J[5]*(1.0);
3747
const double G7_0_0 = 1.0 / (detJ*detJ)*det*J[0]*J[0]*(1.0);
3748
const double G7_0_1 = 1.0 / (detJ*detJ)*det*J[0]*J[1]*(1.0);
3749
const double G7_1_0 = 1.0 / (detJ*detJ)*det*J[1]*J[0]*(1.0);
3750
const double G7_1_1 = 1.0 / (detJ*detJ)*det*J[1]*J[1]*(1.0);
3751
const double G8_0_0 = 1.0 / (detJ*detJ)*det*J[2]*J[2]*(1.0);
3752
const double G8_0_1 = 1.0 / (detJ*detJ)*det*J[2]*J[3]*(1.0);
3753
const double G8_1_0 = 1.0 / (detJ*detJ)*det*J[3]*J[2]*(1.0);
3754
const double G8_1_1 = 1.0 / (detJ*detJ)*det*J[3]*J[3]*(1.0);
3756
// Compute element tensor
3757
A[0] = 0.083333333*G6_0_0 + 0.041666667*G6_0_1 + 0.041666667*G6_1_0 + 0.083333333*G6_1_1 + 0.083333333*G7_0_0 + 0.041666667*G7_0_1 + 0.041666667*G7_1_0 + 0.083333333*G7_1_1 + 0.083333333*G8_0_0 + 0.041666667*G8_0_1 + 0.041666667*G8_1_0 + 0.083333333*G8_1_1;
3758
A[1] = 0.083333333*G6_0_0 - 0.041666667*G6_0_1 + 0.125*G6_1_0 - 0.083333333*G6_1_1 + 0.083333333*G7_0_0 - 0.041666667*G7_0_1 + 0.125*G7_1_0 - 0.083333333*G7_1_1 + 0.083333333*G8_0_0 - 0.041666667*G8_0_1 + 0.125*G8_1_0 - 0.083333333*G8_1_1;
3759
A[2] = 0.083333333*G6_0_0 - 0.125*G6_0_1 + 0.041666667*G6_1_0 - 0.083333333*G6_1_1 + 0.083333333*G7_0_0 - 0.125*G7_0_1 + 0.041666667*G7_1_0 - 0.083333333*G7_1_1 + 0.083333333*G8_0_0 - 0.125*G8_0_1 + 0.041666667*G8_1_0 - 0.083333333*G8_1_1;
3760
A[3] = 0.5*G0_0_0 + 0.5*G0_1_1 + 0.5*G1_0_0 + 0.5*G1_1_1 + 0.5*G2_0_0 + 0.5*G2_1_1;
3761
A[4] = 0.083333333*G6_0_0 + 0.125*G6_0_1 - 0.041666667*G6_1_0 - 0.083333333*G6_1_1 + 0.083333333*G7_0_0 + 0.125*G7_0_1 - 0.041666667*G7_1_0 - 0.083333333*G7_1_1 + 0.083333333*G8_0_0 + 0.125*G8_0_1 - 0.041666667*G8_1_0 - 0.083333333*G8_1_1;
3762
A[5] = 0.25*G6_0_0 - 0.125*G6_0_1 - 0.125*G6_1_0 + 0.083333333*G6_1_1 + 0.25*G7_0_0 - 0.125*G7_0_1 - 0.125*G7_1_0 + 0.083333333*G7_1_1 + 0.25*G8_0_0 - 0.125*G8_0_1 - 0.125*G8_1_0 + 0.083333333*G8_1_1;
3763
A[6] = 0.083333333*G6_0_0 - 0.20833333*G6_0_1 - 0.041666667*G6_1_0 + 0.083333333*G6_1_1 + 0.083333333*G7_0_0 - 0.20833333*G7_0_1 - 0.041666667*G7_1_0 + 0.083333333*G7_1_1 + 0.083333333*G8_0_0 - 0.20833333*G8_0_1 - 0.041666667*G8_1_0 + 0.083333333*G8_1_1;
3764
A[7] = -0.5*G0_0_0 - 0.5*G0_1_1 - 0.5*G1_0_0 - 0.5*G1_1_1 - 0.5*G2_0_0 - 0.5*G2_1_1;
3765
A[8] = 0.083333333*G6_0_0 + 0.041666667*G6_0_1 - 0.125*G6_1_0 - 0.083333333*G6_1_1 + 0.083333333*G7_0_0 + 0.041666667*G7_0_1 - 0.125*G7_1_0 - 0.083333333*G7_1_1 + 0.083333333*G8_0_0 + 0.041666667*G8_0_1 - 0.125*G8_1_0 - 0.083333333*G8_1_1;
3766
A[9] = 0.083333333*G6_0_0 - 0.041666667*G6_0_1 - 0.20833333*G6_1_0 + 0.083333333*G6_1_1 + 0.083333333*G7_0_0 - 0.041666667*G7_0_1 - 0.20833333*G7_1_0 + 0.083333333*G7_1_1 + 0.083333333*G8_0_0 - 0.041666667*G8_0_1 - 0.20833333*G8_1_0 + 0.083333333*G8_1_1;
3767
A[10] = 0.083333333*G6_0_0 - 0.125*G6_0_1 - 0.125*G6_1_0 + 0.25*G6_1_1 + 0.083333333*G7_0_0 - 0.125*G7_0_1 - 0.125*G7_1_0 + 0.25*G7_1_1 + 0.083333333*G8_0_0 - 0.125*G8_0_1 - 0.125*G8_1_0 + 0.25*G8_1_1;
3768
A[11] = 0.5*G0_0_0 + 0.5*G0_1_1 + 0.5*G1_0_0 + 0.5*G1_1_1 + 0.5*G2_0_0 + 0.5*G2_1_1;
3769
A[12] = 0.5*G3_0_0 + 0.5*G3_1_1 + 0.5*G4_0_0 + 0.5*G4_1_1 + 0.5*G5_0_0 + 0.5*G5_1_1;
3770
A[13] = -0.5*G3_0_0 - 0.5*G3_1_1 - 0.5*G4_0_0 - 0.5*G4_1_1 - 0.5*G5_0_0 - 0.5*G5_1_1;
3771
A[14] = 0.5*G3_0_0 + 0.5*G3_1_1 + 0.5*G4_0_0 + 0.5*G4_1_1 + 0.5*G5_0_0 + 0.5*G5_1_1;
3777
/// This class defines the interface for the assembly of the global
3778
/// tensor corresponding to a form with r + n arguments, that is, a
3781
/// a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
3783
/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
3784
/// global tensor A is defined by
3786
/// A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
3788
/// where each argument Vj represents the application to the
3789
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
3790
/// fixed functions (coefficients).
3792
class projectionmanifold_form_0: public ufc::form
3797
projectionmanifold_form_0() : ufc::form()
3803
virtual ~projectionmanifold_form_0()
3808
/// Return a string identifying the form
3809
virtual const char* signature() const
3811
return "e78a2278f8cab860e9ac0865b6fee052698344942c6cbff9e41a848e6f2a49353e4cd71f65096c2facc2f4034f037d7eebb82d42c0e9ecd9a8ef113467c90a0f";
3814
/// Return the rank of the global tensor (r)
3815
virtual std::size_t rank() const
3820
/// Return the number of coefficients (n)
3821
virtual std::size_t num_coefficients() const
3826
/// Return the number of cell domains
3827
virtual std::size_t num_cell_domains() const
3832
/// Return the number of exterior facet domains
3833
virtual std::size_t num_exterior_facet_domains() const
3838
/// Return the number of interior facet domains
3839
virtual std::size_t num_interior_facet_domains() const
3844
/// Return the number of point domains
3845
virtual std::size_t num_point_domains() const
3850
/// Return whether the form has any cell integrals
3851
virtual bool has_cell_integrals() const
3856
/// Return whether the form has any exterior facet integrals
3857
virtual bool has_exterior_facet_integrals() const
3862
/// Return whether the form has any interior facet integrals
3863
virtual bool has_interior_facet_integrals() const
3868
/// Return whether the form has any point integrals
3869
virtual bool has_point_integrals() const
3874
/// Create a new finite element for argument function i
3875
virtual ufc::finite_element* create_finite_element(std::size_t i) const
3881
return new projectionmanifold_finite_element_2();
3886
return new projectionmanifold_finite_element_2();
3894
/// Create a new dofmap for argument function i
3895
virtual ufc::dofmap* create_dofmap(std::size_t i) const
3901
return new projectionmanifold_dofmap_2();
3906
return new projectionmanifold_dofmap_2();
3914
/// Create a new cell integral on sub domain i
3915
virtual ufc::cell_integral* create_cell_integral(std::size_t i) const
3920
/// Create a new exterior facet integral on sub domain i
3921
virtual ufc::exterior_facet_integral* create_exterior_facet_integral(std::size_t i) const
3926
/// Create a new interior facet integral on sub domain i
3927
virtual ufc::interior_facet_integral* create_interior_facet_integral(std::size_t i) const
3932
/// Create a new point integral on sub domain i
3933
virtual ufc::point_integral* create_point_integral(std::size_t i) const
3938
/// Create a new cell integral on everywhere else
3939
virtual ufc::cell_integral* create_default_cell_integral() const
3941
return new projectionmanifold_cell_integral_0_otherwise();
3944
/// Create a new exterior facet integral on everywhere else
3945
virtual ufc::exterior_facet_integral* create_default_exterior_facet_integral() const
3950
/// Create a new interior facet integral on everywhere else
3951
virtual ufc::interior_facet_integral* create_default_interior_facet_integral() const
3956
/// Create a new point integral on everywhere else
3957
virtual ufc::point_integral* create_default_point_integral() const