PyOP2 Laplace demo with a strong boundary condition exemplification

laplace_pyop2.py

# Set up finite element problem. We do this by defining the finite element
# variational forms using the Unified Form Language (UFL)

E = FiniteElement("Lagrange", "triangle", 1)

v = TestFunction(E)
u = TrialFunction(E)
f = Coefficient(E)
g = Coefficient(E)

# LHS form
a = dot(grad(v,),grad(u))*dx
# RHS form
L = v*f*dx

# Generate code for Laplacian and rhs assembly. The compile_form function is
# provided by PyOP2, and it takes a variational form, and returns functions that
# compute the integral over the domain, domain boundary, and cell boundaries.
# Since we only have domain integrals, we can ignore the 2nd and 3rd returned
# values.

laplacian, _, _ = compile_form(a, "laplacian")
rhs, _, _  = compile_form(L, "rhs")

# Assemble matrix and rhs. This is done by invoking the functions returned by
# compile_form with a par_loop. When a par_loop is invoked, the arguments are:
# 1: The function to invoke (laplacian, or rhs) 
# 2: The iteration space. The iteration space consists of a set, and iteration
#    space dimensions. In the laplacian case, we are building a 3x3 local
#    matrix so the dimensions are (3,3). In the RHS case, the local vector has 3
#    elements so the dimensions are (3).
# 3: The local tensor. This takes mappings for the test and trial spaces. In
#    this case, the test and trial spaces are the same, and the mapping is
#    called elem_node - it maps from elements (which the loop iterates over) to
#    the local degrees of freedom. These are indexed by the special indices,
#    op2.i[], which refer to the current value of the i and j indices of the
#    local tensor. The access descriptor op2.INC tells PyOP2 that we're adding
#    into the tensor.
# 4: The coordinates. These are accessed through the elem_node mapping, but
#    don't have an index into them, which tells PyOP2 that it needs to pass the
#    coordinates for all vertices to the local assembly function.
# 5: Each of the functions (referred to as Coefficients in UFL parlance) used by
#    the form. The laplacian uses none, and the right hand side uses the
#    function f.

op2.par_loop(laplacian, elements(3,3),
             mat((elem_node[op2.i[0]], elem_node[op2.i[1]]), op2.INC),
             coords(elem_node, op2.READ))

op2.par_loop(rhs, elements(3),
                     b(elem_node[op2.i[0]], op2.INC),
                     coords(elem_node, op2.READ),
                     f(elem_node, op2.READ))

# We also want to apply a strong boundary condition. Weak boundary conditions
# can be applied by incorporating them into the finite element variational form.
# Strong boundary conditions don't have such a formulation, but this isn't a
# problem because PyOP2 allows the user to write functions in C that will be
# applied at each relevant node.

# This is similar to PETSc's row-zeroing, in that the diagonal term of the
# zeroed rows is set equal to the last parameter (1.0 in this case). The nodes
# on the boundary were 0, 1, 2, 6, 7, and 8.
mat.zero_rows([0, 1, 2, 6, 7, 8], 1.0)
# This is the kernel that sets the value of the target node to the given value.
strongbc_rhs = op2.Kernel("""
void strongbc_rhs(double *val, double *target) { *target = *val; }
""", "strongbc_rhs")
# Call the kernel. The bdry parameter contains the boundary values, which are
# applied to the RHS vector b that was previously assembled.
op2.par_loop(strongbc_rhs, bdry_nodes,
             bdry(op2.IdentityMap, op2.READ),
             b(bdry_node_node[0], op2.WRITE))

# Solve system. After calling op2.solve, x contains the solution to the
# problem.
op2.solve(mat, b, x)