~njansson/dolfin/hpc

return "Mixed finite element: [Lagrange finite element of degree 1 on a tetrahedron, Lagrange finite element of degree 1 on a tetrahedron, Lagrange finite element of degree 1 on a tetrahedron]";

1623

}

1624

1625

/// Return the cell shape

1626

virtual ufc::shape cell_shape() const

1627

{

1628

return ufc::tetrahedron;

1629

}

1630

1631

/// Return the dimension of the finite element function space

1632

virtual unsigned int space_dimension() const

1633

{

1634

return 12;

1635

}

1636

1637

/// Return the rank of the value space

1638

virtual unsigned int value_rank() const

1639

{

1640

return 1;

1641

}

1642

1643

/// Return the dimension of the value space for axis i

1644

virtual unsigned int value_dimension(unsigned int i) const

1645

{

1646

return 3;

1647

}

1648

1649

/// Evaluate basis function i at given point in cell

1650

virtual void evaluate_basis(unsigned int i,

1651

double* values,

1652

const double* coordinates,

1653

const ufc::cell& c) const

1654

{

1655

// Extract vertex coordinates

1656

const double * const * element_coordinates = c.coordinates;

1657

1658

// Compute Jacobian of affine map from reference cell

1659

const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];

1660

const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];

1661

const double J_02 = element_coordinates[3][0] - element_coordinates[0][0];

1662

const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];

1663

const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];

1664

const double J_12 = element_coordinates[3][1] - element_coordinates[0][1];

1665

const double J_20 = element_coordinates[1][2] - element_coordinates[0][2];

1666

const double J_21 = element_coordinates[2][2] - element_coordinates[0][2];

1667

const double J_22 = element_coordinates[3][2] - element_coordinates[0][2];

1668

1669

// Compute sub determinants

1670

const double d00 = J_11*J_22 - J_12*J_21;

1671

const double d01 = J_12*J_20 - J_10*J_22;

1672

const double d02 = J_10*J_21 - J_11*J_20;

1673

1674

const double d10 = J_02*J_21 - J_01*J_22;

1675

const double d11 = J_00*J_22 - J_02*J_20;

1676

const double d12 = J_01*J_20 - J_00*J_21;

1677

1678

const double d20 = J_01*J_12 - J_02*J_11;

1679

const double d21 = J_02*J_10 - J_00*J_12;

1680

const double d22 = J_00*J_11 - J_01*J_10;

1681

1682

// Compute determinant of Jacobian

1683

double detJ = J_00*d00 + J_10*d10 + J_20*d20;

1684

1685

// Compute inverse of Jacobian

1686

1687

// Compute constants

1688

const double C0 = d00*(element_coordinates[0][0] - element_coordinates[2][0] - element_coordinates[3][0]) \

1689

+ d10*(element_coordinates[0][1] - element_coordinates[2][1] - element_coordinates[3][1]) \

1690

+ d20*(element_coordinates[0][2] - element_coordinates[2][2] - element_coordinates[3][2]);

1691

1692

const double C1 = d01*(element_coordinates[0][0] - element_coordinates[1][0] - element_coordinates[3][0]) \

1693

+ d11*(element_coordinates[0][1] - element_coordinates[1][1] - element_coordinates[3][1]) \

1694

+ d21*(element_coordinates[0][2] - element_coordinates[1][2] - element_coordinates[3][2]);

1695

1696

const double C2 = d02*(element_coordinates[0][0] - element_coordinates[1][0] - element_coordinates[2][0]) \

1697

+ d12*(element_coordinates[0][1] - element_coordinates[1][1] - element_coordinates[2][1]) \

1698

+ d22*(element_coordinates[0][2] - element_coordinates[1][2] - element_coordinates[2][2]);

1699

1700

// Get coordinates and map to the UFC reference element

1701

double x = (C0 + d00*coordinates[0] + d10*coordinates[1] + d20*coordinates[2]) / detJ;

1702

double y = (C1 + d01*coordinates[0] + d11*coordinates[1] + d21*coordinates[2]) / detJ;

1703

double z = (C2 + d02*coordinates[0] + d12*coordinates[1] + d22*coordinates[2]) / detJ;

1704

1705

// Map coordinates to the reference cube

1706

if (std::abs(y + z - 1.0) < 1e-14)

1707

x = 1.0;

1708

else

1709

x = -2.0 * x/(y + z - 1.0) - 1.0;

1710

if (std::abs(z - 1.0) < 1e-14)

1711

y = -1.0;

1712

else

1713

y = 2.0 * y/(1.0 - z) - 1.0;

1714

z = 2.0 * z - 1.0;

1715

1716

// Reset values

1717

values[0] = 0;

1718

values[1] = 0;

1719

values[2] = 0;

1720

1721

if (0 <= i && i <= 3)

1722

{

1723

// Map degree of freedom to element degree of freedom

1724

const unsigned int dof = i;

1725

1726

// Generate scalings

1727

const double scalings_y_0 = 1;

1728

const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);

1729

const double scalings_z_0 = 1;

1730

const double scalings_z_1 = scalings_z_0*(0.5 - 0.5*z);

1731

1732

// Compute psitilde_a

1733

const double psitilde_a_0 = 1;

1734

const double psitilde_a_1 = x;

1735

1736

// Compute psitilde_bs

1737

const double psitilde_bs_0_0 = 1;

1738

const double psitilde_bs_0_1 = 1.5*y + 0.5;

1739

const double psitilde_bs_1_0 = 1;

1740

1741

// Compute psitilde_cs

1742

const double psitilde_cs_00_0 = 1;

1743

const double psitilde_cs_00_1 = 2*z + 1;

1744

const double psitilde_cs_01_0 = 1;

1745

const double psitilde_cs_10_0 = 1;

1746

1747

// Compute basisvalues

1748

const double basisvalue0 = 0.866025403784439*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_0;

1749

const double basisvalue1 = 2.73861278752583*psitilde_a_1*scalings_y_1*psitilde_bs_1_0*scalings_z_1*psitilde_cs_10_0;

1750

const double basisvalue2 = 1.58113883008419*psitilde_a_0*scalings_y_0*psitilde_bs_0_1*scalings_z_1*psitilde_cs_01_0;

1751

const double basisvalue3 = 1.11803398874989*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_1;

1752

1753

// Table(s) of coefficients

1754

const static double coefficients0[4][4] = \

1755

{{0.288675134594813, -0.182574185835055, -0.105409255338946, -0.074535599249993},

1756

{0.288675134594813, 0.182574185835055, -0.105409255338946, -0.074535599249993},

1757

{0.288675134594813, 0, 0.210818510677892, -0.074535599249993},

1758

{0.288675134594813, 0, 0, 0.223606797749979}};

1759

1760

// Extract relevant coefficients

1761

const double coeff0_0 = coefficients0[dof][0];

1762

const double coeff0_1 = coefficients0[dof][1];

1763

const double coeff0_2 = coefficients0[dof][2];

1764

const double coeff0_3 = coefficients0[dof][3];

1765

1766

// Compute value(s)

1767

values[0] = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2 + coeff0_3*basisvalue3;

1768

}

1769

1770

if (4 <= i && i <= 7)

1771

{

1772

// Map degree of freedom to element degree of freedom

1773

const unsigned int dof = i - 4;

1774

1775

// Generate scalings

1776

const double scalings_y_0 = 1;

1777

const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);

1778

const double scalings_z_0 = 1;

1779

const double scalings_z_1 = scalings_z_0*(0.5 - 0.5*z);

1780

1781

// Compute psitilde_a

1782

const double psitilde_a_0 = 1;

1783

const double psitilde_a_1 = x;

1784

1785

// Compute psitilde_bs

1786

const double psitilde_bs_0_0 = 1;

1787

const double psitilde_bs_0_1 = 1.5*y + 0.5;

1788

const double psitilde_bs_1_0 = 1;

1789

1790

// Compute psitilde_cs

1791

const double psitilde_cs_00_0 = 1;

1792

const double psitilde_cs_00_1 = 2*z + 1;

1793

const double psitilde_cs_01_0 = 1;

1794

const double psitilde_cs_10_0 = 1;

1795

1796

// Compute basisvalues

1797

const double basisvalue0 = 0.866025403784439*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_0;

1798

const double basisvalue1 = 2.73861278752583*psitilde_a_1*scalings_y_1*psitilde_bs_1_0*scalings_z_1*psitilde_cs_10_0;

1799

const double basisvalue2 = 1.58113883008419*psitilde_a_0*scalings_y_0*psitilde_bs_0_1*scalings_z_1*psitilde_cs_01_0;

