How to compute the boundary normal for a 2D -> 3D mapped geometry

[SteinStoter]

Hi,

I'm getting the error message "The dimension of geometry must equal the sum of the dimensions of the given spaces" when attempting to perform an operation that I think should be valid. I've simplified my case to the below snippet:

from nutils.expression_v2 import Namespace
from numpy import pi, linspace

# Parameters
ns = Namespace()
ns.R = 100
ns.L = 100
ns.δ3 = function.eye(3)
ns.pi = pi

# Reference mesh in R2
verts = linspace(0,1,6)
topo, geom = mesh.rectilinear([verts, verts],periodic=(0,))

# Domain and field
ns.θ = geom
ns.define_for('θ', gradient='∇', normal='n', jacobians=('ds', 'dl'))
ns.basis = topo.basis('spline',degree=2)
ns.u = function.dotarg('udofs',ns.basis)

# Position vector after transformation in R3
ns.X_i = 'R cos( 2 pi θ_0 ) δ3_i0 + R sin( 2 pi θ_0 ) δ3_i1 + L θ_1 δ3_i2'
ns.define_for('X', gradient='∇X', normal='N', jacobians=('dV', 'dS','dL'))

# Covariant and contravariant bases
ns.G_ij  = '∇_j(X_i)'
ns.CG_ij = 'G_ki G_kj' # 2x2 Covariant surface metric coefficients
ns.CGI = function.inverse(ns.CG) # Contravariant surface metric coefficients
ns.Gcontra_ij = 'CGI_jk G_ik'

# Faulty integrals
bdyint = topo.boundary['top'].integral(f'∇_i( u ) Gcontra_ji N_j dL' @ns, degree=4)
# bdyint = topo.boundary['top'].integral(f'∇X_j( u ) N_j dL' @ns, degree=4) # Also gives error about having to invert a non-square Jacobian
bdyint.eval()

In this snipped, I'm just integrating the normal flux (in curvilinear coordinates we get ∇X_j( u ) = ∇_i( u ) Gcontra_ji). The issue appears to arise when the normal N_j is introduced in the integral.

In reality, the integrals are quite a bit more convoluted, so I'd really like to avoid having to do too many manual transformations.

Could someone help me understand what is going wrong here?

Best,
Stein

P.s. Am I correct to assume that the normal in define_for for the X coordinate is the normal to the boundary tangent to the transformed surface?

[gertjanvanzwieten]

The problem here is that Nutils does not know how the normal to the 1D boundary in 3D space is defined. Naturally, you want the one that is tangential to the manifold, but Nutils only sees a 1x3 coordinate transformation that it does not know how to handle. I'm not sure if we have an elegant solution in place for this situation, but the following (untested) workaround comes to mind:

ns.N = ('n_i ∇_i(X_j)' @ ns).normalized()

[gertjanvanzwieten]

I imagine this should work - does it? But in the meantime we looked into the current implementation of normal vectors, and realized it should be straightforward to generalize it in a way that it correctly evaluates to the desired normal in this situation as well, so that the code from your question will work without modification. I'll try to make a draft implementation soon.

[SteinStoter]

It seems like it is doing the trick. I'll mark the question as answered. Let me know when the update is committed!

[gertjanvanzwieten]

Will do!

[gertjanvanzwieten]

I have an implementation in branch "normal" (draft PR #768), could you try it out to see if your original code evaluates correctly?

[SteinStoter]

My code is built for v7.0, and I'm getting some errors trying to make it run on that v8.0 branch (initially about the utils module having been removed, and now something regarding evaluation). I'll email you the link to the log file if you'd like me to debug that.