### Mesh refinement for moving heat source error

75

views

0

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

Here is part of my code including the refinement piece.

```
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)))',
P=eff*2*Power*scale/(pi*radi**2), b=2/(radi**2), xx=xx, yy=yy, zz=zz, dp = 0.00002)
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

Your help is highly appreciated and thank you in advance

Community: FEniCS Project

Please login to add an answer/comment or follow this question.