Skip to content

Commit

Permalink
conformal torus mesh
Browse files Browse the repository at this point in the history
  • Loading branch information
@stla
stla committed Oct 3, 2023
1 parent cf80f91 commit ca2a430
Show file tree
Hide file tree
Showing 2 changed files with 334 additions and 74 deletions.
238 changes: 164 additions & 74 deletions R/someMeshes.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,132 @@ sphereMesh <- function(x = 0, y = 0, z = 0, r = 1, iterations = 3L) {
sphere
}

torusMesh1 <- function(R, r, nu, nv, nrmls) {
nu <- as.integer(nu)
nv <- as.integer(nv)
nunv <- nu*nv
vs <- matrix(NA_real_, nrow = 3L, ncol = nunv)
normals <- if(nrmls) matrix(NA_real_, nrow = nunv, ncol = 3L)
tris1 <- matrix(NA_integer_, nrow = 3L, ncol = nunv)
tris2 <- matrix(NA_integer_, nrow = 3L, ncol = nunv)
u_ <- seq(0, 2*pi, length.out = nu + 1L)[-1L]
cosu_ <- cos(u_)
sinu_ <- sin(u_)
v_ <- seq(0, 2*pi, length.out = nv + 1L)[-1L]
cosv_ <- cos(v_)
sinv_ <- sin(v_)
Rrcosv_ <- R + r*cosv_
rsinv_ <- r*sinv_
jp1_ <- c(2L:nv, 1L)
j_ <- 1L:nv
for(i in 1L:(nu-1L)){
i_nv <- i*nv
k1 <- i_nv - nv
rg <- (k1 + 1L):i_nv
cosu_i <- cosu_[i]
sinu_i <- sinu_[i]
vs[, rg] <- rbind(
cosu_i * Rrcosv_,
sinu_i * Rrcosv_,
rsinv_
)
if(nrmls) {
normals[rg, ] <- cbind(
cosu_i * cosv_,
sinu_i * cosv_,
sinv_
)
}
k_ <- k1 + j_
l_ <- k1 + jp1_
m_ <- i_nv + j_
tris1[, k_] <- rbind(m_, l_, k_)
tris2[, k_] <- rbind(m_, i_nv + jp1_, l_)
}
i_nv <- nunv
k1 <- i_nv - nv
rg <- (k1 + 1L):i_nv
vs[, rg] <- rbind(
Rrcosv_,
0,
rsinv_
)
if(nrmls) {
normals[rg, ] <- cbind(
cosv_,
0,
sinv_
)
}
l_ <- k1 + jp1_
k_ <- k1 + j_
tris1[, k_] <- rbind(j_, l_, k_)
tris2[, k_] <- rbind(j_, jp1_, l_)
tmesh3d(
vertices = vs,
indices = cbind(tris1, tris2),
normals = normals,
homogeneous = FALSE
)
}

torusMesh2 <- function(R, r, nu, nv) {
nu <- as.integer(nu)
nv <- as.integer(nv)
nunv <- nu*nv
vs <- matrix(NA_real_, nrow = 3L, ncol = nunv)
tris1 <- matrix(NA_integer_, nrow = 3L, ncol = nunv)
tris2 <- matrix(NA_integer_, nrow = 3L, ncol = nunv)
u_ <- seq(0, 2*pi, length.out = nu + 1L)[-1L]
cosu_ <- cos(u_)
sinu_ <- sin(u_)
v_ <- seq(0, 2*pi, length.out = nv + 1L)[-1L]
cosv_ <- cos(v_)
sinv_ <- sin(v_)
kxy <- R*R - r*r
kz <- sqrt(kxy) * r
kcosu_ <- kxy * cosu_
ksinu_ <- kxy * sinu_
w_ <- R - r * cosv_
h_ <- kz * sinv_ / w_
jp1_ <- c(2L:nv, 1L)
j_ <- 1L:nv
for(i in 1L:(nu-1L)){
i_nv <- i*nv
k1 <- i_nv - nv
rg <- (k1 + 1L):i_nv
kcosu_i <- kcosu_[i]
ksinu_i <- ksinu_[i]
vs[, rg] <- rbind(
kcosu_i / w_,
ksinu_i / w_,
h_
)
k_ <- k1 + j_
l_ <- k1 + jp1_
m_ <- i_nv + j_
tris1[, k_] <- rbind(m_, l_, k_)
tris2[, k_] <- rbind(m_, i_nv + jp1_, l_)
}
i_nv <- nunv
k1 <- i_nv - nv
rg <- (k1 + 1L):i_nv
vs[, rg] <- rbind(
kxy / w_,
0,
h_
)
l_ <- k1 + jp1_
k_ <- k1 + j_
tris1[, k_] <- rbind(j_, l_, k_)
tris2[, k_] <- rbind(j_, jp1_, l_)
tmesh3d(
vertices = vs,
indices = cbind(tris1, tris2),
homogeneous = FALSE
)
}

