2D Case - Figures 12-13-15

Convergence to reference solution for training performed using Adam and L-BFGS.

We investigate the effect of r-adaptivity and rh-adaptivity on the accuracy of displacement and von Mises stress.

#%% Libraries import
# import HiDeNN library
import sys  
# sys.path.append("../neurom/")
from neurom.HiDeNN_PDE import MeshNN, NeuROM, MeshNN_2D, MeshNN_1D
# Import pre-processing functions
import neurom.src.Pre_processing as pre
# Import torch librairies
import torch
import torch.nn as nn
# Import mechanical functions
from neurom.src.PDE_Library import Strain, Stress,VonMises_plain_strain
# Import Training funcitons
from neurom.src.Training import Training_2D_FEM
#Import post processing libraries
import neurom.Post.Plots as Pplot
import time
import os
import torch._dynamo as dynamo
mps_device = torch.device("mps")
from importlib import reload  # Python 3.4+
import tomllib
import numpy as numpy
import argparse

Default_config_file = 'Configurations/config_2D.toml'

with open(Default_config_file, mode="rb") as f:
    config = tomllib.load(f)

# Prerequisites
mesh_resolution = [44, 144, 484, 1804, 7088 ]

eval_coord_file = "GroundTruth/eval_coordinates.npy"
num_displ_file = "GroundTruth/num_displacement.npy"
num_VM_stress_file = "GroundTruth/num_VM_stress.npy"
eval_coord =  torch.tensor(numpy.load(eval_coord_file), dtype=torch.float64, requires_grad=True)
num_displ = torch.tensor(numpy.load(num_displ_file))
num_VM_stress = torch.tensor(numpy.load(num_VM_stress_file))

# Experiment
element_size = [2.0, 1.0, 0.5, 0.25, 0.125]
config["solver"]["FrozenMesh"] = True
optimizers = ["adam","lbfgs"]

error_u = numpy.zeros((len(element_size),len(optimizers)))
error_v = numpy.zeros((len(element_size),len(optimizers)))
error_stress = numpy.zeros((len(element_size),len(optimizers)))
error_stress_max = numpy.zeros((len(element_size),len(optimizers)))

for e in range(len(element_size)):
    config["interpolation"]["MaxElemSize2D"] = element_size[e]
    for op in range(len(optimizers)):
        config["training"]["optimizer"] = optimizers[op]

        if config["interpolation"]["dimension"] == 1:
            Mat = pre.Material(             flag_lame = True,                               # If True should input lmbda and mu instead of E and nu
                                            coef1     = config["material"]["E"],            # Young Modulus
                                            coef2     = config["geometry"]["A"]             # Section area of the 1D bar
                                )
        elif config["interpolation"]["dimension"] == 2:
            try:
                Mat = pre.Material(         flag_lame = False,                              # If True should input lmbda and mu instead of E and nu
                                            coef1     = config["material"]["E"],            # Young Modulus
                                            coef2     = config["material"]["nu"]            # Poisson's ratio
                                )
            except:
                Mat = pre.Material(         flag_lame = True,                               # If True should input lmbda and mu instead of E and nu
                                            coef1     = config["material"]["lmbda"],        # First Lame's coef
                                            coef2     = config["material"]["mu"]            # Second Lame's coef
                                )

        MaxElemSize = pre.ElementSize(
                                        dimension     = config["interpolation"]["dimension"],
                                        L             = config["geometry"]["L"],
                                        order         = config["interpolation"]["order"],
                                        np            = config["interpolation"]["np"],
                                        MaxElemSize2D = config["interpolation"]["MaxElemSize2D"]
                                    )
        Excluded = []
        Mesh_object = pre.Mesh( 
                                        config["geometry"]["Name"],                 # Create the mesh object
                                        MaxElemSize, 
                                        config["interpolation"]["order"], 
                                        config["interpolation"]["dimension"]
                                )

        Mesh_object.AddBorders(config["Borders"]["Borders"])
        Mesh_object.AddBCs(                                                         # Include Boundary physical domains infos (BCs+volume)
                                        config["geometry"]["Volume_element"],
                                        Excluded,
                                        config["DirichletDictionryList"]
                            )                   

        Mesh_object.MeshGeo()                                                       # Mesh the .geo file if .msh does not exist
        Mesh_object.ReadMesh()                                                      # Parse the .msh file
        Mesh_object.ExportMeshVtk()

        if int(Mesh_object.dim) != int(Mesh_object.dimension):
            raise ValueError("The dimension of the provided geometry does not match the job dimension")

        Model_FEM = MeshNN_2D(Mesh_object, n_components = 2)

        Model_FEM.Freeze_Mesh()
        Model_FEM.UnFreeze_FEM()

        if not config["solver"]["FrozenMesh"]:
            Model_FEM.UnFreeze_Mesh() 

        Model_FEM.RefinementParameters( MaxGeneration = config["training"]["h_adapt_MaxGeneration"], 
                                    Jacobian_threshold = config["training"]["h_adapt_J_thrshld"])

        Model_FEM.TrainingParameters(   loss_decrease_c = config["training"]["loss_decrease_c"], 
                                        Max_epochs = config["training"]["n_epochs"], 
                                        learning_rate = config["training"]["learning_rate"])
        Model_FEM = Training_2D_FEM(Model_FEM, config, Mat)
        
        # evaluation 
        Model_FEM.eval()

        elem_IDs = torch.tensor(Model_FEM.mesh.GetCellIds(eval_coord),dtype=torch.int)
        u = Model_FEM(eval_coord, elem_IDs)
        eps =  Strain(u,eval_coord)
        sigma =  torch.stack(Stress(eps[:,0], eps[:,1], eps[:,2], Mat.lmbda, Mat.mu),dim=1)
        sigma_VM = VonMises_plain_strain(sigma, Mat.lmbda, Mat.mu)

        error_u[e,op] = (torch.linalg.vector_norm(num_displ[:,0] - u[0,:])/torch.linalg.vector_norm(num_displ[:,0])).item()
        error_v[e,op] = (torch.linalg.vector_norm(num_displ[:,1] - u[1,:])/torch.linalg.vector_norm(num_displ[:,1])).item()
        error_stress[e,op] = (torch.linalg.vector_norm(num_VM_stress - sigma_VM)/torch.linalg.vector_norm(num_VM_stress)).item()
        error_stress_max[e,op] = (max(sigma_VM)/max(num_VM_stress)).item()