1800

const double basisvalue3 = 1.11803398874989*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_1;

1801

1802

// Table(s) of coefficients

1803

const static double coefficients0[4][4] = \

1804

{{0.288675134594813, -0.182574185835055, -0.105409255338946, -0.074535599249993},

1805

{0.288675134594813, 0.182574185835055, -0.105409255338946, -0.074535599249993},

1806

{0.288675134594813, 0, 0.210818510677892, -0.074535599249993},

1807

{0.288675134594813, 0, 0, 0.223606797749979}};

1808

1809

// Extract relevant coefficients

1810

const double coeff0_0 = coefficients0[dof][0];

1811

const double coeff0_1 = coefficients0[dof][1];

1812

const double coeff0_2 = coefficients0[dof][2];

1813

const double coeff0_3 = coefficients0[dof][3];

1814

1815

// Compute value(s)

1816

values[1] = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2 + coeff0_3*basisvalue3;

1817

}

1818

1819

if (8 <= i && i <= 11)

1820

{

1821

// Map degree of freedom to element degree of freedom

1822

const unsigned int dof = i - 8;

1823

1824

// Generate scalings

1825

const double scalings_y_0 = 1;

1826

const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);

1827

const double scalings_z_0 = 1;

1828

const double scalings_z_1 = scalings_z_0*(0.5 - 0.5*z);

1829

1830

// Compute psitilde_a

1831

const double psitilde_a_0 = 1;

1832

const double psitilde_a_1 = x;

1833

1834

// Compute psitilde_bs

1835

const double psitilde_bs_0_0 = 1;

1836

const double psitilde_bs_0_1 = 1.5*y + 0.5;

1837

const double psitilde_bs_1_0 = 1;

1838

1839

// Compute psitilde_cs

1840

const double psitilde_cs_00_0 = 1;

1841

const double psitilde_cs_00_1 = 2*z + 1;

1842

const double psitilde_cs_01_0 = 1;

1843

const double psitilde_cs_10_0 = 1;

1844

1845

// Compute basisvalues

1846

const double basisvalue0 = 0.866025403784439*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_0;

1847

const double basisvalue1 = 2.73861278752583*psitilde_a_1*scalings_y_1*psitilde_bs_1_0*scalings_z_1*psitilde_cs_10_0;

1848

const double basisvalue2 = 1.58113883008419*psitilde_a_0*scalings_y_0*psitilde_bs_0_1*scalings_z_1*psitilde_cs_01_0;

1849

const double basisvalue3 = 1.11803398874989*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_1;

1850

1851

// Table(s) of coefficients

1852

const static double coefficients0[4][4] = \

1853

{{0.288675134594813, -0.182574185835055, -0.105409255338946, -0.074535599249993},

1854

{0.288675134594813, 0.182574185835055, -0.105409255338946, -0.074535599249993},

1855

{0.288675134594813, 0, 0.210818510677892, -0.074535599249993},

1856

{0.288675134594813, 0, 0, 0.223606797749979}};

1857

1858

// Extract relevant coefficients

1859

const double coeff0_0 = coefficients0[dof][0];

1860

const double coeff0_1 = coefficients0[dof][1];

1861

const double coeff0_2 = coefficients0[dof][2];

1862

const double coeff0_3 = coefficients0[dof][3];

1863

1864

// Compute value(s)

1865

values[2] = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2 + coeff0_3*basisvalue3;

1866

}

1867

1868

}

1869

1870

/// Evaluate all basis functions at given point in cell

1871

virtual void evaluate_basis_all(double* values,

1872

const double* coordinates,

1873

const ufc::cell& c) const

1874

{

1875

throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");

1876

}

1877

1878

/// Evaluate order n derivatives of basis function i at given point in cell

1879

virtual void evaluate_basis_derivatives(unsigned int i,

1880

unsigned int n,

1881

double* values,

1882

const double* coordinates,

1883

const ufc::cell& c) const

1884

{

1885

// Extract vertex coordinates

1886

const double * const * element_coordinates = c.coordinates;

1887

1888

// Compute Jacobian of affine map from reference cell

1889

const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];

1890

const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];

1891

const double J_02 = element_coordinates[3][0] - element_coordinates[0][0];

1892

const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];

1893

const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];

1894

const double J_12 = element_coordinates[3][1] - element_coordinates[0][1];

1895

const double J_20 = element_coordinates[1][2] - element_coordinates[0][2];

1896

const double J_21 = element_coordinates[2][2] - element_coordinates[0][2];

1897

const double J_22 = element_coordinates[3][2] - element_coordinates[0][2];

1898

1899

// Compute sub determinants

1900

const double d00 = J_11*J_22 - J_12*J_21;

1901

const double d01 = J_12*J_20 - J_10*J_22;

1902

const double d02 = J_10*J_21 - J_11*J_20;

1903

1904

const double d10 = J_02*J_21 - J_01*J_22;

1905

const double d11 = J_00*J_22 - J_02*J_20;

1906

const double d12 = J_01*J_20 - J_00*J_21;

1907

1908

const double d20 = J_01*J_12 - J_02*J_11;

1909

const double d21 = J_02*J_10 - J_00*J_12;

1910

const double d22 = J_00*J_11 - J_01*J_10;

1911

1912

// Compute determinant of Jacobian

1913

double detJ = J_00*d00 + J_10*d10 + J_20*d20;

1914

1915

// Compute inverse of Jacobian

1916

1917

// Compute constants

1918

const double C0 = d00*(element_coordinates[0][0] - element_coordinates[2][0] - element_coordinates[3][0]) \

1919

+ d10*(element_coordinates[0][1] - element_coordinates[2][1] - element_coordinates[3][1]) \

1920

+ d20*(element_coordinates[0][2] - element_coordinates[2][2] - element_coordinates[3][2]);

1921

1922

const double C1 = d01*(element_coordinates[0][0] - element_coordinates[1][0] - element_coordinates[3][0]) \

1923

+ d11*(element_coordinates[0][1] - element_coordinates[1][1] - element_coordinates[3][1]) \

1924

+ d21*(element_coordinates[0][2] - element_coordinates[1][2] - element_coordinates[3][2]);

1925

1926

const double C2 = d02*(element_coordinates[0][0] - element_coordinates[1][0] - element_coordinates[2][0]) \

1927

+ d12*(element_coordinates[0][1] - element_coordinates[1][1] - element_coordinates[2][1]) \

1928

+ d22*(element_coordinates[0][2] - element_coordinates[1][2] - element_coordinates[2][2]);

1929

1930

// Get coordinates and map to the UFC reference element

1931

double x = (C0 + d00*coordinates[0] + d10*coordinates[1] + d20*coordinates[2]) / detJ;

1932

double y = (C1 + d01*coordinates[0] + d11*coordinates[1] + d21*coordinates[2]) / detJ;

1933

double z = (C2 + d02*coordinates[0] + d12*coordinates[1] + d22*coordinates[2]) / detJ;

1934

1935

// Map coordinates to the reference cube

1936

if (std::abs(y + z - 1.0) < 1e-14)

1937

x = 1.0;

1938

else

1939

x = -2.0 * x/(y + z - 1.0) - 1.0;

1940

if (std::abs(z - 1.0) < 1e-14)

1941

y = -1.0;

1942

else

1943

y = 2.0 * y/(1.0 - z) - 1.0;

1944

z = 2.0 * z - 1.0;

1945

1946

// Compute number of derivatives

1947

unsigned int num_derivatives = 1;

1948

1949

for (unsigned int j = 0; j < n; j++)

1950

num_derivatives *= 3;

1951

1952

1953

// Declare pointer to two dimensional array that holds combinations of derivatives and initialise

1954

unsigned int **combinations = new unsigned int *[num_derivatives];

1955

1956

for (unsigned int j = 0; j < num_derivatives; j++)

1957

{

1958

combinations[j] = new unsigned int [n];

1959

for (unsigned int k = 0; k < n; k++)

1960

combinations[j][k] = 0;

1961

}

1962

1963

// Generate combinations of derivatives

1964

for (unsigned int row = 1; row < num_derivatives; row++)

1965

{

1966

for (unsigned int num = 0; num < row; num++)

1967

{

1968

for (unsigned int col = n-1; col+1 > 0; col--)

1969

{

1970

if (combinations[row][col] + 1 > 2)

1971

combinations[row][col] = 0;

1972

else

1973

{

1974

combinations[row][col] += 1;

1975

break;

1976

}

1977

}

1978

}

1979

}

1980

1981

// Compute inverse of Jacobian

1982

