meta_materoal_indentation

autotoc_md58

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

name: python3

Indentation Test

MoFEM eccosystem

MoFEM has complex software eccosystem composed with many libraries. It scalable can be run on laptops with limited resurcers, up to large coumputer clusters.

MoFEM is designed to run complex multi-physics problems. It supports research, it generated impact case for REF2021, and is used by indsutry.

Set problem parameters

Mesh flie geometry can be made in gMesh, Salome or Coreform.

gMesh, Salome are free and open sourve, Coreform is free to use for research.

Note: MoFEM reads Cubit files, not Cube HDF5 files.

from IPython.display import HTML
 
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
 
plt.rcParams['figure.figsize'] = [15, 10]
 
HTML("<style>.container { width:100% !important; }</style>")

# mesh dimension
dimension=2
 
# material parameters 
young_modulus=1
poisson_ratio=0.3
 
# mesh file
mesh_file='anit-chiral2d.cub'
# output
time_depths_factor=1
 
# surface distance function
import importlib
import sdf
importlib.reload(sdf)
print('radius = %f' % sdf.r)

import importlib
import sdf
importlib.reload(sdf)
print('radius = %f' % sdf.r)
 
import pyvista as pv
from pyvirtualdisplay import Display
 
def show_results(params):
    last_file = params.show_file
 
    mesh = pv.read(last_file[:-3] + "vtk")
    mesh.point_data.keys()
    mesh.point_data.remove('GLOBAL_ID')
    mesh.point_data['xyz']=mesh.points
 
    if params.show_edges:
        mesh = mesh.shrink(0.95)
 
    p = pv.Plotter(notebook=True)
 
    cmap = "turbo"
    
    # Add original mesh
    p.add_mesh(
        mesh,
        component=None,
        cmap=cmap,
        smooth_shading=True,
        name="original"
    )
    
    sphere = pv.Sphere(radius=sdf.r, center=(sdf.xc, sdf.yc, sdf.zc), phi_resolution=45, theta_resolution=45)
    p.add_mesh(sphere, opacity=0.50)
 
    p.camera_position = "xy"
    p.enable_parallel_projection()
    p.show()
    
class AttrDict(dict):
    def __getattr__(self, attr):
        if attr in self:
            return self[attr]
        raise AttributeError(f"'AttrDict' object has no attribute '{attr}'")
 
mesh_file_vtk=mesh_file[:-3]+'vtk'
!~/view_mofem_ep/bin/mbconvert {mesh_file} {mesh_file[:-3]}vtk
        
# Post-processing parameters
params = AttrDict() # Attribute dictionary for storing the parameters
params.show_file = mesh_file_vtk        # file name to visualise, intial configuration is always "out_sol_elastic_0" final step depends on simulation
params.show_edges = True
show_results(params)

wd=!pwd
wd=wd[0]
 
path='~/view_mofem_ep/tutorials/adv-1'
 
if dimension == 3:
    code=path+'/contact_3d'
elif dimension == 2:
    code=path+'/contact_2d'
else:
    raise Exception("Only works for 2d or 3d")
    
print(wd)
print(code)
 
# number of processors
nb_proc=12
 
# parrition mesh
!~/view_mofem_ep/bin/mofem_part -file_name {mesh_file} -my_nparts {nb_proc} -dim {dimension} -adj_dim 0 -output_file one.h5m

# time stepping
time=10
time_step=0.1
 
# computational parameters
cn=1e-2/young_modulus
spring_stiffness=0
alpha_damping=1e-5 # should be small value, used for dynamic relaxation
 
log_file='log'
 
!rm out*h5m
!export OMPI_MCA_btl_vader_single_copy_mechanism=none && \
nice -n 5 mpirun --allow-run-as-root -np {nb_proc} {code} -file_name one.h5m \
-ts_type beuler \
-ts_max_time {time} \
-ts_dt {time_step} \
-time_scalar_file scale_displacemets.txt \
-log_no_color \
-snes_linesearch_type l2 \
-cn {cn} \
-spring_stiffness {spring_stiffness} \
-alpha_damping {alpha_damping} \
-young_modulus 1 \
-poisson_ratio {poisson_ratio} \
-max_post_proc_ref_level 4 \
-order 4 2>&1 | tee {log_file}

