Nothing
R/plot.R
In fdaPDE: Physics-Informed Spatial and Functional Data Analysis
Defines functions
R_elementProperties
R_eval_local.FEM
R_image.ORDN.FEM
R_image.ORD1.FEM
image.FEM.time
image.FEM
R_plot_graph
R_plot_volume
R_plot_manifold
R_plot.ORDN.FEM
R_plot.ORD1.FEM
plot.mesh.1.5D
plot.mesh.3D
plot.mesh.2.5D
plot.mesh.2D
plot.FEM.time
plot.FEM
Documented in
image.FEM
image.FEM.time
plot.FEM
plot.FEM.time
plot.mesh.1.5D
plot.mesh.2.5D
plot.mesh.2D
plot.mesh.3D
#' Plot a \code{FEM} object
#'
#' @param x A \code{FEM} object.
#' @param colormap A colormap exploited in the plot. The default value is the heat colormap.
#' @param num_refinements A natural number specifying how many bisections should be applied to each triangular element for
#' plotting purposes. This functionality is useful where a discretization with 2nd order Finite Element is applied.
#' This parameter can be specified only when a FEM object defined over a 2D mesh is plotted.
#' @param ... Arguments representing graphical options to be passed to \link[rgl]{plot3d}.
#' @description Three-dimensional plot of a \code{FEM} object, generated by \code{FEM} or returned by
#' \code{smooth.FEM} or \code{FPCA.FEM}.
#' If the \code{mesh} of the \code{FEMbasis} component is of class \code{mesh.2D} both the 3rd axis and the color represent
#' the value of the coefficients for the Finite Element basis expansion (\code{coeff} component of the \code{FEM} object).
#' If the \code{mesh} is of class \code{mesh.3D}, the color of each triangle or tetrahedron represent the mean value of
#' the coefficients for the Finite Element basis expansion (\code{coeff}).
#' @usage \method{plot}{FEM}(x, colormap = "heat.colors", num_refinements = NULL, ...)
#' @export
#' @seealso \code{\link{FEM}}, \code{\link{image.FEM}}
#' @return No return value
#' @examples
#' library(fdaPDE)
#' ## Upload the horseshoe2D data
#' data(horseshoe2D)
#' boundary_nodes = horseshoe2D$boundary_nodes
#' boundary_segments = horseshoe2D$boundary_segments
#' locations = horseshoe2D$locations
#'
#' ## Create the 2D mesh
#' mesh = create.mesh.2D(nodes = rbind(boundary_nodes, locations), segments = boundary_segments)
#' ## Create the FEM basis
#' FEMbasis = create.FEM.basis(mesh)
#' ## Compute the coeff vector evaluating the desired function at the mesh nodes
#' ## In this case we consider the fs.test() function introduced by Wood et al. 2008
#' coeff = fs.test(mesh$nodes[,1], mesh$nodes[,2])
#' ## Create the FEM object
#' FEMfunction = FEM(coeff, FEMbasis)
#'
#' ## Plot the FEM function
#' plot(FEMfunction)
plot.FEM = function(x, colormap = "heat.colors", num_refinements = NULL, ...)
{
if(is(x$FEMbasis$mesh, "mesh.2D")){
if(x$FEMbasis$order == 1)
{
R_plot.ORD1.FEM(x, colormap, ...)
}else{
R_plot.ORDN.FEM(x, colormap, num_refinements, ...)
}
}else if(is(x$FEMbasis$mesh, "mesh.2.5D")){
R_plot_manifold(x,...)
}else if(is(x$FEMbasis$mesh, "mesh.3D")){
R_plot_volume(x,...)
}else if(is(x$FEMbasis$mesh, "mesh.1.5D")){
R_plot_graph(x, ...)
}
}
#' Plot a \code{FEM.time} object at a given time
#'
#' @param x A \code{FEM.time} object.
#' @param time_locations A vector with the instants in which plot the spatial field
#' @param locations A 2-column matrix when \code{x$FEMbasis$mesh} is of class \code{mesh.2D}
#' or a 3-column matrix otherwise with the spatial locations for which plot the temporal evolution
#' @param lambdaS Index of the space penalization parameter to use for the plot, useful when \code{FEM.time} returned by \code{smooth.FEM.time} using GCV
#' @param lambdaT Index of the time penalization parameter to use for the plot, useful when \code{FEM.time} returned by \code{smooth.FEM.time} using GCV
#' @param num_refinements A natural number specifying how many bisections should be applied to each triangular element for
#' plotting purposes. This functionality is useful where a discretization with 2nd order Finite Element is applied.
#' This parameter can be specified only when a FEM object defined over a 2D mesh is plotted.
#' @param Nt The number of instants to plot when \code{locations} is not NULL
#' @param add Boolean, used only when locations is not NULL, if TURE it performs the graphic over the last plot
#' @param main The title of the plot when \code{locations} is not NULL
#' @param col The color of the plot when \code{locations} is not NULL. May be a single color or a vector of colors
#' @param ... Arguments representing graphical options to be passed to \link[rgl]{plot3d}.
#' @description Plot of a \code{FEM.time} object, generated by \code{FEM.time} or returned by
#' \code{smooth.FEM.time}. \code{time_locations} and \code{locations} must not be both provided.
#' If \code{time_locations} is provided, the spatial field is plotted for the provided temporal instnts.
#' If \code{locations} is provided, the temporal evolution in the provided space locations is plotted.
#' If both \code{time_locations} and \code{locations} are NULL a default plot is provided.
#' If the \code{mesh} of the \code{FEMbasis} component is of class \code{mesh.2D} both the 3rd axis and the color represent
#' the value of the coefficients for the Finite Element basis expansion (\code{coeff} component of the \code{FEM.time} object).
#' If the \code{mesh} is of class \code{mesh.3D}, the color of each triangle or tetrahedron represent the mean value of
#' the coefficients for the Finite Element basis expansion (\code{coeff}).
#' @usage \method{plot}{FEM.time}(x, time_locations = NULL, locations = NULL,
#' lambdaS = NULL, lambdaT = NULL, num_refinements = NULL, Nt = 100,
#' add = FALSE, main = NULL, col = "red", ...)
#' @export
#' @seealso \code{\link{FEM.time}}, \code{\link{image.FEM.time}}
#' @return No return value
#' @examples
#' library(fdaPDE)
#' ## Upload the horseshoe2D data
#' data(horseshoe2D)
#' boundary_nodes = horseshoe2D$boundary_nodes
#' boundary_segments = horseshoe2D$boundary_segments
#' locations = horseshoe2D$locations
#'
#' ## Create the 2D mesh
#' mesh = create.mesh.2D(nodes = rbind(boundary_nodes, locations), segments = boundary_segments)
#' ## Create the FEM basis
#' FEMbasis = create.FEM.basis(mesh)
#' ## Compute the coeff vector evaluating the desired function at the mesh nodes
#' ## In this case we consider the fs.test() function introduced by Wood et al. 2008
#' time = 1:5
#' coeff = rep(fs.test(mesh$nodes[,1], mesh$nodes[,2]),5)*time
#' ## Create the FEM.time object
#' FEM_time_function = FEM.time(coeff=coeff, time_mesh=1:5, FEMbasis=FEMbasis, FLAG_PARABOLIC=TRUE)
#'
#' ## Plot the FEM function
#' plot(FEM_time_function)
#'
#' ## plot spatial field in some instants
#' t = c(1.2,1.5,3.6,2.4,4.5)
#' plot(FEM_time_function, t)
#'
#' ## plot time evolution in some locations
#' plot(FEM_time_function, locations = locations[1:10,])
plot.FEM.time = function(x, time_locations = NULL, locations = NULL,
lambdaS = NULL, lambdaT = NULL, num_refinements = NULL, Nt = 100,
add = FALSE, main = NULL, col = "red", ...)
{
if(!is(x, "FEM.time"))
stop("x is not of class 'FEM.time'")
if (!is.null(time_locations) && !is.null(locations))
stop("time_locations and locations can not be both provided")
if (is.null(time_locations) && is.null(locations)) {
# default plot
# space evolution:
time_locations = mean(x$mesh_time)
if(is.null(lambdaS) || is.null(lambdaT)){
if(exists("x$bestlambda")){
lambdaS = x$bestlambda[1]
lambdaT = x$bestlambda[2]
}else{
lambdaS = 1
lambdaT = 1
}
}
N = nrow(x$FEMbasis$mesh$nodes)
storage.mode(time_locations) <- "double"
storage.mode(N) <- "integer"
storage.mode(x$mesh_time) <- "double"
storage.mode(x$coeff) <- "double"
storage.mode(x$FLAG_PARABOLIC) <- "integer"
solution <- .Call("eval_FEM_time_nodes",N,x$mesh_time,time_locations,x$coeff[,lambdaS,lambdaT],x$FLAG_PARABOLIC, PACKAGE = "fdaPDE")
plot = FEM(solution[1:N],x$FEMbasis)
plot.FEM(x = plot, num_refinements = num_refinements,main = paste0("Spatial function in t = ",time_locations),...)
# time evolution:
locations = x$FEMbasis$mesh$nodes[1,]
# add default title to the plot
if(is.null(main))
main = paste0("Time evolution in (",round(locations,2),")")
# define the evaluation points
t = seq(min(x$mesh_time),max(x$mesh_time),length.out=Nt)
eval_points = cbind(t, t(kronecker(matrix(1,ncol=Nt), locations)))
# evaluate the function
eval_sol = eval.FEM.time(FEM.time = x, space.time.locations = eval_points, lambdaT=lambdaT, lambdaS=lambdaS)
# define variables for the plot
type = "l"
ylab = ""
xlab = "Time"
ylim = range(eval_sol)
# plot the function evaluations
if (add == FALSE) {
plot(t[1:Nt],eval_sol[1:Nt],type=type,xlab=xlab,ylab=ylab,main=main,col=col,ylim=ylim, ...)
}else{
points(t[1:Nt],eval_sol[1:Nt],type=type,col=col, ...)
}
}
if (!is.null(time_locations) && is.null(locations)) {
# plot the spatial function for the provided time_locations
# check that time_locations is a vector
if (!is.vector(time_locations) && nrow(time_locations)!=1 && ncol(time_locations)!=1)
stop("time_locations must be a vector of a column/row matrix")
time_locations = as.vector(time_locations)
# check that time_locations is in the correct range
if(min(time_locations)<x$mesh_time[1] || max(time_locations)>x$mesh_time[length(x$mesh_time)])
stop("time_locations must be within the 'time_mesh' range")
if(is.null(lambdaS) || is.null(lambdaT)){
if(exists("x$bestlambda")){
lambdaS = x$bestlambda[1]
lambdaT = x$bestlambda[2]
}else{
lambdaS = 1
lambdaT = 1
}
}
N = nrow(x$FEMbasis$mesh$nodes)
storage.mode(time_locations) <- "double"
storage.mode(N) <- "integer"
storage.mode(x$mesh_time) <- "double"
storage.mode(x$coeff) <- "double"
storage.mode(x$FLAG_PARABOLIC) <- "integer"
solution <- .Call("eval_FEM_time_nodes",N,x$mesh_time,time_locations,x$coeff[,lambdaS,lambdaT],x$FLAG_PARABOLIC, PACKAGE = "fdaPDE")
for(i in 1:length(t))
{
plot = FEM(solution[(1+(i-1)*N):(N+(i-1)*N)],x$FEMbasis)
plot.FEM(x = plot, num_refinements = num_refinements,...)
}
}
if (is.null(time_locations) && !is.null(locations)) {
# plot temporal evolution in the provided locations
if (is(x$FEMbasis$mesh, "mesh.2D") | is(x$FEMbasis$mesh, "mesh.1.5D")){
# check locations dimension
if (is.vector(locations) && length(locations) != 2)
stop("locations does not have the correct dimension")
if (!is.vector(locations) && ncol(locations) != 2)
stop("locations does not have the correct dimension")
if (is.vector(locations))
locations = matrix(locations, ncol = 2)
# add default title to the plot
if(is.null(main) && nrow(locations)==1)
main = paste0("Time evolution in (",locations[1,1],",",locations[1,2],")")
# define the evaluation points
t = rep(seq(min(x$mesh_time),max(x$mesh_time),length.out=Nt),times=nrow(locations))
eval_points = cbind(t, rep(locations[,1], each = Nt), rep(locations[,2], each = Nt))
}
if (is(x$FEMbasis$mesh, "mesh.2.5D") | is(x$FEMbasis$mesh, "mesh.3D")){
# check locations dimension
if (is.vector(locations) && length(locations) != 3)
stop("locations does not have the correct dimension")
if (!is.vector(locations) && ncol(locations) != 3)
stop("locations does not have the correct dimension")
if (is.vector(locations))
locations = matrix(locations, ncol = 3)
# add default title to the plot
if(is.null(main) && nrow(locations)==1)
main = paste0("Time evolution in (",locations[1,1],",",locations[1,2],",",locations[1,3],")")
# define the evaluation points
t = rep(seq(min(x$mesh_time),max(x$mesh_time),length.out=Nt),times=nrow(locations))
eval_points = cbind(t, rep(locations[,1], each = Nt),
rep(locations[,2], each = Nt), rep(locations[,3], each = Nt))
}
# evaluate the function
eval_sol = eval.FEM.time(FEM.time = x, space.time.locations = eval_points, lambdaT=lambdaT, lambdaS=lambdaS)
# define variables for the plot
if (length(col) < nrow(locations)){
# recicle colors if not provided enough
col = rep(col, nrow(locations))
}
if (!exists("type")) {
type = "l"
}
if (!exists("ylab")) {
ylab = ""
}
if (!exists("xlab")) {
xlab = "Time"
}
if (!exists("ylim")) {
ylim = range(eval_sol)
}
# plot the function evaluations
for (ind in 1:nrow(locations)) {
if (add == FALSE && ind == 1) {
plot(t[1:Nt],eval_sol[(ind-1)*Nt+(1:Nt)],type=type,xlab=xlab,ylab=ylab,main=main,col=col[ind],ylim=ylim, ...)
}else{
points(t[1:Nt],eval_sol[(ind-1)*Nt+(1:Nt)],type=type,col=col[ind], ...)
}
}
}
}
#' Plot a mesh.2D object
#'
#' @param x A \code{mesh.2D} object defining the triangular mesh, as generated by \code{create.mesh.2D}
#' or \code{refine.mesh.2D}.
#' @param ... Arguments representing graphical options to be passed to \link[graphics]{par}.
#' @description Plot a mesh.2D object, generated by \code{create.mesh.2D} or \code{refine.mesh.2D}.
#' @name plot.mesh.2D
#' @usage \method{plot}{mesh.2D}(x, ...)
#' @export
#' @return No return value
#' @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
#' mesh = create.mesh.2D(nodes = rbind(boundary_nodes, locations), segments = boundary_segments)
#'
#' ## Plot the mesh
#' plot(mesh)
plot.mesh.2D<-function(x, ...)
{
plot(x$nodes, xlab="", ylab="", xaxt="n", yaxt="n", bty="n", ...)
segments(x$nodes[x$edges[,1],1], x$nodes[x$edges[,1],2],
x$nodes[x$edges[,2],1], x$nodes[x$edges[,2],2], ...)
segments(x$nodes[x$segments[,1],1], x$nodes[x$segments[,1],2],
x$nodes[x$segments[,2],1], x$nodes[x$segments[,2],2], col="red", ...)
}
#' Plot a mesh.2.5D object
#'
#' @param x A \code{mesh.2.5D} object generated by \code{create.mesh.2.5D}.
#' @param ... Arguments representing graphical options to be passed to \link[graphics]{par}.
#' @description Plot the triangulation of a \code{mesh.2.5D} object, generated by \code{create.mesh.2.5D}
#' @export
#' @name plot.mesh.2.5D
#' @usage \method{plot}{mesh.2.5D}(x, ...)
#' @return No return value
#' @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)
#' plot(mesh)
plot.mesh.2.5D<-function(x,...){
nodes <- x$nodes
edges <- as.vector(t(x$edges))
open3d()
axes3d()
pop3d("lights")
light3d(specular="black")
points3d(nodes[,1], nodes[,2], nodes[,3], col="black", ...)
segments3d(nodes[edges,1], nodes[edges,2], nodes[edges,3], col="black",...)
aspect3d("iso")
view3d(0,-45)
}
#' Plot a mesh.3D object
#'
#' @param x A \code{mesh.3D} object generated by \code{create.mesh.3D}.
#' @param ... Arguments representing graphical options to be passed to \link[graphics]{par}.
#' @description Plot a \code{mesh.3D} object, generated by \code{create.mesh.3D}.
#' @export
#' @name plot.mesh.3D
#' @usage \method{plot}{mesh.3D}(x, ...)
#' @return No return value
#' @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)
#'
#' ##Plot the triangulation of the object
#' plot(mesh)
plot.mesh.3D<-function(x,...){
nodes <- x$nodes
faces <- as.vector(t(x$faces[x$facesmarkers,]))
aux_mesh <- create.mesh.2.5D(nodes=nodes, triangles=x$faces[x$facesmarkers,], order=1)
edges <- as.vector(t(aux_mesh$edges))
open3d()
axes3d()
pop3d("lights")
light3d(specular="black")
triangles3d(nodes[faces,1],nodes[faces,2],nodes[faces,3],col="white",...)
points3d(nodes[,1], nodes[,2], nodes[,3], col="black", ...)
segments3d(nodes[edges,1], nodes[edges,2], nodes[edges,3], col="black",...)
aspect3d("iso")
view3d(0,-45)
}
#' Plot a mesh.1.5D object
#'
#' @param x A \code{mesh.1.5D} object defining the triangular mesh, as generated by \code{create.mesh.1.5D}
#' or \code{refine.mesh.1.5D}.
#' @param ... Arguments representing graphical options to be passed to \link[graphics]{par}.
#' @description Plot a mesh.1.5D object, generated by \code{create.mesh.1.5D} or \code{refine.mesh.1.5D}.
#' @name plot.mesh.1.5D
#' @usage \method{plot}{mesh.1.5D}(x, ...)
#' @return No return value
#' @export
plot.mesh.1.5D<-function(x, ...)
{
if( x$order == 1 ){
plot(x$nodes, xlab="", ylab="", xaxt="n", yaxt="n", bty="n", ...)
segments(x$nodes[x$edges[,1],1], x$nodes[x$edges[,1],2],
x$nodes[x$edges[,2],1], x$nodes[x$edges[,2],2], ...)
}
else{
plot(x$nodes, xlab="", ylab="", xaxt="n", yaxt="n", bty="n", ...)
segments(x$nodes[x$edges[,1],1], x$nodes[x$edges[,1],2],
x$nodes[x$edges[,2],1], x$nodes[x$edges[,2],2], ...)
points(x$nodes[x$edges[,3],1], x$nodes[x$edges[,3],2],col="red", ...)
}
}
R_plot.ORD1.FEM = function(FEM, colormap, ...)
{
# PLOT Plots a FEM object FDOBJ over a rectangular grid defined by
# vectors X and Y;
#
nodes <- FEM$FEMbasis$mesh$nodes
triangles <- as.vector(t(FEM$FEMbasis$mesh$triangles))
colormap <- match.fun(colormap)
heat <- colormap(100)
# How many plots are needed?
nplots <- ncol(FEM$coeff)
for (i in 1:nplots) {
if (i > 1)
readline("Press any key for the next plot...")
open3d()
axes3d()
z <- FEM$coeff[triangles,i]
triangles3d(nodes[triangles,1], nodes[triangles,2], z,
color = heat[round(99*(z-min(z))/diff(range(z)))+1],...)
aspect3d(2,2,1)
view3d(0,-45)
}
}
R_plot.ORDN.FEM = function(FEM, colormap, num_refinements, ...)
{
# num_refinements sets the number od division on each triangle edge to be applied for rifenment
coeff = FEM$coeff
FEMbasis = FEM$FEMbasis
mesh = FEMbasis$mesh
colormap <- match.fun(colormap)
heat <- colormap(100)
coeff = FEM$coeff
if(is.null(num_refinements))
{
num_refinements = 10
}
# For the reference triangles we construct a regular mesh
x = seq(from = 0, to = 1, length.out = num_refinements+1)
y = seq(from = 0, to = 1, length.out = num_refinements+1)
points_ref = expand.grid(x,y)
points_ref = points_ref[points_ref[,1] + points_ref[,2] <= 1,]
meshi = create.mesh.2D(nodes = points_ref, order = 1)
#plot(meshi)
# locations is the matrix with that will contain the coordinate of the points where the function is
# evaluated (1st and 2nd column) and the columns with the evaluation of the ith fucntion on that point
locations = matrix(nrow = nrow(mesh$triangles)*nrow(meshi$nodes), ncol = 2+ncol(coeff))
triangles = matrix(nrow = nrow(mesh$triangles)*nrow(meshi$triangles), ncol = 3)
tot = 0
properties<-R_elementProperties(mesh)
for (i in 1:nrow(mesh$triangles))
{
# For each traingle we define a fine mesh as the transofrmation of the one constructed for the reference
transf<-rbind(cbind(properties$transf_coord$diff1x[i],properties$transf_coord$diff2x[i]),c(properties$transf_coord$diff1y[i],properties$transf_coord$diff2y[i]))
pointsi = t(transf%*%t(meshi$nodes) + mesh$nodes[mesh$triangles[i,1],])
#We evaluate the fine mesh OBS: we know the triangle we are working on no need for point location
z = R_eval_local.FEM(FEM, transf=properties, locations = pointsi, element_index = i)
#We store the results
locations[((i-1)*nrow(pointsi)+1):(i*nrow(pointsi)),] = cbind(pointsi,z)
triangles[((i-1)*nrow(meshi$triangles)+1):(i*nrow(meshi$triangles)),] = meshi$triangles+tot
tot = tot + nrow(meshi$nodes)
}
nsurf = dim(coeff)[[2]]
for (isurf in 1:nsurf)
{
open3d()
axes3d()
pop3d("lights")
light3d(specular="black")
z = locations[as.vector(t(triangles)), 2 + isurf]
triangles3d(x = locations[as.vector(t(triangles)) ,1], y = locations[as.vector(t(triangles)) ,2],
z = z,
color = heat[round(99*(z-min(z))/(max(z)-min(z)))+1],...)
aspect3d(2,2,1)
view3d(0,-45)
if (nsurf > 1)
{readline("Press a button for the next plot...")}
}
}
R_plot_manifold = function(FEM, ...){
nodes <- FEM$FEMbasis$mesh$nodes
triangles <- as.vector(t(FEM$FEMbasis$mesh$triangles))
edges <- as.vector(t(FEM$FEMbasis$mesh$edges))
coeff <- FEM$coeff
p <- jet.col(n=128,alpha=0.8)
palette(p)
ncolor <- length(p)
nplots <- ncol(coeff)
for (i in 1:nplots){
if (i > 1)
readline("Press any key for the next plot...")
open3d()
axes3d()
pop3d("lights")
light3d(specular="black")
col <- coeff[triangles,i]
col <- (ncolor-1)*(col-min(col))/diff(range(col))+1
col <- p[col]
triangles3d(nodes[triangles,1], nodes[triangles,2],
nodes[triangles,3], color = col,...)
segments3d(nodes[edges,1], nodes[edges,2], nodes[edges,3],
color = "black",...)
aspect3d("iso")
view3d(0,-45)
}
}
R_plot_volume = function(FEM,...){
nodes <- FEM$FEMbasis$mesh$nodes
faces <- as.vector(t(FEM$FEMbasis$mesh$faces[as.logical(FEM$FEMbasis$mesh$facesmarkers),]))
# edges <- as.vector(t(FEM$FEMbasis$mesh$edges[as.logical(FEM$FEMbasis$mesh$edgesmarkers),]))
coeff <- FEM$coeff
p <- jet.col(n=128,alpha=0.8)
palette(p)
ncolor <- length(p)
nplots <- ncol(FEM$coeff)
for (i in 1:nplots){
if (i > 1)
readline("Press any key for the next plot...")
open3d()
axes3d()
pop3d("lights")
light3d(specular="black")
col <- coeff[faces,i]
col <- (ncolor-1)*(col-min(col))/diff(range(col))+1
col <- p[col]
triangles3d(nodes[faces,1], nodes[faces,2],
nodes[faces,3], color = col,...)
# segments3d(nodes[edges,1], nodes[edges,2], nodes[edges,3],
# color = "black",...)
aspect3d("iso")
view3d(0,-45)
}
}
R_plot_graph = function(FEM, ...){
nodes <- FEM$FEMbasis$mesh$nodes
if(FEM$FEMbasis$order==1){
edges <- as.vector(t(FEM$FEMbasis$mesh$edges))
}else{
edges <- cbind(FEM$FEMbasis$mesh$edges[,1],
FEM$FEMbasis$mesh$edges[,3],
FEM$FEMbasis$mesh$edges[,3],
FEM$FEMbasis$mesh$edges[,2])
edges <- as.vector(t(edges))
}
coeff <- FEM$coeff
p <- jet.col(n=128,alpha=0.8)
palette(p)
ncolor <- length(p)
nplots <- ncol(coeff)
for (i in 1:nplots){
if (i > 1)
readline("Press any key for the next plot...")
open3d()
pop3d("lights")
light3d(specular="black")
col <- coeff[edges,i]
col <- (ncolor-1)*(col-min(col))/diff(range(col))+1
col <- p[col]
segments3d(nodes[edges,1], nodes[edges,2],rep(0,dim(nodes)[1]),
color = col,lwd=2.5,...)
view3d(0,0,zoom=0.75)
}
}
#' Image Plot of a 2D FEM object
#'
#' @param x A 2D-mesh \code{FEM} object.
#' @param num_refinements A natural number specifying how many bisections should by applied to each triangular element for
#' plotting purposes. This functionality is useful where a discretization with 2nd order Finite Element is applied.
#' @param ... Arguments representing graphical options to be passed to \link[rgl]{plot3d}.
#' @description Image plot of a \code{FEM} object, generated by the function \code{FEM} or returned by
#' \code{smooth.FEM} and \code{FPCA.FEM}. Only FEM objects defined over a 2D mesh can be plotted with this method.
#' @usage \method{image}{FEM}(x, num_refinements, ...)
#' @method image FEM
#' @seealso \code{\link{FEM}} \code{\link{plot.FEM}}
#' @export
#' @return No return value
#' @examples
#' library(fdaPDE)
#' ## Upload the horseshoe2D data
#' data(horseshoe2D)
#' boundary_nodes = horseshoe2D$boundary_nodes
#' boundary_segments = horseshoe2D$boundary_segments
#' locations = horseshoe2D$locations
#'
#' ## Create the 2D mesh
#' mesh = create.mesh.2D(nodes = rbind(boundary_nodes, locations), segments = boundary_segments)
#' ## Create the FEM basis
#' FEMbasis = create.FEM.basis(mesh)
#' ## Compute the coeff vector evaluating the desired function at the mesh nodes
#' ## In this case we consider the fs.test() function introduced by Wood et al. 2008
#' coeff = fs.test(mesh$nodes[,1], mesh$nodes[,2])
#' ## Create the FEM object
#' FEMfunction = FEM(coeff, FEMbasis)
#'
#' ## Plot the FEM function
#' image(FEMfunction)
image.FEM = function(x, num_refinements = NULL, ...)
{
if(!is(x$FEMbasis$mesh, "mesh.2D"))
stop('This function is implemented only for 2D mesh FEM objects')
if(x$FEMbasis$order == 1)
{
R_image.ORD1.FEM(x, ...)
}else{
R_image.ORDN.FEM(x, num_refinements, ...)
}
}
#' Image plot of a 2D FEM.time object at a given time
#'
#' @param x A 2D-mesh \code{FEM.time} object.
#' @param t time at which do the plot
#' @param lambdaS index of the space penalization parameter to use for the plot, useful when \code{FEM.time} returned by \code{smooth.FEM.time} using GCV
#' @param lambdaT index of the time penalization parameter to use for the plot, useful when \code{FEM.time} returned by \code{smooth.FEM.time} using GCV
#' @param num_refinements A natural number specifying how many bisections should by applied to each triangular element for
#' plotting purposes. This functionality is useful where a discretization with 2nd order Finite Element is applied.
#' @param ... Arguments representing graphical options to be passed to \link[rgl]{plot3d}.
#' @description Image plot of a \code{FEM.time} object, generated by the function \code{FEM.time} or returned by
#' \code{smooth.FEM.time}. Only FEM objects defined over a 2D mesh can be plotted with this method.
#' @usage \method{image}{FEM.time}(x,t,lambdaS=1,lambdaT=1,num_refinements=NULL,...)
#' @method image FEM.time
#' @seealso \code{\link{FEM.time}} \code{\link{image.FEM.time}}
#' @return No return value
#' @export
#' @examples
#' library(fdaPDE)
#' ## Upload the horseshoe2D data
#' data(horseshoe2D)
#' boundary_nodes = horseshoe2D$boundary_nodes
#' boundary_segments = horseshoe2D$boundary_segments
#' locations = horseshoe2D$locations
#'
#' ## Create the 2D mesh
#' mesh = create.mesh.2D(nodes = rbind(boundary_nodes, locations), segments = boundary_segments)
#' ## Create the FEM basis
#' FEMbasis = create.FEM.basis(mesh)
#' ## Compute the coeff vector evaluating the desired function at the mesh nodes
#' ## In this case we consider the fs.test() function introduced by Wood et al. 2008
#' time = 1:5
#' coeff = rep(fs.test(mesh$nodes[,1], mesh$nodes[,2]),5)*time
#' ## Create the FEM.time object
#' FEM_time_function = FEM.time(coeff=coeff, time_mesh=1:5,FEMbasis=FEMbasis,FLAG_PARABOLIC=TRUE)
#'
#' ## Plot the FEM function
#' t = c(1.2,1.5,3.6,2.4,4.5)
#' image(FEM_time_function,t)
image.FEM.time = function(x,t,lambdaS=1,lambdaT=1,num_refinements=NULL,...)
{
t <- as.matrix(t)
if(!is(x, "FEM.time"))
stop("x is not of class 'FEM.time'")
if(is.null(t))
stop("time required; is NULL")
if(ncol(t)>1)
stop("t must be a column vector")
if(min(t)<x$mesh_time[1] || max(t)>x$mesh_time[length(x$mesh_time)])
stop("time provided out of the 'time_mesh'")
if(dim(x$coeff)[2]>1 && lambdaS==1)
warning("the first value of lambdaS is being used")
if(dim(x$coeff)[3]>1 && lambdaT==1)
warning("the first value of lambdaT is being used")
if(is(x$FEMbasis$mesh, "mesh.2D"))
N = nrow(x$FEMbasis$mesh$nodes)
else
stop("x$FEMbasis$mesh is not of class 'mesh.2D'")
t <- as.vector(t)
storage.mode(t) <- "double"
storage.mode(N) <- "integer"
storage.mode(x$mesh_time) <- "double"
storage.mode(x$coeff) <- "double"
storage.mode(x$FLAG_PARABOLIC) <- "integer"
solution <- .Call("eval_FEM_time_nodes",N,x$mesh_time,t,x$coeff[,lambdaS,lambdaT],x$FLAG_PARABOLIC, PACKAGE = "fdaPDE")
for(i in 1:length(t))
{
plot = FEM(solution[(1+(i-1)*N):(N+(i-1)*N)],x$FEMbasis)
image.FEM(plot,num_refinements,...)
}
}
R_image.ORD1.FEM = function(FEM)
{
# PLOT Plots a FEM object FDOBJ over a rectangular grid defined by
# vectors X and Y;
#
nodes <- FEM$FEMbasis$mesh$nodes
triangles <- as.vector(t(FEM$FEMbasis$mesh$triangles))
heat <- heat.colors(100)
# How many plots are needed?
nplots <- ncol(FEM$coeff)
for (i in 1:nplots) {
if (i > 1)
readline("Press any key for the next plot...")
open3d()
axes3d()
pop3d("lights")
light3d(specular="black")
z <- FEM$coeff[triangles,i]
triangles3d(nodes[triangles,1], nodes[triangles,2], 0,
color = heat[round(99*(z-min(z))/diff(range(z)))+1])
aspect3d(2,2,1)
view3d(0,0)
}
}
R_image.ORDN.FEM = function(FEM, num_refinements)
{
coeff = FEM$coeff
FEMbasis = FEM$FEMbasis
mesh = FEMbasis$mesh
heat = heat.colors(100)
coeff = FEM$coeff
if(is.null(num_refinements))
{
num_refinements = 10
}
x = seq(from = 0, to = 1, length.out = num_refinements+1)
y = seq(from = 0, to = 1, length.out = num_refinements+1)
points_ref = expand.grid(x,y)
points_ref = points_ref[points_ref[,1] + points_ref[,2] <= 1,]
meshi = create.mesh.2D(nodes = points_ref, order = 1)
#plot(meshi)
locations = matrix(nrow = nrow(mesh$triangles)*nrow(meshi$nodes), ncol = 3*ncol(coeff))
triangles = matrix(nrow = nrow(mesh$triangles)*nrow(meshi$triangles), ncol = 3*ncol(coeff))
tot = 0
properties<-R_elementProperties(mesh)
for (i in 1:nrow(mesh$triangles))
{
# For each traingle we define a fine mesh as the transofrmation of the one constructed for the reference
transf<-rbind(cbind(properties$transf_coord$diff1x[i],properties$transf_coord$diff2x[i]),c(properties$transf_coord$diff1y[i],properties$transf_coord$diff2y[i]))
pointsi = t(transf%*%t(meshi$nodes) + mesh$nodes[mesh$triangles[i,1],])
#We evaluate the fine mesh OBS: we know the triangle we are working on no need for point location
z = R_eval_local.FEM(FEM, transf=properties, locations = pointsi, element_index = i)
#We store the results
locations[((i-1)*nrow(pointsi)+1):(i*nrow(pointsi)),] = cbind(pointsi,z)
triangles[((i-1)*nrow(meshi$triangles)+1):(i*nrow(meshi$triangles)),] = meshi$triangles+tot
tot = tot + nrow(meshi$nodes)
}
heat = heat.colors(100)
nsurf = dim(coeff)[[2]]
for (isurf in 1:nsurf)
{
open3d()
axes3d()
pop3d("lights")
light3d(specular="black")
z = locations[as.vector(t(triangles)), 2 + isurf];
triangles3d(x = locations[as.vector(t(triangles)) ,1], y = locations[as.vector(t(triangles)) ,2],
z=0,
color = heat[round(99*(z- min(z))/(max(z)-min(z)))+1])
aspect3d(2,2,1)
view3d(0,0)
if (nsurf > 1)
{readline("Press a button for the next plot...")}
}
}
R_eval_local.FEM = function(FEM, transf, locations, element_index)
{
N = nrow(locations)
# Augment Xvec and Yvec by ones for computing barycentric coordinates
Pgpts = cbind(matrix(1,N,1),locations[,1],locations[,2])
# Get nodes and index
FEMbasis = FEM$FEMbasis
mesh = FEMbasis$mesh
nodes = mesh$nodes
triangles = mesh$triangles
coeff = FEM$coeff
nsurf = dim(coeff)[2]
order = FEMbasis$order
#nodeindex = params$nodeindex
detJ = transf$detJ
# 1st, 2nd, 3rd vertices of triangles
v1 = nodes[triangles[element_index,1],]
v2 = nodes[triangles[element_index,2],]
v3 = nodes[triangles[element_index,3],]
if(order !=2 && order != 1)
{
stop('ORDER is neither 1 or 2.')
}
# Denominator of change of coordinates chsange matrix
modJ = transf$detJ[element_index]
ones3 = matrix(1,3,1)
#modJMat = modJ %*% t(ones3)
M1 = c(v2[1]*v3[2] - v3[1]*v2[2], v2[2] - v3[2], v3[1] - v2[1])/(modJ)
M2 = c(v3[1]*v1[2] - v1[1]*v3[2], v3[2] - v1[2], v1[1] - v3[1])/(modJ)
M3 = c(v1[1]*v2[2] - v2[1]*v1[2], v1[2] - v2[2], v2[1] - v1[1])/(modJ)
evalmat = matrix(NA, nrow=N, ncol=nsurf)
for (isurf in 1:nsurf)
{
for(i in 1:N)
{
baryc1 = (M1*Pgpts[i,]) %*% ones3
baryc2 = (M2*Pgpts[i,]) %*% ones3
baryc3 = (M3*Pgpts[i,]) %*% ones3
if(order == 2)
{
c1 = coeff[triangles[element_index,1],isurf]
c2 = coeff[triangles[element_index,2],isurf]
c3 = coeff[triangles[element_index,3],isurf]
c4 = coeff[triangles[element_index,6],isurf]
c5 = coeff[triangles[element_index,4],isurf]
c6 = coeff[triangles[element_index,5],isurf]
fval = c1*(2* baryc1^2 - baryc1) +
c2*(2* baryc2^2 - baryc2) +
c3*(2* baryc3^2 - baryc3) +
c4*(4* baryc1 * baryc2) +
c5*(4* baryc2 * baryc3) +
c6*(4* baryc3 * baryc1)
evalmat[i,isurf] = fval
}else{
c1 = coeff[triangles[element_index,1],isurf]
c2 = coeff[triangles[element_index,2],isurf]
c3 = coeff[triangles[element_index,3],isurf]
fval = c1*baryc1 + c2*baryc2 + c3*baryc3
evalmat[i,isurf] = fval
}
}
}
return(evalmat)
}
R_elementProperties=function(mesh)
{
nele = dim(mesh$triangles)[[1]]
nodes = mesh$nodes
triangles = mesh$triangles
#detJ = matrix(0,nele,1) # vector of determinant of transformations
#metric = array(0,c(nele,2,2)) # 3-d array of metric matrices
#transf = array(0,c(nele,2,2))
transf_coord = NULL
transf_coord$diff1x = nodes[triangles[,2],1] - nodes[triangles[,1],1]
transf_coord$diff1y = nodes[triangles[,2],2] - nodes[triangles[,1],2]
transf_coord$diff2x = nodes[triangles[,3],1] - nodes[triangles[,1],1]
transf_coord$diff2y = nodes[triangles[,3],2] - nodes[triangles[,1],2]
# Jacobian or double of the area of triangle
detJ = transf_coord$diff1x*transf_coord$diff2y - transf_coord$diff2x*transf_coord$diff1y
#Too slow, computed only for stiff from diff1x,diff1y,..
# for (i in 1:nele)
# {
# #transf[i,,] = rbind(cbind(diff1x,diff2x),c(diff1y,diff2y))
# # Compute controvariant transformation matrix OSS: This is (tranf)^(-T)
# Ael = matrix(c(diff2y, -diff1y, -diff2x, diff1x),nrow=2,ncol=2,byrow=T)/detJ[i]
#
# # Compute metric matrix
# metric[i,,] = t(Ael)%*%Ael
# }
#FEStruct <- list(detJ=detJ, metric=metric, transf=transf)
FEStruct <- list(detJ=detJ, transf_coord=transf_coord)
return(FEStruct)
}
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 R_elementProperties R_eval_local.FEM R_image.ORDN.FEM R_image.ORD1.FEM image.FEM.time image.FEM R_plot_graph R_plot_volume R_plot_manifold R_plot.ORDN.FEM R_plot.ORD1.FEM plot.mesh.1.5D plot.mesh.3D plot.mesh.2.5D plot.mesh.2D plot.FEM.time plot.FEM
Documented in image.FEM image.FEM.time plot.FEM plot.FEM.time plot.mesh.1.5D plot.mesh.2.5D plot.mesh.2D plot.mesh.3D
#' Plot a \code{FEM} object
#'
#' @param x A \code{FEM} object.
#' @param colormap A colormap exploited in the plot. The default value is the heat colormap.
#' @param num_refinements A natural number specifying how many bisections should be applied to each triangular element for
#' plotting purposes. This functionality is useful where a discretization with 2nd order Finite Element is applied.
#' This parameter can be specified only when a FEM object defined over a 2D mesh is plotted.
#' @param ... Arguments representing graphical options to be passed to \link[rgl]{plot3d}.
#' @description Three-dimensional plot of a \code{FEM} object, generated by \code{FEM} or returned by
#' \code{smooth.FEM} or \code{FPCA.FEM}.
#' If the \code{mesh} of the \code{FEMbasis} component is of class \code{mesh.2D} both the 3rd axis and the color represent
#' the value of the coefficients for the Finite Element basis expansion (\code{coeff} component of the \code{FEM} object).
#' If the \code{mesh} is of class \code{mesh.3D}, the color of each triangle or tetrahedron represent the mean value of
#' the coefficients for the Finite Element basis expansion (\code{coeff}).
#' @usage \method{plot}{FEM}(x, colormap = "heat.colors", num_refinements = NULL, ...)
#' @export
#' @seealso \code{\link{FEM}}, \code{\link{image.FEM}}
#' @return No return value
#' @examples
#' library(fdaPDE)
#' ## Upload the horseshoe2D data
#' data(horseshoe2D)
#' boundary_nodes = horseshoe2D$boundary_nodes
#' boundary_segments = horseshoe2D$boundary_segments
#' locations = horseshoe2D$locations
#'
#' ## Create the 2D mesh
#' mesh = create.mesh.2D(nodes = rbind(boundary_nodes, locations), segments = boundary_segments)
#' ## Create the FEM basis
#' FEMbasis = create.FEM.basis(mesh)
#' ## Compute the coeff vector evaluating the desired function at the mesh nodes
#' ## In this case we consider the fs.test() function introduced by Wood et al. 2008
#' coeff = fs.test(mesh$nodes[,1], mesh$nodes[,2])
#' ## Create the FEM object
#' FEMfunction = FEM(coeff, FEMbasis)
#'
#' ## Plot the FEM function
#' plot(FEMfunction)
plot.FEM = function(x, colormap = "heat.colors", num_refinements = NULL, ...)
{
if(is(x$FEMbasis$mesh, "mesh.2D")){
if(x$FEMbasis$order == 1)
{
R_plot.ORD1.FEM(x, colormap, ...)
}else{
R_plot.ORDN.FEM(x, colormap, num_refinements, ...)
}
}else if(is(x$FEMbasis$mesh, "mesh.2.5D")){
R_plot_manifold(x,...)
}else if(is(x$FEMbasis$mesh, "mesh.3D")){
R_plot_volume(x,...)
}else if(is(x$FEMbasis$mesh, "mesh.1.5D")){
R_plot_graph(x, ...)
}
}
#' Plot a \code{FEM.time} object at a given time
#'
#' @param x A \code{FEM.time} object.
#' @param time_locations A vector with the instants in which plot the spatial field
#' @param locations A 2-column matrix when \code{x$FEMbasis$mesh} is of class \code{mesh.2D}
#' or a 3-column matrix otherwise with the spatial locations for which plot the temporal evolution
#' @param lambdaS Index of the space penalization parameter to use for the plot, useful when \code{FEM.time} returned by \code{smooth.FEM.time} using GCV
#' @param lambdaT Index of the time penalization parameter to use for the plot, useful when \code{FEM.time} returned by \code{smooth.FEM.time} using GCV
#' @param num_refinements A natural number specifying how many bisections should be applied to each triangular element for
#' plotting purposes. This functionality is useful where a discretization with 2nd order Finite Element is applied.
#' This parameter can be specified only when a FEM object defined over a 2D mesh is plotted.
#' @param Nt The number of instants to plot when \code{locations} is not NULL
#' @param add Boolean, used only when locations is not NULL, if TURE it performs the graphic over the last plot
#' @param main The title of the plot when \code{locations} is not NULL
#' @param col The color of the plot when \code{locations} is not NULL. May be a single color or a vector of colors
#' @param ... Arguments representing graphical options to be passed to \link[rgl]{plot3d}.
#' @description Plot of a \code{FEM.time} object, generated by \code{FEM.time} or returned by
#' \code{smooth.FEM.time}. \code{time_locations} and \code{locations} must not be both provided.
#' If \code{time_locations} is provided, the spatial field is plotted for the provided temporal instnts.
#' If \code{locations} is provided, the temporal evolution in the provided space locations is plotted.
#' If both \code{time_locations} and \code{locations} are NULL a default plot is provided.
#' If the \code{mesh} of the \code{FEMbasis} component is of class \code{mesh.2D} both the 3rd axis and the color represent
#' the value of the coefficients for the Finite Element basis expansion (\code{coeff} component of the \code{FEM.time} object).
#' If the \code{mesh} is of class \code{mesh.3D}, the color of each triangle or tetrahedron represent the mean value of
#' the coefficients for the Finite Element basis expansion (\code{coeff}).
#' @usage \method{plot}{FEM.time}(x, time_locations = NULL, locations = NULL,
#' lambdaS = NULL, lambdaT = NULL, num_refinements = NULL, Nt = 100,
#' add = FALSE, main = NULL, col = "red", ...)
#' @export
#' @seealso \code{\link{FEM.time}}, \code{\link{image.FEM.time}}
#' @return No return value
#' @examples
#' library(fdaPDE)
#' ## Upload the horseshoe2D data
#' data(horseshoe2D)
#' boundary_nodes = horseshoe2D$boundary_nodes
#' boundary_segments = horseshoe2D$boundary_segments
#' locations = horseshoe2D$locations
#'
#' ## Create the 2D mesh
#' mesh = create.mesh.2D(nodes = rbind(boundary_nodes, locations), segments = boundary_segments)
#' ## Create the FEM basis
#' FEMbasis = create.FEM.basis(mesh)
#' ## Compute the coeff vector evaluating the desired function at the mesh nodes
#' ## In this case we consider the fs.test() function introduced by Wood et al. 2008
#' time = 1:5
#' coeff = rep(fs.test(mesh$nodes[,1], mesh$nodes[,2]),5)*time
#' ## Create the FEM.time object
#' FEM_time_function = FEM.time(coeff=coeff, time_mesh=1:5, FEMbasis=FEMbasis, FLAG_PARABOLIC=TRUE)
#'
#' ## Plot the FEM function
#' plot(FEM_time_function)
#'
#' ## plot spatial field in some instants
#' t = c(1.2,1.5,3.6,2.4,4.5)
#' plot(FEM_time_function, t)
#'
#' ## plot time evolution in some locations
#' plot(FEM_time_function, locations = locations[1:10,])
plot.FEM.time = function(x, time_locations = NULL, locations = NULL,
lambdaS = NULL, lambdaT = NULL, num_refinements = NULL, Nt = 100,
add = FALSE, main = NULL, col = "red", ...)
{
if(!is(x, "FEM.time"))
stop("x is not of class 'FEM.time'")
if (!is.null(time_locations) && !is.null(locations))
stop("time_locations and locations can not be both provided")
if (is.null(time_locations) && is.null(locations)) {
# default plot
# space evolution:
time_locations = mean(x$mesh_time)
if(is.null(lambdaS) || is.null(lambdaT)){
if(exists("x$bestlambda")){
lambdaS = x$bestlambda[1]
lambdaT = x$bestlambda[2]
}else{
lambdaS = 1
lambdaT = 1
}
}
N = nrow(x$FEMbasis$mesh$nodes)
storage.mode(time_locations) <- "double"
storage.mode(N) <- "integer"
storage.mode(x$mesh_time) <- "double"
storage.mode(x$coeff) <- "double"
storage.mode(x$FLAG_PARABOLIC) <- "integer"
solution <- .Call("eval_FEM_time_nodes",N,x$mesh_time,time_locations,x$coeff[,lambdaS,lambdaT],x$FLAG_PARABOLIC, PACKAGE = "fdaPDE")
plot = FEM(solution[1:N],x$FEMbasis)
plot.FEM(x = plot, num_refinements = num_refinements,main = paste0("Spatial function in t = ",time_locations),...)
# time evolution:
locations = x$FEMbasis$mesh$nodes[1,]
# add default title to the plot
if(is.null(main))
main = paste0("Time evolution in (",round(locations,2),")")
# define the evaluation points
t = seq(min(x$mesh_time),max(x$mesh_time),length.out=Nt)
eval_points = cbind(t, t(kronecker(matrix(1,ncol=Nt), locations)))
# evaluate the function
eval_sol = eval.FEM.time(FEM.time = x, space.time.locations = eval_points, lambdaT=lambdaT, lambdaS=lambdaS)
# define variables for the plot
type = "l"
ylab = ""
xlab = "Time"
ylim = range(eval_sol)
# plot the function evaluations
if (add == FALSE) {
plot(t[1:Nt],eval_sol[1:Nt],type=type,xlab=xlab,ylab=ylab,main=main,col=col,ylim=ylim, ...)
}else{
points(t[1:Nt],eval_sol[1:Nt],type=type,col=col, ...)
}
}
if (!is.null(time_locations) && is.null(locations)) {
# plot the spatial function for the provided time_locations
# check that time_locations is a vector
if (!is.vector(time_locations) && nrow(time_locations)!=1 && ncol(time_locations)!=1)
stop("time_locations must be a vector of a column/row matrix")
time_locations = as.vector(time_locations)
# check that time_locations is in the correct range
if(min(time_locations)<x$mesh_time[1] || max(time_locations)>x$mesh_time[length(x$mesh_time)])
stop("time_locations must be within the 'time_mesh' range")
if(is.null(lambdaS) || is.null(lambdaT)){
if(exists("x$bestlambda")){
lambdaS = x$bestlambda[1]
lambdaT = x$bestlambda[2]
}else{
lambdaS = 1
lambdaT = 1
}
}
N = nrow(x$FEMbasis$mesh$nodes)
storage.mode(time_locations) <- "double"
storage.mode(N) <- "integer"
storage.mode(x$mesh_time) <- "double"
storage.mode(x$coeff) <- "double"
storage.mode(x$FLAG_PARABOLIC) <- "integer"
solution <- .Call("eval_FEM_time_nodes",N,x$mesh_time,time_locations,x$coeff[,lambdaS,lambdaT],x$FLAG_PARABOLIC, PACKAGE = "fdaPDE")
for(i in 1:length(t))
{
plot = FEM(solution[(1+(i-1)*N):(N+(i-1)*N)],x$FEMbasis)
plot.FEM(x = plot, num_refinements = num_refinements,...)
}
}
if (is.null(time_locations) && !is.null(locations)) {
# plot temporal evolution in the provided locations
if (is(x$FEMbasis$mesh, "mesh.2D") | is(x$FEMbasis$mesh, "mesh.1.5D")){
# check locations dimension
if (is.vector(locations) && length(locations) != 2)
stop("locations does not have the correct dimension")
if (!is.vector(locations) && ncol(locations) != 2)
stop("locations does not have the correct dimension")
if (is.vector(locations))
locations = matrix(locations, ncol = 2)
# add default title to the plot
if(is.null(main) && nrow(locations)==1)
main = paste0("Time evolution in (",locations[1,1],",",locations[1,2],")")
# define the evaluation points
t = rep(seq(min(x$mesh_time),max(x$mesh_time),length.out=Nt),times=nrow(locations))
eval_points = cbind(t, rep(locations[,1], each = Nt), rep(locations[,2], each = Nt))
}
if (is(x$FEMbasis$mesh, "mesh.2.5D") | is(x$FEMbasis$mesh, "mesh.3D")){
# check locations dimension
if (is.vector(locations) && length(locations) != 3)
stop("locations does not have the correct dimension")
if (!is.vector(locations) && ncol(locations) != 3)
stop("locations does not have the correct dimension")
if (is.vector(locations))
locations = matrix(locations, ncol = 3)
# add default title to the plot
if(is.null(main) && nrow(locations)==1)
main = paste0("Time evolution in (",locations[1,1],",",locations[1,2],",",locations[1,3],")")
# define the evaluation points
t = rep(seq(min(x$mesh_time),max(x$mesh_time),length.out=Nt),times=nrow(locations))
eval_points = cbind(t, rep(locations[,1], each = Nt),
rep(locations[,2], each = Nt), rep(locations[,3], each = Nt))
}
# evaluate the function
eval_sol = eval.FEM.time(FEM.time = x, space.time.locations = eval_points, lambdaT=lambdaT, lambdaS=lambdaS)
# define variables for the plot
if (length(col) < nrow(locations)){
# recicle colors if not provided enough
col = rep(col, nrow(locations))
}
if (!exists("type")) {
type = "l"
}
if (!exists("ylab")) {
ylab = ""
}
if (!exists("xlab")) {
xlab = "Time"
}
if (!exists("ylim")) {
ylim = range(eval_sol)
}
# plot the function evaluations
for (ind in 1:nrow(locations)) {
if (add == FALSE && ind == 1) {
plot(t[1:Nt],eval_sol[(ind-1)*Nt+(1:Nt)],type=type,xlab=xlab,ylab=ylab,main=main,col=col[ind],ylim=ylim, ...)
}else{
points(t[1:Nt],eval_sol[(ind-1)*Nt+(1:Nt)],type=type,col=col[ind], ...)
}
}
}
}
#' Plot a mesh.2D object
#'
#' @param x A \code{mesh.2D} object defining the triangular mesh, as generated by \code{create.mesh.2D}
#' or \code{refine.mesh.2D}.
#' @param ... Arguments representing graphical options to be passed to \link[graphics]{par}.
#' @description Plot a mesh.2D object, generated by \code{create.mesh.2D} or \code{refine.mesh.2D}.
#' @name plot.mesh.2D
#' @usage \method{plot}{mesh.2D}(x, ...)
#' @export
#' @return No return value
#' @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
#' mesh = create.mesh.2D(nodes = rbind(boundary_nodes, locations), segments = boundary_segments)
#'
#' ## Plot the mesh
#' plot(mesh)
plot.mesh.2D<-function(x, ...)
{
plot(x$nodes, xlab="", ylab="", xaxt="n", yaxt="n", bty="n", ...)
segments(x$nodes[x$edges[,1],1], x$nodes[x$edges[,1],2],
x$nodes[x$edges[,2],1], x$nodes[x$edges[,2],2], ...)
segments(x$nodes[x$segments[,1],1], x$nodes[x$segments[,1],2],
x$nodes[x$segments[,2],1], x$nodes[x$segments[,2],2], col="red", ...)
}
#' Plot a mesh.2.5D object
#'
#' @param x A \code{mesh.2.5D} object generated by \code{create.mesh.2.5D}.
#' @param ... Arguments representing graphical options to be passed to \link[graphics]{par}.
#' @description Plot the triangulation of a \code{mesh.2.5D} object, generated by \code{create.mesh.2.5D}
#' @export
#' @name plot.mesh.2.5D
#' @usage \method{plot}{mesh.2.5D}(x, ...)
#' @return No return value
#' @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)
#' plot(mesh)
plot.mesh.2.5D<-function(x,...){
nodes <- x$nodes
edges <- as.vector(t(x$edges))
open3d()
axes3d()
pop3d("lights")
light3d(specular="black")
points3d(nodes[,1], nodes[,2], nodes[,3], col="black", ...)
segments3d(nodes[edges,1], nodes[edges,2], nodes[edges,3], col="black",...)
aspect3d("iso")
view3d(0,-45)
}
#' Plot a mesh.3D object
#'
#' @param x A \code{mesh.3D} object generated by \code{create.mesh.3D}.
#' @param ... Arguments representing graphical options to be passed to \link[graphics]{par}.
#' @description Plot a \code{mesh.3D} object, generated by \code{create.mesh.3D}.
#' @export
#' @name plot.mesh.3D
#' @usage \method{plot}{mesh.3D}(x, ...)
#' @return No return value
#' @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)
#'
#' ##Plot the triangulation of the object
#' plot(mesh)
plot.mesh.3D<-function(x,...){
nodes <- x$nodes
faces <- as.vector(t(x$faces[x$facesmarkers,]))
aux_mesh <- create.mesh.2.5D(nodes=nodes, triangles=x$faces[x$facesmarkers,], order=1)
edges <- as.vector(t(aux_mesh$edges))
open3d()
axes3d()
pop3d("lights")
light3d(specular="black")
triangles3d(nodes[faces,1],nodes[faces,2],nodes[faces,3],col="white",...)
points3d(nodes[,1], nodes[,2], nodes[,3], col="black", ...)
segments3d(nodes[edges,1], nodes[edges,2], nodes[edges,3], col="black",...)
aspect3d("iso")
view3d(0,-45)
}
#' Plot a mesh.1.5D object
#'
#' @param x A \code{mesh.1.5D} object defining the triangular mesh, as generated by \code{create.mesh.1.5D}
#' or \code{refine.mesh.1.5D}.
#' @param ... Arguments representing graphical options to be passed to \link[graphics]{par}.
#' @description Plot a mesh.1.5D object, generated by \code{create.mesh.1.5D} or \code{refine.mesh.1.5D}.
#' @name plot.mesh.1.5D
#' @usage \method{plot}{mesh.1.5D}(x, ...)
#' @return No return value
#' @export
plot.mesh.1.5D<-function(x, ...)
{
if( x$order == 1 ){
plot(x$nodes, xlab="", ylab="", xaxt="n", yaxt="n", bty="n", ...)
segments(x$nodes[x$edges[,1],1], x$nodes[x$edges[,1],2],
x$nodes[x$edges[,2],1], x$nodes[x$edges[,2],2], ...)
}
else{
plot(x$nodes, xlab="", ylab="", xaxt="n", yaxt="n", bty="n", ...)
segments(x$nodes[x$edges[,1],1], x$nodes[x$edges[,1],2],
x$nodes[x$edges[,2],1], x$nodes[x$edges[,2],2], ...)
points(x$nodes[x$edges[,3],1], x$nodes[x$edges[,3],2],col="red", ...)
}
}
R_plot.ORD1.FEM = function(FEM, colormap, ...)
{
# PLOT Plots a FEM object FDOBJ over a rectangular grid defined by
# vectors X and Y;
#
nodes <- FEM$FEMbasis$mesh$nodes
triangles <- as.vector(t(FEM$FEMbasis$mesh$triangles))
colormap <- match.fun(colormap)
heat <- colormap(100)
# How many plots are needed?
nplots <- ncol(FEM$coeff)
for (i in 1:nplots) {
if (i > 1)
readline("Press any key for the next plot...")
open3d()
axes3d()
z <- FEM$coeff[triangles,i]
triangles3d(nodes[triangles,1], nodes[triangles,2], z,
color = heat[round(99*(z-min(z))/diff(range(z)))+1],...)
aspect3d(2,2,1)
view3d(0,-45)
}
}
R_plot.ORDN.FEM = function(FEM, colormap, num_refinements, ...)
{
# num_refinements sets the number od division on each triangle edge to be applied for rifenment
coeff = FEM$coeff
FEMbasis = FEM$FEMbasis
mesh = FEMbasis$mesh
colormap <- match.fun(colormap)
heat <- colormap(100)
coeff = FEM$coeff
if(is.null(num_refinements))
{
num_refinements = 10
}
# For the reference triangles we construct a regular mesh
x = seq(from = 0, to = 1, length.out = num_refinements+1)
y = seq(from = 0, to = 1, length.out = num_refinements+1)
points_ref = expand.grid(x,y)
points_ref = points_ref[points_ref[,1] + points_ref[,2] <= 1,]
meshi = create.mesh.2D(nodes = points_ref, order = 1)
#plot(meshi)
# locations is the matrix with that will contain the coordinate of the points where the function is
# evaluated (1st and 2nd column) and the columns with the evaluation of the ith fucntion on that point
locations = matrix(nrow = nrow(mesh$triangles)*nrow(meshi$nodes), ncol = 2+ncol(coeff))
triangles = matrix(nrow = nrow(mesh$triangles)*nrow(meshi$triangles), ncol = 3)
tot = 0
properties<-R_elementProperties(mesh)
for (i in 1:nrow(mesh$triangles))
{
# For each traingle we define a fine mesh as the transofrmation of the one constructed for the reference
transf<-rbind(cbind(properties$transf_coord$diff1x[i],properties$transf_coord$diff2x[i]),c(properties$transf_coord$diff1y[i],properties$transf_coord$diff2y[i]))
pointsi = t(transf%*%t(meshi$nodes) + mesh$nodes[mesh$triangles[i,1],])
#We evaluate the fine mesh OBS: we know the triangle we are working on no need for point location
z = R_eval_local.FEM(FEM, transf=properties, locations = pointsi, element_index = i)
#We store the results
locations[((i-1)*nrow(pointsi)+1):(i*nrow(pointsi)),] = cbind(pointsi,z)
triangles[((i-1)*nrow(meshi$triangles)+1):(i*nrow(meshi$triangles)),] = meshi$triangles+tot
tot = tot + nrow(meshi$nodes)
}
nsurf = dim(coeff)[[2]]
for (isurf in 1:nsurf)
{
open3d()
axes3d()
pop3d("lights")
light3d(specular="black")
z = locations[as.vector(t(triangles)), 2 + isurf]
triangles3d(x = locations[as.vector(t(triangles)) ,1], y = locations[as.vector(t(triangles)) ,2],
z = z,
color = heat[round(99*(z-min(z))/(max(z)-min(z)))+1],...)
aspect3d(2,2,1)
view3d(0,-45)
if (nsurf > 1)
{readline("Press a button for the next plot...")}
}
}
R_plot_manifold = function(FEM, ...){
nodes <- FEM$FEMbasis$mesh$nodes
triangles <- as.vector(t(FEM$FEMbasis$mesh$triangles))
edges <- as.vector(t(FEM$FEMbasis$mesh$edges))
coeff <- FEM$coeff
p <- jet.col(n=128,alpha=0.8)
palette(p)
ncolor <- length(p)
nplots <- ncol(coeff)
for (i in 1:nplots){
if (i > 1)
readline("Press any key for the next plot...")
open3d()
axes3d()
pop3d("lights")
light3d(specular="black")
col <- coeff[triangles,i]
col <- (ncolor-1)*(col-min(col))/diff(range(col))+1
col <- p[col]
triangles3d(nodes[triangles,1], nodes[triangles,2],
nodes[triangles,3], color = col,...)
segments3d(nodes[edges,1], nodes[edges,2], nodes[edges,3],
color = "black",...)
aspect3d("iso")
view3d(0,-45)
}
}
R_plot_volume = function(FEM,...){
nodes <- FEM$FEMbasis$mesh$nodes
faces <- as.vector(t(FEM$FEMbasis$mesh$faces[as.logical(FEM$FEMbasis$mesh$facesmarkers),]))
# edges <- as.vector(t(FEM$FEMbasis$mesh$edges[as.logical(FEM$FEMbasis$mesh$edgesmarkers),]))
coeff <- FEM$coeff
p <- jet.col(n=128,alpha=0.8)
palette(p)
ncolor <- length(p)
nplots <- ncol(FEM$coeff)
for (i in 1:nplots){
if (i > 1)
readline("Press any key for the next plot...")
open3d()
axes3d()
pop3d("lights")
light3d(specular="black")
col <- coeff[faces,i]
col <- (ncolor-1)*(col-min(col))/diff(range(col))+1
col <- p[col]
triangles3d(nodes[faces,1], nodes[faces,2],
nodes[faces,3], color = col,...)
# segments3d(nodes[edges,1], nodes[edges,2], nodes[edges,3],
# color = "black",...)
aspect3d("iso")
view3d(0,-45)
}
}
R_plot_graph = function(FEM, ...){
nodes <- FEM$FEMbasis$mesh$nodes
if(FEM$FEMbasis$order==1){
edges <- as.vector(t(FEM$FEMbasis$mesh$edges))
}else{
edges <- cbind(FEM$FEMbasis$mesh$edges[,1],
FEM$FEMbasis$mesh$edges[,3],
FEM$FEMbasis$mesh$edges[,3],
FEM$FEMbasis$mesh$edges[,2])
edges <- as.vector(t(edges))
}
coeff <- FEM$coeff
p <- jet.col(n=128,alpha=0.8)
palette(p)
ncolor <- length(p)
nplots <- ncol(coeff)
for (i in 1:nplots){
if (i > 1)
readline("Press any key for the next plot...")
open3d()
pop3d("lights")
light3d(specular="black")
col <- coeff[edges,i]
col <- (ncolor-1)*(col-min(col))/diff(range(col))+1
col <- p[col]
segments3d(nodes[edges,1], nodes[edges,2],rep(0,dim(nodes)[1]),
color = col,lwd=2.5,...)
view3d(0,0,zoom=0.75)
}
}
#' Image Plot of a 2D FEM object
#'
#' @param x A 2D-mesh \code{FEM} object.
#' @param num_refinements A natural number specifying how many bisections should by applied to each triangular element for
#' plotting purposes. This functionality is useful where a discretization with 2nd order Finite Element is applied.
#' @param ... Arguments representing graphical options to be passed to \link[rgl]{plot3d}.
#' @description Image plot of a \code{FEM} object, generated by the function \code{FEM} or returned by
#' \code{smooth.FEM} and \code{FPCA.FEM}. Only FEM objects defined over a 2D mesh can be plotted with this method.
#' @usage \method{image}{FEM}(x, num_refinements, ...)
#' @method image FEM
#' @seealso \code{\link{FEM}} \code{\link{plot.FEM}}
#' @export
#' @return No return value
#' @examples
#' library(fdaPDE)
#' ## Upload the horseshoe2D data
#' data(horseshoe2D)
#' boundary_nodes = horseshoe2D$boundary_nodes
#' boundary_segments = horseshoe2D$boundary_segments
#' locations = horseshoe2D$locations
#'
#' ## Create the 2D mesh
#' mesh = create.mesh.2D(nodes = rbind(boundary_nodes, locations), segments = boundary_segments)
#' ## Create the FEM basis
#' FEMbasis = create.FEM.basis(mesh)
#' ## Compute the coeff vector evaluating the desired function at the mesh nodes
#' ## In this case we consider the fs.test() function introduced by Wood et al. 2008
#' coeff = fs.test(mesh$nodes[,1], mesh$nodes[,2])
#' ## Create the FEM object
#' FEMfunction = FEM(coeff, FEMbasis)
#'
#' ## Plot the FEM function
#' image(FEMfunction)
image.FEM = function(x, num_refinements = NULL, ...)
{
if(!is(x$FEMbasis$mesh, "mesh.2D"))
stop('This function is implemented only for 2D mesh FEM objects')
if(x$FEMbasis$order == 1)
{
R_image.ORD1.FEM(x, ...)
}else{
R_image.ORDN.FEM(x, num_refinements, ...)
}
}
#' Image plot of a 2D FEM.time object at a given time
#'
#' @param x A 2D-mesh \code{FEM.time} object.
#' @param t time at which do the plot
#' @param lambdaS index of the space penalization parameter to use for the plot, useful when \code{FEM.time} returned by \code{smooth.FEM.time} using GCV
#' @param lambdaT index of the time penalization parameter to use for the plot, useful when \code{FEM.time} returned by \code{smooth.FEM.time} using GCV
#' @param num_refinements A natural number specifying how many bisections should by applied to each triangular element for
#' plotting purposes. This functionality is useful where a discretization with 2nd order Finite Element is applied.
#' @param ... Arguments representing graphical options to be passed to \link[rgl]{plot3d}.
#' @description Image plot of a \code{FEM.time} object, generated by the function \code{FEM.time} or returned by
#' \code{smooth.FEM.time}. Only FEM objects defined over a 2D mesh can be plotted with this method.
#' @usage \method{image}{FEM.time}(x,t,lambdaS=1,lambdaT=1,num_refinements=NULL,...)
#' @method image FEM.time
#' @seealso \code{\link{FEM.time}} \code{\link{image.FEM.time}}
#' @return No return value
#' @export
#' @examples
#' library(fdaPDE)
#' ## Upload the horseshoe2D data
#' data(horseshoe2D)
#' boundary_nodes = horseshoe2D$boundary_nodes
#' boundary_segments = horseshoe2D$boundary_segments
#' locations = horseshoe2D$locations
#'
#' ## Create the 2D mesh
#' mesh = create.mesh.2D(nodes = rbind(boundary_nodes, locations), segments = boundary_segments)
#' ## Create the FEM basis
#' FEMbasis = create.FEM.basis(mesh)
#' ## Compute the coeff vector evaluating the desired function at the mesh nodes
#' ## In this case we consider the fs.test() function introduced by Wood et al. 2008
#' time = 1:5
#' coeff = rep(fs.test(mesh$nodes[,1], mesh$nodes[,2]),5)*time
#' ## Create the FEM.time object
#' FEM_time_function = FEM.time(coeff=coeff, time_mesh=1:5,FEMbasis=FEMbasis,FLAG_PARABOLIC=TRUE)
#'
#' ## Plot the FEM function
#' t = c(1.2,1.5,3.6,2.4,4.5)
#' image(FEM_time_function,t)
image.FEM.time = function(x,t,lambdaS=1,lambdaT=1,num_refinements=NULL,...)
{
t <- as.matrix(t)
if(!is(x, "FEM.time"))
stop("x is not of class 'FEM.time'")
if(is.null(t))
stop("time required; is NULL")
if(ncol(t)>1)
stop("t must be a column vector")
if(min(t)<x$mesh_time[1] || max(t)>x$mesh_time[length(x$mesh_time)])
stop("time provided out of the 'time_mesh'")
if(dim(x$coeff)[2]>1 && lambdaS==1)
warning("the first value of lambdaS is being used")
if(dim(x$coeff)[3]>1 && lambdaT==1)
warning("the first value of lambdaT is being used")
if(is(x$FEMbasis$mesh, "mesh.2D"))
N = nrow(x$FEMbasis$mesh$nodes)
else
stop("x$FEMbasis$mesh is not of class 'mesh.2D'")
t <- as.vector(t)
storage.mode(t) <- "double"
storage.mode(N) <- "integer"
storage.mode(x$mesh_time) <- "double"
storage.mode(x$coeff) <- "double"
storage.mode(x$FLAG_PARABOLIC) <- "integer"
solution <- .Call("eval_FEM_time_nodes",N,x$mesh_time,t,x$coeff[,lambdaS,lambdaT],x$FLAG_PARABOLIC, PACKAGE = "fdaPDE")
for(i in 1:length(t))
{
plot = FEM(solution[(1+(i-1)*N):(N+(i-1)*N)],x$FEMbasis)
image.FEM(plot,num_refinements,...)
}
}
R_image.ORD1.FEM = function(FEM)
{
# PLOT Plots a FEM object FDOBJ over a rectangular grid defined by
# vectors X and Y;
#
nodes <- FEM$FEMbasis$mesh$nodes
triangles <- as.vector(t(FEM$FEMbasis$mesh$triangles))
heat <- heat.colors(100)
# How many plots are needed?
nplots <- ncol(FEM$coeff)
for (i in 1:nplots) {
if (i > 1)
readline("Press any key for the next plot...")
open3d()
axes3d()
pop3d("lights")
light3d(specular="black")
z <- FEM$coeff[triangles,i]
triangles3d(nodes[triangles,1], nodes[triangles,2], 0,
color = heat[round(99*(z-min(z))/diff(range(z)))+1])
aspect3d(2,2,1)
view3d(0,0)
}
}
R_image.ORDN.FEM = function(FEM, num_refinements)
{
coeff = FEM$coeff
FEMbasis = FEM$FEMbasis
mesh = FEMbasis$mesh
heat = heat.colors(100)
coeff = FEM$coeff
if(is.null(num_refinements))
{
num_refinements = 10
}
x = seq(from = 0, to = 1, length.out = num_refinements+1)
y = seq(from = 0, to = 1, length.out = num_refinements+1)
points_ref = expand.grid(x,y)
points_ref = points_ref[points_ref[,1] + points_ref[,2] <= 1,]
meshi = create.mesh.2D(nodes = points_ref, order = 1)
#plot(meshi)
locations = matrix(nrow = nrow(mesh$triangles)*nrow(meshi$nodes), ncol = 3*ncol(coeff))
triangles = matrix(nrow = nrow(mesh$triangles)*nrow(meshi$triangles), ncol = 3*ncol(coeff))
tot = 0
properties<-R_elementProperties(mesh)
for (i in 1:nrow(mesh$triangles))
{
# For each traingle we define a fine mesh as the transofrmation of the one constructed for the reference
transf<-rbind(cbind(properties$transf_coord$diff1x[i],properties$transf_coord$diff2x[i]),c(properties$transf_coord$diff1y[i],properties$transf_coord$diff2y[i]))
pointsi = t(transf%*%t(meshi$nodes) + mesh$nodes[mesh$triangles[i,1],])
#We evaluate the fine mesh OBS: we know the triangle we are working on no need for point location
z = R_eval_local.FEM(FEM, transf=properties, locations = pointsi, element_index = i)
#We store the results
locations[((i-1)*nrow(pointsi)+1):(i*nrow(pointsi)),] = cbind(pointsi,z)
triangles[((i-1)*nrow(meshi$triangles)+1):(i*nrow(meshi$triangles)),] = meshi$triangles+tot
tot = tot + nrow(meshi$nodes)
}
heat = heat.colors(100)
nsurf = dim(coeff)[[2]]
for (isurf in 1:nsurf)
{
open3d()
axes3d()
pop3d("lights")
light3d(specular="black")
z = locations[as.vector(t(triangles)), 2 + isurf];
triangles3d(x = locations[as.vector(t(triangles)) ,1], y = locations[as.vector(t(triangles)) ,2],
z=0,
color = heat[round(99*(z- min(z))/(max(z)-min(z)))+1])
aspect3d(2,2,1)
view3d(0,0)
if (nsurf > 1)
{readline("Press a button for the next plot...")}
}
}
R_eval_local.FEM = function(FEM, transf, locations, element_index)
{
N = nrow(locations)
# Augment Xvec and Yvec by ones for computing barycentric coordinates
Pgpts = cbind(matrix(1,N,1),locations[,1],locations[,2])
# Get nodes and index
FEMbasis = FEM$FEMbasis
mesh = FEMbasis$mesh
nodes = mesh$nodes
triangles = mesh$triangles
coeff = FEM$coeff
nsurf = dim(coeff)[2]
order = FEMbasis$order
#nodeindex = params$nodeindex
detJ = transf$detJ
# 1st, 2nd, 3rd vertices of triangles
v1 = nodes[triangles[element_index,1],]
v2 = nodes[triangles[element_index,2],]
v3 = nodes[triangles[element_index,3],]
if(order !=2 && order != 1)
{
stop('ORDER is neither 1 or 2.')
}
# Denominator of change of coordinates chsange matrix
modJ = transf$detJ[element_index]
ones3 = matrix(1,3,1)
#modJMat = modJ %*% t(ones3)
M1 = c(v2[1]*v3[2] - v3[1]*v2[2], v2[2] - v3[2], v3[1] - v2[1])/(modJ)
M2 = c(v3[1]*v1[2] - v1[1]*v3[2], v3[2] - v1[2], v1[1] - v3[1])/(modJ)
M3 = c(v1[1]*v2[2] - v2[1]*v1[2], v1[2] - v2[2], v2[1] - v1[1])/(modJ)
evalmat = matrix(NA, nrow=N, ncol=nsurf)
for (isurf in 1:nsurf)
{
for(i in 1:N)
{
baryc1 = (M1*Pgpts[i,]) %*% ones3
baryc2 = (M2*Pgpts[i,]) %*% ones3
baryc3 = (M3*Pgpts[i,]) %*% ones3
if(order == 2)
{
c1 = coeff[triangles[element_index,1],isurf]
c2 = coeff[triangles[element_index,2],isurf]
c3 = coeff[triangles[element_index,3],isurf]
c4 = coeff[triangles[element_index,6],isurf]
c5 = coeff[triangles[element_index,4],isurf]
c6 = coeff[triangles[element_index,5],isurf]
fval = c1*(2* baryc1^2 - baryc1) +
c2*(2* baryc2^2 - baryc2) +
c3*(2* baryc3^2 - baryc3) +
c4*(4* baryc1 * baryc2) +
c5*(4* baryc2 * baryc3) +
c6*(4* baryc3 * baryc1)
evalmat[i,isurf] = fval
}else{
c1 = coeff[triangles[element_index,1],isurf]
c2 = coeff[triangles[element_index,2],isurf]
c3 = coeff[triangles[element_index,3],isurf]
fval = c1*baryc1 + c2*baryc2 + c3*baryc3
evalmat[i,isurf] = fval
}
}
}
return(evalmat)
}
R_elementProperties=function(mesh)
{
nele = dim(mesh$triangles)[[1]]
nodes = mesh$nodes
triangles = mesh$triangles
#detJ = matrix(0,nele,1) # vector of determinant of transformations
#metric = array(0,c(nele,2,2)) # 3-d array of metric matrices
#transf = array(0,c(nele,2,2))
transf_coord = NULL
transf_coord$diff1x = nodes[triangles[,2],1] - nodes[triangles[,1],1]
transf_coord$diff1y = nodes[triangles[,2],2] - nodes[triangles[,1],2]
transf_coord$diff2x = nodes[triangles[,3],1] - nodes[triangles[,1],1]
transf_coord$diff2y = nodes[triangles[,3],2] - nodes[triangles[,1],2]
# Jacobian or double of the area of triangle
detJ = transf_coord$diff1x*transf_coord$diff2y - transf_coord$diff2x*transf_coord$diff1y
#Too slow, computed only for stiff from diff1x,diff1y,..
# for (i in 1:nele)
# {
# #transf[i,,] = rbind(cbind(diff1x,diff2x),c(diff1y,diff2y))
# # Compute controvariant transformation matrix OSS: This is (tranf)^(-T)
# Ael = matrix(c(diff2y, -diff1y, -diff2x, diff1x),nrow=2,ncol=2,byrow=T)/detJ[i]
#
# # Compute metric matrix
# metric[i,,] = t(Ael)%*%Ael
# }
#FEStruct <- list(detJ=detJ, metric=metric, transf=transf)
FEStruct <- list(detJ=detJ, transf_coord=transf_coord)
return(FEStruct)
}
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.