const double Jinv[3][3] ={{d00 / detJ, d10 / detJ, d20 / detJ}, {d01 / detJ, d11 / detJ, d21 / detJ}, {d02 / detJ, d12 / detJ, d22 / detJ}};

1983

1984

// Declare transformation matrix

1985

// Declare pointer to two dimensional array and initialise

1986

double **transform = new double *[num_derivatives];

1987

1988

for (unsigned int j = 0; j < num_derivatives; j++)

1989

{

1990

transform[j] = new double [num_derivatives];

1991

for (unsigned int k = 0; k < num_derivatives; k++)

1992

transform[j][k] = 1;

1993

}

1994

1995

// Construct transformation matrix

1996

for (unsigned int row = 0; row < num_derivatives; row++)

1997

{

1998

for (unsigned int col = 0; col < num_derivatives; col++)

1999

{

2000

for (unsigned int k = 0; k < n; k++)

2001

transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];

2002

}

2003

}

2004

2005

// Reset values

2006

for (unsigned int j = 0; j < 3*num_derivatives; j++)

2007

values[j] = 0;

2008

2009

if (0 <= i && i <= 3)

2010

{

2011

// Map degree of freedom to element degree of freedom

2012

const unsigned int dof = i;

2013

2014

// Generate scalings

2015

const double scalings_y_0 = 1;

2016

const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);

2017

const double scalings_z_0 = 1;

2018

const double scalings_z_1 = scalings_z_0*(0.5 - 0.5*z);

2019

2020

// Compute psitilde_a

2021

const double psitilde_a_0 = 1;

2022

const double psitilde_a_1 = x;

2023

2024

// Compute psitilde_bs

2025

const double psitilde_bs_0_0 = 1;

2026

const double psitilde_bs_0_1 = 1.5*y + 0.5;

2027

const double psitilde_bs_1_0 = 1;

2028

2029

// Compute psitilde_cs

2030

const double psitilde_cs_00_0 = 1;

2031

const double psitilde_cs_00_1 = 2*z + 1;

2032

const double psitilde_cs_01_0 = 1;

2033

const double psitilde_cs_10_0 = 1;

2034

2035

// Compute basisvalues

2036

const double basisvalue0 = 0.866025403784439*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_0;

2037

const double basisvalue1 = 2.73861278752583*psitilde_a_1*scalings_y_1*psitilde_bs_1_0*scalings_z_1*psitilde_cs_10_0;

2038

const double basisvalue2 = 1.58113883008419*psitilde_a_0*scalings_y_0*psitilde_bs_0_1*scalings_z_1*psitilde_cs_01_0;

2039

const double basisvalue3 = 1.11803398874989*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_1;

2040

2041

// Table(s) of coefficients

2042

const static double coefficients0[4][4] = \

2043

{{0.288675134594813, -0.182574185835055, -0.105409255338946, -0.074535599249993},

2044

{0.288675134594813, 0.182574185835055, -0.105409255338946, -0.074535599249993},

2045

{0.288675134594813, 0, 0.210818510677892, -0.074535599249993},

2046

{0.288675134594813, 0, 0, 0.223606797749979}};

2047

2048

// Interesting (new) part

2049

// Tables of derivatives of the polynomial base (transpose)

2050

const static double dmats0[4][4] = \

2051

{{0, 0, 0, 0},

2052

{6.32455532033676, 0, 0, 0},

2053

{0, 0, 0, 0},

2054

{0, 0, 0, 0}};

2055

2056

const static double dmats1[4][4] = \

2057

{{0, 0, 0, 0},

2058

{3.16227766016838, 0, 0, 0},

2059

{5.47722557505166, 0, 0, 0},

2060

{0, 0, 0, 0}};

2061

2062

const static double dmats2[4][4] = \

2063

{{0, 0, 0, 0},

2064

{3.16227766016838, 0, 0, 0},

2065

{1.82574185835055, 0, 0, 0},

2066

{5.16397779494322, 0, 0, 0}};

2067

2068

// Compute reference derivatives

2069

// Declare pointer to array of derivatives on FIAT element

2070

double *derivatives = new double [num_derivatives];

2071

2072

// Declare coefficients

2073

double coeff0_0 = 0;

2074

double coeff0_1 = 0;

2075

double coeff0_2 = 0;

2076

double coeff0_3 = 0;

2077

2078

// Declare new coefficients

2079

double new_coeff0_0 = 0;

2080

double new_coeff0_1 = 0;

2081

double new_coeff0_2 = 0;

2082

double new_coeff0_3 = 0;

2083

2084

// Loop possible derivatives

2085

for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)

2086

{

2087

// Get values from coefficients array

2088

new_coeff0_0 = coefficients0[dof][0];

2089

new_coeff0_1 = coefficients0[dof][1];

2090

new_coeff0_2 = coefficients0[dof][2];

2091

new_coeff0_3 = coefficients0[dof][3];

2092

2093

// Loop derivative order

2094

for (unsigned int j = 0; j < n; j++)

2095

{

2096

// Update old coefficients

2097

coeff0_0 = new_coeff0_0;

2098

coeff0_1 = new_coeff0_1;

2099

coeff0_2 = new_coeff0_2;

2100

coeff0_3 = new_coeff0_3;

2101

2102

if(combinations[deriv_num][j] == 0)

2103

{

2104

new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0] + coeff0_3*dmats0[3][0];

2105

new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1] + coeff0_3*dmats0[3][1];

2106

new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2] + coeff0_3*dmats0[3][2];

2107

new_coeff0_3 = coeff0_0*dmats0[0][3] + coeff0_1*dmats0[1][3] + coeff0_2*dmats0[2][3] + coeff0_3*dmats0[3][3];

2108

}

2109

if(combinations[deriv_num][j] == 1)

2110

{

2111

new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0] + coeff0_3*dmats1[3][0];

2112

new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1] + coeff0_3*dmats1[3][1];

2113

new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2] + coeff0_3*dmats1[3][2];

2114

new_coeff0_3 = coeff0_0*dmats1[0][3] + coeff0_1*dmats1[1][3] + coeff0_2*dmats1[2][3] + coeff0_3*dmats1[3][3];

2115

}

2116

if(combinations[deriv_num][j] == 2)

2117

{

2118

new_coeff0_0 = coeff0_0*dmats2[0][0] + coeff0_1*dmats2[1][0] + coeff0_2*dmats2[2][0] + coeff0_3*dmats2[3][0];

2119

new_coeff0_1 = coeff0_0*dmats2[0][1] + coeff0_1*dmats2[1][1] + coeff0_2*dmats2[2][1] + coeff0_3*dmats2[3][1];

2120

new_coeff0_2 = coeff0_0*dmats2[0][2] + coeff0_1*dmats2[1][2] + coeff0_2*dmats2[2][2] + coeff0_3*dmats2[3][2];

2121

new_coeff0_3 = coeff0_0*dmats2[0][3] + coeff0_1*dmats2[1][3] + coeff0_2*dmats2[2][3] + coeff0_3*dmats2[3][3];

2122

}

2123

2124

}

2125

// Compute derivatives on reference element as dot product of coefficients and basisvalues

2126

derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2 + new_coeff0_3*basisvalue3;

2127

}

2128

2129

// Transform derivatives back to physical element

2130

for (unsigned int row = 0; row < num_derivatives; row++)

2131

{

2132

for (unsigned int col = 0; col < num_derivatives; col++)

2133

{

2134

values[row] += transform[row][col]*derivatives[col];

2135

}

2136

}

2137

// Delete pointer to array of derivatives on FIAT element

2138

delete [] derivatives;

2139

2140

// Delete pointer to array of combinations of derivatives and transform

2141

for (unsigned int row = 0; row < num_derivatives; row++)

2142

{

2143

delete [] combinations[row];

2144

delete [] transform[row];

2145

}

2146

2147

delete [] combinations;

2148

delete [] transform;

2149

}

2150

2151

if (4 <= i && i <= 7)

2152

{

2153

// Map degree of freedom to element degree of freedom

2154

const unsigned int dof = i - 4;

2155

2156

// Generate scalings

2157

const double scalings_y_0 = 1;

2158

const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);

2159

const double scalings_z_0 = 1;

2160

const double scalings_z_1 = scalings_z_0*(0.5 - 0.5*z);

2161

2162

// Compute psitilde_a

2163

const double psitilde_a_0 = 1;

2164

const double psitilde_a_1 = x;

2165

2166

// Compute psitilde_bs

2167

const double psitilde_bs_0_0 = 1;

2168

const double psitilde_bs_0_1 = 1.5*y + 0.5;

2169

const double psitilde_bs_1_0 = 1;

2170

