# Setting values to Neumann boundary

Hi everyone. I’m trying to solve the Poisson problem in a spherical shell; the outer sphere is a Dirichlet boundary, and the inner sphere is a Neumann boundary.
I want to specify the du/dn distribution at the inner boundary. I tried spherical harmonics to construct the distribution. The function is (I use sympy because I cannot find legendre in ufl package):

``````def magmap_IB(x,sph_file='SPHs/mrmqc_c2074.dat'):
# read the spherical harmonic coefficients
glm_array = np.zeros((61, 61))
hlm_array = np.zeros((61, 61))
for i in range(sphs_df.shape):
glm_array[sphs_df['l'][i], sphs_df['m'][i]] = sphs_df['glm'][i]
hlm_array[sphs_df['l'][i], sphs_df['m'][i]] = sphs_df['hlm'][I]

import sympy as sp
cos_colat = x / (x ** 2 + x ** 2 + x ** 2) ** (1 / 2)
lon = ufl.acos(x / ufl.sqrt(x ** 2 + x ** 2)) * ufl.sign(x) + ufl.pi - ufl.sign(x) * ufl.pi

l_max = 17

COSCOLAT,LON = sp.symbols('cos_colat,lon')
Br = 0
from ufl import sqrt,sin,cos
Br = sp.symbols('0.')
for i_l in range(l_max + 1):
for i_m in range(i_l + 1):
print('l: ', i_l, '; m:', i_m)
Plm = sp.assoc_legendre(i_l, i_m, COSCOLAT)\
* (2 * math.factorial(i_l - i_m)
/ math.factorial(i_l + i_m)) ** (1 / 2)

Br = Br + Plm * (glm_array[i_l, i_m] * sp.cos(i_m * LON)
+ hlm_array[i_l, i_m] * sp.sin(i_m * LON))
magmap = eval(str(Br))
return magmap
``````

The main code is:

``````    V = fem.FunctionSpace(msh, ('CG', 3))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
x = ufl.SpatialCoordinate(msh)

# 边界自由度
inner_boundary_dof = fem.locate_dofs_topological(V=V, entity_dim=2, entities=inner_boundary)
outer_boundary_dof = fem.locate_dofs_topological(V=V, entity_dim=2, entities=outer_boundary)

phi0 = 0.
B0 = 100. * 1e-9
mu0 = 4 * np.pi * 1e-7

f = fem.Constant(msh, ScalarType(0))
# dirichlet B.C.
bc = fem.dirichletbc(value=ScalarType(phi0), dofs=outer_boundary_dof, V=V)

magmap = magmap_IB(x)

L = dot(f, v) * dx - dot(magmap, v) * ds

# Solve
ksp_type = 'gmres'
problem = fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options={'ksp_type': ksp_type})
uh = problem.solve()
``````

My problems are:
(1) The code cannot work with spherical harmonics over L=17 (the error is as follows). Are there any functions in dolfinx that can generate and handle with legendre polynomials?

``````libffcx_forms_cb3037c64ed6c464f7c84f3d41225a03793c8870.c:59428:26: error: integer literal is too large to be represented in any integer type
sv_e0a = -1 * ((495179769008019818390136611716089140625 * sv_e0a) / 2) + (9738535457157723095006020030416419765625 * pow(sv_e0a, 3)) / 2;
``````

(2) My plan B is to directly interpolate the distribution to the inner boundary. Is there a way in dolfinx to do that?

As you haven’t provided the CSV that you are using or the mesh, it is quite hard to give you any guidance.

The spherical harmonic coefficients (file used in function magmap_IB):

Here’s my full code:

``````import numpy as np
import pyvista
import ufl
from dolfinx import fem, io, plot
from ufl import ds, dx, grad, dot

import pandas as pd

from mpi4py import MPI
from petsc4py.PETSc import ScalarType
import math

def xyz2rlonlat(xyz, for_psi=False):
r_carrington = (xyz ** 2 + xyz ** 2 + xyz ** 2)**(1/2)

lon_carrington = math.asin(xyz / (xyz ** 2 + xyz ** 2)**(1/2))
if xyz < 0:
lon_carrington = np.pi - lon_carrington
if lon_carrington < 0:
lon_carrington += 2 * np.pi

