DOLFIN

Function conditional on subdomain

Asked by Christopher Laing

Hi there,

I'm trying to create a scalar function which is zero on the boundary of my mesh and 1 everywhere else.

I have seen this question (https://answers.launchpad.net/dolfin/+question/213769), and attempted to emulate it:

mesh = UnitCircle(r)
S = FunctionSpace(mesh, "CG", 2)

boundary_parts = MeshFunction("uint", mesh, mesh.topology().dim()-1)
boundary_parts.set_all(5)

class Boundary(SubDomain):
 def inside(self, x, on_boundary):
  tol = 1
  return on_boundary and abs(x[0]**2+x[1]**2-1) < tol

Gamma_0 = Boundary()
Gamma_0.mark(boundary_parts,0)

class const_h(Expression):
 def eval_cell(self, values, x, cell):
  k = boundary_parts.array()[cell.index]
  if k==5:
   values[0] = 1.0
  else:
   values[0] = 0.0

hconst = const_h()
h = interpolate(hconst,S)

This doesn't appear to be producing my desired function, most likely because I don't really understand what I'm doing.

Any help would be appreciated.

Chris

Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Johan Hake
Solved:
Last query:
Last reply:
Revision history for this message
Best Johan Hake (johan-hake) said : 		#1

There are many ways to do it. First read up on how you use SubDomain in
the FEniCS book. It looks like you got it slightly wrong.

from dolfin import *
mesh = UnitCircle(6)
S = FunctionSpace(mesh, "CG", 2)
h = Function(S)
h.vector()[:] = 1
boundary_values = DirichletBC(S, 0., lambda x, on_boundary:on_boundary)
boundary_values.apply(h.vector())
plot(h, interactive=True)

Johan

On 11/12/2012 03:01 AM, Christopher Laing wrote:
> New question #213942 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/213942
>
> Hi there,
>
> I'm trying to create a scalar function which is zero on the boundary of my mesh and 1 everywhere else.
>
> I have seen this question (https://answers.launchpad.net/dolfin/+question/213769), and attempted to emulate it:
>
> mesh = UnitCircle(r)
> S = FunctionSpace(mesh, "CG", 2)
>
> boundary_parts = MeshFunction("uint", mesh, mesh.topology().dim()-1)
> boundary_parts.set_all(5)
>
> class Boundary(SubDomain):
> def inside(self, x, on_boundary):
> tol = 1
> return on_boundary and abs(x[0]**2+x[1]**2-1) < tol
>
> Gamma_0 = Boundary()
> Gamma_0.mark(boundary_parts,0)
>
> class const_h(Expression):
> def eval_cell(self, values, x, cell):
> k = boundary_parts.array()[cell.index]
> if k==5:
> values[0] = 1.0
> else:
> values[0] = 0.0
>
> hconst = const_h()
> h = interpolate(hconst,S)
>
> This doesn't appear to be producing my desired function, most likely because I don't really understand what I'm doing.
>
> Any help would be appreciated.
>
> Chris
>

Revision history for this message
Christopher Laing (9e9o1k-chris) said : 		#2

Johan,

Thanks heaps for that, I'll try it and mark the question as resolved if it works.

I'll also do my utmost to figure out what it's doing, and will read up on subdomains. Sorry if my question was a bit amateur. Before starting my present project, I'd never written a line of code or done a numerical simulation, so my learning curve has been fairly steep, and I'm nowhere near the top of it yet!

Thanks again,

Chris

Revision history for this message
Christopher Laing (9e9o1k-chris) said : 		#3

Thanks Johan Hake, that solved my question.