Applying a Dirac Neumann BC on a variable lying on a boundary part in a multifield problem

[jfbarthelemy]

Here is the description of the problem

I consider a (normalized with no physical unit) electrical ohmic problem on a square with uniform conductivity such that the top line is made of a very high conductive material (infinite conductivity on a thickness tending to 0 such that the product conductivity by thickness is finite ). Thus the top line can be considered as a "lineic material" with its own potential field which is related to the potential of the same point in the square by a linear relationship existing between the local normal flux and the local difference of potentials. In addition I consider a Dirichlet boundary condition at the bottom line and an injection point (Neumann Dirac) at the top left point. At this point there are two choices of injection since two material points coincide on the same node (one of the lineic material and one of the bulk material)

1. If the Dirac is applied on the bulk material there is no problem (except an error message "MultiFieldFEFunction():ERROR: AssertionError: A check failed") which does not seem to hinder the calculation).
2. If the Dirac is applied on the lineic variablen then the following error is raised

ERROR: AssertionError:

Your are trying to integrate a CellField using a CellQuadrature defined on incompatible
triangulations. Verify that either the two objects are defined in the same triangulation
or that the triangulaiton of the CellField is the background triangulation of the CellQuadrature.

This problem is mathematically well-posed and correctly solved by other tools but I have problems to make it work with Gridap.

I attach an example (in .jl.txt to be renamed in .jl).

test_multi_dirac.jl.txt

[JordiManyer]

Hey @jfbarthelemy ,

As promised, here is a solution:

using Gridap
using Gridap.Geometry

D = 3
domain = repeat([0.0,1.0],D)
model = CartesianDiscreteModel(domain,(2,2,3))

face_tags = ["tag_23","tag_26","tag_24","tag_25"]
Γ  = BoundaryTriangulation(model,tags=face_tags)
_Λ  = SkeletonTriangulation(Γ)
glue = get_glue(_Λ,Val(D-1)) # glue from the edges to the attached faces

edge_tags = ["tag_17","tag_18","tag_20","tag_19"]
labeling = get_face_labeling(model)
edge_to_mask = get_face_mask(labeling,edge_tags,D-2)
edge_to_bgedge = findall(edge_to_mask) # ids of the edges we want

topo = get_grid_topology(model)
edge_to_face_map = Geometry.get_faces(topo,D-2,D-1) # edge to face connectivity

# Find local numbering of the edges we want
face_pairs = map(p -> [p...],zip(glue.plus.tface_to_mface,glue.minus.tface_to_mface))
local_edge_ids = map(edge_to_bgedge) do edge
  edge_faces = edge_to_face_map[edge]
  return findfirst(faces -> faces == edge_faces, face_pairs)
end
Λ = view(_Λ,local_edge_ids)

reffe = ReferenceFE(lagrangian,Float64,1)
V = FESpace(Γ,reffe)

dΛ = Measure(Λ,2)
v = interpolate(x -> sum(x),V)

l(v) = ∫(mean(v))*dΛ
b = assemble_vector(l,V)

Let me explain.... The problem is that if we create two separate triangulations directly from the model, Gridap has no way (yet) to change domain between them unless you can jump from one to the other using the model itself. This is currently unavailable for FACES -> CELLS -> EDGES, which is why it tells you triangulations are incompatible (although they technically aren't).

So the workaround is to create an edge triangulation that directly comes from the BoundaryTriangulation where your variable is defined. To do so, we create a SkeletonTriangulation (_Λ) of the boundary (Γ), which contains all the edges on that boundary. Now the problem is filtering out the edges we actually want.
To do so, we create a view of _Λ, which requires finding the local indices (within Γ) of the edges we want. This is the low-level part of the code. With that, we create the view Λ.

Now we can evaluate on Λ things defined on Γ, which is what you need. In my example, I take Γ as the four x-y sides of the a cube and Λ as the four vertical corners. But of course this works with any set of face_tags and edge_tags you want to use without any other modification (I hope).

Note that Λ is a SkeletonTriangulation. This is logical since every edge has two faces around it. If your variable is continuous (which I believe is your case) just use the mean like I did and you'll be fine.

For now you can use this code, I will discuss with the other devs how to integrate this functionality within Gridap in the near future.