1. FEniCS Project Avatar FEniCS Project
  2. Core
  3. DOLFIN
  4. Issues

solving on submesh may lead to wrong results

Issue #88 resolved
Nico Schlömer created an issue

The following code aims at solving Poisson's equation on a subdomain of a given domain, but due to a mixup of meshes produces an incorrect solution. When compiled in Developer mode, Dolfin intercepts the mixup, but doesn't in all other cases.

See discussion on http://fenicsproject.org/qa/544/solving-pde-on-submesh.

from dolfin import *

mesh = UnitSquareMesh(20, 20, 'left/right')

# -------------------------------------------------
class OmegaLeft(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] <= 0.5
omegaL = OmegaLeft()

# Initialize mesh function for interior domains
subdomains = CellFunction('size_t', mesh)
subdomains.set_all(0)
omegaL.mark(subdomains, 1)

dx = Measure('dx')[subdomains]

submesh_omegaL = SubMesh(mesh, subdomains, 1)
# -------------------------------------------------
# define boundaries for the subdomain
class LeftBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and abs(x[0]) < DOLFIN_EPS
left_boundary = LeftBoundary()

omegaL_boundaries = FacetFunction('size_t', submesh_omegaL)
omegaL_boundaries.set_all(0)
left_boundary.mark(omegaL_boundaries, 1)

ds_omegaL = Measure('ds')[omegaL_boundaries]
# -------------------------------------------------
# Solve Poisson's equation.
f = Constant(0.0)

V = FunctionSpace(submesh_omegaL, 'CG', 1)
u = TrialFunction(V)
v = TestFunction(V)
# Boundary conditions of the form
#   n.grad(u) = alpha * (u0 - u)
alpha = 1.0
u0 = 0.421
F = dot(grad(u), grad(v)) * dx(1) \
  - alpha * (u0-u) * v * ds_omegaL(1) \
  - f * v * dx(1)

a = lhs(F)
L = rhs(F)

sol = Function(V)
solve(a == L, sol)

plot(sol)
interactive()
# -------------------------------------------------

Comments (5)

  2. Jan Blechta

    Correct version of this code is

    from dolfin import *

mesh = UnitSquareMesh(20, 20, 'left/right')

# -------------------------------------------------
class OmegaLeft(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] <= 0.5
omegaL = OmegaLeft()

# Initialize mesh function for interior domains
subdomains = CellFunction('size_t', mesh)
subdomains.set_all(0)
omegaL.mark(subdomains, 1)

dx = Measure('dx')[subdomains]

submesh_omegaL = SubMesh(mesh, subdomains, 1)
# -------------------------------------------------
# define boundaries for the subdomain
class LeftBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and abs(x[0]) < DOLFIN_EPS
left_boundary = LeftBoundary()

omegaL_boundaries = FacetFunction('size_t', submesh_omegaL)
omegaL_boundaries.set_all(0)
left_boundary.mark(omegaL_boundaries, 1)

ds_omegaL = Measure('ds')[omegaL_boundaries]

mydomains = CellFunction('size_t', submesh_omegaL)
mydomains.set_all(0)
dx_subdomain = Measure('dx')[mydomains]
# -------------------------------------------------
# Solve Poisson's equation.
f = Constant(1.0)

V = FunctionSpace(submesh_omegaL, 'CG', 1)
u = TrialFunction(V)
v = TestFunction(V)
# Boundary conditions of the form
#   n.grad(u) = alpha * (u0 - u)
alpha = 1.0
u0 = 0.421
F = dot(grad(u), grad(v)) * dx_subdomain(0) \
  - alpha * (u0-u) * v * ds_omegaL(1) \
  - f * v * dx_subdomain(0)

a = lhs(F)
L = rhs(F)

sol = Function(V)
solve(a == L, sol)

plot(sol)
interactive()
# -------------------------------------------------

    This yields constant solution u==0.421. With f = Constant(1.0) non-trivial solution is obtained and u-u0==0.5 is fulfilled on Robin boundary as required by necessary condition of solvability of Neumann-Poisson problem.

    The actual bug is that with incorrect Nico's version of the code user is not warned when assertions are turned off and he is computing rubbish.

  5. Bairav Sabarish

    Hello, I’m trying to solve a poisson equation with non-linear boundary conditions in one subdomain of the domain.

    First, I tried solving the poisson equation in the full domain to check if it is working. The code works fine.

    But when i try to solve it only in a subdomain, i’m not able to get it to work. There is even no change in the geometry of the domains in the first and second case. I’m making some mistake while integrating the poisson solver with the subdomain code shared above.

    I have attached the working code in a single domain first, and then the code to solve in in only one subdomain (not working).

    Please help me with this issue. I’m not able to resolve it.

    Thanks in advance ! @Jan Blechta @Prof Garth Wells 

    from __future__ import print_function
