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, NewtonRaphson, OptimizeError,
18 SecondOrderOptimization,
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}
61
62pub trait SurfaceFiniteElementBlock<C, F, const G: usize, const N: usize, const P: usize, Y>
63where
64 C: Constitutive<Y>,
65 F: SurfaceFiniteElement<C, G, N, P, Y>,
66 Y: Parameters,
67{
68 fn new(
69 constitutive_model_parameters: Y,
70 connectivity: Connectivity<N>,
71 reference_nodal_coordinates: ReferenceNodalCoordinatesBlock,
72 thickness: Scalar,
73 ) -> Self;
74}
75
76impl<C, F, const G: usize, const N: usize> FiniteElementBlockMethods<C, F, G, N>
77 for ElementBlock<F, N>
78where
79 F: FiniteElementMethods<C, G, N>,
80{
81 fn connectivity(&self) -> &Connectivity<N> {
82 &self.connectivity
83 }
84 fn coordinates(&self) -> &ReferenceNodalCoordinatesBlock {
85 &self.coordinates
86 }
87 fn deformation_gradients(
88 &self,
89 nodal_coordinates: &NodalCoordinatesBlock,
90 ) -> Vec<DeformationGradientList<G>> {
91 self.elements()
92 .iter()
93 .zip(self.connectivity().iter())
94 .map(|(element, element_connectivity)| {
95 element.deformation_gradients(
96 &self.nodal_coordinates_element(element_connectivity, nodal_coordinates),
97 )
98 })
99 .collect()
100 }
101 fn elements(&self) -> &[F] {
102 &self.elements
103 }
104 fn nodal_coordinates_element(
105 &self,
106 element_connectivity: &[usize; N],
107 nodal_coordinates: &NodalCoordinatesBlock,
108 ) -> NodalCoordinates<N> {
109 element_connectivity
110 .iter()
111 .map(|node| nodal_coordinates[*node].clone())
112 .collect()
113 }
114}
115
116impl<C, F, const G: usize, const N: usize, Y> FiniteElementBlock<C, F, G, N, Y>
117 for ElementBlock<F, N>
118where
119 C: Constitutive<Y>,
120 F: FiniteElement<C, G, N, Y>,
121 Y: Parameters,
122{
123 fn new(
124 constitutive_model_parameters: Y,
125 connectivity: Connectivity<N>,
126 coordinates: ReferenceNodalCoordinatesBlock,
127 ) -> Self {
128 let elements = connectivity
129 .iter()
130 .map(|element_connectivity| {
131 <F>::new(
132 constitutive_model_parameters,
133 element_connectivity
134 .iter()
135 .map(|&node| coordinates[node].clone())
136 .collect(),
137 )
138 })
139 .collect();
140 Self {
141 connectivity,
142 coordinates,
143 elements,
144 }
145 }
146}
147
148impl<C, F, const G: usize, const N: usize, const P: usize, Y>
149 SurfaceFiniteElementBlock<C, F, G, N, P, Y> for ElementBlock<F, N>
150where
151 C: Constitutive<Y>,
152 F: SurfaceFiniteElement<C, G, N, P, Y>,
153 Y: Parameters,
154{
155 fn new(
156 constitutive_model_parameters: Y,
157 connectivity: Connectivity<N>,
158 coordinates: ReferenceNodalCoordinatesBlock,
159 thickness: Scalar,
160 ) -> Self {
161 let elements = connectivity
162 .iter()
163 .map(|element_connectivity| {
164 <F>::new(
165 constitutive_model_parameters,
166 element_connectivity
167 .iter()
168 .map(|node| coordinates[*node].clone())
169 .collect(),
170 &thickness,
171 )
172 })
173 .collect();
174 Self {
175 connectivity,
176 coordinates,
177 elements,
178 }
179 }
180}
181
182pub trait ElasticFiniteElementBlock<C, F, const G: usize, const N: usize>
183where
184 C: Elastic,
185 F: ElasticFiniteElement<C, G, N>,
186{
187 fn nodal_forces(
188 &self,
189 nodal_coordinates: &NodalCoordinatesBlock,
190 ) -> Result<NodalForcesBlock, ConstitutiveError>;
191 fn nodal_stiffnesses(
192 &self,
193 nodal_coordinates: &NodalCoordinatesBlock,
194 ) -> Result<NodalStiffnessesBlock, ConstitutiveError>;
195 fn root(
196 &self,
197 equality_constraint: EqualityConstraint,
198 ) -> Result<NodalCoordinatesBlock, OptimizeError>;
199}
200
201pub trait HyperelasticFiniteElementBlock<C, F, const G: usize, const N: usize>
202where
203 C: Hyperelastic,
204 F: HyperelasticFiniteElement<C, G, N>,
205 Self: ElasticFiniteElementBlock<C, F, G, N>,
206{
207 fn helmholtz_free_energy(
208 &self,
209 nodal_coordinates: &NodalCoordinatesBlock,
210 ) -> Result<Scalar, ConstitutiveError>;
211 fn minimize(
212 &self,
213 equality_constraint: EqualityConstraint,
214 ) -> Result<NodalCoordinatesBlock, OptimizeError>;
215}
216
217pub trait ViscoelasticFiniteElementBlock<C, F, const G: usize, const N: usize>
218where
219 C: Viscoelastic,
220 F: ViscoelasticFiniteElement<C, G, N>,
221{
222 fn deformation_gradient_rates(
223 &self,
224 nodal_coordinates: &NodalCoordinatesBlock,
225 nodal_velocities: &NodalVelocitiesBlock,
226 ) -> Vec<DeformationGradientRateList<G>>;
227 fn nodal_forces(
228 &self,
229 nodal_coordinates: &NodalCoordinatesBlock,
230 nodal_velocities: &NodalVelocitiesBlock,
231 ) -> Result<NodalForcesBlock, ConstitutiveError>;
232 fn nodal_stiffnesses(
233 &self,
234 nodal_coordinates: &NodalCoordinatesBlock,
235 nodal_velocities: &NodalVelocitiesBlock,
236 ) -> Result<NodalStiffnessesBlock, ConstitutiveError>;
237 fn nodal_velocities_element(
238 &self,
239 element_connectivity: &[usize; N],
240 nodal_velocities: &NodalVelocitiesBlock,
241 ) -> NodalVelocities<N>;
242 fn root(
243 &self,
244 equality_constraint: EqualityConstraint,
245 integrator: impl Explicit<NodalVelocitiesBlock, NodalVelocitiesHistory>,
246 time: &[Scalar],
247 ) -> Result<(Times, NodalCoordinatesHistory, NodalVelocitiesHistory), IntegrationError>;
248 #[doc(hidden)]
249 fn root_inner(
250 &self,
251 equality_constraint: EqualityConstraint,
252 nodal_coordinates: &NodalCoordinatesBlock,
253 ) -> Result<NodalVelocitiesBlock, OptimizeError>;
254}
255
256pub trait ElasticHyperviscousFiniteElementBlock<C, F, const G: usize, const N: usize>
257where
258 C: ElasticHyperviscous,
259 F: ElasticHyperviscousFiniteElement<C, G, N>,
260 Self: ViscoelasticFiniteElementBlock<C, F, G, N>,
261{
262 fn viscous_dissipation(
263 &self,
264 nodal_coordinates: &NodalCoordinatesBlock,
265 nodal_velocities: &NodalVelocitiesBlock,
266 ) -> Result<Scalar, ConstitutiveError>;
267 fn dissipation_potential(
268 &self,
269 nodal_coordinates: &NodalCoordinatesBlock,
270 nodal_velocities: &NodalVelocitiesBlock,
271 ) -> Result<Scalar, ConstitutiveError>;
272 fn minimize(
273 &self,
274 equality_constraint: EqualityConstraint,
275 integrator: impl Explicit<NodalVelocitiesBlock, NodalVelocitiesHistory>,
276 time: &[Scalar],
277 ) -> Result<(Times, NodalCoordinatesHistory, NodalVelocitiesHistory), IntegrationError>;
278 #[doc(hidden)]
279 fn minimize_inner(
280 &self,
281 equality_constraint: EqualityConstraint,
282 nodal_coordinates: &NodalCoordinatesBlock,
283 ) -> Result<NodalVelocitiesBlock, OptimizeError>;
284}
285
286pub trait HyperviscoelasticFiniteElementBlock<C, F, const G: usize, const N: usize>
287where
288 C: Hyperviscoelastic,
289 F: HyperviscoelasticFiniteElement<C, G, N>,
290 Self: ElasticHyperviscousFiniteElementBlock<C, F, G, N>,
291{
292 fn helmholtz_free_energy(
293 &self,
294 nodal_coordinates: &NodalCoordinatesBlock,
295 ) -> Result<Scalar, ConstitutiveError>;
296}
297
298impl<C, F, const G: usize, const N: usize> ElasticFiniteElementBlock<C, F, G, N>
299 for ElementBlock<F, N>
300where
301 C: Elastic,
302 F: ElasticFiniteElement<C, G, N>,
303 Self: FiniteElementBlockMethods<C, F, G, N>,
304{
305 fn nodal_forces(
306 &self,
307 nodal_coordinates: &NodalCoordinatesBlock,
308 ) -> Result<NodalForcesBlock, ConstitutiveError> {
309 let mut nodal_forces = NodalForcesBlock::zero(nodal_coordinates.len());
310 self.elements()
311 .iter()
312 .zip(self.connectivity().iter())
313 .try_for_each(|(element, element_connectivity)| {
314 element
315 .nodal_forces(
316 &self.nodal_coordinates_element(element_connectivity, nodal_coordinates),
317 )?
318 .iter()
319 .zip(element_connectivity.iter())
320 .for_each(|(nodal_force, &node)| nodal_forces[node] += nodal_force);
321 Ok::<(), ConstitutiveError>(())
322 })?;
323 Ok(nodal_forces)
324 }
325 fn nodal_stiffnesses(
326 &self,
327 nodal_coordinates: &NodalCoordinatesBlock,
328 ) -> Result<NodalStiffnessesBlock, ConstitutiveError> {
329 let mut nodal_stiffnesses = NodalStiffnessesBlock::zero(nodal_coordinates.len());
330 self.elements()
331 .iter()
332 .zip(self.connectivity().iter())
333 .try_for_each(|(element, element_connectivity)| {
334 element
335 .nodal_stiffnesses(
336 &self.nodal_coordinates_element(element_connectivity, nodal_coordinates),
337 )?
338 .iter()
339 .zip(element_connectivity.iter())
340 .for_each(|(object, &node_a)| {
341 object.iter().zip(element_connectivity.iter()).for_each(
342 |(nodal_stiffness, &node_b)| {
343 nodal_stiffnesses[node_a][node_b] += nodal_stiffness
344 },
345 )
346 });
347 Ok::<(), ConstitutiveError>(())
348 })?;
349 Ok(nodal_stiffnesses)
350 }
351 fn root(
352 &self,
353 equality_constraint: EqualityConstraint,
354 ) -> Result<NodalCoordinatesBlock, OptimizeError> {
355 let solver = NewtonRaphson {
356 ..Default::default()
357 };
358 solver.root(
359 |nodal_coordinates: &NodalCoordinatesBlock| Ok(self.nodal_forces(nodal_coordinates)?),
360 |nodal_coordinates: &NodalCoordinatesBlock| {
361 Ok(self.nodal_stiffnesses(nodal_coordinates)?)
362 },
363 self.coordinates().clone().into(),
364 equality_constraint,
365 )
366 }
367}
368
369impl<C, F, const G: usize, const N: usize> HyperelasticFiniteElementBlock<C, F, G, N>
370 for ElementBlock<F, N>
371where
372 C: Hyperelastic,
373 F: HyperelasticFiniteElement<C, G, N>,
374 Self: ElasticFiniteElementBlock<C, F, G, N>,
375{
376 fn helmholtz_free_energy(
377 &self,
378 nodal_coordinates: &NodalCoordinatesBlock,
379 ) -> Result<Scalar, ConstitutiveError> {
380 self.elements()
381 .iter()
382 .zip(self.connectivity().iter())
383 .map(|(element, element_connectivity)| {
384 element.helmholtz_free_energy(
385 &self.nodal_coordinates_element(element_connectivity, nodal_coordinates),
386 )
387 })
388 .sum()
389 }
390 fn minimize(
391 &self,
392 equality_constraint: EqualityConstraint,
393 ) -> Result<NodalCoordinatesBlock, OptimizeError> {
394 let banded = band(
395 self.connectivity(),
396 &equality_constraint,
397 self.coordinates().len(),
398 );
399 let solver = NewtonRaphson {
400 ..Default::default()
401 };
402 solver.minimize(
403 |nodal_coordinates: &NodalCoordinatesBlock| {
404 Ok(self.helmholtz_free_energy(nodal_coordinates)?)
405 },
406 |nodal_coordinates: &NodalCoordinatesBlock| Ok(self.nodal_forces(nodal_coordinates)?),
407 |nodal_coordinates: &NodalCoordinatesBlock| {
408 Ok(self.nodal_stiffnesses(nodal_coordinates)?)
409 },
410 self.coordinates().clone().into(),
411 equality_constraint,
412 Some(banded),
413 )
414 }
415}
416
417impl<C, F, const G: usize, const N: usize> ViscoelasticFiniteElementBlock<C, F, G, N>
418 for ElementBlock<F, N>
419where
420 C: Viscoelastic,
421 F: ViscoelasticFiniteElement<C, G, N>,
422 Self: FiniteElementBlockMethods<C, F, G, N>,
423{
424 fn deformation_gradient_rates(
425 &self,
426 nodal_coordinates: &NodalCoordinatesBlock,
427 nodal_velocities: &NodalVelocitiesBlock,
428 ) -> Vec<DeformationGradientRateList<G>> {
429 self.elements()
430 .iter()
431 .zip(self.connectivity().iter())
432 .map(|(element, element_connectivity)| {
433 element.deformation_gradient_rates(
434 &self.nodal_coordinates_element(element_connectivity, nodal_coordinates),
435 &self.nodal_velocities_element(element_connectivity, nodal_velocities),
436 )
437 })
438 .collect()
439 }
440 fn nodal_forces(
441 &self,
442 nodal_coordinates: &NodalCoordinatesBlock,
443 nodal_velocities: &NodalVelocitiesBlock,
444 ) -> Result<NodalForcesBlock, ConstitutiveError> {
445 let mut nodal_forces = NodalForcesBlock::zero(nodal_coordinates.len());
446 self.elements()
447 .iter()
448 .zip(self.connectivity().iter())
449 .try_for_each(|(element, element_connectivity)| {
450 element
451 .nodal_forces(
452 &self.nodal_coordinates_element(element_connectivity, nodal_coordinates),
453 &self.nodal_velocities_element(element_connectivity, nodal_velocities),
454 )?
455 .iter()
456 .zip(element_connectivity.iter())
457 .for_each(|(nodal_force, &node)| nodal_forces[node] += nodal_force);
458 Ok::<(), ConstitutiveError>(())
459 })?;
460 Ok(nodal_forces)
461 }
462 fn nodal_stiffnesses(
463 &self,
464 nodal_coordinates: &NodalCoordinatesBlock,
465 nodal_velocities: &NodalVelocitiesBlock,
466 ) -> Result<NodalStiffnessesBlock, ConstitutiveError> {
467 let mut nodal_stiffnesses = NodalStiffnessesBlock::zero(nodal_coordinates.len());
468 self.elements()
469 .iter()
470 .zip(self.connectivity().iter())
471 .try_for_each(|(element, element_connectivity)| {
472 element
473 .nodal_stiffnesses(
474 &self.nodal_coordinates_element(element_connectivity, nodal_coordinates),
475 &self.nodal_velocities_element(element_connectivity, nodal_velocities),
476 )?
477 .iter()
478 .zip(element_connectivity.iter())
479 .for_each(|(object, &node_a)| {
480 object.iter().zip(element_connectivity.iter()).for_each(
481 |(nodal_stiffness, &node_b)| {
482 nodal_stiffnesses[node_a][node_b] += nodal_stiffness
483 },
484 )
485 });
486 Ok::<(), ConstitutiveError>(())
487 })?;
488 Ok(nodal_stiffnesses)
489 }
490 fn nodal_velocities_element(
491 &self,
492 element_connectivity: &[usize; N],
493 nodal_velocities: &NodalVelocitiesBlock,
494 ) -> NodalVelocities<N> {
495 element_connectivity
496 .iter()
497 .map(|node| nodal_velocities[*node].clone())
498 .collect()
499 }
500 fn root(
501 &self,
502 equality_constraint: EqualityConstraint,
503 integrator: impl Explicit<NodalVelocitiesBlock, NodalVelocitiesHistory>,
504 time: &[Scalar],
505 ) -> Result<(Times, NodalCoordinatesHistory, NodalVelocitiesHistory), IntegrationError> {
506 integrator.integrate(
507 |_: Scalar, nodal_coordinates: &NodalCoordinatesBlock| {
508 Ok(self.root_inner(equality_constraint.clone(), nodal_coordinates)?)
509 },
510 time,
511 self.coordinates().clone().into(),
512 )
513 }
514 #[doc(hidden)]
515 fn root_inner(
516 &self,
517 equality_constraint: EqualityConstraint,
518 nodal_coordinates: &NodalCoordinatesBlock,
519 ) -> Result<NodalVelocitiesBlock, OptimizeError> {
520 let num_coords = nodal_coordinates.len();
521 let solver = NewtonRaphson {
522 ..Default::default()
523 };
524 solver.root(
525 |nodal_velocities: &NodalVelocitiesBlock| {
526 Ok(self.nodal_forces(nodal_coordinates, nodal_velocities)?)
527 },
528 |nodal_velocities: &NodalVelocitiesBlock| {
529 Ok(self.nodal_stiffnesses(nodal_coordinates, nodal_velocities)?)
530 },
531 NodalVelocitiesBlock::zero(num_coords),
532 equality_constraint,
533 )
534 }
535}
536
537impl<C, F, const G: usize, const N: usize> ElasticHyperviscousFiniteElementBlock<C, F, G, N>
538 for ElementBlock<F, N>
539where
540 C: ElasticHyperviscous,
541 F: ElasticHyperviscousFiniteElement<C, G, N>,
542 Self: ViscoelasticFiniteElementBlock<C, F, G, N>,
543{
544 fn viscous_dissipation(
545 &self,
546 nodal_coordinates: &NodalCoordinatesBlock,
547 nodal_velocities: &NodalVelocitiesBlock,
548 ) -> Result<Scalar, ConstitutiveError> {
549 self.elements()
550 .iter()
551 .zip(self.connectivity().iter())
552 .map(|(element, element_connectivity)| {
553 element.viscous_dissipation(
554 &self.nodal_coordinates_element(element_connectivity, nodal_coordinates),
555 &self.nodal_velocities_element(element_connectivity, nodal_velocities),
556 )
557 })
558 .sum()
559 }
560 fn dissipation_potential(
561 &self,
562 nodal_coordinates: &NodalCoordinatesBlock,
563 nodal_velocities: &NodalVelocitiesBlock,
564 ) -> Result<Scalar, ConstitutiveError> {
565 self.elements()
566 .iter()
567 .zip(self.connectivity().iter())
568 .map(|(element, element_connectivity)| {
569 element.dissipation_potential(
570 &self.nodal_coordinates_element(element_connectivity, nodal_coordinates),
571 &self.nodal_velocities_element(element_connectivity, nodal_velocities),
572 )
573 })
574 .sum()
575 }
576 fn minimize(
577 &self,
578 equality_constraint: EqualityConstraint,
579 integrator: impl Explicit<NodalVelocitiesBlock, NodalVelocitiesHistory>,
580 time: &[Scalar],
581 ) -> Result<(Times, NodalCoordinatesHistory, NodalVelocitiesHistory), IntegrationError> {
582 integrator.integrate(
583 |_: Scalar, nodal_coordinates: &NodalCoordinatesBlock| {
584 Ok(self.minimize_inner(equality_constraint.clone(), nodal_coordinates)?)
585 },
586 time,
587 self.coordinates().clone().into(),
588 )
589 }
590 #[doc(hidden)]
591 fn minimize_inner(
592 &self,
593 equality_constraint: EqualityConstraint,
594 nodal_coordinates: &NodalCoordinatesBlock,
595 ) -> Result<NodalVelocitiesBlock, OptimizeError> {
596 let num_coords = nodal_coordinates.len();
597 let banded = band(self.connectivity(), &equality_constraint, num_coords);
598 let solver = NewtonRaphson {
599 ..Default::default()
600 };
601 solver.minimize(
602 |nodal_velocities: &NodalVelocitiesBlock| {
603 Ok(self.dissipation_potential(nodal_coordinates, nodal_velocities)?)
604 },
605 |nodal_velocities: &NodalVelocitiesBlock| {
606 Ok(self.nodal_forces(nodal_coordinates, nodal_velocities)?)
607 },
608 |nodal_velocities: &NodalVelocitiesBlock| {
609 Ok(self.nodal_stiffnesses(nodal_coordinates, nodal_velocities)?)
610 },
611 NodalVelocitiesBlock::zero(num_coords),
612 equality_constraint,
613 Some(banded),
614 )
615 }
616}
617
618impl<C, F, const G: usize, const N: usize> HyperviscoelasticFiniteElementBlock<C, F, G, N>
619 for ElementBlock<F, N>
620where
621 C: Hyperviscoelastic,
622 F: HyperviscoelasticFiniteElement<C, G, N>,
623 Self: ElasticHyperviscousFiniteElementBlock<C, F, G, N>,
624{
625 fn helmholtz_free_energy(
626 &self,
627 nodal_coordinates: &NodalCoordinatesBlock,
628 ) -> Result<Scalar, ConstitutiveError> {
629 self.elements()
630 .iter()
631 .zip(self.connectivity().iter())
632 .map(|(element, element_connectivity)| {
633 element.helmholtz_free_energy(
634 &self.nodal_coordinates_element(element_connectivity, nodal_coordinates),
635 )
636 })
637 .sum()
638 }
639}
640
641fn band<const N: usize>(
642 connectivity: &Connectivity<N>,
643 equality_constraint: &EqualityConstraint,
644 number_of_nodes: usize,
645) -> Banded {
646 match equality_constraint {
647 EqualityConstraint::Linear(matrix, _) => {
648 let neighbors: Vec<Vec<usize>> = invert(connectivity, number_of_nodes)
649 .iter()
650 .map(|elements| {
651 let mut nodes: Vec<usize> = elements
652 .iter()
653 .flat_map(|&element| connectivity[element])
654 .collect();
655 nodes.sort();
656 nodes.dedup();
657 nodes
658 })
659 .collect();
660 let structure: Vec<Vec<bool>> = neighbors
661 .iter()
662 .map(|nodes| (0..number_of_nodes).map(|b| nodes.contains(&b)).collect())
663 .collect();
664 let structure_3d: Vec<Vec<bool>> = structure
665 .iter()
666 .flat_map(|row| {
667 repeat_n(
668 row.iter().flat_map(|entry| repeat_n(*entry, 3)).collect(),
669 3,
670 )
671 })
672 .collect();
673 let num_coords = 3 * number_of_nodes;
674 assert_eq!(matrix.width(), num_coords);
675 let num_dof = matrix.len() + matrix.width();
676 let mut banded = vec![vec![false; num_dof]; num_dof];
677 structure_3d
678 .iter()
679 .zip(banded.iter_mut())
680 .for_each(|(structure_3d_i, banded_i)| {
681 structure_3d_i
682 .iter()
683 .zip(banded_i.iter_mut())
684 .for_each(|(structure_3d_ij, banded_ij)| *banded_ij = *structure_3d_ij)
685 });
686 let mut index = num_coords;
687 matrix.iter().for_each(|matrix_i| {
688 matrix_i.iter().enumerate().for_each(|(j, matrix_ij)| {
689 if matrix_ij != &0.0 {
690 banded[index][j] = true;
691 banded[j][index] = true;
692 index += 1;
693 }
694 })
695 });
696 Banded::from(banded)
697 }
698 _ => unimplemented!(),
699 }
700}
701
702fn invert<const N: usize>(
703 connectivity: &Connectivity<N>,
704 number_of_nodes: usize,
705) -> Vec<Vec<usize>> {
706 let mut inverse_connectivity = vec![vec![]; number_of_nodes];
707 connectivity
708 .iter()
709 .enumerate()
710 .for_each(|(element, nodes)| {
711 nodes
712 .iter()
713 .for_each(|&node| inverse_connectivity[node].push(element))
714 });
715 inverse_connectivity
716}