Reconstruction

import dolfin as df
import math
import matplotlib.pyplot as plt
import numpy as np
import os

from elasticity_problem import elasticity_problem
from noiseless_data_case_rec import noiseless_data_rec

%load_ext autoreload
%autoreload 2

True parameter mu

# mu_1
class Mu1Expression(df.UserExpression):
    def eval(self, value, x):
        value[0]=0
        if (x[0] - 0.5) *(x[0] - 0.5)  + (x[1] - 0.5) *(x[1] - 0.5)  <= 0.2  *0.2 :
            value[0]= 1
        elif  0.2  *0.2 <= (x[0] - 0.5) *(x[0] - 0.5)  + (x[1] - 0.5) *(x[1] - 0.5)  <= 0.3  *0.3 :
            value[0]= (1-(math.sqrt((x[0] - 0.5) *(x[0] - 0.5)  + (x[1] - 0.5) *(x[1] - 0.5)) - 0.2)/0.1 )**2*(1+2*(math.sqrt((x[0] - 0.5)*(x[0] - 0.5)  + (x[1] - 0.5) *(x[1] - 0.5))-0.2) /0.1) 
        value[0] = 1+value[0]

# mu_2
class Mu2Expression(df.UserExpression):
    def eval(self, value, x):
        value[0]=0
        if (x[0] - 0.2) *(x[0] - 0.2)  + (x[1] - 0.8) *(x[1] - 0.8)  <= 0.1  *0.1 :
            value[0]= 1
        elif  0.1  *0.1 <= (x[0] - 0.2) *(x[0] - 0.2)  + (x[1] - 0.8) *(x[1] - 0.8)  <= 0.15  *0.15 :
            value[0]= (1-(math.sqrt((x[0] - 0.2) *(x[0] - 0.2)  + (x[1] - 0.8) *(x[1] - 0.8)) - 0.1)/0.05 )**2*(1+2*(math.sqrt((x[0] - 0.2)*(x[0] - 0.2)  + (x[1] - 0.8) *(x[1] - 0.8))-0.1) /0.05) 
        elif (x[0] - 0.8) *(x[0] - 0.8)  + (x[1] - 0.2) *(x[1] - 0.2)  <= 0.1  *0.1 :
            value[0]= 2
        elif  0.1  *0.1 <= (x[0] - 0.8) *(x[0] - 0.8)  + (x[1] - 0.2) *(x[1] - 0.2)  <= 0.15  *0.15 :
            value[0]=2*(1-(math.sqrt((x[0] - 0.8) *(x[0] - 0.8)  + (x[1] - 0.2) *(x[1] - 0.2)) - 0.1)/0.05 )**2*(1+2*(math.sqrt((x[0] - 0.8)*(x[0] - 0.8)  + (x[1] - 0.2) *(x[1] - 0.2))-0.1) /0.05)     
        value[0] = 1+value[0]

