Unify conformability checks and fix scalar-vector products

src/formal/divergence_2D_s.F90

@@ -16,18 +16,6 @@
 
 contains
 
-  module procedure divergence_2D_conformable_scalar
-    call_julienne_assert(scalar_2D_consistent(scalar_2D))
-    call_julienne_assert(self%tensor_2D_conformable(scalar_2D))
-    conformable = .true.
-  end procedure
-
-  module procedure divergence_2D_conformable_vector
-    call_julienne_assert(vector_2D_consistent(vector_2D))
-    call_julienne_assert(self%tensor_2D_conformable(vector_2D))
-    conformable = .true.
-  end procedure
-
   module procedure divergence_2D_values
     call_julienne_assert(self%consistent())
     divergences = self%points_(1,1,1,1)%values_(:,:)

src/formal/scalar_2D_s.F90

@@ -15,6 +15,7 @@
     ,operator(.isAtLeast.) &
     ,string_t
   use tensors_1D_m, only : cell_centers_extended_1D, scalar_1D_t
+  use interpolator_1D_m, only : centers_to_faces_1D_t
   implicit none
 
 contains
@@ -26,13 +27,6 @@
     self_consistent = .true.
   end procedure
 
-  module procedure scalar_2D_conformable_scalar
-    call_julienne_assert(scalar_2D_consistent(self))
-    call_julienne_assert(scalar_2D_consistent(scalar_2D))
-    call_julienne_assert(self%tensor_2D_conformable(scalar_2D))
-    conformable = .true.
-  end procedure
-
   module procedure construct_2D_scalar_from_components
 
     call_julienne_assert(size(gradient_operator_1D) .equalsExpected. space_dimension)
@@ -170,9 +164,11 @@
 
   module procedure scalar_2D_assign_divergence
 
-     call_julienne_assert(lhs%conformable(rhs))
+     call_julienne_assert(rhs%consistent())
 
+     if (allocated(lhs%points_)) deallocate(lhs%points_)
      allocate(lhs%points_(1,1,1,1))
