grid1d 0.3.6

#![deny(rustdoc::broken_intra_doc_links)]

//! Grid union operations for combining multiple grids with bidirectional mappings.
//!
//! This module provides functionality for creating a unified grid from two grids
//! that share the same domain. The unified grid contains all coordinate points
//! from both input grids, with mappings that allow tracking which intervals in
//! each original grid correspond to each unified interval.
//!
//! ## Key Types
//!
//! | Type | Description |
//! |------|-------------|
//! | [`Grid1DUnion`] | Unified grid with bidirectional mappings to original grids |
//! | [`ErrorsGrid1DUnion`] | Error types for union construction failures |
//!
//! ## Core Concept
//!
//! Given two grids A and B over the same domain, the union operation creates a new
//! grid whose coordinate points are the union of points from both grids:
//!
//! ```text
//! Grid A:    |-------|-------|--------|
//!            0      0.5      1      1.5
//!
//! Grid B:    |---|-----------|--------|
//!            0  0.3          1       1.5
//!
//! Union:     |---|---|-------|--------|
//!            0  0.3 0.5      1       1.5
//! ```
//!
//! Each interval in the unified grid maps back to exactly one interval in each
//! original grid, enabling data transfer and interpolation between grids.
//!
//! ## Basic Usage
//!
//! ```rust
//! use grid1d::{*, intervals::*, scalars::*};
//! use sorted_vec::partial::SortedSet;
//! use try_create::TryNew;
//!
//! let domain = IntervalClosed::new(0.0, 2.0);
//!
//! // Create two grids with different discretizations
//! let grid_a = Grid1D::uniform(domain.clone(), NumIntervals::try_new(4).unwrap());
//! let grid_b = Grid1D::<IntervalClosed<f64>>::try_from_sorted(
//!     SortedSet::from_unsorted(vec![0.0, 0.3, 0.8, 1.5, 2.0])
//! ).unwrap();
//!
//! // Create unified grid
//! let union = Grid1DUnion::try_new(&grid_a, &grid_b).unwrap();
//!
//! // Access mappings
//! for (unified_id, a_id, b_id) in union.iter_interval_mappings() {
//!     println!("Unified interval {} maps to A[{}] and B[{}]",
//!              *unified_id.as_ref(), *a_id.as_ref(), *b_id.as_ref());
//! }
//! ```
//!
//! ## Bidirectional Mappings
//!
//! The union maintains mappings from unified intervals to both original grids:
//!
//! ```rust
//! # use grid1d::{*, intervals::*, scalars::*};
//! # use try_create::TryNew;
//! # let domain = IntervalClosed::new(0.0, 2.0);
//! # let grid_a = Grid1D::uniform(domain.clone(), NumIntervals::try_new(4).unwrap());
//! # let grid_b = Grid1D::uniform(domain, NumIntervals::try_new(2).unwrap());
//! # let union = Grid1DUnion::try_new(&grid_a, &grid_b).unwrap();
//! // For any unified interval, find corresponding original intervals
//! let (a_id, b_id) = union.find_original_intervals(&IntervalId::new(0));
//!
//! // Access the mapping arrays directly
//! let mapping_a = union.mapping_to_grid_a();
//! let mapping_b = union.mapping_to_grid_b();
//! ```
//!
//! ## Use Cases
//!
//! ### Multi-Physics Coupling
//! When different physics simulations use different grid resolutions, the union
//! provides a common discretization for data exchange:
//!
//! ```rust
//! # use grid1d::{*, intervals::*, scalars::*};
//! fn couple_physics(
//!     flow_grid: &Grid1D<IntervalClosed<f64>>,
//!     chemistry_grid: &Grid1D<IntervalClosed<f64>>,
//!     flow_data: &[f64],
//!     chemistry_data: &[f64],
//! ) -> Vec<f64> {
//!     let union = Grid1DUnion::try_new(flow_grid, chemistry_grid).unwrap();
//!     let mut coupled = Vec::with_capacity(*union.num_refined_intervals().as_ref());
//!     
//!     for (_, flow_id, chem_id) in union.iter_interval_mappings() {
//!         // Combine data from both physics
//!         let combined = flow_data[*flow_id.as_ref()] + chemistry_data[*chem_id.as_ref()];
//!         coupled.push(combined);
//!     }
//!     
//!     coupled
//! }
//! ```
//!
//! ### Conservative Data Transfer
//! Transfer integral quantities while preserving conservation:
//!
//! ```rust
//! # use grid1d::{*, intervals::*, scalars::*};
//! # use try_create::TryNew;
//! # let domain = IntervalClosed::new(0.0, 1.0);
//! # let source_grid = Grid1D::uniform(domain.clone(), NumIntervals::try_new(2).unwrap());
//! # let target_grid = Grid1D::uniform(domain, NumIntervals::try_new(4).unwrap());
//! # let source_data = vec![1.0, 2.0];
//! let union = Grid1DUnion::try_new(&source_grid, &target_grid).unwrap();
//!
//! // Query intersections with physical subdomain
//! let subdomain = IntervalClosed::new(0.2, 0.8);
//! let intersections = union.find_intersections_with_mappings(&subdomain);
//! ```
//!
//! ## Requirements
//!
//! - **Domain equality**: Both input grids must have exactly the same domain
//! - **Valid grids**: Both inputs must be valid interval partitions
//!
//! ## Error Handling
//!
//! ```rust
//! use grid1d::{*, intervals::*, scalars::*};
//! use try_create::TryNew;
//!
//! let grid_a = Grid1D::uniform(
//!     IntervalClosed::new(0.0, 1.0),
//!     NumIntervals::try_new(4).unwrap()
//! );
//! let grid_b = Grid1D::uniform(
//!     IntervalClosed::new(0.0, 2.0), // Different domain!
//!     NumIntervals::try_new(4).unwrap()
//! );
//!
//! match Grid1DUnion::try_new(&grid_a, &grid_b) {
//!     Ok(_) => println!("Union created"),
//!     Err(ErrorsGrid1DUnion::DomainsMismatch { .. }) => {
//!         println!("Grids have different domains");
//!     }
//!     Err(e) => println!("Other error: {:?}", e),
//! }
//! ```
//!
//! ## Mathematical Properties
//!
//! - **Domain preservation**: `union.domain() == grid_a.domain() == grid_b.domain()`
//! - **Point superset**: Union coordinates ⊇ grid_a coordinates ∪ grid_b coordinates
//! - **Mapping completeness**: Every unified interval maps to valid intervals in both grids
//! - **Refinement property**: The union is a refinement of both input grids
//!
//! ## Performance
//!
//! | Operation | Complexity | Notes |
//! |-----------|------------|-------|
//! | `try_new` | O(n + m) | Merge-sort style coordinate union |
//! | `find_original_intervals` | O(1) | Direct array lookup |
//! | `iter_interval_mappings` | O(k) | k = unified interval count |
//! | `find_intersections_with_mappings` | O(log k + j) | j = intersection count |
//!
//! ## See Also
//!
//! - [`Grid1DRefinement`](super::refinement::Grid1DRefinement) - For subdividing single grids
//! - [`IntervalPartition`] - Trait implemented by unified grids
//! - [`Grid1D`] - Main grid type

use crate::{
    ErrorsGrid1D, Grid1D, IntervalPartition, SubIntervalInPartition,
    coords::Coords1D,
    intervals::bounded::{IntervalFinitePositiveLength, IntervalFinitePositiveLengthTrait},
    operations::refinement::{Grid1DNonUniformRefinement, Grid1DUniformRefinement},
    scalars::{IntervalId, NumIntervals, PositiveNumPoints1D},
    traits::{BuildIntervalInPartition, HasCoords1D, HasDomain1D},
};
use derive_more::Into;
use num_valid::{Constants, core::errors::capture_backtrace};
use serde::{Deserialize, Serialize};
use std::backtrace::Backtrace;
use thiserror::Error;