# Plot convergence
!sed 's/\[0\] <inform> \[petsc\]//g' {log_file} > log_snes
!grep Function log_snes | sed 's/\[//g' | sed 's/,//g' > snes
newton_data=pd.read_fwf('snes', header=None)
newton_data=newton_data.rename(columns={0: "it", 4: "res", 5: "equlibrium", 6: 'constrain'})
newton_data=newton_data.drop([1, 2, 3], axis=1)
print(newton_data)
 
# newton_data
 
plt.plot(newton_data['res'].to_numpy(),'r^-')
plt.title('Neton method convergence')
plt.ylabel('absolute residial')
plt.xlabel('accumulated iterations')
plt.yscale('log')
plt.grid(True)
plt.show()

# converting analysis files to format readable by post processors
 
!rm -f *.vtk
 
import re, os
list_of_files=!ls -c1 out_*.h5m
def extract_numner(s):
    return int(re.findall(r'\d+',s)[0])
 
size=len(list_of_files)
mod=1 #int(size / 10)
 
list_to_process=[]
for f in list_of_files:
    n=extract_numner(f)
    if n == 0 or n == 1:
        list_to_process.append(f)
    elif n % mod == 0:
        list_to_process.append(f)
        
sorted_list_of_files = sorted(list_to_process, key=extract_numner)
 
out_to_vtk = !ls -c1 out_*h5m
last_file=out_to_vtk[0]
print(last_file)
!~/view_mofem_ep/bin/mbconvert {last_file} {last_file[:-3]}vtk
 
import pyvista as pv
from pyvirtualdisplay import Display
 
 
def show_results(params):
 
    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)
 
    cmap = "turbo"
    #cmap = "viridis"
    #cmap = "Spectral"
    p = pv.Plotter(notebook=True)
 
    # Add original mesh
    p.add_mesh(
        mesh,
        scalars=params.show_field,
        preference="point",
        smooth_shading=False,
        cmap=cmap,
        clim=params.clim if params.clim else None,
        name="original"
    )
 
    p.camera_position = "xy"
    p.enable_parallel_projection()
    p.show()
 
class AttrDict(dict):
    def __getattr__(self, attr):
        if attr in self:
            return self[attr]
        raise AttributeError(f"'AttrDict' object has no attribute '{attr}'")
        
# Post-processing parameters
params = AttrDict() # Attribute dictionary for storing the parameters
params.show_file = last_file        # file name to visualise, intial configuration is always "out_sol_elastic_0" final step depends on simulation
# Available fields - STRESS, U, STRAIN
params.show_field = "STRAIN"
params.warp_field = "U" 
params.warp_factor = 1
params.show_edges = False
params.clim = [0, 0.2]                        # colourmap limits with format [0, 1] use "None" to set auto scaling
show_results(params)

from pathlib import Path
import re
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
 
# expects you already set:
# log_file (path), young_modulus, time_depths_factor
# optional: dimension (2 or 3). If not set, it's inferred.
 
# ----- helpers -----
FLOAT_RE = re.compile(r'[-+]?(?:\d+\.\d*|\.\d+|\d+)(?:[eE][-+]?\d+)?')
 
def strip_tags(line: str) -> str:
    # remove bracket/angle tags anywhere, then compress spaces
    line = re.sub(r'\[[^\]]*\]', ' ', line)   # [Essential], [0], [indent], etc.
    line = re.sub(r'<[^>]*>', ' ', line)      # <inform> etc.
    return ' '.join(line.split())
 
def floats_in(line: str):
    return [float(x) for x in FLOAT_RE.findall(line)]
 
def split_log(log_path: Path):
    reaction, ux, uy, uz = [], [], [], []
    with open(log_path, 'r', encoding='utf-8', errors='replace') as f:
        for raw in f:
            s = strip_tags(raw)
            if 'Total force' in s:
                reaction.append(s)
            if 'Ux' in s:
                ux.append(s)
            if 'Uy' in s:
                uy.append(s)
            if 'Uz' in s:
                uz.append(s)
    return reaction, ux, uy, uz
 
def parse_reaction(lines, dim=None):
    rows = []
    for line in lines:
        nums = floats_in(line)
        if not nums:
            continue
        tx, ty, tz = nums[-3], nums[-2], nums[-1]
        rows.append((tx, ty, tz))
    return pd.DataFrame(rows, columns=['tx', 'ty', 'tz'])
 
