Update comments and explanations in Elasticity examples

examples/Elasticity/plot_1100.py

@@ -1,8 +1,18 @@
 """
 .. _ref_1100:
 
-Displacement control
-^^^^^^^^^^^^^^^^^^^^
+Displacement Control Analysis
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+This script models a linear elastic problem following the theory detailed in (:ref:`theory_elasticity`).
+
+Model Overview
+--------------
+The model represents a square plate where:
+- The bottom edge is fully fixed, preventing both displacement and rotation.
+- The top edge can slide vertically, with a controlled vertical displacement applied.
+
+The geometry, boundary conditions, and mesh are depicted in the diagram below. The plate is discretized into quadrilateral finite elements, appropriate for 2D structural analysis.
 
 .. code-block::
 
@@ -21,6 +31,11 @@
    #      ---X
    #  Z /
 
+
+Material Properties
+-------------------
+The material properties are summarized in the table below:
+
 +----+---------+--------+
 |    | VALUE   | UNITS  |
 +====+=========+========+
@@ -48,14 +63,21 @@
 # -------------------------------
 from phasefieldx.Element.Elasticity.Input import Input
 from phasefieldx.Element.Elasticity.solver.solver import solve
-from phasefieldx.Boundary.boundary_conditions import bc_x, bc_y, bc_xy, get_ds_bound_from_marker
-from phasefieldx.Loading.loading_functions import loading_Txy
+from phasefieldx.Boundary.boundary_conditions import bc_xy, get_ds_bound_from_marker
 from phasefieldx.PostProcessing.ReferenceResult import AllResults
 
 
 ###############################################################################
-# Parameters definition
-# ---------------------
+# Define Simulation Parameters
+# ----------------------------
+# `Data` is an input object containing essential parameters for simulation setup
+# and result storage:
+# - `E`: Young's modulus, set to 210 kN/mm².
+# - `nu`: Poisson's ratio, set to 0.3.
+# - `save_solution_xdmf` and `save_solution_vtu`: Set to `False` and `True`, respectively,
+#   specifying the file format to save displacement results (.vtu here).
+# - `results_folder_name`: Name of the folder for saving results. If it exists,
+#   it will be replaced with a new empty folder.
 Data = Input(E=210.0,
              nu=0.3,
              save_solution_xdmf=False,
@@ -64,31 +86,52 @@
 
 
 ###############################################################################
-# Mesh definition
+# Mesh Definition
 # ---------------
+# The mesh is a structured grid with quadrilateral elements:
+# - `divx`, `divy`: Number of elements along the x and y axes (10 each).
+# - `lx`, `ly`: Physical domain dimensions in x and y (1.0 units each).
+divx, divy = 10, 10
+lx, ly = 1.0, 1.0
 msh = dolfinx.mesh.create_rectangle(mpi4py.MPI.COMM_WORLD,
                                     [np.array([0, 0]),
-                                     np.array([1, 1])],
-                                    [10, 10],
+                                     np.array([lx, ly])],
+                                    [divx, divy],
                                     cell_type=dolfinx.mesh.CellType.quadrilateral)
 
 
+###############################################################################
+# Boundary Identification
+# -----------------------
+# Boundary conditions are applied to specific regions of the domain:
+# - `bottom`: Identifies the $y=0$ boundary.
+# - `top`: Identifies the $y=ly$ boundary.
+# `fdim` is the dimension of boundary facets (1D for a 2D mesh).
 def bottom(x):
     return np.isclose(x[1], 0)
 
-
 def top(x):
-    return np.isclose(x[1], 1)
-
+    return np.isclose(x[1], lx)
 
 fdim = msh.topology.dim - 1
 
+# Using the `bottom` and `top` functions, we locate the facets on the bottom and top sides of the mesh,
+# where $y = 0$ and $y = ly$, respectively. The `locate_entities_boundary` function returns an array of facet
+# indices representing these identified boundary entities.
 bottom_facet_marker = dolfinx.mesh.locate_entities_boundary(msh, fdim, bottom)
 top_facet_marker = dolfinx.mesh.locate_entities_boundary(msh, fdim, top)
 
-ds_bottom = get_ds_bound_from_marker(top_facet_marker, msh, fdim)
+# The `get_ds_bound_from_marker` function generates a measure for applying boundary conditions 
+# specifically to the facets identified by `top_facet_marker` and `bottom_facet_marker`, respectively. 
+# This measure is then assigned to `ds_bottom` and `ds_top`.
+ds_bottom = get_ds_bound_from_marker(bottom_facet_marker, msh, fdim)
 ds_top = get_ds_bound_from_marker(top_facet_marker, msh, fdim)
 
+# `ds_list` is an array that stores boundary condition measures along with names 
+# for each boundary, simplifying result-saving processes. Each entry in `ds_list` 
+# is formatted as `[ds_, "name"]`, where `ds_` represents the boundary condition measure, 
+# and `"name"` is a label used for saving. Here, `ds_bottom` and `ds_top` are labeled 
+# as `"bottom"` and `"top"`, respectively, to ensure clarity when saving results.
 ds_list = np.array([
                    [ds_bottom, "bottom"],
                    [ds_top, "top"]
@@ -106,29 +149,80 @@ def top(x):
 ###############################################################################
 # Boundary Conditions
 # -------------------
-# Apply boundary conditions: bottom nodes fixed in both directions, top nodes can slide vertically.
-bc_bottom = bc_xy(bottom_facet_marker, V_u, fdim)
-bc_top = bc_xy(top_facet_marker, V_u, fdim)
+# Dirichlet boundary conditions are applied as follows:
+# - `bc_bottom`: Fixes x and y displacement to 0 on the bottom boundary.
+# - `bc_top`: Fixes x displacement and allows variable y displacement on the top
+bc_bottom = bc_xy(bottom_facet_marker, V_u, fdim, value_x=0.0, value_y=0.0)
+bc_top = bc_xy(top_facet_marker, V_u, fdim, value_x=0.0, value_y=0.0)
+
+# The bcs_list_u variable is a list that stores all boundary conditions for the displacement
+# field $\boldsymbol u$. This list facilitates easy management of multiple boundary
+# conditions and can be expanded if additional conditions are needed.
 bcs_list_u = [bc_top, bc_bottom]
 
 
+# Function: `update_boundary_conditions`
+# --------------------------------------
+# The `update_boundary_conditions` function dynamically updates the boundary conditions at each
+# time step, enabling a quasi-static analysis by applying incremental displacements to specific
+# degrees of freedom.
+#
+# Parameters:
+# - `bcs`: List of boundary conditions, where each entry corresponds to a boundary condition applied
+#   to a specific facet of the mesh.
+# - `time`: Scalar representing the current time step in the analysis.
+#
+# Inside the function:
+# - `val` is calculated as a linear function of `time`, specifically `val = 0.0003 * time`,
+#   to simulate gradual displacement along the y-axis. This can be modified as needed for different
+#   quasi-static loading schemes.
+# - The value `val` is assigned to the y-component of the displacement field on the boundary,
+#   achieved by updating `bcs[0].g.value[1]`, where `bcs[0]` represents the top boundary condition.
+#
+# Returns:
+# - A tuple `(0, val, 0)` where:
+#   - The first element is zero (indicating no update for the x component in this example).
+#   - The second element is `val`, the calculated y-displacement.
+#   - The third element is zero (indicating no z-component displacement in this 2D example).
+#
+# This function supports the quasi-static analysis by gradually updating the displacement boundary
+# condition over time, allowing for controlled loading in the simulation.
 def update_boundary_conditions(bcs, time):
     val = 0.0003 * time
     bcs[0].g.value[1] = petsc4py.PETSc.ScalarType(val)
     return 0, val, 0
 
-
 T_list_u = None
 update_loading = None
 f = None
 
 
 ###############################################################################
-# Call the solver
-# ---------------
-# We set a final time of t=10 and a time step of \Delta t=1.
-# Note that this is a linear boundary value problem.
-
+# Solver Call for a Static Linear Problem
+# ---------------------------------------
+# We define the parameters for a simple, static linear boundary value problem 
+# with a final time `t = 10.0` and a time step `Δt = 1.0`. Although this setup 
+# includes time parameters, they are primarily used for structural consistency 
+# with a generic solver function and do not affect the result, as the problem 
+# is linear and time-independent.
+#
+# Parameters:
+# - `final_time`: The end time for the simulation, set to 10.0.
+# - `dt`: The time step for the simulation, set to 1.0. In a static context, this
+#   only provides uniformity with dynamic cases but does not change the results.
+# - `path`: Optional path for saving results; set to `None` here to use the default.
+# - `quadrature_degree`: Defines the accuracy of numerical integration; set to 2
+#   for this problem.
+#
+# Function Call:
+# The `solve` function is called with:
+# - `Data`: Simulation data and parameters.
+# - `msh`: Mesh of the domain.
+# - `V_u`: Function space for $\boldsymbol u$.
+# - `bcs_list_u`: List of boundary conditions.
+# - `update_boundary_conditions`, `update_loading`: update the boundary condition for the quasi static analysis
+# - `ds_list`: Boundary measures for integration on specified boundaries.
+# - `dt` and `final_time` to define the static solution time window.
 final_time = 10.0
 dt = 1.0
 
@@ -142,7 +236,9 @@ def update_boundary_conditions(bcs, time):
       T_list_u,
       update_loading,
       ds_list,
-      dt)
+      dt,
+      path=None,
+      quadrature_degree=2)
 
 
 ###############################################################################
@@ -181,7 +277,7 @@ def update_boundary_conditions(bcs, time):
 # Plot the vertical displacement versus the reaction force.
 fig, ax = plt.subplots()
 
-ax.plot(displacement, S.reaction_files['top.reaction']["Ry"], S.color + '.', linewidth=2.0, label=S.label)
+ax.plot(displacement, S.reaction_files['bottom.reaction']["Ry"], S.color + '.', linewidth=2.0, label=S.label)
 
 ax.grid(color='k', linestyle='-', linewidth=0.3)
 ax.set_xlabel('displacement - u $[mm]$')