///  Comprehensive error types for 1D grid union construction and operations.
///
/// [`ErrorsGrid1DUnion`] provides specific error information for all failure modes
/// in grid union operations, enabling robust error handling and debugging when
/// combining multiple 1D grids into a unified partition.
///
/// Grid union operations require strict mathematical consistency between input grids,
/// and this error type captures the specific ways these operations can fail.
///
/// ## Error Philosophy
///
/// The error design follows the library's principle of **mathematical correctness first**:
/// - **Domain compatibility**: Input grids must have identical domains
/// - **Early failure detection**: Errors are caught at construction time, not during use
/// - **Detailed diagnostics**: Each error provides specific information about what went wrong
/// - **Type safety**: Generic over domain types to provide compile-time guarantees
///
/// ## Error Variants
///
/// ### Domain Compatibility Errors
/// - [`ErrorsGrid1DUnion::DomainsMismatch`]: Input grids have incompatible domains
///
/// ### Construction Errors  
/// - [`ErrorsGrid1DUnion::FromGrid1DConstructor`]: Error from underlying unified grid construction
///
/// ## Design Rationale
///
/// ### Why Domain Matching is Required
/// Grid union operations create a **unified partition** that must be mathematically
/// consistent with both input grids. This requires:
/// ```text
/// domain(grid_a) = domain(grid_b) = domain(unified_grid)
/// ```
///
/// If domains don't match, the resulting unified grid would be:
/// - **Mathematically undefined**: Cannot create consistent interval mappings
/// - **Physically meaningless**: Points from different physical domains cannot be combined
/// - **Computationally invalid**: Algorithms expecting domain consistency would fail
///
/// ### Error Propagation Strategy
/// The enum uses `#[from]` delegation for construction errors, ensuring that:
/// - **Root cause preservation**: Original error information is maintained
/// - **Error chain clarity**: Developers can trace errors back to their source
/// - **Consistent error handling**: Same patterns across the library
///
/// ## Usage Examples
///
/// ### Handling Domain Mismatch Errors
/// ```rust
/// use grid1d::{*, intervals::*, scalars::*};
/// use sorted_vec::partial::SortedSet;
///
/// let grid_a = Grid1D::<IntervalClosed<f64>>::try_from_sorted(
///     SortedSet::from_unsorted(vec![0.0, 0.5, 1.0])
/// ).unwrap();
/// let grid_b = Grid1D::<IntervalClosed<f64>>::try_from_sorted(
///     SortedSet::from_unsorted(vec![0.0, 1.0, 2.0])
/// ).unwrap();
///
/// match Grid1DUnion::try_new(&grid_a, &grid_b) {
///     Err(ErrorsGrid1DUnion::DomainsMismatch {
///         domain_grid_a,
///         domain_grid_b,
///         ..
///     }) => {
///         println!("Cannot union grids with different domains:");
///         println!("  Grid A domain: [{}, {}]",
///                  domain_grid_a.lower_bound_value(),
///                  domain_grid_a.upper_bound_value());
///         println!("  Grid B domain: [{}, {}]",
///                  domain_grid_b.lower_bound_value(),
///                  domain_grid_b.upper_bound_value());
///         
///         // Suggest potential fixes
///         println!("Potential solutions:");
///         println!("  1. Ensure both grids cover the same physical domain");
///         println!("  2. Trim grids to their common domain intersection");
///         println!("  3. Extend the smaller grid to match the larger domain");
///     }
///     Ok(union) => {
///         println!("Successfully created grid union");
///     }
///     Err(e) => {
///         println!("Other error: {}", e);
///     }
/// }
/// ```
///
/// ### Handling Propagated Construction Errors
/// ```rust
/// use grid1d::{*, intervals::*, scalars::*};
/// use sorted_vec::partial::SortedSet;
/// use try_create::TryNew;
///
/// fn safe_grid_union<D>(
///     grid_a: &Grid1D<D>,
///     grid_b: &Grid1D<D>,
/// ) -> Result<Grid1DUnion<D>, String>
/// where
///     D: BuildIntervalInPartition + std::fmt::Debug + Clone,
/// {
///     Grid1DUnion::try_new(grid_a, grid_b)
///         .map_err(|e| match e {
///             ErrorsGrid1DUnion::DomainsMismatch { domain_grid_a, domain_grid_b, .. } => {
///                 format!("Domain mismatch: {:?} vs {:?}", domain_grid_a, domain_grid_b)
///             }
///             ErrorsGrid1DUnion::FromGrid1DConstructor { source } => {
///                 format!("Failed to construct unified grid: {}", source)
///             }
///         })
/// }
///
/// // Usage
/// let domain = IntervalClosed::new(0.0, 1.0);
/// let grid_a = Grid1D::uniform(domain.clone(), NumIntervals::try_new(4).unwrap());
/// let grid_b = Grid1D::uniform(domain, NumIntervals::try_new(8).unwrap());
///
/// match safe_grid_union(&grid_a, &grid_b) {
///     Ok(union) => println!("Union successful: {} intervals",
///                          union.num_refined_intervals().as_ref()),
///     Err(msg) => println!("Union failed: {}", msg),
/// }
/// ```
///
/// ### Domain Validation Before Union
/// ```rust
/// use grid1d::{*, intervals::*, scalars::*};
/// use sorted_vec::partial::SortedSet;
///
/// fn validate_grids_for_union<D>(
///     grid_a: &Grid1D<D>,
///     grid_b: &Grid1D<D>,
/// ) -> Result<(), String>
/// where
///     D: BuildIntervalInPartition + PartialEq + std::fmt::Debug,
/// {
///     // Pre-validate domains to provide better error messages
///     if grid_a.domain() != grid_b.domain() {
///         return Err(format!(
///             "Grids have incompatible domains:\n  Grid A: {:?}\n  Grid B: {:?}",
///             grid_a.domain(),
///             grid_b.domain()
///         ));
///     }
///     
///     // Check for degenerate cases
///     if grid_a.num_intervals().as_ref() < &1 || grid_b.num_intervals().as_ref() < &1 {
///         return Err("Both grids must have at least one interval".to_string());
///     }
///     
///     Ok(())
/// }
///
/// fn create_validated_union<D>(
///     grid_a: &Grid1D<D>,
///     grid_b: &Grid1D<D>,
/// ) -> Result<Grid1DUnion<D>, String>
/// where
///     D: BuildIntervalInPartition + PartialEq + std::fmt::Debug + Clone,
/// {
///     // Validate first for better error messages
///     validate_grids_for_union(grid_a, grid_b)?;
///     
///     // Then attempt union
///     Grid1DUnion::try_new(grid_a, grid_b)
///         .map_err(|e| format!("Union failed despite validation: {}", e))
/// }
/// ```
///
/// ### Integration with Error Handling Frameworks
/// ```rust
/// use grid1d::{*, intervals::*, scalars::*};
/// use sorted_vec::partial::SortedSet;
/// use std::collections::HashMap;
///
/// #[derive(Debug)]
/// enum MultiGridError {
///     DomainMismatch(String),
///     ConstructionFailure(String),
///     InvalidInput(String),
/// }
///
/// impl From<ErrorsGrid1DUnion<IntervalClosed<f64>>> for MultiGridError {
///     fn from(error: ErrorsGrid1DUnion<IntervalClosed<f64>>) -> Self {
///         match error {
///             ErrorsGrid1DUnion::DomainsMismatch { domain_grid_a, domain_grid_b, .. } => {
///                 MultiGridError::DomainMismatch(format!(
///                     "Grid domains don't match: [{}, {}] vs [{}, {}]",
///                     domain_grid_a.lower_bound_value(),
///                     domain_grid_a.upper_bound_value(),
///                     domain_grid_b.lower_bound_value(),
///                     domain_grid_b.upper_bound_value()
///                 ))
///             }
///             ErrorsGrid1DUnion::FromGrid1DConstructor { source } => {
///                 MultiGridError::ConstructionFailure(format!(
///                     "Failed to build unified grid: {}", source
///                 ))
///             }
///         }
///     }
/// }
///
/// fn multi_grid_union(
///     grids: &[Grid1D<IntervalClosed<f64>>]
/// ) -> Result<Grid1D<IntervalClosed<f64>>, MultiGridError> {
///     if grids.len() < 2 {
///         return Err(MultiGridError::InvalidInput(
///             "Need at least 2 grids for union".to_string()
///         ));
///     }
///     
///     let (mut result, _, _) = Grid1DUnion::try_new(&grids[0], &grids[1])?.into_parts();
///     
///     for grid in &grids[2..] {
///         let union = Grid1DUnion::try_new(&result, grid)?;
///         (result, _, _) = union.into_parts();
///     }
///     
///     Ok(result)
/// }
/// ```
///
/// ## Error Variant Details
///
/// ### `DomainsMismatch`
///
/// **When it occurs**: When attempting to create a union of grids with different domains.
///
/// **Root cause**: The mathematical requirement that unified grids must have identical domains.
///
/// **Fields**:
/// - `domain_grid_a`: The domain of the first input grid
/// - `domain_grid_b`: The domain of the second input grid  
/// - `backtrace`: Stack trace for debugging (when backtrace capture is enabled)
///
/// **Common scenarios**:
/// - Grids created for different physical regions
/// - Grids with different coordinate systems or units
/// - Programming errors where wrong grids are passed to union
/// - Floating-point precision issues in domain bounds
///
/// **Resolution strategies**:
/// 1. **Domain trimming**: Reduce both grids to their common intersection
/// 2. **Domain extension**: Extend one grid to match the other's domain
/// 3. **Coordinate transformation**: Transform one grid to match the other's coordinate system
/// 4. **Input validation**: Check domain compatibility before attempting union
///
/// ### `FromGrid1DConstructor`
///
/// **When it occurs**: When the unified grid construction fails due to invalid coordinate sets.
///
/// **Root cause**: The merged coordinate set violates [`Grid1D`] construction requirements.
///
/// **Fields**:
/// - `source`: The underlying [`ErrorsGrid1D`] that caused the failure
///
/// **Common scenarios**:
/// - Merged coordinates create an invalid grid structure
/// - Numerical precision issues during coordinate merging
/// - Boundary coordinate mismatches after merging
/// - Memory allocation failures for large unified grids
///
/// **Resolution strategies**:
/// 1. **Input validation**: Verify individual grids before union
/// 2. **Precision management**: Use appropriate scalar types for coordinate precision
/// 3. **Memory management**: Ensure sufficient memory for large unions
/// 4. **Graceful degradation**: Fall back to simpler union strategies
///
/// ## Performance Implications
///
/// ### Error Construction Costs
/// - **Domain comparison**: O(1) for scalar bounds, O(k) for complex domains
/// - **Error message formatting**: Deferred until error is formatted (lazy evaluation)
/// - **Backtrace capture**: Significant overhead when enabled, zero when disabled
/// - **Memory allocation**: Minimal - errors store references and small data
///
/// ### Error Handling Best Practices
/// ```rust
/// use grid1d::{*, intervals::*, scalars::*};
/// use num_valid::core::errors::capture_backtrace;
///
/// // Efficient: Check domains before expensive operations
/// fn efficient_union<D>(
///     grid_a: &Grid1D<D>,
///     grid_b: &Grid1D<D>,
/// ) -> Result<Grid1DUnion<D>, ErrorsGrid1DUnion<D>>
/// where
///     D: BuildIntervalInPartition + PartialEq + Clone,
/// {
///     // Fast domain check before expensive union
///     if grid_a.domain() != grid_b.domain() {
///         return Err(ErrorsGrid1DUnion::DomainsMismatch {
///             domain_grid_a: grid_a.domain().clone(),
///             domain_grid_b: grid_b.domain().clone(),
///             backtrace: capture_backtrace(),
///         });
///     }
///     
///     Grid1DUnion::try_new(grid_a, grid_b)
/// }
///
/// // Memory-efficient: Use references where possible
/// fn analyze_union_error<D>(error: &ErrorsGrid1DUnion<D>) -> String
/// where
///     D: IntervalFinitePositiveLengthTrait,
/// {
///     match error {
///         ErrorsGrid1DUnion::DomainsMismatch { domain_grid_a, domain_grid_b, .. } => {
///             format!("Domain mismatch: {:?} ≠ {:?}", domain_grid_a, domain_grid_b)
///         }
///         ErrorsGrid1DUnion::FromGrid1DConstructor { source } => {
///             format!("Construction error: {}", source)
///         }
///     }
/// }
/// ```
///
/// ## Integration with the other numerical Ecosystem
///
/// ### Error Propagation Patterns
/// ```rust
/// use grid1d::{*, intervals::*, scalars::*};
///
/// // Pattern 1: Direct error propagation
/// fn create_multi_physics_union(
///     fluid_grid: &Grid1D<IntervalClosed<f64>>,
///     thermal_grid: &Grid1D<IntervalClosed<f64>>,
/// ) -> Result<Grid1DUnion<IntervalClosed<f64>>, ErrorsGrid1DUnion<IntervalClosed<f64>>> {
///     Grid1DUnion::try_new(fluid_grid, thermal_grid)
/// }
///
/// // Pattern 2: Error transformation
/// fn create_union_with_context(
///     grid_a: &Grid1D<IntervalClosed<f64>>,
///     grid_b: &Grid1D<IntervalClosed<f64>>,
///     context: &str,
/// ) -> Result<Grid1DUnion<IntervalClosed<f64>>, String> {
///     Grid1DUnion::try_new(grid_a, grid_b)
///         .map_err(|e| format!("Union failed in {}: {}", context, e))
/// }
///
/// // Pattern 3: Error aggregation
/// fn create_union_with_fallback(
///     primary_a: &Grid1D<IntervalClosed<f64>>,
///     primary_b: &Grid1D<IntervalClosed<f64>>,
///     fallback_a: &Grid1D<IntervalClosed<f64>>,
///     fallback_b: &Grid1D<IntervalClosed<f64>>,
/// ) -> Result<Grid1DUnion<IntervalClosed<f64>>, Vec<ErrorsGrid1DUnion<IntervalClosed<f64>>>> {
///     let mut errors = Vec::new();
///     
///     // Try primary union
///     match Grid1DUnion::try_new(primary_a, primary_b) {
///         Ok(union) => return Ok(union),
///         Err(e) => errors.push(e),
///     }
///     
///     // Try fallback union
///     match Grid1DUnion::try_new(fallback_a, fallback_b) {
///         Ok(union) => Ok(union),
///         Err(e) => {
///             errors.push(e);
///             Err(errors)
///         }
///     }
/// }
/// ```
///
/// ## Mathematical Properties and Guarantees
///
/// ### Error Correctness
/// ```text
/// ∀ grids A, B: Grid1DUnion::try_new(A, B) fails ⟺
///   domain(A) ≠ domain(B) ∨ unified_coords_invalid
/// ```
/// Errors occur if and only if there's a genuine mathematical incompatibility.
///
/// ### Error Determinism
/// ```text
/// ∀ grids A, B:
///   Grid1DUnion::try_new(A, B) = Grid1DUnion::try_new(A, B) (same error or success)
/// ```
/// Error conditions are deterministic and reproducible.
///
/// ### Error Completeness
/// ```text
/// ∀ error conditions: ∃ variant ∈ ErrorsGrid1DUnion that captures the condition
/// ```
/// All meaningful failure modes are represented by error variants.
///
/// ## See Also
///
/// - [`Grid1DUnion`]: The struct that uses this error type
/// - [`ErrorsGrid1D`]: Related error type for general grid construction
/// - [`Grid1D`]: The fundamental grid type used in unions
/// - [`IntervalPartition`]: The trait providing partition operations
/// - [`HasDomain1D`]: Core trait for domain access and comparison
/// - Error handling patterns in the [module-level documentation](crate)
#[derive(Debug, Error)]
pub enum ErrorsGrid1DUnion<Domain1D: IntervalFinitePositiveLengthTrait> {
    /// Error indicating that the two grids have incompatible domains.
    #[error("The 1D domains of the two grids do not match!")]
    DomainsMismatch {
        /// The domain of the first grid.
        domain_grid_a: Domain1D,
        /// The domain of the second grid.
        domain_grid_b: Domain1D,
        /// Captured backtrace for debugging.
        backtrace: Backtrace,
    },

