Integral in a weak form

Question

Hi,
I am modeling a waveguide in 3D in frequency domain using Nedelec finite elements.
The problem is that I can't assemble the following form:

$$\int_S v \cdot w \cdot \left(\int_S w \cdot u\,dS \right)\,dS$$

where:

class MyExpression(Expression):
    ...
    def eval_cell(self, values, x, ufc_cell):
        ...
    def value_shape(self):
        return (3,)

V = FunctionSpace(mesh, "N2curl", 1)
u = TrialFunction(V)
v = TestFunction(V)
w= MyExpression()

I've found http://fenicsproject.org/qa/7535/use-integral-coefficient-in-ufl-formulation which seems to describe what I want to achieve, but using UFL's Integral object lead me to:

a0 = ufl.Integral(inner(v,w), ds.integral_type(), ds.domain(), ds.subdomain_id(), ds.metadata(), ds.subdomain_data())
assemble(inner(v,w)*a0*ds)
An integral can only be multiplied by a globally constant scalar expression.
*** UFLException: An integral can only be multiplied by a globally constant scalar expression.

I've been reading about UFL to write this kind of form, and in FEniCS book I read:

Each term of a valid form expression must be a scalar-valued expression integrated exactly once, and they must be linear [...]

which lead me to the idea that I have to calculate $a0 = \int_S w \cdot u\,dS$ first,
and then calculate $\int_S v \cdot w \cdot a0\,dS$.

It should work because $u = \sum_{j=1}^{DOF} N_j \cdot U_j$ and $v = \sum_{i=1}^{DOF} N_i$.

My plan was to do:

a0 = Function(V)
# below I use the fact that basis functions are used as testing functions
# and I can't write assemble(inner(u, w) * ds)
a0.vector()[:] = assemble(inner(v, w) * ds)
P = inner(v, w) * a0 * ds

But it won't work, because I want a0 to be constant over all edges/DOFS (this is also the reason I can't use Discontinous Lagrange function space - because I don't know how to define it over edges of the mesh).

My fallback plan was to work it out on the matrix level:

a0 = assemble(inner(v, w) * ds).array()
p = assemble(inner(v, w) * ds).array() # assume a0 = 1
DOF = len(a0)
a0 = a0.reshape((1, DOF))
p = p.reshape((DOF, 1))
P = np.mat(p) * np.mat(a0)

but it works very slow (because I don't use the fact that the matrices are sparse I guess), and it is not elegant at all.

In general it would be great if I could write something like:

a0 = ... #?
a0.vector()[:] = assemble(inner(v, w) * ds)
assemble(inner(v[i], w) * a0[j] * dS)

which will generate a matrix $K_{i,j}$ where v[i] represents $N_i$ basis function and a0[j] represents a coefficient corresponding to $N_j$ basis function.

Can I solve this problem with UFL/FEniCS?

Answer 1

As you have already noted, the bilinear operator in question is
of rank one with matrix representation $B = b \otimes b$, where

b = assemble(inner(v, w) * ds)

You can form the matrix representation of $B$ explicitly, but a matrix-free implementation may be better. A Krylov subspace solver only requires that action of the operator is implemented, i.e. $$ B x = (b\cdot x)b.$$

An alternative approach may be to introduce an auxilliary variable
$$ \lambda = \int w\cdot u \,dx$$
and solve the constrained problem normally.

See the answers to the related question 7541 for some implementation details.

Comments

Both approaches look great. To go the matrix-free way, I have to implement an operator like CompositeOperator from the linked question, but then I can't use MUMPS or any other non-Krylov subspace solver, right?

I personaly rather use the approach with the auxilliary variable, but should it be a Function, an Argument, VectorConstant? It would be the best to find a type that is evaluated for each DOF, rather than coordinates.

---

The auxiliary variable is handled as an additional unknown, implemented using mixed element method with the additional variable being a scalar:

E = FiniteElement("N1curl", mesh.ufl_cell(), 2)
R = FiniteElement("R", mesh.ufl_cell(), 0)

W = FunctionSpace(mesh, E * R)

u, r = TrialFunctions(W)
v, s = TestFunctions(W)

C = Constant(assemble(1 * ds(domain = mesh))**-1)

a = ... # part of the form that can be handled normally
L = ... 

b = r * inner(w, v) * ds + s * inner(w, u) * ds - C * r * s * ds

solve(a + b == L, ...) 

---

Magne, it will introduce one additional degree of freedom (i.e. E1, E2, E3, ..., En, R1), right?. But for each edge in my mesh I should have different coefficient (R1, R2, R3, ..., Rn).

---

You would need one additional degree of freedom for each term
$$ \int_S w\cdot u\, ds \int_S w\cdot v \, ds, $$
i.e. one for each $S$. If you have one such integral over the boundary, then approaches outlined above are viable. If this is not the case you need a different solution.

Answer 2

Judging from the weak form, it looks like you are trying to do something very similar as I tried here:
http://fenicsproject.org/qa/7541/ufl-unsupported-operand-type-s-for-form-and-form

There are a few possible solutions suggested, including the ones from the answer by Magne Nordaas. But as an easy first proof-of-concept I would recommend the very last one. It is basically assembles the term you are after as a product of an Nx1 and a 1xN matrix.