1#[cfg(test)]
2mod test;
3
4pub mod element;
5
6use self::element::{
7 ElasticFiniteElement, ElasticHyperviscousFiniteElement, FiniteElement, FiniteElementError,
8 FiniteElementMethods, HyperelasticFiniteElement, HyperviscoelasticFiniteElement,
9 SurfaceFiniteElement, ViscoelasticFiniteElement,
10};
11use super::*;
12use crate::{
13 defeat_message,
14 math::{
15 Banded, TestError,
16 integrate::{Explicit, IntegrationError},
17 optimize::{
18 EqualityConstraint, FirstOrderOptimization, FirstOrderRootFinding, OptimizationError,
19 SecondOrderOptimization, ZerothOrderRootFinding,
20 },
21 },
22 mechanics::Times,
23};
24use std::{
25 array::from_fn,
26 fmt::{self, Debug, Display, Formatter},
27 iter::repeat_n,
28};
29
30pub struct ElementBlock<F, const N: usize> {
31 connectivity: Connectivity<N>,
32 coordinates: ReferenceNodalCoordinatesBlock,
33 elements: Vec<F>,
34}
35
36impl<F, const N: usize> Debug for ElementBlock<F, N> {
37 fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
38 let element = match N {
39 3 => "LinearTriangle",
40 4 => "LinearTetrahedron",
41 8 => "LinearHexahedron",
42 10 => "CompositeTetrahedron",
43 _ => panic!(),
44 };
45 write!(
46 f,
47 "ElementBlock {{ elements: [{element}; {}] }}",
48 self.connectivity.len()
49 )
50 }
51}
52
53pub trait FiniteElementBlockMethods<C, F, const G: usize, const N: usize>
54where
55 F: FiniteElementMethods<C, G, N>,
56{
57 fn connectivity(&self) -> &Connectivity<N>;
58 fn coordinates(&self) -> &ReferenceNodalCoordinatesBlock;
59 fn deformation_gradients(
60 &self,
61 nodal_coordinates: &NodalCoordinatesBlock,
62 ) -> Vec<DeformationGradientList<G>>;
63 fn elements(&self) -> &[F];
64 fn nodal_coordinates_element(
65 &self,
66 element_connectivity: &[usize; N],
67 nodal_coordinates: &NodalCoordinatesBlock,
68 ) -> NodalCoordinates<N>;
69}
70
71pub trait FiniteElementBlock<C, F, const G: usize, const N: usize, Y>
72where
73 C: Constitutive<Y>,
74 F: FiniteElement<C, G, N, Y>,
75 Y: Parameters,
76{
77 fn new(
78 constitutive_model_parameters: Y,
79 connectivity: Connectivity<N>,
80 reference_nodal_coordinates: ReferenceNodalCoordinatesBlock,
81 ) -> Self;
82 fn reset(&mut self);
83}
84
85pub trait SurfaceFiniteElementBlock<C, F, const G: usize, const N: usize, const P: usize, Y>
86where
87 C: Constitutive<Y>,
88 F: SurfaceFiniteElement<C, G, N, P, Y>,
89 Y: Parameters,
90{
91 fn new(
92 constitutive_model_parameters: Y,
93 connectivity: Connectivity<N>,
94 reference_nodal_coordinates: ReferenceNodalCoordinatesBlock,
95 thickness: Scalar,
96 ) -> Self;
97}
98
99pub enum FiniteElementBlockError {
100 Upstream(String, String),
101}
102
103impl From<FiniteElementBlockError> for String {
104 fn from(error: FiniteElementBlockError) -> Self {
105 match error {
106 FiniteElementBlockError::Upstream(error, block) => {
107 format!(
108 "{error}\x1b[0;91m\n\
109 In finite element block: {block}."
110 )
111 }
112 }
113 }
114}
115
116impl From<FiniteElementBlockError> for TestError {
117 fn from(error: FiniteElementBlockError) -> Self {
118 Self {
119 message: error.to_string(),
120 }
121 }
122}
123
124impl Debug for FiniteElementBlockError {
125 fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
126 let error = match self {
127 Self::Upstream(error, block) => {
128 format!(
129 "{error}\x1b[0;91m\n\
130 In block: {block}."
131 )
132 }
133 };
134 write!(f, "\n{error}\n\x1b[0;2;31m{}\x1b[0m\n", defeat_message())
135 }
136}
137
138impl Display for FiniteElementBlockError {
139 fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
140 let error = match self {
141 Self::Upstream(error, block) => {
142 format!(
143 "{error}\x1b[0;91m\n\
144 In block: {block}."
145 )
146 }
147 };
148 write!(f, "{error}\x1b[0m")
149 }
150}
151
152impl<C, F, const G: usize, const N: usize> FiniteElementBlockMethods<C, F, G, N>
153 for ElementBlock<F, N>
154where
155 F: FiniteElementMethods<C, G, N>,
156{
157 fn connectivity(&self) -> &Connectivity<N> {
158 &self.connectivity
159 }
160 fn coordinates(&self) -> &ReferenceNodalCoordinatesBlock {
161 &self.coordinates
162 }
163 fn deformation_gradients(
164 &self,
165 nodal_coordinates: &NodalCoordinatesBlock,
166 ) -> Vec<DeformationGradientList<G>> {
167 self.elements()
168 .iter()
169 .zip(self.connectivity().iter())
170 .map(|(element, element_connectivity)| {
171 element.deformation_gradients(
172 &self.nodal_coordinates_element(element_connectivity, nodal_coordinates),
173 )
174 })
175 .collect()
176 }
177 fn elements(&self) -> &[F] {
178 &self.elements
179 }
180 fn nodal_coordinates_element(
181 &self,
182 element_connectivity: &[usize; N],
183 nodal_coordinates: &NodalCoordinatesBlock,
184 ) -> NodalCoordinates<N> {
185 element_connectivity
186 .iter()
187 .map(|node| nodal_coordinates[*node].clone())
188 .collect()
189 }
190}
191
192impl<C, F, const G: usize, const N: usize, Y> FiniteElementBlock<C, F, G, N, Y>
193 for ElementBlock<F, N>
194where
195 C: Constitutive<Y>,
196 F: FiniteElement<C, G, N, Y>,
197 Y: Parameters,
198{
199 fn new(
200 constitutive_model_parameters: Y,
201 connectivity: Connectivity<N>,
202 coordinates: ReferenceNodalCoordinatesBlock,
203 ) -> Self {
204 let elements = connectivity
205 .iter()
206 .map(|element_connectivity| {
207 <F>::new(
208 constitutive_model_parameters,
209 element_connectivity
210 .iter()
211 .map(|&node| coordinates[node].clone())
212 .collect(),
213 )
214 })
215 .collect();
216 Self {
217 connectivity,
218 coordinates,
219 elements,
220 }
221 }
222 fn reset(&mut self) {
223 self.elements.iter_mut().for_each(|element| element.reset())
224 }
225}
226
227impl<C, F, const G: usize, const N: usize, const P: usize, Y>
228 SurfaceFiniteElementBlock<C, F, G, N, P, Y> for ElementBlock<F, N>
229where
230 C: Constitutive<Y>,
231 F: SurfaceFiniteElement<C, G, N, P, Y>,
232 Y: Parameters,
233{
234 fn new(
235 constitutive_model_parameters: Y,
236 connectivity: Connectivity<N>,
237 coordinates: ReferenceNodalCoordinatesBlock,
238 thickness: Scalar,
239 ) -> Self {
240 let elements = connectivity
241 .iter()
242 .map(|element_connectivity| {
243 <F>::new(
244 constitutive_model_parameters,
245 element_connectivity
246 .iter()
247 .map(|node| coordinates[*node].clone())
248 .collect(),
249 &thickness,
250 )
251 })
252 .collect();
253 Self {
254 connectivity,
255 coordinates,
256 elements,
257 }
258 }
259}
260
261pub trait ElasticFiniteElementBlock<C, F, const G: usize, const N: usize>
262where
263 C: Elastic,
264 F: ElasticFiniteElement<C, G, N>,
265{
266 fn nodal_forces(
267 &self,
268 nodal_coordinates: &NodalCoordinatesBlock,
269 ) -> Result<NodalForcesBlock, FiniteElementBlockError>;
270 fn nodal_stiffnesses(
271 &self,
272 nodal_coordinates: &NodalCoordinatesBlock,
273 ) -> Result<NodalStiffnessesBlock, FiniteElementBlockError>;
274}
275
276pub trait ZerothOrderRoot<C, F, const G: usize, const N: usize>
277where
278 C: Elastic,
279 F: ElasticFiniteElement<C, G, N>,
280{
281 fn root(
282 &self,
283 equality_constraint: EqualityConstraint,
284 solver: impl ZerothOrderRootFinding<NodalCoordinatesBlock>,
285 ) -> Result<NodalCoordinatesBlock, OptimizationError>;
286}
287
288pub trait FirstOrderRoot<C, F, const G: usize, const N: usize>
289where
290 C: Elastic,
291 F: ElasticFiniteElement<C, G, N>,
292{
293 fn root(
294 &self,
295 equality_constraint: EqualityConstraint,
296 solver: impl FirstOrderRootFinding<
297 NodalForcesBlock,
298 NodalStiffnessesBlock,
299 NodalCoordinatesBlock,
300 >,
301 ) -> Result<NodalCoordinatesBlock, OptimizationError>;
302}
303
304pub trait HyperelasticFiniteElementBlock<C, F, const G: usize, const N: usize>
305where
306 C: Hyperelastic,
307 F: HyperelasticFiniteElement<C, G, N>,
308 Self: ElasticFiniteElementBlock<C, F, G, N>,
309{
310 fn helmholtz_free_energy(
311 &self,
312 nodal_coordinates: &NodalCoordinatesBlock,
313 ) -> Result<Scalar, FiniteElementBlockError>;
314}
315
316pub trait FirstOrderMinimize<C, F, const G: usize, const N: usize>
317where
318 C: Hyperelastic,
319 F: HyperelasticFiniteElement<C, G, N>,
320 Self: ElasticFiniteElementBlock<C, F, G, N>,
321{
322 fn minimize(
323 &self,
324 equality_constraint: EqualityConstraint,
325 solver: impl FirstOrderOptimization<Scalar, NodalForcesBlock>,
326 ) -> Result<NodalCoordinatesBlock, OptimizationError>;
327}
328
329pub trait SecondOrderMinimize<C, F, const G: usize, const N: usize>
330where
331 C: Hyperelastic,
332 F: HyperelasticFiniteElement<C, G, N>,
333 Self: ElasticFiniteElementBlock<C, F, G, N>,
334{
335 fn minimize(
336 &self,
337 equality_constraint: EqualityConstraint,
338 solver: impl SecondOrderOptimization<
339 Scalar,
340 NodalForcesBlock,
341 NodalStiffnessesBlock,
342 NodalCoordinatesBlock,
343 >,
344 ) -> Result<NodalCoordinatesBlock, OptimizationError>;
345}
346
347pub trait ViscoelasticFiniteElementBlock<C, F, const G: usize, const N: usize>
348where
349 C: Viscoelastic,
350 F: ViscoelasticFiniteElement<C, G, N>,
351{
352 fn deformation_gradient_rates(
353 &self,
354 nodal_coordinates: &NodalCoordinatesBlock,
355 nodal_velocities: &NodalVelocitiesBlock,
356 ) -> Vec<DeformationGradientRateList<G>>;
357 fn nodal_forces(
358 &self,
359 nodal_coordinates: &NodalCoordinatesBlock,
360 nodal_velocities: &NodalVelocitiesBlock,
361 ) -> Result<NodalForcesBlock, FiniteElementBlockError>;
362 fn nodal_stiffnesses(
363 &self,
364 nodal_coordinates: &NodalCoordinatesBlock,
365 nodal_velocities: &NodalVelocitiesBlock,
366 ) -> Result<NodalStiffnessesBlock, FiniteElementBlockError>;
367 fn nodal_velocities_element(
368 &self,
369 element_connectivity: &[usize; N],
370 nodal_velocities: &NodalVelocitiesBlock,
371 ) -> NodalVelocities<N>;
372 fn root(
373 &self,
374 equality_constraint: EqualityConstraint,
375 integrator: impl Explicit<NodalVelocitiesBlock, NodalVelocitiesHistory>,
376 time: &[Scalar],
377 solver: impl FirstOrderRootFinding<
378 NodalForcesBlock,
379 NodalStiffnessesBlock,
380 NodalCoordinatesBlock,
381 >,
382 ) -> Result<(Times, NodalCoordinatesHistory, NodalVelocitiesHistory), IntegrationError>;
383 #[doc(hidden)]
384 fn root_inner(
385 &self,
386 equality_constraint: EqualityConstraint,
387 nodal_coordinates: &NodalCoordinatesBlock,
388 solver: &impl FirstOrderRootFinding<
389 NodalForcesBlock,
390 NodalStiffnessesBlock,
391 NodalCoordinatesBlock,
392 >,
393 initial_guess: &NodalVelocitiesBlock,
394 ) -> Result<NodalVelocitiesBlock, OptimizationError>;
395}
396
397pub trait ElasticHyperviscousFiniteElementBlock<C, F, const G: usize, const N: usize>
398where
399 C: ElasticHyperviscous,
400 F: ElasticHyperviscousFiniteElement<C, G, N>,
401 Self: ViscoelasticFiniteElementBlock<C, F, G, N>,
402{
403 fn viscous_dissipation(
404 &self,
405 nodal_coordinates: &NodalCoordinatesBlock,
406 nodal_velocities: &NodalVelocitiesBlock,
407 ) -> Result<Scalar, FiniteElementBlockError>;
408 fn dissipation_potential(
409 &self,
410 nodal_coordinates: &NodalCoordinatesBlock,
411 nodal_velocities: &NodalVelocitiesBlock,
412 ) -> Result<Scalar, FiniteElementBlockError>;
413 fn minimize(
414 &self,
415 equality_constraint: EqualityConstraint,
416 integrator: impl Explicit<NodalVelocitiesBlock, NodalVelocitiesHistory>,
417 time: &[Scalar],
418 solver: impl SecondOrderOptimization<
419 Scalar,
420 NodalForcesBlock,
421 NodalStiffnessesBlock,
422 NodalCoordinatesBlock,
423 >,
424 ) -> Result<(Times, NodalCoordinatesHistory, NodalVelocitiesHistory), IntegrationError>;
425 #[doc(hidden)]
426 fn minimize_inner(
427 &self,
428 equality_constraint: EqualityConstraint,
429 nodal_coordinates: &NodalCoordinatesBlock,
430 solver: &impl SecondOrderOptimization<
431 Scalar,
432 NodalForcesBlock,
433 NodalStiffnessesBlock,
434 NodalCoordinatesBlock,
435 >,
436 initial_guess: &NodalVelocitiesBlock,
437 ) -> Result<NodalVelocitiesBlock, OptimizationError>;
438}
439
440pub trait HyperviscoelasticFiniteElementBlock<C, F, const G: usize, const N: usize>
441where
442 C: Hyperviscoelastic,
443 F: HyperviscoelasticFiniteElement<C, G, N>,
444 Self: ElasticHyperviscousFiniteElementBlock<C, F, G, N>,
445{
446 fn helmholtz_free_energy(
447 &self,
448 nodal_coordinates: &NodalCoordinatesBlock,
449 ) -> Result<Scalar, FiniteElementBlockError>;
450}
451
452impl<C, F, const G: usize, const N: usize> ElasticFiniteElementBlock<C, F, G, N>
453 for ElementBlock<F, N>
454where
455 C: Elastic,
456 F: ElasticFiniteElement<C, G, N>,
457 Self: FiniteElementBlockMethods<C, F, G, N>,
458{
459 fn nodal_forces(
460 &self,
461 nodal_coordinates: &NodalCoordinatesBlock,
462 ) -> Result<NodalForcesBlock, FiniteElementBlockError> {
463 let mut nodal_forces = NodalForcesBlock::zero(nodal_coordinates.len());
464 match self
465 .elements()
466 .iter()
467 .zip(self.connectivity().iter())
468 .try_for_each(|(element, element_connectivity)| {
469 element
470 .nodal_forces(
471 &self.nodal_coordinates_element(element_connectivity, nodal_coordinates),
472 )?
473 .iter()
474 .zip(element_connectivity.iter())
475 .for_each(|(nodal_force, &node)| nodal_forces[node] += nodal_force);
476 Ok::<(), FiniteElementError>(())
477 }) {
478 Ok(()) => Ok(nodal_forces),
479 Err(error) => Err(FiniteElementBlockError::Upstream(
480 format!("{error}"),
481 format!("{self:?}"),
482 )),
483 }
484 }
485 fn nodal_stiffnesses(
486 &self,
487 nodal_coordinates: &NodalCoordinatesBlock,
488 ) -> Result<NodalStiffnessesBlock, FiniteElementBlockError> {
489 let mut nodal_stiffnesses = NodalStiffnessesBlock::zero(nodal_coordinates.len());
490 match self
491 .elements()
492 .iter()
493 .zip(self.connectivity().iter())
494 .try_for_each(|(element, element_connectivity)| {
495 element
496 .nodal_stiffnesses(
497 &self.nodal_coordinates_element(element_connectivity, nodal_coordinates),
498 )?
499 .iter()
500 .zip(element_connectivity.iter())
501 .for_each(|(object, &node_a)| {
502 object.iter().zip(element_connectivity.iter()).for_each(
503 |(nodal_stiffness, &node_b)| {
504 nodal_stiffnesses[node_a][node_b] += nodal_stiffness
505 },
506 )
507 });
508 Ok::<(), FiniteElementError>(())
509 }) {
510 Ok(()) => Ok(nodal_stiffnesses),
511 Err(error) => Err(FiniteElementBlockError::Upstream(
512 format!("{error}"),
513 format!("{self:?}"),
514 )),
515 }
516 }
517}
518
519impl<C, F, const G: usize, const N: usize> ZerothOrderRoot<C, F, G, N> for ElementBlock<F, N>
520where
521 C: Elastic,
522 F: ElasticFiniteElement<C, G, N>,
523 Self: FiniteElementBlockMethods<C, F, G, N>,
524{
525 fn root(
526 &self,
527 equality_constraint: EqualityConstraint,
528 solver: impl ZerothOrderRootFinding<NodalCoordinatesBlock>,
529 ) -> Result<NodalCoordinatesBlock, OptimizationError> {
530 solver.root(
531 |nodal_coordinates: &NodalCoordinatesBlock| Ok(self.nodal_forces(nodal_coordinates)?),
532 self.coordinates().clone().into(),
533 equality_constraint,
534 )
535 }
536}
537
538impl<C, F, const G: usize, const N: usize> FirstOrderRoot<C, F, G, N> for ElementBlock<F, N>
539where
540 C: Elastic,
541 F: ElasticFiniteElement<C, G, N>,
542 Self: FiniteElementBlockMethods<C, F, G, N>,
543{
544 fn root(
545 &self,
546 equality_constraint: EqualityConstraint,
547 solver: impl FirstOrderRootFinding<
548 NodalForcesBlock,
549 NodalStiffnessesBlock,
550 NodalCoordinatesBlock,
551 >,
552 ) -> Result<NodalCoordinatesBlock, OptimizationError> {
553 solver.root(
554 |nodal_coordinates: &NodalCoordinatesBlock| Ok(self.nodal_forces(nodal_coordinates)?),
555 |nodal_coordinates: &NodalCoordinatesBlock| {
556 Ok(self.nodal_stiffnesses(nodal_coordinates)?)
557 },
558 self.coordinates().clone().into(),
559 equality_constraint,
560 )
561 }
562}
563
564impl<C, F, const G: usize, const N: usize> HyperelasticFiniteElementBlock<C, F, G, N>
565 for ElementBlock<F, N>
566where
567 C: Hyperelastic,
568 F: HyperelasticFiniteElement<C, G, N>,
569 Self: ElasticFiniteElementBlock<C, F, G, N>,
570{
571 fn helmholtz_free_energy(
572 &self,
573 nodal_coordinates: &NodalCoordinatesBlock,
574 ) -> Result<Scalar, FiniteElementBlockError> {
575 match self
576 .elements()
577 .iter()
578 .zip(self.connectivity().iter())
579 .map(|(element, element_connectivity)| {
580 element.helmholtz_free_energy(
581 &self.nodal_coordinates_element(element_connectivity, nodal_coordinates),
582 )
583 })
584 .sum()
585 {
586 Ok(helmholtz_free_energy) => Ok(helmholtz_free_energy),
587 Err(error) => Err(FiniteElementBlockError::Upstream(
588 format!("{error}"),
589 format!("{self:?}"),
590 )),
591 }
592 }
593}
594
595impl<C, F, const G: usize, const N: usize> FirstOrderMinimize<C, F, G, N> for ElementBlock<F, N>
596where
597 C: Hyperelastic,
598 F: HyperelasticFiniteElement<C, G, N>,
599 Self: ElasticFiniteElementBlock<C, F, G, N>,
600{
601 fn minimize(
602 &self,
603 equality_constraint: EqualityConstraint,
604 solver: impl FirstOrderOptimization<Scalar, NodalForcesBlock>,
605 ) -> Result<NodalCoordinatesBlock, OptimizationError> {
606 solver.minimize(
607 |nodal_coordinates: &NodalCoordinatesBlock| {
608 Ok(self.helmholtz_free_energy(nodal_coordinates)?)
609 },
610 |nodal_coordinates: &NodalCoordinatesBlock| Ok(self.nodal_forces(nodal_coordinates)?),
611 self.coordinates().clone().into(),
612 equality_constraint,
613 )
614 }
615}
616
617impl<C, F, const G: usize, const N: usize> SecondOrderMinimize<C, F, G, N> for ElementBlock<F, N>
618where
619 C: Hyperelastic,
620 F: HyperelasticFiniteElement<C, G, N>,
621 Self: ElasticFiniteElementBlock<C, F, G, N>,
622{
623 fn minimize(
624 &self,
625 equality_constraint: EqualityConstraint,
626 solver: impl SecondOrderOptimization<
627 Scalar,
628 NodalForcesBlock,
629 NodalStiffnessesBlock,
630 NodalCoordinatesBlock,
631 >,
632 ) -> Result<NodalCoordinatesBlock, OptimizationError> {
633 let banded = band(
634 self.connectivity(),
635 &equality_constraint,
636 self.coordinates().len(),
637 );
638 solver.minimize(
639 |nodal_coordinates: &NodalCoordinatesBlock| {
640 Ok(self.helmholtz_free_energy(nodal_coordinates)?)
641 },
642 |nodal_coordinates: &NodalCoordinatesBlock| Ok(self.nodal_forces(nodal_coordinates)?),
643 |nodal_coordinates: &NodalCoordinatesBlock| {
644 Ok(self.nodal_stiffnesses(nodal_coordinates)?)
645 },
646 self.coordinates().clone().into(),
647 equality_constraint,
648 Some(banded),
649 )
650 }
651}
652
653impl<C, F, const G: usize, const N: usize> ViscoelasticFiniteElementBlock<C, F, G, N>
654 for ElementBlock<F, N>
655where
656 C: Viscoelastic,
657 F: ViscoelasticFiniteElement<C, G, N>,
658 Self: FiniteElementBlockMethods<C, F, G, N>,
659{
660 fn deformation_gradient_rates(
661 &self,
662 nodal_coordinates: &NodalCoordinatesBlock,
663 nodal_velocities: &NodalVelocitiesBlock,
664 ) -> Vec<DeformationGradientRateList<G>> {
665 self.elements()
666 .iter()
667 .zip(self.connectivity().iter())
668 .map(|(element, element_connectivity)| {
669 element.deformation_gradient_rates(
670 &self.nodal_coordinates_element(element_connectivity, nodal_coordinates),
671 &self.nodal_velocities_element(element_connectivity, nodal_velocities),
672 )
673 })
674 .collect()
675 }
676 fn nodal_forces(
677 &self,
678 nodal_coordinates: &NodalCoordinatesBlock,
679 nodal_velocities: &NodalVelocitiesBlock,
680 ) -> Result<NodalForcesBlock, FiniteElementBlockError> {
681 let mut nodal_forces = NodalForcesBlock::zero(nodal_coordinates.len());
682 match self
683 .elements()
684 .iter()
685 .zip(self.connectivity().iter())
686 .try_for_each(|(element, element_connectivity)| {
687 element
688 .nodal_forces(
689 &self.nodal_coordinates_element(element_connectivity, nodal_coordinates),
690 &self.nodal_velocities_element(element_connectivity, nodal_velocities),
691 )?
692 .iter()
693 .zip(element_connectivity.iter())
694 .for_each(|(nodal_force, &node)| nodal_forces[node] += nodal_force);
695 Ok::<(), FiniteElementError>(())
696 }) {
697 Ok(()) => Ok(nodal_forces),
698 Err(error) => Err(FiniteElementBlockError::Upstream(
699 format!("{error}"),
700 format!("{self:?}"),
701 )),
702 }
703 }
704 fn nodal_stiffnesses(
705 &self,
706 nodal_coordinates: &NodalCoordinatesBlock,
707 nodal_velocities: &NodalVelocitiesBlock,
708 ) -> Result<NodalStiffnessesBlock, FiniteElementBlockError> {
709 let mut nodal_stiffnesses = NodalStiffnessesBlock::zero(nodal_coordinates.len());
710 match self
711 .elements()
712 .iter()
713 .zip(self.connectivity().iter())
714 .try_for_each(|(element, element_connectivity)| {
715 element
716 .nodal_stiffnesses(
717 &self.nodal_coordinates_element(element_connectivity, nodal_coordinates),
718 &self.nodal_velocities_element(element_connectivity, nodal_velocities),
719 )?
720 .iter()
721 .zip(element_connectivity.iter())
722 .for_each(|(object, &node_a)| {
723 object.iter().zip(element_connectivity.iter()).for_each(
724 |(nodal_stiffness, &node_b)| {
725 nodal_stiffnesses[node_a][node_b] += nodal_stiffness
726 },
727 )
728 });
729 Ok::<(), FiniteElementError>(())
730 }) {
731 Ok(()) => Ok(nodal_stiffnesses),
732 Err(error) => Err(FiniteElementBlockError::Upstream(
733 format!("{error}"),
734 format!("{self:?}"),
735 )),
736 }
737 }
738 fn nodal_velocities_element(
739 &self,
740 element_connectivity: &[usize; N],
741 nodal_velocities: &NodalVelocitiesBlock,
742 ) -> NodalVelocities<N> {
743 element_connectivity
744 .iter()
745 .map(|node| nodal_velocities[*node].clone())
746 .collect()
747 }
748 fn root(
749 &self,
750 equality_constraint: EqualityConstraint,
751 integrator: impl Explicit<NodalVelocitiesBlock, NodalVelocitiesHistory>,
752 time: &[Scalar],
753 solver: impl FirstOrderRootFinding<
754 NodalForcesBlock,
755 NodalStiffnessesBlock,
756 NodalCoordinatesBlock,
757 >,
758 ) -> Result<(Times, NodalCoordinatesHistory, NodalVelocitiesHistory), IntegrationError> {
759 let mut solution = NodalVelocitiesBlock::zero(self.coordinates().len());
760 integrator.integrate(
761 |_: Scalar, nodal_coordinates: &NodalCoordinatesBlock| {
762 solution = self.root_inner(
763 equality_constraint.clone(),
764 nodal_coordinates,
765 &solver,
766 &solution,
767 )?;
768 Ok(solution.clone())
769 },
770 time,
771 self.coordinates().clone().into(),
772 )
773 }
774 #[doc(hidden)]
775 fn root_inner(
776 &self,
777 equality_constraint: EqualityConstraint,
778 nodal_coordinates: &NodalCoordinatesBlock,
779 solver: &impl FirstOrderRootFinding<
780 NodalForcesBlock,
781 NodalStiffnessesBlock,
782 NodalCoordinatesBlock,
783 >,
784 initial_guess: &NodalVelocitiesBlock,
785 ) -> Result<NodalVelocitiesBlock, OptimizationError> {
786 solver.root(
787 |nodal_velocities: &NodalVelocitiesBlock| {
788 Ok(self.nodal_forces(nodal_coordinates, nodal_velocities)?)
789 },
790 |nodal_velocities: &NodalVelocitiesBlock| {
791 Ok(self.nodal_stiffnesses(nodal_coordinates, nodal_velocities)?)
792 },
793 initial_guess.clone(),
794 equality_constraint,
795 )
796 }
797}
798
799impl<C, F, const G: usize, const N: usize> ElasticHyperviscousFiniteElementBlock<C, F, G, N>
800 for ElementBlock<F, N>
801where
802 C: ElasticHyperviscous,
803 F: ElasticHyperviscousFiniteElement<C, G, N>,
804 Self: ViscoelasticFiniteElementBlock<C, F, G, N>,
805{
806 fn viscous_dissipation(
807 &self,
808 nodal_coordinates: &NodalCoordinatesBlock,
809 nodal_velocities: &NodalVelocitiesBlock,
810 ) -> Result<Scalar, FiniteElementBlockError> {
811 match self
812 .elements()
813 .iter()
814 .zip(self.connectivity().iter())
815 .map(|(element, element_connectivity)| {
816 element.viscous_dissipation(
817 &self.nodal_coordinates_element(element_connectivity, nodal_coordinates),
818 &self.nodal_velocities_element(element_connectivity, nodal_velocities),
819 )
820 })
821 .sum()
822 {
823 Ok(viscous_dissipation) => Ok(viscous_dissipation),
824 Err(error) => Err(FiniteElementBlockError::Upstream(
825 format!("{error}"),
826 format!("{self:?}"),
827 )),
828 }
829 }
830 fn dissipation_potential(
831 &self,
832 nodal_coordinates: &NodalCoordinatesBlock,
833 nodal_velocities: &NodalVelocitiesBlock,
834 ) -> Result<Scalar, FiniteElementBlockError> {
835 match self
836 .elements()
837 .iter()
838 .zip(self.connectivity().iter())
839 .map(|(element, element_connectivity)| {
840 element.dissipation_potential(
841 &self.nodal_coordinates_element(element_connectivity, nodal_coordinates),
842 &self.nodal_velocities_element(element_connectivity, nodal_velocities),
843 )
844 })
845 .sum()
846 {
847 Ok(dissipation_potential) => Ok(dissipation_potential),
848 Err(error) => Err(FiniteElementBlockError::Upstream(
849 format!("{error}"),
850 format!("{self:?}"),
851 )),
852 }
853 }
854 fn minimize(
855 &self,
856 equality_constraint: EqualityConstraint,
857 integrator: impl Explicit<NodalVelocitiesBlock, NodalVelocitiesHistory>,
858 time: &[Scalar],
859 solver: impl SecondOrderOptimization<
860 Scalar,
861 NodalForcesBlock,
862 NodalStiffnessesBlock,
863 NodalCoordinatesBlock,
864 >,
865 ) -> Result<(Times, NodalCoordinatesHistory, NodalVelocitiesHistory), IntegrationError> {
866 let mut solution = NodalVelocitiesBlock::zero(self.coordinates().len());
867 integrator.integrate(
868 |_: Scalar, nodal_coordinates: &NodalCoordinatesBlock| {
869 solution = self.minimize_inner(
870 equality_constraint.clone(),
871 nodal_coordinates,
872 &solver,
873 &solution,
874 )?;
875 Ok(solution.clone())
876 },
877 time,
878 self.coordinates().clone().into(),
879 )
880 }
881 #[doc(hidden)]
882 fn minimize_inner(
883 &self,
884 equality_constraint: EqualityConstraint,
885 nodal_coordinates: &NodalCoordinatesBlock,
886 solver: &impl SecondOrderOptimization<
887 Scalar,
888 NodalForcesBlock,
889 NodalStiffnessesBlock,
890 NodalCoordinatesBlock,
891 >,
892 initial_guess: &NodalVelocitiesBlock,
893 ) -> Result<NodalVelocitiesBlock, OptimizationError> {
894 let num_coords = nodal_coordinates.len();
895 let banded = band(self.connectivity(), &equality_constraint, num_coords);
896 solver.minimize(
897 |nodal_velocities: &NodalVelocitiesBlock| {
898 Ok(self.dissipation_potential(nodal_coordinates, nodal_velocities)?)
899 },
900 |nodal_velocities: &NodalVelocitiesBlock| {
901 Ok(self.nodal_forces(nodal_coordinates, nodal_velocities)?)
902 },
903 |nodal_velocities: &NodalVelocitiesBlock| {
904 Ok(self.nodal_stiffnesses(nodal_coordinates, nodal_velocities)?)
905 },
906 initial_guess.clone(),
907 equality_constraint,
908 Some(banded),
909 )
910 }
911}
912
913impl<C, F, const G: usize, const N: usize> HyperviscoelasticFiniteElementBlock<C, F, G, N>
914 for ElementBlock<F, N>
915where
916 C: Hyperviscoelastic,
917 F: HyperviscoelasticFiniteElement<C, G, N>,
918 Self: ElasticHyperviscousFiniteElementBlock<C, F, G, N>,
919{
920 fn helmholtz_free_energy(
921 &self,
922 nodal_coordinates: &NodalCoordinatesBlock,
923 ) -> Result<Scalar, FiniteElementBlockError> {
924 match self
925 .elements()
926 .iter()
927 .zip(self.connectivity().iter())
928 .map(|(element, element_connectivity)| {
929 element.helmholtz_free_energy(
930 &self.nodal_coordinates_element(element_connectivity, nodal_coordinates),
931 )
932 })
933 .sum()
934 {
935 Ok(helmholtz_free_energy) => Ok(helmholtz_free_energy),
936 Err(error) => Err(FiniteElementBlockError::Upstream(
937 format!("{error}"),
938 format!("{self:?}"),
939 )),
940 }
941 }
942}
943
944fn band<const N: usize>(
945 connectivity: &Connectivity<N>,
946 equality_constraint: &EqualityConstraint,
947 number_of_nodes: usize,
948) -> Banded {
949 match equality_constraint {
950 EqualityConstraint::Fixed(indices) => {
951 let neighbors: Vec<Vec<usize>> = invert(connectivity, number_of_nodes)
952 .iter()
953 .map(|elements| {
954 let mut nodes: Vec<usize> = elements
955 .iter()
956 .flat_map(|&element| connectivity[element])
957 .collect();
958 nodes.sort();
959 nodes.dedup();
960 nodes
961 })
962 .collect();
963 let structure: Vec<Vec<bool>> = neighbors
964 .iter()
965 .map(|nodes| (0..number_of_nodes).map(|b| nodes.contains(&b)).collect())
966 .collect();
967 let structure_3d: Vec<Vec<bool>> = structure
968 .iter()
969 .flat_map(|row| {
970 repeat_n(
971 row.iter().flat_map(|entry| repeat_n(*entry, 3)).collect(),
972 3,
973 )
974 })
975 .collect();
976 let mut keep = vec![true; structure_3d.len()];
977 indices.iter().for_each(|&index| keep[index] = false);
978 let banded = structure_3d
979 .into_iter()
980 .zip(keep.iter())
981 .filter(|(_, keep)| **keep)
982 .map(|(structure_3d_a, _)| {
983 structure_3d_a
984 .into_iter()
985 .zip(keep.iter())
986 .filter(|(_, keep)| **keep)
987 .map(|(structure_3d_ab, _)| structure_3d_ab)
988 .collect::<Vec<bool>>()
989 })
990 .collect::<Vec<Vec<bool>>>();
991 Banded::from(banded)
992 }
993 EqualityConstraint::Linear(matrix, _) => {
994 let neighbors: Vec<Vec<usize>> = invert(connectivity, number_of_nodes)
995 .iter()
996 .map(|elements| {
997 let mut nodes: Vec<usize> = elements
998 .iter()
999 .flat_map(|&element| connectivity[element])
1000 .collect();
1001 nodes.sort();
1002 nodes.dedup();
1003 nodes
1004 })
1005 .collect();
1006 let structure: Vec<Vec<bool>> = neighbors
1007 .iter()
1008 .map(|nodes| (0..number_of_nodes).map(|b| nodes.contains(&b)).collect())
1009 .collect();
1010 let structure_3d: Vec<Vec<bool>> = structure
1011 .iter()
1012 .flat_map(|row| {
1013 repeat_n(
1014 row.iter().flat_map(|entry| repeat_n(*entry, 3)).collect(),
1015 3,
1016 )
1017 })
1018 .collect();
1019 let num_coords = 3 * number_of_nodes;
1020 assert_eq!(matrix.width(), num_coords);
1021 let num_dof = matrix.len() + matrix.width();
1022 let mut banded = vec![vec![false; num_dof]; num_dof];
1023 structure_3d
1024 .iter()
1025 .zip(banded.iter_mut())
1026 .for_each(|(structure_3d_i, banded_i)| {
1027 structure_3d_i
1028 .iter()
1029 .zip(banded_i.iter_mut())
1030 .for_each(|(structure_3d_ij, banded_ij)| *banded_ij = *structure_3d_ij)
1031 });
1032 let mut index = num_coords;
1033 matrix.iter().for_each(|matrix_i| {
1034 matrix_i.iter().enumerate().for_each(|(j, matrix_ij)| {
1035 if matrix_ij != &0.0 {
1036 banded[index][j] = true;
1037 banded[j][index] = true;
1038 index += 1;
1039 }
1040 })
1041 });
1042 Banded::from(banded)
1043 }
1044 EqualityConstraint::None => {
1045 let neighbors: Vec<Vec<usize>> = invert(connectivity, number_of_nodes)
1046 .iter()
1047 .map(|elements| {
1048 let mut nodes: Vec<usize> = elements
1049 .iter()
1050 .flat_map(|&element| connectivity[element])
1051 .collect();
1052 nodes.sort();
1053 nodes.dedup();
1054 nodes
1055 })
1056 .collect();
1057 let structure: Vec<Vec<bool>> = neighbors
1058 .iter()
1059 .map(|nodes| (0..number_of_nodes).map(|b| nodes.contains(&b)).collect())
1060 .collect();
1061 let structure_3d: Vec<Vec<bool>> = structure
1062 .iter()
1063 .flat_map(|row| {
1064 repeat_n(
1065 row.iter().flat_map(|entry| repeat_n(*entry, 3)).collect(),
1066 3,
1067 )
1068 })
1069 .collect();
1070 Banded::from(structure_3d)
1071 }
1072 }
1073}
1074
1075fn invert<const N: usize>(
1076 connectivity: &Connectivity<N>,
1077 number_of_nodes: usize,
1078) -> Vec<Vec<usize>> {
1079 let mut inverse_connectivity = vec![vec![]; number_of_nodes];
1080 connectivity
1081 .iter()
1082 .enumerate()
1083 .for_each(|(element, nodes)| {
1084 nodes
1085 .iter()
1086 .for_each(|&node| inverse_connectivity[node].push(element))
1087 });
1088 inverse_connectivity
1089}