    /// Error propagated from the underlying [`Grid1D`] construction during union.
    #[error("Error from the Grid1D constructor")]
    FromGrid1DConstructor {
        /// The underlying error from [`Grid1D`] construction.
        #[from]
        source: ErrorsGrid1D<Domain1D>,
    },
}

///  Result of computing the union of two 1D grids with interval mappings.
///
/// This struct contains the refined grid that represents the union of all points
/// from both input grids (defined over the same domain), along with bidirectional mappings that allow tracking
/// how intervals in the refined grid correspond to intervals in the original grids.
///
/// ## Mathematical Definition
///
/// Given two grids `G_A` and `G_B` with point sets `P_A` and `P_B`, the union creates:
/// ```text
/// P_refined = P_A ∪ P_B (sorted and deduplicated)
/// G_refined = Grid(domain, P_refined)
/// ```
///
/// The mapping functions satisfy:
/// ```text
/// map_A: IntervalId(G_refined) → IntervalId(G_A)
/// map_B: IntervalId(G_refined) → IntervalId(G_B)
/// ```
///
/// where `map_A[i]` gives the interval in `G_A` containing the midpoint of refined interval `i`.
///
/// ## Usage Patterns
///
/// ### Basic Access
/// ```rust
/// use grid1d::{*, intervals::*, scalars::*};
/// use sorted_vec::partial::SortedSet;
///
/// let grid_a = Grid1D::<IntervalClosed<f64>>::try_from_sorted(SortedSet::from_unsorted(vec![0.0, 1.0, 2.0])).unwrap();
/// let grid_b = Grid1D::<IntervalClosed<f64>>::try_from_sorted(SortedSet::from_unsorted(vec![0.0, 0.5, 2.0])).unwrap();
///
/// let union = Grid1DUnion::try_new(&grid_a, &grid_b).unwrap();
///
/// // Access the unified grid
/// println!("Unified grid: {:?}", union.unified_grid().coords());
///
/// // Find which original intervals correspond to a refined interval
/// let refined_id = IntervalId::new(1);
/// let original_a_id = union.mapping_to_grid_a()[1];
/// let original_b_id = union.mapping_to_grid_b()[1];
/// println!("Refined interval {} maps to A[{}], B[{}]",
///          refined_id.as_ref(), original_a_id.as_ref(), original_b_id.as_ref());
/// ```
///
/// ### Iteration Patterns
/// ```rust
/// use grid1d::{*, intervals::*, scalars::*};
/// use sorted_vec::partial::SortedSet;
///
/// # let grid_a = Grid1D::<IntervalClosed<f64>>::try_from_sorted(SortedSet::from_unsorted(vec![0.0, 1.0, 2.0])).unwrap();
/// # let grid_b = Grid1D::<IntervalClosed<f64>>::try_from_sorted(SortedSet::from_unsorted(vec![0.0, 0.5, 2.0])).unwrap();
/// # let union = Grid1DUnion::try_new(&grid_a, &grid_b).unwrap();
///
/// // Iterate over all refined intervals with their mappings
/// for (refined_id, a_id, b_id) in union.iter_interval_mappings() {
///     println!("Refined {}: A={}, B={}", refined_id.as_ref(), a_id.as_ref(), b_id.as_ref());
/// }
///
/// // Alternative: using the unified grid directly
/// for (refined_id, _) in union.unified_grid().iter_intervals() {
///     let (a_id, b_id) = union.find_original_intervals(&refined_id);
///     println!("Refined {}: A={}, B={}", refined_id.as_ref(), a_id.as_ref(), b_id.as_ref());
/// }
/// ```
///
/// ### Finite Element Assembly
/// ```rust
/// use grid1d::{*, intervals::*, scalars::*};
/// use sorted_vec::partial::SortedSet;
///
/// # let grid_a = Grid1D::<IntervalClosed<f64>>::try_from_sorted(SortedSet::from_unsorted(vec![0.0, 1.0, 2.0])).unwrap();
/// # let grid_b = Grid1D::<IntervalClosed<f64>>::try_from_sorted(SortedSet::from_unsorted(vec![0.0, 0.5, 2.0])).unwrap();
/// # let union = Grid1DUnion::try_new(&grid_a, &grid_b).unwrap();
/// # struct Element { data: f64 }
/// # let elements_a = vec![Element { data: 1.0 }, Element { data: 2.0 }];
/// # let elements_b = vec![Element { data: 3.0 }, Element { data: 4.0 }];
/// # fn assemble_element_contribution(_: &Element, _: &Element, _: &grid1d::intervals::SubIntervalInPartition<IntervalClosed<f64>>) {}
///
/// // Typical usage in finite element methods
/// for (refined_id, refined_interval) in union.unified_grid().iter_intervals() {
///     let (a_id, b_id) = union.find_original_intervals(&refined_id);
///     let element_a = &elements_a[*a_id.as_ref()];
///     let element_b = &elements_b[*b_id.as_ref()];
///     
///     // Perform element assembly over the refined interval
///     assemble_element_contribution(element_a, element_b, &refined_interval);
/// }
/// ```
///
/// ### Working with Domain Traits
/// ```rust
/// use grid1d::{*, intervals::*, scalars::*};
/// use sorted_vec::partial::SortedSet;
/// use std::ops::Deref;
///
/// # let grid_a = Grid1D::<IntervalClosed<f64>>::try_from_sorted(SortedSet::from_unsorted(vec![0.0, 1.0, 2.0])).unwrap();
/// # let grid_b = Grid1D::<IntervalClosed<f64>>::try_from_sorted(SortedSet::from_unsorted(vec![0.0, 0.5, 2.0])).unwrap();
/// # let union = Grid1DUnion::try_new(&grid_a, &grid_b).unwrap();
///
/// // Grid1DUnion implements HasDomain1D, HasCoords1D, and IntervalPartition
/// assert_eq!(union.domain(), grid_a.domain());
/// assert_eq!(union.domain(), grid_b.domain());
/// assert_eq!(union.coords().deref(), &[0.0, 0.5, 1.0, 2.0]);
/// assert_eq!(union.num_intervals().as_ref(), &3);
///
/// // Can be used with generic functions expecting IntervalPartition
/// fn analyze_partition<P: IntervalPartition>(partition: &P) {
///     println!("Partition has {} intervals", partition.num_intervals().as_ref());
///     for (id, interval) in partition.iter_intervals() {
///         println!("Interval {}: length = {}", id.as_ref(), partition.interval_length(&id));
///     }
/// }
///
/// analyze_partition(&union);
/// ```
///
/// ## Performance Characteristics
///
/// | Operation | Time Complexity | Space Complexity | Notes |
/// |-----------|----------------|------------------|--------|
/// | **Creation** | O(n+m) | O(n+m) | Using merge algorithm |
/// | **Grid access** | O(1) | O(1) | Direct reference |
/// | **Mapping lookup** | O(1) | O(1) | Vector indexing |
/// | **Iteration** | O(k) | O(1) | Where k is result size |
/// | **Point location** | O(log(n+m)) | O(1) | Binary search on unified grid |
///
/// Where `n` and `m` are the sizes of the input grids.
///
/// ## Memory Layout
///
/// - **Unified grid**: Contains merged coordinates and detected structure (uniform/non-uniform)
/// - **Mappings**: Two `Vec<IntervalId>` for O(1) lookup
/// - **Total overhead**: Minimal - only stores the essential mapping information
///
/// ## Integration with other numerical Ecosystem
///
/// `Grid1DUnion` implements the core domain traits:
/// - **[`HasDomain1D`]**: Access to the unified grid's domain
/// - **[`HasCoords1D`]**: Access to the unified grid's coordinates  
/// - **[`IntervalPartition`]**: Full partition interface on the unified grid
///
/// This means `Grid1DUnion` can be used anywhere these traits are expected,
/// providing seamless integration with the broader mathematical ecosystem.
///
/// ## Mathematical Properties
///
/// ### Completeness
/// The unified grid contains all unique points from both input grids:
/// ```text
/// ∀p ∈ P_A ∪ P_B ⟹ p ∈ P_refined
/// ```
///
/// ### Mapping Consistency
/// Each refined interval maps to exactly one interval in each original grid:
/// ```text
/// ∀i ∈ refined_intervals: map_A[i] ∈ [0, |G_A|-1] ∧ map_B[i] ∈ [0, |G_B|-1]
/// ```
///
/// ### Containment Property
/// The midpoint strategy ensures each refined interval is properly contained:
/// ```text
/// ∀i: midpoint(refined_interval[i]) ∈ original_interval[map_A[i]]
/// ```
///
/// ## Error Handling
///
/// [`Grid1DUnion::try_new()`] can fail with:
/// - **[`ErrorsGrid1DUnion::DomainsMismatch`]**: Input grids have different domains
/// - **[`ErrorsGrid1DUnion::FromGrid1DConstructor`]**: Unified grid construction fails
///
/// ### Domain Mismatch
/// ```rust
/// use grid1d::{*, intervals::*, scalars::*};
/// use sorted_vec::partial::SortedSet;
///
/// // Grids on different domains!
/// let grid_a = Grid1D::<IntervalClosed<f64>>::try_from_sorted(SortedSet::from_unsorted(vec![0.0, 1.0])).unwrap();
/// let grid_b = Grid1D::<IntervalClosed<f64>>::try_from_sorted(SortedSet::from_unsorted(vec![0.0, 2.0])).unwrap();
///
/// let err = Grid1DUnion::try_new(&grid_a, &grid_b).unwrap_err();
/// assert!(matches!(err, ErrorsGrid1DUnion::DomainsMismatch { domain_grid_a, domain_grid_b, .. }));
/// ```
///
/// ## See Also
///
/// - [`Grid1D`]: The fundamental 1D grid type
/// - [`IntervalPartition`]: The trait providing partition operations
/// - [`HasDomain1D`], [`HasCoords1D`]: Core domain access traits
#[derive(Debug, Clone, PartialEq, Into, Serialize, Deserialize)]
#[serde(bound(deserialize = "Domain1D: for<'d> serde::Deserialize<'d>, \
                  Domain1D::RealType: for<'d> serde::Deserialize<'d>"))]
