Is this the correct way to use grad() in an H(curl) space?

I have a system of coupled PDEs, and I use the following to generate test and trial functions:

fesm = HCurl(mesh, order=ord, complex=True, dirichlet="outer")
fes = FESpace([fesm,fesm])
# Define test and trial functions:
E, F = fes.TrialFunction()
Et, Ft = fes.TestFunction()   

Now, I want to calculate the following term of a weak form: [tex]\int E \nabla (Et) [/tex].

When I naively try the following for my bilinear form, I get a dimensional mismatch error:

a = BilinearForm(fes, symmetric=True, eliminate_internal=True)
a += SymbolicBFI(InnerProduct(E,grad(Et)))

but, the following seems to work:

a += SymbolicBFI(InnerProduct(E,(grad(Et)[0],grad(Et)[1],grad(Et)[2])))

Is this the correct thing to do? I am still learning NGSolve. Thank you all for your hard work.

I saw a different post [1] saying that VectorH1 spaces will allow all div and curl and grad to be used simultaneously, but I wanted to make sure my attempt makes sense too.

[1] /forum/ngspy-forum/781-fe-space-with-curl-and-div

Hi Thauwa,

the code

a += SymbolicBFI(InnerProduct(E,grad(Et)))

leads to an Exception as E is a vector (with dimension three as I guess that you have a 3D mesh) whereas grad(Et) is a 3x3 matrix and InnerProduct expects that both inputs have the same dimension. Thus, your second code works, but then you have a different integral.

All integrators need at the end a scalar valued function as input, but your integral [tex]\int E\nabla Et[/tex] seems to be vector valued?

Note that if you want the Jacobian you need the “Grad” operator. With “grad” you obtain the transposed Jacobian: grad() = Grad().trans


Hi Michael,

Thank you kindly for your response.
So, my weak formulation appears to be incorrect. In my system of coupled PDEs that focuses on curl of associated vectors (derived from a standard Maxwell system), I have the following term:
\nabla\cdot E

According to 2, I must use the following identity to get the variational formulation:
\int_{\Omega} \nabla\cdot E E_t dx =
-\int_{\Omega} E\cdot \nabla E_t dx +
\int_{\partial\Omega} n \cdot E E_t dS,

where n is the unit outwards normal vector. So, it seems to me that I cannot avoid a dimensional mismatch in the inner product because of the gradient operation. I thought I could just have a derivative taken on each row of my vector and summed up, as usual. That is, for elements v_i of E_t:

[tex](v_1, v_2, v_3) = (\partial_x v_1 + \partial_y v_1 + \partial_z v_1, \partial_x v_2 + \partial_y v_2 + \partial_z v_2, \partial_x v_3 + \partial_y v_3 + \partial_z v_3).[/tex]

But I am not sure whether this is valid anymore.
Do you have any advice on
how I can implement this [tex]E \nabla E_t[/tex] properly in NGSolve?
In the NGSolve forum post 3, they suggest using VectorH1 spaces for situations where grad(), curl() and div() might have to be used simultaneously, but I am concerned it will conflict with HCurl()'s reputation of being suitable for Maxwell-type problems involving curl terms.

Hi Thauwa,

for this integration by parts formula to hold you have (normally) that E is vector valued and Et is scalar valued. Then the term [tex]\int E\cdot\nabla Et[/tex] in NGSolve

a += SymbolicBFI( InnerProduct(E,grad(Et)))


Without knowing the full (coupled) PDE in strong form (+ dimension of the unknowns) I can only guess that one of the unknowns is vector valued (in HCurl) and the other scalar valued (in H1), something like

fes1 = HCurl(mesh, order=1)
fes2 = H1(mesh,order=1)
fes = FESpace( [fes1,fes2] )


Hi Michael,

Thank you very much for your clarifications!


mneunteufel, thanks for explanation.