Plot spatial derivatives of the solution over a part of a boundary

1

I believe this simple question must have been asked at least a million times and yet I still fail to find the answer. What is the simplest way to plot some spatial partial derivative of the solution over a part of the boundary (specified by SubDomain)? The solution is from a mixed function space of degree higher than 1, but it is alright to get the values of the derivative only on the boundary mesh vertices. I am asking this because I want to plot the wall shear stress (which is just a specific type of flux) along the boundary. I want to get ordered vectors of coordinates of the specified part of the boundary, and of the derivative values at those vertices. Below are some of my recent failed attempts.

def LocalFrictionCoefficient(w,subdomain):
    # Inputs
    W = w.function_space()
    mesh = W.mesh()

    # Mark boundary belonging to sub domain
    plate_id = 7 
    ds = BoundaryIntegral(mesh, subdomain, plate_id)

    # Get function space for single velocity component
    V = W.sub(0).sub(0).collapse()
    v = TestFunction(V)
    #v2dof = vertex_to_dof_map(V) # this throws an error, so I replace it with:
    v2dof = V.dofmap().dofs()
    
    # Compute vector of local shear values
    shear = SurfaceShear(w)  # some function involving derivatives and boundary normal
    form = shear *v *ds(7) 
    vector = assemble(form)
    array = vector.get_local()[v2dof]

    return array

I also tried replacing

V = W.sub(0).sub(0).collapse()
v2dof = V.dofmap().dofs()

with

V = FunctionSpace(mesh, "CG", 1)
v2dof = vertex_to_dof_map(V)

but the order of the final array does not seem to match the order of the vertices of the mesh.

Thanks for any hints! I believe most of you are doing this routinely, so there must be some trivial solution. I cannot imagine a more basic task and yet I already spent a crazy amount of time on this.

2

For now, the only solution that works for me is to

  1. project the spatial derivatives of the solution to some function space on the whole domain,
  2. loop through all vertices of interest (those which lie on the boundary segment) and evaluate the derivatives at each vertex,
  3. order the vertices along the boundary segment,
  4. loop through the ordered vertices, compute the tangent and normal vectors manually from the neighboring vertices and use them to evaluate the function of interest.

Here is the code:

def LocalFrictionCoefficient(w,subdomain,mu):
    # Inputs
    W = w.function_space()
    mesh = W.mesh()
    fields = w.split()
    u,v = fields[0].split()

    # Project derivatives on boundary mesh
    bm = BoundaryMesh(mesh,"exterior")
    V = FunctionSpace(mesh, "CG", 1)
    dudx = project(u.dx(0), V)
    dudy = project(u.dx(1), V)
    dvdx = project(v.dx(0), V)
    dvdy = project(v.dx(1), V)
    
    # Get verticies of the boundary segment
    coords = bm.coordinates()
    X,Y,DUDX,DUDY,DVDX,DVDY = [],[],[],[],[],[]
    for iv,xv in enumerate(coords):
        if not subdomain.inside(xv, True): continue
        X.append(xv[0])
        Y.append(xv[1])
        DUDX.append(dudx(xv))
        DUDY.append(dudy(xv))
        DVDX.append(dvdx(xv))
        DVDY.append(dvdy(xv))
        
        # Sanity check
        if   np.isnan(X[-1]): raise Exception('NAN')
        elif np.isnan(Y[-1]): raise Exception('NAN')
        elif np.isnan(DUDX[-1]): raise Exception('NAN')
        elif np.isnan(DUDY[-1]): raise Exception('NAN')
        elif np.isnan(DVDX[-1]): raise Exception('NAN')
        elif np.isnan(DVDY[-1]): raise Exception('NAN')
    
    # Order the vertices
    X, isort = np.unique(X, return_index=True)
    Y        = np.take(Y, isort)
    DUDX, DUDY = np.take(DUDX,isort), np.take(DUDY,isort)
    DVDX, DVDY = np.take(DVDX,isort), np.take(DVDY,isort)
    
    # Get boundary tangent from ordered points
    npoints = len(X)
    shear = np.empty(npoints)
    for i in range(npoints):
        if   i==0        : t = [X[ 1]-X[ 0], Y[ 1]-Y[ 0] ]
        elif i==npoints-1: t = [X[-1]-X[-2], Y[-1]-Y[-2] ]
        else: t = [X[i+1]-X[i-1], Y[i+1]-Y[i-1]]
    
        t = np.divide(t, np.linalg.norm(t) ) # unit tangent vector
        n = [-t[1], t[0]] # boundary normal vector
        
        # Wall shear stress
        shear[i] = t[0]*n[0]*2*DUDX[i]
        shear[i]+= (t[0]*n[1]+t[1]*n[0]) * (DUDY[i]+DVDX[i])
        shear[i]+= t[1]*n[1]*2*DVDY[i]
        shear[i]*= mu
        
        # Sanity check
        if np.isnan(shear[i]): raise Exception('NAN')
    
    return X,Y,shear

It seems incredibly complicated and I guess I am wasting some memory when I store the derivatives on the entire domain. Does anyone have a simpler solution?

3

I guess you can use MeshView as you are using Legacy DOLFINx, and create a projection operator from the parent mesh to the sub mesh.

In DOLFINx, you could use dolfinx.fem.Expression (see for instance Implementation — FEniCSx tutorial) along with sub-meshes, as explained in
Plotting Partial Meshes - #4 by dokken

As interpolation operators in DOLFINx can be supplied with which cells to act on, you can identify which cells are connected to your boundary, and only interpolate over these.