# 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
sphs_df = pd.read_csv(sph_file, header=None, sep='\s+', names=['l', 'm', 'glm', 'hlm'])
glm_array = np.zeros((61, 61))
hlm_array = np.zeros((61, 61))
for i in range(sphs_df.shape[0]):
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[2] / (x[0] ** 2 + x[1] ** 2 + x[2] ** 2) ** (1 / 2)
lon = ufl.acos(x[0] / ufl.sqrt(x[0] ** 2 + x[1] ** 2)) * ufl.sign(x[1]) + ufl.pi - ufl.sign(x[1]) * 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[223] = -1 * ((495179769008019818390136611716089140625 * sv_e0a[2]) / 2) + (9738535457157723095006020030416419765625 * pow(sv_e0a[2], 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.

Thanks a lot for replying! I upload the files to google drive.

The data distribution for the inner boundary (360*180): https://drive.google.com/file/d/1aGFLtr5VzdkmBaeMlFDE7D-JrelwLzJP/view?usp=drive_link

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

Could you please make your script complete, with all import statements and loading of the mesh.

Here’s my full code:

``````import numpy as np
import pyvista
import ufl
from dolfinx import fem, io, plot
from dolfinx.io.gmshio import read_from_msh
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[0] ** 2 + xyz[1] ** 2 + xyz[2] ** 2)**(1/2)

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

lat_carrington = np.pi / 2 - math.acos(xyz[2] / 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'):
sphs_df = pd.read_csv(sph_file, header=None, sep='\s+', names=['l', 'm', 'glm', 'hlm'])
glm_array = np.zeros((61, 61))
hlm_array = np.zeros((61, 61))
for i in range(sphs_df.shape[0]):
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[2] / (x[0] ** 2 + x[1] ** 2 + x[2] ** 2) ** (1 / 2)
lon = ufl.acos(x[0] / ufl.sqrt(x[0] ** 2 + x[1] ** 2)) * ufl.sign(x[1]) + ufl.pi - ufl.sign(x[1]) * 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[2] / (x[0] ** 2 + x[1] ** 2 + x[2] ** 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')