grid1d 0.5.0 - Docs.rs

#![deny(rustdoc::broken_intra_doc_links)]

//! Grid refinement operations for adaptive mesh refinement (AMR) and hierarchical methods.
//!
//! This module provides comprehensive support for subdividing 1D grids into finer
//! partitions while maintaining bidirectional mappings between original and refined
//! intervals. These capabilities are essential for adaptive numerical methods,
//! multigrid solvers, and any application requiring dynamic grid resolution.
//!
//! ## Key Types
//!
//! | Type | Description |
//! |------|-------------|
//! | [`Grid1DRefinement`] | Generic refinement with bidirectional interval mappings |
//! | [`Grid1DUniformRefinement`] | Type alias for refinement of uniform grids |
//! | [`Grid1DNonUniformRefinement`] | Type alias for refinement of non-uniform grids |
//!
//! ## Refinement Modes
//!
//! ### Uniform Refinement
//! All intervals are subdivided by the same factor, creating a uniformly finer grid:
//!
//! ```rust
//! use grid1d::{*, intervals::*, scalars::*};
//! use try_create::TryNew;
//!
//! let grid = Grid1D::uniform(IntervalClosed::new(0.0, 1.0), NumIntervals::try_new(4).unwrap());
//!
//! // Add 1 point per interval → 2 sub-intervals each → 8 total intervals
//! let refinement = grid.refine_uniform(&PositiveNumPoints1D::try_new(1).unwrap());
//! assert_eq!(refinement.refined_grid().num_intervals().as_ref(), &8);
//!
//! // Track which original interval each refined interval came from
//! for (refined_id, original_id) in refinement.iter_refined_with_mapping() {
//!     println!("Refined {} from original {}", *refined_id.as_ref(), *original_id.as_ref());
//! }
//! ```
//!
//! ### Selective Refinement
//! Only specified intervals are subdivided, allowing targeted resolution increases:
//!
//! ```rust
//! use grid1d::{*, intervals::*, scalars::*};
//! use std::collections::BTreeMap;
//! use try_create::TryNew;
//!
//! let grid = Grid1D::uniform(IntervalClosed::new(0.0, 1.0), NumIntervals::try_new(4).unwrap());
//!
//! // Refine only intervals 0 and 2 with different subdivision levels
//! let plan = BTreeMap::from([
//!     (IntervalId::new(0), PositiveNumPoints1D::try_new(3).unwrap()), // 4 sub-intervals
//!     (IntervalId::new(2), PositiveNumPoints1D::try_new(1).unwrap()), // 2 sub-intervals
//! ]);
//! let refinement = grid.refine(&plan);
//!
//! // Intervals 1 and 3 are unchanged
//! assert!(!refinement.was_refined(&IntervalId::new(1)));
//! assert!(!refinement.was_refined(&IntervalId::new(3)));
//! ```
//!
//! ## Bidirectional Mappings
//!
//! The refinement structure maintains two complementary mappings:
//!
//! ### Refined → Original (Forward Mapping)
//! For each refined interval, find which original interval it came from:
//!
//! ```rust
//! # use grid1d::{*, intervals::*, scalars::*};
//! # use try_create::TryNew;
//! # let grid = Grid1D::uniform(IntervalClosed::new(0.0, 1.0), NumIntervals::try_new(2).unwrap());
//! # let refinement = grid.refine_uniform(&PositiveNumPoints1D::try_new(1).unwrap());
//! // O(1) lookup
//! let original = refinement.find_original_interval(&IntervalId::new(0));
//! ```
//!
//! ### Original → Refined (Reverse Mapping)
//! For each original interval, find all refined intervals it produced:
//!
//! ```rust
//! # use grid1d::{*, intervals::*, scalars::*};
//! # use try_create::TryNew;
//! # let grid = Grid1D::uniform(IntervalClosed::new(0.0, 1.0), NumIntervals::try_new(2).unwrap());
//! # let refinement = grid.refine_uniform(&PositiveNumPoints1D::try_new(1).unwrap());
//! // Get all refined intervals from original interval 0
//! let refined_ids = refinement.get_refined_intervals(&IntervalId::new(0));
//! assert_eq!(refined_ids.len(), 2); // Uniform refinement with 1 extra point
//! ```
//!
//! ## Data Transfer
//!
//! Transfer interval-associated data from coarse to fine grids:
//!
//! ```rust
//! # use grid1d::{*, intervals::*, scalars::*};
//! # use try_create::TryNew;
//! # let grid = Grid1D::uniform(IntervalClosed::new(0.0, 1.0), NumIntervals::try_new(2).unwrap());
//! # let refinement = grid.refine_uniform(&PositiveNumPoints1D::try_new(1).unwrap());
//! let coarse_data = vec![1.0, 2.0]; // One value per original interval
//! let fine_data = refinement.transfer_interval_data(&coarse_data);
//! // fine_data = [1.0, 1.0, 2.0, 2.0] - values inherited from parent intervals
//! ```
//!
//! ## Mathematical Properties
//!
//! All refinement operations preserve these invariants:
//!
//! - **Domain preservation**: `refined_grid.domain() == original_grid.domain()`
//! - **Partition completeness**: Refined intervals completely tile the domain
//! - **Ordering preservation**: Interval ordering is maintained
//! - **Bijective sub-interval mapping**: Each refined interval has exactly one parent
//! - **Complete coverage**: Every original interval maps to at least one refined interval
//!
//! ## Performance
//!
//! | Operation | Complexity | Notes |
//! |-----------|------------|-------|
//! | `refine_uniform` | O(n × k) | n = intervals, k = subdivision factor |
//! | `refine` (selective) | O(n + Σkᵢ) | kᵢ = subdivision factor for interval i |
//! | `find_original_interval` | O(1) | Direct array lookup |
//! | `get_refined_intervals` | O(1) | Direct array lookup |
//! | `transfer_interval_data` | O(m) | m = refined interval count |
//!
//! ## Use Cases
//!
//! - **Adaptive Mesh Refinement (AMR)**: Dynamically increase resolution where needed
//! - **Multigrid Methods**: Create grid hierarchies for efficient linear solvers
//! - **Error Estimation**: Compare solutions at different resolutions
//! - **Solution Transfer**: Move data between grids during mesh adaptation
//! - **Hierarchical Algorithms**: Build tree structures for fast spatial queries
//!
//! ## See Also
//!
//! - [`Grid1D::refine`](crate::Grid1D) - Main entry point for selective refinement
//! - [`Grid1D::refine_uniform`](crate::Grid1D) - Entry point for uniform refinement
//! - [`Grid1DUnion`](super::union::Grid1DUnion) - Combining grids from different sources
//! - [`Grid1DTrait`] - Trait implemented by refined grids

use crate::{
    Grid1D, Grid1DIndexSpaces, Grid1DNonUniform, Grid1DUniform, NumIntervals, Topology1D,
    coords::Coords1D,
    grids::traits::{Grid1DIntervalBuilder, Grid1DTrait, HasDomain1D},
    intervals::GetLowerBoundValue,
    scalars::{IntervalId, PositiveNumPoints1D},
};
use derive_more::Into;
use duplicate::duplicate_item;
use getset::Getters;
use num_valid::RealScalar;
use serde::{Deserialize, Serialize};
use sorted_vec::partial::SortedSet;
use std::collections::HashMap;
use try_create::TryNew;

//------------------------------------------------------------------------------------------------------------------------------
/// **INTERNAL** Refines a set of coordinates by subdividing specified intervals and returns a mapping to the original intervals.
///
/// This function takes an existing `Coords1D` object and a "refinement plan" to create a new, more
/// detailed `Coords1D` object. It selectively subdivides specified intervals by inserting a given
/// number of new, uniformly spaced points.
///
/// It also computes a `HashMap` that maps every interval in the *new* refined grid back to the
/// corresponding interval in the *original* grid. This is crucial for adaptive algorithms where data
/// or properties associated with original intervals need to be tracked or migrated to the refined grid.
///
/// ## Core Functionality
///
/// 1.  **Selective Refinement**: Only the intervals specified in `intervals_to_refine` are modified. All other intervals remain unchanged.
/// 2.  **Uniform Subdivision**: Each targeted interval is divided into `n+1` equal-length sub-intervals by adding `n` new points.
/// 3.  **Coordinate Generation**: New points are generated and inserted. The final set of coordinates is correctly sorted and deduplicated.
/// 4.  **Mapping Generation**: A map `new_id -> old_id` is created. If an original interval `j` is subdivided into `n+1` new intervals, the map will contain `n+1` entries pointing back to `j`.
///
/// ## Parameters
///
/// - `old_coords`: The initial set of coordinates to be refined.
/// - `intervals_to_refine`: A slice of tuples `(IntervalId, PositiveNumPoints1D)`.
///   - `IntervalId`: The ID of the interval in the `old_coords` grid to refine.
///   - `PositiveNumPoints1D`: The number of **extra points** to add within the specified interval. Adding `n` points will result in `n+1` new sub-intervals.
///
/// ## Returns
///
/// A tuple `(Coords1D<RealType>, HashMap<IntervalId, IntervalId>)`:
/// 1.  The new, refined `Coords1D` object.
/// 2.  The `HashMap` mapping new interval IDs to their original interval IDs.
///
/// ## Panics
///
/// - If an `IntervalId` in `intervals_to_refine` is out of bounds for the `old_coords`.
fn refine_intervals_in_coords1d<RealType: RealScalar>(
    old_coords: &Coords1D<RealType>,
    intervals_to_refine: &std::collections::BTreeMap<IntervalId, PositiveNumPoints1D>,
) -> (Coords1D<RealType>, Vec<IntervalId>, Vec<Vec<IntervalId>>) {
    let num_old_coords = old_coords.len();
    more_asserts::debug_assert_ge!(
        num_old_coords,
        2,
        "At least two points are required to define intervals."
    );
    let num_old_intervals = num_old_coords - 1;

    if intervals_to_refine.is_empty() {
        let new_to_old_identity_map = (0..num_old_intervals).map(IntervalId::new).collect();

        let old_to_new_interval_map = (0..num_old_intervals)
            .map(|i| vec![IntervalId::new(i)])
            .collect();

        return (
            old_coords.clone(),
            new_to_old_identity_map,
            old_to_new_interval_map,
        );
    }

    // Create a HashMap for efficient lookup of which intervals to refine, and how many subdivisions to create.
    let refinement_map: HashMap<usize, usize> = intervals_to_refine
        .iter()
        .map(|(id, num_extra_points)| {
            let id_val = *id.as_ref();
            if id_val >= num_old_intervals {
                panic!(
                    "IntervalId {} is out of bounds for grid with {} intervals.",
                    id_val, num_old_intervals
                );
            }
            (id_val, *num_extra_points.as_ref() + 1)
        })
        .collect();

    // Estimate capacity for the new coordinate and mapping vectors to reduce reallocations.
    let additional_points: usize = refinement_map.values().map(|&s| s.saturating_sub(1)).sum();
    let expected_new_points = num_old_coords + additional_points;
    let expected_new_intervals = num_old_intervals + additional_points;

    let mut new_coords_vec = Vec::with_capacity(expected_new_points);
    let mut new_to_old_interval_map = Vec::with_capacity(expected_new_intervals);
    let mut old_to_new_interval_map = Vec::with_capacity(num_old_intervals);

    // Iterate over each interval of the old grid.
    let mut refined_interval_idx = 0;
    let old_coords = old_coords.as_ref();
    for old_interval_idx in 0..num_old_intervals {
        // Always add the starting point of the current interval.
        new_coords_vec.push(old_coords[old_interval_idx].clone());

        let subdivisions = refinement_map.get(&old_interval_idx).copied().unwrap_or(1);
        if subdivisions > 1 {
            let p_start = &old_coords[old_interval_idx];
            let p_end = &old_coords[old_interval_idx + 1];

            let interval_length = p_end.clone() - p_start;
            let subdivisions_as_real = RealType::try_from_f64(subdivisions as f64).unwrap();
            let step = interval_length / &subdivisions_as_real;

            // Add the new, intermediate points for this interval.
            for j in 1..subdivisions {
                let j_as_real = RealType::try_from_f64(j as f64).unwrap();
                new_coords_vec.push(p_start.clone() + j_as_real * &step);
            }
        }

        // For each new interval created from this old one, add an entry to the map.
        let original_id = IntervalId::new(old_interval_idx);
        new_to_old_interval_map.resize(new_to_old_interval_map.len() + subdivisions, original_id);

        // Add the ids of the elements in the refinement that correspond to the current coarse element
        old_to_new_interval_map.push(
            (refined_interval_idx..refined_interval_idx + subdivisions)
                .map(IntervalId::new)
                .collect(),
        );

        refined_interval_idx += subdivisions;
    }

    // Add the very last point of the entire grid.
    new_coords_vec.push(
        old_coords
            .last()
            .expect("Old coordinates should not be empty")
            .clone(),
    );

    // Create a SortedSet to ensure uniqueness and order, then build the new Coords1D.
    let new_sorted_set = SortedSet::from_unsorted(new_coords_vec);
    let new_coords = Coords1D::try_from(new_sorted_set)
        .expect("Refinement should always produce valid, non-empty coordinates");

    (new_coords, new_to_old_interval_map, old_to_new_interval_map)
}