print("Adam")
print("u = ", error_u[:,0])
print("v = ", error_v[:,0])
print("s = ", error_stress[:,0])
print("s max = ", error_stress_max[:,0])
print()
print("LBFGS")
print("u = ", error_u[:,1])
print("v = ", error_v[:,1])
print("s = ", error_stress[:,1])
print("s max = ", error_stress_max[:,1])

import matplotlib.pyplot as plt
import matplotlib
plt.rcParams['text.usetex'] = False

# Plot normalized displacement error

fig = matplotlib.pyplot.gcf()
ax = plt.gca()

plt.plot(mesh_resolution, error_u[:,0],'--', color = "darkblue", label = r'$\| e_{u_x}\|_2$'+ ", Adam")
plt.plot(mesh_resolution, error_v[:,0],'--', color = "purple", label = r'$\| e_{u_y}\|_2$'+ ", Adam")
plt.plot(mesh_resolution, error_stress[:,0],'--', color = "green", label = r'$\| e_{\sigma_{VM}}\|_2$'+ ", Adam")

plt.plot(mesh_resolution, error_u[:,1],'-', color = "darkblue", label = r'$\| e_{u_x}\|_2$'+ ", L-BFGS")
plt.plot(mesh_resolution, error_v[:,1],'-', color = "purple", label = r'$\| e_{u_y}\|_2$'+ ", L-BFGS")
plt.plot(mesh_resolution, error_stress[:,1],'-', color = "green", label = r'$\| e_{\sigma_{VM}}\|_2$'+ ", L-BFGS")


ax.set_xscale('log')
ax.set_yscale('log')

plt.xlabel("Number of mesh nodes")
plt.ylabel("Normalized displacement error")
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5),frameon=False )

plt.show()

# Prerequisites
mesh_resolution = [44, 144, 484, 1804, 7088 ]

# eval_coord_file = "../2D_example/eval_coordinates.npy"
# num_displ_file = "../2D_example/num_solution/num_displacement.npy"
# num_VM_stress_file = "../2D_example/num_solution/num_VM_stress.npy"
eval_coord_file = "GroundTruth/eval_coordinates.npy"
num_displ_file = "GroundTruth/num_displacement.npy"
num_VM_stress_file = "GroundTruth/num_VM_stress.npy"

