Creating an animation with d3.js

54 minute read

Introduction

In Part 8 of the series on return algorithms for plasticity we animated the closest-point return algorithm as shown below.

An animated GIF of the plot can be generated using the following R script.

require("ggplot2")
require("animation")
require("latex2exp")

#------------------------------------------------------
# Plot single iteration
#------------------------------------------------------
plotYieldSurface <- function(zrprime_data) {

  zrprime_trial_closest = zrprime_data[which(zrprime_data$Label == "trial" |
                                             zrprime_data$Label == "closest"),]
  plt = ggplot() + 
        geom_path(data = zrprime_data, 
                  aes(x = z*1.0e-6, y = rprime*1.0e-6, group=Label, color = Label),
                      size = 1)+
        geom_point(data = zrprime_data, 
                   aes(x = z*1.0e-6, y = rprime*1.0e-6, group=Label, color = Label)) +
        geom_line(data = zrprime_trial_closest,
                  aes(x = z*1.0e-6, y = rprime*1.0e-6),
                  color = "red", linetype = 1,
                  size = 1) +
        xlab(TeX(paste("$z = I_1/\\sqrt{3}$", "MPa"))) +
        ylab(TeX(paste("$r' = \\beta\\sqrt{3K/2G}\\sqrt{2J_2}$", "MPa"))) +
        coord_fixed() +
        theme_bw() + 
        theme(legend.justification=c(0,1), legend.position=c(0,1),
            plot.title = element_text(size = 14),
            axis.title.x = element_text(size=16),
            axis.title.y = element_text(size=16),
            axis.text.x = element_text(size=14),
            axis.text.y = element_text(size=14),
            legend.text = element_text(size=12))
  print(plt)
}

#------------------------------------------------------
# Animate the plots
#------------------------------------------------------
animateIterations <- function(yieldSurfData, num_iterations, consistency_iter) {
  lapply(seq(1, num_iterations, 1), 
         function(iteration) {
           plotYieldSurface(extractIteration(yieldSurfData, iteration, consistency_iter))
         })
}

#------------------------------------------------------
# Function to create animated gif
#------------------------------------------------------
createGIFIterations <- function(yieldSurfData, consistency_iter) {

  outputGIFFile = paste0("closest_point_iter", ".gif")

  # Compute number of iterations
  num_iterations = length(unique(yieldSurfData$Iteration))

  # Save as animation
  ani.options(ani.height = 600, ani.width = 600)
  saveGIF(animateIterations(yieldSurfData, num_iterations, consistency_iter), 
          interval=1.0, 
          movie.name=outputGIFFile)
}

#-------------------------------------------------------------------------
# Compute the full yield surface and create data frame
#-------------------------------------------------------------------------
ComputeFullYieldSurface <- function(yieldParams, capX, pbar_w, K, G, num_points,
                                    z_r_pt, z_r_closest, z_r_yield_z, z_r_yield_r,
                                    iteration, consistency_iter) {

  # Get yield parameters
  yieldParams = unlist(yieldParams)
  PEAKI1 = as.numeric(yieldParams['PEAKI1'])
  FSLOPE = as.numeric(yieldParams['FSLOPE'])
  STREN  = as.numeric(yieldParams['STREN'])
  YSLOPE = as.numeric(yieldParams['YSLOPE'])
  CR     = as.numeric(yieldParams['CR'])

  # Set up constants
  a1 = STREN
  a2 = (FSLOPE-YSLOPE)/(STREN-YSLOPE*PEAKI1)
  a3 = (STREN-YSLOPE*PEAKI1)*exp(-a2*PEAKI1)
  a4 = YSLOPE

  # Compute kappa
  X_eff = capX + 3*pbar_w
  kappa = PEAKI1 - CR*(PEAKI1 - X_eff)

  # Create an array of I1_eff values
  I1eff_min = X_eff;
  I1eff_max = PEAKI1;

  rad = 0.5*(PEAKI1 - X_eff)
  cen = 0.5*(PEAKI1 + X_eff)
  theta_max = acos(max((I1eff_min - cen)/rad, -1.0))
  theta_min = acos(min((I1eff_max - cen)/rad, 1.0))
  theta_vec = seq(from = theta_min, to = theta_max, length.out = num_points)
  I1_list = cen + rad*cos(theta_vec)
  I1_eff_list = lapply(I1_list, function(val) {max(val, X_eff)})
  sqrtJ2_list = lapply(I1_eff_list,
    function(I1_eff) {
      # Compute F_f
      Ff = a1 - a3*exp(a2*I1_eff) - a4*(I1_eff)

      # Compute Fc
      Fc_sq = 1.0
      if ((I1_eff < kappa) && (X_eff <= I1_eff)) {
        ratio = (kappa - I1_eff)/(kappa - X_eff)
        Fc_sq = 1.0 - ratio^2
      }

      # Compute sqrt(J2)
      sqrtJ2 = Ff*sqrt(Fc_sq)
      return(sqrtJ2)
    })
  zrprime_data = data.frame(z = unlist(I1_eff_list)/sqrt(3), 
                            rprime = unlist(sqrtJ2_list)*sqrt(2)*sqrt(1.5*K/G), 
                            Iteration = as.factor(iteration),
                            CIteration = as.factor(consistency_iter),
                            Label="yield")
  zrprime_data = rbind(zrprime_data,
                       data.frame(z = z_r_pt[1], rprime = z_r_pt[2], 
                                  Iteration = as.factor(iteration),
                                  CIteration = as.factor(consistency_iter),
                                  Label = "trial"))
  zrprime_data = rbind(zrprime_data,
                       data.frame(z = z_r_closest[1], rprime = z_r_closest[2], 
                                  Iteration = as.factor(iteration),
                                  CIteration = as.factor(consistency_iter),
                                  Label = "closest"))
  zrprime_data = rbind(zrprime_data,
                       data.frame(z = z_r_yield_z, rprime = z_r_yield_r, 
                                  Iteration = as.factor(iteration),
                                  CIteration = as.factor(consistency_iter),
                                  Label = "polyline"))
 
  return(zrprime_data)
}