2171

// Compute psitilde_cs

2172

const double psitilde_cs_00_0 = 1;

2173

const double psitilde_cs_00_1 = 2*z + 1;

2174

const double psitilde_cs_01_0 = 1;

2175

const double psitilde_cs_10_0 = 1;

2176

2177

// Compute basisvalues

2178

const double basisvalue0 = 0.866025403784439*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_0;

2179

const double basisvalue1 = 2.73861278752583*psitilde_a_1*scalings_y_1*psitilde_bs_1_0*scalings_z_1*psitilde_cs_10_0;

2180

const double basisvalue2 = 1.58113883008419*psitilde_a_0*scalings_y_0*psitilde_bs_0_1*scalings_z_1*psitilde_cs_01_0;

2181

const double basisvalue3 = 1.11803398874989*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_1;

2182

2183

// Table(s) of coefficients

2184

const static double coefficients0[4][4] = \

2185

{{0.288675134594813, -0.182574185835055, -0.105409255338946, -0.074535599249993},

2186

{0.288675134594813, 0.182574185835055, -0.105409255338946, -0.074535599249993},

2187

{0.288675134594813, 0, 0.210818510677892, -0.074535599249993},

2188

{0.288675134594813, 0, 0, 0.223606797749979}};

2189

2190

// Interesting (new) part

2191

// Tables of derivatives of the polynomial base (transpose)

2192

const static double dmats0[4][4] = \

2193

{{0, 0, 0, 0},

2194

{6.32455532033676, 0, 0, 0},

2195

{0, 0, 0, 0},

2196

{0, 0, 0, 0}};

2197

2198

const static double dmats1[4][4] = \

2199

{{0, 0, 0, 0},

2200

{3.16227766016838, 0, 0, 0},

2201

{5.47722557505166, 0, 0, 0},

2202

{0, 0, 0, 0}};

2203

2204

const static double dmats2[4][4] = \

2205

{{0, 0, 0, 0},

2206

{3.16227766016838, 0, 0, 0},

2207

{1.82574185835055, 0, 0, 0},

2208

{5.16397779494322, 0, 0, 0}};

2209

2210

// Compute reference derivatives

2211

// Declare pointer to array of derivatives on FIAT element

2212

double *derivatives = new double [num_derivatives];

2213

2214

// Declare coefficients

2215

double coeff0_0 = 0;

2216

double coeff0_1 = 0;

2217

double coeff0_2 = 0;

2218

double coeff0_3 = 0;

2219

2220

// Declare new coefficients

2221

double new_coeff0_0 = 0;

2222

double new_coeff0_1 = 0;

2223

double new_coeff0_2 = 0;

2224

double new_coeff0_3 = 0;

2225

2226

// Loop possible derivatives

2227

for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)

2228

{

2229

// Get values from coefficients array

2230

new_coeff0_0 = coefficients0[dof][0];

2231

new_coeff0_1 = coefficients0[dof][1];

2232

new_coeff0_2 = coefficients0[dof][2];

2233

new_coeff0_3 = coefficients0[dof][3];

2234

2235

// Loop derivative order

2236

for (unsigned int j = 0; j < n; j++)

2237

{

2238

// Update old coefficients

2239

coeff0_0 = new_coeff0_0;

2240

coeff0_1 = new_coeff0_1;

2241

coeff0_2 = new_coeff0_2;

2242

coeff0_3 = new_coeff0_3;

2243

2244

if(combinations[deriv_num][j] == 0)

2245

{

2246

new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0] + coeff0_3*dmats0[3][0];

2247

new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1] + coeff0_3*dmats0[3][1];

2248

new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2] + coeff0_3*dmats0[3][2];

2249

new_coeff0_3 = coeff0_0*dmats0[0][3] + coeff0_1*dmats0[1][3] + coeff0_2*dmats0[2][3] + coeff0_3*dmats0[3][3];

2250

}

2251

if(combinations[deriv_num][j] == 1)

2252

{

2253

new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0] + coeff0_3*dmats1[3][0];

2254

new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1] + coeff0_3*dmats1[3][1];

2255

new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2] + coeff0_3*dmats1[3][2];

2256

new_coeff0_3 = coeff0_0*dmats1[0][3] + coeff0_1*dmats1[1][3] + coeff0_2*dmats1[2][3] + coeff0_3*dmats1[3][3];

2257

}

2258

if(combinations[deriv_num][j] == 2)

2259

{

2260

new_coeff0_0 = coeff0_0*dmats2[0][0] + coeff0_1*dmats2[1][0] + coeff0_2*dmats2[2][0] + coeff0_3*dmats2[3][0];

2261

new_coeff0_1 = coeff0_0*dmats2[0][1] + coeff0_1*dmats2[1][1] + coeff0_2*dmats2[2][1] + coeff0_3*dmats2[3][1];

2262

new_coeff0_2 = coeff0_0*dmats2[0][2] + coeff0_1*dmats2[1][2] + coeff0_2*dmats2[2][2] + coeff0_3*dmats2[3][2];

2263

new_coeff0_3 = coeff0_0*dmats2[0][3] + coeff0_1*dmats2[1][3] + coeff0_2*dmats2[2][3] + coeff0_3*dmats2[3][3];

2264

}

2265

2266

}

2267

// Compute derivatives on reference element as dot product of coefficients and basisvalues

2268

derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2 + new_coeff0_3*basisvalue3;

2269

}

2270

2271

// Transform derivatives back to physical element

2272

for (unsigned int row = 0; row < num_derivatives; row++)

2273

{

2274

for (unsigned int col = 0; col < num_derivatives; col++)

2275

{

2276

values[num_derivatives + row] += transform[row][col]*derivatives[col];

2277

}

2278

}

2279

// Delete pointer to array of derivatives on FIAT element

2280

delete [] derivatives;

2281

2282

// Delete pointer to array of combinations of derivatives and transform

2283

for (unsigned int row = 0; row < num_derivatives; row++)

2284

{

2285

delete [] combinations[row];

2286

delete [] transform[row];

2287

}

2288

2289

delete [] combinations;

2290

delete [] transform;

2291

}

2292

2293

if (8 <= i && i <= 11)

2294

{

2295

// Map degree of freedom to element degree of freedom

2296

const unsigned int dof = i - 8;

2297

2298

// Generate scalings

2299

const double scalings_y_0 = 1;

2300

const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);

2301

const double scalings_z_0 = 1;

2302

const double scalings_z_1 = scalings_z_0*(0.5 - 0.5*z);

2303

2304

// Compute psitilde_a

2305

const double psitilde_a_0 = 1;

2306

const double psitilde_a_1 = x;

2307

2308

// Compute psitilde_bs

2309

const double psitilde_bs_0_0 = 1;

2310

const double psitilde_bs_0_1 = 1.5*y + 0.5;

2311

const double psitilde_bs_1_0 = 1;

2312

2313

// Compute psitilde_cs

2314

const double psitilde_cs_00_0 = 1;

2315

const double psitilde_cs_00_1 = 2*z + 1;

2316

const double psitilde_cs_01_0 = 1;

2317

const double psitilde_cs_10_0 = 1;

2318

2319

// Compute basisvalues

2320

const double basisvalue0 = 0.866025403784439*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_0;

2321

const double basisvalue1 = 2.73861278752583*psitilde_a_1*scalings_y_1*psitilde_bs_1_0*scalings_z_1*psitilde_cs_10_0;

2322

const double basisvalue2 = 1.58113883008419*psitilde_a_0*scalings_y_0*psitilde_bs_0_1*scalings_z_1*psitilde_cs_01_0;

2323

const double basisvalue3 = 1.11803398874989*psitilde_a_0*scalings_y_0*psitilde_bs_0_0*scalings_z_0*psitilde_cs_00_1;

2324

2325

// Table(s) of coefficients

2326

const static double coefficients0[4][4] = \

2327

{{0.288675134594813, -0.182574185835055, -0.105409255338946, -0.074535599249993},

2328

{0.288675134594813, 0.182574185835055, -0.105409255338946, -0.074535599249993},

2329

{0.288675134594813, 0, 0.210818510677892, -0.074535599249993},

2330

{0.288675134594813, 0, 0, 0.223606797749979}};

2331

2332

// Interesting (new) part

2333

// Tables of derivatives of the polynomial base (transpose)

2334

const static double dmats0[4][4] = \

2335

{{0, 0, 0, 0},

2336

{6.32455532033676, 0, 0, 0},

2337

{0, 0, 0, 0},

2338

{0, 0, 0, 0}};

2339

2340

const static double dmats1[4][4] = \

2341

