# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     name: python3
# ---

# + [markdown] colab_type="text" id="view-in-github"
# <a href="https://colab.research.google.com/github/kumiori/mec647/blob/main/mec647_Fracture_4.ipynb" target="_parent"><img src="https://colab.research.google.com/assets/colab-badge.svg" alt="Open In Colab"/></a>

# + id="rA4M7CiLCsX_"
# %%capture
import sys

try:
    import google.colab  # noqa: F401
except ImportError:
    import ufl  # noqa: F401
    import dolfinx  # noqa: F401
else:
    try:
        import ufl
        import dolfinx
    except ImportError:
        # !wget "https://fem-on-colab.github.io/releases/fenicsx-install.sh" -O "/tmp/fenicsx-install.sh" && bash "/tmp/fenicsx-install.sh";
        import ufl  # noqa: F401
        import dolfinx  # noqa: F401

try:
    import pyvista
except ImportError:
    # !{sys.executable} -m pip install --upgrade pyvista itkwidgets;
    import pyvista  # noqa: F401
    from pyvista.utilities import xvfb

try:
    import gmsh
except ImportError:
    # !{sys.executable} -m pip install gmsh
    import gmsh

# !sudo apt install libgl1-mesa-glx xvfb;
# !{sys.executable} -m pip install pythreejs;
# !{sys.executable} -m pip install ipygany;
# !{sys.executable} -m pip install --upgrade PyYAML
# # !pip install --ignore-installed PyYAML



# + id="xPTwn6bxC4S5"
# %%capture
# !sudo apt install libgl1-mesa-glx xvfb;
# !{sys.executable} -m pip install pythreejs;
# !{sys.executable} -m pip install ipygany;
# !{sys.executable} -m pip install --upgrade PyYAML
# # !pip install --ignore-installed PyYAML



# + [markdown] id="QhnaGiklC3nQ"
# # Fracture
#
#
# Let $\Omega \subset (0, L)^D$, with $D=1, 2, 3$, $L$ finite, being the (or one) characteristic length of the specimen.
# For any 
# - displacement field $u\in V_t : H^1(\Omega, R^n) + bcs(t)$ with $n=1, 2$ or $3$, and
# - damage field $\alpha \in H^1(\Omega, R)$,
#
# consider the energy $E(u, \alpha)$ defined as
# $$
# E_\ell(u, \alpha)=\frac{1}{2}\int_\Omega a(\alpha) W(u) dx + \underbrace{\frac{G_c}{c_w} \int \left(\frac{1}{\ell}w(\alpha) + \ell |\nabla \alpha|^2 \right)dx}_{\text{Surface energy}}- \int_\Omega f.u dx$$
#
# In practice, $\ell \ll L$.
#
# Above, $W$ is the elastic energy density, reading (in linearised elasticity as) 
# $$ 
# W(u) = Ae(u):e(u)
# $$
# where $A$ is the 4-th order tensor of elasticity, in the isotropic and homogeneous case, it corresponds to a linear combination with two coefficients, say, $A_0$ the stiffness (dimensional), and $\nu$ the Poisson ratio (non-dimensional). 
#
#
# Further, $w(\alpha)$ corresponds to the dissipated energy to damage, homogeneously, the specimen, the gradient term accounts for spatial variations.
#
# **Keypoint:** these two terms are weighted by $\ell$, a parameter that is homogeneous to a length and is understood as a *material* quantity (as opposed to *numerical*).
#
# Define $D(\alpha_0):=\left\{ \alpha \in H^1(\Omega), \alpha \geq \alpha_0 \right\}$, for some $\alpha_0(x)\geq 0$ pointwise. 
#
# We solve two types of problems (by increasing difficulty):
# - **The static problem**: Given a load (boundary conditions) and an initial state of damage $\alpha_0$, what is the equilibrium displacement and repartition of damage?
# In other terms:
#  
# $
# \text{ min loc} \left\{ E_\ell(u, \alpha):
#  u \in V_t, \alpha \in D(\alpha_0) \right\}.
# $
#
# - **The evolution problem**: Given a load **history** (boundary conditions as a function of $t$) and an initial state of damage $\alpha_0$, what is the *evolution* of equilibrium displacement and repartition of damage, i.e. 
# the map $t\mapsto (u_t, \alpha_t)$, such that 
#   - (Irrevers.) $\alpha_t \nearrow t$,
#   - (Stability) $(u_t, \alpha_t) = \operatorname{arg loc min} \left\{ E_\ell (v, \beta), (v, \beta) \in V_t \times D(\alpha_t) \right\}$  -->
#
#

# + [markdown] id="WJD-kxygsGbY"
# ### Parameters
#
# In the energy above:
# - Two elasticity parameters, such as
#   - $A_0$ the stiffness of the sound material
#   - $\nu$ the Poisson ratio
#   - equivalently, $\mu, \lambda$ the Lamé parameters
# - Two fracture/damage parameters:
#   - $\ell$ the internal damage length
#   - $G_c$ the material toughness
#
#
# ### Back of the envelope computation.
#
# 1. Show that the energy above can be written as a function of only two non-dimensional parameters (ex.: $\nu, \tilde \ell)$, by dimensional analysis.

