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