~njansson/dolfin/hpc

return "Mixed finite element: [Discontinuous Lagrange finite element of degree 1 on a tetrahedron, Discontinuous Lagrange finite element of degree 1 on a tetrahedron, Discontinuous 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_23_finite_element_0_0();

2532

break;

2533

case 1:

2534

return new ffc_23_finite_element_0_1();

2535

break;

2536

case 2:

2537

return new ffc_23_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_23_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_23_dof_map_0_0() : ufc::dof_map()

2558

{

2559

__global_dimension = 0;

2560

}

2561

2562

/// Destructor

2563

virtual ~ffc_23_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 Discontinuous 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 false;

2581

break;

2582

case 1:

2583

return false;

2584

break;

2585

case 2:

2586

return false;

2587

break;

2588

case 3:

2589

return true;

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 = 4*m.num_entities[3];

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 0;

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] = 4*c.entity_indices[3][0];

2651

dofs[1] = 4*c.entity_indices[3][0] + 1;

2652

dofs[2] = 4*c.entity_indices[3][0] + 2;

2653

dofs[3] = 4*c.entity_indices[3][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

2664

break;

2665

case 1:

2666

2667

break;

2668

case 2:

2669

2670

break;

2671

case 3:

2672

2673

break;

2674

}

2675

}

2676

2677

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

2678

virtual void tabulate_entity_dofs(unsigned int* dofs,

2679

unsigned int d, unsigned int i) const

2680

{

2681

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

2682

}

2683

2684

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

2685

virtual void tabulate_coordinates(double** coordinates,

2686

const ufc::cell& c) const

2687

{

2688

const double * const * x = c.coordinates;

2689

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

2690

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

2691

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

2692

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

2693

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

2694

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

2695

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

2696

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

2697

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

2698

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

2699

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

2700

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

2701

}

2702

2703

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

2704

virtual unsigned int num_sub_dof_maps() const

2705

{

2706

return 1;

2707

}

2708

2709

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

2710

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

2711

{

2712

return new ffc_23_dof_map_0_0();

2713

}

2714

2715

};

2716

2717

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

2718

/// degrees of freedom (dofs).

2719

2720

class ffc_23_dof_map_0_1: public ufc::dof_map

2721

{

2722

private:

2723

2724

unsigned int __global_dimension;

2725

2726

public:

2727

2728

/// Constructor

2729

ffc_23_dof_map_0_1() : ufc::dof_map()

2730

{

2731

__global_dimension = 0;

2732

}

2733

2734

/// Destructor

2735

virtual ~ffc_23_dof_map_0_1()

2736

{

2737

// Do nothing

2738

}

2739

2740

/// Return a string identifying the dof map

2741

virtual const char* signature() const

2742

{

2743

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

2744

}

2745

2746

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

2747

virtual bool needs_mesh_entities(unsigned int d) const

2748

{

2749

switch ( d )

2750

{

2751

case 0:

2752

return false;

2753

break;

2754

case 1:

2755

return false;

2756

break;

2757

case 2:

2758

return false;

2759

break;

2760

case 3:

2761

return true;

2762

break;

2763

}

2764

return false;

2765

}

2766

2767

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

2768

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

2769

{

2770

__global_dimension = 4*m.num_entities[3];

2771

return false;

2772

}

2773

2774

/// Initialize dof map for given cell

2775

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

2776

const ufc::cell& c)

2777

{

2778

// Do nothing

2779

}

2780

2781

/// Finish initialization of dof map for cells

2782

virtual void init_cell_finalize()

2783

{

2784

// Do nothing

2785

}

2786

2787

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

2788

virtual unsigned int global_dimension() const

2789

{

2790

return __global_dimension;

2791

}

2792

2793

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

2794

virtual unsigned int local_dimension() const

2795

{

2796

return 4;

2797

}

2798

2799

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

2800

virtual unsigned int geometric_dimension() const

2801

{

2802

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

2803

}

2804

2805

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

2806

virtual unsigned int num_facet_dofs() const

2807

{

2808

return 0;

2809

}

2810

2811

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

2812

virtual unsigned int num_entity_dofs(unsigned int d) const

2813

{

2814

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

2815

}

2816

2817

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

2818

virtual void tabulate_dofs(unsigned int* dofs,

2819

const ufc::mesh& m,

2820

const ufc::cell& c) const

2821

{

2822

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

2823

dofs[1] = 4*c.entity_indices[3][0] + 1;

2824

dofs[2] = 4*c.entity_indices[3][0] + 2;

2825

dofs[3] = 4*c.entity_indices[3][0] + 3;

2826

}

2827

2828

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

2829

