<commit_before>use super::Scratchpad;
use crate::StrError;
use russell_lab::mat_mat_mul;

impl Scratchpad {
    /// Calculates the gradient of the interpolation functions
    ///
    /// **Note:** This function works with `geo_ndim == space_ndim` only.
    ///
    /// The gradient is given by:
    ///
    /// ```text
    ///             →
    /// →  →    dNᵐ(ξ)
    /// Gᵐ(ξ) = ——————
    ///            →
    ///           dx
    /// ```
    ///
    /// which can be organized in an (nnode,space_ndim) matrix `G` as follows
    ///
    /// ```text
    /// G = L · J⁻¹
    /// ```
    ///
    /// # Output
    ///
    /// * `deriv` -- interpolation functions (nnode)
    /// * `jacobian` -- Jacobian matrix (space_ndim,geo_ndim)
    /// * `inv_jacobian` -- inverse Jacobian matrix (space_ndim,space_ndim)
    /// * `gradient` -- gradient matrix (nnode,space_ndim)
    /// * Returns the determinant of the Jacobian matrix
    ///
    /// # Input
    ///
    /// * `ksi` -- reference coordinates ξ with len ≥ geo_ndim
    ///
    /// # Example
    ///
    /// ```
    /// use gemlab::shapes::{GeoKind, Scratchpad};
    /// use gemlab::StrError;
    /// use russell_chk::assert_vec_approx_eq;
    /// use russell_lab::Matrix;
    ///
    /// fn main() -> Result<(), StrError> {
    ///     //  3-------------2         ξ₀   ξ₁
    ///     //  |      ξ₁     |  node    r    s
    ///     //  |      |      |     0 -1.0 -1.0
    ///     //  |      +--ξ₀  |     1  1.0 -1.0
    ///     //  |             |     2  1.0  1.0
    ///     //  |             |     3 -1.0  1.0
    ///     //  0-------------1
    ///
    ///     let a = 3.0;
    ///     let space_ndim = 2;
    ///     let mut pad = Scratchpad::new(space_ndim, GeoKind::Qua4)?;
    ///     pad.set_xx(0, 0, 0.0);
    ///     pad.set_xx(0, 1, 0.0);
    ///     pad.set_xx(1, 0, 2.0 * a);
    ///     pad.set_xx(1, 1, 0.0);
    ///     pad.set_xx(2, 0, 2.0 * a);
    ///     pad.set_xx(2, 1, a);
    ///     pad.set_xx(3, 0, 0.0);
    ///     pad.set_xx(3, 1, a);
    ///
    ///     pad.calc_gradient(&[0.0, 0.0])?;
    ///
    ///     let correct_gg = Matrix::from(&[
    ///         [-1.0 / (4.0 * a), -1.0 / (2.0 * a)],
    ///         [1.0 / (4.0 * a), -1.0 / (2.0 * a)],
    ///         [1.0 / (4.0 * a), 1.0 / (2.0 * a)],
    ///         [-1.0 / (4.0 * a), 1.0 / (2.0 * a)],
    ///     ]);
    ///     assert_vec_approx_eq!(pad.gradient.as_data(), correct_gg.as_data(), 1e-15);
    ///     Ok(())
    /// }
    /// ```
    pub fn calc_gradient(&mut self, ksi: &[f64]) -> Result<f64, StrError> {
        // check
        let (space_ndim, geo_ndim) = self.jacobian.dims();
        if geo_ndim != space_ndim {
            return Err("calc_gradient requires that geo_ndim = space_ndim");
        }

        // Jacobian matrix J: dx/dξ
        let det_jac = self.calc_jacobian(ksi)?;

        // gradient: G = L · J⁻¹
        mat_mat_mul(&mut self.gradient, 1.0, &self.deriv, &self.inv_jacobian).unwrap(); // cannot fail because the dims are checked
        Ok(det_jac)
    }
}

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

#[cfg(test)]
mod tests {
    use crate::shapes::scratchpad_testing::aux;
    use crate::shapes::{GeoKind, Scratchpad};
    use crate::util::ONE_BY_3;
    use russell_chk::assert_deriv_approx_eq;
    use russell_lab::{copy_vector, Vector};

    #[test]
    fn calc_gradient_handles_errors() {
        let mut pad = Scratchpad::new(2, GeoKind::Lin2).unwrap();
        assert_eq!(
            pad.calc_gradient(&[0.0, 0.0]).err(),
            Some("calc_gradient requires that geo_ndim = space_ndim")
        );

        let mut pad = Scratchpad::new(2, GeoKind::Tri3).unwrap();
        assert_eq!(
            pad.calc_gradient(&[0.0, 0.0]).err(),
            Some("all components of the coordinates matrix must be set first")
        );
    }

