Hi all,
I am trying to solve the same problem as the demo here by my unit square generated from Gmsh.
However, I got some different results. Following is the result I got from my own geometry,
And the following is the result I got from demo code,
They are similar but not the same, especially the middle peak's location.
Could I ask whether it is just an numerical error or some bugs inside my code?
Following is my code.
from dolfin import *
# Create mesh and define function space
# mesh = UnitSquareMesh(32, 32)
mesh = Mesh("2d.xml")
V = FunctionSpace(mesh, "Lagrange", 1)
# Sub domain
class DitichletBC1(SubDomain):
def inside(self, x, on_boundary):
return x[0] > 1.0 - DOLFIN_EPS and on_boundary
class DitichletBC2(SubDomain):
def inside(self, x, on_boundary):
return x[0] < DOLFIN_EPS and on_boundary
sub_domains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
sub_domains.set_all(2)
dirichletBC1 = DitichletBC1()
dirichletBC1.mark(sub_domains, 0)
dirichletBC2 = DitichletBC2()
dirichletBC2.mark(sub_domains, 1)
file = File("subdomains.xml")
file << sub_domains
boundaries = MeshFunction("size_t", mesh, "subdomains.xml") ## Input your facet file name
# Define boundary condition
u0 = Constant(0.0)
bc1 = DirichletBC(V, u0, boundaries, 0)
bc2 = DirichletBC(V, u0, boundaries, 1)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)")
g = Expression("sin(5*x[0])")
a = inner(grad(u), grad(v))*dx
L = f*v*dx + g*v*ds
# Compute solution
u = Function(V)
solve(a == L, u, [bc1,bc2])
# Plot solution
plot(u, interactive=True)
Thanks a lot!
