diff --git a/aplK.py b/aplK.py
index be84089..0451b68 100644
--- a/aplK.py
+++ b/aplK.py
@@ -1,71 +1,105 @@
-reset()
+r"""
+<Short one-line summary that ends with no period>
+
+<Paragraph description>
+
+EXAMPLES::
+
+<Lots and lots of examples>
+
+AUTHORS:
+
+- YOUR NAME (2005-01-03): initial version
+
+- person (date in ISO year-month-day format): short desc
+
+"""
+
+# ****************************************************************************
+#       Copyright (C) 2021 YOUR NAME <your email>
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 2 of the License, or
+# (at your option) any later version.
+#                  https://www.gnu.org/licenses/
+# ****************************************************************************
 from sage.structure.parent import Parent
 from sage.structure.element import Element
 from sage.structure.unique_representation import UniqueRepresentation
 from sage.categories.rings import Rings
 from sage.structure.richcmp import rich_to_bool, op_EQ, op_NE, op_LT, op_LE, op_GT, op_GE
+from sage.categories.morphism import Morphism
+from sage.rings.rational_field import QQ
+from sage.algebras.commutative_dga import GradedCommutativeAlgebra
+from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
+from sage.functions.other import factorial
 
 #####################################################################
 class Apln(UniqueRepresentation):
-    
-  def __init__(self,n):
-      self.n = n
-      ts = ['t'+str(t) for t in range(n+1)]
-      ys = ['y'+str(t) for t in range(n+1)]
-      degrees = tuple((n+1)*[(0,2)])+tuple((n+1)*[(1,2)])
-      self.gca = GradedCommutativeAlgebra(QQ, names=ts+ys, degrees=degrees)
-      t_relation = -1
-      y_relation = 0
-      #t_relation = sum(self.gca.gens()[:n])-1
-      #y_relation = sum(self.gca.gens()[n:])
-      for i in range(n+1):
-          t_relation += self.gca.gens()[i]
-          y_relation += self.gca.gens()[n+1+i]
-      self.gca_with_relations = self.gca.quotient(self.gca.ideal(t_relation,y_relation),check=False)
-      dfs = {self.gca.gens_dict()[ts[i]]:self.gca.gens_dict()[ys[i]] for i in range(n+1)}
-      self.cdga = self.gca.cdg_algebra(dfs)
-  
-  def _repr_(self):
-    return "Algebra Apln of dimension %s" % self.n
-    
-  def inject_variables(self):
-      self.cdga.inject_variables()
-      
-  def generators(self):
-      return self.cdga.gens()
-      
-  def reduce_with_relations(self,pol):
-      return self.gca_with_relations(pol)
-      
-  def homogeneousPart(self,pol,degree):
-      parts = [pol.homogeneous_parts()[k] for k in pol.homogeneous_parts() if list(k)[0]==degree]
-      hompart = 0
-      for p in parts:
-          hompart += p
-      return hompart
 
-  def face(self,i):
-      if(self.n==0):
-          return 0
-      domain = self.cdga
-      codomain = self.cdga
-      tImages = []
-      yImages = []
-      n = len(domain.gens())//2
-      for k in range(n):
-          if(k<i):
-              tImages.append(codomain.gens()[k])
-              yImages.append(codomain.gens()[n+k])
-          elif(k==i):
-              tImages.append(0)
-              yImages.append(0)
-          elif(k>i):
-              tImages.append(codomain.gens()[k-1])
-              yImages.append(codomain.gens()[n+k-1])
-      images = tImages+yImages
-      return domain.hom(codomain=codomain,im_gens=images)
+    def __init__(self,n):
+        self.n = n
+        ts = ['t'+str(t) for t in range(n+1)]
+        ys = ['y'+str(t) for t in range(n+1)]
+        degrees = tuple((n+1)*[(0,2)])+tuple((n+1)*[(1,2)])
+        self.gca = GradedCommutativeAlgebra(QQ, names=ts+ys, degrees=degrees)
+        t_relation = -1
+        y_relation = 0
+        #t_relation = sum(self.gca.gens()[:n])-1
+        #y_relation = sum(self.gca.gens()[n:])
+        for i in range(n+1):
+            t_relation += self.gca.gens()[i]
+            y_relation += self.gca.gens()[n+1+i]
+        self.gca_with_relations = self.gca.quotient(self.gca.ideal(t_relation,y_relation),check=False)
+        dfs = {self.gca.gens_dict()[ts[i]]:self.gca.gens_dict()[ys[i]] for i in range(n+1)}
+        self.cdga = self.gca.cdg_algebra(dfs)
+        #self.inject_variables()
+  
+    def _repr_(self):
+        return "Algebra Apln of dimension %s" % self.n
+    
+    def inject_variables(self):
+        self.cdga.inject_variables()
       