# + colab={"base_uri": "https://localhost:8080/"} id="FwFJxjZRDD1A" outputId="d1e779e8-742e-429b-9c76-afe215535da2"
branch_name = 'main'

# !rm -rf mec647
try:
  # !git clone -b {branch_name} https://github.com/kumiori/mec647.git
except Exception:
  print('Something went wrong')
  # !rm -rf mec647
  # !git clone https://github.com/kumiori/mec647.git


# + id="B_A6SI2gl7DN"
sys.path.append('mec647/')


# + id="y6YzSDfJ3Vwz"
import matplotlib.pyplot as plt
from petsc4py import PETSc
import numpy as np
from mpi4py import MPI

comm = MPI.COMM_WORLD

from dolfinx.fem import (
    Constant,
    Function,
    FunctionSpace,
    assemble_scalar,
    dirichletbc,
    form,
    locate_dofs_geometrical,
    set_bc,
)

# meshes
import meshes
from meshes import primitives

# + id="xQC3q28P3Xmv"
# Parameters

parameters = {
    'loading': {
        'min': 0,
        'max': 1.5,
        'steps': 10
    },
    'geometry': {
        'geom_type': 'bar',
        'Lx': 1.,
        'Ly': 0.1
    },
    'model': {
        'E': 1.,
        'nu': 0.,
        'w1': 1.,
        'ell': 0.05,
        'k_res': 1.e-8
    },
    'solvers': {
          'elasticity': {        
            'snes': {
                'snes_type': 'newtontr',
                'snes_stol': 1e-8,
                'snes_atol': 1e-8,
                'snes_rtol': 1e-8,
                'snes_max_it': 100,
                'snes_monitor': "",
                'ksp_type': 'preonly',
                'pc_type': 'lu',
                'pc_factor_mat_solver_type': 'mumps'
            }
        },
          'damage': {        
            'snes': {
                'snes_type': 'vinewtonrsls',
                'snes_stol': 1e-5,
                'snes_atol': 1e-5,
                'snes_rtol': 1e-8,
                'snes_max_it': 100,
                'snes_monitor': "",
                'ksp_type': 'preonly',
                'pc_type': 'lu',
                'pc_factor_mat_solver_type': 'mumps'
            },
        },
          'damage_elasticity': {
            "max_it": 100,
            "alpha_rtol": 1.0e-5,
            "criterion": "alpha_H1"
          }
    }
}

# + id="2hjc-VWztAgK"
# Equivalently
import yaml

with open("mec647/test/parameters.yml") as f:
    parameters = yaml.load(f, Loader=yaml.FullLoader)


# + colab={"base_uri": "https://localhost:8080/"} id="wRJtcLKk3iEC" outputId="3258572f-0cbc-4a62-bd87-386952c81018"
parameters.get("model").get("ell")

# + colab={"base_uri": "https://localhost:8080/", "height": 84} id="uGvOIlvU3d0-" outputId="fb961ce3-a5b8-4b20-8a3b-306692f4b62f"
# Mesh

Lx = parameters["geometry"]["Lx"]
Ly = parameters["geometry"]["Ly"]
geom_type = parameters["geometry"]["geom_type"]

gmsh_model, tdim = primitives.mesh_bar_gmshapi(geom_type,
                                    Lx, 
                                    Ly, 
                                    parameters.get("model").get("ell")/3, 
                                    tdim=2)

mesh, mts = meshes.gmsh_model_to_mesh(gmsh_model,
                               cell_data=False,
                               facet_data=True,
                               gdim=2)


from utils.viz import plot_mesh

plt.figure()
ax = plot_mesh(mesh)
fig = ax.get_figure()
fig.savefig(f"mesh.png")


# + id="YbtAhOV538kj"
# Functional Setting

element_u = ufl.VectorElement("Lagrange", mesh.ufl_cell(),
                              degree=1, dim=2)

element_alpha = ufl.FiniteElement("Lagrange", mesh.ufl_cell(),
                              degree=1)

V_u = dolfinx.fem.FunctionSpace(mesh, element_u) 
V_alpha = dolfinx.fem.FunctionSpace(mesh, element_alpha) 

u = dolfinx.fem.Function(V_u, name="Displacement")
u_ = dolfinx.fem.Function(V_u, name="BoundaryDisplacement")


alpha = dolfinx.fem.Function(V_alpha, name="Damage")

# Pack state
state = {"u": u, "alpha": alpha}

# Bounds
alpha_ub = dolfinx.fem.Function(V_alpha, name="UpperBoundDamage")
alpha_lb = dolfinx.fem.Function(V_alpha, name="LowerBoundDamage")



dx = ufl.Measure("dx", domain = mesh)
ds = ufl.Measure("ds", domain = mesh)

# + id="sEJ-jkLR4Wgc"
# Boundary sets

dofs_alpha_left = locate_dofs_geometrical(V_alpha, lambda x: np.isclose(x[0], 0.))
dofs_alpha_right = locate_dofs_geometrical(V_alpha, lambda x: np.isclose(x[0], Lx))