pub struct Grid1DUnion<Domain1D: BuildIntervalInPartition> {
    /// The refined grid containing all unique points from both input grids
    unified_grid: Grid1D<Domain1D>,

    /// Mapping from refined grid intervals to the first input grid intervals
    mapping_to_grid_a: Vec<IntervalId>,

    /// Mapping from refined grid intervals to the second input grid intervals
    mapping_to_grid_b: Vec<IntervalId>,
}

impl<Domain1D: BuildIntervalInPartition> Grid1DUnion<Domain1D> {
    /// Computes the union of two 1D grids with interval mapping.
    ///
    /// This function creates a refined grid containing all unique points from both
    /// input grids, then computes bidirectional mappings between interval IDs.
    ///
    /// # Parameters
    ///
    /// - `grid1d_a`: First input grid
    /// - `grid1d_b`: Second input grid (must have same domain as `grid1d_a`)
    ///
    /// # Returns
    ///
    /// A [`Grid1DUnion`] struct containing the unified grid and interval mappings.
    pub fn try_new(
        grid1d_a: &Grid1D<Domain1D>,
        grid1d_b: &Grid1D<Domain1D>,
    ) -> Result<Self, ErrorsGrid1DUnion<Domain1D>> {
        let domain = grid1d_a.domain().clone();

        if &domain != grid1d_b.domain() {
            return Err(ErrorsGrid1DUnion::DomainsMismatch {
                domain_grid_a: domain,
                domain_grid_b: grid1d_b.domain().clone(),
                backtrace: capture_backtrace(),
            });
        }

        // Get sorted slices
        let coords_a = grid1d_a.coords();
        let coords_b = grid1d_b.coords();

        // Merge the two sorted vectors with deduplication
        let merged_coords = coords_a.union(coords_b);
        let unified_grid = Grid1D::try_from_coords(merged_coords)?;

        let compute_interval_map =
            |unified_grid: &Grid1D<Domain1D>, grid1d_coarse: &Grid1D<Domain1D>| {
                let n_intervals = *unified_grid.num_intervals().as_ref();
                let mut map = Vec::with_capacity(n_intervals);

                let one_div_2 = Domain1D::RealType::one_div_2();

                let v = unified_grid.coords();
                for i in 0..n_intervals {
                    // Find the midpoint of the current interval in the unified grid.
                    let midpoint_interval_unified_grid = (v[i].clone() + &v[i + 1]) * &one_div_2;

                    // Find which interval in the coarse grid contains this midpoint.
                    // Using the panicking version is safe here because the midpoint is guaranteed
                    // to be within the coarse grid's domain.
                    let coarse_interval_id =
                        grid1d_coarse.find_interval_id_of_point(&midpoint_interval_unified_grid);
                    map.push(coarse_interval_id);
                }

                map
            };

        let mapping_to_grid_a = compute_interval_map(&unified_grid, grid1d_a);
        let mapping_to_grid_b = compute_interval_map(&unified_grid, grid1d_b);

        Ok(Self {
            unified_grid,
            mapping_to_grid_a,
            mapping_to_grid_b,
        })
    }

