Saya telah bekerja sedikit pada penyempurnaan adaptif untuk melakukan pekerjaan (lihat kode di bawah). Penskalaan indikator kesalahan dengan ukuran mesh total dan variasi total fungsi mesh tidak sempurna, tetapi Anda dapat menyesuaikannya dengan kebutuhan Anda. Gambar di bawah ini untuk testcase # 4. Jumlah sel meningkat dari 200 menjadi sekitar 24.000, yang mungkin sedikit di atas, tetapi hasilnya cukup bagus. Jala menunjukkan bahwa hanya bagian yang relevan telah disempurnakan. Artefak yang masih bisa Anda lihat, adalah elemen-elemen urutan ketiga itu sendiri yang tidak dapat mewakili cukup akurat.
from dolfin import *
from numpy import abs
def compute_error(expr, mesh):
DG = FunctionSpace(mesh, "DG", 0)
e = project(expr, DG)
err = abs(e.vector().array())
dofmap = DG.dofmap()
return err, dofmap
def refine_by_bool_array(mesh, to_mark, dofmap):
cell_markers = CellFunction("bool", mesh)
cell_markers.set_all(False)
n = 0
for cell in cells(mesh):
index = dofmap.cell_dofs(cell.index())[0]
if to_mark[index]:
cell_markers[cell] = True
n += 1
mesh = refine(mesh, cell_markers)
return mesh, n
def adapt_mesh(f, mesh, max_err=0.001, exp=0):
V = FunctionSpace(mesh, "CG", 1)
while True:
fi = interpolate(f, V)
v = CellVolume(mesh)
expr = v**exp * abs(f-fi)
err, dofmap = compute_error(expr, mesh)
to_mark = (err>max_err)
mesh, n = refine_by_bool_array(mesh, to_mark, dofmap)
if not n:
break
V = FunctionSpace(mesh, "CG", 1)
return fi, mesh
def show_testcase(i, p, N, fac, title1="", title2=""):
funcs = ["sin(60*(x[0]-0.5)*(x[1]-0.5))",
"sin(10*(x[0]-0.5)*(x[1]-0.5))",
"sin(10*(x[0]-0.5))*sin(pow(3*(x[1]-0.05),2))"]
mesh = UnitSquareMesh(N, N)
U = FunctionSpace(mesh, "CG", p)
f = interpolate(Expression(funcs[i]), U)
v0 = (1.0/N) ** 2;
exp = 1
#exp = 0
fac2 = (v0/100)**exp
max_err = fac * fac2
#print v0, fac, exp, fac2, max_err
g, mesh2 = adapt_mesh(f, mesh, max_err=max_err, exp=exp)
plot(mesh, title=title1 + " (mesh)")
plot(f, title=title1)
plot(mesh2, title=title2 + " (mesh)")
plot(g, title=title2)
interactive()
if __name__ == "__main__":
N = 10
fac = 0.01
show_testcase(0, 1, 10, fac, "degree 1 - orig", "degree 1 - refined (no change)")
show_testcase(0, 2, 10, fac, "degree 2 - orig", "degree 2 - refined")
show_testcase(0, 3, 10, fac, "degree 3 - orig", "degree 3 - refined")
show_testcase(0, 3, 10, 0.2*fac, "degree 3 - orig", "degree 3 - more refined")
show_testcase(1, 2, 10, fac, "smooth: degree 2 - orig", "smooth: degree 2 - refined")
show_testcase(1, 3, 10, fac, "smooth: degree 3 - orig", "smooth: degree 3 - refined")
show_testcase(2, 2, 10, fac, "bumps: degree 2 - orig", "bumps: degree 2 - refined")
show_testcase(2, 3, 10, fac, "bumps: degree 3 - orig", "bumps: degree 3 - refined")