Repository: pyccel/pyccel (language: python)

Issue:
Deprecate InhomogeneousTupleVariable
A tuple does not need to have its own type of Variable. This distinction makes it impossible to create tuples in classes. Inhomogeneous tuples can be handled using inline class definitions.
Ambiguous interface doesn't raise an error
### Describe the bug
When printing to Fortran there is no way to differentiate between arrays of the same type but different orders. This is a problem when building an `interface`. Currently the printer filters out one of the competing prototypes arbitrarily. As a result no error is raised, the result is just wrong if you call this function from another function

### To Reproduce
Provide code to reproduce the behaviour:
```python
import numpy as np
from pyccel.decorators import template

@template('T', ['int[:,:]', 'int[:,:](order=F)'])
def print_array(x : 'T'):
    print(x)

if __name__ == '__main__':
    a = np.array(((1,2,3), (4,5,6)))
    b = np.array(((1,2,3), (4,5,6)), order='F')
    print_array(a)
    print_array(b)
```

### Error details
Python output:
```
[[1 2 3]
 [4 5 6]]
[[1 2 3]
 [4 5 6]]
```
Pyccelised output:
```
[[1 2 3]
[4 5 6]]
[[1 4]
[2 5]
[3 6]]
```

### Expected behaviour
An error should be raised, either by Pyccel or by the compiler

### Language
Fortran
Containers: Error required for list of array
### Describe the bug
Containers such as lists and sets need to ensure that each element respects the homogeneity conditions. Each element should have the same:
- datatype
- precision
- rank
- order
- class type

However there is currently no way to determine the order of an element of a container. The type mechanism should be improved to handle this.

~~But in the meantime an error should be raised to prevent the user from creating lists of multi-dimensional arrays (2D or more) as we would not know how to index an element of such a list.~~

### To Reproduce
Provide code to reproduce the behaviour:
```python
a = [np.ones((3,4))]
a.append(np.ones((5,6), order='F'))
```

```python
a : 'float[:,:]' = []
a.append(np.ones((5,6), order='F'))
```

Strings
See old issue #338  
Ambiguous interface doesn't raise an error
### Describe the bug
When printing to Fortran there is no way to differentiate between arrays of the same type but different orders. This is a problem when building an `interface`. Currently the printer filters out one of the competing prototypes arbitrarily. As a result no error is raised, the result is just wrong if you call this function from another function

### To Reproduce
Provide code to reproduce the behaviour:
```python
import numpy as np
from pyccel.decorators import template

@template('T', ['int[:,:]', 'int[:,:](order=F)'])
def print_array(x : 'T'):
    print(x)

if __name__ == '__main__':
    a = np.array(((1,2,3), (4,5,6)))
    b = np.array(((1,2,3), (4,5,6)), order='F')
    print_array(a)
    print_array(b)
```

### Error details
Python output:
```
[[1 2 3]
 [4 5 6]]
[[1 2 3]
 [4 5 6]]
```
Pyccelised output:
```
[[1 2 3]
[4 5 6]]
[[1 4]
[2 5]
[3 6]]
```

### Expected behaviour
An error should be raised, either by Pyccel or by the compiler

### Language
Fortran


Discussion / hints:


cc: @Farouk-Echaref 



Fix patch:
diff --git a/.dict_custom.txt b/.dict_custom.txt
index 2e5398c733..eb2a5fb83f 100644
--- a/.dict_custom.txt
+++ b/.dict_custom.txt
@@ -113,3 +113,5 @@ setter
 bitwise
 datatyping
 datatypes
+indexable
+traceback
diff --git a/CHANGELOG.md b/CHANGELOG.md
index c2ae1c053d..99bbea22a9 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -14,6 +14,8 @@ All notable changes to this project will be documented in this file.
 -   #1750 : Add Python support for set method `remove()`.
 -   #1787 : Ensure `STC` is installed with Pyccel.
 -   #1743 : Add Python support for set method `discard()`.
+-   \[INTERNALS\] Added `container_rank` property to `ast.datatypes.PyccelType` objects.
+-   \[DEVELOPER\] Added an improved traceback to the developer-mode errors for errors in function calls.
 
 ### Fixed
 
@@ -29,6 +31,7 @@ All notable changes to this project will be documented in this file.
 -   #1785 : Add missing cast when creating an array of booleans from non-boolean values.
 -   #1218 : Fix bug when assigning an array to a slice in Fortran.
 -   #1830 : Fix missing allocation when returning an annotated array expression.
+-   #1821 : Ensure an error is raised when creating an ambiguous interface.
 
 ### Changed
 -   #1720 : functions with the `@inline` decorator are no longer exposed to Python in the shared library.
@@ -39,11 +42,15 @@ All notable changes to this project will be documented in this file.
 -   \[INTERNALS\] Build `utilities.metaclasses.ArgumentSingleton` on the fly to ensure correct docstrings.
 -   \[INTERNALS\] Rewrite datatyping system. See #1722.
 -   \[INTERNALS\] Moved precision from `ast.basic.TypedAstNode` to an internal property of `ast.datatypes.FixedSizeNumericType` objects.
+-   \[INTERNALS\] Moved rank from `ast.basic.TypedAstNode` to an internal property of `ast.datatypes.PyccelType` objects.
+-   \[INTERNALS\] Moved order from `ast.basic.TypedAstNode` to an internal property of `ast.datatypes.PyccelType` objects.
 -   \[INTERNALS\] Use cached `__add__` method to determine result type of arithmetic operations.
 -   \[INTERNALS\] Use cached `__and__` method to determine result type of bitwise comparison operations.
 -   \[INTERNALS\] Removed unused `fcode`, `ccode`, `cwrappercode`, `luacode`, and `pycode` functions from printers.
 -   \[INTERNALS\] Removed unused arguments from methods in `pyccel.codegen.codegen.Codegen`.
 -   \[INTERNALS\] Stop storing `FunctionDef`, `ClassDef`, and `Import` objects inside `CodeBlock`s.
+-   \[INTERNALS\] Remove the `order` argument from the `pyccel.ast.core.Allocate` constructor.
+-   \[INTERNALS\] Remove `rank` and `order` arguments from `pyccel.ast.variable.Variable` constructor.
 
 ### Deprecated
 
diff --git a/developer_docs/ast_nodes.md b/developer_docs/ast_nodes.md
index 6f1c3ea4db..785a006868 100644
--- a/developer_docs/ast_nodes.md
+++ b/developer_docs/ast_nodes.md
@@ -54,14 +54,6 @@ The order indicates how an array is laid out in memory. This can either be row-m
 
 The static type is the class type that would be assigned to an object created using an instance of this class as a type annotation.
 
-### Static rank
-
-The static rank is the rank that would be assigned to an object created using an instance of this class as a type annotation.
-
-### Static order
-
-The static order is the order that would be assigned to an object created using an instance of this class as a type annotation.
-
 ## Pyccel Internal Function
 
 The class `pyccel.ast.internals.PyccelInternalFunction` is a super class. This class should never be used directly but provides functionalities which are common to certain AST objects. These AST nodes are those which describe functions which are supported by Pyccel. For example it is used for functions from the `math` library, the `cmath` library, the `numpy` library, etc. `PyccelInternalFunction` inherits from `TypedAstNode`. The type information for the sub-class describes the type of the result of the function.
diff --git a/developer_docs/type_inference.md b/developer_docs/type_inference.md
index bee3e4af61..ad28052d42 100644
--- a/developer_docs/type_inference.md
+++ b/developer_docs/type_inference.md
@@ -68,7 +68,7 @@ The and operator describes what happens when two numeric types are combined in a
 When using these operators on an unknown number of arguments it can be useful to use `NativeGeneric()` as a starting point for the sum.
 
 #### Container Type
-A `ContainerType` is an object which is comprised of `FixedSizeType` objects (e.g. `ndarray`,`list`,`tuple`, custom class). The sub-class `HomogeneousContainerType` describes containers which contain homogeneous data. These objects are characterised by an `element_type`. The elements of a `HomogeneousContainerType` are instances of `PyccelType`, but they can be either `FixedSizeType`s or `ContainerType`s.
+A `ContainerType` is an object which is comprised of `FixedSizeType` objects (e.g. `ndarray`,`list`,`tuple`, custom class). The sub-class `HomogeneousContainerType` describes containers which contain homogeneous data. These objects are characterised by an `element_type`, a `rank` (and associated `container_rank`) and an `order`. The elements of a `HomogeneousContainerType` are instances of `PyccelType`, but they can be either `FixedSizeType`s or `ContainerType`s. The `container_rank` is an integer equal to the number of indices necessary to index the container and get an element of type `element_type`. As these elements may also be indexable the `rank` property allows us to get the number of indices necessary to obtain a scalar element. It is the sum of all the `container_rank`s in the nested types. The `order` specifies the order in which the indices should be used to index the object. This is discussed in detail in the [order docs](./order_docs.md).
 
 `HomogeneousContainerType`s also contain some utility functions. They implement `primitive_type` and `precision` to get the properties of the internal `FixedSizeType` (even if that type is inside another `HomogeneousContainerType`). They also implement `switch_basic_type` which creates a new `HomogeneousContainerType` which is similar to the current `HomogeneousContainerType`. The only difference is that the `FixedSizeType` is exchanged. This is useful when we want to preserve information about the container but need to change the type. For example, when we divide an integer by another we get a floating point type. When we divide a NumPy array or a CuPy array of integers by an integer (or array of integers) we get a NumPy/CuPy array of floating point numbers (with default Python precision). In order to preserve the container type we therefore call `switch_basic_type`. So for the division in the case of NumPy arrays, we want to change the type from `np.ndarray[int]` to `np.ndarray[float]`. This is done in one line:
 ```python
@@ -89,4 +89,8 @@ for container in new_container_types:
 
 The `switch_basic_type` cannot be implemented generally in `PyccelType` as there is no logical interpretation for an inhomogeneous `ContainerType`, however the function is also implemented (as the identity function) for `FixedSizeType`s so `switch_basic_type` can be used without the need for type checks (generally inhomogeneous containers will not be valid arguments to classes which may need to use the `switch_basic_type` function).
 
+`HomogeneousContainerType`s also contain a `switch_rank` function. This function is similar to `switch_basic_type` in that it is used to obtain a type which is similar in all but one characteristic. It is usually used to reduce the rank of an object, for example when calculating the type of a slice, however in the future it can also be used to increase the size of the type (e.g. to implement `np.newaxis`), in this case an order may need to be provided to add additional context. Increasing the rank is only possible for multi-dimensional types (e.g. `NumpyNDArrayType`) however the rank can be decreased for any `ContainerType`. If the rank is reduced by more than the `container_rank`, this function is called recursively.
+
+For multi-dimensional `HomogeneousContainerType`s (e.g. `NumpyNDArrayType`) the function `swap_order` is also implemented. This inverts the ordering, changing from 'C' to 'F' or 'F' to 'C' if the rank is greater than 1.
+
 In order to access the internal `FixedSizeType`, `PyccelType` also implements a `datatype` property. This makes more sense in the case of a `HomogeneousContainerType` however it is also implemented (as the identity function) for `FixedSizeType`s so the low-level type can be obtained without the need for type checks.
diff --git a/docs/ndarrays.md b/docs/ndarrays.md
index 08f2dc15fa..3e392673eb 100644
--- a/docs/ndarrays.md
+++ b/docs/ndarrays.md
@@ -22,9 +22,10 @@ Generally a variable in Pyccel should always keep its initial type, this also tr
 ```Python
 import numpy as np
 
-a = np.array([1, 2, 3], dtype=float)
-#(some code...)
-a = np.array([1, 2, 3], dtype=int)
+if __name__ == '__main__':
+    a = np.array([1, 2, 3], dtype=float)
+    #(some code...)
+    a = np.array([1, 2, 3], dtype=int)
 ```
 
 _OUTPUT_ :
@@ -46,9 +47,10 @@ Pyccel calls its own garbage collector when needed, but has a set of rules to do
     ```Python
     import numpy as np
 
-    a = np.ones((10, 20))
-    #(some code...)
-    a = np.ones(10)
+    if __name__ == '__main__':
+        a = np.ones((10, 20))
+        #(some code...)
+        a = np.ones(10)
     ```
 
     _OUTPUT_ :
@@ -66,9 +68,10 @@ This limitation is due to the fact that the rank of Fortran allocatable objects
     ```Python
     import numpy as np
 
-    a = np.array([1, 2, 3, 4, 5])
-    b = np.array([1, 2, 3])
-    a = b
+    if __name__ == '__main__':
+        a = np.array([1, 2, 3, 4, 5])
+        b = np.array([1, 2, 3])
+        a = b
     ```
 
     _OUTPUT_ :
@@ -139,10 +142,11 @@ This limitation is due to the fact that the rank of Fortran allocatable objects
     ```Python
     import numpy as np
 
-    a = np.ones(10)
-    b = a[:5]
-    #(some code...)
-    a = np.zeros(20)
+    if __name__ == '__main__':
+        a = np.ones(10)
+        b = a[:5]
+        #(some code...)
+        a = np.zeros(20)
     ```
 
     _OUTPUT_ :
@@ -157,7 +161,7 @@ This limitation is set since we need to free the previous data when we reallocat
 
 ### Slicing and indexing ###
 
-The indexing and slicing in Pyccel handles only the basic indexing of [numpy arrays](https://numpy.org/doc/stable/user/basics.indexing.html).
+The indexing and slicing in Pyccel handles only the basic indexing of [numpy arrays](https://numpy.org/doc/stable/user/basics.indexing.html). When multiple indexing expressions are used on the same variable Pyccel squashes them into one object. This means that we do not handle multiple slice indices applied to the same variable (e.g. `a[1::2][2:]`). This is not recommended anyway as it makes code hard to read.
 
 Some examples:
 
@@ -166,8 +170,9 @@ Some examples:
     ```Python
     import numpy as np
 
-    a = np.array([1, 3, 4, 5])
-    a[0] = 0
+    if __name__ == '__main__':
+        a = np.array([1, 3, 4, 5])
+        a[0] = 0
     ```
 
     -   C equivalent:
@@ -211,8 +216,9 @@ Some examples:
     ```Python
     import numpy as np
 
-    a = np.ones((10, 20))
-    b = a[2:, :5]
+    if __name__ == '__main__':
+        a = np.ones((10, 20))
+        b = a[2:, :5]
     ```
 
     -   C equivalent:
@@ -257,10 +263,11 @@ Some examples:
     ```Python
     import numpy as np
 
-    a = np.array([[1, 2, 3, 4], [5, 6, 7, 8]])
-    b = a[1]
-    c = b[2]
-    print(c)
+    if __name__ == '__main__':
+        a = np.array([[1, 2, 3, 4], [5, 6, 7, 8]])
+        b = a[1]
+        c = b[2]
+        print(c)
     ```
 
     -   C equivalent:
@@ -311,6 +318,65 @@ Some examples:
         end program prog_ex
         ```
 