eval_coord =  torch.tensor(numpy.load(eval_coord_file), dtype=torch.float64, requires_grad=True)
num_displ = torch.tensor(numpy.load(num_displ_file))
num_VM_stress = torch.tensor(numpy.load(num_VM_stress_file))

# Experiment
element_size = 2.0
config["solver"]["FrozenMesh"] = False
optimizers = ["adam","lbfgs"]
refinment = [1,2,3,4,5]

config["interpolation"]["MaxElemSize2D"] = element_size

r_adapt_error_u = numpy.zeros((len(refinment),len(optimizers)))
r_adapt_error_v = numpy.zeros((len(refinment),len(optimizers)))
r_adapt_error_stress = numpy.zeros((len(refinment),len(optimizers)))
r_adapt_error_stress_max = numpy.zeros((len(refinment),len(optimizers)))

for e in range(len(refinment)):
    config["training"]["multiscl_max_refinment"] = refinment[e]
    for op in range(len(optimizers)):
        config["training"]["optimizer"] = optimizers[op]

        Mat = pre.Material(             flag_lame = True,                          # If True should input lmbda and mu instead of E and nu
                                        coef1     = config["material"]["lmbda"],        # Young Modulus
                                        coef2     = config["material"]["mu"]        # Poisson's ratio
                            )

        MaxElemSize = pre.ElementSize(
                                        dimension     = config["interpolation"]["dimension"],
                                        L             = config["geometry"]["L"],
                                        order         = config["interpolation"]["order"],
                                        np            = config["interpolation"]["np"],
                                        MaxElemSize2D = config["interpolation"]["MaxElemSize2D"]
                                    )
        Excluded = []
        Mesh_object = pre.Mesh( 
                                        config["geometry"]["Name"],                 # Create the mesh object
                                        MaxElemSize, 
                                        config["interpolation"]["order"], 
                                        config["interpolation"]["dimension"]
                                )

        Mesh_object.AddBorders(config["Borders"]["Borders"])
        Mesh_object.AddBCs(                                                         # Include Boundary physical domains infos (BCs+volume)
                                        config["geometry"]["Volume_element"],
                                        Excluded,
                                        config["DirichletDictionryList"]
                            )                   

        Mesh_object.MeshGeo()                                                       # Mesh the .geo file if .msh does not exist
        Mesh_object.ReadMesh()                                                      # Parse the .msh file
        Mesh_object.ExportMeshVtk()

        if int(Mesh_object.dim) != int(Mesh_object.dimension):
            raise ValueError("The dimension of the provided geometry does not match the job dimension")

        Model_FEM = MeshNN_2D(Mesh_object, n_components = 2)

        Model_FEM.Freeze_Mesh()
        Model_FEM.UnFreeze_FEM()

        if not config["solver"]["FrozenMesh"]:
            Model_FEM.UnFreeze_Mesh() 

        Model_FEM.RefinementParameters( MaxGeneration = config["training"]["h_adapt_MaxGeneration"], 
                                    Jacobian_threshold = config["training"]["h_adapt_J_thrshld"])

        Model_FEM.TrainingParameters(   loss_decrease_c = config["training"]["loss_decrease_c"], 
                                        Max_epochs = config["training"]["n_epochs"], 
                                        learning_rate = config["training"]["learning_rate"])
        Model_FEM = Training_2D_FEM(Model_FEM, config, Mat)
        
        # evaluation 
        Model_FEM.eval()

        # Model_FEM.mesh.Nodes = [[i+1,Model_FEM.coordinates[i][0][0].item(),Model_FEM.coordinates[i][0][1].item(),0] for i in range(len(Model_FEM.coordinates))]
        coordinates_all = torch.ones_like(Model_FEM.coordinates_all)
        coordinates_all[Model_FEM.coord_free] = Model_FEM.coordinates['free']
        coordinates_all[~Model_FEM.coord_free] = Model_FEM.coordinates['imposed']
        Nodes = torch.hstack([torch.linspace(1,coordinates_all.shape[0],coordinates_all.shape[0], dtype = coordinates_all.dtype, device = coordinates_all.device)[:,None],
                                coordinates_all])
        Nodes = torch.hstack([Nodes,torch.zeros(Nodes.shape[0],1, dtype = Nodes.dtype, device = Nodes.device)])
        Model_FEM.mesh.Nodes = Nodes.detach().cpu().numpy()
        Model_FEM.mesh.Connectivity = Model_FEM.connectivity
        Model_FEM.mesh.ExportMeshVtk(flag_update = True)

        elem_IDs = torch.tensor(Model_FEM.mesh.GetCellIds(eval_coord),dtype=torch.int)
        u = Model_FEM(eval_coord, elem_IDs)
        eps =  Strain(u,eval_coord)
        sigma =  torch.stack(Stress(eps[:,0], eps[:,1], eps[:,2], Mat.lmbda, Mat.mu),dim=1)
        sigma_VM = VonMises_plain_strain(sigma, Mat.lmbda, Mat.mu)

        r_adapt_error_u[e,op] = (torch.linalg.vector_norm(num_displ[:,0] - u[0,:])/torch.linalg.vector_norm(num_displ[:,0])).item()
        r_adapt_error_v[e,op] = (torch.linalg.vector_norm(num_displ[:,1] - u[1,:])/torch.linalg.vector_norm(num_displ[:,1])).item()
        r_adapt_error_stress[e,op] = (torch.linalg.vector_norm(num_VM_stress - sigma_VM)/torch.linalg.vector_norm(num_VM_stress)).item()
        r_adapt_error_stress_max[e,op] = (max(sigma_VM)/max(num_VM_stress)).item()

