How to implement a piece wise constant? (curl-curl,electromagnetic Maxwell equations)


215
views
1
5 months ago by
Hi,
If a constant appearing in my equation has different values in different portions of the domain, how do I implement it; I mean, whats the official way?
Scenario: I have a cubic domain 2m^3, ranging from -1 to 1 each of the three dimensions. In the positive z-region, this constant si (sigma) takes the value 1E-9, in the negative z-region (inclusive of 0), the value is 7.3 unless the sigma is inside a 14cm^3 sub-domain. The range of this subdomain appears in the code. Within this small cube in the lower part of the domain, the value is 63000.

(1) If I define the constant as an "expression" which takes different values depending on the mesh coordinates, will this be alright? This is illustrated further in a minimal working code as follows. This code save the system matrix (in Matlab format) as C.dat. I have been following this method for quite sometime, but recently I was pointed to another method, which appears in No. 2 beneath this code. I want to know for certain which of the two methods is the correct FEniCS way for implementing a piece-wise constant.
from dolfin import *
from scipy.sparse import csr_matrix
import scipy.io as io


def main():
    mesh = BoxMesh(Point(-1,-1,-1),Point(1,1,1),32, 32, 32) # unit cube mesh    

    V = FunctionSpace(mesh, "Nedelec 1st kind H(curl)", 1)
    L = FunctionSpace(mesh, "Lagrange", 1)
    E = TrialFunction(V)
    F = TestFunction(V)  

    xl = -0.07
    xr = 0.07
    yl = -0.07
    yr = 0.07
    zr = -0.02
    zl = -0.16    

    # Set the piecewise constant sigma
    si1 = 1.0E-9
    si2 = 7.3
    si3 = 63000
    In_Cube_x = "x[0] >= (xl-DOLFIN_EPS) && x[0] <= (xr+DOLFIN_EPS)"
    In_Cube_y = "x[1] >= (yl-DOLFIN_EPS) && x[1] <= (yr+DOLFIN_EPS)"
    In_Cube_z = "x[2] >= (zl-DOLFIN_EPS) && x[2] <= (zr+DOLFIN_EPS)"
    In_Cube   = In_Cube_x + "&&" + In_Cube_y + "&&" + In_Cube_z
    Val_sigma  = In_Cube + "? si3 : (x[2] <= 0 ? si2 : si1)"
    si   = Expression((Val_sigma),degree=0, si1=si1, si2=si2, si3=si3, xl=xl, xr=xr, yl=yl, yr=yr, zl=zl, zr=zr, domain=mesh)

    curlcurl = inner(curl(E), curl(F))*dx
    mass_si = inner(si*E, F)*dx  
    a = curlcurl + mass_si
    rhs =  inner(Constant((1.0,1.0,1.0)), F)*dx
    
    bc  = DirichletBC(V,Constant((0.0,0.0,0.0)),DirichletBoundary())
    A, b = assemble_system(a, rhs, bc)

    A = as_backend_type(A).mat()
    r, c = A.getSize()
    rptr, cind, vals = A.getValuesCSR()
    A = csr_matrix((vals,cind,rptr), shape=(r,c))
    io.savemat('C.mat',{'C':A})			

# Dirichlet boundary
class DirichletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary
###############################################
main()​

(2) This is the second method. The aim is the same as before, i.e., to let sigma (si) take different values in different parts of the domain, and use it in the forms. This method saves its matrix (in Matlab format) as A.dat.

from dolfin import *
from scipy.sparse import csr_matrix
import scipy.io as io

def main():  
    
    mesh = BoxMesh(Point(-1,-1,-1),Point(1,1,1),32, 32, 32) # unit cube mesh    

    V = FunctionSpace(mesh, "Nedelec 1st kind H(curl)", 1)
    L = FunctionSpace(mesh, "Lagrange", 1)
    E = TrialFunction(V)
    F = TestFunction(V)    
    # subdomains
    Sub1 = sub1()
    Sub2 = sub2()
    Sub3 = sub3()
    #Initialize mesh function for these domains
    regions = CellFunction("size_t", mesh)
    regions.set_all(0)
    Sub1.mark(regions,1)
    Sub2.mark(regions,2)
    Sub3.mark(regions,3)
    dx = Measure('dx', domain = mesh, subdomain_data = regions)

    # Define sigma in the three regions
    si1 = 1E-9
    si2 = 7.3
    si3 = 63000

    curlcurl = inner(curl(E), curl(F))*dx(1) + inner(curl(E), curl(F))*dx(2) + inner(curl(E), curl(F))*dx(3)    
    mass_si = inner(si1*E, F)*dx(1) + inner(si2*E, F)*dx(2) + inner(si3*E, F)*dx(3)  
    a = curlcurl + mass_si
    rhs = inner(Constant((1.0,1.0,1.0)), F)*dx(1) + inner(Constant((1.0,1.0,1.0)), F)*dx(2) + inner(Constant((1.0,1.0,1.0)), F)*dx(3)
    
    # curlcurl = inner(curl(E), curl(F))*dx(0)
    # mass_si = inner(si1*E, F)*dx(1) + inner(si2*E, F)*dx(2) + inner(si3*E, F)*dx(3)  
    # a = curlcurl + mass_si
    # rhs = inner(Constant((1.0,1.0,1.0)), F)*dx(0)

    bc  = DirichletBC(V,Constant((0.0,0.0,0.0)),DirichletBoundary())        
    A, b = assemble_system(a,rhs,bc)


    A = as_backend_type(A).mat()
    r, c = A.getSize()
    rptr, cind, vals = A.getValuesCSR()
    A = csr_matrix((vals,cind,rptr), shape=(r,c))
    io.savemat('A.mat',{'A':A})

################################### Air, brine, graphite cube
class sub1(SubDomain):
    def inside(self,x, on_boundary):
        return x[2]>0

class sub2(SubDomain):
    def inside(self, x, on_boundary):
        return ((x[2]<=0) and not(between(x[0],(-0.07,0.07)) and between(x[1],(-0.07,0.07)) and between(x[2],(-0.16,-0.02))))

class sub3(SubDomain):
    def inside(self,x, on_boundary):
        return (between(x[0],(-0.07,0.07)) and between(x[1],(-0.07,0.07)) and between(x[2],(-0.16,-0.02)))
##################################################################### CLASS OVER-RIDE RETURNING DIRICHLET BC
# Dirichlet boundary
class DirichletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary

main()


When I load the two matrices in Matlab and compare their relative difference in Frobenius norms, it is significantly large. Shouldn't the two methods be equivalent?

Thanks,
Hisham

Community: FEniCS Project

1 Answer


1
5 months ago by
I don't think there is an 'official' way. But I've had good experiences with this approach:
http://home.simula.no/~hpl/homepage/fenics-tutorial/release-1.0-nonabla/webm/materials.html
The basic idea is to make your material parameter a function in a discontinuous piecewise constant function space, i.e. one constant value per cell, and set the values manually.
This is very close to your first version but there may be differences if the criteria in your expressions for identifying the subdomains don't match the cell boundaries.

Your second version should be equivalent in the results. But it's not as tidy and not as close to the mathematical formulation of the problem as the first variant. Imagine you want to include another region with a different parameter. Then you'd have to include another integral, whereas your weak form stays the same in the first approach.

It's just a matter of taste I suppose.
Hi klunkean,
Many thanks. I was unaware of the DG method. i will try this out, and see if this is equivalent to one of the other two. It certainly is more elegant than both of them, and one doesn't require adding an integral for inclusion of another subdomain.

Cheers,
Hisham
written 5 months ago by HbZ Syed  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »