Hi.
I'm trying to solve an adaptive poisson system, based on the demo from the FEniCS website.
I replaced the mesh and the boundary conditions. But it failed to adapt the markers:
*** Error: Unable to adapt markers.
*** Reason: Unable to extract information about parent mesh entities.
*** Where: This error was encountered inside adapt.cpp.
*** Process: 0
*** DOLFIN version: 1.6.0
The code:
from dolfin import *
from mshr import *
# Create mesh and define function space
rec1 = Rectangle(Point(0,0),Point(1,1))
rec2 = Rectangle(Point(0.5,0.5),Point(1,1))
mesh = generate_mesh(rec1-rec2,10)
V = FunctionSpace(mesh, "Lagrange", 1)
# Define boundary condition
u0 = Function(V)
bc = DirichletBC(V, u0, "x[0] < DOLFIN_EPS || x[0] > 1.0 - DOLFIN_EPS")
# Create classes for defining parts of the boundaries
class Mid(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.5)
mid = Mid()
# Initialize mesh function for boundary domains
boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)
mid.mark(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)",
degree=1)
g = Expression("sin(5*x[0])", degree=1)
a = inner(grad(u), grad(v))*dx()
L = f*v*dx() + g*v*ds(1)
# Define function for the solution
u = Function(V)
# Define goal functional (quantity of interest)
M = u*dx()
# Define error tolerance
tol = 1.e-5
# Solve equation a = L with respect to u and the given boundary
# conditions, such that the estimated error (measured in M) is less
# than tol
problem = LinearVariationalProblem(a, L, u, bc)
solver = AdaptiveLinearVariationalSolver(problem, M)
solver.parameters["error_control"]["dual_variational_solver"]["linear_solver"] = "cg"
solver.parameters["error_control"]["dual_variational_solver"]["symmetric"] = True
solver.solve(tol)
solver.summary()
# Plot solution(s)
plot(u.root_node(), title="Solution on initial mesh")
plot(u.leaf_node(), title="Solution on final mesh")
interactive()
I appreciate any help you can provide.
Best, Christian.
