~ubuntu-branches/ubuntu/vivid/ffc/vivid
» Revision
18
: test/regression/references/r_quadrature_O/ProjectionManifold.h
« back to all changes in this revision
Viewing changes to test/regression/references/r_quadrature_O/ProjectionManifold.h
- browse files at revision 18
- compare with another revision
- download diff
- download tarball
- view history from revision 18
- Committer: Package Import Robot
- Author(s): Johannes Ring
- Date: 2014-01-10 13:56:45 UTC
- mfrom: (1.1.14)
- Revision ID: package-import@ubuntu.com-20140110135645-4ozcd71y1oggj44z
Tags: 1.3.0-1
* New upstream release.
* debian/watch: Update URL for move to Bitbucket.
* debian/docs: README -> README.rst and remove TODO.
* debian/control:
- Add python-numpy to Build-Depends.
- Replace python-all with python-all-dev in Build-Depends.
- Add ${shlibs:Depends} to Depends.
- Change to Architecture: any.
- Bump Standards-Version to 3.9.5 (no changes needed).
* debian/rules: Call dh_numpy in override_dh_python2.
* debian/watch: Update URL for move to Bitbucket.
* debian/docs: README -> README.rst and remove TODO.
* debian/control:
- Add python-numpy to Build-Depends.
- Replace python-all with python-all-dev in Build-Depends.
- Add ${shlibs:Depends} to Depends.
- Change to Architecture: any.
- Bump Standards-Version to 3.9.5 (no changes needed).
* debian/rules: Call dh_numpy in override_dh_python2.
- files added:
- .bzrignore
- ffc/dolfin
- ffc/ext
- test/regression/scripts
- test/unit/elements
- files removed:
- README
- doc/sphinx/source/_static
- doc/sphinx/source/_templates
- sandbox
- sandbox/features
- sandbox/jit
- sandbox/multidomain
- sandbox/swig
- test/regression/references
- test/regression/references/output
- test/regression/references/r_auto
- test/regression/references/r_quadrature
- test/regression/references/r_quadrature_O
- files modified:
- ChangeLog
added
removed
test/regression/references/r_quadrature_O/ProjectionManifold.h
1
// This code conforms with the UFC specification version 2.2.0
2
// and was automatically generated by FFC version 1.2.0.
3
//
4
// This code was generated with the following parameters:
5
//
6
// cache_dir: ''
7
// convert_exceptions_to_warnings: True
8
// cpp_optimize: False
9
// cpp_optimize_flags: '-O2'
10
// epsilon: 1e-14
11
// error_control: False
12
// form_postfix: True
13
// format: 'ufc'
14
// log_level: 20
15
// log_prefix: ''
16
// optimize: True
17
// output_dir: '.'
18
// precision: '8'
19
// quadrature_degree: 'auto'
20
// quadrature_rule: 'auto'
21
// representation: 'quadrature'
22
// split: False
23
// swig_binary: 'swig'
24
// swig_path: ''
25
26
#ifndef __PROJECTIONMANIFOLD_H
27
#define __PROJECTIONMANIFOLD_H
28
29
#include <cmath>
30
#include <stdexcept>
31
#include <fstream>
32
#include <ufc.h>
33
34
/// This class defines the interface for a finite element.
35
36
class projectionmanifold_finite_element_0: public ufc::finite_element
37
{
38
public:
39
40
/// Constructor
41
projectionmanifold_finite_element_0() : ufc::finite_element()
42
{
43
// Do nothing
44
}
45
46
/// Destructor
47
virtual ~projectionmanifold_finite_element_0()
48
{
49
// Do nothing
50
}
51
52
/// Return a string identifying the finite element
53
virtual const char* signature() const
54
{
55
return "FiniteElement('Raviart-Thomas', Domain(Cell('triangle', 3), 'triangle_multiverse', 3, 2), 1, None)";
56
}
57
58
/// Return the cell shape
59
virtual ufc::shape cell_shape() const
60
{
61
return ufc::triangle;
62
}
63
64
/// Return the topological dimension of the cell shape
65
virtual std::size_t topological_dimension() const
66
{
67
return 2;
68
}
69
70
/// Return the geometric dimension of the cell shape
71
virtual std::size_t geometric_dimension() const
72
{
73
return 3;
74
}
75
76
/// Return the dimension of the finite element function space
77
virtual std::size_t space_dimension() const
78
{
79
return 3;
80
}
81
82
/// Return the rank of the value space
83
virtual std::size_t value_rank() const
84
{
85
return 1;
86
}
87
88
/// Return the dimension of the value space for axis i
89
virtual std::size_t value_dimension(std::size_t i) const
90
{
91
switch (i)
92
{
93
case 0:
94
{
95
return 3;
96
break;
97
}
98
}
99
100
return 0;
101
}
102
103
/// Evaluate basis function i at given point x in cell
104
virtual void evaluate_basis(std::size_t i,
105
double* values,
106
const double* x,
107
const double* vertex_coordinates,
108
int cell_orientation) const
109
{
110
// Compute Jacobian
111
double J[6];
112
compute_jacobian_triangle_3d(J, vertex_coordinates);
113
114
// Compute Jacobian inverse and determinant
115
double K[6];
116
double detJ;
117
compute_jacobian_inverse_triangle_3d(K, detJ, J);
118
119
120
// Check orientation
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)
125
detJ *= -1;
126
127
128
const double b0 = vertex_coordinates[0];
129
const double b1 = vertex_coordinates[1];
130
const double b2 = vertex_coordinates[2];
131
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;
135
136
137
// Reset values
138
values[0] = 0.0;
139
values[1] = 0.0;
140
switch (i)
141
{
142
case 0:
143
{
144
145
// Array of basisvalues
146
double basisvalues[3] = {0.0, 0.0, 0.0};
147
148
// Declare helper variables
149
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
150
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);
158
159
// Table(s) of coefficients
160
static const double coefficients0[3] = \
161
{0.47140452, 0.28867513, -0.16666667};
162
163
static const double coefficients1[3] = \
164
{0.47140452, 0.0, 0.33333333};
165
166
// Compute value(s)
167
for (unsigned int r = 0; r < 3; r++)
168
{
169
values[0] += coefficients0[r]*basisvalues[r];
170
values[1] += coefficients1[r]*basisvalues[r];
171
}// end loop over 'r'
172
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);
179
break;
180
}
181
case 1:
182
{
183
184
// Array of basisvalues
185
double basisvalues[3] = {0.0, 0.0, 0.0};
186
187
// Declare helper variables
188
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
189
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);
197
198
// Table(s) of coefficients
199
static const double coefficients0[3] = \
200
{0.94280904, -0.28867513, 0.16666667};
201
202
static const double coefficients1[3] = \
203
{-0.47140452, 0.0, -0.33333333};
204
205
// Compute value(s)
206
for (unsigned int r = 0; r < 3; r++)
207
{
208
values[0] += coefficients0[r]*basisvalues[r];
209
values[1] += coefficients1[r]*basisvalues[r];
210
}// end loop over 'r'
211
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);
218
break;
219
}
220
case 2:
221
{
222
223
// Array of basisvalues
224
double basisvalues[3] = {0.0, 0.0, 0.0};
225
226
// Declare helper variables
227
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
228
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);
236
237
// Table(s) of coefficients
238
static const double coefficients0[3] = \
239
{0.47140452, 0.28867513, -0.16666667};
240
241
static const double coefficients1[3] = \
242
{-0.94280904, 0.0, 0.33333333};
243
244
// Compute value(s)
245
for (unsigned int r = 0; r < 3; r++)
246
{
247
values[0] += coefficients0[r]*basisvalues[r];
248
values[1] += coefficients1[r]*basisvalues[r];
249
}// end loop over 'r'
250
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);
257
break;
258
}
259
}
260
261
}
262
263
/// Evaluate all basis functions at given point x in cell
264
virtual void evaluate_basis_all(double* values,
265
const double* x,
266
const double* vertex_coordinates,
267
int cell_orientation) const
268
{
269
// Helper variable to hold values of a single dof.
270
double dof_values[3] = {0.0, 0.0, 0.0};
271
272
// Loop dofs and call evaluate_basis
273
for (unsigned int r = 0; r < 3; r++)
274
{
275
evaluate_basis(r, dof_values, x, vertex_coordinates, cell_orientation);
276
for (unsigned int s = 0; s < 3; s++)
277
{
278
values[r*3 + s] = dof_values[s];
279
}// end loop over 's'
280
}// end loop over 'r'
281
}
282
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,
285
std::size_t n,
286
double* values,
287
const double* x,
288
const double* vertex_coordinates,
289
int cell_orientation) const
290
{
291
// Compute Jacobian
292
double J[6];
293
compute_jacobian_triangle_3d(J, vertex_coordinates);
294
295
// Compute Jacobian inverse and determinant
296
double K[6];
297
double detJ;
298
compute_jacobian_inverse_triangle_3d(K, detJ, J);
299
300
301
// Check orientation
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)
306
detJ *= -1;
307
308
309
const double b0 = vertex_coordinates[0];
310
const double b1 = vertex_coordinates[1];
311
const double b2 = vertex_coordinates[2];
312
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;
316
317
318
// Compute number of derivatives.
319
unsigned int num_derivatives_t = 1;
320
for (unsigned int r = 0; r < n; r++)
321
{
322
num_derivatives_t *= 2;
323
}// end loop over 'r'
324
325
// Compute number of derivatives.
326
unsigned int num_derivatives_g = 1;
327
for (unsigned int r = 0; r < n; r++)
328
{
329
num_derivatives_g *= 3;
330
}// end loop over 'r'
331
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++)
335
{
336
combinations_t[row] = new unsigned int [n];
337
for (unsigned int col = 0; col < n; col++)
338
combinations_t[row][col] = 0;
339
}
340
341
// Generate combinations of derivatives
342
for (unsigned int row = 1; row < num_derivatives_t; row++)
343
{
344
for (unsigned int num = 0; num < row; num++)
345
{
346
for (unsigned int col = n-1; col+1 > 0; col--)
347
{
348
if (combinations_t[row][col] + 1 > 1)
349
combinations_t[row][col] = 0;
350
else
351
{
352
combinations_t[row][col] += 1;
353
break;
354
}
355
}
356
}
357
}
358
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++)
362
{
363
combinations_g[row] = new unsigned int [n];
364
for (unsigned int col = 0; col < n; col++)
365
combinations_g[row][col] = 0;
366
}
367
368
// Generate combinations of derivatives
369
for (unsigned int row = 1; row < num_derivatives_g; row++)
370
{
371
for (unsigned int num = 0; num < row; num++)
372
{
373
for (unsigned int col = n-1; col+1 > 0; col--)
374
{
375
if (combinations_g[row][col] + 1 > 2)
376
combinations_g[row][col] = 0;
377
else
378
{
379
combinations_g[row][col] += 1;
380
break;
381
}
382
}
383
}
384
}
385
386
// Compute inverse of Jacobian
387
const double Jinv[2][3] = {{K[0], K[1], K[2]}, {K[3], K[4], K[5]}};
388
389
// Declare transformation matrix
390
// Declare pointer to two dimensional array and initialise
391
double **transform = new double *[num_derivatives_g];
392
393
for (unsigned int j = 0; j < num_derivatives_g; j++)
394
{
395
transform[j] = new double [num_derivatives_t];
396
for (unsigned int k = 0; k < num_derivatives_t; k++)
397
transform[j][k] = 1;
398
}
399
400
// Construct transformation matrix
401
for (unsigned int row = 0; row < num_derivatives_g; row++)
402
{
403
for (unsigned int col = 0; col < num_derivatives_t; col++)
404
{
405
for (unsigned int k = 0; k < n; k++)
406
transform[row][col] *= Jinv[combinations_t[col][k]][combinations_g[row][k]];
407
}
408
}
409
410
// Reset values. Assuming that values is always an array.
411
for (unsigned int r = 0; r < 3*num_derivatives_g; r++)
412
{
413
values[r] = 0.0;
414
}// end loop over 'r'
415
416
switch (i)
417
{
418
case 0:
419
{
420
421
// Array of basisvalues
422
double basisvalues[3] = {0.0, 0.0, 0.0};
423
424
// Declare helper variables
425
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
426
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);
434
435
// Table(s) of coefficients
436
static const double coefficients0[3] = \
437
{0.47140452, 0.28867513, -0.16666667};
438
439
static const double coefficients1[3] = \
440
{0.47140452, 0.0, 0.33333333};
441
442
// Tables of derivatives of the polynomial base (transpose).
443
static const double dmats0[3][3] = \
444
{{0.0, 0.0, 0.0},
445
{4.8989795, 0.0, 0.0},
446
{0.0, 0.0, 0.0}};
447
448
static const double dmats1[3][3] = \
449
{{0.0, 0.0, 0.0},
450
{2.4494897, 0.0, 0.0},
451
{4.2426407, 0.0, 0.0}};
452
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++)
457
{
458
derivatives[r] = 0.0;
459
}// end loop over 'r'
460
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++)
464
{
465
derivatives_p[r] = 0.0;
466
}// end loop over 'r'
467
468
// Declare derivative matrix (of polynomial basis).
469
double dmats[3][3] = \
470
{{1.0, 0.0, 0.0},
471
{0.0, 1.0, 0.0},
472
{0.0, 0.0, 1.0}};
473
474
// Declare (auxiliary) derivative matrix (of polynomial basis).
475
double dmats_old[3][3] = \
476
{{1.0, 0.0, 0.0},
477
{0.0, 1.0, 0.0},
478
{0.0, 0.0, 1.0}};
479
480
// Loop possible derivatives.
481
for (unsigned int r = 0; r < num_derivatives_t; r++)
482
{
483
// Resetting dmats values to compute next derivative.
484
for (unsigned int t = 0; t < 3; t++)
485
{
486
for (unsigned int u = 0; u < 3; u++)
487
{
488
dmats[t][u] = 0.0;
489
if (t == u)
490
{
491
dmats[t][u] = 1.0;
492
}
493
494
}// end loop over 'u'
495
}// end loop over 't'
496
497
// Looping derivative order to generate dmats.
498
for (unsigned int s = 0; s < n; s++)
499
{
500
// Updating dmats_old with new values and resetting dmats.
501
for (unsigned int t = 0; t < 3; t++)
502
{
503
for (unsigned int u = 0; u < 3; u++)
504
{
505
dmats_old[t][u] = dmats[t][u];
506
dmats[t][u] = 0.0;
507
}// end loop over 'u'
508
}// end loop over 't'
509
510
// Update dmats using an inner product.
511
if (combinations_t[r][s] == 0)
512
{
513
for (unsigned int t = 0; t < 3; t++)
514
{
515
for (unsigned int u = 0; u < 3; u++)
516
{
517
for (unsigned int tu = 0; tu < 3; tu++)
518
{
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'
523
}
524
525
if (combinations_t[r][s] == 1)
526
{
527
for (unsigned int t = 0; t < 3; t++)
528
{
529
for (unsigned int u = 0; u < 3; u++)
530
{
531
for (unsigned int tu = 0; tu < 3; tu++)
532
{
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'
537
}
538
539
}// end loop over 's'
540
for (unsigned int s = 0; s < 3; s++)
541
{
542
for (unsigned int t = 0; t < 3; t++)
543
{
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'
548
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'
556
557
// Transform derivatives back to physical element
558
for (unsigned int r = 0; r < num_derivatives_g; r++)
559
{
560
for (unsigned int s = 0; s < num_derivatives_t; s++)
561
{
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'
567
568
// Delete pointer to array of derivatives on FIAT element
569
delete [] derivatives;
570
571
// Delete pointer to array of reference derivatives on physical element.
572
delete [] derivatives_p;
573
574
// Delete pointer to array of combinations of derivatives and transform
575
for (unsigned int r = 0; r < num_derivatives_t; r++)
576
{
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++)
581
{
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++)
586
{
587
delete [] transform[r];
588
}// end loop over 'r'
589
delete [] transform;
590
break;
591
}
592
case 1:
593
{
594
595
// Array of basisvalues
596
double basisvalues[3] = {0.0, 0.0, 0.0};
597
598
// Declare helper variables
599
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
600
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);
608
609
// Table(s) of coefficients
610
static const double coefficients0[3] = \
611
{0.94280904, -0.28867513, 0.16666667};
612
613
static const double coefficients1[3] = \
614
{-0.47140452, 0.0, -0.33333333};
615
616
// Tables of derivatives of the polynomial base (transpose).
617
static const double dmats0[3][3] = \
618
{{0.0, 0.0, 0.0},
619
{4.8989795, 0.0, 0.0},
620
{0.0, 0.0, 0.0}};
621
622
static const double dmats1[3][3] = \
623
{{0.0, 0.0, 0.0},
624
{2.4494897, 0.0, 0.0},
625
{4.2426407, 0.0, 0.0}};
626
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++)
631
{
632
derivatives[r] = 0.0;
633
}// end loop over 'r'
634
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++)
638
{
639
derivatives_p[r] = 0.0;
640
}// end loop over 'r'
641
642
// Declare derivative matrix (of polynomial basis).
643
double dmats[3][3] = \
644
{{1.0, 0.0, 0.0},
645
{0.0, 1.0, 0.0},
646
{0.0, 0.0, 1.0}};
647
648
// Declare (auxiliary) derivative matrix (of polynomial basis).
649
double dmats_old[3][3] = \
650
{{1.0, 0.0, 0.0},
651
{0.0, 1.0, 0.0},
652
{0.0, 0.0, 1.0}};
653
654
// Loop possible derivatives.
655
for (unsigned int r = 0; r < num_derivatives_t; r++)
656
{
657
// Resetting dmats values to compute next derivative.
658
for (unsigned int t = 0; t < 3; t++)
659
{
660
for (unsigned int u = 0; u < 3; u++)
661
{
662
dmats[t][u] = 0.0;
663
if (t == u)
664
{
665
dmats[t][u] = 1.0;
666
}
667
668
}// end loop over 'u'
669
}// end loop over 't'
670
671
// Looping derivative order to generate dmats.
672
for (unsigned int s = 0; s < n; s++)
673
{
674
// Updating dmats_old with new values and resetting dmats.
675
for (unsigned int t = 0; t < 3; t++)
676
{
677
for (unsigned int u = 0; u < 3; u++)
678
{
679
dmats_old[t][u] = dmats[t][u];
680
dmats[t][u] = 0.0;
681
}// end loop over 'u'
682
}// end loop over 't'
683
684
// Update dmats using an inner product.
685
if (combinations_t[r][s] == 0)
686
{
687
for (unsigned int t = 0; t < 3; t++)
688
{
689
for (unsigned int u = 0; u < 3; u++)
690
{
691
for (unsigned int tu = 0; tu < 3; tu++)
692
{
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'
697
}
698
699
if (combinations_t[r][s] == 1)
700
{
701
for (unsigned int t = 0; t < 3; t++)
702
{
703
for (unsigned int u = 0; u < 3; u++)
704
{
705
for (unsigned int tu = 0; tu < 3; tu++)
706
{
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'
711
}
712
713
}// end loop over 's'
714
for (unsigned int s = 0; s < 3; s++)
715
{
716
for (unsigned int t = 0; t < 3; t++)
717
{
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'
722
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'
730
731
// Transform derivatives back to physical element
732
for (unsigned int r = 0; r < num_derivatives_g; r++)
733
{
734
for (unsigned int s = 0; s < num_derivatives_t; s++)
735
{
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'
741
742
// Delete pointer to array of derivatives on FIAT element
743
delete [] derivatives;
744
745
// Delete pointer to array of reference derivatives on physical element.
746
delete [] derivatives_p;
747
748
// Delete pointer to array of combinations of derivatives and transform
749
for (unsigned int r = 0; r < num_derivatives_t; r++)
750
{
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++)
755
{
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++)
760
{
761
delete [] transform[r];
762
}// end loop over 'r'
763
delete [] transform;
764
break;
765
}
766
case 2:
767
{
768
769
// Array of basisvalues
770
double basisvalues[3] = {0.0, 0.0, 0.0};
771
772
// Declare helper variables
773
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
774
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);
782
783
// Table(s) of coefficients
784
static const double coefficients0[3] = \
785
{0.47140452, 0.28867513, -0.16666667};
786
787
static const double coefficients1[3] = \
788
{-0.94280904, 0.0, 0.33333333};
789
790
// Tables of derivatives of the polynomial base (transpose).
791
static const double dmats0[3][3] = \
792
{{0.0, 0.0, 0.0},
793
{4.8989795, 0.0, 0.0},
794
{0.0, 0.0, 0.0}};
795
796
static const double dmats1[3][3] = \
797
{{0.0, 0.0, 0.0},
798
{2.4494897, 0.0, 0.0},
799
{4.2426407, 0.0, 0.0}};
800
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++)
805
{
806
derivatives[r] = 0.0;
807
}// end loop over 'r'
808
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++)
812
{
813
derivatives_p[r] = 0.0;
814
}// end loop over 'r'
815
816
// Declare derivative matrix (of polynomial basis).
817
double dmats[3][3] = \
818
{{1.0, 0.0, 0.0},
819
{0.0, 1.0, 0.0},
820
{0.0, 0.0, 1.0}};
821
822
// Declare (auxiliary) derivative matrix (of polynomial basis).
823
double dmats_old[3][3] = \
824
{{1.0, 0.0, 0.0},
825
{0.0, 1.0, 0.0},
826
{0.0, 0.0, 1.0}};
827
828
// Loop possible derivatives.
829
for (unsigned int r = 0; r < num_derivatives_t; r++)
830
{
831
// Resetting dmats values to compute next derivative.
832
for (unsigned int t = 0; t < 3; t++)
833
{
834
for (unsigned int u = 0; u < 3; u++)
835
{
836
dmats[t][u] = 0.0;
837
if (t == u)
838
{
839
dmats[t][u] = 1.0;
840
}
841
842
}// end loop over 'u'
843
}// end loop over 't'
844
845
// Looping derivative order to generate dmats.
846
for (unsigned int s = 0; s < n; s++)
847
{
848
// Updating dmats_old with new values and resetting dmats.
849
for (unsigned int t = 0; t < 3; t++)
850
{
851
for (unsigned int u = 0; u < 3; u++)
852
{
853
dmats_old[t][u] = dmats[t][u];
854
dmats[t][u] = 0.0;
855
}// end loop over 'u'
856
}// end loop over 't'
857
858
// Update dmats using an inner product.
859
if (combinations_t[r][s] == 0)
860
{
861
for (unsigned int t = 0; t < 3; t++)
862
{
863
for (unsigned int u = 0; u < 3; u++)
864
{
865
for (unsigned int tu = 0; tu < 3; tu++)
866
{
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'
871
}
872
873
if (combinations_t[r][s] == 1)
874
{
875
for (unsigned int t = 0; t < 3; t++)
876
{
877
for (unsigned int u = 0; u < 3; u++)
878
{
879
for (unsigned int tu = 0; tu < 3; tu++)
880
{
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'
885
}
886
887
}// end loop over 's'
888
for (unsigned int s = 0; s < 3; s++)
889
{
890
for (unsigned int t = 0; t < 3; t++)
891
{
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'
896
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'
904
905
// Transform derivatives back to physical element
906
for (unsigned int r = 0; r < num_derivatives_g; r++)
907
{
908
for (unsigned int s = 0; s < num_derivatives_t; s++)
909
{
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'
915
916
// Delete pointer to array of derivatives on FIAT element
917
delete [] derivatives;
918
919
// Delete pointer to array of reference derivatives on physical element.
920
delete [] derivatives_p;
921
922
// Delete pointer to array of combinations of derivatives and transform
923
for (unsigned int r = 0; r < num_derivatives_t; r++)
924
{
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++)
929
{
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++)
934
{
935
delete [] transform[r];
936
}// end loop over 'r'
937
delete [] transform;
938
break;
939
}
940
}
941
942
}
943
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,
946
double* values,
947
const double* x,
948
const double* vertex_coordinates,
949
int cell_orientation) const
950
{
951
// Compute number of derivatives.
952
unsigned int num_derivatives_g = 1;
953
for (unsigned int r = 0; r < n; r++)
954
{
955
num_derivatives_g *= 3;
956
}// end loop over 'r'
957
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++)
961
{
962
dof_values[r] = 0.0;
963
}// end loop over 'r'
964
965
// Loop dofs and call evaluate_basis_derivatives.
966
for (unsigned int r = 0; r < 3; r++)
967
{
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++)
970
{
971
values[r*3*num_derivatives_g + s] = dof_values[s];
972
}// end loop over 's'
973
}// end loop over 'r'
974
975
// Delete pointer.
976
delete [] dof_values;
977
}
978
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
985
{
986
// Declare variables for result of evaluation
987
double vals[3];
988
989
// Declare variable for physical coordinates
990
double y[3];
991
992
double result;
993
// Compute Jacobian
994
double J[6];
995
compute_jacobian_triangle_3d(J, vertex_coordinates);
996
997
998
// Compute Jacobian inverse and determinant
999
double K[6];
1000
double detJ;
1001
compute_jacobian_inverse_triangle_3d(K, detJ, J);
1002
1003
1004
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)
1010
detJ *= -1;
1011
switch (i)
1012
{
1013
case 0:
1014
{
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]));
1020
return result;
1021
break;
1022
}
1023
case 1:
1024
{
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]));
1030
return result;
1031
break;
1032
}
1033
case 2:
1034
{
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]));
1040
return result;
1041
break;
1042
}
1043
}
1044
1045
return 0.0;
1046
}
1047
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
1054
{
1055
// Declare variables for result of evaluation
1056
double vals[3];
1057
1058
// Declare variable for physical coordinates
1059
double y[3];
1060
1061
double result;
1062
// Compute Jacobian
1063
double J[6];
1064
compute_jacobian_triangle_3d(J, vertex_coordinates);
1065
1066
1067
// Compute Jacobian inverse and determinant
1068
double K[6];
1069
double detJ;
1070
compute_jacobian_inverse_triangle_3d(K, detJ, J);
1071
1072
1073
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)
1079
detJ *= -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]));
1085
values[0] = result;
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]));
1091
values[1] = result;
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]));
1097
values[2] = result;
1098
}
1099
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
1106
{
1107
// Compute Jacobian
1108
double J[6];
1109
compute_jacobian_triangle_3d(J, vertex_coordinates);
1110
1111
1112
// Compute Jacobian inverse and determinant
1113
double K[6];
1114
double detJ;
1115
compute_jacobian_inverse_triangle_3d(K, detJ, J);
1116
1117
1118
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)
1124
detJ *= -1;
1125
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)));
1133
}
1134
1135
/// Map coordinate xhat from reference cell to coordinate x in cell
1136
virtual void map_from_reference_cell(double* x,
1137
const double* xhat,
1138
const ufc::cell& c) const
1139
{
1140
std::cerr << "*** FFC warning: " << "map_from_reference_cell not yet implemented." << std::endl;
1141
}
1142
1143
/// Map from coordinate x in cell to coordinate xhat in reference cell
1144
virtual void map_to_reference_cell(double* xhat,
1145
const double* x,
1146
const ufc::cell& c) const
1147
{
1148
std::cerr << "*** FFC warning: " << "map_to_reference_cell not yet implemented." << std::endl;
1149
}
1150
1151
/// Return the number of sub elements (for a mixed element)
1152
virtual std::size_t num_sub_elements() const
1153
{
1154
return 0;
1155
}
1156
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
1159
{
1160
return 0;
1161
}
1162
1163
/// Create a new class instance
1164
virtual ufc::finite_element* create() const
1165
{
1166
return new projectionmanifold_finite_element_0();
1167
}
1168
1169
};
1170
1171
/// This class defines the interface for a finite element.
1172
1173
class projectionmanifold_finite_element_1: public ufc::finite_element
1174
{
1175
public:
1176
1177
/// Constructor
1178
projectionmanifold_finite_element_1() : ufc::finite_element()
1179
{
1180
// Do nothing
1181
}
1182
1183
/// Destructor
1184
virtual ~projectionmanifold_finite_element_1()
1185
{
1186
// Do nothing
1187
}
1188
1189
/// Return a string identifying the finite element
1190
virtual const char* signature() const
1191
{
1192
return "FiniteElement('Discontinuous Lagrange', Domain(Cell('triangle', 3), 'triangle_multiverse', 3, 2), 0, None)";
1193
}
1194
1195
/// Return the cell shape
1196
virtual ufc::shape cell_shape() const
1197
{
1198
return ufc::triangle;
1199
}
1200
1201
/// Return the topological dimension of the cell shape
1202
virtual std::size_t topological_dimension() const
1203
{
1204
return 2;
1205
}
1206
1207
/// Return the geometric dimension of the cell shape
1208
virtual std::size_t geometric_dimension() const
1209
{
1210
return 3;
1211
}
1212
1213
/// Return the dimension of the finite element function space
1214
virtual std::size_t space_dimension() const
1215
{
1216
return 1;
1217
}
1218
1219
/// Return the rank of the value space
1220
virtual std::size_t value_rank() const
1221
{
1222
return 0;
1223
}
1224
1225
/// Return the dimension of the value space for axis i
1226
virtual std::size_t value_dimension(std::size_t i) const
1227
{
1228
return 1;
1229
}
1230
1231
/// Evaluate basis function i at given point x in cell
1232
virtual void evaluate_basis(std::size_t i,
1233
double* values,
1234
const double* x,
1235
const double* vertex_coordinates,
1236
int cell_orientation) const
1237
{
1238
// Compute Jacobian
1239
double J[6];
1240
compute_jacobian_triangle_3d(J, vertex_coordinates);
1241
1242
// Compute Jacobian inverse and determinant
1243
double K[6];
1244
double detJ;
1245
compute_jacobian_inverse_triangle_3d(K, detJ, J);
1246
1247
1248
1249
// P_FFC = J^dag (p - b), P_FIAT = 2*P_FFC - (1, 1)
1250
1251
1252
// Reset values
1253
*values = 0.0;
1254
1255
// Array of basisvalues
1256
double basisvalues[1] = {0.0};
1257
1258
// Declare helper variables
1259
1260
// Compute basisvalues
1261
basisvalues[0] = 1.0;
1262
1263
// Table(s) of coefficients
1264
static const double coefficients0[1] = \
1265
{1.0};
1266
1267
// Compute value(s)
1268
for (unsigned int r = 0; r < 1; r++)
1269
{
1270
*values += coefficients0[r]*basisvalues[r];
1271
}// end loop over 'r'
1272
}
1273
1274
/// Evaluate all basis functions at given point x in cell
1275
virtual void evaluate_basis_all(double* values,
1276
const double* x,
1277
const double* vertex_coordinates,
1278
int cell_orientation) const
1279
{
1280
// Element is constant, calling evaluate_basis.
1281
evaluate_basis(0, values, x, vertex_coordinates, cell_orientation);
1282
}
1283
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,
1286
std::size_t n,
1287
double* values,
1288
const double* x,
1289
const double* vertex_coordinates,
1290
int cell_orientation) const
1291
{
1292
// Compute Jacobian
1293
double J[6];
1294
compute_jacobian_triangle_3d(J, vertex_coordinates);
1295
1296
// Compute Jacobian inverse and determinant
1297
double K[6];
1298
double detJ;
1299
compute_jacobian_inverse_triangle_3d(K, detJ, J);
1300
1301
1302
1303
// P_FFC = J^dag (p - b), P_FIAT = 2*P_FFC - (1, 1)
1304
1305
1306
// Compute number of derivatives.
1307
unsigned int num_derivatives_t = 1;
1308
for (unsigned int r = 0; r < n; r++)
1309
{
1310
num_derivatives_t *= 2;
1311
}// end loop over 'r'
1312
1313
// Compute number of derivatives.
1314
unsigned int num_derivatives_g = 1;
1315
for (unsigned int r = 0; r < n; r++)
1316
{
1317
num_derivatives_g *= 3;
1318
}// end loop over 'r'
1319
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++)
1323
{
1324
combinations_t[row] = new unsigned int [n];
1325
for (unsigned int col = 0; col < n; col++)
1326
combinations_t[row][col] = 0;
1327
}
1328
1329
// Generate combinations of derivatives
1330
for (unsigned int row = 1; row < num_derivatives_t; row++)
1331
{
1332
for (unsigned int num = 0; num < row; num++)
1333
{
1334
for (unsigned int col = n-1; col+1 > 0; col--)
1335
{
1336
if (combinations_t[row][col] + 1 > 1)
1337
combinations_t[row][col] = 0;
1338
else
1339
{
1340
combinations_t[row][col] += 1;
1341
break;
1342
}
1343
}
1344
}
1345
}
1346
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++)
1350
{
1351
combinations_g[row] = new unsigned int [n];
1352
for (unsigned int col = 0; col < n; col++)
1353
combinations_g[row][col] = 0;
1354
}
1355
1356
// Generate combinations of derivatives
1357
for (unsigned int row = 1; row < num_derivatives_g; row++)
1358
{
1359
for (unsigned int num = 0; num < row; num++)
1360
{
1361
for (unsigned int col = n-1; col+1 > 0; col--)
1362
{
1363
if (combinations_g[row][col] + 1 > 2)
1364
combinations_g[row][col] = 0;
1365
else
1366
{
1367
combinations_g[row][col] += 1;
1368
break;
1369
}
1370
}
1371
}
1372
}
1373
1374
// Compute inverse of Jacobian
1375
const double Jinv[2][3] = {{K[0], K[1], K[2]}, {K[3], K[4], K[5]}};
1376
1377
// Declare transformation matrix
1378
// Declare pointer to two dimensional array and initialise
1379
double **transform = new double *[num_derivatives_g];
1380
1381
for (unsigned int j = 0; j < num_derivatives_g; j++)
1382
{
1383
transform[j] = new double [num_derivatives_t];
1384
for (unsigned int k = 0; k < num_derivatives_t; k++)
1385
transform[j][k] = 1;
1386
}
1387
1388
// Construct transformation matrix
1389
for (unsigned int row = 0; row < num_derivatives_g; row++)
1390
{
1391
for (unsigned int col = 0; col < num_derivatives_t; col++)
1392
{
1393
for (unsigned int k = 0; k < n; k++)
1394
transform[row][col] *= Jinv[combinations_t[col][k]][combinations_g[row][k]];
1395
}
1396
}
1397
1398
// Reset values. Assuming that values is always an array.
1399
for (unsigned int r = 0; r < num_derivatives_g; r++)
1400
{
1401
values[r] = 0.0;
1402
}// end loop over 'r'
1403
1404
1405
// Array of basisvalues
1406
double basisvalues[1] = {0.0};
1407
1408
// Declare helper variables
1409
1410
// Compute basisvalues
1411
basisvalues[0] = 1.0;
1412
1413
// Table(s) of coefficients
1414
static const double coefficients0[1] = \
1415
{1.0};
1416
1417
// Tables of derivatives of the polynomial base (transpose).
1418
static const double dmats0[1][1] = \
1419
{{0.0}};
1420
1421
static const double dmats1[1][1] = \
1422
{{0.0}};
1423
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++)
1428
{
1429
derivatives[r] = 0.0;
1430
}// end loop over 'r'
1431
1432
// Declare derivative matrix (of polynomial basis).
1433
double dmats[1][1] = \
1434
{{1.0}};
1435
1436
// Declare (auxiliary) derivative matrix (of polynomial basis).
1437
double dmats_old[1][1] = \
1438
{{1.0}};
1439
1440
// Loop possible derivatives.
1441
for (unsigned int r = 0; r < num_derivatives_t; r++)
1442
{
1443
// Resetting dmats values to compute next derivative.
1444
for (unsigned int t = 0; t < 1; t++)
1445
{
1446
for (unsigned int u = 0; u < 1; u++)
1447
{
1448
dmats[t][u] = 0.0;
1449
if (t == u)
1450
{
1451
dmats[t][u] = 1.0;
1452
}
1453
1454
}// end loop over 'u'
1455
}// end loop over 't'
1456
1457
// Looping derivative order to generate dmats.
1458
for (unsigned int s = 0; s < n; s++)
1459
{
1460
// Updating dmats_old with new values and resetting dmats.
1461
for (unsigned int t = 0; t < 1; t++)
1462
{
1463
for (unsigned int u = 0; u < 1; u++)
1464
{
1465
dmats_old[t][u] = dmats[t][u];
1466
dmats[t][u] = 0.0;
1467
}// end loop over 'u'
1468
}// end loop over 't'
1469
1470
// Update dmats using an inner product.
1471
if (combinations_t[r][s] == 0)
1472
{
1473
for (unsigned int t = 0; t < 1; t++)
1474
{
1475
for (unsigned int u = 0; u < 1; u++)
1476
{
1477
for (unsigned int tu = 0; tu < 1; tu++)
1478
{
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'
1483
}
1484
1485
if (combinations_t[r][s] == 1)
1486
{
1487
for (unsigned int t = 0; t < 1; t++)
1488
{
1489
for (unsigned int u = 0; u < 1; u++)
1490
{
1491
for (unsigned int tu = 0; tu < 1; tu++)
1492
{
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'
1497
}
1498
1499
}// end loop over 's'
1500
for (unsigned int s = 0; s < 1; s++)
1501
{
1502
for (unsigned int t = 0; t < 1; t++)
1503
{
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'
1508
1509
// Transform derivatives back to physical element
1510
for (unsigned int r = 0; r < num_derivatives_g; r++)
1511
{
1512
for (unsigned int s = 0; s < num_derivatives_t; s++)
1513
{
1514
values[r] += transform[r][s]*derivatives[s];
1515
}// end loop over 's'
1516
}// end loop over 'r'
1517
1518
// Delete pointer to array of derivatives on FIAT element
1519
delete [] derivatives;
1520
1521
// Delete pointer to array of combinations of derivatives and transform
1522
for (unsigned int r = 0; r < num_derivatives_t; r++)
1523
{
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++)
1528
{
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++)
1533
{
1534
delete [] transform[r];
1535
}// end loop over 'r'
1536
delete [] transform;
1537
}
1538
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,
1541
double* values,
1542
const double* x,
1543
const double* vertex_coordinates,
1544
int cell_orientation) const
1545
{
1546
// Element is constant, calling evaluate_basis_derivatives.
1547
evaluate_basis_derivatives(0, n, values, x, vertex_coordinates, cell_orientation);
1548
}
1549
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
1556
{
1557
// Declare variables for result of evaluation
1558
double vals[1];
1559
1560
// Declare variable for physical coordinates
1561
double y[3];
1562
switch (i)
1563
{
1564
case 0:
1565
{
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);
1570
return vals[0];
1571
break;
1572
}
1573
}
1574
1575
return 0.0;
1576
}
1577
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
1584
{
1585
// Declare variables for result of evaluation
1586
double vals[1];
1587
1588
// Declare variable for physical coordinates
1589
double y[3];
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];
1595
}
1596
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
1603
{
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];
1608
}
1609
1610
/// Map coordinate xhat from reference cell to coordinate x in cell
1611
virtual void map_from_reference_cell(double* x,
1612
const double* xhat,
1613
const ufc::cell& c) const
1614
{
1615
std::cerr << "*** FFC warning: " << "map_from_reference_cell not yet implemented." << std::endl;
1616
}
1617
1618
/// Map from coordinate x in cell to coordinate xhat in reference cell
1619
virtual void map_to_reference_cell(double* xhat,
1620
const double* x,
1621
const ufc::cell& c) const
1622
{
1623
std::cerr << "*** FFC warning: " << "map_to_reference_cell not yet implemented." << std::endl;
1624
}
1625
1626
/// Return the number of sub elements (for a mixed element)
1627
virtual std::size_t num_sub_elements() const
1628
{
1629
return 0;
1630
}
1631
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
1634
{
1635
return 0;
1636
}
1637
1638
/// Create a new class instance
1639
virtual ufc::finite_element* create() const
1640
{
1641
return new projectionmanifold_finite_element_1();
1642
}
1643
1644
};
1645
1646
/// This class defines the interface for a finite element.
1647
1648
class projectionmanifold_finite_element_2: public ufc::finite_element
1649
{
1650
public:
1651
1652
/// Constructor
1653
projectionmanifold_finite_element_2() : ufc::finite_element()
1654
{
1655
// Do nothing
1656
}
1657
1658
/// Destructor
1659
virtual ~projectionmanifold_finite_element_2()
1660
{
1661
// Do nothing
1662
}
1663
1664
/// Return a string identifying the finite element
1665
virtual const char* signature() const
1666
{
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,) })";
1668
}
1669
1670
/// Return the cell shape
1671
virtual ufc::shape cell_shape() const
1672
{
1673
return ufc::triangle;
1674
}
1675
1676
/// Return the topological dimension of the cell shape
1677
virtual std::size_t topological_dimension() const
1678
{
1679
return 2;
1680
}
1681
1682
/// Return the geometric dimension of the cell shape
1683
virtual std::size_t geometric_dimension() const
1684
{
1685
return 3;
1686
}
1687
1688
/// Return the dimension of the finite element function space
1689
virtual std::size_t space_dimension() const
1690
{
1691
return 4;
1692
}
1693
1694
/// Return the rank of the value space
1695
virtual std::size_t value_rank() const
1696
{
1697
return 1;
1698
}
1699
1700
/// Return the dimension of the value space for axis i
1701
virtual std::size_t value_dimension(std::size_t i) const
1702
{
1703
switch (i)
1704
{
1705
case 0:
1706
{
1707
return 4;
1708
break;
1709
}
1710
}
1711
1712
return 0;
1713
}
1714
1715
/// Evaluate basis function i at given point x in cell
1716
virtual void evaluate_basis(std::size_t i,
1717
double* values,
1718
const double* x,
1719
const double* vertex_coordinates,
1720
int cell_orientation) const
1721
{
1722
// Compute Jacobian
1723
double J[6];
1724
compute_jacobian_triangle_3d(J, vertex_coordinates);
1725
1726
// Compute Jacobian inverse and determinant
1727
double K[6];
1728
double detJ;
1729
compute_jacobian_inverse_triangle_3d(K, detJ, J);
1730
1731
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)
1737
detJ *= -1;
1738
1739
1740
const double b0 = vertex_coordinates[0];
1741
const double b1 = vertex_coordinates[1];
1742
const double b2 = vertex_coordinates[2];
1743
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;
1747
1748
1749
// Reset values
1750
values[0] = 0.0;
1751
values[1] = 0.0;
1752
values[2] = 0.0;
1753
switch (i)
1754
{
1755
case 0:
1756
{
1757
1758
// Array of basisvalues
1759
double basisvalues[3] = {0.0, 0.0, 0.0};
1760
1761
// Declare helper variables
1762
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1763
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);
1771
1772
// Table(s) of coefficients
1773
static const double coefficients0[3] = \
1774
{0.47140452, 0.28867513, -0.16666667};
1775
1776
static const double coefficients1[3] = \
1777
{0.47140452, 0.0, 0.33333333};
1778
1779
// Compute value(s)
1780
for (unsigned int r = 0; r < 3; r++)
1781
{
1782
values[0] += coefficients0[r]*basisvalues[r];
1783
values[1] += coefficients1[r]*basisvalues[r];
1784
}// end loop over 'r'
1785
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);
1792
break;
1793
}
1794
case 1:
1795
{
1796
1797
// Array of basisvalues
1798
double basisvalues[3] = {0.0, 0.0, 0.0};
1799
1800
// Declare helper variables
1801
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1802
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);
1810
1811
// Table(s) of coefficients
1812
static const double coefficients0[3] = \
1813
{0.94280904, -0.28867513, 0.16666667};
1814
1815
static const double coefficients1[3] = \
1816
{-0.47140452, 0.0, -0.33333333};
1817
1818
// Compute value(s)
1819
for (unsigned int r = 0; r < 3; r++)
1820
{
1821
values[0] += coefficients0[r]*basisvalues[r];
1822
values[1] += coefficients1[r]*basisvalues[r];
1823
}// end loop over 'r'
1824
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);
1831
break;
1832
}
1833
case 2:
1834
{
1835
1836
// Array of basisvalues
1837
double basisvalues[3] = {0.0, 0.0, 0.0};
1838
1839
// Declare helper variables
1840
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1841
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);
1849
1850
// Table(s) of coefficients
1851
static const double coefficients0[3] = \
1852
{0.47140452, 0.28867513, -0.16666667};
1853
1854
static const double coefficients1[3] = \
1855
{-0.94280904, 0.0, 0.33333333};
1856
1857
// Compute value(s)
1858
for (unsigned int r = 0; r < 3; r++)
1859
{
1860
values[0] += coefficients0[r]*basisvalues[r];
1861
values[1] += coefficients1[r]*basisvalues[r];
1862
}// end loop over 'r'
1863
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);
1870
break;
1871
}
1872
case 3:
1873
{
1874
1875
// Array of basisvalues
1876
double basisvalues[1] = {0.0};
1877
1878
// Declare helper variables
1879
1880
// Compute basisvalues
1881
basisvalues[0] = 1.0;
1882
1883
// Table(s) of coefficients
1884
static const double coefficients0[1] = \
1885
{1.0};
1886
1887
// Compute value(s)
1888
for (unsigned int r = 0; r < 1; r++)
1889
{
1890
values[2] += coefficients0[r]*basisvalues[r];
1891
}// end loop over 'r'
1892
break;
1893
}
1894
}
1895
1896
}
1897
1898
/// Evaluate all basis functions at given point x in cell
1899
virtual void evaluate_basis_all(double* values,
1900
const double* x,
1901
const double* vertex_coordinates,
1902
int cell_orientation) const
1903
{
1904
// Helper variable to hold values of a single dof.
1905
double dof_values[4] = {0.0, 0.0, 0.0, 0.0};
1906
1907
// Loop dofs and call evaluate_basis
1908
for (unsigned int r = 0; r < 4; r++)
1909
{
1910
evaluate_basis(r, dof_values, x, vertex_coordinates, cell_orientation);
1911
for (unsigned int s = 0; s < 4; s++)
1912
{
1913
values[r*4 + s] = dof_values[s];
1914
}// end loop over 's'
1915
}// end loop over 'r'
1916
}
1917
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,
1920
std::size_t n,
1921
double* values,
1922
const double* x,
1923
const double* vertex_coordinates,
1924
int cell_orientation) const
1925
{
1926
// Compute Jacobian
1927
double J[6];
1928
compute_jacobian_triangle_3d(J, vertex_coordinates);
1929
1930
// Compute Jacobian inverse and determinant
1931
double K[6];
1932
double detJ;
1933
compute_jacobian_inverse_triangle_3d(K, detJ, J);
1934
1935
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)
1941
detJ *= -1;
1942
1943
1944
const double b0 = vertex_coordinates[0];
1945
const double b1 = vertex_coordinates[1];
1946
const double b2 = vertex_coordinates[2];
1947
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;
1951
1952
1953
// Compute number of derivatives.
1954
unsigned int num_derivatives_t = 1;
1955
for (unsigned int r = 0; r < n; r++)
1956
{
1957
num_derivatives_t *= 2;
1958
}// end loop over 'r'
1959
1960
// Compute number of derivatives.
1961
unsigned int num_derivatives_g = 1;
1962
for (unsigned int r = 0; r < n; r++)
1963
{
1964
num_derivatives_g *= 3;
1965
}// end loop over 'r'
1966
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++)
1970
{
1971
combinations_t[row] = new unsigned int [n];
1972
for (unsigned int col = 0; col < n; col++)
1973
combinations_t[row][col] = 0;
1974
}
1975
1976
// Generate combinations of derivatives
1977
for (unsigned int row = 1; row < num_derivatives_t; row++)
1978
{
1979
for (unsigned int num = 0; num < row; num++)
1980
{
1981
for (unsigned int col = n-1; col+1 > 0; col--)
1982
{
1983
if (combinations_t[row][col] + 1 > 1)
1984
combinations_t[row][col] = 0;
1985
else
1986
{
1987
combinations_t[row][col] += 1;
1988
break;
1989
}
1990
}
1991
}
1992
}
1993
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++)
1997
{
1998
combinations_g[row] = new unsigned int [n];
1999
for (unsigned int col = 0; col < n; col++)
2000
combinations_g[row][col] = 0;
2001
}
2002
2003
// Generate combinations of derivatives
2004
for (unsigned int row = 1; row < num_derivatives_g; row++)
2005
{
2006
for (unsigned int num = 0; num < row; num++)
2007
{
2008
for (unsigned int col = n-1; col+1 > 0; col--)
2009
{
2010
if (combinations_g[row][col] + 1 > 2)
2011
combinations_g[row][col] = 0;
2012
else
2013
{
2014
combinations_g[row][col] += 1;
2015
break;
2016
}
2017
}
2018
}
2019
}
2020
2021
// Compute inverse of Jacobian
2022
const double Jinv[2][3] = {{K[0], K[1], K[2]}, {K[3], K[4], K[5]}};
2023
2024
// Declare transformation matrix
2025
// Declare pointer to two dimensional array and initialise
2026
double **transform = new double *[num_derivatives_g];
2027
2028
for (unsigned int j = 0; j < num_derivatives_g; j++)
2029
{
2030
transform[j] = new double [num_derivatives_t];
2031
for (unsigned int k = 0; k < num_derivatives_t; k++)
2032
transform[j][k] = 1;
2033
}
2034
2035
// Construct transformation matrix
2036
for (unsigned int row = 0; row < num_derivatives_g; row++)
2037
{
2038
for (unsigned int col = 0; col < num_derivatives_t; col++)
2039
{
2040
for (unsigned int k = 0; k < n; k++)
2041
transform[row][col] *= Jinv[combinations_t[col][k]][combinations_g[row][k]];
2042
}
2043
}
2044
2045
// Reset values. Assuming that values is always an array.
2046
for (unsigned int r = 0; r < 4*num_derivatives_g; r++)
2047
{
2048
values[r] = 0.0;
2049
}// end loop over 'r'
2050
2051
switch (i)
2052
{
2053
case 0:
2054
{
2055
2056
// Array of basisvalues
2057
double basisvalues[3] = {0.0, 0.0, 0.0};
2058
2059
// Declare helper variables
2060
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2061
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);
2069
2070
// Table(s) of coefficients
2071
static const double coefficients0[3] = \
2072
{0.47140452, 0.28867513, -0.16666667};
2073
2074
static const double coefficients1[3] = \
2075
{0.47140452, 0.0, 0.33333333};
2076
2077
// Tables of derivatives of the polynomial base (transpose).
2078
static const double dmats0[3][3] = \
2079
{{0.0, 0.0, 0.0},
2080
{4.8989795, 0.0, 0.0},
2081
{0.0, 0.0, 0.0}};
2082
2083
static const double dmats1[3][3] = \
2084
{{0.0, 0.0, 0.0},
2085
{2.4494897, 0.0, 0.0},
2086
{4.2426407, 0.0, 0.0}};
2087
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++)
2092
{
2093
derivatives[r] = 0.0;
2094
}// end loop over 'r'
2095
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++)
2099
{
2100
derivatives_p[r] = 0.0;
2101
}// end loop over 'r'
2102
2103
// Declare derivative matrix (of polynomial basis).
2104
double dmats[3][3] = \
2105
{{1.0, 0.0, 0.0},
2106
{0.0, 1.0, 0.0},
2107
{0.0, 0.0, 1.0}};
2108
2109
// Declare (auxiliary) derivative matrix (of polynomial basis).
2110
double dmats_old[3][3] = \
2111
{{1.0, 0.0, 0.0},
2112
{0.0, 1.0, 0.0},
2113
{0.0, 0.0, 1.0}};
2114
2115
// Loop possible derivatives.
2116
for (unsigned int r = 0; r < num_derivatives_t; r++)
2117
{
2118
// Resetting dmats values to compute next derivative.
2119
for (unsigned int t = 0; t < 3; t++)
2120
{
2121
for (unsigned int u = 0; u < 3; u++)
2122
{
2123
dmats[t][u] = 0.0;
2124
if (t == u)
2125
{
2126
dmats[t][u] = 1.0;
2127
}
2128
2129
}// end loop over 'u'
2130
}// end loop over 't'
2131
2132
// Looping derivative order to generate dmats.
2133
for (unsigned int s = 0; s < n; s++)
2134
{
2135
// Updating dmats_old with new values and resetting dmats.
2136
for (unsigned int t = 0; t < 3; t++)
2137
{
2138
for (unsigned int u = 0; u < 3; u++)
2139
{
2140
dmats_old[t][u] = dmats[t][u];
2141
dmats[t][u] = 0.0;
2142
}// end loop over 'u'
2143
}// end loop over 't'
2144
2145
// Update dmats using an inner product.
2146
if (combinations_t[r][s] == 0)
2147
{
2148
for (unsigned int t = 0; t < 3; t++)
2149
{
2150
for (unsigned int u = 0; u < 3; u++)
2151
{
2152
for (unsigned int tu = 0; tu < 3; tu++)
2153
{
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'
2158
}
2159
2160
if (combinations_t[r][s] == 1)
2161
{
2162
for (unsigned int t = 0; t < 3; t++)
2163
{
2164
for (unsigned int u = 0; u < 3; u++)
2165
{
2166
for (unsigned int tu = 0; tu < 3; tu++)
2167
{
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'
2172
}
2173
2174
}// end loop over 's'
2175
for (unsigned int s = 0; s < 3; s++)
2176
{
2177
for (unsigned int t = 0; t < 3; t++)
2178
{
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'
2183
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'
2191
2192
// Transform derivatives back to physical element
2193
for (unsigned int r = 0; r < num_derivatives_g; r++)
2194
{
2195
for (unsigned int s = 0; s < num_derivatives_t; s++)
2196
{
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'
2202
2203
// Delete pointer to array of derivatives on FIAT element
2204
delete [] derivatives;
2205
2206
// Delete pointer to array of reference derivatives on physical element.
2207
delete [] derivatives_p;
2208
2209
// Delete pointer to array of combinations of derivatives and transform
2210
for (unsigned int r = 0; r < num_derivatives_t; r++)
2211
{
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++)
2216
{
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++)
2221
{
2222
delete [] transform[r];
2223
}// end loop over 'r'
2224
delete [] transform;
2225
break;
2226
}
2227
case 1:
2228
{
2229
2230
// Array of basisvalues
2231
double basisvalues[3] = {0.0, 0.0, 0.0};
2232
2233
// Declare helper variables
2234
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2235
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);
2243
2244
// Table(s) of coefficients
2245
static const double coefficients0[3] = \
2246
{0.94280904, -0.28867513, 0.16666667};
2247
2248
static const double coefficients1[3] = \
2249
{-0.47140452, 0.0, -0.33333333};
2250
2251
// Tables of derivatives of the polynomial base (transpose).
2252
static const double dmats0[3][3] = \
2253
{{0.0, 0.0, 0.0},
2254
{4.8989795, 0.0, 0.0},
2255
{0.0, 0.0, 0.0}};
2256
2257
static const double dmats1[3][3] = \
2258
{{0.0, 0.0, 0.0},
2259
{2.4494897, 0.0, 0.0},
2260
{4.2426407, 0.0, 0.0}};
2261
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++)
2266
{
2267
derivatives[r] = 0.0;
2268
}// end loop over 'r'
2269
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++)
2273
{
2274
derivatives_p[r] = 0.0;
2275
}// end loop over 'r'
2276
2277
// Declare derivative matrix (of polynomial basis).
2278
double dmats[3][3] = \
2279
{{1.0, 0.0, 0.0},
2280
{0.0, 1.0, 0.0},
2281
{0.0, 0.0, 1.0}};
2282
2283
// Declare (auxiliary) derivative matrix (of polynomial basis).
2284
double dmats_old[3][3] = \
2285
{{1.0, 0.0, 0.0},
2286
{0.0, 1.0, 0.0},
2287
{0.0, 0.0, 1.0}};
2288
2289
// Loop possible derivatives.
2290
for (unsigned int r = 0; r < num_derivatives_t; r++)
2291
{
2292
// Resetting dmats values to compute next derivative.
2293
for (unsigned int t = 0; t < 3; t++)
2294
{
2295
for (unsigned int u = 0; u < 3; u++)
2296
{
2297
dmats[t][u] = 0.0;
2298
if (t == u)
2299
{
2300
dmats[t][u] = 1.0;
2301
}
2302
2303
}// end loop over 'u'
2304
}// end loop over 't'
2305
2306
// Looping derivative order to generate dmats.
2307
for (unsigned int s = 0; s < n; s++)
2308
{
2309
// Updating dmats_old with new values and resetting dmats.
2310
for (unsigned int t = 0; t < 3; t++)
2311
{
2312
for (unsigned int u = 0; u < 3; u++)
2313
{
2314
dmats_old[t][u] = dmats[t][u];
2315
dmats[t][u] = 0.0;
2316
}// end loop over 'u'
2317
}// end loop over 't'
2318
2319
// Update dmats using an inner product.
2320
if (combinations_t[r][s] == 0)
2321
{
2322
for (unsigned int t = 0; t < 3; t++)
2323
{
2324
for (unsigned int u = 0; u < 3; u++)
2325
{
2326
for (unsigned int tu = 0; tu < 3; tu++)
2327
{
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'
2332
}
2333
2334
if (combinations_t[r][s] == 1)
2335
{
2336
for (unsigned int t = 0; t < 3; t++)
2337
{
2338
for (unsigned int u = 0; u < 3; u++)
2339
{
2340
for (unsigned int tu = 0; tu < 3; tu++)
2341
{
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'
2346
}
2347
2348
}// end loop over 's'
2349
for (unsigned int s = 0; s < 3; s++)
2350
{
2351
for (unsigned int t = 0; t < 3; t++)
2352
{
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'
2357
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'
2365
2366
// Transform derivatives back to physical element
2367
for (unsigned int r = 0; r < num_derivatives_g; r++)
2368
{
2369
for (unsigned int s = 0; s < num_derivatives_t; s++)
2370
{
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'
2376
2377
// Delete pointer to array of derivatives on FIAT element
2378
delete [] derivatives;
2379
2380
// Delete pointer to array of reference derivatives on physical element.
2381
delete [] derivatives_p;
2382
2383
// Delete pointer to array of combinations of derivatives and transform
2384
for (unsigned int r = 0; r < num_derivatives_t; r++)
2385
{
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++)
2390
{
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++)
2395
{
2396
delete [] transform[r];
2397
}// end loop over 'r'
2398
delete [] transform;
2399
break;
2400
}
2401
case 2:
2402
{
2403
2404
// Array of basisvalues
2405
double basisvalues[3] = {0.0, 0.0, 0.0};
2406
2407
// Declare helper variables
2408
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2409
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);
2417
2418
// Table(s) of coefficients
2419
static const double coefficients0[3] = \
2420
{0.47140452, 0.28867513, -0.16666667};
2421
2422
static const double coefficients1[3] = \
2423
{-0.94280904, 0.0, 0.33333333};
2424
2425
// Tables of derivatives of the polynomial base (transpose).
2426
static const double dmats0[3][3] = \
2427
{{0.0, 0.0, 0.0},
2428
{4.8989795, 0.0, 0.0},
2429
{0.0, 0.0, 0.0}};
2430
2431
static const double dmats1[3][3] = \
2432
{{0.0, 0.0, 0.0},
2433
{2.4494897, 0.0, 0.0},
2434
{4.2426407, 0.0, 0.0}};
2435
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++)
2440
{
2441
derivatives[r] = 0.0;
2442
}// end loop over 'r'
2443
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++)
2447
{
2448
derivatives_p[r] = 0.0;
2449
}// end loop over 'r'
2450
2451
// Declare derivative matrix (of polynomial basis).
2452
double dmats[3][3] = \
2453
{{1.0, 0.0, 0.0},
2454
{0.0, 1.0, 0.0},
2455
{0.0, 0.0, 1.0}};
2456
2457
// Declare (auxiliary) derivative matrix (of polynomial basis).
2458
double dmats_old[3][3] = \
2459
{{1.0, 0.0, 0.0},
2460
{0.0, 1.0, 0.0},
2461
{0.0, 0.0, 1.0}};
2462
2463
// Loop possible derivatives.
2464
for (unsigned int r = 0; r < num_derivatives_t; r++)
2465
{
2466
// Resetting dmats values to compute next derivative.
2467
for (unsigned int t = 0; t < 3; t++)
2468
{
2469
for (unsigned int u = 0; u < 3; u++)
2470
{
2471
dmats[t][u] = 0.0;
2472
if (t == u)
2473
{
2474
dmats[t][u] = 1.0;
2475
}
2476
2477
}// end loop over 'u'
2478
}// end loop over 't'
2479
2480
// Looping derivative order to generate dmats.
2481
for (unsigned int s = 0; s < n; s++)
2482
{
2483
// Updating dmats_old with new values and resetting dmats.
2484
for (unsigned int t = 0; t < 3; t++)
2485
{
2486
for (unsigned int u = 0; u < 3; u++)
2487
{
2488
dmats_old[t][u] = dmats[t][u];
2489
dmats[t][u] = 0.0;
2490
}// end loop over 'u'
2491
}// end loop over 't'
2492
2493
// Update dmats using an inner product.
2494
if (combinations_t[r][s] == 0)
2495
{
2496
for (unsigned int t = 0; t < 3; t++)
2497
{
2498
for (unsigned int u = 0; u < 3; u++)
2499
{
2500
for (unsigned int tu = 0; tu < 3; tu++)
2501
{
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'
2506
}
2507
2508
if (combinations_t[r][s] == 1)
2509
{
2510
for (unsigned int t = 0; t < 3; t++)
2511
{
2512
for (unsigned int u = 0; u < 3; u++)
2513
{
2514
for (unsigned int tu = 0; tu < 3; tu++)
2515
{
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'
2520
}
2521
2522
}// end loop over 's'
2523
for (unsigned int s = 0; s < 3; s++)
2524
{
2525
for (unsigned int t = 0; t < 3; t++)
2526
{
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'
2531
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'
2539
2540
// Transform derivatives back to physical element
2541
for (unsigned int r = 0; r < num_derivatives_g; r++)
2542
{
2543
for (unsigned int s = 0; s < num_derivatives_t; s++)
2544
{
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'
2550
2551
// Delete pointer to array of derivatives on FIAT element
2552
delete [] derivatives;
2553
2554
// Delete pointer to array of reference derivatives on physical element.
2555
delete [] derivatives_p;
2556
2557
// Delete pointer to array of combinations of derivatives and transform
2558
for (unsigned int r = 0; r < num_derivatives_t; r++)
2559
{
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++)
2564
{
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++)
2569
{
2570
delete [] transform[r];
2571
}// end loop over 'r'
2572
delete [] transform;
2573
break;
2574
}
2575
case 3:
2576
{
2577
2578
// Array of basisvalues
2579
double basisvalues[1] = {0.0};
2580
2581
// Declare helper variables
2582
2583
// Compute basisvalues
2584
basisvalues[0] = 1.0;
2585
2586
// Table(s) of coefficients
2587
static const double coefficients0[1] = \
2588
{1.0};
2589
2590
// Tables of derivatives of the polynomial base (transpose).
2591
static const double dmats0[1][1] = \
2592
{{0.0}};
2593
2594
static const double dmats1[1][1] = \
2595
{{0.0}};
2596
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++)
2601
{
2602
derivatives[r] = 0.0;
2603
}// end loop over 'r'
2604
2605
// Declare derivative matrix (of polynomial basis).
2606
double dmats[1][1] = \
2607
{{1.0}};
2608
2609
// Declare (auxiliary) derivative matrix (of polynomial basis).
2610
double dmats_old[1][1] = \
2611
{{1.0}};
2612
2613
// Loop possible derivatives.
2614
for (unsigned int r = 0; r < num_derivatives_t; r++)
2615
{
2616
// Resetting dmats values to compute next derivative.
2617
for (unsigned int t = 0; t < 1; t++)
2618
{
2619
for (unsigned int u = 0; u < 1; u++)
2620
{
2621
dmats[t][u] = 0.0;
2622
if (t == u)
2623
{
2624
dmats[t][u] = 1.0;
2625
}
2626
2627
}// end loop over 'u'
2628
}// end loop over 't'
2629
2630
// Looping derivative order to generate dmats.
2631
for (unsigned int s = 0; s < n; s++)
2632
{
2633
// Updating dmats_old with new values and resetting dmats.
2634
for (unsigned int t = 0; t < 1; t++)
2635
{
2636
for (unsigned int u = 0; u < 1; u++)
2637
{
2638
dmats_old[t][u] = dmats[t][u];
2639
dmats[t][u] = 0.0;
2640
}// end loop over 'u'
2641
}// end loop over 't'
2642
2643
// Update dmats using an inner product.
2644
if (combinations_t[r][s] == 0)
2645
{
2646
for (unsigned int t = 0; t < 1; t++)
2647
{
2648
for (unsigned int u = 0; u < 1; u++)
2649
{
2650
for (unsigned int tu = 0; tu < 1; tu++)
2651
{
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'
2656
}
2657
2658
if (combinations_t[r][s] == 1)
2659
{
2660
for (unsigned int t = 0; t < 1; t++)
2661
{
2662
for (unsigned int u = 0; u < 1; u++)
2663
{
2664
for (unsigned int tu = 0; tu < 1; tu++)
2665
{
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'
2670
}
2671
2672
}// end loop over 's'
2673
for (unsigned int s = 0; s < 1; s++)
2674
{
2675
for (unsigned int t = 0; t < 1; t++)
2676
{
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'
2681
2682
// Transform derivatives back to physical element
2683
for (unsigned int r = 0; r < num_derivatives_g; r++)
2684
{
2685
for (unsigned int s = 0; s < num_derivatives_t; s++)
2686
{
2687
values[2*num_derivatives_g + r] += transform[r][s]*derivatives[s];
2688
}// end loop over 's'
2689
}// end loop over 'r'
2690
2691
// Delete pointer to array of derivatives on FIAT element
2692
delete [] derivatives;
2693
2694
// Delete pointer to array of combinations of derivatives and transform
2695
for (unsigned int r = 0; r < num_derivatives_t; r++)
2696
{
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++)
2701
{
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++)
2706
{
2707
delete [] transform[r];
2708
}// end loop over 'r'
2709
delete [] transform;
2710
break;
2711
}
2712
}
2713
2714
}
2715
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,
2718
double* values,
2719
const double* x,
2720
const double* vertex_coordinates,
2721
int cell_orientation) const
2722
{
2723
// Compute number of derivatives.
2724
unsigned int num_derivatives_g = 1;
2725
for (unsigned int r = 0; r < n; r++)
2726
{
2727
num_derivatives_g *= 3;
2728
}// end loop over 'r'
2729
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++)
2733
{
2734
dof_values[r] = 0.0;
2735
}// end loop over 'r'
2736
2737
// Loop dofs and call evaluate_basis_derivatives.
2738
for (unsigned int r = 0; r < 4; r++)
2739
{
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++)
2742
{
2743
values[r*4*num_derivatives_g + s] = dof_values[s];
2744
}// end loop over 's'
2745
}// end loop over 'r'
2746
2747
// Delete pointer.
2748
delete [] dof_values;
2749
}
2750
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
2757
{
2758
// Declare variables for result of evaluation
2759
double vals[4];
2760
2761
// Declare variable for physical coordinates
2762
double y[3];
2763
2764
double result;
2765
// Compute Jacobian
2766
double J[6];
2767
compute_jacobian_triangle_3d(J, vertex_coordinates);
2768
2769
2770
// Compute Jacobian inverse and determinant
2771
double K[6];
2772
double detJ;
2773
compute_jacobian_inverse_triangle_3d(K, detJ, J);
2774
2775
2776
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)
2782
detJ *= -1;
2783
switch (i)
2784
{
2785
case 0:
2786
{
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]));
2792
return result;
2793
break;
2794
}
2795
case 1:
2796
{
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]));
2802
return result;
2803
break;
2804
}
2805
case 2:
2806
{
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]));
2812
return result;
2813
break;
2814
}
2815
case 3:
2816
{
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);
2821
return vals[3];
2822
break;
2823
}
2824
}
2825
2826
return 0.0;
2827
}
2828
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
2835
{
2836
// Declare variables for result of evaluation
2837
double vals[4];
2838
2839
// Declare variable for physical coordinates
2840
double y[3];
2841
2842
double result;
2843
// Compute Jacobian
2844
double J[6];
2845
compute_jacobian_triangle_3d(J, vertex_coordinates);
2846
2847
2848
// Compute Jacobian inverse and determinant
2849
double K[6];
2850
double detJ;
2851
compute_jacobian_inverse_triangle_3d(K, detJ, J);
2852
2853
2854
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)
2860
detJ *= -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]));
2866
values[0] = result;
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]));
2872
values[1] = result;
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]));
2878
values[2] = result;
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];
2884
}
2885
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
2892
{
2893
// Compute Jacobian
2894
double J[6];
2895
compute_jacobian_triangle_3d(J, vertex_coordinates);
2896
2897
2898
// Compute Jacobian inverse and determinant
2899
double K[6];
2900
double detJ;
2901
compute_jacobian_inverse_triangle_3d(K, detJ, J);
2902
2903
2904
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)
2910
detJ *= -1;
2911
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];
2923
}
2924
2925
/// Map coordinate xhat from reference cell to coordinate x in cell
2926
virtual void map_from_reference_cell(double* x,
2927
const double* xhat,
2928
const ufc::cell& c) const
2929
{
2930
std::cerr << "*** FFC warning: " << "map_from_reference_cell not yet implemented." << std::endl;
2931
}
2932
2933
/// Map from coordinate x in cell to coordinate xhat in reference cell
2934
virtual void map_to_reference_cell(double* xhat,
2935
const double* x,
2936
const ufc::cell& c) const
2937
{
2938
std::cerr << "*** FFC warning: " << "map_to_reference_cell not yet implemented." << std::endl;
2939
}
2940
2941
/// Return the number of sub elements (for a mixed element)
2942
virtual std::size_t num_sub_elements() const
2943
{
2944
return 2;
2945
}
2946
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
2949
{
2950
switch (i)
2951
{
2952
case 0:
2953
{
2954
return new projectionmanifold_finite_element_0();
2955
break;
2956
}
2957
case 1:
2958
{
2959
return new projectionmanifold_finite_element_1();
2960
break;
2961
}
2962
}
2963
2964
return 0;
2965
}
2966
2967
/// Create a new class instance
2968
virtual ufc::finite_element* create() const
2969
{
2970
return new projectionmanifold_finite_element_2();
2971
}
2972
2973
};
2974
2975
/// This class defines the interface for a local-to-global mapping of
2976
/// degrees of freedom (dofs).
2977
2978
class projectionmanifold_dofmap_0: public ufc::dofmap
2979
{
2980
public:
2981
2982
/// Constructor
2983
projectionmanifold_dofmap_0() : ufc::dofmap()
2984
{
2985
// Do nothing
2986
}
2987
2988
/// Destructor
2989
virtual ~projectionmanifold_dofmap_0()
2990
{
2991
// Do nothing
2992
}
2993
2994
/// Return a string identifying the dofmap
2995
virtual const char* signature() const
2996
{
2997
return "FFC dofmap for FiniteElement('Raviart-Thomas', Domain(Cell('triangle', 3), 'triangle_multiverse', 3, 2), 1, None)";
2998
}
2999
3000
/// Return true iff mesh entities of topological dimension d are needed
3001
virtual bool needs_mesh_entities(std::size_t d) const
3002
{
3003
switch (d)
3004
{
3005
case 0:
3006
{
3007
return false;
3008
break;
3009
}
3010
case 1:
3011
{
3012
return true;
3013
break;
3014
}
3015
case 2:
3016
{
3017
return false;
3018
break;
3019
}
3020
}
3021
3022
return false;
3023
}
3024
3025
/// Return the topological dimension of the associated cell shape
3026
virtual std::size_t topological_dimension() const
3027
{
3028
return 2;
3029
}
3030
3031
/// Return the geometric dimension of the associated cell shape
3032
virtual std::size_t geometric_dimension() const
3033
{
3034
return 3;
3035
}
3036
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
3040
{
3041
return num_global_entities[1];
3042
}
3043
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
3046
{
3047
return 3;
3048
}
3049
3050
/// Return the maximum dimension of the local finite element function space
3051
virtual std::size_t max_local_dimension() const
3052
{
3053
return 3;
3054
}
3055
3056
/// Return the number of dofs on each cell facet
3057
virtual std::size_t num_facet_dofs() const
3058
{
3059
return 1;
3060
}
3061
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
3064
{
3065
switch (d)
3066
{
3067
case 0:
3068
{
3069
return 0;
3070
break;
3071
}
3072
case 1:
3073
{
3074
return 1;
3075
break;
3076
}
3077
case 2:
3078
{
3079
return 0;
3080
break;
3081
}
3082
}
3083
3084
return 0;
3085
}
3086
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
3091
{
3092
dofs[0] = c.entity_indices[1][0];
3093
dofs[1] = c.entity_indices[1][1];
3094
dofs[2] = c.entity_indices[1][2];
3095
}
3096
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
3100
{
3101
switch (facet)
3102
{
3103
case 0:
3104
{
3105
dofs[0] = 0;
3106
break;
3107
}
3108
case 1:
3109
{
3110
dofs[0] = 1;
3111
break;
3112
}
3113
case 2:
3114
{
3115
dofs[0] = 2;
3116
break;
3117
}
3118
}
3119
3120
}
3121
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
3125
{
3126
if (d > 2)
3127
{
3128
std::cerr << "*** FFC warning: " << "d is larger than dimension (2)" << std::endl;
3129
}
3130
3131
switch (d)
3132
{
3133
case 0:
3134
{
3135
3136
break;
3137
}
3138
case 1:
3139
{
3140
if (i > 2)
3141
{
3142
std::cerr << "*** FFC warning: " << "i is larger than number of entities (2)" << std::endl;
3143
}
3144
3145
switch (i)
3146
{
3147
case 0:
3148
{
3149
dofs[0] = 0;
3150
break;
3151
}
3152
case 1:
3153
{
3154
dofs[0] = 1;
3155
break;
3156
}
3157
case 2:
3158
{
3159
dofs[0] = 2;
3160
break;
3161
}
3162
}
3163
3164
break;
3165
}
3166
case 2:
3167
{
3168
3169
break;
3170
}
3171
}
3172
3173
}
3174
3175
/// Tabulate the coordinates of all dofs on a cell
3176
virtual void tabulate_coordinates(double** dof_coordinates,
3177
const double* vertex_coordinates) const
3178
{
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];
3188
}
3189
3190
/// Return the number of sub dofmaps (for a mixed element)
3191
virtual std::size_t num_sub_dofmaps() const
3192
{
3193
return 0;
3194
}
3195
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
3198
{
3199
return 0;
3200
}
3201
3202
/// Create a new class instance
3203
virtual ufc::dofmap* create() const
3204
{
3205
return new projectionmanifold_dofmap_0();
3206
}
3207
3208
};
3209
3210
/// This class defines the interface for a local-to-global mapping of
3211
/// degrees of freedom (dofs).
3212
3213
class projectionmanifold_dofmap_1: public ufc::dofmap
3214
{
3215
public:
3216
3217
/// Constructor
3218
projectionmanifold_dofmap_1() : ufc::dofmap()
3219
{
3220
// Do nothing
3221
}
3222
3223
/// Destructor
3224
virtual ~projectionmanifold_dofmap_1()
3225
{
3226
// Do nothing
3227
}
3228
3229
/// Return a string identifying the dofmap
3230
virtual const char* signature() const
3231
{
3232
return "FFC dofmap for FiniteElement('Discontinuous Lagrange', Domain(Cell('triangle', 3), 'triangle_multiverse', 3, 2), 0, None)";
3233
}
3234
3235
/// Return true iff mesh entities of topological dimension d are needed
3236
virtual bool needs_mesh_entities(std::size_t d) const
3237
{
3238
switch (d)
3239
{
3240
case 0:
3241
{
3242
return false;
3243
break;
3244
}
3245
case 1:
3246
{
3247
return false;
3248
break;
3249
}
3250
case 2:
3251
{
3252
return true;
3253
break;
3254
}
3255
}
3256
3257
return false;
3258
}
3259
3260
/// Return the topological dimension of the associated cell shape
3261
virtual std::size_t topological_dimension() const
3262
{
3263
return 2;
3264
}
3265
3266
/// Return the geometric dimension of the associated cell shape
3267
virtual std::size_t geometric_dimension() const
3268
{
3269
return 3;
3270
}
3271
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
3275
{
3276
return num_global_entities[2];
3277
}
3278
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
3281
{
3282
return 1;
3283
}
3284
3285
/// Return the maximum dimension of the local finite element function space
3286
virtual std::size_t max_local_dimension() const
3287
{
3288
return 1;
3289
}
3290
3291
/// Return the number of dofs on each cell facet
3292
virtual std::size_t num_facet_dofs() const
3293
{
3294
return 0;
3295
}
3296
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
3299
{
3300
switch (d)
3301
{
3302
case 0:
3303
{
3304
return 0;
3305
break;
3306
}
3307
case 1:
3308
{
3309
return 0;
3310
break;
3311
}
3312
case 2:
3313
{
3314
return 1;
3315
break;
3316
}
3317
}
3318
3319
return 0;
3320
}
3321
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
3326
{
3327
dofs[0] = c.entity_indices[2][0];
3328
}
3329
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
3333
{
3334
switch (facet)
3335
{
3336
case 0:
3337
{
3338
3339
break;
3340
}
3341
case 1:
3342
{
3343
3344
break;
3345
}
3346
case 2:
3347
{
3348
3349
break;
3350
}
3351
}
3352
3353
}
3354
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
3358
{
3359
if (d > 2)
3360
{
3361
std::cerr << "*** FFC warning: " << "d is larger than dimension (2)" << std::endl;
3362
}
3363
3364
switch (d)
3365
{
3366
case 0:
3367
{
3368
3369
break;
3370
}
3371
case 1:
3372
{
3373
3374
break;
3375
}
3376
case 2:
3377
{
3378
if (i > 0)
3379
{
3380
std::cerr << "*** FFC warning: " << "i is larger than number of entities (0)" << std::endl;
3381
}
3382
3383
dofs[0] = 0;
3384
break;
3385
}
3386
}
3387
3388
}
3389
3390
/// Tabulate the coordinates of all dofs on a cell
3391
virtual void tabulate_coordinates(double** dof_coordinates,
3392
const double* vertex_coordinates) const
3393
{
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];
3397
}
3398
3399
/// Return the number of sub dofmaps (for a mixed element)
3400
virtual std::size_t num_sub_dofmaps() const
3401
{
3402
return 0;
3403
}
3404
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
3407
{
3408
return 0;
3409
}
3410
3411
/// Create a new class instance
3412
virtual ufc::dofmap* create() const
3413
{
3414
return new projectionmanifold_dofmap_1();
3415
}
3416
3417
};
3418
3419
/// This class defines the interface for a local-to-global mapping of
3420
/// degrees of freedom (dofs).
3421
3422
class projectionmanifold_dofmap_2: public ufc::dofmap
3423
{
3424
public:
3425
3426
/// Constructor
3427
projectionmanifold_dofmap_2() : ufc::dofmap()
3428
{
3429
// Do nothing
3430
}
3431
3432
/// Destructor
3433
virtual ~projectionmanifold_dofmap_2()
3434
{
3435
// Do nothing
3436
}
3437
3438
/// Return a string identifying the dofmap
3439
virtual const char* signature() const
3440
{
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,) })";
3442
}
3443
3444
/// Return true iff mesh entities of topological dimension d are needed
3445
virtual bool needs_mesh_entities(std::size_t d) const
3446
{
3447
switch (d)
3448
{
3449
case 0:
3450
{
3451
return false;
3452
break;
3453
}
3454
case 1:
3455
{
3456
return true;
3457
break;
3458
}
3459
case 2:
3460
{
3461
return true;
3462
break;
3463
}
3464
}
3465
3466
return false;
3467
}
3468
3469
/// Return the topological dimension of the associated cell shape
3470
virtual std::size_t topological_dimension() const
3471
{
3472
return 2;
3473
}
3474
3475
/// Return the geometric dimension of the associated cell shape
3476
virtual std::size_t geometric_dimension() const
3477
{
3478
return 3;
3479
}
3480
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
3484
{
3485
return num_global_entities[1] + num_global_entities[2];
3486
}
3487
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
3490
{
3491
return 4;
3492
}
3493
3494
/// Return the maximum dimension of the local finite element function space
3495
virtual std::size_t max_local_dimension() const
3496
{
3497
return 4;
3498
}
3499
3500
/// Return the number of dofs on each cell facet
3501
virtual std::size_t num_facet_dofs() const
3502
{
3503
return 1;
3504
}
3505
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
3508
{
3509
switch (d)
3510
{
3511
case 0:
3512
{
3513
return 0;
3514
break;
3515
}
3516
case 1:
3517
{
3518
return 1;
3519
break;
3520
}
3521
case 2:
3522
{
3523
return 1;
3524
break;
3525
}
3526
}
3527
3528
return 0;
3529
}
3530
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
3535
{
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];
3543
}
3544
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
3548
{
3549
switch (facet)
3550
{
3551
case 0:
3552
{
3553
dofs[0] = 0;
3554
break;
3555
}
3556
case 1:
3557
{
3558
dofs[0] = 1;
3559
break;
3560
}
3561
case 2:
3562
{
3563
dofs[0] = 2;
3564
break;
3565
}
3566
}
3567
3568
}
3569
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
3573
{
3574
if (d > 2)
3575
{
3576
std::cerr << "*** FFC warning: " << "d is larger than dimension (2)" << std::endl;
3577
}
3578
3579
switch (d)
3580
{
3581
case 0:
3582
{
3583
3584
break;
3585
}
3586
case 1:
3587
{
3588
if (i > 2)
3589
{
3590
std::cerr << "*** FFC warning: " << "i is larger than number of entities (2)" << std::endl;
3591
}
3592
3593
switch (i)
3594
{
3595
case 0:
3596
{
3597
dofs[0] = 0;
3598
break;
3599
}
3600
case 1:
3601
{
3602
dofs[0] = 1;
3603
break;
3604
}
3605
case 2:
3606
{
3607
dofs[0] = 2;
3608
break;
3609
}
3610
}
3611
3612
break;
3613
}
3614
case 2:
3615
{
3616
if (i > 0)
3617
{
3618
std::cerr << "*** FFC warning: " << "i is larger than number of entities (0)" << std::endl;
3619
}
3620
3621
dofs[0] = 3;
3622
break;
3623
}
3624
}
3625
3626
}
3627
3628
/// Tabulate the coordinates of all dofs on a cell
3629
virtual void tabulate_coordinates(double** dof_coordinates,
3630
const double* vertex_coordinates) const
3631
{
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];
3644
}
3645
3646
/// Return the number of sub dofmaps (for a mixed element)
3647
virtual std::size_t num_sub_dofmaps() const
3648
{
3649
return 2;
3650
}
3651
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
3654
{
3655
switch (i)
3656
{
3657
case 0:
3658
{
3659
return new projectionmanifold_dofmap_0();
3660
break;
3661
}
3662
case 1:
3663
{
3664
return new projectionmanifold_dofmap_1();
3665
break;
3666
}
3667
}
3668
3669
return 0;
3670
}
3671
3672
/// Create a new class instance
3673
virtual ufc::dofmap* create() const
3674
{
3675
return new projectionmanifold_dofmap_2();
3676
}
3677
3678
};
3679
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.
3683
3684
class projectionmanifold_cell_integral_0_otherwise: public ufc::cell_integral
3685
{
3686
public:
3687
3688
/// Constructor
3689
projectionmanifold_cell_integral_0_otherwise() : ufc::cell_integral()
3690
{
3691
// Do nothing
3692
}
3693
3694
/// Destructor
3695
virtual ~projectionmanifold_cell_integral_0_otherwise()
3696
{
3697
// Do nothing
3698
}
3699
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
3705
{
3706
// Compute Jacobian
3707
double J[6];
3708
compute_jacobian_triangle_3d(J, vertex_coordinates);
3709
3710
// Compute Jacobian inverse and determinant
3711
double K[6];
3712
double detJ;
3713
compute_jacobian_inverse_triangle_3d(K, detJ, J);
3714
3715
// Check orientation
3716
if (cell_orientation == -1)
3717
throw std::runtime_error("cell orientation must be defined (not -1)");
3718
// (If cell_orientation == 1 = down, multiply det(J) by -1)
3719
else if (cell_orientation == 1)
3720
detJ *= -1;
3721
3722
// Set scale factor
3723
const double det = std::abs(detJ);
3724
3725
// Cell volume of triangle in 3D
3726
3727
// Compute circumradius of triangle in 3D
3728
3729
3730
// Facet area
3731
3732
// Array of quadrature weights.
3733
static const double W3[3] = {0.16666667, 0.16666667, 0.16666667};
3734
// Quadrature points on the UFC reference element: (0.16666667, 0.16666667), (0.16666667, 0.66666667), (0.66666667, 0.16666667)
3735
3736
// Value of basis functions at quadrature points.
3737
static const double FE0_C0[3][3] = \
3738
{{0.16666667, 0.83333333, 0.16666667},
3739
{0.16666667, 0.83333333, 0.16666667},
3740
{0.66666667, 0.33333333, 0.66666667}};
3741
3742
// Array of non-zero columns
3743
static const unsigned int nzc0[3] = {0, 1, 2};
3744
3745
static const double FE0_C0_D10[3][3] = \
3746
{{1.0, -1.0, 1.0},
3747
{1.0, -1.0, 1.0},
3748
{1.0, -1.0, 1.0}};
3749
3750
// Array of non-zero columns
3751
static const unsigned int nzc1[3] = {0, 1, 2};
3752
3753
static const double FE0_C1[3][3] = \
3754
{{0.16666667, -0.16666667, -0.83333333},
3755
{0.66666667, -0.66666667, -0.33333333},
3756
{0.16666667, -0.16666667, -0.83333333}};
3757
3758
// Array of non-zero columns
3759
static const unsigned int nzc2[3] = {0, 1, 2};
3760
3761
// Reset values in the element tensor.
3762
for (unsigned int r = 0; r < 16; r++)
3763
{
3764
A[r] = 0.0;
3765
}// end loop over 'r'
3766
// Number of operations to compute geometry constants: 37.
3767
double G[4];
3768
G[0] = det*(J[0]*K[0] + J[1]*K[3] + J[2]*K[1] + J[3]*K[4] + J[4]*K[2] + J[5]*K[5])/detJ;
3769
G[1] = det*(J[0]*J[0] + J[2]*J[2] + J[4]*J[4])/(detJ*detJ);
3770
G[2] = det*(J[0]*J[1] + J[2]*J[3] + J[4]*J[5])/(detJ*detJ);
3771
G[3] = det*(J[1]*J[1] + J[3]*J[3] + J[5]*J[5])/(detJ*detJ);
3772
3773
// Compute element tensor using UFL quadrature representation
3774
// Optimisations: ('eliminate zeros', True), ('ignore ones', True), ('ignore zero tables', True), ('optimisation', 'simplify_expressions'), ('remove zero terms', True)
3775
3776
// Loop quadrature points for integral.
3777
// Number of operations to compute element tensor for following IP loop = 372
3778
for (unsigned int ip = 0; ip < 3; ip++)
3779
{
3780
3781
// Number of operations to compute ip constants: 4
3782
double I[4];
3783
// Number of operations: 1
3784
I[0] = G[0]*W3[ip];
3785
3786
// Number of operations: 1
3787
I[1] = G[1]*W3[ip];
3788
3789
// Number of operations: 1
3790
I[2] = G[2]*W3[ip];
3791
3792
// Number of operations: 1
3793
I[3] = G[3]*W3[ip];
3794
3795
3796
// Number of operations for primary indices: 108
3797
for (unsigned int j = 0; j < 3; j++)
3798
{
3799
for (unsigned int k = 0; k < 3; k++)
3800
{
3801
// Number of operations to compute entry: 3
3802
A[nzc0[j]*4 + nzc0[k]] += FE0_C0[ip][j]*FE0_C0[ip][k]*I[1];
3803
// Number of operations to compute entry: 3
3804
A[nzc0[j]*4 + nzc2[k]] += FE0_C0[ip][j]*FE0_C1[ip][k]*I[2];
3805
// Number of operations to compute entry: 3
3806
A[nzc2[j]*4 + nzc0[k]] += FE0_C0[ip][k]*FE0_C1[ip][j]*I[2];
3807
// Number of operations to compute entry: 3
3808
A[nzc2[j]*4 + nzc2[k]] += FE0_C1[ip][j]*FE0_C1[ip][k]*I[3];
3809
}// end loop over 'k'
3810
}// end loop over 'j'
3811
3812
// Number of operations for primary indices: 6
3813
for (unsigned int j = 0; j < 3; j++)
3814
{
3815
// Number of operations to compute entry: 2
3816
A[nzc1[j]*4 + 3] += FE0_C0_D10[ip][j]*I[0];
3817
}// end loop over 'j'
3818
3819
// Number of operations for primary indices: 6
3820
for (unsigned int k = 0; k < 3; k++)
3821
{
3822
// Number of operations to compute entry: 2
3823
A[3*4 + nzc1[k]] += FE0_C0_D10[ip][k]*I[0];
3824
}// end loop over 'k'
3825
}// end loop over 'ip'
3826
}
3827
3828
};
3829
3830
/// This class defines the interface for the assembly of the global
3831
/// tensor corresponding to a form with r + n arguments, that is, a
3832
/// mapping
3833
///
3834
/// a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
3835
///
3836
/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
3837
/// global tensor A is defined by
3838
///
3839
/// A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
3840
///
3841
/// where each argument Vj represents the application to the
3842
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
3843
/// fixed functions (coefficients).
3844
3845
class projectionmanifold_form_0: public ufc::form
3846
{
3847
public:
3848
3849
/// Constructor
3850
projectionmanifold_form_0() : ufc::form()
3851
{
3852
// Do nothing
3853
}
3854
3855
/// Destructor
3856
virtual ~projectionmanifold_form_0()
3857
{
3858
// Do nothing
3859
}
3860
3861
/// Return a string identifying the form
3862
virtual const char* signature() const
3863
{
3864
return "e78a2278f8cab860e9ac0865b6fee052698344942c6cbff9e41a848e6f2a49353e4cd71f65096c2facc2f4034f037d7eebb82d42c0e9ecd9a8ef113467c90a0f";
3865
}
3866
3867
/// Return the rank of the global tensor (r)
3868
virtual std::size_t rank() const
3869
{
3870
return 2;
3871
}
3872
3873
/// Return the number of coefficients (n)
3874
virtual std::size_t num_coefficients() const
3875
{
3876
return 0;
3877
}
3878
3879
/// Return the number of cell domains
3880
virtual std::size_t num_cell_domains() const
3881
{
3882
return 0;
3883
}
3884
3885
/// Return the number of exterior facet domains
3886
virtual std::size_t num_exterior_facet_domains() const
3887
{
3888
return 0;
3889
}
3890
3891
/// Return the number of interior facet domains
3892
virtual std::size_t num_interior_facet_domains() const
3893
{
3894
return 0;
3895
}
3896
3897
/// Return the number of point domains
3898
virtual std::size_t num_point_domains() const
3899
{
3900
return 0;
3901
}
3902
3903
/// Return whether the form has any cell integrals
3904
virtual bool has_cell_integrals() const
3905
{
3906
return true;
3907
}
3908
3909
/// Return whether the form has any exterior facet integrals
3910
virtual bool has_exterior_facet_integrals() const
3911
{
3912
return false;
3913
}
3914
3915
/// Return whether the form has any interior facet integrals
3916
virtual bool has_interior_facet_integrals() const
3917
{
3918
return false;
3919
}
3920
3921
/// Return whether the form has any point integrals
3922
virtual bool has_point_integrals() const
3923
{
3924
return false;
3925
}
3926
3927
/// Create a new finite element for argument function i
3928
virtual ufc::finite_element* create_finite_element(std::size_t i) const
3929
{
3930
switch (i)
3931
{
3932
case 0:
3933
{
3934
return new projectionmanifold_finite_element_2();
3935
break;
3936
}
3937
case 1:
3938
{
3939
return new projectionmanifold_finite_element_2();
3940
break;
3941
}
3942
}
3943
3944
return 0;
3945
}
3946
3947
/// Create a new dofmap for argument function i
3948
virtual ufc::dofmap* create_dofmap(std::size_t i) const
3949
{
3950
switch (i)
3951
{
3952
case 0:
3953
{
3954
return new projectionmanifold_dofmap_2();
3955
break;
3956
}
3957
case 1:
3958
{
3959
return new projectionmanifold_dofmap_2();
3960
break;
3961
}
3962
}
3963
3964
return 0;
3965
}
3966
3967
/// Create a new cell integral on sub domain i
3968
virtual ufc::cell_integral* create_cell_integral(std::size_t i) const
3969
{
3970
return 0;
3971
}
3972
3973
/// Create a new exterior facet integral on sub domain i
3974
virtual ufc::exterior_facet_integral* create_exterior_facet_integral(std::size_t i) const
3975
{
3976
return 0;
3977
}
3978
3979
/// Create a new interior facet integral on sub domain i
3980
virtual ufc::interior_facet_integral* create_interior_facet_integral(std::size_t i) const
3981
{
3982
return 0;
3983
}
3984
3985
/// Create a new point integral on sub domain i
3986
virtual ufc::point_integral* create_point_integral(std::size_t i) const
3987
{
3988
return 0;
3989
}
3990
3991
/// Create a new cell integral on everywhere else
3992
virtual ufc::cell_integral* create_default_cell_integral() const
3993
{
3994
return new projectionmanifold_cell_integral_0_otherwise();
3995
}
3996
3997
/// Create a new exterior facet integral on everywhere else
3998
virtual ufc::exterior_facet_integral* create_default_exterior_facet_integral() const
3999
{
4000
return 0;
4001
}
4002
4003
/// Create a new interior facet integral on everywhere else
4004
virtual ufc::interior_facet_integral* create_default_interior_facet_integral() const
4005
{
4006
return 0;
4007
}
4008
4009
/// Create a new point integral on everywhere else
4010
virtual ufc::point_integral* create_default_point_integral() const
4011
{
4012
return 0;
4013
}
4014
4015
};
4016
4017
#endif
Loggerhead is a web-based interface for Breezy
Version: 2.0.1