Hi all,
I am perseverating in my effort to port this notebook to Fenicsx. There are still numberless rough edges (e.g. as I understand there is no LocalSolver in Fenics I need to do code the projection myself, I think it is just a matter of a simple weak form but just to get it going I am brutally interpolating for now). Also, I am very unsure on how to translate assemble from Dolfin, and eventually I found some posts on the same error suggesting compliance = fem.assemble_scalar(fem.form(ufl.action(L, uh)) ) . I honestly have not fully understood at all the necessity for fem.form , but at least I do not get a syntax error, the sooner I get it at least running I am hopeful I can try to make more sense of it.
Anyhow, I got close to the end without syntax errors at least, but I now get an error which I cannot get an idea about. This is the code, I have to attach it all including the optimisation loop as the simple elastic problem works solved over one iterations works fine
from dolfinx import *
import matplotlib.pyplot as plt
import numpy as np
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
from mpi4py import MPI
from dolfinx.fem import FunctionSpace, Function
from dolfinx.fem import VectorFunctionSpace
from dolfinx import fem, io
from dolfinx import mesh
from dolfinx.mesh import CellType, create_rectangle
#from dolfinx.la import solver
import ufl
# Algorithmic parameters
niternp = 20 # number of non-penalized iterations
niter = 80 # total number of iterations
pmax = 4 # maximum SIMP exponent
exponent_update_frequency = 4 # minimum number of steps between exponent update
tol_mass = 1e-4 # tolerance on mass when finding Lagrange multiplier
thetamin = 0.001 # minimum density modeling void
# Mesh
L = 5
domain = create_rectangle(MPI.COMM_WORLD, np.array([[0,0],[3*L, L]]), [30,30], cell_type=CellType.triangle)
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(domain.topology)
# Problem parameters
thetamoy = 0.4 # target average material density
E = fem.Constant(domain, ScalarType(1))
nu = fem.Constant(domain, ScalarType(0.3))
lamda = E*nu/(1+nu)/(1-2*nu)
mu = E/(2*(1+nu))
# Traction and body forces
f = fem.Constant(domain, ScalarType([0,-1])) # vertical downwards force
body = fem.Constant(domain, ScalarType((0, 0.1)))
# Function space for density field
V0 = fem.FunctionSpace(domain, ("DG", 0))
# Function space for displacement
V2 = fem.VectorFunctionSpace(domain, ("CG", 2))
# Definition for boundary conditions
def clamped_boundary(x):
return np.isclose(x[0], 0)
fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(domain, fdim, clamped_boundary)
u_D = np.array([0,0], dtype=ScalarType)
bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V2, fdim, boundary_facets), V2)
def left(x):
return np.isclose(x[0], 0)
def right(x):
return np.isclose(x[0], 3*L)
fdim = domain.topology.dim -1
left_facets = mesh.locate_entities_boundary(domain, fdim, left)
right_facets = mesh.locate_entities_boundary(domain, fdim, right)
# Concatenate and sort the arrays based on facet indices. Left facets marked with 1, right facets with two
marked_facets = np.hstack([left_facets, right_facets])
marked_values = np.hstack([np.full_like(left_facets, 1), np.full_like(right_facets, 2)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])
metadata = {"quadrature_degree": 4}
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag, metadata=metadata)
dx = ufl.Measure("dx", domain=domain, metadata=metadata)
# SIMP parameters
p = fem.Constant(domain, ScalarType(1)) # SIMP penalty exponent
exponent_counter = 0 # exponent update counter
lagrange = fem.Constant(domain, ScalarType(1))# Lagrange multiplier for volume constraint
thetaold.x.set(thetamoy)
coeff = thetaold**p
theta = Function(V0)
volume = fem.assemble_scalar(fem.form(thetaold*ufl.dx))
avg_density_0 = fem.assemble_scalar(fem.form(thetaold*ufl.dx))/volume
avg_density = 0.
thetaold = Function(V0, name="Density")
# Functions definitions for weak form
def eps(v):
return ufl.sym(ufl.grad(v))
def sigma(v):
return coeff*(lamda*ufl.nabla_div(v)*ufl.Identity(2)+2*mu*eps(v))
# return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2*mu*epsilon(u)
def energy_density(u, v):
return ufl.inner(sigma(u), eps(v))
# Inhomogeneous elastic variational problem
u_ = ufl.TestFunction(V2)
du = ufl.TrialFunction(V2)
#a = ufl.inner(sigma(u_), eps(du))*ufl.dx
a = ufl.inner(sigma(u_), eps(du)) * ufl.dx
L = ufl.dot(f, u_)*ds(2)
# Functions for handling SIMP density field
def local_project(v, V):
dv = ufl.TrialFunction(V)
v_ = ufl.TestFunction(V)
a_proj = ufl.inner(dv, v_)*dx
b_proj = ufl.inner(v, v_)*dx
u_proj = fem.Function(V)
u_proj.interpolate(v)
# solver.solve_local_rhs(u)
return u_proj
def update_theta():
theta_new = Function(V0)
theta.x.array[:] = local_project((p*coeff*energy_density(u, u)/lagrange)**(1/(p+1)), V0).x.array
# theta.assign(local_project((p*coeff*energy_density(u, u)/lagrange)**(1/(p+1)), V0))
thetav = theta.vector().get_local()
theta.vector().set_local(np.maximum(np.minimum(1, thetav), thetamin))
theta.vector().apply("insert")
# avg_density = fem.assemble(theta*dx)/volume
avg_density = fem.assemble_scalar(fem.form(theta*ufl.dx))/volume
return avg_density
# Update Lagrange multiplier and SIMP penalty exponent
def update_lagrange_multiplier(avg_density):
avg_density1 = avg_density
# Initial bracketing of Lagrange multiplier
if (avg_density1 < avg_density_0):
lagmin = float(lagrange)
while (avg_density < avg_density_0):
lagrange.assign(fem.Constant(domain, ScalarType(lagrange/2)))
avg_density = update_theta()
lagmax = float(lagrange)
elif (avg_density1 > avg_density_0):
lagmax = float(lagrange)
while (avg_density > avg_density_0):
lagrange.assign(Constant(fem.Constant(domain, ScalarType(lagrange*2))))
avg_density = update_theta()
lagmin = float(lagrange)
else:
lagmin = float(lagrange)
lagmax = float(lagrange)
# Dichotomy on Lagrange multiplier
inddico=0
while ((abs(1.-avg_density/avg_density_0)) > tol_mass):
lagrange.assign(Constant((lagmax+lagmin)/2))
avg_density = update_theta()
inddico += 1;
if (avg_density < avg_density_0):
lagmin = float(lagrange)
else:
lagmax = float(lagrange)
print(" Dichotomy iterations:", inddico)
def update_exponent(exponent_counter):
exponent_counter += 1
if (i < niternp):
p.assign(fem.Constant(domain, ScalarType(1)))
elif (i >= niternp):
if i == niternp:
print("\n Starting penalized iterations\n")
if ((abs(compliance-old_compliance) < 0.01*compliance_history[0]) and
(exponent_counter > exponent_update_frequency) ):
# average gray level
gray_level = fem.assemble_scalar((theta-thetamin)*(1.-theta)*dx)*4/volume
# p.assign(Constant(min(float(p)*(1+0.3**(1.+gray_level/2)), pmax)))
p.assign(fem.Constant(domain, ScalarType (min(float(p)*(1+0.3**(1.+gray_level/2)), pmax))))
fem.Constant(domain, ScalarType(0.3))
exponent_counter = 0
print(" Updated SIMP exponent p = ", float(p))
return exponent_counter
#### Main Iteration
u = Function(V2, name="Displacement")
old_compliance = 1e30
compliance_history = []
for i in range(niter):
problem = fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
compliance = fem.assemble_scalar(fem.form(ufl.action(L, uh)) )
##################
compliance_history.append(compliance)
print("Iteration {}: compliance =".format(i), compliance)
avg_density = update_theta()
update_lagrange_multiplier(avg_density)
exponent_counter = update_exponent(exponent_counter)
# Update theta field and compliance
thetaold.assign(theta)
old_compliance = compliance
I get the error message, after solving the first elastic iteration
Couldn't map 'c_4' to a float, returning ufl object without evaluation.
I cannot figure it out what it could refer to, the traceback does not even point to the problematic code line. I got the same error earlier during the day when trying to understand how to interpolate a constant function, and I implemented the solution above, but in this case it is very unclear where it comes from.
Thanks so much