use crate::math::assert::Assert;
use crate::math::assert::perturbation;
use crate::math::{Current, Quantity};
use crate::{
EPSILON,
constitutive::cohesive::elastic::LinearElastic,
fem::block::element::{
ElementNodalCoordinates, ElementNodalReferenceCoordinates, FiniteElement,
cohesive::{
elastic::ElasticCohesiveElement,
linear::wedge::{N, P, Wedge},
},
solid::{ElementNodalForcesSolid, ElementNodalStiffnessesSolid},
},
math::{Rank2, Scalar, Tensor, TensorRank2, assert::AssertionError},
mechanics::test::get_rotation_reference_configuration,
units::{Area, Force, Length, Stress, StressPerLength},
};
const NORMAL_DISPLACEMENT: Scalar = 1.2;
const NORMAL_STIFFNESS: Quantity<StressPerLength> = StressPerLength::pascals_per_meter(3.4);
const TANGENTIAL_STIFFNESS: Quantity<StressPerLength> = StressPerLength::pascals_per_meter(5.6);
const TANGENTIAL_DISPLACEMENT: Scalar = 7.8;
const TANGENTIAL_TRACTION_P: Quantity<Stress> =
Quantity::new(TANGENTIAL_STIFFNESS.value() * TANGENTIAL_DISPLACEMENT / P as Scalar);
const NORMAL_TRACTION_P: Quantity<Stress> =
Quantity::new(NORMAL_STIFFNESS.value() * NORMAL_DISPLACEMENT / P as Scalar);
const COORDINATES: [[Scalar; 3]; N] = [
[-0.47979299, 0.48230032, 0.0],
[2.69165013, 0.37308724, 0.0],
[-0.47600989, 3.2116273, 0.0],
[-0.47979299, 0.48230032, 0.0],
[2.69165013, 0.37308724, 0.0],
[-0.47600989, 3.2116273, 0.0],
];
const MODEL: LinearElastic = LinearElastic {
normal_stiffness: NORMAL_STIFFNESS,
tangential_stiffness: TANGENTIAL_STIFFNESS,
};
#[test]
fn temporary_1() -> Result<(), AssertionError> {
let coordinates = ElementNodalReferenceCoordinates::from(COORDINATES);
let element = Wedge::from(coordinates.clone());
Assert::default().eq_within_tols(
&element.nodal_forces(&MODEL, &coordinates.into())?,
&[[0.0; 3]; N].into(),
)
}
#[test]
fn temporary_2() -> Result<(), AssertionError> {
let mut coordinates = ElementNodalReferenceCoordinates::from(COORDINATES);
let element = Wedge::from(coordinates.clone());
coordinates.iter_mut().skip(P).for_each(|coordinate| {
coordinate[0] += Quantity::new(TANGENTIAL_DISPLACEMENT);
coordinate[2] += Quantity::new(NORMAL_DISPLACEMENT)
});
let area = element
.integration_weights()
.into_iter()
.sum::<Quantity<Area>>();
let tangential_force = (TANGENTIAL_TRACTION_P * area).value_as::<Force>();
let normal_force = (NORMAL_TRACTION_P * area).value_as::<Force>();
Assert::default().eq_within_tols(
&element.nodal_forces(&MODEL, &coordinates.into())?,
&[
[-tangential_force, 0.0, -normal_force],
[-tangential_force, 0.0, -normal_force],
[-tangential_force, 0.0, -normal_force],
[tangential_force, 0.0, normal_force],
[tangential_force, 0.0, normal_force],
[tangential_force, 0.0, normal_force],
]
.into(),
)
}
#[test]
fn temporary_3() -> Result<(), AssertionError> {
let coordinates = ElementNodalReferenceCoordinates::from(COORDINATES)
.iter()
.map(|coordinate| get_rotation_reference_configuration() * coordinate)
.collect::<ElementNodalReferenceCoordinates<N>>();
let element = Wedge::from(coordinates.clone());
Assert::default().eq_within_tols(
&element.nodal_forces(&MODEL, &coordinates.into())?,
&[[0.0; 3]; N].into(),
)
}
#[test]
fn temporary_4() -> Result<(), AssertionError> {
let coordinates_0 = ElementNodalReferenceCoordinates::from(COORDINATES)
.iter()
.map(|coordinate| get_rotation_reference_configuration() * coordinate)
.collect::<ElementNodalReferenceCoordinates<N>>();
let element = Wedge::from(coordinates_0);
let mut coordinates = ElementNodalReferenceCoordinates::from(COORDINATES);
coordinates.iter_mut().skip(P).for_each(|coordinate| {
coordinate[0] += Quantity::new(TANGENTIAL_DISPLACEMENT);
coordinate[2] += Quantity::new(NORMAL_DISPLACEMENT);
});
coordinates = coordinates
.into_iter()
.map(|coordinate| get_rotation_reference_configuration() * coordinate)
.collect();
let area = element
.integration_weights()
.into_iter()
.sum::<Quantity<Area>>();
let tangential_force = (TANGENTIAL_TRACTION_P * area).value_as::<Force>();
let normal_force = (NORMAL_TRACTION_P * area).value_as::<Force>();
let nodal_forces_rotated_back = element
.nodal_forces(&MODEL, &coordinates.into())?
.into_iter()
.map(|nodal_force| {
TensorRank2::<3, Current, Current>::from(
get_rotation_reference_configuration().transpose(),
) * nodal_force
})
.collect::<ElementNodalForcesSolid<N>>();
Assert::default().eq_within_tols(
&nodal_forces_rotated_back,
&[
[-tangential_force, 0.0, -normal_force],
[-tangential_force, 0.0, -normal_force],
[-tangential_force, 0.0, -normal_force],
[tangential_force, 0.0, normal_force],
[tangential_force, 0.0, normal_force],
[tangential_force, 0.0, normal_force],
]
.into(),
)
}
#[test]
fn temporary_5() -> Result<(), AssertionError> {
let coordinates_0 = ElementNodalReferenceCoordinates::from(COORDINATES);
let coordinates = ElementNodalCoordinates::from(coordinates_0.clone());
let element = Wedge::from(coordinates_0);
let mut finite_difference = Quantity::<Force>::new(0.0);
let nodal_stiffnesses_fd = (0..N)
.map(|a| {
(0..N)
.map(|b| {
(0..3)
.map(|i| {
(0..3)
.map(|j| {
let mut nodal_coordinates = coordinates.clone();
nodal_coordinates[b][j] += perturbation(0.5 * EPSILON);
finite_difference =
element.nodal_forces(&MODEL, &nodal_coordinates)?[a][i];
nodal_coordinates[b][j] -= perturbation(EPSILON);
finite_difference -=
element.nodal_forces(&MODEL, &nodal_coordinates)?[a][i];
Ok(finite_difference / perturbation::<Length>(EPSILON))
})
.collect()
})
.collect()
})
.collect()
})
.collect::<Result<ElementNodalStiffnessesSolid<N>, AssertionError>>()?;
Assert::default().eq_within_fd_tol(
&element.nodal_stiffnesses(&MODEL, &coordinates)?,
&nodal_stiffnesses_fd,
)
}
#[test]
fn temporary_6() -> Result<(), AssertionError> {
let coordinates_0 = ElementNodalReferenceCoordinates::from([
[-0.57177033, -0.20395894, 0.23629102],
[1.49477913, 1.72253902, 1.40527015],
[-2.31789525, -0.2546453, 2.40281722],
[-0.57177033, -0.20395894, 0.23629102],
[1.49477913, 1.72253902, 1.40527015],
[-2.31789525, -0.2546453, 2.40281722],
]);
let element = Wedge::from(coordinates_0);
let coordinates = ElementNodalCoordinates::from([
[-0.64542355, -0.31521986, 0.2103109],
[1.50161765, 1.80846799, 1.49664724],
[-2.29750971, -0.08562506, 2.28063606],
[4.72044386, 3.95736046, 4.01544368],
[6.80745386, 6.13361434, 5.46225216],
[3.14323173, 3.98543986, 6.22717385],
]);
let mut finite_difference = Quantity::<Force>::new(0.0);
let nodal_stiffnesses_fd = (0..N)
.map(|a| {
(0..N)
.map(|b| {
(0..3)
.map(|i| {
(0..3)
.map(|j| {
let mut nodal_coordinates = coordinates.clone();
nodal_coordinates[b][j] += perturbation(0.5 * EPSILON);
finite_difference =
element.nodal_forces(&MODEL, &nodal_coordinates)?[a][i];
nodal_coordinates[b][j] -= perturbation(EPSILON);
finite_difference -=
element.nodal_forces(&MODEL, &nodal_coordinates)?[a][i];
Ok(finite_difference / perturbation::<Length>(EPSILON))
})
.collect()
})
.collect()
})
.collect()
})
.collect::<Result<ElementNodalStiffnessesSolid<N>, AssertionError>>()?;
Assert::default().eq_within_fd_tol(
&element.nodal_stiffnesses(&MODEL, &coordinates)?,
&nodal_stiffnesses_fd,
)
}