from dolfin import *
import matplotlib.pyplot as plt
# from __future__ import print_function
from fenics import *
from mshr import *
import numpy as np
from ufl import nabla_div
from ufl import sinh
import random


domain = Rectangle(Point(0,0), Point(100e-6, 200e-6))
mesh = generate_mesh(domain,100)
V = FunctionSpace(mesh, 'P', 1)

f = Constant(0.0)

i0 = 12.6
r = Constant(2*i0)

iApp = -100
g = iApp

F = 96485.33
R = 8.314
T = 300
C = Constant(F/(2*R*T))

kappa = 1


class BottomBoundary(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14
        return on_boundary and abs(x[1]) < tol

class TopBoundary(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14
        return on_boundary and abs(x[1] - 200e-6) < tol

bottom_boundary = BottomBoundary()        
top_boundary = TopBoundary()   

boundary_parts = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_parts.set_all(0)
bottom_boundary.mark(boundary_parts, 1)
top_boundary.mark(boundary_parts, 2)
ds = Measure("ds")[boundary_parts]


u = TrialFunction(V)
v = TestFunction(V)

u = Function(V)

bcs = []

F = kappa*inner(nabla_grad(u), nabla_grad(v))*dx + g*v*ds(2) + r*sinh(C*u)*v*ds(1) - f*v*dx
J = derivative(F, u)
problem = NonlinearVariationalProblem(F, u, bcs, J)
solver  = NonlinearVariationalSolver(problem)
solver.solve()

import matplotlib.pyplot as plt
p = plot(u, mode='color')
plt.colorbar(p)
plt.title(r"$\Phi$",fontsize=26)
plt.show()

    ‌

    from __future__ import print_function
from dolfin import *
import matplotlib.pyplot as plt
# from __future__ import print_function
from fenics import *
from mshr import *
import numpy as np
from ufl import nabla_div
from ufl import sinh
import random


domain = Rectangle(Point(0,0), Point(100e-6, 400e-6))
mesh = generate_mesh(domain,100)

f = Constant(0.0)

i0 = 12.6
r = Constant(2*i0)

iApp = -100
g = iApp

F = 96485.33
R = 8.314
T = 300
C = Constant(F/(2*R*T))

kappa = 1

# -------------------------------------------------
class OmegaMiddle(SubDomain):
    def inside(self, x, on_boundary):
        return x[1] <= 200e-6
omegaM = OmegaMiddle()

# Initialize mesh function for interior domains
subdomains = MeshFunction('size_t', mesh, 2)
subdomains.set_all(0)
omegaM.mark(subdomains, 1)

dx = Measure('dx')[subdomains]

submesh_omegaL = SubMesh(mesh, subdomains, 1)
# -------------------------------------------------
# define boundaries for the subdomain
class BottomBoundary(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14
        return on_boundary and abs(x[1]) < tol
Bottom_boundary = BottomBoundary()

class TopBoundary(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14
        return on_boundary and abs(x[1]-200e-6) < tol
Top_boundary = TopBoundary()

omegaL_boundaries = MeshFunction('size_t', submesh_omegaL,1)
omegaL_boundaries.set_all(0)
Bottom_boundary.mark(omegaL_boundaries, 1)
Top_boundary.mark(omegaL_boundaries, 2)

ds_omegaL = Measure('ds')[omegaL_boundaries]

mydomains = MeshFunction('size_t', submesh_omegaL,2)
mydomains.set_all(0)
dx_subdomain = Measure('dx')[mydomains]
# -------------------------------------------------
# Solve Poisson's equation.

V = FunctionSpace(submesh_omegaL, 'P', 1)
u = TrialFunction(V)
v = TestFunction(V)

u = Function(V)

bcs = []

F = kappa*inner(nabla_grad(u), nabla_grad(v))*dx_subdomain(0) + g*v*ds_omegaL(2) + r*sinh(C*u)*v*ds_omegaL(1) - f*v*dx_subdomain(0)

J = derivative(F, u)
problem = NonlinearVariationalProblem(F, u, bcs, J)
solver  = NonlinearVariationalSolver(problem)
solver.solve()

import matplotlib.pyplot as plt
p = plot(u, mode='color')
plt.colorbar(p)
# plt.title(r"$\Phi$",fontsize=26)
plt.show()

    ‌

  6. Log in to comment
Assignee
–
Type
bug
Priority
minor
Status
resolved
Component
–
Milestone
–
Version
–
Votes
0
Watchers
3
Jira: the preferred issue tracker for Bitbucket. Join the team!