I can't see my subdomains or meshfunction

mesh
1

Hi, can’t see my subdomains or meshfunction, I do not get any error, but when plotting nothing is shown. When I save my results and use ParaView I get the error:

Cannot read cell connectivity from Cells in piece 0 because the “connectivity” array is not long enough.

here is my code

from fenics import *
from mshr import *
import dolfin as df
import numpy as np
from math import sin, cos, pi, sqrt
import matplotlib.pyplot as plt

# My domain is a cylinder
height_cylinder = 150.00
r = 25.0
my_domain = Cylinder(Point(0, 0, 0), Point(0, 0, height_cylinder), r, r)

# I generate the mesh 
mesh = generate_mesh(my_domain, 48) 

# I have one subdomains which is a cylindrical shell inside the domain

ri_shell = 13.00 
ro_shell = 17.00
tol = 5e-1

class Shell(SubDomain):

    def __init__(self): 
        super().__init__()
    def inside(self, x, on_boundary):
        if ((x[0]**2 + x[1]**2 <= (ro_shell + tol)**2) and (x[0]**2 + x[1]**2 >= (ri_shell - tol)**2 ) and between(x[2], (0, height_cylinder)) ):
            return True

shell = Shell()

# I have a second subdomain which is a small cylinder inside the shell in coordinates x = 3 y y = 2.1
r_small = 3

class SmallCylinder(SubDomain):

    def __init__(self): 
        super().__init__()
        
    def inside(self, x, on_boundary): 
        if ( ( (x[0] - 3.0)**2 + (x[1] - 2.1)**2 <= (r_small + tol)**2 ) and between(x[2], (0, height_cylinder)) ):
            return True

smallCylinder = SmallCylinder() 

# After the subdomains have been created I make a MeshFunction to mark the subdomains
materials = MeshFunction("size_t", mesh, 0, mesh.domains())

shell.mark(materials, 1)
smallCylinder.mark(materials, 2)

# I generate the function spaces
V = FunctionSpace(mesh, 'P', 2)
dx = Measure('dx', domain=mesh, subdomain_data=materials)
V0 = FunctionSpace(mesh, 'DG', 0)

# Trying to plot with dolfin and 
%matplotlib inline
plt.figure(figsize=(20,12))
df.plot(materials, linewidth=0.5, alpha=0.35)
plt.show()

materials_file = df.File("materials.pvd")
materials_file << materials

# With the subdomains marked assign permitivities with a userExpression class:
e_cylinder = 1.27E+1
e_shell = 6.60E+1 
e_small_cylinder = 7.29E+1

class Permittivity(UserExpression): 
   
    def __init__(self, materials, **kwargs):
        super().__init__(**kwargs) 
        self.materials = materials
   
    def eval_cell(self, values, x, ufc_cell):

        if self.materials[ufc_cell.index] == 1:
            values[0] = e_shell
        elif self.materials[ufc_cell.index] == 2:
            values[0] = e_small_cylinder
        else:
            values[0] = e_cylinder

epsilon = Permittivity(mesh, degree=1)
2

I was able to do it in vtkplotter! but the issue that remains is why I can not do it with matplotlib plot or with ParaView

Screenshot_20191003_162056
Screenshot_20191003_162056.png642×643 66.6 KB

3

You can see it in Paraview. . I just did some modifications in your code. You can consider the following:

from fenics import *
from mshr import *
import dolfin as df
import numpy as np
from math import sin, cos, pi, sqrt
import matplotlib.pyplot as plt

height_cylinder = 150
r = 25
my_domain = Cylinder(Point(0, 0, 0), Point(0, 0, height_cylinder), r, r)

mesh = generate_mesh(my_domain, 150)
ri_shell = 13
ro_shell = 17
tol = 5E-1

class Shell(SubDomain):
    def inside(self, x, on_boundary):
        return True if x[0] ** 2 + x[1] ** 2 <= (ro_shell ) ** 2. + DOLFIN_EPS  and x[0] ** 2. + x[1] ** 2. >= (ri_shell ) ** 2.- DOLFIN_EPS and between(x[2], (0, height_cylinder)) else False

shell = Shell()
r_small = 3.

class SmallCylinder(SubDomain):
    def inside(self, x, on_boundary):
        return True if (x[0] - 3.0) ** 2. + (x[1] - 2.1) ** 2. <= ((r_small ) ** 2.)+ DOLFIN_EPS and (between(x[2], (0, height_cylinder))) else False

smallCylinder = SmallCylinder()
materials = MeshFunction('size_t', mesh, mesh.topology().dim())

materials.set_all(0)
shell.mark(materials, 1)
smallCylinder.mark(materials, 2)

materials_file = df.File("materials.pvd")
materials_file << materials

Here is what you can see in Paraview:

PHOTO
PHOTO.png712×473 30.7 KB

1 Like