print("Adam")
print("u = ", r_adapt_error_u[:,0])
print("v = ", r_adapt_error_v[:,0])
print("s = ", r_adapt_error_stress[:,0])
print("s max = ", r_adapt_error_stress_max[:,0])
print()
print("LBFGS")
print("u = ", r_adapt_error_u[:,1])
print("v = ", r_adapt_error_v[:,1])
print("s = ", r_adapt_error_stress[:,1])
print("s max = ", r_adapt_error_stress_max[:,1])

# Plot normalized displacement error

fig = matplotlib.pyplot.gcf()
ax = plt.gca()

plt.plot(mesh_resolution, error_u[:,1],'-', color = "darkblue", label = r'$\| e_{u_x}\|_2$'+ ", L-BFGS")
plt.plot(mesh_resolution, error_v[:,1],'-', color = "purple", label = r'$\| e_{u_y}\|_2$'+ ", L-BFGS")
plt.plot(mesh_resolution, error_stress[:,1],'-', color = "green", label = r'$\| e_{\sigma_{VM}}\|_2$'+ ", L-BFGS")

plt.plot(mesh_resolution, r_adapt_error_u[:,1],'--', color = "darkblue", label = r'$\| e_{u_x}\|_2$'+ ", L-BFGS, r-adapt.")
plt.plot(mesh_resolution, r_adapt_error_v[:,1],'--', color = "purple", label = r'$\| e_{u_y}\|_2$'+ ", L-BFGS, r-adapt.")
plt.plot(mesh_resolution, r_adapt_error_stress[:,1],'--', color = "green", label = r'$\| e_{\sigma_{VM}}\|_2$'+ ", L-BFGS, r-adapt.")

plt.plot(mesh_resolution, r_adapt_error_u[:,0],':', color = "darkblue", label = r'$\| e_{u_x}\|_2$'+ ", Adam, r-adapt.")
plt.plot(mesh_resolution, r_adapt_error_v[:,0],':', color = "purple", label = r'$\| e_{u_y}\|_2$'+ ", Adam, r-adapt.")
plt.plot(mesh_resolution, r_adapt_error_stress[:,0],':', color = "green", label = r'$\| e_{\sigma_{VM}}\|_2$'+ ", Adam, r-adapt.")

ax.set_xscale('log')
ax.set_yscale('log')
ax.set_ylim([0.00008, 0.3])

plt.xlabel("Number of mesh nodes")
plt.ylabel("Normalized displacement error")
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5),frameon=False )

plt.show()

# Plot maximal stress

fig = matplotlib.pyplot.gcf()
ax = plt.gca()

plt.plot(mesh_resolution, error_stress_max[:,0],'-', color = "purple", label = "Adam")
plt.plot(mesh_resolution, r_adapt_error_stress_max[:,0],':', color = "purple", label = "Adam, r-adaptivity")
plt.plot(mesh_resolution, error_stress_max[:,1],'-', color = "teal", label = "L-BFGS")
plt.plot(mesh_resolution, r_adapt_error_stress_max[:,1],':', color = "teal", label =  "L-BFGS, r-adaptivity")

ax.set_xscale('log')
ax.set_yscale('log')

