How to compute point-wise quantity along a boundary

Question

I am interested in obtaining point-wise tangential gradient of a FE solution along the boundary. Most FEniCS examples only show how to compute its integral as follows

ds = ds[markers]
n = FacetNormal(mesh)
L_integral = (n[0]*grad(p)[1]- n[1] * grad(p)[0])**2*ds(1)

where p is a FE solution.

For a demo of the integral along the boundary: demo_lift-drag.py (see link https://bitbucket.org/fenics-project/dolfin/src/master/demo/undocumented/lift-drag/python/demo_lift-drag.py?at=master)

How can I obtain a list of point-wise tangential gradient along the boundary? The length of the list should equal to the number of verticies along the boundary. Like this?

L_pointwise = (n[0]*grad(p)[1]- n[1] * grad(p)[0])**2

Answer 1

Let me answer my own question. The following code computes the tangential gradient along the boundary. It gives the tangential gradient for each facet. It assumes the element is a linear triangle. Every time the FE solution changes, only the last 4 lines need to be called to recompute the tangential gradient.

from numpy import *
from dolfin import *
import numpy as np
mesh = UnitSquareMesh(5,5)

def left_boundary(x, on_boundary):
    return on_boundary and \
    ( \
    x[0] < 0.5
    )

# Mark a CG1 Function with ones on part of the boundary
V = FunctionSpace(mesh, 'CG', 1)
bc = DirichletBC(V, 1, left_boundary)
u = Function(V)
u.vector()[:] = 0.
bc.apply(u.vector())

# Get vertices-DOF map
d2v = V.dofmap().dof_to_vertex_map(mesh)
v2d = V.dofmap().vertex_to_dof_map(mesh)

# Mark VertexFunction to check
vf = VertexFunction('size_t', mesh, 0)
vf.array()[:] = u.vector()[d2v]


D = mesh.topology().dim()
mesh.init(D-1,D) # build connection btw facets and cells
coords = mesh.coordinates()

# obtain the number of boundary facet/edges
n_be = 0 # number of boundary edges
for i in range(mesh.num_facets()):
    # for each facet (EDGE in 2D)
    cc = Facet(mesh,i)
    # check if it is boundary edge with one cell (FE element)
    if len(cc.entities(D)) == 2:# facet's adj cells
        continue
    n_be = n_be +1
print 'num of boundary facets', n_be
nB = zeros([n_be,3])
nDOF = arange(3*n_be).reshape(n_be,3)
nVIds =  arange(2*n_be).reshape(n_be,2)

# construct nB vector, fixed for each boundary facet
bei = -1
# for each facet (EDGE in 2D)
for i in range(mesh.num_facets()):  
    cc = Facet(mesh,i)
    # check if it is boundary edge with one cell (FE element)
    if len(cc.entities(D)) == 2:# facet's adj cells
    # interior facets: two adj cells
        continue
    bei = bei +1
    # obtain the boundary facet's cell's verticies
    cellVId =  mesh.cells()[cc.entities(D)]
    # coordinates of cell vertices
    v1= coords[cellVId[0][0]]
    v2= coords[cellVId[0][1]]
    v3= coords[cellVId[0][2]]
    # obtain facet (edge) normal
    n = cc.normal()
    # cell area
    A = ((v2[0]-v1[0])*(v3[1]-v1[1])-(v3[0]-v1[0])*(v2[1]-v1[1]))/2
    # linear triangle element's shape function's derivative
    B = array([[v2[1]-v3[1], v3[1]-v1[1], v1[1]-v2[1]],
        [v3[0]-v2[0], v1[0]-v3[0], v2[0]-v1[0]]])/2/A
    #obtain nB
    nB[bei][:] = n[0]*B[1]-n[1]*B[0]
    # obtain vertices index
    nDOF[bei] = d2v[cellVId[0]]
    nVIds[bei] = cc.entities(0)

# compute tangential gradient for all the boundary facets
for i in range(0, n_be):
    # obtain nodal values for a given boundary facet
    Us = u.vector()[nDOF[i]]
    tgrad = np.dot(nB[i], Us)
    print 'facet:', nVIds[i], 'tgrad:', tgrad


plot(vf)
interactive()

Comments

Hi xpq,

I am also interested in obtaining the nodal values of the tangent gradient of the solution.
However, in my case the solution u is a TestFunction(V) and apparently TestFunctions don't have 'vector' as attribute. Do you have an idea of how I could do it then?
For the moment I'm able to compute the integral of the tangent derivative along a boundary, i.e.

s0 = inner(grad(v('+')),n('+'))*dS(4)  + 1e-15*v*dx
S0 = assemble(s0)

where dS(4) is an internal boundary labeled with number 4 and n = FacetNormal(mesh).

Do you have any idea of how to get the nodal values in this case?

Thank you!

PS in the new version of FEniCS you need to substitute

# Get vertices-DOF map
d2v = V.dofmap().dof_to_vertex_map(mesh)
v2d = V.dofmap().vertex_to_dof_map(mesh)

with

# Get vertices-DOF map
d2v = dof_to_vertex_map(V)
v2d = vertex_to_dof_map(V)