rm(list=ls())
# Load required packages
library(readr)
library(ggplot2)
library(tidyr)
library(plotly)
library(processx)
library(reshape2)
#font
t <- list(
family = "Arial",
size = 14,
color = 'black')
# Read the CSV file
gyroid_u1.2<-read.csv(file='gyroid_geometry_cong_u1.2.csv', header = TRUE, dec='.', sep=',')
# rela_density ~ C: rela_density = a*C + b
lm_model <- lm(rela_density ~ C, data = gyroid_u1.2)
# Extract the coefficients
a <- coef(lm_model)[2] # coefficient for C
b <- coef(lm_model)[1] # intercept term
# Print the coefficients
cat("Coefficient a:", a, "\n")
## Coefficient a: -0.3331213
cat("Coefficient b:", b, "\n")
## Coefficient b: 0.5000124
# Get the model summary
summary(lm_model)
##
## Call:
## lm(formula = rela_density ~ C, data = gyroid_u1.2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0119839 -0.0039080 0.0001466 0.0038832 0.0119592
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.5000124 0.0009114 548.6 <2e-16 ***
## C -0.3331213 0.0010499 -317.3 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.004984 on 28 degrees of freedom
## Multiple R-squared: 0.9997, Adjusted R-squared: 0.9997
## F-statistic: 1.007e+05 on 1 and 28 DF, p-value: < 2.2e-16
# Calculate the correlation coefficient
correlation <- cor(gyroid_u1.2$rela_density, gyroid_u1.2$C)
cat("Correlation coefficient:", correlation, "\n")
## Correlation coefficient: -0.999861
# Create scatter plot
scatter_plot_1 <- ggplot(gyroid_u1.2, aes(x = C, y = rela_density)) +
geom_point() +
labs(x = "C", y = "Relative density [%]") +
scale_x_continuous(breaks = seq(-1.5, 1.5, 0.6)) +
scale_y_continuous(breaks = seq(0, 1, 0.2), labels = paste0(seq(0, 100, 20), "%")) +
theme(
axis.text = element_text(family = "Arial",size = 10,color = "black"),
axis.title = element_text(family = "Arial",size = 10),
axis.ticks = element_line(color = "black", size = 0.5),
axis.ticks.length = unit(-0.1, "mm"),
panel.background = element_rect(fill = "white", color = "black"),
panel.border = element_rect(color = "black",size = 1, fill = NA),
#panel.grid.major.y = element_line(color = "gray80"), # Set the major grid lines color
#panel.grid.major.x = element_line(color = "gray80"),
)+
geom_abline(intercept = b, slope = a, linetype = "dashed", color = "black")
scatter_plot_1
pore_channel, pore_core, strut diameter shows “Hyperbolic paraboloid” with unit cell size and C (relative density). analytical result as red surfaces
# Read the CSV file
gyroid<-read.csv(file='gyroid_geometry_cong_c-1to1.csv', header = TRUE, dec='.', sep=',')
gyroid$unit<-as.numeric(gyroid$unit)
gyroid_dia<-read.csv(file='gyroid_geometry_cong_c-1to1_dia.csv', header = TRUE, dec='.', sep=',')
gyroid_dia$unit<-as.numeric(gyroid_dia$unit)
pore_fig<-plot_ly(gyroid,x = ~unit, y = ~C, z = ~pore_core,
marker = list(color = ~pore_core, colorscale = 'Greens', cmin = -0.2, cmax = 3,showscale = TRUE)) %>%add_markers()
pore_fig <- pore_fig %>% layout(scene = list(xaxis = list(title = 'Unit cell size [mm]',
tickvals =c(1, 2, 3, 4),
ticktext = paste0(c(1, 2, 3, 4))),
yaxis = list(title = 'C',
tickvals =c(-1.0,-0.8,-0.5,-0.3,0.0,0.3,0.5,0.8,1.0),
ticktext = paste0(c(-1.0,"",-0.5,"",0.0,"",0.5,"",1.0))),
zaxis = list(title = 'Pore core size [mm]',
range = c(-0.2, 3.0)),
#camera perspective
camera = list(eye = list(x = -1.25, y = -1.25, z = 0.1))),
font=t,
annotations = list(
x = 1.13,
y = 1.05,
z = 1,
font=t,
text = 'Pore core size',
xref = 'paper',
yref = 'paper',
showarrow = FALSE
))
# Create the plane based on the equation
dummy_unit <- seq(1, 4, length.out = 9)
dummy_c <- seq(-0.9, 1, length.out = 20)
#an analytical surface for pore
ana_pore_core_data<-expand.grid(unit = dummy_unit,C = dummy_c,KEEP.OUT.ATTRS = F)
ana_pore_core_data$pore_core <-ana_pore_core_data$unit+ana_pore_core_data$unit*(atan((1-sqrt(2-ana_pore_core_data$C*ana_pore_core_data$C))/(ana_pore_core_data$C+1)))/pi-ana_pore_core_data$unit*(atan((1+sqrt(2-ana_pore_core_data$C*ana_pore_core_data$C))/(ana_pore_core_data$C+1)))/pi
ana_pore_core_data <- acast(ana_pore_core_data, C~unit , value.var = "pore_core") #y ~ x
# add analytical surface into 3d plot
pore_fig <- add_trace(p=pore_fig, type = "surface",
x = dummy_unit,
y = dummy_c,
z = ana_pore_core_data,
type = "surface",
showscale = FALSE,
colorscale = list(c(0, 1), c("#f0202080", "#f0202080"))
)
#
pore_fig
pore_fig_2<-plot_ly(gyroid, x = ~unit, y = ~C, z = ~pore_channel,
marker = list(color = ~pore_channel, colorscale = 'Greens', cmin = -0.2, cmax = 3,showscale = TRUE)) %>% add_markers()
pore_fig_2 <- pore_fig_2 %>% layout(scene = list(xaxis = list(title = 'Unit cell size [mm]',
tickvals =c(1, 2, 3, 4),
ticktext = paste0(c(1, 2, 3, 4))),
yaxis = list(title = 'C',
tickvals =c(-1.0,-0.8,-0.5,-0.3,0.0,0.3,0.5,0.8,1.0),
ticktext = paste0(c(-1.0,"",-0.5,"",0.0,"",0.5,"",1.0))),
zaxis = list(title = 'Pore channel size [mm]',
range = c(-0.2, 3.0)),
#camera perspective
camera = list(eye = list(x = -1.25, y = -1.25, z = 0.1))),
font=t,
annotations = list(
x = 1.13,
y = 1.05,
z = 1,
font=t,
text = 'Pore channel size',
xref = 'paper',
yref = 'paper',
showarrow = FALSE
))
# Create the plane based on the equation
#axis
range_dummy_unit <- range(gyroid$unit)
range_dummy_c <- range(gyroid$C)
dummy_unit <- seq(1, 4, length.out = 9)
dummy_c <- seq(-0.9, 1, length.out = 20)
#an analytical surface for pore
ana_pore_channel_data<-expand.grid(unit = dummy_unit,C = dummy_c,KEEP.OUT.ATTRS = F)
ana_pore_channel_data$pore_channel <-ana_pore_channel_data$unit-2*ana_pore_channel_data$unit*(atan((1+sqrt(2-ana_pore_channel_data$C*ana_pore_channel_data$C))/(ana_pore_channel_data$C+1)))/pi
ana_pore_channel_data <- acast(ana_pore_channel_data, C~unit , value.var = "pore_channel") #y ~ x
pore_fig_2 <- add_trace(p=pore_fig_2, type = "surface",
x = dummy_unit,
y = dummy_c,
z = ana_pore_channel_data,
type = "surface",
showscale = FALSE,
colorscale = list(c(0, 1), c("#f0202080", "#f0202080"))
)
#
pore_fig_2
pore_fig_3<-plot_ly(gyroid_dia, x = ~unit, y = ~C, z = ~strut_dia,
marker = list(color = ~strut_dia, colorscale = 'Greens', cmin = -0.2, cmax = 3, showscale = TRUE)) %>% add_markers()
pore_fig_3 <- pore_fig_3 %>% layout(scene = list(xaxis = list(title = 'Unit cell size [mm]',
tickvals =c(1, 2, 3, 4),
ticktext = paste0(c(1, 2, 3, 4))),
yaxis = list(title = 'C',
tickvals =c(-1.0,-0.8,-0.5,-0.3,0.0,0.3,0.5,0.8,1.0),
ticktext = paste0(c(-1.0,"",-0.5,"",0.0,"",0.5,"",1.0))),
zaxis = list(title = 'Strut diameter [mm]',
range = c(-0.2, 3.0)),
#camera perspective
camera = list(eye = list(x = -1.25, y = -1.25, z = 0.1))),
font=t,
annotations = list(
x = 1.13,
y = 1.05,
z = 1,
font=t,
text = 'Strut diameter',
xref = 'paper',
yref = 'paper',
showarrow = FALSE
))
# Create the plane based on the equation
#axis
range_dummy_unit <- range(gyroid_dia$unit)
range_dummy_c <- range(gyroid_dia$C)
dummy_unit <- seq(1, 4, length.out = 9)
dummy_c <- seq(-0.9, 1, length.out = 20)
#an analytical surface for pore
ana_strut_dia_data<-expand.grid(unit = dummy_unit,C = dummy_c,KEEP.OUT.ATTRS = F)
ana_strut_dia_data$strut_dia<-sqrt(2)*0.5*ana_strut_dia_data$unit-ana_strut_dia_data$unit+2*ana_strut_dia_data$unit*(atan((1+sqrt(2-ana_strut_dia_data$C*ana_strut_dia_data$C))/(ana_strut_dia_data$C+1)))/pi
ana_strut_dia_data <- acast(ana_strut_dia_data, C~unit , value.var = "strut_dia") #y ~ x
# add plane into 3d plot
pore_fig_3 <- add_trace(p=pore_fig_3, type = "surface",
x = dummy_unit,
y = dummy_c,
z = ana_strut_dia_data,
type = "surface",
showscale = FALSE,
colorscale = list(c(0, 1), c("#f0202080", "#f0202080"))
)
#
pore_fig_3