Problems with integration of boundary terms

[santiagobadia]

Hi @fverdugo

Sorry for raising more problems... but I am trying to integrate a boundary term and it does not work.

If you do the same for interior no problem. I have also considered to replace the integration by another FE function that results from solving the problem, and the problem persists.

The idea is to include a noise term in the projector we are working on. There are probably other ways to do it, but I have just considered a normal distribution for all dofs and create a FE function this way. Surely, I can create a random_fe_function(fespace) method but this is not the point of this discussion...

model = CartesianDiscreteModel(domain=(0.0,1.0,0.0,1.0), partition=(20,20))

labels = FaceLabels(model)
add_tag_from_tags!(labels,"diri_0",[1,3,7])
add_tag_from_tags!(labels,"diri_1",[2,4,8])

order = 2
diritags = ["diri_0", "diri_1"]
fespace = FESpace(reffe=:Lagrangian, conformity=:H1, valuetype=VectorValue{2,Float64},
model=model, labels=labels, order=order, diritags=["diri_0","diri_1"])

A = 0.01
x0 = A*randn(Float64,num_free_dofs(fespace))
x0D = A*randn(Float64,num_diri_dofs(fespace))
noise = FEFunction(U,x0,x0D)

btrian = BoundaryTriangulation(model)
bquad = CellQuadrature(btrian,degree=2*order)

eul2 = sqrt(sum( integrate(inner(noise,noise),btrian,bquad) ))

[santiagobadia]

I guess there is some kind of FE restriction method and it cannot be written automatically... I will look for it, please let me know if it does exist and which is its name

[fverdugo]

The noise fe function is aligned with the volume mesh. Before to integrate it on the boudnary, you have to restrict it, with bnoise = restrict(noise,btrian) or perhaps bnoise = restrict(btrian,nose) I don't remember exactly. Then, use bnoise to compute the integrals.

[santiagobadia]

I see that in all the tests, the right hand side seems to be an analytical function. Have you ever tried to include a term defined in terms of a FE function, e.g.,

function B_Ω(y)
    v, q = y
    inner(v,noise)
  end

This term does not work even at the bulk.

Using bnoise = restrict(btrian,nose) does work. Now I wanted to create a linear form on the RHS with inner(v,noise). Since you do not need to define explictly the restriction for the trial and test spaces, e.g.

function B_∂Ω(y)
    v, q = y
    + (γ/h) * inner(v,ud) - inner(∇(v), outer(nb,ud_cf)) + inner(q*nb,ud)
    # + inner(v,noise)
  end

or

function A_∂Ω(y,x)
    u, p = x
    v, q = y
    (γ/h) * inner(v,u) - inner(outer(nb,v), ∇(u)) - inner(∇(v), outer(nb,u)) + inner(v, p*nb) + inner(q*nb,u)
  end

it seems more consistent not to explicitly require to restrict boundary terms and just write

function B_∂Ω(y)
    v, q = y
   + inner(v,noise)
  end

Which is the expected usage you have in mind?

In any case, let us consider the bulk term first...

[fverdugo]

The term in the bulk should work...unless there is a bug.

For the boundaries, the current usage is to restrict all the cell fields that are captured by the function B_∂Ω(y). I don't see any easy way of doing these extra restrictions internally in a general way. It can be done, but it would require to attach a lot of metadata to the cell fields in order to decide to restrict or not in function of this metadata...

[santiagobadia]

For the boundary, is this the expected usage?

gdofs = interpolate_values(V,u)
gh = FEFunction(V,gdofs[1],gdofs[2])

bgh = restrict(btrian,gh)

b_∂Ω(v) = inner(v,bgh) 

It does not work. I want to create a FEFunction that I can integrate for boundary terms.

[fverdugo]

I am about to be ready with the new version for single-field problems.

It makes more sense to wait and see if the problem persist or not in the new version.

[santiagobadia]

Ok, so I guess that the usage is correct but there is a bug... We can wait

[fverdugo]

These problems should be fixed with v0.7.0. Closing issue.