A simulation of a soap film modeled as a minimal surface satisfying Laplace's equation, \(\nabla^2 u = 0\), on a unit square domain \([0,1] \times [0,1]\). The simulation uses the fdaPDE package in R to implement the Finite Element Method (FEM) and plotly for 3D visualization. The R script solves the PDE with Dirichlet boundary conditions \(u(x,y) = (1-x)(1-y)\) and visualizes the resulting surface and boundary frame. Below, we break down the R script, explaining mesh creation, boundary conditions, FEM solution, and visualization. .
The simulation begins by loading the necessary R libraries: fdaPDE for FEM computations, plotly for interactive 3D visualizations, and RColorBrewer for color schemes. The main function, solve_soap_film_fdaPDE(), orchestrates mesh generation, boundary condition setup, PDE solving, and result preparation.
library(fdaPDE)
library(plotly)
library(RColorBrewer)
# Function to solve soap film using fdaPDE FEM library
solve_soap_film_fdaPDE <- function(n = 50) {
A triangular mesh is created over the unit square \([0,1] \times [0,1]\) to discretize the domain for FEM. The mesh divides the domain into small triangular elements, enabling numerical solution of the PDE.
# Create rectangular domain mesh
x_coords <- seq(0, 1, length.out = n)
y_coords <- seq(0, 1, length.out = n)
# Create grid points
locations <- expand.grid(x = x_coords, y = y_coords)
# Define boundary vertices
boundary_vertices <- rbind(
c(0, 0), # bottom-left
c(1, 0), # bottom-right
c(1, 1), # top-right
c(0, 1) # top-left
)
# Create refined mesh
mesh <- create.mesh.2D(nodes = boundary_vertices,
segments = rbind(c(1,2), c(2,3), c(3,4), c(4,1)),
holes = NULL, triangles = NULL, order = 1, verbosity = 0)
# Refine mesh for better resolution
mesh <- refine.mesh.2D(mesh, maximum_area = 0.001, minimum_angle = 20)
The code generates a grid of points using seq and expand.grid. Boundary vertices are defined at the corners of the square, and create.mesh.2D constructs a triangular mesh. The mesh is refined with refine.mesh.2D to ensure triangles have a maximum area of 0.001 and a minimum angle of 20 degrees, enhancing numerical accuracy and stability.
The soap film is constrained by the Dirichlet boundary condition \(u(x,y) = (1-x)(1-y)\) on the edges of the unit square. This defines the height of the film at the boundary, forming the "frame" of the soap film.
# Get mesh coordinates
mesh_locations <- mesh$nodes
n_nodes <- nrow(mesh_locations)
# Define boundary conditions: u = (1-x)*(1-y) on boundary
boundary_indices <- c()
boundary_values <- c()
for (i in 1:n_nodes) {
x_i <- mesh_locations[i, 1]
y_i <- mesh_locations[i, 2]
# Check if point is on boundary (within tolerance)
tol <- 1e-10
if (abs(x_i) < tol || abs(x_i - 1) < tol || abs(y_i) < tol || abs(y_i - 1) < tol) {
boundary_indices <- c(boundary_indices, i)
boundary_values <- c(boundary_values, (1 - x_i) * (1 - y_i))
}
}
# Set up boundary conditions
BC <- NULL
if (length(boundary_indices) > 0) {
BC <- list(BC_indices = boundary_indices, BC_values = boundary_values)
}
Boundary nodes are identified by checking if their coordinates lie on the domain edges (within a tolerance of \(10^{-10}\)). The condition \(u = (1-x)(1-y)\) is applied to these nodes, and their indices and values are stored for the FEM solver.
The minimal surface of the soap film satisfies Laplace’s equation, \(\nabla^2 u = 0\). The FEM discretizes this PDE over the triangular mesh, solving for \(u(x,y)\) that meets the boundary conditions.
# Create FEM object for Laplace equation
observations <- rep(0, n_nodes) # Zero forcing term
# Create FEMbasis object
FEMbasis <- create.FEM.basis(mesh)
# Create synthetic data at mesh nodes for boundary fitting
data_values <- rep(0, n_nodes)
for (i in 1:n_nodes) {
x_i <- mesh_locations[i, 1]
y_i <- mesh_locations[i, 2]
data_values[i] <- (1 - x_i) * (1 - y_i) # Target function
}
# Solve FEM system
solution <- smooth.FEM.basis(observations = data_values,
FEMbasis = FEMbasis,
lambda = 0,
BC = BC)
# Evaluate solution on regular grid for visualization
x_eval <- seq(0, 1, length.out = 50)
y_eval <- seq(0, 1, length.out = 50)
eval_locations <- expand.grid(x = x_eval, y = y_eval)
# Evaluate FEM solution at regular grid points
eval_result <- eval.FEM(solution$fit.FEM, locations = eval_locations)
# Reshape to matrix for plotting
z_matrix <- matrix(eval_result, nrow = length(x_eval), ncol = length(y_eval), byrow = TRUE)
The create.FEM.basis function sets up a linear FEM basis. The smooth.FEM.basis function solves the PDE with a zero forcing term (homogeneous Laplace equation) and the specified boundary conditions. The parameter lambda = 0 ensures pure interpolation at the boundary. The solution is evaluated on a 50x50 grid for visualization, and results are reshaped into a matrix.
The boundary frame is visualized with alternating pink and black dashed-dotted lines to enhance aesthetics. The create_alternating_segments function generates these segments.
create_alternating_segments <- function(x, y, z, color1 = "#db00fd", color2 = "#b738ca",
dash_len = 3, dot_len = 1) {
segments <- list()
n <- length(x)
i <- 1
while (i < n) {
# Purple dash: longer segment
end_dash <- min(i + dash_len, n)
seg_x <- x[i:end_dash]
seg_y <- y[i:end_dash]
seg_z <- z[i:end_dash]
segments <- append(segments, list(list(
x = seg_x, y = seg_y, z = seg_z,
type = "scatter3d", mode = "lines",
line = list(color = color1, width = 15),
showlegend = FALSE
)), after = length(segments))
i <- end_dash
if (i >= n) break
# Black dot: very short segment
end_dot <- min(i + dot_len, n)
seg_x <- x[i:end_dot]
seg_y <- y[i:end_dot]
seg_z <- z[i:end_dot]
if (length(seg_x) < 2) {
seg_x <- c(seg_x[1], seg_x[1] + 1e-6)
seg_y <- c(seg_y[1], seg_y[1])
seg_z <- c(seg_z[1], seg_z[1])
}
segments <- append(segments, list(list(
x = seg_x, y = seg_y, z = seg_z,
type = "scatter3d", mode = "lines",
line = list(color = color2, width = 10),
showlegend = FALSE
)), after = length(segments))
i <- end_dot
}
return(segments)
}
This function alternates between pink dashes (color #db00fd, width 15) and black dots (color #b738ca, width 10), creating a visually distinct boundary outline.
Two interactive 3D visualizations are produced using plotly, saved as separate HTML files for modularity. These can be viewed by clicking the links below.
The rectangular frame is plotted with a slight lift in the z-direction (\(u(x,y) = (1-x)(1-y) + 0.2\)) using solid black lines to represent the boundary above the soap film surface.
fig1 <- plot_ly() %>%
add_trace(
x = boundary_x, y = boundary_y, z = boundary_z_top,
type = "scatter3d", mode = "lines",
line = list(color = "black", width = 15),
name = "Top Outline"
) %>%
layout(
title = "Top Lifted Outline",
scene = list(
xaxis = list(title = "x"),
yaxis = list(title = "y"),
zaxis = list(title = "z")
),
margin = list(l = 0, r = 0, t = 40, b = 0),
height = 600
)
htmlwidgets::saveWidget(fig1, "structure.html", selfcontained = TRUE)
The soap film is visualized as a 3D surface, with the boundary frame overlaid using alternating pink and black dashed-dotted lines. The surface uses the "Viridis" colorscale to represent height values.
fdaPDE package. The surface smoothly satisfies Laplace's equation with imposed boundary values.Figure 2: Soap film surface computed with FEM using the fdaPDE package, with boundary frame in pink and black dashed-dotted lines.
fig2 <- plot_ly() %>%
add_surface(
x = x_grid, y = y_grid, z = result$z,
colorscale = "Viridis",
opacity = 0.95,
showscale = TRUE,
name = "Soap Film (fdaPDE FEM Solution)"
)
segments <- create_alternating_segments(boundary_x, boundary_y, boundary_z_soap,
"#db00fd", "#b738ca", dash_len = 1, dot_len = 1)
for (seg in segments) {
fig2 <- fig2 %>% add_trace(
x = seg$x, y = seg$y, z = seg$z,
type = seg$type, mode = seg$mode,
line = seg$line,
showlegend = seg$showlegend
)
}
fig2 <- fig2 %>%
layout(
title = "Soap Film Using fdaPDE FEM Modeling with Triangular Mesh",
scene = list(
xaxis = list(title = "x"),
yaxis = list(title = "y"),
zaxis = list(title = "z")
),
margin = list(l = 0, r = 0, t = 40, b = 0),
height = 600
)
htmlwidgets::saveWidget(fig2, "soap_film.html", selfcontained = TRUE)
The top frame uses a scatter3d trace for the boundary outline, while the soap film uses a surface trace. The boundary segments are added to the surface plot using the create_alternating_segments function for visual enhancement.
In case the FEM solver fails, a fallback approach computes the analytical solution \(u(x,y) = (1-x)(1-y)\), which is visualized similarly to the FEM solution.
# Fallback to basic implementation
x <- seq(0, 1, length.out = 50)
y <- seq(0, 1, length.out = 50)
x_grid <- outer(x, rep(1, length(y)))
y_grid <- outer(rep(1, length(x)), y)
z_basic <- (1 - x_grid) * (1 - y_grid)
This fallback directly evaluates the boundary condition across the domain, bypassing FEM computations for robustness. The visualizations are generated similarly and saved to the same HTML files.
This simulation demonstrates:
fdaPDE.