Quad mesh analysis problem

mesh
1

Hi everyone,
I am trying to solve eigenvalue problem by using quad mesh elements in dolfin-x. But the results are not match together.

The proposed code is:

import numpy as np
from mpi4py import MPI
from dolfinx import Function, FunctionSpace, UnitSquareMesh, RectangleMesh
from dolfinx.fem import assemble_matrix, Form
from ufl import TrialFunction, TestFunction, grad, dx, dot
from slepc4py import SLEPc
from dolfinx.cpp.mesh import CellType

nev=1

mesh = RectangleMesh(
    MPI.COMM_WORLD,
    [np.array([0, 0, 0]), np.array([1, 1, 0])], [10, 10],
    cell_type=CellType.triangle)
V = FunctionSpace(mesh, ("CG", 1))

u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u), grad(v))*dx

A = assemble_matrix(a, [])
A.assemble()

eigensolver = SLEPc.EPS().create(MPI.COMM_WORLD)
eigensolver.setOperators(A)

eigensolver.solve()
for i in range (nev):
    l = eigensolver.getEigenpair(i)
    print(l.real)

and the result is:
7.822379233170955

When I change cell_type=CellType.triangle to cell_type=CellType.quadrilateral), the results is completely different:
3.890511304872971

Could you please guide me what is the reason of this error?
Thanks

2

Hi,

I want to update my question:

Previously I tried to compare the results of a simple eigenvalue problem with triangle and quadrilateral mesh elements in dolfin. The results were completely consistent.

Triangle mesh analysis in dolfin:

from dolfin import *

mesh = RectangleMesh(Point(0, 0), Point(1.0, 0.5), 16, 8)

V = FunctionSpace(mesh, "Lagrange", 1)
v = TestFunction(V)
u = TrialFunction(V)

s=inner(grad(u),grad(v))*dx

S = PETScMatrix()
assemble(s, tensor=S)

esolver = SLEPcEigenSolver(S)
esolver.solve()

cutoff = None
for i in range(100):
    (lr, lc) = esolver.get_eigenvalue(i)
        
    print(lr)

the output:

7.828767380299578
7.721888938729298
7.547871587566922
7.440194906508066
7.334226555416579

quadrilateral mesh analysis in dolfin:

from dolfin import *

mesh = RectangleMesh.create([Point(0, 0), Point(1, 1)],[16,8],CellType.Type.quadrilateral)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
ME = FunctionSpace(mesh, P1)
v = TestFunction(ME)
u = TrialFunction(ME)

s=inner(grad(u),grad(v))*dx

S = PETScMatrix()
assemble(s, tensor=S)

esolver = SLEPcEigenSolver(S)
esolver.solve()

cutoff = None
for i in range(100):
    (lr, lc) = esolver.get_eigenvalue(i)
        
    print(lr)

the output:

7.778949022072244
7.582136411376642
7.338799592271945
7.261579663802271
7.156952001449903

but similar comparison in dolfin-x doesn’t give acceptable answers.

Triangle mesh analysis in dolfin-x:

import numpy as np
import matplotlib.pyplot as plt
from dolfinx.plotting import plot
from mpi4py import MPI
from dolfinx import Function, FunctionSpace, UnitSquareMesh, RectangleMesh
from dolfinx.fem import assemble_matrix, Form
from ufl import TrialFunction, TestFunction, grad, dx, dot
from slepc4py import SLEPc
from dolfinx.cpp.mesh import CellType
import dolfinx
from dolfinx.io import XDMFFile

nev=5

mesh = RectangleMesh(
    MPI.COMM_WORLD,
    [np.array([0, 0, 0]), np.array([1, 1, 0])], [10, 10],
    cell_type=CellType.triangle)

V = FunctionSpace(mesh, ("CG", 1))

u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u), grad(v))*dx

A = assemble_matrix(a, [])
A.assemble()

eigensolver = SLEPc.EPS().create(MPI.COMM_WORLD)
eigensolver.setDimensions(nev)
eigensolver.setOperators(A)

