specifying an equation to particular subdomain(s)

Asked by Hossein Nili

Hi

Is there a way I can solve an equation in, say, two of three subdomains and set the solution variable identical to zero everywhere in the third subdomain?

[I am trying to avoid importing of mesh from outside FEniCS, hence the inclusion of a third (in fact redundant!) subdomain]

Many thanks,

Question information

Language:
English Edit question
Status:
Solved
For:
FEniCS Project Edit question
Assignee:
No assignee Edit question
Solved by:
Jan Blechta
Solved:
Last query:
Last reply:
Revision history for this message
Jan Blechta (blechta) said :
#1

On Thu, 02 May 2013 16:16:18 -0000
Hossein Nili <email address hidden> wrote:
> New question #228027 on FEniCS Project:
> https://answers.launchpad.net/fenics/+question/228027
>
> Hi
>
> Is there a way I can solve an equation in, say, two of three
> subdomains and set the solution variable identical to zero everywhere
> in the third subdomain?

I believe that this infrastructure is currently being implemented. So
far you can use a hack. You use DirichletBC setting solution on
(interior) domain where you don't want to define your problem. I used it
succesfully when solving multiphysics problem with flow only in a part
of domain. I did something like

bc_u_zero = DirichletBC(W.sub(0), (0.0, 0.0), zero_velocity_edges, 1)
bc_p_zero = DirichletBC(W.sub(1), 0.0, zero_pressure_points, "pointwise")

Difference is that velocity condition works on (some) facets and pressure
condition on vertices (assuming W.sub(1) is CG1 space to make sense).
Defining zero_pressure_points to make it work in parallel is quite tricky.

Jan

>
> [I am trying to avoid importing of mesh from outside FEniCS, hence
> the inclusion of a third (in fact redundant!) subdomain]
>
> Many thanks,
>

Revision history for this message
Hossein Nili (hossein-nili) said :
#2

Thanks for responding.

I'm afraid I can't follow the trick - is there a fuller version of this (or a similar) code, so that I can better see what is happening? I cannot understand the boundary condition definitions you have made and how they address this issue.

Many thanks,

Hossein

Revision history for this message
Jan Blechta (blechta) said :
#3

On Fri, 03 May 2013 10:56:13 -0000
Hossein Nili <email address hidden> wrote:
> Question #228027 on FEniCS Project changed:
> https://answers.launchpad.net/fenics/+question/228027
>
> Status: Answered => Open
>
> Hossein Nili is still having a problem:
> Thanks for responding.
>
> I'm afraid I can't follow the trick - is there a fuller version of
> this (or a similar) code, so that I can better see what is happening?
> I cannot understand the boundary condition definitions you have made
> and how they address this issue.
>
> Many thanks,
>
> Hossein
>

Well, this was only very problem-specific example, so it is no worth to
explain here how zero_velocity_edges and zero_pressure_points are
computed.

Generally, you need to define Dirichlet subdomain according to your
needs. There's plenty of such methods - check
http://fenicsproject.org/documentation/dolfin/1.2.0/python/programmers-reference/cpp/fem/DirichletBC.html#dolfin.cpp.fem.DirichletBC

Can you post your mesh and tell us where do you want to define
Dirichlet subdomain?

Jan

Revision history for this message
Hossein Nili (hossein-nili) said :
#4

Thanks for replying.

Here's the actual problem: I am solving PDEs that are not terribly complicated (time-consuming) to solve, yet the challenge lies in a boundary that is moving with time. And the moving boundary does not have a simple shape; a simplified assumption for it is to be sinusoidal shaped.

When I attempt to define subdomains with inner boundaries (not moving yet) that are sinusoidal in shape, I get an error message saying the different parameter values I want to attribute to the different subdomains cannot be attributed as requested. I use the recommended syntax of:

m_values = [1, 0, 0] # values of m in the three subdomains

for cell_no in range(len(subdomains.array())):
    subdomain_no = subdomains.array()[cell_no]
    m.vector()[cell_no] = m_values[subdomain_no]

I understand why: the sinusoidal shape leaves cells that belong to none of the subdomains and that crashes the solver.

My question therefore is whether or not I can define subdomains (with different parameters) that are separated by irregularly-shaped boundaries (e.g. sinusoidal, Gaussian, etc.) in FEniCS? Do I have to define the mesh elsewhere, e.g. in Gmsh, and if I do, is there a two-way communication protocol between FEniCS and, say, Gmsh, so that I can update the subdomain boundaries with new shapes as the boundary moves after a time step?

Apologies for the lengthiness and many many thanks for helping,

Hossein

Revision history for this message
Jan Blechta (blechta) said :
#5

