From 33be602b581697a3b34fb1f7c911cebd6686f8bb Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Thu, 15 Jan 2026 11:46:07 -0800
Subject: [PATCH 01/17] fix(vector_1D_t): mk weighted_premultiply private

---
 src/formal/tensors_1D_m.F90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/formal/tensors_1D_m.F90 b/src/formal/tensors_1D_m.F90
index 2615a7c..aa379cf 100644
--- a/src/formal/tensors_1D_m.F90
+++ b/src/formal/tensors_1D_m.F90
@@ -116,7 +116,7 @@ pure module function construct_1D_scalar_from_function(initializer, order, cells
     generic :: weights => gradient_1D_weights
 #endif
     procedure, non_overridable :: dA
-    procedure, non_overridable, pass(vector_1D) :: weighted_premultiply
+    procedure, non_overridable, private, pass(vector_1D) :: weighted_premultiply
     procedure, non_overridable, private :: div
     procedure, non_overridable, private :: dot_surface_normal
     procedure, non_overridable, private :: vector_1D_grid

From 112b667d7c3296f9b8eeb6d1bf17cb556bd60c24 Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Wed, 28 Jan 2026 12:35:47 -0800
Subject: [PATCH 02/17] feat(burgers-1D): define initial condition, order

---
 example/burgers-1D.F90 | 64 ++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 64 insertions(+)
 create mode 100644 example/burgers-1D.F90

diff --git a/example/burgers-1D.F90 b/example/burgers-1D.F90
new file mode 100644
index 0000000..b923d56
--- /dev/null
+++ b/example/burgers-1D.F90
@@ -0,0 +1,64 @@
+! Copyright (c) 2026, The Regents of the University of California
+! Terms of use are as specified in LICENSE.txt
+
+module initial_condition_m
+  implicit none
+
+contains
+
+  pure function initial_condition(x)
+    !! Initial solution to Burgers equation
+    double precision, intent(in) :: x(:)
+    double precision, allocatable :: initial_condition(:)
+    double precision, parameter :: pi = acos(-1D0)
+    initial_condition = sin(2*pi*x)
+    ! To change this function, please edit only the right-hand-side (RHS) expression,
+    ! keeping the rest in place for proper display of the function at runtime.
+  end function
+
+end module
+
+program burgers_1D
+  !! Advance the 1D Burgers partial differential equation over time.
+  use initial_condition_m, only : initial_condition
+  use julienne_m, only :  command_line_t
+  use formal_m, only : vector_1D_t, vector_1D_initializer_i
+  implicit none
+
+  procedure(vector_1D_initializer_i), pointer :: vector_1D_initializer => initial_condition
+  character(len=:), allocatable :: order_string
+  type(command_line_t) command_line
+  type(vector_1D_t) :: u
+  integer order
+
+  if (command_line%argument_present([character(len=len("--help")) :: ("--help"), "-h"])) then
+    stop                  new_line('') // new_line('') &
+      // 'Usage:'                      // new_line('') &
+      // '  fpm run \'                 // new_line('') &
+      // '  --example burgers-1D \'    // new_line('') &
+      // '  --compiler flang-new \'    // new_line('') &
+      // '  --flag "-O3" \'            // new_line('') & 
+      // '  -- [--help|-h] | [--order <integer>]' // new_line('') // new_line('') &
+      // 'where square brackets indicate optional arguments and angular brackets indicate user input values.' // new_line('')
+  end if
+
+  order_string = command_line%flag_value("--order")
+
+  print *,new_line('')
+  print *,"   Initial condition"
+  print *,"   ================="
+  call execute_command_line("grep 'initial_condition =' example/burgers-1D.F90 | grep -v execute_command", wait=.true.)
+
+  if (len(order_string)==0) then 
+    order = 4 
+  else
+    read(order_string,"(i)") order 
+  end if
+
+  print  *,"order = ",order
+
+#ifdef __GFORTRAN__
+    stop
+#endif
+
+end program

From 8ebe23e746c44f14403c8db35a4d0ff540cbba0e Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Thu, 29 Jan 2026 23:12:11 -0800
Subject: [PATCH 03/17] feat(dyad_1D_t): type, operator(/), operator(*)

This commit

1. Adds a dyad_1D_t type,,
2. Adds a vector_1D_t operator(/) with a RHS integer operand,
3. Adds operator(*) for vector_1D_t operands & a dyad_1D_t result,
4. Updates the burgers-1D example to compute u*u/2.
---
 example/burgers-1D.F90      | 12 +++++++++---
 src/formal/dyad_1D_s.F90    | 12 ++++++++++++
 src/formal/tensors_1D_m.F90 | 24 ++++++++++++++++++++++++
 src/formal/vector_1D_s.F90  |  8 ++++++++
 4 files changed, 53 insertions(+), 3 deletions(-)
 create mode 100644 src/formal/dyad_1D_s.F90

diff --git a/example/burgers-1D.F90 b/example/burgers-1D.F90
index b923d56..8914f59 100644
--- a/example/burgers-1D.F90
+++ b/example/burgers-1D.F90
@@ -42,20 +42,26 @@ program burgers_1D
       // 'where square brackets indicate optional arguments and angular brackets indicate user input values.' // new_line('')
   end if
 
-  order_string = command_line%flag_value("--order")
 
-  print *,new_line('')
+  print *, new_line('')
   print *,"   Initial condition"
   print *,"   ================="
+
   call execute_command_line("grep 'initial_condition =' example/burgers-1D.F90 | grep -v execute_command", wait=.true.)
 
+  order_string = command_line%flag_value("--order")
+
   if (len(order_string)==0) then 
     order = 4 
   else
     read(order_string,"(i)") order 
   end if
 
-  print  *,"order = ",order
+  print *, "order = ", order
+
+  u = vector_1D_t(vector_1D_initializer, order, x_min=0D0, x_max=1D0, cells=50)
+  associate(uu_2 => u*u/2)
+  end associate
 
 #ifdef __GFORTRAN__
     stop
diff --git a/src/formal/dyad_1D_s.F90 b/src/formal/dyad_1D_s.F90
new file mode 100644
index 0000000..73718bf
--- /dev/null
+++ b/src/formal/dyad_1D_s.F90
@@ -0,0 +1,12 @@
+! Copyright (c) 2026, The Regents of the University of California
+! Terms of use are as specified in LICENSE.txt
+
+submodule(tensors_1D_m) dyad_1D_s
+  implicit none
+contains
+
+  module procedure dyad_over_integer
+    ratio%tensor_1D_t = tensor_1D_t(self%tensor_1D_t%values_/numerator, self%x_min_, self%x_max_, self%cells_, order = self%order_)
+  end procedure
+
+end submodule dyad_1D_s
diff --git a/src/formal/tensors_1D_m.F90 b/src/formal/tensors_1D_m.F90
index aa379cf..b6e635d 100644
--- a/src/formal/tensors_1D_m.F90
+++ b/src/formal/tensors_1D_m.F90
@@ -16,6 +16,7 @@ module tensors_1D_m
 
   public :: scalar_1D_t
   public :: vector_1D_t
+  public :: dyad_1D_t 
   public :: gradient_1D_t
   public :: laplacian_1D_t
   public :: divergence_1D_t
@@ -110,6 +111,7 @@ pure module function construct_1D_scalar_from_function(initializer, order, cells
     generic :: operator(.x.)   => weighted_premultiply
     generic :: operator(.div.) => div
     generic :: operator(.dot.) => dot_surface_normal
+    generic :: operator(*)     => outer_product
     generic :: grid   => vector_1D_grid
     generic :: values => vector_1D_values
 #ifdef __INTEL_COMPILER
@@ -117,12 +119,19 @@ pure module function construct_1D_scalar_from_function(initializer, order, cells
 #endif
     procedure, non_overridable :: dA
     procedure, non_overridable, private, pass(vector_1D) :: weighted_premultiply
+    procedure, non_overridable, private :: outer_product
     procedure, non_overridable, private :: div
     procedure, non_overridable, private :: dot_surface_normal
     procedure, non_overridable, private :: vector_1D_grid
     procedure, non_overridable, private :: vector_1D_values
   end type
 
+  type, extends(tensor_1D_t) :: dyad_1D_t 
+  contains
+    generic :: operator(/) =>  dyad_over_integer
+    procedure, non_overridable, private :: dyad_over_integer
+  end type
+
   type, extends(tensor_1D_t) :: weighted_product_1D_t
   contains
     generic :: operator(.SS.) => surface_integrate_vector_x_scalar_1D
@@ -202,6 +211,14 @@ pure module function construct_from_components(tensor_1D, divergence_operator_1D
 
   interface
 
+    pure module function dyad_over_integer(self, numerator) result(ratio)
+      !! Result is the dyad divided by the numerator
+      implicit none
+      class(dyad_1D_t), intent(in) :: self
+      integer, intent(in) :: numerator
+      type(dyad_1D_t) ratio
+    end function
+
     pure module function dA(self)
       !! Result is the grid's discrete surface-area differential for use in surface integrals of the form
       !! .SS. (f .x. (v .dot. dA))
@@ -334,6 +351,13 @@ pure module function weighted_premultiply(scalar_1D, vector_1D) result(weighted_
       type(weighted_product_1D_t) weighted_product_1D
     end function
 
+    pure module function outer_product(lhs, rhs) result(dyad_1D)
+      !! Result is the dyadic product of two vector operands
+      implicit none
+      class(vector_1D_t), intent(in) :: lhs, rhs
+      type(dyad_1D_t) dyad_1D
+    end function
+
     pure module function gradient_1D_weights(self) result(weights)
       !! Result is an array of quadrature coefficients that can be used to compute a weighted
       !! inner product  of a vector_1D_t object and a gradient_1D_t object.
diff --git a/src/formal/vector_1D_s.F90 b/src/formal/vector_1D_s.F90
index 5bd892e..3c06737 100644
--- a/src/formal/vector_1D_s.F90
+++ b/src/formal/vector_1D_s.F90
@@ -22,6 +22,14 @@
 
 contains
 
+  module procedure outer_product
+
+     call_julienne_assert(.all. ([lhs%x_min_, lhs%x_max_] .approximates.   [rhs%x_min_, rhs%x_max_] .within. 0D0))
+     call_julienne_assert(.all. ([lhs%cells_, lhs%order_] .equalsExpected. [rhs%cells_, rhs%order_]))
+
+     dyad_1D%tensor_1D_t = tensor_1D_t(lhs%values_ * rhs%values_, lhs%x_min_, lhs%x_max_, lhs%cells_, lhs%order_ ) 
+  end procedure
+
   module procedure dot_surface_normal
      v_dot_dS%tensor_1D_t = tensor_1D_t(vector_1D%values_*dS, vector_1D%x_min_, vector_1D%x_max_, vector_1D%cells_, vector_1D%order_)
      v_dot_dS%divergence_operator_1D_ = vector_1D%divergence_operator_1D_

From fa84b55e9c503ab61eebf8a164c56f0521934cf1 Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Fri, 30 Jan 2026 00:09:45 -0800
Subject: [PATCH 04/17] feat(dyad_1D_t): add operator(.div.)

This commit adds support expressions of the form .div. uu, where
uu is a 1D dyadic tensor and the result is stored at cell centers.

To Do:
1. Embed an interpolation operator into the divergence operator
   to interpolate the result back out to the cell faces.
2. Use the interpolated values to construct the vector_1D_t result.
---
 example/burgers-1D.F90      |  2 +-
 src/formal/dyad_1D_s.F90    | 38 +++++++++++++++++++++++++++++++++++++
 src/formal/tensors_1D_m.F90 | 23 +++++++++++++++++++++-
 src/formal/vector_1D_s.F90  | 10 +++++++---
 4 files changed, 68 insertions(+), 5 deletions(-)

diff --git a/example/burgers-1D.F90 b/example/burgers-1D.F90
index 8914f59..1e2bde9 100644
--- a/example/burgers-1D.F90
+++ b/example/burgers-1D.F90
@@ -60,7 +60,7 @@ program burgers_1D
   print *, "order = ", order
 
   u = vector_1D_t(vector_1D_initializer, order, x_min=0D0, x_max=1D0, cells=50)
-  associate(uu_2 => u*u/2)
+  associate(div_uu_2 => .div. (u*u/2)) ! result is at cell centers; To Do: interpolate to cell faces
   end associate
 
 #ifdef __GFORTRAN__
diff --git a/src/formal/dyad_1D_s.F90 b/src/formal/dyad_1D_s.F90
index 73718bf..7687ab4 100644
--- a/src/formal/dyad_1D_s.F90
+++ b/src/formal/dyad_1D_s.F90
@@ -1,12 +1,50 @@
 ! Copyright (c) 2026, The Regents of the University of California
 ! Terms of use are as specified in LICENSE.txt
 
+#include "julienne-assert-macros.h"
+
 submodule(tensors_1D_m) dyad_1D_s
+  use julienne_m, only : call_julienne_assert_
   implicit none
 contains
 
   module procedure dyad_over_integer
     ratio%tensor_1D_t = tensor_1D_t(self%tensor_1D_t%values_/numerator, self%x_min_, self%x_max_, self%cells_, order = self%order_)
+    ratio%divergence_operator_1D_ = self%divergence_operator_1D_
+  end procedure
+
+  module procedure construct_1D_dyad_from_components
+    dyad_1D%tensor_1D_t = tensor_1D
+    dyad_1D%divergence_operator_1D_ = divergence_operator_1D
+  end procedure
+
+  module procedure div_dyad
+
+    integer center
+   
+#ifdef NAGFOR
+    associate(D => self%divergence_operator_1D_)
+#else
+    associate(D => (self%divergence_operator_1D_))
+#endif
+      call_julienne_assert(size(self%values_))
+      associate(Dv => D .x. self%values_)
+        divergence_1D%tensor_1D_t = tensor_1D_t(Dv(2:size(Dv)-1), self%x_min_, self%x_max_, self%cells_, self%order_)
+        error stop "To Do: interploate the dyad's divergence to the cell faces and construct a vector_1D_t"
+#if ASSERTIONS
+        associate( &
+           q  => divergence_1D%weights() &
+          ,dx => (self%x_max_ - self%x_min_)/self%cells_ &
+          ,b => [-1D0, [(0D0, center = 1, self%cells_-1)], 1D0] &
+        )
+          call_julienne_assert(.all. ([size(Dv), size(q)] .equalsExpected. self%cells_+2))
+          call_julienne_assert((.all. (matmul(transpose(D%assemble()), q) .approximates. b/dx .within. double_equivalence)))
+            ! Check D^T * a = b_{m+1},  Eq. (19), Corbino & Castillo (2020)
+        end associate
+#endif
+      end associate
+    end associate
+
   end procedure
 
 end submodule dyad_1D_s
diff --git a/src/formal/tensors_1D_m.F90 b/src/formal/tensors_1D_m.F90
index b6e635d..27cec5c 100644
--- a/src/formal/tensors_1D_m.F90
+++ b/src/formal/tensors_1D_m.F90
@@ -127,11 +127,25 @@ pure module function construct_1D_scalar_from_function(initializer, order, cells
   end type
 
   type, extends(tensor_1D_t) :: dyad_1D_t 
+    type(divergence_operator_1D_t) divergence_operator_1D_
   contains
-    generic :: operator(/) =>  dyad_over_integer
+    generic :: operator(.div.) =>  div_dyad
+    generic :: operator(/)     =>  dyad_over_integer
     procedure, non_overridable, private :: dyad_over_integer
+    procedure, non_overridable, private :: div_dyad
   end type
 
+  interface dyad_1D_t
+
+    pure module function construct_1D_dyad_from_components(tensor_1D, divergence_operator_1D) result(dyad_1D)
+      !! Result is a 1D dyadic tensor with the provided parent component tensor_1D and the provided divergence operatror
+      type(tensor_1D_t), intent(in) :: tensor_1D
+      type(divergence_operator_1D_t), intent(in) :: divergence_operator_1D
+      type(dyad_1D_t) dyad_1D
+    end function
+
+  end interface
+
   type, extends(tensor_1D_t) :: weighted_product_1D_t
   contains
     generic :: operator(.SS.) => surface_integrate_vector_x_scalar_1D
@@ -297,6 +311,13 @@ pure module function reduced_order_boundary_depth(self) result(num_nodes)
       integer num_nodes
     end function
 
+    pure module function div_dyad(self) result(divergence_1D)
+      !! Result is mimetic divergence of the dyad_1D_t "self"
+      implicit none
+      class(dyad_1D_t), intent(in) :: self
+      type(divergence_1D_t) divergence_1D !! discrete divergence
+    end function
+
     pure module function div(self) result(divergence_1D)
       !! Result is mimetic divergence of the vector_1D_t "self"
       implicit none
diff --git a/src/formal/vector_1D_s.F90 b/src/formal/vector_1D_s.F90
index 3c06737..f4f7159 100644
--- a/src/formal/vector_1D_s.F90
+++ b/src/formal/vector_1D_s.F90
@@ -24,10 +24,12 @@
 
   module procedure outer_product
 
-     call_julienne_assert(.all. ([lhs%x_min_, lhs%x_max_] .approximates.   [rhs%x_min_, rhs%x_max_] .within. 0D0))
-     call_julienne_assert(.all. ([lhs%cells_, lhs%order_] .equalsExpected. [rhs%cells_, rhs%order_]))
+    call_julienne_assert(.all. ([lhs%x_min_, lhs%x_max_] .approximates.   [rhs%x_min_, rhs%x_max_] .within. 0D0))
+    call_julienne_assert(.all. ([lhs%cells_, lhs%order_] .equalsExpected. [rhs%cells_, rhs%order_]))
+
+    dyad_1D%tensor_1D_t = tensor_1D_t(lhs%values_ * rhs%values_, lhs%x_min_, lhs%x_max_, lhs%cells_, lhs%order_ ) 
+    dyad_1D%divergence_operator_1D_ = divergence_operator_1D_t(k=lhs%order_, dx=(lhs%x_max_ - lhs%x_min_)/lhs%cells_, cells=lhs%cells_)
 
-     dyad_1D%tensor_1D_t = tensor_1D_t(lhs%values_ * rhs%values_, lhs%x_min_, lhs%x_max_, lhs%cells_, lhs%order_ ) 
   end procedure
 
   module procedure dot_surface_normal
@@ -84,6 +86,7 @@ pure module function construct_1D_vector_from_function(initializer, order, cells
 #endif
       associate(Dv => D .x. self%values_)
         divergence_1D%tensor_1D_t = tensor_1D_t(Dv(2:size(Dv)-1), self%x_min_, self%x_max_, self%cells_, self%order_)
+#if ASSERTIONS
         associate( &
            q  => divergence_1D%weights() &
           ,dx => (self%x_max_ - self%x_min_)/self%cells_ &
@@ -93,6 +96,7 @@ pure module function construct_1D_vector_from_function(initializer, order, cells
           call_julienne_assert((.all. (matmul(transpose(D%assemble()), q) .approximates. b/dx .within. double_equivalence)))
             ! Check D^T * a = b_{m+1},  Eq. (19), Corbino & Castillo (2020)
         end associate
+#endif
       end associate
     end associate
 

From 3cb46429a6819b2abcbe170020f16c6e38ea8ebb Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Thu, 19 Feb 2026 21:55:43 -0800
Subject: [PATCH 05/17] refac: rename mimetic_* -> differential_operator_*

This commit prepares for disambiguation of two types of mimetic
matrix abstractions: one for differential operators and one for
interpolation operators.
---
 ... => differential_operator_matrix_1D_s.F90} | 10 +++++-----
 ..._m.F90 => differential_operators_1D_m.F90} | 20 +++++++++----------
 src/formal/divergence_operator_1D_s.F90       |  4 ++--
 src/formal/gradient_operator_1D_s.F90         |  4 ++--
 src/formal/tensors_1D_m.F90                   |  2 +-
 src/formal_m.f90                              |  2 +-
 6 files changed, 21 insertions(+), 21 deletions(-)
 rename src/formal/{mimetic_matrix_1D_s.F90 => differential_operator_matrix_1D_s.F90} (76%)
 rename src/formal/{mimetic_operators_1D_m.F90 => differential_operators_1D_m.F90} (91%)

diff --git a/src/formal/mimetic_matrix_1D_s.F90 b/src/formal/differential_operator_matrix_1D_s.F90
similarity index 76%
rename from src/formal/mimetic_matrix_1D_s.F90
rename to src/formal/differential_operator_matrix_1D_s.F90
index fd0b114..b40d573 100644
--- a/src/formal/mimetic_matrix_1D_s.F90
+++ b/src/formal/differential_operator_matrix_1D_s.F90
@@ -1,16 +1,16 @@
 ! Copyright (c) 2026, The Regents of the University of California
 ! Terms of use are as specified in LICENSE.txt
 
-submodule(mimetic_operators_1D_m) mimetic_matrix_1D_s
+submodule(differential_operators_1D_m) differential_operator_matrix_1D_s
   use julienne_m, only : string_t, operator(.csv.)
   implicit none
 
 contains
 
   module procedure construct_matrix_operator
-    mimetic_matrix_1D%upper_ = upper
-    mimetic_matrix_1D%inner_ = inner
-    mimetic_matrix_1D%lower_ = lower
+    differential_operator_matrix_1D%upper_ = upper
+    differential_operator_matrix_1D%inner_ = inner
+    differential_operator_matrix_1D%lower_ = lower
   end procedure
 
   module procedure to_file_t
@@ -33,4 +33,4 @@
 
   end procedure
 
-end submodule mimetic_matrix_1D_s
\ No newline at end of file
+end submodule differential_operator_matrix_1D_s
\ No newline at end of file
diff --git a/src/formal/mimetic_operators_1D_m.F90 b/src/formal/differential_operators_1D_m.F90
similarity index 91%
rename from src/formal/mimetic_operators_1D_m.F90
rename to src/formal/differential_operators_1D_m.F90
index 22a4670..41f7bc2 100644
--- a/src/formal/mimetic_operators_1D_m.F90
+++ b/src/formal/differential_operators_1D_m.F90
@@ -3,7 +3,7 @@
 
 #include "formal-language-support.F90"
 
-module mimetic_operators_1D_m
+module differential_operators_1D_m
   !! Define sparse matrix storage formats and operators tailored to the one-dimensional (1D) mimetic discretizations
   !! detaild by Corbino & Castillo (2020) https://doi.org/10.1016/j.cam.2019.06.042.
   use julienne_m, only : file_t
@@ -13,8 +13,8 @@ module mimetic_operators_1D_m
   public :: gradient_operator_1D_t
   public :: divergence_operator_1D_t
 
-  type mimetic_matrix_1D_t
-    !! Encapsulate a mimetic matrix with a corresponding matrix-vector product operator
+  type differential_operator_matrix_1D_t
+    !! Encapsulate a mimetic matrix
     private
     double precision, allocatable :: upper_(:,:) !! A  submatrix block (cf. Corbino & Castillo, 2020)
     double precision, allocatable :: inner_(:)   !! M  submatrix row   (cf. Corbino & Castillo, 2020)
@@ -23,20 +23,20 @@ module mimetic_operators_1D_m
     procedure, non_overridable :: to_file_t
   end type
 
-  interface mimetic_matrix_1D_t
+  interface differential_operator_matrix_1D_t
 
-    pure module function construct_matrix_operator(upper, inner, lower) result(mimetic_matrix_1D)
+    pure module function construct_matrix_operator(upper, inner, lower) result(differential_operator_matrix_1D)
       !! Construct discrete operator from matrix blocks
       implicit none
       double precision, intent(in) :: upper(:,:) !! A  submatrix block (cf. Corbino & Castillo, 2020)
       double precision, intent(in) :: inner(:)   !! M  submatrix row   (cf. Corbino & Castillo, 2020)
       double precision, intent(in) :: lower(:,:) !! A' submatrix block (cf. Corbino & Castillo, 2020)
-      type(mimetic_matrix_1D_t) mimetic_matrix_1D
+      type(differential_operator_matrix_1D_t) differential_operator_matrix_1D
     end function
 
   end interface
 
-  type, extends(mimetic_matrix_1D_t) :: gradient_operator_1D_t
+  type, extends(differential_operator_matrix_1D_t) :: gradient_operator_1D_t
     !! Encapsulate a 1D mimetic gradient operator
     private
     integer k_ !! order of accuracy
@@ -62,7 +62,7 @@ pure module function construct_1D_gradient_operator(k, dx, cells) result(gradien
 
   end interface
 
-  type, extends(mimetic_matrix_1D_t) :: divergence_operator_1D_t
+  type, extends(differential_operator_matrix_1D_t) :: divergence_operator_1D_t
     !! Encapsulate kth-order mimetic divergence operator on m_ cells of width dx
     private
     integer k_, m_
@@ -129,7 +129,7 @@ pure module function divergence_matrix_multiply(self, vec) result(matvec_product
 
      pure module function to_file_t(self) result(file)
        implicit none
-       class(mimetic_matrix_1D_t), intent(in) :: self
+       class(differential_operator_matrix_1D_t), intent(in) :: self
        type(file_t) file
      end function
 
@@ -164,4 +164,4 @@ pure function negate_and_flip(A) result(Ap)
 
 #endif
 
-end module mimetic_operators_1D_m
+end module differential_operators_1D_m
diff --git a/src/formal/divergence_operator_1D_s.F90 b/src/formal/divergence_operator_1D_s.F90
index 3da38d0..68c23d4 100644
--- a/src/formal/divergence_operator_1D_s.F90
+++ b/src/formal/divergence_operator_1D_s.F90
@@ -4,7 +4,7 @@
 #include "formal-language-support.F90"
 #include "julienne-assert-macros.h"
 
-submodule(mimetic_operators_1D_m) divergence_operator_1D_s
+submodule(differential_operators_1D_m) divergence_operator_1D_s
   use julienne_m, only : call_julienne_assert_, string_t
 #if ASSERTIONS
   use julienne_m, only : operator(.isAtLeast.), operator(.equalsExpected.)
@@ -48,7 +48,7 @@ pure function negate_and_flip(A) result(Ap)
       else
         allocate(Ap, mold = A)
       end if
-      divergence_operator_1D%mimetic_matrix_1D_t = mimetic_matrix_1D_t(A, M(k, dx), Ap)
+      divergence_operator_1D%differential_operator_matrix_1D_t = differential_operator_matrix_1D_t(A, M(k, dx), Ap)
       divergence_operator_1D%k_  = k
       divergence_operator_1D%dx_ = dx
       divergence_operator_1D%m_  = cells
diff --git a/src/formal/gradient_operator_1D_s.F90 b/src/formal/gradient_operator_1D_s.F90
index 9f6ebed..e758b01 100644
--- a/src/formal/gradient_operator_1D_s.F90
+++ b/src/formal/gradient_operator_1D_s.F90
@@ -4,7 +4,7 @@
 #include "julienne-assert-macros.h"
 #include "formal-language-support.F90"
 
-submodule(mimetic_operators_1D_m) gradient_operator_1D_s
+submodule(differential_operators_1D_m) gradient_operator_1D_s
   use julienne_m, only : call_julienne_assert_, string_t
 #if ASSERTIONS
   use julienne_m, only : operator(.isAtLeast.)
@@ -42,7 +42,7 @@ pure function negate_and_flip(A) result(Ap)
     call_julienne_assert(cells .isAtLeast. 2*k)
 
     associate(A => corbino_castillo_A(k, dx), M => corbino_castillo_M(k, dx))
-      gradient_operator_1D%mimetic_matrix_1D_t = mimetic_matrix_1D_t(A, M, negate_and_flip(A))
+      gradient_operator_1D%differential_operator_matrix_1D_t = differential_operator_matrix_1D_t(A, M, negate_and_flip(A))
       gradient_operator_1D%k_  = k
       gradient_operator_1D%dx_ = dx
       gradient_operator_1D%m_  = cells
diff --git a/src/formal/tensors_1D_m.F90 b/src/formal/tensors_1D_m.F90
index 27cec5c..0a3527f 100644
--- a/src/formal/tensors_1D_m.F90
+++ b/src/formal/tensors_1D_m.F90
@@ -8,7 +8,7 @@ module tensors_1D_m
   !! divergence, and Laplacian operators as detailed by Corbino & Castillo (2020)
   !! https://doi.org/10.1016/j.cam.2019.06.042.
   use julienne_m, only : file_t
-  use mimetic_operators_1D_m, only : divergence_operator_1D_t, gradient_operator_1D_t
+  use differential_operators_1D_m, only : divergence_operator_1D_t, gradient_operator_1D_t
     
   implicit none
 
diff --git a/src/formal_m.f90 b/src/formal_m.f90
index 523cfa2..8128923 100644
--- a/src/formal_m.f90
+++ b/src/formal_m.f90
@@ -16,7 +16,7 @@ module formal_m
     ,scalar_1D_initializer_i  & ! scalar_1D_t initializer abstract interface
     ,vector_1D_initializer_i    ! vector_1D_t initializar abstract interface
 
-  use mimetic_operators_1D_m, only : &
+  use differential_operators_1D_m, only : &
      gradient_operator_1D_t & ! matrix operator defining a 1D mimetic gradient
     ,divergence_operator_1D_t ! matrix operator defining a 1D mimetic divergence
 

From 936d70b2f22568586073734b860af7d70d7f0c67 Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Thu, 19 Feb 2026 23:27:26 -0800
Subject: [PATCH 06/17] feat(interpolation): add operators

This commit adds interpolation operator matrix abstractions to
go from cell centers to faces and vice versa based on the
definitions in Dumett & Castillo (2022).

https://doi.org/10.13140/RG.2.2.26630.14400
---
 src/formal/interpolation_operators_1D_m.F90 | 54 ++++++++++++++++++++
 src/formal/interpolation_operators_1D_s.F90 | 55 +++++++++++++++++++++
 2 files changed, 109 insertions(+)
 create mode 100644 src/formal/interpolation_operators_1D_m.F90
 create mode 100644 src/formal/interpolation_operators_1D_s.F90

diff --git a/src/formal/interpolation_operators_1D_m.F90 b/src/formal/interpolation_operators_1D_m.F90
new file mode 100644
index 0000000..df8b156
--- /dev/null
+++ b/src/formal/interpolation_operators_1D_m.F90
@@ -0,0 +1,54 @@
+! Copyright (c) 2026, The Regents of the University of California
+! Terms of use are as specified in LICENSE.txt
+
+module interpolation_operators_1D_m
+  !! Define a sparse matrix storage format tailored to the staggered-grid interpolation matrix
+  !! operators detailed by Dumett & Castillo (2022) https://doi.org/10.13140/RG.2.2.26630.14400
+  implicit none
+
+  private
+  public :: cells_to_faces_t
+  public :: faces_to_cells_t
+
+  type interpolation_operator_1D_t
+    !! Encapsulate a staggered-grid interpolation matrix with a corresponding matrix-vector product operator
+    private
+    integer order
+    double precision                 first_
+    double precision, allocatable :: upper_(:,:) 
+    double precision, allocatable :: inner_(:)
+    double precision, allocatable :: lower_(:,:)
+    double precision                 final_
+  end type
+
+  type, extends(interpolation_operator_1D_t) :: centers_to_faces_1D_t
+  end type
+
+  type, extends(interpolation_operator_1D_t) :: faces_to_centers_1D_t
+  end type
+
+  interface centers_to_faces_1D_t
+
+    pure module function c2f_component_constructor(order, cells, dx) result(centers_to_faces_1D)
+      !! Construct centers-to-faces interpolation operator
+      implicit none
+      integer, intent(in) :: order, cells
+      double precision, intent(in) :: dx
+      type(centers_to_faces_1D_t) centers_to_faces_1D
+    end function
+
+  end interface
+
+  interface faces_to_centers_1D_t
+
+    pure module function f2c_component_constructor(order, cells, dx) result(faces_to_centers_1D)
+      !! Construct faces-to-centers interpolation operator
+      implicit none
+      integer, intent(in) :: order, cells
+      double precision, intent(in) :: dx
+      type(faces_to_centers_1D_t) faces_to_centers_1D
+    end function
+
+  end interface
+
+end module interpolation_operators_1D_m
diff --git a/src/formal/interpolation_operators_1D_s.F90 b/src/formal/interpolation_operators_1D_s.F90
new file mode 100644
index 0000000..05549ab
--- /dev/null
+++ b/src/formal/interpolation_operators_1D_s.F90
@@ -0,0 +1,55 @@
+! Copyright (c) 2026, The Regents of the University of California
+! Terms of use are as specified in LICENSE.txt
+
+submodule(interpolation_operators_1D_m) interpolation_operators_1D_s
+  implicit none
+
+contains
+
+  module procedure c2f_component_constructor
+
+    centers_to_faces_1D%cells_ = cells
+
+    select case(order)
+    case(2)
+      centers_to_faces_1D%first_ = 1D0
+      centers_to_faces_1D%upper_ = reshape([double precision::], [0,3])
+      centers_to_faces_1D%inner_ = [1D0,1D0]/2D0
+      centers_to_faces_1D%lower_ = reshape([double precision::], [0,3])
+      centers_to_faces_1D%final_ = 1D0
+    case(4)
+      centers_to_faces_1D%first_ = 1D0
+      centers_to_faces_1D%upper_ = reshape([-16,  70, 70, -14,   2], [1,5])/112D0
+      centers_to_faces_1D%inner_ =              [ -7,  63, 63,  -7]        /112D0
+      centers_to_faces_1D%lower_ = reshape([  2, -14, 70,  70, -16], [1,5])/112D0
+      centers_to_faces_1D%final_ = 1D0
+    case default
+      error stop "c2f_component_constructor: unsupported order"
+    end select
+
+  end procedure
+
+  module procedure f2c_component_constructor
+
+    faces_to_centers_1D%cells_ = cells
+
+    select case(order)
+    case(2)
+      faces_to_centers_1D%first_ = 1D0
+      faces_to_centers_1D%upper_ = reshape([double precision::], [0,3])
+      faces_to_centers_1D%inner_ = [1D0,1D0]/2D0
+      faces_to_centers_1D%lower_ = reshape([double precision::], [0,3])
+      faces_to_centers_1D%final_ = 1D0
+    case(4)
+      faces_to_centers_1D%first_ = 1D0
+      faces_to_centers_1D%upper_ = reshape([35, 140, -70, 28, -5], [1,5])/128D0
+      faces_to_centers_1D%inner_ =         [-8,  72,  72, -8]            /128D0
+      faces_to_centers_1D%lower_ = reshape([-5, 28, -70, 140, 35], [1,5])/128D0
+      faces_to_centers_1D%final_ = 1D0
+    case default
+      error stop "f2c_component_constructor: unsupported order"
+    end select
+
+  end procedure
+
+end submodule interpolation_operators_1D_s

From 1fe4b878362b85f7030d289ab298a9872e209457 Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Mon, 23 Feb 2026 14:47:50 -0800
Subject: [PATCH 07/17] WIP: start defining divergence of dyad

---
 src/formal/tensors_1D_m.F90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/formal/tensors_1D_m.F90 b/src/formal/tensors_1D_m.F90
index 0a3527f..0249276 100644
--- a/src/formal/tensors_1D_m.F90
+++ b/src/formal/tensors_1D_m.F90
@@ -315,7 +315,7 @@ pure module function div_dyad(self) result(divergence_1D)
       !! Result is mimetic divergence of the dyad_1D_t "self"
       implicit none
       class(dyad_1D_t), intent(in) :: self
-      type(divergence_1D_t) divergence_1D !! discrete divergence
+      type(vector_1D_t) vector_1D
     end function
 
     pure module function div(self) result(divergence_1D)

From fc3a408d4e7edcea2eec5c0c28a4932e517daabe Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Fri, 13 Mar 2026 17:05:43 -0700
Subject: [PATCH 08/17] doc(README): fix test command for fpm version<0.13

---
 README.md | 16 ++++++++++++----
 1 file changed, 12 insertions(+), 4 deletions(-)

diff --git a/README.md b/README.md
index f1721b6..938adbf 100644
--- a/README.md
+++ b/README.md
@@ -31,7 +31,7 @@ $$ \iiint_V (\vec{v} \cdot \nabla f) dV + \iiint_V (f \nabla \cdot \vec{v}) dV =
 
 Running the program as follows
 ```fortran
-fpm run --example extended-gauss-divergence --compiler flang-new --flag -O3
+fpm run --example extended-gauss-divergence --compiler flang --profile release
 ```
 produces output that includes actual program syntax:
 ```fortran
@@ -77,11 +77,19 @@ Building and testing
  GCC      | `gfortran`  | 13                    | `fpm test --compiler gfortran --profile release --flag "-ffree-line-length-none"`
  Intel    | `ifx`       | 2025.1.2              | `FOR_COARRAY_NUM_IMAGES=1 fpm test --compiler ifx --flag "-fpp -O3 -coarray" --profile release`
  LFortran | `lfortran`  | 0.60.0-421-ge2c448c79 | `fpm test --compiler lfortran --flag "--cpp --realloc-lhs-arrays"`
- LLVM     | `flang`     | 20-21                 | `fpm test --compiler flang --flag "-O3"`
- LLVM     | `flang`     | 19                    | `fpm test --compiler flang --flag "-O3 -mmlir -allow-assumed-rank"`
+ LLVM     | `flang`     | 20-21                 | `fpm test --compiler flang --profile release
+ LLVM     | `flang`     | 19                    | `fpm test --compiler flang --profile release --flag "-mmlir -allow-assumed-rank"`
  NAG      | `nagfor`    | 7.2 Build 7242        | `fpm test --compiler nagfor --flag "-O3 -fpp"`
 
-With `fpm` versions before 0.13.0, replace `flang` with `flang-new` above.
+### `fpm` versions before 0.13.0
+With LLVM 20-22, replace the above `flang` command with
+```
+  fpm test --compiler flang --flag "-O3"
+```
+With LLVM 19, replace the above `flang` command with
+```
+ fpm test --compiler flang --flag "-O3 -mmlir -allow-assumed-rank"
+```
 
 Documentation
 -------------

From 6097fc0fa4c170fc9abdf6e5050850dee05deb07 Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Fri, 13 Mar 2026 17:08:44 -0700
Subject: [PATCH 09/17] doc(README): flang-new -> flang with old fpm

---
 README.md | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/README.md b/README.md
index 938adbf..977e80b 100644
--- a/README.md
+++ b/README.md
@@ -84,11 +84,11 @@ Building and testing
 ### `fpm` versions before 0.13.0
 With LLVM 20-22, replace the above `flang` command with
 ```
-  fpm test --compiler flang --flag "-O3"
+  fpm test --compiler flang-new --flag "-O3"
 ```
 With LLVM 19, replace the above `flang` command with
 ```
- fpm test --compiler flang --flag "-O3 -mmlir -allow-assumed-rank"
+ fpm test --compiler flang-new --flag "-O3 -mmlir -allow-assumed-rank"
 ```
 
 Documentation

From e3ab6918b7bb474acf1c9b6c3a0484f15337a2ab Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Sat, 14 Mar 2026 10:12:38 -0700
Subject: [PATCH 10/17] doc(README): merge upstream changes

---
 README.md | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/README.md b/README.md
index 977e80b..fd10780 100644
--- a/README.md
+++ b/README.md
@@ -77,16 +77,16 @@ Building and testing
  GCC      | `gfortran`  | 13                    | `fpm test --compiler gfortran --profile release --flag "-ffree-line-length-none"`
  Intel    | `ifx`       | 2025.1.2              | `FOR_COARRAY_NUM_IMAGES=1 fpm test --compiler ifx --flag "-fpp -O3 -coarray" --profile release`
  LFortran | `lfortran`  | 0.60.0-421-ge2c448c79 | `fpm test --compiler lfortran --flag "--cpp --realloc-lhs-arrays"`
- LLVM     | `flang`     | 20-21                 | `fpm test --compiler flang --profile release
+ LLVM     | `flang`     | 20-22                 | `fpm test --compiler flang --profile release`
  LLVM     | `flang`     | 19                    | `fpm test --compiler flang --profile release --flag "-mmlir -allow-assumed-rank"`
  NAG      | `nagfor`    | 7.2 Build 7242        | `fpm test --compiler nagfor --flag "-O3 -fpp"`
 
 ### `fpm` versions before 0.13.0
-With LLVM 20-22, replace the above `flang` command with
+With LLVM 20-22, replace the corresponding command above with
 ```
   fpm test --compiler flang-new --flag "-O3"
 ```
-With LLVM 19, replace the above `flang` command with
+With LLVM 19, replace the corresponding command above with
 ```
  fpm test --compiler flang-new --flag "-O3 -mmlir -allow-assumed-rank"
 ```

From a5ea34c81b417a6816195eff7c09690e8ebc0097 Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Sat, 14 Mar 2026 10:14:31 -0700
Subject: [PATCH 11/17] doc(README): resolve git-stash conflict

---
 README.md | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/README.md b/README.md
index fd10780..977e80b 100644
--- a/README.md
+++ b/README.md
@@ -77,16 +77,16 @@ Building and testing
  GCC      | `gfortran`  | 13                    | `fpm test --compiler gfortran --profile release --flag "-ffree-line-length-none"`
  Intel    | `ifx`       | 2025.1.2              | `FOR_COARRAY_NUM_IMAGES=1 fpm test --compiler ifx --flag "-fpp -O3 -coarray" --profile release`
  LFortran | `lfortran`  | 0.60.0-421-ge2c448c79 | `fpm test --compiler lfortran --flag "--cpp --realloc-lhs-arrays"`
- LLVM     | `flang`     | 20-22                 | `fpm test --compiler flang --profile release`
+ LLVM     | `flang`     | 20-21                 | `fpm test --compiler flang --profile release
  LLVM     | `flang`     | 19                    | `fpm test --compiler flang --profile release --flag "-mmlir -allow-assumed-rank"`
  NAG      | `nagfor`    | 7.2 Build 7242        | `fpm test --compiler nagfor --flag "-O3 -fpp"`
 
 ### `fpm` versions before 0.13.0
-With LLVM 20-22, replace the corresponding command above with
+With LLVM 20-22, replace the above `flang` command with
 ```
   fpm test --compiler flang-new --flag "-O3"
 ```
-With LLVM 19, replace the corresponding command above with
+With LLVM 19, replace the above `flang` command with
 ```
  fpm test --compiler flang-new --flag "-O3 -mmlir -allow-assumed-rank"
 ```

From f9d08bb2d48d625bc024fe0cf3b64b5bacc45797 Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Mon, 16 Mar 2026 01:47:12 -0700
Subject: [PATCH 12/17] feat(interpolate): .div.dyad maps centers to faces

This commit adds a centers-to-faces interpolator and uses that
interpolator when computing the divergence of a dyad.
---
 example/burgers-1D.F90                        | 19 ++++++----
 src/formal/dyad_1D_s.F90                      | 31 +++++++++------
 ...erators_1D_m.F90 => interpolator_1D_m.F90} | 38 +++++++++++++------
 ...erators_1D_s.F90 => interpolator_1D_s.F90} | 33 ++++++++++++++--
 src/formal/tensors_1D_m.F90                   |  4 +-
 5 files changed, 89 insertions(+), 36 deletions(-)
 rename src/formal/{interpolation_operators_1D_m.F90 => interpolator_1D_m.F90} (52%)
 rename src/formal/{interpolation_operators_1D_s.F90 => interpolator_1D_s.F90} (57%)

diff --git a/example/burgers-1D.F90 b/example/burgers-1D.F90
index 1e2bde9..d7d7085 100644
--- a/example/burgers-1D.F90
+++ b/example/burgers-1D.F90
@@ -35,14 +35,13 @@ program burgers_1D
     stop                  new_line('') // new_line('') &
       // 'Usage:'                      // new_line('') &
       // '  fpm run \'                 // new_line('') &
-      // '  --example burgers-1D \'    // new_line('') &
-      // '  --compiler flang-new \'    // new_line('') &
-      // '  --flag "-O3" \'            // new_line('') & 
-      // '  -- [--help|-h] | [--order <integer>]' // new_line('') // new_line('') &
+      // '    --example burgers-1D \'  // new_line('') &
+      // '    --compiler flang \'      // new_line('') &
+      // '    --flag "-O3" \'          // new_line('') & 
+      // '    [--help|-h] | [--order <integer>]' // new_line('') // new_line('') &
       // 'where square brackets indicate optional arguments and angular brackets indicate user input values.' // new_line('')
   end if
 
-
   print *, new_line('')
   print *,"   Initial condition"
   print *,"   ================="
@@ -60,8 +59,14 @@ program burgers_1D
   print *, "order = ", order
 
   u = vector_1D_t(vector_1D_initializer, order, x_min=0D0, x_max=1D0, cells=50)
-  associate(div_uu_2 => .div. (u*u/2)) ! result is at cell centers; To Do: interpolate to cell faces
-  end associate
+
+  block
+    double precision dt
+    dt = 1D0
+    !associate(dt_grad_u2_2 => dt * (.div. ((u*u))))
+    associate(dt_div_uu => .div. (u*u))
+    end associate
+  end block
 
 #ifdef __GFORTRAN__
     stop
diff --git a/src/formal/dyad_1D_s.F90 b/src/formal/dyad_1D_s.F90
index 7687ab4..3b685e4 100644
--- a/src/formal/dyad_1D_s.F90
+++ b/src/formal/dyad_1D_s.F90
@@ -5,6 +5,7 @@
 
 submodule(tensors_1D_m) dyad_1D_s
   use julienne_m, only : call_julienne_assert_
+  use interpolator_1D_m, only : centers_to_faces_1D_t
   implicit none
 contains
 
@@ -28,18 +29,26 @@
     associate(D => (self%divergence_operator_1D_))
 #endif
       call_julienne_assert(size(self%values_))
-      associate(Dv => D .x. self%values_)
-        divergence_1D%tensor_1D_t = tensor_1D_t(Dv(2:size(Dv)-1), self%x_min_, self%x_max_, self%cells_, self%order_)
-        error stop "To Do: interploate the dyad's divergence to the cell faces and construct a vector_1D_t"
+      associate( &
+         Dv => D .x. self%values_ &
+        ,dx => (self%x_max_ - self%x_min_)/self%cells_ &
+      )
+        associate(interpolator => centers_to_faces_1D_t(order=self%order_, cells=self%cells_, dx=dx))
+          vector_1D = vector_1D_t( &
+              tensor_1D_t(interpolator .center. Dv(2:size(Dv)-1), self%x_min_, self%x_max_, self%cells_, self%order_) &
+            ,divergence_operator_1D_t(k=self%order_, dx=dx, cells=self%cells_) &
+          )
+        end associate
 #if ASSERTIONS
-        associate( &
-           q  => divergence_1D%weights() &
-          ,dx => (self%x_max_ - self%x_min_)/self%cells_ &
-          ,b => [-1D0, [(0D0, center = 1, self%cells_-1)], 1D0] &
-        )
-          call_julienne_assert(.all. ([size(Dv), size(q)] .equalsExpected. self%cells_+2))
-          call_julienne_assert((.all. (matmul(transpose(D%assemble()), q) .approximates. b/dx .within. double_equivalence)))
-            ! Check D^T * a = b_{m+1},  Eq. (19), Corbino & Castillo (2020)
+        associate(divergence_1D => divergence_1D_t(tensor_1D_t(Dv(2:size(Dv)-1), self%x_min_, self%x_max_, self%cells_, self%order_)))
+          associate( &
+             q => divergence_1D%weights() &
+            ,b => [-1D0, [(0D0, center = 1, self%cells_-1)], 1D0] &
+          )
+            call_julienne_assert(.all. ([size(Dv), size(q)] .equalsExpected. self%cells_+2))
+            call_julienne_assert((.all. (matmul(transpose(D%assemble()), q) .approximates. b/dx .within. double_equivalence)))
+              ! Check D^T * a = b_{m+1},  Eq. (19), Corbino & Castillo (2020)
+          end associate
         end associate
 #endif
       end associate
diff --git a/src/formal/interpolation_operators_1D_m.F90 b/src/formal/interpolator_1D_m.F90
similarity index 52%
rename from src/formal/interpolation_operators_1D_m.F90
rename to src/formal/interpolator_1D_m.F90
index df8b156..38d85c0 100644
--- a/src/formal/interpolation_operators_1D_m.F90
+++ b/src/formal/interpolator_1D_m.F90
@@ -1,19 +1,18 @@
 ! Copyright (c) 2026, The Regents of the University of California
 ! Terms of use are as specified in LICENSE.txt
 
-module interpolation_operators_1D_m
+module interpolator_1D_m
   !! Define a sparse matrix storage format tailored to the staggered-grid interpolation matrix
   !! operators detailed by Dumett & Castillo (2022) https://doi.org/10.13140/RG.2.2.26630.14400
   implicit none
 
   private
-  public :: cells_to_faces_t
-  public :: faces_to_cells_t
+  public :: centers_to_faces_1D_t
 
-  type interpolation_operator_1D_t
+  type, abstract :: interpolator_1D_t
     !! Encapsulate a staggered-grid interpolation matrix with a corresponding matrix-vector product operator
     private
-    integer order
+    integer order_, cells_, dx_
     double precision                 first_
     double precision, allocatable :: upper_(:,:) 
     double precision, allocatable :: inner_(:)
@@ -21,15 +20,18 @@ module interpolation_operators_1D_m
     double precision                 final_
   end type
 
-  type, extends(interpolation_operator_1D_t) :: centers_to_faces_1D_t
+  type, extends(interpolator_1D_t) :: centers_to_faces_1D_t
+  contains
+    generic :: operator(.center.) => interpolate_to_faces
+    procedure, non_overridable, private :: interpolate_to_faces
   end type
 
-  type, extends(interpolation_operator_1D_t) :: faces_to_centers_1D_t
+  type, extends(interpolator_1D_t) :: faces_to_centers_1D_t
   end type
 
   interface centers_to_faces_1D_t
 
-    pure module function c2f_component_constructor(order, cells, dx) result(centers_to_faces_1D)
+    pure module function c2f_constructor(order, cells, dx) result(centers_to_faces_1D)
       !! Construct centers-to-faces interpolation operator
       implicit none
       integer, intent(in) :: order, cells
@@ -41,14 +43,26 @@ pure module function c2f_component_constructor(order, cells, dx) result(centers_
 
   interface faces_to_centers_1D_t
 
-    pure module function f2c_component_constructor(order, cells, dx) result(faces_to_centers_1D)
-      !! Construct faces-to-centers interpolation operator
+    pure module function f2c_constructor(order, cells, dx) result(faces_to_centers_1D)
+      !! Construct centers-to-faces interpolation operator
       implicit none
       integer, intent(in) :: order, cells
       double precision, intent(in) :: dx
-      type(faces_to_centers_1D_t) faces_to_centers_1D
+      type(centers_to_faces_1D_t) faces_to_centers_1D
+    end function
+
+  end interface
+
+  interface
+
+    pure module function interpolate_to_faces(self, centers) result(faces)
+      !! Interpolate cell-centered values to face-centered values
+      implicit none
+      class(centers_to_faces_1D_t), intent(in) :: self
+      double precision, intent(in) :: centers(:)
+      double precision, allocatable :: faces(:)
     end function
 
   end interface
 
-end module interpolation_operators_1D_m
+end module interpolator_1D_m
\ No newline at end of file
diff --git a/src/formal/interpolation_operators_1D_s.F90 b/src/formal/interpolator_1D_s.F90
similarity index 57%
rename from src/formal/interpolation_operators_1D_s.F90
rename to src/formal/interpolator_1D_s.F90
index 05549ab..ce5fe2a 100644
--- a/src/formal/interpolation_operators_1D_s.F90
+++ b/src/formal/interpolator_1D_s.F90
@@ -1,14 +1,19 @@
 ! Copyright (c) 2026, The Regents of the University of California
 ! Terms of use are as specified in LICENSE.txt
 
-submodule(interpolation_operators_1D_m) interpolation_operators_1D_s
+#include "julienne-assert-macros.h"
+
+submodule(interpolator_1D_m) interpolator_1D_s
+  use julienne_m, only : call_julienne_assert_, operator(.equalsExpected.)
   implicit none
 
 contains
 
-  module procedure c2f_component_constructor
+  module procedure c2f_constructor
 
+    centers_to_faces_1D%order_ = order
     centers_to_faces_1D%cells_ = cells
+    centers_to_faces_1D%dx_    = dx
 
     select case(order)
     case(2)
@@ -27,11 +32,14 @@
       error stop "c2f_component_constructor: unsupported order"
     end select
 
+    call_julienne_assert(shape(self%lower_) .equalsExpected. shape(self%upper_))
   end procedure
 
-  module procedure f2c_component_constructor
+  module procedure f2c_constructor
 
+    faces_to_centers_1D%order_ = order
     faces_to_centers_1D%cells_ = cells
+    faces_to_centers_1D%dx_    = dx
 
     select case(order)
     case(2)
@@ -50,6 +58,23 @@
       error stop "f2c_component_constructor: unsupported order"
     end select
 
+    call_julienne_assert(shape(self%lower_) .equalsExpected. shape(self%upper_))
+  end procedure
+
+  module procedure interpolate_to_faces
+    integer row
+    associate( &
+       N => size(centers)               , inner_cols => size(self%inner_) &
+      ,upper_rows => size(self%upper_,1), upper_cols => size(self%upper_,2) &
+      ,lower_rows => size(self%lower_,1), lower_cols => size(self%lower_,2) &
+    )
+      faces = [        self%first_ * centers(1) &
+        ,       matmul(self%upper_ , centers(1:upper_cols)) &
+        ,[(dot_product(self%inner_ , centers(row-upper_rows+1 : row-upper_rows+inner_cols)), row = upper_rows+1, N-lower_rows-1)] &
+        ,       matmul(self%lower_ , centers(N-lower_cols+1:N)) &
+        ,              self%final_ * centers(N) &
+      ]
+    end associate
   end procedure
 
-end submodule interpolation_operators_1D_s
+end submodule interpolator_1D_s
\ No newline at end of file
diff --git a/src/formal/tensors_1D_m.F90 b/src/formal/tensors_1D_m.F90
index 0249276..39ccc0e 100644
--- a/src/formal/tensors_1D_m.F90
+++ b/src/formal/tensors_1D_m.F90
@@ -311,8 +311,8 @@ pure module function reduced_order_boundary_depth(self) result(num_nodes)
       integer num_nodes
     end function
 
-    pure module function div_dyad(self) result(divergence_1D)
-      !! Result is mimetic divergence of the dyad_1D_t "self"
+    pure module function div_dyad(self) result(vector_1D)
+      !! Result is the vector divergence of the dyad_1D_t "self"
       implicit none
       class(dyad_1D_t), intent(in) :: self
       type(vector_1D_t) vector_1D

From 8f3258985fc0d6598e16ee60d81e0898c3cdb1cb Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Thu, 2 Apr 2026 23:18:06 -0700
Subject: [PATCH 13/17] feat(interpolation): centers-extended to faces

This commit adds the interpolate_1D_t type with a type-bound
procedure that interpolates using values at cell centers plus
boundaries to estimate values at faces using the 2nd- and
4th-order schemes of Dumett & Castillo (2022).
---
 src/formal/dyad_1D_s.F90         | 13 ++++-
 src/formal/interpolator_1D_m.F90 |  7 +--
 src/formal/interpolator_1D_s.F90 | 38 +++++++------
 src/formal_m.f90                 |  3 +
 test/driver.f90                  |  5 +-
 test/interpolator_1D_test_m.F90  | 97 ++++++++++++++++++++++++++++++++
 6 files changed, 139 insertions(+), 24 deletions(-)
 create mode 100644 test/interpolator_1D_test_m.F90

diff --git a/src/formal/dyad_1D_s.F90 b/src/formal/dyad_1D_s.F90
index 3b685e4..260ef3b 100644
--- a/src/formal/dyad_1D_s.F90
+++ b/src/formal/dyad_1D_s.F90
@@ -4,9 +4,17 @@
 #include "julienne-assert-macros.h"
 
 submodule(tensors_1D_m) dyad_1D_s
-  use julienne_m, only : call_julienne_assert_
+  use julienne_m, only : &
+     call_julienne_assert_ &
+    ,operator(.all.) &
+    ,operator(.approximates.) &
+    ,operator(.equalsExpected.) &
+    ,operator(.within.)
   use interpolator_1D_m, only : centers_to_faces_1D_t
   implicit none
+
+  double precision, parameter :: double_equivalence = 2D-4
+
 contains
 
   module procedure dyad_over_integer
@@ -28,14 +36,13 @@
 #else
     associate(D => (self%divergence_operator_1D_))
 #endif
-      call_julienne_assert(size(self%values_))
       associate( &
          Dv => D .x. self%values_ &
         ,dx => (self%x_max_ - self%x_min_)/self%cells_ &
       )
         associate(interpolator => centers_to_faces_1D_t(order=self%order_, cells=self%cells_, dx=dx))
           vector_1D = vector_1D_t( &
-              tensor_1D_t(interpolator .center. Dv(2:size(Dv)-1), self%x_min_, self%x_max_, self%cells_, self%order_) &
+              tensor_1D_t(interpolator%face_values(Dv(2:size(Dv)-1)), self%x_min_, self%x_max_, self%cells_, self%order_) &
             ,divergence_operator_1D_t(k=self%order_, dx=dx, cells=self%cells_) &
           )
         end associate
diff --git a/src/formal/interpolator_1D_m.F90 b/src/formal/interpolator_1D_m.F90
index 38d85c0..e8a2fb0 100644
--- a/src/formal/interpolator_1D_m.F90
+++ b/src/formal/interpolator_1D_m.F90
@@ -22,8 +22,7 @@ module interpolator_1D_m
 
   type, extends(interpolator_1D_t) :: centers_to_faces_1D_t
   contains
-    generic :: operator(.center.) => interpolate_to_faces
-    procedure, non_overridable, private :: interpolate_to_faces
+    procedure, non_overridable :: face_values
   end type
 
   type, extends(interpolator_1D_t) :: faces_to_centers_1D_t
@@ -55,11 +54,11 @@ pure module function f2c_constructor(order, cells, dx) result(faces_to_centers_1
 
   interface
 
-    pure module function interpolate_to_faces(self, centers) result(faces)
+    pure module function face_values(self, centers_extended) result(faces)
       !! Interpolate cell-centered values to face-centered values
       implicit none
       class(centers_to_faces_1D_t), intent(in) :: self
-      double precision, intent(in) :: centers(:)
+      double precision, intent(in) :: centers_extended(:)
       double precision, allocatable :: faces(:)
     end function
 
diff --git a/src/formal/interpolator_1D_s.F90 b/src/formal/interpolator_1D_s.F90
index ce5fe2a..2cc7559 100644
--- a/src/formal/interpolator_1D_s.F90
+++ b/src/formal/interpolator_1D_s.F90
@@ -4,7 +4,7 @@
 #include "julienne-assert-macros.h"
 
 submodule(interpolator_1D_m) interpolator_1D_s
-  use julienne_m, only : call_julienne_assert_, operator(.equalsExpected.)
+  use julienne_m, only : call_julienne_assert_, operator(.all.), operator(.equalsExpected.)
   implicit none
 
 contains
@@ -17,11 +17,11 @@
 
     select case(order)
     case(2)
-      centers_to_faces_1D%first_ = 1D0
-      centers_to_faces_1D%upper_ = reshape([double precision::], [0,3])
-      centers_to_faces_1D%inner_ = [1D0,1D0]/2D0
-      centers_to_faces_1D%lower_ = reshape([double precision::], [0,3])
-      centers_to_faces_1D%final_ = 1D0
+      centers_to_faces_1D%first_ = (2D0                                 )/2
+      centers_to_faces_1D%upper_ = (reshape([double precision::], [0,3]))/2
+      centers_to_faces_1D%inner_ = ([1D0,1D0]                           )/2
+      centers_to_faces_1D%lower_ = (reshape([double precision::], [0,3]))/2
+      centers_to_faces_1D%final_ = (2D0                                 )/2
     case(4)
       centers_to_faces_1D%first_ = 1D0
       centers_to_faces_1D%upper_ = reshape([-16,  70, 70, -14,   2], [1,5])/112D0
@@ -32,7 +32,7 @@
       error stop "c2f_component_constructor: unsupported order"
     end select
 
-    call_julienne_assert(shape(self%lower_) .equalsExpected. shape(self%upper_))
+    call_julienne_assert(.all. (shape(centers_to_faces_1D%lower_) .equalsExpected. shape(centers_to_faces_1D%upper_)))
   end procedure
 
   module procedure f2c_constructor
@@ -58,22 +58,28 @@
       error stop "f2c_component_constructor: unsupported order"
     end select
 
-    call_julienne_assert(shape(self%lower_) .equalsExpected. shape(self%upper_))
+    call_julienne_assert(.all. (shape(faces_to_centers_1D%lower_) .equalsExpected. shape(faces_to_centers_1D%upper_)))
   end procedure
 
-  module procedure interpolate_to_faces
+  module procedure face_values
     integer row
+    integer, parameter :: end_point = 1
     associate( &
-       N => size(centers)               , inner_cols => size(self%inner_) &
+       N => size(centers_extended)      , inner_cols => size(self%inner_) &
       ,upper_rows => size(self%upper_,1), upper_cols => size(self%upper_,2) &
       ,lower_rows => size(self%lower_,1), lower_cols => size(self%lower_,2) &
     )
-      faces = [        self%first_ * centers(1) &
-        ,       matmul(self%upper_ , centers(1:upper_cols)) &
-        ,[(dot_product(self%inner_ , centers(row-upper_rows+1 : row-upper_rows+inner_cols)), row = upper_rows+1, N-lower_rows-1)] &
-        ,       matmul(self%lower_ , centers(N-lower_cols+1:N)) &
-        ,              self%final_ * centers(N) &
-      ]
+      call_julienne_assert(N .equalsExpected. self%cells_ + 2*end_point)
+      associate(inner_rows => N - (2*end_point + upper_rows + lower_rows))
+        faces = [        self%first_ * centers_extended(1) &
+          ,       matmul(self%upper_ , centers_extended(1:upper_cols)) &
+          ,[(dot_product(self%inner_ , centers_extended(row - upper_rows : row - upper_rows + inner_cols - 1)),  &
+            row = end_point + upper_rows + 1, end_point + upper_rows + inner_rows - 1)] &
+          ,       matmul(self%lower_ , centers_extended(N-lower_cols+1:N)) &
+          ,              self%final_ * centers_extended(N) &
+        ]
+      end associate
+      call_julienne_assert(size(faces) .equalsExpected. self%cells_ + 1)
     end associate
   end procedure
 
diff --git a/src/formal_m.f90 b/src/formal_m.f90
index 8128923..14e481f 100644
--- a/src/formal_m.f90
+++ b/src/formal_m.f90
@@ -20,6 +20,9 @@ module formal_m
      gradient_operator_1D_t & ! matrix operator defining a 1D mimetic gradient
     ,divergence_operator_1D_t ! matrix operator defining a 1D mimetic divergence
 
+  use interpolator_1D_m, only : &
+      centers_to_faces_1D_t ! 1D mimetic estimator producing cell-centered values given facial values
+
   implicit none
 
 end module formal_m
diff --git a/test/driver.f90 b/test/driver.f90
index 5c12c1d..72454f5 100644
--- a/test/driver.f90
+++ b/test/driver.f90
@@ -2,11 +2,13 @@
 ! Terms of use are as specified in LICENSE.txt
 
 program test_suite_driver
+  !! Run the test suite and report results
   use julienne_m, only : test_fixture_t, test_harness_t
   use gradient_operator_1D_test_m, only : gradient_operator_1D_test_t
   use divergence_operator_1D_test_m, only : divergence_operator_1D_test_t
   use laplacian_operator_1D_test_m, only : laplacian_operator_1D_test_t
   use integration_operators_1D_test_m, only : integration_operators_1D_test_t
+  use interpolator_1D_test_m, only : interpolator_1D_test_t
   implicit none
 
   associate(test_harness => test_harness_t([ &
@@ -14,7 +16,8 @@ program test_suite_driver
     ,test_fixture_t(divergence_operator_1D_test_t()) &
     ,test_fixture_t(laplacian_operator_1D_test_t()) &
     ,test_fixture_t(integration_operators_1D_test_t()) &
+    ,test_fixture_t(interpolator_1D_test_t()) &
   ]))
     call test_harness%report_results
   end associate
-end program test_suite_driver
+end program test_suite_driver
\ No newline at end of file
diff --git a/test/interpolator_1D_test_m.F90 b/test/interpolator_1D_test_m.F90
new file mode 100644
index 0000000..bc1360f
--- /dev/null
+++ b/test/interpolator_1D_test_m.F90
@@ -0,0 +1,97 @@
+! Copyright (c) 2026, The Regents of the University of California
+! Terms of use are as specified in LICENSE.txt
+
+module interpolator_1D_test_m
+  use julienne_m, only : &
+     operator(//) &
+    ,operator(.all.) &
+    ,operator(.also.) &
+    ,operator(.approximates.) &
+    ,operator(.within.) &
+    ,passing_test &
+    ,string_t &
+    ,test_description_t &
+    ,test_diagnosis_t &
+    ,test_result_t &
+    ,test_t &
+    ,usher
+  use formal_m, only : centers_to_faces_1D_t, scalar_1D_t, scalar_1D_initializer_i
+
+  implicit none
+
+  type, extends(test_t) :: interpolator_1D_test_t
+  contains
+    procedure, nopass :: subject
+    procedure, nopass :: results
+  end type
+
+  double precision, parameter :: tolerance = 1D-11
+
+contains
+
+  pure function subject() result(test_subject)
+    character(len=:), allocatable :: test_subject
+    test_subject = 'A 1D mimetic interpolator'
+  end function
+
+  function results() result(test_results)
+   type(interpolator_1D_test_t) interpolator_1D_test
+   type(test_result_t), allocatable :: test_results(:)
+
+   test_results = interpolator_1D_test%run([ & 
+      test_description_t('estimating values at cell centers given face values with 2nd- & 4th-order interpolators', usher(check_centers_to_faces)) &
+   ])
+  end function
+
+  pure function line(x) result(y)
+    double precision, intent(in) :: x(:)
+    double precision, allocatable :: y(:)
+    y = 3*x + 5
+  end function
+
+  pure function cubic(x) result(y)
+    double precision, intent(in) :: x(:)
+    double precision, allocatable :: y(:)
+    y = 7*x**3 + 4*x**2 + x + 2
+  end function
+
+  function check_centers_to_faces() result(test_diagnosis)
+    type(test_diagnosis_t) test_diagnosis
+    procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => null()
+    double precision, parameter :: x_min = 0D0, x_max = 20D0
+    integer order, cells
+
+    test_diagnosis = passing_test()
+
+    do order = 2, 4, 2
+
+      select case(order)
+      case(2)
+        cells = 10
+        scalar_1D_initializer => line
+      case(4)
+        cells = 20
+        scalar_1D_initializer => cubic
+      case default
+        error stop "check_centers_to_faces (interpolator_1D_test_m) unsupported order of accuracy"
+      end select
+
+      associate( &
+         scalar_1D => scalar_1D_t(scalar_1D_initializer, order=order, cells=cells, x_min=x_min, x_max=x_max) &
+        ,interpolator => centers_to_faces_1D_t(order=order, cells=cells, dx=(x_max - x_min)/cells) &
+      )
+        associate( &
+           scalar_at_faces => interpolator%face_values(scalar_1D%values()) &
+          ,gradient => .grad. scalar_1D &
+        )
+          associate(face_locations => gradient%grid())
+            test_diagnosis = test_diagnosis .also. .all. (scalar_at_faces .approximates. scalar_1D_initializer(face_locations) .within. tolerance) &
+              // string_t(" for order ") // string_t(order)
+          end associate
+        end associate
+      end associate
+
+    end do
+  end function
+
+end module interpolator_1D_test_m
\ No newline at end of file

From 0d72b08ef9957cc06da1d96c9d30d63965bb0d15 Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Sat, 4 Apr 2026 14:02:36 -0700
Subject: [PATCH 14/17] feat(interpolation): faces to centers-extended

This commit adds the faces_to_centers_t type with a type-bound
procedure that interpolates using values at the faces to estimate
to estimate values at cell centers plus boundaries using the 2nd-
and 4th-order schemes of Dumett & Castillo (2022).
---
 src/formal/interpolator_1D_m.F90 | 13 +++++++-
 src/formal/interpolator_1D_s.F90 | 32 +++++++++++++++---
 src/formal_m.f90                 |  3 +-
 test/interpolator_1D_test_m.F90  | 56 ++++++++++++++++++++++++++++----
 4 files changed, 91 insertions(+), 13 deletions(-)

diff --git a/src/formal/interpolator_1D_m.F90 b/src/formal/interpolator_1D_m.F90
index e8a2fb0..537c6c4 100644
--- a/src/formal/interpolator_1D_m.F90
+++ b/src/formal/interpolator_1D_m.F90
@@ -8,6 +8,7 @@ module interpolator_1D_m
 
   private
   public :: centers_to_faces_1D_t
+  public :: faces_to_centers_1D_t
 
   type, abstract :: interpolator_1D_t
     !! Encapsulate a staggered-grid interpolation matrix with a corresponding matrix-vector product operator
@@ -26,6 +27,8 @@ module interpolator_1D_m
   end type
 
   type, extends(interpolator_1D_t) :: faces_to_centers_1D_t
+  contains
+    procedure, non_overridable :: center_values
   end type
 
   interface centers_to_faces_1D_t
@@ -47,7 +50,7 @@ pure module function f2c_constructor(order, cells, dx) result(faces_to_centers_1
       implicit none
       integer, intent(in) :: order, cells
       double precision, intent(in) :: dx
-      type(centers_to_faces_1D_t) faces_to_centers_1D
+      type(faces_to_centers_1D_t) faces_to_centers_1D
     end function
 
   end interface
@@ -62,6 +65,14 @@ pure module function face_values(self, centers_extended) result(faces)
       double precision, allocatable :: faces(:)
     end function
 
+    pure module function center_values(self, faces) result(centers_extended)
+      !! Interpolate face-centered values to cell-centered values
+      implicit none
+      class(faces_to_centers_1D_t), intent(in) :: self
+      double precision, intent(in) :: faces(:)
+      double precision, allocatable :: centers_extended(:)
+    end function
+
   end interface
 
 end module interpolator_1D_m
\ No newline at end of file
diff --git a/src/formal/interpolator_1D_s.F90 b/src/formal/interpolator_1D_s.F90
index 2cc7559..d7e5592 100644
--- a/src/formal/interpolator_1D_s.F90
+++ b/src/formal/interpolator_1D_s.F90
@@ -43,11 +43,11 @@
 
     select case(order)
     case(2)
-      faces_to_centers_1D%first_ = 1D0
-      faces_to_centers_1D%upper_ = reshape([double precision::], [0,3])
-      faces_to_centers_1D%inner_ = [1D0,1D0]/2D0
-      faces_to_centers_1D%lower_ = reshape([double precision::], [0,3])
-      faces_to_centers_1D%final_ = 1D0
+      faces_to_centers_1D%first_ = (2D0                               )/2
+      faces_to_centers_1D%upper_ = reshape([double precision::], [0,3])/2
+      faces_to_centers_1D%inner_ =                            [1D0,1D0]/2
+      faces_to_centers_1D%lower_ = reshape([double precision::], [0,3])/2
+      faces_to_centers_1D%final_ = (2D0                               )/2
     case(4)
       faces_to_centers_1D%first_ = 1D0
       faces_to_centers_1D%upper_ = reshape([35, 140, -70, 28, -5], [1,5])/128D0
@@ -83,4 +83,26 @@
     end associate
   end procedure
 
+  module procedure center_values
+    integer row
+    integer, parameter :: end_point = 1
+    associate( &
+       N => size(faces)                 , inner_cols => size(self%inner_) &
+      ,upper_rows => size(self%upper_,1), upper_cols => size(self%upper_,2) &
+      ,lower_rows => size(self%lower_,1), lower_cols => size(self%lower_,2) &
+    )
+      call_julienne_assert(N .equalsExpected. self%cells_ + 1)
+      associate(inner_rows => N + 1 - (2*end_point + upper_rows + lower_rows))
+        centers_extended = &
+          [               self%first_ * faces(1) &
+           ,       matmul(self%upper_ , faces(1:upper_cols)) &
+           ,[(dot_product(self%inner_ , faces(row : row + inner_cols - 1)), row = 1, inner_rows )] &
+           ,       matmul(self%lower_ , faces(N-lower_cols+1:N)) &
+           ,              self%final_ * faces(N) &
+          ]
+      end associate
+      call_julienne_assert(size(centers_extended) .equalsExpected. self%cells_ + 2)
+    end associate
+  end procedure
+
 end submodule interpolator_1D_s
\ No newline at end of file
diff --git a/src/formal_m.f90 b/src/formal_m.f90
index 14e481f..7f507d9 100644
--- a/src/formal_m.f90
+++ b/src/formal_m.f90
@@ -21,7 +21,8 @@ module formal_m
     ,divergence_operator_1D_t ! matrix operator defining a 1D mimetic divergence
 
   use interpolator_1D_m, only : &
-      centers_to_faces_1D_t ! 1D mimetic estimator producing cell-centered values given facial values
+     centers_to_faces_1D_t & ! 1D mimetic interpolator producing cell-centered values from face-centered values
+    ,faces_to_centers_1D_t   ! 1D mimetic interpolator producing face-centered values from cell-centered values
 
   implicit none
 
diff --git a/test/interpolator_1D_test_m.F90 b/test/interpolator_1D_test_m.F90
index bc1360f..4cd32bd 100644
--- a/test/interpolator_1D_test_m.F90
+++ b/test/interpolator_1D_test_m.F90
@@ -14,8 +14,11 @@ module interpolator_1D_test_m
     ,test_diagnosis_t &
     ,test_result_t &
     ,test_t &
-    ,usher
-  use formal_m, only : centers_to_faces_1D_t, scalar_1D_t, scalar_1D_initializer_i
+    ,usher &
+    ,csv  
+  use formal_m, only : &
+     centers_to_faces_1D_t, scalar_1D_t, scalar_1D_initializer_i &
+    ,faces_to_centers_1D_t, vector_1D_t, vector_1D_initializer_i
 
   implicit none
 
@@ -39,14 +42,15 @@ function results() result(test_results)
    type(test_result_t), allocatable :: test_results(:)
 
    test_results = interpolator_1D_test%run([ & 
-      test_description_t('estimating values at cell centers given face values with 2nd- & 4th-order interpolators', usher(check_centers_to_faces)) &
+      test_description_t('estimating center values from face values with 2nd- & 4th-order interpolators', usher(check_centers_to_faces)) &
+     ,test_description_t('estimating face values from center values with 2nd- & 4th-order interpolators', usher(check_faces_to_centers)) &
    ])
   end function
 
   pure function line(x) result(y)
     double precision, intent(in) :: x(:)
     double precision, allocatable :: y(:)
-    y = 3*x + 5
+    y = x
   end function
 
   pure function cubic(x) result(y)
@@ -54,7 +58,6 @@ pure function cubic(x) result(y)
     double precision, allocatable :: y(:)
     y = 7*x**3 + 4*x**2 + x + 2
   end function
-
   function check_centers_to_faces() result(test_diagnosis)
     type(test_diagnosis_t) test_diagnosis
     procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => null()
@@ -94,4 +97,45 @@ function check_centers_to_faces() result(test_diagnosis)
     end do
   end function
 
-end module interpolator_1D_test_m
\ No newline at end of file
+  function check_faces_to_centers() result(test_diagnosis)
+    type(test_diagnosis_t) test_diagnosis
+    procedure(vector_1D_initializer_i), pointer :: vector_1D_initializer => null()
+    double precision, parameter :: x_min = 0D0, x_max = 20D0
+    integer order, cells
+
+    test_diagnosis = passing_test()
+
+    do order = 2, 4, 2
+
+      select case(order)
+      case(2)
+        cells = 10
+        vector_1D_initializer => line
+      case(4)
+        cells = 20
+        vector_1D_initializer => cubic
+      case default
+        error stop "check_faces_to_centers (interpolator_1D_test_m) unsupported order of accuracy"
+      end select
+
+      associate( &
+         vector_1D => vector_1D_t(vector_1D_initializer, order=order, cells=cells, x_min=x_min, x_max=x_max) &
+        ,interpolator => faces_to_centers_1D_t(order=order, cells=cells, dx=(x_max - x_min)/cells) &
+      )
+        associate( &
+           vector_at_centers => interpolator%center_values(vector_1D%values()) &
+          ,divergence => .div. vector_1D &
+          ,faces => vector_1D%grid() &
+        )
+          associate( centers_extended => [faces(lbound(faces,1)), divergence%grid(), faces(ubound(faces,1))] )
+            test_diagnosis = test_diagnosis .also. &
+              .all. (vector_at_centers .approximates. vector_1D_initializer(centers_extended) .within. tolerance) &
+              // string_t(" for order ") // string_t(order)
+          end associate
+        end associate
+      end associate
+
+    end do
+  end function
+
+end module interpolator_1D_test_m

From 452c648367095ddc7edae6784656a9528b7338d4 Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Sun, 5 Apr 2026 16:14:57 -0600
Subject: [PATCH 15/17] fix(burgers-1D): work around gfortran issue

---
 example/burgers-1D.F90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/example/burgers-1D.F90 b/example/burgers-1D.F90
index d7d7085..66ac30c 100644
--- a/example/burgers-1D.F90
+++ b/example/burgers-1D.F90
@@ -53,7 +53,7 @@ program burgers_1D
   if (len(order_string)==0) then 
     order = 4 
   else
-    read(order_string,"(i)") order 
+    read(order_string,"(i1)") order 
   end if
 
   print *, "order = ", order

From 1961d9d994c40890ca850b6d113bf9dba4c1b575 Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Sun, 5 Apr 2026 16:29:55 -0600
Subject: [PATCH 16/17] fix(interpolator_test): work around gfortran issue

---
 test/interpolator_1D_test_m.F90 | 12 +++++++-----
 1 file changed, 7 insertions(+), 5 deletions(-)

diff --git a/test/interpolator_1D_test_m.F90 b/test/interpolator_1D_test_m.F90
index 4cd32bd..9262c9f 100644
--- a/test/interpolator_1D_test_m.F90
+++ b/test/interpolator_1D_test_m.F90
@@ -58,9 +58,11 @@ pure function cubic(x) result(y)
     double precision, allocatable :: y(:)
     y = 7*x**3 + 4*x**2 + x + 2
   end function
+
   function check_centers_to_faces() result(test_diagnosis)
     type(test_diagnosis_t) test_diagnosis
     procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => null()
+    procedure(vector_1D_initializer_i), pointer :: vector_1D_initializer => line
     double precision, parameter :: x_min = 0D0, x_max = 20D0
     integer order, cells
 
@@ -81,16 +83,15 @@ function check_centers_to_faces() result(test_diagnosis)
 
       associate( &
          scalar_1D => scalar_1D_t(scalar_1D_initializer, order=order, cells=cells, x_min=x_min, x_max=x_max) &
+        ,vector_1D => vector_1D_t(vector_1D_initializer, order=order, cells=cells, x_min=x_min, x_max=x_max) &
         ,interpolator => centers_to_faces_1D_t(order=order, cells=cells, dx=(x_max - x_min)/cells) &
       )
         associate( &
            scalar_at_faces => interpolator%face_values(scalar_1D%values()) &
-          ,gradient => .grad. scalar_1D &
+          ,face_locations => vector_1D%grid() &
         )
-          associate(face_locations => gradient%grid())
             test_diagnosis = test_diagnosis .also. .all. (scalar_at_faces .approximates. scalar_1D_initializer(face_locations) .within. tolerance) &
               // string_t(" for order ") // string_t(order)
-          end associate
         end associate
       end associate
 
@@ -100,6 +101,7 @@ function check_centers_to_faces() result(test_diagnosis)
   function check_faces_to_centers() result(test_diagnosis)
     type(test_diagnosis_t) test_diagnosis
     procedure(vector_1D_initializer_i), pointer :: vector_1D_initializer => null()
+    procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => line
     double precision, parameter :: x_min = 0D0, x_max = 20D0
     integer order, cells
 
@@ -120,14 +122,14 @@ function check_faces_to_centers() result(test_diagnosis)
 
       associate( &
          vector_1D => vector_1D_t(vector_1D_initializer, order=order, cells=cells, x_min=x_min, x_max=x_max) &
+        ,scalar_1D => scalar_1D_t(scalar_1D_initializer, order=order, cells=cells, x_min=x_min, x_max=x_max) &
         ,interpolator => faces_to_centers_1D_t(order=order, cells=cells, dx=(x_max - x_min)/cells) &
       )
         associate( &
            vector_at_centers => interpolator%center_values(vector_1D%values()) &
-          ,divergence => .div. vector_1D &
           ,faces => vector_1D%grid() &
         )
-          associate( centers_extended => [faces(lbound(faces,1)), divergence%grid(), faces(ubound(faces,1))] )
+          associate( centers_extended => scalar_1D%grid() )
             test_diagnosis = test_diagnosis .also. &
               .all. (vector_at_centers .approximates. vector_1D_initializer(centers_extended) .within. tolerance) &
               // string_t(" for order ") // string_t(order)

From 1da22a54f03e23a7872336b9020cb3ab3ad27f20 Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Sun, 5 Apr 2026 17:09:41 -0600
Subject: [PATCH 17/17] fix(dyad_1D_s): work around lfortran issue

---
 src/formal/dyad_1D_s.F90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/formal/dyad_1D_s.F90 b/src/formal/dyad_1D_s.F90
index 260ef3b..50c3ed5 100644
--- a/src/formal/dyad_1D_s.F90
+++ b/src/formal/dyad_1D_s.F90
@@ -18,7 +18,7 @@
 contains
 
   module procedure dyad_over_integer
-    ratio%tensor_1D_t = tensor_1D_t(self%tensor_1D_t%values_/numerator, self%x_min_, self%x_max_, self%cells_, order = self%order_)
+    ratio%tensor_1D_t = tensor_1D_t(self%values_/numerator, self%x_min_, self%x_max_, self%cells_, order = self%order_)
     ratio%divergence_operator_1D_ = self%divergence_operator_1D_
   end procedure