{{0, 0, 0, 0},

2342

{3.16227766016838, 0, 0, 0},

2343

{5.47722557505166, 0, 0, 0},

2344

{0, 0, 0, 0}};

2345

2346

const static double dmats2[4][4] = \

2347

{{0, 0, 0, 0},

2348

{3.16227766016838, 0, 0, 0},

2349

{1.82574185835055, 0, 0, 0},

2350

{5.16397779494322, 0, 0, 0}};

2351

2352

// Compute reference derivatives

2353

// Declare pointer to array of derivatives on FIAT element

2354

double *derivatives = new double [num_derivatives];

2355

2356

// Declare coefficients

2357

double coeff0_0 = 0;

2358

double coeff0_1 = 0;

2359

double coeff0_2 = 0;

2360

double coeff0_3 = 0;

2361

2362

// Declare new coefficients

2363

double new_coeff0_0 = 0;

2364

double new_coeff0_1 = 0;

2365

double new_coeff0_2 = 0;

2366

double new_coeff0_3 = 0;

2367

2368

// Loop possible derivatives

2369

for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)

2370

{

2371

// Get values from coefficients array

2372

new_coeff0_0 = coefficients0[dof][0];

2373

new_coeff0_1 = coefficients0[dof][1];

2374

new_coeff0_2 = coefficients0[dof][2];

2375

new_coeff0_3 = coefficients0[dof][3];

2376

2377

// Loop derivative order

2378

for (unsigned int j = 0; j < n; j++)

2379

{

2380

// Update old coefficients

2381

coeff0_0 = new_coeff0_0;

2382

coeff0_1 = new_coeff0_1;

2383

coeff0_2 = new_coeff0_2;

2384

coeff0_3 = new_coeff0_3;

2385

2386

if(combinations[deriv_num][j] == 0)

2387

{

2388

new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0] + coeff0_3*dmats0[3][0];

2389

new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1] + coeff0_3*dmats0[3][1];

2390

new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2] + coeff0_3*dmats0[3][2];

2391

new_coeff0_3 = coeff0_0*dmats0[0][3] + coeff0_1*dmats0[1][3] + coeff0_2*dmats0[2][3] + coeff0_3*dmats0[3][3];

2392

}

2393

if(combinations[deriv_num][j] == 1)

2394

{

2395

new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0] + coeff0_3*dmats1[3][0];

2396

new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1] + coeff0_3*dmats1[3][1];

2397

new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2] + coeff0_3*dmats1[3][2];

2398

new_coeff0_3 = coeff0_0*dmats1[0][3] + coeff0_1*dmats1[1][3] + coeff0_2*dmats1[2][3] + coeff0_3*dmats1[3][3];

2399

}

2400

if(combinations[deriv_num][j] == 2)

2401

{

2402

new_coeff0_0 = coeff0_0*dmats2[0][0] + coeff0_1*dmats2[1][0] + coeff0_2*dmats2[2][0] + coeff0_3*dmats2[3][0];

2403

new_coeff0_1 = coeff0_0*dmats2[0][1] + coeff0_1*dmats2[1][1] + coeff0_2*dmats2[2][1] + coeff0_3*dmats2[3][1];

2404

new_coeff0_2 = coeff0_0*dmats2[0][2] + coeff0_1*dmats2[1][2] + coeff0_2*dmats2[2][2] + coeff0_3*dmats2[3][2];

2405

new_coeff0_3 = coeff0_0*dmats2[0][3] + coeff0_1*dmats2[1][3] + coeff0_2*dmats2[2][3] + coeff0_3*dmats2[3][3];

2406

}

2407

2408

}

2409

// Compute derivatives on reference element as dot product of coefficients and basisvalues

2410

derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2 + new_coeff0_3*basisvalue3;

2411

}

2412

2413

// Transform derivatives back to physical element

2414

for (unsigned int row = 0; row < num_derivatives; row++)

2415

{

2416

for (unsigned int col = 0; col < num_derivatives; col++)

2417

{

2418

values[2*num_derivatives + row] += transform[row][col]*derivatives[col];

2419

}

2420

}

2421

// Delete pointer to array of derivatives on FIAT element

2422

delete [] derivatives;

2423

2424

// Delete pointer to array of combinations of derivatives and transform

2425

for (unsigned int row = 0; row < num_derivatives; row++)

2426

{

2427

delete [] combinations[row];

2428

delete [] transform[row];

2429

}

2430

2431

delete [] combinations;

2432

delete [] transform;

2433

}

2434

2435

}

2436

2437

/// Evaluate order n derivatives of all basis functions at given point in cell

2438

virtual void evaluate_basis_derivatives_all(unsigned int n,

2439

double* values,

2440

const double* coordinates,

2441

const ufc::cell& c) const

2442

{

2443

throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");

2444

}

2445

2446

/// Evaluate linear functional for dof i on the function f

2447

virtual double evaluate_dof(unsigned int i,

2448

const ufc::function& f,

2449

const ufc::cell& c) const

2450

{

2451

// The reference points, direction and weights:

2452

const static double X[12][1][3] = {{{0, 0, 0}}, {{1, 0, 0}}, {{0, 1, 0}}, {{0, 0, 1}}, {{0, 0, 0}}, {{1, 0, 0}}, {{0, 1, 0}}, {{0, 0, 1}}, {{0, 0, 0}}, {{1, 0, 0}}, {{0, 1, 0}}, {{0, 0, 1}}};

2453

const static double W[12][1] = {{1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}};

2454

const static double D[12][1][3] = {{{1, 0, 0}}, {{1, 0, 0}}, {{1, 0, 0}}, {{1, 0, 0}}, {{0, 1, 0}}, {{0, 1, 0}}, {{0, 1, 0}}, {{0, 1, 0}}, {{0, 0, 1}}, {{0, 0, 1}}, {{0, 0, 1}}, {{0, 0, 1}}};

2455

2456

const double * const * x = c.coordinates;

2457

double result = 0.0;

2458

// Iterate over the points:

2459

// Evaluate basis functions for affine mapping

2460

const double w0 = 1.0 - X[i][0][0] - X[i][0][1] - X[i][0][2];

2461

const double w1 = X[i][0][0];

2462

const double w2 = X[i][0][1];

2463

const double w3 = X[i][0][2];

2464

2465

// Compute affine mapping y = F(X)

2466

double y[3];

2467

y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0] + w3*x[3][0];

2468

y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1] + w3*x[3][1];

2469

y[2] = w0*x[0][2] + w1*x[1][2] + w2*x[2][2] + w3*x[3][2];

2470

2471

// Evaluate function at physical points

2472

double values[3];

2473

f.evaluate(values, y, c);

2474

2475

// Map function values using appropriate mapping

2476

// Affine map: Do nothing

2477

2478

// Note that we do not map the weights (yet).

2479

2480

// Take directional components

2481

for(int k = 0; k < 3; k++)

2482

result += values[k]*D[i][0][k];

2483

// Multiply by weights

2484

result *= W[i][0];

2485

2486

return result;

2487

}

2488

2489

/// Evaluate linear functionals for all dofs on the function f

2490

virtual void evaluate_dofs(double* values,

2491

const ufc::function& f,

2492

const ufc::cell& c) const

2493

{

2494

throw std::runtime_error("Not implemented (introduced in UFC v1.1).");

2495

}

2496

2497

/// Interpolate vertex values from dof values

2498

virtual void interpolate_vertex_values(double* vertex_values,

2499

const double* dof_values,

2500

const ufc::cell& c) const

2501

{

2502

// Evaluate at vertices and use affine mapping

2503

vertex_values[0] = dof_values[0];

2504

vertex_values[1] = dof_values[1];

2505

vertex_values[2] = dof_values[2];

2506

vertex_values[3] = dof_values[3];

2507

// Evaluate at vertices and use affine mapping

2508

vertex_values[4] = dof_values[4];

2509

vertex_values[5] = dof_values[5];

2510

vertex_values[6] = dof_values[6];

2511

vertex_values[7] = dof_values[7];

2512

// Evaluate at vertices and use affine mapping

2513

vertex_values[8] = dof_values[8];

2514

vertex_values[9] = dof_values[9];

2515

vertex_values[10] = dof_values[10];

2516

vertex_values[11] = dof_values[11];

2517

}

2518

2519

/// Return the number of sub elements (for a mixed element)

2520

virtual unsigned int num_sub_elements() const

2521

{

2522

return 3;

2523

}

2524

2525

/// Create a new finite element for sub element i (for a mixed element)

2526

virtual ufc::finite_element* create_sub_element(unsigned int i) const

2527