On Wed, 08 May 2013 11:01:30 -0000
Hossein Nili <email address hidden> wrote:
> Question #228027 on FEniCS Project changed:
> https://answers.launchpad.net/fenics/+question/228027
>
> Status: Answered => Open
>
> Hossein Nili is still having a problem:
> Thanks for replying.
>
> Here's the actual problem: I am solving PDEs that are not terribly
> complicated (time-consuming) to solve, yet the challenge lies in a
> boundary that is moving with time. And the moving boundary does not
> have a simple shape; a simplified assumption for it is to be
> sinusoidal shaped.
>
> When I attempt to define subdomains with inner boundaries (not moving
> yet) that are sinusoidal in shape, I get an error message saying the
> different parameter values I want to attribute to the different
> subdomains cannot be attributed as requested. I use the recommended
> syntax of:
>
> m_values = [1, 0, 0] # values of m in the three subdomains
>
> for cell_no in range(len(subdomains.array())):
> subdomain_no = subdomains.array()[cell_no]
> m.vector()[cell_no] = m_values[subdomain_no]
>
> I understand why: the sinusoidal shape leaves cells that belong to
> none of the subdomains and that crashes the solver.
>
> My question therefore is whether or not I can define subdomains (with
> different parameters) that are separated by irregularly-shaped
> boundaries (e.g. sinusoidal, Gaussian, etc.) in FEniCS? Do I have to
> define the mesh elsewhere, e.g. in Gmsh, and if I do, is there a
> two-way communication protocol between FEniCS and, say, Gmsh, so that
> I can update the subdomain boundaries with new shapes as the boundary
> moves after a time step?
>
> Apologies for the lengthiness and many many thanks for helping,
>
> Hossein
>

If I understand your problem correctly you want to define regions not
fitting mesh cells. But I don't fully understand the way you're trying
to implement this as you don't tell us which DOLFIN objects you use and
how are they initialized. Note that there are plenty of ways. It could
be helpful if you gave very short complete code.

If you're using SubDomain.mark() method then bug 1119803 may be
relevant.

Regarding mesh generation DOLFIN can't instruct external mesh
generators except CGAL CSG. For example I used to instruct Triangle by
calling it as shell command from python then calling DOLFIN mesh
conversion script. It was a little bit tricky that input file with
geometry for Triangle contained unspecified values which were
substituted on each call, so I was able to generate mesh by DOLFIN
program itself for run-time computed geometry parameter. So you need to
program it by yourself. But I guess your intent is to have static mesh
and update only definitions of some regions - you don't need to
re-execute mesh generator for this.

Jan

Revision history for this message
Jan Blechta (blechta) said :
#6

On Wed, 08 May 2013 11:01:30 -0000
Hossein Nili <email address hidden> wrote:
> Question #228027 on FEniCS Project changed:
> https://answers.launchpad.net/fenics/+question/228027
>
> Status: Answered => Open
>
> Hossein Nili is still having a problem:
> Thanks for replying.
>
> Here's the actual problem: I am solving PDEs that are not terribly
> complicated (time-consuming) to solve, yet the challenge lies in a
> boundary that is moving with time. And the moving boundary does not
> have a simple shape; a simplified assumption for it is to be
> sinusoidal shaped.
>
> When I attempt to define subdomains with inner boundaries (not moving
> yet) that are sinusoidal in shape, I get an error message saying the
> different parameter values I want to attribute to the different
> subdomains cannot be attributed as requested. I use the recommended
> syntax of:
>
> m_values = [1, 0, 0] # values of m in the three subdomains
>
> for cell_no in range(len(subdomains.array())):
> subdomain_no = subdomains.array()[cell_no]
> m.vector()[cell_no] = m_values[subdomain_no]
>
> I understand why: the sinusoidal shape leaves cells that belong to
> none of the subdomains and that crashes the solver.
>
> My question therefore is whether or not I can define subdomains (with
> different parameters) that are separated by irregularly-shaped
> boundaries (e.g. sinusoidal, Gaussian, etc.) in FEniCS? Do I have to
> define the mesh elsewhere, e.g. in Gmsh, and if I do, is there a
> two-way communication protocol between FEniCS and, say, Gmsh, so that
> I can update the subdomain boundaries with new shapes as the boundary
> moves after a time step?
>
> Apologies for the lengthiness and many many thanks for helping,
>
> Hossein
>

https://answers.launchpad.net/dolfin/+question/221327 may also be
relevant

Revision history for this message
Hossein Nili (hossein-nili) said :
#7

Thanks.

That's right, I wish to define regions not fitting mesh cells. At the bottom of this message is the bit of the code where I define the subdomains. I understand there's a way round this, right?

On the moving boundary, that is right, I want a static mesh but need to be able to update the boundary with a self-defined function (expression in terms of spatial coordinates). So if I find a way to define subdomain boundaries that follow self-defined mathematical expressions (e.g. a sine function like below), I have got what I need.