virtual void tabulate_facet_dofs(unsigned int* dofs,

2830

unsigned int facet) const

2831

{

2832

switch ( facet )

2833

{

2834

case 0:

2835

2836

break;

2837

case 1:

2838

2839

break;

2840

case 2:

2841

2842

break;

2843

case 3:

2844

2845

break;

2846

}

2847

}

2848

2849

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

2850

virtual void tabulate_entity_dofs(unsigned int* dofs,

2851

unsigned int d, unsigned int i) const

2852

{

2853

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

2854

}

2855

2856

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

2857

virtual void tabulate_coordinates(double** coordinates,

2858

const ufc::cell& c) const

2859

{

2860

const double * const * x = c.coordinates;

2861

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

2862

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

2863

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

2864

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

2865

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

2866

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

2867

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

2868

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

2869

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

2870

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

2871

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

2872

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

2873

}

2874

2875

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

2876

virtual unsigned int num_sub_dof_maps() const

2877

{

2878

return 1;

2879

}

2880

2881

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

2882

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

2883

{

2884

return new ffc_23_dof_map_0_1();

2885

}

2886

2887

};

2888

2889

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

2890

/// degrees of freedom (dofs).

2891

2892

class ffc_23_dof_map_0_2: public ufc::dof_map

2893

{

2894

private:

2895

2896

unsigned int __global_dimension;

2897

2898

public:

2899

2900

/// Constructor

2901

ffc_23_dof_map_0_2() : ufc::dof_map()

2902

{

2903

__global_dimension = 0;

2904

}

2905

2906

/// Destructor

2907

virtual ~ffc_23_dof_map_0_2()

2908

{

2909

// Do nothing

2910

}

2911

2912

/// Return a string identifying the dof map

2913

virtual const char* signature() const

2914

{

2915

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

2916

}

2917

2918

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

2919

virtual bool needs_mesh_entities(unsigned int d) const

2920

{

2921

switch ( d )

2922

{

2923

case 0:

2924

return false;

2925

break;

2926

case 1:

2927

return false;

2928

break;

2929

case 2:

2930

return false;

2931

break;

2932

case 3:

2933

return true;

2934

break;

2935

}

2936

return false;

2937

}

2938

2939

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

2940

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

2941

{

2942

__global_dimension = 4*m.num_entities[3];

2943

return false;

2944

}

2945

2946

/// Initialize dof map for given cell

2947

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

2948

const ufc::cell& c)

2949

{

2950

// Do nothing

2951

}

2952

2953

/// Finish initialization of dof map for cells

2954

virtual void init_cell_finalize()

2955

{

2956

// Do nothing

2957

}

2958

2959

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

2960

virtual unsigned int global_dimension() const

2961

{

2962

return __global_dimension;

2963

}

2964

2965

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

2966

virtual unsigned int local_dimension() const

2967

{

2968

return 4;

2969

}

2970

2971

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

2972

virtual unsigned int geometric_dimension() const

2973

{

2974

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

2975

}

2976

2977

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

2978

virtual unsigned int num_facet_dofs() const

2979

{

2980

return 0;

2981

}

2982

2983

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

2984

virtual unsigned int num_entity_dofs(unsigned int d) const

2985

{

2986

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

2987

}

2988

2989

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

2990

virtual void tabulate_dofs(unsigned int* dofs,

2991

const ufc::mesh& m,

2992

const ufc::cell& c) const

2993

{

2994

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

2995

dofs[1] = 4*c.entity_indices[3][0] + 1;

2996

dofs[2] = 4*c.entity_indices[3][0] + 2;

2997

dofs[3] = 4*c.entity_indices[3][0] + 3;

2998

}

2999

3000

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

3001

virtual void tabulate_facet_dofs(unsigned int* dofs,

3002

unsigned int facet) const

3003

{

3004

switch ( facet )

3005

{

3006

case 0:

3007

3008

break;

3009

case 1:

3010

3011

break;

3012

case 2:

3013

3014

break;

3015

case 3:

3016

3017

break;

3018

}

3019

}

3020

3021

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

3022

virtual void tabulate_entity_dofs(unsigned int* dofs,

3023

unsigned int d, unsigned int i) const

3024

{

3025

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

3026

}

3027

3028

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

3029

virtual void tabulate_coordinates(double** coordinates,

3030

const ufc::cell& c) const

3031

{

3032

const double * const * x = c.coordinates;

3033

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

3034

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

3035

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

3036

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

3037

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

3038

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

3039

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

3040

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

3041

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

3042

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

3043

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

3044

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

3045

}

3046

3047

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

3048

virtual unsigned int num_sub_dof_maps() const

3049

{

3050

return 1;

3051

}