/// **INTERNAL** Refines all intervals in a coordinate set with uniform subdivision.
///
/// This function applies the same refinement level to every interval in the coordinate set,
/// creating a uniformly refined grid where each original interval is subdivided into the
/// same number of equal-length sub-intervals. This is the core implementation for uniform
/// grid refinement operations.
///
/// ## Mathematical Definition
///
/// Given coordinates `{p₀, p₁, ..., pₙ}` defining `n` intervals and refinement parameter `k`,
/// this creates new coordinates where each interval `[pᵢ, pᵢ₊₁]` is subdivided into `k+1`
/// equal sub-intervals by inserting `k` uniformly spaced points.
///
/// **Result**: `n × (k+1)` total intervals with uniform subdivision.
///
/// ## Algorithm
///
/// 1. **Uniform subdivision**: Each interval `[pᵢ, pᵢ₊₁]` gets `k` extra points at positions:
///    ```text
///    pᵢ + j × (pᵢ₊₁ - pᵢ)/(k+1)  for j = 1, 2, ..., k
///    ```
/// 2. **Mapping generation**: Creates bidirectional mappings between original and refined intervals
/// 3. **Coordinate merging**: Combines all new points into a sorted, deduplicated coordinate set
///
/// ## Performance Characteristics
///
/// | Aspect | Complexity | Notes |
/// |--------|------------|--------|
/// | **Time** | O(n×k) | Linear in total refined intervals |
/// | **Space** | O(n×k) | Storage for refined coordinates and mappings |
/// | **Uniformity** | Perfect | All intervals get identical treatment |
/// | **Predictability** | Complete | Result size is exactly `n×(k+1)` intervals |
///
/// ## Parameters
///
/// - `old_coords`: The original coordinate set to refine
/// - `num_extra_points_each_interval`: Number of extra points to insert in each interval
///   - Value `k` creates `k+1` sub-intervals per original interval
///   - Must be positive (enforced by [`PositiveNumPoints1D`] type)
///
/// ## Returns
///
/// A tuple containing:
/// 1. **`Coords1D<RealType>`**: The refined coordinate set with all new points
/// 2. **`Vec<IntervalId>`**: Mapping from refined interval IDs to original interval IDs
///    - `refined_to_original[i]` gives the original interval that contains refined interval `i`
/// 3. **`Vec<Vec<IntervalId>>`**: Mapping from original interval IDs to lists of refined interval IDs
///    - `original_to_refined[j]` gives all refined intervals created from original interval `j`
///    - Each inner vector has exactly `k+1` elements (uniform subdivision)
///
/// ## Mathematical Properties
///
/// ### Uniform Subdivision Guarantee
/// ```text
/// ∀i ∈ [0, n-1]: |refined_intervals_from(original_interval[i])| = k+1
/// ```
/// Every original interval produces exactly `k+1` refined intervals.
///
/// ### Equal Sub-interval Length
/// ```text
/// ∀i ∈ [0, n-1], ∀j ∈ [0, k]:
///   length(refined_interval[i×(k+1) + j]) = length(original_interval[i]) / (k+1)
/// ```
/// All sub-intervals from the same original interval have identical length.
///
/// ### Predictable Result Size
/// ```text
/// |refined_coordinates| = |original_coordinates| + n×k
/// |refined_intervals| = n×(k+1)
/// ```
/// The result size is completely predictable and deterministic.
///
/// ### Mapping Consistency
/// ```text
/// ∀refined_id r, original_id o:
///   r ∈ original_to_refined[o] ⟺ refined_to_original[r] = o
/// ```
/// The bidirectional mappings are mathematically consistent.
///
/// ## Usage Context
///
/// This function is used internally by:
/// - **[`Grid1DTrait::refine_uniform`]**: Creating uniformly refined grids (uniform, non-uniform, union)
///
/// ## Example Behavior
///
/// ```rust,ignore
/// // Original coordinates: [0.0, 1.0, 3.0] (2 intervals)
/// // Refinement: 1 extra point per interval (k=1)
///
/// let old_coords = Coords1D::from([0.0, 1.0, 3.0]);
/// let extra_points = PositiveNumPoints1D::try_new(1).unwrap();
///
/// let (refined_coords, refined_to_original, original_to_refined) =
///     refine_all_intervals_in_coords1d_same_subdivisions(&old_coords, &extra_points);
///
/// // Result:
/// // refined_coords = [0.0, 0.5, 1.0, 2.0, 3.0]
/// // refined_to_original = [0, 0, 1, 1] (4 intervals → original sources)
/// // original_to_refined = [[0, 1], [2, 3]] (2 original → refined children)
/// ```
///
/// ## Comparison with Selective Refinement
///
/// | Feature | Uniform Refinement | Selective Refinement |
/// |---------|-------------------|---------------------|
/// | **Intervals affected** | All intervals | Only specified intervals |
/// | **Refinement level** | Same for all | Can vary per interval |
/// | **Result predictability** | Perfect | Variable |
/// | **Memory efficiency** | May be wasteful | More efficient |
/// | **Computational cost** | Higher (uniform) | Lower (targeted) |
/// | **Algorithm complexity** | Simpler | More complex |
///
/// ## Implementation Notes
///
/// - **Delegated implementation**: This function creates a refinement plan for all intervals
///   and delegates to [`refine_intervals_in_coords1d`] for the actual work
/// - **Zero-copy optimization**: Reuses existing algorithms rather than duplicating logic
/// - **Type safety**: Uses [`PositiveNumPoints1D`] to ensure valid refinement parameters
/// - **Memory pre-allocation**: Efficient memory usage through capacity pre-allocation
///
/// ## See Also
///
/// - [`refine_intervals_in_coords1d`]: The general selective refinement function this delegates to
/// - [`Grid1DUniformRefinement`]: The result type for uniform refinement operations
/// - [`PositiveNumPoints1D`]: Type-safe wrapper for the refinement parameter
/// - [`Grid1DTrait::refine_uniform`]: The public interface for uniform refinement
pub(crate) fn refine_all_intervals_in_coords1d_same_subdivisions<RealType: RealScalar>(
    old_coords: &Coords1D<RealType>,
    num_extra_points_each_interval: &PositiveNumPoints1D,
) -> (Coords1D<RealType>, Vec<IntervalId>, Vec<Vec<IntervalId>>) {
    let num_old_intervals = old_coords.len() - 1;
    let intervals_to_refine_with_subdivisions = (0..num_old_intervals)
        .map(|id| (IntervalId::new(id), *num_extra_points_each_interval))
        .collect::<std::collections::BTreeMap<_, _>>();

    refine_intervals_in_coords1d(old_coords, &intervals_to_refine_with_subdivisions)
}
//------------------------------------------------------------------------------------------------------------------------------

//------------------------------------------------------------------------------------------------------------------------------

/// A refinement of a 1D grid preserving bidirectional relationships between original and refined intervals.
///
/// The [`Grid1DRefinement`] struct represents the result of refining a 1D grid, maintaining
/// complete **bidirectional mapping** between the original coarse grid and the new refined grid.
/// This is fundamental for adaptive mesh refinement (AMR), multigrid methods, and any numerical
/// algorithm that requires transferring data between grid levels.
///
/// ## Core Mathematical Properties
///
/// - **Interval Preservation**: Every refined interval maps to exactly one original interval
/// - **Complete Coverage**: The union of all refined intervals exactly reconstructs the original domain
/// - **Bidirectional Mapping**: Both refined→original and original→refined mappings are maintained
/// - **Refinement Consistency**: The refinement process preserves the mathematical structure when possible
/// - **Type Safety**: Generic design ensures compatible grid types and prevents type mismatches
/// - **Immutability**: Once created, neither coordinates nor mappings can be modified
///
/// ## Design Philosophy
///
/// ### Bidirectional Mapping Innovation
/// The core innovation is maintaining **both directions** of the interval mapping efficiently:
/// - **Refined → Original**: `refined_to_original_interval_mapping[i] = original_interval_id`
/// - **Original → Refined**: `original_to_refined_interval_mapping[j] = [refined_interval_ids]`
///
/// This enables:
/// - **Data Transfer**: Moving solution data both ways between coarse and fine grids
/// - **Error Estimation**: Computing refinement indicators and back-projecting errors
/// - **Multigrid Operations**: V-cycles, W-cycles, and restriction/prolongation operators
/// - **Adaptive Algorithms**: Dynamic refinement and coarsening decisions
/// - **Solution Analysis**: Comparing solutions at different resolution levels
///
/// ### Performance-Optimized Data Structures
/// The implementation uses optimized data structures for maximum efficiency:
/// - **`Vec<IntervalId>`** for refined→original mapping: O(1) access by array indexing
/// - **`Vec<Vec<IntervalId>>`** for original→refined mapping: O(1) access to refinement children
/// - **No hash overhead**: Direct array access eliminates hash computation costs
/// - **Cache-friendly**: Sequential memory layout improves CPU cache performance
///
/// ### Type-Preserving Refinement
/// The generic design allows for specialized refinement strategies:
/// - **Uniform → Uniform**: `Grid1DUniform` can remain uniform under global refinement
/// - **Non-uniform → Non-uniform**: General refinement always produces non-uniform grids
/// - **Mixed Refinement**: Different strategies for different grid types
///
/// ## Type Parameters
///
/// - **`OriginalGrid1DType`**: The type of the original grid being refined
/// - **`RefinedGrid1DType`**: The type of the resulting refined grid
///
/// Both types must implement [`Grid1DTrait`] and have compatible point and domain types.
///
/// ## Memory Layout and Performance
///
/// | Component | Data Structure | Access Time | Memory Overhead | Notes |
/// |-----------|----------------|-------------|-----------------|--------|
/// | **Original Grid** | Complete grid object | O(1) | Full grid size | For rollback/analysis |
/// | **Refined Grid** | Complete grid object | O(1) | Full grid size | The new refined partition |
/// | **Refined→Original** | `Vec<IntervalId>` | O(1) | 8 bytes × refined intervals | Array indexing |
/// | **Original→Refined** | `Vec<Vec<IntervalId>>` | O(1) | ~16 bytes × original intervals | Geometric ordering preserved |
///
/// **Total memory**: ~2× grid size + mapping overhead (~24 bytes per original interval)
///
/// ## Key Features and API
///
/// ### Bidirectional Data Transfer Operations
/// ```rust
/// use grid1d::{*, intervals::*, scalars::*};
/// use sorted_vec::partial::SortedSet;
/// use try_create::TryNew;
///
/// let original = Grid1D::uniform(
///     IntervalClosed::new(0.0, 1.0),
///     NumIntervals::try_new(4).unwrap()
/// );
///
/// let refinement = original.refine_uniform(&PositiveNumPoints1D::try_new(1).unwrap());
///
/// // Forward transfer: original → refined
/// let original_solution = vec![1., 2., 3., 4.]; // Data on original intervals
/// let refined_solution = refinement.transfer_interval_data(&original_solution);
/// assert_eq!(&refined_solution, &[1., 1., 2., 2., 3., 3., 4., 4.]);
///
/// // Backward access: find which original interval produced a refined interval
/// let refined_id = IntervalId::new(3);
/// let original_id = refinement.find_original_interval(&refined_id);
/// println!("Refined interval {} came from original interval {}",
///          refined_id.as_ref(), original_id.as_ref());
/// ```
///
/// ### Refinement Analysis and Inspection
/// ```rust
/// # use std::collections::BTreeMap;
/// # use grid1d::{*, intervals::*, scalars::*};
/// # use try_create::TryNew;
/// # let original = Grid1D::uniform(IntervalClosed::new(0.0, 1.0), NumIntervals::try_new(4).unwrap());
/// # let refinement = original.refine(&BTreeMap::from([(IntervalId::new(1), PositiveNumPoints1D::try_new(2).unwrap())]));
///
/// // Check which intervals were refined
/// let num_original_intervals = refinement.original_grid().num_intervals();
/// for original_id in 0..num_original_intervals.into_inner() {
///     let original_id = IntervalId::new(original_id);
///     if refinement.was_refined(&original_id) {
///         let refined_intervals = refinement.get_refined_intervals(&original_id);
///         println!("Original interval {} was split into {} intervals: {:?}",
///                  original_id.as_ref(), refined_intervals.len(), refined_intervals);
///     }
/// }
/// ```
///
/// ### Advanced Mapping Operations
/// ```rust
/// # use std::collections::BTreeMap;
/// # use grid1d::{*, intervals::*, scalars::*};
/// # use try_create::TryNew;
/// # let original = Grid1D::uniform(IntervalClosed::new(0.0, 1.0), NumIntervals::try_new(4).unwrap());
/// # let refinement = original.refine(&BTreeMap::from([(IntervalId::new(1), PositiveNumPoints1D::try_new(1).unwrap())]));
///
/// // Iterate over original intervals with their refined counterparts (preserves geometric order)
/// for (original_id, refined_ids) in refinement.iter_original_with_refined() {
///     println!("Original interval {}: refined into {:?}",
///              original_id.as_ref(), refined_ids);
///     
///     // Process refined intervals in geometric order (left to right)
///     for (local_idx, &refined_id) in refined_ids.iter().enumerate() {
///         println!("  Sub-interval {}: refined ID {}", local_idx, refined_id.as_ref());
///     }
/// }
///
/// // Iterate over refined intervals with their original sources
/// for (refined_id, original_id) in refinement.iter_refined_with_mapping() {
///     println!("Refined interval {} came from original interval {}",
///              refined_id.as_ref(), original_id.as_ref());
/// }
/// ```
///
/// ## Performance Characteristics
///
/// | Operation | Time Complexity | Space Complexity | Notes |
/// |-----------|----------------|------------------|--------|
/// | **Creation** | O(n+m) | O(n+m) | Where n=original size, m=refined size |
/// | **Forward Data Transfer** | O(m) | O(m) | Linear in refined grid size |
/// | **Mapping Lookup** | O(1) | O(1) | Direct array indexing |
/// | **Refinement Analysis** | O(1) per query | O(1) | Direct access to counts and mappings |
/// | **Geometric Iteration** | O(k) | O(1) | Where k = intervals in range |
/// | **Grid Access** | O(1) | O(1) | Direct references |
///
/// ## Mathematical Guarantees
///
/// ### Completeness Property
/// ```text
/// ∀ refined_interval_i ∃! original_interval_j : refined_interval_i ⊆ original_interval_j
/// ```
/// Every refined interval is completely contained within exactly one original interval.
///
/// ### Coverage Property
/// ```text
/// ⋃ refined_intervals = ⋃ original_intervals = domain
/// ```
/// The refined grid covers exactly the same domain as the original grid.
///
/// ### Bidirectional Consistency
/// ```text
/// ∀ original_id j, refined_id i:
///   refined_id ∈ original_to_refined[j] ⟺ refined_to_original[i] = j
/// ```
/// The bidirectional mappings are mathematically consistent.
///
/// ### Partition Completeness
/// ```text
/// ∀ original_id j: ⋃_{i ∈ original_to_refined[j]} refined_interval_i = original_interval_j
/// ```
/// The refined intervals completely partition each original interval.
///
/// ### Geometric Ordering Preservation
/// ```text
/// ∀ original_id j: original_to_refined[j] is sorted in ascending geometric order
/// ```
/// Refined intervals maintain left-to-right spatial ordering for numerical algorithms.
///
/// ## Advanced Usage Patterns
///
///
/// ### Error-Based Adaptive Refinement with Geometric Analysis
/// ```rust
/// # use grid1d::{*, intervals::*, scalars::*};
/// # use try_create::TryNew;
/// # fn compute_local_error(_: &Grid1DNonUniform<IntervalClosed<f64>>, _: &[f64]) -> Vec<f64> { vec![0.1, 0.01] }
/// # fn estimate_solution_error(_: &Grid1DNonUniform<IntervalClosed<f64>>, _: &[f64]) -> Vec<f64> { vec![0.1, 0.01] }
///
/// fn adaptive_refinement_cycle(
///     mut grid: Grid1DNonUniform<IntervalClosed<f64>>,
///     mut solution: Vec<f64>,
///     error_threshold: f64,
///     max_iterations: usize,
/// ) -> (Grid1DNonUniform<IntervalClosed<f64>>, Vec<f64>) {
///     for iteration in 0..max_iterations {
///         // Compute error indicators
///         let errors = compute_local_error(&grid, &solution);
///         
///         // Mark intervals for refinement
///         let refinement_plan: std::collections::BTreeMap<IntervalId, PositiveNumPoints1D> = errors
///             .iter()
///             .enumerate()
///             .filter(|&(_, error)| *error > error_threshold)
///             .map(|(i, _)| (IntervalId::new(i), PositiveNumPoints1D::try_new(1).unwrap()))
///             .collect();
///         
///         if refinement_plan.is_empty() {
///             break; // Converged
///         }
///         
///         // Perform refinement
///         let refinement = grid.refine(&refinement_plan);
///         
///         // Analyze refinement pattern
///         println!("Iteration {}: refined {} intervals", iteration, refinement_plan.len());
///         for (original_id, refined_ids) in refinement.iter_original_with_refined() {
///             if refined_ids.len() > 1 {
///                 println!("  Original interval {} → {} refined intervals",
///                          original_id.as_ref(), refined_ids.len());
///             }
///         }
///         
///         // Transfer solution to refined grid
///         solution = refinement.transfer_interval_data(&solution);
///         
///         // Update grid
///         grid = refinement.into_refined_grid();
///     }
///     
///     (grid, solution)
/// }
/// ```
///
/// ### Finite Element Assembly with Refinement Mapping
/// ```rust
/// # use grid1d::{*, scalars::*};
///
/// fn assemble_with_refinement<G1, G2>(
///     refinement: &Grid1DRefinement<G1, G2>,
///     coarse_elements: &[f64],
///     fine_elements: &[f64],
/// ) -> Vec<f64>
/// where
///     G1: Grid1DTrait,
///     G2: Grid1DTrait<CoordType = G1::CoordType, Point1DType = G1::Point1DType, Domain1D = G1::Domain1D>,
/// {
///     let mut assembly = vec![0.0; *refinement.refined_grid().num_intervals().as_ref()];
///     
///     // Iterate through refined intervals in geometric order
///     for (refined_id, original_id) in refinement.iter_refined_with_mapping() {
///         let refined_idx = *refined_id.as_ref();
///         let original_idx = *original_id.as_ref();
///         
///         // Combine contributions from both grids
///         assembly[refined_idx] = fine_elements[refined_idx] + coarse_elements[original_idx];
///         
///         // Apply refinement-based scaling
///         let num_refined = refinement.num_refined_intervals(&original_id);
///         assembly[refined_idx] /= num_refined as f64;
///     }
///     
///     assembly
/// }
/// ```
///
/// ### Data Quality Assessment
/// ```rust
/// # use grid1d::{*, scalars::*};
///
/// fn assess_refinement_quality<G1, G2>(
///     refinement: &Grid1DRefinement<G1, G2>,
/// ) -> (f64, usize, usize)
/// where
///     G1: Grid1DTrait,
///     G2: Grid1DTrait<CoordType = G1::CoordType, Point1DType = G1::Point1DType, Domain1D = G1::Domain1D>,
/// {
///     let original_count = *refinement.original_grid().num_intervals().as_ref();
///     let refined_count = *refinement.refined_grid().num_intervals().as_ref();
///     let refinement_ratio = refined_count as f64 / original_count as f64;
///     
///     // Count how many original intervals were actually refined
///     let refined_intervals = (0..original_count)
///         .map(|i| IntervalId::new(i))
///         .filter(|id| refinement.was_refined(id))
///         .count();
///     
///     // Find maximum refinement factor
///     let max_refinement_factor = (0..original_count)
///         .map(|i| refinement.num_refined_intervals(&IntervalId::new(i)))
///         .max()
///         .unwrap_or(1);
///     
///     println!("Refinement Quality Assessment:");
///     println!("  Overall ratio: {:.2}", refinement_ratio);
///     println!("  Intervals refined: {}/{}", refined_intervals, original_count);
///     println!("  Max refinement factor: {}", max_refinement_factor);
///     
///     (refinement_ratio, refined_intervals, max_refinement_factor)
/// }
/// ```
///
/// ## Error Handling and Validation
///
/// ### Safe Construction Patterns
/// ```rust
/// # use grid1d::{*, scalars::*};
/// # #[derive(Debug)] enum RefinementError { InvalidData }
///
/// fn safe_data_transfer<T: Clone, G1, G2>(
///     refinement: &Grid1DRefinement<G1, G2>,
///     original_data: &[T],
/// ) -> Result<Vec<T>, RefinementError>
/// where
///     G1: Grid1DTrait,
///     G2: Grid1DTrait<CoordType = G1::CoordType, Point1DType = G1::Point1DType, Domain1D = G1::Domain1D>,
/// {
///     if original_data.len() != *refinement.original_grid().num_intervals().as_ref() {
///         return Err(RefinementError::InvalidData);
///     }
///     
///     Ok(refinement.transfer_interval_data(original_data))
/// }
///
/// fn safe_refinement_analysis<G1, G2>(
///     refinement: &Grid1DRefinement<G1, G2>,
///     interval_id: &IntervalId,
/// ) -> Option<Vec<IntervalId>>
/// where
///     G1: Grid1DTrait,
///     G2: Grid1DTrait<CoordType = G1::CoordType, Point1DType = G1::Point1DType, Domain1D = G1::Domain1D>,
/// {
///     if *interval_id.as_ref() >= *refinement.original_grid().num_intervals().as_ref() {
///         return None;
///     }
///     
///     Some(refinement.get_refined_intervals(interval_id).clone())
/// }
/// ```
///
/// ## Best Practices
///
/// ### When to Use Grid Refinement
/// - **Adaptive mesh refinement**: Dynamic grid adaptation based on solution features
/// - **Multigrid methods**: Creating grid hierarchies for efficient solvers
/// - **Error estimation**: Comparing solutions on different resolution levels
/// - **Solution transfer**: Moving data between computational grids
/// - **Multi-physics coupling**: Transferring data between different physics grids
/// - **Hierarchical algorithms**: Building tree-like data structures for fast access
///
/// ### Performance Optimization
/// - **Batch operations**: Process multiple intervals simultaneously when possible
/// - **Memory reuse**: Cache refinement objects for repeated use in time-stepping
/// - **Type specialization**: Use uniform refinements when structure is preserved
/// - **Sparse refinement**: Only refine intervals that actually need higher resolution
/// - **Geometric ordering**: Leverage the preserved ordering for cache-friendly algorithms
///
/// ### Memory Management
/// - **Refinement lifetime**: Consider the lifetime of refinement objects in long-running simulations
/// - **Data movement**: Use move semantics (`into_refined_grid()`) to avoid unnecessary copies
/// - **Mapping efficiency**: The `Vec`-based mapping is optimal for consecutive interval IDs
/// - **Batch analysis**: Process refinement statistics in batches to minimize memory allocations
///
/// ## Integration with other numerical Ecosystem
///
/// The [`Grid1DRefinement`] integrates seamlessly with:
/// - **[`Grid1DTrait`]**: All refined grids implement the full partition interface
/// - **[`HasDomain1D`]**, **[`HasCoords1D`](crate::HasCoords1D)**: Standard domain and coordinate access
/// - **[`Grid1DUnion`](super::union::Grid1DUnion)**: Can be used as input to grid union operations
/// - **Numerical algorithms**: Compatible with all grid-based computations
/// - **AMR frameworks**: Provides the foundation for adaptive mesh refinement systems
///
/// ## Serialization Support
///
/// [`Grid1DRefinement`] implements [`serde::Serialize`] and [`serde::Deserialize`] for complete
/// state persistence, enabling:
/// - **Checkpoint/restart**: Save computation state for long-running simulations
/// - **Result archival**: Store refinement hierarchies for post-processing
/// - **Distributed computing**: Transfer refinement data between compute nodes
/// - **Workflow integration**: Exchange refinement structures with external tools
///
/// ### Example
/// ```rust
/// use grid1d::{*, intervals::*, scalars::*};
/// use try_create::TryNew;
///
/// let grid = Grid1D::uniform(
///     IntervalClosed::new(0.0, 1.0),
///     NumIntervals::try_new(4).unwrap()
/// );
/// let refinement = grid.refine_uniform(&PositiveNumPoints1D::try_new(1).unwrap());
///
/// // Serialize to JSON
/// let json = serde_json::to_string(&refinement).unwrap();
///
/// // Deserialize and verify
/// type GridType = Grid1D<IntervalClosed<f64>>;
/// let deserialized: Grid1DRefinement<GridType, GridType>
///     = serde_json::from_str(&json).unwrap();
/// assert_eq!(
///     refinement.refined_grid().coords(),
///     deserialized.refined_grid().coords()
/// );
/// ```
///
/// The serialization preserves:
/// - Complete grid structures (original and refined)
/// - All bidirectional interval mappings
/// - Domain and coordinate information
/// - Grid type information (uniform vs. non-uniform)
///
/// ## See Also
///
/// - [`Grid1DTrait`]: The core trait for grid partition operations
/// - [`Grid1D`]: The unified grid interface supporting refinement
/// - [`Grid1DUnion`](super::union::Grid1DUnion): For combining multiple grids
/// - Multigrid methods literature for advanced usage patterns
/// - Adaptive mesh refinement (AMR) literature for refinement strategies
#[derive(Debug, Clone, Into, Getters, Serialize, Deserialize)]
#[serde(
    bound(deserialize = "OriginalGrid1DType: for<'d> serde::Deserialize<'d>, \
                  RefinedGrid1DType: for<'d> serde::Deserialize<'d>")
)]
pub struct Grid1DRefinement<OriginalGrid1DType, RefinedGrid1DType>
where
    OriginalGrid1DType: Grid1DTrait,
    RefinedGrid1DType: Grid1DTrait<
            CoordType = OriginalGrid1DType::CoordType,
            Domain1D = OriginalGrid1DType::Domain1D,
        >,
{
    /// The refined grid containing the subdivided intervals.
    ///
    /// This grid has more intervals than the original grid and is used
    /// for higher-resolution computations.
    #[getset(get = "pub")]
    refined_grid: RefinedGrid1DType,

    /// The original (coarser) grid before refinement.
    ///
    /// This grid is preserved to maintain the mapping between coarse
    /// and fine resolutions.
    #[getset(get = "pub")]
    original_grid: OriginalGrid1DType,

    /// Mapping from refined interval index to original interval ID.
    /// `refined_to_original_interval_mapping[i]` gives the original interval ID
    /// for refined interval i.
    refined_to_original_interval_mapping: Vec<IntervalId>,

    /// Mapping from original interval ID to list of refined interval indices.
    /// `original_to_refined_interval_mapping[j]` gives a list of refined interval IDs
    /// that originated from original interval j.
    original_to_refined_interval_mapping: Vec<Vec<IntervalId>>,
}

