R/DS_plot.R

In maxgav13/GMisc: Miscellaneous Functions for Geology (and its areas)

#' @title Draws the results for a direct shear test
#' @description Draws the sigma-tau pairs along with the Mohr-Coulomb failure envelope for a set of direct shear lab test results, and annotates the graph with the values of cohesion and friction angle.
#' @param sig.n A numeric vector of the normal stress values
#' @param tau A numeric vector of the shear stress values
#' @param units A string with the units of the measurements
#' @export
#' @return A ggplot of the direct shear test results
#' @import ggplot2
#' @references Coduto, D. P. (1999). Geotechnical Engineering - Principles and Practices. Prentice Hall.
#' @references Holtz, R. D., Kovacs, W. D. & Sheahan, T. C. (2011). An Introduction to Geotechnical Engineering. Prentice Hall.
#' @references Gonzalez de Vallejo, L. I. (2004). Ingenieria Geologica. Prentice Hall.
#' @examples
#' sig.n = c(80,237,395)
#' tau = c(127,345,475)
#' DS_plot(sig.n, tau)
#'
DS_plot <- function(sig.n, tau, units = 'kPa') {
  # Calculates failure envelope and c and phi parameters
  ds = stats::lm(tau ~ sig.n)
  DS = stats::coef(ds)
  if(DS[[1]] > 0){
    c = DS[[1]]
    phi = degs(atan(DS[[2]]))
  } else {
    DS = stats::coef(stats::lm(tau ~ sig.n - 1))
    c = 0
    phi = degs(atan(DS[[1]]))
  }
  r2 = summary(ds)$r.squared
  r = stats::cor(sig.n,tau)
  (rmse = sqrt(sum((stats::residuals(ds))^2)/stats::df.residual(ds)))

  # Calculates center and radius for each pair
  C = sig.n+(tau/tan(rads(90-phi)))
  R = tau/sin(rads(90-phi))

  # Calculates principal stresses
  sig1 = C+R
  sig3 = C-R

  # MC failure envelope
  x = seq(0,max(sig.n)*1.2)
  y = c+tan(rads(phi))*x

  # Mohr circle calculations
  circ = 0:180
  circ_x = matrix(0,length(circ),length(tau))
  circ_y = matrix(0,length(circ),length(tau))

  for (i in 1:length(tau)) {
    circ_x[,i] = C[i]+R[i]*cos(rads(circ))
  }

  for (i in 1:length(tau)) {
    circ_y[,i] = R[i]*sin(rads(circ))
  }

  Mohr = list()
  for (i in 1:length(tau)) {
    Mohr[[i]] = data.frame(X=circ_x[,i],Y=circ_y[,i],cir=i)
  }
  #MOHR = do.call(rbind,Mohr)
  MOHR = dplyr::bind_rows(Mohr,.id = 'circ')
  MOHR$cir = as.factor(MOHR$cir)

  # Text to annotate with phi and c parameters
  labelphi = sprintf("phi*minute == %0.1f*degree",phi)
  labelc = sprintf("c*minute == %0.1f~%s",c,units)

  # Plot failure envelope and sig.n-tau pairs
  gst = ggplot()+
    geom_point(aes(sig.n,tau),col="blue",size=2)+
    geom_line(aes(x,y),col="red",na.rm = T)+
    coord_fixed()+
    xlim(0,1.1*sig.n[length(sig.n)])+
    ylim(0,1.1*tau[length(tau)])+
    labs(x=paste0('Normal stress (',units,')'),
         y=paste0('Shear stress (',units,')')) +
    # labs(x=TeX(str_glue('$\\sigma$ ({units})')),
    #      y=TeX(str_glue('$\\tau$ ({units})')))+
    annotate("text",x=(max(sig.n)-min(sig.n))/2,y=max(tau),
             label=labelphi,parse=T)+
    annotate("text",x=(max(sig.n)-min(sig.n))/2,y=0.95*max(tau),
             label=labelc,parse=T)+
    theme_bw()
  return(gst)
}