### weighted residuals evaluation

106
views
0
3 months ago by
To compute a posteriori error estimators it is needed to evaluate weighted face/element residuals. Before I used:
h_faces = FacetArea(mesh)
h_elems = CellVolume(mesh)


But the problem is that dolfin-adjoint can not derive gradients of these functions.
Is there any effective built-in functions instead?
Thanks!

Community: FEniCS Project
written 3 months ago by Miguel

Yes, I want to derive gradients from error estimator. Actually, I want to implement the topology optimization problem with regularization by a posteriori error estimator. The code looks as follows (it is just the tutorial on linear elasticity + few strings dolfin-adjoint code to derive gradients):

from fenics import *
import pyipopt

# Scaled variables
L = 3; W = 1
mu = 1
rho = 1
delta = W/L
gamma = 0.4*delta**2
lambda_ = 1.25
g = gamma

N=10
# Create mesh and define function space
mesh = BoxMesh(Point(0, 0, 0), Point(W, L, L), 3*N, N, N)
normals = FacetNormal(mesh)
h_faces = FacetArea(mesh)
h_cells = CellVolume(mesh)

V = VectorFunctionSpace(mesh, 'P', 1) # admissible solutions space
P = FunctionSpace(mesh,'DG', 0) # admissible designs space (piece-wise constant)

# initial design (with respect to which gradient should be computed)
x = interpolate(Constant(1), P)

# Define boundary condition
tol = 1E-14

def clamped_boundary(x, on_boundary):
return on_boundary and x[0] < tol

bc = DirichletBC(V, Constant((0, 0, 0)), clamped_boundary)

# Define strain and stress

def epsilon(u):

def sigma(u):
return lambda_*nabla_div(u)*Identity(d) + 2*mu*epsilon(u)

# Define variational problem
u = TrialFunction(V)
d = u.geometric_dimension()  # space dimension
v = TestFunction(V)
f = Constant((0, 0, -rho*g))
T = Constant((0, 0, 0))
a = inner(x*sigma(u), epsilon(v))*dx
L = dot(f, v)*dx + dot(T, v)*ds

# Compute solution
u = Function(V)
solve(a == L, u, bc)

# Define cost functional: compliance + a posteriori error
J = assemble(dot(u,f)*dx + h_cells**2*(f + div(x*sigma(u)))**2*dx + h_faces('+')*jump(x*sigma(u),normals)**2*dS)

dJdx = compute_gradient(J, Control(x))

It works well for h_cells and it used to work well for h_faces('+'), but now it gives the following error:

Solving linear variational problem.
Traceback (most recent call last):
File "t.py", line 58, in <module>
J = assemble(inner(u,f)*dx + h_cells**2*(f + div(x*sigma(u)))**2*dx + h_faces('+')*jump(x*sigma(u),normals)**2*dS)
Abort trap: 6