dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

I just went through implementing Bernstein polynomials in Python
because Marie and I needed them for our experiments with adaptivity
and the FIAT basis functions could not be evaluated outside the
triangles (the mapping being singular along y = 1).

I've attached the code for the Bernstein polynomials. We will need to
implement this in C++ at some point when things are working.

--
Anders


On Wed, Sep 23, 2009 at 05:53:30PM -0500, Robert Kirby wrote:
> It probably should be added in dolfin/ufc.  Sometimes we want to play around
> with things that aren't supported by ffc yet.  Having a full-featured c++
> library (and python bindings) would help.  An example would be Berstein
> polynomials, which use a different evaluation method than in ffc/fiat.  Another
> example is the transformation business I'm working that needs cell jacobians
> and other geometric information.  I'd like to look at hand-done implementations
> using the dolfin library (preferably PyDolfin) to make sure stuff works before
> going through the labor of teaching ffc how to do it.
>
> Rob
>
> On Wed, Sep 23, 2009 at 5:45 PM, Garth N. Wells <gnw20@xxxxxxxxx> wrote:
>
>
>
>     Robert Kirby wrote:
>
>         Is there a method for getting the Jacobian of a mapping to a reference
>         cell in dolfin, or is this always generated by ffc?
>
>
>
>     It's generated by FFC. Could be a suggestion for an addition to UFC.
>
>     Garth
>
>
>
>
>         ------------------------------------------------------------------------
>
>         _______________________________________________
>         DOLFIN-dev mailing list
>         DOLFIN-dev@xxxxxxxxxx
>         http://www.fenics.org/mailman/listinfo/dolfin-dev
>
>
>
>
>

> _______________________________________________
> DOLFIN-dev mailing list
> DOLFIN-dev@xxxxxxxxxx
> http://www.fenics.org/mailman/listinfo/dolfin-dev

# Copyright and LGPL v3 or later etc, Marie Rognes and Anders Logg

from numpy import array, zeros, dot, sqrt

class Bernstein:

    def __init__(self, vertices, degree):
        "Create Bernstein polynomials of given degree on given simplex."

        # Store points and degree
        self.vertices = vertices
        self.degree = degree

        # Compute order tuples and scale factor
        self.order_tuples = _compute_order_tuples(degree, len(vertices))
        self.scale_factor = _factorial(degree)

    def __call__(self, i, x):
        "Evaluate polynomial i at point x."
        value = self.scale_factor
        for (j, r) in enumerate(self.order_tuples[i]):
            value *= _barycentric(self.vertices, j, x)**r / _factorial(r)
        return value

def _barycentric(vertices, i, x):
    "Evaluate barycentric coordinate i at point x."
    x = _inverse_affine_map(vertices, x)
    dim = len(vertices) - 1
    if i == 0:
        return 1.0 - sum(x)
    return x[i - 1]

# _inverse_affine_map copied from codesnippets.py in FFC