# mu_3
class Mu3Expression(df.UserExpression):
    def eval(self, value, x):
        value[0]=1
        if (x[0] - 0.5) *(x[0] - 0.5)  + (x[1] - 0.5) *(x[1] - 0.5)  <= 0.4   *0.4 :
            value[0]= 2
        elif  0.4  *0.4 <= (x[0] - 0.5) *(x[0] - 0.5)  + (x[1] - 0.5) *(x[1] - 0.5)  <= 0.45  *0.45:
            value[0]= 1+(1-(math.sqrt((x[0] - 0.5) *(x[0] - 0.5)  + (x[1] - 0.5) *(x[1] - 0.5)) - 0.4)/0.05 )**2*(1+2*(math.sqrt((x[0] - 0.5)*(x[0] - 0.5)  + (x[1] - 0.5) *(x[1] - 0.5))-0.4) /0.05) 

        if (x[0] - 0.35) *(x[0] - 0.35)  + (x[1] - 0.65) *(x[1] - 0.65)  <= 0.1   *0.1 :
            value[0]= 3
        elif  0.1  *0.1<= (x[0] - 0.35) *(x[0] - 0.35)  + (x[1] - 0.65) *(x[1] - 0.65)  <= 0.15  *0.15:
            value[0]= 2+(1-(math.sqrt((x[0] - 0.35) *(x[0] - 0.35)  + (x[1] - 0.65) *(x[1] - 0.65)) - 0.1)/0.05 )**2*(1+2*(math.sqrt((x[0] - 0.35)*(x[0] - 0.35)  + (x[1] - 0.65) *(x[1] - 0.65))-0.1) /0.05) 

        if (x[0] - 0.5) *(x[0] - 0.5)  + (x[1] - 0.3) *(x[1] - 0.3)  <= 0.1   *0.1 :
            value[0]= 4  
        elif  0.1  *0.1<= (x[0] - 0.5) *(x[0] - 0.5)  + (x[1] - 0.3) *(x[1] - 0.3) <= 0.15  *0.15:
           value[0]= 2+2*(1-(math.sqrt((x[0] - 0.5) *(x[0] - 0.5)  + (x[1] - 0.3) *(x[1] - 0.3)) - 0.1)/0.05 )**2*(1+2*(math.sqrt((x[0] - 0.5) *(x[0] - 0.5)  + (x[1] - 0.3) *(x[1] - 0.3))-0.1) /0.05)            

        if (x[0] - 0.65) *(x[0] - 0.65)  + (x[1] - 0.65) *(x[1] - 0.65)  <= 0.1   *0.1 :
            value[0]= 3.5       
        elif  0.1  *0.1<= (x[0] - 0.65) *(x[0] - 0.65)  + (x[1] - 0.65) *(x[1] - 0.65)   <= 0.15  *0.15:
            value[0]= 2+1.5*(1-(math.sqrt((x[0] - 0.65) *(x[0] - 0.65)  + (x[1] - 0.65) *(x[1] - 0.65) ) - 0.1)/0.05 )**2*(1+2*(math.sqrt((x[0] - 0.65) *(x[0] - 0.65)  + (x[1] - 0.65) *(x[1] - 0.65) )-0.1) /0.05)     
        value[0] = value[0]

# mu = Mu1Expression()
mu = Mu2Expression()
# mu = Mu3Expression()

Problem parameters

omega = 1
rho = 1
degree = 3
dim = 2
Nx = 200
Ny = 200
u_boundary = df.Expression(('x[1]', 0), degree=1)
f = df.Constant((0,0))
tol = 1e-15

mesh_params = {"Nx":Nx, "Ny":Ny, "degree":degree}
load_params = {"u_boundary":u_boundary, "f":f} 
mat_params  = {"mu":mu, "rho":rho, "omega":omega}

Direct problem

u, eps_12, mu_function = elasticity_problem(
    mesh_params=mesh_params,
    mat_params=mat_params,
    load_params=load_params)

Figures 1 and 2

#Plot u 
u1, u2 = u.split()
plt.subplot(2,2,1)
p = df.plot(u1, title="u_x", cmap="inferno")
plt.colorbar(p) 

plt.subplot(2,2,2)
p = df.plot(u2, title="u_y", cmap="inferno")
plt.colorbar(p) 

plt.subplot(2,2,3)
p = df.plot(mu_function, title="mu", cmap="inferno")
plt.colorbar(p) 


plt.subplot(2,2,4)
p = df.plot(eps_12, title="$eps_12$", cmap="inferno")
plt.colorbar(p) 

Inverse Problem - Noiseless data - Figure 3b

mu_rec, err, prim = noiseless_data_rec(
    u,
    mesh_params=mesh_params,
    mat_params=mat_params,
    load_params=load_params)

plt.subplot(2,2,3)
p = df.plot(mu_rec, title="mu_rec", cmap="inferno", vmin=1.0, vmax=3.0)
plt.colorbar(p) 

Inverse Problem - noiseless data - Figure 4

Nx_coarse_list = [200, 100, 50, 40, 25, 20]
l = len(Nx_coarse_list)

Error_list = np.zeros(l)
h_list = np.zeros(l)

