# ASME V&V 10.1- 2012
#
# An Illustration of the Concepts of Verification and Validation in
# Computational Solid Mechanics
#
# ASME The American Society of Mechanical Engineers
# April 16, 2012

# Walter Ponge-Ferreira
# Cotia, 16.06.2018

setwd("~/Documents/Data/Pasta Tecnica/R/ASME V&V 10-1")

# Validation Approach 1

# Model Results
# X~N(m_mod,s_mod)
m_mod <- -14.2  # mean
s_mod <- 0.71/3 # sd

# Experimental Data
# Y~N(m_exp,s_exp)
m_exp <- -15.0  # mean
s_exp <- 0.75/3 # sd

# PDF
x <- seq(-16.5,-13,0.01)
p_x <- dnorm(x,mean = m_mod, sd = s_mod)
p_y <- dnorm(x,mean = m_exp, sd = s_exp)

plot(x,p_x, type = 'l', col = "blue", lwd = 1.5, xlab = "x / mm", ylab = "p(X)")
abline(v=m_mod+s_mod*c(-3,0,3), lty = c(2,1,2), col = "blue")
lines(x,p_y, type = 'l', col = "red", lwd = 1.5)
abline(v=m_exp+s_exp*c(-3,0,3), lty = c(2,1,2), col = "red")
legend("topleft",c("Model", "Experimental"), col = c("blue","red"), pch = c(15,15))
title("Validation - SRQ Comparison")
text(-16.5, 1.2,pos = 4, expression(hat(mu) == -14.2), col = "blue")
text(-16.5, 1.0,pos = 4, expression(hat(sigma) == 0.24), col = "blue")
text(-16.5, 0.8,pos = 4, expression(hat(mu) == -15.0), col = "red")
text(-16.5, 0.6,pos = 4, expression(hat(sigma) == 0.25), col = "red")

# CDF
P_x <- pnorm(x,mean = m_mod, sd = s_mod)
P_y <- pnorm(x,mean = m_exp, sd = s_exp)

plot(x,P_x, type = 'l', col = "blue", lwd = 1.5, xlab = "x / mm", ylab = "P( X < x )")
lines(x,P_y, type = 'l', col = "red", lwd = 1.5)
title("Area Metric")
legend("topleft",c("Model", "Experimental"), col = c("blue","red"), pch = c(15,15))

# Normalized Area Metric M
dx <- mean(diff(x))
(A <- sum(abs(P_x-P_y))*dx)
(M <- A/abs(m_exp))
plot(x,abs(P_x-P_y), type = 'l', main = "CPF Difference", xlab = "x / mm", ylab = "| Pmod(X) - Pexp(X) |")
text(-16.5,0.8,pos = 4, paste("M = ",round(M*100, digits=1),"%"))
text(-16.5,0.7,pos = 4, paste("A = ",round(A, digits=2)," mm"))
text(-16.5,0.6,pos = 4, paste("Delta = ",round(abs(m_mod-m_exp),digits = 2)," mm"))

# Validation Approach 2

# Measured Beam-Tip Deflections from the Validation Experiments
w <- c(-16.3, -15.5, -16.1, -14.8, -14.4, -15.0, -15.3, -15.4, -15.6, -15.2)
(w_m <- mean(w))
(s_m <- sd(w))

# Empirical Cumulative Distribution Function
Fn <- ecdf(w)
plot(Fn, xlim = c(-16.5,-13.5), xlab = "x / mm")
lines(x,P_x, lty = 1, col = "blue")
lines(x,P_y, lty = 2, col = "red")
grid()

# Test Measurements of the Modulus of Elasticity, E, GPa
E <- c(69.1, 68.8, 74.4, 72.6, 72.9, 67.5, 74.1, 68.3, 71.0, 63.2)
hist(E, main = "Material Modulus", xlab = "E / Gpa")
text(72,3.0, pos = 4, paste("mean = ",round(mean(E),digits = 1), " GPa"))
text(72,2.8, pos = 4, paste("sd = ",round(sd(E),digits = 1), " GPa"))

# Test Estimates of the Support Flexibility, f, rad/Nm x 10e-7
f <- c(8.5, 8.2, 8.1, 8.4, 7.8, 7.7, 8.8, 8.6, 8.1, 8.6, 8.8, 8.5, 8.0, 9.3, 8.1, 9.3, 8.1, 8.3, 8.3)
hist(f, main = "Support Flexibility", xlab = "f / rad/Nm x 10e-7")
text(8.5,5.0, pos = 4, paste("mean = ",round(mean(f),digits=1), " rad/Nm x 10e-7"))
text(8.5,4.5, pos = 4, paste("sd = ",round(sd(f),digits=1), " rad/Nm x 10e-7"))

# The End