+-   Python code:
+
+    ```Python
+    import numpy as np
+
+    if __name__ == '__main__':
+        a = np.array([1, 2, 3, 4, 5, 6, 7, 8])
+        b = a[1::2][2]
+        print(b)
+    ```
+
+    -   C equivalent:
+
+        ```C
+        #include <stdlib.h>
+        #include "ndarrays.h"
+        #include <stdint.h>
+        #include <string.h>
+        #include <stdio.h>
+        #include <inttypes.h>
+        int main()
+        {
+            t_ndarray a = {.shape = NULL};
+            int64_t b;
+            a = array_create(1, (int64_t[]){INT64_C(8)}, nd_int64, false, order_c);
+            int64_t Dummy_0000[] = {INT64_C(1), INT64_C(2), INT64_C(3), INT64_C(4), INT64_C(5), INT64_C(6), INT64_C(7), INT64_C(8)};
+            memcpy(&a.nd_int64[INT64_C(0)], Dummy_0000, 8 * a.type_size);
+            b = GET_ELEMENT(a, nd_int64, INT64_C(5));
+            printf("%"PRId64"\n", b);
+            free_array(&a);
+            return 0;
+        }
+        ```
+
+    -   Fortran equivalent:
+
+        ```Fortran
+        program prog_prog_tmp_index
+
+          use tmp_index
+
+          use, intrinsic :: ISO_C_Binding, only : i64 => C_INT64_T
+          use, intrinsic :: ISO_FORTRAN_ENV, only : stdout => output_unit
+          implicit none
+
+          integer(i64), allocatable :: a(:)
+          integer(i64) :: b
+
+          allocate(a(0:7_i64))
+          a = [1_i64, 2_i64, 3_i64, 4_i64, 5_i64, 6_i64, 7_i64, 8_i64]
+          b = a(5_i64)
+          write(stdout, '(I0)', advance="yes") b
+          if (allocated(a)) then
+            deallocate(a)
+          end if
+
+        end program prog_prog_tmp_index
+        ```
+
 ## NumPy [ndarray](https://numpy.org/doc/stable/reference/generated/numpy.ndarray.html) functions/properties progress in Pyccel ##
 
 -   Supported [types](https://numpy.org/devdocs/user/basics.types.html):
diff --git a/pyccel/ast/basic.py b/pyccel/ast/basic.py
index b105a55c7a..63ecd80b66 100644
--- a/pyccel/ast/basic.py
+++ b/pyccel/ast/basic.py
@@ -534,7 +534,7 @@ def rank(self):
         Number of dimensions of the object. If the object is a scalar then
         this is equal to 0.
         """
-        return self._rank # pylint: disable=no-member
+        return self.class_type.rank
 
     @property
     def dtype(self):
@@ -557,7 +557,7 @@ def order(self):
         ('F') format. This is only relevant if rank > 1. When it is not relevant
         this function returns None.
         """
-        return self._order # pylint: disable=no-member
+        return self.class_type.order
 
     @property
     def class_type(self):
@@ -570,33 +570,6 @@ def class_type(self):
         """
         return self._class_type # pylint: disable=no-member
 
-    @classmethod
-    def static_rank(cls):
-        """
-        Number of dimensions of the object.
-
-        Number of dimensions of the object. If the object is a scalar then
-        this is equal to 0.
-
-        This function is static and will return an AttributeError if the
-        class does not have a predetermined rank.
-        """
-        return cls._rank # pylint: disable=no-member
-
-    @classmethod
-    def static_order(cls):
-        """
-        The data layout ordering in memory.
-
-        Indicates whether the data is stored in row-major ('C') or column-major
-        ('F') format. This is only relevant if rank > 1. When it is not relevant
-        this function returns None.
-
-        This function is static and will return an AttributeError if the
-        class does not have a predetermined order.
-        """
-        return cls._order # pylint: disable=no-member
-
     @classmethod
     def static_type(cls):
         """
@@ -625,8 +598,6 @@ def copy_attributes(self, x):
             The node from which the attributes should be copied.
         """
         self._shape      = x.shape
-        self._rank       = x.rank
-        self._order      = x.order
         self._class_type = x.class_type
 
 
diff --git a/pyccel/ast/bind_c.py b/pyccel/ast/bind_c.py
index d1af786bab..cb6a4bd09b 100644
--- a/pyccel/ast/bind_c.py
+++ b/pyccel/ast/bind_c.py
@@ -341,7 +341,7 @@ def __init__(self, var, original_res_var, scope, **kwargs):
         name = original_res_var.name
         self._shape   = [scope.get_temporary_variable(PythonNativeInt(),
                             name=f'{name}_shape_{i+1}')
-                         for i in range(original_res_var._rank)]
+                         for i in range(original_res_var.rank)]
         self._original_res_var = original_res_var
         super().__init__(var, **kwargs)
 
diff --git a/pyccel/ast/bitwise_operators.py b/pyccel/ast/bitwise_operators.py
index a2ba0e3ad6..5b68c417fd 100644
--- a/pyccel/ast/bitwise_operators.py
+++ b/pyccel/ast/bitwise_operators.py
@@ -86,8 +86,6 @@ class PyccelBitOperator(PyccelOperator):
         The second argument passed to the operator.
     """
     _shape = None
-    _rank  = 0
-    _order = None
     __slots__ = ('_class_type',)
 
     def __init__(self, arg1, arg2):
@@ -130,12 +128,12 @@ def _calculate_type(self, arg1, arg2):
 
         return class_type
 
-    def _set_shape_rank(self):
+    def _set_shape(self):
         """
-        Set the shape and rank of the resulting object.
+        Set the shape of the resulting object.
 
-        Set the shape and rank of the resulting object. For a PyccelBitOperator,
-        the shape and rank are class attributes so nothing needs to be done.
+        Set the shape of the resulting object. For a PyccelBitOperator,
+        the shape is a class attribute so nothing needs to be done.
         """
 
     def __repr__(self):
diff --git a/pyccel/ast/builtin_methods/list_methods.py b/pyccel/ast/builtin_methods/list_methods.py
index 46ac641f1d..6e59e39c5e 100644
--- a/pyccel/ast/builtin_methods/list_methods.py
+++ b/pyccel/ast/builtin_methods/list_methods.py
@@ -78,20 +78,13 @@ class ListAppend(ListMethod):
         The argument passed to append() method.
     """
     __slots__ = ()
-    _dtype = VoidType()
     _shape = None
-    _order = None
-    _rank = 0
     _class_type = VoidType()
     name = 'append'
 
     def __init__(self, list_obj, new_elem) -> None:
         expected_type = list_obj.class_type.element_type
-        is_homogeneous = (
-            new_elem.class_type == expected_type and
-            list_obj.rank - 1 == new_elem.rank
-        )
-        if not is_homogeneous:
+        if new_elem.class_type != expected_type:
             raise TypeError(f"Expecting an argument of the same type as the elements of the list ({expected_type}) but received {new_elem.class_type}")
         super().__init__(list_obj, new_elem)
 
@@ -115,13 +108,11 @@ class ListPop(ListMethod) :
     index_element : TypedAstNode
         The current index value for the element to be popped.
     """
-    __slots__ = ('_class_type', '_rank', '_shape', '_order')
+    __slots__ = ('_class_type', '_shape')
     name = 'pop'
 
     def __init__(self, list_obj, index_element=None) -> None:
-        self._rank = list_obj.rank - 1
         self._shape = (None if len(list_obj.shape) == 1 else tuple(list_obj.shape[1:]))
-        self._order = (None if self._shape is None or len(self._shape) == 1 else list_obj.order)
         self._class_type = list_obj.class_type.element_type
         super().__init__(list_obj, index_element)
 
@@ -145,8 +136,6 @@ class ListClear(ListMethod) :
         The list object which the method is called from.
     """
     __slots__ = ()
-    _rank = 0
-    _order = None
     _shape = None
     _class_type = VoidType()
     name = 'clear'
@@ -182,18 +171,11 @@ class ListInsert(ListMethod):
     """
     __slots__ = ()
     _shape = None
-    _order = None
-    _rank = 0
     _class_type = VoidType()
     name = 'insert'
 
     def __init__(self, list_obj, index, new_elem) -> None:
-        expected_type = list_obj.class_type.element_type
-        is_homogeneous = (
-            new_elem.class_type == expected_type and
-            list_obj.rank - 1 == new_elem.rank
-        )
-        if not is_homogeneous:
+        if new_elem.class_type != list_obj.class_type.element_type:
             raise TypeError("Expecting an argument of the same type as the elements of the list")
         super().__init__(list_obj, index, new_elem)
 
@@ -256,18 +238,11 @@ class ListRemove(ListMethod) :
     """
     __slots__ = ()
     _shape = None
-    _order = None
-    _rank = 0
     _class_type = VoidType()
     name = 'remove'
 
     def __init__(self, list_obj, removed_obj) -> None:
-        expected_type = list_obj.class_type.element_type
-        is_homogeneous = (
-            removed_obj.class_type == expected_type and
-            list_obj.rank - 1 == removed_obj.rank
-        )
-        if not is_homogeneous:
+        if removed_obj.class_type != list_obj.class_type.element_type:
             raise TypeError(f"Can't remove an element of type {removed_obj.class_type} from {list_obj.class_type}")
         super().__init__(list_obj, removed_obj)
 
@@ -297,12 +272,10 @@ class ListCopy(ListMethod) :
     list_obj : TypedAstNode
         The list object which the method is called from.
     """
-    __slots__ = ('_class_type', '_rank', '_shape', '_order')
+    __slots__ = ('_class_type', '_shape')
     name = 'copy'
 
     def __init__(self, list_obj) -> None:
-        self._rank = list_obj.rank
         self._shape = list_obj.shape
-        self._order = list_obj.order
         self._class_type = list_obj.class_type
         super().__init__(list_obj)
diff --git a/pyccel/ast/builtin_methods/set_methods.py b/pyccel/ast/builtin_methods/set_methods.py
index feba49d7ff..7505901e5b 100644
--- a/pyccel/ast/builtin_methods/set_methods.py
+++ b/pyccel/ast/builtin_methods/set_methods.py
@@ -72,17 +72,11 @@ class SetAdd(SetMethod) :
     """
     __slots__ = ()
     _shape = None
-    _order = None
-    _rank = 0
     _class_type = VoidType()
     name = 'add'
 
     def __init__(self, set_variable, new_elem) -> None:
-        is_homogeneous = (
-            set_variable.class_type.element_type == new_elem.class_type and
-            set_variable.rank - 1 == new_elem.rank
-        )
-        if not is_homogeneous:
+        if set_variable.class_type.element_type != new_elem.class_type:
             raise TypeError("Expecting an argument of the same type as the elements of the set")
         super().__init__(set_variable, new_elem)
 
@@ -101,8 +95,6 @@ class SetClear(SetMethod):
     """
     __slots__ = ()
     _shape = None
-    _order = None
-    _rank = 0
     _class_type = VoidType()
     name = 'clear'
 
@@ -122,13 +114,11 @@ class SetCopy(SetMethod):
     set_variable : TypedAstNode
         The set on which the method will operate.
     """
-    __slots__ = ("_shape", "_order", "_rank", "_class_type",)
+    __slots__ = ("_shape", "_class_type",)
     name = 'copy'
 
     def __init__(self, set_variable):
         self._shape = set_variable._shape
-        self._order = set_variable._order
-        self._rank = set_variable._rank
         self._class_type = set_variable._class_type
         super().__init__(set_variable)
 
@@ -149,8 +139,6 @@ class SetPop(SetMethod):
         The name of the set.
     """
     __slots__ = ('_class_type',)
-    _rank = 0
-    _order = None
     _shape = None
     name = 'pop'
 
@@ -176,20 +164,13 @@ class SetRemove(SetMethod):
     """
     __slots__ = ()
     _shape = None
-    _order = None
-    _rank = 0
     _class_type = VoidType()
     name = 'remove'
 
     def __init__(self, set_variable, item) -> None:
         if not isinstance(item, TypedAstNode):
             raise TypeError(f"It is not possible to look for a {type(item).__name__} object in a set of {set_variable.dtype}")
-        expected_type = set_variable.class_type.element_type
-        is_homogeneous = (
-            expected_type == item.class_type and
-            set_variable.rank - 1 == item.rank
-        )
-        if not is_homogeneous:
+        if item.class_type != set_variable.class_type.element_type:
             raise TypeError(f"Can't remove an element of type {item.dtype} from a set of {set_variable.dtype}")
         super().__init__(set_variable, item)
 
@@ -212,17 +193,10 @@ class SetDiscard(SetMethod):
     """
     __slots__ = ()
     _shape = None
-    _order = None
-    _rank = 0
     _class_type = VoidType()
     name = 'discard'
 
     def __init__(self, set_variable, item) -> None:
-        expected_type = set_variable.class_type.element_type
-        is_homogeneous = (
-            expected_type == item.class_type and
-            set_variable.rank - 1 == item.rank
-        )
-        if not is_homogeneous:
+        if set_variable.class_type.element_type != item.class_type:
             raise TypeError("Expecting an argument of the same type as the elements of the set")
         super().__init__(set_variable, item)
diff --git a/pyccel/ast/builtins.py b/pyccel/ast/builtins.py
index e8399de319..0dcb5bb903 100644
--- a/pyccel/ast/builtins.py
+++ b/pyccel/ast/builtins.py
@@ -77,9 +77,7 @@ class PythonComplexProperty(PyccelInternalFunction):
         The object which the property is called from.
     """
     __slots__ = ()
-    _rank  = 0
     _shape = None
-    _order = None
     _class_type = PythonNativeFloat()
 
     def __init__(self, arg):
@@ -165,9 +163,7 @@ class PythonConjugate(PyccelInternalFunction):
         conjugate function.
     """
     __slots__ = ()
-    _rank  = 0
     _shape = None
-    _order = None
     _class_type = PythonNativeComplex()
     name = 'conjugate'
 
@@ -206,9 +202,7 @@ class PythonBool(PyccelInternalFunction):
     __slots__ = ()
     name = 'bool'
     _static_type = PythonNativeBool()
-    _rank  = 0
     _shape = None
-    _order = None
     _class_type = PythonNativeBool()
 
     def __new__(cls, arg):
@@ -251,9 +245,7 @@ class PythonComplex(PyccelInternalFunction):
     name = 'complex'
 
     _static_type = PythonNativeComplex()
-    _rank  = 0
     _shape = None
-    _order = None
     _class_type = PythonNativeComplex()
     _real_cast = PythonReal
     _imag_cast = PythonImag
@@ -413,9 +405,7 @@ class PythonFloat(PyccelInternalFunction):
     __slots__ = ()
     name = 'float'
     _static_type = PythonNativeFloat()
-    _rank  = 0
     _shape = None
-    _order = None
     _class_type = PythonNativeFloat()
 
     def __new__(cls, arg):
@@ -457,9 +447,7 @@ class PythonInt(PyccelInternalFunction):
     __slots__ = ()
     name = 'int'
     _static_type = PythonNativeInt()
-    _rank  = 0
     _shape = None
-    _order = None
     _class_type = PythonNativeInt()
 
     def __new__(cls, arg):
@@ -493,9 +481,8 @@ class PythonTuple(TypedAstNode):
     *args : tuple of TypedAstNode
         The arguments passed to the tuple function.
     """
-    __slots__ = ('_args','_is_homogeneous',
-            '_rank','_shape','_order', '_class_type')
-    _iterable        = True
+    __slots__ = ('_args','_is_homogeneous', '_shape', '_class_type')
+    _iterable = True
     _attribute_nodes = ('_args',)
 
     def __init__(self, *args):
@@ -505,9 +492,7 @@ def __init__(self, *args):
             return
         elif len(args) == 0:
             self._class_type = HomogeneousTupleType(GenericType())
-            self._rank  = 0
-            self._shape = None
-            self._order = None
+            self._shape = (LiteralInteger(0),)
             self._is_homogeneous = False
             return
 
@@ -534,7 +519,6 @@ def __init__(self, *args):
         if is_homogeneous:
             inner_shape = [() if a.rank == 0 else a.shape for a in args]
             self._shape = (LiteralInteger(len(args)), ) + inner_shape[0]
-            self._rank  = len(self._shape)
 
             if contains_pointers:
                 self._class_type = InhomogeneousTupleType(*[a.class_type for a in args])
@@ -542,12 +526,9 @@ def __init__(self, *args):
                 self._class_type = HomogeneousTupleType(args[0].class_type)
 
         else:
-            self._rank       = 1 if not is_homogeneous else max(a.rank for a in args) + 1
             self._class_type = InhomogeneousTupleType(*[a.class_type for a in args])
             self._shape      = (LiteralInteger(len(args)), )
 
-        self._order = None if self._rank < 2 else 'C'
-
     def __getitem__(self,i):
         def is_int(a):
             return isinstance(a, (int, LiteralInteger)) or \
@@ -669,7 +650,7 @@ class PythonList(TypedAstNode):
     FunctionalFor
         The `[]` function when it describes a comprehension.
     """
-    __slots__ = ('_args','_rank','_shape','_order', '_class_type')
+    __slots__ = ('_args', '_shape', '_class_type')
     _attribute_nodes = ('_args',)
 
     def __init__(self, *args):
@@ -678,30 +659,23 @@ def __init__(self, *args):
         if pyccel_stage == 'syntactic':
             return
         elif len(args) == 0:
-            self._rank  = 0
-            self._shape = None
-            self._order = None
+            self._shape = (LiteralInteger(0),)
             self._class_type = HomogeneousListType(GenericType())
             return
         arg0 = args[0]
         is_homogeneous = arg0.class_type is not GenericType() and \
                          all(a.class_type is not GenericType() and \
-                             arg0.class_type == a.class_type and \
-                             arg0.rank  == a.rank  and \
-                             arg0.order == a.order for a in args[1:])
+                             arg0.class_type == a.class_type for a in args[1:])
         if is_homogeneous:
             dtype = arg0.class_type
 
             inner_shape = [() if a.rank == 0 else a.shape for a in args]
-            self._rank = max(a.rank for a in args) + 1
             self._shape = (LiteralInteger(len(args)), ) + inner_shape[0]
-            self._rank  = len(self._shape)
 
         else:
             raise TypeError("Can't create an inhomogeneous list")
 
         self._class_type = HomogeneousListType(dtype)
-        self._order = None if self._rank < 2 else 'C'
 
     def __iter__(self):
         return self._args.__iter__()
@@ -747,7 +721,7 @@ class PythonSet(TypedAstNode):
     *args : tuple of TypedAstNodes
         The arguments passed to the operator.
     """
-    __slots__ = ('_args','_class_type','_rank','_shape','_order')
+    __slots__ = ('_args','_class_type','_shape')
     _attribute_nodes = ('_args',)
 
     def __init__(self, *args):
@@ -758,21 +732,17 @@ def __init__(self, *args):
         arg0 = args[0]
         is_homogeneous = arg0.class_type is not GenericType() and \
                          all(a.class_type is not GenericType() and \
-                             arg0.class_type == a.class_type and \
-                             arg0.rank  == a.rank  and \
-                             arg0.order == a.order for a in args[1:])
+                             arg0.class_type == a.class_type for a in args[1:])
         if is_homogeneous:
-            dtype = arg0.dtype
+            elem_type = arg0.class_type
             inner_shape = [() if a.rank == 0 else a.shape for a in args]
             self._shape = (LiteralInteger(len(args)), ) + inner_shape[0]
-            self._rank  = len(self._shape)
-            if self._rank > 1:
+            if elem_type.rank > 0:
                 raise TypeError("Pyccel can't hash non-scalar types")
         else:
             raise TypeError("Can't create an inhomogeneous set")
 
-        self._class_type = HomogeneousSetType(dtype)
-        self._order = None if self._rank < 2 else 'C'
+        self._class_type = HomogeneousSetType(elem_type)
 
     def __iter__(self):
         return self._args.__iter__()
@@ -998,12 +968,10 @@ class PythonAbs(PyccelInternalFunction):
     x : TypedAstNode
         The argument passed to the function.
     """
-    __slots__ = ('_rank','_shape','_order','_class_type')
+    __slots__ = ('_shape','_class_type')
     name = 'abs'
     def __init__(self, x):
         self._shape     = x.shape
-        self._rank      = x.rank
-        self._order     = x.order
         self._class_type = PythonNativeInt() if x.dtype is PythonNativeInt() else PythonNativeFloat()
         super().__init__(x)
 
@@ -1030,9 +998,7 @@ class PythonSum(PyccelInternalFunction):
     """
     __slots__ = ('_class_type',)
     name   = 'sum'
-    _rank  = 0
     _shape = None
-    _order = None
 
     def __init__(self, arg):
         if not isinstance(arg, TypedAstNode):
@@ -1066,9 +1032,7 @@ class PythonMax(PyccelInternalFunction):
     """
     __slots__ = ('_class_type',)
     name   = 'max'
-    _rank  = 0
     _shape = None
-    _order = None
 
     def __init__(self, *x):
         if len(x)==1:
@@ -1104,9 +1068,8 @@ class PythonMin(PyccelInternalFunction):
     """
     __slots__ = ('_class_type',)
     name   = 'min'
-    _rank  = 0
     _shape = None
-    _order = None
+
     def __init__(self, *x):
         if len(x)==1:
             x = x[0]
diff --git a/pyccel/ast/c_concepts.py b/pyccel/ast/c_concepts.py
index abd684742b..2788e51495 100644
--- a/pyccel/ast/c_concepts.py
+++ b/pyccel/ast/c_concepts.py
@@ -47,6 +47,8 @@ class CStackArray(HomogeneousContainerType, metaclass=ArgumentSingleton):
     """
     __slots__ = ('_element_type',)
     _name = 'c_stackarray'
+    _container_rank = 1
+    _order = None
 
     def __init__(self, element_type):
         assert isinstance(element_type, FixedSizeType)
@@ -76,16 +78,14 @@ class ObjectAddress(TypedAstNode):
     'a'
     """
 
-    __slots__ = ('_obj', '_rank', '_shape', '_order', '_class_type')
+    __slots__ = ('_obj', '_shape', '_class_type')
     _attribute_nodes = ('_obj',)
 
     def __init__(self, obj):
         if not isinstance(obj, TypedAstNode):
             raise TypeError("object must be an instance of TypedAstNode")
         self._obj        = obj
-        self._rank       = obj.rank
         self._shape      = obj.shape
-        self._order      = obj.order
         self._class_type = obj.class_type
         super().__init__()
 
@@ -121,8 +121,7 @@ class PointerCast(TypedAstNode):
     cast_type : TypedAstNode
         A TypedAstNode describing the object resulting from the cast.
     """
-    __slots__ = ('_obj', '_rank', '_shape', '_order',
-            '_class_type', '_cast_type')
+    __slots__ = ('_obj', '_shape', '_class_type', '_cast_type')
     _attribute_nodes = ('_obj',)
 
     def __init__(self, obj, cast_type):
@@ -130,10 +129,8 @@ def __init__(self, obj, cast_type):
             raise TypeError("object must be an instance of TypedAstNode")
         assert getattr(obj, 'is_alias', False)
         self._obj        = obj
-        self._rank       = cast_type.rank
         self._shape      = cast_type.shape
         self._class_type = cast_type.class_type
-        self._order      = cast_type.order
         self._cast_type  = cast_type
         super().__init__()
 
diff --git a/pyccel/ast/cmathext.py b/pyccel/ast/cmathext.py
index af1901d63d..5f821dfeea 100644
--- a/pyccel/ast/cmathext.py
+++ b/pyccel/ast/cmathext.py
@@ -10,7 +10,7 @@
 from pyccel.ast.builtins  import PythonReal, PythonImag
 from pyccel.ast.core      import PyccelFunctionDef, Module
 from pyccel.ast.datatypes import PythonNativeBool, PythonNativeFloat, PythonNativeComplex
-from pyccel.ast.datatypes import PrimitiveComplexType
+from pyccel.ast.datatypes import PrimitiveComplexType, HomogeneousTupleType
 from pyccel.ast.internals import PyccelInternalFunction
 from pyccel.ast.literals  import LiteralInteger
 from pyccel.ast.operators import PyccelAnd, PyccelOr
@@ -75,8 +75,6 @@ class CmathFunctionComplex(MathFunctionBase):
     """
     __slots__ = ()
     _shape = None
-    _rank  = 0
-    _order = None
     _class_type = PythonNativeComplex()
 
     def __init__(self, z : 'TypedAstNode'):
@@ -436,9 +434,7 @@ class CmathPolar(PyccelInternalFunction):
     __slots__ = ()
     name = 'polar'
     _shape = (LiteralInteger(2),)
-    _rank  = 1
-    _order = None
-    _class_type = PythonNativeFloat()
+    _class_type = HomogeneousTupleType(PythonNativeFloat())
 
     def __init__(self, z):
         super().__init__(z)
@@ -459,8 +455,6 @@ class CmathRect(PyccelInternalFunction):
     __slots__ = ()
     name = 'rect'
     _shape = None
-    _rank  = 0
-    _order = None
     _class_type = PythonNativeComplex()
     def __init__(self, r, phi):
         super().__init__(r, phi)
diff --git a/pyccel/ast/core.py b/pyccel/ast/core.py
index 54a36ec8f6..a5e959bec0 100644
--- a/pyccel/ast/core.py
+++ b/pyccel/ast/core.py
@@ -198,13 +198,11 @@ class Duplicate(TypedAstNode):
     length : TypedAstNode
         The number of times the val should appear in the final object.
     """
-    __slots__ = ('_val', '_length','_rank','_shape','_order','_class_type')
+    __slots__ = ('_val', '_length','_shape','_class_type')
     _attribute_nodes = ('_val', '_length')
 
     def __init__(self, val, length):
-        self._rank       = val.rank
         self._shape      = tuple(s if i!= 0 else PyccelMul(s, length, simplify=True) for i,s in enumerate(val.shape))
-        self._order      = val.order
         self._class_type = val.class_type
 
         self._val       = val
@@ -238,14 +236,12 @@ class Concatenate(TypedAstNode):
     arg2 : TypedAstNodes
            The second tuple.
     """
-    __slots__ = ('_args','_rank','_shape','_order','_class_type')
+    __slots__ = ('_args','_shape','_class_type')
     _attribute_nodes = ('_args',)
 
     def __init__(self, arg1, arg2):
-        self._rank       = arg1.rank
         shape_addition   = arg2.shape[0]
         self._shape      = tuple(s if i!= 0 else PyccelAdd(s, shape_addition) for i,s in enumerate(arg1.shape))
-        self._order      = arg1.order
         self._class_type = arg1.class_type
 
         self._args = (arg1, arg2)
@@ -388,10 +384,6 @@ class Allocate(PyccelAstNode):
     shape : int or iterable or None
         Shape of the array after allocation (None for scalars).
 
-    order : str {'C'|'F'}
-        Ordering of multi-dimensional array after allocation
-        ('C' = row-major, 'F' = column-major).
-
     status : str {'allocated'|'unallocated'|'unknown'}
         Variable allocation status at object creation.
 
@@ -409,7 +401,7 @@ class Allocate(PyccelAstNode):
     _attribute_nodes = ('_variable', '_like')
 
     # ...
-    def __init__(self, variable, *, shape, order, status, like = None):
+    def __init__(self, variable, *, shape, status, like = None):
 
         if not isinstance(variable, (Variable, PointerCast)):
             raise TypeError(f"Can only allocate a 'Variable' object, got {type(variable)} instead")
@@ -421,13 +413,10 @@ def __init__(self, variable, *, shape, order, status, like = None):
         if shape and not isinstance(shape, (int, tuple, list)):
             raise TypeError(f"Cannot understand 'shape' parameter of type '{type(shape)}'")
 
-        if variable.rank != len(shape):
+        class_type = variable.class_type
+        if class_type.rank != len(shape):
             raise ValueError("Incompatible rank in variable allocation")
 
-        # rank is None for lambda functions
-        if variable.rank is not None and variable.rank > 1 and variable.order != order:
-            raise ValueError("Incompatible order in variable allocation")
-
         if not isinstance(status, str):
             raise TypeError(f"Cannot understand 'status' parameter of type '{type(status)}'")
 
@@ -436,7 +425,7 @@ def __init__(self, variable, *, shape, order, status, like = None):
 
         self._variable = variable
         self._shape    = shape
-        self._order    = order
+        self._order    = variable.order
         self._status   = status
         self._like = like
         super().__init__()
@@ -1950,7 +1939,7 @@ class FunctionCall(TypedAstNode):
         The function where the call takes place.
     """
     __slots__ = ('_arguments','_funcdef','_interface','_func_name','_interface_name',
-                 '_shape','_rank','_order','_class_type')
+                 '_shape','_class_type')
     _attribute_nodes = ('_arguments','_funcdef','_interface')
 
     def __init__(self, func, args, current_function=None):
@@ -2024,14 +2013,10 @@ def __init__(self, func, args, current_function=None):
         self._func_name  = func.name
         n_results = len(func.results)
         if n_results == 1:
-            self._rank       = func.results[0].var.rank
             self._shape      = func.results[0].var.shape
-            self._order      = func.results[0].var.order
             self._class_type = func.results[0].var.class_type
         elif n_results == 0:
-            self._rank       = 0
             self._shape      = None
-            self._order      = None
             self._class_type = VoidType()
         else:
             dtypes = [r.var.dtype for r in func.results]
@@ -2039,9 +2024,7 @@ def __init__(self, func, args, current_function=None):
                 dtype = HomogeneousTupleType(dtypes[0])
             else:
                 dtype = InhomogeneousTupleType(*dtypes)
-            self._rank       = 1
             self._shape      = (LiteralInteger(n_results),)
-            self._order      = None
             self._class_type = dtype
 
         super().__init__()
@@ -4463,6 +4446,7 @@ def args_var(self):
 
 # ...
 
+
 class Decorator(PyccelAstNode):
     """ Class representing a function decorator.
     For now this is just designed to handle the pyccel decorators
diff --git a/pyccel/ast/cwrapper.py b/pyccel/ast/cwrapper.py
index bc796550e6..6104719597 100644
--- a/pyccel/ast/cwrapper.py
+++ b/pyccel/ast/cwrapper.py
@@ -267,10 +267,8 @@ class PyBuildValueNode(PyccelInternalFunction):
     """
     __slots__ = ('_flags','_result_args')
     _attribute_nodes = ('_result_args',)
-    _rank = 0
     _shape = ()
     _class_type = PyccelPyObject()
-    _order = None
 
     def __init__(self, result_args = ()):
         self._flags = ''
@@ -310,8 +308,6 @@ class PyModule_AddObject(PyccelInternalFunction):
     """
     __slots__ = ('_mod_name','_name','_var')
     _attribute_nodes = ('_name','_var')
-    #_precision = 4
-    _rank = 0
     _shape = None
     _class_type = PythonNativeInt()
 
@@ -361,9 +357,7 @@ class PyModule_Create(PyccelInternalFunction):
     """
     __slots__ = ('_module_def_name',)
     _attribute_nodes = ()
-    _rank = 0
     _shape = ()
-    _order = None
     _class_type = PyccelPyObject()
 
     def __init__(self, module_def_name):
@@ -403,9 +397,7 @@ class PyCapsule_New(PyccelInternalFunction):
     """
     __slots__ = ('_capsule_name', '_API_var')
     _attribute_nodes = ('_API_var',)
-    _rank = 0
     _shape = ()
-    _order = None
     _class_type = PyccelPyObject()
 
     def __init__(self, API_var, module_name):
@@ -453,9 +445,7 @@ class PyCapsule_Import(PyccelInternalFunction):
     """
     __slots__ = ('_capsule_name',)
     _attribute_nodes = ()
-    _rank = 0
     _shape = ()
-    _order = None
     _class_type = BindCPointer()
 
     def __init__(self, module_name):