def parse_u(lines):
    rows = []
    for line in lines:
        s = line
        # try to read labeled min/max/time if present; otherwise fall back to heuristics
        m_min = re.search(r'\bmin\b[^\d\-+]*(' + FLOAT_RE.pattern + r')', s, flags=re.I)
        m_max = re.search(r'\bmax\b[^\d\-+]*(' + FLOAT_RE.pattern + r')', s, flags=re.I)
        m_tim = re.search(r'\btime\b[^\d\-+]*(' + FLOAT_RE.pattern + r')', s, flags=re.I)
        nums = floats_in(s)
        if not nums:
            continue
        time = float(m_tim.group(1)) if m_tim else nums[0]
        if m_min and m_max:
            vmin = float(m_min.group(1))
            vmax = float(m_max.group(1))
        elif len(nums) >= 3:
            # heuristic: last two numbers on the line are min/max (swap if needed)
            vmin, vmax = nums[-2], nums[-1]
            if vmin > vmax:
                vmin, vmax = vmax, vmin
        else:
            continue
        rows.append((time, vmin, vmax))
    return pd.DataFrame(rows, columns=['time', 'min', 'max'])
 
# ----- main -----
log_path = Path(log_file)
if not log_path.exists():
    raise FileNotFoundError(f"Log file not found: {log_path}")
 
reaction_lines, ux_lines, uy_lines, uz_lines = split_log(log_path)
 
data_reaction = parse_reaction(reaction_lines, dim=None)
if data_reaction.empty:
    raise ValueError("No usable 'Force' lines found (after stripping tags).")
 
# infer dimension from presence of tz
try:
    dimension
except NameError:
    dimension = 3 if data_reaction['tz'].notna().any() else 2
 
data_ux = parse_u(ux_lines)
data_uy = parse_u(uy_lines)
data_uz = parse_u(uz_lines) if dimension == 3 else pd.DataFrame()
 
if data_uy.empty:
    raise ValueError("No usable 'Uy' lines with numeric min/max found.")
 
# numeric arrays
yt = float(young_modulus)
tf = float(time_depths_factor)
 
t  = tf * data_ux['time'].to_numpy(dtype=float)
tx = yt * data_reaction['tx'].to_numpy(dtype=float)
ty = yt * data_reaction['ty'].to_numpy(dtype=float)
tz = yt * data_reaction['tz'].fillna(0.0).to_numpy(dtype=float) if dimension == 3 else None
 
# ----- plotting -----
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(7, 8), constrained_layout=True)
 
ax1.plot(t, tx, 'ro-', label='tx')
ax1.plot(t, ty, 'g^-', label='ty')
if dimension == 3 and np.any(~np.isnan(data_reaction['tz'])):
    ax1.plot(t, tz, 'b+-', label='tz')
ax1.set(xlabel='indentation depth', ylabel='force', title='Load path: indentation depth vs force')
ax1.grid(True)
ax1.legend(loc='upper left', shadow=True)
 
uy_min = data_uy['min'].to_numpy(dtype=float)
uy_max = data_uy['max'].to_numpy(dtype=float)
 
n1 = min(len(uy_min), len(tx))
n2 = min(len(uy_max), len(tx))
 
ax2.plot(uy_min[:n1], tx[:n1], 'ro-', label='min(Uy)-tx')
ax2.plot(uy_min[:n1], ty[:n1], 'g^-', label='min(Uy)-ty')
if dimension == 3 and tz is not None:
    ax2.plot(uy_min[:n1], tz[:n1], 'b+-', label='min(Uy)-tz')
 
ax2.plot(uy_max[:n2], tx[:n2], 'r.--', label='max(Uy)-tx')
ax2.plot(uy_max[:n2], ty[:n2], 'g.--', label='max(Uy)-ty')
if dimension == 3 and tz is not None:
    ax2.plot(uy_max[:n2], tz[:n2], 'b.--', label='max(Uy)-tz')
 
ax2.set(xlabel='displacement', ylabel='force', title='Load–displacement path: min[Uy] and max[Uy]')
ax2.grid(True)
ax2.legend(loc='upper left', shadow=True)
plt.show()