eigensolver.solve()
vr, vi = A.createVecs()
for i in range (nev):
    l = eigensolver.getEigenpair(i ,vr, vi)
    print(l.real)

the output:

7.822379233170923
7.563558237348287
7.563558237348272
7.306318516404386
7.156141931378886

quadrilateral mesh analysis in dolfin-x:

import numpy as np
import matplotlib.pyplot as plt
from dolfinx.plotting import plot
from mpi4py import MPI
from dolfinx import Function, FunctionSpace, UnitSquareMesh, RectangleMesh
from dolfinx.fem import assemble_matrix, Form
from ufl import TrialFunction, TestFunction, grad, dx, dot
from slepc4py import SLEPc
from dolfinx.cpp.mesh import CellType
import dolfinx
from dolfinx.io import XDMFFile

nev=5

mesh = RectangleMesh(
    MPI.COMM_WORLD,
    [np.array([0, 0, 0]), np.array([1, 1, 0])], [10, 10],
    cell_type=CellType.quadrilateral)

V = FunctionSpace(mesh, ("CG", 1))

u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u), grad(v))*dx

A = assemble_matrix(a, [])
A.assemble()

eigensolver = SLEPc.EPS().create(MPI.COMM_WORLD)
eigensolver.setDimensions(nev)
#eigensolver.setWhichEigenpairs(2)
eigensolver.setOperators(A)

eigensolver.solve()
vr, vi = A.createVecs()
for i in range (nev):
    l = eigensolver.getEigenpair(i ,vr, vi)
    print(l.real)

the output:

3.890511304873001
3.8904289059861816
3.8074416364789347
3.807441636478916
3.6775942077819552

is the method of quadrilateral mesh analysis in dolfin-x incorrect? Introducing any source will be helpful in this issue.
Thanks in advace.

3

Note that the cause for the different eigenvaluesis that you are using different meshes for your different problems, causing inconsistency.
Using the same meshes for all problems gives you different eigenvalues for quadrilaterals and triangles, as they produce quite different stiffness matrices.

The eigenvalues are consistent between dolfinx and dolfin is you keep a consistent mesh.

You can verify this by using numpy to get the EigenValues:

import numpy as np
import matplotlib.pyplot as plt
from dolfinx.plotting import plot
from mpi4py import MPI
from dolfinx import Function, FunctionSpace, UnitSquareMesh, RectangleMesh
from dolfinx.fem import assemble_matrix, Form
from ufl import TrialFunction, TestFunction, grad, dx, dot
from slepc4py import SLEPc
from dolfinx.cpp.mesh import CellType
import dolfinx
from dolfinx.io import XDMFFile

nev= 5

mesh = RectangleMesh(
    MPI.COMM_WORLD,
    [np.array([0, 0, 0]), np.array([1, 1, 0])], [10 ,10],
    cell_type=CellType.quadrilateral)
V = FunctionSpace(mesh, ("CG", 1))

u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u), grad(v))*dx

A = assemble_matrix(a, [])
A.assemble()
A_numpy = A.convert("dense").getDenseArray()
w, v = np.linalg.eig(A_numpy)

eigensolver = SLEPc.EPS().create(MPI.COMM_WORLD)
eigensolver.setDimensions(nev)
eigensolver.setOperators(A)

eigensolver.solve()
vr, vi = A.createVecs()
eigen_values = np.zeros(nev)
for i in range (nev):
    l = eigensolver.getEigenpair(i ,vr, vi)
    assert(np.isclose(w-l.real, 0).any())

Here I have checked that the eigen solver returns an eigen value that is also returned by numpy (up to machine precision).

3 Likes
4

Perhaps you’re confusing eigenvalues of the operator A rather than that expected from the generalised eigenvalue problem A x = \lambda M x?

5

Hi dokken,
Thanks for your reply. you are right, I was wrong in mesh defining. Now I see the results of dolfin and dolfin-x are consistent. But the concept of differences in eigenvalues when using triangle and quadrilateral is still unclear. This is my first experience using quadrilateral elements and I think that I have to read more. Please let me know if you know of a reference for reading about quadrilateral analysis.

Thanks for your time and kind help.