diff --git a/pyccel/ast/datatypes.py b/pyccel/ast/datatypes.py
index 03f9576f1a..290b1e07c2 100644
--- a/pyccel/ast/datatypes.py
+++ b/pyccel/ast/datatypes.py
@@ -209,6 +209,27 @@ def primitive_type(self):
         """
         return self._primitive_type # pylint: disable=no-member
 
+    @property
+    def rank(self):
+        """
+        Number of dimensions of the object.
+
+        Number of dimensions of the object. If the object is a scalar then
+        this is equal to 0.
+        """
+        return 0
+
+    @property
+    def order(self):
+        """
+        The data layout ordering in memory.
+
+        Indicates whether the data is stored in row-major ('C') or column-major
+        ('F') format. This is only relevant if rank > 1. When it is not relevant
+        this function returns None.
+        """
+        return None
+
     def __reduce__(self):
         """
         Function called during pickling.
@@ -556,6 +577,69 @@ def switch_basic_type(self, new_type):
         cls = type(self)
         return cls(self.element_type.switch_basic_type(new_type))
 
+    def switch_rank(self, new_rank, new_order = None):
+        """
+        Get a type which is identical to this type in all aspects except the rank.
+
+        Get a type which is identical to this type in all aspects except the rank.
+        The order must be provided if the rank is increased from 1. This is never
+        the case for 1D containers.
+
+        Parameters
+        ----------
+        new_rank : int
+            The rank of the new type.
+
+        new_order : str, optional
+            The order of the new type. For 1D containers this should not be provided.
+
+        Returns
+        -------
+        PyccelType
+            The new type.
+        """
+        assert new_order is None
+        rank = self.rank
+        assert new_rank < rank
+
+        if new_rank == rank:
+            return self
+        elif rank - new_rank == self.container_rank:
+            return self.element_type
+        else:
+            return self.element_type.switch_rank(new_rank - self.container_rank)
+
+    @property
+    def container_rank(self):
+        """
+        Number of dimensions of the container.
+
+        Number of dimensions of the object described by the container. This is
+        equal to the number of values required to index an element of this container.
+        """
+        return self._container_rank # pylint: disable=no-member
+
+    @property
+    def rank(self):
+        """
+        Number of dimensions of the object.
+
+        Number of dimensions of the object. If the object is a scalar then
+        this is equal to 0.
+        """
+        return self.container_rank + self.element_type.rank
+
+    @property
+    def order(self):
+        """
+        The data layout ordering in memory.
+
+        Indicates whether the data is stored in row-major ('C') or column-major
+        ('F') format. This is only relevant if rank > 1. When it is not relevant
+        this function returns None.
+        """
+        return self._order # pylint: disable=no-member
+
 class StringType(HomogeneousContainerType, metaclass = Singleton):
     """
     Class representing Python's native string type.
@@ -565,6 +649,8 @@ class StringType(HomogeneousContainerType, metaclass = Singleton):
     __slots__ = ()
     _name = 'str'
     _element_type = PrimitiveCharacterType()
+    _container_rank = 1
+    _order = None
 
     @property
     def datatype(self):
@@ -603,6 +689,16 @@ def primitive_type(self):
         """
         return self
 
+    @property
+    def rank(self):
+        """
+        Number of dimensions of the object.
+
+        Number of dimensions of the object. If the object is a scalar then
+        this is equal to 0.
+        """
+        return 1
+
 class HomogeneousTupleType(HomogeneousContainerType, TupleType, metaclass = ArgumentSingleton):
     """
     Class representing the homogeneous tuple type.
@@ -615,11 +711,13 @@ class HomogeneousTupleType(HomogeneousContainerType, TupleType, metaclass = Argu
     element_type : PyccelType
         The type of the elements of the homogeneous tuple.
     """
-    __slots__ = ('_element_type',)
+    __slots__ = ('_element_type', '_order')
+    _container_rank = 1
 
     def __init__(self, element_type):
         assert isinstance(element_type, PyccelType)
         self._element_type = element_type
+        self._order = 'C' if (element_type.order == 'C' or element_type.rank == 1) else None
         super().__init__()
 
     def __str__(self):
@@ -637,12 +735,14 @@ class HomogeneousListType(HomogeneousContainerType, metaclass = ArgumentSingleto
     element_type : PyccelType
         The type which is stored in the homogeneous list.
     """
-    __slots__ = ('_element_type',)
+    __slots__ = ('_element_type', '_order')
     _name = 'list'
+    _container_rank = 1
 
     def __init__(self, element_type):
         assert isinstance(element_type, PyccelType)
         self._element_type = element_type
+        self._order = 'C' if (element_type.order == 'C' or element_type.rank == 1) else None
         super().__init__()
 
 class HomogeneousSetType(HomogeneousContainerType, metaclass = ArgumentSingleton):
@@ -659,6 +759,8 @@ class HomogeneousSetType(HomogeneousContainerType, metaclass = ArgumentSingleton
     """
     __slots__ = ('_element_type',)
     _name = 'set'
+    _container_rank = 1
+    _order = None
 
     def __init__(self, element_type):
         assert isinstance(element_type, PyccelType)
@@ -701,6 +803,27 @@ def __reduce__(self):
         """
         return (self.__class__, ())
 
+    @property
+    def rank(self):
+        """
+        Number of dimensions of the object.
+
+        Number of dimensions of the object. If the object is a scalar then
+        this is equal to 0.
+        """
+        return 0
+
+    @property
+    def order(self):
+        """
+        The data layout ordering in memory.
+
+        Indicates whether the data is stored in row-major ('C') or column-major
+        ('F') format. This is only relevant if rank > 1. When it is not relevant
+        this function returns None.
+        """
+        return None
+
 class InhomogeneousTupleType(ContainerType, TupleType, metaclass = ArgumentSingleton):
     """
     Class representing the inhomogeneous tuple type.
@@ -714,14 +837,35 @@ class InhomogeneousTupleType(ContainerType, TupleType, metaclass = ArgumentSingl
     *args : tuple of DataTypes
         The datatypes stored in the inhomogeneous tuple.
     """
-    __slots__ = ('_element_types', '_datatype')
+    __slots__ = ('_element_types', '_datatype', '_container_rank', '_order')
 
     def __init__(self, *args):
         self._element_types = args
 
+        # Determine datatype
         possible_types = set(t.datatype for t in self._element_types)
         dtype = possible_types.pop()
         self._datatype = dtype if all(d == dtype for d in possible_types) else self
+
+        # Determine rank
+        elem_ranks = set(elem.rank for elem in self._element_types)
+        if len(elem_ranks) == 1:
+            self._container_rank = elem_ranks.pop() + 1
+        else:
+            self._container_rank = 1
+
+        # Determine order
+        if self._container_rank == 2:
+            self._order = 'C'
+        elif self._container_rank > 2:
+            elem_orders = set(elem.order for elem in self._element_types)
+            if len(elem_orders) == 1 and elem_orders.pop() == 'C':
+                self._order = 'C'
+            else:
+                self._order = None
+        else:
+            self._order = None
+
         super().__init__()
 
     def __str__(self):
@@ -764,6 +908,37 @@ def datatype(self):
         """
         return self._datatype
 
+    @property
+    def container_rank(self):
+        """
+        Number of dimensions of the container.
+
+        Number of dimensions of the object described by the container. This is
+        equal to the number of values required to index an element of this container.
+        """
+        return 1
+
+    @property
+    def rank(self):
+        """
+        Number of dimensions of the object.
+
+        Number of dimensions of the object. If the object is a scalar then
+        this is equal to 0.
+        """
+        return self._container_rank
+
+    @property
+    def order(self):
+        """
+        The data layout ordering in memory.
+
+        Indicates whether the data is stored in row-major ('C') or column-major
+        ('F') format. This is only relevant if rank > 1. When it is not relevant
+        this function returns None.
+        """
+        return self._order
+
 class DictType(ContainerType, metaclass = ArgumentSingleton):
     """
     Class representing the homogeneous dictionary type.
@@ -780,6 +955,8 @@ class DictType(ContainerType, metaclass = ArgumentSingleton):
     """
     __slots__ = ('_index_type', '_value_type')
     _name = 'map'
+    _container_rank = 1
+    _order = None
 
     def __init__(self, index_type, value_type):
         self._index_type = index_type
diff --git a/pyccel/ast/functionalexpr.py b/pyccel/ast/functionalexpr.py
index 4d8e2764f0..5498a92960 100644
--- a/pyccel/ast/functionalexpr.py
+++ b/pyccel/ast/functionalexpr.py
@@ -55,7 +55,7 @@ class FunctionalFor(TypedAstNode):
               Index is `Dummy_0`.
     """
     __slots__ = ('_loops','_expr', '_lhs', '_indices','_index',
-            '_rank','_shape','_order','_class_type')
+            '_shape','_class_type')
     _attribute_nodes = ('_loops','_expr', '_lhs', '_indices','_index')
 
     def __init__(
@@ -74,9 +74,7 @@ def __init__(
         super().__init__()
 
         if pyccel_stage != 'syntactic':
-            self._rank       = lhs.rank
             self._shape      = lhs.shape
-            self._order      = lhs.order
             self._class_type = lhs.class_type
 
     @property
diff --git a/pyccel/ast/internals.py b/pyccel/ast/internals.py
index ead8fc1dba..f502d9abfe 100644
--- a/pyccel/ast/internals.py
+++ b/pyccel/ast/internals.py
@@ -82,9 +82,7 @@ class PyccelArraySize(PyccelInternalFunction):
     __slots__ = ()
     name = 'size'
 
-    _rank  = 0
     _shape = None
-    _order = None
     _class_type = PythonNativeInt()
 
     def __init__(self, arg):
@@ -128,9 +126,7 @@ class PyccelArrayShapeElement(PyccelInternalFunction):
     __slots__ = ()
     name = 'shape'
 
-    _rank  = 0
     _shape = None
-    _order = None
     _class_type = PythonNativeInt()
 
     def __init__(self, arg, index):
diff --git a/pyccel/ast/literals.py b/pyccel/ast/literals.py
index fac53ca364..dd9f662513 100644
--- a/pyccel/ast/literals.py
+++ b/pyccel/ast/literals.py
@@ -40,9 +40,7 @@ class Literal(TypedAstNode):
     """
     __slots__ = ()
     _attribute_nodes  = ()
-    _rank      = 0
     _shape     = None
-    _order     = None
 
     @property
     def python_value(self):
diff --git a/pyccel/ast/macros.py b/pyccel/ast/macros.py
index 013f5c74fb..460fd98cde 100644
--- a/pyccel/ast/macros.py
+++ b/pyccel/ast/macros.py
@@ -12,6 +12,7 @@
 from .basic          import TypedAstNode
 from .datatypes      import PythonNativeInt, GenericType
 from .internals      import PyccelSymbol
+from .numpytypes     import NumpyNDArrayType, NumpyInt
 from .variable       import Variable
 
 pyccel_stage = PyccelStage()
@@ -69,20 +70,20 @@ class MacroShape(Macro):
     index : int | LiteralInteger
         The index of the element of the shape we are accessing.
     """
-    __slots__ = ('_index','_rank','_shape')
+    __slots__ = ('_index','_shape','_class_type')
     _name      = 'shape'
     _order     = None
-    _class_type = PythonNativeInt()
 
     def __init__(self, argument, index=None):
         if index is not None:
-            self._rank = 0
+            self._class_type = PythonNativeInt()
             self._shape = None
         elif pyccel_stage != "syntactic":
-            self._rank      = int(argument.rank>1)
+            rank      = int(argument.rank>1)
             self._shape     = (argument.rank,)
+            self._class_type = NumpyNDArrayType(NumpyInt, rank, None)
         else:
-            self._rank      = 1
+            self._class_type = NumpyNDArrayType(NumpyInt, 0, None)
             self._shape     = ()
         self._index = index
         super().__init__(argument)
@@ -111,9 +112,7 @@ class MacroType(Macro):
     """
     __slots__ = ()
     _name      = 'dtype'
-    _rank      = 0
     _shape     = None
-    _order     = None
     _class_type = GenericType()
 
     def __str__(self):
@@ -134,9 +133,7 @@ class MacroCount(Macro):
     """
     __slots__ = ()
     _name      = 'count'
-    _rank      = 0
     _shape     = None
-    _order     = None
     _class_type = PythonNativeInt()
 
     def __str__(self):
diff --git a/pyccel/ast/mathext.py b/pyccel/ast/mathext.py
index 491288126d..c48c249da9 100644
--- a/pyccel/ast/mathext.py
+++ b/pyccel/ast/mathext.py
@@ -89,8 +89,6 @@ class MathFunctionBase(PyccelInternalFunction):
     """
     __slots__ = ()
     _shape = None
-    _rank  = 0
-    _order = None
 
 
 class MathFunctionFloat(MathFunctionBase):
diff --git a/pyccel/ast/numpy_wrapper.py b/pyccel/ast/numpy_wrapper.py
index afe71b88e6..deabb3a645 100644
--- a/pyccel/ast/numpy_wrapper.py
+++ b/pyccel/ast/numpy_wrapper.py
@@ -9,6 +9,8 @@
 
 import numpy as np
 
+from .bind_c            import BindCPointer
+
 from .datatypes         import PythonNativeBool, GenericType, VoidType, FixedSizeType
 
 from .cwrapper          import PyccelPyObject, check_type_registry, c_to_py_registry, pytype_parse_registry
@@ -93,7 +95,7 @@ def get_numpy_max_acceptable_version_file():
                 name      = 'pyarray_to_ndarray',
                 arguments = [FunctionDefArgument(Variable(PyccelPyObject(), 'a', memory_handling = 'alias'))],
                 body      = [],
-                results   = [FunctionDefResult(Variable(NumpyNDArrayType(GenericType()), 'array'))])
+                results   = [FunctionDefResult(Variable(NumpyNDArrayType(GenericType(), 1, None), 'array'))])
 
 # NumPy array check elements : function definition in pyccel/stdlib/cwrapper/cwrapper_ndarrays.c
 pyarray_check = FunctionDef(
@@ -141,7 +143,7 @@ def get_numpy_max_acceptable_version_file():
 array_get_data  = FunctionDef(name   = 'nd_data',
                            body      = [],
                            arguments = [FunctionDefArgument(Variable(VoidType(), name = 'o', is_optional=True))],
-                           results   = [FunctionDefResult(Variable(VoidType(), name = 'v', memory_handling='alias', rank = 1))])
+                           results   = [FunctionDefResult(Variable(BindCPointer(), name = 'v', memory_handling='alias'))])
 
 PyArray_SetBaseObject = FunctionDef(name   = 'PyArray_SetBaseObject',
                                     body      = [],
diff --git a/pyccel/ast/numpyext.py b/pyccel/ast/numpyext.py
index 13f7a13441..b1a889ba90 100644
--- a/pyccel/ast/numpyext.py
+++ b/pyccel/ast/numpyext.py
@@ -183,15 +183,15 @@ class NumpyFloat(PythonFloat):
     arg : TypedAstNode
         The argument passed to the function.
     """
-    __slots__ = ('_rank','_shape','_order','_class_type')
+    __slots__ = ('_shape','_class_type')
     _static_type = NumpyFloat64Type()
     name = 'float'
 
     def __init__(self, arg):
         self._shape = arg.shape
-        self._rank  = arg.rank
-        self._order = arg.order
-        self._class_type = NumpyNDArrayType(self.static_type()) if self._rank else self.static_type()
+        rank  = arg.rank
+        order = arg.order
+        self._class_type = NumpyNDArrayType(self.static_type(), rank, order) if rank else self.static_type()
         super().__init__(arg)
 
     @property
@@ -246,13 +246,13 @@ class NumpyBool(PythonBool):
     arg : TypedAstNode
         The argument passed to the function.
     """
-    __slots__ = ('_shape','_rank','_order','_class_type')
+    __slots__ = ('_shape','_class_type')
     name = 'bool'
     def __init__(self, arg):
         self._shape = arg.shape
-        self._rank  = arg.rank
-        self._order = arg.order
-        self._class_type = NumpyNDArrayType(self.static_type()) if self._rank else self.static_type()
+        rank  = arg.rank
+        order = arg.order
+        self._class_type = NumpyNDArrayType(self.static_type(), rank, order) if rank else self.static_type()
         super().__init__(arg)
 
     @property
@@ -278,15 +278,15 @@ class NumpyInt(PythonInt):
     arg : TypedAstNode
         The argument passed to the function.
     """
-    __slots__ = ('_shape','_rank','_order','_class_type')
+    __slots__ = ('_shape','_class_type')
     _static_type = numpy_precision_map[(PrimitiveIntegerType(), PythonInt._static_type.precision)]
     name = 'int'
 
     def __init__(self, arg=None, base=10):
         self._shape = arg.shape
-        self._rank  = arg.rank
-        self._order = arg.order
-        self._class_type = NumpyNDArrayType(self.static_type()) if self._rank else self.static_type()
+        rank  = arg.rank
+        order = arg.order
+        self._class_type = NumpyNDArrayType(self.static_type(), rank, order) if rank else self.static_type()
         super().__init__(arg)
 
     @property
@@ -374,7 +374,7 @@ class NumpyReal(PythonReal):
     arg : TypedAstNode
         The argument passed to the function.
     """
-    __slots__ = ('_rank','_shape','_order','_class_type')
+    __slots__ = ('_shape','_class_type')
     name = 'real'
     def __new__(cls, arg):
         if isinstance(arg.dtype, PythonNativeBool):
@@ -387,10 +387,11 @@ def __new__(cls, arg):
 
     def __init__(self, arg):
         super().__init__(arg)
-        self._order = arg.order
-        self._rank  = arg.rank
-        self._shape = process_shape(self._rank == 0, self.internal_var.shape)
-        self._class_type = arg.class_type.switch_basic_type(arg.dtype.element_type)
+        rank  = arg.rank
+        order = arg.order
+        dtype = arg.dtype.element_type
+        self._class_type = NumpyNDArrayType(dtype, rank, order) if rank else dtype
+        self._shape = process_shape(self.rank == 0, self.internal_var.shape)
 
     @property
     def is_elemental(self):
@@ -415,7 +416,7 @@ class NumpyImag(PythonImag):
     arg : TypedAstNode
         The argument passed to the function.
     """
-    __slots__ = ('_rank','_shape','_order','_class_type')
+    __slots__ = ('_shape','_class_type')
     name = 'imag'
     def __new__(cls, arg):
 
@@ -429,10 +430,11 @@ def __new__(cls, arg):
 
     def __init__(self, arg):
         super().__init__(arg)
-        self._order = arg.order
-        self._rank  = arg.rank
-        self._shape = process_shape(self._rank == 0, self.internal_var.shape)
-        self._class_type = arg.class_type.switch_basic_type(arg.dtype.element_type)
+        rank  = arg.rank
+        order = arg.order
+        dtype = arg.dtype.element_type
+        self._class_type = NumpyNDArrayType(dtype, rank, order) if rank else dtype
+        self._shape = process_shape(self.rank == 0, self.internal_var.shape)
 
     @property
     def is_elemental(self):
@@ -458,7 +460,7 @@ class NumpyComplex(PythonComplex):
     """
     _real_cast = NumpyReal
     _imag_cast = NumpyImag
-    __slots__ = ('_rank','_shape','_order','_class_type')
+    __slots__ = ('_shape','_class_type')
     _static_type = NumpyComplex128Type()
     name = 'complex'
 
@@ -466,9 +468,9 @@ def __init__(self, arg0, arg1 = None):
         if arg1 is not None:
             raise NotImplementedError("Use builtin complex function not deprecated np.complex")
         self._shape = arg0.shape
-        self._rank  = arg0.rank
-        self._order = arg0.order
-        self._class_type = NumpyNDArrayType(self.static_type()) if self._rank else self.static_type()
+        rank  = arg0.rank
+        order = arg0.order
+        self._class_type = NumpyNDArrayType(self.static_type(), rank, order) if rank else self.static_type()
         super().__init__(arg0)
 
     @property
@@ -534,9 +536,7 @@ class NumpyResultType(PyccelInternalFunction):
         and no dtypes).
     """
     __slots__ = ('_class_type',)
-    _rank = 0
     _shape = None
-    _order = None
     name = 'result_type'
 
     def __init__(self, *arrays_and_dtypes):
@@ -615,16 +615,17 @@ class NumpyNewArray(PyccelInternalFunction):
     ----------
     *args : tuple of TypedAstNode
         The arguments of the superclass PyccelInternalFunction.
-    dtype : PyccelType
-        The datatype of the new array.
+    class_type : NumpyNDArrayType
+        The type of the new array.
     init_dtype : PythonType, PyccelFunctionDef, LiteralString, str
         The actual dtype passed to the NumPy function.
     """
     __slots__ = ('_init_dtype','_class_type')
 
-    def __init__(self, *args, dtype, init_dtype = None):
+    def __init__(self, *args, class_type, init_dtype = None):
+        assert isinstance(class_type, NumpyNDArrayType)
         self._init_dtype = init_dtype
-        self._class_type = NumpyNDArrayType(dtype) # pylint: disable=no-member
+        self._class_type = class_type # pylint: disable=no-member
 
         super().__init__(*args)
 
@@ -690,7 +691,7 @@ class NumpyArray(NumpyNewArray):
         The minimum number of dimensions that the resulting array should
         have.
     """
-    __slots__ = ('_arg','_shape','_rank','_order')
+    __slots__ = ('_arg','_shape')
     _attribute_nodes = ('_arg',)
     name = 'array'
 
@@ -702,9 +703,9 @@ def __init__(self, arg, dtype=None, order='K', ndmin=None):
         is_homogeneous_tuple = isinstance(arg.class_type, HomogeneousTupleType)
         # Inhomogeneous tuples can contain homogeneous data if it is inhomogeneous due to pointers
         if isinstance(arg.class_type, InhomogeneousTupleType):
+            is_homogeneous_tuple = isinstance(arg.dtype, FixedSizeNumericType) and len(set(a.rank for a in arg))
             if not isinstance(arg, PythonTuple):
                 arg = PythonTuple(*arg)
-            is_homogeneous_tuple = arg.is_homogeneous
 
         # TODO: treat inhomogenous lists and tuples when they have mixed ordering
         if not (is_homogeneous_tuple or isinstance(arg.class_type, HomogeneousContainerType)):
@@ -758,9 +759,7 @@ def __init__(self, arg, dtype=None, order='K', ndmin=None):
 
         self._arg   = arg
         self._shape = shape
-        self._rank  = rank
-        self._order = order
-        super().__init__(dtype = dtype, init_dtype = init_dtype)
+        super().__init__(class_type = NumpyNDArrayType(dtype, rank, order), init_dtype = init_dtype)
 
     def __str__(self):
         return str(self.arg)
@@ -793,8 +792,6 @@ class NumpyArange(NumpyNewArray):
     """
     __slots__ = ('_start','_step','_stop','_shape')
     _attribute_nodes = ('_start','_step','_stop')
-    _rank = 1
-    _order = None
     name = 'arange'
 
     def __init__(self, start, stop = None, step = None, dtype = None):
@@ -814,7 +811,8 @@ def __init__(self, start, stop = None, step = None, dtype = None):
 
         self._shape = (MathCeil(PyccelDiv(PyccelMinus(self._stop, self._start), self._step)))
         self._shape = process_shape(False, self._shape)