for s in range(l):
    Nx_coarse = Nx_coarse_list[s]
    k_0 = Nx/Nx_coarse
    mesh_coarse = df.UnitSquareMesh(Nx_coarse, Nx_coarse)
    dx_coarse = df.Measure("dx", domain=mesh_coarse)
    degree_coarse = degree
    V1_coarse = df.VectorFunctionSpace(mesh_coarse, 'CG', degree_coarse)
    V_S_coarse = df.FunctionSpace(mesh_coarse, 'CG', degree_coarse-1)  
    V2_coarse = df.FunctionSpace(mesh_coarse, 'CG', degree_coarse-1)  

    u_coarse = df.project(u, V1_coarse)
    S_coarse = (df.nabla_grad(u_coarse) + df.nabla_grad(u_coarse).T)/2
    S_11_coarse = S_coarse[0,0]
    S_22_coarse = S_coarse[1,1]
    S_12_coarse = S_coarse[0,1]

    V_rot_coarse = df.VectorFunctionSpace(mesh_coarse, 'CG', degree_coarse-1)
    V_S1_coarse = df.FunctionSpace(mesh_coarse, 'CG', degree_coarse-1) 
    mu1 = df.TrialFunctions(V_rot_coarse)
    nu1 = df.TestFunctions(V_rot_coarse)

    # Normal equation
    a1 = df.dot(df.grad(mu1[0]*(S_11_coarse-S_22_coarse))[0] + df.grad(mu1[0]*S_12_coarse)[1] - df.grad(mu1[1]*S_12_coarse)[0], df.grad(nu1[0]*(S_11_coarse-S_22_coarse))[0] + df.grad(nu1[0]*S_12_coarse)[1] - df.grad(nu1[1]*S_12_coarse)[0]) * dx_coarse + df.dot(df.grad(mu1[1]*S_12_coarse)[1] - df.grad(mu1[0]*S_12_coarse)[0], df.grad(nu1[1]*S_12_coarse)[1] - df.grad(nu1[0]*S_12_coarse)[0]) * dx_coarse
    l1 = df.dot(prim/2, df.grad(nu1[0]*(S_11_coarse-S_22_coarse))[0] + df.grad(nu1[0]*S_12_coarse)[1] - df.grad(nu1[1]*S_12_coarse)[0]) * dx_coarse

    def custom_boundary1(x, on_boundary):
        return on_boundary and (x[0]<tol or x[0]>1-tol or x[1]<tol)
    def custom_boundary2(x, on_boundary):
        return on_boundary  and x[1]<tol 
    bc_mu1 = df.DirichletBC(V_rot_coarse.sub(0), df.Constant(0), custom_boundary1) 
    bc_mu2 = df.DirichletBC(V_rot_coarse.sub(1), df.Constant(0), custom_boundary2)
    bc = [bc_mu1, bc_mu2]

    mu1 = df.Function(V_rot_coarse)
    df.solve(a1 == l1, mu1, bc, solver_parameters={'linear_solver':'mumps'})
    mu_function_rec_coarse = df.project(mu1[0]+1, V_S1_coarse)  
    mu_function_coarse = df.project(mu_function, V_S1_coarse)
    norm1_l2 = df.sqrt(df.assemble((mu_function_coarse - mu_function_rec_coarse)**2 * dx_coarse)) / df.sqrt(df.assemble(mu_function_coarse**2 * dx_coarse))
    h_list[s] = 1/Nx_coarse
    Error_list[s] = norm1_l2
    print('For noiseless data, the L2-norm for the reconstruction:', norm1_l2)
    
print("h_list =", h_list)    
print("Error_list =", Error_list) 

# Calcul convergence entre points consécutifs
slope = np.log(Error_list[:-1] / Error_list[1:]) / np.log(h_list[:-1] / h_list[1:])
print("Convergence orders between points:", slope)

# Régression linéaire en log-log
log_h = np.log(h_list)
log_err = np.log(Error_list)

# Ajustement linéaire : log(Error) ≈ a*log(h) + b
coeffs = np.polyfit(log_h, log_err, 1)
a, b = coeffs
print(f"Régression linéaire: slope = {a:.2f}, intercept = {b:.2f}")

# Erreur prédite par la régression
fit_err = np.exp(a*log_h + b)

# Plot
fig, ax = plt.subplots()

# Courbe des erreurs
ax.plot(h_list, Error_list, "-sr", label="Relative error")

# Régression linéaire (log-log)
ax.plot(h_list, fit_err, "--g", label=f"Fit slope {a:.2f}")

# Références pente 1 et 2
ref_error_2 = 1e3 * h_list**2
ax.plot(h_list, ref_error_2, "--k", label="Slope 2 ($h^2$)")

ref_error_1 = 1 * h_list
ax.plot(h_list, ref_error_1, "--b", label="Slope 1 ($h$)")

# Échelles log
ax.set_xscale("log")
ax.set_yscale("log")

