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

166
views
0
4 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 4 months ago by Adam Janecka