-        super().__init__(dtype = process_dtype(dtype), init_dtype = init_dtype)
+        dtype = process_dtype(dtype)
+        super().__init__(class_type = NumpyNDArrayType(dtype, 1, None), init_dtype = init_dtype)
 
     @property
     def arg(self):
@@ -850,9 +848,7 @@ class NumpySum(PyccelInternalFunction):
     """
     __slots__ = ('_class_type',)
     name = 'sum'
-    _rank  = 0
     _shape = None
-    _order = None
 
     def __init__(self, arg):
         if not isinstance(arg, TypedAstNode):
@@ -883,9 +879,7 @@ class NumpyProduct(PyccelInternalFunction):
     """
     __slots__ = ('_arg','_class_type')
     name = 'product'
-    _rank  = 0
     _shape = None
-    _order = None
 
     def __init__(self, arg):
         if not isinstance(arg, TypedAstNode):
@@ -921,7 +915,7 @@ class NumpyMatmul(PyccelInternalFunction):
     b : TypedAstNode
         The second argument of the matrix multiplication.
     """
-    __slots__ = ('_shape','_rank','_order','_class_type')
+    __slots__ = ('_shape','_class_type')
     name = 'matmul'
 
     def __init__(self, a ,b):
@@ -945,21 +939,21 @@ def __init__(self, a ,b):
             self._shape = (m, n)
 
         if a.rank == 1 and b.rank == 1:
-            self._rank  = 0
+            rank  = 0
             self._shape = None
         elif a.rank == 1 or b.rank == 1:
-            self._rank  = 1
+            rank  = 1
             self._shape = (b.shape[1] if a.rank == 1 else a.shape[0],)
         else:
-            self._rank = 2
+            rank = 2
 
 
         if a.order == b.order:
-            self._order = a.order
+            order = a.order
         else:
-            self._order = None if self._rank < 2 else 'C'
+            order = None if rank < 2 else 'C'
 
-        self._class_type = NumpyNDArrayType(dtype) if self.rank else dtype
+        self._class_type = NumpyNDArrayType(dtype, rank, order) if rank else dtype
 
     @property
     def a(self):
@@ -1028,9 +1022,8 @@ class NumpyLinspace(NumpyNewArray):
            from start and stop, the calculated dtype will never be an integer.
     """
 
-    __slots__ = ('_index','_start','_stop',
-            '_num','_endpoint','_shape', '_rank','_ind','_step',
-            '_py_argument','_order')
+    __slots__ = ('_index','_start','_stop', '_num','_endpoint','_shape', '_ind',
+            '_step', '_py_argument')
     _attribute_nodes = ('_start', '_stop', '_index', '_step', '_num',
             '_endpoint', '_ind')
     name = 'linspace'
@@ -1074,8 +1067,8 @@ def __init__(self, start, stop, num=None, endpoint=True, dtype=None):
         self._shape = (self._num,)
         if shape is not None:
             self._shape += shape
-        self._rank  = len(self._shape)
-        self._order = None if self._rank < 2 else 'C'
+        rank  = len(self._shape)
+        order = None if rank < 2 else 'C'
 
         self._ind = None
 
@@ -1086,7 +1079,9 @@ def __init__(self, start, stop, num=None, endpoint=True, dtype=None):
         else:
             self._step = PyccelDiv(PyccelMinus(self.stop, self.start), PyccelMinus(self.num, PythonInt(self.endpoint)))
 
-        super().__init__(dtype = final_dtype, init_dtype = init_dtype)
+        class_type = NumpyNDArrayType(final_dtype, rank, order)
+
+        super().__init__(class_type = class_type, init_dtype = init_dtype)
 
     @property
     def endpoint(self):
@@ -1153,7 +1148,7 @@ class NumpyWhere(PyccelInternalFunction):
     """
 
     __slots__ = ('_condition', '_value_true', '_value_false',
-                 '_rank', '_shape', '_order', '_class_type')
+                 '_shape', '_class_type')
     _attribute_nodes = ('_condition','_value_true','_value_false')
     name = 'where'
 
@@ -1172,14 +1167,14 @@ def __init__(self, condition, x, y):
 
         args      = (x, y)
         type_info = NumpyResultType(*args)
-        self._class_type = NumpyNDArrayType(process_dtype(type_info.dtype))
 
         shape = broadcast(x.shape, y.shape)
         shape = broadcast(condition.shape, shape)
 
         self._shape = process_shape(False, shape)
-        self._rank  = len(shape)
-        self._order = None if self._rank < 2 else 'C'
+        rank  = len(shape)
+        order = None if rank < 2 else 'C'
+        self._class_type = NumpyNDArrayType(process_dtype(type_info.dtype), rank, order)
         super().__init__(condition, x, y)
 
     @property
@@ -1216,18 +1211,18 @@ class NumpyRand(PyccelInternalFunction):
     *args : tuple of TypedAstNode
         The arguments passed to the function.
     """
-    __slots__ = ('_shape','_rank','_order','_class_type')
+    __slots__ = ('_shape','_class_type')
     name = 'rand'
 
     def __init__(self, *args):
         super().__init__(*args)
-        self._rank  = len(args)
-        self._shape = None if self._rank == 0 else args
-        self._order = None if self._rank < 2 else 'C'
-        if self.rank == 0:
+        rank  = len(args)
+        self._shape = None if rank == 0 else args
+        if rank == 0:
             self._class_type = PythonNativeFloat()
         else:
-            self._class_type = NumpyNDArrayType(NumpyFloat64Type())
+            order = None if rank < 2 else 'C'
+            self._class_type = NumpyNDArrayType(NumpyFloat64Type(), rank, order)
 
 #==============================================================================
 class NumpyRandint(PyccelInternalFunction):
@@ -1247,7 +1242,7 @@ class NumpyRandint(PyccelInternalFunction):
     size : TypedAstNode, optional
         The size of the array that will be generated.
     """
-    __slots__ = ('_rand','_low','_high','_shape','_rank','_order','_class_type')
+    __slots__ = ('_rand','_low','_high','_shape','_class_type')
     name = 'randint'
     _attribute_nodes = ('_low', '_high')
 
@@ -1261,12 +1256,11 @@ def __init__(self, low, high = None, size = None):
 
         self._shape   = size
         if size is None:
-            self._rank = 0
             self._class_type = PythonNativeInt()
         else:
-            self._rank = len(self.shape)
-            self._class_type = NumpyNDArrayType(NumpyInt64Type())
-        self._order   = None if self._rank < 2 else 'C'
+            rank = len(self.shape)
+            order = None if rank < 2 else 'C'
+            self._class_type = NumpyNDArrayType(NumpyInt64Type(), rank, order)
         self._rand    = NumpyRand() if size is None else NumpyRand(*size)
         self._low     = low
         self._high    = high
@@ -1312,7 +1306,7 @@ class NumpyFull(NumpyNewArray):
         Whether to store multidimensional data in C- or Fortran-contiguous
         (row- or column-wise) order in memory.
     """
-    __slots__ = ('_fill_value','_shape','_rank','_order')
+    __slots__ = ('_fill_value','_shape')
     name = 'full'
 
     def __init__(self, shape, fill_value, dtype=None, order='C'):
@@ -1333,10 +1327,12 @@ def __init__(self, shape, fill_value, dtype=None, order='C'):
                 cast_func = DtypePrecisionToCastFunction[dtype]
                 fill_value = cast_func(fill_value)
         self._shape = shape
-        self._rank  = len(self._shape)
-        self._order = NumpyNewArray._process_order(self._rank, order)
+        rank  = len(self._shape)
+        order = NumpyNewArray._process_order(rank, order)
 
-        super().__init__(fill_value, dtype = dtype, init_dtype = init_dtype)
+        class_type = NumpyNDArrayType(dtype, rank, order)
+
+        super().__init__(fill_value, class_type = class_type, init_dtype = init_dtype)
 
     #--------------------------------------------------------------------------
     @property
@@ -1641,7 +1637,7 @@ class NumpyNorm(PyccelInternalFunction):
         The second argument passed to the function, indicating the axis along
         which the norm should be calculated.
     """
-    __slots__ = ('_shape','_rank','_order','_arg','_class_type')
+    __slots__ = ('_shape','_arg','_class_type')
     name = 'norm'
 
     def __init__(self, arg, axis=None):
@@ -1657,13 +1653,12 @@ def __init__(self, arg, axis=None):
             sh = list(arg.shape)
             del sh[self.axis]
             self._shape = tuple(sh)
-            self._rank = len(self._shape)
-            self._order = None if self._rank < 2 else arg.order
+            rank = len(self._shape)
+            order = None if rank < 2 else arg.order
+            self._class_type = NumpyNDArrayType(dtype, rank, order) if rank else dtype
         else:
             self._shape = None
-            self._order = None
-            self._rank  = 0
-        self._class_type = NumpyNDArrayType(dtype) if self.rank else dtype
+            self._class_type = dtype
 
     @property
     def arg(self):
@@ -1701,7 +1696,7 @@ class NumpyUfuncBase(PyccelInternalFunction):
     *args : tuple of TypedAstNode
         The arguments passed to the function.
     """
-    __slots__ = ('_shape','_rank','_order','_class_type')
+    __slots__ = ('_shape','_class_type')
 
     @property
     def is_elemental(self):
@@ -1729,14 +1724,30 @@ class NumpyUfuncUnary(NumpyUfuncBase):
 
     def __init__(self, x):
         dtype = self._get_dtype(x)
-        self._set_shape_rank(x)
-        self._set_order(x)
-        self._class_type = NumpyNDArrayType(dtype) if self.rank else dtype
+        self._shape, rank = self._get_shape_rank(x)
+        order = self._get_order(x, rank)
+        self._class_type = NumpyNDArrayType(dtype, rank, order) if rank else dtype
         super().__init__(x)
 
-    def _set_shape_rank(self, x):
-        self._shape      = x.shape
-        self._rank       = x.rank
+    def _get_shape_rank(self, x):
+        """
+        Get the shape and rank of the result of the function.
+
+        Get the shape and rank of the result of the function.
+
+        Parameters
+        ----------
+        x : TypedAstNode
+            The argument passed to the function.
+
+        Returns
+        -------
+        shape : tuple[TypedAstNode]
+            The shape of the result of the function.
+        rank : int
+            The rank of the result of the function.
+        """
+        return x.shape, x.rank
 
     def _get_dtype(self, x):
         """
@@ -1760,8 +1771,26 @@ def _get_dtype(self, x):
         else:
             return numpy_precision_map[(x_dtype.primitive_type, x_dtype.precision)]
 
-    def _set_order(self, x):
-        self._order      = x.order
+    def _get_order(self, x, rank):
+        """
+        Get the order of the result of the function.
+
+        Get the order of the result of the function.
+
+        Parameters
+        ----------
+        x : TypedAstNode
+            The argument passed to the function.
+
+        rank : int
+            The rank of the result of the function calculated by _get_shape_rank.
+
+        Returns
+        -------
+        str
+            The order of the result of the function.
+        """
+        return x.order
 
     @property
     def arg(self):
@@ -1796,13 +1825,33 @@ class NumpyUfuncBinary(NumpyUfuncBase):
     def __init__(self, x1, x2):
         super().__init__(x1, x2)
         dtype = self._get_dtype(x1, x2)
-        self._set_shape_rank(x1, x2)
-        self._set_order(x1, x2)
-        self._class_type = NumpyNDArrayType(dtype) if self.rank else dtype
+        self._shape, rank = self._get_shape_rank(x1, x2)
+        order = self._get_order(x1, x2, rank)
+        self._class_type = NumpyNDArrayType(dtype, rank, order) if rank else dtype
 
-    def _set_shape_rank(self, x1, x2):
-        self._shape = broadcast(x1.shape, x2.shape)
-        self._rank  = 0 if self._shape is None else len(self._shape)
+    def _get_shape_rank(self, x1, x2):
+        """
+        Get the shape and rank of the result of the function.
+
+        Get the shape and rank of the result of the function.
+
+        Parameters
+        ----------
+        x1 : TypedAstNode
+            The first argument passed to the function.
+        x2 : TypedAstNode
+            The second argument passed to the function.
+
+        Returns
+        -------
+        shape : tuple[TypedAstNode]
+            The shape of the result of the function.
+        rank : int
+            The rank of the result of the function.
+        """
+        shape = broadcast(x1.shape, x2.shape)
+        rank  = 0 if shape is None else len(shape)
+        return shape, rank
 
     def _get_dtype(self, x1, x2):
         """
@@ -1829,11 +1878,30 @@ def _get_dtype(self, x1, x2):
             arg_dtype = x1.dtype + x2.dtype
             return numpy_precision_map[(PrimitiveFloatingPointType(), arg_dtype.precision)]
 
-    def _set_order(self, x1, x2):
+    def _get_order(self, x1, x2, rank):
+        """
+        Get the order of the result of the function.
+
+        Get the order of the result of the function.
+
+        Parameters
+        ----------
+        x1 : TypedAstNode
+            The first argument passed to the function.
+        x2 : TypedAstNode
+            The second argument passed to the function.
+        rank : int
+            The rank of the result of the function calculated by _get_shape_rank.
+
+        Returns
+        -------
+        str
+            The order of the result of the function.
+        """
         if x1.order == x2.order:
-            self._order = x1.order
+            return x1.order
         else:
-            self._order = None if self._rank < 2 else 'C'
+            return None if rank < 2 else 'C'
 
 #------------------------------------------------------------------------------
 # Math operations
@@ -2028,14 +2096,32 @@ def __init__(self, x1, x2):
         x2 = NumpyInt(x2) if isinstance(x2.dtype, PythonNativeBool) else x2
         self._args = (x1, x2)
 
-    def _set_shape_rank(self, x1, x2):
+    def _get_shape_rank(self, x1, x2):
+        """
+        Get the shape and rank of the result of the function.
+
+        Get the shape and rank of the result of the function.
+
+        Parameters
+        ----------
+        x1 : TypedAstNode
+            The first argument passed to the function.
+        x2 : TypedAstNode
+            The second argument passed to the function.
+
+        Returns
+        -------
+        shape : tuple[TypedAstNode]
+            The shape of the result of the function.
+        rank : int
+            The rank of the result of the function.
+        """
         args   = (x1, x2)
         ranks  = [a.rank  for a in args]
         shapes = [a.shape for a in args]
 
         if all(r == 0 for r in ranks):
-            self._rank  = 0
-            self._shape = None
+            return None, 0
         else:
             if len(args) == 1:
                 shape = args[0].shape
@@ -2045,8 +2131,7 @@ def _set_shape_rank(self, x1, x2):
                 for a in args[2:]:
                     shape = broadcast(shape, a.shape)
 
-            self._shape = shape
-            self._rank  = len(shape)
+            return shape, len(shape)
 
     def _get_dtype(self, x1, x2):
         """
@@ -2090,9 +2175,8 @@ class NumpyAmin(PyccelInternalFunction):
     """
     __slots__ = ('_class_type',)
     name = 'amin'
-    _rank = 0
     _shape = None
-    _order = None
+
     def __init__(self, arg):
         super().__init__(arg)
         self._class_type = arg.dtype
@@ -2119,9 +2203,8 @@ class NumpyAmax(PyccelInternalFunction):
     """
     __slots__ = ('_class_type',)
     name = 'amax'
-    _rank = 0
     _shape = None
-    _order = None
+
     def __init__(self, arg):
         super().__init__(arg)
         self._class_type = arg.dtype
@@ -2199,22 +2282,48 @@ def _get_dtype(self, x):
         """
         return process_dtype(x.dtype)
 
-    def _set_shape_rank(self, x):
+    def _get_shape_rank(self, x):
         """
-        Set the shape and rank of the resulting object.
+        Get the shape and rank of the result of the function.
 
-        Set the shape and rank of the resulting object.
+        Get the shape and rank of the result of the function.
 
         Parameters
         ----------
         x : TypedAstNode
             The argument passed to the function.
+
+        Returns
+        -------
+        shape : tuple[TypedAstNode]
+            The shape of the result of the function.
+        rank : int
+            The rank of the result of the function.
+        """
+        shape = tuple(reversed(x.shape))
+        rank  = x.rank
+        return shape, rank
+
+    def _get_order(self, x, rank):
         """
-        self._shape = tuple(reversed(x.shape))
-        self._rank  = x.rank
+        Get the order of the result of the function.
 
-    def _set_order(self, x):
-        self._order = 'C' if x.order=='F' else 'F'
+        Get the order of the result of the function.
+
+        Parameters
+        ----------
+        x : TypedAstNode
+            The argument passed to the function.
+
+        rank : int
+            The rank of the result of the function calculated by _get_shape_rank.
+
+        Returns
+        -------
+        str
+            The order of the result of the function.
+        """
+        return 'C' if x.order=='F' else 'F'
 
     @property
     def is_elemental(self):
@@ -2234,15 +2343,15 @@ class NumpyConjugate(PythonConjugate):
     arg : TypedAstNode
         The argument passed to the function.
     """
-    __slots__ = ('_rank','_shape','_order','_class_type')
+    __slots__ = ('_shape','_class_type')
     name = 'conj'
 
     def __init__(self, arg):
         super().__init__(arg)
-        self._order = arg.order
-        self._rank  = self.internal_var.rank
-        self._shape = process_shape(self._rank == 0, self.internal_var.shape)
-        self._class_type = NumpyNDArrayType(arg.dtype) if self.rank else arg.dtype
+        order = arg.order
+        rank  = arg.rank
+        self._shape = process_shape(rank == 0, arg.shape)
+        self._class_type = NumpyNDArrayType(arg.dtype, rank, order) if rank else arg.dtype
 
     @property
     def is_elemental(self):
@@ -2268,15 +2377,13 @@ class NumpyNonZeroElement(NumpyNewArray):
     __slots__ = ('_arr','_dim','_shape')
     _attribute_nodes = ('_arr',)
     name = 'nonzero'
-    _rank = 1
-    _order = None
 
     def __init__(self, a, dim):
         self._arr = a
         self._dim = dim
 
         self._shape = (NumpyCountNonZero(a),)
-        super().__init__(a, dtype = NumpyInt64Type())
+        super().__init__(a, class_type = NumpyNDArrayType(NumpyInt64Type(), 1, None))
 
     @property
     def array(self):
@@ -2311,9 +2418,7 @@ class NumpyNonZero(PyccelInternalFunction):
     __slots__ = ('_elements','_arr','_shape')
     _attribute_nodes = ('_elements',)
     name = 'nonzero'
-    _rank  = 2
-    _order = 'C'
-    _class_type = HomogeneousTupleType(NumpyNDArrayType(NumpyInt64Type()))
+    _class_type = HomogeneousTupleType(NumpyNDArrayType(NumpyInt64Type(), 1, None))
 
     def __init__(self, a):
         if (a.rank > 1):
@@ -2355,7 +2460,7 @@ class NumpyCountNonZero(PyccelInternalFunction):
         Indicates if output arrays should have the same number of dimensions
         as arg.
     """
