System equivalent to 5 point finite difference stencil
5 months ago by
Consider the code below:
with the different parameters(Option1 and Option 2) available as in Functionspace to get the same stiffness matrix as for the 5 point stencil using the finite difference method ?
import numpy as np from fenics import * mesh = UnitSquareMesh.create(3, 3,CellType.Type_Option1) V = FunctionSpace(mesh, 'Option 2', 1) u = TrialFunction(V) v = TestFunction(V) a = (inner(grad(u), grad(v)))*dx matA = assemble(a)
The MIT matlab code here shows the possibility of obtaining such a matrix using Finite element methods. I would like to get this with Fenics. If I can obtain the same matrix with any of Fenics' capabilities as in the example MIT MATLAB code, that would be great as well.
Community: FEniCS Project
5 months ago by
The long answer is a bit more complicated. This uses a mesh of (n+1)^2 = 16 vertices, so the FE matrix will be 16 x 16. (n-1)^2 = 4 of these vertices are interior to the unit square, and the other 12 are on the boundary. The standard FD matrix is only 4 x 4 in this case: the values at the boundary mesh points are not considered as unknowns. In the FE matrix, 12 rows will be all zeros except for a one on the diagonal. These rows correspond to the boundary nodes and simply impose the boundary condition. If you cross out these rows and also the corresponding columns, you will get the FD matrix. Well, except for one thing: you may have to reorder the vertices (apply a permutation to the rows and columns) in order to make the FE matrix identical to the FD matrix.
from fenics import * n = 3 mesh = UnitSquareMesh(n, n) V = FunctionSpace(mesh, 'Lagrange', 1) u = TrialFunction(V) v = TestFunction(V) a = -inner(grad(u), grad(v))*dx bc = DirichletBC(V, 0., DomainBoundary()) A = assemble(a) bc.apply(A) print(A.array())
Please login to add an answer/comment or follow this question.