### Why do I get this error when computing an auto-adaptive method for the mixed Poisson equation?

347
views
0
10 months ago by
Hello!

I am currently working on mixed formulations for the linear elasticity equations and I wish to implement an adaptive scheme. I tried the adaptive solver implemented in FEniCS and got some errors which I do not understand. I obtain the same errors when solving the mixed Poisson equation. I have constructed a short working example based on the mixed formulation for the Poisson equation given in this example.

The error message that pops up is the following one: "Replacement expressions must have the same shape as what they replace". What is the source of this error?

from dolfin import *

# Create mesh
mesh = UnitSquareMesh(32, 32)

# Define function spaces and mixed (product) space
BDM = FiniteElement("BDM", mesh.ufl_cell(), 1)
DG = FiniteElement("DG", mesh.ufl_cell(), 0)
W = FunctionSpace(mesh, BDM * DG)

# Define trial and test functions
(sigma, u) = TrialFunctions(W)
(tau, v) = TestFunctions(W)
w = Function(W)

# Define source function
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)", degree = 4)

# Define variational form
a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx
L = - f*v*dx

# Define function G such that G \cdot n = g
class BoundarySource(Expression):
def __init__(self, mesh, degree):
self.mesh = mesh
self.degree = degree
def eval_cell(self, values, x, ufc_cell):
cell = Cell(self.mesh, ufc_cell.index)
n = cell.normal(ufc_cell.local_facet)
g = sin(5*x[0])
values[0] = g*n[0]
values[1] = g*n[1]
def value_shape(self):
return (2,)

G = BoundarySource(mesh, degree = 4)

# Define essential boundary
def boundary(x):
return x[1] < DOLFIN_EPS or x[1] > 1.0 - DOLFIN_EPS

bc = DirichletBC(W.sub(0), G, boundary)

# Define goal functional (quantity of interest)
M = w[2]*dx

# Define error tolerance
tol = 1.e-3

# Solve equation a = L with respect to w and the given boundary
# conditions, such that the estimated error (measured in M) is less
# than tol
problem = LinearVariationalProblem(a, L, w, bc)
solver.parameters["error_control"]["dual_variational_solver"]["linear_solver"] = "cg"
solver.parameters["error_control"]["dual_variational_solver"]["symmetric"] = True
solver.solve(tol)

solver.summary()​

Best regards,
Gonzalo

** Dolfin version: 2016.2.0
Community: FEniCS Project
1
Seems like something odd with the function space. Changing it to, e.g.:
BDM = VectorElement('P', mesh.ufl_cell(), 1)
DG  = FiniteElement('DG', mesh.ufl_cell(), 0)​

the error is gone (not that it is converging).

I don't know what is the problem but you might also check out the undocumented auto-adaptive-navier-stokes demo.

written 10 months ago by Adam Janecka

I have noticed that if I choose the piecewise polynomial vector finite element space ''P'', the error no longer appears but the solver does not converge because this combination of finite element spaces does not satisfy the inf-sup stability conditions.

Regarding the auto-adaptive Navier-Stokes demo, in this case they once again use piecewise polynomial finite element spaces. Hence, it seems that the problem is related to the use of Hdiv-conforming finite elements, like 'RT' or 'BDM' (and these are exactly the ones I want to use).
written 10 months ago by Gonzalo González de Diego