-    __slots__ = ('_rank', '_shape', '_order', '_class_type', '_arr',
+    __slots__ = ('_shape', '_class_type', '_arr',
                 '_axis', '_keep_dims')
     _attribute_nodes = ('_arr','_axis')
     name   = 'count_nonzero'
@@ -2368,26 +2473,24 @@ def __init__(self, a, axis = None, *, keepdims = LiteralFalse()):
 
         if keepdims.python_value:
             dtype = NumpyInt64Type()
-            self._rank  = a.rank
-            self._order = a.order
+            rank  = a.rank
+            order = a.order
             if axis is not None:
                 self._shape = list(a.shape)
                 self._shape[axis.python_value] = LiteralInteger(1)
             else:
-                self._shape = (LiteralInteger(1),)*a.rank
-            self._class_type = NumpyNDArrayType(dtype) if self._rank else dtype
+                self._shape = (LiteralInteger(1),)*rank
+            self._class_type = NumpyNDArrayType(dtype, rank, order)
         else:
             if axis is not None:
                 dtype = NumpyInt64Type()
                 self._shape = list(a.shape)
                 self._shape.pop(axis.python_value)
-                self._rank  = a.rank-1
-                self._order = a.order if a.rank>2 else None
-                self._class_type = NumpyNDArrayType(dtype) if a.rank else dtype
+                rank  = a.rank-1
+                order = a.order
+                self._class_type = NumpyNDArrayType(dtype, rank, order)
             else:
-                self._rank  = 0
                 self._shape = None
-                self._order = None
                 self._class_type = PythonNativeInt()
 
         self._arr = a
diff --git a/pyccel/ast/numpytypes.py b/pyccel/ast/numpytypes.py
index 9c9e3e8295..e254471c96 100644
--- a/pyccel/ast/numpytypes.py
+++ b/pyccel/ast/numpytypes.py
@@ -267,14 +267,31 @@ class NumpyNDArrayType(HomogeneousContainerType, metaclass = ArgumentSingleton):
     dtype : NumpyNumericType | PythonNativeBool | GenericType
         The internal datatype of the object (GenericType is allowed for external
         libraries, e.g. MPI).
+    rank : int
+        The rank of the new NumPy array.
+    order : str
+        The order of the memory layout for the new NumPy array.
     """
-    __slots__ = ('_element_type',)
+    __slots__ = ('_element_type', '_container_rank', '_order')
     _name = 'numpy.ndarray'
 
-    def __init__(self, dtype):
+    def __new__(cls, dtype, rank, order):
+        if rank == 0:
+            return dtype
+        else:
+            return super().__new__(cls)
+
+    def __init__(self, dtype, rank, order):
+        assert isinstance(rank, int)
+        assert order in (None, 'C', 'F')
+        assert rank < 2 or order is not None
+
         if pyccel_stage == 'semantic':
             assert isinstance(dtype, (NumpyNumericType, PythonNativeBool, GenericType))
+
         self._element_type = dtype
+        self._container_rank = rank
+        self._order = order
         super().__init__()
 
     @lru_cache
@@ -287,7 +304,14 @@ def __add__(self, other):
         else:
             return NotImplemented
         result_type = original_type_to_pyccel_type[np.result_type(test_type, comparison_type).type]
-        return NumpyNDArrayType(result_type)
+        rank = max(other.rank, self.rank)
+        if rank < 2:
+            order = None
+        else:
+            other_f_contiguous = other.order in (None, 'F')
+            self_f_contiguous = self.order in (None, 'F')
+            order = 'F' if other_f_contiguous and self_f_contiguous else 'C'
+        return NumpyNDArrayType(result_type, rank, order)
 
     @lru_cache
     def __radd__(self, other):
@@ -329,7 +353,77 @@ def switch_basic_type(self, new_type):
         assert isinstance(new_type, FixedSizeNumericType)
         new_type = numpy_precision_map[(new_type.primitive_type, new_type.precision)]
         cls = type(self)
-        return cls(self.element_type.switch_basic_type(new_type))
+        return cls(self.element_type.switch_basic_type(new_type), self._container_rank, self._order)
+
+    def switch_rank(self, new_rank, new_order = None):
+        """
+        Get a type which is identical to this type in all aspects except the rank and/or order.
+
+        Get a type which is identical to this type in all aspects except the rank and/or order.
+        The order must be provided if the rank is increased from 1. Otherwise it defaults to the
+        same order as the current type.
+
+        Parameters
+        ----------
+        new_rank : int
+            The rank of the new type.
+
+        new_order : str, optional
+            The order of the new type. This should be provided if the rank is increased from 1.
+
+        Returns
+        -------
+        PyccelType
+            The new type.
+        """
+        if new_rank == 0:
+            return self.element_type
+        else:
+            new_order = (new_order or self._order) if new_rank > 1 else None
+            return NumpyNDArrayType(self.element_type, new_rank, new_order)
+
+    def swap_order(self):
+        """
+        Get a type which is identical to this type in all aspects except the order.
+
+        Get a type which is identical to this type in all aspects except the order.
+        In the case of a 1D array the final type will be the same as this type. Otherwise
+        if the array is C-ordered the final type will be F-ordered, while if the array
+        is F-ordered the final type will be C-ordered.
+
+        Returns
+        -------
+        PyccelType
+            The new type.
+        """
+        order = None if self._order is None else ('C' if self._order == 'F' else 'F')
+        return NumpyNDArrayType(self.element_type, self._container_rank, order)
+
+    @property
+    def rank(self):
+        """
+        Number of dimensions of the object.
+
+        Number of dimensions of the object. If the object is a scalar then
+        this is equal to 0.
+        """
+        return self._container_rank
+
+    @property
+    def order(self):
+        """
+        The data layout ordering in memory.
+
+        Indicates whether the data is stored in row-major ('C') or column-major
+        ('F') format. This is only relevant if rank > 1. When it is not relevant
+        this function returns None.
+        """
+        return self._order
+
+    def __repr__(self):
+        dims = ','.join(':'*self._container_rank)
+        order_str = f'(order={self._order})' if self._order else ''
+        return f'{self.element_type}[{dims}]{order_str}'
 
 #==============================================================================
 
@@ -368,3 +462,5 @@ def switch_basic_type(self, new_type):
 pyccel_type_to_original_type.update(numpy_type_to_original_type)
 original_type_to_pyccel_type.update({v:k for k,v in numpy_type_to_original_type.items()})
 original_type_to_pyccel_type[np.bool_] = PythonNativeBool()
+
+NumpyInt = numpy_precision_map[PrimitiveIntegerType(), np.dtype(int).alignment]
diff --git a/pyccel/ast/operators.py b/pyccel/ast/operators.py
index 9b93b0b7ea..dcb8702c02 100644
--- a/pyccel/ast/operators.py
+++ b/pyccel/ast/operators.py
@@ -160,10 +160,8 @@ def __init__(self, *args):
         if pyccel_stage == 'syntactic':
             super().__init__()
             return
+        self._set_shape()
         self._set_type()
-        self._set_shape_rank()
-        # rank is None for lambda functions
-        self._set_order()
         super().__init__()
 
     def _set_type(self):
@@ -176,15 +174,15 @@ def _set_type(self):
         """
         self._class_type = self._calculate_type(*self._args)  # pylint: disable=no-member
 
-    def _set_shape_rank(self):
+    def _set_shape(self):
         """
-        Set the shape and rank of the result of the operator.
+        Set the shape of the result of the operator.
 
-        Set the shape and rank of the result of the operator. This function
-        uses the static method `_shape_rank` to set these values. If the
-        values are class parameters in a sub-class, this method must be over-ridden.
+        Set the shape of the result of the operator. This function
+        uses the static method `_shape` to set these values. If the
+        values are class parameters in a sub-class, this method must be overridden.
         """
-        self._shape, self._rank = self._calculate_shape_rank(*self._args)  # pylint: disable=no-member
+        self._shape = self._calculate_shape(*self._args)  # pylint: disable=no-member
 
     @property
     def precedence(self):
@@ -246,24 +244,6 @@ def _handle_precedence(self, args):
     def __str__(self):
         return repr(self)
 
-    def _set_order(self):
-        """
-        Set the order of the result.
-
-        Set the order of the result of the operator.
-        This is chosen to match the arguments if they are in agreement.
-        Otherwise it defaults to 'C'.
-        """
-        if self.rank is not None and self.rank > 1:
-            orders = [a.order for a in self._args if a.order is not None]
-            my_order = orders[0]
-            if all(o == my_order for o in orders):
-                self._order = my_order
-            else:
-                self._order = 'C'
-        else:
-            self._order = None
-
     @property
     def args(self):
         """ Arguments of the operator
@@ -284,7 +264,7 @@ class PyccelUnaryOperator(PyccelOperator):
     arg : TypedAstNode
         The argument passed to the operator.
     """
-    __slots__ = ('_shape','_rank','_order','_class_type')
+    __slots__ = ('_shape','_class_type')
 
     def __init__(self, arg):
         super().__init__(arg)
@@ -311,12 +291,11 @@ def _calculate_type(arg):
         return arg.class_type
 
     @staticmethod
-    def _calculate_shape_rank(arg):
+    def _calculate_shape(arg):
         """
-        Calculate the shape and rank.
+        Calculate the shape.
 
-        Calculate the shape and rank.
-        They are chosen to match the argument
+        Calculate the shape. It is chosen to match the argument.
 
         Parameters
         ----------
@@ -325,14 +304,10 @@ def _calculate_shape_rank(arg):
 
         Returns
         -------
-        shape : tuple[TypedAstNode]
+        tuple[TypedAstNode]
             The shape of the resulting object.
-        rank : int
-            The rank of the resulting object.
         """
-        rank = arg.rank
-        shape = arg.shape
-        return shape, rank
+        return arg.shape
 
 #==============================================================================
 
@@ -413,12 +388,11 @@ def _set_type(self):
         """
 
     @staticmethod
-    def _calculate_shape_rank(arg):
+    def _calculate_shape(arg):
         """
-        Calculate the shape and rank.
+        Calculate the shape.
 
-        Calculate the shape and rank.
-        They are chosen to match the argument
+        Calculate the shape. It is chosen to match the argument.
 
         Parameters
         ----------
@@ -427,14 +401,10 @@ def _calculate_shape_rank(arg):
 
         Returns
         -------
-        shape : tuple[TypedAstNode]
+        tuple[TypedAstNode]
             The shape of the resulting object.
-        rank : int
-            The rank of the resulting object.
         """
-        rank = 0
-        shape = None
-        return shape, rank
+        return None
 
     def __repr__(self):
         return f'not {repr(self.args[0])}'
@@ -477,7 +447,7 @@ class PyccelBinaryOperator(PyccelOperator):
     arg2 : TypedAstNode
         The second argument passed to the operator.
     """
-    __slots__ = ('_shape','_rank','_order','_class_type')
+    __slots__ = ('_shape','_class_type')
 
     def __init__(self, arg1, arg2):
         super().__init__(arg1, arg2)
@@ -518,13 +488,13 @@ def _calculate_type(cls, arg1, arg2):
             raise TypeError(f'Cannot determine the type of ({arg1}, {arg2})') #pylint: disable=raise-missing-from
 
     @staticmethod
-    def _calculate_shape_rank(arg1, arg2):
+    def _calculate_shape(arg1, arg2):
         """
-        Calculate the shape and rank.
+        Calculate the shape.
 
         Strings must be scalars.
 
-        For numeric types the rank and shape is determined according
+        For numeric types the shape is determined according
         to NumPy broadcasting rules where possible.
 
         Parameters
@@ -536,24 +506,20 @@ def _calculate_shape_rank(arg1, arg2):
 
         Returns
         -------
-        shape : tuple[TypedAstNode]
+        tuple[TypedAstNode]
             The shape of the resulting object.
-        rank : int
-            The rank of the resulting object.
         """
         args = (arg1, arg2)
         strs = [a for a in args if isinstance(a.dtype, StringType)]
         if strs:
             other = [a for a in args if isinstance(a.dtype, FixedSizeNumericType)]
             assert len(other) == 0
-            rank  = 0
             shape = None
         else:
             s = broadcast(arg1.shape, arg2.shape)
 
             shape = s
-            rank  = 0 if s is None else len(s)
-        return shape, rank
+        return shape
 
 #==============================================================================
 
@@ -972,7 +938,7 @@ def _calculate_type(cls, arg1, arg2):
         elif len(possible_class_types) == 1:
             class_type = possible_class_types.pop().switch_basic_type(dtype)
         else:
-            description = f"({arg1} {cls.op} {arg2})" # pylint: disable=no-member
+            description = f"({arg1!r} {cls.op} {arg2!r})" # pylint: disable=no-member
             raise NotImplementedError("Can't deduce type for comparison operator"
                                       f" with multiple containers {description}")
         return class_type
@@ -1140,21 +1106,11 @@ class PyccelBooleanOperator(PyccelOperator):
     *args : tuple of TypedAstNode
         The arguments passed to the operator.
     """
-    _rank = 0
     _shape = None
-    _order = None
     _class_type = PythonNativeBool()
 
     __slots__ = ()
 
-    def _set_order(self):
-        """
-        Set the order of the result.
-
-        Set the order of the result of the operator. Nothing needs to
-        be done here as the order is a class variable.
-        """
-
     def _set_type(self):
         """
         Set the type of the result of the operator.
@@ -1163,12 +1119,12 @@ def _set_type(self):
         to be done here as the type is a class variable.
         """
 
-    def _set_shape_rank(self):
+    def _set_shape(self):
         """
-        Set the shape and rank of the result of the operator.
+        Set the shape of the result of the operator.
 
-        Set the shape and rank of the result of the operator. Nothing needs
-        to be done here as the shape and rank are class variables.
+        Set the shape of the result of the operator. Nothing needs
+        to be done here as the shape is a class variable.
         """
 
 #==============================================================================
@@ -1352,7 +1308,7 @@ class IfTernaryOperator(PyccelOperator):
     >>> IfTernaryOperator(PyccelGt(n > 1),  5,  2)
     IfTernaryOperator(PyccelGt(n > 1),  5,  2)
     """
-    __slots__ = ('_shape','_rank','_order','_class_type')
+    __slots__ = ('_shape','_class_type')
     _precedence = 3
 
     def __init__(self, cond, value_true, value_false):
@@ -1395,7 +1351,7 @@ def _calculate_type(cond, value_true, value_false):
         DataType
             The Python type of the object.
         """
-        if value_true.dtype is value_false.dtype and value_true.class_type is value_false.class_type:
+        if value_true.class_type is value_false.class_type:
             return value_true.class_type
 
         try:
@@ -1403,19 +1359,36 @@ def _calculate_type(cond, value_true, value_false):
         except NotImplementedError:
             raise TypeError(f'Cannot determine the type of ({value_true}, {value_false})') #pylint: disable=raise-missing-from
 
+        if value_false.class_type.order != value_true.class_type.order :
+            errors.report('Ternary Operator results should have the same order', severity='fatal')
+
         return class_type
 
     @staticmethod
-    def _calculate_shape_rank(cond, value_true, value_false):
+    def _calculate_shape(cond, value_true, value_false):
         """
-        Sets the shape and rank and the order for IfTernaryOperator
+        Calculate the shape of the result.
+
+        Calculate the shape of the result of the IfTernaryOperator.
+        The shape is equal to the shape of one of the outputs.
+
+        Parameters
+        ----------
+        cond : TypedAstNode
+            The first argument passed to the operator representing the condition.
+        value_true : TypedAstNode
+            The second argument passed to the operator representing the result if the
+            condition is true.
+        value_false : TypedAstNode
+            The third argument passed to the operator representing the result if the
+            condition is false.
+
+        Returns
+        -------
+        tuple[TypedAstNode]
+            The shape of the resulting object.
         """
-        shape = value_true.shape
-        rank  = value_true.rank
-        if rank is not None and rank > 1:
-            if value_false.order != value_true.order :
-                errors.report('Ternary Operator results should have the same order', severity='fatal')
-        return shape, rank
+        return value_true.shape
 
     @property
     def cond(self):
diff --git a/pyccel/ast/sysext.py b/pyccel/ast/sysext.py
index 2a45b48a19..73a017a725 100644
--- a/pyccel/ast/sysext.py
+++ b/pyccel/ast/sysext.py
@@ -32,9 +32,7 @@ class SysExit(PyccelInternalFunction):
     __slots__ = ()
     name      = 'exit'
     _class_type     = VoidType()
-    _rank      = 0
     _shape     = None
-    _order     = None
 
     def __init__(self, status=LiteralInteger(0)):
         super().__init__(status)
diff --git a/pyccel/ast/type_annotations.py b/pyccel/ast/type_annotations.py
index dcb7562e70..7a70c6c1c5 100644
--- a/pyccel/ast/type_annotations.py
+++ b/pyccel/ast/type_annotations.py
@@ -42,23 +42,13 @@ class VariableTypeAnnotation(PyccelAstNode):
     class_type : DataType
         The requested Python type of the variable.
 
-    rank : int
-        The rank of the variable.
-
-    order : str
-        The order of the variable.
-
     is_const : bool, default=False
         True if the variable cannot be modified, false otherwise.
     """
-    __slots__ = ('_class_type', '_rank',
-                 '_order', '_is_const')
+    __slots__ = ('_class_type', '_is_const')
     _attribute_nodes = ()
-    def __init__(self, class_type : 'DataType',
-            rank : int = 0, order : str = None, is_const : bool = False):
+    def __init__(self, class_type : 'DataType', is_const : bool = False):
         self._class_type = class_type
-        self._rank = rank
-        self._order = order
         self._is_const = is_const
 
         super().__init__()
@@ -77,28 +67,12 @@ def class_type(self):
     @property
     def rank(self):
         """
-        Get the rank of the object.
-
-        Get the rank of the object that should be created. The rank indicates the
-        number of dimensions.
-        """
-        return self._rank
-
-    @property
-    def order(self):
-        """
-        Get the order of the object.
+        Number of dimensions of the object.
 
-        Get the order in which the memory will be laid out in the object. For objects
-        with rank > 1 this is either 'C' or 'F'. 
+        Number of dimensions of the object. If the object is a scalar then
+        this is equal to 0.
         """
-        return self._order
-
-    @order.setter
-    def order(self, order):
-        if order not in ('C', 'F', None):
-            raise ValueError("Order must be C, F, or None")
-        self._order = order
+        return self.class_type.rank
 
     @property
     def is_const(self):
@@ -117,24 +91,17 @@ def is_const(self, val):
         self._is_const = val
 
     def __hash__(self):
-        return hash((self.class_type, self.rank, self.order))
+        return hash(self.class_type)
 
     def __eq__(self, other):
         # Needed for set
         if isinstance(other, VariableTypeAnnotation):
-            return self.class_type == other.class_type and \
-                   self.rank == other.rank and \
-                   self.order == other.order
+            return self.class_type == other.class_type
         else:
             return False
 
     def __repr__(self):
-        dtype = str(self._class_type)
-        if self._rank:
-            dtype += '['+','.join(':'*self._rank)+']'
-        if self._order:
-            dtype += '(order = {self._order})'
-        return dtype
+        return repr(self._class_type)
 
 #==============================================================================
 
diff --git a/pyccel/ast/utilities.py b/pyccel/ast/utilities.py
index eb14caeb78..dd7f917245 100644
--- a/pyccel/ast/utilities.py
+++ b/pyccel/ast/utilities.py
@@ -235,6 +235,36 @@ def compatible_operation(*args, language_has_vectors = True):
     else:
         return all(a.rank == 0 for a in args)
 
+#==============================================================================
+def get_deep_indexed_element(expr, indices):
+    """
+    Get the scalar element obtained by indexing the expression with all the indices.
+
+    Get the scalar element obtained by indexed the expression with all the provided
+    indices. This element is constructed by calling IndexedElement multiple times
+    to create a recursive object with one IndexedElement for each container type.
+    This function is used by the functions which unravel vector expressions.
+
+    Parameters
+    ----------
+    expr : TypedAstNode
+        The base object being indexed.
+    indices : list[TypedAstNode]
+        A list of the indices used to obtain the scalar element.
+
+    Returns
+    -------
+    IndexedElement
+        The scalar indexed element.
+    """
+    assert len(indices) == expr.rank
+    result = expr
+    while indices:
+        depth = result.class_type.container_rank
+        result = IndexedElement(result, *indices[:depth])
+        indices = indices[depth:]
+    return result
+
 #==============================================================================
 def insert_index(expr, pos, index_var):
     """
@@ -260,12 +290,13 @@ def insert_index(expr, pos, index_var):
     Examples
     --------
     >>> from pyccel.ast.core import Variable, Assign
+    >>> from pyccel.ast.datatypes import PythonNativeInt
     >>> from pyccel.ast.operators import PyccelAdd
     >>> from pyccel.ast.utilities import insert_index
-    >>> a = Variable('int', 'a', shape=(4,), rank=1)
-    >>> b = Variable('int', 'b', shape=(4,), rank=1)
-    >>> c = Variable('int', 'c', shape=(4,), rank=1)
-    >>> i = Variable('int', 'i')
+    >>> a = Variable(PythonNativeInt(), 'a', shape=(4,))
+    >>> b = Variable(PythonNativeInt(), 'b', shape=(4,))
+    >>> c = Variable(PythonNativeInt(), 'c', shape=(4,))
+    >>> i = Variable(PythonNativeInt(), 'i')
     >>> d = PyccelAdd(a,b)
     >>> expr = Assign(c,d)
     >>> insert_index(expr, 0, i)
@@ -282,7 +313,7 @@ def insert_index(expr, pos, index_var):
 
         # Add index at the required position
         indexes = [Slice(None,None)]*(expr.rank+pos) + [index_var]+[Slice(None,None)]*(-1-pos)
-        return IndexedElement(expr, *indexes)
+        return get_deep_indexed_element(expr, indexes)
 
     elif isinstance(expr, NumpyTranspose):
         if expr.rank==0 or -pos>expr.rank:
@@ -299,15 +330,30 @@ def insert_index(expr, pos, index_var):
 
     elif isinstance(expr, IndexedElement):
         base = expr.base
+        rank = base.rank
+
+        # If pos indexes base then recurse
+        base_container_rank = base.class_type.container_rank
+        if -pos < rank-base_container_rank:
+            return insert_index(base, pos+base_container_rank, index_var)
+
+        # Ensure current indices are fully defined
         indices = list(expr.indices)
         if len(indices) == 1 and isinstance(indices[0], LiteralEllipsis):
-            indices = [Slice(None,None)]*base.rank
-        i = -1
-        while i>=pos and -i<=base.rank:
+            indices = [Slice(None,None)]*base_container_rank
+
+        if len(indices)<rank:
+            indices += [Slice(None,None)]*(rank-base_container_rank)
+
+        # Start from last index in this indexed element
+        i = base_container_rank-rank-1
+        while i>=pos and -i<=base_container_rank:
             if not isinstance(indices[i], Slice):
                 pos -= 1
             i -= 1
-        if -pos>base.rank:
+
+        # if no slices were found then the object is already correctly indexed
+        if -pos > rank:
             return expr
 
         # Add index at the required position
@@ -323,8 +369,11 @@ def insert_index(expr, pos, index_var):
             if indices[pos].start is not None:
                 index_var = PyccelAdd(index_var, indices[pos].start, simplify=True)
 
+        # Update index
         indices[pos] = index_var
-        return IndexedElement(base, *indices)
+
+        # Get new indexed object
+        return get_deep_indexed_element(base, indices)
 
     elif isinstance(expr, PyccelArithmeticOperator):
         return type(expr)(insert_index(expr.args[0], pos, index_var),
@@ -620,21 +669,27 @@ def collect_loops(block, indices, new_index, language_has_vectors = False, resul
 
 def insert_fors(blocks, indices, scope, level = 0):
     """
+    Create For loops as requested by the output of collect_loops.
+
     Run through the output of collect_loops and create For loops of the
-    requested sizes
+    requested sizes.
 
     Parameters
-    ==========
-    block   : list of LoopCollection
-            The result of a call to collect_loops
+    ----------
+    blocks : list of LoopCollection
+        The result of a call to collect_loops.
     indices : list
-            The index variables
-    level   : int
-            The index of the index variable used in the outermost loop
-    Results
-    =======
-    block : list of TypedAstNodes
-            The modified expression
+        The index variables.
+    scope : Scope
+        The scope on which the loop is defined. This is where the scope for
+        the new For loop will be created.
+    level : int, default=0
+        The index of the index variable used in the outermost loop.
+
+    Returns
+    -------
+    list[TypedAstNode]
+        The modified expression.
     """
     if all(not isinstance(b, LoopCollection) for b in blocks.body):
         body = blocks.body
@@ -678,9 +733,9 @@ def expand_inhomog_tuple_assignments(block, language_has_vectors = False):
     >>> from pyccel.ast.literals  import LiteralInteger
     >>> from pyccel.ast.utilities import expand_to_loops
     >>> from pyccel.ast.variable  import Variable
-    >>> a = Variable(PythonNativeInt(), 'a', shape=(,), rank=0)
-    >>> b = Variable(PythonNativeInt(), 'b', shape=(,), rank=0)
-    >>> c = Variable(PythonNativeInt(), 'c', shape=(,), rank=0)
+    >>> a = Variable(PythonNativeInt(), 'a')
+    >>> b = Variable(PythonNativeInt(), 'b')
+    >>> c = Variable(PythonNativeInt(), 'c')
     >>> expr = [Assign(PythonTuple(a,b,c),PythonTuple(LiteralInteger(0),LiteralInteger(1),LiteralInteger(2))]
     >>> expand_inhomog_tuple_assignments(CodeBlock(expr))
     [Assign(a, LiteralInteger(0)), Assign(b, LiteralInteger(1)), Assign(c, LiteralInteger(2))]
@@ -713,33 +768,46 @@ def expand_inhomog_tuple_assignments(block, language_has_vectors = False):
 #==============================================================================
 def expand_to_loops(block, new_index, scope, language_has_vectors = False):
     """
-    Re-write a list of expressions to include explicit loops where necessary
+    Re-write a list of expressions to include explicit loops where necessary.
+
+    Re-write a list of expressions to include explicit loops where necessary.
+    The provided expression is the Pyccel representation of the user code. It
+    is the output of the semantic stage. The result of this function is the
+    equivalent code where any vector expressions are unrolled into explicit
+    loops. The unrolling is done completely for languages such as C which have
+    no support for vector operations and partially for languages such as
+    Fortran which have support for vector operations on objects of the same
+    shape.
 
     Parameters
-    ==========
-    block          : CodeBlock
-                     The expressions to be modified
-    new_index      : function
-                     A function which provides a new variable from a base name,
-                     avoiding name collisions
+    ----------
+    block : CodeBlock
+        The expressions to be modified.
+    new_index : function
+        A function which provides a new variable from a base name, avoiding
+        name collisions.
+    scope : Scope
+        The scope on which the loop is defined. This is where the scope for
+        the new For loop will be created.
     language_has_vectors : bool
-                     Indicates if the language has support for vector
-                     operations of the same shape
+        Indicates if the language has support for vector operations of the
+        same shape.
 
     Returns
-    =======
-    expr        : list of Ast Nodes
-                The expressions with For loops inserted where necessary
+    -------
+    list[PyccelAstNode]
+        The expressions with `For` loops inserted where necessary.
 
     Examples
     --------
     >>> from pyccel.ast.core import Variable, Assign
+    >>> from pyccel.ast.datatypes import PythonNativeInt
     >>> from pyccel.ast.operators import PyccelAdd
     >>> from pyccel.ast.utilities import expand_to_loops
-    >>> a = Variable('int', 'a', shape=(4,), rank=1)
-    >>> b = Variable('int', 'b', shape=(4,), rank=1)
-    >>> c = Variable('int', 'c', shape=(4,), rank=1)
-    >>> i = Variable('int', 'i')
+    >>> a = Variable(PythonNativeInt(), 'a', shape=(4,))
+    >>> b = Variable(PythonNativeInt(), 'b', shape=(4,))
+    >>> c = Variable(PythonNativeInt(), 'c', shape=(4,))
+    >>> i = Variable(PythonNativeInt(), 'i')
     >>> d = PyccelAdd(a,b)
     >>> expr = [Assign(c,d)]
     >>> expand_to_loops(expr, language_has_vectors = False)
diff --git a/pyccel/ast/variable.py b/pyccel/ast/variable.py
index 3810603113..b465e930da 100644
--- a/pyccel/ast/variable.py
+++ b/pyccel/ast/variable.py
@@ -51,9 +51,6 @@ class Variable(TypedAstNode):
         The name of the variable represented. This can be either a string
         or a dotted name, when using a Class attribute.
 
-    rank : int, default: 0
-        The number of dimensions for an array.
-
     memory_handling : str, default: 'stack'
         'heap' is used for arrays, if we need to allocate memory on the heap.
         'stack' if memory should be allocated on the stack, represents stack arrays and scalars.
@@ -78,10 +75,6 @@ class Variable(TypedAstNode):
     cls_base : class, default: None
         Class base if variable is an object or an object member.
 
-    order : str, default: 'C'
-        Used for arrays. Indicates whether the data is stored in C or Fortran format in memory.
-        See order_docs.md in the developer docs for more details.
-
     is_argument : bool, default: False
         Indicates if object is the argument of a function.
 
@@ -95,19 +88,19 @@ class Variable(TypedAstNode):
 
     Examples
     --------
+    >>> from pyccel.ast.datatypes import PythonNativeInt, PythonNativeFloat
     >>> from pyccel.ast.core import Variable
-    >>> Variable('int', 'n')
+    >>> Variable(PythonNativeInt(), 'n')
     n
     >>> n = 4
-    >>> Variable('float', 'x', rank=2, shape=(n,2), memory_handling='heap')
+    >>> Variable(PythonNativeFloat(), 'x', shape=(n,2), memory_handling='heap')
     x
-    >>> Variable('int', DottedName('matrix', 'n_rows'))
+    >>> Variable(PythonNativeInt(), DottedName('matrix', 'n_rows'))
     matrix.n_rows
     """
-    __slots__ = ('_name', '_alloc_shape', '_memory_handling', '_is_const',
-            '_is_target', '_is_optional', '_allows_negative_indexes',
-            '_cls_base', '_is_argument', '_is_temp',
-            '_rank','_shape','_order','_is_private','_class_type')
+    __slots__ = ('_name', '_alloc_shape', '_memory_handling', '_is_const', '_is_target',
+            '_is_optional', '_allows_negative_indexes', '_cls_base', '_is_argument', '_is_temp',
+            '_shape','_is_private','_class_type')
     _attribute_nodes = ()
 
     def __init__(
@@ -115,7 +108,6 @@ def __init__(
         class_type,
         name,
         *,
-        rank=0,
         memory_handling='stack',
         is_const=False,
         is_target=False,
@@ -123,7 +115,6 @@ def __init__(
         is_private=False,
         shape=None,
         cls_base=None,
-        order=None,
         is_argument=False,
         is_temp =False,
         allows_negative_indexes=False
@@ -171,34 +162,22 @@ def __init__(
         self._allows_negative_indexes = allows_negative_indexes
 
         self._cls_base       = cls_base
-        self._order          = order
         self._is_argument    = is_argument
         self._is_temp        = is_temp
 
         # ------------ TypedAstNode Properties ---------------
         assert isinstance(class_type, PyccelType)
-        assert isinstance(rank, int)
+        rank = class_type.rank
 
         if rank == 0:
             assert shape is None
-            assert order is None
 
         elif shape is None:
             shape = tuple(None for i in range(rank))
-        else:
-            assert len(shape) == rank
-
-        if rank == 1:
-            assert order is None
-        elif rank > 1:
-            assert order in ('C', 'F')
 
         self._alloc_shape = shape
         self._class_type = class_type
-        self._rank  = rank
         self._shape = self.process_shape(shape)
-        if self._rank < 2:
-            self._order = None
 
     def process_shape(self, shape):
         """
@@ -338,12 +317,6 @@ def is_stack_array(self):
         """
         return self.on_stack and self.rank > 0
 
-    @property
-    def is_stack_scalar(self):
-        """ Indicates if the variable is located on stack and is a scalar
-        """
-        return self.on_stack and self.rank == 0
-
     @property
     def cls_base(self):
         """ Class from which the Variable inherits
@@ -427,17 +400,15 @@ def is_ndarray(self):
         """
         User friendly method to check if the variable is a numpy.ndarray.
 
-        User friendly method to check if the variable is an ndarray:
-            1. have a rank > 0
-            2. class type is NumpyNDArrayType
+        User friendly method to check if the variable is an ndarray.
         """
-        return isinstance(self.class_type, NumpyNDArrayType) and self.rank > 0
+        return isinstance(self.class_type, NumpyNDArrayType)
 
     def __str__(self):
         return str(self.name)
 
     def __repr__(self):
-        return f'{type(self).__name__}({self.name}, type={self.class_type})'
+        return f'{type(self).__name__}({self.name}, type={repr(self.class_type)})'
 
     def __eq__(self, other):
         if type(self) is type(other):
@@ -458,8 +429,6 @@ def inspect(self):
         print('>>> Variable')
         print(f'  name               = {self.name}')
         print(f'  type               = {self.class_type}')
-        print(f'  rank               = {self.rank}')
-        print(f'  order              = {self.order}')
         print(f'  memory_handling    = {self.memory_handling}')
         print(f'  shape              = {self.shape}')
         print(f'  cls_base           = {self.cls_base}')
@@ -535,10 +504,8 @@ def __reduce_ex__(self, i):
             self.class_type,
             self.name)
         kwargs = {
-            'rank' : self.rank,
             'memory_handling': self.memory_handling,
             'is_optional':self.is_optional,
-            'order':self.order,
             'cls_base':self.cls_base,
             }
 
