R/plot.R

Defines functions R_eval_local.FEM R_image.ORDN.FEM R_image.ORD1.FEM image.FEM R_plot_volume R_plot_manifold R_plot.ORDN.FEM R_plot.ORD1.FEM plot.mesh.3D plot.mesh.2.5D plot.mesh.2D plot.FEM

Documented in image.FEM plot.FEM plot.mesh.2.5D plot.mesh.2D

#' Plot a \code{FEM} object
#' 
#' @param x A \code{FEM} object.
#' @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). 
#' @usage \method{plot}{FEM}(x, num_refinements, ...)  
#' @export
#' @seealso \code{\link{FEM}}, \code{\link{image.FEM}}
#' @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, num_refinements = NULL, ...)  
{
if(class(x$FEMbasis$mesh)=="mesh.2D"){
  if(x$FEMbasis$order == 1)
  {
    R_plot.ORD1.FEM(x, ...)
  }else{
    R_plot.ORDN.FEM(x, num_refinements, ...)
  }
}else if(class(x$FEMbasis$mesh)=="mesh.2.5D"){
	R_plot_manifold(x,...)
}else if(class(x$FEMbasis$mesh)=="mesh.3D"){
	R_plot_volume(x,...)
}
}


#' 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
#' @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, ...)
#' @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,...){

  # if(!require(rgl)){
  #   stop("The plot mesh.2.5D_function(...) requires the R package rgl, please install it and try again!")
  # }
  # 
  mesh<-x
  triangles = c(t(mesh$triangles))
  ntriangles=mesh$ntriangles
  
  order=mesh$order
  
  nodes=mesh$nodes
  
  edges=matrix(rep(0,6*ntriangles),ncol=2)
  for(i in 0:(ntriangles-1)){
  edges[3*i+1,]=c(triangles[3*order*i+1],triangles[3*order*i+2])
  edges[3*i+2,]=c(triangles[3*order*i+1],triangles[3*order*i+3])
  edges[3*i+3,]=c(triangles[3*order*i+2],triangles[3*order*i+3])
  }
  edges=edges[!duplicated(edges),]
  edges<-as.vector(t(edges))
  open3d()
  axes3d()
  rgl.pop("lights") 
  light3d(specular="black") 
  
  rgl.points(x = nodes[ ,1], y = nodes[ ,2], 
                  z=nodes[,3],col="black", ...)
                  
  rgl.lines(x = nodes[edges ,1], y = nodes[edges ,2], 
                  z=nodes[edges,3],col="black",...)
    
  aspect3d("iso")
  rgl.viewpoint(0,-45)

}