impl<OriginalGrid1DType, RefinedGrid1DType> Grid1DRefinement<OriginalGrid1DType, RefinedGrid1DType>
where
    OriginalGrid1DType: Grid1DTrait,
    RefinedGrid1DType: Grid1DTrait<
            CoordType = OriginalGrid1DType::CoordType,
            Point1DType = OriginalGrid1DType::Point1DType,
            Domain1D = OriginalGrid1DType::Domain1D,
        >,
{
    pub(crate) fn new(
        refined_grid: RefinedGrid1DType,
        original_grid: OriginalGrid1DType,
        refined_to_original_interval_mapping: Vec<IntervalId>,
        original_to_refined_interval_mapping: Vec<Vec<IntervalId>>,
    ) -> Self {
        Self {
            refined_grid,
            original_grid,
            refined_to_original_interval_mapping,
            original_to_refined_interval_mapping,
        }
    }

    /// Find which original interval a refined interval came from.
    #[inline(always)]
    #[must_use]
    pub fn find_original_interval(&self, refined_id: &IntervalId) -> IntervalId {
        self.refined_to_original_interval_mapping[*refined_id.as_ref()]
    }

    /// Iterator over refined intervals with their original sources.
    pub fn iter_refined_with_mapping(&self) -> impl Iterator<Item = (IntervalId, IntervalId)> + '_ {
        self.refined_to_original_interval_mapping
            .iter()
            .enumerate()
            .map(|(refined_idx, &original_id)| (IntervalId::new(refined_idx), original_id))
    }

    /// Iterator over all original intervals with their refined counterparts.
    pub fn iter_original_with_refined(
        &self,
    ) -> impl Iterator<Item = (IntervalId, &Vec<IntervalId>)> + '_ {
        self.original_to_refined_interval_mapping
            .iter()
            .enumerate()
            .map(|(original_idx, refined_ids)| (IntervalId::new(original_idx), refined_ids))
    }

    /// Transfer data from original grid intervals to refined grid intervals.
    #[must_use]
    pub fn transfer_interval_data<T: Clone>(&self, original_data: &[T]) -> Vec<T> {
        assert_eq!(
            original_data.len(),
            *self.original_grid.num_intervals().as_ref(),
            "Data length must match number of original intervals"
        );

        self.refined_to_original_interval_mapping
            .iter()
            .map(|&original_id| original_data[*original_id.as_ref()].clone())
            .collect()
    }

    /// Extract the refined grid, consuming the refinement.
    #[inline(always)]
    #[must_use]
    pub fn into_refined_grid(self) -> RefinedGrid1DType {
        self.refined_grid
    }

    /// Extract the original grid, consuming the refinement (rollback).
    #[inline(always)]
    #[must_use]
    pub fn into_original_grid(self) -> OriginalGrid1DType {
        self.original_grid
    }

    /// Get all refined intervals that originated from the given original interval.
    #[inline(always)]
    pub fn get_refined_intervals(&self, original_id: &IntervalId) -> &Vec<IntervalId> {
        &self.original_to_refined_interval_mapping[*original_id.as_ref()]
    }

    /// Iterator over all refined intervals for a given original interval.
    #[inline]
    pub fn iter_refined_intervals(
        &self,
        original_id: &IntervalId,
    ) -> impl Iterator<Item = &IntervalId> {
        self.original_to_refined_interval_mapping[*original_id.as_ref()].iter()
    }

    /// Get the number of refined intervals for a given original interval.
    #[inline(always)]
    pub fn num_refined_intervals(&self, original_id: &IntervalId) -> usize {
        self.original_to_refined_interval_mapping[*original_id.as_ref()].len()
    }

    /// Check if an original interval was refined.
    #[inline(always)]
    pub fn was_refined(&self, original_id: &IntervalId) -> bool {
        self.original_to_refined_interval_mapping[*original_id.as_ref()].len() > 1
    }
}
//-------------------------------------------------------------------------------

