### Mesh refinement for moving heat source error

75
views
0
10 weeks ago by
I found many tutorials and examples for mesh refinement here. I tried running multiple ways to refine my BoxMesh and they seem to work fine when I plot my mesh. However when I add the mesh refinement piece of code to my bigger heat moving source code, I get some errors. The mesh refinment code I want to use is the following

p=Point(0.01,0.01,0.001)
for i in range(2):

cell_markers = MeshFunction("bool", mesh, mesh.topology().dim())
for c in cells(mesh):
if c.midpoint().distance(p) <= 0.0005:
cell_markers[c] = True
else:
cell_markers[c] = False

mesh = refine(mesh, cell_markers)

Here is part of my code including the refinement piece.

def mark_boundaries_in_hypercube(
mesh, d=3, x0=0, x1=1, y0=0, y1=1, z0=0, z1=1):
"""
Return mesh function FacetFunction with each side in a hypercube
in d dimensions. Sides are marked by indicators 0, 1, 2, ..., 5.
Side 0(x=x0), 1(x=x1), 2(y=y0), 3(y=y1), 4(z=z0), 5(z=z1)
"""
# Define boundary subdomains
tol = 1e-14

class BoundaryX0(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 0.0, tol)

class BoundaryX1(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 0.002, tol)

class BoundaryY0(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 0.0, tol)

class BoundaryY1(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 0.002, tol)

class BoundaryZ0(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[2], 0.0, tol)

class BoundaryZ1(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[2], 0.0010, tol)

# Mark boundaries
boundary_markers = FacetFunction('size_t', mesh)
boundary_markers.set_all(9999)
bx0 = BoundaryX0()
bx1 = BoundaryX1()
by0 = BoundaryY0()
by1 = BoundaryY1()
bz0 = BoundaryZ0()
bz1 = BoundaryZ1()
bx0.mark(boundary_markers, 0)
bx1.mark(boundary_markers, 1)
by0.mark(boundary_markers, 2)
by1.mark(boundary_markers, 3)
bz0.mark(boundary_markers, 4)
bz1.mark(boundary_markers, 5)
return boundary_markers

l = 0.002                # length (m): 100 mm
thickness = 0.001       # thickness (m): 1 mm

divisions = (20, 20, 2)
L = (l, l, thickness)
mesh = BoxMesh(Point(0.0,0.0,0.0), Point(*L), *divisions)

boundary_markers = mark_boundaries_in_hypercube(mesh, 3)

da = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
dv = Measure('dx', domain=mesh, subdomain_data=cells)

Space = FunctionSpace(mesh, 'P', 1)
cells = CellFunction('size_t', mesh)

f = Expression('P/dp*exp(-b*(pow((x[0]-xx),2) + pow((x[1]-yy),2) + pow((x[2]-zz)/dp,2)))',

T = TrialFunction(Space)
del_T = TestFunction(Space)

g0 = 1.0E-14  # Neumann
g1 = 4.0E5  # Neumann
#Preheat = Constant(1273.0) # Dirichlet

boundary_conditions = {0: {'Neumann': g0},
1: {'Neumann': g0},
2: {'Neumann': g0},
3: {'Neumann': g0},
4: {'Neumann': g1},
5: {'Neumann': g0}}

bcs = []
for k in boundary_conditions:
if 'Dirichlet' in boundary_conditions[k]:
bc = DirichletBC(Space, boundary_conditions[k]['Dirichlet'], boundary_markers, k)
bcs.append(bc)

# Neumann boundary condition
integrals_N = []
for k in boundary_conditions:
if 'Neumann' in boundary_conditions[k]:
if boundary_conditions[k]['Neumann'] != 0:
g = boundary_conditions[k]['Neumann']
integrals_N.append(g*del_T*ds(k))

F = (rho*c/Dt*(T-T0)*del_T \
+ kappa*T.dx(i)*del_T.dx(i) \
- rho*f*del_T)*dv \
+ h*(T-Ta)*del_T*da

F += sum(integrals_N)
T_M = Function(Space)
F += emis*st_bol*(T_M**4-Ta**4)*del_T*ds(5)

left=lhs(F)
right=rhs(F)

A = assemble(left)
b = None
T = Function(Space)
​

When I add the refinement code right after I define my mesh I get the following error
line 175, in <module>
- rho*f*del_T)*dv \
TypeError: unsupported operand type(s) for *: 'float' and 'Cell'

When I add it after I define the variables space and cells I get
line 113, in <module>
for c in cells(mesh):
TypeError: 'CellFunctionSizet' object is not callable