Defining and plotting subdomains for different materials

1

批注 2020-05-30 223447

I tried to define and plot subdomains for two different materials, but i can’t get correct answer.

from dolfin import *
from fenics import *
import numpy as np
import matplotlib.pyplot as plt
# Create mesh and define function space
nx = ny = 8
mesh = UnitSquareMesh(nx, ny)
V = FunctionSpace(mesh, 'P', 2)
# set subdomains
tol = 1E-10
class Omega_0(SubDomain):
    def inside(self, x , on_boundary):
        return (x[0] >= 0.5 - tol and x[0]<= 0.6+ tol and x[1]>= 0.5-tol and x[1] <= 0.6+tol)
class Omega_1(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] <= 0.5+tol
materials = MeshFunction('size_t', mesh, mesh.topology().dim())
subdomain_0 = Omega_0()
subdomain_1 = Omega_1()
# mark
subdomain_0.mark(materials, 0)
subdomain_1.mark(materials, 1)
# plot subdomain
plot(materials)
plt.savefig('ttheat2/materials.png')
2

This is because the MeshFunction is initialized with the 0 value (as you have not specified any other value in the constructor). Additionally, your mesh is too coarse to resolve this square, consider the following:

from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
# Create mesh and define function space
nx = ny = 20
mesh = UnitSquareMesh(nx, ny)
V = FunctionSpace(mesh, 'P', 2)
# set subdomains
tol = 1E-10
class Omega_0(SubDomain):
    def inside(self, x , on_boundary):
        return (x[0] >= 0.5 - tol and x[0]<= 0.6+ tol and x[1]>= 0.5-tol and x[1] <= 0.6+tol)
class Omega_1(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] <= 0.5+tol
materials = MeshFunction('size_t', mesh, mesh.topology().dim(), 0)
subdomain_0 = Omega_0()
subdomain_1 = Omega_1()
# mark
subdomain_0.mark(materials, 1)
subdomain_1.mark(materials, 2)
# plot subdomain
plot(materials)
plt.savefig('materials.png')
3

Thank you very much, dokken.Your answer is very helpful to me. I modified my code, and the piture is in line with my idea.

from dolfin import *
import numpy as np
import matplotlib.pyplot as plt

Create mesh and define function space

nx = ny = 20
mesh = UnitSquareMesh(nx, ny)

set subdomains

tol = 1E-10
class Omega_0(SubDomain):
def inside(self, x , on_boundary):
return (x[0] >= 0.5 - tol and x[0]<= 0.6+ tol and x[1]>= 0.5-tol and x[1] <= 0.6+tol)
class Omega_1(SubDomain):
def inside(self, x, on_boundary):
return x[0] <= 0.5+tol or x[0] >= 0.6 - tol or x[1] >= 0.6 -tol or x[1] <= 0.5+ tol

mark

materials = MeshFunction(‘size_t’, mesh, mesh.topology().dim(),0)
subdomain_0 = Omega_0()
subdomain_1 = Omega_1()
subdomain_0.mark(materials, 0)
subdomain_1.mark(materials, 1)

plot subdomain

plot(materials)
plt.savefig(‘ttheat2/materials.png’)

4

This code still mark the domain with a zero marker.