~ubuntu-branches/ubuntu/trusty/ffc/trusty
» Revision
4
: test/regression/references/TensorWeightedPoisson.h
« back to all changes in this revision
Viewing changes to test/regression/references/TensorWeightedPoisson.h
- browse files at revision 4
- compare with another revision
- download diff
- download tarball
- view history from revision 4
- Committer: Bazaar Package Importer
- Author(s): Johannes Ring
- Date: 2010-02-03 20:22:35 UTC
- mfrom: (1.1.2 upstream)
- Revision ID: james.westby@ubuntu.com-20100203202235-fe8d0kajuvgy2sqn
Tags: 0.9.0-1
* New upstream release.
* debian/control: Bump Standards-Version (no changes needed).
* Update debian/copyright and debian/copyright_hints.
* debian/control: Bump Standards-Version (no changes needed).
* Update debian/copyright and debian/copyright_hints.
- files added:
- demo/VectorPoisson.ufl
- ffc/quadrature
- ffc/tensor
- test/evaluate_basis
- test/evaluate_basis_derivatives
- test/regression/references
- test/unit/evaluate_basis
- test/unit/misc
- test/unit/symbolics
- files removed:
- ffc/common
- ffc/compiler
- ffc/compiler/quadrature
- ffc/compiler/tensor
- ffc/fem
- ffc/jit
- test/benchmarks
- test/regression/reference
- test/regression/reference/quadrature
- test/regression/reference/tensor
- test/simple_verify_tensors
- test/verify_tensors
- test/verify_tensors/references
- test/verify_tensors/test_cache
- files modified:
- AUTHORS
added
removed
test/regression/references/TensorWeightedPoisson.h
1
// This code conforms with the UFC specification version 1.4
2
// and was automatically generated by FFC version 0.9.0.
3
4
#ifndef __TENSORWEIGHTEDPOISSON_H
5
#define __TENSORWEIGHTEDPOISSON_H
6
7
#include <cmath>
8
#include <stdexcept>
9
#include <fstream>
10
#include <ufc.h>
11
12
/// This class defines the interface for a finite element.
13
14
class tensorweightedpoisson_finite_element_0: public ufc::finite_element
15
{
16
public:
17
18
/// Constructor
19
tensorweightedpoisson_finite_element_0() : ufc::finite_element()
20
{
21
// Do nothing
22
}
23
24
/// Destructor
25
virtual ~tensorweightedpoisson_finite_element_0()
26
{
27
// Do nothing
28
}
29
30
/// Return a string identifying the finite element
31
virtual const char* signature() const
32
{
33
return "FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 0)";
34
}
35
36
/// Return the cell shape
37
virtual ufc::shape cell_shape() const
38
{
39
return ufc::triangle;
40
}
41
42
/// Return the dimension of the finite element function space
43
virtual unsigned int space_dimension() const
44
{
45
return 1;
46
}
47
48
/// Return the rank of the value space
49
virtual unsigned int value_rank() const
50
{
51
return 0;
52
}
53
54
/// Return the dimension of the value space for axis i
55
virtual unsigned int value_dimension(unsigned int i) const
56
{
57
return 1;
58
}
59
60
/// Evaluate basis function i at given point in cell
61
virtual void evaluate_basis(unsigned int i,
62
double* values,
63
const double* coordinates,
64
const ufc::cell& c) const
65
{
66
// Extract vertex coordinates
67
68
// Compute Jacobian of affine map from reference cell
69
70
// Compute determinant of Jacobian
71
72
// Compute inverse of Jacobian
73
74
// Compute constants
75
76
// Get coordinates and map to the reference (FIAT) element
77
78
// Reset values
79
*values = 0.00000000;
80
81
// Map degree of freedom to element degree of freedom
82
const unsigned int dof = i;
83
84
// Array of basisvalues
85
double basisvalues[1] = {0.00000000};
86
87
// Declare helper variables
88
89
// Compute basisvalues
90
basisvalues[0] = 1.00000000;
91
92
// Table(s) of coefficients
93
static const double coefficients0[1][1] = \
94
{{1.00000000}};
95
96
// Compute value(s).
97
for (unsigned int r = 0; r < 1; r++)
98
{
99
*values += coefficients0[dof][r]*basisvalues[r];
100
}// end loop over 'r'
101
}
102
103
/// Evaluate all basis functions at given point in cell
104
virtual void evaluate_basis_all(double* values,
105
const double* coordinates,
106
const ufc::cell& c) const
107
{
108
// Element is constant, calling evaluate_basis.
109
evaluate_basis(0, values, coordinates, c);
110
}
111
112
/// Evaluate order n derivatives of basis function i at given point in cell
113
virtual void evaluate_basis_derivatives(unsigned int i,
114
unsigned int n,
115
double* values,
116
const double* coordinates,
117
const ufc::cell& c) const
118
{
119
// Extract vertex coordinates
120
const double * const * x = c.coordinates;
121
122
// Compute Jacobian of affine map from reference cell
123
const double J_00 = x[1][0] - x[0][0];
124
const double J_01 = x[2][0] - x[0][0];
125
const double J_10 = x[1][1] - x[0][1];
126
const double J_11 = x[2][1] - x[0][1];
127
128
// Compute determinant of Jacobian
129
double detJ = J_00*J_11 - J_01*J_10;
130
131
// Compute inverse of Jacobian
132
const double K_00 = J_11 / detJ;
133
const double K_01 = -J_01 / detJ;
134
const double K_10 = -J_10 / detJ;
135
const double K_11 = J_00 / detJ;
136
137
// Compute constants
138
139
// Get coordinates and map to the reference (FIAT) element
140
141
// Compute number of derivatives.
142
unsigned int num_derivatives = 1;
143
144
for (unsigned int r = 0; r < n; r++)
145
{
146
num_derivatives *= 2;
147
}// end loop over 'r'
148
149
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
150
unsigned int **combinations = new unsigned int *[num_derivatives];
151
for (unsigned int row = 0; row < num_derivatives; row++)
152
{
153
combinations[row] = new unsigned int [n];
154
for (unsigned int col = 0; col < n; col++)
155
combinations[row][col] = 0;
156
}
157
158
// Generate combinations of derivatives
159
for (unsigned int row = 1; row < num_derivatives; row++)
160
{
161
for (unsigned int num = 0; num < row; num++)
162
{
163
for (unsigned int col = n-1; col+1 > 0; col--)
164
{
165
if (combinations[row][col] + 1 > 1)
166
combinations[row][col] = 0;
167
else
168
{
169
combinations[row][col] += 1;
170
break;
171
}
172
}
173
}
174
}
175
176
// Compute inverse of Jacobian
177
const double Jinv[2][2] = {{K_00, K_01}, {K_10, K_11}};
178
179
// Declare transformation matrix
180
// Declare pointer to two dimensional array and initialise
181
double **transform = new double *[num_derivatives];
182
183
for (unsigned int j = 0; j < num_derivatives; j++)
184
{
185
transform[j] = new double [num_derivatives];
186
for (unsigned int k = 0; k < num_derivatives; k++)
187
transform[j][k] = 1;
188
}
189
190
// Construct transformation matrix
191
for (unsigned int row = 0; row < num_derivatives; row++)
192
{
193
for (unsigned int col = 0; col < num_derivatives; col++)
194
{
195
for (unsigned int k = 0; k < n; k++)
196
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
197
}
198
}
199
200
// Reset values. Assuming that values is always an array.
201
for (unsigned int r = 0; r < num_derivatives; r++)
202
{
203
values[r] = 0.00000000;
204
}// end loop over 'r'
205
206
// Map degree of freedom to element degree of freedom
207
const unsigned int dof = i;
208
209
// Array of basisvalues
210
double basisvalues[1] = {0.00000000};
211
212
// Declare helper variables
213
214
// Compute basisvalues
215
basisvalues[0] = 1.00000000;
216
217
// Table(s) of coefficients
218
static const double coefficients0[1][1] = \
219
{{1.00000000}};
220
221
// Tables of derivatives of the polynomial base (transpose).
222
static const double dmats0[1][1] = \
223
{{0.00000000}};
224
225
static const double dmats1[1][1] = \
226
{{0.00000000}};
227
228
// Compute reference derivatives
229
// Declare pointer to array of derivatives on FIAT element
230
double *derivatives = new double [num_derivatives];
231
for (unsigned int r = 0; r < num_derivatives; r++)
232
{
233
derivatives[r] = 0.00000000;
234
}// end loop over 'r'
235
236
// Declare derivative matrix (of polynomial basis).
237
double dmats[1][1] = \
238
{{1.00000000}};
239
240
// Declare (auxiliary) derivative matrix (of polynomial basis).
241
double dmats_old[1][1] = \
242
{{1.00000000}};
243
244
// Loop possible derivatives.
245
for (unsigned int r = 0; r < num_derivatives; r++)
246
{
247
// Resetting dmats values to compute next derivative.
248
for (unsigned int t = 0; t < 1; t++)
249
{
250
for (unsigned int u = 0; u < 1; u++)
251
{
252
dmats[t][u] = 0.00000000;
253
if (t == u)
254
{
255
dmats[t][u] = 1.00000000;
256
}
257
258
}// end loop over 'u'
259
}// end loop over 't'
260
261
// Looping derivative order to generate dmats.
262
for (unsigned int s = 0; s < n; s++)
263
{
264
// Updating dmats_old with new values and resetting dmats.
265
for (unsigned int t = 0; t < 1; t++)
266
{
267
for (unsigned int u = 0; u < 1; u++)
268
{
269
dmats_old[t][u] = dmats[t][u];
270
dmats[t][u] = 0.00000000;
271
}// end loop over 'u'
272
}// end loop over 't'
273
274
// Update dmats using an inner product.
275
if (combinations[r][s] == 0)
276
{
277
for (unsigned int t = 0; t < 1; t++)
278
{
279
for (unsigned int u = 0; u < 1; u++)
280
{
281
for (unsigned int tu = 0; tu < 1; tu++)
282
{
283
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
284
}// end loop over 'tu'
285
}// end loop over 'u'
286
}// end loop over 't'
287
}
288
289
if (combinations[r][s] == 1)
290
{
291
for (unsigned int t = 0; t < 1; t++)
292
{
293
for (unsigned int u = 0; u < 1; u++)
294
{
295
for (unsigned int tu = 0; tu < 1; tu++)
296
{
297
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
298
}// end loop over 'tu'
299
}// end loop over 'u'
300
}// end loop over 't'
301
}
302
303
}// end loop over 's'
304
for (unsigned int s = 0; s < 1; s++)
305
{
306
for (unsigned int t = 0; t < 1; t++)
307
{
308
derivatives[r] += coefficients0[dof][s]*dmats[s][t]*basisvalues[t];
309
}// end loop over 't'
310
}// end loop over 's'
311
}// end loop over 'r'
312
313
// Transform derivatives back to physical element
314
for (unsigned int row = 0; row < num_derivatives; row++)
315
{
316
for (unsigned int col = 0; col < num_derivatives; col++)
317
{
318
values[row] += transform[row][col]*derivatives[col];
319
}
320
}
321
322
// Delete pointer to array of derivatives on FIAT element
323
delete [] derivatives;
324
325
// Delete pointer to array of combinations of derivatives and transform
326
for (unsigned int r = 0; r < num_derivatives; r++)
327
{
328
delete [] combinations[r];
329
delete [] transform[r];
330
}// end loop over 'r'
331
delete [] combinations;
332
delete [] transform;
333
}
334
335
/// Evaluate order n derivatives of all basis functions at given point in cell
336
virtual void evaluate_basis_derivatives_all(unsigned int n,
337
double* values,
338
const double* coordinates,
339
const ufc::cell& c) const
340
{
341
// Element is constant, calling evaluate_basis_derivatives.
342
evaluate_basis_derivatives(0, n, values, coordinates, c);
343
}
344
345
/// Evaluate linear functional for dof i on the function f
346
virtual double evaluate_dof(unsigned int i,
347
const ufc::function& f,
348
const ufc::cell& c) const
349
{
350
// Declare variables for result of evaluation
351
double vals[1];
352
353
// Declare variable for physical coordinates
354
double y[2];
355
356
const double * const * x = c.coordinates;
357
switch (i)
358
{
359
case 0:
360
{
361
y[0] = 0.333333333333*x[0][0] + 0.333333333333*x[1][0] + 0.333333333333*x[2][0];
362
y[1] = 0.333333333333*x[0][1] + 0.333333333333*x[1][1] + 0.333333333333*x[2][1];
363
f.evaluate(vals, y, c);
364
return vals[0];
365
break;
366
}
367
}
368
369
return 0.0;
370
}
371
372
/// Evaluate linear functionals for all dofs on the function f
373
virtual void evaluate_dofs(double* values,
374
const ufc::function& f,
375
const ufc::cell& c) const
376
{
377
// Declare variables for result of evaluation
378
double vals[1];
379
380
// Declare variable for physical coordinates
381
double y[2];
382
383
const double * const * x = c.coordinates;
384
y[0] = 0.333333333333*x[0][0] + 0.333333333333*x[1][0] + 0.333333333333*x[2][0];
385
y[1] = 0.333333333333*x[0][1] + 0.333333333333*x[1][1] + 0.333333333333*x[2][1];
386
f.evaluate(vals, y, c);
387
values[0] = vals[0];
388
}
389
390
/// Interpolate vertex values from dof values
391
virtual void interpolate_vertex_values(double* vertex_values,
392
const double* dof_values,
393
const ufc::cell& c) const
394
{
395
// Evaluate function and change variables
396
vertex_values[0] = dof_values[0];
397
vertex_values[1] = dof_values[0];
398
vertex_values[2] = dof_values[0];
399
}
400
401
/// Return the number of sub elements (for a mixed element)
402
virtual unsigned int num_sub_elements() const
403
{
404
return 0;
405
}
406
407
/// Create a new finite element for sub element i (for a mixed element)
408
virtual ufc::finite_element* create_sub_element(unsigned int i) const
409
{
410
return 0;
411
}
412
413
};
414
415
/// This class defines the interface for a finite element.
416
417
class tensorweightedpoisson_finite_element_1: public ufc::finite_element
418
{
419
public:
420
421
/// Constructor
422
tensorweightedpoisson_finite_element_1() : ufc::finite_element()
423
{
424
// Do nothing
425
}
426
427
/// Destructor
428
virtual ~tensorweightedpoisson_finite_element_1()
429
{
430
// Do nothing
431
}
432
433
/// Return a string identifying the finite element
434
virtual const char* signature() const
435
{
436
return "TensorElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 0, (2, 2), None)";
437
}
438
439
/// Return the cell shape
440
virtual ufc::shape cell_shape() const
441
{
442
return ufc::triangle;
443
}
444
445
/// Return the dimension of the finite element function space
446
virtual unsigned int space_dimension() const
447
{
448
return 4;
449
}
450
451
/// Return the rank of the value space
452
virtual unsigned int value_rank() const
453
{
454
return 2;
455
}
456
457
/// Return the dimension of the value space for axis i
458
virtual unsigned int value_dimension(unsigned int i) const
459
{
460
switch (i)
461
{
462
case 0:
463
{
464
return 2;
465
break;
466
}
467
case 1:
468
{
469
return 2;
470
break;
471
}
472
}
473
474
return 0;
475
}
476
477
/// Evaluate basis function i at given point in cell
478
virtual void evaluate_basis(unsigned int i,
479
double* values,
480
const double* coordinates,
481
const ufc::cell& c) const
482
{
483
// Extract vertex coordinates
484
485
// Compute Jacobian of affine map from reference cell
486
487
// Compute determinant of Jacobian
488
489
// Compute inverse of Jacobian
490
491
// Compute constants
492
493
// Get coordinates and map to the reference (FIAT) element
494
495
// Reset values
496
values[0] = 0.00000000;
497
values[1] = 0.00000000;
498
values[2] = 0.00000000;
499
values[3] = 0.00000000;
500
if (0 <= i && i <= 0)
501
{
502
// Map degree of freedom to element degree of freedom
503
const unsigned int dof = i;
504
505
// Array of basisvalues
506
double basisvalues[1] = {0.00000000};
507
508
// Declare helper variables
509
510
// Compute basisvalues
511
basisvalues[0] = 1.00000000;
512
513
// Table(s) of coefficients
514
static const double coefficients0[1][1] = \
515
{{1.00000000}};
516
517
// Compute value(s).
518
for (unsigned int r = 0; r < 1; r++)
519
{
520
values[0] += coefficients0[dof][r]*basisvalues[r];
521
}// end loop over 'r'
522
}
523
524
if (1 <= i && i <= 1)
525
{
526
// Map degree of freedom to element degree of freedom
527
const unsigned int dof = i - 1;
528
529
// Array of basisvalues
530
double basisvalues[1] = {0.00000000};
531
532
// Declare helper variables
533
534
// Compute basisvalues
535
basisvalues[0] = 1.00000000;
536
537
// Table(s) of coefficients
538
static const double coefficients0[1][1] = \
539
{{1.00000000}};
540
541
// Compute value(s).
542
for (unsigned int r = 0; r < 1; r++)
543
{
544
values[1] += coefficients0[dof][r]*basisvalues[r];
545
}// end loop over 'r'
546
}
547
548
if (2 <= i && i <= 2)
549
{
550
// Map degree of freedom to element degree of freedom
551
const unsigned int dof = i - 2;
552
553
// Array of basisvalues
554
double basisvalues[1] = {0.00000000};
555
556
// Declare helper variables
557
558
// Compute basisvalues
559
basisvalues[0] = 1.00000000;
560
561
// Table(s) of coefficients
562
static const double coefficients0[1][1] = \
563
{{1.00000000}};
564
565
// Compute value(s).
566
for (unsigned int r = 0; r < 1; r++)
567
{
568
values[2] += coefficients0[dof][r]*basisvalues[r];
569
}// end loop over 'r'
570
}
571
572
if (3 <= i && i <= 3)
573
{
574
// Map degree of freedom to element degree of freedom
575
const unsigned int dof = i - 3;
576
577
// Array of basisvalues
578
double basisvalues[1] = {0.00000000};
579
580
// Declare helper variables
581
582
// Compute basisvalues
583
basisvalues[0] = 1.00000000;
584
585
// Table(s) of coefficients
586
static const double coefficients0[1][1] = \
587
{{1.00000000}};
588
589
// Compute value(s).
590
for (unsigned int r = 0; r < 1; r++)
591
{
592
values[3] += coefficients0[dof][r]*basisvalues[r];
593
}// end loop over 'r'
594
}
595
596
}
597
598
/// Evaluate all basis functions at given point in cell
599
virtual void evaluate_basis_all(double* values,
600
const double* coordinates,
601
const ufc::cell& c) const
602
{
603
// Helper variable to hold values of a single dof.
604
double dof_values[4] = {0.00000000, 0.00000000, 0.00000000, 0.00000000};
605
606
// Loop dofs and call evaluate_basis.
607
for (unsigned int r = 0; r < 4; r++)
608
{
609
evaluate_basis(r, dof_values, coordinates, c);
610
for (unsigned int s = 0; s < 4; s++)
611
{
612
values[r*4 + s] = dof_values[s];
613
}// end loop over 's'
614
}// end loop over 'r'
615
}
616
617
/// Evaluate order n derivatives of basis function i at given point in cell
618
virtual void evaluate_basis_derivatives(unsigned int i,
619
unsigned int n,
620
double* values,
621
const double* coordinates,
622
const ufc::cell& c) const
623
{
624
// Extract vertex coordinates
625
const double * const * x = c.coordinates;
626
627
// Compute Jacobian of affine map from reference cell
628
const double J_00 = x[1][0] - x[0][0];
629
const double J_01 = x[2][0] - x[0][0];
630
const double J_10 = x[1][1] - x[0][1];
631
const double J_11 = x[2][1] - x[0][1];
632
633
// Compute determinant of Jacobian
634
double detJ = J_00*J_11 - J_01*J_10;
635
636
// Compute inverse of Jacobian
637
const double K_00 = J_11 / detJ;
638
const double K_01 = -J_01 / detJ;
639
const double K_10 = -J_10 / detJ;
640
const double K_11 = J_00 / detJ;
641
642
// Compute constants
643
644
// Get coordinates and map to the reference (FIAT) element
645
646
// Compute number of derivatives.
647
unsigned int num_derivatives = 1;
648
649
for (unsigned int r = 0; r < n; r++)
650
{
651
num_derivatives *= 2;
652
}// end loop over 'r'
653
654
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
655
unsigned int **combinations = new unsigned int *[num_derivatives];
656
for (unsigned int row = 0; row < num_derivatives; row++)
657
{
658
combinations[row] = new unsigned int [n];
659
for (unsigned int col = 0; col < n; col++)
660
combinations[row][col] = 0;
661
}
662
663
// Generate combinations of derivatives
664
for (unsigned int row = 1; row < num_derivatives; row++)
665
{
666
for (unsigned int num = 0; num < row; num++)
667
{
668
for (unsigned int col = n-1; col+1 > 0; col--)
669
{
670
if (combinations[row][col] + 1 > 1)
671
combinations[row][col] = 0;
672
else
673
{
674
combinations[row][col] += 1;
675
break;
676
}
677
}
678
}
679
}
680
681
// Compute inverse of Jacobian
682
const double Jinv[2][2] = {{K_00, K_01}, {K_10, K_11}};
683
684
// Declare transformation matrix
685
// Declare pointer to two dimensional array and initialise
686
double **transform = new double *[num_derivatives];
687
688
for (unsigned int j = 0; j < num_derivatives; j++)
689
{
690
transform[j] = new double [num_derivatives];
691
for (unsigned int k = 0; k < num_derivatives; k++)
692
transform[j][k] = 1;
693
}
694
695
// Construct transformation matrix
696
for (unsigned int row = 0; row < num_derivatives; row++)
697
{
698
for (unsigned int col = 0; col < num_derivatives; col++)
699
{
700
for (unsigned int k = 0; k < n; k++)
701
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
702
}
703
}
704
705
// Reset values. Assuming that values is always an array.
706
for (unsigned int r = 0; r < 4*num_derivatives; r++)
707
{
708
values[r] = 0.00000000;
709
}// end loop over 'r'
710
711
if (0 <= i && i <= 0)
712
{
713
// Map degree of freedom to element degree of freedom
714
const unsigned int dof = i;
715
716
// Array of basisvalues
717
double basisvalues[1] = {0.00000000};
718
719
// Declare helper variables
720
721
// Compute basisvalues
722
basisvalues[0] = 1.00000000;
723
724
// Table(s) of coefficients
725
static const double coefficients0[1][1] = \
726
{{1.00000000}};
727
728
// Tables of derivatives of the polynomial base (transpose).
729
static const double dmats0[1][1] = \
730
{{0.00000000}};
731
732
static const double dmats1[1][1] = \
733
{{0.00000000}};
734
735
// Compute reference derivatives
736
// Declare pointer to array of derivatives on FIAT element
737
double *derivatives = new double [num_derivatives];
738
for (unsigned int r = 0; r < num_derivatives; r++)
739
{
740
derivatives[r] = 0.00000000;
741
}// end loop over 'r'
742
743
// Declare derivative matrix (of polynomial basis).
744
double dmats[1][1] = \
745
{{1.00000000}};
746
747
// Declare (auxiliary) derivative matrix (of polynomial basis).
748
double dmats_old[1][1] = \
749
{{1.00000000}};
750
751
// Loop possible derivatives.
752
for (unsigned int r = 0; r < num_derivatives; r++)
753
{
754
// Resetting dmats values to compute next derivative.
755
for (unsigned int t = 0; t < 1; t++)
756
{
757
for (unsigned int u = 0; u < 1; u++)
758
{
759
dmats[t][u] = 0.00000000;
760
if (t == u)
761
{
762
dmats[t][u] = 1.00000000;
763
}
764
765
}// end loop over 'u'
766
}// end loop over 't'
767
768
// Looping derivative order to generate dmats.
769
for (unsigned int s = 0; s < n; s++)
770
{
771
// Updating dmats_old with new values and resetting dmats.
772
for (unsigned int t = 0; t < 1; t++)
773
{
774
for (unsigned int u = 0; u < 1; u++)
775
{
776
dmats_old[t][u] = dmats[t][u];
777
dmats[t][u] = 0.00000000;
778
}// end loop over 'u'
779
}// end loop over 't'
780
781
// Update dmats using an inner product.
782
if (combinations[r][s] == 0)
783
{
784
for (unsigned int t = 0; t < 1; t++)
785
{
786
for (unsigned int u = 0; u < 1; u++)
787
{
788
for (unsigned int tu = 0; tu < 1; tu++)
789
{
790
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
791
}// end loop over 'tu'
792
}// end loop over 'u'
793
}// end loop over 't'
794
}
795
796
if (combinations[r][s] == 1)
797
{
798
for (unsigned int t = 0; t < 1; t++)
799
{
800
for (unsigned int u = 0; u < 1; u++)
801
{
802
for (unsigned int tu = 0; tu < 1; tu++)
803
{
804
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
805
}// end loop over 'tu'
806
}// end loop over 'u'
807
}// end loop over 't'
808
}
809
810
}// end loop over 's'
811
for (unsigned int s = 0; s < 1; s++)
812
{
813
for (unsigned int t = 0; t < 1; t++)
814
{
815
derivatives[r] += coefficients0[dof][s]*dmats[s][t]*basisvalues[t];
816
}// end loop over 't'
817
}// end loop over 's'
818
}// end loop over 'r'
819
820
// Transform derivatives back to physical element
821
for (unsigned int row = 0; row < num_derivatives; row++)
822
{
823
for (unsigned int col = 0; col < num_derivatives; col++)
824
{
825
values[row] += transform[row][col]*derivatives[col];
826
}
827
}
828
829
// Delete pointer to array of derivatives on FIAT element
830
delete [] derivatives;
831
832
// Delete pointer to array of combinations of derivatives and transform
833
for (unsigned int r = 0; r < num_derivatives; r++)
834
{
835
delete [] combinations[r];
836
delete [] transform[r];
837
}// end loop over 'r'
838
delete [] combinations;
839
delete [] transform;
840
}
841
842
if (1 <= i && i <= 1)
843
{
844
// Map degree of freedom to element degree of freedom
845
const unsigned int dof = i - 1;
846
847
// Array of basisvalues
848
double basisvalues[1] = {0.00000000};
849
850
// Declare helper variables
851
852
// Compute basisvalues
853
basisvalues[0] = 1.00000000;
854
855
// Table(s) of coefficients
856
static const double coefficients0[1][1] = \
857
{{1.00000000}};
858
859
// Tables of derivatives of the polynomial base (transpose).
860
static const double dmats0[1][1] = \
861
{{0.00000000}};
862
863
static const double dmats1[1][1] = \
864
{{0.00000000}};
865
866
// Compute reference derivatives
867
// Declare pointer to array of derivatives on FIAT element
868
double *derivatives = new double [num_derivatives];
869
for (unsigned int r = 0; r < num_derivatives; r++)
870
{
871
derivatives[r] = 0.00000000;
872
}// end loop over 'r'
873
874
// Declare derivative matrix (of polynomial basis).
875
double dmats[1][1] = \
876
{{1.00000000}};
877
878
// Declare (auxiliary) derivative matrix (of polynomial basis).
879
double dmats_old[1][1] = \
880
{{1.00000000}};
881
882
// Loop possible derivatives.
883
for (unsigned int r = 0; r < num_derivatives; r++)
884
{
885
// Resetting dmats values to compute next derivative.
886
for (unsigned int t = 0; t < 1; t++)
887
{
888
for (unsigned int u = 0; u < 1; u++)
889
{
890
dmats[t][u] = 0.00000000;
891
if (t == u)
892
{
893
dmats[t][u] = 1.00000000;
894
}
895
896
}// end loop over 'u'
897
}// end loop over 't'
898
899
// Looping derivative order to generate dmats.
900
for (unsigned int s = 0; s < n; s++)
901
{
902
// Updating dmats_old with new values and resetting dmats.
903
for (unsigned int t = 0; t < 1; t++)
904
{
905
for (unsigned int u = 0; u < 1; u++)
906
{
907
dmats_old[t][u] = dmats[t][u];
908
dmats[t][u] = 0.00000000;
909
}// end loop over 'u'
910
}// end loop over 't'
911
912
// Update dmats using an inner product.
913
if (combinations[r][s] == 0)
914
{
915
for (unsigned int t = 0; t < 1; t++)
916
{
917
for (unsigned int u = 0; u < 1; u++)
918
{
919
for (unsigned int tu = 0; tu < 1; tu++)
920
{
921
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
922
}// end loop over 'tu'
923
}// end loop over 'u'
924
}// end loop over 't'
925
}
926
927
if (combinations[r][s] == 1)
928
{
929
for (unsigned int t = 0; t < 1; t++)
930
{
931
for (unsigned int u = 0; u < 1; u++)
932
{
933
for (unsigned int tu = 0; tu < 1; tu++)
934
{
935
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
936
}// end loop over 'tu'
937
}// end loop over 'u'
938
}// end loop over 't'
939
}
940
941
}// end loop over 's'
942
for (unsigned int s = 0; s < 1; s++)
943
{
944
for (unsigned int t = 0; t < 1; t++)
945
{
946
derivatives[r] += coefficients0[dof][s]*dmats[s][t]*basisvalues[t];
947
}// end loop over 't'
948
}// end loop over 's'
949
}// end loop over 'r'
950
951
// Transform derivatives back to physical element
952
for (unsigned int row = 0; row < num_derivatives; row++)
953
{
954
for (unsigned int col = 0; col < num_derivatives; col++)
955
{
956
values[num_derivatives + row] += transform[row][col]*derivatives[col];
957
}
958
}
959
960
// Delete pointer to array of derivatives on FIAT element
961
delete [] derivatives;
962
963
// Delete pointer to array of combinations of derivatives and transform
964
for (unsigned int r = 0; r < num_derivatives; r++)
965
{
966
delete [] combinations[r];
967
delete [] transform[r];
968
}// end loop over 'r'
969
delete [] combinations;
970
delete [] transform;
971
}
972
973
if (2 <= i && i <= 2)
974
{
975
// Map degree of freedom to element degree of freedom
976
const unsigned int dof = i - 2;
977
978
// Array of basisvalues
979
double basisvalues[1] = {0.00000000};
980
981
// Declare helper variables
982
983
// Compute basisvalues
984
basisvalues[0] = 1.00000000;
985
986
// Table(s) of coefficients
987
static const double coefficients0[1][1] = \
988
{{1.00000000}};
989
990
// Tables of derivatives of the polynomial base (transpose).
991
static const double dmats0[1][1] = \
992
{{0.00000000}};
993
994
static const double dmats1[1][1] = \
995
{{0.00000000}};
996
997
// Compute reference derivatives
998
// Declare pointer to array of derivatives on FIAT element
999
double *derivatives = new double [num_derivatives];
1000
for (unsigned int r = 0; r < num_derivatives; r++)
1001
{
1002
derivatives[r] = 0.00000000;
1003
}// end loop over 'r'
1004
1005
// Declare derivative matrix (of polynomial basis).
1006
double dmats[1][1] = \
1007
{{1.00000000}};
1008
1009
// Declare (auxiliary) derivative matrix (of polynomial basis).
1010
double dmats_old[1][1] = \
1011
{{1.00000000}};
1012
1013
// Loop possible derivatives.
1014
for (unsigned int r = 0; r < num_derivatives; r++)
1015
{
1016
// Resetting dmats values to compute next derivative.
1017
for (unsigned int t = 0; t < 1; t++)
1018
{
1019
for (unsigned int u = 0; u < 1; u++)
1020
{
1021
dmats[t][u] = 0.00000000;
1022
if (t == u)
1023
{
1024
dmats[t][u] = 1.00000000;
1025
}
1026
1027
}// end loop over 'u'
1028
}// end loop over 't'
1029
1030
// Looping derivative order to generate dmats.
1031
for (unsigned int s = 0; s < n; s++)
1032
{
1033
// Updating dmats_old with new values and resetting dmats.
1034
for (unsigned int t = 0; t < 1; t++)
1035
{
1036
for (unsigned int u = 0; u < 1; u++)
1037
{
1038
dmats_old[t][u] = dmats[t][u];
1039
dmats[t][u] = 0.00000000;
1040
}// end loop over 'u'
1041
}// end loop over 't'
1042
1043
// Update dmats using an inner product.
1044
if (combinations[r][s] == 0)
1045
{
1046
for (unsigned int t = 0; t < 1; t++)
1047
{
1048
for (unsigned int u = 0; u < 1; u++)
1049
{
1050
for (unsigned int tu = 0; tu < 1; tu++)
1051
{
1052
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1053
}// end loop over 'tu'
1054
}// end loop over 'u'
1055
}// end loop over 't'
1056
}
1057
1058
if (combinations[r][s] == 1)
1059
{
1060
for (unsigned int t = 0; t < 1; t++)
1061
{
1062
for (unsigned int u = 0; u < 1; u++)
1063
{
1064
for (unsigned int tu = 0; tu < 1; tu++)
1065
{
1066
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1067
}// end loop over 'tu'
1068
}// end loop over 'u'
1069
}// end loop over 't'
1070
}
1071
1072
}// end loop over 's'
1073
for (unsigned int s = 0; s < 1; s++)
1074
{
1075
for (unsigned int t = 0; t < 1; t++)
1076
{
1077
derivatives[r] += coefficients0[dof][s]*dmats[s][t]*basisvalues[t];
1078
}// end loop over 't'
1079
}// end loop over 's'
1080
}// end loop over 'r'
1081
1082
// Transform derivatives back to physical element
1083
for (unsigned int row = 0; row < num_derivatives; row++)
1084
{
1085
for (unsigned int col = 0; col < num_derivatives; col++)
1086
{
1087
values[2*num_derivatives + row] += transform[row][col]*derivatives[col];
1088
}
1089
}
1090
1091
// Delete pointer to array of derivatives on FIAT element
1092
delete [] derivatives;
1093
1094
// Delete pointer to array of combinations of derivatives and transform
1095
for (unsigned int r = 0; r < num_derivatives; r++)
1096
{
1097
delete [] combinations[r];
1098
delete [] transform[r];
1099
}// end loop over 'r'
1100
delete [] combinations;
1101
delete [] transform;
1102
}
1103
1104
if (3 <= i && i <= 3)
1105
{
1106
// Map degree of freedom to element degree of freedom
1107
const unsigned int dof = i - 3;
1108
1109
// Array of basisvalues
1110
double basisvalues[1] = {0.00000000};
1111
1112
// Declare helper variables
1113
1114
// Compute basisvalues
1115
basisvalues[0] = 1.00000000;
1116
1117
// Table(s) of coefficients
1118
static const double coefficients0[1][1] = \
1119
{{1.00000000}};
1120
1121
// Tables of derivatives of the polynomial base (transpose).
1122
static const double dmats0[1][1] = \
1123
{{0.00000000}};
1124
1125
static const double dmats1[1][1] = \
1126
{{0.00000000}};
1127
1128
// Compute reference derivatives
1129
// Declare pointer to array of derivatives on FIAT element
1130
double *derivatives = new double [num_derivatives];
1131
for (unsigned int r = 0; r < num_derivatives; r++)
1132
{
1133
derivatives[r] = 0.00000000;
1134
}// end loop over 'r'
1135
1136
// Declare derivative matrix (of polynomial basis).
1137
double dmats[1][1] = \
1138
{{1.00000000}};
1139
1140
// Declare (auxiliary) derivative matrix (of polynomial basis).
1141
double dmats_old[1][1] = \
1142
{{1.00000000}};
1143
1144
// Loop possible derivatives.
1145
for (unsigned int r = 0; r < num_derivatives; r++)
1146
{
1147
// Resetting dmats values to compute next derivative.
1148
for (unsigned int t = 0; t < 1; t++)
1149
{
1150
for (unsigned int u = 0; u < 1; u++)
1151
{
1152
dmats[t][u] = 0.00000000;
1153
if (t == u)
1154
{
1155
dmats[t][u] = 1.00000000;
1156
}
1157
1158
}// end loop over 'u'
1159
}// end loop over 't'
1160
1161
// Looping derivative order to generate dmats.
1162
for (unsigned int s = 0; s < n; s++)
1163
{
1164
// Updating dmats_old with new values and resetting dmats.
1165
for (unsigned int t = 0; t < 1; t++)
1166
{
1167
for (unsigned int u = 0; u < 1; u++)
1168
{
1169
dmats_old[t][u] = dmats[t][u];
1170
dmats[t][u] = 0.00000000;
1171
}// end loop over 'u'
1172
}// end loop over 't'
1173
1174
// Update dmats using an inner product.
1175
if (combinations[r][s] == 0)
1176
{
1177
for (unsigned int t = 0; t < 1; t++)
1178
{
1179
for (unsigned int u = 0; u < 1; u++)
1180
{
1181
for (unsigned int tu = 0; tu < 1; tu++)
1182
{
1183
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1184
}// end loop over 'tu'
1185
}// end loop over 'u'
1186
}// end loop over 't'
1187
}
1188
1189
if (combinations[r][s] == 1)
1190
{
1191
for (unsigned int t = 0; t < 1; t++)
1192
{
1193
for (unsigned int u = 0; u < 1; u++)
1194
{
1195
for (unsigned int tu = 0; tu < 1; tu++)
1196
{
1197
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1198
}// end loop over 'tu'
1199
}// end loop over 'u'
1200
}// end loop over 't'
1201
}
1202
1203
}// end loop over 's'
1204
for (unsigned int s = 0; s < 1; s++)
1205
{
1206
for (unsigned int t = 0; t < 1; t++)
1207
{
1208
derivatives[r] += coefficients0[dof][s]*dmats[s][t]*basisvalues[t];
1209
}// end loop over 't'
1210
}// end loop over 's'
1211
}// end loop over 'r'
1212
1213
// Transform derivatives back to physical element
1214
for (unsigned int row = 0; row < num_derivatives; row++)
1215
{
1216
for (unsigned int col = 0; col < num_derivatives; col++)
1217
{
1218
values[3*num_derivatives + row] += transform[row][col]*derivatives[col];
1219
}
1220
}
1221
1222
// Delete pointer to array of derivatives on FIAT element
1223
delete [] derivatives;
1224
1225
// Delete pointer to array of combinations of derivatives and transform
1226
for (unsigned int r = 0; r < num_derivatives; r++)
1227
{
1228
delete [] combinations[r];
1229
delete [] transform[r];
1230
}// end loop over 'r'
1231
delete [] combinations;
1232
delete [] transform;
1233
}
1234
1235
}
1236
1237
/// Evaluate order n derivatives of all basis functions at given point in cell
1238
virtual void evaluate_basis_derivatives_all(unsigned int n,
1239
double* values,
1240
const double* coordinates,
1241
const ufc::cell& c) const
1242
{
1243
// Compute number of derivatives.
1244
unsigned int num_derivatives = 1;
1245
1246
for (unsigned int r = 0; r < n; r++)
1247
{
1248
num_derivatives *= 2;
1249
}// end loop over 'r'
1250
1251
// Helper variable to hold values of a single dof.
1252
double *dof_values = new double [4*num_derivatives];
1253
for (unsigned int r = 0; r < 4*num_derivatives; r++)
1254
{
1255
dof_values[r] = 0.00000000;
1256
}// end loop over 'r'
1257
1258
// Loop dofs and call evaluate_basis_derivatives.
1259
for (unsigned int r = 0; r < 4; r++)
1260
{
1261
evaluate_basis_derivatives(r, n, dof_values, coordinates, c);
1262
for (unsigned int s = 0; s < 4*num_derivatives; s++)
1263
{
1264
values[r*4*num_derivatives + s] = dof_values[s];
1265
}// end loop over 's'
1266
}// end loop over 'r'
1267
1268
// Delete pointer.
1269
delete [] dof_values;
1270
}
1271
1272
/// Evaluate linear functional for dof i on the function f
1273
virtual double evaluate_dof(unsigned int i,
1274
const ufc::function& f,
1275
const ufc::cell& c) const
1276
{
1277
// Declare variables for result of evaluation
1278
double vals[4];
1279
1280
// Declare variable for physical coordinates
1281
double y[2];
1282
1283
const double * const * x = c.coordinates;
1284
switch (i)
1285
{
1286
case 0:
1287
{
1288
y[0] = 0.333333333333*x[0][0] + 0.333333333333*x[1][0] + 0.333333333333*x[2][0];
1289
y[1] = 0.333333333333*x[0][1] + 0.333333333333*x[1][1] + 0.333333333333*x[2][1];
1290
f.evaluate(vals, y, c);
1291
return vals[0];
1292
break;
1293
}
1294
case 1:
1295
{
1296
y[0] = 0.333333333333*x[0][0] + 0.333333333333*x[1][0] + 0.333333333333*x[2][0];
1297
y[1] = 0.333333333333*x[0][1] + 0.333333333333*x[1][1] + 0.333333333333*x[2][1];
1298
f.evaluate(vals, y, c);
1299
return vals[1];
1300
break;
1301
}
1302
case 2:
1303
{
1304
y[0] = 0.333333333333*x[0][0] + 0.333333333333*x[1][0] + 0.333333333333*x[2][0];
1305
y[1] = 0.333333333333*x[0][1] + 0.333333333333*x[1][1] + 0.333333333333*x[2][1];
1306
f.evaluate(vals, y, c);
1307
return vals[2];
1308
break;
1309
}
1310
case 3:
1311
{
1312
y[0] = 0.333333333333*x[0][0] + 0.333333333333*x[1][0] + 0.333333333333*x[2][0];
1313
y[1] = 0.333333333333*x[0][1] + 0.333333333333*x[1][1] + 0.333333333333*x[2][1];
1314
f.evaluate(vals, y, c);
1315
return vals[3];
1316
break;
1317
}
1318
}
1319
1320
return 0.0;
1321
}
1322
1323
/// Evaluate linear functionals for all dofs on the function f
1324
virtual void evaluate_dofs(double* values,
1325
const ufc::function& f,
1326
const ufc::cell& c) const
1327
{
1328
// Declare variables for result of evaluation
1329
double vals[4];
1330
1331
// Declare variable for physical coordinates
1332
double y[2];
1333
1334
const double * const * x = c.coordinates;
1335
y[0] = 0.333333333333*x[0][0] + 0.333333333333*x[1][0] + 0.333333333333*x[2][0];
1336
y[1] = 0.333333333333*x[0][1] + 0.333333333333*x[1][1] + 0.333333333333*x[2][1];
1337
f.evaluate(vals, y, c);
1338
values[0] = vals[0];
1339
y[0] = 0.333333333333*x[0][0] + 0.333333333333*x[1][0] + 0.333333333333*x[2][0];
1340
y[1] = 0.333333333333*x[0][1] + 0.333333333333*x[1][1] + 0.333333333333*x[2][1];
1341
f.evaluate(vals, y, c);
1342
values[1] = vals[1];
1343
y[0] = 0.333333333333*x[0][0] + 0.333333333333*x[1][0] + 0.333333333333*x[2][0];
1344
y[1] = 0.333333333333*x[0][1] + 0.333333333333*x[1][1] + 0.333333333333*x[2][1];
1345
f.evaluate(vals, y, c);
1346
values[2] = vals[2];
1347
y[0] = 0.333333333333*x[0][0] + 0.333333333333*x[1][0] + 0.333333333333*x[2][0];
1348
y[1] = 0.333333333333*x[0][1] + 0.333333333333*x[1][1] + 0.333333333333*x[2][1];
1349
f.evaluate(vals, y, c);
1350
values[3] = vals[3];
1351
}
1352
1353
/// Interpolate vertex values from dof values
1354
virtual void interpolate_vertex_values(double* vertex_values,
1355
const double* dof_values,
1356
const ufc::cell& c) const
1357
{
1358
// Evaluate function and change variables
1359
vertex_values[0] = dof_values[0];
1360
vertex_values[4] = dof_values[0];
1361
vertex_values[8] = dof_values[0];
1362
// Evaluate function and change variables
1363
vertex_values[1] = dof_values[1];
1364
vertex_values[5] = dof_values[1];
1365
vertex_values[9] = dof_values[1];
1366
// Evaluate function and change variables
1367
vertex_values[2] = dof_values[2];
1368
vertex_values[6] = dof_values[2];
1369
vertex_values[10] = dof_values[2];
1370
// Evaluate function and change variables
1371
vertex_values[3] = dof_values[3];
1372
vertex_values[7] = dof_values[3];
1373
vertex_values[11] = dof_values[3];
1374
}
1375
1376
/// Return the number of sub elements (for a mixed element)
1377
virtual unsigned int num_sub_elements() const
1378
{
1379
return 4;
1380
}
1381
1382
/// Create a new finite element for sub element i (for a mixed element)
1383
virtual ufc::finite_element* create_sub_element(unsigned int i) const
1384
{
1385
switch (i)
1386
{
1387
case 0:
1388
{
1389
return new tensorweightedpoisson_finite_element_0();
1390
break;
1391
}
1392
case 1:
1393
{
1394
return new tensorweightedpoisson_finite_element_0();
1395
break;
1396
}
1397
case 2:
1398
{
1399
return new tensorweightedpoisson_finite_element_0();
1400
break;
1401
}
1402
case 3:
1403
{
1404
return new tensorweightedpoisson_finite_element_0();
1405
break;
1406
}
1407
}
1408
1409
return 0;
1410
}
1411
1412
};
1413
1414
/// This class defines the interface for a finite element.
1415
1416
class tensorweightedpoisson_finite_element_2: public ufc::finite_element
1417
{
1418
public:
1419
1420
/// Constructor
1421
tensorweightedpoisson_finite_element_2() : ufc::finite_element()
1422
{
1423
// Do nothing
1424
}
1425
1426
/// Destructor
1427
virtual ~tensorweightedpoisson_finite_element_2()
1428
{
1429
// Do nothing
1430
}
1431
1432
/// Return a string identifying the finite element
1433
virtual const char* signature() const
1434
{
1435
return "FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1)";
1436
}
1437
1438
/// Return the cell shape
1439
virtual ufc::shape cell_shape() const
1440
{
1441
return ufc::triangle;
1442
}
1443
1444
/// Return the dimension of the finite element function space
1445
virtual unsigned int space_dimension() const
1446
{
1447
return 3;
1448
}
1449
1450
/// Return the rank of the value space
1451
virtual unsigned int value_rank() const
1452
{
1453
return 0;
1454
}
1455
1456
/// Return the dimension of the value space for axis i
1457
virtual unsigned int value_dimension(unsigned int i) const
1458
{
1459
return 1;
1460
}
1461
1462
/// Evaluate basis function i at given point in cell
1463
virtual void evaluate_basis(unsigned int i,
1464
double* values,
1465
const double* coordinates,
1466
const ufc::cell& c) const
1467
{
1468
// Extract vertex coordinates
1469
const double * const * x = c.coordinates;
1470
1471
// Compute Jacobian of affine map from reference cell
1472
const double J_00 = x[1][0] - x[0][0];
1473
const double J_01 = x[2][0] - x[0][0];
1474
const double J_10 = x[1][1] - x[0][1];
1475
const double J_11 = x[2][1] - x[0][1];
1476
1477
// Compute determinant of Jacobian
1478
double detJ = J_00*J_11 - J_01*J_10;
1479
1480
// Compute inverse of Jacobian
1481
1482
// Compute constants
1483
const double C0 = x[1][0] + x[2][0];
1484
const double C1 = x[1][1] + x[2][1];
1485
1486
// Get coordinates and map to the reference (FIAT) element
1487
double X = (J_01*(C1 - 2.0*coordinates[1]) + J_11*(2.0*coordinates[0] - C0)) / detJ;
1488
double Y = (J_00*(2.0*coordinates[1] - C1) + J_10*(C0 - 2.0*coordinates[0])) / detJ;
1489
1490
// Reset values
1491
*values = 0.00000000;
1492
1493
// Map degree of freedom to element degree of freedom
1494
const unsigned int dof = i;
1495
1496
// Array of basisvalues
1497
double basisvalues[3] = {0.00000000, 0.00000000, 0.00000000};
1498
1499
// Declare helper variables
1500
unsigned int rr = 0;
1501
unsigned int ss = 0;
1502
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
1503
1504
// Compute basisvalues
1505
basisvalues[0] = 1.00000000;
1506
basisvalues[1] = tmp0;
1507
for (unsigned int r = 0; r < 1; r++)
1508
{
1509
rr = (r + 1)*(r + 1 + 1)/2 + 1;
1510
ss = r*(r + 1)/2;
1511
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
1512
}// end loop over 'r'
1513
for (unsigned int r = 0; r < 2; r++)
1514
{
1515
for (unsigned int s = 0; s < 2 - r; s++)
1516
{
1517
rr = (r + s)*(r + s + 1)/2 + s;
1518
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
1519
}// end loop over 's'
1520
}// end loop over 'r'
1521
1522
// Table(s) of coefficients
1523
static const double coefficients0[3][3] = \
1524
{{0.47140452, -0.28867513, -0.16666667},
1525
{0.47140452, 0.28867513, -0.16666667},
1526
{0.47140452, 0.00000000, 0.33333333}};
1527
1528
// Compute value(s).
1529
for (unsigned int r = 0; r < 3; r++)
1530
{
1531
*values += coefficients0[dof][r]*basisvalues[r];
1532
}// end loop over 'r'
1533
}
1534
1535
/// Evaluate all basis functions at given point in cell
1536
virtual void evaluate_basis_all(double* values,
1537
const double* coordinates,
1538
const ufc::cell& c) const
1539
{
1540
// Helper variable to hold values of a single dof.
1541
double dof_values = 0.00000000;
1542
1543
// Loop dofs and call evaluate_basis.
1544
for (unsigned int r = 0; r < 3; r++)
1545
{
1546
evaluate_basis(r, &dof_values, coordinates, c);
1547
values[r] = dof_values;
1548
}// end loop over 'r'
1549
}
1550
1551
/// Evaluate order n derivatives of basis function i at given point in cell
1552
virtual void evaluate_basis_derivatives(unsigned int i,
1553
unsigned int n,
1554
double* values,
1555
const double* coordinates,
1556
const ufc::cell& c) const
1557
{
1558
// Extract vertex coordinates
1559
const double * const * x = c.coordinates;
1560
1561
// Compute Jacobian of affine map from reference cell
1562
const double J_00 = x[1][0] - x[0][0];
1563
const double J_01 = x[2][0] - x[0][0];
1564
const double J_10 = x[1][1] - x[0][1];
1565
const double J_11 = x[2][1] - x[0][1];
1566
1567
// Compute determinant of Jacobian
1568
double detJ = J_00*J_11 - J_01*J_10;
1569
1570
// Compute inverse of Jacobian
1571
const double K_00 = J_11 / detJ;
1572
const double K_01 = -J_01 / detJ;
1573
const double K_10 = -J_10 / detJ;
1574
const double K_11 = J_00 / detJ;
1575
1576
// Compute constants
1577
const double C0 = x[1][0] + x[2][0];
1578
const double C1 = x[1][1] + x[2][1];
1579
1580
// Get coordinates and map to the reference (FIAT) element
1581
double X = (J_01*(C1 - 2.0*coordinates[1]) + J_11*(2.0*coordinates[0] - C0)) / detJ;
1582
double Y = (J_00*(2.0*coordinates[1] - C1) + J_10*(C0 - 2.0*coordinates[0])) / detJ;
1583
1584
// Compute number of derivatives.
1585
unsigned int num_derivatives = 1;
1586
1587
for (unsigned int r = 0; r < n; r++)
1588
{
1589
num_derivatives *= 2;
1590
}// end loop over 'r'
1591
1592
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
1593
unsigned int **combinations = new unsigned int *[num_derivatives];
1594
for (unsigned int row = 0; row < num_derivatives; row++)
1595
{
1596
combinations[row] = new unsigned int [n];
1597
for (unsigned int col = 0; col < n; col++)
1598
combinations[row][col] = 0;
1599
}
1600
1601
// Generate combinations of derivatives
1602
for (unsigned int row = 1; row < num_derivatives; row++)
1603
{
1604
for (unsigned int num = 0; num < row; num++)
1605
{
1606
for (unsigned int col = n-1; col+1 > 0; col--)
1607
{
1608
if (combinations[row][col] + 1 > 1)
1609
combinations[row][col] = 0;
1610
else
1611
{
1612
combinations[row][col] += 1;
1613
break;
1614
}
1615
}
1616
}
1617
}
1618
1619
// Compute inverse of Jacobian
1620
const double Jinv[2][2] = {{K_00, K_01}, {K_10, K_11}};
1621
1622
// Declare transformation matrix
1623
// Declare pointer to two dimensional array and initialise
1624
double **transform = new double *[num_derivatives];
1625
1626
for (unsigned int j = 0; j < num_derivatives; j++)
1627
{
1628
transform[j] = new double [num_derivatives];
1629
for (unsigned int k = 0; k < num_derivatives; k++)
1630
transform[j][k] = 1;
1631
}
1632
1633
// Construct transformation matrix
1634
for (unsigned int row = 0; row < num_derivatives; row++)
1635
{
1636
for (unsigned int col = 0; col < num_derivatives; col++)
1637
{
1638
for (unsigned int k = 0; k < n; k++)
1639
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
1640
}
1641
}
1642
1643
// Reset values. Assuming that values is always an array.
1644
for (unsigned int r = 0; r < num_derivatives; r++)
1645
{
1646
values[r] = 0.00000000;
1647
}// end loop over 'r'
1648
1649
// Map degree of freedom to element degree of freedom
1650
const unsigned int dof = i;
1651
1652
// Array of basisvalues
1653
double basisvalues[3] = {0.00000000, 0.00000000, 0.00000000};
1654
1655
// Declare helper variables
1656
unsigned int rr = 0;
1657
unsigned int ss = 0;
1658
double tmp0 = (1.00000000 + Y + 2.00000000*X)/2.00000000;
1659
1660
// Compute basisvalues
1661
basisvalues[0] = 1.00000000;
1662
basisvalues[1] = tmp0;
1663
for (unsigned int r = 0; r < 1; r++)
1664
{
1665
rr = (r + 1)*(r + 1 + 1)/2 + 1;
1666
ss = r*(r + 1)/2;
1667
basisvalues[rr] = basisvalues[ss]*(0.50000000 + r + Y*(1.50000000 + r));
1668
}// end loop over 'r'
1669
for (unsigned int r = 0; r < 2; r++)
1670
{
1671
for (unsigned int s = 0; s < 2 - r; s++)
1672
{
1673
rr = (r + s)*(r + s + 1)/2 + s;
1674
basisvalues[rr] *= std::sqrt((0.50000000 + r)*(1.00000000 + r + s));
1675
}// end loop over 's'
1676
}// end loop over 'r'
1677
1678
// Table(s) of coefficients
1679
static const double coefficients0[3][3] = \
1680
{{0.47140452, -0.28867513, -0.16666667},
1681
{0.47140452, 0.28867513, -0.16666667},
1682
{0.47140452, 0.00000000, 0.33333333}};
1683
1684
// Tables of derivatives of the polynomial base (transpose).
1685
static const double dmats0[3][3] = \
1686
{{0.00000000, 0.00000000, 0.00000000},
1687
{4.89897949, 0.00000000, 0.00000000},
1688
{0.00000000, 0.00000000, 0.00000000}};
1689
1690
static const double dmats1[3][3] = \
1691
{{0.00000000, 0.00000000, 0.00000000},
1692
{2.44948974, 0.00000000, 0.00000000},
1693
{4.24264069, 0.00000000, 0.00000000}};
1694
1695
// Compute reference derivatives
1696
// Declare pointer to array of derivatives on FIAT element
1697
double *derivatives = new double [num_derivatives];
1698
for (unsigned int r = 0; r < num_derivatives; r++)
1699
{
1700
derivatives[r] = 0.00000000;
1701
}// end loop over 'r'
1702
1703
// Declare derivative matrix (of polynomial basis).
1704
double dmats[3][3] = \
1705
{{1.00000000, 0.00000000, 0.00000000},
1706
{0.00000000, 1.00000000, 0.00000000},
1707
{0.00000000, 0.00000000, 1.00000000}};
1708
1709
// Declare (auxiliary) derivative matrix (of polynomial basis).
1710
double dmats_old[3][3] = \
1711
{{1.00000000, 0.00000000, 0.00000000},
1712
{0.00000000, 1.00000000, 0.00000000},
1713
{0.00000000, 0.00000000, 1.00000000}};
1714
1715
// Loop possible derivatives.
1716
for (unsigned int r = 0; r < num_derivatives; r++)
1717
{
1718
// Resetting dmats values to compute next derivative.
1719
for (unsigned int t = 0; t < 3; t++)
1720
{
1721
for (unsigned int u = 0; u < 3; u++)
1722
{
1723
dmats[t][u] = 0.00000000;
1724
if (t == u)
1725
{
1726
dmats[t][u] = 1.00000000;
1727
}
1728
1729
}// end loop over 'u'
1730
}// end loop over 't'
1731
1732
// Looping derivative order to generate dmats.
1733
for (unsigned int s = 0; s < n; s++)
1734
{
1735
// Updating dmats_old with new values and resetting dmats.
1736
for (unsigned int t = 0; t < 3; t++)
1737
{
1738
for (unsigned int u = 0; u < 3; u++)
1739
{
1740
dmats_old[t][u] = dmats[t][u];
1741
dmats[t][u] = 0.00000000;
1742
}// end loop over 'u'
1743
}// end loop over 't'
1744
1745
// Update dmats using an inner product.
1746
if (combinations[r][s] == 0)
1747
{
1748
for (unsigned int t = 0; t < 3; t++)
1749
{
1750
for (unsigned int u = 0; u < 3; u++)
1751
{
1752
for (unsigned int tu = 0; tu < 3; tu++)
1753
{
1754
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1755
}// end loop over 'tu'
1756
}// end loop over 'u'
1757
}// end loop over 't'
1758
}
1759
1760
if (combinations[r][s] == 1)
1761
{
1762
for (unsigned int t = 0; t < 3; t++)
1763
{
1764
for (unsigned int u = 0; u < 3; u++)
1765
{
1766
for (unsigned int tu = 0; tu < 3; tu++)
1767
{
1768
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1769
}// end loop over 'tu'
1770
}// end loop over 'u'
1771
}// end loop over 't'
1772
}
1773
1774
}// end loop over 's'
1775
for (unsigned int s = 0; s < 3; s++)
1776
{
1777
for (unsigned int t = 0; t < 3; t++)
1778
{
1779
derivatives[r] += coefficients0[dof][s]*dmats[s][t]*basisvalues[t];
1780
}// end loop over 't'
1781
}// end loop over 's'
1782
}// end loop over 'r'
1783
1784
// Transform derivatives back to physical element
1785
for (unsigned int row = 0; row < num_derivatives; row++)
1786
{
1787
for (unsigned int col = 0; col < num_derivatives; col++)
1788
{
1789
values[row] += transform[row][col]*derivatives[col];
1790
}
1791
}
1792
1793
// Delete pointer to array of derivatives on FIAT element
1794
delete [] derivatives;
1795
1796
// Delete pointer to array of combinations of derivatives and transform
1797
for (unsigned int r = 0; r < num_derivatives; r++)
1798
{
1799
delete [] combinations[r];
1800
delete [] transform[r];
1801
}// end loop over 'r'
1802
delete [] combinations;
1803
delete [] transform;
1804
}
1805
1806
/// Evaluate order n derivatives of all basis functions at given point in cell
1807
virtual void evaluate_basis_derivatives_all(unsigned int n,
1808
double* values,
1809
const double* coordinates,
1810
const ufc::cell& c) const
1811
{
1812
// Compute number of derivatives.
1813
unsigned int num_derivatives = 1;
1814
1815
for (unsigned int r = 0; r < n; r++)
1816
{
1817
num_derivatives *= 2;
1818
}// end loop over 'r'
1819
1820
// Helper variable to hold values of a single dof.
1821
double *dof_values = new double [num_derivatives];
1822
for (unsigned int r = 0; r < num_derivatives; r++)
1823
{
1824
dof_values[r] = 0.00000000;
1825
}// end loop over 'r'
1826
1827
// Loop dofs and call evaluate_basis_derivatives.
1828
for (unsigned int r = 0; r < 3; r++)
1829
{
1830
evaluate_basis_derivatives(r, n, dof_values, coordinates, c);
1831
for (unsigned int s = 0; s < num_derivatives; s++)
1832
{
1833
values[r*num_derivatives + s] = dof_values[s];
1834
}// end loop over 's'
1835
}// end loop over 'r'
1836
1837
// Delete pointer.
1838
delete [] dof_values;
1839
}
1840
1841
/// Evaluate linear functional for dof i on the function f
1842
virtual double evaluate_dof(unsigned int i,
1843
const ufc::function& f,
1844
const ufc::cell& c) const
1845
{
1846
// Declare variables for result of evaluation
1847
double vals[1];
1848
1849
// Declare variable for physical coordinates
1850
double y[2];
1851
1852
const double * const * x = c.coordinates;
1853
switch (i)
1854
{
1855
case 0:
1856
{
1857
y[0] = x[0][0];
1858
y[1] = x[0][1];
1859
f.evaluate(vals, y, c);
1860
return vals[0];
1861
break;
1862
}
1863
case 1:
1864
{
1865
y[0] = x[1][0];
1866
y[1] = x[1][1];
1867
f.evaluate(vals, y, c);
1868
return vals[0];
1869
break;
1870
}
1871
case 2:
1872
{
1873
y[0] = x[2][0];
1874
y[1] = x[2][1];
1875
f.evaluate(vals, y, c);
1876
return vals[0];
1877
break;
1878
}
1879
}
1880
1881
return 0.0;
1882
}
1883
1884
/// Evaluate linear functionals for all dofs on the function f
1885
virtual void evaluate_dofs(double* values,
1886
const ufc::function& f,
1887
const ufc::cell& c) const
1888
{
1889
// Declare variables for result of evaluation
1890
double vals[1];
1891
1892
// Declare variable for physical coordinates
1893
double y[2];
1894
1895
const double * const * x = c.coordinates;
1896
y[0] = x[0][0];
1897
y[1] = x[0][1];
1898
f.evaluate(vals, y, c);
1899
values[0] = vals[0];
1900
y[0] = x[1][0];
1901
y[1] = x[1][1];
1902
f.evaluate(vals, y, c);
1903
values[1] = vals[0];
1904
y[0] = x[2][0];
1905
y[1] = x[2][1];
1906
f.evaluate(vals, y, c);
1907
values[2] = vals[0];
1908
}
1909
1910
/// Interpolate vertex values from dof values
1911
virtual void interpolate_vertex_values(double* vertex_values,
1912
const double* dof_values,
1913
const ufc::cell& c) const
1914
{
1915
// Evaluate function and change variables
1916
vertex_values[0] = dof_values[0];
1917
vertex_values[1] = dof_values[1];
1918
vertex_values[2] = dof_values[2];
1919
}
1920
1921
/// Return the number of sub elements (for a mixed element)
1922
virtual unsigned int num_sub_elements() const
1923
{
1924
return 0;
1925
}
1926
1927
/// Create a new finite element for sub element i (for a mixed element)
1928
virtual ufc::finite_element* create_sub_element(unsigned int i) const
1929
{
1930
return 0;
1931
}
1932
1933
};
1934
1935
/// This class defines the interface for a local-to-global mapping of
1936
/// degrees of freedom (dofs).
1937
1938
class tensorweightedpoisson_dof_map_0: public ufc::dof_map
1939
{
1940
private:
1941
1942
unsigned int _global_dimension;
1943
1944
public:
1945
1946
/// Constructor
1947
tensorweightedpoisson_dof_map_0() : ufc::dof_map()
1948
{
1949
_global_dimension = 0;
1950
}
1951
1952
/// Destructor
1953
virtual ~tensorweightedpoisson_dof_map_0()
1954
{
1955
// Do nothing
1956
}
1957
1958
/// Return a string identifying the dof map
1959
virtual const char* signature() const
1960
{
1961
return "FFC dofmap for FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 0)";
1962
}
1963
1964
/// Return true iff mesh entities of topological dimension d are needed
1965
virtual bool needs_mesh_entities(unsigned int d) const
1966
{
1967
switch (d)
1968
{
1969
case 0:
1970
{
1971
return false;
1972
break;
1973
}
1974
case 1:
1975
{
1976
return false;
1977
break;
1978
}
1979
case 2:
1980
{
1981
return true;
1982
break;
1983
}
1984
}
1985
1986
return false;
1987
}
1988
1989
/// Initialize dof map for mesh (return true iff init_cell() is needed)
1990
virtual bool init_mesh(const ufc::mesh& m)
1991
{
1992
_global_dimension = m.num_entities[2];
1993
return false;
1994
}
1995
1996
/// Initialize dof map for given cell
1997
virtual void init_cell(const ufc::mesh& m,
1998
const ufc::cell& c)
1999
{
2000
// Do nothing
2001
}
2002
2003
/// Finish initialization of dof map for cells
2004
virtual void init_cell_finalize()
2005
{
2006
// Do nothing
2007
}
2008
2009
/// Return the dimension of the global finite element function space
2010
virtual unsigned int global_dimension() const
2011
{
2012
return _global_dimension;
2013
}
2014
2015
/// Return the dimension of the local finite element function space for a cell
2016
virtual unsigned int local_dimension(const ufc::cell& c) const
2017
{
2018
return 1;
2019
}
2020
2021
/// Return the maximum dimension of the local finite element function space
2022
virtual unsigned int max_local_dimension() const
2023
{
2024
return 1;
2025
}
2026
2027
// Return the geometric dimension of the coordinates this dof map provides
2028
virtual unsigned int geometric_dimension() const
2029
{
2030
return 2;
2031
}
2032
2033
/// Return the number of dofs on each cell facet
2034
virtual unsigned int num_facet_dofs() const
2035
{
2036
return 0;
2037
}
2038
2039
/// Return the number of dofs associated with each cell entity of dimension d
2040
virtual unsigned int num_entity_dofs(unsigned int d) const
2041
{
2042
switch (d)
2043
{
2044
case 0:
2045
{
2046
return 0;
2047
break;
2048
}
2049
case 1:
2050
{
2051
return 0;
2052
break;
2053
}
2054
case 2:
2055
{
2056
return 1;
2057
break;
2058
}
2059
}
2060
2061
return 0;
2062
}
2063
2064
/// Tabulate the local-to-global mapping of dofs on a cell
2065
virtual void tabulate_dofs(unsigned int* dofs,
2066
const ufc::mesh& m,
2067
const ufc::cell& c) const
2068
{
2069
dofs[0] = c.entity_indices[2][0];
2070
}
2071
2072
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
2073
virtual void tabulate_facet_dofs(unsigned int* dofs,
2074
unsigned int facet) const
2075
{
2076
switch (facet)
2077
{
2078
case 0:
2079
{
2080
2081
break;
2082
}
2083
case 1:
2084
{
2085
2086
break;
2087
}
2088
case 2:
2089
{
2090
2091
break;
2092
}
2093
}
2094
2095
}
2096
2097
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
2098
virtual void tabulate_entity_dofs(unsigned int* dofs,
2099
unsigned int d, unsigned int i) const
2100
{
2101
if (d > 2)
2102
{
2103
std::cerr << "*** FFC warning: " << "d is larger than dimension (2)" << std::endl;
2104
}
2105
2106
switch (d)
2107
{
2108
case 0:
2109
{
2110
2111
break;
2112
}
2113
case 1:
2114
{
2115
2116
break;
2117
}
2118
case 2:
2119
{
2120
if (i > 0)
2121
{
2122
std::cerr << "*** FFC warning: " << "i is larger than number of entities (0)" << std::endl;
2123
}
2124
2125
dofs[0] = 0;
2126
break;
2127
}
2128
}
2129
2130
}
2131
2132
/// Tabulate the coordinates of all dofs on a cell
2133
virtual void tabulate_coordinates(double** coordinates,
2134
const ufc::cell& c) const
2135
{
2136
const double * const * x = c.coordinates;
2137
2138
coordinates[0][0] = 0.333333333333*x[0][0] + 0.333333333333*x[1][0] + 0.333333333333*x[2][0];
2139
coordinates[0][1] = 0.333333333333*x[0][1] + 0.333333333333*x[1][1] + 0.333333333333*x[2][1];
2140
}
2141
2142
/// Return the number of sub dof maps (for a mixed element)
2143
virtual unsigned int num_sub_dof_maps() const
2144
{
2145
return 0;
2146
}
2147
2148
/// Create a new dof_map for sub dof map i (for a mixed element)
2149
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
2150
{
2151
return 0;
2152
}
2153
2154
};
2155
2156
/// This class defines the interface for a local-to-global mapping of
2157
/// degrees of freedom (dofs).
2158
2159
class tensorweightedpoisson_dof_map_1: public ufc::dof_map
2160
{
2161
private:
2162
2163
unsigned int _global_dimension;
2164
2165
public:
2166
2167
/// Constructor
2168
tensorweightedpoisson_dof_map_1() : ufc::dof_map()
2169
{
2170
_global_dimension = 0;
2171
}
2172
2173
/// Destructor
2174
virtual ~tensorweightedpoisson_dof_map_1()
2175
{
2176
// Do nothing
2177
}
2178
2179
/// Return a string identifying the dof map
2180
virtual const char* signature() const
2181
{
2182
return "FFC dofmap for TensorElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 0, (2, 2), None)";
2183
}
2184
2185
/// Return true iff mesh entities of topological dimension d are needed
2186
virtual bool needs_mesh_entities(unsigned int d) const
2187
{
2188
switch (d)
2189
{
2190
case 0:
2191
{
2192
return false;
2193
break;
2194
}
2195
case 1:
2196
{
2197
return false;
2198
break;
2199
}
2200
case 2:
2201
{
2202
return true;
2203
break;
2204
}
2205
}
2206
2207
return false;
2208
}
2209
2210
/// Initialize dof map for mesh (return true iff init_cell() is needed)
2211
virtual bool init_mesh(const ufc::mesh& m)
2212
{
2213
_global_dimension = 4*m.num_entities[2];
2214
return false;
2215
}
2216
2217
/// Initialize dof map for given cell
2218
virtual void init_cell(const ufc::mesh& m,
2219
const ufc::cell& c)
2220
{
2221
// Do nothing
2222
}
2223
2224
/// Finish initialization of dof map for cells
2225
virtual void init_cell_finalize()
2226
{
2227
// Do nothing
2228
}
2229
2230
/// Return the dimension of the global finite element function space
2231
virtual unsigned int global_dimension() const
2232
{
2233
return _global_dimension;
2234
}
2235
2236
/// Return the dimension of the local finite element function space for a cell
2237
virtual unsigned int local_dimension(const ufc::cell& c) const
2238
{
2239
return 4;
2240
}
2241
2242
/// Return the maximum dimension of the local finite element function space
2243
virtual unsigned int max_local_dimension() const
2244
{
2245
return 4;
2246
}
2247
2248
// Return the geometric dimension of the coordinates this dof map provides
2249
virtual unsigned int geometric_dimension() const
2250
{
2251
return 2;
2252
}
2253
2254
/// Return the number of dofs on each cell facet
2255
virtual unsigned int num_facet_dofs() const
2256
{
2257
return 0;
2258
}
2259
2260
/// Return the number of dofs associated with each cell entity of dimension d
2261
virtual unsigned int num_entity_dofs(unsigned int d) const
2262
{
2263
switch (d)
2264
{
2265
case 0:
2266
{
2267
return 0;
2268
break;
2269
}
2270
case 1:
2271
{
2272
return 0;
2273
break;
2274
}
2275
case 2:
2276
{
2277
return 4;
2278
break;
2279
}
2280
}
2281
2282
return 0;
2283
}
2284
2285
/// Tabulate the local-to-global mapping of dofs on a cell
2286
virtual void tabulate_dofs(unsigned int* dofs,
2287
const ufc::mesh& m,
2288
const ufc::cell& c) const
2289
{
2290
unsigned int offset = 0;
2291
2292
dofs[0] = offset + c.entity_indices[2][0];
2293
offset += m.num_entities[2];
2294
dofs[1] = offset + c.entity_indices[2][0];
2295
offset += m.num_entities[2];
2296
dofs[2] = offset + c.entity_indices[2][0];
2297
offset += m.num_entities[2];
2298
dofs[3] = offset + c.entity_indices[2][0];
2299
offset += m.num_entities[2];
2300
}
2301
2302
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
2303
virtual void tabulate_facet_dofs(unsigned int* dofs,
2304
unsigned int facet) const
2305
{
2306
switch (facet)
2307
{
2308
case 0:
2309
{
2310
2311
break;
2312
}
2313
case 1:
2314
{
2315
2316
break;
2317
}
2318
case 2:
2319
{
2320
2321
break;
2322
}
2323
}
2324
2325
}
2326
2327
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
2328
virtual void tabulate_entity_dofs(unsigned int* dofs,
2329
unsigned int d, unsigned int i) const
2330
{
2331
if (d > 2)
2332
{
2333
std::cerr << "*** FFC warning: " << "d is larger than dimension (2)" << std::endl;
2334
}
2335
2336
switch (d)
2337
{
2338
case 0:
2339
{
2340
2341
break;
2342
}
2343
case 1:
2344
{
2345
2346
break;
2347
}
2348
case 2:
2349
{
2350
if (i > 0)
2351
{
2352
std::cerr << "*** FFC warning: " << "i is larger than number of entities (0)" << std::endl;
2353
}
2354
2355
dofs[0] = 0;
2356
dofs[1] = 1;
2357
dofs[2] = 2;
2358
dofs[3] = 3;
2359
break;
2360
}
2361
}
2362
2363
}
2364
2365
/// Tabulate the coordinates of all dofs on a cell
2366
virtual void tabulate_coordinates(double** coordinates,
2367
const ufc::cell& c) const
2368
{
2369
const double * const * x = c.coordinates;
2370
2371
coordinates[0][0] = 0.333333333333*x[0][0] + 0.333333333333*x[1][0] + 0.333333333333*x[2][0];
2372
coordinates[0][1] = 0.333333333333*x[0][1] + 0.333333333333*x[1][1] + 0.333333333333*x[2][1];
2373
coordinates[1][0] = 0.333333333333*x[0][0] + 0.333333333333*x[1][0] + 0.333333333333*x[2][0];
2374
coordinates[1][1] = 0.333333333333*x[0][1] + 0.333333333333*x[1][1] + 0.333333333333*x[2][1];
2375
coordinates[2][0] = 0.333333333333*x[0][0] + 0.333333333333*x[1][0] + 0.333333333333*x[2][0];
2376
coordinates[2][1] = 0.333333333333*x[0][1] + 0.333333333333*x[1][1] + 0.333333333333*x[2][1];
2377
coordinates[3][0] = 0.333333333333*x[0][0] + 0.333333333333*x[1][0] + 0.333333333333*x[2][0];
2378
coordinates[3][1] = 0.333333333333*x[0][1] + 0.333333333333*x[1][1] + 0.333333333333*x[2][1];
2379
}
2380
2381
/// Return the number of sub dof maps (for a mixed element)
2382
virtual unsigned int num_sub_dof_maps() const
2383
{
2384
return 4;
2385
}
2386
2387
/// Create a new dof_map for sub dof map i (for a mixed element)
2388
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
2389
{
2390
switch (i)
2391
{
2392
case 0:
2393
{
2394
return new tensorweightedpoisson_dof_map_0();
2395
break;
2396
}
2397
case 1:
2398
{
2399
return new tensorweightedpoisson_dof_map_0();
2400
break;
2401
}
2402
case 2:
2403
{
2404
return new tensorweightedpoisson_dof_map_0();
2405
break;
2406
}
2407
case 3:
2408
{
2409
return new tensorweightedpoisson_dof_map_0();
2410
break;
2411
}
2412
}
2413
2414
return 0;
2415
}
2416
2417
};
2418
2419
/// This class defines the interface for a local-to-global mapping of
2420
/// degrees of freedom (dofs).
2421
2422
class tensorweightedpoisson_dof_map_2: public ufc::dof_map
2423
{
2424
private:
2425
2426
unsigned int _global_dimension;
2427
2428
public:
2429
2430
/// Constructor
2431
tensorweightedpoisson_dof_map_2() : ufc::dof_map()
2432
{
2433
_global_dimension = 0;
2434
}
2435
2436
/// Destructor
2437
virtual ~tensorweightedpoisson_dof_map_2()
2438
{
2439
// Do nothing
2440
}
2441
2442
/// Return a string identifying the dof map
2443
virtual const char* signature() const
2444
{
2445
return "FFC dofmap for FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1)";
2446
}
2447
2448
/// Return true iff mesh entities of topological dimension d are needed
2449
virtual bool needs_mesh_entities(unsigned int d) const
2450
{
2451
switch (d)
2452
{
2453
case 0:
2454
{
2455
return true;
2456
break;
2457
}
2458
case 1:
2459
{
2460
return false;
2461
break;
2462
}
2463
case 2:
2464
{
2465
return false;
2466
break;
2467
}
2468
}
2469
2470
return false;
2471
}
2472
2473
/// Initialize dof map for mesh (return true iff init_cell() is needed)
2474
virtual bool init_mesh(const ufc::mesh& m)
2475
{
2476
_global_dimension = m.num_entities[0];
2477
return false;
2478
}
2479
2480
/// Initialize dof map for given cell
2481
virtual void init_cell(const ufc::mesh& m,
2482
const ufc::cell& c)
2483
{
2484
// Do nothing
2485
}
2486
2487
/// Finish initialization of dof map for cells
2488
virtual void init_cell_finalize()
2489
{
2490
// Do nothing
2491
}
2492
2493
/// Return the dimension of the global finite element function space
2494
virtual unsigned int global_dimension() const
2495
{
2496
return _global_dimension;
2497
}
2498
2499
/// Return the dimension of the local finite element function space for a cell
2500
virtual unsigned int local_dimension(const ufc::cell& c) const
2501
{
2502
return 3;
2503
}
2504
2505
/// Return the maximum dimension of the local finite element function space
2506
virtual unsigned int max_local_dimension() const
2507
{
2508
return 3;
2509
}
2510
2511
// Return the geometric dimension of the coordinates this dof map provides
2512
virtual unsigned int geometric_dimension() const
2513
{
2514
return 2;
2515
}
2516
2517
/// Return the number of dofs on each cell facet
2518
virtual unsigned int num_facet_dofs() const
2519
{
2520
return 2;
2521
}
2522
2523
/// Return the number of dofs associated with each cell entity of dimension d
2524
virtual unsigned int num_entity_dofs(unsigned int d) const
2525
{
2526
switch (d)
2527
{
2528
case 0:
2529
{
2530
return 1;
2531
break;
2532
}
2533
case 1:
2534
{
2535
return 0;
2536
break;
2537
}
2538
case 2:
2539
{
2540
return 0;
2541
break;
2542
}
2543
}
2544
2545
return 0;
2546
}
2547
2548
/// Tabulate the local-to-global mapping of dofs on a cell
2549
virtual void tabulate_dofs(unsigned int* dofs,
2550
const ufc::mesh& m,
2551
const ufc::cell& c) const
2552
{
2553
dofs[0] = c.entity_indices[0][0];
2554
dofs[1] = c.entity_indices[0][1];
2555
dofs[2] = c.entity_indices[0][2];
2556
}
2557
2558
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
2559
virtual void tabulate_facet_dofs(unsigned int* dofs,
2560
unsigned int facet) const
2561
{
2562
switch (facet)
2563
{
2564
case 0:
2565
{
2566
dofs[0] = 1;
2567
dofs[1] = 2;
2568
break;
2569
}
2570
case 1:
2571
{
2572
dofs[0] = 0;
2573
dofs[1] = 2;
2574
break;
2575
}
2576
case 2:
2577
{
2578
dofs[0] = 0;
2579
dofs[1] = 1;
2580
break;
2581
}
2582
}
2583
2584
}
2585
2586
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
2587
virtual void tabulate_entity_dofs(unsigned int* dofs,
2588
unsigned int d, unsigned int i) const
2589
{
2590
if (d > 2)
2591
{
2592
std::cerr << "*** FFC warning: " << "d is larger than dimension (2)" << std::endl;
2593
}
2594
2595
switch (d)
2596
{
2597
case 0:
2598
{
2599
if (i > 2)
2600
{
2601
std::cerr << "*** FFC warning: " << "i is larger than number of entities (2)" << std::endl;
2602
}
2603
2604
switch (i)
2605
{
2606
case 0:
2607
{
2608
dofs[0] = 0;
2609
break;
2610
}
2611
case 1:
2612
{
2613
dofs[0] = 1;
2614
break;
2615
}
2616
case 2:
2617
{
2618
dofs[0] = 2;
2619
break;
2620
}
2621
}
2622
2623
break;
2624
}
2625
case 1:
2626
{
2627
2628
break;
2629
}
2630
case 2:
2631
{
2632
2633
break;
2634
}
2635
}
2636
2637
}
2638
2639
/// Tabulate the coordinates of all dofs on a cell
2640
virtual void tabulate_coordinates(double** coordinates,
2641
const ufc::cell& c) const
2642
{
2643
const double * const * x = c.coordinates;
2644
2645
coordinates[0][0] = x[0][0];
2646
coordinates[0][1] = x[0][1];
2647
coordinates[1][0] = x[1][0];
2648
coordinates[1][1] = x[1][1];
2649
coordinates[2][0] = x[2][0];
2650
coordinates[2][1] = x[2][1];
2651
}
2652
2653
/// Return the number of sub dof maps (for a mixed element)
2654
virtual unsigned int num_sub_dof_maps() const
2655
{
2656
return 0;
2657
}
2658
2659
/// Create a new dof_map for sub dof map i (for a mixed element)
2660
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
2661
{
2662
return 0;
2663
}
2664
2665
};
2666
2667
/// This class defines the interface for the tabulation of the cell
2668
/// tensor corresponding to the local contribution to a form from
2669
/// the integral over a cell.
2670
2671
class tensorweightedpoisson_cell_integral_0_0: public ufc::cell_integral
2672
{
2673
public:
2674
2675
/// Constructor
2676
tensorweightedpoisson_cell_integral_0_0() : ufc::cell_integral()
2677
{
2678
// Do nothing
2679
}
2680
2681
/// Destructor
2682
virtual ~tensorweightedpoisson_cell_integral_0_0()
2683
{
2684
// Do nothing
2685
}
2686
2687
/// Tabulate the tensor for the contribution from a local cell
2688
virtual void tabulate_tensor(double* A,
2689
const double * const * w,
2690
const ufc::cell& c) const
2691
{
2692
// Extract vertex coordinates
2693
const double * const * x = c.coordinates;
2694
2695
// Compute Jacobian of affine map from reference cell
2696
const double J_00 = x[1][0] - x[0][0];
2697
const double J_01 = x[2][0] - x[0][0];
2698
const double J_10 = x[1][1] - x[0][1];
2699
const double J_11 = x[2][1] - x[0][1];
2700
2701
// Compute determinant of Jacobian
2702
double detJ = J_00*J_11 - J_01*J_10;
2703
2704
// Compute inverse of Jacobian
2705
const double K_00 = J_11 / detJ;
2706
const double K_01 = -J_01 / detJ;
2707
const double K_10 = -J_10 / detJ;
2708
const double K_11 = J_00 / detJ;
2709
2710
// Set scale factor
2711
const double det = std::abs(detJ);
2712
2713
// Array of quadrature weights
2714
static const double W1 = 0.50000000;
2715
// Quadrature points on the UFC reference element: (0.33333333, 0.33333333)
2716
2717
// Value of basis functions at quadrature points.
2718
static const double FE0_D01[1][3] = \
2719
{{-1.00000000, 0.00000000, 1.00000000}};
2720
2721
static const double FE0_D10[1][3] = \
2722
{{-1.00000000, 1.00000000, 0.00000000}};
2723
2724
static const double FE1_C0[1][4] = \
2725
{{1.00000000, 0.00000000, 0.00000000, 0.00000000}};
2726
2727
static const double FE1_C1[1][4] = \
2728
{{0.00000000, 1.00000000, 0.00000000, 0.00000000}};
2729
2730
static const double FE1_C2[1][4] = \
2731
{{0.00000000, 0.00000000, 1.00000000, 0.00000000}};
2732
2733
static const double FE1_C3[1][4] = \
2734
{{0.00000000, 0.00000000, 0.00000000, 1.00000000}};
2735
2736
for (unsigned int r = 0; r < 9; r++)
2737
{
2738
A[r] = 0.00000000;
2739
}// end loop over 'r'
2740
2741
// Compute element tensor using UFL quadrature representation
2742
// Optimisations: ('simplify expressions', False), ('ignore zero tables', False), ('non zero columns', False), ('remove zero terms', False), ('ignore ones', False)
2743
2744
// Loop quadrature points for integral
2745
// Number of operations to compute element tensor for following IP loop = 302
2746
// Only 1 integration point, omitting IP loop.
2747
2748
// Coefficient declarations
2749
double F0 = 0.00000000;
2750
double F1 = 0.00000000;
2751
double F2 = 0.00000000;
2752
double F3 = 0.00000000;
2753
2754
// Total number of operations to compute function values = 32
2755
for (unsigned int r = 0; r < 4; r++)
2756
{
2757
F0 += FE1_C0[0][r]*w[0][r];
2758
F1 += FE1_C1[0][r]*w[0][r];
2759
F2 += FE1_C2[0][r]*w[0][r];
2760
F3 += FE1_C3[0][r]*w[0][r];
2761
}// end loop over 'r'
2762
2763
// Number of operations for primary indices: 270
2764
for (unsigned int j = 0; j < 3; j++)
2765
{
2766
for (unsigned int k = 0; k < 3; k++)
2767
{
2768
// Number of operations to compute entry: 30
2769
A[j*3 + k] += ((((((K_00*FE0_D10[0][k] + K_10*FE0_D01[0][k]))*F0 + ((K_01*FE0_D10[0][k] + K_11*FE0_D01[0][k]))*F1))*((K_00*FE0_D10[0][j] + K_10*FE0_D01[0][j])) + ((((K_01*FE0_D10[0][k] + K_11*FE0_D01[0][k]))*F3 + ((K_00*FE0_D10[0][k] + K_10*FE0_D01[0][k]))*F2))*((K_01*FE0_D10[0][j] + K_11*FE0_D01[0][j]))))*W1*det;
2770
}// end loop over 'k'
2771
}// end loop over 'j'
2772
}
2773
2774
};
2775
2776
/// This class defines the interface for the assembly of the global
2777
/// tensor corresponding to a form with r + n arguments, that is, a
2778
/// mapping
2779
///
2780
/// a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
2781
///
2782
/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
2783
/// global tensor A is defined by
2784
///
2785
/// A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
2786
///
2787
/// where each argument Vj represents the application to the
2788
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
2789
/// fixed functions (coefficients).
2790
2791
class tensorweightedpoisson_form_0: public ufc::form
2792
{
2793
public:
2794
2795
/// Constructor
2796
tensorweightedpoisson_form_0() : ufc::form()
2797
{
2798
// Do nothing
2799
}
2800
2801
/// Destructor
2802
virtual ~tensorweightedpoisson_form_0()
2803
{
2804
// Do nothing
2805
}
2806
2807
/// Return a string identifying the form
2808
virtual const char* signature() const
2809
{
2810
return "Form([Integral(IndexSum(Product(Indexed(ComponentTensor(IndexSum(Product(Indexed(Coefficient(TensorElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 0, (2, 2), None), 0), MultiIndex((Index(0), Index(1)), {Index(0): 2, Index(1): 2})), Indexed(ComponentTensor(SpatialDerivative(Argument(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), 1), MultiIndex((Index(2),), {Index(2): 2})), MultiIndex((Index(2),), {Index(2): 2})), MultiIndex((Index(1),), {Index(1): 2}))), MultiIndex((Index(1),), {Index(1): 2})), MultiIndex((Index(0),), {Index(0): 2})), MultiIndex((Index(3),), {Index(3): 2})), Indexed(ComponentTensor(SpatialDerivative(Argument(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), 0), MultiIndex((Index(4),), {Index(4): 2})), MultiIndex((Index(4),), {Index(4): 2})), MultiIndex((Index(3),), {Index(3): 2}))), MultiIndex((Index(3),), {Index(3): 2})), Measure('cell', 0, None))])";
2811
}
2812
2813
/// Return the rank of the global tensor (r)
2814
virtual unsigned int rank() const
2815
{
2816
return 2;
2817
}
2818
2819
/// Return the number of coefficients (n)
2820
virtual unsigned int num_coefficients() const
2821
{
2822
return 1;
2823
}
2824
2825
/// Return the number of cell integrals
2826
virtual unsigned int num_cell_integrals() const
2827
{
2828
return 1;
2829
}
2830
2831
/// Return the number of exterior facet integrals
2832
virtual unsigned int num_exterior_facet_integrals() const
2833
{
2834
return 0;
2835
}
2836
2837
/// Return the number of interior facet integrals
2838
virtual unsigned int num_interior_facet_integrals() const
2839
{
2840
return 0;
2841
}
2842
2843
/// Create a new finite element for argument function i
2844
virtual ufc::finite_element* create_finite_element(unsigned int i) const
2845
{
2846
switch (i)
2847
{
2848
case 0:
2849
{
2850
return new tensorweightedpoisson_finite_element_2();
2851
break;
2852
}
2853
case 1:
2854
{
2855
return new tensorweightedpoisson_finite_element_2();
2856
break;
2857
}
2858
case 2:
2859
{
2860
return new tensorweightedpoisson_finite_element_1();
2861
break;
2862
}
2863
}
2864
2865
return 0;
2866
}
2867
2868
/// Create a new dof map for argument function i
2869
virtual ufc::dof_map* create_dof_map(unsigned int i) const
2870
{
2871
switch (i)
2872
{
2873
case 0:
2874
{
2875
return new tensorweightedpoisson_dof_map_2();
2876
break;
2877
}
2878
case 1:
2879
{
2880
return new tensorweightedpoisson_dof_map_2();
2881
break;
2882
}
2883
case 2:
2884
{
2885
return new tensorweightedpoisson_dof_map_1();
2886
break;
2887
}
2888
}
2889
2890
return 0;
2891
}
2892
2893
/// Create a new cell integral on sub domain i
2894
virtual ufc::cell_integral* create_cell_integral(unsigned int i) const
2895
{
2896
switch (i)
2897
{
2898
case 0:
2899
{
2900
return new tensorweightedpoisson_cell_integral_0_0();
2901
break;
2902
}
2903
}
2904
2905
return 0;
2906
}
2907
2908
/// Create a new exterior facet integral on sub domain i
2909
virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const
2910
{
2911
return 0;
2912
}
2913
2914
/// Create a new interior facet integral on sub domain i
2915
virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const
2916
{
2917
return 0;
2918
}
2919
2920
};
2921
2922
#endif
Loggerhead is a web-based interface for Breezy
Version: 2.0.1