{

2528

switch ( i )

2529

{

2530

case 0:

2531

return new ffc_16_finite_element_0_0();

2532

break;

2533

case 1:

2534

return new ffc_16_finite_element_0_1();

2535

break;

2536

case 2:

2537

return new ffc_16_finite_element_0_2();

2538

break;

2539

}

2540

return 0;

2541

}

2542

2543

};

2544

2545

/// This class defines the interface for a local-to-global mapping of

2546

/// degrees of freedom (dofs).

2547

2548

class ffc_16_dof_map_0_0: public ufc::dof_map

2549

{

2550

private:

2551

2552

unsigned int __global_dimension;

2553

2554

public:

2555

2556

/// Constructor

2557

ffc_16_dof_map_0_0() : ufc::dof_map()

2558

{

2559

__global_dimension = 0;

2560

}

2561

2562

/// Destructor

2563

virtual ~ffc_16_dof_map_0_0()

2564

{

2565

// Do nothing

2566

}

2567

2568

/// Return a string identifying the dof map

2569

virtual const char* signature() const

2570

{

2571

return "FFC dof map for Lagrange finite element of degree 1 on a tetrahedron";

2572

}

2573

2574

/// Return true iff mesh entities of topological dimension d are needed

2575

virtual bool needs_mesh_entities(unsigned int d) const

2576

{

2577

switch ( d )

2578

{

2579

case 0:

2580

return true;

2581

break;

2582

case 1:

2583

return false;

2584

break;

2585

case 2:

2586

return false;

2587

break;

2588

case 3:

2589

return false;

2590

break;

2591

}

2592

return false;

2593

}

2594

2595

/// Initialize dof map for mesh (return true iff init_cell() is needed)

2596

virtual bool init_mesh(const ufc::mesh& m)

2597

{

2598

__global_dimension = m.num_entities[0];

2599

return false;

2600

}

2601

2602

/// Initialize dof map for given cell

2603

virtual void init_cell(const ufc::mesh& m,

2604

const ufc::cell& c)

2605

{

2606

// Do nothing

2607

}

2608

2609

/// Finish initialization of dof map for cells

2610

virtual void init_cell_finalize()

2611

{

2612

// Do nothing

2613

}

2614

2615

/// Return the dimension of the global finite element function space

2616

virtual unsigned int global_dimension() const

2617

{

2618

return __global_dimension;

2619

}

2620

2621

/// Return the dimension of the local finite element function space

2622

virtual unsigned int local_dimension() const

2623

{

2624

return 4;

2625

}

2626

2627

// Return the geometric dimension of the coordinates this dof map provides

2628

virtual unsigned int geometric_dimension() const

2629

{

2630

throw std::runtime_error("Not implemented (introduced in UFC v1.1).");

2631

}

2632

2633

/// Return the number of dofs on each cell facet

2634

virtual unsigned int num_facet_dofs() const

2635

{

2636

return 3;

2637

}

2638

2639

/// Return the number of dofs associated with each cell entity of dimension d

2640

virtual unsigned int num_entity_dofs(unsigned int d) const

2641

{

2642

throw std::runtime_error("Not implemented (introduced in UFC v1.1).");

2643

}

2644

2645

/// Tabulate the local-to-global mapping of dofs on a cell

2646

virtual void tabulate_dofs(unsigned int* dofs,

2647

const ufc::mesh& m,

2648

const ufc::cell& c) const

2649

{

2650

dofs[0] = c.entity_indices[0][0];

2651

dofs[1] = c.entity_indices[0][1];

2652

dofs[2] = c.entity_indices[0][2];

2653

dofs[3] = c.entity_indices[0][3];

2654

}

2655

2656

/// Tabulate the local-to-local mapping from facet dofs to cell dofs

2657

virtual void tabulate_facet_dofs(unsigned int* dofs,

2658

unsigned int facet) const

2659

{

2660

switch ( facet )

2661

{

2662

case 0:

2663

dofs[0] = 1;

2664

dofs[1] = 2;

2665

dofs[2] = 3;

2666

break;

2667

case 1:

2668

dofs[0] = 0;

2669

dofs[1] = 2;

2670

dofs[2] = 3;

2671

break;

2672

case 2:

2673

dofs[0] = 0;

2674

dofs[1] = 1;

2675

dofs[2] = 3;

2676

break;

2677

case 3:

2678

dofs[0] = 0;

2679

dofs[1] = 1;

2680

dofs[2] = 2;

2681

break;

2682

}

2683

}

2684

2685

/// Tabulate the local-to-local mapping of dofs on entity (d, i)

2686

virtual void tabulate_entity_dofs(unsigned int* dofs,

2687

unsigned int d, unsigned int i) const

2688

{

2689

throw std::runtime_error("Not implemented (introduced in UFC v1.1).");

2690

}

2691

2692

/// Tabulate the coordinates of all dofs on a cell

2693

virtual void tabulate_coordinates(double** coordinates,

2694

const ufc::cell& c) const

2695

{

2696

const double * const * x = c.coordinates;

2697

coordinates[0][0] = x[0][0];

2698

coordinates[0][1] = x[0][1];

2699

coordinates[0][2] = x[0][2];

2700

coordinates[1][0] = x[1][0];

2701

coordinates[1][1] = x[1][1];

2702

coordinates[1][2] = x[1][2];

2703

coordinates[2][0] = x[2][0];

2704

coordinates[2][1] = x[2][1];

2705

coordinates[2][2] = x[2][2];

2706

coordinates[3][0] = x[3][0];

2707

coordinates[3][1] = x[3][1];

2708

coordinates[3][2] = x[3][2];

2709

}

2710

2711

/// Return the number of sub dof maps (for a mixed element)

2712

virtual unsigned int num_sub_dof_maps() const

2713

{

2714

return 1;

2715

}

2716

2717

/// Create a new dof_map for sub dof map i (for a mixed element)

2718

virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const

2719

{

2720

return new ffc_16_dof_map_0_0();

2721

}

2722

2723

};

2724

2725

/// This class defines the interface for a local-to-global mapping of

2726

/// degrees of freedom (dofs).

2727

2728

class ffc_16_dof_map_0_1: public ufc::dof_map

2729