    /// Returns a reference to the unified grid.
    ///
    /// The unified grid contains all unique points from both input grids,
    /// creating a refined partition that allows consistent operations
    /// across both original grids.
    ///
    /// # Example
    ///
    /// ```rust
    /// use grid1d::{*, intervals::*, scalars::*};
    /// use sorted_vec::partial::SortedSet;
    /// use std::ops::Deref;
    ///
    /// let grid_a = Grid1D::<IntervalClosed<f64>>::try_from_sorted(SortedSet::from_unsorted(vec![0.0, 1.0, 2.0])).unwrap();
    /// let grid_b = Grid1D::<IntervalClosed<f64>>::try_from_sorted(SortedSet::from_unsorted(vec![0.0, 0.5, 2.0])).unwrap();
    ///
    /// let union = Grid1DUnion::try_new(&grid_a, &grid_b).unwrap();
    /// let refined_coords = union.unified_grid().coords();
    /// assert_eq!(refined_coords.deref(), &[0.0, 0.5, 1.0, 2.0]);
    /// ```
    pub fn unified_grid(&self) -> &Grid1D<Domain1D> {
        &self.unified_grid
    }

    /// Returns a reference to the mapping from refined grid intervals to the first input grid.
    ///
    /// This mapping allows you to determine which interval in the first input grid
    /// contains the midpoint of each refined interval.
    ///
    /// # Example
    ///
    /// ```rust
    /// use grid1d::{*, intervals::*, scalars::*};
    /// use sorted_vec::partial::SortedSet;
    ///
    /// let grid_a = Grid1D::<IntervalClosed<f64>>::try_from_sorted(SortedSet::from_unsorted(vec![0.0, 1.0, 2.0])).unwrap();
    /// let grid_b = Grid1D::<IntervalClosed<f64>>::try_from_sorted(SortedSet::from_unsorted(vec![0.0, 0.5, 2.0])).unwrap();
    ///
    /// let union = Grid1DUnion::try_new(&grid_a, &grid_b).unwrap();
    /// let mapping_a = union.mapping_to_grid_a();
    ///
    /// // Refined interval [0.0, 0.5] maps to grid_a interval [0.0, 1.0]
    /// assert_eq!(mapping_a[0], IntervalId::new(0));
    /// ```
    pub fn mapping_to_grid_a(&self) -> &Vec<IntervalId> {
        &self.mapping_to_grid_a
    }

    /// Returns a reference to the mapping from refined grid intervals to the second input grid.
    ///
    /// This mapping allows you to determine which interval in the second input grid
    /// contains the midpoint of each refined interval.
    ///
    /// # Example
    ///
    /// ```rust
    /// use grid1d::{*, intervals::*, scalars::*};
    /// use sorted_vec::partial::SortedSet;
    ///
    /// let grid_a = Grid1D::<IntervalClosed<f64>>::try_from_sorted(SortedSet::from_unsorted(vec![0.0, 1.0, 2.0])).unwrap();
    /// let grid_b = Grid1D::<IntervalClosed<f64>>::try_from_sorted(SortedSet::from_unsorted(vec![0.0, 0.5, 2.0])).unwrap();
    ///
    /// let union = Grid1DUnion::try_new(&grid_a, &grid_b).unwrap();
    /// let mapping_b = union.mapping_to_grid_b();
    ///
    /// // Refined interval [0.0, 0.5] maps to grid_b interval [0.0, 0.5]
    /// assert_eq!(mapping_b[0], IntervalId::new(0));
    /// ```
    pub fn mapping_to_grid_b(&self) -> &Vec<IntervalId> {
        &self.mapping_to_grid_b
    }

    /// Returns the number of intervals in the unified grid.
    ///
    /// # Example
    ///
    /// ```rust
    /// use grid1d::{*, intervals::*, scalars::*};
    /// use sorted_vec::partial::SortedSet;
    ///
    /// let grid_a = Grid1D::<IntervalClosed<f64>>::try_from_sorted(SortedSet::from_unsorted(vec![0.0, 1.0, 2.0])).unwrap();
    /// let grid_b = Grid1D::<IntervalClosed<f64>>::try_from_sorted(SortedSet::from_unsorted(vec![0.0, 0.5, 2.0])).unwrap();
    ///
    /// let union = Grid1DUnion::try_new(&grid_a, &grid_b).unwrap();
    /// assert_eq!(union.num_refined_intervals().as_ref(), &3);
    /// ```
    pub fn num_refined_intervals(&self) -> NumIntervals {
        self.unified_grid.num_intervals()
    }

    /// Finds the corresponding interval IDs in both original grids for a given refined interval.
    ///
    /// This is a convenience method that combines both mapping lookups.
    ///
    /// # Parameters
    ///
    /// - `refined_id`: The interval ID in the unified grid
    ///
    /// # Returns
    ///
    /// A tuple `(grid_a_id, grid_b_id)` containing the corresponding interval IDs
    /// in both original grids.
    ///
    /// # Panics
    ///
    /// Panics if the `refined_id` is not found in the mappings (which should never
    /// happen if the struct is constructed correctly).
    ///
    /// # Example
    ///
    /// ```rust
    /// use grid1d::{*, intervals::*, scalars::*};
    /// use sorted_vec::partial::SortedSet;
    ///
    /// let grid_a = Grid1D::<IntervalClosed<f64>>::try_from_sorted(SortedSet::from_unsorted(vec![0.0, 1.0, 2.0])).unwrap();
    /// let grid_b = Grid1D::<IntervalClosed<f64>>::try_from_sorted(SortedSet::from_unsorted(vec![0.0, 0.5, 2.0])).unwrap();
    ///
    /// let union = Grid1DUnion::try_new(&grid_a, &grid_b).unwrap();
    /// let (a_id, b_id) = union.find_original_intervals(&IntervalId::new(1));
    ///
    /// println!("Refined interval 1 maps to: grid_a[{}], grid_b[{}]", a_id.as_ref(), b_id.as_ref());
    /// ```
    pub fn find_original_intervals(&self, refined_id: &IntervalId) -> (IntervalId, IntervalId) {
        let refined_id = *refined_id.as_ref();
        let grid_a_id = self.mapping_to_grid_a[refined_id];
        let grid_b_id = self.mapping_to_grid_b[refined_id];
        (grid_a_id, grid_b_id)
    }

    /// Returns an iterator over all refined intervals with their corresponding original interval IDs.
    ///
    /// This provides a convenient way to iterate over all mappings simultaneously.
    ///
    /// # Returns
    ///
    /// An iterator yielding tuples of `(refined_id, grid_a_id, grid_b_id)`.
    ///
    /// # Example
    ///
    /// ```rust
    /// use grid1d::{*, intervals::*, scalars::*};
    /// use sorted_vec::partial::SortedSet;
    ///
    /// let grid_a = Grid1D::<IntervalClosed<f64>>::try_from_sorted(SortedSet::from_unsorted(vec![0.0, 1.0, 2.0])).unwrap();
    /// let grid_b = Grid1D::<IntervalClosed<f64>>::try_from_sorted(SortedSet::from_unsorted(vec![0.0, 0.5, 2.0])).unwrap();
    ///
    /// let union = Grid1DUnion::try_new(&grid_a, &grid_b).unwrap();
    ///
    /// for (refined_id, a_id, b_id) in union.iter_interval_mappings() {
    ///     println!("Refined {}: A={}, B={}", *refined_id.as_ref(), *a_id.as_ref(), *b_id.as_ref());
    /// }
    /// ```
    pub fn iter_interval_mappings(
        &self,
    ) -> impl Iterator<Item = (IntervalId, IntervalId, IntervalId)> + '_ {
        (0..*self.num_refined_intervals().as_ref()).map(move |i| {
            let refined_id = IntervalId::new(i);
            let (a_id, b_id) = self.find_original_intervals(&refined_id);
            (refined_id, a_id, b_id)
        })
    }

    /// Returns an iterator over refined intervals with their physical domains.
    ///
    /// This provides a convenient way to iterate over refined intervals
    /// along with their actual geometric domains.
    ///
    /// # Example
    ///
    /// ```rust
    /// use grid1d::{*, intervals::*, scalars::*};
    /// use sorted_vec::partial::SortedSet;
    ///
    /// let grid_a = Grid1D::<IntervalClosed<f64>>::try_from_sorted(SortedSet::from_unsorted(vec![0.0, 1.0, 2.0])).unwrap();
    /// let grid_b = Grid1D::<IntervalClosed<f64>>::try_from_sorted(SortedSet::from_unsorted(vec![0.0, 0.5, 2.0])).unwrap();
    ///
    /// let union = Grid1DUnion::try_new(&grid_a, &grid_b).unwrap();
    ///
    /// for (refined_id, interval, a_id, b_id) in union.iter_intervals_with_mappings() {
    ///     println!("Refined interval {}: {:?} -> A[{}], B[{}]",
    ///              refined_id.as_ref(), interval, a_id.as_ref(), b_id.as_ref());
    /// }
    /// ```
    pub fn iter_intervals_with_mappings(
        &self,
    ) -> impl Iterator<
        Item = (
            IntervalId,
            SubIntervalInPartition<Domain1D>,
            IntervalId,
            IntervalId,
        ),
    > + '_ {
        self.unified_grid()
            .iter_intervals()
            .map(move |(refined_id, interval)| {
                let (a_id, b_id) = self.find_original_intervals(&refined_id);
                (refined_id, interval, a_id, b_id)
            })
    }

    /// Finds all refined intervals that intersect with the given domain.
    ///
    /// This combines the functionality of `intervals_in_intersection` with
    /// the mapping information, providing a complete view of how the intersection
    /// relates to both original grids.
    ///
    /// # Example
    ///
    /// ```rust
    /// use grid1d::{*, intervals::*, scalars::*};
    /// use sorted_vec::partial::SortedSet;
    ///
    /// let grid_a = Grid1D::<IntervalClosed<f64>>::try_from_sorted(SortedSet::from_unsorted(vec![0.0, 1.0, 2.0])).unwrap();
    /// let grid_b = Grid1D::<IntervalClosed<f64>>::try_from_sorted(SortedSet::from_unsorted(vec![0.0, 0.5, 2.0])).unwrap();
    ///
    /// let union = Grid1DUnion::try_new(&grid_a, &grid_b).unwrap();
    /// let subdomain = IntervalClosed::new(0.3, 1.2);
    ///
    /// let intersections = union.find_intersections_with_mappings(&subdomain);
    /// for (refined_id, intersection, a_id, b_id) in intersections {
    ///     println!("Intersection in refined[{}]: {:?} -> A[{}], B[{}]",
    ///              refined_id.as_ref(), intersection, a_id.as_ref(), b_id.as_ref());
    /// }
    /// ```
    pub fn find_intersections_with_mappings<
        IntervalType: IntervalFinitePositiveLengthTrait<RealType = Domain1D::RealType>,
    >(
        &self,
        domain_in: &IntervalType,
    ) -> Vec<(
        IntervalId,
        IntervalFinitePositiveLength<Domain1D::RealType>,
        IntervalId,
        IntervalId,
    )> {
        self.unified_grid()
            .intervals_in_intersection(domain_in)
            .into_iter()
            .map(|(refined_id, intersection)| {
                let (a_id, b_id) = self.find_original_intervals(&refined_id);
                (refined_id, intersection, a_id, b_id)
            })
            .collect()
    }

    /// Deconstructs the `Grid1DUnion` into its component parts.
    ///
    /// This method consumes the struct and returns the individual components,
    /// which can be useful when you need to work with the parts separately
    /// or integrate with legacy code that expects the tuple format.
    ///
    /// # Returns
    ///
    /// A tuple containing:
    /// 1. The unified grid
    /// 2. Mapping from refined grid to first input grid
    /// 3. Mapping from refined grid to second input grid
    ///
    /// # Example
    ///
    /// ```rust
    /// use grid1d::{*, intervals::*, scalars::*};
    /// use sorted_vec::partial::SortedSet;
    ///
    /// let grid_a = Grid1D::<IntervalClosed<f64>>::try_from_sorted(SortedSet::from_unsorted(vec![0.0, 1.0, 2.0])).unwrap();
    /// let grid_b = Grid1D::<IntervalClosed<f64>>::try_from_sorted(SortedSet::from_unsorted(vec![0.0, 0.5, 2.0])).unwrap();
    ///
    /// let union = Grid1DUnion::try_new(&grid_a, &grid_b).unwrap();
    /// let (unified_grid, map_a, map_b) = union.into_parts();
    ///
    /// // Now you can use the individual components
    /// println!("Unified grid has {} intervals", unified_grid.num_intervals().as_ref());
    /// ```
    pub fn into_parts(self) -> (Grid1D<Domain1D>, Vec<IntervalId>, Vec<IntervalId>) {
        (
            self.unified_grid,
            self.mapping_to_grid_a,
            self.mapping_to_grid_b,
        )
    }
}

