carapace_core/model/elements/
shell3.rs

1//! `Shell3`: the 3-node flat shell (OpenSees `ShellDKGT`): a DKT plate-bending triangle (Kirchhoff, so no transverse
2//! shear) combined with an Allman-type membrane triangle that carries the drilling rotation, six DOFs per node
3//! (18 total), the 4-point rule `(1/3,1/3,1/3)` with weight `-9/16` and `(1/5,3/5,1/5)` permutations with `25/48`.
4//!
5//! Like [`Shell4`](super::Shell4) it reproduces OpenSees term by term: the local frame is `e1 = x1 - x0`, `e2` its
6//! Gram-Schmidt completion with `x2 - x0`, `e3 = e1 x e2`; the section ordering, sign convention and bending-row
7//! negation are the same, and the mass is a row-sum lumped translational mass (OpenSees' is consistent). The DKGT has
8//! no drilling penalty and no shear strains: its `gamma_xz`, `gamma_yz` rows are zero.
9
10use nalgebra::{SMatrix, SVector, Vector3};
11
12use super::{DofMask, DofRef, ElementForce, NodeView, ShellResponse, TangentSink, VectorSink};
13use crate::model::continuum::validate_triangle;
14use crate::model::{ElementLoad3, Node3Id, ShellSection, ShellVector, SPATIAL_NDF, SPATIAL_NDIM};
15
16type View<'a> = NodeView<'a, SPATIAL_NDIM, SPATIAL_NDF, Node3Id>;
17type Vec18 = SVector<f64, 18>;
18type Mat18 = SMatrix<f64, 18, 18>;
19type Bmat = SMatrix<f64, 8, 18>;
20
21/// Area coordinates `(L1, L2, L3)` of the four Gauss points and their weights.
22const GAUSS: [([f64; 3], f64); 4] = [
23    ([1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0], -9.0 / 16.0),
24    ([0.2, 0.6, 0.2], 25.0 / 48.0),
25    ([0.6, 0.2, 0.2], 25.0 / 48.0),
26    ([0.2, 0.2, 0.6], 25.0 / 48.0),
27];
28
29#[derive(Debug, Clone)]
30pub struct Shell3 {
31    pub nodes: [Node3Id; 3],
32    sections: [ShellSection; 4],
33    cache: Option<Cache>,
34}
35
36#[derive(Debug, Clone, Copy, PartialEq)]
37struct Cache {
38    /// Local axes `e1`, `e2`, `e3` in global components.
39    axes: [Vector3<f64>; 3],
40    points: [PointGeometry; 4],
41}
42
43#[derive(Debug, Clone, Copy, PartialEq)]
44struct PointGeometry {
45    /// Strain-displacement matrix of the generalized strains, in global DOFs.
46    b: Bmat,
47    /// The same with the bending-bending block negated, for the internal force.
48    b_force: Bmat,
49    /// `0.5 w xsj` (`w` summing to 1, so the weights sum to the area).
50    dvol: f64,
51    /// Area coordinates.
52    n: [f64; 3],
53}
54
55type Frame = ([Vector3<f64>; 3], [[f64; 3]; 2]);
56
57fn frame(coords: &[[f64; 3]; 3]) -> Result<Frame, &'static str> {
58    let x: [Vector3<f64>; 3] = std::array::from_fn(|i| Vector3::from(coords[i]));
59    let mut v1 = x[1] - x[0];
60    let mut v2 = x[2] - x[0];
61    let n1 = v1.norm();
62    if !n1.is_finite() || n1 == 0.0 {
63        return Err("degenerate shell geometry");
64    }
65    v1 /= n1;
66    v2 -= v2.dot(&v1) * v1;
67    let n2 = v2.norm();
68    if !n2.is_finite() || n2 == 0.0 {
69        return Err("degenerate shell geometry");
70    }
71    v2 /= n2;
72    let v3 = v1.cross(&v2);
73    Ok((
74        [v1, v2, v3],
75        [
76            std::array::from_fn(|i| x[i].dot(&v1)),
77            std::array::from_fn(|i| x[i].dot(&v2)),
78        ],
79    ))
80}
81
82/// The DKT bending shape-function derivatives, OpenSees' `shapeBend`: rows `Hx,x`, `Hx,y`, `Hy,x`, `Hy,y`
83/// over the nine DOFs `(w, theta1, theta2)` per node.
84fn shape_bend(l: [f64; 3], x: &[[f64; 3]; 2], area: f64) -> [[f64; 9]; 4] {
85    let (x12, x23, x31) = (x[0][0] - x[0][1], x[0][1] - x[0][2], x[0][2] - x[0][0]);
86    let (y12, y23, y31) = (x[1][0] - x[1][1], x[1][1] - x[1][2], x[1][2] - x[1][0]);
87    let (l12, l23, l31) = (
88        x12 * x12 + y12 * y12,
89        x23 * x23 + y23 * y23,
90        x31 * x31 + y31 * y31,
91    );
92    // Per mid-side 4 (23), 5 (31), 6 (12): a, b, c, d, e.
93    let side = |xs: f64, ys: f64, len: f64| {
94        (
95            -xs / len,
96            0.25 * 3.0 * xs * ys / len,
97            0.25 * (xs * xs - 2.0 * ys * ys) / len,
98            -ys / len,
99            0.25 * (ys * ys - 2.0 * xs * xs) / len,
100        )
101    };
102    let (a4, b4, _c4, d4, _e4) = side(x23, y23, l23);
103    let (a5, b5, _c5, d5, _e5) = side(x31, y31, l31);
104    let (a6, b6, _c6, d6, _e6) = side(x12, y12, l12);
105    let (p4, p5, p6) = (6.0 * a4, 6.0 * a5, 6.0 * a6);
106    let (t4, t5, t6) = (6.0 * d4, 6.0 * d5, 6.0 * d6);
107    let (q4, q5, q6) = (4.0 * b4, 4.0 * b5, 4.0 * b6);
108    let (r4, r5, r6) = (
109        3.0 * y23 * y23 / l23,
110        3.0 * y31 * y31 / l31,
111        3.0 * y12 * y12 / l12,
112    );
113
114    let (l1, l2) = (l[1], l[2]);
115    let hx_xi = [
116        p6 * (1.0 - 2.0 * l1) + (p5 - p6) * l2,
117        q6 * (1.0 - 2.0 * l1) - (q5 + q6) * l2,
118        -4.0 + 6.0 * (l1 + l2) + r6 * (1.0 - 2.0 * l1) - l2 * (r5 + r6),
119        -p6 * (1.0 - 2.0 * l1) + l2 * (p4 + p6),
120        q6 * (1.0 - 2.0 * l1) - l2 * (q6 - q4),
121        -2.0 + 6.0 * l1 + r6 * (1.0 - 2.0 * l1) + l2 * (r4 - r6),
122        -l2 * (p5 + p4),
123        l2 * (q4 - q5),
124        -l2 * (r5 - r4),
125    ];
126    let hx_eta = [
127        -p5 * (1.0 - 2.0 * l2) - l1 * (p6 - p5),
128        q5 * (1.0 - 2.0 * l2) - l1 * (q5 + q6),
129        -4.0 + 6.0 * (l1 + l2) + r5 * (1.0 - 2.0 * l2) - l1 * (r5 + r6),
130        l1 * (p4 + p6),
131        l1 * (q4 - q6),
132        -l1 * (r6 - r4),
133        p5 * (1.0 - 2.0 * l2) - l1 * (p4 + p5),
134        q5 * (1.0 - 2.0 * l2) + l1 * (q4 - q5),
135        -2.0 + 6.0 * l2 + r5 * (1.0 - 2.0 * l2) + l1 * (r4 - r5),
136    ];
137    let hy_xi = [
138        t6 * (1.0 - 2.0 * l1) + l2 * (t5 - t6),
139        1.0 + r6 * (1.0 - 2.0 * l1) - l2 * (r5 + r6),
140        -q6 * (1.0 - 2.0 * l1) + l2 * (q5 + q6),
141        -t6 * (1.0 - 2.0 * l1) + l2 * (t4 + t6),
142        -1.0 + r6 * (1.0 - 2.0 * l1) + l2 * (r4 - r6),
143        -q6 * (1.0 - 2.0 * l1) - l2 * (q4 - q6),
144        -l2 * (t4 + t5),
145        l2 * (r4 - r5),
146        -l2 * (q4 - q5),
147    ];
148    let hy_eta = [
149        -t5 * (1.0 - 2.0 * l2) - l1 * (t6 - t5),
150        1.0 + r5 * (1.0 - 2.0 * l2) - l1 * (r5 + r6),
151        -q5 * (1.0 - 2.0 * l2) + l1 * (q5 + q6),
152        l1 * (t4 + t6),
153        l1 * (r4 - r6),
154        -l1 * (q4 - q6),
155        t5 * (1.0 - 2.0 * l2) - l1 * (t4 + t5),
156        -1.0 + r5 * (1.0 - 2.0 * l2) + l1 * (r4 - r5),
157        -q5 * (1.0 - 2.0 * l2) - l1 * (q4 - q5),
158    ];
159    let s = 1.0 / (2.0 * area);
160    let mut out = [[0.0; 9]; 4];
161    for i in 0..9 {
162        out[0][i] = (hx_xi[i] * y31 + hx_eta[i] * y12) * s;
163        out[1][i] = (hx_xi[i] * -x31 + hx_eta[i] * -x12) * s;
164        out[2][i] = (hy_xi[i] * y31 + hy_eta[i] * y12) * s;
165        out[3][i] = (hy_xi[i] * -x31 + hy_eta[i] * -x12) * s;
166    }
167    out
168}
169
170fn point_geometry(axes: &[Vector3<f64>; 3], x: &[[f64; 3]; 2]) -> [PointGeometry; 4] {
171    let area = 0.5
172        * (x[0][0] * x[1][1] + x[0][1] * x[1][2] + x[0][2] * x[1][0]
173            - x[0][0] * x[1][2]
174            - x[0][1] * x[1][0]
175            - x[0][2] * x[1][1]);
176    let b: [f64; 3] = [x[1][1] - x[1][2], x[1][2] - x[1][0], x[1][0] - x[1][1]];
177    let c: [f64; 3] = [-x[0][1] + x[0][2], -x[0][2] + x[0][0], -x[0][0] + x[0][1]];
178    let xsj = (x[0][1] - x[0][0]) * (x[1][2] - x[1][0]) - (x[1][1] - x[1][0]) * (x[0][2] - x[0][0]);
179
180    std::array::from_fn(|g| {
181        let (l, w) = GAUSS[g];
182        let bend = shape_bend(l, x, area);
183        let mut local = Bmat::zeros();
184        for j in 0..3 {
185            let (nx, ny) = (b[j] / (2.0 * area), c[j] / (2.0 * area));
186            let d = drill_row(j, l, &b, &c, area);
187            let col = 6 * j;
188            // Membrane rows: [N,x 0 Nu,x; 0 N,y Nv,y; N,y N,x Nv,x + Nu,y] on (u, v, theta_z).
189            local[(0, col)] = nx;
190            local[(0, col + 5)] = d[0];
191            local[(1, col + 1)] = ny;
192            local[(1, col + 5)] = d[3];
193            local[(2, col)] = ny;
194            local[(2, col + 1)] = nx;
195            local[(2, col + 5)] = d[1] + d[2];
196            // Bending rows on (w, theta_x, theta_y) = local columns 2, 3, 4; OpenSees negates Bbend.
197            for k in 0..3 {
198                local[(3, col + 2 + k)] = -bend[0][3 * j + k];
199                local[(4, col + 2 + k)] = -bend[3][3 * j + k];
200                local[(5, col + 2 + k)] = -(bend[1][3 * j + k] + bend[2][3 * j + k]);
201            }
202        }
203        // Local -> global DOFs: each node's translations and rotations rotate with the axes.
204        let mut t = SMatrix::<f64, 18, 18>::zeros();
205        for j in 0..3 {
206            for block in 0..2 {
207                for (r, axis) in axes.iter().enumerate() {
208                    for k in 0..3 {
209                        t[(6 * j + 3 * block + r, 6 * j + 3 * block + k)] = axis[k];
210                    }
211                }
212            }
213        }
214        let b_global = local * t;
215        let mut b_force = b_global;
216        // Negate the bending-bending block of the *local* B before rotating: rows 3..6, local columns 2..5.
217        let mut local_force = local;
218        for j in 0..3 {
219            for r in 3..6 {
220                for k in 2..5 {
221                    local_force[(r, 6 * j + k)] *= -1.0;
222                }
223            }
224        }
225        b_force.copy_from(&(local_force * t));
226        PointGeometry {
227            b: b_global,
228            b_force,
229            dvol: 0.5 * w * xsj,
230            n: l,
231        }
232    })
233}
234
235/// `[shpDrill[0..4][i]]` of OpenSees at area coordinates `l`.
236fn drill_row(i: usize, l: [f64; 3], b: &[f64; 3], c: &[f64; 3], area: f64) -> [f64; 4] {
237    let s = 1.0 / (4.0 * area);
238    match i {
239        0 => {
240            let pb = b[2] * l[1] - b[1] * l[2];
241            let pc = c[2] * l[1] - c[1] * l[2];
242            [b[0] * pb * s, c[0] * pb * s, b[0] * pc * s, c[0] * pc * s]
243        }
244        1 => {
245            let pb = b[0] * l[2] - b[2] * l[0];
246            let pc = c[0] * l[2] - c[2] * l[0];
247            [b[1] * pb * s, c[1] * pb * s, b[1] * pc * s, c[1] * pc * s]
248        }
249        _ => {
250            let pb = b[1] * l[0] - b[0] * l[1];
251            let pc = c[1] * l[0] - c[0] * l[1];
252            [b[2] * pb * s, c[2] * pb * s, b[2] * pc * s, c[2] * pc * s]
253        }
254    }
255}
256
257fn dofs(nodes: &[Node3Id; 3]) -> [DofRef<Node3Id>; 18] {
258    std::array::from_fn(|a| (nodes[a / 6], (a % 6) as u8))
259}
260
261impl Shell3 {
262    pub fn new(nodes: [Node3Id; 3], section: ShellSection) -> Self {
263        Shell3 {
264            nodes,
265            sections: std::array::from_fn(|_| section.clone()),
266            cache: None,
267        }
268    }
269
270    /// The Gauss points' sections, in OpenSees' Gauss-point order.
271    pub fn sections(&self) -> &[ShellSection; 4] {
272        &self.sections
273    }
274
275    fn coords(&self, view: &View<'_>) -> [[f64; 3]; 3] {
276        std::array::from_fn(|i| view.get(self.nodes[i]).coords)
277    }
278
279    fn compute_cache(&self, view: &View<'_>) -> Result<Cache, &'static str> {
280        let (axes, xl) = frame(&self.coords(view))?;
281        Ok(Cache {
282            axes,
283            points: point_geometry(&axes, &xl),
284        })
285    }
286
287    fn cache_or_compute(&self, view: &View<'_>) -> Cache {
288        self.cache.unwrap_or_else(|| {
289            self.compute_cache(view)
290                .expect("unvalidated Shell3 geometry")
291        })
292    }
293
294    pub(super) fn dof_mask(&self) -> DofMask {
295        DofMask::all(SPATIAL_NDF)
296    }
297
298    pub(super) fn validate(&self, view: &View<'_>) -> Result<(), &'static str> {
299        self.sections[0].validate().map_err(|e| e.message())?;
300        let (_, x) = frame(&self.coords(view))?;
301        validate_triangle(&[[x[0][0], x[1][0]], [x[0][1], x[1][1]], [x[0][2], x[1][2]]])
302    }
303
304    pub(super) fn prepare(&mut self, view: &View<'_>) {
305        self.cache = self.compute_cache(view).ok();
306    }
307
308    fn displacement(&self, view: &View<'_>) -> Vec18 {
309        Vec18::from_fn(|a, _| view.get(self.nodes[a / 6]).displacement[a % 6])
310    }
311
312    pub(super) fn assemble_tangent<S: TangentSink<Node3Id>>(&self, view: &View<'_>, sink: &mut S) {
313        let cache = self.cache_or_compute(view);
314        let u = self.displacement(view);
315        let mut k = Mat18::zeros();
316        let mut r = Vec18::zeros();
317        for (g, section) in cache.points.iter().zip(&self.sections) {
318            let (stress, d) = section.trial_resultant(&(g.b * u));
319            k += g.dvol * g.b_force.transpose() * d * g.b;
320            r += g.dvol * g.b_force.transpose() * stress;
321        }
322        sink.add(&dofs(&self.nodes), &k, &r);
323    }
324
325    /// `integral(N_i dA)` per node (`area / 3`, exactly): the pressure and body-force share and the lumped mass.
326    fn nodal_weights(&self, cache: &Cache) -> [f64; 3] {
327        let mut w = [0.0; 3];
328        for g in &cache.points {
329            for (i, wi) in w.iter_mut().enumerate() {
330                *wi += g.n[i] * g.dvol;
331            }
332        }
333        w
334    }
335
336    pub(super) fn assemble_mass<S: VectorSink<Node3Id>>(&self, view: &View<'_>, sink: &mut S) {
337        let w = self.nodal_weights(&self.cache_or_compute(view));
338        let rho_h = self.sections[0].rho_h();
339        let m = Vec18::from_fn(|a, _| if a % 6 < 3 { rho_h * w[a / 6] } else { 0.0 });
340        sink.add(&dofs(&self.nodes), &m);
341    }
342
343    /// Equivalent nodal forces of `load` in global axes (body acceleration times `rho h`, pressure along `e3`).
344    pub(super) fn load_vector(&self, view: &View<'_>, load: &ElementLoad3) -> Vec18 {
345        let cache = self.cache_or_compute(view);
346        let w = self.nodal_weights(&cache);
347        let rho_h = self.sections[0].rho_h();
348        let normal = cache.axes[2];
349        Vec18::from_fn(|a, _| {
350            let (node, slot) = (a / 6, a % 6);
351            if slot < 3 {
352                w[node] * (rho_h * load.body[slot] + load.pressure * normal[slot])
353            } else {
354                0.0
355            }
356        })
357    }
358
359    pub(super) fn assemble_load<S: VectorSink<Node3Id>>(
360        &self,
361        view: &View<'_>,
362        load: &ElementLoad3,
363        sink: &mut S,
364    ) {
365        sink.add(&dofs(&self.nodes), &self.load_vector(view, load));
366    }
367
368    pub(super) fn commit(&mut self, view: &View<'_>) {
369        let cache = self.cache_or_compute(view);
370        let u = self.displacement(view);
371        for (section, g) in self.sections.iter_mut().zip(&cache.points) {
372            let strain: ShellVector = g.b * u;
373            *section = section.commit(&strain);
374        }
375    }
376
377    /// The element's resisting force in global axes (18).
378    pub(super) fn local_force(&self, view: &View<'_>) -> ElementForce {
379        let cache = self.cache_or_compute(view);
380        let u = self.displacement(view);
381        let mut r = Vec18::zeros();
382        for (g, section) in cache.points.iter().zip(&self.sections) {
383            let (stress, _) = section.trial_resultant(&(g.b * u));
384            r += g.dvol * g.b_force.transpose() * stress;
385        }
386        r.into()
387    }
388
389    /// Generalized strain and resultants at each Gauss point (local axes, section ordering).
390    pub(super) fn shell_responses(&self, view: &View<'_>) -> Vec<ShellResponse> {
391        let cache = self.cache_or_compute(view);
392        let u = self.displacement(view);
393        cache
394            .points
395            .iter()
396            .zip(&self.sections)
397            .map(|(g, section)| {
398                let strain = g.b * u;
399                ShellResponse {
400                    strain: strain.into(),
401                    resultant: section.trial_resultant(&strain).0.into(),
402                }
403            })
404            .collect()
405    }
406}
carapace_core/model/elements/shell3.rs

carapace_core/model/elements/
shell3.rs
