Hello,
I am using fenics to solve a linear elastic problem with dirichlet boundary conditions.
I want to mimimize the functionnal :
integral_omega (grad(u) -grad(ud):B:grad(u) -grad(ud))
under the condition div (Bgrad(u)) =0 et u=0 on the bottom and u=u0 on the top.
The solution u obtained does not satisfy the top dirichlet boundary condition.
When I solve the problem :
div (Bgrad(u)) =0 et u=0 on the bottom and u=u0 on the top.
The solution ud seems good.
And I need your help.....
Thanks
1. Import dependencies
from dolfin import *
from petsc4py import *
import numpy as np
import scipy.optimize as opt
import time
import logging
import matplotlib.pyplot as plt
import nb
import math
start = time.clock()
logging.getLogger('FFC').setLevel(logging.WARNING)
logging.getLogger('UFL').setLevel(logging.WARNING)
set_log_active(False)
np.random.seed(seed=1)
2. Model set up:
load mesh and define function spaces for u
mesh_f = Mesh("../Composite_Gmsh/composite_g.xml.gz")
subdomains_f = MeshFunction("size_t", mesh_f, "../Composite_Gmsh/composite_g_physical_region.xml")
boundaries_f = MeshFunction("size_t", mesh_f, "../Composite_Gmsh/composite_g_facet_region.xml")
#
Vu = VectorFunctionSpace(mesh_f, 'Lagrange', 1)
Vu2 = Vu*Vu
Vu2 = MixedFunctionSpace([Vu,Vu])
Va = FunctionSpace(mesh_f, 'DG', 0)
Stress space
Project and write stress field to post-processing file
VS = TensorFunctionSpace(mesh_f, 'DG', 0)
load mesh and define function spaces for a
mesh_g = Mesh("../Composite_Gmsh/composite_f.xml.gz")
subdomains_g = MeshFunction("size_t", mesh_g, "../Composite_Gmsh/composite_f_physical_region.xml")
boundaries_g = MeshFunction("size_t", mesh_g, "../Composite_Gmsh/composite_f_facet_region.xml")
#
Wu = VectorFunctionSpace(mesh_g, 'Lagrange', 1)
Wa = FunctionSpace(mesh_g, 'DG', 0)
Stress space
Project and write stress field to post-processing file
WS = TensorFunctionSpace(mesh_g, 'DG', 0)
The true and inverted parameter
Elasticity parameters
E1 = 1.0e9
nu1 = 0.3
G1 = E1/(2.0(1.0 + nu1))
K1 = E1/(3.0(1.0 - 2.0nu1))
E2 = 5.0e9
nu2 = 0.3
G2 = E2/(2.0(1.0 + nu2))
K2 = E2/(3.0(1.0 - 2.0nu2))
Ktrue = Function(Va)
Gtrue = Function(Va)
K_values = [K1, K2] # values of k in the two subdomain
help = np.asarray(subdomains_f.array(), dtype=np.int32)
Ktrue.vector()[:] = np.choose(help, K_values)
K = interpolate(Expression("2.0e9"),Va)
K = Ktrue
K_g = project(K,Wa)
#
G_values = [G1, G2] # values of k in the two subdomain
help = np.asarray(subdomains_f.array(), dtype=np.int32)
Gtrue.vector()[:] = np.choose(help, G_values)
G = interpolate(Expression("2.0e9"),Va)
G = Gtrue
G_g = project(G,Wa)
Stress computation
def sigma(v):
gdim = v.geometric_dimension()
return (3.0K-2.0G)/3.0 * tr(sym(grad(v)))Identity(gdim) + 2.0G * sym(grad(v))
Stress computation
def sigma_g(v):
gdim = v.geometric_dimension()
return (3.0K_g-2.0G_g)/3.0 * tr(sym(grad(v)))Identity(gdim) + 2.0G_g * sym(grad(v))
Stress computation
def sigma_true(v):
gdim = v.geometric_dimension()
return (3.0Ktrue-2.0Gtrue)/3.0 * tr(sym(grad(v)))Identity(gdim) + 2.0Gtrue * sym(grad(v))
define function for state and adjoint
u = Function(Vu)
p = Function(Vu)
utrue = Function(Vu)
ptrue = Function(Vu)
define Trial and Test Functions
t_trial = TrialFunction(Vu2)
(u_trial,p_trial)=split(t_trial)
t_test = TestFunction(Vu2)
(u_test,p_test)=split(t_test)
K_trial = TrialFunction(Va)
G_trial = TrialFunction(Va)
initialize input functions
f = Constant((0.0, 0.0))
u0 = Constant((0.0, 0.3))
u00 = Constant((0.0, 0.0))
plot
plt.figure(figsize=(15,5))
nb.plot(mesh_g,subplot_loc=121, mytitle="Mesh", show_axis='on')
nb.plot(Ktrue,subplot_loc=122, mytitle="True parameter field")
set up dirichlet boundary conditions
bc_state_1 = DirichletBC(Vu2.sub(0), u0, boundaries_f, 2)
bc_state_2 = DirichletBC(Vu2.sub(0), u00, boundaries_f, 1)
bc_adj_1 = DirichletBC(Vu2.sub(1), u00, boundaries_f, 2)
bc_adj_2 = DirichletBC(Vu2.sub(1), u00, boundaries_f, 1)
bcs = [bc_state_1, bc_state_2, bc_adj_1, bc_adj_2]
bcs = [bc_state_2, bc_adj_1, bc_adj_2]
3. Set up synthetic observations:
ut_trial = TrialFunction(Vu)
ut_test = TestFunction(Vu)
bc_1 = DirichletBC(Vu, u0, boundaries_f, 2)
bc_2 = DirichletBC(Vu, u00, boundaries_f, 1)
weak form for setting up the synthetic observations
at_goal = inner( sigma_true(ut_trial), grad(ut_test)) * dx
Lt_goal = inner(f, ut_test) * dx
solve the forward/state problem to generate synthetic observations
goalt_A, goalt_b = assemble_system(at_goal, Lt_goal, [bc_1, bc_2])
utrue = Function(Vu)
solve(goalt_A, utrue.vector(), goalt_b)
ud = Function(Vu)
ud.assign(utrue)
ud_g = project(ud, Wu)
plot
plt.figure(figsize=(15,5))
nb.plot(Ktrue,subplot_loc=121, mytitle="K_ini", vmin=Ktrue.vector().min(), vmax=Ktrue.vector().max())
nb.plot(utrue,subplot_loc=122, mytitle="u(a_ini)")
4. Solving the global system
define parameters for the optimization
tol = 1e-4
maxiter = 500
initialize iter counters
iter = 1
converged = False
initializations
K_prev = Function(Va)
G_prev = Function(Va)
K_prev_g = Function(Wa)
G_prev_g = Function(Wa)
K_prev.assign(K)
K_prev_g.assign(K_g)
G_prev.assign(G)
G_prev_g.assign(G_g)
u_prev = Function(Vu)
a1 = Function(Wa)
a2 = Function(Wa)
grad_u_g = grad(ud_g)
print "Nit Norme(a - aprec)"
Preconditionner
m = inner(grad(u_trial), grad(u_test))dx + inner(grad(p_trial), grad(p_test))dx
L_goal = dot(f, p_test) * dx + inner( sigma(ud), grad(u_test)) * dx
(M, _) = assemble_system(m, L_goal, bcs)
4. 1st step on inverse problem
weak form for setting up the synthetic observations
A_goal = inner( sigma(u_trial), sym(grad(p_test))) * dx + inner( sigma(u_trial), sym(grad(u_test))) * dx + inner( sigma(p_trial), sym(grad(u_test))) * dx
L_goal = dot(f, p_test) * dx + inner( sigma(ud), sym(grad(u_test))) * dx
solve the forward/state problem to generate synthetic observations
goal_A, goal_L = assemble_system(A_goal, L_goal, bcs)
usol = Function(Vu2)
fenics iterative solve
solver = KrylovSolver("tfqmr", "amg")
solver.set_operators(goal_A, M)
solver.solve(usol.vector(), goal_L)
(u, p) = usol.split(deepcopy=True)
solve minimization on a
u_g = project(u, Wu)
plot(u.sub(1))
interactive()