impl<Domain1D: BuildIntervalInPartition> HasDomain1D for Grid1DUnion<Domain1D> {
    type Domain1D = Domain1D;

    /// Returns a reference to the domain of the unified grid.
    ///
    /// The domain is shared among both original grids and the unified grid,
    /// as grid union operations require identical domains.
    ///
    /// # Example
    /// ```rust
    /// use grid1d::{*, intervals::*, scalars::*};
    /// use sorted_vec::partial::SortedSet;
    ///
    /// let grid_a = Grid1D::<IntervalClosed<f64>>::try_from_sorted(SortedSet::from_unsorted(vec![0.0, 1.0, 2.0])).unwrap();
    /// let grid_b = Grid1D::<IntervalClosed<f64>>::try_from_sorted(SortedSet::from_unsorted(vec![0.0, 0.5, 2.0])).unwrap();
    ///
    /// let union = Grid1DUnion::try_new(&grid_a, &grid_b).unwrap();
    /// assert_eq!(union.domain(), grid_a.domain());
    /// assert_eq!(union.domain(), grid_b.domain());
    /// ```
    fn domain(&self) -> &Self::Domain1D {
        self.unified_grid.domain()
    }
}

impl<Domain1D: BuildIntervalInPartition> HasCoords1D for Grid1DUnion<Domain1D> {
    type Point1DType = Domain1D::RealType;

    /// Returns a reference to the coordinates of the unified grid.
    ///
    /// These coordinates represent the union of all unique points from both
    /// original grids, sorted in ascending order.
    ///
    /// # Example
    /// ```rust
    /// use grid1d::{*, intervals::*, scalars::*};
    /// use sorted_vec::partial::SortedSet;
    /// use std::ops::Deref;
    ///
    /// let grid_a = Grid1D::<IntervalClosed<f64>>::try_from_sorted(SortedSet::from_unsorted(vec![0.0, 1.0, 2.0])).unwrap();
    /// let grid_b = Grid1D::<IntervalClosed<f64>>::try_from_sorted(SortedSet::from_unsorted(vec![0.0, 0.5, 2.0])).unwrap();
    ///
    /// let union = Grid1DUnion::try_new(&grid_a, &grid_b).unwrap();
    /// assert_eq!(union.coords().deref(), &[0.0, 0.5, 1.0, 2.0]);
    /// ```
    fn coords(&self) -> &Coords1D<Self::Point1DType> {
        self.unified_grid.coords()
    }
}

impl<Domain1D: BuildIntervalInPartition> IntervalPartition for Grid1DUnion<Domain1D> {
    type UniformlyRefinedGrid1DType = Grid1D<Domain1D>;

    fn refine_uniform(
        self,
        num_extra_points_each_interval: &PositiveNumPoints1D,
    ) -> Grid1DUniformRefinement<Self> {
        let unified_grid = self.unified_grid;
        let mapping_to_grid_a = self.mapping_to_grid_a;
        let mapping_to_grid_b = self.mapping_to_grid_b;

        let (
            refined_grid,
            original_unified_grid,
            refined_to_original_interval_mapping,
            original_to_refined_interval_mapping,
        ) = unified_grid
            .refine_uniform(num_extra_points_each_interval)
            .into();

        let original_grid = Grid1DUnion {
            unified_grid: original_unified_grid,
            mapping_to_grid_a,
            mapping_to_grid_b,
        };

        Grid1DUniformRefinement::new(
            refined_grid,
            original_grid,
            refined_to_original_interval_mapping,
            original_to_refined_interval_mapping,
        )
    }

    fn refine(
        self,
        intervals_to_refine: &std::collections::BTreeMap<IntervalId, PositiveNumPoints1D>,
    ) -> Grid1DNonUniformRefinement<Self> {
        let unified_grid = self.unified_grid;
        let mapping_to_grid_a = self.mapping_to_grid_a;
        let mapping_to_grid_b = self.mapping_to_grid_b;

        let (
            refined_grid,
            original_unified_grid,
            refined_to_original_interval_mapping,
            original_to_refined_interval_mapping,
        ) = unified_grid.refine(intervals_to_refine).into();

        let original_grid = Grid1DUnion {
            unified_grid: original_unified_grid,
            mapping_to_grid_a,
            mapping_to_grid_b,
        };

        Grid1DNonUniformRefinement::new(
            refined_grid,
            original_grid,
            refined_to_original_interval_mapping,
            original_to_refined_interval_mapping,
        )
    }
}
//===============================================================================
// TESTS
//===============================================================================

#[cfg(test)]
mod tests {
    use super::*;
    use crate::{
        intervals::IntervalClosed,
        scalars::NumIntervals,
        traits::{HasCoords1D, HasDomain1D, IntervalPartition},
    };
    use num_valid::RealNative64StrictFinite;
    use sorted_vec::partial::SortedSet;
    use std::ops::Deref;
    use try_create::TryNew;

    type Real = RealNative64StrictFinite;

    /// Helper function to create a Real value from f64
    fn real(x: f64) -> Real {
        Real::try_new(x).unwrap()
    }

    mod construction {
        use super::*;

        #[test]
        fn union_of_identical_grids() {
            let grid = Grid1D::uniform(
                IntervalClosed::new(real(0.0), real(2.0)),
                NumIntervals::try_new(2).unwrap(),
            );
            let union = Grid1DUnion::try_new(&grid, &grid).unwrap();

            // Union of identical grids should be the same grid
            assert_eq!(union.unified_grid().coords(), grid.coords());
            assert_eq!(union.num_refined_intervals(), grid.num_intervals());

            // All mappings should point to corresponding intervals
            for i in 0..*grid.num_intervals().as_ref() {
                assert_eq!(union.mapping_to_grid_a()[i], IntervalId::new(i));
                assert_eq!(union.mapping_to_grid_b()[i], IntervalId::new(i));
            }
        }

