Nothing
R/mesh.R
In fdaPDE: Physics-Informed Spatial and Functional Data Analysis
Defines functions
refine.mesh.1.5D
projection.points.1.5D
create.mesh.1.5D
create.mesh.3D
projection.points.2.5D
create.mesh.2.5D
refine.mesh.2D
create.mesh.2D
triangulate_native
Documented in
create.mesh.1.5D
create.mesh.2.5D
create.mesh.2D
create.mesh.3D
projection.points.1.5D
projection.points.2.5D
refine.mesh.1.5D
refine.mesh.2D
triangulate_native <- function(P, PB, PA, S, SB,H, TR, flags) {
## It is necessary to check for NAs and NaNs, as the triangulate C
## code crashes if fed with them
P <- as.matrix(P)
PB <- as.integer(PB)
PA <- as.matrix(PA)
S <- as.matrix(S)
SB <- as.integer(SB)
H <- as.matrix(H)
TR <- as.matrix(TR)
storage.mode(P) <- "double"
storage.mode(PA) <- "double"
#storage.mode(PB) <- "integer"
storage.mode(S) <- "integer"
#storage.mode(SB) <- "integer"
storage.mode(H) <- "double"
storage.mode(TR) <- "integer"
storage.mode(flags) <- 'character'
## Call the main routine
out <- .Call("R_triangulate_native",
t(P),
PB,
PA,
t(S),
SB,
H,
t(TR),
flags,
PACKAGE="fdaPDE")
names(out) <- c("P", "PB", "PA", "T", "S", "SB", "E", "EB","TN", "VP", "VE", "VN", "VA")
class(out) <- "triangulation"
return(out)
}
#' Create a 2D triangular mesh
#'
#' @param nodes A #nodes-by-2 matrix containing the x and y coordinates of the mesh nodes.
#' @param nodesattributes A matrix with #nodes rows containing nodes' attributes.
#' These are passed unchanged to the output. If a node is added during the triangulation process or mesh refinement, its attributes are computed
#' by linear interpolation using the attributes of neighboring nodes. This functionality is for instance used to compute the value
#' of a Dirichlet boundary condition at boundary nodes added during the triangulation process.
#' @param segments A #segments-by-2 matrix. Each row contains the row's indices in \code{nodes} of the vertices where the segment starts from and ends to.
#' Segments are edges that are not splitted during the triangulation process. These are for instance used to define the boundaries
#' of the domain. If this is input is NULL, it generates a triangulation over the
#' convex hull of the points specified in \code{nodes}.
#' @param holes A #holes-by-2 matrix containing the x and y coordinates of a point internal to each hole of the mesh. These points are used to carve holes
#' in the triangulation, when the domain has holes.
#' @param triangles A #triangles-by-3 (when \code{order} = 1) or #triangles-by-6 (when \code{order} = 2) matrix.
#' This option is used when a triangulation is already available. It specifies the triangles giving the row's indices in \code{nodes} of the triangles' vertices and (when \code{nodes} = 2) also if the triangles' edges midpoints. The triangles' vertices and midpoints are ordered as described
#' at \cr https://www.cs.cmu.edu/~quake/triangle.highorder.html.
#' In this case the function \code{create.mesh.2D} is used to produce a complete mesh.2D object.
#' @param order Either '1' or '2'. It specifies wether each mesh triangle should be represented by 3 nodes (the triangle' vertices) or by 6 nodes (the triangle's vertices and midpoints).
#' These are
#' respectively used for linear (order = 1) and quadratic (order = 2) Finite Elements. Default is \code{order} = 1.
#' @param verbosity This can be '0', '1' or '2'. It indicates the level of verbosity in the triangulation process. When \code{verbosity} = 0 no message is returned
#' during the triangulation. When \code{verbosity} = 2 the triangulation process is described step by step by displayed messages.
#' Default is \code{verbosity} = 0.
#' @description This function is a wrapper of the Triangle library (http://www.cs.cmu.edu/~quake/triangle.html). It can be used
#' to create a triangulation of the domain of interest starting from a list of points, to be used as triangles' vertices, and a list of segments, that define the domain boundary. The resulting
#' mesh is a Constrained Delaunay triangulation. This is constructed in a way to preserve segments provided in the input \code{segments} without splitting them. This imput can be used to define the boundaries
#' of the domain. If this imput is NULL, it generates a triangulation over the
#' convex hull of the points.
#' It is also possible to create a mesh.2D from the nodes locations and the connectivity matrix.
#' @usage create.mesh.2D(nodes, nodesattributes = NA, segments = NA, holes = NA,
#' triangles = NA, order = 1, verbosity = 0)
#' @seealso \code{\link{refine.mesh.2D}}, \code{\link{create.FEM.basis}}
#' @return An object of the class mesh.2D with the following output:
#' \itemize{
#' \item{\code{nodes}}{A #nodes-by-2 matrix containing the x and y coordinates of the mesh nodes.}
#' \item{\code{nodesmarkers}}{A vector of length #nodes, with entries either '1' or '0'. An entry '1' indicates that the corresponding node is a boundary node; an entry '0' indicates that the corresponding node is not a boundary node.}
#' \item{\code{nodesattributes}}{A matrix with #nodes rows containing nodes' attributes.
#' These are passed unchanged from the input.}
#' \item{\code{triangles}}{A #triangles-by-3 (when \code{order} = 1) or #triangles-by-6 (when \code{order} = 2) matrix.
#' This option is used when a triangulation is already available. It specifies the triangles giving the indices in \code{nodes} of the triangles' vertices and (when \code{nodes} = 2) also if the triangles' edges midpoints. The triangles' vertices and midpoints are ordered as described
#' at \cr https://www.cs.cmu.edu/~quake/triangle.highorder.html.}
#' \item{\code{segmentsmarker}}{A vector of length #segments with entries either '1' or '0'. An entry '1' indicates that the corresponding element in \code{segments} is a boundary segment;
#' an entry '0' indicates that the corresponding segment is not a boundary segment.}
#' \item{\code{edges}}{A #edges-by-2 matrix containing all the edges of the triangles in the output triangulation. Each row contains the row's indices in \code{nodes}, indicating the nodes where the edge starts from and ends to.}
#' \item{\code{edgesmarkers}}{A vector of lenght #edges with entries either '1' or '0'. An entry '1' indicates that the corresponding element in \code{edge} is a boundary edge;
#' an entry '0' indicates that the corresponding edge is not a boundary edge.}
#' \item{\code{neighbors}}{A #triangles-by-3 matrix. Each row contains the indices of the three neighbouring triangles. An entry '-1' indicates that
#' one edge of the triangle is a boundary edge.}
#' \item{\code{holes}}{A #holes-by-2 matrix containing the x and y coordinates of a point internal to each hole of the mesh. These points are used to carve holes
#' in the triangulation, when the domain has holes.}
#' \item{\code{order}}{Either '1' or '2'. It specifies wether each mesh triangle should be represented by 3 nodes (the triangle' vertices) or by 6 nodes (the triangle's vertices and midpoints).
#' These are respectively used for linear (order = 1) and quadratic (order = 2) Finite Elements.}
#' }
#' @export
#' @examples
#' library(fdaPDE)
#'
#' ## Upload the quasicirle2D data
#' data(quasicircle2D)
#' boundary_nodes = quasicircle2D$boundary_nodes
#' boundary_segments = quasicircle2D$boundary_segments
#' locations = quasicircle2D$locations
#' data = quasicircle2D$data
#'
#' ## Create mesh from boundary
#' ## if the domain is convex it is sufficient to call:
#' mesh = create.mesh.2D(nodes = rbind(boundary_nodes, locations))
#' plot(mesh)
#'
#' ## if the domain is not convex, pass in addition the segments the compose the boundary:
#' mesh = create.mesh.2D(nodes = rbind(boundary_nodes, locations), segments = boundary_segments)
#'
#' ## Create mesh from data locations (without knowing the boundary)
#' mesh = create.mesh.2D(nodes = locations)
#' plot(mesh)
#' ## In this case the domain is the convex hull of the data locations.
#' ## Do this only if you do not have any information about the shape of the domain of interest.
create.mesh.2D <- function(nodes, nodesattributes = NA, segments = NA, holes = NA, triangles = NA, order = 1, verbosity = 0)
{
##########################
### Input checking ###
##########################
# Triangle finds out which are on the border (see https://www.cs.cmu.edu/~quake/triangle.help.html)
nodesmarkers = vector(mode = "integer", 0)
segmentsmarkers = vector(mode = "integer", 0)
nodes = as.matrix(nodes)
if (ncol(nodes) != 2)
stop("Matrix of nodes should have 2 columns")
if (anyDuplicated(nodes))
stop("Duplicated nodes")
## If attributes not specified, set them to a matrix with zero columns
if (any(is.na(nodesattributes))) {
nodesattributes <- matrix(0, nrow(nodes), 0)
}else{
nodesattributes <- as.matrix(nodesattributes)
if (nrow(nodesattributes) != nrow(nodes))
stop("Point attribute matrix \'nodesattributes\' does not have same number of rows the point matrix \'nodes\'")
}
## If boundary nodes not specified, set them to 0
# if (any(is.na(nodesmarkers))) {
# nodesmarkers <- vector(mode = "integer", 0)
# }else{
# nodesmarkers = as.vector(nodesmarkers)
# }
## Deal with segments
if (any(is.na(segments))) {
segments <- matrix(0, 0, 2)
} else {
segments <- as.matrix(segments)
if (ncol(segments) != 2) {
stop("Matrix of segments should have 2 columns")
}
}
## If boundary segments not specified, set them to 0
# if (any(is.na(segmentsmarkers))) {
# segmentsmarkers <- vector(mode = "integer", 0)
# }else{
# segmentsmarkers = as.vector(segmentsmarkers)
# }
## If hole not specified, set it to empty matrix
if (any(is.na(holes)))
holes <- matrix(0, 0, 2)
holes = as.matrix(holes)
## If triangles are not already specified
if(any(is.na(triangles)))
triangles = matrix(0,nrow = 0, ncol = 3)
triangles = as.matrix(triangles)
## Set meshing parameters ##
flags="ven"
if(nrow(segments) == 0){
flags = paste(flags,"c",sep = '')
}
if(nrow(segments)>0){
flags = paste(flags,"p",sep = '')
}
#If order=2 add flag for second order nodes
if(order == 2){
flags = paste(flags,"o2",sep = '')
}
if(order < 1 || order >2){
stop('Order must be 1 or 2')
}
if(nrow(triangles) > 0){
flags = paste(flags,"r",sep = '')
}
if (verbosity == 0) {
flags = paste(flags,"Q",sep = '')
}
if (verbosity == 1) {
flags = paste(flags,"V",sep = '')
}
if (verbosity == 2) {
flags = paste(flags,"VV",sep = '')
}
out<-NULL
#If triangles is null it makes the trianglulation
#If triangle is not null it makes a refinement with no parameter, to compose the mesh object
out <- triangulate_native(
nodes,
nodesmarkers,
nodesattributes,
segments,
segmentsmarkers,
t(holes),
triangles,
flags
)
names(out)[1]<-"nodes"
names(out)[2]<-"nodesmarkers"
names(out)[3]<-"nodesattributes"
names(out)[4]<-"triangles"
names(out)[5]<-"segments"
names(out)[6]<-"segmentsmarkers"
names(out)[7]<-"edges"
names(out)[8]<-"edgesmarkers"
names(out)[9]<-"neighbors"
out[13]<-NULL
out[12]<-NULL
out[11]<-NULL
out[10]<-NULL
out[[10]] = holes
names(out)[10]<-"holes"
out[[11]] = order
names(out)[11]<-"order"
class(out)<-"mesh.2D"
return(out)
}
#' Refine a 2D triangular mesh
#'
#' @param mesh A mesh.2D object representing the triangular mesh, created by \link{create.mesh.2D}.
#' @param minimum_angle A scalar specifying a minimun value for the triangles angles.
#' @param maximum_area A scalar specifying a maximum value for the triangles areas.
#' @param delaunay A boolean parameter indicating whether or not the output mesh should satisfy the Delaunay condition.
#' @param verbosity This can be '0', '1' or '2'. It indicates the level of verbosity in the triangulation process.
#' @description This function refines a Constrained Delaunay triangulation into a Conforming Delaunay triangulation. This is a wrapper of the Triangle library (http://www.cs.cmu.edu/~quake/triangle.html). It can be used to
#' refine a mesh previously created with \link{create.mesh.2D}. The algorithm can add Steiner points (points through which the \code{segments} are splitted)
#' in order to meet the imposed refinement conditions.
#' @usage refine.mesh.2D(mesh, minimum_angle, maximum_area, delaunay, verbosity)
#' @seealso \code{\link{create.mesh.2D}}, \code{\link{create.FEM.basis}}
#' @return A mesh.2D object representing the refined triangular mesh, with the following output:
#' \itemize{
#' \item{\code{nodes}}{A #nodes-by-2 matrix containing the x and y coordinates of the mesh nodes.}
#' \item{\code{nodesmarkers}}{A vector of length #nodes, with entries either '1' or '0'. An entry '1' indicates that the corresponding node is a boundary node; an entry '0' indicates that the corresponding node is not a boundary node.}
#' \item{\code{nodesattributes}}{nodesattributes A matrix with #nodes rows containing nodes' attributes.
#' These are passed unchanged to the output. If a node is added during the triangulation process or mesh refinement, its attributes are computed
#' by linear interpolation using the attributes of neighboring nodes. This functionality is for instance used to compute the value
#' of a Dirichlet boundary condition at boundary nodes added during the triangulation process.}
#' \item{\code{triangles}}{A #triangles-by-3 (when \code{order} = 1) or #triangles-by-6 (when \code{order} = 2) matrix.}
#' \item{\code{edges}}{A #edges-by-2 matrix. Each row contains the row's indices of the nodes where the edge starts from and ends to.}
#' \item{\code{edgesmarkers}}{A vector of lenght #edges with entries either '1' or '0'. An entry '1' indicates that the corresponding element in \code{edge} is a boundary edge;
#' an entry '0' indicates that the corresponding edge is not a boundary edge.}
#' \item{\code{neighbors}}{A #triangles-by-3 matrix. Each row contains the indices of the three neighbouring triangles. An entry '-1' indicates that
#' one edge of the triangle is a boundary edge.}
#' \item{\code{holes}}{A #holes-by-2 matrix containing the x and y coordinates of a point internal to each hole of the mesh. These points are used to carve holes
#' in the triangulation, when the domain has holes.}
#' \item{\code{order}}{Either '1' or '2'. It specifies wether each mesh triangle should be represented by 3 nodes (the triangle' vertices) or by 6 nodes (the triangle's vertices and midpoints).
#' These are respectively used for linear (order = 1) and quadratic (order = 2) Finite Elements.}
#' }
#' @examples
#' library(fdaPDE)
#'
#' ## Upload the quasicircle2D data
#' data(quasicircle2D)
#' boundary_nodes = quasicircle2D$boundary_nodes
#' boundary_segments = quasicircle2D$boundary_segments
#' locations = quasicircle2D$locations
#' data = quasicircle2D$data
#'
#' ## Create mesh from boundary:
#' mesh = create.mesh.2D(nodes = boundary_nodes, segments = boundary_segments)
#' plot(mesh)
#' ## Refine the mesh with the maximum area criterion:
#' finemesh = refine.mesh.2D(mesh = mesh, maximum_area = 0.1)
#' plot(finemesh)
#' ## Refine the mesh with the minimum angle criterion:
#' finemesh2 = refine.mesh.2D(mesh = mesh, minimum_angle = 30)
#' plot(finemesh2)
#' @export
refine.mesh.2D<-function(mesh, minimum_angle = NA, maximum_area = NA, delaunay = FALSE, verbosity = 0)
{
if(!is(mesh, "mesh.2D"))
stop("Sorry, this function is implemented just for mesh.2D class ")
flags="rpven"
if(!is.na(minimum_angle)){
flags <- paste(flags, "q", sprintf("%.12f", minimum_angle), sep='')
}
if(!is.na(maximum_area)){
flags <- paste(flags, "a", sprintf("%.12f", maximum_area), sep='')
}
if(delaunay){
flags <- paste(flags, "D", sep='')
}
if(mesh$order==2){
flags <- paste(flags, "o2", sep='')
}
if (verbosity == 0) {
flags = paste(flags,"Q",sep = '')
}
if (verbosity == 1) {
flags = paste(flags,"V",sep = '')
}
if (verbosity == 2) {
flags = paste(flags,"VV",sep = '')
}
# Triangle finds out which are on the border (see https://www.cs.cmu.edu/~quake/triangle.help.html)
mesh$nodesmarkers = vector(mode = "integer", 0)
mesh$segmentsmarkers = vector(mode = "integer", 0)
out<-NULL
#If triangles is null it makes the trianglulation
#If triangle is not null it makes a refinement with no parameter, to compose the mesh object
out <- triangulate_native(
mesh$nodes,
mesh$nodesmarkers,
mesh$nodesattributes,
mesh$segments,
mesh$segmentmarkers,
t(mesh$holes),
mesh$triangles,
flags
)
names(out)[1]<-"nodes"
names(out)[2]<-"nodesmarkers"
names(out)[3]<-"nodesattributes"
names(out)[4]<-"triangles"
names(out)[5]<-"segments"
names(out)[6]<-"segmentsmarkers"
names(out)[7]<-"edges"
names(out)[8]<-"edgesmarkers"
names(out)[9]<-"neighbors"
out[13]<-NULL
out[12]<-NULL
out[11]<-NULL
out[10]<-NULL
out[[10]] = mesh$holes
names(out)[10]<-"holes"
out[[11]] = mesh$order
names(out)[11]<-"order"
class(out)<-"mesh.2D"
return(out)
}
#' Create a \code{mesh.2.5D} object from the nodes locations and the connectivity matrix
#'
#' @param nodes A #nodes-by-3 matrix containing the x, y, z coordinates of the mesh nodes.
#' @param triangles A #triangles-by-3 (when \code{order} = 1) or #triangles-by-6 (when \code{order} = 2) matrix.
#' It specifies the triangles giving the row's indices in \code{nodes} of the triangles' vertices and (when \code{nodes} = 2) also if the triangles' edges midpoints. The triangles' vertices and midpoints are ordered as described
#' at \cr https://www.cs.cmu.edu/~quake/triangle.highorder.html.
#' @param order Either '1' or '2'. It specifies wether each mesh triangle should be represented by 3 nodes (the triangle' vertices) or by 6 nodes (the triangle's vertices and midpoints).
#' These are
#' respectively used for linear (order = 1) and quadratic (order = 2) Finite Elements. Default is \code{order} = 1.
#' @param nodesattributes A matrix with #nodes rows containing nodes' attributes.
#' These are passed unchanged to the output. This has been added for consistency with the function \code{create.mesh.2D}.
#' @param segments A #segments-by-2 matrix. Each row contains the row's indices in \code{nodes} of the vertices where the segment starts from and ends to.
#' Segments are edges that are not splitted during the triangulation process. These are for instance used to define the boundaries
#' of the domain. This has been added for consistency with the function \code{create.mesh.2D}.
#' @param holes A #holes-by-3 matrix containing the x, y, z coordinates of a point internal to each hole of the mesh. These points are used to carve holes
#' in the triangulation, when the domain has holes. This has been added for consistency with the function \code{create.mesh.2D}.
#' @return An object of the class mesh.2.5D with the following output:
#' \itemize{
#' \item{\code{nodes}}{A #nodes-by-3 matrix containing the x, y, z coordinates of the mesh nodes.}
#' \item{\code{nodesmarkers}}{A vector of length #nodes, with entries either '1' or '0'. An entry '1' indicates that the corresponding node is a boundary node; an entry '0' indicates that the corresponding node is not a boundary node.}
#' \item{\code{nodesattributes}}{A matrix with #nodes rows containing nodes' attributes.
#' These are passed unchanged from the input.}
#' \item{\code{triangles}}{A #triangles-by-3 (when \code{order} = 1) or #triangles-by-6 (when \code{order} = 2) matrix.
#' It specifies the triangles giving the indices in \code{nodes} of the triangles' vertices and (when \code{nodes} = 2) also if the triangles' edges midpoints. The triangles' vertices and midpoints are ordered as described
#' at \cr https://www.cs.cmu.edu/~quake/triangle.highorder.html.}
#' \item{\code{segmentsmarker}}{A vector of length #segments with entries either '1' or '0'. An entry '1' indicates that the corresponding element in \code{segments} is a boundary segment;
#' an entry '0' indicates that the corresponding segment is not a boundary segment.}
#' \item{\code{edges}}{A #edges-by-2 matrix containing all the edges of the triangles in the output triangulation. Each row contains the row's indices in \code{nodes}, indicating the nodes where the edge starts from and ends to.}
#' \item{\code{edgesmarkers}}{A vector of lenght #edges with entries either '1' or '0'. An entry '1' indicates that the corresponding element in \code{edge} is a boundary edge;
#' an entry '0' indicates that the corresponding edge is not a boundary edge.}
#' \item{\code{neighbors}}{A #triangles-by-3 matrix. Each row contains the indices of the three neighbouring triangles. An entry '-1' indicates that
#' one edge of the triangle is a boundary edge.}
#' \item{\code{holes}}{A #holes-by-3 matrix containing the x, y, z coordinates of a point internal to each hole of the mesh. These points are used to carve holes
#' in the triangulation, when the domain has holes. These are passed unchanged from the input.}
#' \item{\code{order}}{Either '1' or '2'. It specifies wether each mesh triangle should be represented by 3 nodes (the triangle' vertices) or by 6 nodes (the triangle's vertices and midpoints).
#' These are respectively used for linear (order = 1) and quadratic (order = 2) Finite Elements.}
#' }
#' @export
#' @examples
#' library(fdaPDE)
#'
#' ## Upload the hub2.5D the data
#' data(hub2.5D)
#' hub2.5D.nodes = hub2.5D$hub2.5D.nodes
#' hub2.5D.triangles = hub2.5D$hub2.5D.triangles
#'
#' ## Create mesh from nodes and connectivity matrix:
#' mesh = create.mesh.2.5D(nodes = hub2.5D.nodes, triangles = hub2.5D.triangles)
#' plot(mesh)
create.mesh.2.5D<- function(nodes, triangles = NULL, order = 1, nodesattributes = NULL, segments = NULL, holes = NULL)
{
##########################
### Input checking ###
##########################
segmentsmarkers <- vector(mode = "integer", 0)
nodes <- as.matrix(nodes)
if (ncol(nodes) != 3)
stop("Matrix of nodes should have 3 columns")
if (anyDuplicated(nodes))
stop("Duplicated nodes")
## If attributes not specified, set them to a matrix with zero columns
if (any(is.null(nodesattributes))) {
nodesattributes <- matrix(0, nrow(nodes), 0)
}else{
nodesattributes <- as.matrix(nodesattributes)
if (nrow(nodesattributes) != nrow(nodes))
stop("Point attribute matrix \'nodesattributes\' does not have same number of rows the point matrix \'nodes\'")
}
## Deal with segments
if (any(is.null(segments))) {
segments <- matrix(0, 0, 2)
} else {
segments <- as.matrix(segments)
if (ncol(segments) != 2) {
stop("Matrix of segments should have 2 columns")
}
}
if (any(is.null(holes)))
holes <- matrix(0, 0, 2)
holes <- as.matrix(holes)
## If triangles are not already specified
if(any(is.null(triangles))){
stop("Missing triangles argument is needed!")
# triangles = matrix(0,nrow = 0, ncol = 3)
} else {
triangles = as.matrix(triangles)
}
# Indexes in C++ starts from 0, in R from 1, needed transformations!
triangles = triangles - 1
## Set proper type for correct C++ reading
storage.mode(triangles) <- "integer"
storage.mode(nodes) <- "double"
out<-NULL
if(order==1 && ncol(triangles) == 3){
outCPP <- .Call("CPP_SurfaceMeshHelper", triangles, nodes, PACKAGE = "fdaPDE")
out <- list(nodes=nodes, nodesmarkers=outCPP[[3]], nodesattributes=nodesattributes,
triangles=triangles+1, segments=segments, segmentsmarkers=segmentsmarkers,
edges=outCPP[[1]], edgesmarkers=outCPP[[2]], neighbors=outCPP[[4]], holes=holes, order=order)
}
else if(order==2 && ncol(triangles) == 6){
outCPP <- .Call("CPP_SurfaceMeshHelper", triangles[,1:3], nodes, PACKAGE = "fdaPDE")
out <- list(nodes=nodes, nodesmarkers=outCPP[[3]], nodesattributes=nodesattributes,
triangles=triangles+1, segments=segments, segmentsmarkers=segmentsmarkers,
edges=outCPP[[1]], edgesmarkers=outCPP[[2]], neighbors=outCPP[[4]], holes=holes, order=order)
}
else if(order==2 && ncol(triangles) == 3){
message("You set order=2 but passed a matrix of triangles with just 3 columns. The midpoints for each edge will be computed.")
outCPP <- .Call("CPP_SurfaceMeshOrder2", triangles, nodes, PACKAGE = "fdaPDE")
out <- list(nodes=rbind(nodes, outCPP[[5]]), nodesmarkers=c(outCPP[[3]], outCPP[[2]]), nodesattributes=nodesattributes,
triangles=cbind(triangles+1, outCPP[[6]]), segments=segments, segmentsmarkers=segmentsmarkers,
edges=outCPP[[1]], edgesmarkers=outCPP[[2]], neighbors=outCPP[[4]], holes=holes, order=order)
if (ncol(out$nodesattributes)==0) {
out$nodesattributes <- matrix(0, nrow(out$nodes), 0)
}
}
else{
stop("The number of columns of triangles matrix is not consistent with the order parameter")
}
class(out)<-"mesh.2.5D"
return(out)
}
#' Project 3D points onto 2D 2.5D triangular mesh
#'
#' @param mesh A mesh.2.5D object representing the triangular mesh, created by \link{create.mesh.2.5D}.
#' @param locations 3D points to be projected onto 2.5D triangular mesh.
#' @description This function projects any 3D points onto 2.5D triangular mesh.
#' @return 3D points projected onto 2.5D triangluar mesh.
#' @export
#' @examples
#' library(fdaPDE)
#'
#' ## Upload the hub2.5D the data
#' data(hub2.5D)
#' hub2.5D.nodes = hub2.5D$hub2.5D.nodes
#' hub2.5D.triangles = hub2.5D$hub2.5D.triangles
#'
#' ## Create mesh
#' mesh = create.mesh.2.5D(nodes = hub2.5D.nodes, triangles = hub2.5D.triangles)
#'
#' ## Create 3D points to be projected
#' x <- cos(seq(0,2*pi, length.out = 9))
#' y <- sin(seq(0,2*pi, length.out = 9))
#' z <- rep(0.5, 9)
#' locations = cbind(x,y,z)
#'
#' ## Project the points on the mesh
#' loc = projection.points.2.5D(mesh, locations)
projection.points.2.5D<-function(mesh, locations) {
if(!is(mesh, "mesh.2.5D"))
stop("Data projection is only available for 2.5D mesh ")
mesh$triangles = mesh$triangles - 1
mesh$edges = mesh$edges - 1
mesh$neighbors[mesh$neighbors != -1] = mesh$neighbors[mesh$neighbors != -1] - 1
mydim=2
ndim=3
# Imposing types, this is necessary for correct reading from C++
## Set proper type for correct C++ reading
locations <- as.matrix(locations)
storage.mode(locations) <- "double"
storage.mode(mesh$nodes) <- "double"
storage.mode(mesh$triangles) <- "integer"
storage.mode(mesh$edges) <- "integer"
storage.mode(mesh$neighbors) <- "integer"
storage.mode(mesh$order) <- "integer"
storage.mode(mydim) <- "integer"
storage.mode(ndim) <-"integer"
## Call C++ function
evalmat <- .Call("points_projection", mesh, locations, mydim, ndim, PACKAGE = "fdaPDE")
#Returning the evaluation matrix
return(evalmat)
}
#' Create a \code{mesh.3D} object from the connectivity matrix and nodes locations
#'
#' @param nodes A #nodes-by-3 matrix containing the x, y, z coordinates of the mesh nodes.
#' @param tetrahedrons A #tetrahedrons-by-4 (when \code{order} = 1) or #tetrahedrons-by-10 (when \code{order} = 2) matrix.
#' It specifies the tetrahedrons giving the row's indices in \code{nodes} of the tetrahedrons' vertices and (when \code{nodes} = 2) also if the tetrahedrons' edges midpoints. The tetrahedrons' vertices and midpoints are ordered as described
#' in "The Finite Element Method its Basis and Fundamentals" by O. C. Zienkiewicz, R. L. Taylor and J.Z. Zhu
#' @param order Either '1' or '2'. It specifies wether each mesh tetrahedron should be represented by 4 nodes (the tetrahedron's vertices) or by 10 nodes (the tetrahedron's vertices and edge midpoints).
#' These are
#' respectively used for linear (order = 1) and quadratic (order = 2) Finite Elements. Default is \code{order} = 1.
#' @param nodesattributes A matrix with #nodes rows containing nodes' attributes.
#' These are passed unchanged to the output. This has been added for consistency with the function \code{create.mesh.2D}.
#' @param segments A #segments-by-2 matrix. Each row contains the row's indices in \code{nodes} of the vertices where the segment starts from and ends to.
#' Segments are edges that are not splitted during the triangulation process. These are for instance used to define the boundaries
#' of the domain. This has been added for consistency with the function \code{create.mesh.2D}.
#' @param holes A #holes-by-3 matrix containing the x, y, z coordinates of a point internal to each hole of the mesh. These points are used to carve holes
#' in the triangulation, when the domain has holes. This has been added for consistency with the function \code{create.mesh.2D}.
#' @return An object of the class mesh.3D with the following output:
#' \itemize{
#' \item{\code{nodes}}{A #nodes-by-3 matrix containing the x, y, z coordinates of the mesh nodes.}
#' \item{\code{nodesmarkers}}{A vector of length #nodes, with entries either '1' or '0'. An entry '1' indicates that the corresponding node is a boundary node; an entry '0' indicates that the corresponding node is not a boundary node.}
#' \item{\code{nodesattributes}}{A matrix with #nodes rows containing nodes' attributes.
#' These are passed unchanged from the input.}
#' \item{\code{tetrahedrons}}{A #tetrahedrons-by-4 (when \code{order} = 1) or #tetrahedrons-by-10 (when \code{order} = 2) matrix.
#' It specifies the tetrahedrons giving the indices in \code{nodes} of the tetrahedrons' vertices and (when \code{nodes} = 2) also if the tetrahedrons' edges midpoints.}
#' \item{\code{segmentsmarker}}{A vector of length #segments with entries either '1' or '0'. An entry '1' indicates that the corresponding element in \code{segments} is a boundary segment;
#' an entry '0' indicates that the corresponding segment is not a boundary segment.}
#' \item{\code{faces}}{A #faces-by-3 matrix containing all the faces of the tetrahedrons in the output triangulation. Each row contains the row's indices in \code{nodes}, indicating the nodes where the face starts from and ends to.}
#' \item{\code{facesmarkers}}{A vector of lenght #faces with entries either '1' or '0'. An entry '1' indicates that the corresponding element in \code{faces} is a boundary face;
#' an entry '0' indicates that the corresponding edge is not a boundary face.}
#' \item{\code{neighbors}}{A #triangles-by-4 matrix. Each row contains the indices of the four neighbouring tetrahedrons An entry '-1' indicates that
#' one face of the tetrahedrons is a boundary face.}
#' \item{\code{holes}}{A #holes-by-3 matrix containing the x, y, z coordinates of a point internal to each hole of the mesh. These points are used to carve holes
#' in the triangulation, when the domain has holes. These are passed unchanged from the input.}
#' \item{\code{order}}{Either '1' or '2'. It specifies wether each mesh tetrahedron should be represented by 3 nodes (the tetrahedron's vertices) or by 6 nodes (the tetrahedron's vertices and midpoints).
#' These are respectively used for linear (order = 1) and quadratic (order = 2) Finite Elements.}
#' }
#' @export
#' @examples
#' library(fdaPDE)
#'
#' ##Load the matrix nodes and tetrahedrons
#' data(sphere3Ddata)
#'
#' nodes=sphere3Ddata$nodes
#' tetrahedrons=sphere3Ddata$tetrahedrons
#'
#' ##Create the triangulated mesh from the connectivity matrix and nodes locations
#' mesh=create.mesh.3D(nodes,tetrahedrons)
create.mesh.3D<- function(nodes, tetrahedrons, order = 1, nodesattributes = NULL, segments = NULL, holes = NULL)
{
segmentsmarkers <- vector(mode = "integer", 0)
##########################
### Input checking ###
##########################
nodes = as.matrix(nodes)
if (ncol(nodes) != 3)
stop("Matrix of nodes should have 3 columns")
if (anyDuplicated(nodes))
stop("Duplicated nodes")
## If attributes not specified, set them to a matrix with zero columns
if (any(is.null(nodesattributes))) {
nodesattributes <- matrix(0, nrow(nodes), 0)
}else{
nodesattributes <- as.matrix(nodesattributes)
if (nrow(nodesattributes) != nrow(nodes))
stop("Point attribute matrix \'nodesattributes\' does not have same number of rows the point matrix \'nodes\'")
}
## Deal with segments
if (any(is.null(segments))) {
segments <- matrix(0, 0, 2)
} else {
segments <- as.matrix(segments)
if (ncol(segments) != 2) {
stop("Matrix of segments should have 2 columns")
}
}
if (any(is.null(holes)))
holes <- matrix(0, 0, 2)
holes <- as.matrix(holes)
## If tetrahedrons are not already specified
if(any(is.null(tetrahedrons))){
stop("Per il momento in questo caso serve tetrahedrons")
# triangles = matrix(0,nrow = 0, ncol = 3)
} else {
tetrahedrons <- as.matrix(tetrahedrons)
}
# Indexes in C++ starts from 0, in R from 1, needed transformations!
tetrahedrons = tetrahedrons - 1
## Set proper type for correct C++ reading
storage.mode(tetrahedrons) <- "integer"
storage.mode(nodes) <- "double"
out<-NULL
if(order==1 && ncol(tetrahedrons) == 4){
outCPP <- .Call("CPP_VolumeMeshHelper", tetrahedrons, nodes, PACKAGE = "fdaPDE")
out <- list(nodes=nodes, nodesmarkers=outCPP[[3]], nodesattributes=nodesattributes,
tetrahedrons=tetrahedrons+1, segments=segments, segmentsmarkers=segmentsmarkers,
faces=outCPP[[1]], facesmarkers=outCPP[[2]], neighbors=outCPP[[4]],
holes=holes, order=order)
}
else if(order==2 && ncol(tetrahedrons) == 10){
outCPP <- .Call("CPP_VolumeMeshHelper", tetrahedrons[,1:4], nodes, PACKAGE = "fdaPDE")
out <- list(nodes=nodes, nodesmarkers=outCPP[[3]], nodesattributes=nodesattributes,
tetrahedrons=tetrahedrons+1, segments=segments, segmentsmarkers=segmentsmarkers,
faces=outCPP[[1]], facesmarkers=outCPP[[2]], neighbors=outCPP[[4]],
holes=holes, order=order)
}
else if(order==2 && ncol(tetrahedrons) == 4){
message("You set order=2 but passed a matrix of tetrahedrons with just 4 columns. The midpoints for each edge will be computed.")
outCPP <- .Call("CPP_VolumeMeshOrder2", tetrahedrons, nodes, PACKAGE = "fdaPDE")
out <- list(nodes=rbind(nodes, outCPP[[5]]), nodesmarkers=c(outCPP[[3]], rep(FALSE, nrow(outCPP[[5]]))), nodesattributes=nodesattributes,
tetrahedrons=cbind(tetrahedrons+1, outCPP[[6]]), segments=segments, segmentsmarkers=segmentsmarkers,
faces=outCPP[[1]], facesmarkers=outCPP[[2]], neighbors=outCPP[[4]],
holes=holes, order=order)
if (ncol(out$nodesattributes)==0) {
out$nodesattributes <- matrix(0, nrow(out$nodes), 0)
}
}
else{
stop("The number of columns of tetrahedrons matrix is not consistent with the order parameter")
}
class(out)<-"mesh.3D"
return(out)
}
#' Create a \code{mesh.2D} object by splitting each triangle of a given mesh into four subtriangles.
#'
#' @param mesh a \code{mesh.2D} object to split
#' @return An object of class mesh.2D with splitted triangles
#' @export
refine.by.splitting.mesh.2D <- function (mesh=NULL){
if(is.null(mesh))
stop("No mesh passed as input!")
if(!is(mesh, "mesh.2D"))
stop("Wrong mesh class! Should be mesh.2D")
# Indexes in C++ starts from 0, in R from 1, needed transformations!
mesh$triangles = mesh$triangles - 1
storage.mode(mesh$triangles) <- "integer"
storage.mode(mesh$nodes) <- "double"
if(mesh$order==1){
outCPP <- .Call("CPP_TriangleMeshSplit", mesh$triangles[,1:3], mesh$nodes)
splittedmesh<-create.mesh.2D(nodes=rbind(mesh$nodes, outCPP[[2]]), triangles=outCPP[[1]])
}
else if(mesh$order==2){
nnodes<-nrow(mesh$nodes)-nrow(mesh$edges)
outCPP <- .Call("CPP_TriangleMeshSplitOrder2", mesh$triangles[,1:3], mesh$nodes[1:nnodes,])
splittedmesh<-create.mesh.2D(nodes=mesh$nodes, triangles=outCPP[[1]], order=2)
}
return(splittedmesh)
}
#' Create a \code{mesh.2.5D} object by splitting each triangle of a given mesh into four subtriangles.
#'
#' @param mesh a \code{mesh.2.5D} object to split
#' @return An object of class mesh.2.5D with splitted triangles
#' @export
refine.by.splitting.mesh.2.5D <- function (mesh=NULL){
if(is.null(mesh))
stop("No mesh passed as input!")
if(!is(mesh, "mesh.2.5D"))
stop("Wrong mesh class! Should be mesh.2.5D")
# Indexes in C++ starts from 0, in R from 1, needed transformations!
mesh$triangles = mesh$triangles - 1
storage.mode(mesh$triangles) <- "integer"
storage.mode(mesh$nodes) <- "double"
if(mesh$order==1){
outCPP <- .Call("CPP_TriangleMeshSplit", mesh$triangles[,1:3], mesh$nodes)
splittedmesh<-create.mesh.2.5D(nodes=rbind(mesh$nodes, outCPP[[2]]), triangles=outCPP[[1]])
}
else if(mesh$order==2){
nnodes<-nrow(mesh$nodes)-nrow(mesh$edges)
outCPP <- .Call("CPP_TriangleMeshSplitOrder2", mesh$triangles[,1:3], mesh$nodes[1:nnodes,])
splittedmesh<-create.mesh.2.5D(nodes=mesh$nodes, triangles=outCPP[[1]], order=2)
}
return(splittedmesh)
}
#' Create a \code{mesh.3D} object by splitting each tetrahedron of a given mesh into eight subtetrahedrons.
#'
#' @param mesh a \code{mesh.3D} object to split
#' @return An object of class mesh.3D with splitted tetrahedrons
#' @export
refine.by.splitting.mesh.3D <- function (mesh=NULL){
if(is.null(mesh))
stop("No mesh passed as input!")
if(!is(mesh, "mesh.3D"))
stop("Wrong mesh class! Should be mesh.3D")
# Indexes in C++ starts from 0, in R from 1, needed transformations!
mesh$tetrahedrons = mesh$tetrahedrons - 1
storage.mode(mesh$tetrahedrons) <- "integer"
storage.mode(mesh$nodes) <- "double"
if(mesh$order==1){
outCPP <- .Call("CPP_TetraMeshSplit", mesh$tetrahedrons[,1:4], mesh$nodes)
splittedmesh<-create.mesh.3D(nodes=rbind(mesh$nodes, outCPP[[2]]), tetrahedrons=outCPP[[1]])
}
else if(mesh$order==2){
nnodes<-max(mesh$tetrahedrons[,1:4])+1
outCPP <- .Call("CPP_TetraMeshSplitOrder2", mesh$tetrahedrons[,1:4], mesh$nodes[1:nnodes,])
splittedmesh<-create.mesh.3D(nodes=mesh$nodes, tetrahedrons=outCPP[[1]], order=2)
}
return(splittedmesh)
}
#' Create a 1.5D linear network mesh
#'
#' @param nodes A #nodes-by-2 matrix containing the x and y coordinates of the mesh nodes.
#' @param nodesattributes A matrix with #nodes rows containing nodes' attributes.
#' These are passed unchanged to the output. If a node is added during the triangulation process or mesh refinement, its attributes are computed
#' by linear interpolation using the attributes of neighboring nodes. This functionality is for instance used to compute the value
#' of a Dirichlet boundary condition at boundary nodes added during the triangulation process.
#' @param edges A #edges-by-2 (when \code{order} = 1) or #triangles-by-3 (when \code{order} = 2) matrix.
#' This option is used when a triangulation is already available. It specifies the edges giving the row's indices in \code{nodes} of the edges' vertices and (when \code{nodes} = 2) also if the triangles' edges midpoints. The triangles' vertices and midpoints are ordered as
#' 1---3---2
#' In this case the function \code{create.mesh.1.5D} is used to produce a complete mesh.1.5D object.
#' @param order Either '1' or '2'. It specifies wether each mesh should be represented by 2 nodes (the edges vertices) or by 3 nodes (the edges's vertices and midpoint).
#' These are
#' respectively used for linear (order = 1) and quadratic (order = 2) Finite Elements. Default is \code{order} = 1.
#' @usage create.mesh.1.5D(nodes, edges = NULL, order = 1, nodesattributes = NULL)
#' @return An object of the class mesh.1.5D with the following output:
#' \itemize{
#' \item{\code{nodes}}{A #nodes-by-2 matrix containing the x and y coordinates of the mesh nodes.}
#' \item{\code{nodesmarkers}}{A vector of length #nodes, with entries either '1' or '0'. An entry '1' indicates that the corresponding node is a boundary node; an entry '0' indicates that the corresponding node is not a boundary node.}
#' \item{\code{nodesattributes}}{A matrix with #nodes rows containing nodes' attributes.
#' These are passed unchanged from the input.}
#' \item{\code{edges}}{A #edges-by-2 matrix containing all the edges of the triangles in the output triangulation. Each row contains the row's indices in \code{nodes}, indicating the nodes where the edge starts from and ends to.}
#' \item{\code{neighbors}}{A #edges-by-2 matrix of list. Each row contains the indices of the neighbouring edges. An empty entry indicates that one node of the edge is a boundary node.}
#' \item{\code{order}}{Either '1' or '2'. It specifies wether each mesh triangle should be represented by 3 nodes (the triangle' vertices) or by 6 nodes (the triangle's vertices and midpoints).
#' These are respectively used for linear (order = 1) and quadratic (order = 2) Finite Elements.}
#' }
#' @export
create.mesh.1.5D <- function(nodes, edges = NULL, order = 1, nodesattributes = NULL)
{
nodes <- as.matrix(nodes)
if(ncol(nodes) != 2)
stop("Matrix of nodes should have 2 columns")
if(anyDuplicated(nodes))
stop("Duplicated nodes")
if(any(is.null(edges)))
stop("Missing edges argument is needed!")
else
edges = as.matrix(edges)
# Indexes in C++ starts from 0, in R from 1, needed transformations!
edges = edges - 1
storage.mode(edges) <- "integer"
storage.mode(nodes) <- "double"
out <- NULL
# length(out) == 11, according to create.mesh.2.5D / create.mesh.3D
if(order == 1 && ncol(edges) == 2){
outCPP <- .Call("CPP_EdgeMeshHelper", edges, nodes, PACKAGE = "fdaPDE")
#out <- list(nodes = nodes, nodesmarkers=outCPP[[2]], nodesattributes=nodesattributes,
# edges=edges+1,
# order=order, neighbors=outCPP[[4]], num_neighbors=outCPP[[3]])
out <- list(nodes = nodes, nodesmarkers=outCPP[[2]], nodesattributes=nodesattributes,
edges=edges+1)
out[[9]]<- outCPP[[4]]
names(out)[9] = "neighbors"
out[[11]]<- order
names(out)[11] <- "order"
}
else if(order==2 && ncol(edges) == 3){
edges_adj = matrix(edges[,1:2], nrow=dim(edges)[1],ncol=2)
storage.mode(edges_adj) <-"integer"
outCPP <- .Call("CPP_EdgeMeshHelper", edges_adj, nodes, PACKAGE = "fdaPDE")
#out <- list(nodes=nodes, nodesmarkers=outCPP[[2]], nodesattributes=nodesattributes,
# edges=edges+1,
# order=order,neighbors=outCPP[[4]], num_neighbors=outCPP[[3]])
out <- list(nodes = nodes, nodesmarkers=outCPP[[2]], nodesattributes=nodesattributes,
edges=edges+1)
out[[9]]<- outCPP[[4]]
names(out)[9] = "neighbors"
out[[11]]<- order
names(out)[11] <- "order"
}
else if( order==2 && ncol(edges)==2){
message("You set order=2 but passed a matrix of edges with just 2 columns. The midpoints for each edge will be computed.")
outCPP <- .Call("CPP_EdgeMeshOrder2", edges, nodes, PACKAGE = "fdaPDE")
edges = cbind( edges , outCPP[[6]])
nodes=rbind(nodes, outCPP[[5]])
out <- list(nodes=nodes, nodesmarkers=outCPP[[2]], nodesattributes=nodesattributes,
edges=edges+1)
out[[9]]<- outCPP[[4]]
names(out)[9] = "neighbors"
out[[11]]<- order
names(out)[11] <- "order"
}
class(out) <- "mesh.1.5D"
return(out)
}
#' Project 2D points onto 1.5D linear network mesh
#'
#' @param mesh A mesh.1.5D object representing the graph mesh, created by \link{create.mesh.1.5D}.
#' @param locations 2D points to be projected onto 1.5D mesh.
#' @description This function projects any 2D points onto 1.5D linear network mesh.
#' @return 2D points projected onto 1.5D linear network mesh.
#' @export
#' @examples
#' library(fdaPDE)
#'##Create Mesh
#'
#'nodes=matrix(c(0.25,0.25,0.5,0.25,0.75,0.5,0.75,0.), nrow = 4, byrow=TRUE)
#'edges=matrix(c(1,2,2,3,2,4),nrow = 3,byrow = TRUE)
#'mesh_ = create.mesh.1.5D(nodes,edges,order=1)
#'
#' ## Create 2D points to be projected
#'locations=matrix(nrow=5,ncol=2)
#'locations[,1] = runif(5,min=0.25,max=0.75)
#'locations[,2] = runif(5,min=0.25,max=0.5)
#'
#' ## Project the points on the mesh
#' loc = projection.points.1.5D(mesh_, locations)
projection.points.1.5D<-function(mesh, locations) {
if(!is(mesh, "mesh.1.5D"))
stop("Data projection is only available for 1.5D mesh ")
mesh$edges = mesh$edges - 1
mydim=1
ndim=2
# Imposing types, this is necessary for correct reading from C++
## Set proper type for correct C++ reading
locations <- as.matrix(locations)
storage.mode(locations) <- "double"
storage.mode(mesh$nodes) <- "double"
storage.mode(mesh$edges) <- "integer"
storage.mode(mesh$order) <- "integer"
storage.mode(mydim) <- "integer"
storage.mode(ndim) <-"integer"
## Call C++ function
evalmat <- .Call("points_projection", mesh, locations, mydim, ndim, PACKAGE = "fdaPDE")
#Returning the evaluation matrix
return(evalmat)
}
#' Refine 1.5D mesh
#'
#' @param mesh a \code{mesh.1.5D} object to refine
#' @param delta the maximum allowed length
#' @return An object of class mesh.1.5D with refined edges
#' @export
refine.mesh.1.5D <-function(mesh,delta){
if(is.null(mesh))
stop("No mesh passed as input!")
if(!is(mesh, "mesh.1.5D"))
stop("Wrong mesh class! Should be mesh.1.5D")
if( delta <= 0)
stop("Wrong delta value! Should be a positive number")
# Indexes in C++ starts from 0, in R from 1, needed transformations!
mesh$edges = mesh$edges - 1
edges <- as.matrix(mesh$edges[,1:2])
storage.mode(edges) <- "integer"
storage.mode(mesh$nodes) <- "double"
# outCPP[2] -> edges
# outCPP[1] -> new_nodes
nnodes<-max(mesh$edges[,1:2])+1
outCPP <- .Call("refine1D",mesh$nodes[1:nnodes,],edges,delta)
nodes <- rbind( mesh$nodes[1:nnodes,], outCPP[[1]])
edges <- outCPP[[2]]
ref_mesh <- create.mesh.1.5D(nodes = nodes, edges = edges, order = mesh$order)
return(ref_mesh)
}
#' Create a \code{mesh.1.5D} object by splitting each edge of a given mesh into two subedges.
#'
#' @param mesh a \code{mesh.1.5D} object to split
#' @return An object of class mesh.1.5D with splitted edges
#' @export
refine.by.splitting.mesh.1.5D <- function (mesh=NULL){
if(is.null(mesh))
stop("No mesh passed as input!")
if(!is(mesh, "mesh.1.5D"))
stop("Wrong mesh class! Should be mesh.1.5D")
# Indexes in C++ starts from 0, in R from 1, needed transformations!
mesh$edges = mesh$edges - 1
storage.mode(mesh$edges) <- "integer"
storage.mode(mesh$nodes) <- "double"
if(mesh$order==1){
outCPP <- .Call("CPP_EdgeMeshSplit", mesh$edges, mesh$nodes)
splittedmesh<-create.mesh.1.5D(nodes=rbind(mesh$nodes, outCPP[[2]]), edges=outCPP[[1]])
}
else if(mesh$order==2){
nnodes<-max(mesh$edges[,1:2])+1
edges_adj = as.matrix(mesh$edges)
storage.mode(edges_adj)<- "integer"
outCPP <- .Call("CPP_EdgeMeshSplit", edges_adj, mesh$nodes[1:nnodes,])
splittedmesh <- create.mesh.1.5D(nodes= rbind(mesh$nodes[1:nnodes,], outCPP[[2]]), edges = outCPP[[1]], order = 2)
}
return(splittedmesh)
}
Try the fdaPDE package in your browser
Any scripts or data that you put into this service are public.
fdaPDE documentation built on March 7, 2023, 5:28 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions refine.mesh.1.5D projection.points.1.5D create.mesh.1.5D create.mesh.3D projection.points.2.5D create.mesh.2.5D refine.mesh.2D create.mesh.2D triangulate_native
Documented in create.mesh.1.5D create.mesh.2.5D create.mesh.2D create.mesh.3D projection.points.1.5D projection.points.2.5D refine.mesh.1.5D refine.mesh.2D
triangulate_native <- function(P, PB, PA, S, SB,H, TR, flags) {
## It is necessary to check for NAs and NaNs, as the triangulate C
## code crashes if fed with them
P <- as.matrix(P)
PB <- as.integer(PB)
PA <- as.matrix(PA)
S <- as.matrix(S)
SB <- as.integer(SB)
H <- as.matrix(H)
TR <- as.matrix(TR)
storage.mode(P) <- "double"
storage.mode(PA) <- "double"
#storage.mode(PB) <- "integer"
storage.mode(S) <- "integer"
#storage.mode(SB) <- "integer"
storage.mode(H) <- "double"
storage.mode(TR) <- "integer"
storage.mode(flags) <- 'character'
## Call the main routine
out <- .Call("R_triangulate_native",
t(P),
PB,
PA,
t(S),
SB,
H,
t(TR),
flags,
PACKAGE="fdaPDE")
names(out) <- c("P", "PB", "PA", "T", "S", "SB", "E", "EB","TN", "VP", "VE", "VN", "VA")
class(out) <- "triangulation"
return(out)
}
#' Create a 2D triangular mesh
#'
#' @param nodes A #nodes-by-2 matrix containing the x and y coordinates of the mesh nodes.
#' @param nodesattributes A matrix with #nodes rows containing nodes' attributes.
#' These are passed unchanged to the output. If a node is added during the triangulation process or mesh refinement, its attributes are computed
#' by linear interpolation using the attributes of neighboring nodes. This functionality is for instance used to compute the value
#' of a Dirichlet boundary condition at boundary nodes added during the triangulation process.
#' @param segments A #segments-by-2 matrix. Each row contains the row's indices in \code{nodes} of the vertices where the segment starts from and ends to.
#' Segments are edges that are not splitted during the triangulation process. These are for instance used to define the boundaries
#' of the domain. If this is input is NULL, it generates a triangulation over the
#' convex hull of the points specified in \code{nodes}.
#' @param holes A #holes-by-2 matrix containing the x and y coordinates of a point internal to each hole of the mesh. These points are used to carve holes
#' in the triangulation, when the domain has holes.
#' @param triangles A #triangles-by-3 (when \code{order} = 1) or #triangles-by-6 (when \code{order} = 2) matrix.
#' This option is used when a triangulation is already available. It specifies the triangles giving the row's indices in \code{nodes} of the triangles' vertices and (when \code{nodes} = 2) also if the triangles' edges midpoints. The triangles' vertices and midpoints are ordered as described
#' at \cr https://www.cs.cmu.edu/~quake/triangle.highorder.html.
#' In this case the function \code{create.mesh.2D} is used to produce a complete mesh.2D object.
#' @param order Either '1' or '2'. It specifies wether each mesh triangle should be represented by 3 nodes (the triangle' vertices) or by 6 nodes (the triangle's vertices and midpoints).
#' These are
#' respectively used for linear (order = 1) and quadratic (order = 2) Finite Elements. Default is \code{order} = 1.
#' @param verbosity This can be '0', '1' or '2'. It indicates the level of verbosity in the triangulation process. When \code{verbosity} = 0 no message is returned
#' during the triangulation. When \code{verbosity} = 2 the triangulation process is described step by step by displayed messages.
#' Default is \code{verbosity} = 0.
#' @description This function is a wrapper of the Triangle library (http://www.cs.cmu.edu/~quake/triangle.html). It can be used
#' to create a triangulation of the domain of interest starting from a list of points, to be used as triangles' vertices, and a list of segments, that define the domain boundary. The resulting
#' mesh is a Constrained Delaunay triangulation. This is constructed in a way to preserve segments provided in the input \code{segments} without splitting them. This imput can be used to define the boundaries
#' of the domain. If this imput is NULL, it generates a triangulation over the
#' convex hull of the points.
#' It is also possible to create a mesh.2D from the nodes locations and the connectivity matrix.
#' @usage create.mesh.2D(nodes, nodesattributes = NA, segments = NA, holes = NA,
#' triangles = NA, order = 1, verbosity = 0)
#' @seealso \code{\link{refine.mesh.2D}}, \code{\link{create.FEM.basis}}
#' @return An object of the class mesh.2D with the following output:
#' \itemize{
#' \item{\code{nodes}}{A #nodes-by-2 matrix containing the x and y coordinates of the mesh nodes.}
#' \item{\code{nodesmarkers}}{A vector of length #nodes, with entries either '1' or '0'. An entry '1' indicates that the corresponding node is a boundary node; an entry '0' indicates that the corresponding node is not a boundary node.}
#' \item{\code{nodesattributes}}{A matrix with #nodes rows containing nodes' attributes.
#' These are passed unchanged from the input.}
#' \item{\code{triangles}}{A #triangles-by-3 (when \code{order} = 1) or #triangles-by-6 (when \code{order} = 2) matrix.
#' This option is used when a triangulation is already available. It specifies the triangles giving the indices in \code{nodes} of the triangles' vertices and (when \code{nodes} = 2) also if the triangles' edges midpoints. The triangles' vertices and midpoints are ordered as described
#' at \cr https://www.cs.cmu.edu/~quake/triangle.highorder.html.}
#' \item{\code{segmentsmarker}}{A vector of length #segments with entries either '1' or '0'. An entry '1' indicates that the corresponding element in \code{segments} is a boundary segment;
#' an entry '0' indicates that the corresponding segment is not a boundary segment.}
#' \item{\code{edges}}{A #edges-by-2 matrix containing all the edges of the triangles in the output triangulation. Each row contains the row's indices in \code{nodes}, indicating the nodes where the edge starts from and ends to.}
#' \item{\code{edgesmarkers}}{A vector of lenght #edges with entries either '1' or '0'. An entry '1' indicates that the corresponding element in \code{edge} is a boundary edge;
#' an entry '0' indicates that the corresponding edge is not a boundary edge.}
#' \item{\code{neighbors}}{A #triangles-by-3 matrix. Each row contains the indices of the three neighbouring triangles. An entry '-1' indicates that
#' one edge of the triangle is a boundary edge.}
#' \item{\code{holes}}{A #holes-by-2 matrix containing the x and y coordinates of a point internal to each hole of the mesh. These points are used to carve holes
#' in the triangulation, when the domain has holes.}
#' \item{\code{order}}{Either '1' or '2'. It specifies wether each mesh triangle should be represented by 3 nodes (the triangle' vertices) or by 6 nodes (the triangle's vertices and midpoints).
#' These are respectively used for linear (order = 1) and quadratic (order = 2) Finite Elements.}
#' }
#' @export
#' @examples
#' library(fdaPDE)
#'
#' ## Upload the quasicirle2D data
#' data(quasicircle2D)
#' boundary_nodes = quasicircle2D$boundary_nodes
#' boundary_segments = quasicircle2D$boundary_segments
#' locations = quasicircle2D$locations
#' data = quasicircle2D$data
#'
#' ## Create mesh from boundary
#' ## if the domain is convex it is sufficient to call:
#' mesh = create.mesh.2D(nodes = rbind(boundary_nodes, locations))
#' plot(mesh)
#'
#' ## if the domain is not convex, pass in addition the segments the compose the boundary:
#' mesh = create.mesh.2D(nodes = rbind(boundary_nodes, locations), segments = boundary_segments)
#'
#' ## Create mesh from data locations (without knowing the boundary)
#' mesh = create.mesh.2D(nodes = locations)
#' plot(mesh)
#' ## In this case the domain is the convex hull of the data locations.
#' ## Do this only if you do not have any information about the shape of the domain of interest.
create.mesh.2D <- function(nodes, nodesattributes = NA, segments = NA, holes = NA, triangles = NA, order = 1, verbosity = 0)
{
##########################
### Input checking ###
##########################
# Triangle finds out which are on the border (see https://www.cs.cmu.edu/~quake/triangle.help.html)
nodesmarkers = vector(mode = "integer", 0)
segmentsmarkers = vector(mode = "integer", 0)
nodes = as.matrix(nodes)
if (ncol(nodes) != 2)
stop("Matrix of nodes should have 2 columns")
if (anyDuplicated(nodes))
stop("Duplicated nodes")
## If attributes not specified, set them to a matrix with zero columns
if (any(is.na(nodesattributes))) {
nodesattributes <- matrix(0, nrow(nodes), 0)
}else{
nodesattributes <- as.matrix(nodesattributes)
if (nrow(nodesattributes) != nrow(nodes))
stop("Point attribute matrix \'nodesattributes\' does not have same number of rows the point matrix \'nodes\'")
}
## If boundary nodes not specified, set them to 0
# if (any(is.na(nodesmarkers))) {
# nodesmarkers <- vector(mode = "integer", 0)
# }else{
# nodesmarkers = as.vector(nodesmarkers)
# }
## Deal with segments
if (any(is.na(segments))) {
segments <- matrix(0, 0, 2)
} else {
segments <- as.matrix(segments)
if (ncol(segments) != 2) {
stop("Matrix of segments should have 2 columns")
}
}
## If boundary segments not specified, set them to 0
# if (any(is.na(segmentsmarkers))) {
# segmentsmarkers <- vector(mode = "integer", 0)
# }else{
# segmentsmarkers = as.vector(segmentsmarkers)
# }
## If hole not specified, set it to empty matrix
if (any(is.na(holes)))
holes <- matrix(0, 0, 2)
holes = as.matrix(holes)
## If triangles are not already specified
if(any(is.na(triangles)))
triangles = matrix(0,nrow = 0, ncol = 3)
triangles = as.matrix(triangles)
## Set meshing parameters ##
flags="ven"
if(nrow(segments) == 0){
flags = paste(flags,"c",sep = '')
}
if(nrow(segments)>0){
flags = paste(flags,"p",sep = '')
}
#If order=2 add flag for second order nodes
if(order == 2){
flags = paste(flags,"o2",sep = '')
}
if(order < 1 || order >2){
stop('Order must be 1 or 2')
}
if(nrow(triangles) > 0){
flags = paste(flags,"r",sep = '')
}
if (verbosity == 0) {
flags = paste(flags,"Q",sep = '')
}
if (verbosity == 1) {
flags = paste(flags,"V",sep = '')
}
if (verbosity == 2) {
flags = paste(flags,"VV",sep = '')
}
out<-NULL
#If triangles is null it makes the trianglulation
#If triangle is not null it makes a refinement with no parameter, to compose the mesh object
out <- triangulate_native(
nodes,
nodesmarkers,
nodesattributes,
segments,
segmentsmarkers,
t(holes),
triangles,
flags
)
names(out)[1]<-"nodes"
names(out)[2]<-"nodesmarkers"
names(out)[3]<-"nodesattributes"
names(out)[4]<-"triangles"
names(out)[5]<-"segments"
names(out)[6]<-"segmentsmarkers"
names(out)[7]<-"edges"
names(out)[8]<-"edgesmarkers"
names(out)[9]<-"neighbors"
out[13]<-NULL
out[12]<-NULL
out[11]<-NULL
out[10]<-NULL
out[[10]] = holes
names(out)[10]<-"holes"
out[[11]] = order
names(out)[11]<-"order"
class(out)<-"mesh.2D"
return(out)
}
#' Refine a 2D triangular mesh
#'
#' @param mesh A mesh.2D object representing the triangular mesh, created by \link{create.mesh.2D}.
#' @param minimum_angle A scalar specifying a minimun value for the triangles angles.
#' @param maximum_area A scalar specifying a maximum value for the triangles areas.
#' @param delaunay A boolean parameter indicating whether or not the output mesh should satisfy the Delaunay condition.
#' @param verbosity This can be '0', '1' or '2'. It indicates the level of verbosity in the triangulation process.
#' @description This function refines a Constrained Delaunay triangulation into a Conforming Delaunay triangulation. This is a wrapper of the Triangle library (http://www.cs.cmu.edu/~quake/triangle.html). It can be used to
#' refine a mesh previously created with \link{create.mesh.2D}. The algorithm can add Steiner points (points through which the \code{segments} are splitted)
#' in order to meet the imposed refinement conditions.
#' @usage refine.mesh.2D(mesh, minimum_angle, maximum_area, delaunay, verbosity)
#' @seealso \code{\link{create.mesh.2D}}, \code{\link{create.FEM.basis}}
#' @return A mesh.2D object representing the refined triangular mesh, with the following output:
#' \itemize{
#' \item{\code{nodes}}{A #nodes-by-2 matrix containing the x and y coordinates of the mesh nodes.}
#' \item{\code{nodesmarkers}}{A vector of length #nodes, with entries either '1' or '0'. An entry '1' indicates that the corresponding node is a boundary node; an entry '0' indicates that the corresponding node is not a boundary node.}
#' \item{\code{nodesattributes}}{nodesattributes A matrix with #nodes rows containing nodes' attributes.
#' These are passed unchanged to the output. If a node is added during the triangulation process or mesh refinement, its attributes are computed
#' by linear interpolation using the attributes of neighboring nodes. This functionality is for instance used to compute the value
#' of a Dirichlet boundary condition at boundary nodes added during the triangulation process.}
#' \item{\code{triangles}}{A #triangles-by-3 (when \code{order} = 1) or #triangles-by-6 (when \code{order} = 2) matrix.}
#' \item{\code{edges}}{A #edges-by-2 matrix. Each row contains the row's indices of the nodes where the edge starts from and ends to.}
#' \item{\code{edgesmarkers}}{A vector of lenght #edges with entries either '1' or '0'. An entry '1' indicates that the corresponding element in \code{edge} is a boundary edge;
#' an entry '0' indicates that the corresponding edge is not a boundary edge.}
#' \item{\code{neighbors}}{A #triangles-by-3 matrix. Each row contains the indices of the three neighbouring triangles. An entry '-1' indicates that
#' one edge of the triangle is a boundary edge.}
#' \item{\code{holes}}{A #holes-by-2 matrix containing the x and y coordinates of a point internal to each hole of the mesh. These points are used to carve holes
#' in the triangulation, when the domain has holes.}
#' \item{\code{order}}{Either '1' or '2'. It specifies wether each mesh triangle should be represented by 3 nodes (the triangle' vertices) or by 6 nodes (the triangle's vertices and midpoints).
#' These are respectively used for linear (order = 1) and quadratic (order = 2) Finite Elements.}
#' }
#' @examples
#' library(fdaPDE)
#'
#' ## Upload the quasicircle2D data
#' data(quasicircle2D)
#' boundary_nodes = quasicircle2D$boundary_nodes
#' boundary_segments = quasicircle2D$boundary_segments
#' locations = quasicircle2D$locations
#' data = quasicircle2D$data
#'
#' ## Create mesh from boundary:
#' mesh = create.mesh.2D(nodes = boundary_nodes, segments = boundary_segments)
#' plot(mesh)
#' ## Refine the mesh with the maximum area criterion:
#' finemesh = refine.mesh.2D(mesh = mesh, maximum_area = 0.1)
#' plot(finemesh)
#' ## Refine the mesh with the minimum angle criterion:
#' finemesh2 = refine.mesh.2D(mesh = mesh, minimum_angle = 30)
#' plot(finemesh2)
#' @export
refine.mesh.2D<-function(mesh, minimum_angle = NA, maximum_area = NA, delaunay = FALSE, verbosity = 0)
{
if(!is(mesh, "mesh.2D"))
stop("Sorry, this function is implemented just for mesh.2D class ")
flags="rpven"
if(!is.na(minimum_angle)){
flags <- paste(flags, "q", sprintf("%.12f", minimum_angle), sep='')
}
if(!is.na(maximum_area)){
flags <- paste(flags, "a", sprintf("%.12f", maximum_area), sep='')
}
if(delaunay){
flags <- paste(flags, "D", sep='')
}
if(mesh$order==2){
flags <- paste(flags, "o2", sep='')
}
if (verbosity == 0) {
flags = paste(flags,"Q",sep = '')
}
if (verbosity == 1) {
flags = paste(flags,"V",sep = '')
}
if (verbosity == 2) {
flags = paste(flags,"VV",sep = '')
}
# Triangle finds out which are on the border (see https://www.cs.cmu.edu/~quake/triangle.help.html)
mesh$nodesmarkers = vector(mode = "integer", 0)
mesh$segmentsmarkers = vector(mode = "integer", 0)
out<-NULL
#If triangles is null it makes the trianglulation
#If triangle is not null it makes a refinement with no parameter, to compose the mesh object
out <- triangulate_native(
mesh$nodes,
mesh$nodesmarkers,
mesh$nodesattributes,
mesh$segments,
mesh$segmentmarkers,
t(mesh$holes),
mesh$triangles,
flags
)
names(out)[1]<-"nodes"
names(out)[2]<-"nodesmarkers"
names(out)[3]<-"nodesattributes"
names(out)[4]<-"triangles"
names(out)[5]<-"segments"
names(out)[6]<-"segmentsmarkers"
names(out)[7]<-"edges"
names(out)[8]<-"edgesmarkers"
names(out)[9]<-"neighbors"
out[13]<-NULL
out[12]<-NULL
out[11]<-NULL
out[10]<-NULL
out[[10]] = mesh$holes
names(out)[10]<-"holes"
out[[11]] = mesh$order
names(out)[11]<-"order"
class(out)<-"mesh.2D"
return(out)
}
#' Create a \code{mesh.2.5D} object from the nodes locations and the connectivity matrix
#'
#' @param nodes A #nodes-by-3 matrix containing the x, y, z coordinates of the mesh nodes.
#' @param triangles A #triangles-by-3 (when \code{order} = 1) or #triangles-by-6 (when \code{order} = 2) matrix.
#' It specifies the triangles giving the row's indices in \code{nodes} of the triangles' vertices and (when \code{nodes} = 2) also if the triangles' edges midpoints. The triangles' vertices and midpoints are ordered as described
#' at \cr https://www.cs.cmu.edu/~quake/triangle.highorder.html.
#' @param order Either '1' or '2'. It specifies wether each mesh triangle should be represented by 3 nodes (the triangle' vertices) or by 6 nodes (the triangle's vertices and midpoints).
#' These are
#' respectively used for linear (order = 1) and quadratic (order = 2) Finite Elements. Default is \code{order} = 1.
#' @param nodesattributes A matrix with #nodes rows containing nodes' attributes.
#' These are passed unchanged to the output. This has been added for consistency with the function \code{create.mesh.2D}.
#' @param segments A #segments-by-2 matrix. Each row contains the row's indices in \code{nodes} of the vertices where the segment starts from and ends to.
#' Segments are edges that are not splitted during the triangulation process. These are for instance used to define the boundaries
#' of the domain. This has been added for consistency with the function \code{create.mesh.2D}.
#' @param holes A #holes-by-3 matrix containing the x, y, z coordinates of a point internal to each hole of the mesh. These points are used to carve holes
#' in the triangulation, when the domain has holes. This has been added for consistency with the function \code{create.mesh.2D}.
#' @return An object of the class mesh.2.5D with the following output:
#' \itemize{
#' \item{\code{nodes}}{A #nodes-by-3 matrix containing the x, y, z coordinates of the mesh nodes.}
#' \item{\code{nodesmarkers}}{A vector of length #nodes, with entries either '1' or '0'. An entry '1' indicates that the corresponding node is a boundary node; an entry '0' indicates that the corresponding node is not a boundary node.}
#' \item{\code{nodesattributes}}{A matrix with #nodes rows containing nodes' attributes.
#' These are passed unchanged from the input.}
#' \item{\code{triangles}}{A #triangles-by-3 (when \code{order} = 1) or #triangles-by-6 (when \code{order} = 2) matrix.
#' It specifies the triangles giving the indices in \code{nodes} of the triangles' vertices and (when \code{nodes} = 2) also if the triangles' edges midpoints. The triangles' vertices and midpoints are ordered as described
#' at \cr https://www.cs.cmu.edu/~quake/triangle.highorder.html.}
#' \item{\code{segmentsmarker}}{A vector of length #segments with entries either '1' or '0'. An entry '1' indicates that the corresponding element in \code{segments} is a boundary segment;
#' an entry '0' indicates that the corresponding segment is not a boundary segment.}
#' \item{\code{edges}}{A #edges-by-2 matrix containing all the edges of the triangles in the output triangulation. Each row contains the row's indices in \code{nodes}, indicating the nodes where the edge starts from and ends to.}
#' \item{\code{edgesmarkers}}{A vector of lenght #edges with entries either '1' or '0'. An entry '1' indicates that the corresponding element in \code{edge} is a boundary edge;
#' an entry '0' indicates that the corresponding edge is not a boundary edge.}
#' \item{\code{neighbors}}{A #triangles-by-3 matrix. Each row contains the indices of the three neighbouring triangles. An entry '-1' indicates that
#' one edge of the triangle is a boundary edge.}
#' \item{\code{holes}}{A #holes-by-3 matrix containing the x, y, z coordinates of a point internal to each hole of the mesh. These points are used to carve holes
#' in the triangulation, when the domain has holes. These are passed unchanged from the input.}
#' \item{\code{order}}{Either '1' or '2'. It specifies wether each mesh triangle should be represented by 3 nodes (the triangle' vertices) or by 6 nodes (the triangle's vertices and midpoints).
#' These are respectively used for linear (order = 1) and quadratic (order = 2) Finite Elements.}
#' }
#' @export
#' @examples
#' library(fdaPDE)
#'
#' ## Upload the hub2.5D the data
#' data(hub2.5D)
#' hub2.5D.nodes = hub2.5D$hub2.5D.nodes
#' hub2.5D.triangles = hub2.5D$hub2.5D.triangles
#'
#' ## Create mesh from nodes and connectivity matrix:
#' mesh = create.mesh.2.5D(nodes = hub2.5D.nodes, triangles = hub2.5D.triangles)
#' plot(mesh)
create.mesh.2.5D<- function(nodes, triangles = NULL, order = 1, nodesattributes = NULL, segments = NULL, holes = NULL)
{
##########################
### Input checking ###
##########################
segmentsmarkers <- vector(mode = "integer", 0)
nodes <- as.matrix(nodes)
if (ncol(nodes) != 3)
stop("Matrix of nodes should have 3 columns")
if (anyDuplicated(nodes))
stop("Duplicated nodes")
## If attributes not specified, set them to a matrix with zero columns
if (any(is.null(nodesattributes))) {
nodesattributes <- matrix(0, nrow(nodes), 0)
}else{
nodesattributes <- as.matrix(nodesattributes)
if (nrow(nodesattributes) != nrow(nodes))
stop("Point attribute matrix \'nodesattributes\' does not have same number of rows the point matrix \'nodes\'")
}
## Deal with segments
if (any(is.null(segments))) {
segments <- matrix(0, 0, 2)
} else {
segments <- as.matrix(segments)
if (ncol(segments) != 2) {
stop("Matrix of segments should have 2 columns")
}
}
if (any(is.null(holes)))
holes <- matrix(0, 0, 2)
holes <- as.matrix(holes)
## If triangles are not already specified
if(any(is.null(triangles))){
stop("Missing triangles argument is needed!")
# triangles = matrix(0,nrow = 0, ncol = 3)
} else {
triangles = as.matrix(triangles)
}
# Indexes in C++ starts from 0, in R from 1, needed transformations!
triangles = triangles - 1
## Set proper type for correct C++ reading
storage.mode(triangles) <- "integer"
storage.mode(nodes) <- "double"
out<-NULL
if(order==1 && ncol(triangles) == 3){
outCPP <- .Call("CPP_SurfaceMeshHelper", triangles, nodes, PACKAGE = "fdaPDE")
out <- list(nodes=nodes, nodesmarkers=outCPP[[3]], nodesattributes=nodesattributes,
triangles=triangles+1, segments=segments, segmentsmarkers=segmentsmarkers,
edges=outCPP[[1]], edgesmarkers=outCPP[[2]], neighbors=outCPP[[4]], holes=holes, order=order)
}
else if(order==2 && ncol(triangles) == 6){
outCPP <- .Call("CPP_SurfaceMeshHelper", triangles[,1:3], nodes, PACKAGE = "fdaPDE")
out <- list(nodes=nodes, nodesmarkers=outCPP[[3]], nodesattributes=nodesattributes,
triangles=triangles+1, segments=segments, segmentsmarkers=segmentsmarkers,
edges=outCPP[[1]], edgesmarkers=outCPP[[2]], neighbors=outCPP[[4]], holes=holes, order=order)
}
else if(order==2 && ncol(triangles) == 3){
message("You set order=2 but passed a matrix of triangles with just 3 columns. The midpoints for each edge will be computed.")
outCPP <- .Call("CPP_SurfaceMeshOrder2", triangles, nodes, PACKAGE = "fdaPDE")
out <- list(nodes=rbind(nodes, outCPP[[5]]), nodesmarkers=c(outCPP[[3]], outCPP[[2]]), nodesattributes=nodesattributes,
triangles=cbind(triangles+1, outCPP[[6]]), segments=segments, segmentsmarkers=segmentsmarkers,
edges=outCPP[[1]], edgesmarkers=outCPP[[2]], neighbors=outCPP[[4]], holes=holes, order=order)
if (ncol(out$nodesattributes)==0) {
out$nodesattributes <- matrix(0, nrow(out$nodes), 0)
}
}
else{
stop("The number of columns of triangles matrix is not consistent with the order parameter")
}
class(out)<-"mesh.2.5D"
return(out)
}
#' Project 3D points onto 2D 2.5D triangular mesh
#'
#' @param mesh A mesh.2.5D object representing the triangular mesh, created by \link{create.mesh.2.5D}.
#' @param locations 3D points to be projected onto 2.5D triangular mesh.
#' @description This function projects any 3D points onto 2.5D triangular mesh.
#' @return 3D points projected onto 2.5D triangluar mesh.
#' @export
#' @examples
#' library(fdaPDE)
#'
#' ## Upload the hub2.5D the data
#' data(hub2.5D)
#' hub2.5D.nodes = hub2.5D$hub2.5D.nodes
#' hub2.5D.triangles = hub2.5D$hub2.5D.triangles
#'
#' ## Create mesh
#' mesh = create.mesh.2.5D(nodes = hub2.5D.nodes, triangles = hub2.5D.triangles)
#'
#' ## Create 3D points to be projected
#' x <- cos(seq(0,2*pi, length.out = 9))
#' y <- sin(seq(0,2*pi, length.out = 9))
#' z <- rep(0.5, 9)
#' locations = cbind(x,y,z)
#'
#' ## Project the points on the mesh
#' loc = projection.points.2.5D(mesh, locations)
projection.points.2.5D<-function(mesh, locations) {
if(!is(mesh, "mesh.2.5D"))
stop("Data projection is only available for 2.5D mesh ")
mesh$triangles = mesh$triangles - 1
mesh$edges = mesh$edges - 1
mesh$neighbors[mesh$neighbors != -1] = mesh$neighbors[mesh$neighbors != -1] - 1
mydim=2
ndim=3
# Imposing types, this is necessary for correct reading from C++
## Set proper type for correct C++ reading
locations <- as.matrix(locations)
storage.mode(locations) <- "double"
storage.mode(mesh$nodes) <- "double"
storage.mode(mesh$triangles) <- "integer"
storage.mode(mesh$edges) <- "integer"
storage.mode(mesh$neighbors) <- "integer"
storage.mode(mesh$order) <- "integer"
storage.mode(mydim) <- "integer"
storage.mode(ndim) <-"integer"
## Call C++ function
evalmat <- .Call("points_projection", mesh, locations, mydim, ndim, PACKAGE = "fdaPDE")
#Returning the evaluation matrix
return(evalmat)
}
#' Create a \code{mesh.3D} object from the connectivity matrix and nodes locations
#'
#' @param nodes A #nodes-by-3 matrix containing the x, y, z coordinates of the mesh nodes.
#' @param tetrahedrons A #tetrahedrons-by-4 (when \code{order} = 1) or #tetrahedrons-by-10 (when \code{order} = 2) matrix.
#' It specifies the tetrahedrons giving the row's indices in \code{nodes} of the tetrahedrons' vertices and (when \code{nodes} = 2) also if the tetrahedrons' edges midpoints. The tetrahedrons' vertices and midpoints are ordered as described
#' in "The Finite Element Method its Basis and Fundamentals" by O. C. Zienkiewicz, R. L. Taylor and J.Z. Zhu
#' @param order Either '1' or '2'. It specifies wether each mesh tetrahedron should be represented by 4 nodes (the tetrahedron's vertices) or by 10 nodes (the tetrahedron's vertices and edge midpoints).
#' These are
#' respectively used for linear (order = 1) and quadratic (order = 2) Finite Elements. Default is \code{order} = 1.
#' @param nodesattributes A matrix with #nodes rows containing nodes' attributes.
#' These are passed unchanged to the output. This has been added for consistency with the function \code{create.mesh.2D}.
#' @param segments A #segments-by-2 matrix. Each row contains the row's indices in \code{nodes} of the vertices where the segment starts from and ends to.
#' Segments are edges that are not splitted during the triangulation process. These are for instance used to define the boundaries
#' of the domain. This has been added for consistency with the function \code{create.mesh.2D}.
#' @param holes A #holes-by-3 matrix containing the x, y, z coordinates of a point internal to each hole of the mesh. These points are used to carve holes
#' in the triangulation, when the domain has holes. This has been added for consistency with the function \code{create.mesh.2D}.
#' @return An object of the class mesh.3D with the following output:
#' \itemize{
#' \item{\code{nodes}}{A #nodes-by-3 matrix containing the x, y, z coordinates of the mesh nodes.}
#' \item{\code{nodesmarkers}}{A vector of length #nodes, with entries either '1' or '0'. An entry '1' indicates that the corresponding node is a boundary node; an entry '0' indicates that the corresponding node is not a boundary node.}
#' \item{\code{nodesattributes}}{A matrix with #nodes rows containing nodes' attributes.
#' These are passed unchanged from the input.}
#' \item{\code{tetrahedrons}}{A #tetrahedrons-by-4 (when \code{order} = 1) or #tetrahedrons-by-10 (when \code{order} = 2) matrix.
#' It specifies the tetrahedrons giving the indices in \code{nodes} of the tetrahedrons' vertices and (when \code{nodes} = 2) also if the tetrahedrons' edges midpoints.}
#' \item{\code{segmentsmarker}}{A vector of length #segments with entries either '1' or '0'. An entry '1' indicates that the corresponding element in \code{segments} is a boundary segment;
#' an entry '0' indicates that the corresponding segment is not a boundary segment.}
#' \item{\code{faces}}{A #faces-by-3 matrix containing all the faces of the tetrahedrons in the output triangulation. Each row contains the row's indices in \code{nodes}, indicating the nodes where the face starts from and ends to.}
#' \item{\code{facesmarkers}}{A vector of lenght #faces with entries either '1' or '0'. An entry '1' indicates that the corresponding element in \code{faces} is a boundary face;
#' an entry '0' indicates that the corresponding edge is not a boundary face.}
#' \item{\code{neighbors}}{A #triangles-by-4 matrix. Each row contains the indices of the four neighbouring tetrahedrons An entry '-1' indicates that
#' one face of the tetrahedrons is a boundary face.}
#' \item{\code{holes}}{A #holes-by-3 matrix containing the x, y, z coordinates of a point internal to each hole of the mesh. These points are used to carve holes
#' in the triangulation, when the domain has holes. These are passed unchanged from the input.}
#' \item{\code{order}}{Either '1' or '2'. It specifies wether each mesh tetrahedron should be represented by 3 nodes (the tetrahedron's vertices) or by 6 nodes (the tetrahedron's vertices and midpoints).
#' These are respectively used for linear (order = 1) and quadratic (order = 2) Finite Elements.}
#' }
#' @export
#' @examples
#' library(fdaPDE)
#'
#' ##Load the matrix nodes and tetrahedrons
#' data(sphere3Ddata)
#'
#' nodes=sphere3Ddata$nodes
#' tetrahedrons=sphere3Ddata$tetrahedrons
#'
#' ##Create the triangulated mesh from the connectivity matrix and nodes locations
#' mesh=create.mesh.3D(nodes,tetrahedrons)
create.mesh.3D<- function(nodes, tetrahedrons, order = 1, nodesattributes = NULL, segments = NULL, holes = NULL)
{
segmentsmarkers <- vector(mode = "integer", 0)
##########################
### Input checking ###
##########################
nodes = as.matrix(nodes)
if (ncol(nodes) != 3)
stop("Matrix of nodes should have 3 columns")
if (anyDuplicated(nodes))
stop("Duplicated nodes")
## If attributes not specified, set them to a matrix with zero columns
if (any(is.null(nodesattributes))) {
nodesattributes <- matrix(0, nrow(nodes), 0)
}else{
nodesattributes <- as.matrix(nodesattributes)
if (nrow(nodesattributes) != nrow(nodes))
stop("Point attribute matrix \'nodesattributes\' does not have same number of rows the point matrix \'nodes\'")
}
## Deal with segments
if (any(is.null(segments))) {
segments <- matrix(0, 0, 2)
} else {
segments <- as.matrix(segments)
if (ncol(segments) != 2) {
stop("Matrix of segments should have 2 columns")
}
}
if (any(is.null(holes)))
holes <- matrix(0, 0, 2)
holes <- as.matrix(holes)
## If tetrahedrons are not already specified
if(any(is.null(tetrahedrons))){
stop("Per il momento in questo caso serve tetrahedrons")
# triangles = matrix(0,nrow = 0, ncol = 3)
} else {
tetrahedrons <- as.matrix(tetrahedrons)
}
# Indexes in C++ starts from 0, in R from 1, needed transformations!
tetrahedrons = tetrahedrons - 1
## Set proper type for correct C++ reading
storage.mode(tetrahedrons) <- "integer"
storage.mode(nodes) <- "double"
out<-NULL
if(order==1 && ncol(tetrahedrons) == 4){
outCPP <- .Call("CPP_VolumeMeshHelper", tetrahedrons, nodes, PACKAGE = "fdaPDE")
out <- list(nodes=nodes, nodesmarkers=outCPP[[3]], nodesattributes=nodesattributes,
tetrahedrons=tetrahedrons+1, segments=segments, segmentsmarkers=segmentsmarkers,
faces=outCPP[[1]], facesmarkers=outCPP[[2]], neighbors=outCPP[[4]],
holes=holes, order=order)
}
else if(order==2 && ncol(tetrahedrons) == 10){
outCPP <- .Call("CPP_VolumeMeshHelper", tetrahedrons[,1:4], nodes, PACKAGE = "fdaPDE")
out <- list(nodes=nodes, nodesmarkers=outCPP[[3]], nodesattributes=nodesattributes,
tetrahedrons=tetrahedrons+1, segments=segments, segmentsmarkers=segmentsmarkers,
faces=outCPP[[1]], facesmarkers=outCPP[[2]], neighbors=outCPP[[4]],
holes=holes, order=order)
}
else if(order==2 && ncol(tetrahedrons) == 4){
message("You set order=2 but passed a matrix of tetrahedrons with just 4 columns. The midpoints for each edge will be computed.")
outCPP <- .Call("CPP_VolumeMeshOrder2", tetrahedrons, nodes, PACKAGE = "fdaPDE")
out <- list(nodes=rbind(nodes, outCPP[[5]]), nodesmarkers=c(outCPP[[3]], rep(FALSE, nrow(outCPP[[5]]))), nodesattributes=nodesattributes,
tetrahedrons=cbind(tetrahedrons+1, outCPP[[6]]), segments=segments, segmentsmarkers=segmentsmarkers,
faces=outCPP[[1]], facesmarkers=outCPP[[2]], neighbors=outCPP[[4]],
holes=holes, order=order)
if (ncol(out$nodesattributes)==0) {
out$nodesattributes <- matrix(0, nrow(out$nodes), 0)
}
}
else{
stop("The number of columns of tetrahedrons matrix is not consistent with the order parameter")
}
class(out)<-"mesh.3D"
return(out)
}
#' Create a \code{mesh.2D} object by splitting each triangle of a given mesh into four subtriangles.
#'
#' @param mesh a \code{mesh.2D} object to split
#' @return An object of class mesh.2D with splitted triangles
#' @export
refine.by.splitting.mesh.2D <- function (mesh=NULL){
if(is.null(mesh))
stop("No mesh passed as input!")
if(!is(mesh, "mesh.2D"))
stop("Wrong mesh class! Should be mesh.2D")
# Indexes in C++ starts from 0, in R from 1, needed transformations!
mesh$triangles = mesh$triangles - 1
storage.mode(mesh$triangles) <- "integer"
storage.mode(mesh$nodes) <- "double"
if(mesh$order==1){
outCPP <- .Call("CPP_TriangleMeshSplit", mesh$triangles[,1:3], mesh$nodes)
splittedmesh<-create.mesh.2D(nodes=rbind(mesh$nodes, outCPP[[2]]), triangles=outCPP[[1]])
}
else if(mesh$order==2){
nnodes<-nrow(mesh$nodes)-nrow(mesh$edges)
outCPP <- .Call("CPP_TriangleMeshSplitOrder2", mesh$triangles[,1:3], mesh$nodes[1:nnodes,])
splittedmesh<-create.mesh.2D(nodes=mesh$nodes, triangles=outCPP[[1]], order=2)
}
return(splittedmesh)
}
#' Create a \code{mesh.2.5D} object by splitting each triangle of a given mesh into four subtriangles.
#'
#' @param mesh a \code{mesh.2.5D} object to split
#' @return An object of class mesh.2.5D with splitted triangles
#' @export
refine.by.splitting.mesh.2.5D <- function (mesh=NULL){
if(is.null(mesh))
stop("No mesh passed as input!")
if(!is(mesh, "mesh.2.5D"))
stop("Wrong mesh class! Should be mesh.2.5D")
# Indexes in C++ starts from 0, in R from 1, needed transformations!
mesh$triangles = mesh$triangles - 1
storage.mode(mesh$triangles) <- "integer"
storage.mode(mesh$nodes) <- "double"
if(mesh$order==1){
outCPP <- .Call("CPP_TriangleMeshSplit", mesh$triangles[,1:3], mesh$nodes)
splittedmesh<-create.mesh.2.5D(nodes=rbind(mesh$nodes, outCPP[[2]]), triangles=outCPP[[1]])
}
else if(mesh$order==2){
nnodes<-nrow(mesh$nodes)-nrow(mesh$edges)
outCPP <- .Call("CPP_TriangleMeshSplitOrder2", mesh$triangles[,1:3], mesh$nodes[1:nnodes,])
splittedmesh<-create.mesh.2.5D(nodes=mesh$nodes, triangles=outCPP[[1]], order=2)
}
return(splittedmesh)
}
#' Create a \code{mesh.3D} object by splitting each tetrahedron of a given mesh into eight subtetrahedrons.
#'
#' @param mesh a \code{mesh.3D} object to split
#' @return An object of class mesh.3D with splitted tetrahedrons
#' @export
refine.by.splitting.mesh.3D <- function (mesh=NULL){
if(is.null(mesh))
stop("No mesh passed as input!")
if(!is(mesh, "mesh.3D"))
stop("Wrong mesh class! Should be mesh.3D")
# Indexes in C++ starts from 0, in R from 1, needed transformations!
mesh$tetrahedrons = mesh$tetrahedrons - 1
storage.mode(mesh$tetrahedrons) <- "integer"
storage.mode(mesh$nodes) <- "double"
if(mesh$order==1){
outCPP <- .Call("CPP_TetraMeshSplit", mesh$tetrahedrons[,1:4], mesh$nodes)
splittedmesh<-create.mesh.3D(nodes=rbind(mesh$nodes, outCPP[[2]]), tetrahedrons=outCPP[[1]])
}
else if(mesh$order==2){
nnodes<-max(mesh$tetrahedrons[,1:4])+1
outCPP <- .Call("CPP_TetraMeshSplitOrder2", mesh$tetrahedrons[,1:4], mesh$nodes[1:nnodes,])
splittedmesh<-create.mesh.3D(nodes=mesh$nodes, tetrahedrons=outCPP[[1]], order=2)
}
return(splittedmesh)
}
#' Create a 1.5D linear network mesh
#'
#' @param nodes A #nodes-by-2 matrix containing the x and y coordinates of the mesh nodes.
#' @param nodesattributes A matrix with #nodes rows containing nodes' attributes.
#' These are passed unchanged to the output. If a node is added during the triangulation process or mesh refinement, its attributes are computed
#' by linear interpolation using the attributes of neighboring nodes. This functionality is for instance used to compute the value
#' of a Dirichlet boundary condition at boundary nodes added during the triangulation process.
#' @param edges A #edges-by-2 (when \code{order} = 1) or #triangles-by-3 (when \code{order} = 2) matrix.
#' This option is used when a triangulation is already available. It specifies the edges giving the row's indices in \code{nodes} of the edges' vertices and (when \code{nodes} = 2) also if the triangles' edges midpoints. The triangles' vertices and midpoints are ordered as
#' 1---3---2
#' In this case the function \code{create.mesh.1.5D} is used to produce a complete mesh.1.5D object.
#' @param order Either '1' or '2'. It specifies wether each mesh should be represented by 2 nodes (the edges vertices) or by 3 nodes (the edges's vertices and midpoint).
#' These are
#' respectively used for linear (order = 1) and quadratic (order = 2) Finite Elements. Default is \code{order} = 1.
#' @usage create.mesh.1.5D(nodes, edges = NULL, order = 1, nodesattributes = NULL)
#' @return An object of the class mesh.1.5D with the following output:
#' \itemize{
#' \item{\code{nodes}}{A #nodes-by-2 matrix containing the x and y coordinates of the mesh nodes.}
#' \item{\code{nodesmarkers}}{A vector of length #nodes, with entries either '1' or '0'. An entry '1' indicates that the corresponding node is a boundary node; an entry '0' indicates that the corresponding node is not a boundary node.}
#' \item{\code{nodesattributes}}{A matrix with #nodes rows containing nodes' attributes.
#' These are passed unchanged from the input.}
#' \item{\code{edges}}{A #edges-by-2 matrix containing all the edges of the triangles in the output triangulation. Each row contains the row's indices in \code{nodes}, indicating the nodes where the edge starts from and ends to.}
#' \item{\code{neighbors}}{A #edges-by-2 matrix of list. Each row contains the indices of the neighbouring edges. An empty entry indicates that one node of the edge is a boundary node.}
#' \item{\code{order}}{Either '1' or '2'. It specifies wether each mesh triangle should be represented by 3 nodes (the triangle' vertices) or by 6 nodes (the triangle's vertices and midpoints).
#' These are respectively used for linear (order = 1) and quadratic (order = 2) Finite Elements.}
#' }
#' @export
create.mesh.1.5D <- function(nodes, edges = NULL, order = 1, nodesattributes = NULL)
{
nodes <- as.matrix(nodes)
if(ncol(nodes) != 2)
stop("Matrix of nodes should have 2 columns")
if(anyDuplicated(nodes))
stop("Duplicated nodes")
if(any(is.null(edges)))
stop("Missing edges argument is needed!")
else
edges = as.matrix(edges)
# Indexes in C++ starts from 0, in R from 1, needed transformations!
edges = edges - 1
storage.mode(edges) <- "integer"
storage.mode(nodes) <- "double"
out <- NULL
# length(out) == 11, according to create.mesh.2.5D / create.mesh.3D
if(order == 1 && ncol(edges) == 2){
outCPP <- .Call("CPP_EdgeMeshHelper", edges, nodes, PACKAGE = "fdaPDE")
#out <- list(nodes = nodes, nodesmarkers=outCPP[[2]], nodesattributes=nodesattributes,
# edges=edges+1,
# order=order, neighbors=outCPP[[4]], num_neighbors=outCPP[[3]])
out <- list(nodes = nodes, nodesmarkers=outCPP[[2]], nodesattributes=nodesattributes,
edges=edges+1)
out[[9]]<- outCPP[[4]]
names(out)[9] = "neighbors"
out[[11]]<- order
names(out)[11] <- "order"
}
else if(order==2 && ncol(edges) == 3){
edges_adj = matrix(edges[,1:2], nrow=dim(edges)[1],ncol=2)
storage.mode(edges_adj) <-"integer"
outCPP <- .Call("CPP_EdgeMeshHelper", edges_adj, nodes, PACKAGE = "fdaPDE")
#out <- list(nodes=nodes, nodesmarkers=outCPP[[2]], nodesattributes=nodesattributes,
# edges=edges+1,
# order=order,neighbors=outCPP[[4]], num_neighbors=outCPP[[3]])
out <- list(nodes = nodes, nodesmarkers=outCPP[[2]], nodesattributes=nodesattributes,
edges=edges+1)
out[[9]]<- outCPP[[4]]
names(out)[9] = "neighbors"
out[[11]]<- order
names(out)[11] <- "order"
}
else if( order==2 && ncol(edges)==2){
message("You set order=2 but passed a matrix of edges with just 2 columns. The midpoints for each edge will be computed.")
outCPP <- .Call("CPP_EdgeMeshOrder2", edges, nodes, PACKAGE = "fdaPDE")
edges = cbind( edges , outCPP[[6]])
nodes=rbind(nodes, outCPP[[5]])
out <- list(nodes=nodes, nodesmarkers=outCPP[[2]], nodesattributes=nodesattributes,
edges=edges+1)
out[[9]]<- outCPP[[4]]
names(out)[9] = "neighbors"
out[[11]]<- order
names(out)[11] <- "order"
}
class(out) <- "mesh.1.5D"
return(out)
}
#' Project 2D points onto 1.5D linear network mesh
#'
#' @param mesh A mesh.1.5D object representing the graph mesh, created by \link{create.mesh.1.5D}.
#' @param locations 2D points to be projected onto 1.5D mesh.
#' @description This function projects any 2D points onto 1.5D linear network mesh.
#' @return 2D points projected onto 1.5D linear network mesh.
#' @export
#' @examples
#' library(fdaPDE)
#'##Create Mesh
#'
#'nodes=matrix(c(0.25,0.25,0.5,0.25,0.75,0.5,0.75,0.), nrow = 4, byrow=TRUE)
#'edges=matrix(c(1,2,2,3,2,4),nrow = 3,byrow = TRUE)
#'mesh_ = create.mesh.1.5D(nodes,edges,order=1)
#'
#' ## Create 2D points to be projected
#'locations=matrix(nrow=5,ncol=2)
#'locations[,1] = runif(5,min=0.25,max=0.75)
#'locations[,2] = runif(5,min=0.25,max=0.5)
#'
#' ## Project the points on the mesh
#' loc = projection.points.1.5D(mesh_, locations)
projection.points.1.5D<-function(mesh, locations) {
if(!is(mesh, "mesh.1.5D"))
stop("Data projection is only available for 1.5D mesh ")
mesh$edges = mesh$edges - 1
mydim=1
ndim=2
# Imposing types, this is necessary for correct reading from C++
## Set proper type for correct C++ reading
locations <- as.matrix(locations)
storage.mode(locations) <- "double"
storage.mode(mesh$nodes) <- "double"
storage.mode(mesh$edges) <- "integer"
storage.mode(mesh$order) <- "integer"
storage.mode(mydim) <- "integer"
storage.mode(ndim) <-"integer"
## Call C++ function
evalmat <- .Call("points_projection", mesh, locations, mydim, ndim, PACKAGE = "fdaPDE")
#Returning the evaluation matrix
return(evalmat)
}
#' Refine 1.5D mesh
#'
#' @param mesh a \code{mesh.1.5D} object to refine
#' @param delta the maximum allowed length
#' @return An object of class mesh.1.5D with refined edges
#' @export
refine.mesh.1.5D <-function(mesh,delta){
if(is.null(mesh))
stop("No mesh passed as input!")
if(!is(mesh, "mesh.1.5D"))
stop("Wrong mesh class! Should be mesh.1.5D")
if( delta <= 0)
stop("Wrong delta value! Should be a positive number")
# Indexes in C++ starts from 0, in R from 1, needed transformations!
mesh$edges = mesh$edges - 1
edges <- as.matrix(mesh$edges[,1:2])
storage.mode(edges) <- "integer"
storage.mode(mesh$nodes) <- "double"
# outCPP[2] -> edges
# outCPP[1] -> new_nodes
nnodes<-max(mesh$edges[,1:2])+1
outCPP <- .Call("refine1D",mesh$nodes[1:nnodes,],edges,delta)
nodes <- rbind( mesh$nodes[1:nnodes,], outCPP[[1]])
edges <- outCPP[[2]]
ref_mesh <- create.mesh.1.5D(nodes = nodes, edges = edges, order = mesh$order)
return(ref_mesh)
}
#' Create a \code{mesh.1.5D} object by splitting each edge of a given mesh into two subedges.
#'
#' @param mesh a \code{mesh.1.5D} object to split
#' @return An object of class mesh.1.5D with splitted edges
#' @export
refine.by.splitting.mesh.1.5D <- function (mesh=NULL){
if(is.null(mesh))
stop("No mesh passed as input!")
if(!is(mesh, "mesh.1.5D"))
stop("Wrong mesh class! Should be mesh.1.5D")
# Indexes in C++ starts from 0, in R from 1, needed transformations!
mesh$edges = mesh$edges - 1
storage.mode(mesh$edges) <- "integer"
storage.mode(mesh$nodes) <- "double"
if(mesh$order==1){
outCPP <- .Call("CPP_EdgeMeshSplit", mesh$edges, mesh$nodes)
splittedmesh<-create.mesh.1.5D(nodes=rbind(mesh$nodes, outCPP[[2]]), edges=outCPP[[1]])
}
else if(mesh$order==2){
nnodes<-max(mesh$edges[,1:2])+1
edges_adj = as.matrix(mesh$edges)
storage.mode(edges_adj)<- "integer"
outCPP <- .Call("CPP_EdgeMeshSplit", edges_adj, mesh$nodes[1:nnodes,])
splittedmesh <- create.mesh.1.5D(nodes= rbind(mesh$nodes[1:nnodes,], outCPP[[2]]), edges = outCPP[[1]], order = 2)
}
return(splittedmesh)
}
Try the fdaPDE package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.