r"""
.. _ref_ex_advanced2:

Mesh Refinement
---------------
Perform a mesh refinement study.

In this example the convergence of the torsion constant is investigated through an analysis of an
I Section. The mesh is refined both by modifying the mesh size and by specifying the number of
points making up the root radius. The figure below the example code shows that mesh refinement
adjacent to the root radius is a far more efficient method in obtaining fast convergence when
compared to reducing the mesh area size for the entire section.
"""

# sphinx_gallery_thumbnail_number = 1

import numpy as np
import matplotlib.pyplot as plt
import sectionproperties.pre.library.steel_sections as steel_sections
from sectionproperties.analysis.section import Section

# %%
# Define mesh sizes
mesh_size_list = [200, 100, 50, 20, 10, 5]
nr_list = [4, 8, 12, 16, 20, 24, 32]

# %%
# Initialise result lists
mesh_results = []
mesh_elements = []
nr_results = []
nr_elements = []

# %%
# Calculate reference solution
geometry = steel_sections.i_section(d=203, b=133, t_f=7.8, t_w=5.8, r=8.9, n_r=32)
geometry.create_mesh(mesh_sizes=[5])  # create mesh
section = Section(geometry)  # create a Section object
section.calculate_geometric_properties()
section.calculate_warping_properties()
j_reference = section.get_j()  # get the torsion constant

# %%
# Run through mesh_sizes with n_r = 8
for mesh_size in mesh_size_list:
    geometry = steel_sections.i_section(d=203, b=133, t_f=7.8, t_w=5.8, r=8.9, n_r=8)
    geometry.create_mesh(mesh_sizes=[mesh_size])  # create mesh
    section = Section(geometry)  # create a Section object
    section.calculate_geometric_properties()
    section.calculate_warping_properties()

    mesh_elements.append(len(section.elements))
    mesh_results.append(section.get_j())

# %%
# Run through n_r with mesh_size = 10
for n_r in nr_list:
    geometry = steel_sections.i_section(d=203, b=133, t_f=7.8, t_w=5.8, r=8.9, n_r=n_r)
    geometry.create_mesh(mesh_sizes=[10])  # create mesh
    section = Section(geometry)  # create a Section object
    section.calculate_geometric_properties()
    section.calculate_warping_properties()

    nr_elements.append(len(section.elements))
    nr_results.append(section.get_j())

# %%
# Convert results to a numpy array and compute the error
mesh_results = np.array(mesh_results)
nr_results = np.array(nr_results)
mesh_error_vals = (mesh_results - j_reference) / mesh_results * 100
nr_error_vals = (nr_results - j_reference) / nr_results * 100

# %%
# Plot the results
(fig, ax) = plt.subplots()
ax.loglog(mesh_elements, mesh_error_vals, "kx-", label="Mesh Size Refinement")
ax.loglog(nr_elements, nr_error_vals, "rx-", label="Root Radius Refinement")
plt.xlabel("Number of Elements")
plt.ylabel("Torsion Constant Error [%]")
plt.legend(loc="center left", bbox_to_anchor=(1, 0.5))
plt.tight_layout()
plt.show()