        #[test]
        fn union_with_disjoint_points() {
            let grid_a =
                Grid1D::<IntervalClosed<Real>>::try_from_sorted(SortedSet::from_unsorted(vec![
                    real(0.0),
                    real(2.0),
                    real(4.0),
                ]))
                .unwrap();
            let grid_b =
                Grid1D::<IntervalClosed<Real>>::try_from_sorted(SortedSet::from_unsorted(vec![
                    real(0.0),
                    real(1.0),
                    real(3.0),
                    real(4.0),
                ]))
                .unwrap();

            let union = Grid1DUnion::try_new(&grid_a, &grid_b).unwrap();

            // Unified grid should have all unique points
            assert_eq!(
                union.unified_grid().coords().deref(),
                &[real(0.0), real(1.0), real(2.0), real(3.0), real(4.0)]
            );
            assert_eq!(union.num_refined_intervals().as_ref(), &4);
        }

        #[test]
        fn union_preserves_domain() {
            let domain = IntervalClosed::new(real(-5.0), real(10.0));
            let grid_a = Grid1D::uniform(domain.clone(), NumIntervals::try_new(3).unwrap());
            let grid_b = Grid1D::uniform(domain.clone(), NumIntervals::try_new(5).unwrap());

            let union = Grid1DUnion::try_new(&grid_a, &grid_b).unwrap();

            assert_eq!(union.domain(), &domain);
        }

        #[test]
        fn union_with_uniform_and_non_uniform() {
            let domain = IntervalClosed::new(real(0.0), real(1.0));
            let uniform_grid = Grid1D::uniform(domain, NumIntervals::try_new(4).unwrap());
            let non_uniform_grid =
                Grid1D::<IntervalClosed<Real>>::try_from_sorted(SortedSet::from_unsorted(vec![
                    real(0.0),
                    real(0.1),
                    real(0.9),
                    real(1.0),
                ]))
                .unwrap();

            let union = Grid1DUnion::try_new(&uniform_grid, &non_uniform_grid).unwrap();

            // Should contain points from both grids
            assert!(union.unified_grid().coords().len() >= 5); // At least 5 unique points
        }
    }

    mod domain_mismatch {
        use super::*;

        #[test]
        fn different_upper_bounds() {
            let grid_a = Grid1D::uniform(
                IntervalClosed::new(real(0.0), real(1.0)),
                NumIntervals::try_new(2).unwrap(),
            );
            let grid_b = Grid1D::uniform(
                IntervalClosed::new(real(0.0), real(2.0)),
                NumIntervals::try_new(2).unwrap(),
            );

            let result = Grid1DUnion::try_new(&grid_a, &grid_b);
            assert!(matches!(
                result,
                Err(ErrorsGrid1DUnion::DomainsMismatch { .. })
            ));
        }

        #[test]
        fn different_lower_bounds() {
            let grid_a = Grid1D::uniform(
                IntervalClosed::new(real(-1.0), real(1.0)),
                NumIntervals::try_new(2).unwrap(),
            );
            let grid_b = Grid1D::uniform(
                IntervalClosed::new(real(0.0), real(1.0)),
                NumIntervals::try_new(2).unwrap(),
            );

            let result = Grid1DUnion::try_new(&grid_a, &grid_b);
            assert!(matches!(
                result,
                Err(ErrorsGrid1DUnion::DomainsMismatch { .. })
            ));
        }

        #[test]
        fn completely_different_domains() {
            let grid_a = Grid1D::uniform(
                IntervalClosed::new(real(0.0), real(1.0)),
                NumIntervals::try_new(2).unwrap(),
            );
            let grid_b = Grid1D::uniform(
                IntervalClosed::new(real(10.0), real(20.0)),
                NumIntervals::try_new(2).unwrap(),
            );

            let result = Grid1DUnion::try_new(&grid_a, &grid_b);
            assert!(matches!(
                result,
                Err(ErrorsGrid1DUnion::DomainsMismatch { .. })
            ));
        }
    }

    mod mappings {
        use super::*;

        #[test]
        fn find_original_intervals() {
            let grid_a =
                Grid1D::<IntervalClosed<Real>>::try_from_sorted(SortedSet::from_unsorted(vec![
                    real(0.0),
                    real(1.0),
                    real(2.0),
                ]))
                .unwrap();
            let grid_b =
                Grid1D::<IntervalClosed<Real>>::try_from_sorted(SortedSet::from_unsorted(vec![
                    real(0.0),
                    real(0.5),
                    real(2.0),
                ]))
                .unwrap();

            let union = Grid1DUnion::try_new(&grid_a, &grid_b).unwrap();

            // Verify find_original_intervals matches direct mapping access
            for i in 0..*union.num_refined_intervals().as_ref() {
                let refined_id = IntervalId::new(i);
                let (a_id, b_id) = union.find_original_intervals(&refined_id);
                assert_eq!(a_id, union.mapping_to_grid_a()[i]);
                assert_eq!(b_id, union.mapping_to_grid_b()[i]);
            }
        }

        #[test]
        fn mapping_to_grid_a() {
            let domain = IntervalClosed::new(real(0.0), real(2.0));
            let grid_a = Grid1D::uniform(domain.clone(), NumIntervals::try_new(2).unwrap());
            let grid_b = Grid1D::uniform(domain, NumIntervals::try_new(4).unwrap());

            let union = Grid1DUnion::try_new(&grid_a, &grid_b).unwrap();

            // Grid B has more intervals, so unified grid should match grid_b
            // Grid A intervals: [0, 1], [1, 2]
            // Grid B intervals: [0, 0.5], [0.5, 1], [1, 1.5], [1.5, 2]
            // Mapping to A: [0, 0, 1, 1]
            assert_eq!(union.mapping_to_grid_a()[0], IntervalId::new(0));
            assert_eq!(union.mapping_to_grid_a()[1], IntervalId::new(0));
            assert_eq!(union.mapping_to_grid_a()[2], IntervalId::new(1));
            assert_eq!(union.mapping_to_grid_a()[3], IntervalId::new(1));
        }
    }

    mod iteration {
        use super::*;

        #[test]
        fn iter_interval_mappings() {
            let grid_a =
                Grid1D::<IntervalClosed<Real>>::try_from_sorted(SortedSet::from_unsorted(vec![
                    real(0.0),
                    real(1.0),
                    real(2.0),
                ]))
                .unwrap();
            let grid_b =
                Grid1D::<IntervalClosed<Real>>::try_from_sorted(SortedSet::from_unsorted(vec![
                    real(0.0),
                    real(0.5),
                    real(2.0),
                ]))
                .unwrap();

            let union = Grid1DUnion::try_new(&grid_a, &grid_b).unwrap();

            let mappings: Vec<(IntervalId, IntervalId, IntervalId)> =
                union.iter_interval_mappings().collect();

            assert_eq!(mappings.len(), *union.num_refined_intervals().as_ref());

            // Verify consistency with find_original_intervals
            for (refined_id, expected_a, expected_b) in &mappings {
                let (actual_a, actual_b) = union.find_original_intervals(refined_id);
                assert_eq!(actual_a, *expected_a);
                assert_eq!(actual_b, *expected_b);
            }
        }

        #[test]
        fn iter_intervals_with_mappings() {
            let domain = IntervalClosed::new(real(0.0), real(2.0));
            let grid_a = Grid1D::uniform(domain.clone(), NumIntervals::try_new(2).unwrap());
            let grid_b = Grid1D::uniform(domain, NumIntervals::try_new(4).unwrap());

            let union = Grid1DUnion::try_new(&grid_a, &grid_b).unwrap();

            let interval_mappings: Vec<_> = union.iter_intervals_with_mappings().collect();

            assert_eq!(
                interval_mappings.len(),
                *union.num_refined_intervals().as_ref()
            );

            // Each tuple should have (refined_id, interval, a_id, b_id)
            for (refined_id, _interval, a_id, b_id) in &interval_mappings {
                let (check_a, check_b) = union.find_original_intervals(refined_id);
                assert_eq!(&check_a, a_id);
                assert_eq!(&check_b, b_id);
            }
        }
    }

    mod intersection_queries {
        use super::*;

        #[test]
        fn find_intersections_with_mappings() {
            let domain = IntervalClosed::new(real(0.0), real(3.0));
            let grid_a = Grid1D::uniform(domain, NumIntervals::try_new(3).unwrap());
            let grid_b =
                Grid1D::<IntervalClosed<Real>>::try_from_sorted(SortedSet::from_unsorted(vec![
                    real(0.0),
                    real(0.5),
                    real(1.5),
                    real(3.0),
                ]))
                .unwrap();

            let union = Grid1DUnion::try_new(&grid_a, &grid_b).unwrap();
            let subdomain = IntervalClosed::new(real(0.7), real(2.2));

            let intersections = union.find_intersections_with_mappings(&subdomain);
            assert!(!intersections.is_empty());

            // Verify that each intersection has valid mapping information
            for (refined_id, _intersection, a_id, b_id) in &intersections {
                let (check_a, check_b) = union.find_original_intervals(refined_id);
                assert_eq!(check_a, *a_id);
                assert_eq!(check_b, *b_id);
            }
        }

        #[test]
        fn intersection_subdomain_at_boundary() {
            let domain = IntervalClosed::new(real(0.0), real(2.0));
            let grid_a = Grid1D::uniform(domain.clone(), NumIntervals::try_new(4).unwrap());
            let grid_b = Grid1D::uniform(domain.clone(), NumIntervals::try_new(4).unwrap());

            let union = Grid1DUnion::try_new(&grid_a, &grid_b).unwrap();

            // Query subdomain that starts at a grid point
            let subdomain = IntervalClosed::new(real(0.5), real(1.5));
            let intersections = union.find_intersections_with_mappings(&subdomain);

            // Should find intervals in the range [0.5, 1.5]
            assert!(!intersections.is_empty());
        }
    }

