runmat_analysis_fea/assembly/
mod.rs

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

runmat_analysis_fea/assembly/
mod.rs