plot.mesh.3D<-function(x,...){

  # if(!require(rgl)){
  #   stop("The plot mesh.2.5D_function(...) requires the R package rgl, please install it and try again!")
  # }
  #if(!require(geometry)){
   # stop("The plot mesh.3D_function(...) requires the R package geometry, please install it and try again!")
  # }
  mesh<-x
  tetrahedrons = mesh$tetrahedrons
  ntetrahedrons=mesh$ntetrahedrons
  
  nodes=mesh$nodes
 
  tet=t(rbind(tetrahedrons[,-1],tetrahedrons[,-2],tetrahedrons[,-3],tetrahedrons[,-4]))
  
  open3d()
  axes3d()
  rgl.pop("lights") 
  light3d(specular="black") 
  
  rgl.triangles(nodes[tet,1],nodes[tet,2],nodes[tet,3],col="white",...)
    
  aspect3d("iso")
  rgl.viewpoint(0,-45)

}
 
 
 R_plot.ORD1.FEM = function(FEM, ...)  
 {
   # PLOT  Plots a FEM object FDOBJ over a rectangular grid defined by 
   # vectors X and Y;
   #
   
   
   
   #   if (!is.fd(fdobj))
   #   {
   #     stop('FDOBJ is not an FD object')
   #   }
   
   nodes = FEM$FEMbasis$mesh$nodes
   triangles = FEM$FEMbasis$mesh$triangles
   
   coeff = FEM$coeff
   
   FEMbasis = FEM$basis
   
   mesh = FEMbasis$mesh
   
   heat = heat.colors(100)
   
   nsurf = dim(coeff)[[2]]
   for (isurf in 1:nsurf)
   {
     open3d()
     axes3d()
     
     z = coeff[as.vector(t(triangles)),isurf]
     rgl.triangles(x = nodes[as.vector(t(triangles)) ,1], y = nodes[as.vector(t(triangles)) ,2], 
                   z=coeff[as.vector(t(triangles)),isurf], 
                   color = heat[round(99*(z- min(z))/(max(z)-min(z)))+1],...)
     aspect3d(2,2,1)
     rgl.viewpoint(0,-45)
     if (nsurf > 1)
     {readline("Press a button for the next plot...")}
   }
 }
 
 R_plot.ORDN.FEM = function(FEM, num_refinements, ...)  
 {
   coeff = FEM$coeff
   
   FEMbasis = FEM$FEMbasis
   
   mesh = FEMbasis$mesh
   
   heat = heat.colors(100)
   
   coeff = FEM$coeff
   
   # num_refinements sets the number od division on each triangle edge to be applied for rifenment
   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
   
   
   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(FEMbasis$transf_coord$diff1x[i],FEMbasis$transf_coord$diff2x[i]),c(FEMbasis$transf_coord$diff1y[i],FEMbasis$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, 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()
     rgl.pop("lights") 
     light3d(specular="black") 
     z = locations[as.vector(t(triangles)), 2 + isurf]
     rgl.triangles(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)
     rgl.viewpoint(0,-45)
     if (nsurf > 1)
     {readline("Press a button for the next plot...")}
   }
 }
 
 
 R_plot_manifold = function(FEM, ...){
   triangles = c(t(FEM$FEMbasis$mesh$triangles))
   ntriangles = FEM$FEMbasis$mesh$ntriangles
   order=FEM$FEMbasis$mesh$order
   nodes=FEM$FEMbasis$mesh$nodes
   edges=matrix(rep(0,6*ntriangles),ncol=2)
   for(i in 0:(ntriangles-1)){
     edges[3*i+1,]=c(triangles[3*order*i+1],triangles[3*order*i+2])
     edges[3*i+2,]=c(triangles[3*order*i+1],triangles[3*order*i+3])
     edges[3*i+3,]=c(triangles[3*order*i+2],triangles[3*order*i+3])
   }
   edges=edges[!duplicated(edges),]
   edges<-as.vector(t(edges))
   
   coeff = FEM$coeff
   
   FEMbasis = FEM$FEMbasis
   
   mesh = FEMbasis$mesh
   
   p=jet.col(n=128,alpha=0.8)
   #p <- colorRampPalette(c("#0E1E44","#3E6DD8","#68D061","#ECAF53", "#EB5F5F","#E11F1C"))(128)
   palette(p)
   
   ncolor=length(p)
   
   nsurf = dim(coeff)[[2]]
   for (isurf in 1:nsurf)
   {
     open3d()
     axes3d()
     rgl.pop("lights") 
     light3d(specular="black") 
     
     diffrange = max(coeff[,isurf])-min(coeff[,isurf])
     
     col = coeff[triangles,isurf]
     col= (col - min(coeff[,isurf]))/diffrange*(ncolor-1)+1
     
     rgl.triangles(x = nodes[triangles ,1], y = nodes[triangles ,2], 
                   z=nodes[triangles,3], 
                   color = col,...)
     rgl.lines(x = nodes[edges ,1], y = nodes[edges ,2], 
               z=nodes[edges,3], 
               color = "black",...)
     aspect3d("iso")
     rgl.viewpoint(0,-45)
     if (nsurf > 1 && isurf<nsurf)
     {readline("Press a button for the next plot...")}
   }
 }
 
 
 
 R_plot_volume = function(FEM,...){
   
   # if(!require(rgl)){
   #   stop("The plot mesh.3D_function(...) requires the R package rgl, please install it and try again!")
   # }
   # # if(!require(geometry)){
   #   stop("The plot mesh.3D_function(...) requires the R package geometry, please install it and try again!")
   # }
   
   tetrahedrons = FEM$FEMbasis$mesh$tetrahedrons
   
   nodes=FEM$FEMbasis$mesh$nodes
   
   ntetrahedrons = FEM$FEMbasis$mesh$ntetrahedrons
   
   tet=t(rbind(tetrahedrons[,-1],tetrahedrons[,-2],tetrahedrons[,-3],tetrahedrons[,-4]))
   
   coeff = FEM$coeff
   
   nsurf = dim(coeff)[[2]]
   
   p=jet.col(n=128,alpha=0.8)
   #p <- colorRampPalette(c("#0E1E44","#3E6DD8","#68D061","#ECAF53", "#EB5F5F","#E11F1C"))(128)
   ncolors=length(p)
   for (isurf in 1:nsurf)
   {	col=rep(0,ntetrahedrons)
   for(j in 1:ntetrahedrons)
     col[j]=mean(c(coeff[tetrahedrons[j,1],isurf],coeff[tetrahedrons[j,2],isurf],
                   coeff[tetrahedrons[j,3],isurf],coeff[tetrahedrons[j,4],isurf]))
   
   
   col=rep(1,3) %o% col
   
   
   diffrange = max(col)-min(col)
   
   col= (col - min(col))/diffrange*(ncolors-1)+1
   col=p[col]
   
   open3d()
   axes3d()
   rgl.pop("lights") 
   light3d(specular="black") 
   
   rgl.triangles(nodes[tet,1],nodes[tet,2],nodes[tet,3],col=col,alpha=0.7,...)
   
   #tetramesh(tetrahedrons,nodes,col=col,alpha=0.7)
   
   aspect3d("iso")
   rgl.viewpoint(0,-45)
   if (nsurf > 1 && isurf<nsurf)
   {readline("Press a button for the next plot...")}
   }
   
 }
 
  
 #' 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, ...)  
 #' @seealso \code{\link{FEM}} \code{\link{plot.FEM}}
 #' @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
 #' 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(class(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, ...)
   }
 }
 
 
 R_image.ORD1.FEM = function(FEM)  
 {
   # PLOT  Plots a FEM object FDOBJ over a rectangular grid defined by 
   # vectors X and Y;
   #
   
   #   if (!is.fd(fdobj))
   #   {
   #     stop('FDOBJ is not an FD object')
   #   }
   
   nodes = FEM$FEMbasis$mesh$nodes
   triangles = FEM$FEMbasis$mesh$triangles
   
   coeff = FEM$coeff
   
   FEMbasis = FEM$FEMbasis
   
   mesh = FEMbasis$mesh
   
   heat = heat.colors(100)
   
   nsurf = dim(coeff)[[2]]
   for (isurf in 1:nsurf)
   {
     open3d()
     axes3d()
     rgl.pop("lights") 
     light3d(specular="black") 
     z = coeff[as.vector(t(triangles)),isurf]
     rgl.triangles(x = nodes[as.vector(t(triangles)) ,1], y = nodes[as.vector(t(triangles)) ,2], 
                   z=0, 
                   color = heat[round(99*(z- min(z))/(max(z)-min(z)))+1])
     aspect3d(2,2,1)
     rgl.viewpoint(0,0)
     if (nsurf > 1)
     {readline("Press a button for the next plot...")}
   }
 }
 
 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
   
   for (i in 1:nrow(mesh$triangles))
   {
     transf<-rbind(cbind(FEMbasis$transf_coord$diff1x[i],FEMbasis$transf_coord$diff2x[i]),c(FEMbasis$transf_coord$diff1y[i],FEMbasis$transf_coord$diff2y[i]))
     pointsi = t(transf%*%t(meshi$nodes) + mesh$nodes[mesh$triangles[i,1],])
     z = R_eval_local.FEM(FEM, locations = pointsi, element_index = i)
     
     #mesh3 <- addNormals(subdivision3d(tmesh3d(vertices = t(locations), indices = t(triangles), homogeneous = FALSE)),deform = TRUE)
     
     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()
     rgl.pop("lights") 
     light3d(specular="black") 
     z = locations[as.vector(t(triangles)), 2 + isurf];
     rgl.triangles(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)
     rgl.viewpoint(0,0)
     if (nsurf > 1)
     {readline("Press a button for the next plot...")}
   }
 }
 
 
 R_eval_local.FEM = function(FEM, 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 = FEMbasis$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 = FEMbasis$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)
 }
 

Try the fdaPDE package in your browser

Any scripts or data that you put into this service are public.

fdaPDE documentation built on July 2, 2020, 2:22 a.m.