    mod decomposition {
        use super::*;

        #[test]
        fn into_parts() {
            let domain = IntervalClosed::new(real(0.0), real(1.0));
            let grid_a = Grid1D::uniform(domain.clone(), NumIntervals::try_new(2).unwrap());
            let grid_b = Grid1D::uniform(domain, NumIntervals::try_new(4).unwrap());

            let union = Grid1DUnion::try_new(&grid_a, &grid_b).unwrap();
            let num_intervals = *union.num_refined_intervals().as_ref();

            let (unified_grid, map_a, map_b) = union.into_parts();

            assert_eq!(unified_grid.num_intervals().as_ref(), &num_intervals);
            assert_eq!(map_a.len(), num_intervals);
            assert_eq!(map_b.len(), num_intervals);
        }

        #[test]
        fn decomposition_preserves_data() {
            let domain = IntervalClosed::new(real(0.0), real(2.0));
            let grid_a = Grid1D::uniform(domain.clone(), NumIntervals::try_new(2).unwrap());
            let grid_b = Grid1D::uniform(domain, NumIntervals::try_new(4).unwrap());

            let union = Grid1DUnion::try_new(&grid_a, &grid_b).unwrap();

            // Store expected values before decomposition
            let expected_coords = union.unified_grid().coords().deref().to_vec();
            let expected_map_a: Vec<_> = union.mapping_to_grid_a().to_vec();
            let expected_map_b: Vec<_> = union.mapping_to_grid_b().to_vec();

            let (unified_grid, map_a, map_b) = union.into_parts();

            assert_eq!(unified_grid.coords().deref(), &expected_coords);
            assert_eq!(map_a, expected_map_a);
            assert_eq!(map_b, expected_map_b);
        }
    }

    mod refinement_of_union {
        use super::*;
        use std::collections::BTreeMap;

        #[test]
        fn uniform_refinement_of_union() {
            let domain = IntervalClosed::new(real(0.0), real(2.0));
            let grid_a = Grid1D::uniform(domain.clone(), NumIntervals::try_new(2).unwrap());
            let grid_b = Grid1D::uniform(domain, NumIntervals::try_new(4).unwrap());

            let union = Grid1DUnion::try_new(&grid_a, &grid_b).unwrap();
            let original_intervals = *union.num_refined_intervals().as_ref();

            // Uniform refinement doubles the number of intervals
            let uniform_refinement =
                union.refine_uniform(&PositiveNumPoints1D::try_new(1).unwrap());

            assert_eq!(
                uniform_refinement.refined_grid().num_intervals().as_ref(),
                &(original_intervals * 2)
            );
        }

        #[test]
        fn selective_refinement_of_union() {
            let domain = IntervalClosed::new(real(0.0), real(2.0));
            let grid_a = Grid1D::uniform(domain.clone(), NumIntervals::try_new(2).unwrap());
            let grid_b = Grid1D::uniform(domain, NumIntervals::try_new(4).unwrap());

            let union = Grid1DUnion::try_new(&grid_a, &grid_b).unwrap();

            let refinement_plan =
                BTreeMap::from([(IntervalId::new(1), PositiveNumPoints1D::try_new(1).unwrap())]);
            let selective_refinement = union.refine(&refinement_plan);

            // Original 4 intervals, refine interval 1 into 2 -> 5 intervals
            assert_eq!(
                selective_refinement.refined_grid().num_intervals().as_ref(),
                &5
            );
        }

        #[test]
        fn chained_union_refinements() {
            let domain = IntervalClosed::new(real(0.0), real(2.0));
            let grid_a = Grid1D::uniform(domain.clone(), NumIntervals::try_new(2).unwrap());
            let grid_b = Grid1D::uniform(domain, NumIntervals::try_new(4).unwrap());

            let union = Grid1DUnion::try_new(&grid_a, &grid_b).unwrap();

            // First uniform refinement
            let first_refinement = union.refine_uniform(&PositiveNumPoints1D::try_new(1).unwrap());
            assert_eq!(first_refinement.refined_grid().num_intervals().as_ref(), &8);

            // Second selective refinement on the result
            let refinement_plan =
                BTreeMap::from([(IntervalId::new(1), PositiveNumPoints1D::try_new(1).unwrap())]);
            let second_refinement = first_refinement
                .into_original_grid()
                .refine(&refinement_plan);
            assert_eq!(
                second_refinement.refined_grid().num_intervals().as_ref(),
                &5
            );
        }
    }

    mod trait_implementations {
        use super::*;

        #[test]
        fn has_domain_1d() {
            let domain = IntervalClosed::new(real(0.0), real(2.0));
            let grid_a = Grid1D::uniform(domain.clone(), NumIntervals::try_new(2).unwrap());
            let grid_b = Grid1D::uniform(domain.clone(), NumIntervals::try_new(4).unwrap());

            let union = Grid1DUnion::try_new(&grid_a, &grid_b).unwrap();

            assert_eq!(union.domain(), &domain);
        }

        #[test]
        fn has_coords_1d() {
            let grid_a =
                Grid1D::<IntervalClosed<Real>>::try_from_sorted(SortedSet::from_unsorted(vec![
                    real(0.0),
                    real(1.0),
                    real(2.0),
                ]))
                .unwrap();
            let grid_b =
                Grid1D::<IntervalClosed<Real>>::try_from_sorted(SortedSet::from_unsorted(vec![
                    real(0.0),
                    real(0.5),
                    real(2.0),
                ]))
                .unwrap();

            let union = Grid1DUnion::try_new(&grid_a, &grid_b).unwrap();

            // Coords should be the union of all unique points
            assert_eq!(
                union.coords().deref(),
                &[real(0.0), real(0.5), real(1.0), real(2.0)]
            );
        }

        #[test]
        fn interval_partition_num_intervals() {
            let domain = IntervalClosed::new(real(0.0), real(2.0));
            let grid_a = Grid1D::uniform(domain.clone(), NumIntervals::try_new(2).unwrap());
            let grid_b = Grid1D::uniform(domain, NumIntervals::try_new(4).unwrap());

            let union = Grid1DUnion::try_new(&grid_a, &grid_b).unwrap();

            // The unified grid should have the resolution of the finer grid
            assert_eq!(union.num_intervals().as_ref(), &4);
        }
    }

    mod edge_cases {
        use super::*;

        #[test]
        fn union_with_single_interval_grids() {
            let domain = IntervalClosed::new(real(0.0), real(1.0));
            let grid_a = Grid1D::uniform(domain.clone(), NumIntervals::try_new(1).unwrap());
            let grid_b = Grid1D::uniform(domain, NumIntervals::try_new(1).unwrap());

            let union = Grid1DUnion::try_new(&grid_a, &grid_b).unwrap();

            assert_eq!(union.num_refined_intervals().as_ref(), &1);
        }

        #[test]
        fn union_with_very_different_resolutions() {
            let domain = IntervalClosed::new(real(0.0), real(1.0));
            let coarse_grid = Grid1D::uniform(domain.clone(), NumIntervals::try_new(2).unwrap());
            let fine_grid = Grid1D::uniform(domain, NumIntervals::try_new(100).unwrap());

            let union = Grid1DUnion::try_new(&coarse_grid, &fine_grid).unwrap();

            // Should have at least as many intervals as the fine grid
            assert!(*union.num_refined_intervals().as_ref() >= 100);
        }

        #[test]
        fn union_of_non_uniform_grids() {
            let grid_a =
                Grid1D::<IntervalClosed<Real>>::try_from_sorted(SortedSet::from_unsorted(vec![
                    real(0.0),
                    real(0.1),
                    real(0.2),
                    real(1.0),
                ]))
                .unwrap();
            let grid_b =
                Grid1D::<IntervalClosed<Real>>::try_from_sorted(SortedSet::from_unsorted(vec![
                    real(0.0),
                    real(0.8),
                    real(0.9),
                    real(1.0),
                ]))
                .unwrap();

            let union = Grid1DUnion::try_new(&grid_a, &grid_b).unwrap();

            // Unified grid should have all unique points from both grids
            let expected_points = vec![
                real(0.0),
                real(0.1),
                real(0.2),
                real(0.8),
                real(0.9),
                real(1.0),
            ];
            assert_eq!(union.unified_grid().coords().deref(), &expected_points);
        }

        #[test]
        fn union_with_overlapping_points() {
            // Both grids share some points
            let grid_a =
                Grid1D::<IntervalClosed<Real>>::try_from_sorted(SortedSet::from_unsorted(vec![
                    real(0.0),
                    real(0.5),
                    real(1.0),
                    real(2.0),
                ]))
                .unwrap();
            let grid_b =
                Grid1D::<IntervalClosed<Real>>::try_from_sorted(SortedSet::from_unsorted(vec![
                    real(0.0),
                    real(1.0),
                    real(1.5),
                    real(2.0),
                ]))
                .unwrap();

            let union = Grid1DUnion::try_new(&grid_a, &grid_b).unwrap();

            // Unified should have unique points: 0, 0.5, 1, 1.5, 2
            assert_eq!(
                union.unified_grid().coords().deref(),
                &[real(0.0), real(0.5), real(1.0), real(1.5), real(2.0)]
            );
            assert_eq!(union.num_refined_intervals().as_ref(), &4);
        }
    }
}