''' 
Virtual Fields Method Verification
==================================
In this example, we verify the virtual 
power residual calculation for a 
boundary value problem (BVP) with an 
analytical solution. To conform 
to our VFM implementation constraints, 
we chose a thin plate of length :math:`H`
loaded in uniaxial tension in the Y 
direction with a distributed load with 
an over all magnitude of
:math:`L`. 
Assuming linear elasticity and 
small deformation, the stress
in the Y direction in the plate is 

.. math::

    \sigma_{yy} =\\frac{L}{W T}

where :math:`W` is the plate 
width and :math:`T` is the 
plate thickness. All other 
stresses are zero. 
The in-plane displacements in the plate 
are given by 

.. math::
    
    u=-\\nu\\frac{\sigma_{yy}}{E}x

and 

.. math::
    
    v=\\frac{\sigma_{yy}}{E}y

where :math:`E` is the material elastic modulus, 
:math:`\\nu` is the material Poisson's ratio, and
:math:`x` and :math:`y` are the locations on the plate. 

The internal virtual power becomes

.. math::
    
    P_{int}=\\frac{L}{Wt}\int_{z=-T/2}^{T/2}\int_{y=-H/2}^{H/2}\int_{x=-W/2}^{W/2}\dot{\mathbf{F}}^*_{yy}dxdydz

where :math:`\dot{\mathbf{F}}^*_{yy}` is the :math:`\\frac{\partial\mathbf{v}^*}{\partial y}`. 
The virtual velocity function used is

.. math::

    \mathbf{v}^*_y=\\frac{2y+H}{2H}


With this virtual velocity function, the virtual 
velocity gradient component of interest is 

.. math::

    \\frac{\partial\mathbf{v}^*}{\partial y} = 1/h.

As a result, the internal virtual power is given as 


.. math::
    
    P_{int}=L

The external virtual power, as expected, 
results an equivalent value. 

.. math::

    P_{ext} = L\mathbf{v}^*_y = L\\frac{2H/2+H}{2H} = L

With an analytical solution for the 
virtual internal and external powers 
available for the BVP
of interest, we will create MatCal VFM 
objects with the proper input and verify 
that they return the expect values. 

We will begin by creating synthetic data 
that represents the experimental data 
for this BVP. To do so, we import 
MatCal's tools and create the variables needed
for the analytical solutions above. 
We will need:

#.  The materials elastics constants :math:`E`
    and :math:`\\nu`. We choose values that 
    are similar to steel with :math:`E=200` GPa
    and :math:`\\nu=0.25`.
#.  The plate dimensions :math:`H=15.2` mm, :math:`W=7.6`, 
    and :math:`T=0.1` mm with a discretization in Y and Y.
#.  The load magnitude 
    as a function of time :math:`L(t)`. 
    We will choose :math:`L(t)` such that 
    the maximum stress in the plate is 100 MPa.
#.  The number of time steps and the time values.
'''
from matcal import *
import numpy as np

E = 200e9
nu = 0.25

H = 15.2e-3
W = 7.6e-3
T = 0.1e-3

measured_num_points = 10
measured_xs = np.linspace(0, W, measured_num_points)
measured_ys = np.linspace(0, H, measured_num_points)
measured_x_grid, measured_y_grid = np.meshgrid(measured_xs, measured_ys)


n_time_steps = 10
L = np.linspace(0, 100e6*W*T,n_time_steps).reshape(n_time_steps,1)
time = np.linspace(0,1,n_time_steps)

# %%
# With all the required inputs initialized, we can 
# create the :class:`~matcal.full_field.data.FieldData` objects
# that are need to have MatCal evaluate the virtual powers
# for the problem. The :class:`~matcal.full_field.data.FieldData`
# object requires several fields:
#
# #. time
# #. load
# #. The X locations of the field data
# #. The Y locations of the data
# #. The X displacements for each point (U)
# #. The Y displacements for each point (V)