Here is the code part. I understand why I get error running this, yet I have included it to show what I am after.

Thank you,

from dolfin import *
import math, numpy, sys
# choose solver type (automated approach for solving nonlinear PDE)
answer = sys.argv[1] # g (GMRES) or l (sparse LU) solver
iterative_solver = True if answer == 'g' else False

# define mesh
h=0.1; hb=0.1
mesh = Rectangle(0, 0, 1, htopv, 100, 100)

# define a meshfunction for numbering subdomains
subdomains = MeshFunction("uint", mesh, 2)

tol = 1E-14
# define the subdomains
class Biomass(SubDomain):
 def inside(self, x, on_boundary):
         return True if x[1] <= 0.1 + 0.05*sin(pi*x[0]) else False
class Interface(SubDomain):
 def inside(self, x, on_boundary):
  return True if x[1] >= 0.1 + 0.05*sin(pi*x[0]) and x[1] <= 0.2 + 0.05*sin(pi*x[0]) else False
class Bulk(SubDomain):
 def inside(self, x, on_boundary):
  return True if x[1] >= 0.2 + 0.05*sin(pi*x[0]) else False

# mark the subdomains with numbers
subdomain0 = Biomass()
subdomain0.mark(subdomains, 0)
subdomain1 = Interface()
subdomain1.mark(subdomains, 1)
subdomain2 = Bulk()
subdomain2.mark(subdomains, 2)

# define functionspace for m
V0 = FunctionSpace(mesh, "DG", 0)
m = Function(V0)

# attribute m values to appropriate subdomains
m_values = [1, 0, 0] # values of m in the two subdomains
for cell_no in range(len(subdomains.array())):
    subdomain_no = subdomains.array()[cell_no]
    m.vector()[cell_no] = m_values[subdomain_no]

Revision history for this message
Jan Blechta (blechta) said :
#8

On Wed, 08 May 2013 14:41:43 -0000
Hossein Nili <email address hidden> wrote:
> Question #228027 on FEniCS Project changed:
> https://answers.launchpad.net/fenics/+question/228027
>
> Status: Answered => Open
>
> Hossein Nili is still having a problem:
> Thanks.
>
> That's right, I wish to define regions not fitting mesh cells. At the
> bottom of this message is the bit of the code where I define the
> subdomains. I understand there's a way round this, right?

Maybe:)

>
> On the moving boundary, that is right, I want a static mesh but need
> to be able to update the boundary with a self-defined function
> (expression in terms of spatial coordinates). So if I find a way to
> define subdomain boundaries that follow self-defined mathematical
> expressions (e.g. a sine function like below), I have got what I need.
>
> Here is the code part. I understand why I get error running this, yet
> I have included it to show what I am after.

I can't reproduce an error - could you report which version of FEniCS do
you use and post an error message?

>
> Thank you,
>
> from dolfin import *
> import math, numpy, sys
> # choose solver type (automated approach for solving nonlinear PDE)
> answer = sys.argv[1] # g (GMRES) or l (sparse LU) solver
> iterative_solver = True if answer == 'g' else False
>
> # define mesh
> h=0.1; hb=0.1
> mesh = Rectangle(0, 0, 1, htopv, 100, 100)
>
> # define a meshfunction for numbering subdomains
> subdomains = MeshFunction("uint", mesh, 2)
>
> tol = 1E-14
> # define the subdomains
> class Biomass(SubDomain):
> def inside(self, x, on_boundary):
> return True if x[1] <= 0.1 + 0.05*sin(pi*x[0]) else
> False class Interface(SubDomain):
> def inside(self, x, on_boundary):
> return True if x[1] >= 0.1 + 0.05*sin(pi*x[0]) and
> x[1] <= 0.2 + 0.05*sin(pi*x[0]) else False class Bulk(SubDomain):
> def inside(self, x, on_boundary):
> return True if x[1] >= 0.2 + 0.05*sin(pi*x[0]) else
> False
>
>
> # mark the subdomains with numbers
> subdomain0 = Biomass()
> subdomain0.mark(subdomains, 0)
> subdomain1 = Interface()
> subdomain1.mark(subdomains, 1)
> subdomain2 = Bulk()
> subdomain2.mark(subdomains, 2)

You're actually using SubDomain.mark() method! Did you read
https://bugs.launchpad.net/dolfin/+bug/1119803 and
https://answers.launchpad.net/dolfin/+question/221327 ? One workaround
to this bug could be to initialize

subdomains.set_all(default_value)

where default_value is the most convenient value for cells which don't
belong to any subdomain.

If one default value is not acceptable you need to create more
overlapping subdomains which would entirely cover whole mesh and apply
them in needed order.

Jan

