Thank you for the idea: element_boundary would indeed work for the example I posted. Unfortunately, I lied about the boundary term (sorry
). We actually need
Code:
phi = 0.5*(v*(1-u) + v.Other(0)*(1-u.Other(0)))*n
phi += 0.5*(v-v.Other(0))
phiR = IfPos(n, n*IfPos(u-0.5, 0.25, u*(1-u)), 0)
a += SymbolicBFI(-v*(1-u)*grad(w))
a += SymbolicBFI(phi*(w - w.Other()), skeleton=True)
a += SymbolicBFI(phiR*w, BND, skeleton=True)
Am I wrong, or is there no way to restrict element_boundary integrals to inner facets? I guess I could subtract the wrong boundary contribution from element_boundary in a skeleton integrator....
Code:
phi = 0.5*(v*(1-u))*n
phi += 0.25*(v-v.Other(0))
phiR = IfPos(n, n*IfPos(u-0.5, 0.25, u*(1-u)), 0)
a += SymbolicBFI(-v*(1-u)*grad(w))
a += SymbolicBFI(phi*(w - w.Other()), element_boundary=True)
a += SymbolicBFI(-1*phi*w, BND, skeleton=True)
a += SymbolicBFI(phiR*w, BND, skeleton=True)
...but that seems very hacky.