绘制左心室网格

# triangulation
from fenics import *
from mshr import *
big = Ellipsoid(Point(0,0,0),10,10,20)
small = Ellipsoid(Point(0,0,0),7,7,17)
gaizi = Box(Point(-11,-11,5),Point(11,11,21))
domain = big-small-gaizi
mesh = generate_mesh(domain,50)
# File("mesh.xml") << mesh

# mesh = Mesh("mesh.xml")
# mark the boundary
mesh_function = MeshFunction("size_t", mesh, 2, value = 0)

class VentricleBase(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[2],5.0)


class VentricleWall(SubDomain):
    def inside(self, x, on_boundary):
        temp = (x[0]*x[0] + x[1]*x[1])/7.0/7.0 + x[2]*x[2]/17.0/17.0
        return on_boundary and temp < 1.1


venbase = VentricleBase()
venbase.mark(mesh_function, 1)
venwall = VentricleWall()
venwall.mark(mesh_function, 2)

# file = File("base.xml")
# file << mesh_function

# visulize mesh_function
V = FunctionSpace(mesh, "P", 1)

bcs = [DirichletBC(V, Constant(100), mesh_function, 1),
      DirichletBC(V, Constant(10000), mesh_function, 2)]

u = Function(V)
for bc in bcs:
    bc.apply(u.vector())