Problem with interpolating tangent vectors defined on boundaries to surface

[steppenwolf]

I have a rectangle defined by mesh = RectangleMesh(Point(0.0, 0.0), Point(5.0, 0.1), 200, 200, "right/left"). The top boundary is marked as 3 and the bottom boundary is marked as 4. The problem is to find the unit tangent vectors to the top and the bottom boundaries and then interpolate those unit tangent vectors to the interior nodes which we will call the “interior tangent vectors”.

In the below MWE, I am able to successfully calculate the unit tangent vectors to the top and bottom boundaries. However, when I am trying to interpolate those to the interior nodes, there are two problems:

1. The direction of the interior tangent vectors are random. That is, some are pointing in the positive X direction while some are pointing in the negative X direction. But the tangent vectors at the top and the bottom boundaries are all oriented in the negative X direction.
2. The magnitude of the tangent vectors on the top and the bottom boundaries before interpolation are 1. But after interpolation, the magnitudes are no longer 1 on the surface as well as on the boundaries.

MWE:

from dolfin import *

# Creating Rectangular mesh #

mesh = RectangleMesh(Point(0.0, 0.0), Point(5.0, 0.1), 200, 200, "right/left")

# Defining the boundaries of the rectangle #

class bottom(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14
        return on_boundary and near(x[1], 0, tol)
    
class top(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14
        return on_boundary and near(x[1], 0.1, tol)
    
class left(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14
        return on_boundary and near(x[0], 0, tol)
    
class right(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14
        return on_boundary and near(x[0], 5.0, tol)
    
# Defining boundary markers, surface and boundary measures #
    
boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
surface_markers = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)

top().mark(boundary_markers, 3)
bottom().mark(boundary_markers, 4)
left().mark(boundary_markers, 1)
right().mark(boundary_markers, 2)

ds_custom = Measure("ds", domain=mesh, subdomain_data=boundary_markers) 
dx_custom = Measure("dx", domain=mesh, subdomain_data=surface_markers)

# Finding unit tangent vectors to the top and bottom boundaries #

n = FacetNormal(mesh)
t = as_vector([-n[1], n[0]])

V = VectorFunctionSpace(mesh, "CG", 1)

u = TrialFunction(V)
v = TestFunction(V)
a = inner(u,v)*ds_custom(3) + inner(u,v)*ds_custom(4) 
l = inner(t, v)*ds_custom(3) + inner(-t,v)*ds_custom(4)
A = assemble(a, keep_diagonal=True)
L = assemble(l)

A.ident_zeros()
tangent_boundary = Function(V)

solve(A, tangent_boundary.vector(), L)
tangent_boundary.rename('tangent_boundary','tangent_boundary')
File("tangent_boundary.pvd") << tangent_boundary

# Interpolating the unit tangent vectors on the top and bottom boundaries to the surface #

u1 = TrialFunction(V)
v1 = TestFunction(V)
a1 = inner(u1,v1)*dx_custom()
l1 = inner(tangent_boundary, v1)*ds_custom(3) + inner(tangent_boundary, v1)*ds_custom(4)
A1 = assemble(a1, keep_diagonal=True)
L1 = assemble(l)

A1.ident_zeros()
tangent_surface = Function(V)

solve(A1, tangent_surface.vector(), L1)
tangent_surface.rename('tangent_surface','tangent_surface')
File("tangent_surface.pvd") << tangent_surface

Results from Paraview:

1. Unit tangent vectors on the top and bottom boundaries (just showing a few glyphs):
   
   
2. Non-unit tangent vectors with different orientations when interpolated to the surface:
   
   

I want unit tangent vectors on the surface with all vectors oriented in the direction of the unit tangent vectors defined on the top and the bottom boundaries as shown in Result 1.

I am not sure what I am missing or doing wrong here.
Please advise.
Thanks in advance!

[hherlyng]

Hi steppenwolf,

In this post, I answered a question regarding computation of facet vectors, including tangential facet vectors. There I have linked to this gist which has dolfinx code that computes tangential facet vectors in such a way that all the tangent vectors are “oriented in the same way”, e.g. such that all have a positive x component.

That code also includes normalization of the tangent vectors, such that their length equal unity.

Let me know if this solves your problem.

Cheers, H

[steppenwolf]

@hherlyng
Thanks for your response.
The code in the link you provided calculates the tangent vectors to the facets, if I am not wrong. Can it be also used to interpolate the tangential vectors from the boundary nodes/vertices to the interior vertices? I need the tangents at the nodes and not the facets.
Please advise.

[hherlyng]

Yes, you’re right that the code I linked to approximates the tangent vectors of the facets.

Could you elaborate on why you need the tangent vectors at the vertices? In a point, such as a vertex or node on a finite element mesh, a tangent vector is not uniquely defined.

[steppenwolf]

@hherlyng

My actual problem involves a hyperelastic framework that requires a growth tensor to be defined at vertices of a mesh. As the body deforms, I need to apply a rotation matrix (defined pointwise on vertices) to the growth tensor (defined pointwise on vertices) to orient the growth tensor along the tangents defined at vertices so that the body grows in the direction of the deformation. So, I need the tangents defined at the interior vertices to be in the same direction as the tangents defined at nodes on the boundaries. So, I am trying to find the unit tangents on the top and bottom boundary nodes (which I am able to do) and interpolate those to the interior nodes such that the tangents at the interior vertices have the same direction/orientation as those on the top and bottom boundary nodes.

[hherlyng]

I’m still not sure I understand. In this drawing (please don’t judge), let’s say the node where the three arrows originate is the node in question. Asssume this is an interior node and that the mesh extends further in all directions. How would one choose any of the three (green, blue and red) tangent vectors as the “tangent vector of the node”?

Maybe you could provide an illustration of your problem?

Also, when you say your rotation matrix and growth tensor are defined pointwise on vertices, does that mean you actually have arrays with their values at the vertices? Or do you have finite element functions representing them, such that you have the degree of freedom/coefficient values at the vertices?

[steppenwolf]

@hherlyng
Here is a schematic diagram

Consider two vertices (red circles) on the boundaries 3 and 4. We find the tangent vectors at those vertices. Next, consider a set of vertices of the mesh that lies on a straight line connecting those two points. What I am aiming for is to define unit vectors at those vertices lying on the straight line joining the two opposite vertices on the boundaries such that vectors lying on the vertices on that line is kind of like a smooth interpolation of the tangents at the boundaries.

I think the word “tangent at interior nodes” is a bit misleading. I should have stated that we need to assign vectors to the interior nodes lying on a straight line connecting two vertices on the opposite boundaries such that the interior vectors provide a smooth interpolation along the straight line connecting the tangent vectors defined on the vertices of the boundaries.

I have matrices defined pointwise on the vertices of the mesh and not any finite element functions representing them.

Let me know if you have any questions.

[hherlyng]

All right, this was very helpful, I think I understand now. One possible solution would be the following.

If you have a mesh of this, approximate facet vectors \mathbf{v} of your mesh with the code I linked to. The vectors you seek in those nodes along the line connecting the boundary nodes, let’s denote those vectors \mathbf{v}_{a}. Let \mathbf{v}_{\perp} be the orthogonal projection of \mathbf{v} onto the line connecting the boundary nodes. Then \mathbf{v}_a = \mathbf{v} - \mathbf{v}_{\perp} should be the vectors you are seeking. This projection is pretty straightforward to perform in 2D. A way to do this is to define a vector \mathbf{u} that goes from one boundary node to the other, and then performing Gram-Schmidt orthogonalization of \mathbf{v} with respect to \mathbf{u}. You would presumably need to normalize the vector \mathbf{v}_a again after performing this projection, if you want them to be of unit length.

In case you have a 3D mesh, you would have to project the tangent vectors into the plane spanned by either of the boundary tangent vectors and the straight line connecting the two boundary nodes.

Let me know if this was of help.

[steppenwolf]

@hherlyng
Thanks for the suggestion. That will work I guess.

@hherlyng

Let \mathbf{v}{\perp}v⊥\mathbf{v}{\perp} be the orthogonal projection of \mathbf{v}v\mathbf{v} onto the line connecting the boundary nodes.

Also, I don’t have a physical line connecting the vertices on the boundaries.

I was wondering if there is any way to modify the MWE I provided in my initial post to get it done. I don’t understand why that MWE isn’t working. In the MWE, I already have the unit tangents at the nodes on the boundaries. I am not sure why the interpolation section of that MWE isn’t working.

[nami]

@steppenwolf, I believe your picture might be misleading. I understand that your growth must happen along the length of the medium. If that is correct, then the two points in your picture must not align like that, rather the bottom point should be slightly moved to the right so that the inner vectors are aligned along the longitudinal axis of the domain.
That being said, you can adopt @hherlyng method or do a weighted average of the two boundary vectors. The weights should be the shortest distance of the inner point from each boundary.

[hherlyng]

@steppenwolf I made a few changes to your MEW that produces normalized tangent vectors in the interior of domain, here is a clip of the plot:

Here is the second part of your code with the changes I made:

v1 = TestFunction(V)
a1 = inner(u1('+'), v1('+'))*dS_custom() + inner(u1('-'), v1('-'))*dS_custom()
l1 = inner(t, v)*ds_custom(3) + inner(-t,v)*ds_custom(4)
A1 = assemble(a1, keep_diagonal=True)
L1 = assemble(l1)

A1.ident_zeros()
tangent_surface = Function(V)

solve(A1, tangent_surface.vector(), L1)
tangent_norm = sqrt(inner(tangent_surface, tangent_surface))
normed_vec   = as_vector((abs(tangent_surface[0])/tangent_norm, tangent_surface[1]/tangent_norm))
tangents_normed = project(normed_vec, V)
File("tangents_normed.pvd") << tangents_normed

The integral measure in a1 is defined as dS_custom = Measure("dS", domain=mesh, subdomain_data=surface_markers). I made this change as I thought it made more sense to solve for u1 on the interior facets, and not the interior cell volumes (which was the case in your code since you used a dx measure). However, it seems like they produce the same results for the current element choice.

Also, I just used the tangent vector t directly instead of the vectors tangent_boundary approximated in the first variational problem of your code.

When normalizing the vectors to have unit length, I have also added taking the absolute value of the x component to make all the vectors align in the same way. If you want your tangent vectors to be aligned in some other direction, you would have to change this.

Hope this is of help:)

[steppenwolf]

@hherlyng
Thanks for the explanation and the code. Now I understand what I did wrong.
Really appreciate it!