plt.xlabel("Number of mesh nodes")
plt.ylabel(r'$\sigma^{max}_{VM}$')
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5),frameon=False )

plt.show()

# Prerequisites
mesh_resolution = [44, 144, 484, 1804, 7088 ]

# eval_coord_file = "../2D_example/eval_coordinates.npy"
# num_displ_file = "../2D_example/num_solution/num_displacement.npy"
# num_VM_stress_file = "../2D_example/num_solution/num_VM_stress.npy"

eval_coord_file = "GroundTruth/eval_coordinates.npy"
num_displ_file = "GroundTruth/num_displacement.npy"
num_VM_stress_file = "GroundTruth/num_VM_stress.npy"

eval_coord =  torch.tensor(numpy.load(eval_coord_file), dtype=torch.float64, requires_grad=True)
num_displ = torch.tensor(numpy.load(num_displ_file))
num_VM_stress = torch.tensor(numpy.load(num_VM_stress_file))

# Experiment
element_size = 2.0
config["solver"]["FrozenMesh"] = False
optimizer = "lbfgs"
refinment = [1,2,3,4,5]

config["interpolation"]["MaxElemSize2D"] = element_size
config["training"]["optimizer"] = optimizer
config["training"]["h_adapt_MaxGeneration"] = 2

rh_adapt_error_u = numpy.zeros((len(refinment)))
rh_adapt_error_v = numpy.zeros((len(refinment)))
rh_adapt_error_stress = numpy.zeros((len(refinment)))
rh_adapt_error_stress_max = numpy.zeros((len(refinment)))

for e in range(len(refinment)):
    config["training"]["multiscl_max_refinment"] = refinment[e]

    Mat = pre.Material(             flag_lame = True,                          # If True should input lmbda and mu instead of E and nu
                                    coef1     = config["material"]["lmbda"],        # Young Modulus
                                    coef2     = config["material"]["mu"]        # Poisson's ratio
                        )

    MaxElemSize = pre.ElementSize(
                                    dimension     = config["interpolation"]["dimension"],
                                    L             = config["geometry"]["L"],
                                    order         = config["interpolation"]["order"],
                                    np            = config["interpolation"]["np"],
                                    MaxElemSize2D = config["interpolation"]["MaxElemSize2D"]
                                )
    Excluded = []
    Mesh_object = pre.Mesh( 
                                    config["geometry"]["Name"],                 # Create the mesh object
                                    MaxElemSize, 
                                    config["interpolation"]["order"], 
                                    config["interpolation"]["dimension"]
                            )

    Mesh_object.AddBorders(config["Borders"]["Borders"])
    Mesh_object.AddBCs(                                                         # Include Boundary physical domains infos (BCs+volume)
                                    config["geometry"]["Volume_element"],
                                    Excluded,
                                    config["DirichletDictionryList"]
                        )                   

    Mesh_object.MeshGeo()                                                       # Mesh the .geo file if .msh does not exist
    Mesh_object.ReadMesh()                                                      # Parse the .msh file
    Mesh_object.ExportMeshVtk()

    if int(Mesh_object.dim) != int(Mesh_object.dimension):
        raise ValueError("The dimension of the provided geometry does not match the job dimension")

    Model_FEM = MeshNN_2D(Mesh_object, n_components = 2)

    Model_FEM.Freeze_Mesh()
    Model_FEM.UnFreeze_FEM()

    if not config["solver"]["FrozenMesh"]:
        Model_FEM.UnFreeze_Mesh() 

    Model_FEM.RefinementParameters( MaxGeneration = config["training"]["h_adapt_MaxGeneration"], 
                                Jacobian_threshold = config["training"]["h_adapt_J_thrshld"])

    Model_FEM.TrainingParameters(   loss_decrease_c = config["training"]["loss_decrease_c"], 
                                    Max_epochs = config["training"]["n_epochs"], 
                                    learning_rate = config["training"]["learning_rate"])
    Model_FEM = Training_2D_FEM(Model_FEM, config, Mat)
    
    # evaluation 
    Model_FEM.eval()

    # Model_FEM.mesh.Nodes = [[i+1,Model_FEM.coordinates[i][0][0].item(),Model_FEM.coordinates[i][0][1].item(),0] for i in range(len(Model_FEM.coordinates))]
    coordinates_all = torch.ones_like(Model_FEM.coordinates_all)
    coordinates_all[Model_FEM.coord_free] = Model_FEM.coordinates['free']
    coordinates_all[~Model_FEM.coord_free] = Model_FEM.coordinates['imposed']
    Nodes = torch.hstack([torch.linspace(1,coordinates_all.shape[0],coordinates_all.shape[0], dtype = coordinates_all.dtype, device = coordinates_all.device)[:,None],
                            coordinates_all])
    Nodes = torch.hstack([Nodes,torch.zeros(Nodes.shape[0],1, dtype = Nodes.dtype, device = Nodes.device)])
    Model_FEM.mesh.Nodes = Nodes.detach().cpu().numpy()
    Model_FEM.mesh.Connectivity = Model_FEM.connectivity
    Model_FEM.mesh.ExportMeshVtk(flag_update = True)

    elem_IDs = torch.tensor(Model_FEM.mesh.GetCellIds(eval_coord),dtype=torch.int)
    u = Model_FEM(eval_coord, elem_IDs)
    eps =  Strain(u,eval_coord)
    sigma =  torch.stack(Stress(eps[:,0], eps[:,1], eps[:,2], Mat.lmbda, Mat.mu),dim=1)
    sigma_VM = VonMises_plain_strain(sigma, Mat.lmbda, Mat.mu)

    rh_adapt_error_u[e] = (torch.linalg.vector_norm(num_displ[:,0] - u[0,:])/torch.linalg.vector_norm(num_displ[:,0])).item()
    rh_adapt_error_v[e] = (torch.linalg.vector_norm(num_displ[:,1] - u[1,:])/torch.linalg.vector_norm(num_displ[:,1])).item()
    rh_adapt_error_stress[e] = (torch.linalg.vector_norm(num_VM_stress - sigma_VM)/torch.linalg.vector_norm(num_VM_stress)).item()
    rh_adapt_error_stress_max[e] = (max(sigma_VM)/max(num_VM_stress)).item()

