Instantly share code, notes, and snippets.

# Balaje/eigen-func-4.png

Last active October 22, 2021 03:17
Show Gist options
• Save Balaje/102485bb14ec6daf677f938fbd8f3ebb to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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ₕ)])

## 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:

### ThomasHaine commented Sep 22, 2021

Works great! Just what I was looking for.