1#[cfg(test)]
2mod test;
3
4pub mod element;
5
6use self::element::{
7 ElasticFiniteElement, ElasticHyperviscousFiniteElement, FiniteElement, FiniteElementMethods,
8 HyperelasticFiniteElement, HyperviscoelasticFiniteElement, SurfaceFiniteElement,
9 ViscoelasticFiniteElement,
10};
11use super::*;
12use crate::{
13 math::{
14 Banded,
15 integrate::{Explicit, IntegrationError},
16 optimize::{
17 EqualityConstraint, FirstOrderRootFinding, OptimizeError, SecondOrderOptimization,
18 ZerothOrderRootFinding,
19 },
20 },
21 mechanics::Times,
22};
23use std::{array::from_fn, iter::repeat_n};
24
25pub struct ElementBlock<F, const N: usize> {
26 connectivity: Connectivity<N>,
27 coordinates: ReferenceNodalCoordinatesBlock,
28 elements: Vec<F>,
29}
30
31pub trait FiniteElementBlockMethods<C, F, const G: usize, const N: usize>
32where
33 F: FiniteElementMethods<C, G, N>,
34{
35 fn connectivity(&self) -> &Connectivity<N>;
36 fn coordinates(&self) -> &ReferenceNodalCoordinatesBlock;
37 fn deformation_gradients(
38 &self,
39 nodal_coordinates: &NodalCoordinatesBlock,
40 ) -> Vec<DeformationGradientList<G>>;
41 fn elements(&self) -> &[F];
42 fn nodal_coordinates_element(
43 &self,
44 element_connectivity: &[usize; N],
45 nodal_coordinates: &NodalCoordinatesBlock,
46 ) -> NodalCoordinates<N>;
47}
48
49pub trait FiniteElementBlock<C, F, const G: usize, const N: usize, Y>
50where
51 C: Constitutive<Y>,
52 F: FiniteElement<C, G, N, Y>,
53 Y: Parameters,
54{
55 fn new(
56 constitutive_model_parameters: Y,
57 connectivity: Connectivity<N>,
58 reference_nodal_coordinates: ReferenceNodalCoordinatesBlock,
59 ) -> Self;
60 fn reset(&mut self);
61}
62
63pub trait SurfaceFiniteElementBlock<C, F, const G: usize, const N: usize, const P: usize, Y>
64where
65 C: Constitutive<Y>,
66 F: SurfaceFiniteElement<C, G, N, P, Y>,
67 Y: Parameters,
68{
69 fn new(
70 constitutive_model_parameters: Y,
71 connectivity: Connectivity<N>,
72 reference_nodal_coordinates: ReferenceNodalCoordinatesBlock,
73 thickness: Scalar,
74 ) -> Self;
75}
76
77impl<C, F, const G: usize, const N: usize> FiniteElementBlockMethods<C, F, G, N>
78 for ElementBlock<F, N>
79where
80 F: FiniteElementMethods<C, G, N>,
81{
82 fn connectivity(&self) -> &Connectivity<N> {
83 &self.connectivity
84 }
85 fn coordinates(&self) -> &ReferenceNodalCoordinatesBlock {
86 &self.coordinates
87 }
88 fn deformation_gradients(
89 &self,
90 nodal_coordinates: &NodalCoordinatesBlock,
91 ) -> Vec<DeformationGradientList<G>> {
92 self.elements()
93 .iter()
94 .zip(self.connectivity().iter())
95 .map(|(element, element_connectivity)| {
96 element.deformation_gradients(
97 &self.nodal_coordinates_element(element_connectivity, nodal_coordinates),
98 )
99 })
100 .collect()
101 }
102 fn elements(&self) -> &[F] {
103 &self.elements
104 }
105 fn nodal_coordinates_element(
106 &self,
107 element_connectivity: &[usize; N],
108 nodal_coordinates: &NodalCoordinatesBlock,
109 ) -> NodalCoordinates<N> {
110 element_connectivity
111 .iter()
112 .map(|node| nodal_coordinates[*node].clone())
113 .collect()
114 }
115}
116
117impl<C, F, const G: usize, const N: usize, Y> FiniteElementBlock<C, F, G, N, Y>
118 for ElementBlock<F, N>
119where
120 C: Constitutive<Y>,
121 F: FiniteElement<C, G, N, Y>,
122 Y: Parameters,
123{
124 fn new(
125 constitutive_model_parameters: Y,
126 connectivity: Connectivity<N>,
127 coordinates: ReferenceNodalCoordinatesBlock,
128 ) -> Self {
129 let elements = connectivity
130 .iter()
131 .map(|element_connectivity| {
132 <F>::new(
133 constitutive_model_parameters,
134 element_connectivity
135 .iter()
136 .map(|&node| coordinates[node].clone())
137 .collect(),
138 )
139 })
140 .collect();
141 Self {
142 connectivity,
143 coordinates,
144 elements,
145 }
146 }
147 fn reset(&mut self) {
148 self.elements.iter_mut().for_each(|element| element.reset())
149 }
150}
151
152impl<C, F, const G: usize, const N: usize, const P: usize, Y>
153 SurfaceFiniteElementBlock<C, F, G, N, P, Y> for ElementBlock<F, N>
154where
155 C: Constitutive<Y>,
156 F: SurfaceFiniteElement<C, G, N, P, Y>,
157 Y: Parameters,
158{
159 fn new(
160 constitutive_model_parameters: Y,
161 connectivity: Connectivity<N>,
162 coordinates: ReferenceNodalCoordinatesBlock,
163 thickness: Scalar,
164 ) -> Self {
165 let elements = connectivity
166 .iter()
167 .map(|element_connectivity| {
168 <F>::new(
169 constitutive_model_parameters,
170 element_connectivity
171 .iter()
172 .map(|node| coordinates[*node].clone())
173 .collect(),
174 &thickness,
175 )
176 })
177 .collect();
178 Self {
179 connectivity,
180 coordinates,
181 elements,
182 }
183 }
184}
185
186pub trait ElasticFiniteElementBlock<C, F, const G: usize, const N: usize>
187where
188 C: Elastic,
189 F: ElasticFiniteElement<C, G, N>,
190{
191 fn nodal_forces(
192 &self,
193 nodal_coordinates: &NodalCoordinatesBlock,
194 ) -> Result<NodalForcesBlock, ConstitutiveError>;
195 fn nodal_stiffnesses(
196 &self,
197 nodal_coordinates: &NodalCoordinatesBlock,
198 ) -> Result<NodalStiffnessesBlock, ConstitutiveError>;
199}
200
201pub trait ZerothOrderRoot<C, F, const G: usize, const N: usize>
202where
203 C: Elastic,
204 F: ElasticFiniteElement<C, G, N>,
205{
206 fn root(
207 &self,
208 equality_constraint: EqualityConstraint,
209 solver: impl ZerothOrderRootFinding<NodalCoordinatesBlock>,
210 ) -> Result<NodalCoordinatesBlock, OptimizeError>;
211}
212
213pub trait FirstOrderRoot<C, F, const G: usize, const N: usize>
214where
215 C: Elastic,
216 F: ElasticFiniteElement<C, G, N>,
217{
218 fn root(
219 &self,
220 equality_constraint: EqualityConstraint,
221 solver: impl FirstOrderRootFinding<
222 NodalForcesBlock,
223 NodalStiffnessesBlock,
224 NodalCoordinatesBlock,
225 >,
226 ) -> Result<NodalCoordinatesBlock, OptimizeError>;
227}
228
229pub trait HyperelasticFiniteElementBlock<C, F, const G: usize, const N: usize>
230where
231 C: Hyperelastic,
232 F: HyperelasticFiniteElement<C, G, N>,
233 Self: ElasticFiniteElementBlock<C, F, G, N>,
234{
235 fn helmholtz_free_energy(
236 &self,
237 nodal_coordinates: &NodalCoordinatesBlock,
238 ) -> Result<Scalar, ConstitutiveError>;
239}
240
241pub trait SecondOrderMinimize<C, F, const G: usize, const N: usize>
242where
243 C: Hyperelastic,
244 F: HyperelasticFiniteElement<C, G, N>,
245 Self: ElasticFiniteElementBlock<C, F, G, N>,
246{
247 fn minimize(
248 &self,
249 equality_constraint: EqualityConstraint,
250 solver: impl SecondOrderOptimization<
251 Scalar,
252 NodalForcesBlock,
253 NodalStiffnessesBlock,
254 NodalCoordinatesBlock,
255 >,
256 ) -> Result<NodalCoordinatesBlock, OptimizeError>;
257}
258
259pub trait ViscoelasticFiniteElementBlock<C, F, const G: usize, const N: usize>
260where
261 C: Viscoelastic,
262 F: ViscoelasticFiniteElement<C, G, N>,
263{
264 fn deformation_gradient_rates(
265 &self,
266 nodal_coordinates: &NodalCoordinatesBlock,
267 nodal_velocities: &NodalVelocitiesBlock,
268 ) -> Vec<DeformationGradientRateList<G>>;
269 fn nodal_forces(
270 &self,
271 nodal_coordinates: &NodalCoordinatesBlock,
272 nodal_velocities: &NodalVelocitiesBlock,
273 ) -> Result<NodalForcesBlock, ConstitutiveError>;
274 fn nodal_stiffnesses(
275 &self,
276 nodal_coordinates: &NodalCoordinatesBlock,
277 nodal_velocities: &NodalVelocitiesBlock,
278 ) -> Result<NodalStiffnessesBlock, ConstitutiveError>;
279 fn nodal_velocities_element(
280 &self,
281 element_connectivity: &[usize; N],
282 nodal_velocities: &NodalVelocitiesBlock,
283 ) -> NodalVelocities<N>;
284 fn root(
285 &self,
286 equality_constraint: EqualityConstraint,
287 integrator: impl Explicit<NodalVelocitiesBlock, NodalVelocitiesHistory>,
288 time: &[Scalar],
289 solver: impl FirstOrderRootFinding<
290 NodalForcesBlock,
291 NodalStiffnessesBlock,
292 NodalCoordinatesBlock,
293 >,
294 ) -> Result<(Times, NodalCoordinatesHistory, NodalVelocitiesHistory), IntegrationError>;
295 #[doc(hidden)]
296 fn root_inner(
297 &self,
298 equality_constraint: EqualityConstraint,
299 nodal_coordinates: &NodalCoordinatesBlock,
300 solver: &impl FirstOrderRootFinding<
301 NodalForcesBlock,
302 NodalStiffnessesBlock,
303 NodalCoordinatesBlock,
304 >,
305 initial_guess: &NodalVelocitiesBlock,
306 ) -> Result<NodalVelocitiesBlock, OptimizeError>;
307}
308
309pub trait ElasticHyperviscousFiniteElementBlock<C, F, const G: usize, const N: usize>
310where
311 C: ElasticHyperviscous,
312 F: ElasticHyperviscousFiniteElement<C, G, N>,
313 Self: ViscoelasticFiniteElementBlock<C, F, G, N>,
314{
315 fn viscous_dissipation(
316 &self,
317 nodal_coordinates: &NodalCoordinatesBlock,
318 nodal_velocities: &NodalVelocitiesBlock,
319 ) -> Result<Scalar, ConstitutiveError>;
320 fn dissipation_potential(
321 &self,
322 nodal_coordinates: &NodalCoordinatesBlock,
323 nodal_velocities: &NodalVelocitiesBlock,
324 ) -> Result<Scalar, ConstitutiveError>;
325 fn minimize(
326 &self,
327 equality_constraint: EqualityConstraint,
328 integrator: impl Explicit<NodalVelocitiesBlock, NodalVelocitiesHistory>,
329 time: &[Scalar],
330 solver: impl SecondOrderOptimization<
331 Scalar,
332 NodalForcesBlock,
333 NodalStiffnessesBlock,
334 NodalCoordinatesBlock,
335 >,
336 ) -> Result<(Times, NodalCoordinatesHistory, NodalVelocitiesHistory), IntegrationError>;
337 #[doc(hidden)]
338 fn minimize_inner(
339 &self,
340 equality_constraint: EqualityConstraint,
341 nodal_coordinates: &NodalCoordinatesBlock,
342 solver: &impl SecondOrderOptimization<
343 Scalar,
344 NodalForcesBlock,
345 NodalStiffnessesBlock,
346 NodalCoordinatesBlock,
347 >,
348 initial_guess: &NodalVelocitiesBlock,
349 ) -> Result<NodalVelocitiesBlock, OptimizeError>;
350}
351
352pub trait HyperviscoelasticFiniteElementBlock<C, F, const G: usize, const N: usize>
353where
354 C: Hyperviscoelastic,
355 F: HyperviscoelasticFiniteElement<C, G, N>,
356 Self: ElasticHyperviscousFiniteElementBlock<C, F, G, N>,
357{
358 fn helmholtz_free_energy(
359 &self,
360 nodal_coordinates: &NodalCoordinatesBlock,
361 ) -> Result<Scalar, ConstitutiveError>;
362}
363
364impl<C, F, const G: usize, const N: usize> ElasticFiniteElementBlock<C, F, G, N>
365 for ElementBlock<F, N>
366where
367 C: Elastic,
368 F: ElasticFiniteElement<C, G, N>,
369 Self: FiniteElementBlockMethods<C, F, G, N>,
370{
371 fn nodal_forces(
372 &self,
373 nodal_coordinates: &NodalCoordinatesBlock,
374 ) -> Result<NodalForcesBlock, ConstitutiveError> {
375 let mut nodal_forces = NodalForcesBlock::zero(nodal_coordinates.len());
376 self.elements()
377 .iter()
378 .zip(self.connectivity().iter())
379 .try_for_each(|(element, element_connectivity)| {
380 element
381 .nodal_forces(
382 &self.nodal_coordinates_element(element_connectivity, nodal_coordinates),
383 )?
384 .iter()
385 .zip(element_connectivity.iter())
386 .for_each(|(nodal_force, &node)| nodal_forces[node] += nodal_force);
387 Ok::<(), ConstitutiveError>(())
388 })?;
389 Ok(nodal_forces)
390 }
391 fn nodal_stiffnesses(
392 &self,
393 nodal_coordinates: &NodalCoordinatesBlock,
394 ) -> Result<NodalStiffnessesBlock, ConstitutiveError> {
395 let mut nodal_stiffnesses = NodalStiffnessesBlock::zero(nodal_coordinates.len());
396 self.elements()
397 .iter()
398 .zip(self.connectivity().iter())
399 .try_for_each(|(element, element_connectivity)| {
400 element
401 .nodal_stiffnesses(
402 &self.nodal_coordinates_element(element_connectivity, nodal_coordinates),
403 )?
404 .iter()
405 .zip(element_connectivity.iter())
406 .for_each(|(object, &node_a)| {
407 object.iter().zip(element_connectivity.iter()).for_each(
408 |(nodal_stiffness, &node_b)| {
409 nodal_stiffnesses[node_a][node_b] += nodal_stiffness
410 },
411 )
412 });
413 Ok::<(), ConstitutiveError>(())
414 })?;
415 Ok(nodal_stiffnesses)
416 }
417}
418
419impl<C, F, const G: usize, const N: usize> FirstOrderRoot<C, F, G, N> for ElementBlock<F, N>
439where
440 C: Elastic,
441 F: ElasticFiniteElement<C, G, N>,
442 Self: FiniteElementBlockMethods<C, F, G, N>,
443{
444 fn root(
445 &self,
446 equality_constraint: EqualityConstraint,
447 solver: impl FirstOrderRootFinding<
448 NodalForcesBlock,
449 NodalStiffnessesBlock,
450 NodalCoordinatesBlock,
451 >,
452 ) -> Result<NodalCoordinatesBlock, OptimizeError> {
453 solver.root(
454 |nodal_coordinates: &NodalCoordinatesBlock| Ok(self.nodal_forces(nodal_coordinates)?),
455 |nodal_coordinates: &NodalCoordinatesBlock| {
456 Ok(self.nodal_stiffnesses(nodal_coordinates)?)
457 },
458 self.coordinates().clone().into(),
459 equality_constraint,
460 )
461 }
462}
463
464impl<C, F, const G: usize, const N: usize> HyperelasticFiniteElementBlock<C, F, G, N>
465 for ElementBlock<F, N>
466where
467 C: Hyperelastic,
468 F: HyperelasticFiniteElement<C, G, N>,
469 Self: ElasticFiniteElementBlock<C, F, G, N>,
470{
471 fn helmholtz_free_energy(
472 &self,
473 nodal_coordinates: &NodalCoordinatesBlock,
474 ) -> Result<Scalar, ConstitutiveError> {
475 self.elements()
476 .iter()
477 .zip(self.connectivity().iter())
478 .map(|(element, element_connectivity)| {
479 element.helmholtz_free_energy(
480 &self.nodal_coordinates_element(element_connectivity, nodal_coordinates),
481 )
482 })
483 .sum()
484 }
485}
486
487impl<C, F, const G: usize, const N: usize> SecondOrderMinimize<C, F, G, N> for ElementBlock<F, N>
488where
489 C: Hyperelastic,
490 F: HyperelasticFiniteElement<C, G, N>,
491 Self: ElasticFiniteElementBlock<C, F, G, N>,
492{
493 fn minimize(
494 &self,
495 equality_constraint: EqualityConstraint,
496 solver: impl SecondOrderOptimization<
497 Scalar,
498 NodalForcesBlock,
499 NodalStiffnessesBlock,
500 NodalCoordinatesBlock,
501 >,
502 ) -> Result<NodalCoordinatesBlock, OptimizeError> {
503 let banded = band(
504 self.connectivity(),
505 &equality_constraint,
506 self.coordinates().len(),
507 );
508 solver.minimize(
509 |nodal_coordinates: &NodalCoordinatesBlock| {
510 Ok(self.helmholtz_free_energy(nodal_coordinates)?)
511 },
512 |nodal_coordinates: &NodalCoordinatesBlock| Ok(self.nodal_forces(nodal_coordinates)?),
513 |nodal_coordinates: &NodalCoordinatesBlock| {
514 Ok(self.nodal_stiffnesses(nodal_coordinates)?)
515 },
516 self.coordinates().clone().into(),
517 equality_constraint,
518 Some(banded),
519 )
520 }
521}
522
523impl<C, F, const G: usize, const N: usize> ViscoelasticFiniteElementBlock<C, F, G, N>
524 for ElementBlock<F, N>
525where
526 C: Viscoelastic,
527 F: ViscoelasticFiniteElement<C, G, N>,
528 Self: FiniteElementBlockMethods<C, F, G, N>,
529{
530 fn deformation_gradient_rates(
531 &self,
532 nodal_coordinates: &NodalCoordinatesBlock,
533 nodal_velocities: &NodalVelocitiesBlock,
534 ) -> Vec<DeformationGradientRateList<G>> {
535 self.elements()
536 .iter()
537 .zip(self.connectivity().iter())
538 .map(|(element, element_connectivity)| {
539 element.deformation_gradient_rates(
540 &self.nodal_coordinates_element(element_connectivity, nodal_coordinates),
541 &self.nodal_velocities_element(element_connectivity, nodal_velocities),
542 )
543 })
544 .collect()
545 }
546 fn nodal_forces(
547 &self,
548 nodal_coordinates: &NodalCoordinatesBlock,
549 nodal_velocities: &NodalVelocitiesBlock,
550 ) -> Result<NodalForcesBlock, ConstitutiveError> {
551 let mut nodal_forces = NodalForcesBlock::zero(nodal_coordinates.len());
552 self.elements()
553 .iter()
554 .zip(self.connectivity().iter())
555 .try_for_each(|(element, element_connectivity)| {
556 element
557 .nodal_forces(
558 &self.nodal_coordinates_element(element_connectivity, nodal_coordinates),
559 &self.nodal_velocities_element(element_connectivity, nodal_velocities),
560 )?
561 .iter()
562 .zip(element_connectivity.iter())
563 .for_each(|(nodal_force, &node)| nodal_forces[node] += nodal_force);
564 Ok::<(), ConstitutiveError>(())
565 })?;
566 Ok(nodal_forces)
567 }
568 fn nodal_stiffnesses(
569 &self,
570 nodal_coordinates: &NodalCoordinatesBlock,
571 nodal_velocities: &NodalVelocitiesBlock,
572 ) -> Result<NodalStiffnessesBlock, ConstitutiveError> {
573 let mut nodal_stiffnesses = NodalStiffnessesBlock::zero(nodal_coordinates.len());
574 self.elements()
575 .iter()
576 .zip(self.connectivity().iter())
577 .try_for_each(|(element, element_connectivity)| {
578 element
579 .nodal_stiffnesses(
580 &self.nodal_coordinates_element(element_connectivity, nodal_coordinates),
581 &self.nodal_velocities_element(element_connectivity, nodal_velocities),
582 )?
583 .iter()
584 .zip(element_connectivity.iter())
585 .for_each(|(object, &node_a)| {
586 object.iter().zip(element_connectivity.iter()).for_each(
587 |(nodal_stiffness, &node_b)| {
588 nodal_stiffnesses[node_a][node_b] += nodal_stiffness
589 },
590 )
591 });
592 Ok::<(), ConstitutiveError>(())
593 })?;
594 Ok(nodal_stiffnesses)
595 }
596 fn nodal_velocities_element(
597 &self,
598 element_connectivity: &[usize; N],
599 nodal_velocities: &NodalVelocitiesBlock,
600 ) -> NodalVelocities<N> {
601 element_connectivity
602 .iter()
603 .map(|node| nodal_velocities[*node].clone())
604 .collect()
605 }
606 fn root(
607 &self,
608 equality_constraint: EqualityConstraint,
609 integrator: impl Explicit<NodalVelocitiesBlock, NodalVelocitiesHistory>,
610 time: &[Scalar],
611 solver: impl FirstOrderRootFinding<
612 NodalForcesBlock,
613 NodalStiffnessesBlock,
614 NodalCoordinatesBlock,
615 >,
616 ) -> Result<(Times, NodalCoordinatesHistory, NodalVelocitiesHistory), IntegrationError> {
617 let mut solution = NodalVelocitiesBlock::zero(self.coordinates().len());
618 integrator.integrate(
619 |_: Scalar, nodal_coordinates: &NodalCoordinatesBlock| {
620 solution = self.root_inner(
621 equality_constraint.clone(),
622 nodal_coordinates,
623 &solver,
624 &solution,
625 )?;
626 Ok(solution.clone())
627 },
628 time,
629 self.coordinates().clone().into(),
630 )
631 }
632 #[doc(hidden)]
633 fn root_inner(
634 &self,
635 equality_constraint: EqualityConstraint,
636 nodal_coordinates: &NodalCoordinatesBlock,
637 solver: &impl FirstOrderRootFinding<
638 NodalForcesBlock,
639 NodalStiffnessesBlock,
640 NodalCoordinatesBlock,
641 >,
642 initial_guess: &NodalVelocitiesBlock,
643 ) -> Result<NodalVelocitiesBlock, OptimizeError> {
644 solver.root(
645 |nodal_velocities: &NodalVelocitiesBlock| {
646 Ok(self.nodal_forces(nodal_coordinates, nodal_velocities)?)
647 },
648 |nodal_velocities: &NodalVelocitiesBlock| {
649 Ok(self.nodal_stiffnesses(nodal_coordinates, nodal_velocities)?)
650 },
651 initial_guess.clone(),
652 equality_constraint,
653 )
654 }
655}
656
657impl<C, F, const G: usize, const N: usize> ElasticHyperviscousFiniteElementBlock<C, F, G, N>
658 for ElementBlock<F, N>
659where
660 C: ElasticHyperviscous,
661 F: ElasticHyperviscousFiniteElement<C, G, N>,
662 Self: ViscoelasticFiniteElementBlock<C, F, G, N>,
663{
664 fn viscous_dissipation(
665 &self,
666 nodal_coordinates: &NodalCoordinatesBlock,
667 nodal_velocities: &NodalVelocitiesBlock,
668 ) -> Result<Scalar, ConstitutiveError> {
669 self.elements()
670 .iter()
671 .zip(self.connectivity().iter())
672 .map(|(element, element_connectivity)| {
673 element.viscous_dissipation(
674 &self.nodal_coordinates_element(element_connectivity, nodal_coordinates),
675 &self.nodal_velocities_element(element_connectivity, nodal_velocities),
676 )
677 })
678 .sum()
679 }
680 fn dissipation_potential(
681 &self,
682 nodal_coordinates: &NodalCoordinatesBlock,
683 nodal_velocities: &NodalVelocitiesBlock,
684 ) -> Result<Scalar, ConstitutiveError> {
685 self.elements()
686 .iter()
687 .zip(self.connectivity().iter())
688 .map(|(element, element_connectivity)| {
689 element.dissipation_potential(
690 &self.nodal_coordinates_element(element_connectivity, nodal_coordinates),
691 &self.nodal_velocities_element(element_connectivity, nodal_velocities),
692 )
693 })
694 .sum()
695 }
696 fn minimize(
697 &self,
698 equality_constraint: EqualityConstraint,
699 integrator: impl Explicit<NodalVelocitiesBlock, NodalVelocitiesHistory>,
700 time: &[Scalar],
701 solver: impl SecondOrderOptimization<
702 Scalar,
703 NodalForcesBlock,
704 NodalStiffnessesBlock,
705 NodalCoordinatesBlock,
706 >,
707 ) -> Result<(Times, NodalCoordinatesHistory, NodalVelocitiesHistory), IntegrationError> {
708 let mut solution = NodalVelocitiesBlock::zero(self.coordinates().len());
709 integrator.integrate(
710 |_: Scalar, nodal_coordinates: &NodalCoordinatesBlock| {
711 solution = self.minimize_inner(
712 equality_constraint.clone(),
713 nodal_coordinates,
714 &solver,
715 &solution,
716 )?;
717 Ok(solution.clone())
718 },
719 time,
720 self.coordinates().clone().into(),
721 )
722 }
723 #[doc(hidden)]
724 fn minimize_inner(
725 &self,
726 equality_constraint: EqualityConstraint,
727 nodal_coordinates: &NodalCoordinatesBlock,
728 solver: &impl SecondOrderOptimization<
729 Scalar,
730 NodalForcesBlock,
731 NodalStiffnessesBlock,
732 NodalCoordinatesBlock,
733 >,
734 initial_guess: &NodalVelocitiesBlock,
735 ) -> Result<NodalVelocitiesBlock, OptimizeError> {
736 let num_coords = nodal_coordinates.len();
737 let banded = band(self.connectivity(), &equality_constraint, num_coords);
738 solver.minimize(
739 |nodal_velocities: &NodalVelocitiesBlock| {
740 Ok(self.dissipation_potential(nodal_coordinates, nodal_velocities)?)
741 },
742 |nodal_velocities: &NodalVelocitiesBlock| {
743 Ok(self.nodal_forces(nodal_coordinates, nodal_velocities)?)
744 },
745 |nodal_velocities: &NodalVelocitiesBlock| {
746 Ok(self.nodal_stiffnesses(nodal_coordinates, nodal_velocities)?)
747 },
748 initial_guess.clone(),
749 equality_constraint,
750 Some(banded),
751 )
752 }
753}
754
755impl<C, F, const G: usize, const N: usize> HyperviscoelasticFiniteElementBlock<C, F, G, N>
756 for ElementBlock<F, N>
757where
758 C: Hyperviscoelastic,
759 F: HyperviscoelasticFiniteElement<C, G, N>,
760 Self: ElasticHyperviscousFiniteElementBlock<C, F, G, N>,
761{
762 fn helmholtz_free_energy(
763 &self,
764 nodal_coordinates: &NodalCoordinatesBlock,
765 ) -> Result<Scalar, ConstitutiveError> {
766 self.elements()
767 .iter()
768 .zip(self.connectivity().iter())
769 .map(|(element, element_connectivity)| {
770 element.helmholtz_free_energy(
771 &self.nodal_coordinates_element(element_connectivity, nodal_coordinates),
772 )
773 })
774 .sum()
775 }
776}
777
778fn band<const N: usize>(
779 connectivity: &Connectivity<N>,
780 equality_constraint: &EqualityConstraint,
781 number_of_nodes: usize,
782) -> Banded {
783 match equality_constraint {
784 EqualityConstraint::Linear(matrix, _) => {
785 let neighbors: Vec<Vec<usize>> = invert(connectivity, number_of_nodes)
786 .iter()
787 .map(|elements| {
788 let mut nodes: Vec<usize> = elements
789 .iter()
790 .flat_map(|&element| connectivity[element])
791 .collect();
792 nodes.sort();
793 nodes.dedup();
794 nodes
795 })
796 .collect();
797 let structure: Vec<Vec<bool>> = neighbors
798 .iter()
799 .map(|nodes| (0..number_of_nodes).map(|b| nodes.contains(&b)).collect())
800 .collect();
801 let structure_3d: Vec<Vec<bool>> = structure
802 .iter()
803 .flat_map(|row| {
804 repeat_n(
805 row.iter().flat_map(|entry| repeat_n(*entry, 3)).collect(),
806 3,
807 )
808 })
809 .collect();
810 let num_coords = 3 * number_of_nodes;
811 assert_eq!(matrix.width(), num_coords);
812 let num_dof = matrix.len() + matrix.width();
813 let mut banded = vec![vec![false; num_dof]; num_dof];
814 structure_3d
815 .iter()
816 .zip(banded.iter_mut())
817 .for_each(|(structure_3d_i, banded_i)| {
818 structure_3d_i
819 .iter()
820 .zip(banded_i.iter_mut())
821 .for_each(|(structure_3d_ij, banded_ij)| *banded_ij = *structure_3d_ij)
822 });
823 let mut index = num_coords;
824 matrix.iter().for_each(|matrix_i| {
825 matrix_i.iter().enumerate().for_each(|(j, matrix_ij)| {
826 if matrix_ij != &0.0 {
827 banded[index][j] = true;
828 banded[j][index] = true;
829 index += 1;
830 }
831 })
832 });
833 Banded::from(banded)
834 }
835 _ => unimplemented!(),
836 }
837}
838
839fn invert<const N: usize>(
840 connectivity: &Connectivity<N>,
841 number_of_nodes: usize,
842) -> Vec<Vec<usize>> {
843 let mut inverse_connectivity = vec![vec![]; number_of_nodes];
844 connectivity
845 .iter()
846 .enumerate()
847 .for_each(|(element, nodes)| {
848 nodes
849 .iter()
850 .for_each(|&node| inverse_connectivity[node].push(element))
851 });
852 inverse_connectivity
853}