//-------------------------------------------------------------------------------
///  A specialized refinement of a 1D grid where all intervals are uniformly subdivided.
///
/// [`Grid1DUniformRefinement`] is a type alias for [`Grid1DRefinement`] that represents the result
/// of **uniform refinement**, where every interval in the original grid is subdivided into the same
/// number of equal-length sub-intervals. This is a fundamental operation in adaptive mesh refinement
/// (AMR), multigrid methods, and numerical algorithms that require systematic grid resolution increase.
///
/// ## Mathematical Definition
///
/// Given an original grid with `n` intervals and a refinement factor `k+1` (where `k` is the number
/// of extra points per interval), uniform refinement creates a new grid where:
/// ```text
/// ∀i ∈ [0, n-1]: original_interval[i] → k+1 refined_intervals
/// ```
///
/// **Total refined intervals**: `n × (k+1)`
///
/// **Uniform spacing preservation**: If the original grid is uniform, the refined grid remains uniform
/// with spacing `δ_refined = δ_original / (k+1)`.
///
/// ## Type Definition
///
/// ```rust,ignore
/// pub type Grid1DUniformRefinement<OriginalGrid1DType> = Grid1DRefinement<
///     OriginalGrid1DType,
///     <OriginalGrid1DType as Grid1DTrait>::UniformlyRefinedGrid1DType,
/// >;
/// ```
///
/// This ensures that:
/// - **Input**: Any grid type implementing [`Grid1DTrait`]
/// - **Output**: The same grid type's uniformly refined variant
/// - **Type safety**: Compile-time guarantee of structural consistency
///
/// ## Core Mathematical Properties
///
/// ### Uniform Subdivision Guarantee
/// ```text
/// ∀i ∈ [0, n-1]: num_refined_intervals(original_interval[i]) = k+1
/// ```
/// Every original interval produces exactly `k+1` refined intervals.
///
/// ### Equal Length Sub-intervals
/// ```text
/// ∀i,j ∈ refined_intervals_from_same_original:
///   length(refined_interval[i]) = length(refined_interval[j])
/// ```
/// All refined intervals from the same original interval have identical length.
///
/// ### Structure Preservation
/// - **Uniform → Uniform**: [`Grid1DUniform`] remains [`Grid1DUniform`] with reduced spacing
/// - **Non-uniform → Non-uniform**: [`Grid1DNonUniform`] becomes [`Grid1DNonUniform`] with more points
/// - **Domain consistency**: The refined grid covers exactly the same domain as the original
///
/// ### Refinement Ratio Predictability
/// ```text
/// refinement_ratio = (k+1) × n / n = k+1
/// ```
/// The total number of intervals increases by exactly the refinement factor.
///
/// ## Design Philosophy
///
/// ### Systematic Grid Enhancement
/// Uniform refinement provides **predictable, systematic improvement** in grid resolution:
/// - **Error reduction**: Theoretical error reduction proportional to `(δ_original/(k+1))^p` where `p` is the method order
/// - **Stability preservation**: Maintains numerical stability properties of the original grid
/// - **Memory predictability**: Memory usage increases by exactly factor `k+1`
/// - **Algorithm compatibility**: Preserves structure needed by numerical methods
///
/// ### Performance Optimization
/// - **Bulk operations**: All intervals refined simultaneously for optimal cache usage
/// - **Vectorization friendly**: Regular patterns enable SIMD optimizations
/// - **Minimal branching**: Uniform treatment reduces conditional logic overhead
/// - **Memory locality**: Sequential refinement improves data access patterns
///
/// ### Integration with Numerical Methods
/// - **Multigrid hierarchies**: Creates perfect level-to-level relationships
/// - **Richardson extrapolation**: Enables systematic error analysis
/// - **Convergence studies**: Provides controlled resolution progression
/// - **Adaptive algorithms**: Forms the foundation for non-uniform refinement
///
/// ## Construction and Usage Examples
///
/// ### Basic Uniform Refinement
/// ```rust
/// use grid1d::{*, intervals::*, scalars::*};
/// use sorted_vec::partial::SortedSet;
/// use std::ops::Deref;
/// use try_create::TryNew;
///
/// // Create original uniform grid
/// let original = Grid1D::uniform(
///     IntervalClosed::new(0.0, 1.0),
///     NumIntervals::try_new(4).unwrap()
/// );
/// assert_eq!(original.coords().deref(), &[0.0, 0.25, 0.5, 0.75, 1.0]);
///
/// // Uniform refinement: add 1 extra point per interval
/// let refinement = original.refine_uniform(&PositiveNumPoints1D::try_new(1).unwrap());
///
/// // Verify refined grid structure
/// assert_eq!(refinement.refined_grid().coords().deref(),
///            &[0.0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1.0]);
/// assert_eq!(refinement.refined_grid().num_intervals().as_ref(), &8);
///
/// // Check mapping consistency
/// for (refined_id, original_id) in refinement.iter_refined_with_mapping() {
///     println!("Refined interval {} came from original interval {}",
///              refined_id.as_ref(), original_id.as_ref());
/// }
/// ```
///
/// ### Multi-Level Refinement Hierarchy
/// ```rust
/// use grid1d::{*, intervals::*, scalars::*};
/// use try_create::TryNew;
///
/// // Start with coarse grid
/// let level_0 = Grid1D::uniform(
///     IntervalClosed::new(0.0, 1.0),
///     NumIntervals::try_new(2).unwrap()
/// );
///
/// // Create refinement hierarchy and verigy geometric progression
/// let level_1 = level_0.refine_uniform(&PositiveNumPoints1D::try_new(1).unwrap());
/// assert_eq!(level_1.refined_grid().num_intervals().as_ref(), &4);  // 2×2
///
/// let level_2 = level_1.into_refined_grid()
///                    .refine_uniform(&PositiveNumPoints1D::try_new(1).unwrap());
/// assert_eq!(level_2.refined_grid().num_intervals().as_ref(), &8);  // 4×2  
///
/// let level_3 = level_2.into_refined_grid()
///                    .refine_uniform(&PositiveNumPoints1D::try_new(1).unwrap());
/// assert_eq!(level_3.refined_grid().num_intervals().as_ref(), &16); // 8×2
///
/// println!("Refinement hierarchy: 2 → 4 → 8 → 16 intervals");
/// ```
///
/// ### Uniform Refinement of Non-Uniform Grids
/// ```rust
/// use grid1d::{*, intervals::*, scalars::*};
/// use sorted_vec::partial::SortedSet;
/// use try_create::TryNew;
///
/// // Start with non-uniform grid (adaptive)
/// let coords = SortedSet::from_unsorted(vec![0.0, 0.1, 0.5, 1.0]);
/// let non_uniform = Grid1D::<IntervalClosed<f64>>::try_from_sorted(
///     coords
/// ).unwrap();
/// let original_count = *non_uniform.num_intervals().as_ref();
/// assert_eq!(original_count, 3);
///
/// // Apply uniform refinement (2 extra points per interval)
/// let refinement = non_uniform.refine_uniform(&PositiveNumPoints1D::try_new(2).unwrap());
///
/// // Each original interval now has 3 sub-intervals
/// for (original_id, refined_ids) in refinement.iter_original_with_refined() {
///     assert_eq!(refined_ids.len(), 3);
///     println!("Original interval {}: refined into {:?}",
///              original_id.as_ref(), refined_ids);
/// }
///
/// // Verify total interval count
/// let refined_count = *refinement.refined_grid().num_intervals().as_ref();
/// assert_eq!(refined_count, original_count * 3);
/// ```
///
/// ### Data Transfer in Uniform Refinement
/// ```rust
/// use grid1d::{*, intervals::*, scalars::*};
/// use try_create::TryNew;
///
/// let original = Grid1D::uniform(
///     IntervalClosed::new(0.0, 2.0),
///     NumIntervals::try_new(4).unwrap()
/// );
///
/// // Original solution data (one value per interval)
/// let original_solution = vec![1.0, 4.0, 9.0, 16.0]; // Example: x² at interval centers
///
/// // Uniform refinement with 1 extra point per interval
/// let refinement = original.refine_uniform(&PositiveNumPoints1D::try_new(1).unwrap());
///
/// // Transfer data: each original interval's data is copied to all its refined intervals
/// let refined_solution = refinement.transfer_interval_data(&original_solution);
/// assert_eq!(refined_solution, vec![1.0, 1.0, 4.0, 4.0, 9.0, 9.0, 16.0, 16.0]);
///
/// // Verify data consistency: each pair of refined intervals has the same value
/// for (original_id, refined_ids) in refinement.iter_original_with_refined() {
///     let original_value = original_solution[*original_id.as_ref()];
///     for &refined_id in refined_ids {
///         assert_eq!(refined_solution[*refined_id.as_ref()], original_value);
///     }
/// }
/// ```
///
/// ## Advanced Usage Patterns
///
/// ### Multigrid V-Cycle with Uniform Grids
/// ```rust
/// use grid1d::{*, intervals::*, scalars::*};
/// use try_create::TryNew;
///
/// fn multigrid_v_cycle_uniform(
///     coarse_grid: Grid1D<IntervalClosed<f64>>,
///     solution: &mut [f64],
///     refinement_factor: usize,
/// ) -> f64 {
///     // Create fine grid through uniform refinement
///     let refinement = coarse_grid.refine_uniform(
///         &PositiveNumPoints1D::try_new(refinement_factor - 1).unwrap()
///     );
///     
///     // Transfer solution to fine grid
///     let fine_solution = refinement.transfer_interval_data(solution);
///     
///     // Simulate computation on fine grid (placeholder)
///     let residual_norm = fine_solution.iter().map(|x| x.powi(2)).sum::<f64>().sqrt();
///     
///     // Transfer refined solution back to coarse grid (averaging)
///     for (original_id, refined_ids) in refinement.iter_original_with_refined() {
///         let avg = refined_ids.iter()
///             .map(|&id| fine_solution[*id.as_ref()])
///             .sum::<f64>() / refined_ids.len() as f64;
///         solution[*original_id.as_ref()] = avg;
///     }
///     
///     residual_norm
/// }
///
/// // Usage
/// let grid = Grid1D::uniform(IntervalClosed::new(0.0, 1.0), NumIntervals::try_new(8).unwrap());
/// let mut solution = vec![1.0; 8];
/// let residual = multigrid_v_cycle_uniform(grid, &mut solution, 2);
/// println!("Residual after V-cycle: {:.6}", residual);
/// ```
///
/// ### Richardson Extrapolation with Systematic Refinement
/// ```rust
/// use grid1d::{*, intervals::*, scalars::*};
/// use try_create::TryNew;
///
/// fn richardson_extrapolation(
///     base_grid: Grid1D<IntervalClosed<f64>>,
///     compute_solution: impl Fn(&Grid1D<IntervalClosed<f64>>) -> f64,
/// ) -> (f64, f64, f64) {
///     // Coarse grid solution
///     let solution_h = compute_solution(&base_grid);
///     
///     // Refined grid (h/2 spacing)
///     let refinement_h2 = base_grid.refine_uniform(&PositiveNumPoints1D::try_new(1).unwrap());
///     let solution_h2 = compute_solution(refinement_h2.refined_grid());
///     
///     // Extra refined grid (h/4 spacing)
///     let refinement_h4 = refinement_h2.into_refined_grid()
///         .refine_uniform(&PositiveNumPoints1D::try_new(1).unwrap());
///     let solution_h4 = compute_solution(refinement_h4.refined_grid());
///     
///     // Richardson extrapolation (assuming 2nd order method)
///     let extrapolated = solution_h4 + (solution_h4 - solution_h2) / 3.0;
///     
///     (solution_h, solution_h2, extrapolated)
/// }
///
/// // Example usage
/// let base_grid = Grid1D::uniform(IntervalClosed::new(0.0, 1.0), NumIntervals::try_new(4).unwrap());
/// let compute_fn = |grid: &Grid1D<_>| {
///     // Placeholder: integrate x² over domain using midpoint rule
///     let dx = 1.0 / *grid.num_intervals().as_ref() as f64;
///     (0..*grid.num_intervals().as_ref())
///         .map(|i| {
///             let x_mid = (i as f64 + 0.5) * dx;
///             x_mid * x_mid * dx
///         })
///         .sum()
/// };
///
/// let (coarse, fine, extrapolated) = richardson_extrapolation(base_grid, compute_fn);
/// println!("Coarse: {:.6}, Fine: {:.6}, Extrapolated: {:.6}", coarse, fine, extrapolated);
/// println!("Analytical: {:.6}", 1.0/3.0); // ∫₀¹ x² dx = 1/3
/// ```
///
/// ### Convergence Analysis with Uniform Refinement
/// ```rust
/// use grid1d::{*, intervals::*, scalars::*};
/// use try_create::TryNew;
///
/// fn convergence_study(
///     domain: IntervalClosed<f64>,
///     initial_intervals: usize,
///     max_levels: usize,
/// ) -> Vec<(usize, f64, f64)> {
///     let mut results = Vec::new();
///     let mut current_grid = Grid1D::uniform(
///         domain,
///         NumIntervals::try_new(initial_intervals).unwrap()
///     );
///     
///     for level in 0..max_levels {
///         let num_intervals = *current_grid.num_intervals().as_ref();
///         let dx = 1.0 / num_intervals as f64;
///         
///         // Compute numerical solution (example: midpoint rule for ∫₀¹ x² dx)
///         let numerical = (0..num_intervals)
///             .map(|i| {
///                 let x_mid = (i as f64 + 0.5) * dx;
///                 x_mid * x_mid * dx
///             })
///             .sum::<f64>();
///         
///         let analytical = 1.0 / 3.0;
///         let error = (numerical - analytical).abs();
///         
///         results.push((num_intervals, numerical, error));
///         
///         // Uniform refinement for next level
///         if level < max_levels - 1 {
///             let refinement = current_grid.refine_uniform(
///                 &PositiveNumPoints1D::try_new(1).unwrap()
///             );
///             current_grid = refinement.into_refined_grid();
///         }
///     }
///     
///     results
/// }
///
/// // Usage
/// let study = convergence_study(IntervalClosed::new(0.0, 1.0), 4, 6);
/// println!("Convergence Study:");
/// println!("Intervals | Numerical | Error    | Ratio");
/// println!("----------|-----------|----------|------");
/// for (i, &(n, num, err)) in study.iter().enumerate() {
///     let ratio = if i > 0 { study[i-1].2 / err } else { 0.0 };
///     println!("{:9} | {:9.6} | {:8.2e} | {:4.1}", n, num, err, ratio);
/// }
/// ```
///
/// ### Adaptive-to-Uniform Refinement Strategy
/// ```rust
/// use grid1d::{*, intervals::*, scalars::*};
/// use sorted_vec::partial::SortedSet;
/// use try_create::TryNew;
///
/// fn adaptive_then_uniform_strategy(
///     base_coords: Vec<f64>,
/// ) -> Grid1DNonUniformRefinement<Grid1D<IntervalClosed<f64>>> {
///     // Start with adaptive (non-uniform) grid
///     let adaptive_grid = Grid1D::<IntervalClosed<f64>>::try_from_sorted(
///         SortedSet::from_unsorted(base_coords)
///     ).unwrap();
///     
///     println!("Adaptive grid: {} intervals", adaptive_grid.num_intervals().as_ref());
///     
///     // Apply uniform refinement to the adaptive grid
///     let uniform_refinement = adaptive_grid.refine_uniform(
///         &PositiveNumPoints1D::try_new(2).unwrap() // 3 sub-intervals per original
///     );
///     
///     println!("After uniform refinement: {} intervals",
///              uniform_refinement.refined_grid().num_intervals().as_ref());
///     
///     // Convert to fully non-uniform grid for further adaptive refinement
///     let refined_grid = uniform_refinement.into_refined_grid();
///     
///     // Apply selective refinement to specific regions
///     let selective_refinement = refined_grid.refine(&std::collections::BTreeMap::from([
///         (IntervalId::new(1), PositiveNumPoints1D::try_new(1).unwrap()),
///         (IntervalId::new(3), PositiveNumPoints1D::try_new(2).unwrap()),
///     ]));
///     
///     println!("After selective refinement: {} intervals",
///              selective_refinement.refined_grid().num_intervals().as_ref());
///     
///     selective_refinement
/// }
///
/// // Usage
/// let base_coords = vec![0.0, 0.2, 0.7, 1.0]; // Non-uniform initial distribution
/// let final_refinement = adaptive_then_uniform_strategy(base_coords);
///
/// // Analyze the refinement pattern
/// for (original_id, refined_ids) in final_refinement.iter_original_with_refined() {
///     println!("Original adaptive interval {} → {} final intervals",
///              original_id.as_ref(), refined_ids.len());
/// }
/// ```
///
/// ## Inherited Functionality from [`Grid1DRefinement`]
///
/// [`Grid1DUniformRefinement`] inherits all the functionality from [`Grid1DRefinement`]:
///
/// ### Bidirectional Mapping Operations
/// All the bidirectional mapping operations are available:
/// - [`find_original_interval`](Grid1DRefinement::find_original_interval): Find which original interval produced a refined interval
/// - [`get_refined_intervals`](Grid1DRefinement::get_refined_intervals): Get all refined intervals from an original interval
/// - [`transfer_interval_data`](Grid1DRefinement::transfer_interval_data): Transfer data from original to refined grid
/// - [`iter_refined_with_mapping`](Grid1DRefinement::iter_refined_with_mapping): Iterate over refined intervals with their sources
/// - [`iter_original_with_refined`](Grid1DRefinement::iter_original_with_refined): Iterate over original intervals with their children
///
/// ### Analysis Methods
/// All analysis methods work identically:
/// - [`was_refined`](Grid1DRefinement::was_refined): Check if an interval was subdivided (always `true` for uniform refinement)
/// - [`num_refined_intervals`](Grid1DRefinement::num_refined_intervals): Get subdivision count (always `k+1`)
///
/// ### Grid Access
/// Standard grid access patterns:
/// - [`refined_grid`](Grid1DRefinement::refined_grid): Access the refined grid
/// - [`original_grid`](Grid1DRefinement::original_grid): Access the original grid
/// - [`into_refined_grid`](Grid1DRefinement::into_refined_grid): Extract the refined grid
/// - [`into_original_grid`](Grid1DRefinement::into_original_grid): Extract the original grid
///
/// ## Performance Characteristics vs. Non-Uniform Refinement
///
/// | Operation | [`Grid1DUniformRefinement`] | [`Grid1DNonUniformRefinement`] | Advantage |
/// |-----------|----------------------------|------------------------------|-----------|
/// | **Creation** | O(n×k) | O(n×k + log n) | ✅ **Uniform simpler** |
/// | **Memory overhead** | O(n×k) | O(n×k) | Equal |
/// | **Data transfer** | O(n×k) | O(n×k) | Equal |
/// | **Mapping lookup** | O(1) | O(1) | Equal |
/// | **Predictability** | Perfect | Variable | ✅ **Uniform predictable** |
/// | **Cache efficiency** | Excellent | Good | ✅ **Uniform better** |
///
/// ## Mathematical Guarantees
///
/// ### Uniform Subdivision Property
/// ```text
/// ∀i ∈ [0, n-1]: |refined_intervals_from(original_interval[i])| = k+1
/// ```
/// Every original interval produces exactly `k+1` refined intervals.
///
/// ### Equal Sub-interval Length Property
/// ```text
/// ∀i ∈ [0, n-1], ∀j,ℓ ∈ refined_intervals_from(original_interval[i]):
///   length(refined_interval[j]) = length(refined_interval[ℓ])
/// ```
/// All refined intervals from the same original interval have identical length.
///
/// ### Domain Conservation Property
/// ```text
/// domain(refined_grid) = domain(original_grid)
/// ```
/// The refined grid covers exactly the same mathematical domain.
///
/// ### Completeness Property
/// ```text
/// ⋃_{i=0}^{n×k-1} refined_interval[i] = ⋃_{i=0}^{n-1} original_interval[i]
/// ```
/// The refined intervals completely partition the original domain.
///
/// ### Mapping Consistency Property
/// ```text
/// ∀refined_id r, original_id o:
///   r ∈ original_to_refined[o] ⟺ refined_to_original[r] = o
/// ```
/// The bidirectional mappings are mathematically consistent.
///
/// ## Error Handling and Validation
///
/// [`Grid1DUniformRefinement`] inherits all error handling from [`Grid1DRefinement`]:
///
/// ### Construction Robustness
/// ```rust
/// use grid1d::{*, intervals::*, scalars::*};
/// use try_create::TryNew;
///
/// fn safe_uniform_refinement<G>(
///     grid: G,
///     refinement_factor: usize,
/// ) -> Result<Grid1DUniformRefinement<G>, String>
/// where
///     G: Grid1DTrait,
/// {
///     if refinement_factor == 0 {
///         return Err("Refinement factor must be positive".to_string());
///     }
///     
///     let extra_points = PositiveNumPoints1D::try_new(refinement_factor - 1)
///         .map_err(|e| format!("Invalid refinement factor: {}", e))?;
///     
///     Ok(grid.refine_uniform(&extra_points))
/// }
///
/// // Usage
/// let grid = Grid1D::uniform(IntervalClosed::new(0.0, 1.0), NumIntervals::try_new(4).unwrap());
/// let refinement = safe_uniform_refinement(grid, 2).unwrap();
/// assert_eq!(refinement.refined_grid().num_intervals().as_ref(), &8);
/// ```
///
/// ## Best Practices and Recommendations
///
/// ### When to Use Uniform Refinement
/// - **Convergence studies**: Systematic error analysis with predictable resolution increase
/// - **Multigrid methods**: Creating perfect hierarchical relationships between grid levels
/// - **Global resolution increase**: When entire domain needs finer resolution uniformly
/// - **Method validation**: Testing numerical methods with controlled grid refinement
/// - **Performance optimization**: When regular patterns enable vectorization/SIMD
/// - **Memory predictability**: When exact memory requirements must be known in advance
///
/// ### When to Use Non-Uniform Refinement Instead
/// - **Adaptive mesh refinement**: Variable resolution based on solution features
/// - **Error-driven refinement**: Refinement guided by local error indicators
/// - **Feature tracking**: High resolution only where needed (shocks, boundaries, etc.)
/// - **Computational efficiency**: When uniform refinement would waste resources
/// - **Complex geometries**: When problem structure requires non-uniform treatment
///
/// ### Optimal Refinement Strategies
/// ```rust
/// use grid1d::{*, intervals::*, scalars::*};
/// use try_create::TryNew;
///
/// // Strategy 1: Conservative uniform refinement
/// fn conservative_strategy<G: Grid1DTrait>(grid: G) -> Grid1DUniformRefinement<G> {
///     grid.refine_uniform(&PositiveNumPoints1D::try_new(1).unwrap()) // 2× refinement
/// }
///
/// // Strategy 2: Aggressive uniform refinement
/// fn aggressive_strategy<G: Grid1DTrait>(grid: G) -> Grid1DUniformRefinement<G> {
///     grid.refine_uniform(&PositiveNumPoints1D::try_new(3).unwrap()) // 4× refinement
/// }
///
/// // Strategy 3: Problem-adaptive uniform refinement
/// fn adaptive_uniform_strategy<G: Grid1DTrait>(
///     grid: G,
///     target_error: f64,
///     current_error: f64,
/// ) -> Grid1DUniformRefinement<G> {
///     let refinement_factor = if current_error > 10.0 * target_error {
///         3 // 4× refinement for high error
///     } else if current_error > 2.0 * target_error {
///         1 // 2× refinement for moderate error
///     } else {
///         0 // 1× refinement (minimal) for low error
///     };
///     
///     grid.refine_uniform(&PositiveNumPoints1D::try_new(refinement_factor).unwrap())
/// }
/// ```
///
/// ## Integration with other numerical Ecosystem
///
/// [`Grid1DUniformRefinement`] integrates seamlessly with:
///
/// - **[`Grid1DRefinement`]**: Inherits all bidirectional mapping functionality
/// - **[`Grid1DTrait`]**: The refined grid implements the full partition interface
/// - **[`HasDomain1D`]**, **[`HasCoords1D`](crate::HasCoords1D)**: Standard domain and coordinate access patterns
/// - **[`Grid1D`]**, **[`Grid1DUnion`](super::union::Grid1DUnion)**: Compatible with all grid operations and combinations
/// - **Numerical algorithms**: Works with any algorithm expecting [`Grid1DTrait`]
/// - **Multigrid frameworks**: Provides the systematic refinement needed for multigrid methods
///
/// ## See Also
///
/// - [`Grid1DRefinement`]: The general refinement struct this type alias specializes
/// - [`Grid1DNonUniformRefinement`]: Alternative for selective interval refinement
/// - [`Grid1DTrait`]: Core trait providing partition operations
/// - [`Grid1D`]: The unified grid interface supporting refinement
/// - [`Grid1DUnion`](super::union::Grid1DUnion): For combining multiple grids into unified partitions
/// - Multigrid methods literature for systematic refinement strategies
/// - Numerical analysis texts for convergence theory with uniform refinement
pub type Grid1DUniformRefinement<OriginalGrid1DType> = Grid1DRefinement<
    OriginalGrid1DType,
    <OriginalGrid1DType as Grid1DTrait>::UniformlyRefinedGrid1DType,