>
> # define functionspace for m
> V0 = FunctionSpace(mesh, "DG", 0)
> m = Function(V0)
>
> # attribute m values to appropriate subdomains
> m_values = [1, 0, 0] # values of m in the two subdomains
> for cell_no in range(len(subdomains.array())):
> subdomain_no = subdomains.array()[cell_no]
> m.vector()[cell_no] = m_values[subdomain_no]
>

Revision history for this message
Hossein Nili (hossein-nili) said :
#9

Thank you. Yes, sorry, my bad in seeing the problem cause:

The error message I receive relates to the boundary condition definition, which goes back to my previous problem emerging again I guess (https://answers.launchpad.net/fenics/+question/228021). I have copied the problematic bit. As you see, I am trying to impose a Dirichlet BC on the subdomain boundary (i.e. interior boundary). When I don't use "pointwise" I get the 'Found no facets matching domain for boundary condition' message. When I do use "pointwise", the error message disappears, I get a solution for my PDE, but it's equal to 0.000 everywhere in the problem domain. I'm happy to copy-paste the variational problem definition part of the code too, but I thought that's not the issue.

Can you please help with this? Is this to be dealt with using the subdomains.set_all(default_value) command?

Many thanks,

# define functionspace for variational problem
VS = FunctionSpace(mesh, "Lagrange", 1)

# define Dirichlet BC
tol = 1E-14 # tolerance for coordinate comparisons
class TopBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[1] - (0.2 + 0.05*sin(pi*x[0]))) < tol

Domain_top_BCS = DirichletBC(VS, Constant(1.0), TopBoundary(), method="pointwise")

Revision history for this message
Best Jan Blechta (blechta) said :
#10

On Wed, 08 May 2013 16:26:20 -0000
Hossein Nili <email address hidden> wrote:
> Question #228027 on FEniCS Project changed:
> https://answers.launchpad.net/fenics/+question/228027
>
> Status: Answered => Open
>
> Hossein Nili is still having a problem:
> Thank you. Yes, sorry, my bad in seeing the problem cause:
>
> The error message I receive relates to the boundary condition
> definition, which goes back to my previous problem emerging again I
> guess (https://answers.launchpad.net/fenics/+question/228021). I have
> copied the problematic bit. As you see, I am trying to impose a
> Dirichlet BC on the subdomain boundary (i.e. interior boundary). When
> I don't use "pointwise" I get the 'Found no facets matching domain for
> boundary condition' message. When I do use "pointwise", the error
> message disappears, I get a solution for my PDE, but it's equal to
> 0.000 everywhere in the problem domain. I'm happy to copy-paste the
> variational problem definition part of the code too, but I thought
> that's not the issue.

I would bet that "no facets found" is the same issue as with marking
cells - i.e. when whole facet (i.e. all its vertices and possibly
barycenter) doesn't fullfill the inside condition (which is quite
tight) facet is not marked. So you could also try to increase tolerance
many orders of magnitude - say to
  tol = 1.e-3
But check that subdomains CellFunction is also set as expected (you
can do plot(subdomains)).

>
> Can you please help with this? Is this to be dealt with using the
> subdomains.set_all(default_value) command?

That probably does not make sense as you have three different regions.

You also could most conveniently write your own marking functions:

subdomains = CellFunction(mesh)
for c in cells(mesh):
    subdomains[c] = (here your decision criterion)

facet_domains = FacetFunction(mesh)
for f in facets(mesh):
    domains = []
    for c in cells(f):
        domains.append(subdomains[c])
    if max(domains) != min(domains):
        print "this is boundary between 2 regions"
        facet_domains[f] = max(domains)*10 + min(domains)
    else:
        print "this interior facet"
        facet_domains[f] = 0

bc = DiricheltBC(V, 0.0, facet_domains, index_you_need)

This ensures that all entries in subdomains and facet_domains are set
and algorithm for that is in your hand.

Jan

>
> Many thanks,
>
>
> # define functionspace for variational problem
> VS = FunctionSpace(mesh, "Lagrange", 1)
>
> # define Dirichlet BC
> tol = 1E-14 # tolerance for coordinate comparisons
> class TopBoundary(SubDomain):
> def inside(self, x, on_boundary):
> return abs(x[1] - (0.2 + 0.05*sin(pi*x[0]))) < tol
>
> Domain_top_BCS = DirichletBC(VS, Constant(1.0), TopBoundary(),
> method="pointwise")
>

Revision history for this message
Hossein Nili (hossein-nili) said :
#11

Thanks very much!

The changing of tol to 1E-3 solved my problem for the time being...as I get along updating the boundary with new functions, I'm quite sure I will get to a point where I need to define my own marking functions which should provide a good way out of the issue. I didn't know of the syntax for doing so; thanks for posting it.

Thanks a lot for your help,

Hossein

Revision history for this message
Hossein Nili (hossein-nili) said :
#12

Thanks Jan Blechta, that solved my question.