dofs_u_left = locate_dofs_geometrical(V_u, lambda x: np.isclose(x[0], 0.))
dofs_u_right = locate_dofs_geometrical(V_u, lambda x: np.isclose(x[0], Lx))

# + id="9opnlp5F6QRV"
# Boundary data

u_.interpolate(lambda x: (np.ones_like(x[0]), np.zeros_like(x[1])))

# Bounds (nontrivial)

alpha_lb.interpolate(lambda x: np.zeros_like(x[0]))
alpha_ub.interpolate(lambda x: np.ones_like(x[0]))

# + id="O0K_yA_G4glE"
# Boundary conditions

bcs_u = [
         dirichletbc(np.array([0., 0.], dtype=PETSc.ScalarType),
                      dofs_u_left,
                      V_u),
         dirichletbc(u_, dofs_u_right)
         ]

bcs_alpha = [
             dirichletbc(np.array(0., dtype = PETSc.ScalarType),
                         np.concatenate([dofs_alpha_left, dofs_alpha_right]),
                         V_alpha)
]


# + id="nldNjjoF7h1X"
import models
from models import DamageElasticityModel as Brittle

# + id="qi88nswf4iY1"
bcs = {"bcs_u": bcs_u, "bcs_alpha": bcs_alpha}

model = Brittle(parameters["model"])


total_energy = model.total_energy_density(state) * dx

# + id="XgSMuKMq8SlY"
import algorithms
from algorithms import am

# + id="DKcj0EB48dz4"
solver = am.AlternateMinimisation(total_energy,
                         state,
                         bcs,
                         parameters.get("solvers"),
                         bounds=(alpha_lb, alpha_ub)
                         )

# + colab={"base_uri": "https://localhost:8080/"} id="n6gLZCY9-TSC" outputId="111c1bc4-069e-4b55-89fd-212a00bae095"
# Loop for evolution


loads = np.linspace(parameters.get("loading").get("min"),
                    parameters.get("loading").get("max"),
                    parameters.get("loading").get("steps"))

data = {
    'elastic': [],
    'surface': [],
    'total': [],
    'load': []
}
for (i_t, t) in enumerate(loads):
  # update boundary conditions

  u_.interpolate(lambda x: (t * np.ones_like(x[0]), np.zeros_like(x[1])))
  u_.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT,
                        mode=PETSc.ScatterMode.FORWARD)

  # update lower bound for damage
  alpha.vector.copy(alpha_lb.vector)
  alpha.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT,
                        mode=PETSc.ScatterMode.FORWARD)

  # solve for current load step
  print(f"Solving timestep {i_t}, load: {t}")

  solver.solve()

  # postprocessing
  # global

  surface_energy = comm.allreduce(
          assemble_scalar(
            dolfinx.fem.form(model.damage_dissipation_density(state) * dx)),
          op=MPI.SUM
          )

  elastic_energy = comm.allreduce(
      assemble_scalar(
        dolfinx.fem.form(model.elastic_energy_density(state) * dx)),
          op=MPI.SUM
        )
  
  data.get('elastic').append(elastic_energy)
  data.get('surface').append(surface_energy)
  data.get('total').append(surface_energy+elastic_energy)
  data.get('load').append(t)

  print(f"Solved timestep {i_t}, load: {t}")
  print(f"Elastic Energy {elastic_energy:.3g}, Surface energy: {surface_energy:.3g}")
  print("\n\n")

  # # savings?








# + colab={"base_uri": "https://localhost:8080/", "height": 333} id="9ePMBvwzmzZ2" outputId="bf8bff5e-324f-4ff2-e92c-7f25ccaf5a37"
plt.plot(data.get('load'), data.get('surface'), label='surface')
plt.plot(data.get('load'), data.get('elastic'), label='elastic')
plt.plot(data.get('load'), [1./2. * t**2*Ly for t in data.get('load')], label='anal elast', ls=':', c='k')

plt.title('Traction bar energetics')
plt.legend()
plt.yticks([0, 1/20], [0, '$1/2.\sigma_c^2/E_0$'])
plt.xticks([0, 1], [0, 1])

# + id="w4p1xwzGoL1a"
# experiments
params0, data0 = parameters, data


# +
from test_viz import plot_vector, plot_scalar, plot_profile

_ell = parameters.get('model').get('ell')
tol = 1e-3
xs = np.linspace(0 + tol, Lx - tol, 101)
points = np.zeros((3, 101))
points[0] = xs

_plt, data = plot_profile(
    alpha,
    points,
    subplot=(0, 0),
    lineproperties={
        "c": "k",
        "label": f"$alpha$ with $\ell$ = {_ell:.2f}"
    },
)
ax = _plt.gca()
ax.axvline(0.0, c="k")
ax.axvline(2 * _ell, c="k", label='D=$2\ell$')
_plt.legend()
_plt.fill_between(data[0], data[1].reshape(len(data[1])))
_plt.title("Traction bar damage profile")

# -



StarcoderdataJupyterScriptsDedupFiltered