{

2730

private:

2731

2732

unsigned int __global_dimension;

2733

2734

public:

2735

2736

/// Constructor

2737

ffc_16_dof_map_0_1() : ufc::dof_map()

2738

{

2739

__global_dimension = 0;

2740

}

2741

2742

/// Destructor

2743

virtual ~ffc_16_dof_map_0_1()

2744

{

2745

// Do nothing

2746

}

2747

2748

/// Return a string identifying the dof map

2749

virtual const char* signature() const

2750

{

2751

return "FFC dof map for Lagrange finite element of degree 1 on a tetrahedron";

2752

}

2753

2754

/// Return true iff mesh entities of topological dimension d are needed

2755

virtual bool needs_mesh_entities(unsigned int d) const

2756

{

2757

switch ( d )

2758

{

2759

case 0:

2760

return true;

2761

break;

2762

case 1:

2763

return false;

2764

break;

2765

case 2:

2766

return false;

2767

break;

2768

case 3:

2769

return false;

2770

break;

2771

}

2772

return false;

2773

}

2774

2775

/// Initialize dof map for mesh (return true iff init_cell() is needed)

2776

virtual bool init_mesh(const ufc::mesh& m)

2777

{

2778

__global_dimension = m.num_entities[0];

2779

return false;

2780

}

2781

2782

/// Initialize dof map for given cell

2783

virtual void init_cell(const ufc::mesh& m,

2784

const ufc::cell& c)

2785

{

2786

// Do nothing

2787

}

2788

2789

/// Finish initialization of dof map for cells

2790

virtual void init_cell_finalize()

2791

{

2792

// Do nothing

2793

}

2794

2795

/// Return the dimension of the global finite element function space

2796

virtual unsigned int global_dimension() const

2797

{

2798

return __global_dimension;

2799

}

2800

2801

/// Return the dimension of the local finite element function space

2802

virtual unsigned int local_dimension() const

2803

{

2804

return 4;

2805

}

2806

2807

// Return the geometric dimension of the coordinates this dof map provides

2808

virtual unsigned int geometric_dimension() const

2809

{

2810

throw std::runtime_error("Not implemented (introduced in UFC v1.1).");

2811

}

2812

2813

/// Return the number of dofs on each cell facet

2814

virtual unsigned int num_facet_dofs() const

2815

{

2816

return 3;

2817

}

2818

2819

/// Return the number of dofs associated with each cell entity of dimension d

2820

virtual unsigned int num_entity_dofs(unsigned int d) const

2821

{

2822

throw std::runtime_error("Not implemented (introduced in UFC v1.1).");

2823

}

2824

2825

/// Tabulate the local-to-global mapping of dofs on a cell

2826

virtual void tabulate_dofs(unsigned int* dofs,

2827

const ufc::mesh& m,

2828

const ufc::cell& c) const

2829

{

2830

dofs[0] = c.entity_indices[0][0];

2831

dofs[1] = c.entity_indices[0][1];

2832

dofs[2] = c.entity_indices[0][2];

2833

dofs[3] = c.entity_indices[0][3];

2834

}

2835

2836

/// Tabulate the local-to-local mapping from facet dofs to cell dofs

2837

virtual void tabulate_facet_dofs(unsigned int* dofs,

2838

unsigned int facet) const

2839

{

2840

switch ( facet )

2841

{

2842

case 0:

2843

dofs[0] = 1;

2844

dofs[1] = 2;

2845

dofs[2] = 3;

2846

break;

2847

case 1:

2848

dofs[0] = 0;

2849

dofs[1] = 2;

2850

dofs[2] = 3;

2851

break;

2852

case 2:

2853

dofs[0] = 0;

2854

dofs[1] = 1;

2855

dofs[2] = 3;

2856

break;

2857

case 3:

2858

dofs[0] = 0;

2859

dofs[1] = 1;

2860

dofs[2] = 2;

2861

break;

2862

}

2863

}

2864

2865

/// Tabulate the local-to-local mapping of dofs on entity (d, i)

2866

virtual void tabulate_entity_dofs(unsigned int* dofs,

2867

unsigned int d, unsigned int i) const

2868

{

2869

throw std::runtime_error("Not implemented (introduced in UFC v1.1).");

2870

}

2871

2872

/// Tabulate the coordinates of all dofs on a cell

2873

virtual void tabulate_coordinates(double** coordinates,

2874

const ufc::cell& c) const

2875

{

2876

const double * const * x = c.coordinates;

2877

coordinates[0][0] = x[0][0];

2878

coordinates[0][1] = x[0][1];

2879

coordinates[0][2] = x[0][2];

2880

coordinates[1][0] = x[1][0];

2881

coordinates[1][1] = x[1][1];

2882

coordinates[1][2] = x[1][2];

2883

coordinates[2][0] = x[2][0];

2884

coordinates[2][1] = x[2][1];

2885

coordinates[2][2] = x[2][2];

2886

coordinates[3][0] = x[3][0];

2887

coordinates[3][1] = x[3][1];

2888

coordinates[3][2] = x[3][2];

2889

}

2890

2891

/// Return the number of sub dof maps (for a mixed element)

2892

virtual unsigned int num_sub_dof_maps() const

2893

{

2894

return 1;

2895

}

2896

2897

/// Create a new dof_map for sub dof map i (for a mixed element)

2898

virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const

2899

{

2900

return new ffc_16_dof_map_0_1();

2901

}

2902

2903

};

2904

2905

/// This class defines the interface for a local-to-global mapping of

2906

/// degrees of freedom (dofs).

2907

2908

class ffc_16_dof_map_0_2: public ufc::dof_map

2909

{

2910

private:

2911

2912

unsigned int __global_dimension;

2913

2914

public:

2915

2916

/// Constructor

2917

ffc_16_dof_map_0_2() : ufc::dof_map()

2918

{

2919

__global_dimension = 0;

2920

}

2921

2922

/// Destructor

2923

virtual ~ffc_16_dof_map_0_2()

2924

{

2925

// Do nothing

2926

}

2927

2928

/// Return a string identifying the dof map

2929

virtual const char* signature() const

2930

{

2931

return "FFC dof map for Lagrange finite element of degree 1 on a tetrahedron";

2932

}

2933

2934

/// Return true iff mesh entities of topological dimension d are needed

2935

virtual bool needs_mesh_entities(unsigned int d) const

2936

{

2937

switch ( d )

2938

{

2939

case 0:

2940

return true;

2941

break;

2942

case 1:

2943

return false;

2944

break;

2945

case 2:

2946

return false;

2947

break;

2948

case 3:

2949

return false;

2950

break;

2951

}

2952

return false;

2953

}

2954

2955

/// Initialize dof map for mesh (return true iff init_cell() is needed)

2956

virtual bool init_mesh(const ufc::mesh& m)

2957

{

2958

__global_dimension = m.num_entities[0];

2959

return false;

2960

}

2961

2962

/// Initialize dof map for given cell

2963

virtual void init_cell(const ufc::mesh& m,

2964

const ufc::cell& c)

2965

{

2966

// Do nothing

2967

}

2968

2969

/// Finish initialization of dof map for cells

2970

virtual void init_cell_finalize()

2971

{

2972

// Do nothing

2973

}

2974

2975

/// Return the dimension of the global finite element function space

2976

virtual unsigned int global_dimension() const

2977

{

2978

return __global_dimension;

2979

}

2980

2981

/// Return the dimension of the local finite element function space

2982

virtual unsigned int local_dimension() const

2983

{

2984

return 4;

2985

}

2986

2987

// Return the geometric dimension of the coordinates this dof map provides

2988

virtual unsigned int geometric_dimension() const

2989

{

2990

throw std::runtime_error("Not implemented (introduced in UFC v1.1).");

2991

}

2992

2993

/// Return the number of dofs on each cell facet

2994

virtual unsigned int num_facet_dofs() const

2995

{

2996

return 3;

2997

}

2998

2999

/// Return the number of dofs associated with each cell entity of dimension d

3000

virtual unsigned int num_entity_dofs(unsigned int d) const

3001

{

3002

throw std::runtime_error("Not implemented (introduced in UFC v1.1).");

3003

}

3004

3005

/// Tabulate the local-to-global mapping of dofs on a cell

3006

virtual void tabulate_dofs(unsigned int* dofs,

3007

const ufc::mesh& m,

3008

const ufc::cell& c) const

3009

{

3010

dofs[0] = c.entity_indices[0][0];

3011

dofs[1] = c.entity_indices[0][1];

3012

dofs[2] = c.entity_indices[0][2];

3013

dofs[3] = c.entity_indices[0][3];

3014

}

3015

3016

/// Tabulate the local-to-local mapping from facet dofs to cell dofs

3017

virtual void tabulate_facet_dofs(unsigned int* dofs,

3018

unsigned int facet) const

3019

{

3020

switch ( facet )

3021

{

3022

case 0:

3023

dofs[0] = 1;

3024

dofs[1] = 2;

3025

dofs[2] = 3;

3026

break;

3027

case 1:

3028

dofs[0] = 0;

3029

dofs[1] = 2;

3030

dofs[2] = 3;

3031

break;

3032

case 2:

3033

dofs[0] = 0;

3034

dofs[1] = 1;

3035

dofs[2] = 3;

3036

break;

3037

case 3:

3038

dofs[0] = 0;

3039

dofs[1] = 1;

3040

dofs[2] = 2;

3041

break;

3042

}

3043

}

3044

3045

/// Tabulate the local-to-local mapping of dofs on entity (d, i)

3046

virtual void tabulate_entity_dofs(unsigned int* dofs,

3047

unsigned int d, unsigned int i) const

3048

{

3049

throw std::runtime_error("Not implemented (introduced in UFC v1.1).");

3050

}

3051

3052

/// Tabulate the coordinates of all dofs on a cell

3053

virtual void tabulate_coordinates(double** coordinates,

3054

const ufc::cell& c) const

3055

{

3056

const double * const * x = c.coordinates;

3057

coordinates[0][0] = x[0][0];

3058

coordinates[0][1] = x[0][1];

3059

coordinates[0][2] = x[0][2];

3060

coordinates[1][0] = x[1][0];

3061

coordinates[1][1] = x[1][1];

3062

coordinates[1][2] = x[1][2];

3063

coordinates[2][0] = x[2][0];

3064

coordinates[2][1] = x[2][1];

3065

coordinates[2][2] = x[2][2];

3066

coordinates[3][0] = x[3][0];

3067

coordinates[3][1] = x[3][1];

3068

coordinates[3][2] = x[3][2];

3069

}

3070

3071

/// Return the number of sub dof maps (for a mixed element)

3072

virtual unsigned int num_sub_dof_maps() const

3073

{

3074

return 1;

3075

}

3076

3077

/// Create a new dof_map for sub dof map i (for a mixed element)

3078

virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const

3079

{

3080

return new ffc_16_dof_map_0_2();

3081

}

3082

3083

};

3084

3085

/// This class defines the interface for a local-to-global mapping of

3086

