@Stein and @dokken thank you for your reply. While I look into the literature that you cited, I gave a try at computing assemble((dot(v, n(v))**2)*ds). The value is indeed very small and, most importantly, it decreases as beta increases, so this is working.
What is more striking is "Nitsche’s’’ method seems to work even with the nonlinear constraint: In order to make the facet_normal markedly differ from n(v), I changed the matrix B in the code. Then I impose n(v).v=0 with “Nitsche’s” method, and I get that the BC is being enforced correctly: the solution is perpendicular to n(v) but not to facet_normal as it should:
from __future__ import print_function
from fenics import *
import ufl as ufl
from mshr import *
# Create mesh
disk = Circle( Point( 0, 0 ), 1 )
domain = disk
mesh = generate_mesh(domain, 64)
#radius of the smallest cell in the mesh
r_mesh = mesh.hmin()
#nitsche's parameter
beta = 1e2
xdmffile_v = XDMFFile( "v.xdmf" )
xdmffile_facet_normal = XDMFFile( "facet_normal.xdmf" )
xdmffile_n = XDMFFile( "n.xdmf" )
xdmffile_n_dot_v = XDMFFile( "n_dot_v.xdmf" )
i, j = ufl.indices( 2 )
Q = FunctionSpace( mesh, 'P', 2 )
V = VectorFunctionSpace(mesh, 'P', 2)
facet_normal = FacetNormal( mesh )
class grad_v0_expression( UserExpression ):
def eval(self, values, x):
values[0] = 2.0*x[0]
values[1] = 4.0*x[1]
def value_shape(self):
return (4,)
class laplacian_u_expression(UserExpression):
def eval(self, values, x):
values[0] = 6.0
values[1] = -2.0
def value_shape(self):
return (2,)
# Define variational problem
v = Function( V )
n_dot_v = Function( Q )
nu = TestFunction( V )
f = Function(V)
grad_v0 = Function( V )
f.interpolate(laplacian_u_expression(element=V.ufl_element()))
grad_v0.interpolate( grad_v0_expression( element=V.ufl_element() ) )
def calc_normal_cg2(mesh):
n = FacetNormal(mesh)
V = VectorFunctionSpace(mesh, "CG", 2)
u = TrialFunction(V)
v = TestFunction(V)
a = inner(u, v) * ds
l = inner(n, v) * ds
A = assemble(a, keep_diagonal=True)
L = assemble(l)
A.ident_zeros()
nh = Function(V)
solve(A, nh.vector(), L)
return nh
def facet_normal_smooth():
u = calc_normal_cg2(mesh)
return as_tensor(u[i], (i))
def B(v):
return as_tensor( [[1 + (v[0]) ** 2, (v[0]) * (v[1])], [(v[0]) * (v[1]), 2 + (v[1]) ** 2]] )
# return as_tensor( [[1 , 0], [0, 1]] )
def A(v):
return ufl.inv( B( v ) )
def n(u):
return as_tensor( A( u )[i, j] * facet_normal[j], (i) )
def n_smooth(u):
return as_tensor( A( u )[i, j] * facet_normal_smooth()[j], (i) )
F_0 = (((v[i]).dx( j )) * ((nu[i]).dx( j )) + f[i] * nu[i]) * dx - ((facet_normal[i] * grad_v0[i]) * nu[0] + facet_normal[i] * ((v[1]).dx( i )) * nu[1]) * ds
F_N = beta / r_mesh * (n( v )[j] * v[j]) * (n( v )[i] * nu[i]) * ds
F = F_0 + F_N
solve( F == 0, v )
xdmffile_v.write( v, 0 )
n_dot_v.assign( project( n_smooth( v )[i] * v[i], Q ) )
xdmffile_n_dot_v.write( n_dot_v, 0 )
xdmffile_facet_normal.write( facet_normal_smooth(), 0 )
xdmffile_n.write( project(n_smooth(v), V), 0 )
print("\int_{\partial \Omega} [n(v).v]^2 dS = ", assemble((dot(v, n(v))**2)*ds))
print("\int_{\partial \Omega} [facet_normal.v]^2 dS = ", assemble((dot(v, facet_normal)**2)*ds))
$ python3 example.py
No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 7.263e-01 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
[...]
Newton iteration 4: r (abs) = 3.661e-10 (tol = 1.000e-10) r (rel) = 5.041e-10 (tol = 1.000e-09)
Newton solver finished in 4 iterations and 4 linear solver iterations.
\int_{\partial \Omega} [n(v).v]^2 dS = 2.2367125330045466e-05
\int_{\partial \Omega} [facet_normal.v]^2 dS = 0.17642189719899645
This is also clear visually if I plot v, facet_normal and n(v), I clearly see that v is orthogonal to n(v) but not to facet_normal.
So it looks like “Nitsche’s” method is working also with a non-linear constraint, while you said that it is conceived for linear constraints only. Do you know why is it working ?
Thanks
