# Setting Binary initial condition

**URL:** <https://fenicsproject.discourse.group/t/setting-binary-initial-condition/3443>\
**Category:** mesh\
**Created:** [June 5, 2020, 5:30am UTC](https://fenicsproject.discourse.group/t/setting-binary-initial-condition/3443 "2020-06-05T05:30:24Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![Ash](https://avatars.discourse-cdn.com/v4/letter/a/f19dbf/32.png) [@Ash](https://fenicsproject.discourse.group/u/Ash)\
**Post date:** [June 5, 2020, 5:30am UTC](https://fenicsproject.discourse.group/t/setting-binary-initial-condition/3443/1 "2020-06-05T05:30:24Z")

</div>

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](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/1/1da5b2b74d1a4e057d0214a9bf5db64708e6e862.gif)

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
```

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [June 5, 2020, 6:26am UTC](https://fenicsproject.discourse.group/t/setting-binary-initial-condition/3443/2 "2020-06-05T06:26:22Z")

</div>

Consider the following minimal example:

```auto
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

```

---

<div class="post-metadata">

**Author:** ![Ash](https://avatars.discourse-cdn.com/v4/letter/a/f19dbf/32.png) [@Ash](https://fenicsproject.discourse.group/u/Ash)\
**Post date:** [June 5, 2020, 6:53am UTC](https://fenicsproject.discourse.group/t/setting-binary-initial-condition/3443/3 "2020-06-05T06:53:33Z")

</div>

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