Specifying inconsisten initial and boundary conditions

I’d like to solve heat problems with inhomogeneous Dirichlet conditions and a non-zero initial condition that need not be consistent with the boundary conditions.

As a simple example, we might have a square whose initial temperature distribution is given by u(x,y,0)=1-x^2. At time zero, the bottom is set to temp 0 and the top is set to temp 1 while the sides are insulated. How can I model this problem?

Following the Parabolic Model tutorial, I see that I can use syntax such as

``````gfu = GridFunction(fes)
gfu.Set(1-x**2)``````

to specify the initial condition and then iteratively redefine gfu via an implicit Euler method that respects any stated Dirichlet conditions. I’m also aware that I can specify specify inhomogendous Dirichlet conditions via code like

``gfu.Set(values_list, definedon=mesh.Boundaries(dirichlet_boundaries))``

where dirichlet_boundaries is a regular expression specifying the boundaries and values_lists contains the values.

I can’t seem to find a way to accomplish both of these however. Perhaps there’s a way to set gfu to values_list on the boundary but 1-x**2 in the interior? Or, perhaps there’s another way to accomplish this?

My code to this point (which is mostly based on the parabolic tutorial) is below:

``````from netgen import gui
from math import pi
from ngsolve import *
from netgen.geom2d import SplineGeometry

time = 0.0
dt = 0.001
tstep = 4

geo = SplineGeometry()
bcs = ("bottom", "right", "top", "left"))
mesh = Mesh( geo.GenerateMesh(maxh=0.1))
fes = H1(mesh, order=3, dirichlet="bottom|top")
gfu = GridFunction(fes)

## This simple block for a single initial condition
#gfu.Set(1-x**2)

##  This block for inhomogenous boundary conditions
boundary_conditions = {'bottom':'0', 'top':'1'}
dirichlet_boundaries = '|'.join(boundary_conditions.keys())
for k in boundary_conditions.keys():
boundary_conditions[k] = eval(boundary_conditions[k])
values_list = [
boundary_conditions[bc] if bc in boundary_conditions.keys() else 0
for bc in mesh.GetBoundaries()
]
gfu.Set(values_list, definedon=mesh.Boundaries(dirichlet_boundaries))

res = gfu.vec.CreateVector()
Draw(gfu,mesh,"init")

u,v = fes.TnT()
a = BilinearForm(fes, symmetric=False)
a.Assemble()
m = BilinearForm(fes, symmetric=False)
m += u*v*dx
m.Assemble()
mstar = m.mat.CreateMatrix()
mstar.AsVector().data = m.mat.AsVector() + dt * a.mat.AsVector()
invmstar = mstar.Inverse(freedofs=fes.FreeDofs())``````

Then,

``````%%time

t_intermediate=0 # time counter within one block-run
while t_intermediate < tstep - 0.5 * dt:
res.data = - dt * a.mat * gfu.vec
gfu.vec.data += invmstar * res
t_intermediate += dt
print("\r",time+t_intermediate,end="")
Redraw(blocking=True)
print("")
time+=t_intermediate``````

Hi,

Set zeros-out the GridFunction where nothing is set, i.e., when only setting on a boundary. The way to get around this ist to use two GridFunctions, Set one to the boundary condition and one to the initial condition. You can then loop over the interior/free degrees of freedoms and copy the initial values into the GridFunction with the correct boundary condition:

``````gfu, gfu_int = GridFunction(fes), GridFunction(fes)
gfu.Set(values_list, definedon=mesh.Boundaries(dirichlet_boundaries))
gfu_int.Set(1-x**2)

free_dofs = fes.FreeDofs()

for dof in range(fes.ndof):
if free_dofs[dof]:
gfu.vec.data[dof] = gfu_int.vec[dof]``````

Best wishes,
Henry

This works perfectly, Henry - Thank you very much!