Forum Message

 

 

We have moved the forum to https://forum.ngsolve.org . This is an archived version of the topics until 05/05/23. All the topics were moved to the new forum and conversations can be continued there. This forum is just kept as legacy to not invalidate old links. If you want to continue a conversation just look for the topic in the new forum.

Notice

The forum is in read only mode.

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

More
3 years 6 months ago #3180 by Thauwa
I have a system of coupled PDEs, and I use the following to generate test and trial functions:
Code:
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:
Code:
a = BilinearForm(fes, symmetric=True, eliminate_internal=True) a += SymbolicBFI(InnerProduct(E,grad(Et)))

but, the following seems to work:
Code:
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
More
3 years 6 months ago #3182 by mneunteufel
Hi Thauwa,

the code
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

Best
Michael
More
3 years 6 months ago #3183 by Thauwa
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:
[tex]
\nabla\cdot E
[/tex]

According to [2], I must use the following identity to get the variational formulation:
[tex]
\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,
[/tex]

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.


[2]: math.uci.edu/~chenlong/226/FEMMaxwell.pdf
[3]: forum/ngspy-forum/781-fe-space-with-curl-and-div
More
3 years 6 months ago #3184 by mneunteufel
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
Code:
a += SymbolicBFI( InnerProduct(E,grad(Et)))
works.

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
Code:
fes1 = HCurl(mesh, order=1) fes2 = H1(mesh,order=1) fes = FESpace( [fes1,fes2] )

Best
Michael
More
3 years 6 months ago #3194 by Thauwa
Hi Michael,

Thank you very much for your clarifications!

Best,
Thauwa
More
3 years 6 months ago - 3 years 6 months ago #3206 by todd11MR
mneunteufel, thanks for explanation.
Last edit: 3 years 6 months ago by todd11MR.
Time to create page: 0.153 seconds