plot_yield_surfaces/
plot_yield_surfaces.rs1use plotpy::{linspace, Canvas, Curve, Legend, Plot, StrError, Surface};
2use std::f64::consts::PI;
3
4fn main() -> Result<(), StrError> {
5 let elev = 20.0;
6 let azim = 10.0;
7 plot_von_mises(elev, azim)?;
8 plot_mohr_coulomb(elev, azim)?;
9 Ok(())
10}
11
12fn plot_von_mises(elev: f64, azim: f64) -> Result<(), StrError> {
13 let radius = 3.0;
14 let height = 15.0;
15
16 let n_z = 50;
17 let n_theta = 50;
18
19 let z_cyl = linspace(0.0, height, n_z);
20 let theta = linspace(0.0, 2.0 * PI, n_theta);
21
22 let v = [1.0 / 3.0f64.sqrt(), 1.0 / 3.0f64.sqrt(), 1.0 / 3.0f64.sqrt()];
23 let z_axis = [0.0, 0.0, 1.0];
24
25 let mut axis = cross(&z_axis, &v);
26 let n_axis = norm(&axis);
27 for i in 0..3 {
28 axis[i] /= n_axis;
29 }
30 let angle = dot(&z_axis, &v).acos();
31
32 let r = rodrigues_rotation(&axis, angle);
33
34 let mut x_mat = vec![vec![0.0; n_theta]; n_z];
35 let mut y_mat = vec![vec![0.0; n_theta]; n_z];
36 let mut z_mat = vec![vec![0.0; n_theta]; n_z];
37
38 for i in 0..n_z {
39 for j in 0..n_theta {
40 let cx = radius * theta[j].cos();
41 let cy = radius * theta[j].sin();
42 let cz = z_cyl[i];
43
44 x_mat[i][j] = r[0][0] * cx + r[0][1] * cy + r[0][2] * cz;
45 y_mat[i][j] = r[1][0] * cx + r[1][1] * cy + r[1][2] * cz;
46 z_mat[i][j] = r[2][0] * cx + r[2][1] * cy + r[2][2] * cz;
47 }
48 }
49
50 let mut surf = Surface::new();
51 surf.set_surf_color("#7e7ed0bd")
52 .set_with_wireframe(false)
53 .draw(&x_mat, &y_mat, &z_mat);
54
55 let mut axis_line = Curve::new();
56 axis_line
57 .set_line_color("red")
58 .set_line_style("--")
59 .set_line_width(5.0)
60 .set_label("Hydrostatic Axis");
61 axis_line.draw_3d(&[0.0, 15.0], &[0.0, 15.0], &[0.0, 15.0]);
62
63 let mut legend = Legend::new();
64 legend.set_handle_len(7.0).set_location("lower left").draw();
65
66 let mut plot = Plot::new();
67 plot.add(&surf)
68 .add(&axis_line)
69 .add(&legend)
70 .set_camera(elev, azim)
71 .set_equal_axes(true)
72 .set_zoom_3d(2.1)
73 .set_hide_axes(true)
74 .set_title("von Mises Yield Surface (Cylinder)")
75 .set_figure_size_points(800.0, 800.0);
76
77 add_glyph_and_legend(&mut plot);
78
79 plot.save("/tmp/plotpy/von_mises_surface.png")?;
80
81 Ok(())
82}
83
84fn plot_mohr_coulomb(elev: f64, azim: f64) -> Result<(), StrError> {
85 let apex = [-2.0, -2.0, -2.0];
86
87 let r_c = 6.0;
88 let r_e = 3.5;
89
90 let angles = [0.0, 60.0, 120.0, 180.0, 240.0, 300.0];
91 let radii = [r_c, r_e, r_c, r_e, r_c, r_e];
92
93 let base_center = [12.0, 12.0, 12.0];
94
95 let n = [1.0 / 3.0f64.sqrt(), 1.0 / 3.0f64.sqrt(), 1.0 / 3.0f64.sqrt()];
96 let u = [1.0 / 2.0f64.sqrt(), -1.0 / 2.0f64.sqrt(), 0.0];
97 let v_vec = cross(&n, &u);
98
99 let mut xx = vec![apex[0]];
101 let mut yy = vec![apex[1]];
102 let mut zz = vec![apex[2]];
103
104 for i in 0..6 {
105 let ang = angles[i] * PI / 180.0;
106 let r = radii[i];
107 let px = base_center[0] + r * (ang.cos() * u[0] + ang.sin() * v_vec[0]);
108 let py = base_center[1] + r * (ang.cos() * u[1] + ang.sin() * v_vec[1]);
109 let pz = base_center[2] + r * (ang.cos() * u[2] + ang.sin() * v_vec[2]);
110 xx.push(px);
111 yy.push(py);
112 zz.push(pz);
113 }
114
115 let mut connectivity = vec![vec![0; 3]; 6];
117 for i in 0..6 {
118 connectivity[i][0] = 0;
119 connectivity[i][1] = i + 1;
120 connectivity[i][2] = if i == 5 { 1 } else { i + 2 };
121 }
122
123 let mut canvas = Canvas::new();
124 canvas
125 .set_face_color("#e8b25550")
126 .set_edge_color("black")
127 .set_line_width(2.0)
128 .draw_triangles_3d(&xx, &yy, &zz, &connectivity);
129
130 let mut axis_line = Curve::new();
131 axis_line
132 .set_line_color("red")
133 .set_line_style("--")
134 .set_line_width(5.0)
135 .set_label("Hydrostatic Axis");
136 axis_line.draw_3d(&[-3.0, 17.0], &[-3.0, 17.0], &[-3.0, 17.0]);
137
138 let mut legend = Legend::new();
139 legend.set_handle_len(7.0).set_location("lower left").draw();
140
141 let mut plot = Plot::new();
142 plot.add(&canvas)
143 .add(&axis_line)
144 .add(&legend)
145 .set_equal_axes(true)
146 .set_camera(elev, azim)
147 .set_zoom_3d(2.0)
148 .set_hide_axes(true)
149 .set_title("Mohr-Coulomb Yield Surface (Hexagonal Pyramid)")
150 .set_figure_size_points(800.0, 800.0);
151
152 add_glyph_and_legend(&mut plot);
153
154 plot.save("/tmp/plotpy/mohr_coulomb_surface.png")?;
155
156 Ok(())
157}
158
159fn cross(a: &[f64; 3], b: &[f64; 3]) -> [f64; 3] {
162 [
163 a[1] * b[2] - a[2] * b[1],
164 a[2] * b[0] - a[0] * b[2],
165 a[0] * b[1] - a[1] * b[0],
166 ]
167}
168
169fn dot(a: &[f64; 3], b: &[f64; 3]) -> f64 {
170 a[0] * b[0] + a[1] * b[1] + a[2] * b[2]
171}
172
173fn norm(a: &[f64; 3]) -> f64 {
174 dot(a, a).sqrt()
175}
176
177fn rodrigues_rotation(axis: &[f64; 3], angle: f64) -> [[f64; 3]; 3] {
178 let [kx, ky, kz] = *axis;
179 let k = [[0.0, -kz, ky], [kz, 0.0, -kx], [-ky, kx, 0.0]];
180 let mut k2 = [[0.0; 3]; 3];
181 for i in 0..3 {
182 for j in 0..3 {
183 for m in 0..3 {
184 k2[i][j] += k[i][m] * k[m][j];
185 }
186 }
187 }
188
189 let cos_a = angle.cos();
190 let sin_a = angle.sin();
191
192 let mut r = [[0.0; 3]; 3];
193 for i in 0..3 {
194 for j in 0..3 {
195 r[i][j] = (if i == j { 1.0 } else { 0.0 }) + sin_a * k[i][j] + (1.0 - cos_a) * k2[i][j];
196 }
197 }
198 r
199}
200
201fn add_glyph_and_legend(plot: &mut Plot) {
202 let mut glyph = Canvas::new();
203 glyph
204 .set_glyph_label_x(r"$-\\hat{\\sigma}_1$")
205 .set_glyph_label_y(r"$-\\hat{\\sigma}_2$")
206 .set_glyph_label_z(r"$-\\hat{\\sigma}_3$")
207 .set_glyph_label_color("")
208 .set_glyph_color_x("#000000")
209 .set_glyph_color_y("#000000")
210 .set_glyph_color_z("#000000")
211 .set_glyph_line_width(4.0)
212 .set_glyph_label_fontsize(16.0)
213 .set_glyph_bbox("boxstyle='circle,pad=0.2',facecolor='white',edgecolor='black'")
214 .set_glyph_size(10.0)
215 .draw_glyph_3d(0.0, 0.0, 0.0);
216
217 plot.add(&glyph);
218}