3052

3053

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

3054

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

3055

{

3056

return new ffc_23_dof_map_0_2();

3057

}

3058

3059

};

3060

3061

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

3062

/// degrees of freedom (dofs).

3063

3064

class ffc_23_dof_map_0: public ufc::dof_map

3065

{

3066

private:

3067

3068

unsigned int __global_dimension;

3069

3070

public:

3071

3072

/// Constructor

3073

ffc_23_dof_map_0() : ufc::dof_map()

3074

{

3075

__global_dimension = 0;

3076

}

3077

3078

/// Destructor

3079

virtual ~ffc_23_dof_map_0()

3080

{

3081

// Do nothing

3082

}

3083

3084

/// Return a string identifying the dof map

3085

virtual const char* signature() const

3086

{

3087

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

3088

}

3089

3090

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

3091

virtual bool needs_mesh_entities(unsigned int d) const

3092

{

3093

switch ( d )

3094

{

3095

case 0:

3096

return false;

3097

break;

3098

case 1:

3099

return false;

3100

break;

3101

case 2:

3102

return false;

3103

break;

3104

case 3:

3105

return true;

3106

break;

3107

}

3108

return false;

3109

}

3110

3111

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

3112

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

3113

{

3114

__global_dimension = 12*m.num_entities[3];

3115

return false;

3116

}

3117

3118

/// Initialize dof map for given cell

3119

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

3120

const ufc::cell& c)

3121

{

3122

// Do nothing

3123

}

3124

3125

/// Finish initialization of dof map for cells

3126

virtual void init_cell_finalize()

3127

{

3128

// Do nothing

3129

}

3130

3131

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

3132

virtual unsigned int global_dimension() const

3133

{

3134

return __global_dimension;

3135

}

3136

3137

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

3138

virtual unsigned int local_dimension() const

3139

{

3140

return 12;

3141

}

3142

3143

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

3144

virtual unsigned int geometric_dimension() const

3145

{

3146

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

3147

}

3148

3149

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

3150

virtual unsigned int num_facet_dofs() const

3151

{

3152

return 0;

3153

}

3154

3155

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

3156

virtual unsigned int num_entity_dofs(unsigned int d) const

3157

{

3158

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

3159

}

3160

3161

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

3162

virtual void tabulate_dofs(unsigned int* dofs,

3163

const ufc::mesh& m,

3164

const ufc::cell& c) const

3165

{

3166

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

3167

dofs[1] = 4*c.entity_indices[3][0] + 1;

3168

dofs[2] = 4*c.entity_indices[3][0] + 2;

3169

dofs[3] = 4*c.entity_indices[3][0] + 3;

3170

unsigned int offset = 4*m.num_entities[3];

3171

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

3172

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

3173

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

3174

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

3175

offset = offset + 4*m.num_entities[3];

3176

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

3177

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

3178

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

3179

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

3180

}

3181

3182

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

3183

virtual void tabulate_facet_dofs(unsigned int* dofs,

3184

unsigned int facet) const

3185

{

3186

switch ( facet )

3187

{

3188

case 0:

3189

3190

break;

3191

case 1:

3192

3193

break;

3194

case 2:

3195

3196

break;

3197

case 3:

3198

3199

break;

3200

}

3201

}

3202

3203

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

3204

virtual void tabulate_entity_dofs(unsigned int* dofs,

3205

unsigned int d, unsigned int i) const

3206

{

3207

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

3208

}

3209

3210

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

3211

virtual void tabulate_coordinates(double** coordinates,

3212

const ufc::cell& c) const

3213

{

3214

const double * const * x = c.coordinates;

3215

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

3216

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

3217

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

3218

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

3219

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

3220

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

3221

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

3222

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

3223

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

3224

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

3225

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

3226

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

3227

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

3228

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

3229

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

3230

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

3231

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

3232

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

3233

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

3234

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

3235

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

3236

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

3237

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

3238

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

3239

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

3240

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

3241

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

3242

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

3243

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

3244

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

3245

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

3246

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

3247

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

3248

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

3249

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

3250

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

3251

}

3252

3253

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

3254

virtual unsigned int num_sub_dof_maps() const

3255

{

3256

return 3;

3257

}

3258

3259

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

3260

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

3261

{

3262

switch ( i )

3263

{

3264

case 0:

3265

return new ffc_23_dof_map_0_0();

3266

break;

3267

case 1:

3268

return new ffc_23_dof_map_0_1();

3269

break;

3270

case 2:

3271

return new ffc_23_dof_map_0_2();

3272

break;

3273

}

3274

return 0;

3275

}

3276

3277

};

3278

3279

#endif

Older »