    // Holds arguments for numerical differentiation of N with respect to x => G (gradient) matrix
    struct ArgsNumGrad {
        pad: Scratchpad, // scratchpad to send to calc_coords
        at_x: Vector,    // at x coord value
        x: Vector,       // temporary x coord
        ksi: Vec<f64>,   // temporary reference coord
        m: usize,        // node index from 0 to nnode
        j: usize,        // dimension index from 0 to space_ndim
    }

    // Computes Nᵐ(ξ(x)) with variable v := xⱼ
    fn nn_given_x(v: f64, args: &mut ArgsNumGrad) -> f64 {
        copy_vector(&mut args.x, &args.at_x).unwrap();
        args.x[args.j] = v;
        args.pad.approximate_ksi(&mut args.ksi, &args.x, 10, 1e-14).unwrap();
        (args.pad.fn_interp)(&mut args.pad.interp, &args.ksi);
        args.pad.interp[args.m]
    }

    #[test]
    fn calc_gradient_works() {
        // kind (except Lin) and tolerances
        let problem = vec![
            // Tri
            (GeoKind::Tri3, 1e-12),
            (GeoKind::Tri6, 1e-10),
            (GeoKind::Tri10, 1e-9),
            (GeoKind::Tri15, 1e-9),
            // Qua
            (GeoKind::Qua4, 1e-11),
            (GeoKind::Qua8, 1e-10),
            (GeoKind::Qua9, 1e-10),
            (GeoKind::Qua12, 1e-10),
            (GeoKind::Qua16, 1e-10),
            (GeoKind::Qua17, 1e-9),
            // Tet
            (GeoKind::Tet4, 1e-12),
            (GeoKind::Tet10, 1e-9),
            (GeoKind::Tet20, 1e-9),
            // Hex
            (GeoKind::Hex8, 1e-11),
            (GeoKind::Hex20, 1e-10),
            (GeoKind::Hex32, 1e-9),
        ];

        // loop over shapes
        for (kind, tol) in problem {
            // scratchpad with coordinates
            let geo_ndim = kind.ndim();
            let space_ndim = usize::max(2, geo_ndim);
            let mut pad = aux::gen_scratchpad_with_coords(space_ndim, kind);
<commit_msg>Increase tol of gradient test for notebook<commit_after>use super::Scratchpad;
use crate::StrError;
use russell_lab::mat_mat_mul;

impl Scratchpad {
    /// Calculates the gradient of the interpolation functions
    ///
    /// **Note:** This function works with `geo_ndim == space_ndim` only.
    ///
    /// The gradient is given by:
    ///
    /// ```text
    ///             →
    /// →  →    dNᵐ(ξ)
    /// Gᵐ(ξ) = ——————
    ///            →
    ///           dx
    /// ```
    ///
    /// which can be organized in an (nnode,space_ndim) matrix `G` as follows
    ///
    /// ```text
    /// G = L · J⁻¹
    /// ```
    ///
    /// # Output
    ///
    /// * `deriv` -- interpolation functions (nnode)
    /// * `jacobian` -- Jacobian matrix (space_ndim,geo_ndim)
    /// * `inv_jacobian` -- inverse Jacobian matrix (space_ndim,space_ndim)
    /// * `gradient` -- gradient matrix (nnode,space_ndim)
    /// * Returns the determinant of the Jacobian matrix
    ///
    /// # Input
    ///
    /// * `ksi` -- reference coordinates ξ with len ≥ geo_ndim
    ///
    /// # Example
    ///
    /// ```
    /// use gemlab::shapes::{GeoKind, Scratchpad};
    /// use gemlab::StrError;
    /// use russell_chk::assert_vec_approx_eq;
    /// use russell_lab::Matrix;
    ///
    /// fn main() -> Result<(), StrError> {
    ///     //  3-------------2         ξ₀   ξ₁
    ///     //  |      ξ₁     |  node    r    s
    ///     //  |      |      |     0 -1.0 -1.0
    ///     //  |      +--ξ₀  |     1  1.0 -1.0
    ///     //  |             |     2  1.0  1.0
    ///     //  |             |     3 -1.0  1.0
    ///     //  0-------------1
    ///
    ///     let a = 3.0;
    ///     let space_ndim = 2;
    ///     let mut pad = Scratchpad::new(space_ndim, GeoKind::Qua4)?;
    ///     pad.set_xx(0, 0, 0.0);
    ///     pad.set_xx(0, 1, 0.0);
    ///     pad.set_xx(1, 0, 2.0 * a);
    ///     pad.set_xx(1, 1, 0.0);
    ///     pad.set_xx(2, 0, 2.0 * a);
    ///     pad.set_xx(2, 1, a);
    ///     pad.set_xx(3, 0, 0.0);
    ///     pad.set_xx(3, 1, a);
    ///
    ///     pad.calc_gradient(&[0.0, 0.0])?;
    ///
    ///     let correct_gg = Matrix::from(&[
    ///         [-1.0 / (4.0 * a), -1.0 / (2.0 * a)],
    ///         [1.0 / (4.0 * a), -1.0 / (2.0 * a)],
    ///         [1.0 / (4.0 * a), 1.0 / (2.0 * a)],
    ///         [-1.0 / (4.0 * a), 1.0 / (2.0 * a)],
    ///     ]);
    ///     assert_vec_approx_eq!(pad.gradient.as_data(), correct_gg.as_data(), 1e-15);
    ///     Ok(())
    /// }
    /// ```
    pub fn calc_gradient(&mut self, ksi: &[f64]) -> Result<f64, StrError> {
        // check
        let (space_ndim, geo_ndim) = self.jacobian.dims();
        if geo_ndim != space_ndim {
            return Err("calc_gradient requires that geo_ndim = space_ndim");
        }

        // Jacobian matrix J: dx/dξ
        let det_jac = self.calc_jacobian(ksi)?;

        // gradient: G = L · J⁻¹
        mat_mat_mul(&mut self.gradient, 1.0, &self.deriv, &self.inv_jacobian).unwrap(); // cannot fail because the dims are checked
        Ok(det_jac)
    }
}

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

#[cfg(test)]
mod tests {
    use crate::shapes::scratchpad_testing::aux;
    use crate::shapes::{GeoKind, Scratchpad};
    use crate::util::ONE_BY_3;
    use russell_chk::assert_deriv_approx_eq;
    use russell_lab::{copy_vector, Vector};

    #[test]
    fn calc_gradient_handles_errors() {
        let mut pad = Scratchpad::new(2, GeoKind::Lin2).unwrap();
        assert_eq!(
            pad.calc_gradient(&[0.0, 0.0]).err(),
            Some("calc_gradient requires that geo_ndim = space_ndim")
        );

        let mut pad = Scratchpad::new(2, GeoKind::Tri3).unwrap();
        assert_eq!(
            pad.calc_gradient(&[0.0, 0.0]).err(),
            Some("all components of the coordinates matrix must be set first")
        );
    }

    // Holds arguments for numerical differentiation of N with respect to x => G (gradient) matrix
    struct ArgsNumGrad {
        pad: Scratchpad, // scratchpad to send to calc_coords
        at_x: Vector,    // at x coord value
        x: Vector,       // temporary x coord
        ksi: Vec<f64>,   // temporary reference coord
        m: usize,        // node index from 0 to nnode
        j: usize,        // dimension index from 0 to space_ndim
    }

    // Computes Nᵐ(ξ(x)) with variable v := xⱼ
    fn nn_given_x(v: f64, args: &mut ArgsNumGrad) -> f64 {
        copy_vector(&mut args.x, &args.at_x).unwrap();
        args.x[args.j] = v;
        args.pad.approximate_ksi(&mut args.ksi, &args.x, 10, 1e-14).unwrap();
        (args.pad.fn_interp)(&mut args.pad.interp, &args.ksi);
        args.pad.interp[args.m]
    }

    #[test]
    fn calc_gradient_works() {
        // kind (except Lin) and tolerances
        let problem = vec![
            // Tri
            (GeoKind::Tri3, 1e-12),
            (GeoKind::Tri6, 1e-10),
            (GeoKind::Tri10, 1e-9),
            (GeoKind::Tri15, 1e-9),
            // Qua
            (GeoKind::Qua4, 1e-11),
            (GeoKind::Qua8, 1e-10),
            (GeoKind::Qua9, 1e-10),
            (GeoKind::Qua12, 1e-10),
            (GeoKind::Qua16, 1e-9),
            (GeoKind::Qua17, 1e-9),
            // Tet
            (GeoKind::Tet4, 1e-12),
            (GeoKind::Tet10, 1e-9),
            (GeoKind::Tet20, 1e-9),
            // Hex
            (GeoKind::Hex8, 1e-11),
            (GeoKind::Hex20, 1e-10),
            (GeoKind::Hex32, 1e-9),
        ];

        // loop over shapes
        for (kind, tol) in problem {
            // scratchpad with coordinates
            let geo_ndim = kind.ndim();
            let space_ndim = usize::max(2, geo_ndim);
            let mut pad = aux::gen_scratchpad_with_coords(space_ndim, kind);

{'fraction_non_alphanumeric': 0.2000343997248022, 'fraction_numerical': 0.033453732370141036, 'mean_word_length': 1.3864149394623435, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 17, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4425935', 'n_tokens_mistral': 4371, 'n_tokens_neox': 4176, 'n_words': 888}