-  def _reduce_polynomial(self,p,n):
+    def generators(self):
+        return self.cdga.gens()
+      
+    def reduce_with_relations(self,pol):
+        return self.gca_with_relations(pol)
+      
+    def homogeneousPart(self,pol,degree):
+
+        parts = [pol.homogeneous_parts()[k] for k in pol.homogeneous_parts() if list(k)[0]==degree]
+        hompart = 0
+        for p in parts:
+            hompart += p
+        return hompart
+
+    def face(self,i):
+
+        if(self.n==0):
+            return 0
+        domain = self.cdga
+        codomain = self.cdga
+        tImages = []
+        yImages = []
+        n = len(domain.gens())//2
+        for k in range(n):
+            if(k<i):
+                tImages.append(codomain.gens()[k])
+                yImages.append(codomain.gens()[n+k])
+            elif(k==i):
+                tImages.append(0)
+                yImages.append(0)
+            elif(k>i):
+                tImages.append(codomain.gens()[k-1])
+                yImages.append(codomain.gens()[n+k-1])
+        images = tImages+yImages
+        return domain.hom(codomain=codomain,im_gens=images)
+      
+    def _reduce_polynomial(self,p,n):
+
     ## Input:
     #  - p: polynomial in variables t_0,..t_n,y_0,...,y_n
     #  - n: dimension of the algebra (apl)_n where p is defined
@@ -73,25 +107,25 @@ class Apln(UniqueRepresentation):
     #  - The polynomial after applying the corresponding relations to reduce the variables t_n, y_n
     #  - For practical reasons, we change the nomenclature in this way: t_i ==> xs_i, y_i ==> zs_i
     
-    xs = ['x'+str(x) for x in range(n+1)]
-    zs = ['z'+str(z) for z in range(n+1)]
-    R = PolynomialRing(QQ,names=xs+zs)
-    r_gens = R.gens()
-    gca_gens = self.gca.gens()
-    pol = R(0)
-    for monomial_exponents in list(p.dict().keys()):
-      monomial = R(1)
-      coef = p.dict()[monomial_exponents]
-      for i in range(n+1):
-        monomial *= r_gens[i]**monomial_exponents[i]
-        monomial *= r_gens[n+1+i]**monomial_exponents[self.n+1+i]
-      pol += coef*monomial
-    xn_relation = 1-sum(r_gens[:n])
-    zn_relation = -sum(r_gens[n+1:-1])
-    pol = pol.subs({r_gens[n]:xn_relation,r_gens[2*n+1]:zn_relation})
-    return pol
+        xs = ['x'+str(x) for x in range(n+1)]
+        zs = ['z'+str(z) for z in range(n+1)]
+        R = PolynomialRing(QQ,names=xs+zs)
+        r_gens = R.gens()
+        gca_gens = self.gca.gens()
+        pol = R(0)
+        for monomial_exponents in list(p.dict().keys()):
+            monomial = R(1)
+            coef = p.dict()[monomial_exponents]
+            for i in range(n+1):
+                monomial *= r_gens[i]**monomial_exponents[i]
+                monomial *= r_gens[n+1+i]**monomial_exponents[self.n+1+i]
+            pol += coef*monomial
+        xn_relation = 1-sum(r_gens[:n])
+        zn_relation = -sum(r_gens[n+1:-1])
+        pol = pol.subs({r_gens[n]:xn_relation,r_gens[2*n+1]:zn_relation})
+        return pol
     
-  def _convert_integrable_polynomial(self,p,n):
+    def _convert_integrable_polynomial(self,p,n):
     ## Input:
     #  - p: polynomial in variables t_0,..t_n,y_0,...,y_n
     #  - n: dimension of the algebra (apl)_n where p is defined
@@ -100,284 +134,260 @@ class Apln(UniqueRepresentation):
     #    (After reducing the variables t_{n},y_{n} by the corresponding relations, 
     #     the algorithm removes the variables y_{0},..,y_{n-1}). 
 
