diff --git a/Cargo.toml b/Cargo.toml
index 19e0626..c6eb2db 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -8,4 +8,6 @@ edition = "2021"
 [dependencies]
 num = "0.4"
 bempp-traits = { git = "https://github.com/bempp/bempp-rs.git" }
+bempp-tools = { git = "https://github.com/bempp/bempp-rs.git" }
 bempp-element = { git = "https://github.com/bempp/bempp-rs.git" }
+rlst-dense = { git = "https://github.com/linalg-rs/rlst.git" }
diff --git a/src/grid.rs b/src/grid.rs
index db133b1..36dde20 100644
--- a/src/grid.rs
+++ b/src/grid.rs
@@ -1,5 +1,6 @@
 //! Grid implementation
 
+pub mod topology;
 pub mod traits;
 
 pub use crate::grid::traits::Geometry;
diff --git a/src/grid/topology.rs b/src/grid/topology.rs
new file mode 100644
index 0000000..5caec61
--- /dev/null
+++ b/src/grid/topology.rs
@@ -0,0 +1,339 @@
+//! Implementation of grid topology
+
+use crate::grid::traits::{Ownership, Topology};
+use crate::reference_cell::{entity_counts, entity_types};
+use crate::traits::cell::ReferenceCellType;
+use bempp_element::cell;
+use bempp_tools::arrays::AdjacencyList;
+use bempp_traits::arrays::AdjacencyListAccess;
+use bempp_traits::cell::ReferenceCell;
+
+/// Topology of a serial grid
+pub struct SerialTopology {
+    dim: usize,
+    connectivity: Vec<Vec<AdjacencyList<usize>>>,
+    index_map: Vec<usize>,
+    starts: Vec<usize>,
+    cell_types: Vec<ReferenceCellType>,
+    entity_types: Vec<Vec<ReferenceCellType>>,
+}
+
+fn get_reference_cell(cell_type: ReferenceCellType) -> Box<dyn ReferenceCell> {
+    match cell_type {
+        ReferenceCellType::Interval => Box::new(cell::Interval),
+        ReferenceCellType::Triangle => Box::new(cell::Triangle),
+        ReferenceCellType::Quadrilateral => Box::new(cell::Quadrilateral),
+        _ => {
+            panic!("Unsupported cell type (for now)");
+        }
+    }
+}
+
+unsafe impl Sync for SerialTopology {}
+
+impl SerialTopology {
+    pub fn new(cells: &AdjacencyList<usize>, cell_types: &[ReferenceCellType]) -> Self {
+        let mut index_map = vec![];
+        let mut vertices = vec![];
+        let mut starts = vec![];
+        let mut cell_types_new = vec![];
+        let dim = get_reference_cell(cell_types[0]).dim();
+
+        let mut connectivity = vec![];
+        for i in 0..dim + 1 {
+            connectivity.push(vec![]);
+            for _j in 0..dim + 1 {
+                connectivity[i].push(AdjacencyList::<usize>::new());
+            }
+        }
+
+        // dim0 = dim, dim1 = 0
+        for c in cell_types {
+            if dim != get_reference_cell(*c).dim() {
+                panic!("Grids with cells of mixed topological dimension not supported.");
+            }
+            if !cell_types_new.contains(c) {
+                starts.push(connectivity[dim][0].num_rows());
+                cell_types_new.push(*c);
+                let n = get_reference_cell(*c).vertex_count();
+                for (i, cell) in cells.iter_rows().enumerate() {
+                    if cell_types[i] == *c {
+                        index_map.push(i);
+                        // Note: this hard codes that the first n points are at the vertices
+                        let mut row = vec![];
+                        for v in &cell[..n] {
+                            if !vertices.contains(v) {
+                                vertices.push(*v);
+                            }
+                            row.push(vertices.iter().position(|&r| r == *v).unwrap());
+                        }
+                        connectivity[dim][0].add_row(&row);
+                    }
+                }
+            }
+        }
+
+        // dim1 == 0
+        for dim0 in 1..dim {
+            let mut cty = AdjacencyList::<usize>::new();
+            let cells = &connectivity[dim][0];
+            for (i, cell_type) in cell_types_new.iter().enumerate() {
+                let ref_cell = get_reference_cell(*cell_type);
+                let ref_entities = (0..entity_counts(*cell_type)[dim0])
+                    .map(|x| ref_cell.connectivity(dim0, x, 0).unwrap())
+                    .collect::<Vec<Vec<usize>>>();
+
+                let cstart = starts[i];
+                let cend = if i == starts.len() - 1 {
+                    connectivity[2][0].num_rows()
+                } else {
+                    starts[i + 1]
+                };
+                for c in cstart..cend {
+                    let cell = cells.row(c).unwrap();
+                    for e in &ref_entities {
+                        let vertices = e.iter().map(|x| cell[*x]).collect::<Vec<usize>>();
+                        let mut found = false;
+                        for entity in cty.iter_rows() {
+                            if all_equal(entity, &vertices) {
+                                found = true;
+                                break;
+                            }
+                        }
+                        if !found {
+                            cty.add_row(&vertices);
+                        }
+                    }
+                }
+            }
+            connectivity[dim0][0] = cty;
+        }
+
+        // dim0 == dim1 == 0
+        let mut nvertices = 0;
+        let mut cty = AdjacencyList::<usize>::new();
+        let cells = &connectivity[dim][0];
+        for cell in cells.iter_rows() {
+            for j in cell {
+                if *j >= nvertices {
+                    nvertices = *j + 1;
+                }
+            }
+        }
+        for i in 0..nvertices {
+            cty.add_row(&[i]);
+        }
+        connectivity[0][0] = cty;
+
+        // dim0 == dim1
+        for (dim0, c) in connectivity.iter_mut().enumerate().skip(1) {
+            for i in 0..c[0].num_rows() {
+                c[dim0].add_row(&[i]);
+            }
+        }
+
+        // dim0 == dim
+        for dim1 in 1..dim + 1 {
+            let mut cty = AdjacencyList::<usize>::new();
+            let entities0 = &connectivity[dim][0];
+            let entities1 = &connectivity[dim1][0];
+
+            let mut sub_cell_types = vec![ReferenceCellType::Point; entities0.num_rows()];
+            for (i, cell_type) in cell_types_new.iter().enumerate() {
+                let etypes = &entity_types(*cell_type)[dim];
+
+                let cstart = starts[i];
+                let cend = if i == starts.len() - 1 {
+                    connectivity[2][0].num_rows()
+                } else {
+                    starts[i + 1]
+                };
+                for t in sub_cell_types.iter_mut().skip(cstart).take(cend) {
+                    *t = etypes[0];
+                }
+            }
+            for (ei, entity0) in entities0.iter_rows().enumerate() {
+                let entity = get_reference_cell(sub_cell_types[ei]);
+                let mut row = vec![];
+                for i in 0..entity.entity_count(dim1) {
+                    let vertices = entity
+                        .connectivity(dim1, i, 0)
+                        .unwrap()
+                        .iter()
+                        .map(|x| entity0[*x])
+                        .collect::<Vec<usize>>();
+                    for (j, entity1) in entities1.iter_rows().enumerate() {
+                        if all_equal(&vertices, entity1) {
+                            row.push(j);
+                            break;
+                        }
+                    }
+                }
+                cty.add_row(&row);
+            }
+            connectivity[dim][dim1] = cty
+        }
+
+        // dim1 < dim0
+        for dim1 in 1..dim + 1 {
+            for dim0 in dim1 + 1..dim {
+                let mut cty = AdjacencyList::<usize>::new();
+                let entities0 = &connectivity[dim0][0];
+                let entities1 = &connectivity[dim1][0];
+                let cell_to_entities0 = &connectivity[dim][dim0];
+
+                let mut sub_cell_types = vec![ReferenceCellType::Point; entities0.num_rows()];
+                for (i, cell_type) in cell_types_new.iter().enumerate() {
+                    let etypes = &entity_types(*cell_type)[dim0];
+
+                    let cstart = starts[i];
+                    let cend = if i == starts.len() - 1 {
+                        connectivity[2][0].num_rows()
+                    } else {
+                        starts[i + 1]
+                    };
+                    for c in cstart..cend {
+                        for (e, t) in cell_to_entities0.row(c).unwrap().iter().zip(etypes) {
+                            sub_cell_types[*e] = *t;
+                        }
+                    }
+                }
+                for (ei, entity0) in entities0.iter_rows().enumerate() {
+                    let entity = get_reference_cell(sub_cell_types[ei]);
+                    let mut row = vec![];
+                    for i in 0..entity.entity_count(dim1) {
+                        let vertices = entity
+                            .connectivity(dim1, i, 0)
+                            .unwrap()
+                            .iter()
+                            .map(|x| entity0[*x])
+                            .collect::<Vec<usize>>();
+                        for (j, entity1) in entities1.iter_rows().enumerate() {
+                            if all_equal(&vertices, entity1) {
+                                row.push(j);
+                                break;
+                            }
+                        }
+                    }
+                    cty.add_row(&row);
+                }
+                connectivity[dim0][dim1] = cty;
+            }
+        }
+
+        // dim1 > dim0
+        for dim1 in 1..dim + 1 {
+            for dim0 in 0..dim1 {
+                let mut data = vec![vec![]; connectivity[dim0][0].num_rows()];
+                for (i, row) in connectivity[dim1][dim0].iter_rows().enumerate() {
+                    for v in row {
+                        data[*v].push(i);
+                    }
+                }
+                for row in data {
+                    connectivity[dim0][dim1].add_row(&row);
+                }
+            }
+        }
+
+        let mut all_entity_types = vec![vec![], vec![], vec![], vec![]];
+        for c in &cell_types_new {
+            let et = entity_types(*c);
+            for (dim, t) in et.iter().enumerate() {
+                for e in t {
+                    if !all_entity_types[dim].contains(e) {
+                        all_entity_types[dim].push(*e);
+                    }
+                }
+            }
+        }
+
+        Self {
+            dim,
+            connectivity,
+            index_map,
+            starts,
+            cell_types: cell_types_new,
+            entity_types: all_entity_types,
+        }
+    }
+}
+
+fn all_equal(a: &[usize], b: &[usize]) -> bool {
+    if a.len() != b.len() {
+        false
+    } else {
+        all_in(a, b)
+    }
+}
+
+fn all_in(a: &[usize], b: &[usize]) -> bool {
+    for i in a {
+        if !b.contains(i) {
+            return false;
+        }
+    }
+    true
+}
+
+impl Topology for SerialTopology {
+    //type Connectivity = AdjacencyList<usize>;
+
+    fn dim(&self) -> usize {
+        self.dim
+    }
+    fn index_map(&self) -> &[usize] {
+        &self.index_map
+    }
+    fn entity_count(&self, dim: usize) -> usize {
+        self.connectivity[dim][0].num_rows()
+    }
+    fn cell(&self, index: usize) -> Option<&[usize]> {
+        if index < self.entity_count(self.dim) {
+            Some(self.connectivity[self.dim][0].row(index).unwrap())
+        } else {
+            None
+        }
+    }
+    fn cell_type(&self, index: usize) -> Option<ReferenceCellType> {
+        for (i, start) in self.starts.iter().enumerate() {
+            let end = if i == self.starts.len() - 1 {
+                self.connectivity[2][0].num_rows()
+            } else {
+                self.starts[i + 1]
+            };
+            if *start <= index && index < end {
+                return Some(self.cell_types[i]);
+            }
+        }
+        None
+    }
+
+    fn entity_types(&self, dim: usize) -> &[ReferenceCellType] {
+        &self.entity_types[dim]
+    }
+
+    /// Get the indices of entities of type `etype` that are connected to each cell of type `cell_type`
+    fn cell_entities(
+        &self,
+        _cell_type: ReferenceCellType,
+        _etype: ReferenceCellType,
+    ) -> Option<&[usize]> {
+        // TODO
+        None
+    }
+
+    /// Get the indices of entities of dimension `dim` that are connected to the entity of type `etype` with index `index`
+    fn connectivity(
+        &self,
+        _etype: ReferenceCellType,
+        _index: usize,
+        _dim: usize,
+    ) -> Option<&[usize]> {
+        None
+    }
+
+    fn entity_ownership(&self, _dim: usize, _index: usize) -> Ownership {
+        Ownership::Owned
+    }
+}
diff --git a/src/reference_cell.rs b/src/reference_cell.rs
index d735e98..95320a3 100644
--- a/src/reference_cell.rs
+++ b/src/reference_cell.rs
@@ -55,3 +55,55 @@ pub fn vertices<T: Float>(cell: ReferenceCellType) -> Vec<T> {
         ],
     }
 }