lat_carrington = np.pi / 2 - math.acos(xyz / r_carrington)
if for_psi:
lat_carrington = np.pi / 2 - lat_carrington
return [r_carrington, lon_carrington, lat_carrington]

def magmap_IB(x,sph_file='SPHs/mrmqc_c2158.dat'):
glm_array = np.zeros((61, 61))
hlm_array = np.zeros((61, 61))
for i in range(sphs_df.shape):
glm_array[sphs_df['l'][i], sphs_df['m'][i]] = sphs_df['glm'][i]
hlm_array[sphs_df['l'][i], sphs_df['m'][i]] = sphs_df['hlm'][i]
import sympy as sp

cos_colat = x / (x ** 2 + x ** 2 + x ** 2) ** (1 / 2)
lon = ufl.acos(x / ufl.sqrt(x ** 2 + x ** 2)) * ufl.sign(x) + ufl.pi - ufl.sign(x) * ufl.pi

l_max = 60

COSCOLAT,LON = sp.symbols('cos_colat,lon')
Br = 0
from ufl import sqrt,sin,cos
Br = sp.symbols('0.')
for i_l in range(l_max + 1):
for i_m in range(i_l + 1):
print('l: ', i_l, '; m:', i_m)
Plm = sp.assoc_legendre(i_l, i_m, COSCOLAT)\
* (2 * math.factorial(i_l - i_m)
/ math.factorial(i_l + i_m)) ** (1 / 2)

Br = Br + Plm * (glm_array[i_l, i_m] * sp.cos(i_m * LON)
+ hlm_array[i_l, i_m] * sp.sin(i_m * LON))
magmap = eval(str(Br))
return magmap

def fem_laplace_in_shell(path,filename,inner_boundary_marker,outer_boundary_marker):
msh, cell_tags, facet_tags = read_from_msh(path + filename + '.msh', MPI.COMM_WORLD, 0, gdim=3)
inner_boundary = facet_tags.find(inner_boundary_marker)
outer_boundary = facet_tags.find(outer_boundary_marker)

# %% 创建函数空间，
V = fem.FunctionSpace(msh, ('CG', 3))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

x = ufl.SpatialCoordinate(msh)

inner_boundary_dof = fem.locate_dofs_topological(V=V, entity_dim=2, entities=inner_boundary)
outer_boundary_dof = fem.locate_dofs_topological(V=V, entity_dim=2, entities=outer_boundary)

phi0 = 0.
B0 = 100. * 1e-9
mu0 = 4 * np.pi * 1e-7

f = fem.Constant(msh, ScalarType(0))

bc = fem.dirichletbc(value=ScalarType(phi0), dofs=outer_boundary_dof, V=V)

# dipole field
# u_dipole = -B0 / mu0 * x / (x ** 2 + x ** 2 + x ** 2) ** (1 / 2)

magmap = magmap_IB(x)

L = dot(f, v) * dx - dot(magmap, v) * ds

# Solve
ksp_type = 'gmres'
from datetime import datetime
d1 = datetime.now()
problem = fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options={'ksp_type': ksp_type})
# problem.solver.view()
uh = problem.solve()
print('Solved')
d2 = datetime.now()
print(d2 - d1)

# SAVE
# with io.XDMFFile(msh.comm, "out_potentialfield/second_dipole_" + ksp_type + '_' + filename + ".xdmf", "w") as file:
#     file.write_mesh(msh)
#     file.write_function(uh)

cells, types, x = plot.create_vtk_mesh(V)
result = pyvista.UnstructuredGrid(cells, types, x)
result.point_data["u"] = uh.x.array.real
return result

if __name__ == '__main__':
path = 'MESH/PHYSICAL/'
filename = '(SphO1-SphI1)phy'

# Extract boundaries
inner_boundary_marker = 2
outer_boundary_marker = 1

result = fem_laplace_in_shell(path,filename,inner_boundary_marker,outer_boundary_marker)
# %% Visualization
result.set_active_scalars("u")
plotter = pyvista.Plotter()
plotter.add_mesh(result, show_edges=True, opacity=0.5, clim=[-0.1, 0.1], cmap='seismic')