@@ -827,7 +794,7 @@ class IndexedElement(TypedAstNode):
     >>> IndexedElement(A, i, j) == A[i, j]
     True
     """
-    __slots__ = ('_label', '_indices','_shape','_rank','_order','_class_type')
+    __slots__ = ('_label', '_indices','_shape','_class_type', '_is_slice')
     _attribute_nodes = ('_label', '_indices', '_shape')
 
     def __init__(self, base, *indices):
@@ -843,7 +810,8 @@ def __init__(self, base, *indices):
             return
 
         shape = base.shape
-        rank  = base.rank
+        rank  = base.class_type.container_rank
+        assert len(indices) <= rank
 
         if any(not isinstance(a, (int, TypedAstNode, Slice, LiteralEllipsis)) for a in indices):
             errors.report("Index is not of valid type",
@@ -881,17 +849,23 @@ def __init__(self, base, *indices):
 
                     _shape = MathCeil(PyccelDiv(_shape, step, simplify=True))
                 new_shape.append(_shape)
-        self._rank  = len(new_shape)
-        self._shape = None if self._rank == 0 else tuple(new_shape)
-
-        base_type = base.class_type
-        rank = base.rank
-        for _ in range(base.rank-self._rank):
-            rank -= 1
-            if not (rank and isinstance(base_type, NumpyNDArrayType)):
-                base_type = base_type.element_type
-        self._class_type = base_type
-        self._order = None if self.rank < 2 else base.order
+        new_rank = len(new_shape)
+
+        if new_rank == 0:
+            self._class_type = base.class_type.element_type
+            self._is_slice = False
+            if self._class_type.rank:
+                self._shape = base.shape[rank:]
+            else:
+                self._shape = None
+        elif new_rank != rank:
+            self._class_type = base.class_type.switch_rank(new_rank)
+            self._is_slice = True
+            self._shape = tuple(new_shape) + base.shape[rank:]
+        else:
+            self._class_type = base.class_type
+            self._is_slice = True
+            self._shape = tuple(new_shape) + base.shape[rank:]
 
         super().__init__()
 
@@ -917,27 +891,34 @@ def __repr__(self):
 
     def __getitem__(self, *args):
 
+        if self.class_type.container_rank < len(args):
+            raise IndexError('Rank mismatch.')
+
         if len(args) == 1 and isinstance(args[0], (tuple, list)):
             args = args[0]
 
-        if self.rank < len(args):
-            raise IndexError('Rank mismatch.')
-
-        new_indexes = []
-        j = 0
-        base = self.base
-        for i in self.indices:
-            if isinstance(i, Slice) and j<len(args):
-                if i.step == 1 or i.step is None:
-                    incr = args[j]
-                else:
-                    incr = PyccelMul(i.step, args[j], simplify = True)
-                if i.start != 0 and i.start is not None:
-                    incr = PyccelAdd(i.start, incr, simplify = True)
-                i = incr
-                j += 1
-            new_indexes.append(i)
-        return IndexedElement(base, *new_indexes)
+        if self._is_slice:
+            new_indexes = []
+            j = 0
+            base = self.base
+            for i in self.indices:
+                if isinstance(i, Slice) and j<len(args):
+                    current_arg = args[j]
+                    if isinstance(current_arg, Slice):
+                        raise NotImplementedError("Can't extract a slice from a slice")
+                    else:
+                        if i.step == 1 or i.step is None:
+                            incr = current_arg
+                        else:
+                            incr = PyccelMul(i.step, current_arg, simplify = True)
+                        if i.start != 0 and i.start is not None:
+                            incr = PyccelAdd(i.start, incr, simplify = True)
+                    i = incr
+                    j += 1
+                new_indexes.append(i)
+            return IndexedElement(base, *new_indexes)
+        else:
+            return IndexedElement(self, *args)
 
     @property
     def is_const(self):
@@ -949,6 +930,15 @@ def is_const(self):
         """
         return self.base.is_const
 