+
+/// The edges of the reference cell
+pub fn edges(cell: ReferenceCellType) -> Vec<Vec<usize>> {
+    match cell {
+        ReferenceCellType::Point => vec![],
+        ReferenceCellType::Interval => vec![vec![0, 1]],
+        ReferenceCellType::Triangle => vec![vec![1, 2], vec![0, 2], vec![0, 1]],
+        ReferenceCellType::Quadrilateral => vec![vec![0, 1], vec![0, 2], vec![1, 3], vec![2, 3]],
+        _ => {
+            panic!("Not implemented yet");
+        }
+    }
+}
+
+pub fn entity_types(cell: ReferenceCellType) -> Vec<Vec<ReferenceCellType>> {
+    match cell {
+        ReferenceCellType::Point => vec![vec![ReferenceCellType::Point], vec![], vec![], vec![]],
+        ReferenceCellType::Interval => vec![
+            vec![ReferenceCellType::Point; 2],
+            vec![ReferenceCellType::Interval],
+            vec![],
+            vec![],
+        ],
+        ReferenceCellType::Triangle => vec![
+            vec![ReferenceCellType::Point; 3],
+            vec![ReferenceCellType::Interval; 3],
+            vec![ReferenceCellType::Triangle],
+            vec![],
+        ],
+        ReferenceCellType::Quadrilateral => vec![
+            vec![ReferenceCellType::Point; 4],
+            vec![ReferenceCellType::Interval; 4],
+            vec![ReferenceCellType::Quadrilateral],
+            vec![],
+        ],
+        _ => {
+            panic!("Not implemented yet");
+        }
+    }
+}
+
+pub fn entity_counts(cell: ReferenceCellType) -> Vec<usize> {
+    match cell {
+        ReferenceCellType::Point => vec![1, 0, 0, 0],
+        ReferenceCellType::Interval => vec![2, 1, 0, 0],
+        ReferenceCellType::Triangle => vec![3, 3, 1, 0],
+        ReferenceCellType::Quadrilateral => vec![4, 4, 1, 0],
+        _ => {
+            panic!("Not implemented yet");
+        }
+    }
+}