#-------------------------------------------------------------------------
# Read data and plot for animation
#-------------------------------------------------------------------------
ReadAndPlotClosestPoint <- function() {)

  num_points = 50
  consistency_iter = 1

  iteration = 1
  K = 3.98447e+07
  G = 1.32816e+07
  X = -1.87891e+07
  pbar_w = 5.24774e+06
  yieldParams = list(BETA = 1, CR = 0.5, FSLOPE = 0.355309, PEAKI1 = 996.118, STREN = 1.76382e+07, YSLOPE = 0.353622)
  z_r_pt = c(52172.3,3.60195e+06)
  z_r_closest = c(-1.15079e+06,2.02165e+06)
  z_r_yield_z = c(575.109,-167406,-607187,-1.15079e+06,-1.59057e+06,-1.75855e+06)
  z_r_yield_r = c(2.45257e-09,310133,1.12207e+06,2.02165e+06,1.72669e+06,0)
  zr_df = 
  ComputeFullYieldSurface(yieldParams, X, pbar_w, K, G, num_points,
                          z_r_pt, z_r_closest, z_r_yield_z, z_r_yield_r,
                          iteration, consistency_iter)
  iteration = 2
  K = 3.98447e+07
  G = 1.32816e+07
  X = -1.87891e+07
  pbar_w = 5.24774e+06
  yieldParams = list(BETA = 1, CR = 0.5, FSLOPE = 0.355309, PEAKI1 = 996.118, STREN = 1.76382e+07, YSLOPE = 0.353622)
  z_r_pt = c(52172.3,3.60195e+06)
  z_r_closest = c(-1.15079e+06,2.02165e+06)
  z_r_yield_z = c(-878986,-1.15079e+06,-1.39598e+06,-1.59057e+06,-1.7155e+06,-1.75855e+06)
  z_r_yield_r = c(1.62388e+06,2.02165e+06,2.08595e+06,1.72669e+06,979052,0)
  zr_df = rbind(zr_df,
  ComputeFullYieldSurface(yieldParams, X, pbar_w, K, G, num_points,
                          z_r_pt, z_r_closest, z_r_yield_z, z_r_yield_r,
                          iteration, consistency_iter))
  iteration = 3
  K = 3.98447e+07
  G = 1.32816e+07
  X = -1.87891e+07
  pbar_w = 5.24774e+06
  yieldParams = list(BETA = 1, CR = 0.5, FSLOPE = 0.355309, PEAKI1 = 996.118, STREN = 1.76382e+07, YSLOPE = 0.353622)
  z_r_pt = c(52172.3,3.60195e+06)
  z_r_closest = c(-1.15079e+06,2.02165e+06)
  z_r_yield_z = c(-878986,-970925,-1.06186e+06,-1.15079e+06,-1.23674e+06,-1.31877e+06)
  z_r_yield_r = c(1.62388e+06,1.7838e+06,1.91864e+06,2.02165e+06,2.08688e+06,2.10948e+06)
  zr_df = rbind(zr_df,
  ComputeFullYieldSurface(yieldParams, X, pbar_w, K, G, num_points,
                          z_r_pt, z_r_closest, z_r_yield_z, z_r_yield_r,
                          iteration, consistency_iter))
  iteration = 4
  K = 3.98447e+07
  G = 1.32816e+07
  X = -1.87891e+07
  pbar_w = 5.24774e+06
  yieldParams = list(BETA = 1, CR = 0.5, FSLOPE = 0.355309, PEAKI1 = 996.118, STREN = 1.76382e+07, YSLOPE = 0.353622)
  z_r_pt = c(52172.3,3.60195e+06)
  z_r_closest = c(-1.18969e+06,2.05584e+06)
  z_r_yield_z = c(-1.09888e+06,-1.14468e+06,-1.18969e+06,-1.2338e+06,-1.27687e+06,-1.31877e+06)
  z_r_yield_r = c(1.96539e+06,2.01563e+06,2.05584e+06,2.0853e+06,2.10336e+06,2.10948e+06)
  zr_df = rbind(zr_df,
  ComputeFullYieldSurface(yieldParams, X, pbar_w, K, G, num_points,
                          z_r_pt, z_r_closest, z_r_yield_z, z_r_yield_r,
                          iteration, consistency_iter))
  iteration = 5
  K = 3.98447e+07
  G = 1.32816e+07
  X = -1.87891e+07
  pbar_w = 5.24774e+06
  yieldParams = list(BETA = 1, CR = 0.5, FSLOPE = 0.355309, PEAKI1 = 996.118, STREN = 1.76382e+07, YSLOPE = 0.353622)
  z_r_pt = c(52172.3,3.60195e+06)
  z_r_closest = c(-1.18723e+06,2.05389e+06)
  z_r_yield_z = c(-1.09888e+06,-1.12123e+06,-1.14342e+06,-1.16542e+06,-1.18723e+06,-1.20882e+06)
  z_r_yield_r = c(1.96539e+06,1.99102e+06,2.01438e+06,2.03536e+06,2.05389e+06,2.06989e+06)
  zr_df = rbind(zr_df,
  ComputeFullYieldSurface(yieldParams, X, pbar_w, K, G, num_points,
                          z_r_pt, z_r_closest, z_r_yield_z, z_r_yield_r,
                          iteration, consistency_iter))
  iteration = 6
  K = 3.98447e+07
  G = 1.32816e+07
  X = -1.87891e+07
  pbar_w = 5.24774e+06
  yieldParams = list(BETA = 1, CR = 0.5, FSLOPE = 0.355309, PEAKI1 = 996.118, STREN = 1.76382e+07, YSLOPE = 0.353622)
  z_r_pt = c(52172.3,3.60195e+06)
  z_r_closest = c(-1.18699e+06,2.05371e+06)
  z_r_yield_z = c(-1.15385e+06,-1.16495e+06,-1.176e+06,-1.18699e+06,-1.19794e+06,-1.20882e+06)
  z_r_yield_r = c(2.0246e+06,2.03493e+06,2.04464e+06,2.05371e+06,2.06213e+06,2.06989e+06)
  zr_df = rbind(zr_df,
  ComputeFullYieldSurface(yieldParams, X, pbar_w, K, G, num_points,
                          z_r_pt, z_r_closest, z_r_yield_z, z_r_yield_r,
                          iteration, consistency_iter))
  iteration = 7
  K = 3.98447e+07
  G = 1.32816e+07
  X = -1.87891e+07
  pbar_w = 5.24774e+06
  yieldParams = list(BETA = 1, CR = 0.5, FSLOPE = 0.355309, PEAKI1 = 996.118, STREN = 1.76382e+07, YSLOPE = 0.353622)
  z_r_pt = c(52172.3,3.60195e+06)
  z_r_closest = c(-1.18686e+06,2.0536e+06)
  z_r_yield_z = c(-1.18133e+06,-1.18686e+06,-1.19237e+06,-1.19787e+06,-1.20335e+06,-1.20882e+06)
  z_r_yield_r = c(2.04911e+06,2.0536e+06,2.05792e+06,2.06208e+06,2.06607e+06,2.06989e+06)
  zr_df = rbind(zr_df,
  ComputeFullYieldSurface(yieldParams, X, pbar_w, K, G, num_points,
                          z_r_pt, z_r_closest, z_r_yield_z, z_r_yield_r,
                          iteration, consistency_iter))
  iteration = 8
  K = 3.98447e+07
  G = 1.32816e+07
  X = -1.87891e+07
  pbar_w = 5.24774e+06
  yieldParams = list(BETA = 1, CR = 0.5, FSLOPE = 0.355309, PEAKI1 = 996.118, STREN = 1.76382e+07, YSLOPE = 0.353622)
  z_r_pt = c(52172.3,3.60195e+06)
  z_r_closest = c(-1.18684e+06,2.05359e+06)
  z_r_yield_z = c(-1.18133e+06,-1.18409e+06,-1.18684e+06,-1.18959e+06,-1.19234e+06,-1.19508e+06)
  z_r_yield_r = c(2.04911e+06,2.05137e+06,2.05359e+06,2.05576e+06,2.05789e+06,2.05999e+06)
  zr_df = rbind(zr_df,
  ComputeFullYieldSurface(yieldParams, X, pbar_w, K, G, num_points,
                          z_r_pt, z_r_closest, z_r_yield_z, z_r_yield_r,
                          iteration, consistency_iter))
  iteration = 9
  K = 3.98447e+07
  G = 1.32816e+07
  X = -1.87891e+07
  pbar_w = 5.24774e+06
  yieldParams = list(BETA = 1, CR = 0.5, FSLOPE = 0.355309, PEAKI1 = 996.118, STREN = 1.76382e+07, YSLOPE = 0.353622)
  z_r_pt = c(52172.3,3.60195e+06)
  z_r_closest = c(-1.18683e+06,2.05358e+06)
  z_r_yield_z = c(-1.18133e+06,-1.18271e+06,-1.18409e+06,-1.18546e+06,-1.18683e+06,-1.18821e+06)
  z_r_yield_r = c(2.04911e+06,2.05025e+06,2.05137e+06,2.05248e+06,2.05358e+06,2.05467e+06)
  zr_df = rbind(zr_df,
  ComputeFullYieldSurface(yieldParams, X, pbar_w, K, G, num_points,
                          z_r_pt, z_r_closest, z_r_yield_z, z_r_yield_r,
                          iteration, consistency_iter))
  iteration = 10
  K = 3.98447e+07
  G = 1.32816e+07
  X = -1.87891e+07
  pbar_w = 5.24774e+06
  yieldParams = list(BETA = 1, CR = 0.5, FSLOPE = 0.355309, PEAKI1 = 996.118, STREN = 1.76382e+07, YSLOPE = 0.353622)
  z_r_pt = c(52172.3,3.60195e+06)
  z_r_closest = c(-1.18658e+06,2.05338e+06)
  z_r_yield_z = c(-1.18477e+06,-1.18546e+06,-1.18615e+06,-1.18683e+06,-1.18752e+06,-1.18821e+06)
  z_r_yield_r = c(2.05192e+06,2.05248e+06,2.05303e+06,2.05358e+06,2.05412e+06,2.05467e+06)
  zr_df = rbind(zr_df,
  ComputeFullYieldSurface(yieldParams, X, pbar_w, K, G, num_points,
                          z_r_pt, z_r_closest, z_r_yield_z, z_r_yield_r,
                          iteration, consistency_iter))
  iteration = 11
  K = 3.98447e+07
  G = 1.32816e+07
  X = -1.87891e+07
  pbar_w = 5.24774e+06
  yieldParams = list(BETA = 1, CR = 0.5, FSLOPE = 0.355309, PEAKI1 = 996.118, STREN = 1.76382e+07, YSLOPE = 0.353622)
  z_r_pt = c(52172.3,3.60195e+06)
  z_r_closest = c(-1.18649e+06,2.0533e+06)
  z_r_yield_z = c(-1.18649e+06,-1.18683e+06,-1.18718e+06,-1.18752e+06,-1.18786e+06,-1.18821e+06)
  z_r_yield_r = c(2.0533e+06,2.05358e+06,2.05385e+06,2.05412e+06,2.0544e+06,2.05467e+06)
  zr_df = rbind(zr_df,
  ComputeFullYieldSurface(yieldParams, X, pbar_w, K, G, num_points,
                          z_r_pt, z_r_closest, z_r_yield_z, z_r_yield_r,
                          iteration, consistency_iter))
  iteration = 12
  K = 3.98447e+07
  G = 1.32816e+07
  X = -1.87891e+07
  pbar_w = 5.24774e+06
  yieldParams = list(BETA = 1, CR = 0.5, FSLOPE = 0.355309, PEAKI1 = 996.118, STREN = 1.76382e+07, YSLOPE = 0.353622)
  z_r_pt = c(52172.3,3.60195e+06)
  z_r_closest = c(-1.18649e+06,2.0533e+06)
  z_r_yield_z = c(-1.18649e+06,-1.18666e+06,-1.18683e+06,-1.187e+06,-1.18718e+06,-1.18735e+06)
  z_r_yield_r = c(2.0533e+06,2.05344e+06,2.05358e+06,2.05371e+06,2.05385e+06,2.05399e+06)
  zr_df = rbind(zr_df,
  ComputeFullYieldSurface(yieldParams, X, pbar_w, K, G, num_points,
                          z_r_pt, z_r_closest, z_r_yield_z, z_r_yield_r,
                          iteration, consistency_iter))
  createGIFIterations(zr_df, consistency_iter)
}