+    @property
+    def allows_negative_indexes(self):
+        """
+        Indicate whether variables used to index this Variable can be negative.
+
+        Indicate whether variables used to index this Variable can be negative.
+        """
+        return self.base.allows_negative_indexes
+
 class DottedVariable(Variable):
     """
     Class representing a dotted variable.
diff --git a/pyccel/codegen/printing/ccode.py b/pyccel/codegen/printing/ccode.py
index a8b277e1fc..bc29c4118b 100644
--- a/pyccel/codegen/printing/ccode.py
+++ b/pyccel/codegen/printing/ccode.py
@@ -420,8 +420,9 @@ def copy_NumpyArray_Data(self, expr):
         if order == "F":
             # If the order is F then the data should be copied non-contiguously so a temporary
             # variable is required to pass to array_copy_data
-            temp_var = self.scope.get_temporary_variable(lhs, order='C')
-            operations += self._print(Allocate(temp_var, shape=lhs.shape, order="C", status="unallocated"))
+            new_dtype = lhs.class_type.swap_order()
+            temp_var = self.scope.get_temporary_variable(lhs, class_type=new_dtype)
+            operations += self._print(Allocate(temp_var, shape=lhs.shape, status="unallocated"))
             copy_to = temp_var
         else:
             copy_to = lhs
@@ -1098,7 +1099,7 @@ def formatted_args_to_printf(args_format, args, end):
                 args_format.append(CStringExpression('(', tmp_arg_format_list, ')'))
                 assign = Assign(tmp_list, f)
                 self._additional_code += self._print(assign)
-            elif f.rank > 0:
+            elif f.rank > 0 and not isinstance(f.class_type, StringType):
                 if args_format:
                     code += formatted_args_to_printf(args_format, args, sep)
                     args_format = []
@@ -1363,7 +1364,20 @@ def _print_IndexedElement(self, expr):
         base = expr.base
         inds = list(expr.indices)
         base_shape = base.shape
-        allow_negative_indexes = True if isinstance(base, PythonTuple) else base.allows_negative_indexes
+        allow_negative_indexes = expr.allows_negative_indexes
+        if isinstance(base.class_type, NumpyNDArrayType):
+            #set dtype to the C struct types
+            dtype = self.find_in_ndarray_type_registry(expr.dtype)
+        elif isinstance(base.class_type, HomogeneousContainerType):
+            dtype = self.find_in_ndarray_type_registry(numpy_precision_map[(expr.dtype.primitive_type, expr.dtype.precision)])
+        else:
+            raise NotImplementedError(f"Don't know how to index {expr.class_type} type")
+
+        if isinstance(base, IndexedElement):
+            while isinstance(base, IndexedElement) and isinstance(base.class_type, HomogeneousContainerType):
+                inds = list(base.indices) + inds
+                base = base.base
+
         for i, ind in enumerate(inds):
             if isinstance(ind, PyccelUnarySub) and isinstance(ind.args[0], LiteralInteger):
                 inds[i] = PyccelMinus(base_shape[i], ind.args[0], simplify = True)
@@ -1375,30 +1389,21 @@ def _print_IndexedElement(self, expr):
                         not isinstance(ind, LiteralInteger) and not isinstance(ind, Slice):
                     inds[i] = IfTernaryOperator(PyccelLt(ind, LiteralInteger(0)),
                         PyccelAdd(base_shape[i], ind, simplify = True), ind)
-        if isinstance(base.class_type, NumpyNDArrayType):
-            #set dtype to the C struct types
-            dtype = self.find_in_ndarray_type_registry(expr.dtype)
-        elif isinstance(base.class_type, HomogeneousContainerType):
-            dtype = self.find_in_ndarray_type_registry(numpy_precision_map[(expr.dtype.primitive_type, expr.dtype.precision)])
-        else:
-            raise NotImplementedError(f"Don't know how to index {expr.class_type} type")
+
         base_name = self._print(base)
-        if getattr(base, 'is_ndarray', False) or isinstance(base.class_type, HomogeneousContainerType):
-            if expr.rank > 0:
-                #managing the Slice input
-                for i , ind in enumerate(inds):
-                    if isinstance(ind, Slice):
-                        inds[i] = self._new_slice_with_processed_arguments(ind, PyccelArrayShapeElement(base, i),
-                            allow_negative_indexes)
-                    else:
-                        inds[i] = Slice(ind, PyccelAdd(ind, LiteralInteger(1), simplify = True), LiteralInteger(1),
-                            Slice.Element)
-                inds = [self._print(i) for i in inds]
-                return "array_slicing(%s, %s, %s)" % (base_name, expr.rank, ", ".join(inds))
-            inds = [self._cast_to(i, NumpyInt64Type()).format(self._print(i)) for i in inds]
-        else:
-            raise NotImplementedError(expr)
-        return "GET_ELEMENT(%s, %s, %s)" % (base_name, dtype, ", ".join(inds))
+        if expr.rank > 0:
+            #managing the Slice input
+            for i , ind in enumerate(inds):
+                if isinstance(ind, Slice):
+                    inds[i] = self._new_slice_with_processed_arguments(ind, PyccelArrayShapeElement(base, i),
+                        allow_negative_indexes)
+                else:
+                    inds[i] = Slice(ind, PyccelAdd(ind, LiteralInteger(1), simplify = True), LiteralInteger(1),
+                        Slice.Element)
+            indices = ", ".join(self._print(i) for i in inds)
+            return f"array_slicing({base_name}, {expr.rank}, {indices})"
+        indices = ", ".join(self._cast_to(i, NumpyInt64Type()).format(self._print(i)) for i in inds)
+        return f"GET_ELEMENT({base_name}, {dtype}, {indices})"
 
 
     def _cast_to(self, expr, dtype):
@@ -1539,7 +1544,7 @@ def _print_Allocate(self, expr):
             is_view = 'false' if variable.on_heap else 'true'
             order = "order_f" if expr.order == "F" else "order_c"
             alloc_code = f"{self._print(variable)} = array_create({variable.rank}, {shape_Assign}, {dtype}, {is_view}, {order});\n"
-            return '{}{}'.format(free_code, alloc_code)
+            return f'{free_code}{alloc_code}'
         elif variable.is_alias:
             var_code = self._print(ObjectAddress(variable))
             if expr.like:
diff --git a/pyccel/codegen/printing/fcode.py b/pyccel/codegen/printing/fcode.py
index dc86daf7db..af9d9fd598 100644
--- a/pyccel/codegen/printing/fcode.py
+++ b/pyccel/codegen/printing/fcode.py
@@ -21,7 +21,7 @@
 
 from pyccel.ast.builtins import PythonInt, PythonType, PythonPrint, PythonRange
 from pyccel.ast.builtins import PythonTuple, DtypePrecisionToCastFunction
-from pyccel.ast.builtins import PythonBool, PythonAbs
+from pyccel.ast.builtins import PythonBool
 
 from pyccel.ast.core import FunctionDef
 from pyccel.ast.core import SeparatorComment, Comment
@@ -34,9 +34,9 @@
 from pyccel.ast.core import FunctionCall, PyccelFunctionDef
 
 from pyccel.ast.datatypes import PrimitiveBooleanType, PrimitiveIntegerType, PrimitiveFloatingPointType, PrimitiveComplexType
-from pyccel.ast.datatypes import SymbolicType, StringType, FixedSizeNumericType
+from pyccel.ast.datatypes import SymbolicType, StringType, FixedSizeNumericType, HomogeneousContainerType
 from pyccel.ast.datatypes import PythonNativeInt
-from pyccel.ast.datatypes import CustomDataType, InhomogeneousTupleType
+from pyccel.ast.datatypes import CustomDataType, InhomogeneousTupleType, TupleType
 from pyccel.ast.datatypes import pyccel_type_to_original_type
 
 from pyccel.ast.internals import Slice, PrecomputedCode, PyccelArrayShapeElement
@@ -52,8 +52,8 @@
 from pyccel.ast.numpyext import NumpyEmpty, NumpyInt32
 from pyccel.ast.numpyext import NumpyFloat, NumpyBool
 from pyccel.ast.numpyext import NumpyReal, NumpyImag
-from pyccel.ast.numpyext import NumpyRand
-from pyccel.ast.numpyext import NumpyNewArray
+from pyccel.ast.numpyext import NumpyRand, NumpyAbs
+from pyccel.ast.numpyext import NumpyNewArray, NumpyArray
 from pyccel.ast.numpyext import NumpyNonZero
 from pyccel.ast.numpyext import NumpySign
 from pyccel.ast.numpyext import NumpyIsFinite, NumpyIsNan
@@ -733,7 +733,8 @@ def _print_PythonPrint(self, expr):
             elif isinstance(f, PythonType):
                 args_format.append('A')
                 args.append(self._print(f.print_string))
-            elif f.rank > 0 and not isinstance(f, FunctionCall):
+            elif isinstance(f.class_type, (TupleType, HomogeneousContainerType)) and not isinstance(f.class_type, StringType) \
+                    and not isinstance(f, FunctionCall):
                 if args_format:
                     code += self._formatted_args_to_print(args_format, args, sep, separator, expr)
                     args_format = []
@@ -1043,7 +1044,7 @@ def _print_NumpyEmpty(self, expr):
 
     def _print_NumpyNorm(self, expr):
         """Fortran print."""
-        arg = PythonAbs(expr.arg) if isinstance(expr.arg.dtype.primitive_type, PrimitiveComplexType) else expr.arg
+        arg = NumpyAbs(expr.arg) if isinstance(expr.arg.dtype.primitive_type, PrimitiveComplexType) else expr.arg
         if expr.axis:
             axis = expr.axis
             if arg.order != 'F':
@@ -1196,10 +1197,10 @@ def _print_NumpyArray(self, expr):
             cast_func = DtypePrecisionToCastFunction[expr.dtype]
         except KeyError:
             errors.report(PYCCEL_RESTRICTION_TODO, severity='fatal')
-        arg = expr.arg if expr.arg.dtype == expr.dtype else cast_func(expr.arg)
         # If Numpy array is stored with column-major ordering, transpose values
         # use reshape with order for rank > 2
         if expr.rank <= 2:
+            arg = expr.arg if expr.arg.dtype == expr.dtype else cast_func(expr.arg)
             rhs_code = self._print(arg)
             if expr.arg.order and expr.arg.order != expr.order:
                 rhs_code = f'transpose({rhs_code})'
@@ -1216,8 +1217,15 @@ def _print_NumpyArray(self, expr):
             inv_order = 'C' if order == 'F' else 'F'
             for a in expr_args:
                 ac = self._print(a)
+
+                # Pack list/tuple of array/list/tuple into array
+                if a.order is None and a.rank > 1:
+                    a = NumpyArray(a)
+                    ac = self._print(a)
+
+                # Reshape array element if out of order
                 if a.order == inv_order:
-                    shape = a.shape if a.order == 'C' else a.shape[::-1]
+                    shape = a.shape[::-1] if a.order == 'F' else a.shape
                     shape_code = ', '.join(self._print(i) for i in shape)
                     order_code = ', '.join(self._print(LiteralInteger(i)) for i in range(a.rank, 0, -1))
                     ac = f'reshape({ac}, [{shape_code}], order=[{order_code}])'
@@ -1341,8 +1349,7 @@ def _print_NumpyRand(self, expr):
         if (not self._additional_code):
             self._additional_code = ''
         var = self.scope.get_temporary_variable(expr.dtype, memory_handling = 'stack',
-                shape = expr.shape,
-                order = expr.order, rank = expr.rank)
+                shape = expr.shape)
 
         self._additional_code = self._additional_code + self._print(Assign(var,expr)) + '\n'
         return self._print(var)
@@ -1537,9 +1544,9 @@ def _print_Declare(self, expr):
         vstr = self._print(expr.variable.name)
 
         # arrays are 0-based in pyccel, to avoid ambiguity with range
-        s = self._print(LiteralInteger(0))
+        start_val = self._print(LiteralInteger(0))
         if not(is_static) and (on_heap or (var.shape is None)):
-            s = ''
+            start_val = ''
 
         # Default empty strings
         intentstr      = ''
@@ -1548,6 +1555,7 @@ def _print_Declare(self, expr):
         privatestr     = ''
         rankstr        = ''
         externalstr    = ''
+        is_string = isinstance(var.class_type, StringType)
 
         # Compute intent string
         if intent:
@@ -1557,10 +1565,10 @@ def _print_Declare(self, expr):
                 if is_const:
                     intentstr += ', intent(in)'
             else:
-                intentstr = ', intent({})'.format(intent)
+                intentstr = f', intent({intent})'
 
         # Compute allocatable string
-        if not is_static:
+        if not is_static and not is_string:
             if is_alias:
                 allocatablestr = ', pointer'
 
@@ -1569,7 +1577,7 @@ def _print_Declare(self, expr):
 
             # ISSUES #177: var is allocatable and target
             if is_target:
-                allocatablestr = '{}, target'.format(allocatablestr)
+                allocatablestr = f'{allocatablestr}, target'
 
         # Compute optional string
         if is_optional:
@@ -1585,29 +1593,29 @@ def _print_Declare(self, expr):
 
         # Compute rank string
         # TODO: improve
-        if ((rank == 1) and (isinstance(shape, (int, TypedAstNode))) and (is_static or on_stack)):
-            rankstr = '({0}:{1})'.format(self._print(s), self._print(PyccelMinus(shape, LiteralInteger(1), simplify = True)))
+        if not is_string:
+            if ((rank == 1) and (isinstance(shape, (int, TypedAstNode))) and (is_static or on_stack)):
+                ubound = PyccelMinus(shape, LiteralInteger(1), simplify = True)
+                rankstr = f'({self._print(start_val)}:{self._print(ubound)})'
 
-        elif ((rank > 0) and (isinstance(shape, (PythonTuple, tuple))) and (is_static or on_stack)):
-            #TODO fix bug when we include shape of type list
+            elif ((rank > 0) and (isinstance(shape, (PythonTuple, tuple))) and (is_static or on_stack)):
+                #TODO fix bug when we include shape of type list
 
-            if var.order == 'C':
-                rankstr = ','.join('{0}:{1}'.format(self._print(s),
-                                                      self._print(PyccelMinus(i, LiteralInteger(1), simplify = True))) for i in shape[::-1])
-            else:
-                rankstr =  ','.join('{0}:{1}'.format(self._print(s),
-                                                     self._print(PyccelMinus(i, LiteralInteger(1), simplify = True))) for i in shape)
-            rankstr = '({rank})'.format(rank=rankstr)
+                ordered_shape = shape[::-1] if var.order == 'C' else shape
+                ubounds = [PyccelMinus(s, LiteralInteger(1), simplify = True) for s in ordered_shape]
+                rankstr = ','.join(f'{self._print(start_val)}:{self._print(u)}' for u in ubounds)
+                rankstr = f'({rankstr})'
 
-        elif (rank > 0) and on_heap and intent_in:
-            rankstr = '({})'.format(','.join(['0:'] * rank))
+            elif (rank > 0) and on_heap and intent_in:
+                rankstr = ','.join(['0:'] * rank)
+                rankstr = f'({rankstr})'
 
-        elif (rank > 0) and (on_heap or is_alias):
-            rankstr = '({})'.format(','.join( [':'] * rank))
+            elif (rank > 0) and (on_heap or is_alias):
+                rankstr = "(:" + ",:" * (rank-1) + ")"
 
-#        else:
-#            errors.report(PYCCEL_RESTRICTION_TODO, symbol=expr,
-#                severity='fatal')
+    #        else:
+    #            errors.report(PYCCEL_RESTRICTION_TODO, symbol=expr,
+    #                severity='fatal')
 
         mod_str = ''
         if expr.module_variable and not is_private and isinstance(expr.variable.class_type, FixedSizeNumericType):
@@ -1616,7 +1624,7 @@ def _print_Declare(self, expr):
         # Construct declaration
         left  = dtype + allocatablestr + optionalstr + privatestr + externalstr + mod_str + intentstr
         right = vstr + rankstr + code_value
-        return '{} :: {}\n'.format(left, right)
+        return f'{left} :: {right}\n'
 
     def _print_AliasAssign(self, expr):
         code = ''
@@ -2887,13 +2895,25 @@ def _print_LiteralComplex(self, expr):
 
     def _print_IndexedElement(self, expr):
         base = expr.base
-        base_code = self._print(base)
 
         inds = list(expr.indices)
         if len(inds) == 1 and isinstance(inds[0], LiteralEllipsis):
             inds = [Slice(None,None)]*expr.rank
 
-        if expr.base.order == 'C':
+        # Condense all indices on homogeneous objects into one IndexedElement for printing
+        # This should be removed when support for lists is added
+        if isinstance(base, IndexedElement):
+            while isinstance(base, IndexedElement) and isinstance(base.class_type, HomogeneousContainerType):
+                inds = list(base.indices) + inds
+                base = base.base
+
+        rank = base.rank
+        if len(inds)<rank:
+            inds += [Slice(None,None)]*(rank-base.class_type.container_rank)
+
+        base_code = self._print(base)
+
+        if base.order != 'F':
             inds = inds[::-1]
         allow_negative_indexes = base.allows_negative_indexes
 
diff --git a/pyccel/codegen/printing/pycode.py b/pyccel/codegen/printing/pycode.py
index 628ef21e05..79729b398e 100644
--- a/pyccel/codegen/printing/pycode.py
+++ b/pyccel/codegen/printing/pycode.py
@@ -245,18 +245,7 @@ def _print_FunctionDefArgument(self, expr):
             type_annotation = f"'{self._print(expr.annotation)}'"
         else:
             var = expr.var
-            def get_type_annotation(var):
-                if isinstance(var, Variable):
-                    type_annotation = str(var.class_type)
-                    if var.rank:
-                        type_annotation += '[' + ','.join(':' for _ in range(var.rank)) + ']'
-                    return f"'{type_annotation}'"
-                elif isinstance(var, FunctionAddress):
-                    results = ', '.join(get_type_annotation(r.var) for r in var.results)
-                    arguments = ', '.join(get_type_annotation(a.var) for a in var.arguments)
-                    return f'"({results})({arguments})"'
-
-            type_annotation = get_type_annotation(var)
+            type_annotation = f"'{var.class_type}'"
 
         if expr.has_default:
             if isinstance(expr.value, FunctionDef):
@@ -282,11 +271,7 @@ def _print_IndexedElement(self, expr):
             if len(indices) == 1 and isinstance(indices[0], (tuple, list)):
                 indices = indices[0]
 
-            indices = [self._print(i) for i in indices]
-            if expr.pyccel_staging != 'syntactic' and isinstance(expr.base.class_type, (HomogeneousTupleType, HomogeneousListType)):
-                indices = ']['.join(i for i in indices)
-            else:
-                indices = ','.join(i for i in indices)
+            indices = ','.join(self._print(i) for i in indices)
         else:
             errors.report(PYCCEL_RESTRICTION_TODO, symbol=expr,
                 severity='fatal')
diff --git a/pyccel/codegen/wrapper/c_to_python_wrapper.py b/pyccel/codegen/wrapper/c_to_python_wrapper.py
index 56dc3822e9..460420f686 100644
--- a/pyccel/codegen/wrapper/c_to_python_wrapper.py
+++ b/pyccel/codegen/wrapper/c_to_python_wrapper.py
@@ -303,7 +303,7 @@ def _get_check_function(self, py_obj, arg, raise_error):
             func_call = FunctionCall(check_func, [py_obj, type_ref, LiteralInteger(rank), flag])
 
         if raise_error:
-            message = LiteralString(f"Expected an argument of type {dtype} for argument {arg.name}")
+            message = LiteralString(f"Expected an argument of type {arg.class_type} for argument {arg.name}")
             python_error = FunctionCall(PyErr_SetString, [PyTypeError, message])
             error_code = (python_error,)
 
@@ -611,7 +611,7 @@ def _build_module_init_function(self, expr, imports):
         # Create necessary variables
         module_var = self.get_new_PyObject("mod")
         API_var_name = self.scope.get_new_name(f'Py{mod_name}_API')
-        API_var = Variable(CStackArray(BindCPointer()), API_var_name, rank=1, shape = (n_classes,),
+        API_var = Variable(CStackArray(BindCPointer()), API_var_name, shape = (n_classes,),
                                     cls_base = StackArrayClass)
         self.scope.insert_variable(API_var)
         capsule_obj = self.get_new_PyObject(self.scope.get_new_name('c_api_object'))
@@ -708,7 +708,7 @@ def _build_module_import_function(self, expr):
         func_name = self.scope.get_new_name(f'{mod_name}_import')
 
         API_var_name = self.scope.get_new_name(f'Py{mod_name}_API')
-        API_var = Variable(CStackArray(BindCPointer()), API_var_name, rank=1, shape = (None,),
+        API_var = Variable(CStackArray(BindCPointer()), API_var_name, shape = (None,),
                                     cls_base = StackArrayClass,
                                     memory_handling = 'alias')
         self.scope.insert_variable(API_var)
@@ -772,7 +772,7 @@ def _allocate_class_instance(self, class_var, scope, is_alias):
 
         alias_val = LiteralTrue() if is_alias else LiteralFalse()
 
-        return [Allocate(class_var, shape=(), order=None, status='unallocated'),
+        return [Allocate(class_var, shape=(), status='unallocated'),
                 AliasAssign(ref_list, FunctionCall(PyList_New, ())),
                 Assign(alias_bool, alias_val)]
 
@@ -824,7 +824,7 @@ def _get_class_allocator(self, class_dtype, func = None):
         else:
             result_name = self.scope.get_new_name('result')
             result = Variable(class_dtype, result_name)
-            body.append(Allocate(c_res, shape=(), order=None, status='unallocated',
+            body.append(Allocate(c_res, shape=(), status='unallocated',
                          like = result))
 
         body.append(Return([ObjectAddress(PointerCast(python_result_var, func_results[0].var))]))
@@ -1299,7 +1299,7 @@ def _wrap_FunctionDef(self, expr):
             if isinstance(r, DottedVariable):
                 self.scope.remove_variable(r, name=n)
                 if not o_r.var.is_alias:
-                    body.append(Allocate(r, shape=(), order=None, status='unallocated', like=o_r.var))
+                    body.append(Allocate(r, shape=(), status='unallocated', like=o_r.var))
         c_results = [ObjectAddress(r) if r.dtype is BindCPointer() else r for r in c_results]
         c_results = [PointerCast(r, cast_type = o_r.var) if isinstance(r, DottedVariable) else r for r,o_r in zip(c_results, original_c_results)]
 
@@ -1655,7 +1655,7 @@ def _wrap_BindCFunctionDefResult(self, expr):
             # Save so we can find by iterating over func.bind_c_results
             self.scope.insert_variable(c_res, orig_var_name)
 
-            body.append(Allocate(c_res, shape = shape_vars, order = orig_var.order, status='unallocated'))
+            body.append(Allocate(c_res, shape = shape_vars, status='unallocated'))
             body.append(AliasAssign(DottedVariable(VoidType(), 'raw_data', memory_handling = 'alias', lhs=c_res), arg_var))
         elif isinstance(orig_var.dtype, CustomDataType):
             c_res = expr.var.clone(self.scope.get_expected_name(var_name), is_argument = False, memory_handling='alias')
@@ -1752,7 +1752,7 @@ def _wrap_BindCArrayVariable(self, expr):
         # Create ndarray to store array data
         nd_var = self.scope.get_temporary_variable(dtype_or_var = v,
                 name = v.name, memory_handling = 'alias')
-        alloc = Allocate(nd_var, shape=shape, order=nd_var.order, status='unallocated')
+        alloc = Allocate(nd_var, shape=shape, status='unallocated')
         # Save raw_data into ndarray to obtain useable pointer
         set_data = AliasAssign(DottedVariable(VoidType(), 'raw_data',
                 memory_handling = 'alias', lhs=nd_var), var)
@@ -1817,7 +1817,7 @@ def _wrap_DottedVariable(self, expr):
             body = [AliasAssign(new_res_val, attrib), *res_wrapper]
         elif isinstance(expr.dtype, CustomDataType):
             self.scope.remove_variable(new_res_val, name = expr.name)
-            body = [Allocate(getter_result, shape=(), order=None, status='unallocated'),
+            body = [Allocate(getter_result, shape=(), status='unallocated'),
                     AliasAssign(new_res_val, attrib),
                     *res_wrapper]
         else:
@@ -1953,7 +1953,7 @@ def _wrap_BindCClassProperty(self, expr):
             call = Assign(c_results, FunctionCall(expr.getter, (class_obj,)))
 
         if isinstance(getter_result.dtype, CustomDataType):
-            arg_code.append(Allocate(getter_result, shape=(), order=None, status='unallocated'))
+            arg_code.append(Allocate(getter_result, shape=(), status='unallocated'))
 
         wrapped_var = expr.getter.original_function
         res_wrapper.extend(self._incref_return_pointer(getter_args[0], getter_result, wrapped_var))
diff --git a/pyccel/codegen/wrapper/fortran_to_c_wrapper.py b/pyccel/codegen/wrapper/fortran_to_c_wrapper.py
index 125796efd1..bd8050a274 100644
--- a/pyccel/codegen/wrapper/fortran_to_c_wrapper.py
+++ b/pyccel/codegen/wrapper/fortran_to_c_wrapper.py
@@ -414,8 +414,7 @@ def _wrap_FunctionDefResult(self, expr):
                 scope.insert_variable(ptr_var)
 
                 # Define the additional steps necessary to define and fill ptr_var
-                alloc = Allocate(ptr_var, shape=result.shape,
-                                 order=var.order, status='unallocated')
+                alloc = Allocate(ptr_var, shape=result.shape, status='unallocated')
                 copy = Assign(ptr_var, local_var)
                 self._additional_exprs.extend([alloc, copy])
             else:
@@ -608,7 +607,7 @@ def _wrap_ClassDef(self, expr):
         result = BindCFunctionDefResult(bind_var, local_var, func_scope)
 
         # Define the additional steps necessary to define and fill ptr_var
-        alloc = Allocate(local_var, shape=(), order=None, status='unallocated')
+        alloc = Allocate(local_var, shape=(), status='unallocated')
         c_loc = CLocFunc(local_var, bind_var)
         body = [alloc, c_loc]
 
diff --git a/pyccel/errors/messages.py b/pyccel/errors/messages.py
index 944b128751..5be9064878 100644
--- a/pyccel/errors/messages.py
+++ b/pyccel/errors/messages.py
@@ -77,7 +77,7 @@
 INCOMPATIBLE_TYPEVAR_TO_FUNC = 'TypeError: ufunc "{}" not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule \'safe\''
 UNSUPPORTED_ARGUMENT_2_FOR_SUPER = 'Unsupported argument 2 for "super"'
 WRONG_NUMBER_OUTPUT_ARGS = 'Number of output arguments does not match number of provided variables'
-INDEXED_TUPLE = 'Tuples must be indexed with constant integers for the type inference to work'
+INDEXED_TUPLE = 'Inhomogeneous tuples must be indexed with constant integers for the type inference to work'
 LIST_OF_TUPLES = 'Cannot create list of non-homogeneous tuples'
 
 UNUSED_DECORATORS = 'Decorator(s) not used'
diff --git a/pyccel/parser/semantic.py b/pyccel/parser/semantic.py
index 688d61a01f..de71c02823 100644
--- a/pyccel/parser/semantic.py
+++ b/pyccel/parser/semantic.py
@@ -10,8 +10,6 @@
 from itertools import chain, product
 import warnings
 
-from math import prod
-
 from sympy.utilities.iterables import iterable as sympy_iterable
 
 from sympy import Sum as Summation
@@ -123,11 +121,11 @@
 
 from pyccel.errors.messages import (PYCCEL_RESTRICTION_TODO, UNDERSCORE_NOT_A_THROWAWAY,
         UNDEFINED_VARIABLE, IMPORTING_EXISTING_IDENTIFIED, INDEXED_TUPLE, LIST_OF_TUPLES,
-        INVALID_INDICES, INCOMPATIBLE_ARGUMENT, INCOMPATIBLE_ORDERING,
+        INVALID_INDICES, INCOMPATIBLE_ARGUMENT,
         UNRECOGNISED_FUNCTION_CALL, STACK_ARRAY_SHAPE_UNPURE_FUNC, STACK_ARRAY_UNKNOWN_SHAPE,
         ARRAY_DEFINITION_IN_LOOP, STACK_ARRAY_DEFINITION_IN_LOOP, MISSING_TYPE_ANNOTATIONS,
         INCOMPATIBLE_TYPES_IN_ASSIGNMENT, ARRAY_ALREADY_IN_USE, ASSIGN_ARRAYS_ONE_ANOTHER,
-        INVALID_POINTER_REASSIGN, INCOMPATIBLE_REDEFINITION, ARRAY_IS_ARG,
+        INVALID_POINTER_REASSIGN, ARRAY_IS_ARG,
         INCOMPATIBLE_REDEFINITION_STACK_ARRAY, ARRAY_REALLOCATION, RECURSIVE_RESULTS_REQUIRED,
         PYCCEL_RESTRICTION_INHOMOG_LIST, UNDEFINED_IMPORT_OBJECT, UNDEFINED_LAMBDA_VARIABLE,
         UNDEFINED_LAMBDA_FUNCTION, UNDEFINED_INIT_METHOD, UNDEFINED_FUNCTION,
@@ -135,7 +133,7 @@
         PYCCEL_RESTRICTION_LIST_COMPREHENSION_LIMITS, PYCCEL_RESTRICTION_LIST_COMPREHENSION_SIZE,
         UNUSED_DECORATORS, UNSUPPORTED_POINTER_RETURN_VALUE, PYCCEL_RESTRICTION_OPTIONAL_NONE,
         PYCCEL_RESTRICTION_PRIMITIVE_IMMUTABLE, PYCCEL_RESTRICTION_IS_ISNOT,
-        FOUND_DUPLICATED_IMPORT, UNDEFINED_WITH_ACCESS, MACRO_MISSING_HEADER_OR_FUNC, PYCCEL_RESTRICTION_INHOMOG_SET)
+        FOUND_DUPLICATED_IMPORT, UNDEFINED_WITH_ACCESS, MACRO_MISSING_HEADER_OR_FUNC)
 
 from pyccel.parser.base      import BasicParser
 from pyccel.parser.syntactic import SyntaxParser
@@ -669,10 +667,8 @@ def _infer_type(self, expr):
 
         Create a dictionary describing all the type information that can be
         inferred about the expression `expr`. This includes information about:
-        - `datatype`
-        - `rank`
+        - `class_type`
         - `shape`
-        - `order`
         - `memory_handling`
         - `cls_base`
         - `is_target`
@@ -690,8 +686,6 @@ def _infer_type(self, expr):
         """
         d_var = {
                 'shape'    : expr.shape,
-                'rank'     : expr.rank,
-                'order'    : expr.order,
                 'class_type' : expr.class_type,
             }
 
@@ -702,6 +696,10 @@ def _infer_type(self, expr):
             d_var['is_target'      ] = expr.is_target
             return d_var
 
+        elif isinstance(expr, LiteralString):
+            d_var['memory_handling'] = 'stack'
+            return d_var
+
         elif isinstance(expr, PythonTuple):
             d_var['cls_base'       ] = TupleClass
             d_var['memory_handling'] = 'heap'
@@ -977,31 +975,6 @@ def _handle_function_args(self, arguments):
                 args.append(a)
         return args
 