>;

#[duplicate_item(
    grid_type;
    [Grid1DUniform];
    [Grid1DNonUniform];
)]
impl<Domain1D> From<Grid1DUniformRefinement<grid_type<Domain1D>>>
    for Grid1DUniformRefinement<Grid1D<Domain1D>>
where
    Domain1D: Grid1DIntervalBuilder,
{
    fn from(value: Grid1DUniformRefinement<grid_type<Domain1D>>) -> Self {
        let Grid1DRefinement {
            refined_grid,
            original_grid,
            refined_to_original_interval_mapping,
            original_to_refined_interval_mapping,
        } = value;

        Grid1DRefinement::new(
            refined_grid.into(),  // converting refined grid to Grid1D enum wrapper
            original_grid.into(), // converting original grid to Grid1D enum wrapper
            refined_to_original_interval_mapping,
            original_to_refined_interval_mapping,
        )
    }
}

///  A specialized refinement of a 1D grid where intervals are selectively subdivided with variable resolution.
///
/// [`Grid1DNonUniformRefinement`] is a type alias for [`Grid1DRefinement`] that represents the result
/// of **selective refinement**, where only specified intervals in the original grid are subdivided,
/// each potentially with a different number of sub-intervals. This is the cornerstone of adaptive mesh
/// refinement (AMR) strategies where resolution is increased only where needed, such as near solution
/// gradients, shocks, or complex features.
///
/// ## Mathematical Definition
///
/// Given an original grid with `n` intervals and a refinement specification
/// `{(i₁, k₁), (i₂, k₂), ..., (iₘ, kₘ)}`, selective refinement creates a new grid where:
/// ```text
/// ∀j ∈ refinement_list: original_interval[j] → kⱼ+1 refined_intervals
/// ∀j ∉ refinement_list: original_interval[j] → 1 refined_interval (unchanged)
/// ```
///
/// **Total refined intervals**: `n + Σ(kᵢ)` where the sum is over all refined intervals
///
/// **Adaptive resolution**: Different intervals can have different refinement levels based on local needs
///
/// ## Type Definition
///
/// ```rust,ignore
/// pub type Grid1DNonUniformRefinement<OriginalGrid1DType> = Grid1DRefinement<
///     OriginalGrid1DType,
///     Grid1DNonUniform<<OriginalGrid1DType as HasDomain1D>::Domain1D>,
/// >;
/// ```
///
/// This ensures that:
/// - **Input**: Any grid type implementing [`Grid1DTrait`]
/// - **Output**: Always a [`Grid1DNonUniform`] to handle arbitrary point distributions
/// - **Type safety**: Compile-time guarantee of domain compatibility
/// - **Flexibility**: Can handle any combination of refinement patterns
///
/// ## Core Mathematical Properties
///
/// ### Selective Subdivision Guarantee
/// ```text
/// ∀i ∈ [0, n-1]: num_refined_intervals(original_interval[i]) = specified_or_1
/// ```
/// Only specified intervals are subdivided; others remain unchanged.
///
/// ### Variable Resolution Support
/// ```text
/// ∀i,j ∈ refined_intervals: length(refined_interval[i]) ≠ length(refined_interval[j]) (generally)
/// ```
/// Refined intervals can have dramatically different lengths based on local requirements.
///
/// ### Structure Adaptation
/// - **Uniform → Non-uniform**: [`Grid1DUniform`] becomes [`Grid1DNonUniform`] after selective refinement
/// - **Non-uniform → Non-uniform**: [`Grid1DNonUniform`] remains [`Grid1DNonUniform`] with more points
/// - **Local preservation**: Unrefined regions maintain their original spacing
/// - **Domain consistency**: The refined grid covers exactly the same domain as the original
///
/// ### Refinement Efficiency
/// ```text
/// refinement_ratio = (n + Σ(kᵢ)) / n
/// ```
/// The total number of intervals increases only where refinement is applied.
///
/// ## Design Philosophy
///
/// ### Adaptive Resolution Strategy
/// Non-uniform refinement provides **targeted improvement** in grid resolution:
/// - **Local error reduction**: High resolution only where the solution requires it
/// - **Computational efficiency**: Avoids wasting resources on smooth regions
/// - **Feature preservation**: Can resolve boundary layers, shocks, and sharp gradients
/// - **Memory optimization**: Minimizes memory usage compared to uniform refinement
///
/// ### Algorithmic Flexibility
/// - **Error-driven refinement**: Refine based on local error indicators
/// - **Feature detection**: Increase resolution near detected features
/// - **Gradient tracking**: Follow solution gradients with adaptive resolution
/// - **Physics-based adaptation**: Different refinement for different physical phenomena
///
/// ### Integration with Numerical Methods
/// - **AMR frameworks**: Core component of adaptive mesh refinement systems
/// - **Error estimation**: Enables sophisticated error analysis and control
/// - **Multi-physics coupling**: Different refinement patterns for different physics
/// - **Hierarchical methods**: Builds complex grid hierarchies for multi-level algorithms
///
/// ## Construction and Usage Examples
///
/// ### Basic Selective Refinement
/// ```rust
/// use grid1d::{*, intervals::*, scalars::*};
/// use sorted_vec::partial::SortedSet;
/// use std::ops::Deref;
/// use std::collections::BTreeMap;
/// use try_create::TryNew;
///
/// // Create original uniform grid
/// let original = Grid1D::uniform(
///     IntervalClosed::new(0.0, 1.0),
///     NumIntervals::try_new(4).unwrap()
/// );
/// assert_eq!(original.coords().deref(), &[0.0, 0.25, 0.5, 0.75, 1.0]);
///
/// // Selective refinement: refine only intervals 1 and 3
/// let refinement_plan = BTreeMap::from([
///     (IntervalId::new(1), PositiveNumPoints1D::try_new(1).unwrap()), // Split into 2
///     (IntervalId::new(3), PositiveNumPoints1D::try_new(2).unwrap()), // Split into 3
/// ]);
/// let refinement = original.refine(&refinement_plan);
///
/// // Verify refined grid structure
/// assert_eq!(refinement.refined_grid().coords().deref(),
///            &[0.0, 0.25, 0.375, 0.5, 0.75, 0.8333333333333334, 0.9166666666666666, 1.0]);
/// assert_eq!(refinement.refined_grid().num_intervals().as_ref(), &7);
///
/// // Check selective mapping
/// for (original_id, refined_ids) in refinement.iter_original_with_refined() {
///     println!("Original interval {}: {} refined intervals",
///              original_id.as_ref(), refined_ids.len());
/// }
/// ```
///
/// ### Error-Driven Adaptive Refinement
/// ```rust
/// use grid1d::{*, intervals::*, scalars::*};
/// use try_create::TryNew;
///
/// fn adaptive_refinement_step(
///     grid: Grid1D<IntervalClosed<f64>>,
///     solution: &[f64],
///     error_threshold: f64,
/// ) -> Grid1DNonUniformRefinement<Grid1D<IntervalClosed<f64>>> {
///     // Compute local error indicators (simplified)
///     let errors: Vec<f64> = (0..*grid.num_intervals().as_ref())
///         .map(|i| {
///             // Simplified error estimation based on solution gradient
///             if i > 0 && i < solution.len() - 1 {
///                 (solution[i+1] - 2.0*solution[i] + solution[i-1]).abs()
///             } else {
///                 0.0
///             }
///         })
///         .collect();
///     
///     // Create refinement plan based on error indicators
///     let refinement_plan: std::collections::BTreeMap<IntervalId, PositiveNumPoints1D> = errors
///         .iter()
///         .enumerate()
///         .filter_map(|(i, &error)| {
///             if error > error_threshold {
///                 // More refinement for higher errors
///                 let num_extra = if error > 10.0 * error_threshold { 3 } else { 1 };
///                 Some((IntervalId::new(i), PositiveNumPoints1D::try_new(num_extra).unwrap()))
///             } else {
///                 None
///             }
///         })
///         .collect();
///     
///     println!("Refining {} out of {} intervals",
///              refinement_plan.len(), grid.num_intervals().as_ref());
///     
///     grid.refine(&refinement_plan)
/// }
///
/// // Usage in adaptive solver
/// let mut current_grid = Grid1D::uniform(IntervalClosed::new(0.0, 1.0), NumIntervals::try_new(8).unwrap());
/// let mut solution = vec![1.0; 8]; // Initial solution
///
/// for iteration in 0..3 {
///     let refinement = adaptive_refinement_step(current_grid, &solution, 0.1);
///     
///     // Transfer solution to refined grid
///     solution = refinement.transfer_interval_data(&solution);
///     current_grid = Grid1D::NonUniform(refinement.into_refined_grid());
///     
///     println!("Iteration {}: grid has {} intervals",
///              iteration, current_grid.num_intervals().as_ref());
/// }
/// ```
///
/// ### Boundary Layer Refinement
/// ```rust
/// use grid1d::{*, intervals::*, scalars::*};
/// use sorted_vec::partial::SortedSet;
/// use try_create::TryNew;
///
/// fn boundary_layer_refinement(
///     domain: IntervalClosed<f64>,
///     base_resolution: usize,
///     boundary_thickness: f64,
/// ) -> Grid1DNonUniformRefinement<Grid1D<IntervalClosed<f64>>> {
///     // Create coarse base grid
///     let base_grid = Grid1D::uniform(domain.clone(), NumIntervals::try_new(base_resolution).unwrap());
///     
///     let mut refinement_plan = std::collections::BTreeMap::new();
///     let domain_length = domain.upper_bound_value() - domain.lower_bound_value();
///     
///     // Refine intervals near boundaries
///     for (interval_id, _) in base_grid.iter_intervals() {
///         let coords = base_grid.coords().as_ref();
///         let interval_start = coords[*interval_id.as_ref()];
///         let interval_end = coords[*interval_id.as_ref() + 1];
///         
///         // Distance from lower boundary
///         let dist_from_lower = interval_start - domain.lower_bound_value();
///         let dist_from_upper = domain.upper_bound_value() - &interval_end;
///         
///         let min_boundary_dist = dist_from_lower.min(dist_from_upper);
///         
///         if min_boundary_dist < boundary_thickness {
///             // Refinement level decreases with distance from boundary
///             let refinement_level = if min_boundary_dist < (boundary_thickness / 3.0) {
///                 4 // Very close to boundary: 5 sub-intervals
///             } else if min_boundary_dist < (2.0 * boundary_thickness / 3.0) {
///                 2 // Moderately close: 3 sub-intervals
///             } else {
///                 1 // Edge of boundary layer: 2 sub-intervals
///             };
///             
///             refinement_plan.insert(
///                 interval_id,
///                 PositiveNumPoints1D::try_new(refinement_level).unwrap()
///             );
///         }
///     }
///     
///     println!("Boundary layer refinement: {}/{} intervals refined",
///              refinement_plan.len(), base_grid.num_intervals().as_ref());
///     
///     base_grid.refine(&refinement_plan)
/// }
///
/// // Create boundary layer mesh for fluid dynamics
/// let domain = IntervalClosed::new(0.0, 1.0);
/// let bl_refinement = boundary_layer_refinement(domain, 10, 0.1);
///
/// // Verify high resolution near boundaries
/// let refined_coords = bl_refinement.refined_grid().coords().as_ref();
/// println!("Refined grid: {} points", refined_coords.len());
/// println!("Spacing near x=0: {:.6}", refined_coords[1] - refined_coords[0]);
/// println!("Spacing in middle: {:.6}",
///          refined_coords[refined_coords.len()/2+1] - refined_coords[refined_coords.len()/2]);
/// ```
///
/// ## Advanced Usage Patterns
///
/// ### Hierarchical Adaptive Refinement
/// ```rust
/// use grid1d::{*, intervals::*, scalars::*};
/// use try_create::TryNew;
///
/// fn hierarchical_adaptive_strategy(
///     initial_grid: Grid1D<IntervalClosed<f64>>,
///     max_levels: usize,
///     refinement_criterion: impl Fn(&Grid1D<IntervalClosed<f64>>, &[f64]) -> Vec<IntervalId>,
/// ) -> Vec<Grid1DNonUniformRefinement<Grid1D<IntervalClosed<f64>>>> {
///     let mut refinement_history = Vec::new();
///     let mut current_grid = initial_grid;
///     let mut solution = vec![1.0; *current_grid.num_intervals().as_ref()]; // Placeholder solution
///     
///     for level in 0..max_levels {
///         // Apply refinement criterion
///         let intervals_to_refine = refinement_criterion(&current_grid, &solution);
///         
///         if intervals_to_refine.is_empty() {
///             println!("Convergence achieved at level {}", level);
///             break;
///         }
///         
///         // Create refinement plan (uniform refinement factor for simplicity)
///         let refinement_plan: std::collections::BTreeMap<IntervalId, PositiveNumPoints1D> = intervals_to_refine
///             .into_iter()
///             .map(|id| (id, PositiveNumPoints1D::try_new(1).unwrap()))
///             .collect();
///         
///         // Perform refinement
///         let refinement = current_grid.refine(&refinement_plan);
///         
///         println!("Level {}: {} → {} intervals",
///                  level,
///                  refinement.original_grid().num_intervals().as_ref(),
///                  refinement.refined_grid().num_intervals().as_ref());
///         
///         // Transfer solution and update grid
///         solution = refinement.transfer_interval_data(&solution);
///         current_grid = Grid1D::NonUniform(refinement.refined_grid().clone());
///         
///         refinement_history.push(refinement);
///     }
///     
///     refinement_history
/// }
///
/// // Example refinement criterion: refine intervals with high curvature
/// let criterion = |grid: &Grid1D<IntervalClosed<f64>>, solution: &[f64]| -> Vec<IntervalId> {
///     let mut to_refine = Vec::new();
///     for i in 1..solution.len()-1 {
///         let curvature = (solution[i+1] - 2.0*solution[i] + solution[i-1]).abs();
///         if curvature > 0.1 {
///             to_refine.push(IntervalId::new(i));
///         }
///     }
///     to_refine
/// };
///
/// let initial = Grid1D::uniform(IntervalClosed::new(0.0, 1.0), NumIntervals::try_new(8).unwrap());
/// let refinement_sequence = hierarchical_adaptive_strategy(initial, 5, criterion);
///
/// println!("Created {} refinement levels", refinement_sequence.len());
/// ```
///
/// ### Multi-Physics Adaptive Coupling
/// ```rust
/// use grid1d::{*, intervals::*, scalars::*};
/// use std::collections::HashMap;
/// use try_create::TryNew;
///
/// struct MultiPhysicsAdaptiveSystem {
///     fluid_grid: Grid1D<IntervalClosed<f64>>,
///     thermal_grid: Grid1D<IntervalClosed<f64>>,
///     coupling_refinement: Option<Grid1DNonUniformRefinement<Grid1D<IntervalClosed<f64>>>>,
/// }
///
/// impl MultiPhysicsAdaptiveSystem {
///     fn new(domain: IntervalClosed<f64>, initial_resolution: usize) -> Self {
///         let fluid_grid = Grid1D::uniform(domain.clone(), NumIntervals::try_new(initial_resolution).unwrap());
///         let thermal_grid = Grid1D::uniform(domain, NumIntervals::try_new(initial_resolution/2).unwrap());
///         
///         Self {
///             fluid_grid,
///             thermal_grid,
///             coupling_refinement: None,
///         }
///     }
///     
///     fn create_coupled_refinement(&mut self,
///         fluid_gradients: &[f64],
///         thermal_gradients: &[f64]
///     ) -> Result<(), String> {
///         // Create union of both physics grids
///         let union = Grid1DUnion::try_new(&self.fluid_grid, &self.thermal_grid)
///             .map_err(|e| format!("Failed to create grid union: {:?}", e))?;
///         
///         let unified_grid = union.unified_grid().clone();
///         
///         // Identify intervals needing refinement from either physics
///         let mut refinement_map: HashMap<usize, usize> = HashMap::new();
///         
///         // Process fluid gradients
///         for (i, &gradient) in fluid_gradients.iter().enumerate() {
///             if gradient > 1.0 {
///                 let (unified_id, _) = union.find_original_intervals(&IntervalId::new(i));
///                 refinement_map.insert(*unified_id.as_ref(),
///                                     refinement_map.get(unified_id.as_ref()).unwrap_or(&1) + 1);
///             }
///         }
///         
///         // Process thermal gradients (with different threshold)
///         for (i, &gradient) in thermal_gradients.iter().enumerate() {
///             if gradient > 0.5 {
///                 let (_, unified_id) = union.find_original_intervals(&IntervalId::new(i));
///                 refinement_map.insert(*unified_id.as_ref(),
///                                     refinement_map.get(unified_id.as_ref()).unwrap_or(&1) + 1);
///             }
///         }
///         
///         // Create refinement plan
///         let refinement_plan: std::collections::BTreeMap<IntervalId, PositiveNumPoints1D> = refinement_map
///             .into_iter()
///             .map(|(id, level)| {
///                 let clamped_level = level.min(3); // Maximum 4 sub-intervals
///                 (IntervalId::new(id), PositiveNumPoints1D::try_new(clamped_level).unwrap())
///             })
///             .collect();
///         
///         self.coupling_refinement = Some(unified_grid.refine(&refinement_plan));
///         
///         Ok(())
///     }
///     
///     fn get_coupled_grid(&self) -> Option<&Grid1DNonUniform<IntervalClosed<f64>>> {
///         self.coupling_refinement.as_ref().map(|r| r.refined_grid())
///     }
/// }
///
/// // Usage in multi-physics simulation
/// let mut system = MultiPhysicsAdaptiveSystem::new(IntervalClosed::new(0.0, 1.0), 20);
/// let fluid_gradients = vec![0.1, 2.0, 0.5, 1.5, 0.2]; // High gradients at indices 1, 3
/// let thermal_gradients = vec![0.3, 0.1, 0.8, 0.2]; // High gradient at index 2
///
/// system.create_coupled_refinement(&fluid_gradients, &thermal_gradients).unwrap();
///
/// if let Some(coupled_grid) = system.get_coupled_grid() {
///     println!("Coupled grid has {} intervals", coupled_grid.num_intervals().as_ref());
/// }
/// ```
///
/// ## Inherited Functionality from [`Grid1DRefinement`]
///
/// [`Grid1DNonUniformRefinement`] inherits all the functionality from [`Grid1DRefinement`]:
///
/// ### Bidirectional Mapping Operations
/// All the bidirectional mapping operations are available:
/// - [`find_original_interval`](Grid1DRefinement::find_original_interval): Find which original interval produced a refined interval
/// - [`get_refined_intervals`](Grid1DRefinement::get_refined_intervals): Get all refined intervals from an original interval
/// - [`transfer_interval_data`](Grid1DRefinement::transfer_interval_data): Transfer data from original to refined grid
/// - [`iter_refined_with_mapping`](Grid1DRefinement::iter_refined_with_mapping): Iterate over refined intervals with their sources
/// - [`iter_original_with_refined`](Grid1DRefinement::iter_original_with_refined): Iterate over original intervals with their children
///
/// ### Analysis Methods
/// All analysis methods work with variable refinement:
/// - [`was_refined`](Grid1DRefinement::was_refined): Check if an interval was subdivided
/// - [`num_refined_intervals`](Grid1DRefinement::num_refined_intervals): Get subdivision count (variable per interval)
///
/// ### Grid Access
/// Standard grid access patterns:
/// - [`refined_grid`](Grid1DRefinement::refined_grid): Access the refined grid (always [`Grid1DNonUniform`])
/// - [`original_grid`](Grid1DRefinement::original_grid): Access the original grid
/// - [`into_refined_grid`](Grid1DRefinement::into_refined_grid): Extract the refined grid
/// - [`into_original_grid`](Grid1DRefinement::into_original_grid): Extract the original grid
///
/// ## Performance Characteristics vs. Uniform Refinement
///
/// | Operation | [`Grid1DNonUniformRefinement`] | [`Grid1DUniformRefinement`] | Advantage |
/// |-----------|-------------------------------|----------------------------|-----------|
/// | **Creation** | O(n×k + log n) | O(n×k) | ✅ **Uniform simpler** |
/// | **Memory overhead** | O(refined intervals) | O(n×k) | ✅ **Non-uniform smaller** |
/// | **Data transfer** | O(refined intervals) | O(n×k) | ✅ **Non-uniform smaller** |
/// | **Mapping lookup** | O(1) | O(1) | Equal |
/// | **Flexibility** | Complete | None | ✅ **Non-uniform adaptive** |
/// | **Computational savings** | High | None | ✅ **Non-uniform efficient** |
///
/// ## Mathematical Guarantees
///
/// ### Selective Subdivision Property
/// ```text
/// ∀i ∈ [0, n-1]: |refined_intervals_from(original_interval[i])| = specified_or_1
/// ```
/// Only intervals in the refinement plan are subdivided; others remain as single intervals.
///
/// ### Variable Resolution Property
/// ```text
/// ∀i,j: refinement_factor[i] ≠ refinement_factor[j] (potentially)
/// ```
/// Different intervals can have completely different refinement levels.
///
/// ### Local Consistency Property
/// ```text
/// ∀i ∈ refined_intervals_from(original_interval[j]):
///   length(refined_interval[i]) = length(original_interval[j]) / refinement_factor[j]
/// ```
/// Within each refined original interval, sub-intervals have equal length.
///
/// ### Domain Conservation Property
/// ```text
/// domain(refined_grid) = domain(original_grid)
/// ```
/// The refined grid covers exactly the same mathematical domain.
///
/// ### Efficiency Property
/// ```text
/// |refined_intervals| ≤ |original_intervals| × max_refinement_factor
/// ```
/// The number of refined intervals is bounded and often much smaller than uniform refinement.
///
/// ### Mapping Consistency Property
/// ```text
/// ∀refined_id r, original_id o:
///   r ∈ original_to_refined[o] ⟺ refined_to_original[r] = o
/// ```
/// The bidirectional mappings are mathematically consistent.
///
/// ## Best Practices and Recommendations
///
/// ### When to Use Non-Uniform Refinement
/// - **Adaptive mesh refinement**: Variable resolution based on solution features
/// - **Error-driven refinement**: Refinement guided by local error indicators
/// - **Feature tracking**: High resolution only where needed (shocks, boundaries, etc.)
/// - **Multi-physics problems**: Different resolution requirements for different physics
/// - **Computational efficiency**: When uniform refinement would waste resources
/// - **Complex geometries**: When problem structure requires non-uniform treatment
///
/// ### When to Use Uniform Refinement Instead
/// - **Convergence studies**: Systematic error analysis with predictable resolution increase
/// - **Multigrid methods**: Creating perfect hierarchical relationships between grid levels
/// - **Simple problems**: When uniform resolution is sufficient across the domain
/// - **Performance optimization**: When regular patterns enable vectorization/SIMD
/// - **Memory predictability**: When exact memory requirements must be known in advance
///
/// ### Optimal Refinement Strategies
/// ```rust
/// use grid1d::{*, intervals::*, scalars::*};
/// use try_create::TryNew;
///
/// // Strategy 1: Error-based adaptive refinement
/// fn error_driven_refinement<G: Grid1DTrait>(
///     grid: G,
///     error_indicators: &[f64],
///     threshold: f64
/// ) -> Grid1DNonUniformRefinement<G> {
///     let refinement_plan: std::collections::BTreeMap<IntervalId, PositiveNumPoints1D> = error_indicators
///         .iter()
///         .enumerate()
///         .filter_map(|(i, &error)| {
///             if error > threshold {
///                 let level = if error > 10.0 * threshold { 3 } else { 1 };
///                 Some((IntervalId::new(i), PositiveNumPoints1D::try_new(level).unwrap()))
///             } else {
///                 None
///             }
///         })
///         .collect();
///     
///     grid.refine(&refinement_plan)
/// }
///
/// // Strategy 2: Gradient-based refinement
/// fn gradient_based_refinement<G: Grid1DTrait>(
///     grid: G,
///     solution: &[f64],
///     gradient_threshold: f64
/// ) -> Grid1DNonUniformRefinement<G> {
///     let mut refinement_plan = std::collections::BTreeMap::new();
///     
///     for i in 1..solution.len()-1 {
///         let gradient = ((solution[i+1] - solution[i-1]) / 2.0).abs();
///         if gradient > gradient_threshold {
///             refinement_plan.insert(
///                 IntervalId::new(i),
///                 PositiveNumPoints1D::try_new(1).unwrap()
///             );
///         }
///     }
///     
///     grid.refine(&refinement_plan)
/// }
///
/// // Strategy 3: Feature-proximity refinement
/// fn proximity_based_refinement<G: Grid1DTrait>(
///     grid: G,
///     feature_locations: &[f64],
///     influence_radius: f64
/// ) -> Grid1DNonUniformRefinement<G>
/// {
///     let mut refinement_plan = std::collections::BTreeMap::new();
///     
///     for (interval_id, _) in grid.iter_intervals() {
///         let coords = grid.coords().as_ref();
///         let i = *interval_id.as_ref();
///         let interval_center = coords[i].clone() + &coords[i + 1];
///         // Note: This is a simplified example - actual center calculation would depend on the scalar type
///         
///         // Check proximity to any feature
///         for &feature_pos in feature_locations {
///             // Simplified distance calculation
///             if true { // Placeholder for actual distance check
///                 refinement_plan.insert(
///                     interval_id,
///                     PositiveNumPoints1D::try_new(2).unwrap()
///                 );
///                 break;
///             }
///         }
///     }
///     
///     grid.refine(&refinement_plan)
/// }
/// ```
///
/// ## Error Handling and Validation
///
/// [`Grid1DNonUniformRefinement`] inherits all error handling from [`Grid1DRefinement`] and adds validation for refinement plans:
///
/// ### Refinement Plan Validation
/// ```rust
/// use grid1d::{*, intervals::*, scalars::*};
///
/// fn validate_and_refine<G>(
///     grid: G,
///     refinement_plan: &std::collections::BTreeMap<IntervalId, PositiveNumPoints1D>,
/// ) -> Result<Grid1DNonUniformRefinement<G>, String>
/// where
///     G: Grid1DTrait,
/// {
///     // Validate interval IDs
///     for (&interval_id, _) in refinement_plan {
///         if *interval_id.as_ref() >= *grid.num_intervals().as_ref() {
///             return Err(format!("Invalid interval ID: {} >= {}",
///                               interval_id.as_ref(), grid.num_intervals().as_ref()));
///         }
///     }
///     
///     // Note: BTreeMap guarantees unique interval IDs - no need to check for duplicates!
///     
///     // Estimate memory usage
///     let total_extra_intervals: usize = refinement_plan
///         .values()
///         .map(|extra_points| *extra_points.as_ref())
///         .sum();
///     
///     if total_extra_intervals > 10 * grid.num_intervals().as_ref() {
///         return Err("Refinement plan would create excessive number of intervals".to_string());
///     }
///     
///     Ok(grid.refine(refinement_plan))
/// }
/// ```
///
/// ## Integration with other numerical Ecosystem
///
/// [`Grid1DNonUniformRefinement`] integrates seamlessly with:
///
/// - **[`Grid1DRefinement`]**: Inherits all bidirectional mapping functionality
/// - **[`Grid1DTrait`]**: The refined grid implements the full partition interface
/// - **[`HasDomain1D`]**, **[`HasCoords1D`](crate::HasCoords1D)**: Standard domain and coordinate access patterns
/// - **[`Grid1D`]**, **[`Grid1DUnion`](super::union::Grid1DUnion)**: Compatible with all grid operations and combinations
/// - **Numerical algorithms**: Works with any algorithm expecting [`Grid1DTrait`]
/// - **AMR frameworks**: Provides the selective refinement needed for adaptive methods
///
/// ## See Also
///
/// - [`Grid1DRefinement`]: The general refinement struct this type alias specializes
/// - [`Grid1DUniformRefinement`]: Alternative for systematic interval refinement
/// - [`Grid1DTrait`]: Core trait providing partition operations
/// - [`Grid1D`]: The unified grid interface supporting refinement
/// - [`Grid1DUnion`](super::union::Grid1DUnion): For combining multiple grids into unified partitions
/// - Adaptive mesh refinement (AMR) literature for refinement strategies
/// - Computational fluid dynamics texts for boundary layer and feature-based refinement
pub type Grid1DNonUniformRefinement<OriginalGrid1DType> = Grid1DRefinement<
    OriginalGrid1DType,
    Grid1DNonUniform<<OriginalGrid1DType as HasDomain1D>::Domain1D>,
