delta functions?

[stevengj]

It would be nice to be able to represent Dirac delta distributions in order to solve a PDE with a delta-function right-hand-side.

Basically each delta-function term δ(x-x₀) would reduce the dimensionality of the integration (for assembly) by 1, but the implementation requires more understanding of the CellField mechanisms than we currently have.

cc @wenjie

[fverdugo]

Hi @stevengj!

yes, this is in our TO-DO list.

I am working in an enhanced notation for weak forms (see issue #265). Once this is available one could implement something like:

f(x) = x[1]
x0 = Point(1.4, 3.0)
l(v)  = ∫( f*v )*dΩ + 3.5*v(x0)

or

l(v)  = ∫( f*v  +  3.5*δ(x0) )*dΩ

[stevengj]

That sounds great. Ideally, we'd also want surface delta functions, or at least the simplest case thereof which is a delta function of just one or two of the cartesian coordinates. e.g.

l(v) = ∫ (f*v)*δ(x[1]-x0[1]) *dΩ

(For example, to implement mode-launching sources in wave equations.)

[fverdugo]

@santiagobadia @stevengj

I have added (PR #485) DiracDelta to impose point, line, or surface loads. The usage is pretty straightforward.

degree = 2
δ_point = DiracDelta{0}(model,tags=["mypoint1","mypoint2"])
δ_line = DiracDelta{1}(model,degree,tags="myline")
δ_surf = DiracDelta{2}(model,degree,tags="mysurf")

For point loads degree is optional. For line/surface loads degree is used to choose the quadrature rule. The implementation is efficient, since we relay on physical tags defined in the model (we do not need to look for which cell includes a point)

Once defined, it is used in the linear form together with other terms:

l(v) = δ_point(v) + # other terms

Example: linear elasticity with point + body forces:

using Gridap

domain = (0,1,0,1)
cells = (30,30)
model = CartesianDiscreteModel(domain,cells)

Ω = Triangulation(model)

order = 2
reffe = ReferenceFE(:Lagrangian,VectorValue{2,Float64},order)
V = FESpace(model,reffe,dirichlet_tags=[1,2,5])

dΩ = LebesgueMeasure(Ω,2*order)

δ = DiracDelta{0}(model,tags=4)

E = 1
ν = 0.33
λ = (E*ν)/((1+ν)*(1-2*ν))
μ = E/(2*(1+ν))
σ(ε) = λ*tr(ε)*one(ε) + 2*μ*ε

f = VectorValue(0,-10)
g = VectorValue(-1,0)

a(u,v) = ∫( ε(v) ⊙ (σ∘ε(u)) )dΩ
l(v) = ∫( f⋅v )dΩ + δ( g⋅v )

op = AffineFEOperator(a,l,V,V)
uh = solve(op)

writevtk(Ω,"results",cellfields=["uh"=>uh])

[stevengj]

Great! Note that there are also derivatives of delta functions ("dipoles", "quadrupoles", etcetera), which correspond to evaluating derivatives of the test function.

In principle, I guess you could implement delta-function derivatives using the same scheme…

[stevengj]

Note also that it would be nice to be able to use these in the bilinear form, not just for loads, e.g. to include a term δ(surface)*u*v in the bilinear form.

This shows up, for example, in modeling electromagnetic waves in graphene, which is treated as a delta-function susceptibility (coefficient).

[fverdugo]

a(u,v)=delta(u*v) should also work

[stevengj]

I also put a couple of other questions in #485:

1. Can you multiply by an amplitude, including variable amplitudes f(x) via f * δ?

2. What happens if the test function is discontinuous, e.g. for piecewise-constant spaces?

[fverdugo]

Hi @stevengj

1. You can define f(x) and then use it in l(v) = δ(f*v) (always inside δ)
2. DirachDelta assumes continuous functions. If you have discontinuous funcions, you can use the machinery associated with SkeletonTriangulation to impose surface loads, e.g., l(v) = ∫( mean(v) )dΓ where dΓ is the Lebesgue measure of a SkeletonTriangulation. (No current way to deal with point / line loads for discontinuous funcions. Do they make sense for discontinuous functions at all?)

[fverdugo]

BTW, you can even use the surface normal for surface loads:

n = get_normal_vector(get_triangulation(δ))
l(v) = δ(n·v)

get_tangent_vector is not implemented but should be easy to add it.

[fverdugo]

Internally get_triangulation(δ) returns a BoundaryTriangulation that maps to the first cell around the point/line/surface (that is why it only works for continuous functions)

[stevengj]

δ(f*v) is a bit of a weird syntax, mathematically, since it would normally imply a delta function where f*v == 0. It's fine, of course, to have a different syntax for computational purposes, but is there a reason why you can't support (f*v)*δ

Do they make sense for discontinuous functions at all?

Ordinarily one doesn't allow δ for discontinuous functions, which is why I was wondering what happens — is there an error, or does it pick the value from one side or the other, or …?

[santiagobadia]

Hi @fverdugo @stevengj

Let us focus on the point Dirac delta for the moment. I like the construction of the Dirac delta with your point entity

δ = DiracDelta{0}(model,tags=4)

Here, you are defining where you are positioning your Diract delta. The Dirac delta is a functional (over smooth functions), so I like

l(v) = δ( g⋅v )

In this case g⋅v is the function on which you are applying the functional. I like it more than the standard "integral" form.

However, I consider that we need another constructor that receives as input a point, as you said here. There are many problems in which you have sources/sinks within the domain that are not mesh nodes. E.g., think about the position of electrodes and measurements in electrical resistivity tomography. We consider adaptive octree meshes and this points are not nodes of our mesh.

With regard to applying Dirac Deltas to discontinuous functions, it is true that it does not make sense. But it does not make sense for H1 functions in d>=2 and you have done it in your test :-) . I think that, at the end of the day, you want a "cell-wise delta function". When the point is on a n-face of a cell, pick one, and only consider the delta function in the selected cell. This way, the result is the same if the solution is C0 and it also "works" for discontinuous functions. I guess you are already doing this, somehow, but I didn't have the time to check the code yet.

[fverdugo]

The Dirac delta is evaluated by taking the value of the given function in the first cell around the point/line/surface.

Here the user is at her own risk when using a DiracDelta. If discontinuous funcions are used, the code will provide a result and the user is responsable to interprete this result for her particular setup.

I believe that current implementation accounts for many practical cases and it is a good starting point. When using Gmsh, one can easily enforce that some desired points are represented by the mesh.

[stevengj]

I agree that you seem to have handled the main issue here, so I'm going to close this. I'll open up a new issue to track delta-function derivatives.