-    xs = ['x'+str(x) for x in range(n)]
-    zs = ['z'+str(z) for z in range(n)]
-    R = PolynomialRing(QQ,names=xs+zs)
-    pol = R(self._reduce_polynomial(p,n))
-    r_gens = R.gens()
-    relations = {}
-    for i in range(n):
-      relations[r_gens[n+i]]=1
-    return pol.subs(relations)
+
+        xs = ['x'+str(x) for x in range(n)]
+        zs = ['z'+str(z) for z in range(n)]
+        R = PolynomialRing(QQ,names=xs+zs)
+        pol = R(self._reduce_polynomial(p,n))
+        r_gens = R.gens()
+        relations = {}
+        for i in range(n):
+            relations[r_gens[n+i]]=1
+        return pol.subs(relations)
         
-  def integral(self,p,n):
+    def integral(self,p,n):
     ## Input:
     #  - p: polynomial in variables t_0,..t_n,y_0,...,y_n
     #  - n: dimension of the algebra (apl)_n where p is defined
     ## Output:
     #  - the integral of the polynomial after reducing the variables t_n, y_n by the corresponding relations 
     #
-    xs = ['x'+str(x) for x in range(n)]
-    R = PolynomialRing(QQ,names=xs)
-    r_gens = R.gens()
-    pol = R(self._convert_integrable_polynomial(p,n))
-    result = pol
-    for i in range(n):
-      result = result.integral(r_gens[n-1-i])
-      result = result.subs({r_gens[n-1-i]:1-sum(r_gens[:n-1-i])})
-    return R(((-1)**(n*(n-1)/2)))*result
+
+        xs = ['x'+str(x) for x in range(n)]
+        R = PolynomialRing(QQ,names=xs)
+        r_gens = R.gens()
+        pol = R(self._convert_integrable_polynomial(p,n))
+        result = pol
+        for i in range(n):
+            result = result.integral(r_gens[n-1-i])
+            result = result.subs({r_gens[n-1-i]:1-sum(r_gens[:n-1-i])})
+        return R(((-1)**(n*(n-1)/2)))*result
+        
+    def integral2(self,p,n):
+    ## Input:
+    #  - p: polynomial in variables t_0,..t_n,y_0,...,y_n
+    #  - n: dimension of the algebra (apl)_n where p is defined
+    ## Output:
+    #  - the integral of the polynomial after reducing the variables t_n, y_n by the corresponding relations 
+    #
+
+        xs = ['x'+str(x) for x in range(n)]
+        R = PolynomialRing(QQ,names=xs)
+        r_gens = R.gens()
+        pol = R(self._convert_integrable_polynomial(p,n))
+        result = 0
+        for k in pol.dict().keys(): # k is a monomial
+            numerator = 1
+            denominator = n
+            for e in k: # e is an exponents
+                numerator *= factorial(e)
+                denominator += e
+            coef = pol.dict()[k]
+            result += coef*numerator/factorial(denominator)
+        return R(((-1)**(n*(n-1)/2)))*result
             
 #####################################################################
+
+class MorphismFromQQToAplK(Morphism):
+
+    def __init__(self, codomain):
+        self.simplicialComplex= codomain.K
+        Morphism.__init__(self, QQ, codomain)
+
+    def _call_(self, x):
+        S = self.simplicialComplex
+        params1 = {}
+        for d in range(S.dimension()+1):
+            for c in S._n_cells_sorted(d):
+                params1[tuple(c)]=x
+        A = AplK(S)
+        return AplKElement(A,params1)
+    
+#####################################################################
+
 class AplK(UniqueRepresentation, Parent):
     
-  def __init__(self,K, category=None):
-    self.K = K
-    self.n = K.dimension()
-    Parent.__init__(self,category=category )
+    def __init__(self,K, category=None):
+        self.K = K
+        self.n = K.dimension()
+        Parent.__init__(self,category=category )
+        f = MorphismFromQQToAplK(self)
+        self.register_coercion(f)
    
-  def _repr_(self):
-    return "Algebra AplK over %s" % self.K
+    def _repr_(self):
+        return "Algebra AplK over %s" % self.K
     
-  def _element_constructor_(self,values):
-    return self.element_class(self,values)
+    def _element_constructor_(self,values):
+        return self.element_class(self,values)
     
-  def dimension(self):
-    return self.K.dimension()
+    def dimension(self):
+        return self.K.dimension()
     
-  def simplicialComplex(self):
-    return self.K
+    def simplicialComplex(self):
+        return self.K
     
