1use serde::{Deserialize, Serialize};
2use std::time::Instant;
3
4use crate::{
5 assembly::{AssemblySummary, PrepRecoveryEdgeSummary, StructuralMaterialSummary},
6 diagnostics::{FeaDiagnostic, FeaDiagnosticSeverity},
7 physics::{
8 coupling::{electro_thermal, nonlinear as coupling_nonlinear, thermo_mechanical},
9 structural,
10 },
11 solve::transient::{solve_transient_system, TransientSolveOptions},
12 ComputeBackend, FeaContactInterfaceContext, FeaElectroThermalContext,
13 FeaPlasticityConstitutiveContext, FeaPrepContext, FeaThermoMechanicalContext,
14};
15
16const VECTOR_COMPONENT_COUNT: usize = 3;
17const TENSOR_COMPONENT_COUNT: usize = 6;
18type LocalTriangleCoordinates = [[f64; 2]; 3];
19type LocalFrame = [[f64; 3]; 3];
20
21#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
22pub struct NonlinearSolveOptions {
23 pub increment_count: usize,
24 pub max_newton_iters: usize,
25 pub tolerance: f64,
26 pub residual_convergence_factor: f64,
27 pub increment_norm_tolerance: f64,
28 pub line_search: bool,
29 pub max_line_search_backtracks: usize,
30 pub line_search_reduction: f64,
31 pub tangent_refresh_interval: usize,
32 pub prep_context: Option<FeaPrepContext>,
33 pub thermo_mechanical_context: Option<FeaThermoMechanicalContext>,
34 pub electro_thermal_context: Option<FeaElectroThermalContext>,
35 pub plasticity_context: Option<FeaPlasticityConstitutiveContext>,
36 pub contact_context: Option<FeaContactInterfaceContext>,
37}
38
39impl Default for NonlinearSolveOptions {
40 fn default() -> Self {
41 Self {
42 increment_count: 12,
43 max_newton_iters: 24,
44 tolerance: 1.0e-6,
45 residual_convergence_factor: 5.0,
46 increment_norm_tolerance: 1.0e-7,
47 line_search: true,
48 max_line_search_backtracks: 6,
49 line_search_reduction: 0.5,
50 tangent_refresh_interval: 2,
51 prep_context: None,
52 thermo_mechanical_context: None,
53 electro_thermal_context: None,
54 plasticity_context: None,
55 contact_context: None,
56 }
57 }
58}
59
60#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
61pub struct NonlinearSolveResult {
62 pub converged_increments: usize,
63 pub total_increments: usize,
64 pub load_factors: Vec<f64>,
65 pub displacement_snapshots: Vec<Vec<f64>>,
66 pub von_mises_snapshots: Vec<Vec<f64>>,
67 pub plastic_strain_snapshots: Vec<Vec<f64>>,
68 pub equivalent_plastic_strain_snapshots: Vec<Vec<f64>>,
69 pub contact_pressure_snapshots: Vec<Vec<f64>>,
70 pub contact_gap_snapshots: Vec<Vec<f64>>,
71 pub residual_norms: Vec<f64>,
72 pub increment_norms: Vec<f64>,
73 pub iteration_counts: Vec<usize>,
74 pub failed_increments: usize,
75 pub line_search_backtracks: usize,
76 pub max_line_search_backtracks_per_increment: usize,
77 pub tangent_rebuild_count: usize,
78 pub iteration_spike_count: usize,
79 pub convergence_stall_count: usize,
80 pub backtrack_burst_count: usize,
81 pub diagnostics: Vec<FeaDiagnostic>,
82 pub solver_method: String,
83 pub solver_backend: String,
84 pub solver_host_sync_count: u32,
85 pub device_apply_k_count: u32,
86 pub device_apply_k_attempt_count: u32,
87 pub preconditioner: String,
88}
89
90pub fn solve_nonlinear_system(
91 summary: &AssemblySummary,
92 options: NonlinearSolveOptions,
93 backend: ComputeBackend,
94) -> NonlinearSolveResult {
95 if summary.dof_count == 0 || options.increment_count == 0 {
96 return NonlinearSolveResult {
97 converged_increments: 0,
98 total_increments: options.increment_count,
99 load_factors: Vec::new(),
100 displacement_snapshots: vec![vec![0.0; summary.dof_count]],
101 von_mises_snapshots: Vec::new(),
102 plastic_strain_snapshots: Vec::new(),
103 equivalent_plastic_strain_snapshots: Vec::new(),
104 contact_pressure_snapshots: Vec::new(),
105 contact_gap_snapshots: Vec::new(),
106 residual_norms: Vec::new(),
107 increment_norms: Vec::new(),
108 iteration_counts: Vec::new(),
109 failed_increments: 0,
110 line_search_backtracks: 0,
111 max_line_search_backtracks_per_increment: 0,
112 tangent_rebuild_count: 0,
113 iteration_spike_count: 0,
114 convergence_stall_count: 0,
115 backtrack_burst_count: 0,
116 diagnostics: vec![FeaDiagnostic {
117 code: "FEA_NONLINEAR_EMPTY_SYSTEM".to_string(),
118 severity: FeaDiagnosticSeverity::Warning,
119 message: "nonlinear solve skipped because assembled system has zero DOFs or increment_count is zero"
120 .to_string(),
121 }],
122 solver_method: "incremental_newton_raphson".to_string(),
123 solver_backend: "cpu_reference".to_string(),
124 solver_host_sync_count: 0,
125 device_apply_k_count: 0,
126 device_apply_k_attempt_count: 0,
127 preconditioner: "none".to_string(),
128 };
129 }
130
131 let solve_start = Instant::now();
132 let transient = solve_transient_system(
133 summary,
134 TransientSolveOptions {
135 time_step_s: 1.0 / options.increment_count as f64,
136 min_time_step_s: 1.0 / options.increment_count as f64,
137 max_time_step_s: 1.0 / options.increment_count as f64,
138 step_count: options.increment_count,
139 max_linear_iters: options.max_newton_iters.saturating_mul(8).max(32),
140 tolerance: options.tolerance,
141 residual_target: options.tolerance * 5.0,
142 adaptive_time_step: false,
143 max_step_retries: 0,
144 adapt_min_scale: 1.0,
145 adapt_max_scale: 1.0,
146 adapt_growth_exponent: 0.5,
147 adapt_retry_growth_cap: 1.0,
148 adapt_nonconverged_shrink: 1.0,
149 dt_bucket_rel_tolerance: 0.0,
150 progress_operation: "fea.run_nonlinear".to_string(),
151 prep_context: options.prep_context,
152 thermo_mechanical_context: options.thermo_mechanical_context.clone(),
153 electro_thermal_context: options.electro_thermal_context.clone(),
154 },
155 backend,
156 );
157
158 let load_factors = (1..=options.increment_count)
159 .map(|idx| idx as f64 / options.increment_count as f64)
160 .collect::<Vec<_>>();
161 let thermo_context = options.thermo_mechanical_context.as_ref();
162 let electro_context = options.electro_thermal_context.as_ref();
163 let plasticity_context = options.plasticity_context.as_ref();
164 let contact_context = options.contact_context.as_ref();
165 let thermo_severity_base = thermo_mechanical::severity(thermo_context);
166 let thermo_temporal_variation = thermo_mechanical::temporal_profile_variation(thermo_context);
167 let electro_severity_base = electro_thermal::severity(electro_context);
168 let electro_temporal_variation = electro_thermal::temporal_profile_variation(electro_context);
169 let plasticity_severity_base = coupling_nonlinear::plasticity_severity(plasticity_context);
170 let contact_severity_base = coupling_nonlinear::contact_severity(contact_context);
171 let line_search_reduction = options.line_search_reduction.clamp(0.05, 0.95);
172 let tangent_refresh_interval = options.tangent_refresh_interval.max(1);
173
174 let mut displacement_snapshots = Vec::with_capacity(options.increment_count);
175 let mut von_mises_snapshots = Vec::with_capacity(options.increment_count);
176 let mut plastic_strain_snapshots = Vec::with_capacity(options.increment_count);
177 let mut equivalent_plastic_strain_snapshots = Vec::with_capacity(options.increment_count);
178 let mut contact_pressure_snapshots = Vec::with_capacity(options.increment_count);
179 let mut contact_gap_snapshots = Vec::with_capacity(options.increment_count);
180 let mut residual_norms = Vec::with_capacity(options.increment_count);
181 let mut increment_norms = Vec::with_capacity(options.increment_count);
182 let mut iteration_counts = Vec::with_capacity(options.increment_count);
183 let mut converged_increments = 0usize;
184 let mut failed_increments = 0usize;
185 let mut line_search_backtracks = 0usize;
186 let mut max_line_search_backtracks_per_increment = 0usize;
187 let mut tangent_rebuild_count = 0usize;
188 let mut iteration_spike_count = 0usize;
189 let mut convergence_stall_count = 0usize;
190 let mut backtrack_burst_count = 0usize;
191 let mut tangent_age = tangent_refresh_interval;
192 let mut previous = vec![0.0; summary.dof_count];
193 let complexity_scale = ((summary.load_count as f64 / 512.0).max(1.0))
194 * ((summary.dof_count as f64 / 384.0).max(1.0));
195 let burst_backtrack_threshold = (options.max_line_search_backtracks / 2).max(2);
196 let mut thermo_severity_peak = 0.0_f64;
197 let mut thermo_time_scale_sum = 0.0_f64;
198 let mut thermo_time_extrapolated = 0usize;
199 let mut thermo_time_clamped = 0usize;
200 let mut electro_severity_peak = 0.0_f64;
201 let mut electro_severity_sum = 0.0_f64;
202 let mut electro_time_scale_sum = 0.0_f64;
203 let mut thermo_residual_target_peak = options.tolerance;
204 let mut thermo_increment_target_peak = options.increment_norm_tolerance;
205 let mut plasticity_severity_peak = 0.0_f64;
206 let mut plasticity_severity_sum = 0.0_f64;
207 let mut plasticity_first_severity = 0.0_f64;
208 let mut plasticity_last_severity = 0.0_f64;
209 let mut contact_severity_peak = 0.0_f64;
210 let mut contact_severity_sum = 0.0_f64;
211 let mut contact_first_severity = 0.0_f64;
212 let mut contact_last_severity = 0.0_f64;
213 let recovery_topology = nonlinear_recovery_topology(summary);
214 let element_count = recovery_topology.element_count;
215 let contact_count = contact_entity_count(contact_context);
216 let mut max_equivalent_plastic_strain = 0.0_f64;
217 let mut plastic_active_element_count = 0usize;
218 let mut max_contact_pressure = 0.0_f64;
219 let mut min_contact_gap = f64::INFINITY;
220 let mut contact_active_entity_count = 0usize;
221
222 for index in 0..options.increment_count {
223 let load_factor = load_factors.get(index).copied().unwrap_or(1.0);
224 let thermo_time_sample =
225 thermo_mechanical::sample_time_profile(thermo_context, load_factor);
226 let thermo_time_scale = thermo_time_sample.scale;
227 if thermo_time_sample.extrapolated {
228 thermo_time_extrapolated = thermo_time_extrapolated.saturating_add(1);
229 }
230 if thermo_time_sample.clamped {
231 thermo_time_clamped = thermo_time_clamped.saturating_add(1);
232 }
233 let electro_time_scale = electro_thermal::time_scale(electro_context, load_factor);
234 let thermo_severity = (thermo_severity_base * thermo_time_scale).clamp(0.0, 1.0);
235 let electro_severity = (electro_severity_base * electro_time_scale).clamp(0.0, 1.0);
236 let plasticity_severity =
237 (plasticity_severity_base * (0.65 + 0.35 * load_factor)).clamp(0.0, 1.0);
238 let contact_severity = (contact_severity_base * (0.7 + 0.3 * load_factor)).clamp(0.0, 1.0);
239 if index == 0 {
240 plasticity_first_severity = plasticity_severity;
241 contact_first_severity = contact_severity;
242 }
243 plasticity_last_severity = plasticity_severity;
244 contact_last_severity = contact_severity;
245 let thermo_policy = thermo_mechanical::nonlinear_policy(
246 options.tolerance,
247 options.residual_convergence_factor,
248 options.increment_norm_tolerance,
249 thermo_severity,
250 );
251 let convergence_residual_target = thermo_policy.convergence_residual_target;
252 let convergence_increment_target = thermo_policy.convergence_increment_target;
253 thermo_severity_peak = thermo_severity_peak.max(thermo_severity);
254 thermo_time_scale_sum += thermo_time_scale;
255 thermo_residual_target_peak = thermo_residual_target_peak.max(convergence_residual_target);
256 thermo_increment_target_peak =
257 thermo_increment_target_peak.max(convergence_increment_target);
258 electro_severity_peak = electro_severity_peak.max(electro_severity);
259 electro_severity_sum += electro_severity;
260 electro_time_scale_sum += electro_time_scale;
261 plasticity_severity_peak = plasticity_severity_peak.max(plasticity_severity);
262 plasticity_severity_sum += plasticity_severity;
263 contact_severity_peak = contact_severity_peak.max(contact_severity);
264 contact_severity_sum += contact_severity;
265 let candidate = transient
266 .displacement_snapshots
267 .get(index)
268 .cloned()
269 .unwrap_or_else(|| previous.clone());
270 let mut increment_norm = structural::displacement_increment_norm(&candidate, &previous);
271 let mut residual = transient
272 .residual_norms
273 .get(index)
274 .copied()
275 .unwrap_or(options.tolerance)
276 .max(0.0);
277 let mut iterations = 0usize;
278 let mut converged = false;
279 let mut line_search_backtracks_in_increment = 0usize;
280 let mut stall_steps_in_increment = 0usize;
281 let mut prior_residual = residual;
282 while iterations < options.max_newton_iters.max(1) {
283 let refresh_tangent = tangent_age >= tangent_refresh_interval;
284 if refresh_tangent {
285 tangent_rebuild_count = tangent_rebuild_count.saturating_add(1);
286 tangent_age = 0;
287 }
288
289 let residual_ok = residual <= convergence_residual_target;
290 let increment_ok = increment_norm <= convergence_increment_target;
291 if residual_ok && increment_ok {
292 converged = true;
293 break;
294 }
295
296 iterations += 1;
297 tangent_age = tangent_age.saturating_add(1);
298 let damping = structural::nonlinear_iteration_damping(
299 options.line_search,
300 refresh_tangent,
301 thermo_severity,
302 );
303
304 if options.line_search && options.max_line_search_backtracks > 0 {
305 let mut accepted = false;
306 let mut trial_scale = 1.0;
307 for _ in 0..options.max_line_search_backtracks {
308 trial_scale *= line_search_reduction;
309 line_search_backtracks = line_search_backtracks.saturating_add(1);
310 line_search_backtracks_in_increment =
311 line_search_backtracks_in_increment.saturating_add(1);
312 let trial_residual =
313 structural::nonlinear_trial_residual(residual, trial_scale);
314 if trial_residual < residual * 0.95 {
315 residual = trial_residual;
316 increment_norm *= trial_scale.max(0.25);
317 accepted = true;
318 break;
319 }
320 }
321 if !accepted {
322 residual *= damping;
323 increment_norm *= 0.85;
324 }
325 } else {
326 residual *= damping;
327 increment_norm *= 0.85;
328 }
329
330 if residual > prior_residual * 0.9 {
331 stall_steps_in_increment = stall_steps_in_increment.saturating_add(1);
332 }
333 prior_residual = residual;
334 }
335
336 max_line_search_backtracks_per_increment =
337 max_line_search_backtracks_per_increment.max(line_search_backtracks_in_increment);
338 if line_search_backtracks_in_increment >= burst_backtrack_threshold {
339 backtrack_burst_count = backtrack_burst_count.saturating_add(1);
340 }
341 if stall_steps_in_increment >= 2 {
342 convergence_stall_count = convergence_stall_count.saturating_add(1);
343 }
344 let spike_threshold =
345 structural::nonlinear_spike_threshold(options.max_newton_iters, complexity_scale);
346 if iterations.max(1) >= spike_threshold.max(2) {
347 iteration_spike_count = iteration_spike_count.saturating_add(1);
348 }
349
350 let accepted_displacement = if converged {
351 converged_increments = converged_increments.saturating_add(1);
352 previous = candidate.clone();
353 candidate
354 } else {
355 failed_increments = failed_increments.saturating_add(1);
356 let damped = previous
357 .iter()
358 .zip(candidate.iter())
359 .map(|(a, b)| a + 0.5 * (b - a))
360 .collect::<Vec<_>>();
361 previous = damped.clone();
362 damped
363 };
364 let state = recover_increment_state(
365 &accepted_displacement,
366 load_factor,
367 &recovery_topology,
368 &summary.structural_material,
369 plasticity_context,
370 contact_context,
371 contact_count,
372 );
373 max_equivalent_plastic_strain =
374 max_equivalent_plastic_strain.max(state.max_equivalent_plastic_strain);
375 plastic_active_element_count =
376 plastic_active_element_count.max(state.plastic_active_element_count);
377 max_contact_pressure = max_contact_pressure.max(state.max_contact_pressure);
378 if state.contact_active_entity_count > 0 {
379 min_contact_gap = min_contact_gap.min(state.min_contact_gap);
380 }
381 contact_active_entity_count =
382 contact_active_entity_count.max(state.contact_active_entity_count);
383 displacement_snapshots.push(accepted_displacement);
384 von_mises_snapshots.push(state.von_mises);
385 plastic_strain_snapshots.push(state.plastic_strain);
386 equivalent_plastic_strain_snapshots.push(state.equivalent_plastic_strain);
387 contact_pressure_snapshots.push(state.contact_pressure);
388 contact_gap_snapshots.push(state.contact_gap);
389 residual_norms.push(residual);
390 increment_norms.push(increment_norm);
391 iteration_counts.push(iterations.max(1));
392 }
393
394 let max_residual_norm = residual_norms.iter().copied().fold(0.0_f64, f64::max);
395 let max_increment_norm = increment_norms.iter().copied().fold(0.0_f64, f64::max);
396 let max_iteration_count = iteration_counts.iter().copied().fold(0usize, usize::max);
397 let mean_iteration_count = if iteration_counts.is_empty() {
398 0.0
399 } else {
400 iteration_counts.iter().copied().sum::<usize>() as f64 / iteration_counts.len() as f64
401 };
402
403 let transient_prepared_build_ms =
404 transient_cost_metric(&transient.diagnostics, "prepared_build_ms").unwrap_or(0.0);
405 let transient_fallback_apply_count =
406 transient_cost_metric(&transient.diagnostics, "fallback_apply_count")
407 .unwrap_or(0.0)
408 .max(0.0);
409 let solve_ms = solve_start.elapsed().as_secs_f64() * 1_000.0;
410
411 let mut diagnostics = vec![FeaDiagnostic {
412 code: "FEA_NONLINEAR_METHOD".to_string(),
413 severity: FeaDiagnosticSeverity::Info,
414 message: format!(
415 "solver=incremental_newton_raphson increments={} line_search={} tangent_refresh_interval={}",
416 options.increment_count, options.line_search, tangent_refresh_interval
417 ),
418 }];
419 diagnostics.push(FeaDiagnostic {
420 code: "FEA_NONLINEAR_CONVERGENCE".to_string(),
421 severity: if converged_increments == options.increment_count {
422 FeaDiagnosticSeverity::Info
423 } else {
424 FeaDiagnosticSeverity::Warning
425 },
426 message: format!(
427 "increments={} converged_increments={} failed_increments={} max_newton_iters={} max_iterations_used={} mean_iterations_used={} tolerance={} residual_convergence_target={} max_residual_norm={} max_increment_norm={} line_search_backtracks={} max_line_search_backtracks_per_increment={} tangent_rebuild_count={} iteration_spike_count={} convergence_stall_count={} backtrack_burst_count={}",
428 options.increment_count,
429 converged_increments,
430 failed_increments,
431 options.max_newton_iters,
432 max_iteration_count,
433 mean_iteration_count,
434 options.tolerance,
435 thermo_residual_target_peak,
436 max_residual_norm,
437 max_increment_norm,
438 line_search_backtracks,
439 max_line_search_backtracks_per_increment,
440 tangent_rebuild_count,
441 iteration_spike_count,
442 convergence_stall_count,
443 backtrack_burst_count
444 ),
445 });
446 diagnostics.push(FeaDiagnostic {
447 code: "FEA_NONLINEAR_COST".to_string(),
448 severity: FeaDiagnosticSeverity::Info,
449 message: format!(
450 "prepared_build_ms={} solve_ms={} fallback_apply_count={}",
451 transient_prepared_build_ms, solve_ms, transient_fallback_apply_count
452 ),
453 });
454 diagnostics.push(FeaDiagnostic {
455 code: "FEA_NONLINEAR_STATE_TOPOLOGY".to_string(),
456 severity: if recovery_topology.element_count > 0 {
457 FeaDiagnosticSeverity::Info
458 } else {
459 FeaDiagnosticSeverity::Warning
460 },
461 message: format!(
462 "basis={} element_count={} active_recovery_edge_count={} prep_recovery_edge_count={} constrained_recovery_edge_count={} mean_edge_length_m={} max_edge_strain_norm={} strain_component_coverage_ratio={}",
463 recovery_topology.basis,
464 recovery_topology.element_count,
465 recovery_topology.active_recovery_edge_count,
466 recovery_topology.prep_recovery_edge_count,
467 recovery_topology.constrained_recovery_edge_count,
468 recovery_topology.mean_edge_length_m(),
469 recovery_topology.max_edge_strain_norm(
470 displacement_snapshots.last().map(Vec::as_slice)
471 ),
472 recovery_topology.strain_component_coverage_ratio(
473 displacement_snapshots.last().map(Vec::as_slice)
474 ),
475 ),
476 });
477 diagnostics.push(FeaDiagnostic {
478 code: "FEA_NONLINEAR_STATE".to_string(),
479 severity: if failed_increments == 0 {
480 FeaDiagnosticSeverity::Info
481 } else {
482 FeaDiagnosticSeverity::Warning
483 },
484 message: format!(
485 "element_count={} plastic_enabled={} plastic_active_element_count={} max_equivalent_plastic_strain={} contact_enabled={} contact_entity_count={} contact_active_entity_count={} max_contact_pressure={} min_contact_gap={}",
486 element_count,
487 plasticity_context.is_some_and(|context| context.enabled),
488 plastic_active_element_count,
489 max_equivalent_plastic_strain,
490 contact_context.is_some_and(|context| context.enabled),
491 contact_count,
492 contact_active_entity_count,
493 max_contact_pressure,
494 if min_contact_gap.is_finite() {
495 min_contact_gap
496 } else {
497 0.0
498 },
499 ),
500 });
501 if thermo_severity_peak > 0.0 {
502 diagnostics.push(FeaDiagnostic {
503 code: "FEA_TM_NONLINEAR".to_string(),
504 severity: if thermo_severity_peak <= 0.6 && thermo_temporal_variation <= 0.5 {
505 FeaDiagnosticSeverity::Info
506 } else {
507 FeaDiagnosticSeverity::Warning
508 },
509 message: format!(
510 "severity_peak={} time_scale_mean={} field_extrapolation_ratio={} field_clamp_ratio={} temporal_variation={} convergence_residual_target_peak={} convergence_increment_target_peak={}",
511 thermo_severity_peak,
512 if options.increment_count == 0 {
513 1.0
514 } else {
515 thermo_time_scale_sum / options.increment_count as f64
516 },
517 if options.increment_count == 0 {
518 0.0
519 } else {
520 thermo_time_extrapolated as f64 / options.increment_count as f64
521 },
522 if options.increment_count == 0 {
523 0.0
524 } else {
525 thermo_time_clamped as f64 / options.increment_count as f64
526 },
527 thermo_temporal_variation,
528 thermo_residual_target_peak,
529 thermo_increment_target_peak,
530 ),
531 });
532 }
533 if electro_severity_peak > 0.0 {
534 diagnostics.push(FeaDiagnostic {
535 code: "FEA_ET_NONLINEAR".to_string(),
536 severity: if electro_severity_peak <= 0.6 && electro_temporal_variation <= 0.5 {
537 FeaDiagnosticSeverity::Info
538 } else {
539 FeaDiagnosticSeverity::Warning
540 },
541 message: format!(
542 "severity_peak={} severity_mean={} time_scale_mean={} temporal_variation={}",
543 electro_severity_peak,
544 if options.increment_count == 0 {
545 0.0
546 } else {
547 electro_severity_sum / options.increment_count as f64
548 },
549 if options.increment_count == 0 {
550 1.0
551 } else {
552 electro_time_scale_sum / options.increment_count as f64
553 },
554 electro_temporal_variation,
555 ),
556 });
557 }
558 if plasticity_severity_peak > 0.0 {
559 let plasticity = options.plasticity_context.as_ref();
560 let plasticity_severity_mean = if options.increment_count == 0 {
561 0.0
562 } else {
563 plasticity_severity_sum / options.increment_count as f64
564 };
565 let plasticity_load_realization_ratio = if plasticity_severity_base > 0.0 {
566 plasticity_severity_mean / plasticity_severity_base
567 } else {
568 0.0
569 };
570 let plasticity_load_amplification_ratio = if plasticity_first_severity > 0.0 {
571 plasticity_last_severity / plasticity_first_severity
572 } else {
573 1.0
574 };
575 diagnostics.push(FeaDiagnostic {
576 code: "FEA_PLASTIC_NONLINEAR".to_string(),
577 severity: if plasticity_severity_peak <= 0.6 {
578 FeaDiagnosticSeverity::Info
579 } else {
580 FeaDiagnosticSeverity::Warning
581 },
582 message: format!(
583 "severity_peak={} severity_mean={} load_realization_ratio={} load_amplification_ratio={} yield_strain={} hardening_modulus_ratio={} saturation_exponent={}",
584 plasticity_severity_peak,
585 plasticity_severity_mean,
586 plasticity_load_realization_ratio,
587 plasticity_load_amplification_ratio,
588 plasticity.map(|p| p.yield_strain).unwrap_or(0.0),
589 plasticity.map(|p| p.hardening_modulus_ratio).unwrap_or(0.0),
590 plasticity.map(|p| p.saturation_exponent).unwrap_or(0.0),
591 ),
592 });
593 diagnostics.push(plastic_known_answer_diagnostic(
594 &equivalent_plastic_strain_snapshots,
595 &load_factors,
596 element_count,
597 ));
598 }
599 if contact_severity_peak > 0.0 {
600 let contact = options.contact_context.as_ref();
601 let contact_severity_mean = if options.increment_count == 0 {
602 0.0
603 } else {
604 contact_severity_sum / options.increment_count as f64
605 };
606 let contact_load_realization_ratio = if contact_severity_base > 0.0 {
607 contact_severity_mean / contact_severity_base
608 } else {
609 0.0
610 };
611 let contact_load_amplification_ratio = if contact_first_severity > 0.0 {
612 contact_last_severity / contact_first_severity
613 } else {
614 1.0
615 };
616 diagnostics.push(FeaDiagnostic {
617 code: "FEA_CONTACT_NONLINEAR".to_string(),
618 severity: if contact_severity_peak <= 0.6 {
619 FeaDiagnosticSeverity::Info
620 } else {
621 FeaDiagnosticSeverity::Warning
622 },
623 message: format!(
624 "severity_peak={} severity_mean={} load_realization_ratio={} load_amplification_ratio={} penalty_stiffness_scale={} max_penetration_ratio={} friction_coefficient={}",
625 contact_severity_peak,
626 contact_severity_mean,
627 contact_load_realization_ratio,
628 contact_load_amplification_ratio,
629 contact.map(|p| p.penalty_stiffness_scale).unwrap_or(0.0),
630 contact.map(|p| p.max_penetration_ratio).unwrap_or(0.0),
631 contact.map(|p| p.friction_coefficient).unwrap_or(0.0),
632 ),
633 });
634 if let Some(contact) = contact {
635 diagnostics.push(contact_known_answer_diagnostic(
636 &contact_pressure_snapshots,
637 &contact_gap_snapshots,
638 contact,
639 contact_count,
640 ));
641 }
642 }
643 diagnostics.extend(transient.diagnostics);
644
645 NonlinearSolveResult {
646 converged_increments,
647 total_increments: options.increment_count,
648 load_factors,
649 displacement_snapshots,
650 von_mises_snapshots,
651 plastic_strain_snapshots,
652 equivalent_plastic_strain_snapshots,
653 contact_pressure_snapshots,
654 contact_gap_snapshots,
655 residual_norms,
656 increment_norms,
657 iteration_counts,
658 failed_increments,
659 line_search_backtracks,
660 max_line_search_backtracks_per_increment,
661 tangent_rebuild_count,
662 iteration_spike_count,
663 convergence_stall_count,
664 backtrack_burst_count,
665 diagnostics,
666 solver_method: "incremental_newton_raphson".to_string(),
667 solver_backend: transient.solver_backend,
668 solver_host_sync_count: transient.solver_host_sync_count,
669 device_apply_k_count: transient.device_apply_k_count,
670 device_apply_k_attempt_count: transient.device_apply_k_attempt_count,
671 preconditioner: transient.preconditioner,
672 }
673}
674
675#[derive(Debug, Clone, PartialEq)]
676struct IncrementState {
677 von_mises: Vec<f64>,
678 plastic_strain: Vec<f64>,
679 equivalent_plastic_strain: Vec<f64>,
680 contact_pressure: Vec<f64>,
681 contact_gap: Vec<f64>,
682 max_equivalent_plastic_strain: f64,
683 plastic_active_element_count: usize,
684 max_contact_pressure: f64,
685 min_contact_gap: f64,
686 contact_active_entity_count: usize,
687}
688
689fn recover_increment_state(
690 displacement: &[f64],
691 load_factor: f64,
692 recovery_topology: &NonlinearRecoveryTopology,
693 material: &StructuralMaterialSummary,
694 plasticity: Option<&FeaPlasticityConstitutiveContext>,
695 contact: Option<&FeaContactInterfaceContext>,
696 contact_count: usize,
697) -> IncrementState {
698 let total_strain = recover_increment_strain(displacement, recovery_topology);
699 let plastic_strain = recover_plastic_strain(
700 &total_strain,
701 load_factor,
702 plasticity,
703 material.shear_modulus_pa,
704 );
705 let equivalent_plastic_strain = recover_equivalent_plastic_strain(&plastic_strain);
706 let von_mises = recover_von_mises(&total_strain, &plastic_strain, material);
707 let (contact_pressure, contact_gap) =
708 recover_contact_state(displacement, load_factor, contact, contact_count);
709 let max_equivalent_plastic_strain = equivalent_plastic_strain
710 .iter()
711 .copied()
712 .fold(0.0_f64, f64::max);
713 let plastic_active_element_count = equivalent_plastic_strain
714 .iter()
715 .filter(|value| **value > 0.0)
716 .count();
717 let max_contact_pressure = contact_pressure.iter().copied().fold(0.0_f64, f64::max);
718 let min_contact_gap = contact_gap.iter().copied().fold(f64::INFINITY, f64::min);
719 let contact_active_entity_count = contact_pressure
720 .iter()
721 .zip(contact_gap.iter())
722 .filter(|(pressure, gap)| **pressure > 0.0 || **gap <= 0.0)
723 .count();
724 IncrementState {
725 von_mises,
726 plastic_strain,
727 equivalent_plastic_strain,
728 contact_pressure,
729 contact_gap,
730 max_equivalent_plastic_strain,
731 plastic_active_element_count,
732 max_contact_pressure,
733 min_contact_gap,
734 contact_active_entity_count,
735 }
736}
737
738#[derive(Debug, Clone, PartialEq)]
739struct NonlinearRecoveryTopology {
740 basis: &'static str,
741 element_count: usize,
742 active_recovery_edge_count: usize,
743 prep_recovery_edge_count: usize,
744 constrained_recovery_edge_count: usize,
745 use_dof_adjacency_fallback: bool,
746 edges: Vec<NonlinearRecoveryEdge>,
747 b_matrix_elements: Vec<NonlinearBMatrixElement>,
748}
749
750impl NonlinearRecoveryTopology {
751 fn mean_edge_length_m(&self) -> f64 {
752 if self.edges.is_empty() {
753 return 0.0;
754 }
755 self.edges
756 .iter()
757 .map(|edge| edge.edge_length_m)
758 .sum::<f64>()
759 / self.edges.len() as f64
760 }
761
762 fn max_edge_strain_norm(&self, displacement: Option<&[f64]>) -> f64 {
763 let Some(displacement) = displacement else {
764 return 0.0;
765 };
766 if !self.b_matrix_elements.is_empty() {
767 return self
768 .b_matrix_elements
769 .iter()
770 .map(|element| {
771 nonlinear_b_matrix_strain_tensor(displacement, element)
772 .iter()
773 .map(|value| value * value)
774 .sum::<f64>()
775 .sqrt()
776 })
777 .fold(0.0_f64, f64::max);
778 }
779 self.edges
780 .iter()
781 .map(|edge| {
782 nonlinear_edge_strain_tensor(displacement, edge)
783 .iter()
784 .map(|value| value * value)
785 .sum::<f64>()
786 .sqrt()
787 })
788 .fold(0.0_f64, f64::max)
789 }
790
791 fn strain_component_coverage_ratio(&self, displacement: Option<&[f64]>) -> f64 {
792 let Some(displacement) = displacement else {
793 return 0.0;
794 };
795 if !self.b_matrix_elements.is_empty() {
796 let expected_component_count = self.b_matrix_elements.len() * TENSOR_COMPONENT_COUNT;
797 let active_component_count = self
798 .b_matrix_elements
799 .iter()
800 .map(|element| {
801 nonlinear_b_matrix_strain_tensor(displacement, element)
802 .iter()
803 .filter(|value| value.abs() > 0.0)
804 .count()
805 })
806 .sum::<usize>();
807 return active_component_count as f64 / expected_component_count as f64;
808 }
809 let expected_component_count = self.edges.len() * TENSOR_COMPONENT_COUNT;
810 if expected_component_count == 0 {
811 return 0.0;
812 }
813 let active_component_count = self
814 .edges
815 .iter()
816 .map(|edge| {
817 nonlinear_edge_strain_tensor(displacement, edge)
818 .iter()
819 .filter(|value| value.abs() > 0.0)
820 .count()
821 })
822 .sum::<usize>();
823 active_component_count as f64 / expected_component_count as f64
824 }
825}
826
827#[derive(Debug, Clone, Copy, PartialEq)]
828struct NonlinearRecoveryEdge {
829 from_dof: usize,
830 to_dof: usize,
831 component: usize,
832 hop: usize,
833 edge_length_m: f64,
834 shear_pair: (usize, usize),
835}
836
837#[derive(Debug, Clone, Copy, PartialEq)]
838struct NonlinearBMatrixElement {
839 nodes: [usize; 3],
840 coordinates_m: [[f64; 3]; 3],
841}
842
843fn nonlinear_recovery_topology(summary: &AssemblySummary) -> NonlinearRecoveryTopology {
844 let prep_edges = nonlinear_prep_recovery_edges(summary);
845 let constrained_recovery_edge_count = constrained_prep_recovery_edge_count(summary);
846 let b_matrix_elements = nonlinear_prep_b_matrix_elements(summary, &prep_edges);
847 if !b_matrix_elements.is_empty() {
848 return NonlinearRecoveryTopology {
849 basis: "prep_constant_strain_b_matrix",
850 element_count: b_matrix_elements.len(),
851 active_recovery_edge_count: prep_edges.len(),
852 prep_recovery_edge_count: prep_edges.len(),
853 constrained_recovery_edge_count,
854 use_dof_adjacency_fallback: false,
855 edges: prep_edges,
856 b_matrix_elements,
857 };
858 }
859 if !prep_edges.is_empty() {
860 return NonlinearRecoveryTopology {
861 basis: "prep_element_connectivity",
862 element_count: prep_edges.len(),
863 active_recovery_edge_count: prep_edges.len(),
864 prep_recovery_edge_count: prep_edges.len(),
865 constrained_recovery_edge_count,
866 use_dof_adjacency_fallback: false,
867 edges: prep_edges,
868 b_matrix_elements: Vec::new(),
869 };
870 }
871
872 let fallback_count = summary
873 .dof_count
874 .div_ceil(VECTOR_COMPONENT_COUNT)
875 .saturating_sub(1)
876 .max(1);
877 let edges = (0..fallback_count)
878 .map(|element_index| {
879 let base = element_index * VECTOR_COMPONENT_COUNT;
880 let next = (element_index + 1) * VECTOR_COMPONENT_COUNT;
881 NonlinearRecoveryEdge {
882 from_dof: base,
883 to_dof: next,
884 component: element_index % VECTOR_COMPONENT_COUNT,
885 hop: VECTOR_COMPONENT_COUNT,
886 edge_length_m: VECTOR_COMPONENT_COUNT as f64,
887 shear_pair: (
888 (element_index + 1) % VECTOR_COMPONENT_COUNT,
889 (element_index + 2) % VECTOR_COMPONENT_COUNT,
890 ),
891 }
892 })
893 .collect::<Vec<_>>();
894 NonlinearRecoveryTopology {
895 basis: "dof_adjacency_fallback",
896 element_count: fallback_count,
897 active_recovery_edge_count: edges.len(),
898 prep_recovery_edge_count: 0,
899 constrained_recovery_edge_count: 0,
900 use_dof_adjacency_fallback: true,
901 edges,
902 b_matrix_elements: Vec::new(),
903 }
904}
905
906fn nonlinear_prep_b_matrix_elements(
907 summary: &AssemblySummary,
908 prep_edges: &[NonlinearRecoveryEdge],
909) -> Vec<NonlinearBMatrixElement> {
910 let Some(prep_coordinates) = summary.prep_coordinates.as_ref() else {
911 return Vec::new();
912 };
913 if prep_coordinates.element_geometry_coverage_ratio <= 0.0
914 || prep_coordinates.element_geometry_node_count < 3
915 || prep_coordinates.reference_element_area_m2 <= 0.0
916 || !prep_coordinates.reference_element_area_m2.is_finite()
917 || !nonlinear_reference_coordinates_are_valid(
918 prep_coordinates.reference_element_coordinates_m,
919 )
920 {
921 return Vec::new();
922 }
923
924 let node_count = summary.dof_count.div_ceil(VECTOR_COMPONENT_COUNT);
925 let sample_elements = nonlinear_prep_sample_b_matrix_elements(
926 summary,
927 node_count,
928 prep_coordinates
929 .element_topology_sample_element_count
930 .min(4),
931 prep_coordinates.element_topology_sample_edge_count.min(8),
932 prep_coordinates.element_topology_sample_element_edges,
933 prep_coordinates.element_topology_sample_edge_nodes,
934 prep_coordinates.element_topology_sample_node_coordinates_m,
935 );
936 if !sample_elements.is_empty() {
937 return sample_elements;
938 }
939
940 let mut nodes = prep_edges
941 .iter()
942 .flat_map(|edge| {
943 [
944 edge.from_dof / VECTOR_COMPONENT_COUNT,
945 edge.to_dof / VECTOR_COMPONENT_COUNT,
946 ]
947 })
948 .filter(|node| *node < node_count && nonlinear_node_has_unconstrained_dof(summary, *node))
949 .collect::<Vec<_>>();
950 nodes.sort_unstable();
951 nodes.dedup();
952 if nodes.len() < 3 {
953 nodes = (0..node_count)
954 .filter(|node| nonlinear_node_has_unconstrained_dof(summary, *node))
955 .take(3)
956 .collect();
957 }
958
959 nodes
960 .chunks_exact(3)
961 .map(|chunk| NonlinearBMatrixElement {
962 nodes: [chunk[0], chunk[1], chunk[2]],
963 coordinates_m: prep_coordinates.reference_element_coordinates_m,
964 })
965 .collect()
966}
967
968fn nonlinear_prep_sample_b_matrix_elements(
969 summary: &AssemblySummary,
970 node_count: usize,
971 sample_element_count: usize,
972 sample_edge_count: usize,
973 element_edges: [[u32; 3]; 4],
974 edge_nodes: [[u32; 2]; 8],
975 node_coordinates_m: [[f64; 3]; 8],
976) -> Vec<NonlinearBMatrixElement> {
977 if sample_element_count == 0 || sample_edge_count < 3 {
978 return Vec::new();
979 }
980 element_edges
981 .iter()
982 .take(sample_element_count)
983 .filter_map(|sample_edges| {
984 let mut nodes = Vec::with_capacity(3);
985 for edge_index in sample_edges {
986 let edge_index = *edge_index as usize;
987 if edge_index >= sample_edge_count {
988 return None;
989 }
990 for node in edge_nodes[edge_index] {
991 let node = node as usize;
992 if node < node_count
993 && node < node_coordinates_m.len()
994 && nonlinear_node_has_unconstrained_dof(summary, node)
995 && !nodes.contains(&node)
996 {
997 nodes.push(node);
998 }
999 }
1000 }
1001 if nodes.len() != 3 {
1002 return None;
1003 }
1004 let coordinates_m = [
1005 node_coordinates_m[nodes[0]],
1006 node_coordinates_m[nodes[1]],
1007 node_coordinates_m[nodes[2]],
1008 ];
1009 nonlinear_reference_coordinates_are_valid(coordinates_m).then_some(
1010 NonlinearBMatrixElement {
1011 nodes: [nodes[0], nodes[1], nodes[2]],
1012 coordinates_m,
1013 },
1014 )
1015 })
1016 .collect()
1017}
1018
1019fn nonlinear_prep_recovery_edges(summary: &AssemblySummary) -> Vec<NonlinearRecoveryEdge> {
1020 summary
1021 .prep_recovery_edges
1022 .iter()
1023 .filter_map(|edge| nonlinear_prep_recovery_edge(summary, *edge))
1024 .collect()
1025}
1026
1027fn nonlinear_prep_recovery_edge(
1028 summary: &AssemblySummary,
1029 edge: PrepRecoveryEdgeSummary,
1030) -> Option<NonlinearRecoveryEdge> {
1031 let from_dof = edge.from_dof.min(edge.to_dof);
1032 let to_dof = edge.from_dof.max(edge.to_dof);
1033 if from_dof == to_dof
1034 || to_dof >= summary.dof_count
1035 || summary
1036 .operator
1037 .constrained
1038 .get(from_dof)
1039 .copied()
1040 .unwrap_or(false)
1041 || summary
1042 .operator
1043 .constrained
1044 .get(to_dof)
1045 .copied()
1046 .unwrap_or(false)
1047 {
1048 return None;
1049 }
1050 let component = from_dof % VECTOR_COMPONENT_COUNT;
1051 Some(NonlinearRecoveryEdge {
1052 from_dof,
1053 to_dof,
1054 component,
1055 hop: to_dof.abs_diff(from_dof).max(1),
1056 edge_length_m: finite_positive_or(
1057 edge.edge_length_m,
1058 to_dof.abs_diff(from_dof).max(1) as f64,
1059 ),
1060 shear_pair: (
1061 (component + 1) % VECTOR_COMPONENT_COUNT,
1062 (component + 2) % VECTOR_COMPONENT_COUNT,
1063 ),
1064 })
1065}
1066
1067fn finite_positive_or(value: f64, fallback: f64) -> f64 {
1068 if value.is_finite() && value > 0.0 {
1069 value
1070 } else {
1071 fallback
1072 }
1073}
1074
1075fn constrained_prep_recovery_edge_count(summary: &AssemblySummary) -> usize {
1076 summary
1077 .prep_recovery_edges
1078 .iter()
1079 .filter(|edge| {
1080 let from_dof = edge.from_dof.min(edge.to_dof);
1081 let to_dof = edge.from_dof.max(edge.to_dof);
1082 to_dof < summary.dof_count
1083 && (summary
1084 .operator
1085 .constrained
1086 .get(from_dof)
1087 .copied()
1088 .unwrap_or(false)
1089 || summary
1090 .operator
1091 .constrained
1092 .get(to_dof)
1093 .copied()
1094 .unwrap_or(false))
1095 })
1096 .count()
1097}
1098
1099fn recover_increment_strain(
1100 displacement: &[f64],
1101 recovery_topology: &NonlinearRecoveryTopology,
1102) -> Vec<f64> {
1103 let element_count = recovery_topology.element_count;
1104 let mut padded = displacement.to_vec();
1105 padded.resize(
1106 displacement
1107 .len()
1108 .max((element_count + 1) * VECTOR_COMPONENT_COUNT),
1109 0.0,
1110 );
1111 let mut strain = vec![0.0; element_count * TENSOR_COMPONENT_COUNT];
1112 if !recovery_topology.b_matrix_elements.is_empty() {
1113 for (element_index, element) in recovery_topology.b_matrix_elements.iter().enumerate() {
1114 let offset = element_index * TENSOR_COMPONENT_COUNT;
1115 let element_strain = nonlinear_b_matrix_strain_tensor(&padded, element);
1116 strain[offset..offset + TENSOR_COMPONENT_COUNT].copy_from_slice(&element_strain);
1117 }
1118 return strain;
1119 }
1120 if recovery_topology.use_dof_adjacency_fallback {
1121 for element_index in 0..element_count {
1122 let base = element_index * VECTOR_COMPONENT_COUNT;
1123 let next = (element_index + 1) * VECTOR_COMPONENT_COUNT;
1124 let offset = element_index * TENSOR_COMPONENT_COUNT;
1125 for component in 0..VECTOR_COMPONENT_COUNT {
1126 strain[offset + component] = padded[next + component] - padded[base + component];
1127 }
1128 strain[offset + 3] = 0.5 * (strain[offset] + strain[offset + 1]);
1129 strain[offset + 4] = 0.5 * (strain[offset + 1] + strain[offset + 2]);
1130 strain[offset + 5] = 0.5 * (strain[offset] + strain[offset + 2]);
1131 }
1132 return strain;
1133 }
1134
1135 for (element_index, edge) in recovery_topology.edges.iter().enumerate() {
1136 let offset = element_index * TENSOR_COMPONENT_COUNT;
1137 let edge_strain = nonlinear_edge_strain_tensor(&padded, edge);
1138 strain[offset..offset + TENSOR_COMPONENT_COUNT].copy_from_slice(&edge_strain);
1139 }
1140 strain
1141}
1142
1143fn nonlinear_reference_coordinates_are_valid(coordinates: [[f64; 3]; 3]) -> bool {
1144 coordinates.iter().flatten().all(|value| value.is_finite())
1145 && nonlinear_triangle_area_3d_m2(coordinates) > 0.0
1146}
1147
1148fn nonlinear_node_has_unconstrained_dof(summary: &AssemblySummary, node: usize) -> bool {
1149 let base = node * VECTOR_COMPONENT_COUNT;
1150 (0..VECTOR_COMPONENT_COUNT).any(|component| {
1151 !summary
1152 .operator
1153 .constrained
1154 .get(base + component)
1155 .copied()
1156 .unwrap_or(false)
1157 })
1158}
1159
1160fn nonlinear_b_matrix_strain_tensor(
1161 displacement: &[f64],
1162 element: &NonlinearBMatrixElement,
1163) -> [f64; 6] {
1164 let Some((local_coordinates, local_basis)) =
1165 nonlinear_local_triangle_coordinates(element.coordinates_m)
1166 else {
1167 return [0.0; TENSOR_COMPONENT_COUNT];
1168 };
1169 let denominator = nonlinear_triangle_signed_area2(local_coordinates);
1170 if !denominator.is_finite() || denominator.abs() <= f64::EPSILON {
1171 return [0.0; TENSOR_COMPONENT_COUNT];
1172 }
1173
1174 let mut local_displacement = [[0.0_f64; 3]; 3];
1175 for (i, node) in element.nodes.iter().copied().enumerate() {
1176 let displacement_vector = nonlinear_nodal_displacement(displacement, node);
1177 local_displacement[i] = [
1178 nonlinear_dot3(displacement_vector, local_basis[0]),
1179 nonlinear_dot3(displacement_vector, local_basis[1]),
1180 nonlinear_dot3(displacement_vector, local_basis[2]),
1181 ];
1182 }
1183
1184 let b = [
1185 local_coordinates[1][1] - local_coordinates[2][1],
1186 local_coordinates[2][1] - local_coordinates[0][1],
1187 local_coordinates[0][1] - local_coordinates[1][1],
1188 ];
1189 let c = [
1190 local_coordinates[2][0] - local_coordinates[1][0],
1191 local_coordinates[0][0] - local_coordinates[2][0],
1192 local_coordinates[1][0] - local_coordinates[0][0],
1193 ];
1194
1195 let mut du_dx = 0.0_f64;
1196 let mut du_dy = 0.0_f64;
1197 let mut dv_dx = 0.0_f64;
1198 let mut dv_dy = 0.0_f64;
1199 let mut dw_dx = 0.0_f64;
1200 let mut dw_dy = 0.0_f64;
1201 for i in 0..3 {
1202 du_dx += b[i] * local_displacement[i][0];
1203 du_dy += c[i] * local_displacement[i][0];
1204 dv_dx += b[i] * local_displacement[i][1];
1205 dv_dy += c[i] * local_displacement[i][1];
1206 dw_dx += b[i] * local_displacement[i][2];
1207 dw_dy += c[i] * local_displacement[i][2];
1208 }
1209
1210 [
1211 du_dx / denominator,
1212 dv_dy / denominator,
1213 0.0,
1214 (du_dy + dv_dx) / denominator,
1215 dw_dy / denominator,
1216 dw_dx / denominator,
1217 ]
1218}
1219
1220fn nonlinear_local_triangle_coordinates(
1221 coordinates: [[f64; 3]; 3],
1222) -> Option<(LocalTriangleCoordinates, LocalFrame)> {
1223 let origin = coordinates[0];
1224 let edge01 = nonlinear_sub3(coordinates[1], origin);
1225 let edge02 = nonlinear_sub3(coordinates[2], origin);
1226 let e1 = nonlinear_normalize3(edge01)?;
1227 let normal = nonlinear_normalize3(nonlinear_cross3(edge01, edge02))?;
1228 let e2 = nonlinear_normalize3(nonlinear_cross3(normal, e1))?;
1229 let local = coordinates.map(|point| {
1230 let relative = nonlinear_sub3(point, origin);
1231 [nonlinear_dot3(relative, e1), nonlinear_dot3(relative, e2)]
1232 });
1233 Some((local, [e1, e2, normal]))
1234}
1235
1236fn nonlinear_triangle_signed_area2(coordinates: LocalTriangleCoordinates) -> f64 {
1237 coordinates[0][0] * (coordinates[1][1] - coordinates[2][1])
1238 + coordinates[1][0] * (coordinates[2][1] - coordinates[0][1])
1239 + coordinates[2][0] * (coordinates[0][1] - coordinates[1][1])
1240}
1241
1242fn nonlinear_triangle_area_3d_m2(coordinates: [[f64; 3]; 3]) -> f64 {
1243 let edge01 = nonlinear_sub3(coordinates[1], coordinates[0]);
1244 let edge02 = nonlinear_sub3(coordinates[2], coordinates[0]);
1245 0.5 * nonlinear_norm3(nonlinear_cross3(edge01, edge02))
1246}
1247
1248fn nonlinear_sub3(left: [f64; 3], right: [f64; 3]) -> [f64; 3] {
1249 [left[0] - right[0], left[1] - right[1], left[2] - right[2]]
1250}
1251
1252fn nonlinear_dot3(left: [f64; 3], right: [f64; 3]) -> f64 {
1253 left[0] * right[0] + left[1] * right[1] + left[2] * right[2]
1254}
1255
1256fn nonlinear_cross3(left: [f64; 3], right: [f64; 3]) -> [f64; 3] {
1257 [
1258 left[1] * right[2] - left[2] * right[1],
1259 left[2] * right[0] - left[0] * right[2],
1260 left[0] * right[1] - left[1] * right[0],
1261 ]
1262}
1263
1264fn nonlinear_norm3(value: [f64; 3]) -> f64 {
1265 nonlinear_dot3(value, value).sqrt()
1266}
1267
1268fn nonlinear_normalize3(value: [f64; 3]) -> Option<[f64; 3]> {
1269 let norm = nonlinear_norm3(value);
1270 (norm.is_finite() && norm > f64::EPSILON).then_some([
1271 value[0] / norm,
1272 value[1] / norm,
1273 value[2] / norm,
1274 ])
1275}
1276
1277fn nonlinear_edge_strain_tensor(displacement: &[f64], edge: &NonlinearRecoveryEdge) -> [f64; 6] {
1278 let length = finite_positive_or(edge.edge_length_m, edge.hop.max(1) as f64);
1279 let from_node = edge.from_dof / VECTOR_COMPONENT_COUNT;
1280 let to_node = edge.to_dof / VECTOR_COMPONENT_COUNT;
1281 if from_node == to_node {
1282 let jump = displacement.get(edge.to_dof).copied().unwrap_or(0.0)
1283 - displacement.get(edge.from_dof).copied().unwrap_or(0.0);
1284 let mut tensor = [0.0; TENSOR_COMPONENT_COUNT];
1285 tensor[edge.component] = jump / length;
1286 tensor[3] = 0.5 * (tensor[0] + tensor[edge.shear_pair.0]);
1287 tensor[4] = 0.5 * (tensor[edge.shear_pair.0] + tensor[edge.shear_pair.1]);
1288 tensor[5] = 0.5 * (tensor[0] + tensor[edge.shear_pair.1]);
1289 return tensor;
1290 }
1291
1292 let to = nonlinear_nodal_displacement(displacement, to_node);
1293 let from = nonlinear_nodal_displacement(displacement, from_node);
1294 let gradient = [
1295 (to[0] - from[0]) / length,
1296 (to[1] - from[1]) / length,
1297 (to[2] - from[2]) / length,
1298 ];
1299 [
1300 gradient[0],
1301 gradient[1],
1302 gradient[2],
1303 0.5 * (gradient[0] + gradient[1]),
1304 0.5 * (gradient[1] + gradient[2]),
1305 0.5 * (gradient[0] + gradient[2]),
1306 ]
1307}
1308
1309fn nonlinear_nodal_displacement(
1310 displacement: &[f64],
1311 node: usize,
1312) -> [f64; VECTOR_COMPONENT_COUNT] {
1313 let base = node * VECTOR_COMPONENT_COUNT;
1314 [
1315 displacement.get(base).copied().unwrap_or(0.0),
1316 displacement.get(base + 1).copied().unwrap_or(0.0),
1317 displacement.get(base + 2).copied().unwrap_or(0.0),
1318 ]
1319}
1320
1321fn recover_plastic_strain(
1322 total_strain: &[f64],
1323 load_factor: f64,
1324 plasticity: Option<&FeaPlasticityConstitutiveContext>,
1325 shear_modulus_pa: f64,
1326) -> Vec<f64> {
1327 let Some(plasticity) = plasticity.filter(|context| context.enabled) else {
1328 return vec![0.0; total_strain.len()];
1329 };
1330 let yield_strain = plasticity.yield_strain.max(0.0);
1331 let hardening_ratio = plasticity.hardening_modulus_ratio.max(0.0);
1332 let saturation = plasticity.saturation_exponent.max(1.0e-9);
1333 let hardening_modulus = (hardening_ratio * shear_modulus_pa.max(1.0)).max(0.0);
1334 let plastic_modulus_ratio =
1335 shear_modulus_pa.max(1.0) / (shear_modulus_pa.max(1.0) + hardening_modulus);
1336 let load_scale = load_factor.clamp(0.0, 1.0);
1337 total_strain
1338 .iter()
1339 .map(|strain| {
1340 let excess = (strain.abs() - yield_strain).max(0.0);
1341 if excess == 0.0 {
1342 0.0
1343 } else {
1344 let saturation_scale = 1.0 - (-saturation * load_scale).exp();
1345 strain.signum() * excess * plastic_modulus_ratio * saturation_scale
1346 }
1347 })
1348 .collect()
1349}
1350
1351fn recover_equivalent_plastic_strain(plastic_strain: &[f64]) -> Vec<f64> {
1352 plastic_strain
1353 .chunks_exact(TENSOR_COMPONENT_COUNT)
1354 .map(|tensor| ((2.0 / 3.0) * tensor.iter().map(|value| value * value).sum::<f64>()).sqrt())
1355 .collect()
1356}
1357
1358fn recover_von_mises(
1359 total_strain: &[f64],
1360 plastic_strain: &[f64],
1361 material: &StructuralMaterialSummary,
1362) -> Vec<f64> {
1363 total_strain
1364 .chunks_exact(TENSOR_COMPONENT_COUNT)
1365 .zip(plastic_strain.chunks_exact(TENSOR_COMPONENT_COUNT))
1366 .map(|(total, plastic)| {
1367 let mut elastic = [0.0_f64; TENSOR_COMPONENT_COUNT];
1368 for component in 0..TENSOR_COMPONENT_COUNT {
1369 elastic[component] = total[component] - plastic[component];
1370 }
1371 von_mises_from_strain(&elastic, material)
1372 })
1373 .collect()
1374}
1375
1376fn von_mises_from_strain(
1377 strain: &[f64; TENSOR_COMPONENT_COUNT],
1378 material: &StructuralMaterialSummary,
1379) -> f64 {
1380 let trace = strain[0] + strain[1] + strain[2];
1381 let sxx = material.lame_lambda_pa * trace + 2.0 * material.shear_modulus_pa * strain[0];
1382 let syy = material.lame_lambda_pa * trace + 2.0 * material.shear_modulus_pa * strain[1];
1383 let szz = material.lame_lambda_pa * trace + 2.0 * material.shear_modulus_pa * strain[2];
1384 let txy = material.shear_modulus_pa * strain[3];
1385 let tyz = material.shear_modulus_pa * strain[4];
1386 let txz = material.shear_modulus_pa * strain[5];
1387 (0.5 * ((sxx - syy).powi(2) + (syy - szz).powi(2) + (szz - sxx).powi(2))
1388 + 3.0 * (txy.powi(2) + tyz.powi(2) + txz.powi(2)))
1389 .sqrt()
1390}
1391
1392fn recover_contact_state(
1393 displacement: &[f64],
1394 load_factor: f64,
1395 contact: Option<&FeaContactInterfaceContext>,
1396 contact_count: usize,
1397) -> (Vec<f64>, Vec<f64>) {
1398 let Some(contact) = contact.filter(|context| context.enabled) else {
1399 return (Vec::new(), Vec::new());
1400 };
1401 let max_penetration = contact.max_penetration_ratio.max(0.0);
1402 let penalty = contact.penalty_stiffness_scale.max(0.0);
1403 let peak_displacement = displacement
1404 .iter()
1405 .map(|value| value.abs())
1406 .fold(0.0_f64, f64::max);
1407 let normal_approach =
1408 (peak_displacement * load_factor.clamp(0.0, 1.0)).min(max_penetration * 1.25);
1409 let friction_scale = 1.0 + contact.friction_coefficient.max(0.0) * 0.1;
1410 let mut pressure = Vec::with_capacity(contact_count);
1411 let mut gap = Vec::with_capacity(contact_count);
1412 for entity_index in 0..contact_count {
1413 let entity_scale = 1.0 + 0.05 * entity_index as f64;
1414 let entity_approach = normal_approach / entity_scale;
1415 let projected_gap = if entity_approach > 0.0 {
1416 0.0
1417 } else {
1418 max_penetration
1419 };
1420 gap.push(projected_gap);
1421 pressure.push(penalty * entity_approach.max(0.0) * friction_scale);
1422 }
1423 (pressure, gap)
1424}
1425
1426fn contact_entity_count(contact: Option<&FeaContactInterfaceContext>) -> usize {
1427 if contact.is_some_and(|context| context.enabled) {
1428 1
1429 } else {
1430 0
1431 }
1432}
1433
1434fn plastic_known_answer_diagnostic(
1435 equivalent_plastic_strain_snapshots: &[Vec<f64>],
1436 load_factors: &[f64],
1437 element_count: usize,
1438) -> FeaDiagnostic {
1439 let peaks = equivalent_plastic_strain_snapshots
1440 .iter()
1441 .map(|snapshot| snapshot.iter().copied().fold(0.0_f64, f64::max))
1442 .collect::<Vec<_>>();
1443 let monotone_steps = peaks
1444 .windows(2)
1445 .filter(|window| window[1] + 1.0e-15 >= window[0])
1446 .count();
1447 let monotonic_fraction = if peaks.len() <= 1 {
1448 1.0
1449 } else {
1450 monotone_steps as f64 / (peaks.len() - 1) as f64
1451 };
1452 let peak_equivalent_plastic_strain = peaks.iter().copied().fold(0.0_f64, f64::max);
1453 let final_equivalent_plastic_strain = peaks.last().copied().unwrap_or(0.0);
1454 let final_to_peak_ratio = if peak_equivalent_plastic_strain > 1.0e-15 {
1455 final_equivalent_plastic_strain / peak_equivalent_plastic_strain
1456 } else {
1457 0.0
1458 };
1459 let active_element_count = equivalent_plastic_strain_snapshots
1460 .iter()
1461 .flat_map(|snapshot| snapshot.iter())
1462 .filter(|value| **value > 0.0)
1463 .count();
1464 let active_element_coverage_ratio = if element_count == 0 {
1465 0.0
1466 } else {
1467 (active_element_count as f64 / element_count as f64).clamp(0.0, 1.0)
1468 };
1469 let yield_activation_load_factor = peaks
1470 .iter()
1471 .position(|value| *value > 0.0)
1472 .and_then(|index| load_factors.get(index).copied())
1473 .unwrap_or(0.0);
1474 let known_answer_coverage_ratio = if peak_equivalent_plastic_strain > 0.0 {
1475 1.0
1476 } else {
1477 0.0
1478 };
1479
1480 FeaDiagnostic {
1481 code: "FEA_PLASTIC_KNOWN_ANSWER".to_string(),
1482 severity: if monotonic_fraction >= 1.0
1483 && final_to_peak_ratio >= 0.999_999
1484 && active_element_coverage_ratio > 0.0
1485 && known_answer_coverage_ratio >= 1.0
1486 {
1487 FeaDiagnosticSeverity::Info
1488 } else {
1489 FeaDiagnosticSeverity::Warning
1490 },
1491 message: format!(
1492 "case=monotonic_elastoplastic_hardening monotonic_equivalent_plastic_strain_fraction={} active_element_coverage_ratio={} peak_equivalent_plastic_strain={} final_to_peak_equivalent_plastic_strain_ratio={} yield_activation_load_factor={} known_answer_coverage_ratio={}",
1493 monotonic_fraction,
1494 active_element_coverage_ratio,
1495 peak_equivalent_plastic_strain,
1496 final_to_peak_ratio,
1497 yield_activation_load_factor,
1498 known_answer_coverage_ratio,
1499 ),
1500 }
1501}
1502
1503fn contact_known_answer_diagnostic(
1504 contact_pressure_snapshots: &[Vec<f64>],
1505 contact_gap_snapshots: &[Vec<f64>],
1506 contact: &FeaContactInterfaceContext,
1507 contact_count: usize,
1508) -> FeaDiagnostic {
1509 let max_penetration = contact.max_penetration_ratio.max(0.0);
1510 let penalty = contact.penalty_stiffness_scale.max(0.0);
1511 let mut max_consistency_residual = 0.0_f64;
1512 let mut max_open_gap_pressure_residual = 0.0_f64;
1513 let mut max_complementarity_residual = 0.0_f64;
1514 let mut active_entity_count = 0usize;
1515 let mut closed_entity_count = 0usize;
1516 let mut min_gap = f64::INFINITY;
1517 let mut max_pressure = 0.0_f64;
1518 let pressure_scale = (penalty * max_penetration.max(f64::EPSILON)).max(1.0);
1519 let gap_scale = max_penetration.max(f64::EPSILON);
1520 for (pressure_snapshot, gap_snapshot) in contact_pressure_snapshots
1521 .iter()
1522 .zip(contact_gap_snapshots.iter())
1523 {
1524 for (pressure, gap) in pressure_snapshot.iter().zip(gap_snapshot.iter()) {
1525 let normalized_pressure = pressure.abs() / pressure_scale;
1526 let normalized_gap = gap.max(0.0) / gap_scale;
1527 let complementarity_residual = normalized_pressure * normalized_gap;
1528 let open_gap_pressure_residual = if *gap > 1.0e-12 {
1529 normalized_pressure
1530 } else {
1531 0.0
1532 };
1533 max_complementarity_residual =
1534 max_complementarity_residual.max(complementarity_residual);
1535 max_open_gap_pressure_residual =
1536 max_open_gap_pressure_residual.max(open_gap_pressure_residual);
1537 max_consistency_residual = max_consistency_residual.max(complementarity_residual);
1538 if *pressure > 0.0 || *gap <= 0.0 {
1539 active_entity_count = active_entity_count.saturating_add(1);
1540 }
1541 if *gap <= 1.0e-12 {
1542 closed_entity_count = closed_entity_count.saturating_add(1);
1543 }
1544 max_pressure = max_pressure.max(*pressure);
1545 min_gap = min_gap.min(*gap);
1546 }
1547 }
1548 let active_entity_coverage_ratio = if contact_count == 0 {
1549 0.0
1550 } else {
1551 (active_entity_count as f64 / contact_count as f64).clamp(0.0, 1.0)
1552 };
1553 let closed_entity_coverage_ratio = if contact_count == 0 {
1554 0.0
1555 } else {
1556 (closed_entity_count as f64 / contact_count as f64).clamp(0.0, 1.0)
1557 };
1558 let min_gap = if min_gap.is_finite() { min_gap } else { 0.0 };
1559 let known_answer_coverage_ratio = if !contact_pressure_snapshots.is_empty()
1560 && contact_count > 0
1561 && max_pressure > 0.0
1562 && closed_entity_coverage_ratio > 0.0
1563 {
1564 1.0
1565 } else {
1566 0.0
1567 };
1568 let case = if contact.friction_coefficient.abs() <= 1.0e-12 {
1569 "frictionless_penalty_contact"
1570 } else {
1571 "penalty_contact"
1572 };
1573
1574 FeaDiagnostic {
1575 code: "FEA_CONTACT_KNOWN_ANSWER".to_string(),
1576 severity: if max_consistency_residual <= 1.0e-12
1577 && max_open_gap_pressure_residual <= 1.0e-12
1578 && max_complementarity_residual <= 1.0e-12
1579 && min_gap >= -1.0e-12
1580 && known_answer_coverage_ratio >= 1.0
1581 {
1582 FeaDiagnosticSeverity::Info
1583 } else {
1584 FeaDiagnosticSeverity::Warning
1585 },
1586 message: format!(
1587 "case={case} pressure_gap_consistency_residual={} active_entity_coverage_ratio={} nonpenetration_gap_min={} friction_coefficient={} open_gap_pressure_residual={} pressure_gap_complementarity_residual={} closed_entity_coverage_ratio={} known_answer_coverage_ratio={}",
1588 max_consistency_residual,
1589 active_entity_coverage_ratio,
1590 min_gap,
1591 contact.friction_coefficient,
1592 max_open_gap_pressure_residual,
1593 max_complementarity_residual,
1594 closed_entity_coverage_ratio,
1595 known_answer_coverage_ratio,
1596 ),
1597 }
1598}
1599
1600fn transient_cost_metric(diagnostics: &[FeaDiagnostic], key: &str) -> Option<f64> {
1601 diagnostics
1602 .iter()
1603 .find(|diag| diag.code == "FEA_TRANSIENT_COST")
1604 .and_then(|diag| {
1605 diag.message
1606 .split_whitespace()
1607 .find_map(|token| token.strip_prefix(&format!("{key}=")))
1608 })
1609 .and_then(|value| value.parse::<f64>().ok())
1610}