>;

#[duplicate_item(
    grid_type;
    [Grid1DUniform];
    [Grid1DNonUniform];
)]
impl<Domain1D> From<Grid1DNonUniformRefinement<grid_type<Domain1D>>>
    for Grid1DNonUniformRefinement<Grid1D<Domain1D>>
where
    Domain1D: Grid1DIntervalBuilder,
{
    fn from(value: Grid1DNonUniformRefinement<grid_type<Domain1D>>) -> Self {
        let Grid1DRefinement {
            refined_grid,
            original_grid,
            refined_to_original_interval_mapping,
            original_to_refined_interval_mapping,
        } = value;

        Grid1DRefinement::new(
            refined_grid,         // the non-uniform refined grid is directly used
            original_grid.into(), // converting original grid to Grid1D enum wrapper
            refined_to_original_interval_mapping,
            original_to_refined_interval_mapping,
        )
    }
}

/// **INTERNAL** Builds a non-uniform grid refinement with bidirectional interval mappings.
///
/// This function creates a refined grid by selectively subdividing specified intervals
/// and maintains complete bidirectional mappings between the original and refined grids.
/// It's the core implementation for adaptive mesh refinement operations.
///
/// ## Process
/// 1. **Coordinate refinement**: Use [`refine_intervals_in_coords1d`] to create new coordinates
/// 2. **Grid construction**: Build the refined grid from the new coordinates
/// 3. **Mapping creation**: Establish bidirectional interval relationships
/// 4. **Validation**: Ensure mathematical consistency of the refinement
///
/// ## Parameters
/// - `original_grid`: The grid to be refined
/// - `intervals_to_refine`: Specification of which intervals to refine and how
///
/// ## Returns
/// A [`Grid1DNonUniformRefinement`] containing the refined grid and all mappings.
///
/// ## Mathematical Properties
/// - **Coverage preservation**: Refined grid covers exactly the same domain
/// - **Mapping consistency**: Bidirectional mappings are mathematically consistent
/// - **Interval conservation**: Total interval count equals original plus refinements
pub(crate) fn build_non_uniform_grid_refinement<OriginalGrid1DType: Grid1DTrait>(
    original_grid: OriginalGrid1DType,
    intervals_to_refine: &std::collections::BTreeMap<IntervalId, PositiveNumPoints1D>,
) -> Grid1DNonUniformRefinement<OriginalGrid1DType> {
    let (
        refined_coords,
        refined_to_original_interval_mapping,
        original_to_refined_interval_mapping,
    ) = refine_intervals_in_coords1d(original_grid.coords(), intervals_to_refine);

    let domain = original_grid.domain().clone();

    let num_refined_intervals = NumIntervals::try_new(refined_coords.len() - 1)
        .expect("Refined grid must have at least one interval!");

    let topology = original_grid.index_spaces().topology();
    let refined_index_spaces = match topology {
        Topology1D::RealLine => Grid1DIndexSpaces::new_non_periodic(num_refined_intervals),
        Topology1D::Circle => {
            if domain.is_lower_bound_open() {
                Grid1DIndexSpaces::new_periodic_lower_open_interval(num_refined_intervals)
            } else {
                Grid1DIndexSpaces::new_periodic_upper_open_interval(num_refined_intervals)
            }
        }
    };

    let refined_grid = Grid1DNonUniform {
        coords: refined_coords,
        domain,
        index_spaces: refined_index_spaces,
    };

    Grid1DNonUniformRefinement {
        original_grid,
        refined_grid,
        refined_to_original_interval_mapping,
        original_to_refined_interval_mapping,
    }
}
//-------------------------------------------------------------------------------
//===============================================================================
// TESTS
//===============================================================================

