-
Notifications
You must be signed in to change notification settings - Fork 97
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.
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
delta functions? #408
Comments
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.) |
I have added (PR #485) 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 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])
|
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… |
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 This shows up, for example, in modeling electromagnetic waves in graphene, which is treated as a delta-function susceptibility (coefficient). |
a(u,v)=delta(u*v) should also work |
I also put a couple of other questions in #485:
|
Hi @stevengj
|
BTW, you can even use the surface normal for surface loads: n = get_normal_vector(get_triangulation(δ))
l(v) = δ(n·v)
|
Internally |
Ordinarily one doesn't allow |
Let us focus on the point Dirac delta for the moment. I like the construction of the Dirac delta with your point entity
Here, you are defining where you are positioning your Diract delta. The Dirac delta is a functional (over smooth functions), so I like
In this case 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. |
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. |
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. |
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
The text was updated successfully, but these errors were encountered: