Setting Binary initial condition

mesh
1

Hi, I am a novice to FEniCS and having trouble setting an initial condition (3D) to solve a heat equation. What I want to realize is,

CodeCogsEqn

I made an example code below. As I check the vtk file, materials look as I expected (0 on boundary and 1 inside domain), but u_n shows 1 almost everywhere except some part on the boundary. The issue may arise where I define the eval_cell method, but I have no idea how I can fix this.

from fenics import *
import numpy as np

# Create mesh and define function space
mesh = Mesh("meshfile.xml")
V = FunctionSpace(mesh, 'P', 1)

# Define boundary condition
class Boundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary
boundary_D = Boundary()
bc = DirichletBC(V, Constant(0.0), boundary_D)

# Set MeshFunction
materials = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
materials.set_all(1)
boundary_D.mark(materials, 0)

# Set Initial condition
class InitCond(UserExpression):
    def __init__(self, materials, k_0, k_1, **kwargs):
        try:    super(InitCond, self).__init__(**kwargs)
        except: pass
        self.materials = materials
        self.k_0 = k_0
        self.k_1 = k_1
    def eval_cell(self, values, x, cell):
        if self.materials[cell.index] == 0:
            values[0] = self.k_0
        else:
            values[0] = self.k_1
    def value_shape(self):
        return()

u_0 = InitCond(materials, 0.0, 1.0, degree=1)
u_n = interpolate(u_0, V)


# Create VTK file for saving solution
vtkfile = File('u_n.pvd')
vtkfile << u_n
materials_file = File("materials.pvd")
materials_file << materials
2

Consider the following minimal example:

from dolfin import *

mesh = UnitSquareMesh(10,10)

V = FunctionSpace(mesh, "CG", 1)
u_n = Function(V)
u_n.vector()[:] = 1
bc = DirichletBC(V, Constant(0), "on_boundary")
bc.apply(u_n.vector())
File("u.pvd") << u_n
1 Like
3

Thank you for such a quick reply!
It worked with your example.
Thank you!!!