#' @title Torus mesh
#' @description Triangle mesh of a torus.
#'
Expand All @@ -36,6 +162,9 @@ sphereMesh <- function(x = 0, y = 0, z = 0, r = 1, iterations = 3L) {
#' \code{NULL}, the torus has equatorial plane z=0 and the
#' z-axis as revolution axis
#' @param nu,nv numbers of subdivisions, integers (at least 3)
#' @param normals a Boolean value, whether to compute the normals of the mesh
#' @param conformal a Boolean value, whether to use a conformal
#' parameterization of the torus
#'
#' @return A triangle \strong{rgl} mesh (class \code{mesh3d}).
#' @export
Expand All @@ -48,8 +177,8 @@ sphereMesh <- function(x = 0, y = 0, z = 0, r = 1, iterations = 3L) {
#' mesh <- torusMesh(R = 3, r = 1)
#' open3d(windowRect = 50 + c(0, 0, 512, 512))
#' view3d(0, 0, zoom = 0.75)
#' shade3d(mesh, color = "green")
#' wire3d(mesh)
#' shade3d(mesh, color = "green", polygon_offset = 1)
#' wire3d(mesh, color = "black")
#'
#' # Villarceau circles ####
#' Villarceau <- function(beta, theta0, phi) {
Expand Down Expand Up @@ -77,7 +206,10 @@ sphereMesh <- function(x = 0, y = 0, z = 0, r = 1, iterations = 3L) {
#' rmesh <- torusMesh(r = 0.05, p1 = p1, p2 = p2, p3 = p3)
#' shade3d(rmesh, color = colors[i])
#' }}
torusMesh <- function(R, r, p1 = NULL, p2 = NULL, p3 = NULL, nu = 50, nv = 30) {
torusMesh <- function(
R, r, p1 = NULL, p2 = NULL, p3 = NULL, nu = 50, nv = 30,
normals = TRUE, conformal = FALSE
) {
transformation <- !is.null(p1) && !is.null(p2) && !is.null(p3)
if(transformation) {
ccircle <- circumcircle(p1, p2, p3)
Expand All @@ -89,75 +221,31 @@ torusMesh <- function(R, r, p1 = NULL, p2 = NULL, p3 = NULL, nu = 50, nv = 30) {
stopifnot(isPositiveNumber(R), isPositiveNumber(r))
stopifnot(R > r)
stopifnot(nu >= 3, nv >= 3)
nu <- as.integer(nu)
nv <- as.integer(nv)
nunv <- nu*nv
vs <- matrix(NA_real_, nrow = 3L, ncol = nunv)
normals <- matrix(NA_real_, nrow = nunv, ncol = 3L)
tris1 <- matrix(NA_integer_, nrow = 3L, ncol = nunv)
tris2 <- matrix(NA_integer_, nrow = 3L, ncol = nunv)
u_ <- seq(0, 2*pi, length.out = nu + 1L)[-1L]
cosu_ <- cos(u_)
sinu_ <- sin(u_)
v_ <- seq(0, 2*pi, length.out = nv + 1L)[-1L]
cosv_ <- cos(v_)
sinv_ <- sin(v_)
Rrcosv_ <- R + r*cosv_
rsinv_ <- r*sinv_
jp1_ <- c(2L:nv, 1L)
j_ <- 1L:nv
for(i in 1L:(nu-1L)){
i_nv <- i*nv
k1 <- i_nv - nv
rg <- (k1 + 1L):i_nv
cosu_i <- cosu_[i]
sinu_i <- sinu_[i]
vs[, rg] <- rbind(
cosu_i * Rrcosv_,
sinu_i * Rrcosv_,
rsinv_
)
normals[rg, ] <- cbind(
cosu_i * cosv_,
sinu_i * cosv_,
sinv_
)
k_ <- k1 + j_
l_ <- k1 + jp1_
m_ <- i_nv + j_
tris1[, k_] <- rbind(m_, l_, k_)
tris2[, k_] <- rbind(m_, i_nv + jp1_, l_)
}
i_nv <- nunv
k1 <- i_nv - nv
rg <- (k1 + 1L):i_nv
vs[, rg] <- rbind(
Rrcosv_,
0,
rsinv_
)
normals[rg, ] <- cbind(
cosv_,
0,
sinv_
)
l_ <- k1 + jp1_
k_ <- k1 + j_
tris1[, k_] <- rbind(j_, l_, k_)
tris2[, k_] <- rbind(j_, jp1_, l_)
rmesh <- tmesh3d(
vertices = vs,
indices = cbind(tris1, tris2),
normals = normals,
homogeneous = FALSE
)
if(transformation) {
rmesh <- translate3d(
rotate3d(
rmesh, matrix = rotMatrix
),
x = center[1L], y = center[2L], z = center[3L]
)
stopifnot(isBoolean(normals))
stopifnot(isBoolean(conformal))
if(conformal) {
rmesh <- torusMesh2(R, r, nu, nv)
if(transformation) {
rmesh <- translate3d(
rotate3d(
rmesh, matrix = rotMatrix
),
x = center[1L], y = center[2L], z = center[3L]
)
}
if(normals) {
rmesh <- addNormals(rmesh)
}
} else {
rmesh <- torusMesh1(R, r, nu, nv, normals)
if(transformation) {
rmesh <- translate3d(
rotate3d(
rmesh, matrix = rotMatrix
),
x = center[1L], y = center[2L], z = center[3L]
)
}
}
rmesh
}
Expand All @@ -168,6 +256,8 @@ torusMesh <- function(R, r, p1 = NULL, p2 = NULL, p3 = NULL, nu = 50, nv = 30) {
#' @param a,c,mu cyclide parameters, positive numbers such that
#' \code{c < mu < a}
#' @param nu,nv numbers of subdivisions, integers (at least 3)
#' @param conformal a Boolean value, whether to use a conformal
#' parameterization of the torus
#'
#' @return A triangle \strong{rgl} mesh (class \code{mesh3d}).
#'
Expand All @@ -194,7 +284,7 @@ torusMesh <- function(R, r, p1 = NULL, p2 = NULL, p3 = NULL, nu = 50, nv = 30) {
#' wire3d(mesh)
#' shade3d(sphere, color = "red")
#' wire3d(sphere)
cyclideMesh <- function(a, c, mu, nu = 90L, nv = 40L){
cyclideMesh <- function(a, c, mu, nu = 90L, nv = 40L, conformal = FALSE){
stopifnot(c > 0, a > mu, mu > c)
stopifnot(nu >= 3, nv >= 3)
nu <- as.integer(nu)
Expand All @@ -220,7 +310,7 @@ cyclideMesh <- function(a, c, mu, nu = 90L, nv = 40L){
b2 <- (a*mu*(c+mu)*(a-c) + bb2 - c*c + bb*(c*(a-mu-c) + 2*a*mu))/denb2
omegaT <- (b1 + b2)/2
OmegaT <- c(omegaT, 0, 0)
tormesh <- torusMesh(R, r, nu = nu, nv = nv)
tormesh <- torusMesh(R, r, nu = nu, nv = nv, conformal = conformal)
rtnormals <- r * tormesh[["normals"]][1L:3L, ]
xvertices <- tormesh[["vb"]][1L:3L, ] + OmegaT
for(i in 1L:nu){
Expand Down
Loading

0 comments on commit ca2a430

Please sign in to comment.