"""
Mesh quality improvement
========================

See also: :mod:`~mymesh.improvement`, :mod:`~mymesh.quality`, :ref:`theory_improvement`

High quality meshes are essential for accuracy in simulations (e.g. finite element, finite volume methods) and other mesh-based computation. 
There are different approaches to improving mesh quality, including smoothing and modifications to the mesh connectivity (see also: :ref:`theory_improvement`).
Often, the best approach is to use a combination of these methods :cite:p:`Klingner2008`.

The :func:`mymesh.improvement.Improve` function combines edge contraction for coarsening (:func:`~mymesh.improvement.Contract`), edge splitting for refinement (:func:`~mymesh.improvement.Split`), edge flipping (:func:`~mymesh.improvement.Flip`), and smoothing to improve mesh quality while also achieving a target edge length for triangular or tetrahedral meshes.

:func:`~mymesh.improvement.Improve` repeats a schedule of operations, which is by default splitting ('s'), contraction ('c'), flipping ('f'), and smoothing ('S') defined by a string (``'scfS'``), however the schedule can be customized to remove or rearrange the order of operations (e.g. ``'Sfc'``).

"""

#%%
# 2D Surface Improvement
# ----------------------
# :func:`~mymesh.improvement.Improve` can be used on a two dimensional 
# triangular surface while preserving the original boundaries of the mesh 

from mymesh import implicit, improvement, primitives, quality, visualize
import numpy as np

C = primitives.Circle([0,0,0], 1, ElemType='tri', 
                        theta_resolution=40, radial_resolution=20)
C.ElemData['skew'] = quality.Skewness(*C)

target_edge_length = 0.1
C2 = improvement.Improve(C, target_edge_length, schedule = 'scfS', 
                         repeat=10, verbose=False)
C2.ElemData['skew'] = quality.Skewness(*C2)

means = [np.mean(C.ElemData['skew']), np.mean(C2.ElemData['skew'])]

C.plot(scalars='skew', clim=(0,1), show_edges=True, view='xy', 
       title=f'Original\nMean Skewness={means[0]:.2f}')
C2.plot(scalars='skew', clim=(0,1), show_edges=True, view='xy',
        title=f'Improved\nMean Skewness={means[1]:.2f}')

#%%
# 3D Surface Improvement
# ----------------------
# Similarly, three dimensional triangular surfaces can be improved. 
# Implicit or image-based surfaces often contain some highly skewed elements
# as a result of the contouring process, which can be improved by 
# :func:`~mymesh.improvement.Improve`.

S = implicit.SurfaceMesh(implicit.sphere([0,0,0], 1), [-1,1,-1,1,-1,1], .1)
S.ElemData['skew'] = quality.Skewness(*S)

target_edge_length = 0.2
S2 = improvement.Improve(S, target_edge_length, schedule='scfS', 
                         repeat=10, verbose=False)
S2.ElemData['skew'] = quality.Skewness(*S2)

means = [np.mean(S.ElemData['skew']), np.mean(S2.ElemData['skew'])]
S.plot(scalars='skew', clim=(0,1), show_edges=True, view='trimetric', 
       title=f'Original\nMean Skewness={means[0]:.2f}')
S2.plot(scalars='skew', clim=(0,1), show_edges=True, view='trimetric',
        title=f'Improved\nMean Skewness={means[1]:.2f}')
#%%
# 3D Volume Improvement
# ---------------------
# The tetrahedral meshes generated by :func:`~mymesh.contour.MarchingTetrahedra`
# (:func:`mymesh.implicit.TetMesh`, :func:`mymesh.image.TetMesh`) generally 
# structured interiors with some low quality elements on the surface, however,
# even the structured interiors have elements of only 'average' quality. 

S = implicit.TetMesh(implicit.sphere([0,0,0], 1), [-1,1,-1,1,-1,1], .1)
S.ElemData['skew'] = quality.Skewness(*S)

target_edge_length = 0.25
S2 = improvement.Improve(S, target_edge_length, schedule='scfS', 
                         repeat=20, verbose=False)
S2.ElemData['skew'] = quality.Skewness(*S2)

means = [np.mean(S.ElemData['skew']), np.mean(S2.ElemData['skew'])]
S.Clip().plot(scalars='skew', clim=(0,1), show_edges=True, view='trimetric', 
       title=f'Original\nMean Skewness={means[0]:.2f}')
S2.Clip().plot(scalars='skew', clim=(0,1), show_edges=True, view='trimetric',
        title=f'Improved\nMean Skewness={means[1]:.2f}')

# sphinx_gallery_thumbnail_number = 2
# %%