Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Difficulty with complex vector-valued interior facet weak forms #931

Closed
kroppert opened this issue Aug 9, 2023 · 2 comments
Closed

Difficulty with complex vector-valued interior facet weak forms #931

kroppert opened this issue Aug 9, 2023 · 2 comments

Comments

@kroppert
Copy link

kroppert commented Aug 9, 2023

Dear all,

I am trying to define an IPDG (interior penalty discontinuous Galerkin) version of time-harmonic Maxwell's equations using Nedelec reference elements. The setup I wanted to test it on is just a cube with a (tangential vector) Dirichlet excitation on one side of the cube (I know it doesn't make much sense but I wanted to keep it simple).

This is the setup of the geometry, the definition of the measures and the volume bilinear form:

using Gridap
using Gridap.Geometry

L = 1.0
domain = (0.0, 1.0, 0.0, 1.0, 0.0, 1.0)
partition = (4,4,4)
model = CartesianDiscreteModel(domain,partition)
labels = get_face_labeling(model)
add_tag_from_tags!(labels,"S_x_pos",[1,3,5,7,13,15,17,19,25])

order = 2
V = TestFESpace(model, ReferenceFE(nedelec, Float64, order), conformity=:L2, dirichlet_tags="S_x_pos", vector_type=Vector{ComplexF64})
U = TrialFESpace(V, VectorValue(0.0+0.0im, 0.0+0.0im, 1.0+0.0im))

Ω = Triangulation(model)
Γ = BoundaryTriangulation(model)
Λ = SkeletonTriangulation(model)
degree = 2*order
dΩ = Measure(Ω,degree)
dΓ = Measure(Γ,degree)
dΛ = Measure(Λ,degree)
n_Γ = get_normal_vector(Γ)
n_Λ = get_normal_vector(Λ)

a_Ω(u,v) = ∫( curl(v) ⋅ curl(u) )dΩ - ∫( v ⋅ u )dΩ

However, I am running into a few problems, which I hope, are just due to my inability to define them properly :)

  1. When using the mean and jump operators (which are involving a cross product in this space) , defined as
# integrals on the interior facets
 a_Λ(u,v) = ∫( - mean(curl(u)) ⋅ jump(n_Λ × v)
              - mean(curl(v)) ⋅ jump(n_Λ × u)
              + jump(n_Λ × u) ⋅ jump(n_Λ × v) )dΛ
a(u,v) = a_Ω(u,v) + a_Λ(u,v) 
l(v) = 0
op = AffineFEOperator(a, l, U, V)
uh = solve(op)

I get the error message:

ERROR: MethodError: no method matching cross(::SkeletonPair{Gridap.CellData.GenericCellField{ReferenceDomain}, Gridap.CellData.GenericCellField{ReferenceDomain}},
::Gridap.FESpaces.SingleFieldFEBasis{Gridap.FESpaces.TestBasis, ReferenceDomain})

Then I found this issue #793 and I guess I can define the operators on my own - so the cross product should not be an issue. Then I defined a very simple integral on the interior facets to see if it works in principle

a_Λ(u,v) = ∫( u.plus ⋅ v.plus)dΛ              
a(u,v) = a_Ω(u,v) + a_Λ(u,v) 
l(v) = 0
op = AffineFEOperator(a, l, U, V)
uh = solve(op)

but then I get the following error:

ERROR: MethodError: Cannot `convert` an object of type 
   Gridap.Fields.ArrayBlock{Array{ComplexF64,1},1} to an object of type 
   Gridap.Fields.ArrayBlock{Array{Float64,1},1}

Can anyone see an issue in my definitions or is this functionality currently not implemented?

  1. When skipping the interior facets bilinear form and just solving a continuous time-harmonic curl-curl problem, I get a solution without running into an error but then writing it to a vtk file is cumbersome.
writevtk(Ω,"results",cellfields=["uh"=>uh])
# results in
# ERROR: ArgumentError: data type not supported by VTK: ComplexF64

and

writevtk(Ω,"results",cellfields=["uh_r"=>real(uh)])
# results in 
# ERROR: MethodError: no method matching real(::VectorValue{3, ComplexF64})

Is there any way to write a complex vector-valued result to a vtk file?

Sorry for the long question but I would be very grateful if someone here could help me out. I have also attached the script as a txt file, maybe this is easier to handle than the code snippets above.
IPDG.txt
@HelgeGehring
Copy link
Contributor

Hey kroppert,

I recently had the same problem with writing complex VectorValues to .vtk-files. I made a pull request (#934) to fix this.

Until it's merged you can add

Base.real(x::VectorValue{D,ComplexF64}) where {D} = VectorValue(real.(x.data))
Base.imag(x::VectorValue{D,ComplexF64}) where {D} = VectorValue(imag.(x.data))

to your script :)

@santiagobadia
Copy link
Member

Hi @kroppert @HelgeGehring

This PR has been merged.

I close the issue.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
3 participants
@santiagobadia @HelgeGehring @kroppert