FEFLOW.R

library(xts)
library(lubridate)

# Write time-series for FEFLOW
# userunit = m^3/d, l/min, m, ^0C
writeTS <- function(fPath, x,  
         option = c("linear", "cyclic"), 
         type = c("Polylined", "Constant"),
         unitclass = c("PUMP_RATE", "TEMPERATURE", "LENGTH", "VELOCITY", "PUMP_HEAT_RATE"),
         userunit,
         timeunit = "d",
         username = "NAME",
         TSID0 = 1,
         t0 = 0){
  
  type <- match.arg(type,  c("Polylined", "Constant"))
  option <- match.arg(option, c("linear", "cyclic"))
  unitclass <- match.arg(unitclass, c("PUMP_RATE", "TEMPERATURE", "LENGTH", "VELOCITY", "PUMP_HEAT_RATE"))
  
  # first clean the file if it exists
  dsn <- file(fPath, open = "wb")
  close(dsn)
  
  # switch in appending mode
  dsn <- file(fPath, open = "ab")
  
  if(is.list(x)){
    n <- length(x)
  }else if(is.zoo(x)){
    n <- ncol(x)
  }else if(is.matrix(x)){
    n <- ncol(x) - 1
  }else{
    stop("lkjlkjklj")
  }
  if(length(username) != n) username <- paste(username, 1:n)
  if(length(TSID0) != n) TSID0 <- TSID0[1] + (1:n) - 1
  for(i in 1:n){
    writeChar(paste0("# ", TSID0[i], "\r\n"), con = dsn, eos = NULL)
    writeChar(paste0("! ", username[i], "\r\n"),  con = dsn, eos = NULL)
    writeChar(paste0("! [type=", 
                     type,
                     ";option=", 
                     option,
                     ";timeunit=",
                     timeunit,
                     ";unitclass=",
                     unitclass,
                     ";userunit=",
                     userunit,
                     "]\r\n"), con = dsn, eos = NULL)
    if(is.list(x)){
      if(is.zoo(x[[i]])){
        ttt <- (as.numeric(time(x[[i]])) - as.numeric(t0))/3600/24
        XX <- paste0(as.character(ttt), " ", as.character(x[[i]]), "\r\n")
      }else if(is.matrix(x[[i]])){
        XX <- paste(as.character(x[[i]][,1] - as.numeric(t0)), as.character(x[[i]][,2]), "\r\n")
      }else{
        stop("lkjlkjlj")
      }
    }else if(is.zoo(x)){
      ttt <- (as.numeric(time(x)) - as.numeric(t0))/3600/24
      XX <- paste(as.character(ttt), as.character(x[,i]), "\r\n")
    }else if(is.matrix(x)){
      XX <- paste(as.character(x[,1] - as.numeric(t0)), as.character(x[,i+1]), "\r\n")
    }else{
      stop("lkjlkjklj")
    }
    writeChar(XX, con = dsn, eos = NULL)
    writeChar("END\r\n",  con = dsn, eos = NULL)
  }
  writeChar("END\r\n",  con = dsn, eos = NULL)
  
  close(dsn)  
  
}



#' Write mesh nodes at an ideal distance from a well
#' 
#' @param xy [\code{numeric(2)}] Well coordinates
#' @param r [\code{numeric(1)}] Well radius
#' @param n [\code{integer}] Number of points
#' @param fPath [\code{character(1)}] File names for exporting the points
# example
# xy <- c(651645, 225560)
# r <- 0.108
# n <- 6
# 
# xy_nodes <- wellMeshNodes(xy = c(651645, 225560), r = 0.108, n = 6)
# 
# plot(xy_nodes, col = "blue", asp = 1)
# points(t(xy), col = "red", pch = 20)
# 
# writeWellMeshNodes("test.pnt", xy = c(651645, 225560), r = 0.108, n = 6)

wellMeshNodes <- function(xy = c(0,0), r = 1 , n = 6){
  dist_factor <- exp(2*pi / (n * tan(pi/n)))
  node_dist <- dist_factor * r
  
  theta <- ((1:n) - 1 ) * 2 * pi / n
  
  xy_nodes <- matrix(nrow = n, ncol = 2)
  xy_nodes[, 1] <- xy[1] + node_dist * cos(theta)  
  xy_nodes[, 2] <- xy[2] + node_dist * sin(theta)  
  return(xy_nodes)
}

# @param x sf feature with attribute DN (Durchmesser in meters!!)
# Brunnen Durchmesser (m) angeben falls nicht vorhanden
# wells$DN <- 0.125 # m
# # oder
# wells$DN <- c(0.125, 0.2, 0.2, 1)

# sfWellMeshNodes <- function(x, n = 6){
#  xy <- sf::st_coordinates(x)
#  wnodes <- apply(xy, 1, wellMeshNodes, r = x$DN[1]/2, n = n)
#  wells_nodes <- cbind(as.vector(wnodes[1:n, ]), as.vector(wnodes[(1:n) + n, ]))
#  well_nodes <- sf::st_as_sf(as.data.frame(wells_nodes),  
#                             coords = 1:2, 
#                             crs = sf::st_crs(x))
#  return(well_nodes)
#}

sfWellMeshNodes <- function(x, n = 6){
  xy <- sf::st_coordinates(x)
  wells_nodes <- matrix(NA, nrow = n * nrow(xy), ncol = 2)
  for(i in 1:nrow(xy)){
    wells_nodes[1:n + (i - 1)*n,] <- wellMeshNodes(xy[i, 1:2], r = x$DN[i]/2, n = n)
  }
  well_nodes <- sf::st_as_sf(as.data.frame(wells_nodes),  
                             coords = 1:2, 
                             crs = sf::st_crs(x))
  return(well_nodes)
}


writeWellMeshNodes <- function(fPath, xy = c(0,0), r = 1 , n = 6){
  xy_nodes <- wellMeshNodes(xy = xy, r = r, n = n)
  cat(paste(1:n, xy_nodes[, 1], xy_nodes[, 2], 
            paste("\"nodes", 1:n, "\""), 
            collapse = "\n", sep = "\t"), 
      file = fPath, append = FALSE ,sep="\n")
  cat("END", file = fPath, append = TRUE, sep="\n" )
}

#-----------------------------------------------------------------

# get the extremum values (either min or max)
getExt <- function(x){
  x[which.max(abs(x))[1]]
}

createXTS <- function(x, d0){
  xts::xts(x$T, d0 + x$t*60*60*24)
}

# @param X -> read with 'read.table()'
# @param d0 -> first day of simulation (generally 1st of January)
#              d0 <- lubridate::ymd_hms("2020-10-01 00:00:00")
# reformat the time-series exported from FEFLOW
feflowHistToXts <- function(X, d0){
  colnames(X) <- c("ID", "t", "T", "ID2")
  
  u <- rle(X$ID)
  u$values <- seq_along(u$values)
  
  ID2 <- inverse.rle(u)

  X <- split(X, ID2)
  XS <- lapply(X, createXTS, d0 = d0)
  XS2 <- na.approx(do.call(merge.xts, XS))
  TP_ID2 <- sapply(X, function(x) x$ID2[1])
  colnames(XS2) <- TP_ID2
  
  return(XS2)
}