- edited description
solving on submesh may lead to wrong results
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)
-
reporter -
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
. Withf = Constant(1.0)
non-trivial solution is obtained andu-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.
-
Fix in pull request 48.
-
- changed status to resolved
-
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()
- Log in to comment