Skip to content
Snippets Groups Projects
Commit 18318bbe authored by Ruben Heinrich's avatar Ruben Heinrich
Browse files

[feature] WIP working on .vtu handling

parent c1ad501e
No related branches found
No related tags found
1 merge request!6Merge branch 7 fixed functionality into master
#===== read_in_pvd =====
#'read_in_pvd
#'@description Function to read in a .pvd file
#'@param pvd_path string: Path to .pvd file that should be read in
#'@return character: Vector containing .vtu paths
#'@export
read_in_pvd <- function(pvd_path) {
xml_doc <- validate_read_in_xml(pvd_path)
xpath_expr <- "/VTKFile"
root_node <- xml2::xml_find_first(xml_doc, xpath_expr)
dataset_nodes <- xml2::xml_find_all(root_node, "//Collection/DataSet")
vtu_refs <- character()
for(i in seq_len(length(dataset_nodes))){
vtu_refs <- c(vtu_refs, xml2::xml_attrs(dataset_nodes[[i]])[["file"]])
}
return(invisible(vtu_refs))
}
#===== read_in_vtu =====
#'read_in_vtu
#'@description Wrapper function to read in a whole .vtu file as a OGS6_vtu
#' class object
#'@param vtu_path string: Path to .vtu file that should be read in
#'@return OGS6_vtu: Mesh object
#'@export
read_in_vtu <- function(vtu_path) {
xml_doc <- validate_read_in_xml(vtu_path)
xpath_expr <- "/VTKFile"
root_node <- xml2::xml_find_first(xml_doc, xpath_expr)
vtu_obj <-
node_to_r2ogs6_class_object(xml_node = root_node,
xpath_expr = xpath_expr,
subclasses_names =
get_subclass_names("OGS6_vtu"))
return(invisible(vtu_obj))
}
#===== read_in_PointData_DataArray =====
#'read_in_PointData_DataArray
#'@description Wrapper function to read in a `PointData` `DataArray` element
#' from a .vtu file
#'@param vtu_path string: Path to .vtu file that should be read in
#'@param Name string: `Name` attribute of `DataArray` element
#'@export
read_in_PointData_DataArray <- function(vtu_path,
Name) {
assertthat::assert_that(assertthat::is.string(vtu_path))
assertthat::assert_that(assertthat::is.string(Name))
# load vtu
vtk_xml_ugr <- vtk$vtkXMLUnstructuredGridReader()
vtk_xml_ugr$SetFileName(vtu_path)
vtk_xml_ugr$Update()
# extract data
wrapped_data <- vtk_dsa$WrapDataObject(vtk_xml_ugr$GetOutput())
wrapped_data_arr <- wrapped_data$PointData[[Name]]
return(invisible(wrapped_data_arr))
}
#===== Decoding and decompressing functionality =====
#'decode_appended_data
#'@description Decodes AppendedData
#'@param appended_data character: `AppendedData` parameter of `OGS6_vtu`
#'@param data_arrays list: Content of lists specified in
#' `get_valid_vtu_categories()`
#'@param compressor string: Optional: How the data was compressed, this is the
#' `compressor` parameter of `OGS6_vtu`
#'@return list: DataArrays with `data` element which is the decoded data
decode_appended_data <- function(appended_data,
data_arrays,
compressor = "") {
assertthat::assert_that(is.character(appended_data))
assertthat::assert_that(length(appended_data) == 2)
assertthat::assert_that(assertthat::is.string(compressor))
encoding <- appended_data[["encoding"]]
appended_data <- substring(appended_data[["xml_text"]], 2)
for(i in seq_len(length(data_arrays))){
offset <- data_arrays[[i]][["offset"]]
next_offset <- ifelse(i < length(data_arrays),
data_arrays[[i+1]][["offset"]],
(nchar(appended_data) + 1))
encoded_data <- substring(appended_data, offset, (next_offset - 1))
data <- decode_data_array_data(encoded_data,
encoding = encoding,
compressor = compressor)
data_arrays[[i]] <- c(data_arrays[[i]], data = data)
}
return(invisible(data_arrays))
}
decode_data_array_data <- function(data_array_data,
encoding = "",
compressor = ""){
assertthat::assert_that(assertthat::is.string(data_array_data))
assertthat::assert_that(assertthat::is.string(encoding))
assertthat::assert_that(assertthat::is.string(compressor))
decoded_data <- data_array_data
# Decode
if(encoding == "base64"){
decoded_data <- base64enc::base64decode(decoded_data)
}else{
stop(paste("Encoding of AppendedData is not `base64`."), call. = FALSE)
}
switch(
compressor,
vtkZLibDataCompressor = {
py_env <- reticulate::py_run_string(
paste(
"import zlib",
"def decompress(x):",
"\t",
"\treturn zlib.decompress(x)",
sep = "\n"
),
convert = TRUE
)
decoded_data <- py_env$decompress(decoded_data)
}
)
return(invisible(decoded_data))
}
#===== get_vtu_data_arrays_from_file =====
#'get_vtu_data_arrays_from_file
#'@description Reads DataArray elements from a .vtu file
#'@param vtu_path string: .vtu file path
#'@param Names character: Optional: Select `DataArray` elements by `Name`
#' attribute
#'@param categories character: Optional: One or more of `FieldData`, `PointData`,
#' `CellData`, `Points` or `Cells`. If left empty, will get `DataArray`
#' elements from whole XML document.
get_vtu_data_arrays_from_file <- function(vtu_path,
Names = character(),
categories = character()){
assertthat::assert_that(assertthat::is.string(vtu_path))
assertthat::assert_that(is.character(Names))
assertthat::assert_that(is.character(categories))
valid_categories <- get_valid_vtu_categories()
lapply(categories, function(x){
assertthat::assert_that(x %in% valid_categories)
})
xml_doc <- validate_read_in_xml(vtu_path)
xpath_expr <- "/VTKFile/UnstructuredGrid"
data_array_nodes <- list()
if(length(categories) != 0){
for(i in seq_len(length(categories))){
categories_expr <- ifelse(categories[[i]] != "FieldData",
"/Piece",
"")
categories_expr <- paste0(xpath_expr,
categories_expr,
"/", categories[[i]], "/DataArray")
data_array_nodes <- c(data_array_nodes,
list(xml2::xml_find_all(xml_doc,
categories_expr)))
}
}else{
data_array_nodes <- xml2::xml_find_all(xml_doc, "//DataArray")
}
if(length(Names) != 0){
data_array_nodes <-
data_array_nodes[xml2::xml_attr(data_array_nodes,
"Name") %in% Names]
}
data_arrays <- lapply(data_array_nodes, function(x){
node_to_object(x)
})
return(invisible(data_arrays))
}
#===== OGS6_vtu ===== #===== OGS6_pvd =====
#'OGS6_vtu #'OGS6_pvd
#'@description Constructor for the OGS6_vtu base class #'@description Constructor for the OGS6_pvd base class
#'@export #'@export
OGS6_vtu <- R6::R6Class( OGS6_pvd <- R6::R6Class(
"OGS6_vtu", "OGS6_pvd",
public = list( public = list(
#'@description #'@description
#'Creates new OGS6_vtu object #'Creates new OGS6_pvd object
#'@param type string: #'@param pvd_path string: Path to .pvd file
#'@param version string: initialize = function(pvd_path) {
#'@param byte_order string:
#'@param UnstructuredGrid OGS6_UnstructuredGrid:
#'@param header_type string: Optional:
#'@param compressor string: Optional:
#'@param AppendedData string: Optional:
initialize = function(type,
version,
byte_order,
UnstructuredGrid,
header_type = NULL,
compressor = NULL,
AppendedData = NULL) {
self$type <- type
self$version <- version
self$byte_order <- byte_order
self$UnstructuredGrid <- UnstructuredGrid
if(!is.null(header_type)){
self$header_type <- header_type
}
if(!is.null(compressor)){ xml_doc <- validate_read_in_xml(pvd_path)
self$compressor <- compressor dataset_nodes <- xml2::xml_find_all(xml_doc,
} "/VTKFile/Collection/DataSet")
if(!is.null(AppendedData)){ private$.pvd_path <- pvd_path
self$AppendedData <- AppendedData private$.datasets <- lapply(dataset_nodes,
} node_to_object)
} private$.OGS6_vtus <- lapply(self$abs_vtu_paths,
), read_in_vtu)
active = list(
#'@field type
#'Access to private parameter '.type'
type = function(value) {
if(missing(value)) {
private$.type
}else{
assertthat::assert_that(assertthat::is.string(value))
private$.type <- value
}
}, },
#'@field version #'@description
#'Access to private parameter '.version' #'Returns .vtu path for specified timestep
version = function(value) { #'@param timestep string: Timestep
if(missing(value)) { get_vtu_path_by_timestep = function(timestep){
private$.version
}else{
assertthat::assert_that(assertthat::is.string(value))
private$.version <- value
}
},
#'@field byte_order assertthat::assert_that(assertthat::is.string(timestep))
#'Access to private parameter '.byte_order'
byte_order = function(value) {
if(missing(value)) {
private$.byte_order
}else{
assertthat::assert_that(assertthat::is.string(value))
private$.byte_order <- value
}
},
#'@field UnstructuredGrid for(i in seq_len(length(private$.datasets))){
#'Access to private parameter '.UnstructuredGrid' if(private$.datasets[[i]][["timestep"]] == timestep){
UnstructuredGrid = function(value) { return(private$.datasets[[i]][["file"]])
if(missing(value)) { }
private$.UnstructuredGrid
}else{
assertthat::assert_that("OGS6_UnstructuredGrid" %in%
class(value))
private$.UnstructuredGrid <- value
} }
},
#'@field header_type warning(paste("No .vtu path found for timestep", timestep),
#'Access to private parameter '.header_type' call. = FALSE)
header_type = function(value) {
if(missing(value)) {
private$.header_type
}else{
assertthat::assert_that(assertthat::is.string(value))
private$.header_type <- value
}
}, },
#'@field compressor #'@description
#'Access to private parameter '.compressor' #'Returns timestep for specified .vtu path
compressor = function(value) { #'@param vtu_path string: .vtu path
if(missing(value)) { get_timestep_by_vtu_path = function(vtu_path){
private$.compressor
}else{
assertthat::assert_that(assertthat::is.string(value))
private$.compressor <- value
}
},
#'@field AppendedData assertthat::assert_that(assertthat::is.string(vtu_path))
#'Access to private parameter '.AppendedData '
AppendedData = function(value) { for(i in seq_len(length(private$.datasets))){
if(missing(value)) { if(private$.datasets[[i]][["file"]] == vtu_path){
private$.AppendedData return(private$.datasets[[i]][["timestep"]])
}else{ }
assertthat::assert_that(assertthat::is.string(value))
private$.AppendedData <- value
} }
},
#'@field is_subclass warning(paste("No timestep found for .vtu path", vtu_path),
#'Access to private parameter '.is_subclass' call. = FALSE)
is_subclass = function() {
private$.is_subclass
}, },
#'@field attr_names #'@description
#'Access to private parameter '.attr_names' #'Creates a tibble object from PointData
attr_names = function() { #'@param Names character: `Name` attributes of `DataArray` elements
private$.attr_names get_PointData_time_tibble = function(Names){
},
#'@field flatten_on_exp assertthat::assert_that(is.character(Names))
#'Access to private parameter '.flatten_on_exp'
flatten_on_exp = function() {
private$.flatten_on_exp
}
),
private = list( time_list <- list()
.type = NULL,
.version = NULL,
.byte_order = NULL,
.UnstructuredGrid = NULL,
.header_type = NULL,
.compressor = NULL,
.AppendedData = NULL,
.is_subclass = FALSE,
.attr_names = c("type",
"version",
"byte_order",
"header_type",
"compressor"),
.flatten_on_exp = character()
)
)
# For each .vtu file referenced in pvd_path...
for(i in seq_len(length(self$OGS6_vtus))){
#===== OGS6_UnstructuredGrid ===== new_row <- list()
# ... get row of PointData by Name
for(j in seq_len(length(Names))){
point_data <-
self$OGS6_vtus[[i]]$get_PointData_DataArray(Names[[j]])
new_row <- c(new_row, list(list(point_data)))
names(new_row)[[length(new_row)]] <- Names[[j]]
}
#'OGS6_UnstructuredGrid time_list <- c(time_list,
#'@description Constructor for the OGS6_UnstructuredGrid base class list(tibble::as_tibble_row(new_row)))
#'@export }
OGS6_UnstructuredGrid <- R6::R6Class(
"OGS6_UnstructuredGrid", # Combine into tibble
public = list( time_tibble <- dplyr::bind_rows(time_list)
return(time_tibble)
},
get_PointData_timeline = function(point_id,
Name,
starting_from_timestep,
ending_on_timestep){
# ...
},
#'@description #'@description
#'Creates new OGS6_UnstructuredGrid object #'Gets PointData at specified timestep. Calls `get_PointData_timeline`
#'@param Piece OGS6_Piece: #' internally with `starting_from_timestep` and `ending_on_timestep`
#'@param FieldData character, length == 2: Optional: #' both being `timestep`
initialize = function(Piece, #'@param point_id number: Point ID
FieldData = NULL) { #'@param Name string: `Name` attribute of `DataArray` element
self$Piece <- Piece #'@param timestep string: Timestep
get_PointData_at_timestep = function(point_id,
if(!is.null(FieldData)){ Name,
self$FieldData <- FieldData timestep){
}
self$get_PointData_timeline(point_id = point_id,
Name = Name,
starting_from_timestep = timestep,
ending_on_timestep = timestep)
} }
), ),
active = list( active = list(
#'@field Piece
#'Access to private parameter '.Piece' #'@field pvd_path
Piece = function(value) { #'Getter for private parameter '.pvd_path'
if (missing(value)) { pvd_path = function() {
private$.Piece private$.pvd_path
}else{
assertthat::assert_that("OGS6_Piece" %in% class(value))
private$.Piece <- value
}
}, },
#'@field FieldData #'@field datasets
#'Access to private parameter '.FieldData' #'Getter for private parameter '.datasets'
FieldData = function(value) { datasets = function() {
if (missing(value)) { private$.datasets
private$.FieldData
} else{
private$.FieldData <- value
}
}, },
#'@field is_subclass #'@field vtu_paths
#'Access to private parameter '.is_subclass' #'Getter for `datasets` `file`
is_subclass = function() { #'@return character: .vtu paths as referenced in `pvd_path`, for
private$.is_subclass #' absolute paths use `abs_vtu_paths`
vtu_paths = function() {
vtu_paths <- lapply(private$.datasets, function(x){
x[["file"]]
})
}, },
#'@field attr_names #'@field abs_vtu_paths
#'Access to private parameter '.attr_names' #'Gets absolute .vtu paths, e.g. `dirname(pvd_path)` + `datasets` `file`
attr_names = function() { #'@return character: Absolute .vtu paths
private$.attr_names abs_vtu_paths = function() {
abs_vtu_paths <- lapply(self$vtu_paths, function(x){
abs_vtu_path <- paste0(dirname(self$pvd_path), "/", x)
return(invisible(abs_vtu_path))
})
}, },
#'@field flatten_on_exp #'@field timesteps
#'Access to private parameter '.flatten_on_exp' #'Gets timesteps from private parameter '.datasets'
flatten_on_exp = function() { timesteps = function() {
private$.flatten_on_exp
timesteps <- lapply(private$.datasets, function(x){
x[["timestep"]]
})
},
#'@field OGS6_vtus
#'Getter for private parameter '.OGS6_vtus'
OGS6_vtus = function() {
private$.OGS6_vtus
} }
), ),
private = list( private = list(
.Piece = NULL, .pvd_path = NULL,
.FieldData = NULL, .datasets = NULL,
.is_subclass = TRUE, .OGS6_vtus = NULL
.attr_names = character(),
.flatten_on_exp = character()
) )
) )
#===== OGS6_Piece ===== #===== OGS6_vtu =====
#'OGS6_Piece #'OGS6_vtu
#'@description Constructor for the OGS6_Piece base class #'@description Constructor for the OGS6_vtu base class
#'@export #'@export
OGS6_Piece <- R6::R6Class( OGS6_vtu <- R6::R6Class(
"OGS6_Piece", "OGS6_vtu",
public = list( public = list(
#'@description #'@description
#'Creates new OGS6_Piece object #'Creates new OGS6_vtu object
#'@param NumberOfPoints string | number: #'@param vtkUnstructuredGrid
#'@param NumberOfCells string | number: initialize = function(vtu_path,
#'@param PointData list, : vtkUnstructuredGrid) {
#'@param Points list, : self$vtkUnstructuredGrid <- vtkUnstructuredGrid
#'@param Cells list, :
#'@param CellData list, : Optional:
initialize = function(NumberOfPoints,
NumberOfCells,
PointData,
Points,
Cells,
CellData = NULL) {
self$NumberOfPoints <- NumberOfPoints
self$NumberOfCells <- NumberOfCells
self$PointData <- PointData
self$Points <- Points
self$Cells <- Cells
if(!is.null(CellData)){
self$CellData <- CellData
}
}
),
active = list(
#'@field NumberOfPoints
#'Access to private parameter '.NumberOfPoints'
NumberOfPoints = function(value) {
if (missing(value)) {
private$.NumberOfPoints
} else{
private$.NumberOfPoints <- value
}
},
#'@field NumberOfCells
#'Access to private parameter '.NumberOfCells'
NumberOfCells = function(value) {
if (missing(value)) {
private$.NumberOfCells
} else{
private$.NumberOfCells <- value
}
}, },
#'@field PointData get_PointData_DataArray = function(Name){
#'Access to private parameter '.PointData'
PointData = function(value) {
if (missing(value)) {
private$.PointData
} else{
private$.PointData <- value
}
},
#'@field Points if(missing(Name)){
#'Access to private parameter '.Points' return(self$dsa_wrapped_vtkUnstructuredGrid$PointData)
Points = function(value) {
if (missing(value)) {
private$.Points
} else{
private$.Points <- value
} }
},
#'@field Cells assertthat::assert_that(assertthat::is.string(Name))
#'Access to private parameter '.Cells' return(self$dsa_wrapped_vtkUnstructuredGrid$PointData[[Name]])
Cells = function(value) { }
if (missing(value)) {
private$.Cells
} else{
private$.Cells <- value
}
},
#'@field CellData ),
#'Access to private parameter '.CellData'
CellData = function(value) {
if (missing(value)) {
private$.CellData
} else{
private$.CellData <- value
}
},
#'@field is_subclass active = list(
#'Access to private parameter '.is_subclass'
is_subclass = function() {
private$.is_subclass
},
#'@field subclasses_names #'@field vtkUnstructuredGrid
#'Access to private parameter '.subclasses_names' #'Access to private parameter '.vtkUnstructuredGrid'
subclasses_names = function() { vtkUnstructuredGrid = function(value) {
private$.subclasses_names if(missing(value)) {
private$.vtkUnstructuredGrid
}else{
# Check class
private$.vtkUnstructuredGrid <- value
private$.dsa_wrapped_vtkUnstructuredGrid <-
vtk_dsa$WrapDataObject(value)
}
}, },
#'@field attr_names #'@field dsa_wrapped_vtkUnstructuredGrid
#'Access to private parameter '.attr_names' #'Getter for private parameter '.dsa_wrapped_vtkUnstructuredGrid'
attr_names = function() { dsa_wrapped_vtkUnstructuredGrid = function() {
private$.attr_names private$.dsa_wrapped_vtkUnstructuredGrid
} }
), ),
private = list( private = list(
.NumberOfPoints = NULL, .vtkUnstructuredGrid = NULL,
.NumberOfCells = NULL, .dsa_wrapped_vtkUnstructuredGrid = NULL
.PointData = NULL,
.Points = NULL,
.Cells = NULL,
.CellData = NULL,
.is_subclass = TRUE,
.attr_names = c("NumberOfPoints",
"NumberOfCells"),
.flatten_on_exp = character()
) )
) )
#'read_in_vtu
#'@description Reads in .vtu file via `vtkXMLUnstructuredGridReader` from the
#' python `vtk` library
#'@param vtu_path string: Path to .vtu file
#'@return vtkUnstructuredGrid*: Unstructured Grid
#'@export
read_in_vtu <- function(vtu_path) {
vtk_xml_ugr <- vtk$vtkXMLUnstructuredGridReader()
vtk_xml_ugr$SetFileName(vtu_path)
vtk_xml_ugr$Update()
return(invisible(OGS6_vtu$new(vtu_path,
vtk_xml_ugr$GetOutput())))
}
#===== generate_structured_mesh ===== #===== generate_structured_mesh =====
...@@ -378,39 +262,27 @@ OGS6_Piece <- R6::R6Class( ...@@ -378,39 +262,27 @@ OGS6_Piece <- R6::R6Class(
#'@description Wrapper function to call generateStructuredMesh.exe #'@description Wrapper function to call generateStructuredMesh.exe
#' (VTK mesh generator). For full documentation see #' (VTK mesh generator). For full documentation see
#'https://www.opengeosys.org/docs/tools/meshing/structured-mesh-generation/ #'https://www.opengeosys.org/docs/tools/meshing/structured-mesh-generation/
#'@param ogs6_obj OGS6: Simulation object #'@param args_str string: The arguments the script will be called with
#'@param call_str string: The arguments the script will be called with
#' (EXCEPT -o output_file_name, this will be generated automatically!)
#'@param read_in_vtu flag: Should .vtu file just be copied or read in too?
#'@return string: .vtu file path #'@return string: .vtu file path
#'@export #'@export
generate_structured_mesh = function(ogs6_obj, generate_structured_mesh = function(ogs_bin_path,
call_str, args_str) {
read_in_vtu) {
assertthat::assert_that(inherits(ogs6_obj, "OGS6"))
assertthat::assert_that(assertthat::is.string(call_str))
mesh_number <- 1
is_first <- (length(ogs6_obj$meshes) == 0) if(missing(ogs_bin_path)){
ogs_bin_path <- unlist(options("r2ogs6.default_ogs_bin_path"))
if(!is_first){
mesh_number <- length(ogs6_obj$meshes) + 1
} }
vtu_dir_path <- tempdir() assertthat::assert_that(assertthat::is.string(ogs_bin_path))
vtu_output_filename <- paste0(ogs6_obj$sim_name, "_", mesh_number, ".vtu") assertthat::assert_that(assertthat::is.string(args_str))
vtu_path <- paste0(vtu_dir_path, vtu_output_filename)
system(command = paste0(ogs6_obj$ogs_bin_path, "generateStructuredMesh.exe",
" -o ", vtu_path, " ", call_str))
if(read_in_vtu){ # Get .vtu path from args_str
read_in_vtu(ogs6_obj, vtu_path) vtu_path <- stringr::str_extract(args_str, "-o [^ ]*")
}else{ vtu_path <- stringr::str_remove(vtu_path, "-o ")
ogs6_obj$add_mesh(vtu_path)
} system(command = paste0(ogs_bin_path,
"generateStructuredMesh.exe ",
args_str))
return(invisible(vtu_path)) return(invisible(vtu_path))
} }
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment