Hi,
I’m having troubles to plot the density field on the deformed mesh. The code comes from the exemple in fenics documentation: Topology optimization using the SIMP method — Numerical tours of continuum mechanics using FEniCS master documentation
Blockquote
from dolfin import *
import matplotlib.pyplot as plt
import numpy as np
import timeAlgorithmic 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 voidProblem parameters
thetamoy = 0.4 # target average material density
E = Constant(30)
nu = Constant(0.2)
lamda = Enu/(1+nu)/(1-2nu)
mu = E/(2*(1+nu))
f_grav = Constant((0, -1e-3)) # vertical downwards force
f_surface=Constant((0,-1e-2))Mesh
mesh = RectangleMesh(Point(-2, 0), Point(2, 1), 25, 6, “crossed”) #200,48
Boundaries
def cote(x, on_boundary):
return near (x[0],2) or near(x[0], -2) and on_boundarydef load(x, on_boundary):
return near(x[1], 1)facets = MeshFunction(“size_t”, mesh, 1)
AutoSubDomain(load).mark(facets, 1)
ds = Measure(“ds”, subdomain_data=facets)
#plt.figure(1)
#plot(mesh,linewidth=0.1)
#plt.pause(10)Function space for density field
V0 = FunctionSpace(mesh, “DG”, 0)
Function space for displacement
V2 = VectorFunctionSpace(mesh, “CG”, 2)
Fixed boundary condtions
bc = DirichletBC(V2, Constant((0, 0)), cote)
p = Constant(1) # SIMP penalty exponent
exponent_counter = 0 # exponent update counter
lagrange = Constant(1) # Lagrange multiplier for volume constraintthetaold = Function(V0, name=“Density”)
thetaold.interpolate(Constant(thetamoy))
coeff = thetaold**p
theta = Function(V0)volume = assemble(Constant(1.)dx(domain=mesh))
avg_density_0 = assemble(thetaolddx)/volume # initial average density
avg_density = 0.def eps(v):
return sym(grad(v))
def sigma(v):
return coeff*(lamdadiv(v)Identity(2)+2mueps(v))
def energy_density(u, v):
return inner(sigma(u), eps(v))Inhomogeneous elastic variational problem
u_ = TestFunction(V2)
du = TrialFunction(V2)
a = inner(sigma(u_), eps(du))*dx
#L = dot(f_surface, u_)*ds(1)
L=inner(f_grav,u_)*dx +dot(f_surface, u_)*ds(1)def local_project(v, V):
dv = TrialFunction(V)
v_ = TestFunction(V)
a_proj = inner(dv, v_)dx
b_proj = inner(v, v_)dx
solver = LocalSolver(a_proj, b_proj)
solver.factorize()
u = Function(V)
solver.solve_local_rhs(u)
return u
def update_theta():
theta.assign(local_project((pcoeffenergy_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 = assemble(theta*dx)/volume
return avg_densitydef 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(Constant(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(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(Constant(1))
elif (i >= niternp):
if i == niternp:
print(“\n Starting penalized iterations\n”)
if ((abs(compliance-old_compliance) < 0.01compliance_history[0]) and
(exponent_counter > exponent_update_frequency) ):
# average gray level
gray_level = assemble((theta-thetamin)(1.-theta)dx)4/volume
p.assign(Constant(min(float(p)(1+0.3*(1.+gray_level/2)), pmax)))
exponent_counter = 0
print(" Updated SIMP exponent p = ", float(p))
return exponent_counteru = Function(V2, name=“Displacement”)
old_compliance = 1e30
ffile = XDMFFile(“topology_optimization.xdmf”)
ffile.parameters[“flush_output”]=True
ffile.parameters[“functions_share_mesh”]=True
compliance_history =tini=time.time()
for i in range(niter):
start=time.time()
solve(a == L, u, bc, solver_parameters={“linear_solver”: “cg”, “preconditioner”: “hypre_amg”})
ffile.write(thetaold, i)
ffile.write(u, i)compliance = assemble(action(L, u)) 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 finish=time.time() print("temps écoulé",finish-start) plt.figure(2) plot(theta, cmap="bone_r") plt.title("Density") plt.show(block=False) plt.pause(0.1)tfin=time.time()
plt.figure(4)
plot(u,mode=“displacement”,cmap=“bone_r”)
plt.title(“déplacements”)
the plot at the end represent the amount of deformation on each element in the deformed mesh but not the density field (which is theta). I tried u*theta but it isn’t working. How can I get the kind of plot shown in the fenics documentation?
Thanks for your help!