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)))',
                   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.

Similar posts:
Search »