#-------------------------------------------------------------------------
# Call function to read data and plot
#-------------------------------------------------------------------------
ReadAndPlotClosestPoint()

Let us explore how the animation at the top of the docs was produced using d3.js.

Input data

The input data that can be seen in the R function ReadAndPlotClosestPoint were generated during a simulation using print statements. These have to be converted into a form that can be read easily by Javascript. For our simulation we converted these into JSON and created the file yieldSurfData.json that contains:

{
  "num_points":50,
  "consistency_iter":1,
  "data":[
    {
      "iteration":1,
      "K":3.98447e+07,
      "G":1.32816e+07,
      "X":-1.87891e+07,
      "pbar_w":5.24774e+06,
      "BETA":1,
      "CR":0.5,
      "FSLOPE":0.355309,
      "PEAKI1":996.118,
      "STREN":1.76382e+07,
      "YSLOPE":0.353622,
      "z_r_pt":[
        52172.3,
        3.60195e+06
      ],
      "z_r_closest":[
        -1.15079e+06,
        2.02165e+06
      ],
      "z_r_yield_z":[
        575.109,
        -167406,
        -607187,
        -1.15079e+06,
        -1.59057e+06,
        -1.75855e+06
      ],
      "z_r_yield_r":[
        2.45257e-09,
        310133,
        1.12207e+06,
        2.02165e+06,
        1.72669e+06,
        0
      ]
    },
    .....
  ]
}