def _inverse_affine_map(vertices, x):
    "Map x to reference simplex."

    dim = len(vertices) - 1

    if dim == 1:

        # Compute Jacobian of affine map from reference cell
        J_00 = vertices[1][0] - vertices[0][0];

        # Get coordinates and map to the reference (UFC) element
        xx0 = (x[0] - vertices[0][0]) / J_00

        return array((xx0,))

    elif dim == 2:

        # Compute Jacobian of affine map from reference cell
        J_00 = vertices[1][0] - vertices[0][0]
        J_01 = vertices[2][0] - vertices[0][0]
        J_10 = vertices[1][1] - vertices[0][1]
        J_11 = vertices[2][1] - vertices[0][1]

        # Compute determinant of Jacobian
        detJ = J_00*J_11 - J_01*J_10

        # Compute inverse of Jacobian
        Jinv_00 =  J_11 / detJ
        Jinv_01 = -J_01 / detJ
        Jinv_10 = -J_10 / detJ
        Jinv_11 =  J_00 / detJ

        # Get coordinates and map to the reference (UFC) element
        xx1 = (vertices[0][1]*vertices[2][0] - vertices[0][0]*vertices[2][1] + J_11*x[0] - J_01*x[1]) / detJ
        xx0 = (vertices[1][1]*vertices[0][0] - vertices[1][0]*vertices[0][1] - J_10*x[0] + J_00*x[1]) / detJ

        return array((xx0, xx1))

    elif dim == 3:

        J_00 = vertices[1][0] - vertices[0][0]
        J_01 = vertices[2][0] - vertices[0][0]
        J_02 = vertices[3][0] - vertices[0][0]
        J_10 = vertices[1][1] - vertices[0][1]
        J_11 = vertices[2][1] - vertices[0][1]
        J_12 = vertices[3][1] - vertices[0][1]
        J_20 = vertices[1][2] - vertices[0][2]
        J_21 = vertices[2][2] - vertices[0][2]
        J_22 = vertices[3][2] - vertices[0][2]

        # Compute sub determinants
        d00 = J_11*J_22 - J_12*J_21
        d01 = J_12*J_20 - J_10*J_22
        d02 = J_10*J_21 - J_11*J_20
        d10 = J_02*J_21 - J_01*J_22
        d11 = J_00*J_22 - J_02*J_20
        d12 = J_01*J_20 - J_00*J_21
        d20 = J_01*J_12 - J_02*J_11
        d21 = J_02*J_10 - J_00*J_12
        d22 = J_00*J_11 - J_01*J_10

        # Compute determinant of Jacobian
        detJ = J_00*d00 + J_10*d10 + J_20*d20

        # Compute inverse of Jacobian
        Jinv_00 = d00 / detJ
        Jinv_01 = d10 / detJ
        Jinv_02 = d20 / detJ
        Jinv_10 = d01 / detJ
        Jinv_11 = d11 / detJ
        Jinv_12 = d21 / detJ
        Jinv_20 = d02 / detJ
        Jinv_21 = d12 / detJ
        Jinv_22 = d22 / detJ

        # Compute constants
        C0 = d00*(vertices[0][0] - vertices[2][0] - vertices[3][0]) \
           + d10*(vertices[0][1] - vertices[2][1] - vertices[3][1]) \
           + d20*(vertices[0][2] - vertices[2][2] - vertices[3][2])

        C1 = d01*(vertices[0][0] - vertices[1][0] - vertices[3][0]) \
           + d11*(vertices[0][1] - vertices[1][1] - vertices[3][1]) \
           + d21*(vertices[0][2] - vertices[1][2] - vertices[3][2])

        C2 = d02*(vertices[0][0] - vertices[1][0] - vertices[2][0]) \
           + d12*(vertices[0][1] - vertices[1][1] - vertices[2][1]) \
           + d22*(vertices[0][2] - vertices[1][2] - vertices[2][2])

        # Get coordinates and map to the UFC reference element
        xx0 = (C0 + d00*x[0] + d10*x[1] + d20*x[2]) / detJ
        xx1 = (C1 + d01*x[0] + d11*x[1] + d21*x[2]) / detJ
        xx2 = (C2 + d02*x[0] + d12*x[1] + d22*x[2]) / detJ

        return array((xx0, xx1, xx2))

    raise RuntimeError, "Unable to map point to reference element, only implemented in 1D, 2D, 3D."

# _compute_order_tuples and _factorial copied from Exterior

def _compute_order_tuples(k, n):
    "Compute all tuples of n integers such that the sum is k."
    if n == 1:
        return [(k,)]
    order_tuples = []
    for i in range(k + 1):
        for order_tuple in _compute_order_tuples(k - i, n - 1):
            order_tuples += [order_tuple + (i,)]
    return order_tuples

def _factorial(n):
    "Compute n!"
    if n == 0:
        return 1
    else:
        return n*_factorial(n - 1)

if __name__ == "__main__":

    # Compute some order tuples
    print "Order tuples:", _compute_order_tuples(2, 3)

    # Evaluate some barycentric coordinates
    vertices = [(0, 0), (1, 0), (0, 1)]
    print "lambda_0(0.5, 0.5) =", _barycentric(vertices, 0, (0.5, 0.5))
    print "lambda_1(0.5, 0.5) =", _barycentric(vertices, 1, (0.5, 0.5))
    print "lambda_2(0.5, 0.5) =", _barycentric(vertices, 2, (0.5, 0.5))

    # Create quadratic Bernstein polynomials on reference triangle
    q = 2
    b = Bernstein(vertices, q)

    # Evaluate at some points
    print "b_0(0.5, 0.5) =", b(0, (0.5, 0.5))
    print "b_1(0.5, 0,5) =", b(1, (0.5, 0.5))
    print "b_2(0.5, 0,5) =", b(2, (0.5, 0.5))
    print "b_3(0.5, 0.5) =", b(3, (0.5, 0.5))
    print "b_4(0.5, 0,5) =", b(4, (0.5, 0.5))
    print "b_5(0.5, 0,5) =", b(5, (0.5, 0.5))

Attachment: signature.asc
Description: Digital signature