adolc_plasticity

autotoc_md337

jupyter: jupytext: formats: ipynb,md text_representation: extension: .md format_name: markdown format_version: '1.3' jupytext_version: 1.16.0 kernelspec: display_name: Python 3 (ipykernel) language: python

name: python3

Importing necessary modules, setting parameters and paths

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.pyplot import cm
from scipy import optimize
import time
import os
import os.path
import zipfile
import pandas as pd
from scipy.optimize import curve_fit, least_squares
import sys
import gmsh
import math
 
import matplotlib.image as mpimg
import re
 
import pyvista as pv
import ipywidgets as widgets
pv.set_plot_theme("document")
 
plt.rcParams['figure.figsize'] = [12, 9]
plt.rcParams['figure.dpi'] = 300
plt.rcParams['font.family'] = "Serif"
plt.rcParams['font.size'] = 15
 
from pyvirtualdisplay import Display
display = Display(backend="xvfb", visible=False, size=(1024, 768))
display.start()
    
user_name=!whoami # get user name
user_name=user_name[0]
um_view = "$HOME/um_view_release"
 
# new executables
exe2d = um_view + "/adolc_plasticity/adolc_plasticity_2d"
exe3d = um_view + "/adolc_plasticity/adolc_plasticity_3d"
mofem_part = um_view + "/bin/mofem_part"

Define utility functions (mesh generation, running the analysis, results post-processing)

class AttrDict(dict):
    def __getattr__(self, attr):
        if attr in self:
            return self[attr]
        raise AttributeError(f"'AttrDict' object has no attribute '{attr}'")
        
        
# looking through the output for norms and other phrases 
def parse_log_file(params):
    res = AttrDict()
    res.force_x = []
    res.time = []
    res.min_x = []
    res.max_x = []
    res.scale = []
    with open(params.log_file, "r") as log_file:
        for line in log_file:
            line = line.strip()
            if "nb global dofs" in line:
                res.dof_num = line.split()[13]
            if "norm_u: " in line:
                res.norm_u = float(line.split()[4])
            if "TS dt" in line:
                res.time.append(float(line.split()[8]))
            if "Block DISPLACEMENTSET1 Force:" in line:
                res.force_x.append(float(line.split()[6]))
            if "Min displacement" in line:
                res.min_x.append(float(line.split()[5]))
            if "Max displacement" in line:
                res.max_x.append(float(line.split()[5]))
            if "scale:" in line:
                res.scale.append(float(line.split()[6]))
 
    #print(res.time)
    #print(res.scale)
    #print(res.force_x)
    
    return res
 
# running MoFEM executable
def run_mofem(exe, params):
    if params.order > 5:
        raise Exception("Approximation order greater than 5 should not be used in this tutorial")
    if params.nproc > 32:
        raise Exception("Not more than 4 processors should be used in this tutorial")
    
    # define name of the mesh file
    params.part_file = params.mesh_file[:-4] + "_" + str(params.nproc) + "p.h5m"
    
    # partition the mesh
    !{mofem_part} \
    -my_file {params.mesh_file} \
    -my_nparts {params.nproc} -output_file {params.part_file} 
    
    # run the chosen MoFEM executable
    !export OMPI_MCA_btl_vader_single_copy_mechanism=none && \
    nice -n 10 mpirun --oversubscribe --allow-run-as-root \
    -np {params.nproc} {exe} \
    -snes_linesearch_type {params.line_search} \
    -file_name {params.part_file} \
    -order {params.order} \
    -b_bar {params.bbar} \
    -young_modulus {params.young_modulus} \
    -poisson_ratio {params.poisson_ratio} \
    -sigma_y {params.sigma_y} \
    -H {params.H} \
    -K {params.K} \
    -K {params.phi} \
    -Ht {params.Ht} \
    -Hc {params.Hc} \
    -sigma_yt {params.sigma_yt} \
    -sigma_yc {params.sigma_yc} \
    -nup {params.nup} \
    -material {params.material} \
    -ts_dt {params.dt} \
    -ts_max_time {params.final_time} \
    -ts_max_steps {params.max_steps} \
    -ts_mofem_adapt_desired_it {params.desired} \
    -ts_mofem_adapt_off {params.adapt_off} \
    -time_scalar_file {params.scaling_file} \
    -reaction_print_block_name \
    -log_no_color \
    2>&1 | tee {params.log_file} 
    
    # convert results to vtk if desired
    if params.convert_result:
        !convert.py -np {params.nproc} out*
    
    # return results found by parse_log_file
    return parse_log_file(params)

