# ----------------------------------------------------------------------------
# Example 1: smoothing data over a circular region
# ----------------------------------------------------------------------------
# Set up a FEM object with a approximatedly circular boundary with 12 sides,
# and two rings of 12 points plus a point at the centre.
angle = seq(0,11*pi/6,len=12)
# define the 24 point locations
rad = 2
ctr = c(0,0)
pts = matrix(0,25,2);
pts[ 1:12,1] = rad*cos(angle) + ctr[1]
pts[ 1:12,2] = rad*sin(angle) + ctr[2]
pts[13:24,1] = rad*cos(angle)/2 + ctr[1]
pts[13:24,2] = rad*sin(angle)/2 + ctr[2]
# define the edge indices
edg <- matrix(0,12,2)
for (i in 1:11) {
edg[i,1] <- i
edg[i,2] <- i+1
}
edg[12,] <- c(12,1)
# use delaunay triangulation to define the triangles
tri = delaunayn(pts)
plotFEM.mesh(pts, tri, shift = 0.02)
title("A 25-point curcular mesh")
# Create an order 1 basis
hexFEMbasis <- create.FEM.basis(pts, edg, tri, 1)
# locations of data
obspts <- pts
npts <- dim(obspts)[1]
points(obspts[,1], obspts[,2], col=2, cex=2, lwd=2)
hexfd = fd(matrix(1-c(obspts[,1]^2+obspts[,2]^2),npts,1),hexFEMbasis)
# heights of a hemisphere at these locations
hexhgts <- round(eval.FEM.fd(obspts,hexfd),2)
Run the code above in your browser using DataLab