-    def get_type_description(self, var, include_rank = True):
-        """
-        Get a text description of the type of a variable.
-
-        Get a text description of the type of a variable.
-        This is notably useful for error messages.
-
-        Parameters
-        ----------
-        var : Variable
-            The variable to describe.
-        include_rank : bool, default=True
-            Indicates whether rank information should be included.
-
-        Returns
-        -------
-        str
-            A description of the variable type.
-        """
-        descr = str(var.class_type)
-        if include_rank and var.rank>0:
-            dims = ','.join(':'*var.rank)
-            descr += f'[{dims}]'
-        return descr
-
     def _check_argument_compatibility(self, input_args, func_args, func, elemental, raise_error=True, error_type='error'):
         """
         Check that the provided arguments match the expected types.
@@ -1033,8 +1006,7 @@ def incompatible(i_arg, f_arg):
                 return i_arg.class_type.datatype != f_arg.class_type.datatype
         else:
             def incompatible(i_arg, f_arg):
-                return (i_arg.class_type != f_arg.class_type or
-                        i_arg.rank != f_arg.rank)
+                return i_arg.class_type != f_arg.class_type
 
         err_msgs = []
         # Compare each set of arguments
@@ -1049,14 +1021,11 @@ def incompatible(i_arg, f_arg):
 
             # Check for compatibility
             if incompatible(i_arg, f_arg):
-                expected  = self.get_type_description(f_arg, not elemental)
-                type_name = self.get_type_description(i_arg, not elemental)
+                expected  = str(f_arg.class_type)
+                type_name = str(i_arg.class_type)
                 received  = f'{i_arg} ({type_name})'
                 err_msgs += [INCOMPATIBLE_ARGUMENT.format(idx+1, received, func, expected)]
 
-            if f_arg.rank > 1 and i_arg.order != f_arg.order:
-                err_msgs += [INCOMPATIBLE_ORDERING.format(idx=idx+1, arg=i_arg, func=func, order=f_arg.order)]
-
         if err_msgs:
             if raise_error:
                 bounding_box=(self.current_ast_node.lineno, self.current_ast_node.col_offset)
@@ -1109,8 +1078,9 @@ def _handle_function(self, expr, func, args, is_method = False):
                 if not message:
                     message = UNRECOGNISED_FUNCTION_CALL
                 errors.report(message,
-                                symbol = expr,
-                                severity = 'fatal')
+                              symbol = expr,
+                              traceback = e.__traceback__,
+                              severity = 'fatal')
 
             return new_expr
         else:
@@ -1409,7 +1379,7 @@ def _assign_lhs_variable(self, lhs, d_var, rhs, new_expressions, is_augassign,ar
             name = lhs
             if lhs == '_':
                 name = self.scope.get_new_name()
-            dtype = d_var.pop('class_type')
+            class_type = d_var.pop('class_type')
 
             d_lhs = d_var.copy()
             # ISSUES #177: lhs must be a pointer when rhs is heap array
@@ -1443,7 +1413,7 @@ def _assign_lhs_variable(self, lhs, d_var, rhs, new_expressions, is_augassign,ar
                         attribute_name = lhs.name[-1]
                         new_name = class_def.scope.get_expected_name(attribute_name)
                         # Create the attribute
-                        member = self._create_variable(new_name, dtype, rhs, d_lhs)
+                        member = self._create_variable(new_name, class_type, rhs, d_lhs)
 
                         # Insert the attribute to the class scope
                         # Passing the original name ensures that the attribute can be found under this name
@@ -1469,7 +1439,7 @@ def _assign_lhs_variable(self, lhs, d_var, rhs, new_expressions, is_augassign,ar
 
                 # We cannot allow the definition of a stack array from a shape which
                 # is unknown at the declaration
-                if d_lhs['rank'] > 0 and d_lhs.get('memory_handling', None) == 'stack':
+                if not isinstance(class_type, StringType) and class_type.rank > 0 and d_lhs.get('memory_handling', None) == 'stack':
                     for a in d_lhs['shape']:
                         if (isinstance(a, FunctionCall) and not a.funcdef.is_pure) or \
                                 any(not f.funcdef.is_pure for f in a.get_attribute_nodes(FunctionCall)):
@@ -1487,7 +1457,7 @@ def _assign_lhs_variable(self, lhs, d_var, rhs, new_expressions, is_augassign,ar
                 if not isinstance(lhs, DottedVariable):
                     new_name = self.scope.get_expected_name(name)
                     # Create new variable
-                    lhs = self._create_variable(new_name, dtype, rhs, d_lhs, arr_in_multirets=arr_in_multirets)
+                    lhs = self._create_variable(new_name, class_type, rhs, d_lhs, arr_in_multirets=arr_in_multirets)
 
                     # Add variable to scope
                     self.scope.insert_variable(lhs, name)
@@ -1496,7 +1466,7 @@ def _assign_lhs_variable(self, lhs, d_var, rhs, new_expressions, is_augassign,ar
                 # Add memory allocation if needed
                 array_declared_in_function = (isinstance(rhs, FunctionCall) and not isinstance(rhs.funcdef, PyccelFunctionDef) \
                                             and not getattr(rhs.funcdef, 'is_elemental', False) and not isinstance(lhs.class_type, HomogeneousTupleType)) or arr_in_multirets
-                if lhs.on_heap and not array_declared_in_function:
+                if not isinstance(lhs.class_type, StringType) and lhs.on_heap and not array_declared_in_function:
                     if self.scope.is_loop:
                         # Array defined in a loop may need reallocation at every cycle
                         errors.report(ARRAY_DEFINITION_IN_LOOP, symbol=name,
@@ -1518,15 +1488,15 @@ def _assign_lhs_variable(self, lhs, d_var, rhs, new_expressions, is_augassign,ar
                                     new_args.extend(v for v in a.get_vars() if v.rank>0)
                                 else:
                                     new_expressions.append(Allocate(a,
-                                        shape=a.alloc_shape, order=a.order, status=status))
+                                        shape=a.alloc_shape, status=status))
                             args = new_args
                     else:
-                        new_expressions.append(Allocate(lhs, shape=lhs.alloc_shape, order=lhs.order, status=status))
+                        new_expressions.append(Allocate(lhs, shape=lhs.alloc_shape, status=status))
                 # ...
 
                 # ...
                 # Add memory deallocation
-                if isinstance(lhs.class_type, CustomDataType) or not lhs.on_stack:
+                if isinstance(lhs.class_type, CustomDataType) or (not lhs.on_stack and not isinstance(lhs.class_type, StringType)):
                     if isinstance(lhs, InhomogeneousTupleVariable):
                         args = [v for v in lhs.get_vars() if v.rank>0]
                         new_args = []
@@ -1550,7 +1520,8 @@ def _assign_lhs_variable(self, lhs, d_var, rhs, new_expressions, is_augassign,ar
 
                 # Not yet supported for arrays: x=y+z, x=b[:]
                 # Because we cannot infer shape of right-hand side yet
-                know_lhs_shape = (lhs.rank == 0) or all(sh is not None for sh in lhs.alloc_shape)
+                know_lhs_shape = (lhs.rank == 0) or all(sh is not None for sh in lhs.alloc_shape) \
+                        or isinstance(lhs.class_type, StringType)
 
                 if not know_lhs_shape:
                     msg = f"Cannot infer shape of right-hand side for expression {lhs} = {rhs}"
@@ -1561,9 +1532,9 @@ def _assign_lhs_variable(self, lhs, d_var, rhs, new_expressions, is_augassign,ar
             # Variable already exists
             else:
 
-                self._ensure_inferred_type_matches_existing(dtype, d_var, var, is_augassign, new_expressions, rhs)
+                self._ensure_inferred_type_matches_existing(class_type, d_var, var, is_augassign, new_expressions, rhs)
 
-                # in the case of elemental, lhs is not of the same dtype as
+                # in the case of elemental, lhs is not of the same class_type as
                 # var.
                 # TODO d_lhs must be consistent with var!
                 # the following is a small fix, since lhs must be already
@@ -1642,25 +1613,20 @@ def _ensure_inferred_type_matches_existing(self, class_type, d_var, var, is_auga
                 tmp_result = PyccelAdd(var, rhs)
                 result_type = tmp_result.class_type
                 raise_error = var.class_type != result_type
-            elif class_type == var.class_type and var.rank == 0:
-                # Don't complain about non-numpy and numpy scalars
-                raise_error = False
             else:
                 raise_error = True
 
             if raise_error:
                 name = var.name
                 rhs_str = str(rhs)
-                errors.report(INCOMPATIBLE_TYPES_IN_ASSIGNMENT.format(var.class_type, class_type),
+                errors.report(INCOMPATIBLE_TYPES_IN_ASSIGNMENT.format(repr(var.class_type), repr(class_type)),
                     symbol=f'{name}={rhs_str}',
                     bounding_box=(self.current_ast_node.lineno, self.current_ast_node.col_offset),
                     severity='error')
 
-        elif not is_augassign:
+        elif not is_augassign and not isinstance(var.class_type, StringType):
 
-            rank  = getattr(var, 'rank' , 'None')
-            order = getattr(var, 'order', 'None')
-            shape = getattr(var, 'shape', 'None')
+            shape = var.shape
 
             # Get previous allocation calls
             previous_allocations = var.get_direct_user_nodes(lambda p: isinstance(p, Allocate))
@@ -1668,18 +1634,7 @@ def _ensure_inferred_type_matches_existing(self, class_type, d_var, var, is_auga
             if len(previous_allocations) == 0:
                 var.set_init_shape(d_var['shape'])
 
-            if (d_var['rank'] != rank) or (rank > 1 and d_var['order'] != order):
-
-                txt = '|{name}| {dtype}{old} <-> {dtype}{new}'
-                def format_shape(s):
-                    return "" if s is None else s
-                txt = txt.format(name=var.name, dtype=dtype, old=format_shape(var.shape),
-                    new=format_shape(d_var['shape']))
-                errors.report(INCOMPATIBLE_REDEFINITION, symbol=txt,
-                    bounding_box=(self.current_ast_node.lineno, self.current_ast_node.col_offset),
-                    severity='error')
-
-            elif d_var['shape'] != shape:
+            if d_var['shape'] != shape:
 
                 if var.is_argument:
                     errors.report(ARRAY_IS_ARG, symbol=var,
@@ -1719,9 +1674,7 @@ def format_shape(s):
                     else:
                         status = 'unallocated'
 
-                    new_expressions.append(Allocate(var,
-                        shape=d_var['shape'], order=d_var['order'],
-                        status=status))
+                    new_expressions.append(Allocate(var, shape=d_var['shape'], status=status))
 
                     if status != 'unallocated':
                         errors.report(ARRAY_REALLOCATION, symbol=var.name,
@@ -1733,9 +1686,7 @@ def format_shape(s):
                 # If previously allocated in If still under construction
                 status = previous_allocations[-1].status
 
-                new_expressions.append(Allocate(var,
-                    shape=d_var['shape'], order=d_var['order'],
-                    status=status))
+                new_expressions.append(Allocate(var, shape=d_var['shape'], status=status))
             elif isinstance(var.class_type, CustomDataType) and not var.is_alias:
                 new_expressions.append(Deallocate(var))
 
@@ -1931,18 +1882,18 @@ def _get_indexed_type(self, base, args, expr):
             return annotation
 
         if all(isinstance(a, Slice) for a in args):
+            rank = len(args)
+            order = None if rank < 2 else 'C'
             if isinstance(base, VariableTypeAnnotation):
                 dtype = base.class_type
-                class_type = NumpyNDArrayType(numpy_process_dtype(dtype))
-                if base.rank != 0:
-                    raise errors.report("Can't index a vector type",
-                            severity='fatal', symbol=expr)
+                if dtype.rank != 0:
+                    raise errors.report("NumPy element must be a scalar type", severity='fatal', symbol=expr)
+                class_type = NumpyNDArrayType(numpy_process_dtype(dtype), rank, order)
             elif isinstance(base, PyccelFunctionDef):
                 dtype_cls = base.cls_name
                 dtype = numpy_process_dtype(dtype_cls.static_type())
-                class_type = NumpyNDArrayType(dtype)
-            rank = len(args)
-            return VariableTypeAnnotation(class_type, rank)
+                class_type = NumpyNDArrayType(dtype, rank, order)
+            return VariableTypeAnnotation(class_type)
 
         if not any(isinstance(a, Slice) for a in args):
             if isinstance(base, PyccelFunctionDef):
@@ -1965,10 +1916,7 @@ def _get_indexed_type(self, base, args, expr):
                     raise errors.report(f"Unknown annotation base {base}\n"+PYCCEL_RESTRICTION_TODO,
                             severity='fatal', symbol=expr)
                 for u in internal_datatypes.type_list:
-                    rank = u.rank+1
-                    order = None if rank == 1 else 'C'
-                    type_annotations.append(VariableTypeAnnotation(class_type(u.class_type),
-                        rank, order, u.is_const))
+                    type_annotations.append(VariableTypeAnnotation(class_type(u.class_type), u.is_const))
                 return UnionTypeAnnotation(*type_annotations)
             else:
                 raise errors.report("Cannot handle non-homogenous type index\n"+PYCCEL_RESTRICTION_TODO,
@@ -2185,15 +2133,13 @@ def _visit_Module(self, expr):
                         for v in headers:
                             types = [self._visit(d).type_list[0] for d in v.dtypes]
                             args = [Variable(t.class_type, PyccelSymbol(f'anon_{i}'),
-                                shape = None, rank = t.rank, order = t.order,
-                                is_const = t.is_const, is_optional = False,
-                                cls_base = t.class_type if t.rank == 0 else NumpyNDArrayType(numpy_process_dtype(t.class_type.element_type)),
+                                shape = None, is_const = t.is_const, is_optional = False,
+                                cls_base = t.class_type,
                                 memory_handling = 'heap' if t.rank > 0 else 'stack') for i,t in enumerate(types)]
 
                             types = [self._visit(d).type_list[0] for d in v.results]
-                            results = [Variable(t.class_type, PyccelSymbol(f'result_{i}'),
-                                shape = None, rank = t.rank, order = t.order,
-                                cls_base = t.class_type if t.rank == 0 else NumpyNDArrayType(numpy_process_dtype(t.class_type.element_type)),
+                            results = [Variable(t.class_type, PyccelSymbol(f'result_{i}'), shape = None,
+                                cls_base = t.class_type,
                                 is_const = t.is_const, is_optional = False,
                                 memory_handling = 'heap' if t.rank > 0 else 'stack') for i,t in enumerate(types)]
 
@@ -2564,13 +2510,12 @@ def _visit_AnnotatedPyccelSymbol(self, expr):
                 address = FunctionAddress(name, args, results)
                 possible_args.append(address)
             elif isinstance(t, VariableTypeAnnotation):
-                rank  = t.rank
                 class_type = t.class_type
                 cls_base = get_cls_base(class_type) or self.scope.find(class_type.name, 'classes')
                 v = var_class(class_type, name, cls_base = cls_base,
-                        shape = None, rank = rank, order = t.order,
+                        shape = None,
                         is_const = t.is_const, is_optional = False,
-                        memory_handling = array_memory_handling if rank > 0 else 'stack',
+                        memory_handling = array_memory_handling if class_type.rank > 0 else 'stack',
                         **kwargs)
                 possible_args.append(v)
             else:
@@ -2595,8 +2540,8 @@ def _visit_SyntacticTypeAnnotation(self, expr):
             class_type = dtype_cls.static_type()
             return UnionTypeAnnotation(VariableTypeAnnotation(class_type))
         elif isinstance(visited_dtype, VariableTypeAnnotation):
-            if visited_dtype.rank > 1:
-                visited_dtype.order = order or visited_dtype.order or 'C'
+            if order and order != visited_dtype.class_type.order:
+                visited_dtype = VariableTypeAnnotation(visited_dtype.class_type.swap_order())
             return UnionTypeAnnotation(visited_dtype)
         elif isinstance(visited_dtype, UnionTypeAnnotation):
             return visited_dtype
@@ -3047,7 +2992,6 @@ def _visit_FunctionCall(self, expr):
             d_var = {'class_type' : dtype,
                     'memory_handling':'stack',
                     'shape' : None,
-                    'rank' : 0,
                     'is_target' : False,
                     'cls_base' : cls_def,
                     }
@@ -3239,9 +3183,7 @@ def _visit_Assign(self, expr):
                     assert(len(c_ranks) == 1)
                     arg = call_args[0].value
                     d_var['shape'          ] = arg.shape
-                    d_var['rank'           ] = arg.rank
                     d_var['memory_handling'] = arg.memory_handling
-                    d_var['order'          ] = arg.order
                     d_var['class_type'     ] = arg.class_type
                     d_var['cls_base'       ] = arg.cls_base
 
@@ -3575,30 +3517,31 @@ def _visit_FunctionalFor(self, expr):
                 stop  = a.stop
                 start = a.start
                 step  = a.step
+
             elif isinstance(a, (PythonZip, PythonEnumerate)):
                 dvar  = self._infer_type(a.element)
                 class_type = dvar.pop('class_type')
-                if dvar['rank'] > 0:
-                    dvar['rank' ] -= 1
+                if class_type.rank > 0:
+                    class_type = class_type.switch_rank(class_type.rank-1)
                     dvar['shape'] = (dvar['shape'])[1:]
-                if dvar['rank'] == 0:
+                if class_type.rank == 0:
+                    dvar['shape'] = None
                     dvar['memory_handling'] = 'stack'
                 var  = Variable(class_type, var, **dvar)
                 stop = a.element.shape[0]
+
             elif isinstance(a, Variable):
                 dvar  = self._infer_type(a)
-                class_type = dvar.pop('class_type').element_type
-                if dvar['rank'] == 1:
-                    dvar['rank']  = 0
-                    dvar['shape'] = None
-                if dvar['rank'] > 1:
-                    dvar['rank'] -= 1
+                class_type = dvar.pop('class_type')
+                if class_type.rank > 0:
+                    class_type = class_type.switch_rank(class_type.rank-1)
                     dvar['shape'] = (dvar['shape'])[1:]
-                if dvar['rank'] == 0:
+                if class_type.rank == 0:
+                    dvar['shape'] = None
                     dvar['memory_handling'] = 'stack'
-
                 var  = Variable(class_type, var, **dvar)
                 stop = a.shape[0]
+
             else:
                 errors.report(PYCCEL_RESTRICTION_TODO,
                               bounding_box=(self.current_ast_node.lineno, self.current_ast_node.col_offset),
@@ -3700,16 +3643,12 @@ def _visit_FunctionalFor(self, expr):
                           severity='fatal')
 
         d_var['memory_handling'] = 'heap'
-        d_var['rank'] += 1
-        shape = [dim]
-        if d_var['rank'] != 1:
-            d_var['order'] = 'C'
-            shape += list(d_var['shape'])
-        else:
-            d_var['order'] = None
-        d_var['shape'] = shape
         class_type = HomogeneousListType(class_type)
         d_var['class_type'] = class_type
+        shape = [dim]
+        if d_var['shape']:
+            shape.extend(d_var['shape'])
+        d_var['shape'] = shape
         d_var['cls_base'] = get_cls_base(class_type)
 
         # ...
diff --git a/pyccel/parser/syntactic.py b/pyccel/parser/syntactic.py
index c39c85ea09..e37e287a4a 100644
--- a/pyccel/parser/syntactic.py
+++ b/pyccel/parser/syntactic.py
@@ -1043,17 +1043,10 @@ def _visit_ClassDef(self, stmt):
 
     def _visit_Subscript(self, stmt):
 
-        ch = stmt
-        args = []
-        while isinstance(ch, ast.Subscript):
-            val = self._visit(ch.slice)
-            if isinstance(val, (PythonTuple, list)):
-                args += val
-            else:
-                args.insert(0, val)
-            ch = ch.value
-        args = tuple(args)
-        var = self._visit(ch)
+        args = self._visit(stmt.slice)
+        if not isinstance(args, (PythonTuple, list)):
+            args = (args,)
+        var = self._visit(stmt.value)
         var = IndexedElement(var, *args)
         return var
 


Test patch:
diff --git a/tests/epyccel/modules/arrays.py b/tests/epyccel/modules/arrays.py
index 9152252b89..77e394af5f 100644
--- a/tests/epyccel/modules/arrays.py
+++ b/tests/epyccel/modules/arrays.py
@@ -362,6 +362,11 @@ def array_float_nested_F_array_initialization_3(a : 'float[:,:,:]', e : 'float[:
                         ((1., 2., 3.), (1., 2., 3.)))), dtype='float', order="F")
     x[:,:,:,:] = tmp[:,:,:,:]
 
+def array_float_nested_F_array_initialization_mixed(x : 'float[:,:,:](order=F)', y : 'float[:,:](order=F)', z : 'float[:,:](order=F)', a : 'float[:,:,:,:](order=F)'):
+    from numpy import array
+    tmp      = array((x, (y, z, z), x), dtype='float', order="F")
+    a[:,:,:,:] = tmp[:,:,:,:]
+
 ##==============================================================================
 ## TEST ARRAY VIEW STEPS ARRAY INITIALIZATION ORDER C 1D
 ##==============================================================================
@@ -2098,3 +2103,10 @@ def unpack_array_2D_of_known_size():
     arr = array([[1,2,3], [4,5,6], [7,8,9]], dtype='float64')
     x, y, z = arr[:]
     return x.sum(), y.sum(), z.sum()
+
+#==============================================================================
+# Indexing
+#==============================================================================
+
+def multi_layer_index(x : 'int[:]', start : int, stop : int, step : int, idx : int):
+    return x[start:stop:step][idx]
diff --git a/tests/epyccel/modules/tuples.py b/tests/epyccel/modules/tuples.py
index 32c79c9cc9..11d90aea03 100644
--- a/tests/epyccel/modules/tuples.py
+++ b/tests/epyccel/modules/tuples.py
@@ -7,6 +7,7 @@
         'homogenous_tuple_float',
         'homogenous_tuple_string',
         'homogenous_tuple_math',
+        'homogeneous_tuple_of_arrays',
         'inhomogenous_tuple_1',
         'inhomogenous_tuple_2',
         'inhomogenous_tuple_3',
@@ -452,3 +453,17 @@ def test_tuple_inhomogeneous():
 def tuple_different_ranks():
     a = (1,(2,3))
     return a[0], a[1][0], a[1][1]
+
+def homogeneous_tuple_of_arrays():
+    from numpy import array, empty
+    x = array(((1,2), (3,4)), order='F')
+    y = array(((5,4), (7,8)), order='F')
+    z = array(((9,10), (11,12)), order='F')
+    a = (x, y, z)
+    b = empty((3,2,2))
+    for j in range(2):
+        for k in range(2):
+            b[0,j,k] = a[0][j,k]
+            b[1,j,k] = a[1][j,k]
+            b[2,j,k] = a[2][j,k]
+    return b
diff --git a/tests/epyccel/recognised_functions/test_numpy_funcs.py b/tests/epyccel/recognised_functions/test_numpy_funcs.py
index e4757c250a..f214e61141 100644
--- a/tests/epyccel/recognised_functions/test_numpy_funcs.py
+++ b/tests/epyccel/recognised_functions/test_numpy_funcs.py
@@ -4330,6 +4330,22 @@ def test_int(min_int, max_int, dtype):
     assert epyccel_func(fl32) == get_mod(fl32)
     assert epyccel_func(fl64) == get_mod(fl64)
 
+@pytest.mark.xfail(os.environ.get('PYCCEL_DEFAULT_COMPILER', None) == 'intel', reason='Rounding errors. See #1669')
+def test_numpy_mod_mixed_order(language):
+
+    def get_mod(arr1 : 'float[:,:]', arr2 : 'float[:,:](order=F)'):
+        from numpy import mod, shape
+        a = mod(arr1, arr2)
+        s = shape(a)
+        return len(s), s[0], s[1], a[0,1], a[1,0]
+
+    epyccel_func = epyccel(get_mod, language=language)
+
+    fl1 = uniform(min_float / 2, max_float / 2, size = (2,5))
+    fl2 = uniform(min_float / 2, max_float / 2, size = (5,2)).T
+
+    assert epyccel_func(fl1, fl2) == get_mod(fl1, fl2)
+
 @pytest.mark.parametrize( 'language', (
         pytest.param("fortran", marks = [pytest.mark.fortran]),
         pytest.param("c", marks = [
diff --git a/tests/epyccel/test_arrays.py b/tests/epyccel/test_arrays.py
index af63560927..37bab0f153 100644
--- a/tests/epyccel/test_arrays.py
+++ b/tests/epyccel/test_arrays.py
@@ -4105,6 +4105,23 @@ def test_array_float_nested_F_array_initialization_3(language):
 
     assert np.array_equal(x1, x2)
 
+def test_array_float_nested_F_array_initialization_mixed(language):
+    f1 = arrays.array_float_nested_F_array_initialization_mixed
+    f2 = epyccel(f1, language = language)
+
+    x  = np.array(np.random.random((3,2,4)), order="F")
+    y  = np.array(np.random.random((2,4)), order="F")
+    z  = np.array(np.random.random((2,4)), order="F")
+    a  = np.array([x, [y, z, z], x], order="F")
+
+    x1 = np.zeros_like(a)
+    x2 = np.zeros_like(a)
+
+    f1(x, y, z, x1)
+    f2(x, y, z, x2)
+
+    assert np.array_equal(x1, x2)
+
 ##==============================================================================
 ## TEST SIMPLE ARRAY SLICING WITH ORDER C 1D
 ##==============================================================================
@@ -4202,7 +4219,10 @@ def test_array_view_steps_F_2D_3(language):
 #==============================================================================
 
 @pytest.mark.parametrize( 'language', (
-        pytest.param("fortran", marks = pytest.mark.fortran),
+        pytest.param("fortran", marks = [
+            pytest.mark.skip(reason=("Template makes interface ambiguous")),
+            pytest.mark.fortran]
+        ),
         pytest.param("c", marks = pytest.mark.c),
         pytest.param("python", marks = [
             pytest.mark.skip(reason=("Template results in wrong ordered arrays")),
@@ -4228,7 +4248,10 @@ def test_array_ndmin_1(language):
 
 
 @pytest.mark.parametrize( 'language', (
-        pytest.param("fortran", marks = pytest.mark.fortran),
+        pytest.param("fortran", marks = [
+            pytest.mark.skip(reason=("Template makes interface ambiguous")),
+            pytest.mark.fortran]
+        ),
         pytest.param("c", marks = pytest.mark.c),
         pytest.param("python", marks = [
             pytest.mark.skip(reason=("Template results in wrong ordered arrays")),
@@ -4254,7 +4277,10 @@ def test_array_ndmin_2(language):
 
 
 @pytest.mark.parametrize( 'language', (
-        pytest.param("fortran", marks = pytest.mark.fortran),
+        pytest.param("fortran", marks = [
+            pytest.mark.skip(reason=("Template makes interface ambiguous")),
+            pytest.mark.fortran]
+        ),
         pytest.param("c", marks = pytest.mark.c),
         pytest.param("python", marks = [
             pytest.mark.skip(reason=("Template results in wrong ordered arrays")),
@@ -4280,7 +4306,10 @@ def test_array_ndmin_4(language):
 
 
 @pytest.mark.parametrize( 'language', (
-        pytest.param("fortran", marks = pytest.mark.fortran),
+        pytest.param("fortran", marks = [
+            pytest.mark.skip(reason=("Template makes interface ambiguous")),
+            pytest.mark.fortran]
+        ),
         pytest.param("c", marks = pytest.mark.c),
         pytest.param("python", marks = [
             pytest.mark.skip(reason=("Template results in wrong ordered arrays")),
@@ -6363,3 +6392,12 @@ def test_unpacking_2D_of_known_size(language):
     f1 = arrays.unpack_array_2D_of_known_size
     f2 = epyccel(f1, language = language)
     assert f1() == f2()
+
+##==============================================================================
+## TEST INDEXING
+##==============================================================================
+
+def test_multi_layer_index(language):
+    f1 = arrays.multi_layer_index
+    f2 = epyccel(f1, language = language)
+    assert f1(arrays.a_1d, 3, 18, 5, 2) == f2(arrays.a_1d, 3, 18, 5, 2)
diff --git a/tests/epyccel/test_tuples.py b/tests/epyccel/test_tuples.py
index 8c1103f918..1ecec77c43 100644
--- a/tests/epyccel/test_tuples.py
+++ b/tests/epyccel/test_tuples.py
@@ -42,7 +42,10 @@ def compare_python_pyccel( p_output, f_output ):
 
     for pth, pycc in zip(p_output, f_output):
 
-        if isinstance(pth, bool):
+        if isinstance(pth, np.ndarray):
+            assert np.allclose(pth,pycc)
+
+        elif isinstance(pth, bool):
             pycc_bool = (pycc == 1)
             assert(pth == pycc_bool)