# showing results of vtks
def show_results(params):
    out_to_vtk = !ls -c1 {params.show_file}*vtk
    last_file=out_to_vtk[0]
 
    mesh = pv.read(last_file[:-3] + "vtk")
    if params.warp_field:
        mesh = mesh.warp_by_vector(vectors=params.warp_field, factor=params.warp_factor)
 
    if params.show_edges:
        mesh=mesh.shrink(0.95)
    
    jupyter_backend='ipygany'
    cmap = "turbo"
 
    p = pv.Plotter(notebook=True)
    p.add_mesh(mesh, scalars=params.show_field, component=None, smooth_shading=True, cmap=cmap)
    p.camera_position = "xy"
    
    p.enable_parallel_projection()
    p.enable_image_style()
    
    p.show(jupyter_backend=jupyter_backend)
    
    
def plot_data(params_list):
    # Plot data
    colors_list = ['r', 'b', 'g', 'p']
    fig, ax = plt.subplots(2)
    for params,c in zip(params_list, colors_list):
        ax[0].plot(params.time, -np.array(params.force_x), 
                   c+'o-', label='force vs time')
    ax[0].axis(xmin=0)
    ax[0].legend(loc='best')
    ax[0].grid(True)
    for params,c in zip(params_list, colors_list):
        ax[1].plot(params.scale, -np.array(params.force_x), 
                   c+'o-', label='force vs displacement')
    #ax[1].axis(xmin=0)
    ax[1].legend(loc='best')
    ax[1].grid(True)
    
import os
import zipfile    
    
def zip_vtks(file_name):
    zipf = zipfile.ZipFile(file_name, 'w', zipfile.ZIP_DEFLATED)
    path = "."
    for root, dirs, files in os.walk(path):
        for file in files:
            if file.endswith('.vtk'):        
                zipf.write(os.path.join(root, file), 
                           os.path.relpath(os.path.join(root, file), os.path.join(path, '..')))
    zipf.close()

params = AttrDict() # Attribute dictionary for storing the parameters
 
# material parameters
params.young_modulus = 20
params.poisson_ratio = 0.25
params.sigma_y = 1 
params.H = 1e-2
params.K = 0
params.phi = 0
params.Ht = 1e-2
params.Hc = 1e-2
params.sigma_yt = 1
params.sigma_yc = 2
params.nup = 0.
 
# solution parameters
params.log_file = "log" # log file name 
params.nproc = 32
params.bbar = 1 # only meake sens when orde = 1
params.order = 1 # approximation order for displacements
 
# post-processing defaults
params.warp_factor = 1
 
# general params
params.desired = 5
params.adapt_off = 1
params.max_steps = 125
params.final_time = 1.25
params.dt = 0.01
params.line_search = "cp"
 
# scaling file 
params.scaling_file = "hist.txt"

 magic_args="-s \"$params.scaling_file\""
cat << EOF > $1
0 0
0.25 0.25
0.5 0.0
0.75 -0.25
1.0 0.
1.25 0.25 
EOF

!rm -rf out*                      # remove previous analysis output files
params.convert_result = True      # convert the results to vtk for visualisation
params.mesh_file = "plate_2d.cub"
params.material = "vonMisses"
res_misses_isotropic = run_mofem(exe2d, params)   # run analysis
zip_vtks("misses_isotropic_vtk.zip")

plot_data([res_misses_isotropic])

# Post-processing parameters
params.show_file = "out_step_100"
params.show_field = "PLASTIC_STRAIN"
params.warp_field = "U"     # warp with displacement
params.warp_factor = 0.2
params.show_edges = True
params.bbar = 0
show_results(params)

!rm -rf out*                      # remove previous analysis output files
params.convert_result = True      # convert the results to vtk for visualisation
params.mesh_file = "plate_2d.cub"
params.material = "Paraboloidal"
params.nup = 0.0
res_paraboloidal=run_mofem(exe2d, params)   # run analysis
zip_vtks("paraboloidal_vtk.zip")

plot_data([res_paraboloidal])

# Post-processing parameters
params.show_file = "out_step_125"
params.show_field = "PLASTIC_STRAIN"
params.warp_field = "U"     # warp with displacement
params.warp_factor = 0.2
params.show_edges = True
show_results(params)

plot_data([res_misses_isotropic, res_paraboloidal])