carapace_core/model/elements/
shell4.rs

1//! `Shell4`: the 4-node MITC4 flat shell (OpenSees `ShellMITC4`), small-displacement, 2x2 Gauss
2//! integration, six DOFs per node (24 total).
3//!
4//! The formulation follows OpenSees so results agree to solver precision:
5//!
6//! * Local frame: `e1` from the edge midpoints (`0.5 (x2 + x1 - x3 - x0)`), `e2` its
7//!   Gram-Schmidt completion with `0.5 (x3 + x2 - x1 - x0)`, `e3 = e1 x e2`. The shell is
8//!   treated as flat in that frame (nodes are projected on `e1`, `e2`); warped quads are
9//!   tolerated the same way.
10//! * Membrane: bilinear displacement field. Bending: Reissner-Mindlin rotations, with the
11//!   transverse shear strains interpolated from the four edge midpoints (MITC4) to avoid
12//!   shear locking. Drilling: a penalty on the difference between the in-plane rotation
13//!   and the skew part of the displacement gradient, stiffness the smallest eigenvalue of the
14//!   section's membrane tangent.
15//! * Mass: OpenSees uses a consistent translational mass with no rotational inertia; here
16//!   `assemble_mass` is a row-sum lumped translational mass (same total, same centroid).
17//!
18//! Generalized strains and resultants use the section's ordering and sign convention
19//! (`ShellSection`): `[eps_x, eps_y, gamma_xy, kappa_x, kappa_y, 2 kappa_xy, gamma_xz,
20//! gamma_yz]`. The bending rows of the `B` matrix used for the internal force are negated
21//! because the section's bending tangent is negative (the OpenSees convention).
22
23use nalgebra::{Matrix2, SMatrix, SVector, Vector3};
24
25use super::{DofMask, DofRef, ElementForce, NodeView, ShellResponse, TangentSink, VectorSink};
26use crate::model::continuum::validate_quad;
27use crate::model::{ElementLoad3, Node3Id, ShellSection, ShellVector, SPATIAL_NDF, SPATIAL_NDIM};
28
29type View<'a> = NodeView<'a, SPATIAL_NDIM, SPATIAL_NDF, Node3Id>;
30type Vec24 = SVector<f64, 24>;
31type Mat24 = SMatrix<f64, 24, 24>;
32type Bmat = SMatrix<f64, 8, 24>;
33
34const ONE_OVER_ROOT3: f64 = 0.577_350_269_189_625_7;
35/// Gauss points in OpenSees' order: `(-,-), (+,-), (+,+), (-,+)`.
36const SG: [f64; 4] = [
37    -ONE_OVER_ROOT3,
38    ONE_OVER_ROOT3,
39    ONE_OVER_ROOT3,
40    -ONE_OVER_ROOT3,
41];
42const TG: [f64; 4] = [
43    -ONE_OVER_ROOT3,
44    -ONE_OVER_ROOT3,
45    ONE_OVER_ROOT3,
46    ONE_OVER_ROOT3,
47];
48
49#[derive(Debug, Clone)]
50pub struct Shell4 {
51    pub nodes: [Node3Id; 4],
52    sections: [ShellSection; 4],
53    cache: Option<Cache>,
54}
55
56/// Geometry-only quantities, computed once by `prepare`.
57#[derive(Debug, Clone, Copy, PartialEq)]
58struct Cache {
59    /// Local axes `e1`, `e2`, `e3` in global components.
60    axes: [Vector3<f64>; 3],
61    points: [PointGeometry; 4],
62    /// Drilling stiffness penalty.
63    ktt: f64,
64}
65
66#[derive(Debug, Clone, Copy, PartialEq)]
67struct PointGeometry {
68    /// Strain-displacement matrix of the generalized strains.
69    b: Bmat,
70    /// The same with the bending-bending block negated, for the internal force.
71    b_force: Bmat,
72    /// Drilling strain-displacement row.
73    drill: Vec24,
74    /// `det J * w`.
75    dvol: f64,
76    /// Shape-function values.
77    n: [f64; 4],
78}
79
80type Frame = ([Vector3<f64>; 3], [[f64; 2]; 4]);
81
82/// The local frame and the nodes' coordinates in it.
83fn frame(coords: &[[f64; 3]; 4]) -> Result<Frame, &'static str> {
84    let x: [Vector3<f64>; 4] = std::array::from_fn(|i| Vector3::from(coords[i]));
85    let mut v1 = 0.5 * (x[2] + x[1] - x[3] - x[0]);
86    let mut v2 = 0.5 * (x[3] + x[2] - x[1] - x[0]);
87    let n1 = v1.norm();
88    if !n1.is_finite() || n1 == 0.0 {
89        return Err("degenerate shell geometry");
90    }
91    v1 /= n1;
92    v2 -= v2.dot(&v1) * v1;
93    let n2 = v2.norm();
94    if !n2.is_finite() || n2 == 0.0 {
95        return Err("degenerate shell geometry");
96    }
97    v2 /= n2;
98    let v3 = v1.cross(&v2);
99    Ok((
100        [v1, v2, v3],
101        std::array::from_fn(|i| [x[i].dot(&v1), x[i].dot(&v2)]),
102    ))
103}
104
105/// OpenSees' `shape2d`: `(dN/dx, dN/dy, N)` per node at `(s, t)` and the Jacobian determinant.
106fn shape2d(s: f64, t: f64, xl: &[[f64; 2]; 4]) -> ([[f64; 4]; 3], f64) {
107    const SS: [f64; 4] = [-0.5, 0.5, 0.5, -0.5];
108    const TT: [f64; 4] = [-0.5, -0.5, 0.5, 0.5];
109    let mut shp = [[0.0; 4]; 3];
110    for (i, (&si, &ti)) in SS.iter().zip(&TT).enumerate() {
111        shp[2][i] = (0.5 + si * s) * (0.5 + ti * t);
112        shp[0][i] = si * (0.5 + ti * t);
113        shp[1][i] = ti * (0.5 + si * s);
114    }
115    let mut xs = [[0.0; 2]; 2];
116    for (a, row) in xs.iter_mut().enumerate() {
117        for (b, v) in row.iter_mut().enumerate() {
118            *v = (0..4).map(|k| xl[k][a] * shp[b][k]).sum();
119        }
120    }
121    let det = xs[0][0] * xs[1][1] - xs[0][1] * xs[1][0];
122    let inv = 1.0 / det;
123    let sx = [
124        [xs[1][1] * inv, -xs[0][1] * inv],
125        [-xs[1][0] * inv, xs[0][0] * inv],
126    ];
127    #[allow(clippy::needless_range_loop)] // reads and writes two rows of `shp` at column `i`
128    for i in 0..4 {
129        let dx = shp[0][i] * sx[0][0] + shp[1][i] * sx[1][0];
130        shp[1][i] = shp[0][i] * sx[0][1] + shp[1][i] * sx[1][1];
131        shp[0][i] = dx;
132    }
133    (shp, det)
134}
135
136fn point_geometry(axes: &[Vector3<f64>; 3], xl: &[[f64; 2]; 4]) -> [PointGeometry; 4] {
137    let [g1, g2, g3] = axes;
138    let (dx34, dy34) = (xl[2][0] - xl[3][0], xl[2][1] - xl[3][1]);
139    let (dx21, dy21) = (xl[1][0] - xl[0][0], xl[1][1] - xl[0][1]);
140    let (dx32, dy32) = (xl[2][0] - xl[1][0], xl[2][1] - xl[1][1]);
141    let (dx41, dy41) = (xl[3][0] - xl[0][0], xl[3][1] - xl[0][1]);
142    let q = 0.25;
143    let mut g = SMatrix::<f64, 4, 12>::zeros();
144    for (row, entries) in [
145        [
146            (0, -0.5),
147            (1, -dy41 * q),
148            (2, dx41 * q),
149            (9, 0.5),
150            (10, -dy41 * q),
151            (11, dx41 * q),
152        ],
153        [
154            (0, -0.5),
155            (1, -dy21 * q),
156            (2, dx21 * q),
157            (3, 0.5),
158            (4, -dy21 * q),
159            (5, dx21 * q),
160        ],
161        [
162            (3, -0.5),
163            (4, -dy32 * q),
164            (5, dx32 * q),
165            (6, 0.5),
166            (7, -dy32 * q),
167            (8, dx32 * q),
168        ],
169        [
170            (6, 0.5),
171            (7, -dy34 * q),
172            (8, dx34 * q),
173            (9, -0.5),
174            (10, -dy34 * q),
175            (11, dx34 * q),
176        ],
177    ]
178    .into_iter()
179    .enumerate()
180    {
181        for (col, v) in entries {
182            g[(row, col)] = v;
183        }
184    }
185
186    let ax = -xl[0][0] + xl[1][0] + xl[2][0] - xl[3][0];
187    let bx = xl[0][0] - xl[1][0] + xl[2][0] - xl[3][0];
188    let cx = -xl[0][0] - xl[1][0] + xl[2][0] + xl[3][0];
189    let ay = -xl[0][1] + xl[1][1] + xl[2][1] - xl[3][1];
190    let by = xl[0][1] - xl[1][1] + xl[2][1] - xl[3][1];
191    let cy = -xl[0][1] - xl[1][1] + xl[2][1] + xl[3][1];
192    let alpha = (ay / ax).atan();
193    let beta = std::f64::consts::FRAC_PI_2 - (cx / cy).atan();
194    let rot = Matrix2::new(beta.sin(), -alpha.sin(), -beta.cos(), alpha.cos());
195
196    std::array::from_fn(|i| {
197        let (s, t) = (SG[i], TG[i]);
198        let r1 = (cx + s * bx).hypot(cy + s * by);
199        let r2 = (ax + t * bx).hypot(ay + t * by);
200        let (shp, det) = shape2d(s, t, xl);
201
202        let mut ms = SMatrix::<f64, 2, 4>::zeros();
203        ms[(1, 0)] = 1.0 - s;
204        ms[(0, 1)] = 1.0 - t;
205        ms[(1, 2)] = 1.0 + s;
206        ms[(0, 3)] = 1.0 + t;
207        let mut bsv = ms * g;
208        for j in 0..12 {
209            bsv[(0, j)] *= r1 / (8.0 * det);
210            bsv[(1, j)] *= r2 / (8.0 * det);
211        }
212        let bs = rot * bsv;
213
214        let mut b = Bmat::zeros();
215        let mut drill = Vec24::zeros();
216        for j in 0..4 {
217            let (nx, ny, n) = (shp[0][j], shp[1][j], shp[2][j]);
218            let c = 6 * j;
219            // Membrane: [Nx 0; 0 Ny; Ny Nx] times [g1; g2] on the translations.
220            // Bending: [0 -Nx; Ny 0; Nx -Ny] times [g1; g2] on the rotations.
221            for k in 0..3 {
222                b[(0, c + k)] = nx * g1[k];
223                b[(1, c + k)] = ny * g2[k];
224                b[(2, c + k)] = ny * g1[k] + nx * g2[k];
225                b[(3, c + 3 + k)] = -nx * g2[k];
226                b[(4, c + 3 + k)] = ny * g1[k];
227                b[(5, c + 3 + k)] = nx * g1[k] - ny * g2[k];
228                // Shear: [w, theta1, theta2] -> [g3.u, g1.theta, g2.theta].
229                for r in 0..2 {
230                    b[(6 + r, c + k)] = bs[(r, 3 * j)] * g3[k];
231                    b[(6 + r, c + 3 + k)] = bs[(r, 3 * j + 1)] * g1[k] + bs[(r, 3 * j + 2)] * g2[k];
232                }
233                let (b1, b2, b6) = (-0.5 * ny, 0.5 * nx, -n);
234                drill[c + k] = b1 * g1[k] + b2 * g2[k];
235                drill[c + 3 + k] = b6 * g3[k];
236            }
237        }
238        let mut b_force = b;
239        for j in 0..4 {
240            for r in 3..6 {
241                for k in 0..3 {
242                    b_force[(r, 6 * j + 3 + k)] *= -1.0;
243                }
244            }
245        }
246        PointGeometry {
247            b,
248            b_force,
249            drill,
250            dvol: det,
251            n: shp[2],
252        }
253    })
254}
255
256fn dofs(nodes: &[Node3Id; 4]) -> [DofRef<Node3Id>; 24] {
257    std::array::from_fn(|a| (nodes[a / 6], (a % 6) as u8))
258}
259
260impl Shell4 {
261    pub fn new(nodes: [Node3Id; 4], section: ShellSection) -> Self {
262        Shell4 {
263            nodes,
264            sections: std::array::from_fn(|_| section.clone()),
265            cache: None,
266        }
267    }
268
269    /// The Gauss points' sections, in OpenSees' Gauss-point order.
270    pub fn sections(&self) -> &[ShellSection; 4] {
271        &self.sections
272    }
273
274    fn coords(&self, view: &View<'_>) -> [[f64; 3]; 4] {
275        std::array::from_fn(|i| view.get(self.nodes[i]).coords)
276    }
277
278    fn compute_cache(&self, view: &View<'_>) -> Result<Cache, &'static str> {
279        let (axes, xl) = frame(&self.coords(view))?;
280        Ok(Cache {
281            axes,
282            points: point_geometry(&axes, &xl),
283            ktt: self.sections[0].drilling_stiffness(),
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 Shell4 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 (_, xl) = frame(&self.coords(view))?;
301        validate_quad(&xl)
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<'_>) -> Vec24 {
309        Vec24::from_fn(|a, _| view.get(self.nodes[a / 6]).displacement[a % 6])
310    }
311
312    /// Generalized strain at each Gauss point.
313    fn strains(&self, cache: &Cache, u: &Vec24) -> [ShellVector; 4] {
314        std::array::from_fn(|g| cache.points[g].b * u)
315    }
316
317    pub(super) fn assemble_tangent<S: TangentSink<Node3Id>>(&self, view: &View<'_>, sink: &mut S) {
318        let cache = self.cache_or_compute(view);
319        let u = self.displacement(view);
320        let mut k = Mat24::zeros();
321        let mut r = Vec24::zeros();
322        for (g, section) in cache.points.iter().zip(&self.sections) {
323            let (stress, d) = section.trial_resultant(&(g.b * u));
324            k += g.dvol * g.b_force.transpose() * d * g.b;
325            k += cache.ktt * g.dvol * g.drill * g.drill.transpose();
326            r += g.dvol * g.b_force.transpose() * stress;
327            r += cache.ktt * g.dvol * g.drill.dot(&u) * g.drill;
328        }
329        sink.add(&dofs(&self.nodes), &k, &r);
330    }
331
332    /// `integral(N_i dA)` per node: the pressure and body-force share, and (times `rho h`)
333    /// the row-sum lumped mass.
334    fn nodal_weights(&self, cache: &Cache) -> [f64; 4] {
335        let mut w = [0.0; 4];
336        for g in &cache.points {
337            for (i, wi) in w.iter_mut().enumerate() {
338                *wi += g.n[i] * g.dvol;
339            }
340        }
341        w
342    }
343
344    pub(super) fn assemble_mass<S: VectorSink<Node3Id>>(&self, view: &View<'_>, sink: &mut S) {
345        let w = self.nodal_weights(&self.cache_or_compute(view));
346        let rho_h = self.sections[0].rho_h();
347        let m = Vec24::from_fn(|a, _| if a % 6 < 3 { rho_h * w[a / 6] } else { 0.0 });
348        sink.add(&dofs(&self.nodes), &m);
349    }
350
351    /// Equivalent nodal forces of `load` in global axes: the body force (a body
352    /// acceleration scaled by `rho h`) and the pressure along `e3`.
353    pub(super) fn load_vector(&self, view: &View<'_>, load: &ElementLoad3) -> Vec24 {
354        let cache = self.cache_or_compute(view);
355        let w = self.nodal_weights(&cache);
356        let rho_h = self.sections[0].rho_h();
357        let normal = cache.axes[2];
358        Vec24::from_fn(|a, _| {
359            let (node, slot) = (a / 6, a % 6);
360            if slot < 3 {
361                w[node] * (rho_h * load.body[slot] + load.pressure * normal[slot])
362            } else {
363                0.0
364            }
365        })
366    }
367
368    pub(super) fn assemble_load<S: VectorSink<Node3Id>>(
369        &self,
370        view: &View<'_>,
371        load: &ElementLoad3,
372        sink: &mut S,
373    ) {
374        sink.add(&dofs(&self.nodes), &self.load_vector(view, load));
375    }
376
377    pub(super) fn commit(&mut self, view: &View<'_>) {
378        let u = self.displacement(view);
379        let strains = self.strains(&self.cache_or_compute(view), &u);
380        for (section, strain) in self.sections.iter_mut().zip(&strains) {
381            *section = section.commit(strain);
382        }
383    }
384
385    /// The element's resisting force in global axes (24).
386    pub(super) fn local_force(&self, view: &View<'_>) -> ElementForce {
387        let cache = self.cache_or_compute(view);
388        let u = self.displacement(view);
389        let mut r = Vec24::zeros();
390        for (g, section) in cache.points.iter().zip(&self.sections) {
391            let (stress, _) = section.trial_resultant(&(g.b * u));
392            r += g.dvol * g.b_force.transpose() * stress;
393            r += cache.ktt * g.dvol * g.drill.dot(&u) * g.drill;
394        }
395        r.into()
396    }
397
398    /// Generalized strain and resultants at each Gauss point (local axes, section ordering).
399    pub(super) fn shell_responses(&self, view: &View<'_>) -> Vec<ShellResponse> {
400        let cache = self.cache_or_compute(view);
401        let u = self.displacement(view);
402        cache
403            .points
404            .iter()
405            .zip(&self.sections)
406            .map(|(g, section)| {
407                let strain = g.b * u;
408                ShellResponse {
409                    strain: strain.into(),
410                    resultant: section.trial_resultant(&strain).0.into(),
411                }
412            })
413            .collect()
414    }
415}
carapace_core/model/elements/shell4.rs

carapace_core/model/elements/
shell4.rs
