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