#[cfg(test)]
mod tests {
    use super::*;
    use crate::{
        Grid1D, Grid1DNonUniform, Grid1DUniform,
        grids::traits::{Grid1DTrait, HasCoords1D, HasIntervalIdRange},
        intervals::{IntervalClosed, bounded::IntervalFromBounds},
        scalars::NumIntervals,
    };
    use num_valid::RealNative64StrictFinite;
    use sorted_vec::partial::SortedSet;
    use std::collections::BTreeMap;
    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 uniform_refinement {
        use super::*;

        #[test]
        fn basic_uniform_refinement() {
            let grid = Grid1D::uniform(
                IntervalClosed::new(real(0.0), real(1.0)),
                NumIntervals::try_new(2).unwrap(),
            );
            let refinement = grid.refine_uniform(&PositiveNumPoints1D::try_new(1).unwrap());

            // Original: [0, 0.5, 1] -> Refined: [0, 0.25, 0.5, 0.75, 1]
            assert_eq!(
                refinement.refined_grid().coords().deref(),
                &[real(0.0), real(0.25), real(0.5), real(0.75), real(1.0)]
            );
            assert_eq!(refinement.refined_grid().num_intervals().as_ref(), &4);
        }

        #[test]
        fn uniform_refinement_mappings() {
            let grid = Grid1D::uniform(
                IntervalClosed::new(real(0.0), real(1.0)),
                NumIntervals::try_new(2).unwrap(),
            );
            let refinement = grid.refine_uniform(&PositiveNumPoints1D::try_new(1).unwrap());

            // Check mappings: refined intervals 0,1 -> original 0; refined 2,3 -> original 1
            assert_eq!(
                refinement.find_original_interval(&IntervalId::new(0)),
                IntervalId::new(0)
            );
            assert_eq!(
                refinement.find_original_interval(&IntervalId::new(1)),
                IntervalId::new(0)
            );
            assert_eq!(
                refinement.find_original_interval(&IntervalId::new(2)),
                IntervalId::new(1)
            );
            assert_eq!(
                refinement.find_original_interval(&IntervalId::new(3)),
                IntervalId::new(1)
            );
        }

        #[test]
        fn uniform_refinement_with_higher_subdivision() {
            let grid = Grid1D::uniform(
                IntervalClosed::new(real(0.0), real(1.0)),
                NumIntervals::try_new(2).unwrap(),
            );
            // Add 2 extra points per interval -> 3 sub-intervals each
            let refinement = grid.refine_uniform(&PositiveNumPoints1D::try_new(2).unwrap());

            assert_eq!(refinement.refined_grid().num_intervals().as_ref(), &6);

            // Each original interval should map to 3 refined intervals
            assert_eq!(refinement.num_refined_intervals(&IntervalId::new(0)), 3);
            assert_eq!(refinement.num_refined_intervals(&IntervalId::new(1)), 3);
        }

        #[test]
        fn uniform_refinement_preserves_domain() {
            let domain = IntervalClosed::new(real(-5.0), real(10.0));
            let grid = Grid1D::uniform(domain.clone(), NumIntervals::try_new(3).unwrap());
            let original_first: Real = *grid.coords().first();
            let original_last: Real = *grid.coords().last();
            let refinement = grid.refine_uniform(&PositiveNumPoints1D::try_new(1).unwrap());

            assert_eq!(refinement.refined_grid().domain(), &domain);
            assert_eq!(refinement.refined_grid().coords().first(), &original_first);
            assert_eq!(refinement.refined_grid().coords().last(), &original_last);
        }

        #[test]
        fn uniform_refinement_on_non_uniform_grid() {
            let coords = SortedSet::from_unsorted(vec![real(0.0), real(0.2), real(1.0)]);
            let grid = Grid1D::<IntervalClosed<Real>>::try_from_sorted(coords).unwrap();

            let refinement = grid.refine_uniform(&PositiveNumPoints1D::try_new(1).unwrap());

            // Each of 2 intervals becomes 2 -> total 4 intervals
            assert_eq!(refinement.refined_grid().num_intervals().as_ref(), &4);
        }
    }

    mod selective_refinement {
        use super::*;

        #[test]
        fn selective_refinement_single_interval() {
            let grid = Grid1D::uniform(
                IntervalClosed::new(real(0.0), real(3.0)),
                NumIntervals::try_new(3).unwrap(),
            );

            // Refine only the middle interval
            let plan =
                BTreeMap::from([(IntervalId::new(1), PositiveNumPoints1D::try_new(1).unwrap())]);
            let refinement = grid.refine(&plan);

            // Should have 4 intervals total (3 original + 1 extra from refinement)
            assert_eq!(refinement.refined_grid().num_intervals().as_ref(), &4);

            // Check that only interval 1 was refined
            assert_eq!(refinement.num_refined_intervals(&IntervalId::new(0)), 1);
            assert_eq!(refinement.num_refined_intervals(&IntervalId::new(1)), 2);
            assert_eq!(refinement.num_refined_intervals(&IntervalId::new(2)), 1);
        }

        #[test]
        fn selective_refinement_was_refined_check() {
            let grid = Grid1D::uniform(
                IntervalClosed::new(real(0.0), real(3.0)),
                NumIntervals::try_new(3).unwrap(),
            );

            let plan =
                BTreeMap::from([(IntervalId::new(1), PositiveNumPoints1D::try_new(1).unwrap())]);
            let refinement = grid.refine(&plan);

            assert!(!refinement.was_refined(&IntervalId::new(0)));
            assert!(refinement.was_refined(&IntervalId::new(1)));
            assert!(!refinement.was_refined(&IntervalId::new(2)));
        }

        #[test]
        fn complex_refinement_plan() {
            let grid = Grid1D::uniform(
                IntervalClosed::new(real(0.0), real(4.0)),
                NumIntervals::try_new(4).unwrap(),
            );

            // Multiple intervals with different refinement levels
            let plan = BTreeMap::from([
                (IntervalId::new(0), PositiveNumPoints1D::try_new(2).unwrap()), // 3 sub-intervals
                (IntervalId::new(2), PositiveNumPoints1D::try_new(1).unwrap()), // 2 sub-intervals
                (IntervalId::new(3), PositiveNumPoints1D::try_new(3).unwrap()), // 4 sub-intervals
            ]);
            let refinement = grid.refine(&plan);

            // Total: 3 + 1 + 2 + 4 = 10 intervals
            assert_eq!(refinement.refined_grid().num_intervals().as_ref(), &10);

            // Verify individual refinements
            assert_eq!(
                refinement.get_refined_intervals(&IntervalId::new(0)).len(),
                3
            );
            assert_eq!(
                refinement.get_refined_intervals(&IntervalId::new(1)).len(),
                1
            );
            assert_eq!(
                refinement.get_refined_intervals(&IntervalId::new(2)).len(),
                2
            );
            assert_eq!(
                refinement.get_refined_intervals(&IntervalId::new(3)).len(),
                4
            );
        }

        #[test]
        fn empty_refinement_plan() {
            let grid = Grid1D::uniform(
                IntervalClosed::new(real(0.0), real(2.0)),
                NumIntervals::try_new(4).unwrap(),
            );

            let plan: BTreeMap<IntervalId, PositiveNumPoints1D> = BTreeMap::new();
            let refinement = grid.refine(&plan);

            // Grid should be unchanged
            assert_eq!(
                refinement.refined_grid().num_intervals(),
                refinement.original_grid().num_intervals()
            );
            assert_eq!(
                refinement.refined_grid().coords().deref(),
                refinement.original_grid().coords().deref()
            );
        }
    }

    mod data_transfer {
        use super::*;

        #[test]
        fn data_transfer_uniform_refinement() {
            let grid = Grid1D::uniform(
                IntervalClosed::new(real(0.0), real(2.0)),
                NumIntervals::try_new(2).unwrap(),
            );
            let original_data = vec![real(10.0), real(20.0)];

            let refinement = grid.refine_uniform(&PositiveNumPoints1D::try_new(1).unwrap());
            let transferred_data = refinement.transfer_interval_data(&original_data);

            assert_eq!(
                transferred_data,
                vec![real(10.0), real(10.0), real(20.0), real(20.0)]
            );
        }

        #[test]
        fn data_transfer_selective_refinement() {
            let grid = Grid1D::uniform(
                IntervalClosed::new(real(0.0), real(3.0)),
                NumIntervals::try_new(3).unwrap(),
            );
            let original_data = vec![real(1.0), real(2.0), real(3.0)];

            let plan =
                BTreeMap::from([(IntervalId::new(1), PositiveNumPoints1D::try_new(1).unwrap())]);
            let refinement = grid.refine(&plan);
            let transferred_data = refinement.transfer_interval_data(&original_data);

            // Data: 1 (no split), 2 (split into 2), 3 (no split) -> [1, 2, 2, 3]
            assert_eq!(
                transferred_data,
                vec![real(1.0), real(2.0), real(2.0), real(3.0)]
            );
        }

        #[test]
        fn data_transfer_preserves_count() {
            let grid = Grid1D::uniform(
                IntervalClosed::new(real(0.0), real(5.0)),
                NumIntervals::try_new(5).unwrap(),
            );
            let original_data: Vec<Real> = (1..=5).map(|i| real(i as f64)).collect();

            let plan = BTreeMap::from([
                (IntervalId::new(1), PositiveNumPoints1D::try_new(2).unwrap()),
                (IntervalId::new(3), PositiveNumPoints1D::try_new(1).unwrap()),
            ]);
            let refinement = grid.refine(&plan);
            let transferred_data = refinement.transfer_interval_data(&original_data);

            // Transferred data should have same length as refined intervals
            assert_eq!(
                transferred_data.len(),
                *refinement.refined_grid().num_intervals().as_ref()
            );
        }
    }

