How should I implement Nitsche's method in this simple example?

General
1

Dear all,
I am struggling with implementing Nitsche’s method. Given the Poisson equation

\nabla^2 u = f \textrm{ in } \Omega
u = u_D \textrm{ in } \partial \Omega

where \Omega is the unit square, I can succesfully solve this by implementing the boundary condition with Nitsche’s method, here is a minimal working example:

from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt
import argparse
import numpy as np
from mshr import *
import matplotlib.pyplot as plt
import meshio
import ufl as ufl

i, j = ufl.indices(2)

parser = argparse.ArgumentParser()
args = parser.parse_args()

xdmffile_u = XDMFFile("u.xdmf")

# Create mesh
channel = Rectangle(Point(0, 0), Point(1.0, 1.0))
domain = channel
alpha = Constant(10.0)
h = Constant(2.0)
mesh = generate_mesh(domain, 16)
n = FacetNormal(mesh)

V = FunctionSpace(mesh, 'P', 8)
O = VectorFunctionSpace(mesh, 'P', 8, dim=2)

u = Function(V)
u_D = Function(V)
v = TestFunction(V)
f = Function(V)

class u_expression(UserExpression):
    def eval(self, values, x):
        values[0] =cos(x[0])**2 + 2*sin(x[0]+x[1])**4
        
    def value_shape(self):
        return (1,)
    
class laplacian_u_expression(UserExpression):
    def eval(self, values, x):
        # values[0] = 6.0
        values[0] = - 2 * (cos(2 * x[0]) - 4 * cos(2 * (x[0]+x[1])) + 4 * cos(4 * (x[0]+x[1])))
    def value_shape(self):
        return (1,)


f = interpolate(laplacian_u_expression(element=V.ufl_element()), V)
u_D = interpolate(u_expression(element=V.ufl_element()), V)

# Define variational problem
#this is the ordinary variational functional
F_0 = dot(grad(u), grad(v))*dx + f*v*dx + ( - n[i]*(u.dx(i))  * v ) * ds
#this is the term that enforces the BCs with Nitsche's method
F_BC =  ( ( + alpha/h*(u - u_D ) ) * v  + n[i]*(v.dx(i))*(u-u_D) ) * ds
F = F_0 + F_BC

solve(F == 0, u)

xdmffile_u.write(u, 0)

Notice that I added to the variational functional the terms \alpha/h (u - u_D) \, \nu \, ds and n^i \partial_i \nu (u - u_D) \, ds, where \alpha is the large-enough parameter in Nitsche’s method, h the radius of the smallest circle containing \Omega, and \nu the test function.

Now I would like to implement, with the same method, the BC

n^i \partial_i u = u_D \textrm{ on } \partial \Omega.

What are the terms that I should add to the variational functional to achieve this? And why?

Thank you !

2

That is a Neumann condition, and is typically imposed weakly by adding a term on the right hand side. There no need of using Nitsche for that.

3

I did not say there is one, this is a minimal example to understand how Nitsche’s would be imposed.

4

You wouldn’t use Nitsche for that kind of condition as it comes “natural” from the variational form. See for instance: Setting multiple Dirichlet, Neumann, and Robin conditions — FEniCSx tutorial

5

I understand, here is another example where I would like to impose a BC with Nitsche’s method, and (as far as I can tell) such BC cannot be imposed as a natural BC:

\sigma_i - \partial_i u = 0,
\sum_i (\sigma_i \,\partial_i \sigma_ i) = f
in \Omega, where I wrote the sum over the three repeated indices i explicitly for clarity.
With boundary conditions (BCs)
\sigma_i n^i = g \textrm{ on } \partial \Omega.

Here is a minimal (non) working example:

from dolfin import *
import ufl as ufl


# Create mesh
mesh = UnitSquareMesh(32, 32)
n = FacetNormal(mesh)

i = ufl.indices(1)

# Define function spaces
P_U = FiniteElement('P', triangle, 2)
P_SIGMA = VectorElement('P', triangle, 2)
element = MixedElement([P_U, P_SIGMA])
U_SIGMA = FunctionSpace(mesh, element)
U = U_SIGMA.sub(0).collapse()
SIGMA = U_SIGMA.sub(1).collapse()

# Define trial and test functions
u_sigma = Function(U_SIGMA)
u, sigma = split(u_sigma)
nu, rho = TestFunctions(U_SIGMA)
f = Function(U)

#analytical expression for a general scalar function
class f_Expression(UserExpression):
    def eval(self, values, x):
        values[0] = 1.0
    def value_shape(self):
        return (1,)
    
#define source function
f = interpolate(f_Expression(element=U.ufl_element()), U)

# Define variational form
F_sigma = (sigma[i]*rho[i] + u*((rho[i]).dx(i)) ) * dx - n[i]*rho[i]*u * ds
F_u = ( (sigma[0]*((sigma[0]).dx(0)) + sigma[1]*((sigma[1]).dx(1) )) * nu - f*nu ) * dx
F = F_sigma + F_u

# Compute solution
solve(F == 0, u_sigma)

If I run it, the Newton solver yields some Nans, which is normal because no BC has been imposed:

$ 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) = 1.795e-02 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
  Newton iteration 1: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)

My question is simple: what are the correct terms to add to F in order to impose the BC \sigma_i n^i = g \textrm{ on } \partial \Omega with Nitsche’s method?

Thank you