至于事情的背景：

• 我正在使用专用 FEM 求解器的 FEniCS

• 我要解决的问题是控制 PDE 的静磁问题

  $$\bf{\nabla} \times \frac{1}{\mu} \bf{\nabla} \times \textbf{A} = \textbf{J}$$ 由于域是二维的，它简化为泊松方程： $$\frac{1}{\mu} \bf{\nabla}^{2} \it{A_z} = J_z$$ . 该域包含不同磁概率的子域：
  1. 空气和铜线（带电流源$J_z$) - 相对渗透率为 1
  2. 铁磁芯 - 10e3 的相对磁导率

此外，还有边界条件：$A_z = 0$在封装一切的空气子域的外边缘。

计算成功，我得到了磁矢量势的分布$A_z$. 接下来，我将磁通量分布评估为：

$$\bf{B} = \bf{\nabla} \times \bf{A}$$ .

现在，如果我们尝试将解决方案可视化为表面，这显然需要在单元格内进行插值，我们将获得以下结果：

我在可视化数据时遇到了这个问题。我知道解决方案与节点相关联。在位于边缘的节点上存在一些非零解，边缘是低渗透率域和高渗透率域之间的边界。例如，下一个节点的值为 0（零），但这些值正在单元内插值。我们知道，磁通量不会以这种方式表现，并且这些域之间的过渡应该是陡峭的。

节点处的解决方案似乎是正确的，但是当您尝试插入数据时会出现问题：

• 用于可视化目的
• 在任何点读取数据并将解决方案投影到另一个网格时（模拟多物理场时）

如何解决这个问题？我是否需要在铁磁材料（或与此相关的任何材料）的边缘添加一些额外的边界条件？我不能假设渗透率总是那么高，有时它可能是$\mu_r = 10$.

编辑：这是 Python 代码

from dolfin import *
import matplotlib.pyplot as plt
import numpy as np
import mshr
import copy
from scipy import constants
#############################################################################
import time
#############################################################################
from ProcessSubDomains import *
#############################################################################
mesh = Mesh('Mesh.xml')
# ---------------------------------------------------------------------------
markSubdomains('subdomains.xml', subdom_filename) # custom function that marks each subdomain with a number
materials = MeshFunction('size_t', mesh, mesh.topology().dim())
File(subdom_filename+'.xml') >> materials

#############################################################################
# MATERIAL PROPERTIES
class MaterialProperty(UserExpression):
    def __init__(self, materials, property, material_subdomains, property_index, **kwargs):
        super(MaterialProperty, self).__init__(**kwargs)    
        self.materials = materials
        self.property = property
        self.material_subdomains = material_subdomains
        self.property_index = property_index

    def eval_cell(self, values, x, cell):
        label = self.materials[cell.index]
        for key in self.material_subdomains:
            if label in self.material_subdomains[key]:
                values[0] = self.property[ self.property_index[key] ]

print("Assigning material properties...")
material_subdomains = {'air': [0, 2], 'iron': [1], 'wire': [3, 4, 5, 6]} # numbers associeted with subdomains
property_index = {'air': 0, 'iron': 1, 'wire': 2} # property index in permeability array
permeability = constants.mu_0*np.array([1, 35e3, 0.9991])
mu = MaterialProperty(materials, permeability, material_subdomains, property_index, degree=0)

#############################################################################
# DEFINE FUNCTION SPACE
V = FunctionSpace(mesh, 'P', 1)

#############################################################################
# BOUNDRY CONDITIONS
tol = 1e-14
def boundry(x, on_boundry):
    return on_boundry and ( near(abs(x[0]), 0.1, tol) or near(abs(x[1]), 0.1, tol) ) # checking when on Dirichlet BC
u_D = Constant(0) # value on Dirichlet BC
bc = DirichletBC(V, u_D, boundry)

#############################################################################
# REDEFINE INTEGRATION MEASURES
dx = Measure('dx', domain = mesh, subdomain_data = materials)

#############################################################################
# SOURCES
I = 400000.0
J_A = Constant(I)
J_B = Constant(I)

#############################################################################
# TRIAL AND TEST FUNCTIONS
A_z = TrialFunction(V)
v = TestFunction(V)

#############################################################################
# SOLVE VARIATIONAL PROBLEM
print("Solving variational form...")
a = (1/mu)*inner( grad(A_z), grad(v) )*dx
L = J_A*v*dx(3) - J_A*v*dx(4) + J_B*v*dx(5) - J_B*v*dx(6)