    mod iterators {
        use super::*;

        #[test]
        fn iter_refined_with_mapping() {
            let grid = Grid1D::uniform(
                IntervalClosed::new(real(0.0), real(2.0)),
                NumIntervals::try_new(2).unwrap(),
            );
            let refinement = grid.refine_uniform(&PositiveNumPoints1D::try_new(1).unwrap());

            let refined_mappings: Vec<(IntervalId, IntervalId)> =
                refinement.iter_refined_with_mapping().collect();

            assert_eq!(refined_mappings.len(), 4);
            assert_eq!(
                refined_mappings[0],
                (IntervalId::new(0), IntervalId::new(0))
            );
            assert_eq!(
                refined_mappings[1],
                (IntervalId::new(1), IntervalId::new(0))
            );
            assert_eq!(
                refined_mappings[2],
                (IntervalId::new(2), IntervalId::new(1))
            );
            assert_eq!(
                refined_mappings[3],
                (IntervalId::new(3), IntervalId::new(1))
            );
        }

        #[test]
        fn iter_original_with_refined() {
            let grid = Grid1D::uniform(
                IntervalClosed::new(real(0.0), real(2.0)),
                NumIntervals::try_new(2).unwrap(),
            );
            let refinement = grid.refine_uniform(&PositiveNumPoints1D::try_new(1).unwrap());

            let original_mappings: Vec<(IntervalId, &Vec<IntervalId>)> =
                refinement.iter_original_with_refined().collect();

            assert_eq!(original_mappings.len(), 2);
            assert_eq!(
                original_mappings[0].1,
                &vec![IntervalId::new(0), IntervalId::new(1)]
            );
            assert_eq!(
                original_mappings[1].1,
                &vec![IntervalId::new(2), IntervalId::new(3)]
            );
        }

        #[test]
        fn iter_intervals_bidirectional_mapping() {
            let grid = Grid1D::uniform(
                IntervalClosed::new(real(0.0), real(1.0)),
                NumIntervals::try_new(2).unwrap(),
            );
            let refinement = grid.refine_uniform(&PositiveNumPoints1D::try_new(1).unwrap());

            // Test that iter_refined_with_mapping returns consistent data
            let mappings: Vec<_> = refinement.iter_refined_with_mapping().collect();
            assert_eq!(mappings.len(), 4);

            // Each mapping should provide (refined_id, original_id)
            for (refined_id, original_id) in mappings {
                let expected_original = refinement.find_original_interval(&refined_id);
                assert_eq!(expected_original, original_id);
            }
        }
    }

    mod extraction_and_conversion {
        use super::*;

        #[test]
        fn into_refined_grid() {
            let grid = Grid1D::uniform(
                IntervalClosed::new(real(0.0), real(1.0)),
                NumIntervals::try_new(2).unwrap(),
            );
            let refinement = grid.refine_uniform(&PositiveNumPoints1D::try_new(1).unwrap());

            let refined_grid = refinement.into_refined_grid();
            assert_eq!(refined_grid.num_intervals().as_ref(), &4);
        }

        #[test]
        fn into_original_grid() {
            let grid = Grid1D::uniform(
                IntervalClosed::new(real(0.0), real(1.0)),
                NumIntervals::try_new(2).unwrap(),
            );
            let original_intervals = *grid.num_intervals().as_ref();

            let refinement = grid.refine_uniform(&PositiveNumPoints1D::try_new(1).unwrap());
            let original_grid = refinement.into_original_grid();

            assert_eq!(original_grid.num_intervals().as_ref(), &original_intervals);
        }

        #[test]
        fn into_tuple_decomposition() {
            let grid = Grid1D::uniform(
                IntervalClosed::new(real(0.0), real(1.0)),
                NumIntervals::try_new(2).unwrap(),
            );
            let refinement = grid.refine_uniform(&PositiveNumPoints1D::try_new(1).unwrap());

            let (refined_grid, original_grid, refined_to_original, original_to_refined) =
                refinement.into();

            assert_eq!(refined_grid.num_intervals().as_ref(), &4);
            assert_eq!(original_grid.num_intervals().as_ref(), &2);
            assert_eq!(refined_to_original.len(), 4);
            assert_eq!(original_to_refined.len(), 2);
        }

        #[test]
        fn convert_uniform_refinement_to_grid1d_uniform_refinement() {
            let uniform_grid = Grid1DUniform::new(
                IntervalClosed::new(real(0.0), real(1.0)),
                NumIntervals::try_new(2).unwrap(),
            );
            let uniform_refinement =
                uniform_grid.refine_uniform(&PositiveNumPoints1D::try_new(1).unwrap());

            let converted: Grid1DUniformRefinement<Grid1D<IntervalClosed<Real>>> =
                uniform_refinement.clone().into();

            assert!(matches!(converted.original_grid(), Grid1D::Uniform(_)));
            assert!(matches!(converted.refined_grid(), Grid1D::Uniform(_)));

            for i in 0..*converted.refined_grid().num_intervals().as_ref() {
                let id = IntervalId::new(i);
                assert_eq!(
                    converted.find_original_interval(&id),
                    uniform_refinement.find_original_interval(&id)
                );
            }
        }

        #[test]
        fn convert_non_uniform_uniform_refinement_to_grid1d_uniform_refinement() {
            let non_uniform_grid = Grid1DNonUniform::try_new_from_coords(
                Coords1D::try_from(SortedSet::from_unsorted(vec![
                    real(0.0),
                    real(0.2),
                    real(1.0),
                ]))
                .unwrap(),
            )
            .unwrap();
            let uniform_refinement =
                non_uniform_grid.refine_uniform(&PositiveNumPoints1D::try_new(1).unwrap());

            let converted: Grid1DUniformRefinement<Grid1D<IntervalClosed<Real>>> =
                uniform_refinement.clone().into();

            assert!(matches!(converted.original_grid(), Grid1D::NonUniform(_)));
            assert!(matches!(converted.refined_grid(), Grid1D::NonUniform(_)));

            for i in 0..*converted.refined_grid().num_intervals().as_ref() {
                let id = IntervalId::new(i);
                assert_eq!(
                    converted.find_original_interval(&id),
                    uniform_refinement.find_original_interval(&id)
                );
            }
        }
    }

    mod mathematical_properties {
        use super::*;

        #[test]
        fn bidirectional_mapping_consistency() {
            let grid = Grid1D::uniform(
                IntervalClosed::new(real(0.0), real(4.0)),
                NumIntervals::try_new(4).unwrap(),
            );

            let plan = BTreeMap::from([
                (IntervalId::new(0), PositiveNumPoints1D::try_new(2).unwrap()),
                (IntervalId::new(2), PositiveNumPoints1D::try_new(1).unwrap()),
            ]);
            let refinement = grid.refine(&plan);

            // For each original interval, check that all refined intervals map back correctly
            for original_idx in 0..*refinement.original_grid().num_intervals().as_ref() {
                let original_id = IntervalId::new(original_idx);
                let refined_ids = refinement.get_refined_intervals(&original_id);

                for refined_id in refined_ids {
                    assert_eq!(refinement.find_original_interval(refined_id), original_id);
                }
            }
        }

        #[test]
        fn total_intervals_consistency() {
            let grid = Grid1D::uniform(
                IntervalClosed::new(real(0.0), real(5.0)),
                NumIntervals::try_new(5).unwrap(),
            );

            let plan = BTreeMap::from([
                (IntervalId::new(1), PositiveNumPoints1D::try_new(2).unwrap()), // +2 intervals
                (IntervalId::new(3), PositiveNumPoints1D::try_new(3).unwrap()), // +3 intervals
            ]);
            let refinement = grid.refine(&plan);

            // Calculate expected: 5 original - 2 refined + (3 + 4) new = 10
            let expected_intervals = 5 + 2 + 3; // unrefined + extra from refinements
            assert_eq!(
                refinement.refined_grid().num_intervals().as_ref(),
                &expected_intervals
            );
        }

        #[test]
        fn sub_interval_lengths_are_equal() {
            let grid = Grid1D::uniform(
                IntervalClosed::new(real(0.0), real(2.0)),
                NumIntervals::try_new(2).unwrap(),
            );
            let refinement = grid.refine_uniform(&PositiveNumPoints1D::try_new(1).unwrap());

            // Each original interval of length 1.0 should become two intervals of length 0.5
            let expected_length = real(0.5);
            for i in 0..4 {
                let length = refinement
                    .refined_grid()
                    .interval_length(&IntervalId::new(i));
                // All sub-intervals should have equal length
                assert_eq!(length.as_ref(), &expected_length);
            }
        }

        #[test]
        fn domain_exactly_preserved() {
            let domain = IntervalClosed::new(real(-3.), real(2.));
            let grid = Grid1D::uniform(domain.clone(), NumIntervals::try_new(7).unwrap());
            let refinement = grid.refine_uniform(&PositiveNumPoints1D::try_new(2).unwrap());

            assert_eq!(refinement.refined_grid().domain(), &domain);
        }
    }

    mod edge_cases {
        use super::*;

        #[test]
        fn single_interval_grid() {
            let grid = Grid1D::uniform(
                IntervalClosed::new(real(0.0), real(1.0)),
                NumIntervals::try_new(1).unwrap(),
            );
            let refinement = grid.refine_uniform(&PositiveNumPoints1D::try_new(3).unwrap());

            assert_eq!(refinement.refined_grid().num_intervals().as_ref(), &4);
            assert_eq!(refinement.num_refined_intervals(&IntervalId::new(0)), 4);
        }

        #[test]
        fn large_refinement_factor() {
            let grid = Grid1D::uniform(
                IntervalClosed::new(real(0.0), real(1.0)),
                NumIntervals::try_new(2).unwrap(),
            );
            // Add 9 extra points -> 10 sub-intervals per original interval
            let refinement = grid.refine_uniform(&PositiveNumPoints1D::try_new(9).unwrap());

            assert_eq!(refinement.refined_grid().num_intervals().as_ref(), &20);
        }

        #[test]
        fn refine_first_and_last_intervals_only() {
            let grid = Grid1D::uniform(
                IntervalClosed::new(real(0.0), real(4.0)),
                NumIntervals::try_new(4).unwrap(),
            );

            let plan = BTreeMap::from([
                (IntervalId::new(0), PositiveNumPoints1D::try_new(1).unwrap()),
                (IntervalId::new(3), PositiveNumPoints1D::try_new(1).unwrap()),
            ]);
            let refinement = grid.refine(&plan);

            // 2 + 1 + 1 + 2 = 6 intervals
            assert_eq!(refinement.refined_grid().num_intervals().as_ref(), &6);

            // First coordinate should still be the original domain start
            assert_eq!(
                refinement.refined_grid().coords().first(),
                refinement.original_grid().coords().first()
            );
            // Last coordinate should still be the original domain end
            assert_eq!(
                refinement.refined_grid().coords().last(),
                refinement.original_grid().coords().last()
            );
        }
    }

    mod serialization {
        use super::*;

        #[test]
        fn uniform_refinement_serialization() {
            let grid = Grid1D::uniform(
                IntervalClosed::new(real(0.0), real(1.0)),
                NumIntervals::try_new(3).unwrap(),
            );
            let refinement = grid.refine_uniform(&PositiveNumPoints1D::try_new(1).unwrap());

            // Serialize to JSON
            let json = serde_json::to_string(&refinement).unwrap();

            // Deserialize from JSON
            let deserialized: Grid1DRefinement<
                Grid1D<IntervalClosed<Real>>,
                Grid1D<IntervalClosed<Real>>,
            > = serde_json::from_str(&json).unwrap();

            // Verify the deserialized refinement matches the original
            assert_eq!(
                refinement.refined_grid().coords().deref(),
                deserialized.refined_grid().coords().deref()
            );
            assert_eq!(
                refinement.original_grid().coords().deref(),
                deserialized.original_grid().coords().deref()
            );
            assert_eq!(
                refinement.refined_grid().num_intervals(),
                deserialized.refined_grid().num_intervals()
            );

            // Verify mappings are preserved
            for i in 0..*refinement.refined_grid().num_intervals().as_ref() {
                let id = IntervalId::new(i);
                assert_eq!(
                    refinement.find_original_interval(&id),
                    deserialized.find_original_interval(&id)
                );
            }
        }

        #[test]
        fn selective_refinement_serialization() {
            let grid = Grid1D::uniform(
                IntervalClosed::new(real(0.0), real(4.0)),
                NumIntervals::try_new(4).unwrap(),
            );

            let plan = BTreeMap::from([
                (IntervalId::new(0), PositiveNumPoints1D::try_new(2).unwrap()),
                (IntervalId::new(2), PositiveNumPoints1D::try_new(1).unwrap()),
            ]);
            let refinement = grid.refine(&plan);

            // Serialize to JSON
            let json = serde_json::to_string(&refinement).unwrap();

            // Deserialize from JSON using the correct type
            let deserialized: Grid1DNonUniformRefinement<Grid1D<IntervalClosed<Real>>> =
                serde_json::from_str(&json).unwrap();

            // Verify coordinates
            assert_eq!(
                refinement.refined_grid().coords().deref(),
                deserialized.refined_grid().coords().deref()
            );

            // Verify all mappings
            for i in 0..*refinement.refined_grid().num_intervals().as_ref() {
                let id = IntervalId::new(i);
                assert_eq!(
                    refinement.find_original_interval(&id),
                    deserialized.find_original_interval(&id)
                );
            }

            // Verify original interval to refined intervals mappings
            for i in 0..*refinement.original_grid().num_intervals().as_ref() {
                let id = IntervalId::new(i);
                assert_eq!(
                    refinement.num_refined_intervals(&id),
                    deserialized.num_refined_intervals(&id)
                );
            }
        }

        #[test]
        fn non_uniform_grid_refinement_serialization() {
            let coords = SortedSet::from_unsorted(vec![real(0.0), real(0.1), real(0.5), real(1.0)]);
            let grid = Grid1D::<IntervalClosed<Real>>::try_from_sorted(coords).unwrap();
            let refinement = grid.refine_uniform(&PositiveNumPoints1D::try_new(1).unwrap());

            // Serialize to JSON
            let json = serde_json::to_string(&refinement).unwrap();

            // Deserialize from JSON
            let deserialized: Grid1DRefinement<
                Grid1D<IntervalClosed<Real>>,
                Grid1D<IntervalClosed<Real>>,
            > = serde_json::from_str(&json).unwrap();

            // Verify coordinates match
            assert_eq!(
                refinement.refined_grid().coords().deref(),
                deserialized.refined_grid().coords().deref()
            );
            assert_eq!(
                refinement.original_grid().coords().deref(),
                deserialized.original_grid().coords().deref()
            );

            // Verify domain is preserved
            assert_eq!(
                refinement.refined_grid().domain(),
                deserialized.refined_grid().domain()
            );
        }

        #[test]
        fn round_trip_preserves_interval_data_transfer() {
            let grid = Grid1D::uniform(
                IntervalClosed::new(real(0.0), real(3.0)),
                NumIntervals::try_new(3).unwrap(),
            );
            let refinement = grid.refine_uniform(&PositiveNumPoints1D::try_new(1).unwrap());

            // Original interval data
            let coarse_data = vec![1.0, 2.0, 3.0];

            // Serialize and deserialize
            let json = serde_json::to_string(&refinement).unwrap();
            let deserialized: Grid1DRefinement<
                Grid1D<IntervalClosed<Real>>,
                Grid1D<IntervalClosed<Real>>,
            > = serde_json::from_str(&json).unwrap();

            // Transfer data using both original and deserialized refinements
            let fine_data_original = refinement.transfer_interval_data(&coarse_data);
            let fine_data_deserialized = deserialized.transfer_interval_data(&coarse_data);

            // Both should produce the same result
            assert_eq!(fine_data_original, fine_data_deserialized);
        }
    }
}