How to apply boundary conditions in a mixed-space setting (Nedelec & nodal Galerkin).

Question

Hi FEniCS users,

I have a system of two PDEs which I want to discretize, and then port the relevant matrices to another solution platform. I have so far used assemble_system with success. These are Maxwell's equations in an electromagnetic geophysics application. In the weak form, I have an integral over the inner product of a test function belonging to the Nedelec 1st kind H(curl) space (TA), and the gradient of a trial function belonging to a single degree nodal Galerkin space. I don't quite know how to implement strong Dirichlet BC (homogeneous in this case) to the matrix that represents this particular form.

I am working on the following code; which has 4 assemble_system statements, and I don't know how to apply the boundary conditions in the 3rd usage (which is currently commented out).

    parameters["linear_algebra_backend"] = "PETSc";
    mesh     = UnitCubeMesh(5,5,5)    
    mesh     = refine(mesh)
    VN        = FunctionSpace(mesh, "Nedelec 1st kind H(curl)", 1)
    VG       = FunctionSpace(mesh, "Lagrange", 1)   
    A = TrialFunction(VN)
    phi = TrialFunction(VG)
    TA = TestFunction(VN)
    Tphi = TestFunction(VG)  

    om = 1e-4    
    mu = 1         
    si = 1           
    ep = 1          

    zero   = Constant((0.0,0.0,0.0))   
    BCN     = DirichletBC(VN,zero,DirichletBoundary())
    BCG      = DirichletBC(VG, Constant((0.0)),DirichletBoundary())  (Nedelec)
    c = inner(curl(TA), curl(A))*dx     
    m1 = inner(TA,A)*dx                  
    m2 = inner(TA,grad(phi))*dx       
    #m3 = inner(grad(Tphi),A)*dx     
    m4 = inner(grad(Tphi),grad(phi))*dx  

    dbdt = Constant((0.0, 0.0, 1.0))   
    L1 = inner(TA, dbdt)*dx               
    L2 = inner(grad(Tphi),dbdt)*dx     

    C = PETScMatrix()                      
    M1 = PETScMatrix()
    M2= PETScMatrix()
    M4 = PETScMatrix()

    bc = PETScVector()
    bm1 = PETScVector()
    bm2= PETScVector()
    bm4 = PETScVector()

    assemble_system(c,L1,BCN, x0=None, form_compiler_parameters=None, add_values=False, finalize_tensor=True, keep_diagonal=False, A_tensor=C, b_tensor=bc, backend=PETSc)
    assemble_system(m1,L1,BCN, x0=None, form_compiler_parameters=None, add_values=False, finalize_tensor=True, keep_diagonal=False, A_tensor=M1, b_tensor=bm1, backend=PETSc)
    #assemble_system(m2,L1,???, x0=None, form_compiler_parameters=None, add_values=False, finalize_tensor=True, keep_diagonal=False, A_tensor=M2, b_tensor=bm2, backend=PETSc)    
    assemble_system(m4,L2,BCG, x0=None, form_compiler_parameters=None, add_values=False, finalize_tensor=True, keep_diagonal=False, A_tensor=M4, b_tensor=bm4, backend=PETSc)    

Any help is deeply appreciated; criticism of my usage, and alternate ways of handling this kind of mixed problem are also warmly welcome with thanks.

Cheers,
Hisham

Comments

The code would be easier to read if one didn't have to scroll left and right all the time. I'd scratch def main, hardcode Nx and friends to something, and put the comments _above_ the line they comment on, not right behind.

Apart from that, it's a little unclear to me what it is you're trying to achieve. Can it not be broken down to something more elemental? Like, like, a 2x2 block matrix of Laplaces?

---

Hi Nico,

Done. Comments have been removed, they were not useful.
Nx, Ny, Nz, are now just 5,5,5.

I have a system of 2 PDEs which after discretization I want to develop a block preconditioning technique for.

Two variables appear in the weak form, A (which belongs to the Nedelec first order element space), and phi (which belongs to nodal Lagrangian space of degree 1). I want to assemble the bilinear form m2 into a PETSc matrix. The only hurdle seems to be the application of strong Dirichlet boundary conditions throughout the unit cube domain. The way it appears implemented for c, m1, and m4 matrices works fine. In these three bilinear forms, the space to which the two components of the inner product belong is the same. Nedelec in c, and m1, and Lagrange in m4. However, in m2, one component belongs to the Nedelec space, and the other to Lagrangian. So I just want to figure out, what to put in place where there are currently question marks in the second last line of the code.

Scott suggested to use [BCN, BCG], and although the code then runs without errors, I just need to confirm if this is indeed the right way to bring this about.

Does, this explanation help any further?
Cheers,
Hisham

---

I forgot to comment the original question, despite I would still suggest to use implicit boundary conditions:

"Scott suggested to use [BCN, BCG], and although the code then runs without errors, I just need to confirm if this is indeed the right way to bring this about. "

YES, this is the correct way.

Cheers, Raphael

---

Thanks Raphael.

Answer 1

hi, can I ask you a question?

When you set the boundary condition, what's the meaning of
zero = Constant((0.0,0.0,0.0)) # Expression for zero Dirichlet BC on the Nedelec space
BCN = DirichletBC(VN,zero,DirichletBoundary()) # Zero Dirichlet Boundary conditions (Nedelec) ?

Your boundary condition is A=0 or A*n=0 on boundary?

Comments

Hi,
My boundary condition is A=0 on the boundary. A represents the vector potential (part) of the electric field, and the boundaries are far from the actual domain of interest; so the field dissipates and is negligible.

Cheers,
Hisham

Answer 2

Hi Hisham,

if you want to solve geophysical scattering problems which I assume, just don't set explicit boundary conditions and everthing will be fine! Setting explicit Zero-Dirichlet-BC, disregarding if I use the field or potential approach, made my solutions only more inaccurate at the domain boundaries. Nevertheless, I'm actually not sure if the implementation, as already mentioned in the following question, is correct on Nedelec spaces:

https://fenicsproject.org/qa/11942/how-to-implement-the-boundary-condition-of-maxwell-equation