diff --git a/Mesh.cpp b/Mesh.cpp
index 7fd709d..a0b4a24 100644
--- a/Mesh.cpp
+++ b/Mesh.cpp
@@ -6,11 +6,22 @@ Mesh::Mesh()
 }
 
 size_t Mesh::add_vertice(const Point& pt) {
-    vertices.push_back(pt);
-    return vertices.size() - 1;
+    return inner_add_vertice(pt, NormalVect());
+}
+
+size_t Mesh::add_vertice(const Point& pt, const Point& normal) {
+    return inner_add_vertice(pt, NormalVect(normal, true));
 }
 
 void Mesh::add_face(const Face& face) {
+    for(int f_id=0; f_id < 3; ++f_id) {
+        if((unsigned)face[f_id] >= normals.size())
+            throw std::out_of_range("Vertice out of range for face");
+        if(!normals[face[f_id]].manual)
+            normals[face[f_id]].dirty = true;
+
+        faces_with_vert[face[f_id]].push_back(faces.size());
+    }
     faces.push_back(face);
 }
 void Mesh::add_face(size_t f1, size_t f2, size_t f3) {
@@ -50,3 +61,37 @@ const std::vector<Point>& Mesh::get_vertices() const {
 const std::vector<Face>& Mesh::get_faces() const {
     return faces;
 }
+
+const Point& Mesh::get_normal(size_t vert_id) {
+    if(vert_id >= normals.size())
+        throw std::out_of_range("Normal out of range");
+
+    NormalVect& norm = normals[vert_id];
+    if(norm.dirty) {
+        compute_normal(vert_id);
+        norm.dirty = false;
+    }
+    return norm.vect;
+}
+
+size_t Mesh::inner_add_vertice(const Point& pt, const Mesh::NormalVect& normal)
+{
+    vertices.push_back(pt);
+    normals.push_back(normal);
+    faces_with_vert.push_back(std::vector<size_t>());
+    return vertices.size() - 1;
+}
+
+void Mesh::compute_normal(size_t vert) {
+    Point normal(0., 0., 0.);
+    for(size_t f_id: faces_with_vert[vert]) {
+        const Face& face = faces[f_id];
+        Point e1 = vertices[face[1]] - vertices[face[0]],
+              e2 = vertices[face[2]] - vertices[face[0]];
+        normal += e1.cross_prod(e2);
+    }
+
+    normal = (1. / ((double) faces_with_vert[vert].size())) * normal;
+    normals[vert].dirty = false;
+    normals[vert].vect = normal;
+}
diff --git a/Mesh.hpp b/Mesh.hpp
index 43d668f..fe06831 100644
--- a/Mesh.hpp
+++ b/Mesh.hpp
@@ -15,6 +15,7 @@ class Mesh {
 
         /// Adds a fresh vertice at `pt`, and returns its ID for further use
         size_t add_vertice(const Point& pt);
+        size_t add_vertice(const Point& pt, const Point& normal);
 
         /// Creates a new face out of the three given point IDs
         void add_face(const Face& face);
@@ -31,15 +32,35 @@ class Mesh {
          */
         void normalize_center(bool keep_position=false);
 
+        /// Get the normal vector for vertice `vert_id`
+        const Point& get_normal(size_t vert_id);
+
         /// Gets the various vertices
         const std::vector<Point>& get_vertices() const;
 
         /// Gets the various faces
         const std::vector<Face>& get_faces() const;
 
+    private: //struct
+        struct NormalVect {
+            NormalVect(const Point& vect, bool manual=false)
+                : vect(vect), dirty(false), manual(manual) {}
+            NormalVect()
+                : vect(Point(0., 0., 0.)), dirty(true), manual(false) {}
+            Point vect;
+            bool dirty;
+            bool manual;
+        };
+
+    private: //meth
+        size_t inner_add_vertice(const Point& pt, const NormalVect& normal);
+        void compute_normal(size_t vert);
+
     private:
         std::vector<Point> vertices;
+        std::vector<NormalVect> normals;
         std::vector<Face> faces;
+        std::vector<std::vector<size_t> > faces_with_vert;
 
         Point center;
 };
diff --git a/common_structures.hpp b/common_structures.hpp
index 347470c..df3453a 100644
--- a/common_structures.hpp
+++ b/common_structures.hpp
@@ -46,6 +46,10 @@ struct Point {
         return (*this += Point(-pt.x, -pt.y, -pt.z));
     }
 
+    Point operator-(const Point& pt) const {
+        return (*this) + (-1 * pt);
+    }
+
     friend Point operator*(double scalar, const Point& pt) {
         return Point(
                 scalar * pt.x,
@@ -53,6 +57,13 @@ struct Point {
                 scalar * pt.z);
     }
 
+    Point cross_prod(const Point& pt) const {
+        return Point(
+                (*this)[1] * pt[2] - (*this)[2] * pt[1],
+                (*this)[2] * pt[0] - (*this)[0] * pt[2],
+                (*this)[0] * pt[1] - (*this)[1] * pt[0]);
+    }
+
     bool operator<(const Point& pt) const {
         /// Lexicographic order on (x, y, z)
         if(x == pt.x) {