Nothing
##' A simple mesh generator for non-convex regions in n-D space
##'
##' An unstructured simplex requires a choice of meshpoints (vertex nodes) and
##' a triangulation. This is a simple and short algorithm that improves the
##' quality of a mesh by relocating the meshpoints according to a relaxation
##' scheme of forces in a truss structure. The topology of the truss is reset
##' using Delaunay triangulation. A (sufficiently smooth) user supplied signed
##' distance function (\code{fd}) indicates if a given node is inside or
##' outside the region. Points outside the region are projected back to the
##' boundary.
##'
##' This is an implementation of original Matlab software of Per-Olof Persson.
##'
##' Excerpt (modified) from the reference below:
##'
##' \sQuote{The algorithm is based on a mechanical analogy between a triangular
##' mesh and a n-D truss structure. In the physical model, the edges of the
##' Delaunay triangles of a set of points correspond to bars of a truss. Each
##' bar has a force-displacement relationship \eqn{f(\ell, \ell_{0})}{F(L,L0)}
##' depending on its current length \eqn{\ell}{L} and its unextended length
##' \eqn{\ell_{0}}{L0}.}
##'
##' \sQuote{External forces on the structure come at the boundaries, on which
##' external forces have normal orientations. These external forces are just
##' large enough to prevent nodes from moving outside the boundary. The
##' position of the nodes are the unknowns, and are found by solving for a
##' static force equilibrium. The hope is that (when \code{fh = function(p)
##' return(rep(1,nrow(p)))}), the lengths of all the bars at equilibrium will
##' be nearly equal, giving a well-shaped triangular mesh.}
##'
##' See the references below for all details. Also, see the comments in the
##' source file of \code{distmesh2d}.
##'
##' @param fdist Vectorized signed distance function, for example
##' \code{\link{mesh.dsphere}}, accepting an \code{m}-by-\code{n}
##' matrix, where \code{m} is arbitrary, as the first argument.
##' @param fh Vectorized function, for example \code{\link{mesh.hunif}},
##' that returns desired edge length as a function of position.
##' Accepts an \code{m}-by-\code{n} matrix, where \code{n} is
##' arbitrary, as its first argument.
##' @param h Initial distance between mesh nodes.
##' @param box \code{2}-by-\code{n} matrix that specifies the bounding box.
##' (See \link{distmesh2d} for an example.)
##' @param pfix \code{nfix}-by-2 matrix with fixed node positions.
##' @param \dots parameters that are passed to \code{fdist} and \code{fh}
##' @param ptol Algorithm stops when all node movements are smaller than
##' \code{dptol}
##' @param ttol Controls how far the points can move (relatively) before a
##' retriangulation with \code{\link{delaunayn}}.
##' @param deltat Size of the time step in Eulers method.
##' @param geps Tolerance in the geometry evaluations.
##' @param deps Stepsize \eqn{\Delta x} in numerical derivative computation for
##' distance function.
##' @return \code{m}-by-\code{n} matrix with node positions.
##' @section Wishlist : \itemize{ \item*Implement in C/Fortran \item*Translate
##' other functions of the matlab package }
##' @author Raoul Grasman; translated from original Matlab sources of Per-Olof
##' Persson.
##' @seealso \code{\link{distmesh2d}}, \code{\link[tripack]{tri.mesh}},
##' \code{\link{delaunayn}}, \code{\link{mesh.dsphere}},
##' \code{\link{mesh.hunif}},\cr \code{\link{mesh.diff}},
##' \code{\link{mesh.union}}, \code{\link{mesh.intersect}}
##' @references \url{http://persson.berkeley.edu/distmesh/}
##'
##' \cite{P.-O. Persson, G. Strang, A Simple Mesh Generator in MATLAB. SIAM
##' Review, Volume 46 (2), pp. 329-345, June 2004}
##' @keywords math optimize dplot graphs
##' @examples
##'
##' \dontrun{
##' # examples distmeshnd
##' require(rgl)
##'
##' fd = function(p, ...) sqrt((p^2)%*%c(1,1,1)) - 1
##' # also predefined as `mesh.dsphere'
##' fh = function(p,...) rep(1,nrow(p))
##' # also predefined as `mesh.hunif'
##' bbox = matrix(c(-1,1),2,3)
##' p = distmeshnd(fd,fh,0.2,bbox, maxiter=100)
##' # this may take a while:
##' # press Esc to get result of current iteration
##' }
##'
##' @export
"distmeshnd" <-
function (fdist, fh, h, box, pfix = array(dim = c(0, ncol(box))),
..., ptol = 0.001, ttol = 0.1, deltat = 0.1, geps = 0.1 *
h, deps = sqrt(.Machine$double.eps) * h)
{
# %DISTMESHND N-D Mesh Generator using Distance Functions.
dim = ncol(as.matrix(box))
L0mult = 1 + 0.4/2^(dim - 1)
rownorm2 = function(x) drop(sqrt((x^2) %*% rep(1, ncol(x))))
# %1. Create initial distribution in bounding box
if (dim == 1) {
p = seq(box[1], box[2], by = h)
}
else {
cbox = lapply(1:dim, function(ii) seq(box[1, ii], box[2,
ii], by = h))
p = do.call("expand.grid", cbox)
p = as.matrix(p)
}
# %2. Remove points outside the region, apply the rejection method
p = p[fdist(p, ...) < geps, ]
r0 = fh(p, ...)
p = rbind(pfix, p[stats::runif(nrow(p)) < min(r0)^dim/r0^dim, ])
N = nrow(p)
if (N <= dim + 1)
stop("Not enough starting points inside boundary (is h0 too large?).")
on.exit(return(invisible(p)))
cat("Press esc if the mesh seems fine but the algorithm hasn't converged.\n")
utils::flush.console()
count = 0
p0 = 1/.Machine$double.eps
# mimick Matlab call ``localpairs=nchoosek(1:dim+1,2)'':
localpairs = as.matrix(expand.grid(1:(dim + 1), 1:(dim + 1)))
localpairs = localpairs[lower.tri(matrix(TRUE, dim + 1, dim + 1)), 2:1]
while (TRUE) {
if (max(rownorm2(p - p0)) > ttol * h) {
# %3. Retriangulation by Delaunay:
p0 = p
t = delaunayn(p)
pmid = matrix(0, nrow(t), dim)
for (ii in 1:(dim + 1)) pmid = pmid + p[t[, ii],
]/(dim + 1)
t = t[fdist(pmid, ...) < (-geps), ]
pair = array(dim = c(0, 2))
for (ii in 1:nrow(localpairs)) {
pair = rbind(pair, t[, localpairs[ii, ]])
}
# %4. Describe each edge by a unique pair of nodes
pair = Unique(pair, TRUE); # base-function `unique' is way too slow
if (dim == 2) {
trimesh(t, p[, 1:2])
}
else if (dim == 3) {
if (count%%5 == 0) {
tetramesh(t, p)
}
}
else {
cat("Retriangulation #", 15, "\n")
utils::flush.console()
}
count = count + 1
}
bars = p[pair[, 1], ] - p[pair[, 2], ]
L = rownorm2(bars)
L0 = fh((p[pair[, 1], ] + p[pair[, 2], ])/2, ...)
L0 = L0 * L0mult * (sum(L^dim)/sum(L0^dim))^(1/dim)
F = L0 - L
F[F < 0] = 0
Fbar = cbind(bars, -bars) * matrix(F/L, nrow = nrow(bars),
ncol = 2 * dim)
ii = pair[, t(matrix(1:2, 2, dim))]
jj = rep(1, nrow(pair)) %o% c(1:dim, 1:dim)
s = c(Fbar)
ns = length(s)
dp = matrix(0, N, dim)
dp[1:(dim * N)] = rowsum(s, ii[1:ns] + ns * (jj[1:ns] -
1))
if (nrow(pfix) > 0)
dp[1:nrow(pfix), ] = 0
p = p + deltat * dp
d = fdist(p, ...)
ix = d > 0
gradd = matrix(0, sum(ix), dim)
for (ii in 1:dim) {
a = rep(0, dim)
a[ii] = deps
d1x = fdist(p[ix, ] + rep(1, sum(ix)) %o% a, ...)
gradd[, ii] = (d1x - d[ix])/deps
}
p[ix, ] = p[ix, ] - (d[ix] %o% rep(1, dim)) * gradd
maxdp = max(deltat * rownorm2(dp[d < (-geps), ]))
if (maxdp < ptol * h)
break
}
}
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.