sigma_yy = L/(W*T)

def displacements(x,y):
    return  -nu*sigma_yy/E*x, sigma_yy/E*y

all_x_1D = measured_x_grid.reshape(measured_num_points**2)
all_y_1D = measured_y_grid.reshape(measured_num_points**2)
u,v = displacements(all_x_1D, all_y_1D)

# %%
# After creating the values for the 
# fields that we need, they can be combined
# into a dictionary and converted into 
# a MatCal :class:`~matcal.full_field.data.FieldData`
# object using the 
# :func:`~matcal.full_field.data.convert_dictionary_to_field_data`
# function.

data_dict = {'time':time, 'load':L, 
             "X":all_x_1D, 
             "Y":all_y_1D, 
             "U":u, 
             "V":v}

field_data = convert_dictionary_to_field_data(data_dict, 
                                              coordinate_names=["X", "Y"])


# %%
# Next, a discretization for 
# the VFM models must be created. 
# MatCal has a simple tool
# for creating rectangular volumes
# that we will use, the 
# :func:`~matcal.full_field.TwoDimensionalFieldGrid.auto_generate_two_dimensional_field_grid`
# function. It requires the number of X and Y discretization points and the field data 
# that the mesh is intended to encapsulate as inputs. 
# It will return MatCal's two dimensional mesh class that can be
# used as a mesh input to MatCal's VFM model. It currently 
# can only be used for rectangular shapes without holes.

from matcal.full_field.TwoDimensionalFieldGrid import auto_generate_two_dimensional_field_grid
auto_mesh = auto_generate_two_dimensional_field_grid(measured_num_points*5, 
                                                     measured_num_points*5, field_data)

#%%
# The next item needed for the creation of the 
# MatCal VFM models is a SierraSM material file. 
# We will create a simple elastic material model file
# using Python file tools and then create a 
# MatCal :class:`~matcal.sierra.material.Material`
# object so the VFM models can use that material.

material_file_string = \
"""
  begin material steel
    begin parameters for model elastic
      youngs modulus                =  {E}
      poissons ratio                =   {nu}
    end
  end
"""

with open("elastic_material.inc", "w") as mf:
    mf.write(material_file_string)

material = Material("steel", "elastic_material.inc", "elastic")

#%%
# We can now create the VFM models that we want 
# to evaluate. We will look at two VFM models,
# The :class:`~matcal.sierra.models.VFMUniaxialTensionHexModel`
# and the :class:`~matcal.sierra.models.VFMUniaxialTensionConnectedHexModel`
# will both be evaluated in this verification example. 
# The models require a material, a mesh and thickness 
# to be initialized and field data with a two-dimensional displacement field 
# to be added as the boundary condition data. We optionally set the names
# so their results can be more easily pulled from the results 
# data object returned from our study.

default_hex_model = VFMUniaxialTensionHexModel(material, auto_mesh, T)
default_hex_model.set_name("default_VFM_hex_model")
default_hex_model.add_boundary_condition_data(field_data)
default_hex_model.set_number_of_time_steps(20)

connected_hex_model = VFMUniaxialTensionConnectedHexModel(material, auto_mesh, T)
connected_hex_model.set_name("VFM_connected_hex_model")
connected_hex_model.add_boundary_condition_data(field_data)
connected_hex_model.set_number_of_time_steps(20)


#%%
# The final step is to create the parameters for the study, 
# the objective to be evaluated 
# and the study itself. For our elastic model, the 
# only parameters we need are the elastic modulus and 
# the Poisson's ratio. The study we will use to evaluate
# the models is the :class:`~matcal.core.parameter_studies.ParameterStudy`
# since just need the model results at our predetermined 
# steel-like parameter values. Since this study is 
# evaluating a VFM model, we create a VFM model 
# with no inputs because our global data fields 
# of "time" and "load" match the default ones 
# expected by the objective.
E_param = Parameter('E', 100e9, 300e9)
nu_param = Parameter('nu', 0.1, 0.5)

