1pub mod dofs;
2pub mod elements;
3pub mod solid;
4mod solid_boundary;
5mod solid_matrix;
6pub mod solver_solid;
7
8use std::{collections::BTreeMap, fmt};
9
10use runmat_analysis_core::{
11 AnalysisModel, BeamSectionModel, BoundaryConditionKind, LoadKind, ShellSectionModel,
12 StructuralElementKind, StructuralModel,
13};
14use runmat_meshing_core::AnalysisMeshArtifact;
15use serde::{Deserialize, Serialize};
16
17use self::elements::solid::SolidMaterial;
18use self::{
19 dofs::{StructuralDofKind, StructuralDofLayout, StructuralNodeDofSet},
20 elements::beam::{
21 global_stiffness_matrix as beam_global_stiffness_matrix, transformation_matrix,
22 BeamElementGeometry, BeamMaterial, BeamSection, BeamTransform12, BEAM_ELEMENT_DOF_COUNT,
23 },
24 elements::shell::{
25 global_stiffness_matrix as shell_global_stiffness_matrix, ShellElementGeometry,
26 ShellMaterial, ShellSection, SHELL_ELEMENT_DOF_COUNT, SHELL_NODE_DOF_COUNT,
27 },
28 solid::{
29 assemble_solid_stiffness_csr_with_materials, solid_topology_from_analysis_mesh,
30 SolidAssemblyError,
31 },
32 solid_boundary::apply_analysis_mesh_structural_regions,
33};
34
35use crate::operator::{CsrMatrix, OperatorSystem};
36use crate::physics::coupling::thermo_mechanical;
37use crate::{
38 FeaElectroThermalContext, FeaPrepCalibrationProfile, FeaPrepContext, FeaThermoMechanicalContext,
39};
40
41#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
42pub struct AssemblySummary {
43 pub dof_count: usize,
44 #[serde(default)]
45 pub structural_node_count: usize,
46 #[serde(default)]
47 pub structural_translational_dof_count: usize,
48 #[serde(default)]
49 pub structural_rotational_dof_count: usize,
50 #[serde(default)]
51 pub structural_rotation_node_count: usize,
52 #[serde(default)]
53 pub structural_moment_load_count: usize,
54 #[serde(default)]
55 pub structural_direct_rotational_moment_load_count: usize,
56 #[serde(default)]
57 pub structural_wrench_lowering: Vec<WrenchLoweringSummary>,
58 #[serde(default)]
59 pub structural_rotational_constraint_count: usize,
60 #[serde(default)]
61 pub structural_beam_element_count: usize,
62 #[serde(default)]
63 pub structural_shell_element_count: usize,
64 #[serde(default)]
65 pub structural_solid_element_count: usize,
66 #[serde(default)]
67 pub structural_solid_recovery: Vec<SolidRecoveryElementSummary>,
68 #[serde(default)]
69 pub structural_dof_layout: StructuralDofLayout,
70 #[serde(default)]
71 pub structural_beam_recovery: Vec<BeamRecoveryElementSummary>,
72 #[serde(default)]
73 pub structural_shell_recovery: Vec<ShellRecoveryElementSummary>,
74 pub constrained_dof_count: usize,
75 pub load_count: usize,
76 pub structural_material: StructuralMaterialSummary,
77 pub prep_assembly: Option<PrepAssemblySummary>,
78 pub prep_operator_topology: Option<PrepOperatorTopologySummary>,
79 pub prep_region_topology: Option<PrepRegionTopologySummary>,
80 pub prep_element_assembly: Option<PrepElementAssemblySummary>,
81 pub prep_element_connectivity: Option<PrepElementConnectivitySummary>,
82 pub prep_graph_assembly: Option<PrepGraphAssemblySummary>,
83 pub prep_recovery_edges: Vec<PrepRecoveryEdgeSummary>,
84 pub prep_calibration: Option<PrepCalibrationSummary>,
85 pub prep_acceptance: Option<PrepAcceptanceSummary>,
86 pub prep_coordinates: Option<PrepCoordinateSummary>,
87 pub thermo_mechanical: Option<ThermoMechanicalAssemblySummary>,
88 pub electro_thermal: Option<ElectroThermalAssemblySummary>,
89 pub operator: OperatorSystem,
90}
91
92#[derive(Debug, Clone, PartialEq, Eq)]
93pub enum LinearAssemblyError {
94 SolidStiffness(SolidAssemblyError),
95 AnalysisMeshRegionMapping(AnalysisMeshRegionMappingError),
96}
97
98impl fmt::Display for LinearAssemblyError {
99 fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
100 match self {
101 LinearAssemblyError::SolidStiffness(err) => {
102 write!(f, "solid stiffness assembly failed: {err:?}")
103 }
104 LinearAssemblyError::AnalysisMeshRegionMapping(err) => {
105 write!(f, "{err}")
106 }
107 }
108 }
109}
110
111impl std::error::Error for LinearAssemblyError {}
112
113#[derive(Debug, Clone, PartialEq, Eq)]
114pub enum AnalysisMeshRegionMappingError {
115 UnmappedLoadRegion {
116 load_id: String,
117 region_id: String,
118 load_kind: &'static str,
119 },
120 UnmappedBoundaryConditionRegion {
121 bc_id: String,
122 region_id: String,
123 boundary_condition_kind: &'static str,
124 },
125}
126
127impl fmt::Display for AnalysisMeshRegionMappingError {
128 fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
129 match self {
130 AnalysisMeshRegionMappingError::UnmappedLoadRegion {
131 load_id,
132 region_id,
133 load_kind,
134 } => write!(
135 f,
136 "analysis mesh load region did not resolve to solver boundary entities: load_id={load_id} region_id={region_id} load_kind={load_kind}"
137 ),
138 AnalysisMeshRegionMappingError::UnmappedBoundaryConditionRegion {
139 bc_id,
140 region_id,
141 boundary_condition_kind,
142 } => write!(
143 f,
144 "analysis mesh boundary condition region did not resolve to solver boundary entities: bc_id={bc_id} region_id={region_id} boundary_condition_kind={boundary_condition_kind}"
145 ),
146 }
147 }
148}
149
150#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
151pub struct StructuralMaterialSummary {
152 pub youngs_modulus_pa: f64,
153 pub poisson_ratio: f64,
154 #[serde(default)]
155 pub density_kg_per_m3: f64,
156 pub lame_lambda_pa: f64,
157 pub shear_modulus_pa: f64,
158}
159
160#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
161pub struct BeamRecoveryElementSummary {
162 pub element_id: String,
163 pub region_id: String,
164 pub node_i_index: usize,
165 pub node_j_index: usize,
166 pub length_m: f64,
167 pub section: BeamSection,
168 pub material: BeamMaterial,
169 pub transform_global_to_local: BeamTransform12,
170}
171
172#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
173pub struct ShellRecoveryElementSummary {
174 pub element_id: String,
175 pub region_id: String,
176 pub node_indices: [usize; 3],
177 pub area_m2: f64,
178 pub section: ShellSection,
179 pub material: ShellMaterial,
180}
181
182#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
183pub struct SolidRecoveryElementSummary {
184 pub element_id: String,
185 pub region_id: String,
186 pub node_indices: [usize; 4],
187 pub coordinates_m: [[f64; 3]; 4],
188}
189
190#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
191pub struct WrenchLoweringSummary {
192 pub load_id: String,
193 pub region_id: String,
194 pub target_node_count: usize,
195 pub applied_force: [f64; 3],
196 pub applied_moment_at_point: [f64; 3],
197 pub force_residual: [f64; 3],
198 pub moment_residual: [f64; 3],
199 pub moment_couple_applied: bool,
200}
201
202#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
203pub struct PrepAssemblySummary {
204 pub active_region_count: usize,
205 pub mapped_load_count: usize,
206 pub mapped_bc_count: usize,
207 pub mapped_load_ratio: f64,
208 pub constrained_prep_ratio: f64,
209 pub layout_seed: u64,
210}
211
212#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
213pub struct PrepOperatorTopologySummary {
214 pub stiffness_scale: f64,
215 pub mass_scale: f64,
216 pub damping_scale: f64,
217 pub rhs_scale: f64,
218 pub coupling_nonzero_ratio: f64,
219 pub stiffness_spread_ratio: f64,
220 pub topology_fingerprint: u64,
221}
222
223#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
224pub struct PrepRegionTopologySummary {
225 pub region_block_count: usize,
226 pub inter_block_edge_count: usize,
227 pub coupling_nonzero_ratio: f64,
228 pub block_size_min: usize,
229 pub block_size_max: usize,
230 pub block_size_mean: f64,
231 pub region_topology_fingerprint: u64,
232}
233
234#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
235pub struct PrepElementAssemblySummary {
236 pub assembled_element_count: usize,
237 pub triangle_element_count: usize,
238 pub quad_element_count: usize,
239 pub tetrahedron_element_count: usize,
240 pub hex_element_count: usize,
241 pub mixed_element_count: usize,
242 pub scatter_nnz_count: usize,
243 pub assembly_fingerprint: u64,
244}
245
246#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
247pub struct PrepElementConnectivitySummary {
248 pub assembled_element_count: usize,
249 pub stiffness_offdiag_nnz_count: usize,
250 pub mass_offdiag_nnz_count: usize,
251 pub damping_offdiag_nnz_count: usize,
252 pub triangle_contrib_share: f64,
253 pub quad_contrib_share: f64,
254 pub tetrahedron_contrib_share: f64,
255 pub hex_contrib_share: f64,
256 pub mixed_contrib_share: f64,
257 pub mean_connectivity_hop: f64,
258 pub connectivity_fingerprint: u64,
259}
260
261#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
262pub struct PrepGraphAssemblySummary {
263 pub node_count: usize,
264 pub edge_count: usize,
265 pub degree_min: usize,
266 pub degree_max: usize,
267 pub degree_mean: f64,
268 pub degree_p95: f64,
269 pub fill_ratio: f64,
270 pub connected_component_count: usize,
271 pub ordering_bandwidth_before: usize,
272 pub ordering_bandwidth_after: usize,
273 pub ordering_reduction_ratio: f64,
274 pub ordering_fingerprint: u64,
275 pub recommend_ilu0: bool,
276 pub graph_fingerprint: u64,
277}
278
279#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
280pub struct PrepRecoveryEdgeSummary {
281 pub from_dof: usize,
282 pub to_dof: usize,
283 pub element_family_index: usize,
284 pub edge_length_m: f64,
285}
286
287#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
288pub struct PrepCalibrationSummary {
289 pub profile: String,
290 pub triangle_weight: f64,
291 pub quad_weight: f64,
292 pub tetrahedron_weight: f64,
293 pub hex_weight: f64,
294 pub mixed_weight: f64,
295 pub stiffness_calibration_scale: f64,
296 pub mass_calibration_scale: f64,
297 pub damping_calibration_scale: f64,
298 pub calibration_fingerprint: u64,
299}
300
301#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
302pub struct PrepAcceptanceSummary {
303 pub profile: String,
304 pub accepted: bool,
305 pub bounded_displacement_scale: bool,
306 pub bounded_stress_scale: bool,
307 pub bounded_connectivity_fill: bool,
308 pub acceptance_score: f64,
309 pub acceptance_fingerprint: u64,
310}
311
312#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
313pub struct PrepCoordinateSummary {
314 pub span_m: [f64; 3],
315 pub active_dimension_count: usize,
316 pub characteristic_length_m: f64,
317 pub element_geometry_node_count: usize,
318 pub element_geometry_edge_count: usize,
319 pub mean_element_edge_length_m: f64,
320 pub mean_element_area_m2: f64,
321 pub element_geometry_coverage_ratio: f64,
322 pub reference_element_coordinates_m: [[f64; 3]; 3],
323 pub reference_element_area_m2: f64,
324 pub element_topology_sample_element_count: usize,
325 pub element_topology_sample_edge_count: usize,
326 pub element_topology_sample_edge_nodes: [[u32; 2]; 8],
327 pub element_topology_sample_node_coordinates_m: [[f64; 3]; 8],
328 pub element_topology_sample_element_edges: [[u32; 3]; 4],
329 pub element_topology_sample_element_orientations: [[i8; 3]; 4],
330 pub element_topology_sample_element_areas_m2: [f64; 4],
331 pub element_topology_node_coordinates_m: Vec<[f64; 3]>,
332 pub element_topology_edge_nodes: Vec<[u32; 2]>,
333 pub element_topology_element_edges: Vec<[u32; 3]>,
334 pub element_topology_element_orientations: Vec<[i8; 3]>,
335 pub element_topology_element_areas_m2: Vec<f64>,
336}
337
338#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
339pub struct ThermoMechanicalAssemblySummary {
340 pub enabled: bool,
341 pub reference_temperature_k: f64,
342 pub applied_temperature_delta_k: f64,
343 pub thermal_expansion_coefficient: f64,
344 pub thermal_strain_scale: f64,
345 pub thermal_load_scale: f64,
346 pub constitutive_temperature_factor: f64,
347 pub constitutive_poisson_coupling: f64,
348 pub effective_modulus_scale: f64,
349 pub constitutive_material_spread_ratio: f64,
350 pub assignment_heterogeneity_index: f64,
351 pub spatial_gradient_index: f64,
352 pub spatial_coverage_ratio: f64,
353 pub temporal_profile_variation: f64,
354 pub region_delta_count: usize,
355 pub coupling_fingerprint: u64,
356}
357
358#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
359pub struct ElectroThermalAssemblySummary {
360 pub enabled: bool,
361 pub reference_temperature_k: f64,
362 pub applied_voltage_v: f64,
363 pub base_electrical_conductivity_s_per_m: f64,
364 pub resistive_heating_coefficient: f64,
365 pub joule_heating_scale: f64,
366 pub conductivity_spread_ratio: f64,
367 pub temporal_profile_variation: f64,
368 pub region_scale_count: usize,
369 pub coupling_fingerprint: u64,
370 #[serde(default)]
371 pub prep_recovery_edges: Vec<PrepRecoveryEdgeSummary>,
372 #[serde(default)]
373 pub prep_coordinates: Option<PrepCoordinateSummary>,
374}
375
376pub fn assemble_linear_system(
377 model: &AnalysisModel,
378 prep_context: Option<FeaPrepContext>,
379 analysis_mesh: Option<AnalysisMeshArtifact>,
380 thermo_mechanical_context: Option<FeaThermoMechanicalContext>,
381 electro_thermal_context: Option<FeaElectroThermalContext>,
382) -> AssemblySummary {
383 assemble_linear_system_impl(
384 model,
385 prep_context,
386 analysis_mesh,
387 thermo_mechanical_context,
388 electro_thermal_context,
389 false,
390 )
391 .expect("non-strict assembly should build operator topology")
392}
393
394pub fn try_assemble_linear_system(
395 model: &AnalysisModel,
396 prep_context: Option<FeaPrepContext>,
397 analysis_mesh: Option<AnalysisMeshArtifact>,
398 thermo_mechanical_context: Option<FeaThermoMechanicalContext>,
399 electro_thermal_context: Option<FeaElectroThermalContext>,
400) -> Result<AssemblySummary, LinearAssemblyError> {
401 assemble_linear_system_impl(
402 model,
403 prep_context,
404 analysis_mesh,
405 thermo_mechanical_context,
406 electro_thermal_context,
407 true,
408 )
409}
410
411fn assemble_linear_system_impl(
412 model: &AnalysisModel,
413 prep_context: Option<FeaPrepContext>,
414 analysis_mesh: Option<AnalysisMeshArtifact>,
415 thermo_mechanical_context: Option<FeaThermoMechanicalContext>,
416 electro_thermal_context: Option<FeaElectroThermalContext>,
417 strict_analysis_mesh_stiffness: bool,
418) -> Result<AssemblySummary, LinearAssemblyError> {
419 if analysis_mesh.is_none() {
420 if let Some(summary) = assemble_beam_system(model) {
421 return Ok(summary);
422 }
423 }
424
425 let base_dof_count = (model.loads.len() * 3).max(3);
426 let prep_context_ref = prep_context.as_ref();
427 let solid_topology = analysis_mesh
428 .as_ref()
429 .and_then(|mesh| solid_topology_from_analysis_mesh(mesh, base_dof_count).ok());
430 let dof_count = solid_topology
431 .as_ref()
432 .map(|topology| topology.dof_count)
433 .unwrap_or(base_dof_count);
434 let structural_solid_recovery = analysis_mesh
435 .as_ref()
436 .map(solid_recovery_from_analysis_mesh)
437 .unwrap_or_default();
438
439 let avg_youngs_modulus = if model.materials.is_empty() {
440 1.0e9
441 } else {
442 model
443 .materials
444 .iter()
445 .map(|material| material.mechanical.youngs_modulus_pa.max(1.0))
446 .sum::<f64>()
447 / model.materials.len() as f64
448 };
449 let avg_poisson_ratio = if model.materials.is_empty() {
450 0.3
451 } else {
452 model
453 .materials
454 .iter()
455 .map(|material| material.mechanical.poisson_ratio.clamp(0.0, 0.49))
456 .sum::<f64>()
457 / model.materials.len() as f64
458 };
459 let avg_reference_temperature_k = if model.materials.is_empty() {
460 293.15
461 } else {
462 model
463 .materials
464 .iter()
465 .map(|material| material.thermal.reference_temperature_k)
466 .sum::<f64>()
467 / model.materials.len() as f64
468 };
469 let shear_modulus_pa = avg_youngs_modulus / (2.0 * (1.0 + avg_poisson_ratio)).max(1.0e-9);
470 let avg_density_kg_per_m3 = if model.materials.is_empty() {
471 7850.0
472 } else {
473 model
474 .materials
475 .iter()
476 .map(|material| material.mechanical.density_kg_per_m3.max(1.0))
477 .sum::<f64>()
478 / model.materials.len() as f64
479 };
480 let lame_lambda_pa = avg_youngs_modulus * avg_poisson_ratio
481 / ((1.0 + avg_poisson_ratio) * (1.0 - 2.0 * avg_poisson_ratio)).max(1.0e-9);
482 let structural_material = StructuralMaterialSummary {
483 youngs_modulus_pa: avg_youngs_modulus,
484 poisson_ratio: avg_poisson_ratio,
485 density_kg_per_m3: avg_density_kg_per_m3,
486 lame_lambda_pa,
487 shear_modulus_pa,
488 };
489 let stiffness_base = (avg_youngs_modulus / 2.0e3).max(1.0e5);
490
491 let mut stiffness_diag = vec![0.0; dof_count];
492 let mut stiffness_upper = vec![0.0; dof_count.saturating_sub(1)];
493 let mut mass_diag = vec![0.0; dof_count];
494 let mut damping_diag = vec![0.0; dof_count];
495 for i in 0..dof_count {
496 let factor = 1.0 + (i as f64) * 0.05;
497 stiffness_diag[i] = stiffness_base * factor;
498 mass_diag[i] = 1.0 + (i as f64) * 0.01;
499 damping_diag[i] = 0.05 * factor;
500 }
501
502 let mut rhs = vec![0.0; dof_count];
503 let load_base_index = |i: usize, prep: Option<&FeaPrepContext>| -> usize {
504 if let Some(prep) = prep {
505 let stride = (1 + prep.mapped_load_count.max(1)).min(dof_count.max(1));
506 let offset = (prep.layout_seed as usize) % dof_count.max(1);
507 (offset + i.saturating_mul(stride)) % dof_count.max(1)
508 } else {
509 (i * 3) % dof_count
510 }
511 };
512 for (i, load) in model.loads.iter().enumerate() {
513 let base = load_base_index(i, prep_context_ref);
514 match &load.kind {
515 runmat_analysis_core::LoadKind::Force { fx, fy, fz } => {
516 rhs[base] += *fx;
517 if base + 1 < dof_count {
518 rhs[base + 1] += *fy;
519 }
520 if base + 2 < dof_count {
521 rhs[base + 2] += *fz;
522 }
523 }
524 runmat_analysis_core::LoadKind::Moment { .. } => {}
525 runmat_analysis_core::LoadKind::Wrench { fx, fy, fz, .. } => {
526 rhs[base] += *fx;
527 if base + 1 < dof_count {
528 rhs[base + 1] += *fy;
529 }
530 if base + 2 < dof_count {
531 rhs[base + 2] += *fz;
532 }
533 }
534 runmat_analysis_core::LoadKind::Pressure { magnitude_pa } => {
535 rhs[base] += magnitude_pa * 1.0e-3;
536 if base + 1 < dof_count {
537 rhs[base + 1] -= magnitude_pa * 1.0e-3;
538 }
539 }
540 runmat_analysis_core::LoadKind::BodyForce { gx, gy, gz } => {
541 rhs[base] += *gx;
542 if base + 1 < dof_count {
543 rhs[base + 1] += *gy;
544 }
545 if base + 2 < dof_count {
546 rhs[base + 2] += *gz;
547 }
548 }
549 runmat_analysis_core::LoadKind::CurrentDensity { jx, jy, jz, .. } => {
550 rhs[base] += *jx * 1.0e-3;
551 if base + 1 < dof_count {
552 rhs[base + 1] += *jy * 1.0e-3;
553 }
554 if base + 2 < dof_count {
555 rhs[base + 2] += *jz * 1.0e-3;
556 }
557 }
558 runmat_analysis_core::LoadKind::CoilCurrent { current_a, .. } => {
559 rhs[base] += *current_a * 1.0e-2;
560 }
561 runmat_analysis_core::LoadKind::HeatSource { .. } => {}
562 }
563 }
564
565 let structural_dof_layout = StructuralDofLayout::legacy_translational_rows(dof_count);
566 let structural_moment_load_count = model
567 .loads
568 .iter()
569 .filter(|load| matches!(load.kind, runmat_analysis_core::LoadKind::Moment { .. }))
570 .count();
571 let structural_direct_rotational_moment_load_count =
572 if structural_dof_layout.has_rotational_dofs() {
573 structural_moment_load_count
574 } else {
575 0
576 };
577
578 let legacy_constrained_dof_count = model.boundary_conditions.len().min(dof_count);
579 let mut constrained = vec![false; dof_count];
580 let constraint_offset = prep_context_ref
581 .map(|prep| (prep.layout_seed as usize) % dof_count.max(1))
582 .unwrap_or(0);
583 for idx in 0..legacy_constrained_dof_count {
584 let dof = (constraint_offset + idx * 2) % dof_count.max(1);
585 constrained[dof] = true;
586 rhs[dof] = 0.0;
587 }
588 let mut structural_wrench_lowering = Vec::new();
589 if let (Some(mesh), Some(_)) = (analysis_mesh.as_ref(), solid_topology.as_ref()) {
590 structural_wrench_lowering = apply_analysis_mesh_structural_regions(
591 model,
592 mesh,
593 &structural_dof_layout,
594 &mut constrained,
595 &mut rhs,
596 strict_analysis_mesh_stiffness,
597 )
598 .map_err(LinearAssemblyError::AnalysisMeshRegionMapping)?;
599 }
600
601 let mut prep_load_bonus = 0usize;
602 let mut prep_assembly = None;
603 let mut prep_operator_topology = None;
604 let mut prep_region_topology = None;
605 let mut prep_element_assembly = None;
606 let mut prep_element_connectivity = None;
607 let mut prep_graph_assembly = None;
608 let mut prep_recovery_edges = Vec::new();
609 let mut prep_calibration = None;
610 let mut prep_acceptance = None;
611 let mut prep_coordinates = None;
612 let mut thermo_mechanical = None;
613 let mut electro_thermal = None;
614 let mut topology_stiffness_scale = 1.0;
615 let mut topology_mass_scale = 1.0;
616 let mut topology_damping_scale = 1.0;
617 let mut topology_rhs_scale = 1.0;
618 let mut topology_coupling_scale = 1.0;
619 let mut topology_coupling_anisotropy = 1.0;
620 let mut region_block_sizes = Vec::new();
621 let mut region_boundary_positions = Vec::new();
622 let mut region_coupling_weight = 1.0;
623 if let Some(prep) = prep_context_ref {
624 let mesh_scale = 1.0 + (prep.prepared_mesh_count.min(32) as f64) * 0.01;
625 let density = if prep.prepared_node_count == 0 {
626 1.0
627 } else {
628 prep.prepared_element_count as f64 / prep.prepared_node_count as f64
629 };
630 let quality_scale = (prep.min_scaled_jacobian.clamp(0.5, 1.0)
631 / prep.mean_aspect_ratio.clamp(1.0, 4.0))
632 .clamp(0.25, 1.0);
633 let stiffness_scale =
634 (mesh_scale * (1.0 + 0.05 * density.min(2.0)) * quality_scale).clamp(0.5, 1.5);
635 let rhs_scale = (1.0 / quality_scale).clamp(1.0, 2.0);
636 let region_span_norm = (prep.topology_region_span_mean / 12.0).clamp(0.0, 1.0);
637 topology_stiffness_scale = (1.0
638 + 0.12 * prep.topology_volume_core_ratio
639 + 0.06 * prep.topology_surface_patch_ratio
640 + 0.05 * prep.mapped_region_participation_ratio
641 - 0.08 * prep.topology_mixed_family_ratio)
642 .clamp(0.85, 1.25);
643 topology_mass_scale =
644 (1.0 + 0.05 * prep.topology_surface_patch_ratio + 0.04 * region_span_norm
645 - 0.04 * prep.topology_mixed_family_ratio)
646 .clamp(0.9, 1.2);
647 topology_damping_scale = (1.0
648 + 0.03 * prep.topology_volume_core_ratio
649 + 0.05 * prep.topology_mixed_family_ratio
650 + 0.02 * prep.mapped_region_participation_ratio)
651 .clamp(0.9, 1.2);
652 topology_rhs_scale = (1.0
653 + 0.04 * prep.topology_surface_patch_ratio
654 + 0.03 * prep.topology_region_span_mean.clamp(1.0, 24.0) / 24.0)
655 .clamp(1.0, 1.18);
656 topology_coupling_scale = (1.0
657 + 0.12 * prep.topology_volume_core_ratio
658 + 0.08 * prep.mapped_region_participation_ratio
659 - 0.06 * prep.topology_mixed_family_ratio)
660 .clamp(0.8, 1.3);
661 topology_coupling_anisotropy =
662 (1.0 - 0.18 * prep.topology_mixed_family_ratio).clamp(0.78, 1.0);
663 region_coupling_weight = (1.0
664 + 0.10 * prep.mapped_region_participation_ratio
665 + 0.05 * prep.topology_region_mesh_variance.clamp(0.0, 6.0) / 6.0
666 - 0.04 * prep.topology_mixed_family_ratio)
667 .clamp(0.85, 1.2);
668 let region_block_count = prep.topology_region_block_count.clamp(1, dof_count.max(1));
669 region_block_sizes = build_region_block_sizes(
670 dof_count,
671 region_block_count,
672 prep.layout_seed,
673 prep.topology_region_mesh_mean,
674 prep.topology_region_mesh_variance,
675 prep.mapped_region_participation_ratio,
676 );
677 let region_block_offsets = block_offsets(®ion_block_sizes);
678 region_boundary_positions = region_block_offsets
679 .iter()
680 .skip(1)
681 .map(|offset| offset.saturating_sub(1))
682 .filter(|index| *index < dof_count.saturating_sub(1))
683 .collect::<Vec<_>>();
684
685 for (block_index, offset) in region_block_offsets.iter().enumerate() {
686 let block_size = region_block_sizes[block_index];
687 let block_end = offset.saturating_add(block_size).min(dof_count);
688 let block_bias = block_bias(prep.layout_seed, block_index)
689 * (0.04 + 0.02 * prep.topology_region_mesh_mean.clamp(1.0, 6.0) / 6.0);
690 let stiffness_block_scale = (1.0 + block_bias).clamp(0.9, 1.1);
691 let mass_block_scale = (1.0 + 0.6 * block_bias).clamp(0.92, 1.08);
692 let damping_block_scale = (1.0 + 0.8 * block_bias).clamp(0.9, 1.1);
693 for idx in *offset..block_end {
694 stiffness_diag[idx] *= stiffness_block_scale;
695 mass_diag[idx] *= mass_block_scale;
696 damping_diag[idx] *= damping_block_scale;
697 }
698 }
699
700 for value in &mut stiffness_diag {
701 *value *= stiffness_scale * topology_stiffness_scale;
702 }
703 for value in &mut mass_diag {
704 *value *= (1.0 + 0.02 * density.min(2.0)).clamp(1.0, 1.04) * topology_mass_scale;
705 }
706 for value in &mut damping_diag {
707 *value *= (1.0 + 0.03 * prep.mean_aspect_ratio.min(3.0)).clamp(1.0, 1.09)
708 * topology_damping_scale;
709 }
710 for value in &mut rhs {
711 *value *= rhs_scale * topology_rhs_scale;
712 }
713 prep_element_assembly = Some(apply_prep_native_element_assembly(
714 prep,
715 dof_count,
716 &constrained,
717 &mut stiffness_diag,
718 &mut mass_diag,
719 &mut damping_diag,
720 &mut rhs,
721 ));
722 prep_load_bonus = prep
723 .mapped_region_count
724 .saturating_add(prep.inverted_element_count.min(8));
725 prep_coordinates = Some(PrepCoordinateSummary {
726 span_m: [
727 prep.coordinate_span_x_m,
728 prep.coordinate_span_y_m,
729 prep.coordinate_span_z_m,
730 ],
731 active_dimension_count: prep.coordinate_active_dimension_count.max(1),
732 characteristic_length_m: prep.coordinate_characteristic_length_m,
733 element_geometry_node_count: prep.element_geometry_node_count,
734 element_geometry_edge_count: prep.element_geometry_edge_count,
735 mean_element_edge_length_m: prep.mean_element_edge_length_m,
736 mean_element_area_m2: prep.mean_element_area_m2,
737 element_geometry_coverage_ratio: prep.element_geometry_coverage_ratio,
738 reference_element_coordinates_m: prep.reference_element_coordinates_m,
739 reference_element_area_m2: prep.reference_element_area_m2,
740 element_topology_sample_element_count: prep.element_topology_sample_element_count,
741 element_topology_sample_edge_count: prep.element_topology_sample_edge_count,
742 element_topology_sample_edge_nodes: prep.element_topology_sample_edge_nodes,
743 element_topology_sample_node_coordinates_m: prep
744 .element_topology_sample_node_coordinates_m,
745 element_topology_sample_element_edges: prep.element_topology_sample_element_edges,
746 element_topology_sample_element_orientations: prep
747 .element_topology_sample_element_orientations,
748 element_topology_sample_element_areas_m2: prep.element_topology_sample_element_areas_m2,
749 element_topology_node_coordinates_m: prep.element_topology_node_coordinates_m.clone(),
750 element_topology_edge_nodes: prep.element_topology_edge_nodes.clone(),
751 element_topology_element_edges: prep.element_topology_element_edges.clone(),
752 element_topology_element_orientations: prep
753 .element_topology_element_orientations
754 .clone(),
755 element_topology_element_areas_m2: prep.element_topology_element_areas_m2.clone(),
756 });
757
758 prep_assembly = Some(PrepAssemblySummary {
759 active_region_count: prep.mapped_region_count,
760 mapped_load_count: prep.mapped_load_count,
761 mapped_bc_count: prep.mapped_bc_count,
762 mapped_load_ratio: if model.loads.is_empty() {
763 0.0
764 } else {
765 prep.mapped_load_count as f64 / model.loads.len() as f64
766 },
767 constrained_prep_ratio: if model.boundary_conditions.is_empty() {
768 0.0
769 } else {
770 prep.mapped_bc_count as f64 / model.boundary_conditions.len() as f64
771 },
772 layout_seed: prep.layout_seed,
773 });
774 }
775
776 let bandwidth_stride = prep_context_ref
777 .map(|prep| prep.topology_bandwidth_estimate.max(1) as usize)
778 .unwrap_or(1);
779 for i in 0..stiffness_upper.len() {
780 let at_region_boundary = region_boundary_positions.binary_search(&i).is_ok();
781 let local_coupling_scale = if i % 2 == 0 {
782 topology_coupling_scale
783 } else {
784 topology_coupling_scale * topology_coupling_anisotropy
785 };
786 let region_boundary_scale = if at_region_boundary {
787 region_coupling_weight * 0.75
788 } else {
789 region_coupling_weight
790 };
791 let coupling = 0.05
792 * stiffness_diag[i].min(stiffness_diag[i + 1])
793 * local_coupling_scale
794 * region_boundary_scale;
795 let in_sparse_band = bandwidth_stride <= 1 || (i % bandwidth_stride != 0);
796 stiffness_upper[i] = if constrained[i] || constrained[i + 1] || !in_sparse_band {
797 0.0
798 } else {
799 coupling
800 };
801 }
802
803 if let Some(prep) = prep_context_ref {
804 if let Some(element_summary) = prep_element_assembly.as_ref() {
805 let (connectivity_summary, graph_summary, recovery_edges) =
806 apply_prep_element_connectivity_scatter(
807 prep,
808 &constrained,
809 &mut stiffness_upper,
810 &mut mass_diag,
811 &mut damping_diag,
812 element_summary,
813 );
814 prep_element_connectivity = Some(connectivity_summary);
815 prep_graph_assembly = Some(graph_summary);
816 prep_recovery_edges = recovery_edges;
817 }
818 if let Some(calibration) = apply_prep_calibration(
819 prep,
820 avg_youngs_modulus,
821 prep_graph_assembly.as_ref(),
822 &mut stiffness_diag,
823 &mut mass_diag,
824 &mut damping_diag,
825 &mut rhs,
826 ) {
827 let acceptance = evaluate_prep_acceptance(
828 prep,
829 &calibration,
830 prep_graph_assembly.as_ref(),
831 &stiffness_diag,
832 );
833 prep_calibration = Some(calibration);
834 prep_acceptance = Some(acceptance);
835 }
836 let coupling_nonzero_ratio = if stiffness_upper.is_empty() {
837 0.0
838 } else {
839 stiffness_upper
840 .iter()
841 .filter(|value| value.abs() > 0.0)
842 .count() as f64
843 / stiffness_upper.len() as f64
844 };
845 let max_stiffness = stiffness_diag.iter().copied().fold(0.0_f64, f64::max);
846 let min_stiffness = stiffness_diag
847 .iter()
848 .copied()
849 .filter(|value| *value > 0.0)
850 .fold(f64::INFINITY, f64::min);
851 let stiffness_spread_ratio = if min_stiffness.is_finite() && min_stiffness > 0.0 {
852 max_stiffness / min_stiffness
853 } else {
854 0.0
855 };
856 prep_operator_topology = Some(PrepOperatorTopologySummary {
857 stiffness_scale: topology_stiffness_scale,
858 mass_scale: topology_mass_scale,
859 damping_scale: topology_damping_scale,
860 rhs_scale: topology_rhs_scale,
861 coupling_nonzero_ratio,
862 stiffness_spread_ratio,
863 topology_fingerprint: topology_fingerprint(
864 prep,
865 topology_stiffness_scale,
866 topology_mass_scale,
867 topology_damping_scale,
868 topology_rhs_scale,
869 coupling_nonzero_ratio,
870 stiffness_spread_ratio,
871 ),
872 });
873
874 if !region_block_sizes.is_empty() {
875 let block_size_min = region_block_sizes.iter().copied().min().unwrap_or(0);
876 let block_size_max = region_block_sizes.iter().copied().max().unwrap_or(0);
877 let block_size_mean =
878 region_block_sizes.iter().sum::<usize>() as f64 / region_block_sizes.len() as f64;
879 let inter_block_edge_count = region_boundary_positions.len();
880 prep_region_topology = Some(PrepRegionTopologySummary {
881 region_block_count: region_block_sizes.len(),
882 inter_block_edge_count,
883 coupling_nonzero_ratio,
884 block_size_min,
885 block_size_max,
886 block_size_mean,
887 region_topology_fingerprint: region_topology_fingerprint(
888 prep,
889 ®ion_block_sizes,
890 inter_block_edge_count,
891 coupling_nonzero_ratio,
892 ),
893 });
894 }
895 }
896
897 if let Some(context) = thermo_mechanical_context {
898 if context.enabled {
899 let thermal_strain_scale = (context.thermal_expansion_coefficient
900 * context.applied_temperature_delta_k.abs())
901 .clamp(0.0, 0.05);
902 let thermal_load_scale = (context.applied_temperature_delta_k / 50.0).clamp(-2.0, 2.0);
903 let constitutive_temperature_factor = if model.materials.is_empty() {
904 (-(2.5e-4) * context.applied_temperature_delta_k).clamp(-0.25, 0.25)
905 } else {
906 let response = model
907 .materials
908 .iter()
909 .map(|material| {
910 let adjusted_delta = context.applied_temperature_delta_k
911 + (context.reference_temperature_k
912 - material.thermal.reference_temperature_k)
913 + (avg_reference_temperature_k
914 - material.thermal.reference_temperature_k)
915 * 0.1;
916 material.thermal.modulus_temp_coeff_per_k * adjusted_delta
917 })
918 .sum::<f64>()
919 / model.materials.len() as f64;
920 response.clamp(-0.25, 0.25)
921 };
922 let constitutive_poisson_coupling =
923 (0.6 + avg_poisson_ratio.clamp(0.0, 0.49)).clamp(0.6, 1.2);
924 let modulus_temperature_scale = (1.0
925 + constitutive_temperature_factor * constitutive_poisson_coupling)
926 .clamp(0.72, 1.15);
927 let thermal_stiffening_scale = (1.0 + 0.35 * thermal_strain_scale).clamp(1.0, 1.06);
928 let effective_modulus_scale =
929 (modulus_temperature_scale * thermal_stiffening_scale).clamp(0.75, 1.2);
930 let mut dof_adjustments = vec![0.0_f64; dof_count];
931 let assignment_heterogeneity_index = apply_thermo_material_heterogeneity(
932 model,
933 dof_count,
934 constitutive_temperature_factor,
935 context.reference_temperature_k,
936 context.applied_temperature_delta_k,
937 &mut dof_adjustments,
938 );
939 let spatial_field =
940 apply_thermo_spatial_field(&context, dof_count, &mut dof_adjustments);
941 let temporal_profile_variation =
942 thermo_mechanical::temporal_profile_variation(Some(&context));
943 let mut local_modulus_scales = vec![effective_modulus_scale; dof_count];
944 for i in 0..dof_count {
945 let thermal_bias = 1.0 + thermal_strain_scale * (1.0 + (i % 3) as f64 * 0.1);
946 let local_scale =
947 (effective_modulus_scale * (1.0 + dof_adjustments[i])).clamp(0.75, 1.2);
948 local_modulus_scales[i] = local_scale;
949 stiffness_diag[i] *= thermal_bias * local_scale;
950 if !constrained[i] {
951 rhs[i] += thermal_load_scale * (1.0 + (i % 5) as f64 * 0.05);
952 }
953 }
954 for i in 0..stiffness_upper.len() {
955 let edge_scale = 0.5 * (local_modulus_scales[i] + local_modulus_scales[i + 1]);
956 stiffness_upper[i] *= edge_scale;
957 }
958 let min_modulus_scale = local_modulus_scales
959 .iter()
960 .copied()
961 .fold(f64::INFINITY, f64::min);
962 let max_modulus_scale = local_modulus_scales.iter().copied().fold(0.0_f64, f64::max);
963 let constitutive_material_spread_ratio =
964 if min_modulus_scale.is_finite() && min_modulus_scale > 0.0 {
965 max_modulus_scale / min_modulus_scale
966 } else {
967 1.0
968 };
969 thermo_mechanical = Some(ThermoMechanicalAssemblySummary {
970 enabled: true,
971 reference_temperature_k: context.reference_temperature_k,
972 applied_temperature_delta_k: context.applied_temperature_delta_k,
973 thermal_expansion_coefficient: context.thermal_expansion_coefficient,
974 thermal_strain_scale,
975 thermal_load_scale,
976 constitutive_temperature_factor,
977 constitutive_poisson_coupling,
978 effective_modulus_scale,
979 constitutive_material_spread_ratio,
980 assignment_heterogeneity_index,
981 spatial_gradient_index: spatial_field.gradient_index,
982 spatial_coverage_ratio: spatial_field.coverage_ratio,
983 temporal_profile_variation,
984 region_delta_count: context.region_temperature_deltas.len(),
985 coupling_fingerprint: thermo_mechanical_fingerprint(
986 &context,
987 ThermoMechanicalFingerprintInputs {
988 dof_count,
989 constitutive_temperature_factor,
990 constitutive_poisson_coupling,
991 effective_modulus_scale,
992 constitutive_material_spread_ratio,
993 assignment_heterogeneity_index,
994 spatial_gradient_index: spatial_field.gradient_index,
995 temporal_profile_variation,
996 },
997 ),
998 });
999 }
1000 }
1001
1002 if let Some(context) = electro_thermal_context {
1003 if context.enabled {
1004 let temporal_variation = if context.time_profile.len() < 2 {
1005 0.0
1006 } else {
1007 let min_scale = context
1008 .time_profile
1009 .iter()
1010 .map(|point| point.current_scale)
1011 .fold(f64::INFINITY, f64::min);
1012 let max_scale = context
1013 .time_profile
1014 .iter()
1015 .map(|point| point.current_scale)
1016 .fold(-f64::INFINITY, f64::max);
1017 ((max_scale - min_scale).abs() / 2.0).clamp(0.0, 1.0)
1018 };
1019 let mut conductivity_scales = vec![1.0_f64; dof_count];
1020 for (idx, scale) in context.region_conductivity_scales.iter().enumerate() {
1021 let cursor = (idx * 5 + scale.region_id.len()) % dof_count.max(1);
1022 conductivity_scales[cursor] = scale.conductivity_scale.clamp(0.2, 2.5);
1023 }
1024 let min_scale = conductivity_scales
1025 .iter()
1026 .copied()
1027 .fold(f64::INFINITY, f64::min)
1028 .max(1.0e-6);
1029 let max_scale = conductivity_scales.iter().copied().fold(0.0_f64, f64::max);
1030 let conductivity_spread_ratio = (max_scale / min_scale).clamp(1.0, 8.0);
1031 let joule_heating_scale = (context.applied_voltage_v.powi(2)
1032 * context.base_electrical_conductivity_s_per_m.max(1.0e-9)
1033 * context.resistive_heating_coefficient.max(0.0)
1034 / 1.0e6)
1035 .clamp(0.0, 10.0);
1036
1037 for i in 0..dof_count {
1038 let local = conductivity_scales[i];
1039 damping_diag[i] *= (1.0 + 0.02 * local).clamp(1.0, 1.1);
1040 if !constrained[i] {
1041 rhs[i] += joule_heating_scale * local * (1.0 + (i % 7) as f64 * 0.01);
1042 }
1043 }
1044
1045 electro_thermal = Some(ElectroThermalAssemblySummary {
1046 enabled: true,
1047 reference_temperature_k: context.reference_temperature_k,
1048 applied_voltage_v: context.applied_voltage_v,
1049 base_electrical_conductivity_s_per_m: context.base_electrical_conductivity_s_per_m,
1050 resistive_heating_coefficient: context.resistive_heating_coefficient,
1051 joule_heating_scale,
1052 conductivity_spread_ratio,
1053 temporal_profile_variation: temporal_variation,
1054 region_scale_count: context.region_conductivity_scales.len(),
1055 coupling_fingerprint: electro_thermal_fingerprint(
1056 &context,
1057 dof_count,
1058 joule_heating_scale,
1059 conductivity_spread_ratio,
1060 temporal_variation,
1061 ),
1062 prep_recovery_edges: prep_recovery_edges.clone(),
1063 prep_coordinates: prep_coordinates.clone(),
1064 });
1065 }
1066 }
1067
1068 let stiffness_csr = match analysis_mesh.as_ref() {
1069 Some(mesh) => match assemble_solid_stiffness_csr_with_materials(
1070 mesh,
1071 SolidMaterial {
1072 youngs_modulus_pa: structural_material.youngs_modulus_pa,
1073 poisson_ratio: structural_material.poisson_ratio,
1074 },
1075 &solid_materials_by_region(model),
1076 base_dof_count,
1077 ) {
1078 Ok(dense) => Some(dense),
1079 Err(err) if strict_analysis_mesh_stiffness => {
1080 return Err(LinearAssemblyError::SolidStiffness(err));
1081 }
1082 Err(_) => None,
1083 },
1084 None => None,
1085 };
1086 if let Some(csr) = stiffness_csr.as_ref() {
1087 apply_csr_constraints(csr, &constrained, &mut rhs, dof_count);
1088 for (i, diagonal) in stiffness_diag.iter_mut().enumerate().take(dof_count) {
1089 let start = csr.row_offsets[i];
1090 let end = csr.row_offsets[i + 1];
1091 *diagonal = csr.column_indices[start..end]
1092 .iter()
1093 .zip(csr.values[start..end].iter())
1094 .find_map(|(&column, &value)| (column == i).then_some(value.abs()))
1095 .unwrap_or(1.0e-12)
1096 .max(1.0e-12);
1097 }
1098 stiffness_upper.fill(0.0);
1099 }
1100
1101 let constrained_dof_count = constrained.iter().filter(|value| **value).count();
1102
1103 Ok(AssemblySummary {
1104 dof_count,
1105 structural_node_count: structural_dof_layout.node_count(),
1106 structural_translational_dof_count: structural_dof_layout.translational_dof_count(),
1107 structural_rotational_dof_count: structural_dof_layout.rotational_dof_count(),
1108 structural_rotation_node_count: structural_dof_layout.rotation_node_count(),
1109 structural_moment_load_count,
1110 structural_direct_rotational_moment_load_count,
1111 structural_wrench_lowering,
1112 structural_rotational_constraint_count: 0,
1113 structural_beam_element_count: 0,
1114 structural_shell_element_count: 0,
1115 structural_solid_element_count: solid_topology
1116 .as_ref()
1117 .map(|topology| topology.volume_element_count)
1118 .unwrap_or(0),
1119 structural_solid_recovery,
1120 structural_dof_layout,
1121 structural_beam_recovery: Vec::new(),
1122 structural_shell_recovery: Vec::new(),
1123 constrained_dof_count,
1124 load_count: model.loads.len().saturating_add(prep_load_bonus),
1125 structural_material,
1126 prep_assembly,
1127 prep_operator_topology,
1128 prep_region_topology,
1129 prep_element_assembly,
1130 prep_element_connectivity,
1131 prep_graph_assembly,
1132 prep_recovery_edges,
1133 prep_calibration,
1134 prep_acceptance,
1135 prep_coordinates,
1136 thermo_mechanical,
1137 electro_thermal,
1138 operator: OperatorSystem {
1139 dof_count,
1140 constrained,
1141 stiffness_dense: None,
1142 stiffness_csr,
1143 stiffness_diag,
1144 stiffness_upper,
1145 mass_diag,
1146 damping_diag,
1147 rhs,
1148 },
1149 })
1150}
1151
1152fn solid_recovery_from_analysis_mesh(
1153 mesh: &AnalysisMeshArtifact,
1154) -> Vec<SolidRecoveryElementSummary> {
1155 let node_indices = mesh
1156 .nodes
1157 .iter()
1158 .enumerate()
1159 .map(|(index, node)| (node.node_id, index))
1160 .collect::<BTreeMap<_, _>>();
1161
1162 mesh.volume_elements
1163 .iter()
1164 .filter(|element| element.node_ids.len() == 4)
1165 .filter_map(|element| {
1166 let mut indices = [0_usize; 4];
1167 let mut coordinates_m = [[0.0_f64; 3]; 4];
1168 for (local, node_id) in element.node_ids.iter().copied().enumerate() {
1169 let node_index = *node_indices.get(&node_id)?;
1170 indices[local] = node_index;
1171 coordinates_m[local] = mesh.nodes.get(node_index)?.coordinates_m;
1172 }
1173 Some(SolidRecoveryElementSummary {
1174 element_id: element.element_id.clone(),
1175 region_id: element.material_region_id.clone(),
1176 node_indices: indices,
1177 coordinates_m,
1178 })
1179 })
1180 .collect()
1181}
1182
1183fn electro_thermal_fingerprint(
1184 context: &FeaElectroThermalContext,
1185 dof_count: usize,
1186 joule_heating_scale: f64,
1187 conductivity_spread_ratio: f64,
1188 temporal_profile_variation: f64,
1189) -> u64 {
1190 let mut hash = 1469598103934665603_u64;
1191 for value in [
1192 context.reference_temperature_k.to_bits(),
1193 context.applied_voltage_v.to_bits(),
1194 context.base_electrical_conductivity_s_per_m.to_bits(),
1195 context.resistive_heating_coefficient.to_bits(),
1196 joule_heating_scale.to_bits(),
1197 conductivity_spread_ratio.to_bits(),
1198 temporal_profile_variation.to_bits(),
1199 dof_count as u64,
1200 context.region_conductivity_scales.len() as u64,
1201 context.time_profile.len() as u64,
1202 ] {
1203 hash ^= value;
1204 hash = hash.wrapping_mul(1099511628211);
1205 }
1206 hash
1207}
1208
1209fn topology_fingerprint(
1210 prep: &FeaPrepContext,
1211 stiffness_scale: f64,
1212 mass_scale: f64,
1213 damping_scale: f64,
1214 rhs_scale: f64,
1215 coupling_nonzero_ratio: f64,
1216 stiffness_spread_ratio: f64,
1217) -> u64 {
1218 let mut hash = 1469598103934665603_u64;
1219 for value in [
1220 stiffness_scale.to_bits(),
1221 mass_scale.to_bits(),
1222 damping_scale.to_bits(),
1223 rhs_scale.to_bits(),
1224 coupling_nonzero_ratio.to_bits(),
1225 stiffness_spread_ratio.to_bits(),
1226 prep.topology_surface_patch_ratio.to_bits(),
1227 prep.topology_volume_core_ratio.to_bits(),
1228 prep.topology_mixed_family_ratio.to_bits(),
1229 prep.topology_region_span_mean.to_bits(),
1230 prep.mapped_region_participation_ratio.to_bits(),
1231 prep.topology_dof_multiplier.to_bits(),
1232 prep.topology_bandwidth_estimate as u64,
1233 prep.layout_seed,
1234 ] {
1235 hash ^= value;
1236 hash = hash.wrapping_mul(1099511628211_u64);
1237 }
1238 hash
1239}
1240
1241fn assemble_beam_system(model: &AnalysisModel) -> Option<AssemblySummary> {
1242 let structural = model.structural.as_ref()?;
1243 let beam_elements = structural
1244 .elements
1245 .iter()
1246 .filter(|element| matches!(element.kind, StructuralElementKind::Beam(_)))
1247 .collect::<Vec<_>>();
1248 let shell_elements = structural
1249 .elements
1250 .iter()
1251 .filter(|element| matches!(element.kind, StructuralElementKind::Shell(_)))
1252 .collect::<Vec<_>>();
1253 if (beam_elements.is_empty() && shell_elements.is_empty()) || structural.nodes.is_empty() {
1254 return None;
1255 }
1256
1257 let structural_material = structural_material_summary(model);
1258 let beam_material = BeamMaterial {
1259 youngs_modulus_pa: structural_material.youngs_modulus_pa,
1260 shear_modulus_pa: structural_material.shear_modulus_pa,
1261 };
1262 let shell_material = ShellMaterial {
1263 youngs_modulus_pa: structural_material.youngs_modulus_pa,
1264 poisson_ratio: structural_material.poisson_ratio,
1265 shear_modulus_pa: structural_material.shear_modulus_pa,
1266 };
1267 let node_count = structural.nodes.len();
1268 let structural_dof_layout = StructuralDofLayout::from_node_sets(vec![
1269 StructuralNodeDofSet::translational_rotational();
1270 node_count
1271 ]);
1272 let dof_count = structural_dof_layout.total_dof_count();
1273 let mut dense = vec![0.0_f64; dof_count * dof_count];
1274 let mut rhs = vec![0.0_f64; dof_count];
1275 let mut mass_diag = vec![0.0_f64; dof_count];
1276 let mut damping_diag = vec![0.0_f64; dof_count];
1277 let mut structural_beam_recovery = Vec::new();
1278 let mut structural_shell_recovery = Vec::new();
1279
1280 for element in &beam_elements {
1281 let StructuralElementKind::Beam(beam) = &element.kind else {
1282 continue;
1283 };
1284 let node_i_index = structural_node_index(structural, beam.node_ids[0])?;
1285 let node_j_index = structural_node_index(structural, beam.node_ids[1])?;
1286 let section = structural_beam_section(structural, &beam.section_id)?;
1287 let geometry = BeamElementGeometry {
1288 node_i_m: structural.nodes[node_i_index].coordinates_m,
1289 node_j_m: structural.nodes[node_j_index].coordinates_m,
1290 reference_axis: beam.reference_axis,
1291 };
1292 let frame = geometry.local_frame().ok()?;
1293 let transform_global_to_local = transformation_matrix(frame);
1294 let stiffness = beam_global_stiffness_matrix(section, beam_material, geometry).ok()?;
1295 let element_dofs =
1296 beam_element_dof_indices(&structural_dof_layout, node_i_index, node_j_index)?;
1297 for (local_row, &global_row) in element_dofs.iter().enumerate() {
1298 for (local_col, &global_col) in element_dofs.iter().enumerate() {
1299 dense[global_row * dof_count + global_col] += stiffness[local_row][local_col];
1300 }
1301 }
1302 structural_beam_recovery.push(BeamRecoveryElementSummary {
1303 element_id: element.element_id.clone(),
1304 region_id: element.region_id.clone(),
1305 node_i_index,
1306 node_j_index,
1307 length_m: frame.length_m,
1308 section,
1309 material: beam_material,
1310 transform_global_to_local,
1311 });
1312
1313 add_lumped_beam_mass_and_damping(
1314 &mut mass_diag,
1315 &mut damping_diag,
1316 &element_dofs,
1317 section,
1318 structural_material.density_kg_per_m3,
1319 frame.length_m,
1320 &transform_global_to_local,
1321 );
1322 }
1323
1324 for element in &shell_elements {
1325 let StructuralElementKind::Shell(shell) = &element.kind else {
1326 continue;
1327 };
1328 let node_indices = [
1329 structural_node_index(structural, shell.node_ids[0])?,
1330 structural_node_index(structural, shell.node_ids[1])?,
1331 structural_node_index(structural, shell.node_ids[2])?,
1332 ];
1333 let section = structural_shell_section(structural, &shell.section_id)?;
1334 let geometry = ShellElementGeometry {
1335 nodes_m: [
1336 structural.nodes[node_indices[0]].coordinates_m,
1337 structural.nodes[node_indices[1]].coordinates_m,
1338 structural.nodes[node_indices[2]].coordinates_m,
1339 ],
1340 reference_axis: shell.reference_axis,
1341 };
1342 let frame = geometry.local_frame().ok()?;
1343 let stiffness = shell_global_stiffness_matrix(section, shell_material, geometry).ok()?;
1344 let element_dofs = shell_element_dof_indices(&structural_dof_layout, node_indices)?;
1345 for (local_row, &global_row) in element_dofs.iter().enumerate() {
1346 for (local_col, &global_col) in element_dofs.iter().enumerate() {
1347 dense[global_row * dof_count + global_col] += stiffness[local_row][local_col];
1348 }
1349 }
1350 add_lumped_shell_mass_and_damping(
1351 &mut mass_diag,
1352 &mut damping_diag,
1353 &element_dofs,
1354 section,
1355 structural_material.density_kg_per_m3,
1356 frame.area_m2,
1357 );
1358 structural_shell_recovery.push(ShellRecoveryElementSummary {
1359 element_id: element.element_id.clone(),
1360 region_id: element.region_id.clone(),
1361 node_indices,
1362 area_m2: frame.area_m2,
1363 section,
1364 material: shell_material,
1365 });
1366 }
1367
1368 let mut direct_rotational_moment_load_count = 0usize;
1369 let mut structural_wrench_lowering = Vec::new();
1370 for load in &model.loads {
1371 let target_nodes = structural_target_nodes(structural, &load.region_id);
1372 if target_nodes.is_empty() {
1373 continue;
1374 }
1375 let scale = 1.0 / target_nodes.len() as f64;
1376 match load.kind {
1377 LoadKind::Force { fx, fy, fz } => {
1378 for node_index in target_nodes {
1379 add_rhs(
1380 &structural_dof_layout,
1381 &mut rhs,
1382 node_index,
1383 StructuralDofKind::Ux,
1384 fx * scale,
1385 );
1386 add_rhs(
1387 &structural_dof_layout,
1388 &mut rhs,
1389 node_index,
1390 StructuralDofKind::Uy,
1391 fy * scale,
1392 );
1393 add_rhs(
1394 &structural_dof_layout,
1395 &mut rhs,
1396 node_index,
1397 StructuralDofKind::Uz,
1398 fz * scale,
1399 );
1400 }
1401 }
1402 LoadKind::Moment { mx, my, mz } => {
1403 direct_rotational_moment_load_count += 1;
1404 for node_index in target_nodes {
1405 add_rhs(
1406 &structural_dof_layout,
1407 &mut rhs,
1408 node_index,
1409 StructuralDofKind::Rx,
1410 mx * scale,
1411 );
1412 add_rhs(
1413 &structural_dof_layout,
1414 &mut rhs,
1415 node_index,
1416 StructuralDofKind::Ry,
1417 my * scale,
1418 );
1419 add_rhs(
1420 &structural_dof_layout,
1421 &mut rhs,
1422 node_index,
1423 StructuralDofKind::Rz,
1424 mz * scale,
1425 );
1426 }
1427 }
1428 LoadKind::Wrench {
1429 fx,
1430 fy,
1431 fz,
1432 mx,
1433 my,
1434 mz,
1435 px,
1436 py,
1437 pz,
1438 } => {
1439 let summary = add_wrench_rhs(
1440 structural,
1441 &structural_dof_layout,
1442 &mut rhs,
1443 &target_nodes,
1444 [fx, fy, fz],
1445 [mx, my, mz],
1446 [px, py, pz],
1447 );
1448 structural_wrench_lowering.push(WrenchLoweringSummary {
1449 load_id: load.load_id.clone(),
1450 region_id: load.region_id.clone(),
1451 ..summary
1452 });
1453 }
1454 _ => {}
1455 }
1456 }
1457
1458 let mut constrained = vec![false; dof_count];
1459 for bc in &model.boundary_conditions {
1460 let target_nodes = structural_target_nodes(structural, &bc.region_id);
1461 for node_index in target_nodes {
1462 match bc.kind {
1463 BoundaryConditionKind::Fixed => {
1464 for kind in StructuralDofKind::ORDER {
1465 constrain_dof(
1466 &structural_dof_layout,
1467 &mut constrained,
1468 &mut rhs,
1469 node_index,
1470 kind,
1471 0.0,
1472 );
1473 }
1474 }
1475 BoundaryConditionKind::PrescribedDisplacement => {
1476 for kind in [
1477 StructuralDofKind::Ux,
1478 StructuralDofKind::Uy,
1479 StructuralDofKind::Uz,
1480 ] {
1481 constrain_dof(
1482 &structural_dof_layout,
1483 &mut constrained,
1484 &mut rhs,
1485 node_index,
1486 kind,
1487 0.0,
1488 );
1489 }
1490 }
1491 BoundaryConditionKind::PrescribedRotation { rx, ry, rz } => {
1492 constrain_dof(
1493 &structural_dof_layout,
1494 &mut constrained,
1495 &mut rhs,
1496 node_index,
1497 StructuralDofKind::Rx,
1498 rx,
1499 );
1500 constrain_dof(
1501 &structural_dof_layout,
1502 &mut constrained,
1503 &mut rhs,
1504 node_index,
1505 StructuralDofKind::Ry,
1506 ry,
1507 );
1508 constrain_dof(
1509 &structural_dof_layout,
1510 &mut constrained,
1511 &mut rhs,
1512 node_index,
1513 StructuralDofKind::Rz,
1514 rz,
1515 );
1516 }
1517 _ => {}
1518 }
1519 }
1520 }
1521 apply_dense_constraints(&dense, &constrained, &mut rhs, dof_count);
1522
1523 let stiffness_diag = (0..dof_count)
1524 .map(|idx| {
1525 if constrained[idx] {
1526 1.0
1527 } else {
1528 dense[idx * dof_count + idx].abs().max(1.0e-12)
1529 }
1530 })
1531 .collect::<Vec<_>>();
1532 let stiffness_upper = vec![0.0; dof_count.saturating_sub(1)];
1533 let constrained_dof_count = constrained.iter().filter(|value| **value).count();
1534 let rotational_constraint_count = constrained
1535 .iter()
1536 .enumerate()
1537 .filter(|(_, is_constrained)| **is_constrained)
1538 .filter(|(row, _)| {
1539 structural_dof_layout
1540 .address(*row)
1541 .is_some_and(|address| address.kind.is_rotational())
1542 })
1543 .count();
1544
1545 Some(AssemblySummary {
1546 dof_count,
1547 structural_node_count: structural_dof_layout.node_count(),
1548 structural_translational_dof_count: structural_dof_layout.translational_dof_count(),
1549 structural_rotational_dof_count: structural_dof_layout.rotational_dof_count(),
1550 structural_rotation_node_count: structural_dof_layout.rotation_node_count(),
1551 structural_moment_load_count: model
1552 .loads
1553 .iter()
1554 .filter(|load| matches!(load.kind, LoadKind::Moment { .. }))
1555 .count(),
1556 structural_direct_rotational_moment_load_count: direct_rotational_moment_load_count,
1557 structural_wrench_lowering,
1558 structural_rotational_constraint_count: rotational_constraint_count,
1559 structural_beam_element_count: beam_elements.len(),
1560 structural_shell_element_count: shell_elements.len(),
1561 structural_solid_element_count: 0,
1562 structural_solid_recovery: Vec::new(),
1563 structural_dof_layout,
1564 structural_beam_recovery,
1565 structural_shell_recovery,
1566 constrained_dof_count,
1567 load_count: model.loads.len(),
1568 structural_material,
1569 prep_assembly: None,
1570 prep_operator_topology: None,
1571 prep_region_topology: None,
1572 prep_element_assembly: None,
1573 prep_element_connectivity: None,
1574 prep_graph_assembly: None,
1575 prep_recovery_edges: Vec::new(),
1576 prep_calibration: None,
1577 prep_acceptance: None,
1578 prep_coordinates: None,
1579 thermo_mechanical: None,
1580 electro_thermal: None,
1581 operator: OperatorSystem {
1582 dof_count,
1583 constrained,
1584 stiffness_dense: Some(dense),
1585 stiffness_csr: None,
1586 stiffness_diag,
1587 stiffness_upper,
1588 mass_diag,
1589 damping_diag,
1590 rhs,
1591 },
1592 })
1593}
1594
1595fn structural_material_summary(model: &AnalysisModel) -> StructuralMaterialSummary {
1596 let avg_youngs_modulus = if model.materials.is_empty() {
1597 1.0e9
1598 } else {
1599 model
1600 .materials
1601 .iter()
1602 .map(|material| material.mechanical.youngs_modulus_pa.max(1.0))
1603 .sum::<f64>()
1604 / model.materials.len() as f64
1605 };
1606 let avg_poisson_ratio = if model.materials.is_empty() {
1607 0.3
1608 } else {
1609 model
1610 .materials
1611 .iter()
1612 .map(|material| material.mechanical.poisson_ratio.clamp(0.0, 0.49))
1613 .sum::<f64>()
1614 / model.materials.len() as f64
1615 };
1616 let shear_modulus_pa = avg_youngs_modulus / (2.0 * (1.0 + avg_poisson_ratio)).max(1.0e-9);
1617 let avg_density_kg_per_m3 = if model.materials.is_empty() {
1618 7850.0
1619 } else {
1620 model
1621 .materials
1622 .iter()
1623 .map(|material| material.mechanical.density_kg_per_m3.max(1.0))
1624 .sum::<f64>()
1625 / model.materials.len() as f64
1626 };
1627 let lame_lambda_pa = avg_youngs_modulus * avg_poisson_ratio
1628 / ((1.0 + avg_poisson_ratio) * (1.0 - 2.0 * avg_poisson_ratio)).max(1.0e-9);
1629 StructuralMaterialSummary {
1630 youngs_modulus_pa: avg_youngs_modulus,
1631 poisson_ratio: avg_poisson_ratio,
1632 density_kg_per_m3: avg_density_kg_per_m3,
1633 lame_lambda_pa,
1634 shear_modulus_pa,
1635 }
1636}
1637
1638fn solid_materials_by_region(model: &AnalysisModel) -> BTreeMap<String, SolidMaterial> {
1639 let materials_by_id = model
1640 .materials
1641 .iter()
1642 .map(|material| (material.material_id.as_str(), material))
1643 .collect::<BTreeMap<_, _>>();
1644
1645 model
1646 .material_assignments
1647 .iter()
1648 .filter_map(|assignment| {
1649 let material = materials_by_id.get(assignment.assigned_material_id.as_str())?;
1650 Some((
1651 assignment.region_id.clone(),
1652 SolidMaterial {
1653 youngs_modulus_pa: material.mechanical.youngs_modulus_pa.max(1.0),
1654 poisson_ratio: material.mechanical.poisson_ratio.clamp(0.0, 0.49),
1655 },
1656 ))
1657 })
1658 .collect()
1659}
1660
1661fn structural_node_index(structural: &StructuralModel, node_id: u32) -> Option<usize> {
1662 structural
1663 .nodes
1664 .iter()
1665 .position(|node| node.node_id == node_id)
1666}
1667
1668fn structural_beam_section(structural: &StructuralModel, section_id: &str) -> Option<BeamSection> {
1669 structural
1670 .beam_sections
1671 .iter()
1672 .find(|section| section.section_id == section_id)
1673 .map(section_from_model)
1674}
1675
1676fn structural_shell_section(
1677 structural: &StructuralModel,
1678 section_id: &str,
1679) -> Option<ShellSection> {
1680 structural
1681 .shell_sections
1682 .iter()
1683 .find(|section| section.section_id == section_id)
1684 .map(shell_section_from_model)
1685}
1686
1687fn section_from_model(section: &BeamSectionModel) -> BeamSection {
1688 BeamSection {
1689 area_m2: section.area_m2,
1690 iy_m4: section.iy_m4,
1691 iz_m4: section.iz_m4,
1692 torsion_j_m4: section.torsion_j_m4,
1693 outer_fiber_y_m: section.outer_fiber_y_m,
1694 outer_fiber_z_m: section.outer_fiber_z_m,
1695 torsion_outer_radius_m: section.torsion_outer_radius_m,
1696 }
1697}
1698
1699fn shell_section_from_model(section: &ShellSectionModel) -> ShellSection {
1700 ShellSection {
1701 thickness_m: section.thickness_m,
1702 shear_correction: section.shear_correction,
1703 drilling_stiffness_scale: section.drilling_stiffness_scale,
1704 }
1705}
1706
1707fn beam_element_dof_indices(
1708 layout: &StructuralDofLayout,
1709 node_i_index: usize,
1710 node_j_index: usize,
1711) -> Option<[usize; BEAM_ELEMENT_DOF_COUNT]> {
1712 let mut indices = [0usize; BEAM_ELEMENT_DOF_COUNT];
1713 for (local, kind) in StructuralDofKind::ORDER.iter().copied().enumerate() {
1714 indices[local] = layout.index(node_i_index, kind)?;
1715 indices[local + 6] = layout.index(node_j_index, kind)?;
1716 }
1717 Some(indices)
1718}
1719
1720fn shell_element_dof_indices(
1721 layout: &StructuralDofLayout,
1722 node_indices: [usize; 3],
1723) -> Option<[usize; SHELL_ELEMENT_DOF_COUNT]> {
1724 let mut indices = [0usize; SHELL_ELEMENT_DOF_COUNT];
1725 for (node_offset, node_index) in node_indices.iter().enumerate() {
1726 for (component, kind) in StructuralDofKind::ORDER.iter().enumerate() {
1727 indices[node_offset * SHELL_NODE_DOF_COUNT + component] =
1728 layout.index(*node_index, *kind)?;
1729 }
1730 }
1731 Some(indices)
1732}
1733
1734fn add_lumped_beam_mass_and_damping(
1735 mass_diag: &mut [f64],
1736 damping_diag: &mut [f64],
1737 element_dofs: &[usize; BEAM_ELEMENT_DOF_COUNT],
1738 section: BeamSection,
1739 density_kg_per_m3: f64,
1740 length_m: f64,
1741 transform_global_to_local: &BeamTransform12,
1742) {
1743 let density = density_kg_per_m3.max(1.0);
1744 let length = length_m.max(1.0e-12);
1745 let nodal_mass = density * section.area_m2.max(1.0e-18) * length / 2.0;
1746 let local_rotational_inertia = [
1747 density * section.torsion_j_m4.max(1.0e-24) * length / 2.0,
1748 density * section.iy_m4.max(1.0e-24) * length / 2.0,
1749 density * section.iz_m4.max(1.0e-24) * length / 2.0,
1750 ];
1751
1752 for node_offset in [0usize, 6] {
1753 for component in 0..3 {
1754 let dof = element_dofs[node_offset + component];
1755 mass_diag[dof] += nodal_mass;
1756 damping_diag[dof] += 0.01;
1757 }
1758 for global_component in 0..3 {
1759 let local_col = node_offset + 3 + global_component;
1760 let inertia = local_rotational_inertia
1761 .iter()
1762 .enumerate()
1763 .map(|(local_component, local_inertia)| {
1764 let row = node_offset + 3 + local_component;
1765 local_inertia * transform_global_to_local[row][local_col].powi(2)
1766 })
1767 .sum::<f64>();
1768 let dof = element_dofs[node_offset + 3 + global_component];
1769 mass_diag[dof] += inertia.max(1.0e-18);
1770 damping_diag[dof] += 0.01;
1771 }
1772 }
1773}
1774
1775fn add_lumped_shell_mass_and_damping(
1776 mass_diag: &mut [f64],
1777 damping_diag: &mut [f64],
1778 element_dofs: &[usize; SHELL_ELEMENT_DOF_COUNT],
1779 section: ShellSection,
1780 density_kg_per_m3: f64,
1781 area_m2: f64,
1782) {
1783 let density = density_kg_per_m3.max(1.0);
1784 let thickness = section.thickness_m.max(1.0e-12);
1785 let area = area_m2.max(1.0e-18);
1786 let nodal_mass = density * thickness * area / 3.0;
1787 let nodal_rotary_inertia = nodal_mass * thickness.powi(2) / 12.0;
1788 for node_offset in [0usize, 6, 12] {
1789 for component in 0..3 {
1790 let dof = element_dofs[node_offset + component];
1791 mass_diag[dof] += nodal_mass;
1792 damping_diag[dof] += 0.01;
1793 }
1794 for component in 3..6 {
1795 let dof = element_dofs[node_offset + component];
1796 mass_diag[dof] += nodal_rotary_inertia.max(1.0e-18);
1797 damping_diag[dof] += 0.01;
1798 }
1799 }
1800}
1801
1802fn structural_target_nodes(structural: &StructuralModel, region_id: &str) -> Vec<usize> {
1803 if let Some(node_id) = structural_node_selector(region_id) {
1804 return structural_node_index(structural, node_id)
1805 .into_iter()
1806 .collect();
1807 }
1808 let mut nodes = Vec::new();
1809 for element in &structural.elements {
1810 if element.region_id != region_id {
1811 continue;
1812 }
1813 match &element.kind {
1814 StructuralElementKind::Beam(beam) => {
1815 for node_id in beam.node_ids {
1816 if let Some(index) = structural_node_index(structural, node_id) {
1817 if !nodes.contains(&index) {
1818 nodes.push(index);
1819 }
1820 }
1821 }
1822 }
1823 StructuralElementKind::Shell(shell) => {
1824 for node_id in shell.node_ids {
1825 if let Some(index) = structural_node_index(structural, node_id) {
1826 if !nodes.contains(&index) {
1827 nodes.push(index);
1828 }
1829 }
1830 }
1831 }
1832 }
1833 }
1834 nodes
1835}
1836
1837fn structural_node_selector(region_id: &str) -> Option<u32> {
1838 region_id
1839 .strip_prefix("node:")
1840 .unwrap_or(region_id)
1841 .parse::<u32>()
1842 .ok()
1843}
1844
1845fn add_rhs(
1846 layout: &StructuralDofLayout,
1847 rhs: &mut [f64],
1848 node_index: usize,
1849 kind: StructuralDofKind,
1850 value: f64,
1851) {
1852 if let Some(dof) = layout.index(node_index, kind) {
1853 rhs[dof] += value;
1854 }
1855}
1856
1857fn add_wrench_rhs(
1858 structural: &StructuralModel,
1859 layout: &StructuralDofLayout,
1860 rhs: &mut [f64],
1861 target_nodes: &[usize],
1862 force: [f64; 3],
1863 moment_at_point: [f64; 3],
1864 point_m: [f64; 3],
1865) -> WrenchLoweringSummary {
1866 if target_nodes.is_empty() {
1867 return WrenchLoweringSummary {
1868 load_id: String::new(),
1869 region_id: String::new(),
1870 target_node_count: 0,
1871 applied_force: [0.0; 3],
1872 applied_moment_at_point: [0.0; 3],
1873 force_residual: force,
1874 moment_residual: moment_at_point,
1875 moment_couple_applied: false,
1876 };
1877 }
1878
1879 let centroid = target_centroid(structural, target_nodes);
1880 let scale = 1.0 / target_nodes.len() as f64;
1881 let mut nodal_forces = Vec::with_capacity(target_nodes.len());
1882 for &node_index in target_nodes {
1883 let nodal_force = scale_vec(force, scale);
1884 add_translational_rhs(layout, rhs, node_index, nodal_force);
1885 nodal_forces.push(nodal_force);
1886 }
1887
1888 let force_arm = [
1889 centroid[0] - point_m[0],
1890 centroid[1] - point_m[1],
1891 centroid[2] - point_m[2],
1892 ];
1893 let force_moment = cross(force_arm, force);
1894 let couple = [
1895 moment_at_point[0] - force_moment[0],
1896 moment_at_point[1] - force_moment[1],
1897 moment_at_point[2] - force_moment[2],
1898 ];
1899 let mut moment_couple_applied = false;
1900
1901 if !couple
1902 .iter()
1903 .all(|component| component.abs() <= f64::EPSILON)
1904 {
1905 let mut coupling = [[0.0_f64; 3]; 3];
1906 let offsets = target_nodes
1907 .iter()
1908 .map(|&node_index| {
1909 let node = structural.nodes[node_index].coordinates_m;
1910 [
1911 node[0] - centroid[0],
1912 node[1] - centroid[1],
1913 node[2] - centroid[2],
1914 ]
1915 })
1916 .collect::<Vec<_>>();
1917 for offset in &offsets {
1918 let r2 = dot(*offset, *offset);
1919 for row in 0..3 {
1920 coupling[row][row] += r2;
1921 for col in 0..3 {
1922 coupling[row][col] -= offset[row] * offset[col];
1923 }
1924 }
1925 }
1926
1927 if let Some(inv) = invert_3x3(coupling) {
1928 let lambda = mat_vec(inv, couple);
1929 for ((&node_index, offset), nodal_force) in target_nodes
1930 .iter()
1931 .zip(offsets.iter())
1932 .zip(nodal_forces.iter_mut())
1933 {
1934 let couple_force = cross(lambda, *offset);
1935 add_translational_rhs(layout, rhs, node_index, couple_force);
1936 nodal_force[0] += couple_force[0];
1937 nodal_force[1] += couple_force[1];
1938 nodal_force[2] += couple_force[2];
1939 }
1940 moment_couple_applied = true;
1941 }
1942 }
1943
1944 let (applied_force, applied_moment_at_point) =
1945 wrench_resultants(structural, target_nodes, &nodal_forces, point_m);
1946 WrenchLoweringSummary {
1947 load_id: String::new(),
1948 region_id: String::new(),
1949 target_node_count: target_nodes.len(),
1950 applied_force,
1951 applied_moment_at_point,
1952 force_residual: [
1953 force[0] - applied_force[0],
1954 force[1] - applied_force[1],
1955 force[2] - applied_force[2],
1956 ],
1957 moment_residual: [
1958 moment_at_point[0] - applied_moment_at_point[0],
1959 moment_at_point[1] - applied_moment_at_point[1],
1960 moment_at_point[2] - applied_moment_at_point[2],
1961 ],
1962 moment_couple_applied,
1963 }
1964}
1965
1966fn add_translational_rhs(
1967 layout: &StructuralDofLayout,
1968 rhs: &mut [f64],
1969 node_index: usize,
1970 force: [f64; 3],
1971) {
1972 add_rhs(layout, rhs, node_index, StructuralDofKind::Ux, force[0]);
1973 add_rhs(layout, rhs, node_index, StructuralDofKind::Uy, force[1]);
1974 add_rhs(layout, rhs, node_index, StructuralDofKind::Uz, force[2]);
1975}
1976
1977fn target_centroid(structural: &StructuralModel, target_nodes: &[usize]) -> [f64; 3] {
1978 let mut centroid = [0.0_f64; 3];
1979 for &node_index in target_nodes {
1980 let node = structural.nodes[node_index].coordinates_m;
1981 centroid[0] += node[0];
1982 centroid[1] += node[1];
1983 centroid[2] += node[2];
1984 }
1985 scale_vec(centroid, 1.0 / target_nodes.len() as f64)
1986}
1987
1988fn wrench_resultants(
1989 structural: &StructuralModel,
1990 target_nodes: &[usize],
1991 nodal_forces: &[[f64; 3]],
1992 point_m: [f64; 3],
1993) -> ([f64; 3], [f64; 3]) {
1994 let mut applied_force = [0.0_f64; 3];
1995 let mut applied_moment = [0.0_f64; 3];
1996 for (&node_index, &force) in target_nodes.iter().zip(nodal_forces.iter()) {
1997 applied_force[0] += force[0];
1998 applied_force[1] += force[1];
1999 applied_force[2] += force[2];
2000 let node = structural.nodes[node_index].coordinates_m;
2001 let arm = [
2002 node[0] - point_m[0],
2003 node[1] - point_m[1],
2004 node[2] - point_m[2],
2005 ];
2006 let moment = cross(arm, force);
2007 applied_moment[0] += moment[0];
2008 applied_moment[1] += moment[1];
2009 applied_moment[2] += moment[2];
2010 }
2011 (applied_force, applied_moment)
2012}
2013
2014fn dot(a: [f64; 3], b: [f64; 3]) -> f64 {
2015 a[0] * b[0] + a[1] * b[1] + a[2] * b[2]
2016}
2017
2018fn cross(a: [f64; 3], b: [f64; 3]) -> [f64; 3] {
2019 [
2020 a[1] * b[2] - a[2] * b[1],
2021 a[2] * b[0] - a[0] * b[2],
2022 a[0] * b[1] - a[1] * b[0],
2023 ]
2024}
2025
2026fn scale_vec(value: [f64; 3], scale: f64) -> [f64; 3] {
2027 [value[0] * scale, value[1] * scale, value[2] * scale]
2028}
2029
2030fn mat_vec(matrix: [[f64; 3]; 3], value: [f64; 3]) -> [f64; 3] {
2031 [
2032 dot(matrix[0], value),
2033 dot(matrix[1], value),
2034 dot(matrix[2], value),
2035 ]
2036}
2037
2038fn invert_3x3(matrix: [[f64; 3]; 3]) -> Option<[[f64; 3]; 3]> {
2039 let m = matrix;
2040 let c00 = m[1][1] * m[2][2] - m[1][2] * m[2][1];
2041 let c01 = -(m[1][0] * m[2][2] - m[1][2] * m[2][0]);
2042 let c02 = m[1][0] * m[2][1] - m[1][1] * m[2][0];
2043 let c10 = -(m[0][1] * m[2][2] - m[0][2] * m[2][1]);
2044 let c11 = m[0][0] * m[2][2] - m[0][2] * m[2][0];
2045 let c12 = -(m[0][0] * m[2][1] - m[0][1] * m[2][0]);
2046 let c20 = m[0][1] * m[1][2] - m[0][2] * m[1][1];
2047 let c21 = -(m[0][0] * m[1][2] - m[0][2] * m[1][0]);
2048 let c22 = m[0][0] * m[1][1] - m[0][1] * m[1][0];
2049 let det = m[0][0] * c00 + m[0][1] * c01 + m[0][2] * c02;
2050 if det.abs() <= 1.0e-18 {
2051 return None;
2052 }
2053 let inv_det = 1.0 / det;
2054 Some([
2055 [c00 * inv_det, c10 * inv_det, c20 * inv_det],
2056 [c01 * inv_det, c11 * inv_det, c21 * inv_det],
2057 [c02 * inv_det, c12 * inv_det, c22 * inv_det],
2058 ])
2059}
2060
2061fn constrain_dof(
2062 layout: &StructuralDofLayout,
2063 constrained: &mut [bool],
2064 rhs: &mut [f64],
2065 node_index: usize,
2066 kind: StructuralDofKind,
2067 value: f64,
2068) {
2069 if let Some(dof) = layout.index(node_index, kind) {
2070 constrained[dof] = true;
2071 rhs[dof] = value;
2072 }
2073}
2074
2075fn apply_dense_constraints(dense: &[f64], constrained: &[bool], rhs: &mut [f64], dof_count: usize) {
2076 for dof in 0..dof_count {
2077 if !constrained[dof] {
2078 continue;
2079 }
2080 for row in 0..dof_count {
2081 if !constrained[row] {
2082 rhs[row] -= dense[row * dof_count + dof] * rhs[dof];
2083 }
2084 }
2085 }
2086}
2087
2088fn apply_csr_constraints(csr: &CsrMatrix, constrained: &[bool], rhs: &mut [f64], dof_count: usize) {
2089 for row in 0..dof_count {
2090 if constrained[row] {
2091 continue;
2092 }
2093 let start = csr.row_offsets[row];
2094 let end = csr.row_offsets[row + 1];
2095 for entry in start..end {
2096 let column = csr.column_indices[entry];
2097 if constrained[column] {
2098 rhs[row] -= csr.values[entry] * rhs[column];
2099 }
2100 }
2101 }
2102}
2103
2104fn apply_prep_native_element_assembly(
2105 prep: &FeaPrepContext,
2106 dof_count: usize,
2107 constrained: &[bool],
2108 stiffness_diag: &mut [f64],
2109 mass_diag: &mut [f64],
2110 damping_diag: &mut [f64],
2111 rhs: &mut [f64],
2112) -> PrepElementAssemblySummary {
2113 let element_count = prep
2114 .prepared_element_count
2115 .max(prep.prepared_mesh_count)
2116 .max(1);
2117 let triangle_count = ((element_count as f64)
2118 * prep.topology_triangle_family_ratio.clamp(0.0, 1.0))
2119 .round() as usize;
2120 let quad_count =
2121 ((element_count as f64) * prep.topology_quad_family_ratio.clamp(0.0, 1.0)).round() as usize;
2122 let tetrahedron_count = ((element_count as f64)
2123 * prep.topology_tetrahedron_family_ratio.clamp(0.0, 1.0))
2124 .round() as usize;
2125 let mut hex_count =
2126 ((element_count as f64) * prep.topology_hex_family_ratio.clamp(0.0, 1.0)).round() as usize;
2127 let assigned = triangle_count + quad_count + tetrahedron_count + hex_count;
2128 if assigned > element_count {
2129 let overflow = assigned - element_count;
2130 hex_count = hex_count.saturating_sub(overflow);
2131 }
2132 let mixed_count =
2133 element_count.saturating_sub(triangle_count + quad_count + tetrahedron_count + hex_count);
2134
2135 let mut touched_diag = vec![false; dof_count];
2136 let mut touched_rhs = vec![false; dof_count];
2137 let mut element_cursor = 0usize;
2138 let stride = (prep.topology_bandwidth_estimate.max(1) as usize + 1)
2139 .saturating_add(prep.topology_region_block_count.saturating_sub(1));
2140 let span_scale = 1.0 + 0.08 * prep.topology_region_span_mean.clamp(1.0, 24.0) / 24.0;
2141
2142 let mut apply_family = |count: usize, stiffness_factor: f64, mass_factor: f64| {
2143 for _ in 0..count {
2144 let base = ((prep.layout_seed as usize)
2145 .wrapping_add(element_cursor.saturating_mul(stride.max(1))))
2146 % dof_count.max(1);
2147 let wave = 1.0 + ((element_cursor % 17) as f64) / 80.0;
2148 let assembly_scale = stiffness_factor * span_scale * wave;
2149 let mass_scale = mass_factor * (1.0 + 0.04 * prep.topology_surface_patch_ratio);
2150 let damping_scale = (0.9 + 0.2 * prep.mapped_region_participation_ratio)
2151 * (1.0 + 0.05 * prep.topology_mixed_family_ratio);
2152
2153 stiffness_diag[base] += 7.5e4 * assembly_scale;
2154 mass_diag[base] += 0.25 * mass_scale;
2155 damping_diag[base] += 0.012 * damping_scale;
2156 touched_diag[base] = true;
2157
2158 if !constrained[base] {
2159 rhs[base] += 0.5 * assembly_scale;
2160 touched_rhs[base] = true;
2161 }
2162 if base + 1 < dof_count {
2163 stiffness_diag[base + 1] += 2.0e4 * assembly_scale;
2164 mass_diag[base + 1] += 0.08 * mass_scale;
2165 damping_diag[base + 1] += 0.004 * damping_scale;
2166 touched_diag[base + 1] = true;
2167 }
2168 element_cursor = element_cursor.saturating_add(1);
2169 }
2170 };
2171
2172 apply_family(triangle_count, 0.92, 0.95);
2173 apply_family(quad_count, 1.00, 1.00);
2174 apply_family(tetrahedron_count, 1.07, 1.05);
2175 apply_family(hex_count, 1.15, 1.12);
2176 apply_family(mixed_count, 0.98, 1.02);
2177
2178 let scatter_nnz_count = touched_diag.iter().filter(|&&hit| hit).count()
2179 + touched_rhs.iter().filter(|&&hit| hit).count();
2180 PrepElementAssemblySummary {
2181 assembled_element_count: element_count,
2182 triangle_element_count: triangle_count,
2183 quad_element_count: quad_count,
2184 tetrahedron_element_count: tetrahedron_count,
2185 hex_element_count: hex_count,
2186 mixed_element_count: mixed_count,
2187 scatter_nnz_count,
2188 assembly_fingerprint: element_assembly_fingerprint(
2189 prep,
2190 ElementAssemblyFingerprintInputs {
2191 element_count,
2192 triangle_count,
2193 quad_count,
2194 tetrahedron_count,
2195 hex_count,
2196 mixed_count,
2197 scatter_nnz_count,
2198 },
2199 ),
2200 }
2201}
2202
2203fn apply_prep_element_connectivity_scatter(
2204 prep: &FeaPrepContext,
2205 constrained: &[bool],
2206 stiffness_upper: &mut [f64],
2207 mass_diag: &mut [f64],
2208 damping_diag: &mut [f64],
2209 element_summary: &PrepElementAssemblySummary,
2210) -> (
2211 PrepElementConnectivitySummary,
2212 PrepGraphAssemblySummary,
2213 Vec<PrepRecoveryEdgeSummary>,
2214) {
2215 if stiffness_upper.is_empty() {
2216 let connectivity_summary = PrepElementConnectivitySummary {
2217 assembled_element_count: element_summary.assembled_element_count,
2218 stiffness_offdiag_nnz_count: 0,
2219 mass_offdiag_nnz_count: 0,
2220 damping_offdiag_nnz_count: 0,
2221 triangle_contrib_share: 0.0,
2222 quad_contrib_share: 0.0,
2223 tetrahedron_contrib_share: 0.0,
2224 hex_contrib_share: 0.0,
2225 mixed_contrib_share: 0.0,
2226 mean_connectivity_hop: 0.0,
2227 connectivity_fingerprint: element_connectivity_fingerprint(
2228 prep,
2229 ElementConnectivityFingerprintInputs {
2230 element_summary,
2231 stiffness_offdiag_nnz_count: 0,
2232 mass_offdiag_nnz_count: 0,
2233 damping_offdiag_nnz_count: 0,
2234 mean_connectivity_hop: 0.0,
2235 shares: [0.0; 5],
2236 graph_fingerprint: 0,
2237 },
2238 ),
2239 };
2240 let graph_summary = PrepGraphAssemblySummary {
2241 node_count: constrained.len(),
2242 edge_count: 0,
2243 degree_min: 0,
2244 degree_max: 0,
2245 degree_mean: 0.0,
2246 degree_p95: 0.0,
2247 fill_ratio: 0.0,
2248 connected_component_count: constrained.len().max(1),
2249 ordering_bandwidth_before: 0,
2250 ordering_bandwidth_after: 0,
2251 ordering_reduction_ratio: 0.0,
2252 ordering_fingerprint: 0,
2253 recommend_ilu0: false,
2254 graph_fingerprint: graph_fingerprint(prep, 0, constrained.len().max(1), 0.0, 0.0),
2255 };
2256 return (connectivity_summary, graph_summary, Vec::new());
2257 }
2258
2259 let node_count = constrained.len().max(1);
2260 let edges = build_prep_graph_edges(prep, node_count, element_summary);
2261 let recovery_edges = edges
2262 .iter()
2263 .map(|(left, right, family_index)| PrepRecoveryEdgeSummary {
2264 from_dof: *left,
2265 to_dof: *right,
2266 element_family_index: *family_index,
2267 edge_length_m: prep_recovery_edge_length_m(prep, (*left).abs_diff(*right)),
2268 })
2269 .collect::<Vec<_>>();
2270 let (degree_min, degree_max, degree_mean, degree_p95, component_count) =
2271 graph_degree_stats(node_count, &edges);
2272 let max_edges = node_count.saturating_mul(node_count.saturating_sub(1)) / 2;
2273 let fill_ratio = if max_edges == 0 {
2274 0.0
2275 } else {
2276 edges.len() as f64 / max_edges as f64
2277 };
2278
2279 let mut touched_stiffness = vec![false; stiffness_upper.len()];
2280 let mut touched_mass = vec![false; mass_diag.len()];
2281 let mut touched_damping = vec![false; damping_diag.len()];
2282 let mut family_contrib = [0.0_f64; 5];
2283 let mut hops = Vec::new();
2284 let family_stiffness = [0.85_f64, 0.95_f64, 1.05_f64, 1.15_f64, 0.9_f64];
2285 let family_mass = [0.10_f64, 0.11_f64, 0.12_f64, 0.14_f64, 0.105_f64];
2286 let family_damping = [0.004_f64, 0.0045_f64, 0.005_f64, 0.0055_f64, 0.0042_f64];
2287 let region_bias = 1.0 + 0.05 * prep.topology_region_mesh_mean.clamp(1.0, 8.0) / 8.0;
2288
2289 for (edge_cursor, (left, right, family_index)) in edges.iter().copied().enumerate() {
2290 if constrained[left] || constrained[right] {
2291 continue;
2292 }
2293 let hop = right.abs_diff(left).max(1);
2294 hops.push(hop as f64);
2295 let wave = 1.0 + ((edge_cursor % 23) as f64) / 100.0;
2296 let stiffness_add = 0.012
2297 * family_stiffness[family_index]
2298 * region_bias
2299 * wave
2300 * prep.topology_dof_multiplier.clamp(1.0, 4.0)
2301 / hop as f64;
2302 let lo = left.min(right);
2303 let hi = left.max(right);
2304 for band in lo..hi {
2305 if band >= stiffness_upper.len() {
2306 continue;
2307 }
2308 let attenuation = 1.0 / (1.0 + (band - lo) as f64);
2309 stiffness_upper[band] +=
2310 stiffness_add * attenuation * (stiffness_upper[band].abs() + 1.0);
2311 touched_stiffness[band] = true;
2312 }
2313 family_contrib[family_index] += stiffness_add.abs();
2314
2315 let mass_add = family_mass[family_index] * wave;
2316 mass_diag[left] += mass_add;
2317 mass_diag[right] += mass_add * 0.8;
2318 touched_mass[left] = true;
2319 touched_mass[right] = true;
2320
2321 let damping_add =
2322 family_damping[family_index] * (1.0 + 0.5 * prep.topology_mixed_family_ratio);
2323 damping_diag[left] += damping_add;
2324 damping_diag[right] += damping_add * 0.85;
2325 touched_damping[left] = true;
2326 touched_damping[right] = true;
2327 }
2328
2329 let stiffness_offdiag_nnz_count = touched_stiffness.iter().filter(|&&hit| hit).count();
2330 let mass_offdiag_nnz_count = touched_mass.iter().filter(|&&hit| hit).count();
2331 let damping_offdiag_nnz_count = touched_damping.iter().filter(|&&hit| hit).count();
2332 let total_contrib = family_contrib.iter().sum::<f64>().max(1.0e-12);
2333 let shares = [
2334 family_contrib[0] / total_contrib,
2335 family_contrib[1] / total_contrib,
2336 family_contrib[2] / total_contrib,
2337 family_contrib[3] / total_contrib,
2338 family_contrib[4] / total_contrib,
2339 ];
2340 let mean_connectivity_hop = if hops.is_empty() {
2341 0.0
2342 } else {
2343 hops.iter().sum::<f64>() / hops.len() as f64
2344 };
2345
2346 let graph_fingerprint_value =
2347 graph_fingerprint(prep, edges.len(), component_count, degree_mean, degree_p95);
2348 let ordering_permutation = graph_ordering_permutation(node_count, &edges);
2349 let ordering_bandwidth_before = graph_bandwidth(&edges, None);
2350 let ordering_bandwidth_after = graph_bandwidth(&edges, Some(&ordering_permutation));
2351 let ordering_reduction_ratio = if ordering_bandwidth_before == 0 {
2352 0.0
2353 } else {
2354 1.0 - (ordering_bandwidth_after as f64 / ordering_bandwidth_before as f64)
2355 };
2356 let ordering_fingerprint = graph_ordering_fingerprint(prep, &ordering_permutation);
2357 let recommend_ilu0 = degree_p95 >= 4.0 || fill_ratio >= 0.02 || component_count <= 3;
2358 let graph_summary = PrepGraphAssemblySummary {
2359 node_count,
2360 edge_count: edges.len(),
2361 degree_min,
2362 degree_max,
2363 degree_mean,
2364 degree_p95,
2365 fill_ratio,
2366 connected_component_count: component_count,
2367 ordering_bandwidth_before,
2368 ordering_bandwidth_after,
2369 ordering_reduction_ratio,
2370 ordering_fingerprint,
2371 recommend_ilu0,
2372 graph_fingerprint: graph_fingerprint_value,
2373 };
2374
2375 let connectivity_summary = PrepElementConnectivitySummary {
2376 assembled_element_count: element_summary.assembled_element_count,
2377 stiffness_offdiag_nnz_count,
2378 mass_offdiag_nnz_count,
2379 damping_offdiag_nnz_count,
2380 triangle_contrib_share: shares[0],
2381 quad_contrib_share: shares[1],
2382 tetrahedron_contrib_share: shares[2],
2383 hex_contrib_share: shares[3],
2384 mixed_contrib_share: shares[4],
2385 mean_connectivity_hop,
2386 connectivity_fingerprint: element_connectivity_fingerprint(
2387 prep,
2388 ElementConnectivityFingerprintInputs {
2389 element_summary,
2390 stiffness_offdiag_nnz_count,
2391 mass_offdiag_nnz_count,
2392 damping_offdiag_nnz_count,
2393 mean_connectivity_hop,
2394 shares,
2395 graph_fingerprint: graph_fingerprint_value,
2396 },
2397 ),
2398 };
2399
2400 (connectivity_summary, graph_summary, recovery_edges)
2401}
2402
2403fn prep_recovery_edge_length_m(prep: &FeaPrepContext, hop: usize) -> f64 {
2404 let measured_edge_length = if prep.element_geometry_coverage_ratio > 0.0 {
2405 prep.mean_element_edge_length_m
2406 } else {
2407 0.0
2408 };
2409 let characteristic = if measured_edge_length.is_finite() && measured_edge_length > 0.0 {
2410 measured_edge_length
2411 } else {
2412 prep.coordinate_characteristic_length_m
2413 };
2414 let length = characteristic * hop.max(1) as f64;
2415 if length.is_finite() && length > 0.0 {
2416 length
2417 } else {
2418 hop.max(1) as f64
2419 }
2420}
2421
2422fn apply_prep_calibration(
2423 prep: &FeaPrepContext,
2424 avg_youngs_modulus: f64,
2425 graph_summary: Option<&PrepGraphAssemblySummary>,
2426 stiffness_diag: &mut [f64],
2427 mass_diag: &mut [f64],
2428 damping_diag: &mut [f64],
2429 rhs: &mut [f64],
2430) -> Option<PrepCalibrationSummary> {
2431 if stiffness_diag.is_empty() {
2432 return None;
2433 }
2434 let profile = select_calibration_profile(prep, avg_youngs_modulus, graph_summary);
2435 let (profile_gain, profile_name) = match profile {
2436 CalibrationProfile::Fast => (0.92, "fast"),
2437 CalibrationProfile::Balanced => (1.0, "balanced"),
2438 CalibrationProfile::Conservative => (1.08, "conservative"),
2439 };
2440
2441 let triangle_weight = (0.95 + 0.08 * prep.topology_triangle_family_ratio) * profile_gain;
2442 let quad_weight = (1.0 + 0.06 * prep.topology_quad_family_ratio) * profile_gain;
2443 let tetrahedron_weight = (1.04 + 0.10 * prep.topology_tetrahedron_family_ratio) * profile_gain;
2444 let hex_weight = (1.08 + 0.12 * prep.topology_hex_family_ratio) * profile_gain;
2445 let mixed_weight = (0.9 + 0.05 * prep.topology_mixed_family_ratio) * profile_gain;
2446
2447 let stiffness_calibration_scale = (triangle_weight * prep.topology_triangle_family_ratio
2448 + quad_weight * prep.topology_quad_family_ratio
2449 + tetrahedron_weight * prep.topology_tetrahedron_family_ratio
2450 + hex_weight * prep.topology_hex_family_ratio
2451 + mixed_weight * prep.topology_mixed_family_ratio.max(0.01))
2452 .clamp(0.8, 1.3);
2453 let mass_calibration_scale = (0.96
2454 + 0.03 * prep.topology_surface_patch_ratio
2455 + 0.04 * prep.topology_region_mesh_mean.clamp(1.0, 6.0) / 6.0)
2456 .clamp(0.9, 1.2)
2457 * profile_gain;
2458 let damping_calibration_scale = (0.94
2459 + 0.05 * prep.topology_mixed_family_ratio
2460 + 0.03 * prep.mapped_region_participation_ratio)
2461 .clamp(0.9, 1.2)
2462 * profile_gain;
2463
2464 for value in stiffness_diag.iter_mut() {
2465 *value *= stiffness_calibration_scale;
2466 }
2467 for value in mass_diag.iter_mut() {
2468 *value *= mass_calibration_scale;
2469 }
2470 for value in damping_diag.iter_mut() {
2471 *value *= damping_calibration_scale;
2472 }
2473 for value in rhs.iter_mut() {
2474 *value *= (2.0 - stiffness_calibration_scale).clamp(0.8, 1.2);
2475 }
2476
2477 Some(PrepCalibrationSummary {
2478 profile: profile_name.to_string(),
2479 triangle_weight,
2480 quad_weight,
2481 tetrahedron_weight,
2482 hex_weight,
2483 mixed_weight,
2484 stiffness_calibration_scale,
2485 mass_calibration_scale,
2486 damping_calibration_scale,
2487 calibration_fingerprint: calibration_fingerprint(
2488 prep,
2489 profile_name,
2490 stiffness_calibration_scale,
2491 mass_calibration_scale,
2492 damping_calibration_scale,
2493 ),
2494 })
2495}
2496
2497fn evaluate_prep_acceptance(
2498 prep: &FeaPrepContext,
2499 calibration: &PrepCalibrationSummary,
2500 graph_summary: Option<&PrepGraphAssemblySummary>,
2501 stiffness_diag: &[f64],
2502) -> PrepAcceptanceSummary {
2503 let bounded_displacement_scale = (0.8..=1.3).contains(&calibration.stiffness_calibration_scale);
2504 let bounded_stress_scale = (0.9..=1.25).contains(&calibration.damping_calibration_scale);
2505 let bounded_connectivity_fill = graph_summary
2506 .map(|graph| graph.fill_ratio <= 0.25 && graph.connected_component_count <= 64)
2507 .unwrap_or(true);
2508 let stiffness_max = stiffness_diag.iter().copied().fold(0.0_f64, f64::max);
2509 let stiffness_min = stiffness_diag
2510 .iter()
2511 .copied()
2512 .filter(|value| *value > 0.0)
2513 .fold(f64::INFINITY, f64::min);
2514 let spread = if stiffness_min.is_finite() && stiffness_min > 0.0 {
2515 stiffness_max / stiffness_min
2516 } else {
2517 0.0
2518 };
2519 let spread_penalty = (spread / 100.0).clamp(0.0, 1.0);
2520 let mut acceptance_score = 1.0 - spread_penalty;
2521 if !bounded_displacement_scale {
2522 acceptance_score -= 0.2;
2523 }
2524 if !bounded_stress_scale {
2525 acceptance_score -= 0.2;
2526 }
2527 if !bounded_connectivity_fill {
2528 acceptance_score -= 0.3;
2529 }
2530 acceptance_score = acceptance_score.clamp(0.0, 1.0);
2531 let accepted = bounded_displacement_scale
2532 && bounded_stress_scale
2533 && bounded_connectivity_fill
2534 && acceptance_score >= 0.4
2535 && prep.min_scaled_jacobian >= 0.45;
2536
2537 PrepAcceptanceSummary {
2538 profile: calibration.profile.clone(),
2539 accepted,
2540 bounded_displacement_scale,
2541 bounded_stress_scale,
2542 bounded_connectivity_fill,
2543 acceptance_score,
2544 acceptance_fingerprint: acceptance_fingerprint(
2545 prep,
2546 &calibration.profile,
2547 accepted,
2548 acceptance_score,
2549 spread,
2550 ),
2551 }
2552}
2553
2554#[derive(Debug, Clone, Copy, PartialEq, Eq)]
2555enum CalibrationProfile {
2556 Fast,
2557 Balanced,
2558 Conservative,
2559}
2560
2561fn select_calibration_profile(
2562 prep: &FeaPrepContext,
2563 avg_youngs_modulus: f64,
2564 graph_summary: Option<&PrepGraphAssemblySummary>,
2565) -> CalibrationProfile {
2566 if let Some(profile) = prep.calibration_profile_override {
2567 return match profile {
2568 FeaPrepCalibrationProfile::Fast => CalibrationProfile::Fast,
2569 FeaPrepCalibrationProfile::Balanced => CalibrationProfile::Balanced,
2570 FeaPrepCalibrationProfile::Conservative => CalibrationProfile::Conservative,
2571 };
2572 }
2573 let stiffness_regime = avg_youngs_modulus;
2574 let ordering_gain = graph_summary
2575 .map(|graph| graph.ordering_reduction_ratio)
2576 .unwrap_or(0.0);
2577 if prep.min_scaled_jacobian < 0.65 || prep.topology_mixed_family_ratio > 0.25 {
2578 CalibrationProfile::Conservative
2579 } else if stiffness_regime < 5.0e10 || ordering_gain > 0.2 {
2580 CalibrationProfile::Fast
2581 } else {
2582 CalibrationProfile::Balanced
2583 }
2584}
2585
2586fn build_prep_graph_edges(
2587 prep: &FeaPrepContext,
2588 node_count: usize,
2589 element_summary: &PrepElementAssemblySummary,
2590) -> Vec<(usize, usize, usize)> {
2591 use std::collections::BTreeSet;
2592
2593 let family_counts = [
2594 element_summary.triangle_element_count,
2595 element_summary.quad_element_count,
2596 element_summary.tetrahedron_element_count,
2597 element_summary.hex_element_count,
2598 element_summary.mixed_element_count,
2599 ];
2600 let family_valence = [2usize, 3, 4, 5, 3];
2601 let stride = (prep.topology_bandwidth_estimate.max(1) as usize)
2602 .saturating_add(prep.topology_region_block_count.max(1));
2603 let max_hop = prep
2604 .topology_region_span_mean
2605 .round()
2606 .clamp(1.0, node_count as f64) as usize;
2607 let mut edges = BTreeSet::new();
2608 let mut cursor = 0usize;
2609 for family_index in 0..family_counts.len() {
2610 for _ in 0..family_counts[family_index] {
2611 let base = ((prep.layout_seed as usize)
2612 .wrapping_add(family_index.saturating_mul(31))
2613 .wrapping_add(cursor.saturating_mul(stride.max(1))))
2614 % node_count.max(1);
2615 for k in 0..family_valence[family_index] {
2616 let hop = 1
2617 + ((prep.layout_seed as usize)
2618 .wrapping_add(k)
2619 .wrapping_add(cursor)
2620 .wrapping_add(family_index.saturating_mul(7))
2621 % max_hop.max(1));
2622 let target = (base + hop) % node_count.max(1);
2623 if base == target {
2624 continue;
2625 }
2626 let lo = base.min(target);
2627 let hi = base.max(target);
2628 edges.insert((lo, hi, family_index));
2629 }
2630 cursor = cursor.saturating_add(1);
2631 }
2632 }
2633 edges.into_iter().collect()
2634}
2635
2636fn graph_degree_stats(
2637 node_count: usize,
2638 edges: &[(usize, usize, usize)],
2639) -> (usize, usize, f64, f64, usize) {
2640 if node_count == 0 {
2641 return (0, 0, 0.0, 0.0, 0);
2642 }
2643 let mut degree = vec![0usize; node_count];
2644 let mut parent = (0..node_count).collect::<Vec<_>>();
2645
2646 fn find(parent: &mut [usize], x: usize) -> usize {
2647 if parent[x] != x {
2648 let root = find(parent, parent[x]);
2649 parent[x] = root;
2650 }
2651 parent[x]
2652 }
2653 fn union(parent: &mut [usize], a: usize, b: usize) {
2654 let ra = find(parent, a);
2655 let rb = find(parent, b);
2656 if ra != rb {
2657 parent[rb] = ra;
2658 }
2659 }
2660
2661 for (a, b, _) in edges {
2662 degree[*a] = degree[*a].saturating_add(1);
2663 degree[*b] = degree[*b].saturating_add(1);
2664 union(&mut parent, *a, *b);
2665 }
2666 let degree_min = degree.iter().copied().min().unwrap_or(0);
2667 let degree_max = degree.iter().copied().max().unwrap_or(0);
2668 let degree_mean = degree.iter().sum::<usize>() as f64 / degree.len() as f64;
2669 let mut sorted_degree = degree.iter().map(|d| *d as f64).collect::<Vec<_>>();
2670 sorted_degree.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
2671 let degree_p95 = if sorted_degree.is_empty() {
2672 0.0
2673 } else {
2674 let index = ((sorted_degree.len() - 1) as f64 * 0.95).round() as usize;
2675 sorted_degree[index]
2676 };
2677
2678 let mut roots = std::collections::BTreeSet::new();
2679 for idx in 0..node_count {
2680 roots.insert(find(&mut parent, idx));
2681 }
2682 (degree_min, degree_max, degree_mean, degree_p95, roots.len())
2683}
2684
2685fn graph_ordering_permutation(node_count: usize, edges: &[(usize, usize, usize)]) -> Vec<usize> {
2686 if node_count == 0 {
2687 return Vec::new();
2688 }
2689 let mut degree = vec![0usize; node_count];
2690 for (a, b, _) in edges {
2691 degree[*a] = degree[*a].saturating_add(1);
2692 degree[*b] = degree[*b].saturating_add(1);
2693 }
2694 let mut order = (0..node_count).collect::<Vec<_>>();
2695 order.sort_by(|a, b| degree[*a].cmp(°ree[*b]).then_with(|| a.cmp(b)));
2696 let mut permutation = vec![0usize; node_count];
2697 for (new_idx, old_idx) in order.iter().copied().enumerate() {
2698 permutation[old_idx] = new_idx;
2699 }
2700 permutation
2701}
2702
2703fn graph_bandwidth(edges: &[(usize, usize, usize)], permutation: Option<&[usize]>) -> usize {
2704 let mut max_bw = 0usize;
2705 for (a, b, _) in edges {
2706 let lhs = permutation.map(|perm| perm[*a]).unwrap_or(*a);
2707 let rhs = permutation.map(|perm| perm[*b]).unwrap_or(*b);
2708 max_bw = max_bw.max(lhs.abs_diff(rhs));
2709 }
2710 max_bw
2711}
2712
2713fn graph_ordering_fingerprint(prep: &FeaPrepContext, permutation: &[usize]) -> u64 {
2714 let mut hash = 1469598103934665603_u64;
2715 hash ^= prep.layout_seed;
2716 hash = hash.wrapping_mul(1099511628211_u64);
2717 for value in permutation.iter().take(256) {
2718 hash ^= *value as u64;
2719 hash = hash.wrapping_mul(1099511628211_u64);
2720 }
2721 hash
2722}
2723
2724fn build_region_block_sizes(
2725 dof_count: usize,
2726 block_count: usize,
2727 layout_seed: u64,
2728 region_mesh_mean: f64,
2729 region_mesh_variance: f64,
2730 mapped_region_participation_ratio: f64,
2731) -> Vec<usize> {
2732 let mut sizes = vec![dof_count / block_count; block_count];
2733 for size in &mut sizes {
2734 if *size == 0 {
2735 *size = 1;
2736 }
2737 }
2738 let assigned = sizes.iter().sum::<usize>();
2739 let mut remainder = dof_count.saturating_sub(assigned);
2740 let seed_bias = ((layout_seed % 13) as usize).max(1);
2741 let participation_bias =
2742 (mapped_region_participation_ratio.clamp(0.0, 1.0) * 7.0).round() as usize;
2743 let variance_bias = region_mesh_variance.clamp(0.0, 16.0).round() as usize;
2744 let stride = (seed_bias + participation_bias + variance_bias).max(1);
2745 let mut cursor = (region_mesh_mean.round() as usize + seed_bias) % block_count.max(1);
2746 while remainder > 0 {
2747 sizes[cursor % block_count] = sizes[cursor % block_count].saturating_add(1);
2748 cursor = cursor.saturating_add(stride);
2749 remainder -= 1;
2750 }
2751 sizes
2752}
2753
2754fn block_offsets(sizes: &[usize]) -> Vec<usize> {
2755 let mut offsets = Vec::with_capacity(sizes.len());
2756 let mut current = 0usize;
2757 for size in sizes {
2758 offsets.push(current);
2759 current = current.saturating_add(*size);
2760 }
2761 offsets
2762}
2763
2764fn block_bias(layout_seed: u64, block_index: usize) -> f64 {
2765 let mut hash = layout_seed ^ ((block_index as u64).wrapping_mul(0x9E37_79B9_7F4A_7C15));
2766 hash ^= hash >> 33;
2767 hash = hash.wrapping_mul(0xff51afd7ed558ccd);
2768 hash ^= hash >> 33;
2769 let normalized = (hash % 1000) as f64 / 999.0;
2770 normalized * 2.0 - 1.0
2771}
2772
2773fn region_topology_fingerprint(
2774 prep: &FeaPrepContext,
2775 block_sizes: &[usize],
2776 inter_block_edge_count: usize,
2777 coupling_nonzero_ratio: f64,
2778) -> u64 {
2779 let mut hash = 1469598103934665603_u64;
2780 for value in [
2781 prep.layout_seed,
2782 prep.topology_region_block_count as u64,
2783 prep.topology_region_mesh_mean.to_bits(),
2784 prep.topology_region_mesh_variance.to_bits(),
2785 prep.topology_region_span_mean.to_bits(),
2786 prep.mapped_region_participation_ratio.to_bits(),
2787 inter_block_edge_count as u64,
2788 coupling_nonzero_ratio.to_bits(),
2789 ] {
2790 hash ^= value;
2791 hash = hash.wrapping_mul(1099511628211_u64);
2792 }
2793 for size in block_sizes {
2794 hash ^= *size as u64;
2795 hash = hash.wrapping_mul(1099511628211_u64);
2796 }
2797 hash
2798}
2799
2800#[derive(Debug, Clone, Copy)]
2801struct ElementAssemblyFingerprintInputs {
2802 element_count: usize,
2803 triangle_count: usize,
2804 quad_count: usize,
2805 tetrahedron_count: usize,
2806 hex_count: usize,
2807 mixed_count: usize,
2808 scatter_nnz_count: usize,
2809}
2810
2811fn element_assembly_fingerprint(
2812 prep: &FeaPrepContext,
2813 inputs: ElementAssemblyFingerprintInputs,
2814) -> u64 {
2815 let mut hash = 1469598103934665603_u64;
2816 for value in [
2817 prep.layout_seed,
2818 prep.prepared_element_count as u64,
2819 inputs.element_count as u64,
2820 inputs.triangle_count as u64,
2821 inputs.quad_count as u64,
2822 inputs.tetrahedron_count as u64,
2823 inputs.hex_count as u64,
2824 inputs.mixed_count as u64,
2825 inputs.scatter_nnz_count as u64,
2826 prep.topology_triangle_family_ratio.to_bits(),
2827 prep.topology_quad_family_ratio.to_bits(),
2828 prep.topology_tetrahedron_family_ratio.to_bits(),
2829 prep.topology_hex_family_ratio.to_bits(),
2830 prep.topology_mixed_family_ratio.to_bits(),
2831 ] {
2832 hash ^= value;
2833 hash = hash.wrapping_mul(1099511628211_u64);
2834 }
2835 hash
2836}
2837
2838#[derive(Debug, Clone, Copy)]
2839struct ElementConnectivityFingerprintInputs<'a> {
2840 element_summary: &'a PrepElementAssemblySummary,
2841 stiffness_offdiag_nnz_count: usize,
2842 mass_offdiag_nnz_count: usize,
2843 damping_offdiag_nnz_count: usize,
2844 mean_connectivity_hop: f64,
2845 shares: [f64; 5],
2846 graph_fingerprint: u64,
2847}
2848
2849fn element_connectivity_fingerprint(
2850 prep: &FeaPrepContext,
2851 inputs: ElementConnectivityFingerprintInputs<'_>,
2852) -> u64 {
2853 let mut hash = 1469598103934665603_u64;
2854 for value in [
2855 prep.layout_seed,
2856 inputs.element_summary.assembly_fingerprint,
2857 inputs.element_summary.assembled_element_count as u64,
2858 inputs.stiffness_offdiag_nnz_count as u64,
2859 inputs.mass_offdiag_nnz_count as u64,
2860 inputs.damping_offdiag_nnz_count as u64,
2861 inputs.mean_connectivity_hop.to_bits(),
2862 inputs.shares[0].to_bits(),
2863 inputs.shares[1].to_bits(),
2864 inputs.shares[2].to_bits(),
2865 inputs.shares[3].to_bits(),
2866 inputs.shares[4].to_bits(),
2867 prep.topology_bandwidth_estimate as u64,
2868 inputs.graph_fingerprint,
2869 ] {
2870 hash ^= value;
2871 hash = hash.wrapping_mul(1099511628211_u64);
2872 }
2873 hash
2874}
2875
2876fn graph_fingerprint(
2877 prep: &FeaPrepContext,
2878 edge_count: usize,
2879 connected_component_count: usize,
2880 degree_mean: f64,
2881 degree_p95: f64,
2882) -> u64 {
2883 let mut hash = 1469598103934665603_u64;
2884 for value in [
2885 prep.layout_seed,
2886 prep.topology_bandwidth_estimate as u64,
2887 prep.topology_region_block_count as u64,
2888 edge_count as u64,
2889 connected_component_count as u64,
2890 degree_mean.to_bits(),
2891 degree_p95.to_bits(),
2892 prep.topology_region_span_mean.to_bits(),
2893 prep.topology_mixed_family_ratio.to_bits(),
2894 ] {
2895 hash ^= value;
2896 hash = hash.wrapping_mul(1099511628211_u64);
2897 }
2898 hash
2899}
2900
2901fn calibration_fingerprint(
2902 prep: &FeaPrepContext,
2903 profile: &str,
2904 stiffness_scale: f64,
2905 mass_scale: f64,
2906 damping_scale: f64,
2907) -> u64 {
2908 let mut hash = 1469598103934665603_u64;
2909 for byte in profile.as_bytes() {
2910 hash ^= *byte as u64;
2911 hash = hash.wrapping_mul(1099511628211_u64);
2912 }
2913 for value in [
2914 prep.layout_seed,
2915 prep.prepared_element_count as u64,
2916 stiffness_scale.to_bits(),
2917 mass_scale.to_bits(),
2918 damping_scale.to_bits(),
2919 prep.topology_triangle_family_ratio.to_bits(),
2920 prep.topology_quad_family_ratio.to_bits(),
2921 prep.topology_tetrahedron_family_ratio.to_bits(),
2922 prep.topology_hex_family_ratio.to_bits(),
2923 ] {
2924 hash ^= value;
2925 hash = hash.wrapping_mul(1099511628211_u64);
2926 }
2927 hash
2928}
2929
2930fn acceptance_fingerprint(
2931 prep: &FeaPrepContext,
2932 profile: &str,
2933 accepted: bool,
2934 score: f64,
2935 spread: f64,
2936) -> u64 {
2937 let mut hash = calibration_fingerprint(
2938 prep,
2939 profile,
2940 score,
2941 spread,
2942 if accepted { 1.0 } else { 0.0 },
2943 );
2944 hash ^= accepted as u64;
2945 hash = hash.wrapping_mul(1099511628211_u64);
2946 hash
2947}
2948
2949#[derive(Debug, Clone, Copy)]
2950struct ThermoMechanicalFingerprintInputs {
2951 dof_count: usize,
2952 constitutive_temperature_factor: f64,
2953 constitutive_poisson_coupling: f64,
2954 effective_modulus_scale: f64,
2955 constitutive_material_spread_ratio: f64,
2956 assignment_heterogeneity_index: f64,
2957 spatial_gradient_index: f64,
2958 temporal_profile_variation: f64,
2959}
2960
2961fn thermo_mechanical_fingerprint(
2962 context: &FeaThermoMechanicalContext,
2963 inputs: ThermoMechanicalFingerprintInputs,
2964) -> u64 {
2965 let mut hash = 1469598103934665603_u64;
2966 for value in [
2967 inputs.dof_count as u64,
2968 context.reference_temperature_k.to_bits(),
2969 context.applied_temperature_delta_k.to_bits(),
2970 context.thermal_expansion_coefficient.to_bits(),
2971 inputs.constitutive_temperature_factor.to_bits(),
2972 inputs.constitutive_poisson_coupling.to_bits(),
2973 inputs.effective_modulus_scale.to_bits(),
2974 inputs.constitutive_material_spread_ratio.to_bits(),
2975 inputs.assignment_heterogeneity_index.to_bits(),
2976 inputs.spatial_gradient_index.to_bits(),
2977 inputs.temporal_profile_variation.to_bits(),
2978 ] {
2979 hash ^= value;
2980 hash = hash.wrapping_mul(1099511628211_u64);
2981 }
2982 hash
2983}
2984
2985#[derive(Debug, Clone, Copy, PartialEq)]
2986struct ThermoSpatialFieldSummary {
2987 gradient_index: f64,
2988 coverage_ratio: f64,
2989}
2990
2991fn apply_thermo_spatial_field(
2992 context: &FeaThermoMechanicalContext,
2993 dof_count: usize,
2994 dof_adjustments: &mut [f64],
2995) -> ThermoSpatialFieldSummary {
2996 if dof_count == 0 || context.region_temperature_deltas.is_empty() {
2997 return ThermoSpatialFieldSummary {
2998 gradient_index: 0.0,
2999 coverage_ratio: 0.0,
3000 };
3001 }
3002 let mut touched = vec![false; dof_count];
3003 let mut min_delta = f64::INFINITY;
3004 let mut max_delta = -f64::INFINITY;
3005 for (idx, region_delta) in context.region_temperature_deltas.iter().enumerate() {
3006 min_delta = min_delta.min(region_delta.temperature_delta_k);
3007 max_delta = max_delta.max(region_delta.temperature_delta_k);
3008 let normalized = ((region_delta.temperature_delta_k - context.applied_temperature_delta_k)
3009 / 240.0)
3010 .clamp(-0.45, 0.45);
3011 let start =
3012 ((region_hash(®ion_delta.region_id) as usize).wrapping_add(idx * 5)) % dof_count;
3013 let stride = context
3014 .region_temperature_deltas
3015 .len()
3016 .saturating_add(3)
3017 .max(2);
3018 let mut cursor = start;
3019 for hop in 0..dof_count {
3020 if hop > 0 && cursor == start {
3021 break;
3022 }
3023 let wave = 1.0 + ((hop + idx) % 7) as f64 * 0.02;
3024 dof_adjustments[cursor] += normalized * wave;
3025 touched[cursor] = true;
3026 cursor = (cursor + stride) % dof_count;
3027 }
3028 }
3029 if !min_delta.is_finite() || !max_delta.is_finite() {
3030 return ThermoSpatialFieldSummary {
3031 gradient_index: 0.0,
3032 coverage_ratio: 0.0,
3033 };
3034 }
3035 let touched_count = touched.iter().filter(|entry| **entry).count() as f64;
3036 ThermoSpatialFieldSummary {
3037 gradient_index: ((max_delta - min_delta).abs() / 240.0).clamp(0.0, 1.0),
3038 coverage_ratio: (touched_count / dof_count as f64).clamp(0.0, 1.0),
3039 }
3040}
3041
3042fn apply_thermo_material_heterogeneity(
3043 model: &AnalysisModel,
3044 dof_count: usize,
3045 constitutive_temperature_factor: f64,
3046 reference_temperature_k: f64,
3047 applied_temperature_delta_k: f64,
3048 dof_adjustments: &mut [f64],
3049) -> f64 {
3050 if dof_count == 0 || model.material_assignments.is_empty() {
3051 return 0.0;
3052 }
3053 let base_amplitude = (constitutive_temperature_factor.abs() * 0.8).clamp(0.0, 0.15);
3054 if base_amplitude <= 0.0 {
3055 return 0.0;
3056 }
3057 let mut weighted_activity = 0.0_f64;
3058 let mut weight_sum = 0.0_f64;
3059 for (idx, assignment) in model.material_assignments.iter().enumerate() {
3060 let confidence_weight = match assignment.confidence {
3061 runmat_analysis_core::EvidenceConfidence::Verified => 1.0,
3062 runmat_analysis_core::EvidenceConfidence::Probable => 0.65,
3063 runmat_analysis_core::EvidenceConfidence::Inferred => 0.4,
3064 };
3065 let expected_modulus = model
3066 .materials
3067 .iter()
3068 .find(|material| material.material_id == assignment.expected_material_id)
3069 .map(|material| material.mechanical.youngs_modulus_pa)
3070 .unwrap_or(1.0e9)
3071 .max(1.0);
3072 let assigned_modulus = model
3073 .materials
3074 .iter()
3075 .find(|material| material.material_id == assignment.assigned_material_id)
3076 .map(|material| material.mechanical.youngs_modulus_pa)
3077 .unwrap_or(expected_modulus)
3078 .max(1.0);
3079 let modulus_delta_ratio =
3080 ((assigned_modulus - expected_modulus) / expected_modulus).clamp(-0.6, 0.6);
3081 let expected_temp_response = model
3082 .materials
3083 .iter()
3084 .find(|material| material.material_id == assignment.expected_material_id)
3085 .map(|material| {
3086 material.thermal.modulus_temp_coeff_per_k
3087 * (applied_temperature_delta_k
3088 + (reference_temperature_k - material.thermal.reference_temperature_k))
3089 })
3090 .unwrap_or(constitutive_temperature_factor)
3091 .clamp(-0.4, 0.2);
3092 let assigned_temp_response = model
3093 .materials
3094 .iter()
3095 .find(|material| material.material_id == assignment.assigned_material_id)
3096 .map(|material| {
3097 material.thermal.modulus_temp_coeff_per_k
3098 * (applied_temperature_delta_k
3099 + (reference_temperature_k - material.thermal.reference_temperature_k))
3100 })
3101 .unwrap_or(expected_temp_response)
3102 .clamp(-0.4, 0.2);
3103 let response_delta = (assigned_temp_response - expected_temp_response).clamp(-0.35, 0.35);
3104 let region_phase = ((region_hash(&assignment.region_id) % 11) as f64) / 10.0;
3105 let activity =
3106 (0.7 * modulus_delta_ratio.abs() + 0.3 * response_delta.abs()).clamp(0.0, 1.0);
3107 let signed_bias = base_amplitude
3108 * confidence_weight
3109 * (0.45 * modulus_delta_ratio
3110 + 0.35 * response_delta
3111 + 0.2 * modulus_delta_ratio.signum() * region_phase);
3112 let stride = model.material_assignments.len().saturating_add(1).max(2);
3113 let start = ((region_hash(&assignment.region_id) as usize).wrapping_add(idx * 3))
3114 % dof_count.max(1);
3115 let mut cursor = start;
3116 for hop in 0..dof_count {
3117 if hop > 0 && cursor == start {
3118 break;
3119 }
3120 let wave = 1.0 + ((hop + idx) % 5) as f64 * 0.03;
3121 dof_adjustments[cursor] += signed_bias * wave;
3122 cursor = (cursor + stride) % dof_count.max(1);
3123 }
3124 weighted_activity += activity * confidence_weight;
3125 weight_sum += confidence_weight;
3126 }
3127 for value in dof_adjustments.iter_mut() {
3128 *value = value.clamp(-0.18, 0.18);
3129 }
3130 if weight_sum > 0.0 {
3131 (weighted_activity / weight_sum).clamp(0.0, 1.0)
3132 } else {
3133 0.0
3134 }
3135}
3136
3137fn region_hash(region_id: &str) -> u64 {
3138 let mut hash = 1469598103934665603_u64;
3139 for byte in region_id.as_bytes() {
3140 hash ^= *byte as u64;
3141 hash = hash.wrapping_mul(1099511628211_u64);
3142 }
3143 hash
3144}
3145
3146#[cfg(test)]
3147mod tests {
3148 use super::*;
3149 use crate::fixtures::{fixture_model, FixtureId};
3150 use runmat_meshing_core::{
3151 contracts::artifact::ANALYSIS_MESH_SCHEMA_VERSION, AnalysisBoundaryFace,
3152 AnalysisMeshArtifact, AnalysisMeshNode, AnalysisMeshProvenance, AnalysisMeshQualityReport,
3153 AnalysisVolumeElement, BoundaryElementKind, MeshSizingField, VolumeElementKind,
3154 };
3155
3156 #[test]
3157 fn analysis_mesh_populates_sparse_solid_stiffness_operator() {
3158 let model = fixture_model(FixtureId::CantileverLinearStatic);
3159 let summary = assemble_linear_system(&model, None, Some(tetrahedron4_mesh()), None, None);
3160
3161 assert_eq!(summary.dof_count, 12);
3162 assert_eq!(summary.structural_solid_element_count, 1);
3163 assert_eq!(summary.structural_solid_recovery.len(), 1);
3164 assert_eq!(
3165 summary.structural_solid_recovery[0].node_indices,
3166 [0, 1, 2, 3]
3167 );
3168 assert_eq!(
3169 summary.structural_solid_recovery[0].coordinates_m,
3170 [
3171 [0.0, 0.0, 0.0],
3172 [1.0, 0.0, 0.0],
3173 [0.0, 1.0, 0.0],
3174 [0.0, 0.0, 1.0],
3175 ]
3176 );
3177 assert!(summary.operator.stiffness_dense.is_none());
3178 let csr = summary
3179 .operator
3180 .stiffness_csr
3181 .as_ref()
3182 .expect("analysis mesh should assemble a sparse solid stiffness matrix");
3183 assert_eq!(csr.row_offsets.len(), summary.dof_count + 1);
3184 assert_eq!(csr.row_offsets.last().copied(), Some(csr.values.len()));
3185 assert_eq!(csr.column_indices.len(), csr.values.len());
3186 assert!(summary
3187 .operator
3188 .stiffness_diag
3189 .iter()
3190 .all(|value| *value > 0.0));
3191 for row in 0..summary.dof_count {
3192 let start = csr.row_offsets[row];
3193 let end = csr.row_offsets[row + 1];
3194 let diagonal = csr.column_indices[start..end]
3195 .iter()
3196 .zip(csr.values[start..end].iter())
3197 .find_map(|(&column, &value)| (column == row).then_some(value.abs()))
3198 .expect("csr row should include diagonal");
3199 assert!((summary.operator.stiffness_diag[row] - diagonal) <= 1.0e-8);
3200 }
3201 }
3202
3203 #[test]
3204 fn analysis_mesh_preempts_explicit_beam_topology() {
3205 let mut model = fixture_model(FixtureId::CantileverLinearStatic);
3206 model.structural = Some(runmat_analysis_core::StructuralModel {
3207 nodes: vec![
3208 runmat_analysis_core::StructuralNode {
3209 node_id: 1,
3210 coordinates_m: [0.0, 0.0, 0.0],
3211 },
3212 runmat_analysis_core::StructuralNode {
3213 node_id: 2,
3214 coordinates_m: [1.0, 0.0, 0.0],
3215 },
3216 ],
3217 elements: vec![runmat_analysis_core::StructuralElement {
3218 element_id: "beam_1".to_string(),
3219 region_id: "span".to_string(),
3220 kind: runmat_analysis_core::StructuralElementKind::Beam(
3221 runmat_analysis_core::BeamElementModel {
3222 node_ids: [1, 2],
3223 section_id: "rect".to_string(),
3224 reference_axis: [0.0, 1.0, 0.0],
3225 },
3226 ),
3227 }],
3228 beam_sections: vec![runmat_analysis_core::BeamSectionModel {
3229 section_id: "rect".to_string(),
3230 area_m2: 1.0e-4,
3231 iy_m4: 1.0e-9,
3232 iz_m4: 1.0e-9,
3233 torsion_j_m4: 1.0e-9,
3234 outer_fiber_y_m: 0.01,
3235 outer_fiber_z_m: 0.01,
3236 torsion_outer_radius_m: 0.01,
3237 }],
3238 shell_sections: Vec::new(),
3239 });
3240
3241 let summary = assemble_linear_system(&model, None, Some(tetrahedron4_mesh()), None, None);
3242
3243 assert_eq!(summary.structural_solid_element_count, 1);
3244 assert_eq!(summary.structural_solid_recovery.len(), 1);
3245 assert_eq!(summary.structural_beam_element_count, 0);
3246 assert!(summary.operator.stiffness_csr.is_some());
3247 assert!(summary.operator.stiffness_dense.is_none());
3248 }
3249
3250 #[test]
3251 fn analysis_mesh_material_regions_select_assigned_solid_materials() {
3252 let mut soft_model = fixture_model(FixtureId::CantileverLinearStatic);
3253 let mut hard = soft_model.materials[0].clone();
3254 hard.material_id = "mat_hard".to_string();
3255 hard.mechanical.youngs_modulus_pa = 200.0e9;
3256 hard.mechanical.poisson_ratio = 0.3;
3257 let mut soft = hard.clone();
3258 soft.material_id = "mat_soft".to_string();
3259 soft.mechanical.youngs_modulus_pa = 20.0e9;
3260 soft_model.materials = vec![hard.clone(), soft.clone()];
3261 soft_model.material_assignments = vec![runmat_analysis_core::MaterialAssignment {
3262 region_id: "soft_region".to_string(),
3263 expected_material_id: "mat_hard".to_string(),
3264 assigned_material_id: "mat_soft".to_string(),
3265 confidence: runmat_analysis_core::EvidenceConfidence::Verified,
3266 }];
3267 let mut soft_mesh = tetrahedron4_mesh();
3268 soft_mesh.volume_elements[0].material_region_id = "soft_region".to_string();
3269
3270 let mut hard_model = soft_model.clone();
3271 hard_model.material_assignments = vec![runmat_analysis_core::MaterialAssignment {
3272 region_id: "soft_region".to_string(),
3273 expected_material_id: "mat_hard".to_string(),
3274 assigned_material_id: "mat_hard".to_string(),
3275 confidence: runmat_analysis_core::EvidenceConfidence::Verified,
3276 }];
3277 let hard_mesh = soft_mesh.clone();
3278
3279 let soft_summary = assemble_linear_system(&soft_model, None, Some(soft_mesh), None, None);
3280 let hard_summary = assemble_linear_system(&hard_model, None, Some(hard_mesh), None, None);
3281
3282 assert!(
3283 first_csr_diagonal(&soft_summary) < first_csr_diagonal(&hard_summary) * 0.2,
3284 "soft material assignment should lower solid element stiffness"
3285 );
3286 assert_eq!(
3287 soft_summary.structural_solid_recovery[0].region_id,
3288 "soft_region"
3289 );
3290 }
3291
3292 #[test]
3293 fn strict_analysis_mesh_assembly_rejects_invalid_tetrahedron4_stiffness() {
3294 let model = fixture_model(FixtureId::CantileverLinearStatic);
3295 let mut mesh = tetrahedron4_mesh();
3296 mesh.volume_elements[0].node_ids = vec![1, 3, 2, 4];
3297 mesh.boundary_faces = vec![
3298 boundary_face("root_face", vec![1, 2, 3], &["root"]),
3299 boundary_face("tip_face", vec![1, 2, 4], &["tip"]),
3300 ];
3301
3302 let err = try_assemble_linear_system(&model, None, Some(mesh), None, None).expect_err(
3303 "strict analysis mesh assembly should reject inverted Tetrahedron4 stiffness",
3304 );
3305
3306 assert!(matches!(
3307 err,
3308 LinearAssemblyError::SolidStiffness(SolidAssemblyError::ElementStiffness { .. })
3309 ));
3310 }
3311
3312 #[test]
3313 fn analysis_mesh_boundary_regions_drive_solid_loads_and_constraints() {
3314 let mut model = fixture_model(FixtureId::CantileverLinearStatic);
3315 model.boundary_conditions = vec![runmat_analysis_core::BoundaryCondition {
3316 bc_id: "fixed_root".to_string(),
3317 region_id: "root".to_string(),
3318 kind: BoundaryConditionKind::Fixed,
3319 }];
3320 model.loads = vec![runmat_analysis_core::LoadCase {
3321 load_id: "load_tip".to_string(),
3322 region_id: "tip".to_string(),
3323 kind: LoadKind::Force {
3324 fx: 0.0,
3325 fy: -12.0,
3326 fz: 0.0,
3327 },
3328 }];
3329 let mut mesh = tetrahedron4_mesh();
3330 mesh.boundary_faces = vec![
3331 boundary_face("root_face", vec![1, 2, 3], &["root"]),
3332 boundary_face("tip_face", vec![1, 2, 4], &["tip"]),
3333 ];
3334
3335 let summary = assemble_linear_system(&model, None, Some(mesh), None, None);
3336
3337 assert_eq!(summary.constrained_dof_count, 9);
3338 assert!(summary.operator.constrained[0]);
3339 assert!(summary.operator.constrained[1]);
3340 assert!(summary.operator.constrained[2]);
3341 assert_eq!(summary.operator.rhs[10], -4.0);
3342 }
3343
3344 #[test]
3345 fn strict_analysis_mesh_assembly_rejects_unmapped_load_region() {
3346 let mut model = fixture_model(FixtureId::CantileverLinearStatic);
3347 model.boundary_conditions = vec![runmat_analysis_core::BoundaryCondition {
3348 bc_id: "fixed_root".to_string(),
3349 region_id: "root".to_string(),
3350 kind: BoundaryConditionKind::Fixed,
3351 }];
3352 model.loads = vec![runmat_analysis_core::LoadCase {
3353 load_id: "load_tip".to_string(),
3354 region_id: "missing_tip".to_string(),
3355 kind: LoadKind::Force {
3356 fx: 0.0,
3357 fy: -12.0,
3358 fz: 0.0,
3359 },
3360 }];
3361 let mut mesh = tetrahedron4_mesh();
3362 mesh.boundary_faces = vec![boundary_face("root_face", vec![1, 2, 3], &["root"])];
3363
3364 let err = try_assemble_linear_system(&model, None, Some(mesh), None, None)
3365 .expect_err("strict analysis mesh assembly should reject unmapped loads");
3366
3367 assert!(matches!(
3368 err,
3369 LinearAssemblyError::AnalysisMeshRegionMapping(
3370 AnalysisMeshRegionMappingError::UnmappedLoadRegion { .. }
3371 )
3372 ));
3373 let message = err.to_string();
3374 assert!(message.contains("load_id=load_tip"));
3375 assert!(message.contains("region_id=missing_tip"));
3376 }
3377
3378 #[test]
3379 fn strict_analysis_mesh_assembly_rejects_unmapped_constraint_region() {
3380 let mut model = fixture_model(FixtureId::CantileverLinearStatic);
3381 model.boundary_conditions = vec![runmat_analysis_core::BoundaryCondition {
3382 bc_id: "fixed_root".to_string(),
3383 region_id: "missing_root".to_string(),
3384 kind: BoundaryConditionKind::Fixed,
3385 }];
3386 model.loads = vec![runmat_analysis_core::LoadCase {
3387 load_id: "load_tip".to_string(),
3388 region_id: "tip".to_string(),
3389 kind: LoadKind::Force {
3390 fx: 0.0,
3391 fy: -12.0,
3392 fz: 0.0,
3393 },
3394 }];
3395 let mut mesh = tetrahedron4_mesh();
3396 mesh.boundary_faces = vec![boundary_face("tip_face", vec![1, 2, 4], &["tip"])];
3397
3398 let err = try_assemble_linear_system(&model, None, Some(mesh), None, None)
3399 .expect_err("strict analysis mesh assembly should reject unmapped constraints");
3400
3401 assert!(matches!(
3402 err,
3403 LinearAssemblyError::AnalysisMeshRegionMapping(
3404 AnalysisMeshRegionMappingError::UnmappedBoundaryConditionRegion { .. }
3405 )
3406 ));
3407 let message = err.to_string();
3408 assert!(message.contains("bc_id=fixed_root"));
3409 assert!(message.contains("region_id=missing_root"));
3410 }
3411
3412 #[test]
3413 fn analysis_mesh_boundary_regions_integrate_pressure_loads() {
3414 let mut model = fixture_model(FixtureId::CantileverLinearStatic);
3415 model.boundary_conditions = Vec::new();
3416 model.loads = vec![runmat_analysis_core::LoadCase {
3417 load_id: "pressure_tip".to_string(),
3418 region_id: "tip".to_string(),
3419 kind: LoadKind::Pressure { magnitude_pa: 12.0 },
3420 }];
3421 let mut mesh = tetrahedron4_mesh();
3422 mesh.boundary_faces = vec![boundary_face("tip_face", vec![1, 2, 4], &["tip"])];
3423
3424 let summary = assemble_linear_system(&model, None, Some(mesh), None, None);
3425
3426 assert_close(summary.operator.rhs[1], -2.0);
3427 assert_close(summary.operator.rhs[4], -2.0);
3428 assert_close(summary.operator.rhs[10], -2.0);
3429 }
3430
3431 #[test]
3432 fn analysis_mesh_boundary_regions_lower_wrench_moments() {
3433 let mut model = fixture_model(FixtureId::CantileverLinearStatic);
3434 model.boundary_conditions = Vec::new();
3435 model.loads = vec![runmat_analysis_core::LoadCase {
3436 load_id: "wrench_tip".to_string(),
3437 region_id: "tip".to_string(),
3438 kind: LoadKind::Wrench {
3439 fx: 0.0,
3440 fy: 0.0,
3441 fz: 0.0,
3442 mx: 0.0,
3443 my: 6.0,
3444 mz: 0.0,
3445 px: 0.0,
3446 py: 0.0,
3447 pz: 0.0,
3448 },
3449 }];
3450 let mut mesh = tetrahedron4_mesh();
3451 mesh.boundary_faces = vec![boundary_face("tip_face", vec![1, 2, 4], &["tip"])];
3452
3453 let summary = assemble_linear_system(&model, None, Some(mesh), None, None);
3454
3455 assert_eq!(summary.structural_wrench_lowering.len(), 1);
3456 let lowering = &summary.structural_wrench_lowering[0];
3457 assert_eq!(lowering.load_id, "wrench_tip");
3458 assert_eq!(lowering.region_id, "tip");
3459 assert!(lowering.moment_couple_applied);
3460 assert_close(lowering.applied_moment_at_point[1], 6.0);
3461 assert_close(lowering.moment_residual[1], 0.0);
3462 assert!(summary.operator.rhs.iter().any(|value| value.abs() > 0.0));
3463 }
3464
3465 fn tetrahedron4_mesh() -> AnalysisMeshArtifact {
3466 let mut mesh = AnalysisMeshArtifact {
3467 schema_version: ANALYSIS_MESH_SCHEMA_VERSION.to_string(),
3468 mesh_id: "unit_tetrahedron".to_string(),
3469 nodes: vec![
3470 node(1, [0.0, 0.0, 0.0]),
3471 node(2, [1.0, 0.0, 0.0]),
3472 node(3, [0.0, 1.0, 0.0]),
3473 node(4, [0.0, 0.0, 1.0]),
3474 ],
3475 volume_elements: vec![AnalysisVolumeElement {
3476 element_id: "tetrahedron_1".to_string(),
3477 kind: VolumeElementKind::Tetrahedron4,
3478 node_ids: vec![1, 2, 3, 4],
3479 material_region_id: "solid".to_string(),
3480 provenance: Vec::new(),
3481 }],
3482 boundary_faces: Vec::new(),
3483 boundary_edges: Vec::new(),
3484 quality: AnalysisMeshQualityReport::default(),
3485 sizing: MeshSizingField::default(),
3486 field_topology: Vec::new(),
3487 backend: Default::default(),
3488 adaptive_iterations: Vec::new(),
3489 provenance: AnalysisMeshProvenance {
3490 algorithm: "test".to_string(),
3491 source_geometry_id: "geo:test".to_string(),
3492 source_geometry_revision: 1,
3493 source_geometry_sha256: None,
3494 },
3495 };
3496 mesh.refresh_field_topology();
3497 mesh
3498 }
3499
3500 fn boundary_face(
3501 face_id: &str,
3502 node_ids: Vec<u32>,
3503 region_ids: &[&str],
3504 ) -> AnalysisBoundaryFace {
3505 AnalysisBoundaryFace {
3506 face_id: face_id.to_string(),
3507 kind: BoundaryElementKind::Tri3,
3508 node_ids,
3509 adjacent_volume_element_ids: Vec::new(),
3510 region_ids: region_ids
3511 .iter()
3512 .map(|region| (*region).to_string())
3513 .collect(),
3514 provenance: Vec::new(),
3515 }
3516 }
3517
3518 fn assert_close(actual: f64, expected: f64) {
3519 assert!(
3520 (actual - expected).abs() <= 1.0e-8,
3521 "expected {expected}, got {actual}"
3522 );
3523 }
3524
3525 fn first_csr_diagonal(summary: &AssemblySummary) -> f64 {
3526 let csr = summary
3527 .operator
3528 .stiffness_csr
3529 .as_ref()
3530 .expect("analysis mesh should assemble CSR stiffness");
3531 csr.column_indices[csr.row_offsets[0]..csr.row_offsets[1]]
3532 .iter()
3533 .zip(csr.values[csr.row_offsets[0]..csr.row_offsets[1]].iter())
3534 .find_map(|(&column, &value)| (column == 0).then_some(value.abs()))
3535 .expect("first row should contain diagonal")
3536 }
3537
3538 fn node(node_id: u32, coordinates_m: [f64; 3]) -> AnalysisMeshNode {
3539 AnalysisMeshNode {
3540 node_id,
3541 coordinates_m,
3542 provenance: Vec::new(),
3543 }
3544 }
3545}