+     if (allocated(lhs%points_(1,1,1,1)%values_)) deallocate(lhs%points_(1,1,1,1)%values_)
      allocate(lhs%points_(1,1,1,1)%values_(rhs%cells_(x_dir)+2, rhs%cells_(y_dir)+2))
 
      associate( &
@@ -192,6 +188,7 @@
      lhs%order_ = rhs%order_
 
      call_julienne_assert(lhs%consistent())
+     call_julienne_assert(lhs%conformable(rhs))
 
   end procedure
 
@@ -235,4 +232,39 @@
     file = file_t(lines)
   end procedure
 
+  module procedure scalar_2D_to_faces
+
+    call_julienne_assert(self%consistent())
+
+    construct_interpolator_array: &
+    associate(interpolator => centers_to_faces_1D_t(order=self%order_, cells=self%cells_, dx=(self%x_max_ - self%x_min_)/self%cells_))
+
+      select case(direction)
+
+      case(x_dir)
+
+        allocate(scalars(self%cells_(x_dir)+1, self%cells_(y_dir)+2))
+
+        interpolate_centers_to_x_faces: &
+        do concurrent(integer :: j = 1:size(self%points_(1,1,1,1)%values_,y_dir))
+          scalars(:,j) = interpolator(x_dir)%face_values(self%points_(1,1,1,1)%values_(:,j))
+        end do interpolate_centers_to_x_faces
+
+      case(y_dir)
+
+        allocate(scalars(self%cells_(x_dir)+2, self%cells_(y_dir)+1))
+
+        interpolate_centers_to_y_faces: &
+        do concurrent(integer :: i = 1:size(self%points_(1,1,1,1)%values_,x_dir))
+          scalars(i,:) = interpolator(y_dir)%face_values(self%points_(1,1,1,1)%values_(i,:))
+        end do interpolate_centers_to_y_faces
+
+      case default
+        error stop "scalar_2D_to_faces in scalar_2D_s: invalid direction"
+      end select
+
+    end associate construct_interpolator_array
+
+  end procedure scalar_2D_to_faces
+
 end submodule scalar_2D_s

src/formal/tensors_2D_m.F90

@@ -17,8 +17,10 @@ module tensors_2D_m
   public :: scalar_2D_initializer_i
   public :: vector_2D_initializer_i
   public :: divergence_2D_initializer_i
+  public :: x_dir
+  public :: y_dir
 
-  integer, parameter :: space_dimension = 2, max_tensor_rank = 4, x_dir = 1, y_dir = 2, z_dir = 3
+  integer, parameter :: space_dimension = 2, max_tensor_rank = 4, x_dir = 1, y_dir = 2
 
   abstract interface
 
@@ -61,6 +63,7 @@ pure function vector_2D_initializer_i(x,y) result(v)
     integer cells_(space_dimension) !! number of grid cells spanning the domain
     integer order_ !! order of accuracy of mimetic discretization
   contains
+    generic :: conformable => tensor_2D_conformable
     procedure, non_overridable, private :: tensor_rank
     procedure, non_overridable, private :: tensor_2D_consistent
     procedure, non_overridable, private :: tensor_2D_conformable
@@ -92,15 +95,15 @@ pure module function construct_2D_tensor_from_components(points, cells, x_min, x
     generic :: values => scalar_2D_values
     generic :: grid => scalar_2D_grid
     generic :: consistent => scalar_2D_consistent
-    generic :: conformable => scalar_2D_conformable_scalar
+    generic :: to_faces => scalar_2D_to_faces
     generic :: to_file => scalar_2D_to_file
     procedure, non_overridable, private :: scalar_2D_assign_divergence
     procedure, non_overridable, private :: scalar_2D_to_file
+    procedure, non_overridable, private :: scalar_2D_to_faces
     procedure, non_overridable, private :: scalar_2D_gradient
     procedure, non_overridable, private :: scalar_2D_values
     procedure, non_overridable, private :: scalar_2D_grid
     procedure, non_overridable, private :: scalar_2D_consistent
-    procedure, non_overridable, private :: scalar_2D_conformable_scalar
     procedure, non_overridable, private :: scalar_2D_postmultiply_double
     procedure, non_overridable, private :: scalar_2D_plus_scalar
     procedure, non_overridable, private, pass(rhs) :: scalar_2D_premultiply_double
@@ -143,19 +146,18 @@ pure module function construct_2D_scalar_from_components(tensor_2D, gradient_ope
     type(divergence_operator_1D_t) divergence_operator_1D_(space_dimension)
   contains
     generic :: grid => vector_2D_grid
+    generic :: values => vector_2D_values
     generic :: consistent => vector_2D_consistent
-    generic :: conformable => vector_2D_conformable_vector, vector_2D_conformable_scalar
     generic :: to_centers_extended => vector_2D_to_centers_extended
     generic :: operator(.div.) => vector_2D_divergence
     generic :: operator(.dot.) => vector_2D_dot_vector
     generic :: operator(*) => vector_2D_postmultiply_scalar, vector_2D_premultiply_scalar
     generic :: to_file => vector_2D_to_file
     procedure, non_overridable, private :: vector_2D_to_file
     procedure, non_overridable, private :: vector_2D_grid
+    procedure, non_overridable, private :: vector_2D_values
     procedure, non_overridable, private :: vector_2D_divergence
     procedure, non_overridable, private :: vector_2D_consistent
-    procedure, non_overridable, private :: vector_2D_conformable_vector
-    procedure, non_overridable, private :: vector_2D_conformable_scalar
     procedure, non_overridable, private :: vector_2D_to_centers_extended
     procedure, non_overridable, private :: vector_2D_dot_vector
     procedure, non_overridable, private :: vector_2D_postmultiply_scalar
@@ -230,16 +232,13 @@ pure module function construct_2D_gradient_from_components(tensor_2D, divergence
     generic :: values => divergence_2D_values
     generic :: grid => divergence_2D_grid
     generic :: consistent => tensor_2D_consistent
-    generic :: conformable => divergence_2D_conformable_scalar, divergence_2D_conformable_vector
     generic :: operator(*) => divergence_2D_premultiply_constant, divergence_2D_postmultiply_constant
     generic :: operator(-) => divergence_2D_minus_scalar, divergence_2D_minus_divergence
     generic :: to_file => divergence_2D_to_file
     procedure, non_overridable, private :: divergence_2D_to_file
     procedure, private, non_overridable :: divergence_2D_values
     procedure, private, non_overridable :: divergence_2D_grid
     procedure, private, non_overridable :: divergence_2D_minus_scalar
-    procedure, private, non_overridable :: divergence_2D_conformable_vector
-    procedure, private, non_overridable :: divergence_2D_conformable_scalar
     procedure, private, non_overridable :: divergence_2D_postmultiply_constant
     procedure, private, non_overridable :: divergence_2D_minus_divergence
     procedure, private, non_overridable, pass(rhs) :: divergence_2D_premultiply_constant
@@ -298,49 +297,12 @@ pure module function scalar_2D_consistent(self) result(self_consistent)
       logical self_consistent
     end function
 
-    pure module function scalar_2D_conformable_scalar(self, scalar_2D) result(conformable)
-      !! Assert the arguments' components are conformable, self-consistent, and consistent with each other
+    pure module function scalar_2D_to_faces(self, direction) result(scalars)
+      !! Result is the scalar values interpolated form centers-extended to faces along the requested direction
       implicit none
-      class(scalar_2D_t), intent(in) :: self, scalar_2D
-      logical conformable
-    end function
-
-    pure module function divergence_2D_conformable_scalar(self, scalar_2D) result(conformable)
-      !! Assert the arguments' components are conformable, self-consistent, and consistent with each other
-      implicit none
-      class(divergence_2D_t), intent(in) :: self
-      class(scalar_2D_t), intent(in) :: scalar_2D
-      logical conformable
-    end function
-
-    pure module function divergence_2D_conformable_vector(self, vector_2D) result(conformable)
-      !! Assert the arguments' components are conformable, self-consistent, and consistent with each other
-      implicit none
-      class(divergence_2D_t), intent(in) :: self
-      class(vector_2D_t), intent(in) :: vector_2D
-      logical conformable
-    end function
-
-    pure module function vector_2D_consistent(self) result(self_consistent)
-      !! Assert components allocated and self-consistent, including sufficient accuracy for divergence operator
-      implicit none
-      class(vector_2D_t), intent(in) :: self
-      logical self_consistent
-    end function
-
-    pure module function vector_2D_conformable_vector(self, vector_2D) result(conformable)
-      !! Assert the arguments' components are conformable, self-consistent, and consistent with each other
-      implicit none
-      class(vector_2D_t), intent(in) :: self, vector_2D
-      logical conformable
-    end function
-
-    pure module function vector_2D_conformable_scalar(self, scalar_2D) result(conformable)
-      !! Assert components allocated and self-consistent
-      implicit none
-      class(vector_2D_t), intent(in) :: self
-      class(scalar_2D_t), intent(in) :: scalar_2D
-      logical conformable
+      class(scalar_2D_t), intent(in) :: self
+      integer, intent(in) :: direction
+      double precision, allocatable :: scalars(:,:)
     end function
 
     pure module function scalar_2D_values(self) result(values)
@@ -389,6 +351,14 @@ pure module function vector_2D_grid(self, component, coordinate) result(vector_g
       double precision, allocatable :: vector_grid_1D(:) !! grid points along the requested coordinate direction
     end function
 
+    pure module function vector_2D_values(self, direction) result(vector_values)
+      !! Result contains the vector values for the component designated by "direction"
+      implicit none
+      class(vector_2D_t), intent(in) :: self
+      integer, intent(in) :: direction
+      double precision, allocatable :: vector_values(:,:) 
+    end function
+
     pure module function divergence_2D_grid(self, direction) result(divergence_grid_1D)
       !! Result array contains divergence grid locations along the requested spatial direction
       implicit none
@@ -404,6 +374,13 @@ pure module function vector_2D_to_centers_extended(self) result(vectors)
       double precision, allocatable :: vectors(:,:,:)
     end function
 
+    pure module function vector_2D_consistent(self) result(self_consistent)
+      !! Assert internal self consistency, including sufficient accuracy for divergence operator
+      implicit none
+      class(vector_2D_t), intent(in) :: self
+      logical self_consistent
+    end function
+
     pure module function divergence_2D_values(self) result(divergences)
       !! Vector values getter
       implicit none

src/formal/vector_2D_s.F90

@@ -41,18 +41,8 @@
     self_consistent = .true.
   end procedure
 
-  module procedure vector_2D_conformable_vector
-    call_julienne_assert(vector_2D_consistent(self))
-    call_julienne_assert(vector_2D_consistent(vector_2D))
-    call_julienne_assert(self%tensor_2D_conformable(vector_2D))
-    conformable = .true.
-  end procedure
-
-  module procedure vector_2D_conformable_scalar
-    call_julienne_assert(vector_2D_consistent(self))
-    call_julienne_assert(scalar_2D_consistent(scalar_2D))
-    call_julienne_assert(self%tensor_2D_conformable(scalar_2D))
-    conformable = .true.
+  module procedure vector_2D_values
+    vector_values = self%points_(direction,1,1,1)%values_
   end procedure
 
   module procedure construct_2D_vector_from_function
@@ -183,15 +173,18 @@ pure function description(coordinate, component) result(point_cloud)
 
     call_julienne_assert(vector_2D%conformable(scalar_2D))
 
-    associate(vector => vector_2D%to_centers_extended())
-
-      call_julienne_assert(.all. ([size(vector,x_dir), size(vector,y_dir)] .equalsExpected. shape(scalar_2D%points(1,1,1,1)%values_)))
+    associate( &
+       scalar_x => scalar_2D%to_faces(x_dir) &
+      ,scalar_y => scalar_2D%to_faces(y_dir) &
+    )
+      call_julienne_assert(.all. (shape(scalar_x) .equalsExpected. shape(vector_2D%points_(x_dir,1,1,1)%values_)))
+      call_julienne_assert(.all. (shape(scalar_y) .equalsExpected. shape(vector_2D%points_(y_dir,1,1,1)%values_)))
 
-      vector_x_scalar = construct_2D_vector_from_components( &
+      vector_x_scalar = vector_2D_t( &
          tensor_2D_t( &
             points = reshape( &
-               source = [ points_2D_t(vector(:,:,x_dir) * scalar_2D%points_(1,1,1,1)%values_) &
-                         ,points_2D_t(vector(:,:,y_dir) * scalar_2D%points_(1,1,1,1)%values_)] &
+               source = [ points_2D_t(vector_2D%points_(x_dir,1,1,1)%values_ * scalar_x) &
+                         ,points_2D_t(vector_2D%points_(y_dir,1,1,1)%values_ * scalar_y)] &
               ,shape  = [space_dimension,1,1,1] &
             ) &
            ,cells = vector_2D%cells_ &

src/formal_m.f90

@@ -25,7 +25,8 @@ module formal_m
     ,gradient_2D_t & ! result of `.grad. s` for a scalar_2D_t s
     ,scalar_2D_initializer_i & ! scalar_2D_t initializer abstract interface
     ,vector_2D_initializer_i & ! vector_2D_t initializar abstract interface
-    ,divergence_2D_initializer_i ! divergence_2D_t initializar abstract interface
+    ,divergence_2D_initializer_i & ! divergence_2D_t initializar abstract interface
+    ,x_dir, y_dir
 
   use differential_operators_1D_m, only : &
      gradient_operator_1D_t & ! matrix operator defining a 1D mimetic gradient