Hi all.
I have doubts on how to impose Neumann boundary conditions in my problem. I have an energy functional of \mathbf m(r,\theta) given by
F[\mathbf m]=\int_{\Omega} dx\, J(r,\theta) f_B(\mathbf m) + \int_{\partial\Omega_1} ds_1\, J(r,\theta) f_S(\mathbf m,\mathbf m_1) + \int_{\partial\Omega_2} ds_2\, J(r,\theta) f_S(\mathbf m,\mathbf m_2) \,,
where J(r,\theta)=r^2 \sin\theta is the Jacobian of the transformation to spherical coordinates, \mathbf m_1 and \mathbf m_2 are simply constant vectors and the surface integrals are performed in (different) regions 1 and 2 of the boundary.
As can be seen in the MWE below, I do dF = ufl.derivative(F, m, phi) to find the weak form and then define the non-linear problem fem.petsc.NonlinearProblem(dF, m, J=JF) which automatically includes Robin BCs in boundaries 1 and 2.
In this context, I would like to impose also Neumann BCs in the region 0 of the boundary, such that \partial_r\mathbf{m}=0 and \partial_{\theta}\mathbf{m}=0 for r=0. How can I incorporate this to my formulation of the problem?
Thanks in advance!
import numpy as np
import ufl
from dolfinx import mesh, fem, nls, log, plot, geometry
from petsc4py.PETSc import ScalarType
from mpi4py import MPI
mesh_comm = MPI.COMM_WORLD
OPdim = 3 # dimension of the OP
gdim = 2 # dimension of the geometry
Rbase = 1 # radius of the base (sphere)
def boundaryMeasure():
boundaries = [(0, lambda x: np.isclose(x[0], 0)), # origin
(1, lambda x: np.isclose(x[0], Rbase)), # base
(2, lambda x: np.isclose(x[1], minq)), # flat cap
(3, lambda x: np.isclose(x[1], np.pi))]
facet_indices, facet_markers = [], []
tdim = domain.topology.dim; fdim = tdim - 1
for (marker, locator) in boundaries:
facets = mesh.locate_entities(domain, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = mesh.meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
return ufl.measure.Measure("ds", domain=domain, subdomain_data=facet_tag)
def Jacob(x):
xr, xq = x
return xr**2*ufl.sin(xq)
def curl_sph(A, x):
"""Curl of a vector field depending only on r and theta in spherical coordinates."""
xr, xq = x; jr, jq = 0, 1
Ar, Aq, Af = A
curl_r = (ufl.Dx(Af*ufl.sin(xq),jq))/(xr*ufl.sin(xq))
curl_q = (-ufl.Dx(xr*Af,jr))/xr
curl_f = (ufl.Dx(xr*Aq,jr)-ufl.Dx(Ar,jq))/xr
return ufl.as_vector([curl_r, curl_q, curl_f])
def div_sph(A, x):
"""Divergence of a vector field depending only on r and theta in spherical coordinates."""
xr, xq = x; jr, jq = 0, 1
Ar, Aq, Af = A
div_r = ufl.Dx(xr**2*Ar,jr)/xr**2
div_q = ufl.Dx(Aq*ufl.sin(xq),jq)/(xr*ufl.sin(xq))
div_f = 0
return div_r + div_q + div_f
def fB(u, x):
fT = a/2.*ufl.dot(u,u) + c/4.*ufl.dot(u,u)**2
fD = d*( (div_sph(u,x))**2 + (ufl.dot(u,curl_sph(u,x)))**2 + (ufl.cross(u,curl_sph(u,x)))**2 )
return fT + fD
def fS(u, u0):
u0 = ufl.as_vector(u0)
return w/2*ufl.dot(u-u0,u-u0)
def ET(u, x):
u0base, u0cap = [1,0,0], [1,0,0]
ds = boundaryMeasure() # custom integration measure over the boundaries
return fB(u,x)*Jacob(x)*ufl.dx + fS(u,u0base)*Jacob(x)*ds(1) + fS(u,u0cap)*Jacob(x)*ds(2)
def solveProblem(prob):
solver = nls.petsc.NewtonSolver(mesh_comm, prob)
solver.convergence_criterion = "incremental"
#solver.ksp_type = "preonly"
solver.rtol = 1e-12
solver.max_it = 100
solver.report = True
minit = [1,0,0] # initial guess for the unknown function
minit = ufl.as_vector([fem.Constant(domain, ScalarType(minit[0])),fem.Constant(domain, ScalarType(minit[1])),
fem.Constant(domain, ScalarType(minit[2]))])
expr = fem.Expression(minit, V.element.interpolation_points())
m.interpolate(expr)
solver.nonzero_initial_guess = True
#log.set_log_level(log.LogLevel.INFO)
n, converged = solver.solve(m)
assert(converged)
print(f"Number of iterations: {n:d}")
#### physical parameters
a = 0.1; c = 0.01
d = 0.01
w = 1
#### define the domain of integration
minq = np.pi/2; bins = 10
domain = mesh.create_rectangle(mesh_comm, np.array([[0,minq],[Rbase,np.pi]]),[bins,int(bins*(np.pi-minq))])
x = ufl.SpatialCoordinate(domain)
#### define the space of funtions and the functions that we will use
element = ufl.VectorElement("CG", domain.ufl_cell(), 1, OPdim)
V = fem.FunctionSpace(domain, element) # space of vector functions with gdim components
m = fem.Function(V, name='m') # OPdim-dimensional vector OP
dm = ufl.TrialFunction(V) # variation of the OP
phi = ufl.TestFunction(V) # test function
#### derive the weak formulation directly from the energy functional
F = ET(m, x) # total energy
dF = ufl.derivative(F, m, phi) # variation of the energy
JF = ufl.derivative(dF, m, dm) # Jacobian of the energy
#### define and solve the problem
problem = fem.petsc.NonlinearProblem(dF, m, J=JF)
solveProblem(problem)