print("LBFGS")
print("u = ", rh_adapt_error_u[:])
print("v = ", rh_adapt_error_v[:])
print("s = ", rh_adapt_error_stress[:])
print("s max = ", rh_adapt_error_stress_max[:])

# Plot normalized displacement error

fig = matplotlib.pyplot.gcf()
ax = plt.gca()

print(rh_adapt_error_u.shape)

plt.plot(mesh_resolution, r_adapt_error_u[:,1],'--', color = "darkblue", label = r'$\| e_{u_x}\|_2$'+ ", r-adaptivity")
plt.plot(mesh_resolution, r_adapt_error_v[:,1],'--', color = "purple", label = r'$\| e_{u_y}\|_2$'+ ", r-adaptivity")
plt.plot(mesh_resolution, r_adapt_error_stress[:,1],'--', color = "green", label = r'$\| e_{\sigma_{VM}}\|_2$'+ ", r-adaptivity")

plt.plot(mesh_resolution, rh_adapt_error_u,':', color = "darkblue", label = r'$\| e_{u_x}\|_2$'+ ", rh-adaptivity")
plt.plot(mesh_resolution, rh_adapt_error_v,':', color = "purple", label = r'$\| e_{u_y}\|_2$'+ ", rh-adaptivity")
plt.plot(mesh_resolution, rh_adapt_error_stress,':', color = "green", label = r'$\| e_{\sigma_{VM}}\|_2$'+ ", rh-adaptivity")


ax.set_xscale('log')
ax.set_yscale('log')
ax.set_ylim([0.00006, 0.3])

plt.xlabel("Number of mesh nodes")
plt.ylabel("Normalized displacement error")
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5),frameon=False )

plt.show()

# Plot maximal stress

fig = matplotlib.pyplot.gcf()
ax = plt.gca()

plt.plot(mesh_resolution, error_stress_max[:,1],'-', color = "lightblue", label = "fixed mesh")
plt.plot(mesh_resolution, r_adapt_error_stress_max[:,1],'--', color = "lightblue", label =  "r-adaptivity")
plt.plot(mesh_resolution, rh_adapt_error_stress_max,':', color = "lightblue", label =  "rh-adaptivity")

ax.set_xscale('log')
ax.set_yscale('log')

plt.xlabel("Number of mesh nodes")
# plt.ylabel(r'$\sigma^{max}_{VM}$')
plt.ylabel(r'$\frac{\sigma_{VM}^{max}}{\sigma_{VM, ref}^{max}}$')
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5),frameon=False )

plt.show()