gemlab 2.0.0

Geometry and meshes laboratory for finite element analyses
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232

use super::CommonArgs;
use crate::StrError;
use russell_lab::{Matrix, Vector};

/// Implements the shape(N) times vector(V) dot gradient(B) integration case 09 (e.g., variable density matrix)
///
/// Callback function: `f(v, p, N, B)`
///
/// Variable density coefficients:
///
/// ```text
///       ⌠    →   →
/// Kᵐⁿ = │ Nᵐ v ⊗ Bⁿ α dΩ
/// ▔     ⌡
///       Ωₑ
/// ```
///
/// The numerical integration is:
///
/// ```text
///         nip-1    →   →  →   →   →       →
/// Kᵐⁿᵢⱼ ≈   Σ   Nᵐ(ιᵖ) vᵢ(ιᵖ) Bⁿⱼ(ιᵖ) |J|(ιᵖ) wᵖ α
///          p=0
/// ```
///
/// # Results
///
/// ```text
///     ┌                                               ┐
///     | K⁰⁰₀₀ K⁰⁰₀₁ K⁰¹₀₀ K⁰¹₀₁ K⁰²₀₀ K⁰²₀₁ ··· K⁰ⁿ₀ⱼ |  ⟸  ii0
///     | K⁰⁰₁₀ K⁰⁰₁₁ K⁰¹₁₀ K⁰¹₁₁ K⁰²₁₀ K⁰²₁₁ ··· K⁰ⁿ₁ⱼ |
///     | K¹⁰₀₀ K¹⁰₀₁ K¹¹₀₀ K¹¹₀₁ K¹²₀₀ K¹²₀₁ ··· K¹ⁿ₀ⱼ |
/// K = | K¹⁰₁₀ K¹⁰₁₁ K¹¹₁₀ K¹¹₁₁ K¹²₁₀ K¹²₁₁ ··· K¹ⁿ₁ⱼ |
///     | K²⁰₀₀ K²⁰₀₁ K²¹₀₀ K²¹₀₁ K²²₀₀ K²²₀₁ ··· K²ⁿ₀ⱼ |
///     | K²⁰₁₀ K²⁰₁₁ K²¹₁₀ K²¹₁₁ K²²₁₀ K²²₁₁ ··· K²ⁿ₁ⱼ |
///     |  ···   ···   ···   ···   ···   ···  ···  ···  |
///     | Kᵐ⁰ᵢ₀ Kᵐ⁰ᵢ₁ Kᵐ¹ᵢ₀ Kᵐ¹ᵢ₁ Kᵐ²ᵢ₀ Kᵐ²ᵢ₁ ··· Kᵐⁿᵢⱼ |  ⟸  ii := i + m ⋅ space_ndim
///     └                                               ┘
///        ⇑                                        ⇑
///       jj0                                       jj := j + n ⋅ space_ndim
///
/// m = ii / space_ndim    n = jj / space_ndim
/// i = ii % space_ndim    j = jj % space_ndim
/// ```
///
/// # Arguments
///
/// * `kk` -- A matrix containing all `Kᵐⁿᵢⱼ` values, one after another, and sequentially placed as shown
///   above (in 2D). `m` and `n` are the indices of the node and `i` and `j` correspond to `space_ndim`.
///   The dimensions must be `nrow(K) ≥ ii0 + nnode ⋅ space_ndim` and `ncol(K) ≥ jj0 + nnode ⋅ space_ndim`.
/// * `args` --- Common arguments
/// * `fn_v` -- Function `f(v,p,N,B)` that computes `v(x(ιᵖ))`, given `0 ≤ p ≤ ngauss`,
///   shape functions N(ιᵖ), and gradients B(ιᵖ). `v.dim() = space_ndim`.
pub fn mat_09_nvb<F>(kk: &mut Matrix, args: &mut CommonArgs, mut fn_v: F) -> Result<(), StrError>
where
    F: FnMut(&mut Vector, usize, &Vector, &Matrix) -> Result<(), StrError>,
{
    // check
    let (space_ndim, nnode) = args.pad.xxt.dims();
    let (nrow_kk, ncol_kk) = kk.dims();
    let (ii0, jj0) = (args.ii0, args.jj0);
    if nrow_kk < ii0 + nnode * space_ndim {
        return Err("nrow(K) must be ≥ ii0 + nnode ⋅ space_ndim");
    }
    if ncol_kk < jj0 + nnode * space_ndim {
        return Err("ncol(K) must be ≥ jj0 + nnode ⋅ space_ndim");
    }

    // allocate auxiliary vector
    let mut v = Vector::new(space_ndim);

    // clear output matrix
    if args.clear {
        kk.fill(0.0);
    }

    // loop over integration points
    for p in 0..args.gauss.npoint() {
        // ksi coordinates and weight
        let iota = args.gauss.coords(p);
        let weight = args.gauss.weight(p);

        // calculate interpolation functions, Jacobian and gradient
        (args.pad.fn_interp)(&mut args.pad.interp, iota); // N
        let det_jac = args.pad.calc_gradient(iota)?; // B

        // calculate v
        let nn = &args.pad.interp;
        let bb = &args.pad.gradient;
        fn_v(&mut v, p, nn, bb)?;

        // calculate coefficient
        let c = if args.axisymmetric {
            let mut r = 0.0; // radius @ x(ιᵖ)
            for m in 0..nnode {
                r += nn[m] * args.pad.xxt.get(0, m);
            }
            det_jac * weight * args.alpha * r
        } else {
            det_jac * weight * args.alpha
        };

        // add contribution to K matrix
        if space_ndim == 2 {
            for m in 0..nnode {
                for n in 0..nnode {
                    kk.add(ii0 + 0 + m * 2, jj0 + 0 + n * 2, c * nn[m] * v[0] * bb.get(n, 0));
                    kk.add(ii0 + 0 + m * 2, jj0 + 1 + n * 2, c * nn[m] * v[0] * bb.get(n, 1));

                    kk.add(ii0 + 1 + m * 2, jj0 + 0 + n * 2, c * nn[m] * v[1] * bb.get(n, 0));
                    kk.add(ii0 + 1 + m * 2, jj0 + 1 + n * 2, c * nn[m] * v[1] * bb.get(n, 1));
                }
            }
        } else {
            for m in 0..nnode {
                for n in 0..nnode {
                    kk.add(ii0 + 0 + m * 3, jj0 + 0 + n * 3, c * nn[m] * v[0] * bb.get(n, 0));
                    kk.add(ii0 + 0 + m * 3, jj0 + 1 + n * 3, c * nn[m] * v[0] * bb.get(n, 1));
                    kk.add(ii0 + 0 + m * 3, jj0 + 2 + n * 3, c * nn[m] * v[0] * bb.get(n, 2));

                    kk.add(ii0 + 1 + m * 3, jj0 + 0 + n * 3, c * nn[m] * v[1] * bb.get(n, 0));
                    kk.add(ii0 + 1 + m * 3, jj0 + 1 + n * 3, c * nn[m] * v[1] * bb.get(n, 1));
                    kk.add(ii0 + 1 + m * 3, jj0 + 2 + n * 3, c * nn[m] * v[1] * bb.get(n, 2));

                    kk.add(ii0 + 2 + m * 3, jj0 + 0 + n * 3, c * nn[m] * v[2] * bb.get(n, 0));
                    kk.add(ii0 + 2 + m * 3, jj0 + 1 + n * 3, c * nn[m] * v[2] * bb.get(n, 1));
                    kk.add(ii0 + 2 + m * 3, jj0 + 2 + n * 3, c * nn[m] * v[2] * bb.get(n, 2));
                }
            }
        }
    }
    Ok(())
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

#[cfg(test)]
mod tests {
    use crate::integ::testing::aux;
    use crate::integ::{self, AnalyticalTet4, AnalyticalTri3, CommonArgs, Gauss};
    use russell_lab::{mat_approx_eq, Matrix, Vector};

    #[test]
    fn capture_some_errors() {
        let mut pad = aux::gen_pad_lin2(1.0);
        let mut kk = Matrix::new(4, 4);
        let mut vv = Vector::new(0);
        let nn = Vector::new(0);
        let bb = Matrix::new(0, 0);
        let f = |_: &mut Vector, _: usize, _: &Vector, _: &Matrix| Ok(());
        f(&mut vv, 0, &nn, &bb).unwrap();
        let gauss = Gauss::new(pad.kind);
        let mut args = CommonArgs::new(&mut pad, &gauss);
        args.ii0 = 1;
        assert_eq!(
            integ::mat_09_nvb(&mut kk, &mut args, f).err(),
            Some("nrow(K) must be ≥ ii0 + nnode ⋅ space_ndim")
        );
        args.ii0 = 0;
        args.jj0 = 1;
        assert_eq!(
            integ::mat_09_nvb(&mut kk, &mut args, f).err(),
            Some("ncol(K) must be ≥ jj0 + nnode ⋅ space_ndim")
        );
        args.jj0 = 0;
        // more errors
        assert_eq!(
            integ::mat_09_nvb(&mut kk, &mut args, f).err(),
            Some("calc_gradient requires that geo_ndim = space_ndim")
        );
        let mut pad = aux::gen_pad_tri3();
        let mut kk = Matrix::new(6, 6);
        let mut args = CommonArgs::new(&mut pad, &gauss);
        assert_eq!(
            integ::mat_09_nvb(&mut kk, &mut args, |_, _, _, _| Err("stop")).err(),
            Some("stop")
        );
    }

    #[test]
    fn mat_09_nvb_tri3_works() {
        let mut pad = aux::gen_pad_tri3();
        let mut kk = Matrix::new(3 * 2, 3 * 2);
        let ana = AnalyticalTri3::new(&pad);
        // constant
        let (v0, v1) = (2.0, 3.0);
        let kk_correct = ana.mat_09_nvb(v0, v1);
        // println!("{}", kk_correct);
        let class = pad.kind.class();
        let tolerances = [1e-15];
        let selection: Vec<_> = [3].iter().map(|n| Gauss::new_sized(class, *n).unwrap()).collect();
        selection.iter().zip(tolerances).for_each(|(ips, tol)| {
            // println!("nip={}, tol={:.e}", ips.len(), tol);
            let mut args = CommonArgs::new(&mut pad, ips);
            integ::mat_09_nvb(&mut kk, &mut args, |v, _, _, _| {
                v[0] = v0;
                v[1] = v1;
                Ok(())
            })
            .unwrap();
            // println!("{}", kk);
            mat_approx_eq(&kk, &kk_correct, tol);
        });
    }

    #[test]
    fn mat_09_nvb_tet4_works() {
        let mut pad = aux::gen_pad_tet4();
        let mut kk = Matrix::new(4 * 3, 4 * 3);
        let ana = AnalyticalTet4::new(&pad);
        // constant
        let (v0, v1, v2) = (2.0, 3.0, 4.0);
        let kk_correct = ana.mat_09_nvb(v0, v1, v2);
        // println!("{}", kk_correct);
        let class = pad.kind.class();
        let tolerances = [1e-15];
        let selection: Vec<_> = [4].iter().map(|n| Gauss::new_sized(class, *n).unwrap()).collect();
        selection.iter().zip(tolerances).for_each(|(ips, tol)| {
            // println!("nip={}, tol={:.e}", ips.len(), tol);
            let mut args = CommonArgs::new(&mut pad, ips);
            integ::mat_09_nvb(&mut kk, &mut args, |v, _, _, _| {
                v[0] = v0;
                v[1] = v1;
                v[2] = v2;
                Ok(())
            })
            .unwrap();
            // println!("{}", kk);
            mat_approx_eq(&kk, &kk_correct, tol);
        });
    }
}