-  def _n_cells_sorted(self,n):
-    return self.K._n_cells_sorted(n)
+    def _n_cells_sorted(self,n):
+        return self.K._n_cells_sorted(n)
     
 #####################################################################
 class AplKElement(Element):
 
-  def __init__(self, parent,values):
-    self.elements = values
-    self.apl = Apln(parent.dimension())
-    Element.__init__(self,parent)
-    self.check()
+    def __init__(self, parent,values):
+        self.elements = values
+        self.apl = Apln(parent.dimension())
+        Element.__init__(self,parent)
+        self.check()
     
-  def _repr_(self):
-    return str(self.elements)
+    def _repr_(self):
+        return str(self.elements)
     
-  def dimension(self):
-    return self.parent().dimension()
+    def dimension(self):
+        return self.parent().dimension()
     
-  def nullElement(self):
-    C = self.__class__
-    values = {}
-    for k in self.elements.keys():
-      values[k]=0
-    return C(self.parent(),values)    
+    def nullElement(self):
+        C = self.__class__
+        values = {}
+        for k in self.elements.keys():
+            values[k]=0
+        return C(self.parent(),values)    
     
-  def _add_(self,other):
-    C = self.__class__
-    if(self.elements.keys()!=other.elements.keys()):
-      raise ValueError("error: keys do not match")
-    newValues = {}
-    for k in self.elements.keys():
-      newValues[k] = self.elements[k]+other.elements[k]
-    return C(self.parent(),newValues)
+    def _add_(self,other):
+        r"""
+        Return sum of the elements
+
+        INPUT:
+
+        - ``other`` -- An AplKEelement
+
+
+        OUTPUT: An AplKElement equal to the sum of self and the element passed as argument
+
+        EXAMPLES:
+
+        This example illustrates ... ::
+
+            sage: S = SimplicialComplex([[0,1],[0,2],[1,2]],is_mutable=False) 
+            sage: apl = Apln(S.dimension())
+            sage: (t0, t1, y0, y1) = apl.generators()  
+            sage: params1 = {(0,): t0, (1,): t0, (2,):3*t0, (0,1): t0+t1,(0,2): t0+3*t1, (1,2): t0+3*t1}
+            sage: params2 = {(0,): 3*t0, (1,): 4*t0, (2,):5*t0, (0,1): 3*t0+4*t1,(0,2): 3*t0+5*t1, (1,2): 4*t0+5*t1}
+            sage: A = AplK(S)
+            sage: e1 = AplKElement(A,params1)
+            sage: e2 = AplKElement(A,params2)  
+            sage: e1+e2
+            {(0,): 4*t0, (1,): 5*t0, (2,): 8*t0, (0, 1): 4*t0 + 5*t1, (0, 2): 4*t0 + 8*t1, (1, 2): 5*t0 + 8*t1}
+          
+        """
+        C = self.__class__
+        if(self.elements.keys()!=other.elements.keys()):
+            raise ValueError("error: keys do not match")
+        newValues = {}
+        for k in self.elements.keys():
+            newValues[k] = self.elements[k]+other.elements[k]
+        return C(self.parent(),newValues)
       
-  def _sub_(self,other):
-    C = self.__class__
-    if(self.elements.keys()!=other.elements.keys()):
-      raise ValueError("error: keys do not match")
-    newValues = {}
-    for k in self.elements.keys():
-      newValues[k] = self.elements[k]-other.elements[k]
-    return C(self.parent(),newValues)
+    def _sub_(self,other):
+
+        C = self.__class__
+        if(self.elements.keys()!=other.elements.keys()):
+            raise ValueError("error: keys do not match")
+        newValues = {}
+        for k in self.elements.keys():
+            newValues[k] = self.elements[k]-other.elements[k]
+        return C(self.parent(),newValues)
       
-  def _mul_(self,other):
-    C = self.__class__
-    if(self.elements.keys()!=other.elements.keys()):
-      raise ValueError("error: keys do not match")
-    newValues = {}
-    for k in self.elements.keys():
-      newValues[k] = self.elements[k]*other.elements[k]
-    return C(self.parent(),newValues)
+    def _mul_(self,other):
+
+        C = self.__class__
+        if(self.elements.keys()!=other.elements.keys()):
+            raise ValueError("error: keys do not match")
+        newValues = {}
+        for k in self.elements.keys():
+            newValues[k] = self.elements[k]*other.elements[k]
+        return C(self.parent(),newValues)
     
