2

我是stackoverflow和fenics的新手。我最近使用 env 来计算另一个文件的输入代码。

我想计算有限元并保存组装矩阵，比如在矩形域D=(-1,1)*(-2,2)上。x我在和的方向上设置离散点的数量y，设置nx,ny=9，即2*10*10三角形。

当我计算矩阵并组装它们时，它们的大小(nx+1)*(ny+1)*(nx+1)*(ny+1)而不是(nx+1)*(ny+1).

MWE：

from dolfin import *
import numpy as np

# Mesh and function space
# mesh (generates 2*nx*ny number of triangles) and function space
# number of FEM = (nx+1)*(ny+1) = 100
# matrixSize = (nx+1)*(ny+1) * (nx+1)*(ny+1) = 100*100
nx = 9
ny = 9
# mesh dims
x = [-1, 1]
y = [-2, 2]
mesh = RectangleMesh(Point(x[0], y[0]), Point(x[1], y[1]), nx, ny)
V = FunctionSpace(mesh, "CG", 1)

# Time variables
dt = Constant(0.3)

g_expr = 'alpha'
g = Expression(g_expr , alpha=0.0, degree=2)

u0 = interpolate(g, V)

# Variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Expression('-exp(-(pow(x[0], 2)+pow(x[1], 2)))', degree=2)

a = u*v*dx + dt*inner(grad(u), grad(v))*dx
L = u0*v*dx + dt*f*v*dx
bc = DirichletBC(V, g, "on_boundary")

A = assemble(a)
print(A.array())
# print the sizes of row and col space of assembled matrix A (both = 100)
print((len(A.array()), len(A.array()[0])))

我会很感激任何帮助，谢谢！

4

0 回答 0