vfm_objective = MechanicalVFMObjective()
vfm_objective.set_name('vfm_objective')

#%%
# After creating the parameters, objectives, synthetic data
# and models, we can create the study and add the evaluation 
# sets and parameter values that we want to evaluate. 
# The study is then launched and the results are stored in 
# an object that we can analyze when the study completes.

study = ParameterStudy(E_param, nu_param)
study.add_parameter_evaluation(E=E, nu=nu)
study.add_evaluation_set(default_hex_model, vfm_objective, field_data)
study.add_evaluation_set(connected_hex_model, vfm_objective, field_data)
study.set_core_limit(3)
results = study.launch()

#%%
# In this verification problem, our 
# goal is to verify that our MatCal 
# VFM objective is accurately evaluating 
# the internal and external virtual powers. 
# As a result, we extract these virtual powers 
# from the study results and compare them to 
# the analytically determined value for the internal power 
# :math:`L\left(t\right)`. We compare them by 
# evaluating their percent error with the following equation
#
# .. math::
#    Pe =  100\frac{P_{eval} - L}{\max\left(L\right)}
#
# where :math:`P_{eval}` is the evaluated virtual power 
# (either the internal or external from either model).
# Since we will evaluate this percent error several times,
# we make a function to perform the calculation.

def percent_error(prediction, actual):
    return (prediction - actual)/np.max(actual)*100

#%%
# We then extract the results from each of the models
# and print or plot their errors.

matcal_internal_power_default = results.best_simulation_qois(default_hex_model, 
                                                             vfm_objective, 
                                                             field_data.state, 
                                                             0)["virtual_power"]
matcal_external_power_default = results.get_experiment_qois(default_hex_model, 
                                                             vfm_objective, 
                                                             field_data.state, 
                                                             0)["virtual_power"]

matcal_internal_power_connected = results.best_simulation_qois(connected_hex_model, 
                                                             vfm_objective, 
                                                             field_data.state, 
                                                             0)["virtual_power"]
matcal_external_power_connected = results.get_experiment_qois(connected_hex_model, 
                                                             vfm_objective, 
                                                             field_data.state, 
                                                             0)["virtual_power"]

print(percent_error(matcal_external_power_default, L))
print(percent_error(matcal_external_power_connected, L))

#%%
# Both of the external powers have
# no error. This is expected as our 
# experimental data for this study was 
# exact. This result demonstrates in part that 
# all of the MatCal's VFM objective calculations and virtual 
# fields are correctly implemented. We have unit tests
# that test each of the components to contribute to 
# this individually. 
# The next evaluation compares the internal 
# virtual power to the expected virtual power. 
# This is significantly more involved as it requires 
# the simulation of the VFM model using SierraSM. 

internal_power_error_default = percent_error(matcal_internal_power_default.reshape(10,1), L)
internal_power_error_connected = percent_error(matcal_internal_power_connected.reshape(10,1), L)
print(internal_power_error_default)
print(internal_power_error_connected)

#%%
# When we print the percent errors, it is clear that the 
# some part of the process has introduced errors. The 
# magnitude of the maximum error is 0.05%. To
# further investigate, we plot the internal 
# power errors as a function time. 

import matplotlib.pyplot as plt

plt.plot(time, internal_power_error_default, label="default VFM model", marker='o')
plt.plot(time, internal_power_error_connected, label="connected hex VFM model")

plt.xlabel("time (s)")
plt.ylabel("Internal Power Error (%)")
plt.legend()
plt.show()

#%%
# The error is increasing quadratically as
# the load is increased and the model 
# becomes more deformed. This error is due to 
# our initial assumption of small deformation. 
# SierraSM is formulated as a large deformation code. 
# As a result, our small deformation assumption results 
# in this small, but noticeable error. In later verification 
# examples, we evaluate our methods with synthetic data 
# using large deformation solutions and are able to 
# obtain results with much smaller errors.