/// degrees of freedom (dofs).

3087

3088

class ffc_16_dof_map_0: public ufc::dof_map

3089

{

3090

private:

3091

3092

unsigned int __global_dimension;

3093

3094

public:

3095

3096

/// Constructor

3097

ffc_16_dof_map_0() : ufc::dof_map()

3098

{

3099

__global_dimension = 0;

3100

}

3101

3102

/// Destructor

3103

virtual ~ffc_16_dof_map_0()

3104

{

3105

// Do nothing

3106

}

3107

3108

/// Return a string identifying the dof map

3109

virtual const char* signature() const

3110

{

3111

return "FFC dof map for Mixed finite element: [Lagrange finite element of degree 1 on a tetrahedron, Lagrange finite element of degree 1 on a tetrahedron, Lagrange finite element of degree 1 on a tetrahedron]";

3112

}

3113

3114

/// Return true iff mesh entities of topological dimension d are needed

3115

virtual bool needs_mesh_entities(unsigned int d) const

3116

{

3117

switch ( d )

3118

{

3119

case 0:

3120

return true;

3121

break;

3122

case 1:

3123

return false;

3124

break;

3125

case 2:

3126

return false;

3127

break;

3128

case 3:

3129

return false;

3130

break;

3131

}

3132

return false;

3133

}

3134

3135

/// Initialize dof map for mesh (return true iff init_cell() is needed)

3136

virtual bool init_mesh(const ufc::mesh& m)

3137

{

3138

__global_dimension = 3*m.num_entities[0];

3139

return false;

3140

}

3141

3142

/// Initialize dof map for given cell

3143

virtual void init_cell(const ufc::mesh& m,

3144

const ufc::cell& c)

3145

{

3146

// Do nothing

3147

}

3148

3149

/// Finish initialization of dof map for cells

3150

virtual void init_cell_finalize()

3151

{

3152

// Do nothing

3153

}

3154

3155

/// Return the dimension of the global finite element function space

3156

virtual unsigned int global_dimension() const

3157

{

3158

return __global_dimension;

3159

}

3160

3161

/// Return the dimension of the local finite element function space

3162

virtual unsigned int local_dimension() const

3163

{

3164

return 12;

3165

}

3166

3167

// Return the geometric dimension of the coordinates this dof map provides

3168

virtual unsigned int geometric_dimension() const

3169

{

3170

throw std::runtime_error("Not implemented (introduced in UFC v1.1).");

3171

}

3172

3173

/// Return the number of dofs on each cell facet

3174

virtual unsigned int num_facet_dofs() const

3175

{

3176

return 9;

3177

}

3178

3179

/// Return the number of dofs associated with each cell entity of dimension d

3180

virtual unsigned int num_entity_dofs(unsigned int d) const

3181

{

3182

throw std::runtime_error("Not implemented (introduced in UFC v1.1).");

3183

}

3184

3185

/// Tabulate the local-to-global mapping of dofs on a cell

3186

virtual void tabulate_dofs(unsigned int* dofs,

3187

const ufc::mesh& m,

3188

const ufc::cell& c) const

3189

{

3190

dofs[0] = c.entity_indices[0][0];

3191

dofs[1] = c.entity_indices[0][1];

3192

dofs[2] = c.entity_indices[0][2];

3193

dofs[3] = c.entity_indices[0][3];

3194

unsigned int offset = m.num_entities[0];

3195

dofs[4] = offset + c.entity_indices[0][0];

3196

dofs[5] = offset + c.entity_indices[0][1];

3197

dofs[6] = offset + c.entity_indices[0][2];

3198

dofs[7] = offset + c.entity_indices[0][3];

3199

offset = offset + m.num_entities[0];

3200

dofs[8] = offset + c.entity_indices[0][0];

3201

dofs[9] = offset + c.entity_indices[0][1];

3202

dofs[10] = offset + c.entity_indices[0][2];

3203

dofs[11] = offset + c.entity_indices[0][3];

3204

}

3205

3206

/// Tabulate the local-to-local mapping from facet dofs to cell dofs

3207

virtual void tabulate_facet_dofs(unsigned int* dofs,

3208

unsigned int facet) const

3209

{

3210

switch ( facet )

3211

{

3212

case 0:

3213

dofs[0] = 1;

3214

dofs[1] = 2;

3215

dofs[2] = 3;

3216

dofs[3] = 5;

3217

dofs[4] = 6;

3218

dofs[5] = 7;

3219

dofs[6] = 9;

3220

dofs[7] = 10;

3221

dofs[8] = 11;

3222

break;

3223

case 1:

3224

dofs[0] = 0;

3225

dofs[1] = 2;

3226

dofs[2] = 3;

3227

dofs[3] = 4;

3228

dofs[4] = 6;

3229

dofs[5] = 7;

3230

dofs[6] = 8;

3231

dofs[7] = 10;

3232

dofs[8] = 11;

3233

break;

3234

case 2:

3235

dofs[0] = 0;

3236

dofs[1] = 1;

3237

dofs[2] = 3;

3238

dofs[3] = 4;

3239

dofs[4] = 5;

3240

dofs[5] = 7;

3241

dofs[6] = 8;

3242

dofs[7] = 9;

3243

dofs[8] = 11;

3244

break;

3245

case 3:

3246

dofs[0] = 0;

3247

dofs[1] = 1;

3248

dofs[2] = 2;

3249

dofs[3] = 4;

3250

dofs[4] = 5;

3251

dofs[5] = 6;

3252

dofs[6] = 8;

3253

dofs[7] = 9;

3254

dofs[8] = 10;

3255

break;

3256

}

3257

}

3258

3259

/// Tabulate the local-to-local mapping of dofs on entity (d, i)

3260

virtual void tabulate_entity_dofs(unsigned int* dofs,

3261

unsigned int d, unsigned int i) const

3262

{

3263

throw std::runtime_error("Not implemented (introduced in UFC v1.1).");

3264

}

3265

3266

/// Tabulate the coordinates of all dofs on a cell

3267

virtual void tabulate_coordinates(double** coordinates,

3268

const ufc::cell& c) const

3269

{

3270

const double * const * x = c.coordinates;

3271

coordinates[0][0] = x[0][0];

3272

coordinates[0][1] = x[0][1];

3273

coordinates[0][2] = x[0][2];

3274

coordinates[1][0] = x[1][0];

3275

coordinates[1][1] = x[1][1];

3276

coordinates[1][2] = x[1][2];

3277

coordinates[2][0] = x[2][0];

3278

coordinates[2][1] = x[2][1];

3279

coordinates[2][2] = x[2][2];

3280

coordinates[3][0] = x[3][0];

3281

coordinates[3][1] = x[3][1];

3282

coordinates[3][2] = x[3][2];

3283

coordinates[4][0] = x[0][0];

3284

coordinates[4][1] = x[0][1];

3285

coordinates[4][2] = x[0][2];

3286

coordinates[5][0] = x[1][0];

3287

coordinates[5][1] = x[1][1];

3288

coordinates[5][2] = x[1][2];

3289

coordinates[6][0] = x[2][0];

3290

coordinates[6][1] = x[2][1];

3291

coordinates[6][2] = x[2][2];

3292

coordinates[7][0] = x[3][0];

3293

coordinates[7][1] = x[3][1];

3294

coordinates[7][2] = x[3][2];

3295

coordinates[8][0] = x[0][0];

3296

coordinates[8][1] = x[0][1];

3297

coordinates[8][2] = x[0][2];

3298

coordinates[9][0] = x[1][0];

3299

coordinates[9][1] = x[1][1];

3300

coordinates[9][2] = x[1][2];

3301

coordinates[10][0] = x[2][0];

3302

coordinates[10][1] = x[2][1];

3303

coordinates[10][2] = x[2][2];

3304

coordinates[11][0] = x[3][0];

3305

coordinates[11][1] = x[3][1];

3306

coordinates[11][2] = x[3][2];

3307

}

3308

3309

/// Return the number of sub dof maps (for a mixed element)

3310

virtual unsigned int num_sub_dof_maps() const

3311

{

3312

return 3;

3313

}

3314

3315

/// Create a new dof_map for sub dof map i (for a mixed element)

3316

virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const

3317

{

3318

switch ( i )

3319

{

3320

case 0:

3321

return new ffc_16_dof_map_0_0();

3322

break;

3323

case 1:

3324

return new ffc_16_dof_map_0_1();

3325

break;

3326

case 2:

3327

return new ffc_16_dof_map_0_2();

3328

break;

3329

}

3330

return 0;

3331

}

3332

3333

};

3334

3335

#endif

Older »