-  def _richcmp_(self, other, op):
-    if(op==op_EQ):
-      if(self.elements.keys()!=other.elements.keys()):
-        return False
-      for k in self.elements.keys():
-        if(self.elements[k]!=other.elements[k]):
-          return False
-      return True
-    elif(op==op_NE):
-      return not self._richcmp_(other,op_EQ)
-    else:
-      raise NotImplementedError("operation not implemented for this element class")
+    def _richcmp_(self, other, op):
+
+        if(op==op_EQ):
+            if(self.elements.keys()!=other.elements.keys()):
+                return False
+            for k in self.elements.keys():
+                if(self.elements[k]!=other.elements[k]):
+                    return False
+            return True
+        elif(op==op_NE):
+            return not self._richcmp_(other,op_EQ)
+        else:
+            raise NotImplementedError("operation not implemented for this element class")
       
-  def differential(self):
-    C = self.__class__
-    newValues = {}
-    for k in self.elements.keys():
-      newValues[k] = self.apl.cdga(self.elements[k]).differential()
-    return C(self.parent(),newValues)
+    def differential(self):
+
+        C = self.__class__
+        newValues = {}
+        for k in self.elements.keys():
+            newValues[k] = self.apl.cdga(self.elements[k]).differential()
+        return C(self.parent(),newValues)
       
-  def check(self):
-    S = self.parent().simplicialComplex()
-    dim = S.dimension()
-    for n in range(1,dim+1):
-      print("Checking dimension " + str(n))
-      for simplex in S._n_cells_sorted(n):
-        for f in range(n+1):
-          sFace = tuple(S.face(simplex,f))
-          aFace = self.apl.face(f)(self.elements[tuple(simplex)])
-          if(self.elements[sFace] != aFace):
-            raise ValueError("the assignation is not valid")
+    def check(self):
+
+        S = self.parent().simplicialComplex()
+        dim = S.dimension()
+        for n in range(1,dim+1):
+            #print("Checking dimension " + str(n))
+            for simplex in S._n_cells_sorted(n):
+                for f in range(n+1):
+                    sFace = tuple(S.face(simplex,f))
+                    aFace = self.apl.face(f)(self.elements[tuple(simplex)])
+                    if(self.elements[sFace] != aFace):
+                        raise ValueError("the assignation is not valid")
             
-  def homogeneous_component(self,degree):
-    C = self.__class__
-    component = {}
-    for k in self.elements.keys():
-      component[k] = self.apl.homogeneousPart(self.elements[k],degree)
-    return C(self.parent(),component)
+    def homogeneous_component(self,degree):
+
+        C = self.__class__
+        component = {}
+        for k in self.elements.keys():
+            component[k] = self.apl.homogeneousPart(self.elements[k],degree)
+        return C(self.parent(),component)
     
-  def homogeneous_components(self):
-    components = {}
-    for n in range(self.dimension()+1):
-      components[n] = self.homogeneous_component(n)
-    return components
+    def homogeneous_components(self):
+
+        components = {}
+        for n in range(self.dimension()+1):
+            components[n] = self.homogeneous_component(n)
+        return components
                
-  def is_homogeneous(self):
-    isHomogeneous = True
-    for k in self.elements.keys():
-      isHomogeneous = isHomogeneous and self.apl.cdga(self.elements[k]).is_homogeneous()
-    return isHomogeneous
+    def is_homogeneous(self):
+
+        isHomogeneous = True
+        for k in self.elements.keys():
+            isHomogeneous = isHomogeneous and self.apl.cdga(self.elements[k]).is_homogeneous()
+        return isHomogeneous
   
-  def is_in_kernel(self):
-    diff = self.differential()
-    if(diff == self.nullElement()):
-      return True
-    else:
-      return False
+    def is_in_kernel(self):
+
+        diff = self.differential()
+        if(diff == self.nullElement()):
+            return True
+        else:
+            return False
       
-  def _integrable_component(self):
-    C = self.__class__
-    component = {}
-    for k in self.elements.keys():
-      degree = len(k)-1
-      component[k] = self.apl.homogeneousPart(self.elements[k],degree)
-      #component[k] = self.homogeneous_component(dim).elements[k]
-    return component
+    def _integrable_component(self):
+        C = self.__class__
+        component = {}
+        for k in self.elements.keys():
+            degree = len(k)-1
+            component[k] = self.apl.homogeneousPart(self.elements[k],degree)
+            #component[k] = self.homogeneous_component(dim).elements[k]
+        return component
     