# Mise en forme
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.set_xlabel('Mesh size (h)')
ax.set_ylabel('Error')
ax.legend()
plt.show()

Inverse Problem - noisy data

ureal = u

Noise_list = np.array([1e-7, 1e-6])
Nx_coarse_list = [50] #[200, 100, 50, 40, 25, 20, 10]
l = len(Nx_coarse_list)
tol = 1e-15

k_list = []
size_noise = np.size(Noise_list)
Error_list_with_projection = np.zeros((l, size_noise))
Error_list_without_projection = np.zeros(size_noise)
h = ureal.vector().get_local()
a = h.shape 
V = u.function_space()
mesh = u.function_space().mesh()
dx = df.Measure("dx", domain=mesh)

Inverse Problem - noisy data - Figure 3

for k in range(size_noise) :
    print("Noise = ", Noise_list[k])
    noise_Expression_2 = Function(V)

    for i in range(a[0]) :
        if k==1 :
           noise_Expression_2.vector()[i] = np.random.normal(0,Noise_list[k])
        else :
            noise_Expression_2.vector()[i] = np.random.normal(Noise_list[k-1],Noise_list[k])
    end 
    print( noise_Expression_2.vector())
 
    u_noised = Function(V)
    u_noised.vector()[:] = ureal.vector()[:] + noise_Expression_2.vector()[:] 


    #Without Projection
    mu_rec, err, prim= noiseless_data_rec( u_noised, 
                        mesh_params= mesh_params,
                        mat_params=mat_params,
                        load_params=load_params)
    Error_list_without_projection[k]= err
 
    plt.figure(figsize=(5,4))
    p = plot(mu_rec, title=f"mu_rec (no reg) — noise={Noise_list[k]}", cmap="inferno", vmin=1.0, vmax=3.0)
    plt.colorbar(p)
    plt.tight_layout()
    plt.show(block=False)
    plt.pause(0.5) 


    #With regularization
    for s in range(l) :
        Nx_coarse = Nx_coarse_list[s]
        k_0=Nx/Nx_coarse
        k_list.append(k_0)
        print("k_0= h_0/h_1=", k_0)
        mesh_coarse = UnitSquareMesh(Nx_coarse, Nx_coarse)
        dx_coarse = Measure("dx", domain=mesh_coarse)
        degree_coarse = degree
        V1_coarse = VectorFunctionSpace(mesh_coarse, 'DG', degree_coarse)
        V_tensor_coarse = TensorFunctionSpace(mesh_coarse, 'DG', degree=degree_coarse-1, shape=(dim, dim))
        
        u=TrialFunction(V1_coarse)
        v=TestFunction(V1_coarse)
        a_proj = dot(u,v)*dx_coarse 
        l_proj = dot(u_noised,v)*dx_coarse 
        u_proj= Function(V1_coarse)
        solve(a_proj== l_proj, u_proj,  solver_parameters={'linear_solver': 'mumps'})
        
        V1_DG= VectorFunctionSpace(mesh,  'DG',degree_coarse)
        epsilon_udelta =  0.5*( nabla_grad(interpolate(u_proj,V1_DG)) + nabla_grad(interpolate(u_proj,V1_DG) ).T)
        
        S_11delta = epsilon_udelta [0, 0]
        S_22delta = epsilon_udelta [1, 1]
        S_12delta = epsilon_udelta [0, 1]

        V_mu= mu_function.function_space()
        V_S = FunctionSpace(mesh, 'CG',degree)  

        #By taking into account that the components of eps are in Wˆ{1,\infty}
        S_11delta_cont = project(S_11delta,V_mu )
        S_22delta_cont = project(S_22delta, V_mu )
        S_12delta_cont = project(S_12delta, V_mu ) 

        f1=f[0]
        f2=f[1]

        u1= interpolate(u_proj,V1_DG)[0]
        u2= interpolate(u_proj,V1_DG)[1]

        #Compute the antiderivative
        F=- rho * omega**2* (grad(u1 )[1] - grad( u2 )[0] )  - 2*grad(grad(S_11delta - S_22delta)[1])[0] - 2*grad(grad(S_12delta)[1])[1] + 2*grad(grad(S_12delta)[0])[0]

        def boundary(x, on_boundary):
           return on_boundary and x[0]<tol 
        bc = DirichletBC(V_S, Constant(0), boundary)
        w1 = TrialFunction(V_S)
        w1_test = TestFunction(V_S)
        a2 = grad(w1)[0]*w1_test*dx   
        l2 = F* w1_test*dx
        prim = Function(V_S)
        solve(a2 == l2, prim, bc)

        V_rot = VectorFunctionSpace(mesh, 'CG', degree-1)
        V_S1 = FunctionSpace(mesh, 'CG',degree-1) 
        mu1 = TrialFunctions(V_rot)
        nu1 = TestFunctions(V_rot)

        a1=dot(grad(S_12delta_cont*mu1[0])[1] + 0.5*grad((S_11delta-S_22delta)/S_12delta)[0] * mu1[0]*S_12delta_cont + 0.5*(S_11delta_cont-S_22delta_cont)/S_12delta_cont*grad(mu1[0]*S_12delta_cont)[0] - grad(mu1[1]*S_12delta_cont)[0] ,nu1[0] )*dx  +  dot(grad(S_12delta_cont*mu1[1])[1] ,nu1[1] )*dx - 0.5*dot((S_11delta_cont-S_22delta_cont)*mu1[0],grad(nu1[0])[0])*dx  + dot(S_12delta_cont*mu1[0], grad(nu1[1])[0])*dx    
    
        l1=dot (prim /2,nu1[1])* dx

        def custom_boundary1(x, on_boundary):
            return on_boundary and (x[0]<tol or x[0] > 1-tol or x[1]<tol)
        def custom_boundary2(x, on_boundary):
            return on_boundary  and x[1]<tol
        bc_mu1 = DirichletBC(V_rot.sub(0), Constant(0), custom_boundary1)    
        bc_mu2 = DirichletBC(V_rot.sub(1), Constant(0), custom_boundary2)
        bc = [bc_mu1, bc_mu2]


        #Normal equation
        A = assemble(a1)
        b = assemble(l1)
        [bc.apply(A, b) for bc in bc]
    
        mu1 = Function(V_rot )
        [bc.apply(A, b) for bc in bc]

        # Calculer A^T A et A^T b
        from petsc4py import PETSc
        from dolfin import as_backend_type

        A_petsc = as_backend_type(A).mat()
        #AT_petsc = A_petsc.transpose()
        b_petsc = as_backend_type(b).vec()


        AT_petsc = PETSc.Mat().createTranspose(A_petsc)

        # Multiplier AT par b pour obtenir AT_b
        AT_b_petsc = AT_petsc * b_petsc
    
        # Créer la matrice normale ATA
        ATA_petsc = AT_petsc * A_petsc

        # Définir le vecteur solution
        mu1_petsc = b_petsc.duplicate()

        # Configurer le solveur PETSc
        ksp = PETSc.KSP().create()
        ksp.setOperators(ATA_petsc)
        ksp.setFromOptions()
        ksp.setType('cg')
        ksp.getPC().setType('none')

        # Résoudre le système linéaire
        ksp.solve(AT_b_petsc, mu1_petsc)

        converged_reason = ksp.getConvergedReason()
        #print(converged_reason)

        if converged_reason > 0:
            print("The solver has converged successfully.")
        else:
            print("Solver diverged due to high residual, but let's check the result.")
            print(f"Final residual norm: {ksp.getResidualNorm()}")
       
        mu1 = Function(V_rot)
        mu1.vector().set_local(mu1_petsc.getArray())

        mu_function_rec_bruit  = project(mu1[0]+ 1 , V_S1)  

        plt.figure(figsize=(5,4))
        p = plot(mu_function_rec_bruit, title=f"mu_rec (reg) with k={k_0} — noise={Noise_list[k]}", cmap="inferno", vmin=1.0, vmax=3.0)
        plt.colorbar(p)
        plt.tight_layout()
        plt.show(block=False)
        plt.pause(0.5) 

        norm1_l2_normal = sqrt(assemble((mu_function - mu_function_rec_bruit)**2 * dx)) / sqrt(assemble(mu_function**2 * dx))
        print('erreur en norme L2 pour l approche avec bruit avec projection, equation normale:', norm1_l2_normal)
        #####################
        Error_list_with_projection[s,k] = norm1_l2_normal 
 

print("k_list=", k_list)
print("Noise=", Noise_list)

print("Error_with_projection_normal_eqt = ", Error_list_with_projection)
print("Error_without_projection = ", Error_list_without_projection)