Решить трехмерное уравнение Пуассона со смешанным граничным условием в фенике питона

Question

Я новичок в FEniCS. Я использую приложение Ubuntu на Win10, затем установил python3.6 и fenics.

Я пытаюсь решить трехмерное уравнение Пуассона, такое как Лаплас (u)=f, с граничным условием Дирихле (u=0) для боковых границ по x и y и граничным условием Неймана (∂ u⁄∂z=ug) для верха и низа границы, где данные rhs и граничного условия Неймана этого уравнения (f и ug) являются сеточными данными численного эксперимента. В этой числовой модели моя фактическая сетка имеет размер 301x301x22 с горизонтальным разрешением 2 км и неоднородным вертикальным разрешением в диапазоне от 0,08 км до 1 км.

Меня всегда смущала связь между сеткой в ​​фенике и сеткой в ​​моей настоящей проблеме. Разве сетка и интервал, которые я установил в фенике, должны быть такими же, как у моей настоящей проблемы? Если да, то как мне настроить неравномерную сетку в fenics (однажды я попробовал сопоставление сетки, как в моем примере)?

Упростите мою проблему, например, у меня в коде есть данные сетки f и ug, разрешение по горизонтали которых составляет 2 км, а по вертикали - Z, я пытаюсь создать сетку следующим образом:

from dolfin import *
import numpy as np

# Create mesh and define function space
mesh = BoxMesh(Point(0.0, 0.0, 0.0), Point(10.0, 10.0, 5.0), 5, 5, 5)
Z = [0.02666,  0.10259,  0.22883,  0.41162,  0.65375,  0.95918] #vertical levels of my actual grid data(unit:km)
zz=[]
for vertex in vertices(mesh):
    z = vertex.x(2)
    zz.append(Z[int(z)])
mesh.coordinates()[:,2] = np.array(zz)
V = FunctionSpace(mesh, 'Lagrange', 1)                                                                                                           

#Define boundary conditions
u0 = Constant(0.0)
tol = 1E-14   # tolerance for coordinate comparisons
def Dirichlet_boundary(x):
    return abs(x[0]) < tol or abs(x[0] - 10) < tol or abs(x[1]) < tol or abs(x[1] - 10) < tol
bc = DirichletBC(V, u0, Dirichlet_boundary)

#my input data----RHS and neumann boundary condition of u
f=np.arange(6*6*6).reshape(6,6,6)
ug=np.arange(24,6*6*6+24).reshape(6,6,6)

#Expand 3D data into 1D according to mesh coordinate
f_ve = Function(V) 
d2v = dof_to_vertex_map(V)   
f_value=f.flatten()
f_ve.vector()[:]= f_value[d2v]
g_ve = Function(V)   
ug_value=ug.flatten()
g_ve.vector()[:]= ug_value[d2v]

# Define variational problem and Compute solution
u = TrialFunction(V)
v = TestFunction(V)
a = inner(nabla_grad(u), nabla_grad(v))*dx          
L = g_ve*v*ds - f_ve*v*dx  
u = Function(V)
solve(a == L, u, bc, solver_parameters={"linear_solver": "cg","preconditioner": "amg"})
print(u(2,3,Z[3]))

Этот код может работать без ошибок. Но результат, похоже, не верен! Так есть ли проблемы с моим поколением мешей или другими?

Я буду очень признателен, если вы сможете мне помочь.

Жан Чан