-  def integral(self):
-    component = self._integrable_component()
-    cochain = {}
-    #integration of each polynomial
-    for k in component.keys():
-      print("Integrating simplex " + str(k))
-      n = len(k)
-      pol = self.apl.cdga(component[k])
-      cochain[k] = self.apl.integral(pol,n)
-    return cochain
+    def integral(self):
+        component = self._integrable_component()
+        cochain = {}
+        #integration of each polynomial
+        for k in component.keys():
+            print("Integrating simplex " + str(k))
+            n = len(k)-1
+            pol = self.apl.cdga(component[k])
+            cochain[k] = self.apl.integral(pol,n)
+        return cochain
 
-  
-      
-#####################################################################
-#S = SimplicialComplex([[0,1,2]],is_mutable=False)  
-S = SimplicialComplex([[0,1],[0,2],[1,2]],is_mutable=False)                                                                                 
-apl = Apln(S.dimension())
-apl.inject_variables()
-
-params1 = {}
-params1[(0,)] = t0
-params1[(1,)] = t0
-params1[(2,)] = 3*t0
-params1[(0,1)] = t0+t1
-params1[(0,2)] = t0+3*t1
-params1[(1,2)] = t0+3*t1
-
-params2 = {}
-params2[(0,)] = 3*t0
-params2[(1,)] = 4*t0
-params2[(2,)] = 5*t0
-params2[(0,1)] = 3*t0+4*t1
-params2[(0,2)] = 3*t0+5*t1
-params2[(1,2)] = 4*t0+5*t1
-
-params3 = {}
-for d in range(S.dimension()+1):
-  for c in S._n_cells_sorted(d):
-    params3[tuple(c)]=randint(1,100)*t0
-    
-params4 = {}
-params4[(0,)] = 3*t0+y0
-params4[(1,)] = 4*t0
-params4[(2,)] = 5*t0
-params4[(0,1)] = 3*t0+4*t1+y0
-params4[(0,2)] = 3*t0+5*t1+y0
-params4[(1,2)] = 4*t0+5*t1
-
-A = AplK(S)
-e1 = AplKElement(A,params1)
-e2 = AplKElement(A,params2)
-#e3 = AplKElement(A,params3)
-e4 = AplKElement(A,params4)
-
-
-class morphism_from_QQ_to_AplK(Morphism):
-
-  def __init__(self, codomain):
-    self.simplicialComplex= codomain.K
-    Morphism.__init__(self, QQ, codomain)
-
-  def _call_(self, x):
-    S = self.simplicialComplex
-    params1 = {}
-    for d in range(S.dimension()+1):
-      for c in S._n_cells_sorted(d):
-        params1[tuple(c)]=x
-    A = AplK(S)
-    return AplKElement(A,params1)
-    
-f = morphism_from_QQ_to_AplK(A)
-A.register_coercion(f)
-
-#p = t0*t1**2*y0+t0*y1
-
-
-S = SimplicialComplex([[0,1,2,3,4,5,6]],is_mutable=False)                                                                                 
-apl = Apln(S.dimension())
-apl.inject_variables()
-p = y0*y1*y2*y3*y4*y5+y0*y1*y2*y3*y6+t4*t5
-
-initial_assignation = {}
-initial_assignation[(0,1,2,3,4,5,6)] = p
-
-def generateElement(S,initial_assignation):
-  apl = Apln(S.dimension())
-  aplk_structure = initial_assignation.copy()
-  
-  for n in reversed(range(1,S.dimension()+1)):
-    print("Generating dimension "+ str(n))
-    for simplex in S._n_cells_sorted(n):
-      p = aplk_structure[tuple(simplex)]
-      for f in range(n+1):
-        sFace = tuple(S.face(simplex,f))
-        aFace = apl.face(f)(p)
-        if(sFace not in aplk_structure.keys()):
-          aplk_structure[tuple(sFace)] = aFace
-        elif(aplk_structure[tuple(sFace)] != aFace):
-          raise ValueError("Error during the checking. The assignation is not valid")
-  A = AplK(S)
-  return AplKElement(A,aplk_structure)
-          
-          
-g = generateElement(S,initial_assignation)
-#g.check()
-#cochain = g.integral()
-  
-#n = 1
-#a = Apln(n)
-#a.inject_variables()
-#p = t0*y1
-##print(p)
-#print(a._reduce_polynomial(p,n))
-#print(a._convert_integrable_polynomial(p,n))
-#print(a.integral(p,n))