eigen-func-4.png

eigen.jl

using Gridap
using Gridap.Arrays
using Gridap.FESpaces
using Gridap.CellData
using Gridap.Algebra
using Gridap.CellData
using Arpack

struct EigenOperator{K<:AbstractMatrix,M<:AbstractMatrix} <: NonlinearOperator
  stima::K
  massma::M
end

struct EigenProblem <: FEOperator
  trial::FESpace
  test::FESpace
  op::EigenOperator
  nev::Int64
end

function EigenProblem(weakformₖ::Function, weakformₘ::Function,
                      test::FESpace, trial::FESpace, nev::Int64)
  L(v) = 0
  opK = AffineFEOperator(weakformₖ, L, test, trial)
  opM = AffineFEOperator(weakformₘ, L, test, trial)
  K = opK.op.matrix
  M = opM.op.matrix
  op = EigenOperator(K, M)
  EigenProblem(trial, test, op, nev)
end

function solve(prob::EigenProblem)
  K = prob.op.stima
  M = prob.op.massma
  ξ,Vec = eigs(K,M; nev=prob.nev,which=:SM,sigma=nothing)
  fₕs = Vector{CellField}(undef, prob.nev)
  for m=1:prob.nev
    fₕ = FEFunction(prob.trial, Vec[:,m])
    fₕs[m] = fₕ
  end
  ξ, fₕs
end

eigen_test.jl

include("eigen.jl")

using Gridap: ∇

  order = 1
domain = (0,π,0,π)
partition = (10,10)
model = CartesianDiscreteModel(domain, partition)
Ω = Triangulation(model)
dΩ = Measure(Ω, 4*order)

Vₕ = FESpace(model, ReferenceFE(lagrangian, Float64, order), conformity=:H1,
             vector_type=Vector{ComplexF64},
             dirichlet_tags="boundary")

Vₕ⁰ = TrialFESpace(Vₕ, 0)

a(u,v) = ∫(∇(u)⋅∇(v))dΩ
m(u,v) = ∫(u*v)dΩ

nev = 10
prob = EigenProblem(a, m, Vₕ⁰, Vₕ, nev)
ξ, uₕs = solve(prob)

# Visualize any eigenfunction
uₕ = uₕs[4]
writevtk(get_triangulation(uₕ), "sol", cellfields=["uh"=>real(uₕ)])

feature-list.md

Eigenvalue problems in Gridap

• New method EigenProblem to create a new standard generalized eigenvalue problem. K⋅x = ξ(M⋅x)
• Extend solve to return the eigenvalues ξ and an array of FEFunction corresponding to the eigenvectors.
• Dependencies: Arpack.jl for Eigenvalue problems.

Result

Run eigen-test.jl to solve:

-Δu = ξu   in Ω
u = 0  on ∂Ω

ξ=
10-element Vector{ComplexF64}:
  2.016502905927483 + 4.762186701117926e-17im
  5.141503239816802 - 1.0475743045319342e-16im
  5.141503239816826 - 5.058208907594383e-17im
  8.266503573706123 + 1.7471556284344306e-16im
 10.692073427500683 + 1.487357557157156e-16im
 10.692073427500699 - 4.462073229727142e-16im
 13.817073761390013 - 7.885029881682365e-16im
  13.81707376139005 - 1.5611691573305622e-17im
 19.200724573055854 + 1.5985650768400135e-16im
  19.20072457305593 + 9.243304317297491e-17im

4th eigenmode:

Works great! Just what I was looking for.