The HTML file

In the HTML file, I added the following to allow the animation to be added to the DOM.

  <body>
     ...
     <!-- Create div where the svg canvas will be added -->
     <div class="yield-surf-canvas"> </div>
     <!-- Load D3.js -->
     <script src="https://d3js.org/d3.v4.min.js"></script>
     <!-- Load javascript for animating yield surface -->
     <script src="/ParSim/assets/js/yieldsurface.js"></script>
     <script>
       // Read JSON file and then call javascript function to animate yield surface 
       d3.json("/ParSim/assets/json/yieldSurfData.json", drawYieldSurface);
     </script>
  </body>

The Javascript code

Let us now explore the Javascript code used to generate the animation.

Generating points on the yield surface

We first generate points on the yield surface in \(z-r'\)-space. These are displayed in blue in the figure. The function computeYieldSurface takes an input the data for a single iteration (read from the JSON file) and the number of points to be used to discretize the yield surface.

function computeYieldSurface(iterData, numPts) {

  // Get the moduli
  let K = iterData['K'];
  let G = iterData['G'];

  // Get yield parameters
  let PEAKI1 = iterData['PEAKI1'];
  let FSLOPE = iterData['FSLOPE'];
  let STREN  = iterData['STREN'];
  let YSLOPE = iterData['YSLOPE'];
  let CR     = iterData['CR'];

  // Set up constants
  let a1 = STREN;
  let a2 = (FSLOPE-YSLOPE)/(STREN-YSLOPE*PEAKI1);
  let a3 = (STREN-YSLOPE*PEAKI1)*Math.exp(-a2*PEAKI1);
  let a4 = YSLOPE;

  // Compute kappa
  let X_eff = iterData['X'] + 3.0*iterData['pbar_w'];
  let kappa = PEAKI1 - CR*(PEAKI1 - X_eff);

  // Create an array of I1_eff values
  let I1eff_min = X_eff;
  let I1eff_max = PEAKI1;

  let rad = 0.5*(PEAKI1 - X_eff);
  let cen = 0.5*(PEAKI1 + X_eff);
  let theta_max = Math.acos(Math.max((I1eff_min - cen)/rad, -1.0));
  let theta_min = Math.acos(Math.min((I1eff_max - cen)/rad, 1.0));
  let theta_inc = (theta_max-theta_min)/numPts;
  let theta_vec = d3.range(theta_min, theta_max+theta_inc, theta_inc);
  let I1_list = theta_vec.map(function(theta) {return cen + rad*Math.cos(theta);});
  let I1_eff_list = I1_list.map(function(I1) {return Math.max(I1, X_eff);});

  // Create an array of sqrtJ2 values
  let sqrtJ2_list = I1_eff_list.map(
          function(I1_eff) {
            // Compute F_f
            let Ff = a1 - a3*Math.exp(a2*I1_eff) - a4*(I1_eff);

            // Compute Fc
            let Fc_sq = 1.0;
            if ((I1_eff < kappa) && (X_eff <= I1_eff)) {
              let ratio = (kappa - I1_eff)/(kappa - X_eff);
              Fc_sq = 1.0 - ratio*ratio;
            }

            // Compute sqrt(J2)
            let sqrtJ2 = Ff*Math.sqrt(Fc_sq);
            return(sqrtJ2);
          });

  // Compute z and r'
  let z_list = I1_eff_list.map(function(I1_eff) {return I1_eff/Math.sqrt(3.0);});
  let rprime_list = sqrtJ2_list.map(function(sqrtJ2) {return sqrtJ2*Math.sqrt(3.0*K/G);});

  // Create zipped array
  let zrprime_data = d3.zip(z_list, rprime_list);
  return zrprime_data;
}

The procedure is identical to that used in the R script, except for two differences:

  • the d3.range function is used to create the vector of \(\theta\) values. This function differs from the R function seq in that the end point of the range is not included in the sequence that is produced. Therefore, the set of points produced by d3.js is not identical to that produced by R.
  • the function returns a zipped array created with d3.zip instead of a R dataframe.
Drawing the yield surface and closest-point projection

We use the drawYieldSurface function to plot and animate the yield surface and the closest-point projection. Let us look at the full code first and then explore some of the details.

function drawYieldSurface(data) {

  // Create svg
  var svgsize = {x: 600, y: 600};
  var margin = {top: 50, right: 50, bottom: 50, left: 50};
  var width = svgsize.x - margin.left - margin.right;
  var height = svgsize.y - margin.top - margin.bottom;
  var svg = d3.select(".yield-surf-canvas")
              .append("svg")
                 .attr("width", svgsize.x)
                 .attr("height", svgsize.y);

  // Create chart area
  var chart = svg.append("g")
                    .attr("class", "chart")
                    .attr("transform", 
                          "translate(" + margin.left + "," + margin.top + ")");

  // Process the input data
  let numPts = data['num_points'];
  let numIter = data.data.length;

  // Compute yield surface in z and r' from first element of data
  // (Note: the yield surface does not change at each iteration)
  let iterData = data.data[0];
  let zrprime = computeYieldSurface(iterData, numPts);

  // Get the coordinates of the trial stress
  let ztrial = iterData.z_r_pt[0];
  let rprimetrial = iterData.z_r_pt[1];

  // Compute min max of domain
  let zmin = d3.min(zrprime, function(d) {return d[0];})
  let zmax = d3.max(zrprime, function(d) {return d[0];})
  let rprimemin = d3.min(zrprime, function(d) {return d[1];})
  let rprimemax = d3.max(zrprime, function(d) {return d[1];})
  zmin = Math.min(zmin, ztrial);
  zmax = Math.max(zmax, ztrial);
  rprimemin = Math.min(rprimemin, rprimetrial);
  rprimemax = Math.max(rprimemax, rprimetrial);

  // Create scaling functions
  var xscale = d3.scaleLinear().domain([zmin, zmax]).rangeRound([0, width]);
  var yscale = d3.scaleLinear().domain([rprimemin, rprimemax]).rangeRound([height, 0]);
  var xscaleAxis = d3.scaleLinear().domain([zmin*1.0e-6, zmax*1.0e-6]).rangeRound([0, width]);
  var yscaleAxis = d3.scaleLinear().domain([rprimemin*1.0e-6, rprimemax*1.0e-6]).rangeRound([height, 0]);

  // Create polyline generation function
  var line = d3.line()
                 .x(function(d) { return xscale(d[0]); })
                 .y(function(d) { return yscale(d[1]); });

  // Create bottom axis
  chart.append("g")
          .attr("transform", "translate(0," + height + ")")
        .call(d3.axisBottom(xscaleAxis))
        .append("text")
           .attr("fill", "#000")
           .attr("transform",
                 "translate(" + (width/2) + " ," + 
                                (margin.top - 20) + ")")
           .attr("text-anchor", "end")
           .text("z (MPa)"); 

  // Create left axis
  chart.append("g")
       .call(d3.axisLeft(yscaleAxis))
       .append("text")
          .attr("fill", "#000")
          .attr("transform", "rotate(-90)")
          .attr("y", 6)
          .attr("dy", "0.71em")
          .attr("text-anchor", "end")
          .text("r' (MPa)"); 

  // Plot yield surface
  chart.append("path")
          .attr("class", "yield_surf")
          .datum(zrprime)
          .attr("fill", "none")    
          .attr("stroke", "steelblue")
          .attr("stroke-linejoin", "round")
          .attr("stroke-linecap", "round")
          .attr("stroke-width", 1.5)
          .attr("d", line);

  // Loop through the iterations
  for (let iter = 0; iter < numIter; iter++) {

    iterData = data.data[iter];

    // Create group for each iteration
    let iterGroup = chart.append("g")
                            .attr("class", "iteration"+iter)
                            .attr("opacity", 0);

    // Create group for the yield surface points + circles
    let yieldSurfGroup = iterGroup.append("g") 
                                     .attr("class", "approx_yield_surf"+iter);

    // Get the coordinates of polyline approximation
    let zpoly = iterData.z_r_yield_z;
    let rprimepoly = iterData.z_r_yield_r;
    let zrprimepoly = d3.zip(zpoly, rprimepoly);

    // Add the yield surface polyline to the svg
    yieldSurfGroup.append("path")
                  .datum(zrprimepoly)
                     .attr("fill", "none")    
                     .attr("stroke", "red")
                     .attr("stroke-linejoin", "round")
                     .attr("stroke-linecap", "round")
                     .attr("stroke-width", 1.5)
                     .attr("d", line);

    // Add circles to the yield surface polyline
    yieldSurfGroup.selectAll(".approx_yield_surf_circle"+iter)
                  .data(zrprimepoly)
                  .enter()
                  .append("circle")
                     .attr("r", 2)
                     .attr("cx", function(d) {return xscale(d[0]);})
                     .attr("cy", function(d) {return yscale(d[1]);})
                     .attr("fill", "blue")
                     .attr("stroke", "black");

    // Create group for the projection line and points
    let projLineGroup = iterGroup.append("g") 
                                    .attr("class", "projection_line"+iter);

    // Get the coordinates of the trial stress
    let ztrial = iterData.z_r_pt[0];
    let rprimetrial = iterData.z_r_pt[1];

    // Get the coordinates of the closest point stress
    let zclosest = iterData.z_r_closest[0];
    let rprimeclosest = iterData.z_r_closest[1];

    // Set up projection line
    let zrprimeproj = [[ztrial, rprimetrial],[zclosest, rprimeclosest]];

    // Add the projection line to the svg
    projLineGroup.append("path")
                 .datum(zrprimeproj)
                   .attr("fill", "none")    
                   .attr("stroke", "green")
                   .attr("stroke-linejoin", "round")
                   .attr("stroke-linecap", "round")
                   .attr("stroke-width", 1.5)
                   .attr("d", line);

    // Add the triangles at the end points of the projection line
    projLineGroup.selectAll(".projection_line_point"+iter)
                 .data(zrprimeproj)
                 .enter()
                 .append("path")
                    .attr("class", "projection_line_point"+iter)
                    .attr("d", d3.symbol().type(d3.symbolTriangle))
                    .attr("transform", 
                      function(d) {return "translate(" + xscale(d[0]) + "," + yscale(d[1]) + ")";});

  } // End for loop over iterations

  // Function to do the animation
  function drawAnimation() {
    var iterationID = 0;

    // Function to draw the surface for the next iteration
    function drawNextSurf() {

      d3.active(this)
          .attr("opacity", 1)
        .transition()
          .duration(1000)
          .attr("opacity", 0);
      iterationID++;
      if (iterationID === numIter) {
        iterationID = 0;
      }
      chart.select(".iteration"+iterationID)
           .transition()
             .delay(1000)
             .duration(1000)
           .on("start", drawNextSurf);            

    } // end function drawNextSurf

    // Start the animation
    chart.select(".iteration0")
         .transition()
           .duration(1000)
         .on("start", drawNextSurf);

  } // end function drawAnimation

  // Do the animation
  drawAnimation();

} 

Let us now look at some of the details of the code.

Creating the SVG

First we select the <div> in our HTML file using d3.select and add a SVG canvas to it:

var svg = d3.select(".yield-surf-canvas")
            .append("svg")
              .attr("width", svgsize.x)
              .attr("height", svgsize.y);
Creating the canvas group

Then we create a subset of the SVG canvas as the region where the plot will be produced and identify it using a group g. This area has to be translated so that the top left hand corner is at the right position.

var chart = svg.append("g")
                 .attr("class", "chart")
                 .attr("transform", 
                       "translate(" + margin.left + "," + margin.top + ")");
Creating map from real to canvas coordinates

We then compute the true yield surface and find the extents of the domain of the plot based on the yield surface and the trial stress state. We use these extents to create maps from real coordinates to svg coordinates. Separate maps are also created so that the axis labels can be converted into MPa.

var xscale = d3.scaleLinear()
                 .domain([zmin, zmax])
                 .rangeRound([0, width]);
var yscale = d3.scaleLinear()
                 .domain([rprimemin, rprimemax])
                 .rangeRound([height, 0]);
var xscaleAxis = d3.scaleLinear()
                     .domain([zmin*1.0e-6, zmax*1.0e-6])
                     .rangeRound([0, width]);
var yscaleAxis = d3.scaleLinear()
                     .domain([rprimemin*1.0e-6, rprimemax*1.0e-6])
                     .rangeRound([height, 0]);
Creating SVG polyline generator

To plot the polylines that represent the yield surface, its approximation, and the closest point projection, we need a function that can convert a set of points to the SVG representation of a polyline. We generate that function using d3.line:

var line = d3.line()
               .x(function(d) { return xscale(d[0]); })
               .y(function(d) { return yscale(d[1]); });
Creating the axes and the yield surface

Next, we create the axes. The generation of axes using d3.js is straightforward, but we have to keep in mind that a translation may be required depending on the position of the axis.

After the axes have been created, we add the polyline representing the yield surface to the SVG. This is done using the zipped data returned from the computeYieldSurface function. Because the yield surface does not change during closest-point search iterations, we use the data from the first iteration to plot this curve.

chart.append("path")
       .attr("class", "yield_surf")
     .datum(zrprime)
       .attr("fill", "none")    
       .attr("stroke", "steelblue")
       .attr("stroke-width", 1.5)
       .attr("d", line);

Notice that we use the .datum method in this case and that the attribute d is set using the line function we had defined earlier.

Creating groups for each iteration

Now that the fixed content of the plot has been created, we are ready to create the content that changes between iterations. To do that, we loop through the iterations and add a group for each iteration:

let iterGroup = chart.append("g")
                       .attr("class", "iteration"+iter)
                       .attr("opacity", 0);

We set the opacity of the group to 0 so that although all the objects in the plot have been created, none are visible when the animation starts.

Adding circles to the approximate surface

Now we are ready to add in the polyline approximation to the yield surface in red. The procedure is identical to the plot of the full yield surface. We also add circles to each vertex of the polyline using

yieldSurfGroup.selectAll(".approx_yield_surf_circle"+iter)
                .data(zrprimepoly)
                .enter()
                .append("circle")
                  .attr("class", "approx_yield_surf_circle"+iter)
                  .attr("r", 2)
                  .attr("cx", function(d) {return xscale(d[0]);})
                  .attr("cy", function(d) {return yscale(d[1]);})
                  .attr("fill", "blue")
                  .attr("stroke", "black");
Adding triangles to the projection line

Next we add the line that represents the closest-point projection (in green) and add triangles to the end points for clarity. The triangles are added using

projLineGroup.selectAll(".projection_line_point"+iter)
               .data(zrprimeproj)
               .enter()
               .append("path")
                 .attr("class", "projection_line_point"+iter)
                 .attr("d", d3.symbol().type(d3.symbolTriangle))
                 .attr("transform", 
                      function(d) {return "translate(" + xscale(d[0]) + "," + yscale(d[1]) + ")";});
Setting up the animation

Now we can terminate the loop of the iterations and proceed to setting up the animation. We start the animation by selecting the group that has class iteration0 and use a transition that takes 1 second to move from one frame to the next. We call the drawNextSurf method when the animation starts.

chart.select(".iteration0")
     .transition()
       .duration(1000)
     .on("start", drawNextSurf);

The drawNextSurf method changes the opacity of the active group to 1, waits for 1 second, and then changes back the opacity to 0.

d3.active(this)
    .attr("opacity", 1)
  .transition()
    .duration(1000)
    .attr("opacity", 0);

After that, it increments to iteration count, selects the next group, and recurses after a short delay.

chart.select(".iteration"+iterationID)
     .transition()
       .delay(1000)
       .duration(1000)
     .on("start", drawNextSurf);

To allow for the animation to loop, we set the iteration index to 0 when the total number of iterations in the data is reached.

Remarks

This d3.js animation took me 2 days to complete, starting from absolutely no knowledge of the library, and with significant breaks in the workflow. I’d suggest that, since a novice can learn and use the library so quickly, d3.js should be an essential tool in the computational engineers toolbox to make data available in the form of interactive plots instead of the static (and hard to grasp) plots that are the norm in engineering literature.

In the next article, we will go back to the work on XML particle data that we had started earlier and discuss ways to reading XML files using C++ code.