Nothing
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:
#' \describe{
#' \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:
#' \describe{
#' \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:
#' \describe{
#' \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:
#' \describe{
#' \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:
#' \describe{
#' \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)
}
CPP_search_points<-function(mesh, locations)
{
if(is(mesh, "mesh.2D")){
ndim = 2
mydim = 2
}else if(is(mesh, "mesh.1.5D")){
ndim = 2
mydim = 1
}else if(is(mesh, "mesh.2.5D")){
ndim = 3
mydim = 2
}else if(is(mesh, "mesh.3D")){
ndim = 3
mydim = 3
}else{
stop('Unknown mesh class')
}
## Set propr type for correct C++ reading
if( (ndim==2 && mydim==2) || (ndim==3 && mydim==2) ){
# Indexes in C++ starts from 0, in R from 1, opporGCV.inflation.factor transformation
mesh$triangles = mesh$triangles - 1
mesh$edges = mesh$edges - 1
mesh$neighbors[mesh$neighbors != -1] = mesh$neighbors[mesh$neighbors != -1] - 1
## Set propr type for correct C++ reading
storage.mode(mesh$triangles) <- "integer"
storage.mode(mesh$edges) <- "integer"
storage.mode(mesh$neighbors) <- "integer"
}else if( ndim==2 && mydim==1){
# Indexes in C++ starts from 0, in R from 1, opporGCV.inflation.factor transformation
mesh$edges = mesh$edges - 1
## Set propr type for correct C++ reading
storage.mode(mesh$edges) <- "integer"
}else if( ndim==3 && mydim==3){
# Indexes in C++ starts from 0, in R from 1, opporGCV.inflation.factor transformation
mesh$tetrahedrons = mesh$tetrahedrons - 1
mesh$faces = mesh$faces - 1
mesh$neighbors[mesh$neighbors != -1] = mesh$neighbors[mesh$neighbors != -1] - 1
## Set propr type for correct C++ reading
storage.mode(mesh$faces) <- "integer"
storage.mode(mesh$neighbors) <- "integer"
storage.mode(mesh$tetrahedrons) <- "integer"
}
storage.mode(mesh$order) <- "integer"
storage.mode(mesh$nodes) <- "double"
storage.mode(ndim)<-"integer"
storage.mode(mydim)<-"integer"
## Call C++ function
idx_elems <- .Call("points_search", mesh, locations, mydim, ndim,
PACKAGE = "fdaPDE")
return(idx_elems)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.