A_z = Function(V)
solve( a == L, A_z, bc)

#############################################################################
# POSTPROCESSING
W = VectorFunctionSpace(mesh, 'P', 1)

Bx = A_z.dx(1)
By = -A_z.dx(0)
B = project( as_vector(( Bx, By )), W )

B_abs = np.power( Bx**2 + By**2, 0.5 )

plot(A_z)
plt.show()

plot(B, mode='glyphs')
plt.show()

plot(B_abs)
plt.show()

---

在找到标量泊松方程的解之后 - 值$A_z$最好是创建相应的向量场$\textbf{A} = [0, 0, A_z]$，但以下操作失败：

A_vec = project( as_vector([0, 0, A_z]), W )

然后执行：

B_vec = project(curl(A_vec), W)

但是由于几何的二维性质（可能），这种操作似乎是不可能的。

编辑（解决方案）：

所以我找到了解决方案！在这里：FEniCS course 这需要再次求解 PDE（实际上它是特殊类型的映射）。假设您已经解决了方程并且有$A_z$.

必须在不连续伽辽金函数空间中求解的新公式是：

$$\int_{\Omega} w_h v \mathrm{d}x = \int_{\Omega} f v \mathrm{d}x$$

在哪里$f$是您要映射的先前解决方案，即$f = |\textbf{B}|$

# POSTPROCESSING
# calculate derivatives
Bx = A_z.dx(1)
By = -A_z.dx(0)

B_abs = np.power( Bx**2 + By**2, 0.5 ) # compute length of vector

# define new function space as Discontinuous Galerkin
abs_B = FunctionSpace(mesh, 'DG', 0)
f = B_abs # obtained solution is "source" for solving another PDE

# make new weak formulation
w_h = TrialFunction(abs_B)
v = TestFunction(abs_B)

a = w_h*v*dx
L = f*v*dx

w_h = Function(abs_B)
solve(a == L, w_h)

# plot the solution
plot(w_h)
plt.show()

---

补充说明：

@AntonMenshov：插值是在 ParaView 中自行完成的，也可以通过读取脚本中任意点的值来进行插值。

是的，可视化在商业软件 (Ansys Maxwell/Electronics) 中。不幸的是，除了解决了哪个 PDE 之外，文档没有说明它是如何完成的。第二个问题是，基本上不可能将数据（与节点关联的网格和解决方案）导出到软件之外。

网格被保存为.msh文件（Ansys 扩展）和解决方案为二进制文件。在查看.msh文件时，似乎该程序使用的是二阶元素而不是一阶元素，这是某种线索，为什么它会以这种方式运行。也许可以通过读取每个节点的位置、重建网格并读取这些节点上的值来完成。但这需要一些额外的劳动力，很多。

软件的另一个提示是它在每个子域之间设置默认的 Neumann 边界条件（如果没有另外说明），因此$\bf{B}$和正切分量$\bf{H}$在边界处是连续的。

---

@Paul：是的，我在这里使用的是 Continuous Galerkin 元素（拉格朗日）。我考虑在靠近边界的地方添加额外的元素层，这些元素的值会很快达到较小的值，这将减少后期处理中的这种发光。我还没有尝试过使用不连续 Galerkin 元素，这可能是一种选择，但需要调整变分公式。

---

当尝试沿某条线插入值以解决问题时$A_z$python/FEniCS 没有问题。但是当我试图插入甚至是$\bf{B}$这当然是作为导数计算的，即

B_x = A_z.dx(1)

当我试图绘制组件时，我得到：

Traceback (most recent call last):
...
some lines
...
/anaconda3/lib/python3.6/site-packages/dolfin/function/function.py", line 257, in ufl_evaluate
    assert derivatives == ()   # TODO: Handle derivatives

所以，对于 FEniCS 版本 2018.1.0 来说，这样的操作是不可能的。在 ParaView 或内置 FEniCSplot()功能中对其进行可视化时，这也会影响结果。

但是有一条出路，至少对于可视化部分。正确的解决方案是磁矢量势的分布$\bf{A}$. 然后我们可以计算$\nabla \times \bf{A}$要得到$\bf{B}$这要归功于 ParaView 中的计算器。

python 脚本中的插值仍然缺少此功能。也许可以通过添加自定义代码来执行类似于上面显示的操作。