diff --git a/README.md b/README.md
index a842aa3..dcdd536 100644
--- a/README.md
+++ b/README.md
@@ -7,7 +7,7 @@
 [![Python Versions][pyversions-img]][pyversions-url]
 [![License: MIT][license-img]][license-url]
 
-![compmec-section logo](docs/source/img/logo.svg)
+![compmec-section logo](docs/source/img/logo.png)
 
 A python package for analysis of beam sections by using Boundary Element Method.
 
diff --git a/docs/source/img/logo.png b/docs/source/img/logo.png
new file mode 100644
index 0000000..7747f99
Binary files /dev/null and b/docs/source/img/logo.png differ
diff --git a/pyproject.toml b/pyproject.toml
index adf89ec..4f954e1 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,6 +1,6 @@
 [tool.poetry]
 name = "compmec-section"
-version = "0.3.0"
+version = "0.4.0"
 description = "Analysis of beams cross-section using the boundary element method"
 authors = ["Carlos Adir <carlos.adir.leite@gmail.com>"]
 readme = "README.md"
diff --git a/src/compmec/section/__init__.py b/src/compmec/section/__init__.py
index 8b72ac4..19ed414 100644
--- a/src/compmec/section/__init__.py
+++ b/src/compmec/section/__init__.py
@@ -3,10 +3,13 @@
 It uses mainly curves as boundary to compute the elements.
 """
 
+from .curve import Curve
+from .dataio import JsonIO
+from .geometry import Geometry
 from .material import Isotropic, Material
 from .section import Section
 
-__version__ = "0.3.0"
+__version__ = "0.4.0"
 
 if __name__ == "__main__":
     pass
diff --git a/src/compmec/section/abcs.py b/src/compmec/section/abcs.py
new file mode 100644
index 0000000..0cbb38a
--- /dev/null
+++ b/src/compmec/section/abcs.py
@@ -0,0 +1,344 @@
+"""
+File that contains base classes
+"""
+
+from __future__ import annotations
+
+from abc import ABC, abstractmethod
+from typing import Any, Dict, Optional, Tuple, Union
+
+
+class Tracker(ABC):  # pylint: disable=too-few-public-methods
+    """
+    Parent class to track all the instances
+    """
+
+    instances = {}
+
+    @classmethod
+    def clear(cls, keys: Optional[Tuple[int]] = None):
+        """
+        Remove node labels
+        """
+        if keys is None:
+            cls.instances.clear()
+            return
+        for key in keys:
+            if key in cls.instances:
+                cls.instances.pop(key)
+
+    @classmethod
+    @abstractmethod
+    def _next_available_key(cls) -> Any:
+        """
+        Protected method that gives next key
+        """
+        raise NotImplementedError
+
+    def _get_key(self) -> Any:
+        """
+        Protected method that gives the key of the instance by iterating
+        on the instances dictionary
+
+        Parameters
+        ----------
+
+        :return: The instances key from the dictionary
+        :rtype: Any
+
+        """
+        for key, instance in self.instances.items():
+            if id(instance) == id(self):
+                return key
+        return None
+
+    def _set_key(self, new_key: Any):
+        """
+        Protected method that sets a new key for the instance.
+        Iterates over the dictionary and changes the name.
+
+        If `None` is given, it will attribute the next available key
+
+        Parameters
+        ----------
+
+        :param new_key: The new instance key"
+        :type new_key: Any
+
+        """
+        cur_key = self._get_key()
+        if new_key is None:
+            new_key = self._next_available_key()
+        elif cur_key == new_key:
+            return
+        elif new_key in self.instances:
+            msg = f"Key {new_key} is already used"
+            raise ValueError(msg)
+        elif cur_key is not None:
+            self.instances.pop(cur_key)
+        self.instances[new_key] = self
+
+
+class LabeledTracker(Tracker):
+    """
+    Labeled Base Tracker, to track instances of Curve and Node
+    """
+
+    @classmethod
+    def _next_available_key(cls) -> str:
+        """
+        Protected method that gives next key
+        """
+        index = 1
+        while index in cls.instances:
+            index += 1
+        return index
+
+    @property
+    def label(self) -> int:
+        """
+        Gives the instance label
+
+        :getter: Returns the instance's label
+        :setter: Attribuates a new label for instance
+        :type: str
+
+        """
+        return int(self._get_key())
+
+    @label.setter
+    def label(self, new_label: Union[int, None]):
+        if new_label is not None:
+            new_label = int(new_label)
+        self._set_key(new_label)
+
+
+class NamedTracker(Tracker):
+    """
+    Named Base Tracker, to track instances of Material, Geometry and Section
+    """
+
+    @classmethod
+    def _next_available_key(cls) -> str:
+        """
+        Protected method that gives next key
+        """
+        index = 1
+        while f"instance-{index}" in cls.instances:
+            index += 1
+        return f"instance-{index}"
+
+    @property
+    def name(self) -> str:
+        """
+        Gives the instance name
+
+        :getter: Returns the instance's name
+        :setter: Attribuates a new name for instance
+        :type: str
+
+        """
+        return str(self._get_key())
+
+    @name.setter
+    def name(self, new_name: Union[str, None]):
+        if new_name is not None:
+            new_name = str(new_name)
+        self._set_key(new_name)
+
+
+class ICurve(ABC):
+    """
+    Interface abstract curve to be parent of others.
+
+    This class serves as interface between the curves from others packaged
+    like nurbs.Curve, to this package, expliciting the minimum requirements
+    of a curve must have.
+    With this, it's possible to implement your own type of parametric curve
+    """
+
+    @classmethod
+    @abstractmethod
+    def from_dict(cls, dictionary: Dict) -> ICurve:
+        """
+        Gives a curve instance based on the given parameters
+        """
+        raise NotImplementedError
+
+    @abstractmethod
+    def to_dict(self) -> Dict:
+        """
+        Transforms a curve instance into a dictionary
+        """
+        raise NotImplementedError
+
+    @property
+    @abstractmethod
+    def limits(self) -> Tuple[float]:
+        """
+        Gives the curve's parametric interval
+
+        :getter: Returns the pair [a, b] in which curve is parametric defined
+        :type: Tuple[float]
+
+        """
+        raise NotImplementedError
+
+    @property
+    @abstractmethod
+    def knots(self) -> Tuple[float]:
+        """
+        Gives the curve's knots, in which the parametric interval is divided
+
+        :getter: Returns the knots that divides the curve's interval
+        :type: Tuple[float]
+
+        """
+        raise NotImplementedError
+
+    @abstractmethod
+    def eval(self, parameters: Tuple[float]) -> Tuple[Tuple[float]]:
+        """
+        Evaluates the curves at given parameters.
+
+        :param parameters: A vector-like of lenght n
+        :type parameters: Tuple[float]
+        :return: A matrix of shape (n, 2)
+        :rtype: Tuple[Tuple[float]]
+        """
+        raise NotImplementedError
+
+    @abstractmethod
+    def winding(self, point: Tuple[float]) -> float:
+        """
+        Computes the winding number of the given curve
+
+        The possible results are one in the interval [0, 1]
+
+        * 0 means the point is outside the internal curve
+        * 1 means the point is inside the internal curve
+        * in (0, 1) means it's on the boundary
+
+        :param point: A 2D-point
+        :type point: Tuple[float]
+        :return: The winding value with respect to the curve
+        :rtype: float
+
+        """
+        raise NotImplementedError
+
+
+class IMaterial(ABC):
+    """
+    Material abstract class, the parent of other more specific materials
+    """
+
+    @abstractmethod
+    def to_dict(self) -> Dict:
+        """
+        Converts the material to a dictionary
+        """
+        raise NotImplementedError
+
+    @classmethod
+    @abstractmethod
+    def from_dict(cls, dictionary: Dict) -> IMaterial:
+        """
+        Converts the dictionary in a material instance
+        """
+        raise NotImplementedError
+
+
+class IFileIO(ABC):
+    """
+    Abstract Reader class that serves as basic interface to read/save
+    informations from/to given file.
+    The file format is defined by child class.
+    """
+
+    @abstractmethod
+    def is_open(self) -> bool:
+        """
+        Tells if the reader is open
+        """
+        raise NotImplementedError
+
+    @abstractmethod
+    def open(self, mode: str = "r"):
+        """
+        Opens the file with given mode
+        """
+        raise NotImplementedError
+
+    @abstractmethod
+    def close(self):
+        """
+        Closes the file
+        """
+        raise NotImplementedError
+
+    @abstractmethod
+    def load_nodes(self):
+        """
+        Saves all the nodes from file into Node class
+        """
+        raise NotImplementedError
+
+    @abstractmethod
+    def load_curves(self):
+        """
+        Creates all the curves instances from file
+        """
+        raise NotImplementedError
+
+    @abstractmethod
+    def load_materials(self):
+        """
+        Creates all the materials instances from file
+        """
+        raise NotImplementedError
+
+    @abstractmethod
+    def load_sections(self):
+        """
+        Creates all the sections instances from file
+        """
+        raise NotImplementedError
+
+
+class IField(ABC):
+    """
+    This is a base abstract class parent of others
+
+    It's responsible to decide if a given point is
+    inside/outside of given section.
+    """
+
+    @property
+    @abstractmethod
+    def ndata(self) -> int:
+        """
+        Gives the quantity of output numbers for one point
+
+        :getter: Returns the lenght of output data
+        :type: int
+        """
+        raise NotImplementedError
+
+    @abstractmethod
+    def eval(self, points: Tuple[Tuple[float]]) -> Tuple[Tuple[float]]:
+        """
+        Evaluate the field at given points
+
+        :param points: The wanted points, a matrix of shape (n, 2)
+        :type points: Tuple[Tuple[float]]
+        :return: The results in a matrix of shape (n, ndata)
+        :rtype: Tuple[Tuple[float]]
+        """
+        raise NotImplementedError
+
+
+class ISection(ABC):  # pylint: disable=too-few-public-methods
+    """
+    Section abstract class to serve as interface
+    """
diff --git a/src/compmec/section/curve.py b/src/compmec/section/curve.py
index 36df7cc..4c7fb87 100644
--- a/src/compmec/section/curve.py
+++ b/src/compmec/section/curve.py
@@ -6,18 +6,19 @@
 
 from __future__ import annotations
 
-from abc import ABC, abstractmethod
 from collections import OrderedDict
 from typing import Dict, Optional, Tuple
 
 import numpy as np
 from compmec.shape import JordanCurve
-from compmec.shape.shape import DefinedShape
 
 from compmec import nurbs
 
+from . import integral
+from .abcs import ICurve, LabeledTracker
 
-class Node:
+
+class Node(LabeledTracker):
     """
     Class that stores all nodes
     """
@@ -27,18 +28,6 @@ class Node:
     def __new__(cls, label: int) -> Tuple[float]:
         return Node.instances[label]
 
-    @staticmethod
-    def clear(labels: Optional[Tuple[int]] = None):
-        """
-        Remove node labels
-        """
-        if labels is None:
-            Node.instances.clear()
-            return
-        for label in labels:
-            if label in Node.instances:
-                Node.instances.pop(label)
-
     @staticmethod
     def insert_matrix(matrix: Tuple[Tuple[int, float, float]]):
         """
@@ -81,30 +70,13 @@ def from_labels(labels: Tuple[int]) -> Tuple[Tuple[float]]:
         return tuple(Node.instances[label] for label in labels)
 
 
-class Curve(ABC):
+class Curve(LabeledTracker, ICurve):
     """
-    Abstract curve to be parent of others.
-
-    This class serves as interface between the curves from others packaged
-    like nurbs.Curve, to this package, expliciting the minimum requirements
-    of a curve must have.
-    With this, it's possible to implement your only type of parametric curve
+    Base class that tracks the instances
     """
 
     instances = OrderedDict()
 
-    @staticmethod
-    def clear(labels: Optional[Tuple[int]] = None):
-        """
-        Removes all given instances of Curve
-        """
-        if labels is None:
-            Curve.instances.clear()
-            return
-        for label in labels:
-            if label in Curve.instances:
-                Curve.instances.pop(label)
-
     @staticmethod
     def new_instance(tipo: str, dictionary: Dict) -> Curve:
         """
@@ -121,105 +93,6 @@ def new_instance(tipo: str, dictionary: Dict) -> Curve:
             raise NotImplementedError
         return tipo2class[tipo].from_dict(dictionary)
 
-    @staticmethod
-    def __next_available_label() -> str:
-        index = 1
-        while index in Curve.instances:
-            index += 1
-        return index
-
-    def __new__(cls, label: Optional[int] = None) -> Curve:
-        if label is None:
-            label = Curve.__next_available_label()
-        elif label in Curve.instances:
-            msg = f"Cannot create curve '{label}'! There's already one!"
-            raise ValueError(msg)
-        instance = super().__new__(cls)
-        instance.label = label
-        Curve.instances[label] = instance
-        return instance
-
-    @classmethod
-    @abstractmethod
-    def from_dict(cls, dictionary: Dict) -> Curve:
-        """
-        Gives a curve instance based on the given parameters
-        """
-        raise NotImplementedError
-
-    @abstractmethod
-    def to_dict(self) -> Dict:
-        """
-        Transforms a curve instance into a dictionary
-        """
-        raise NotImplementedError
-
-    @property
-    def label(self) -> int:
-        """
-        Gives the curve label
-
-        :getter: Returns the curve's label
-        :setter: Attribuates a new label for curve
-        :type: int
-
-        """
-        return self.__label
-
-    @label.setter
-    def label(self, new_label: int):
-        try:
-            cur_label = self.label
-        except AttributeError:
-            self.__label = new_label
-            return
-        if cur_label == new_label:
-            return
-        if new_label in self.instances:
-            msg = f"Cannot set the label '{new_label}' for the curve "
-            msg += f"'{cur_label}' cause there's already a curve with "
-            msg += "the same label!\n"
-            msg += f"Cur: '{self.instances[new_label]}'\n"
-            msg += f"New: '{self}'"
-            raise ValueError(msg)
-        self.instances[new_label] = self.instances.pop(cur_label)
-
-    @property
-    @abstractmethod
-    def limits(self) -> Tuple[float]:
-        """
-        Gives the curve's parametric interval
-
-        :getter: Returns the pair [a, b] in which curve is parametric defined
-        :type: Tuple[float]
-
-        """
-        raise NotImplementedError
-
-    @property
-    @abstractmethod
-    def knots(self) -> Tuple[float]:
-        """
-        Gives the curve's knots, in which the parametric interval is divided
-
-        :getter: Returns the knots that divides the curve's interval
-        :type: Tuple[float]
-
-        """
-        raise NotImplementedError
-
-    @abstractmethod
-    def eval(self, parameters: Tuple[float]) -> Tuple[Tuple[float]]:
-        """
-        Evaluates the curves at given parameters.
-
-        :param parameters: A vector-like of lenght n
-        :type parameters: Tuple[float]
-        :return: A matrix of shape (n, 2)
-        :rtype: Tuple[Tuple[float]]
-        """
-        raise NotImplementedError
-
 
 class NurbsCurve(Curve):
     """
@@ -285,12 +158,14 @@ def __new__(cls, nurbs_curve: nurbs.Curve, label: Optional[int] = None):
         if not isinstance(nurbs_curve, nurbs.Curve):
             msg = "Invalid internal curve"
             raise TypeError(msg)
-        instance = super().__new__(cls, label)
+        instance = super().__new__(cls)
+        instance.label = label
         instance.internal = nurbs_curve
         return instance
 
     def eval(self, parameters: Tuple[float]) -> Tuple[Tuple[float]]:
-        return self.internal.eval(parameters)
+        values = self.internal.eval(parameters)
+        return np.array(values, dtype="float64")
 
     @property
     def knots(self) -> Tuple[float]:
@@ -318,27 +193,21 @@ def internal(self, new_curve: nurbs.Curve):
             raise TypeError(msg)
         self.__internal = new_curve
 
-
-def shapes_to_curves(shapes: Tuple[DefinedShape]) -> Tuple[Tuple[int]]:
-    """
-    Transform shapes instances into the pair (curves, geom_labels)
-    curves contains all the parametric curves used to define the shapes
-    geom_labels relates which curves are used to describe each shape
-
-    :param shapes: The group of shapes used in section
-    :type shapes: Tuple[DefinedShape]
-    :return: The pair (curves, geom_labels)
-    :rtype: Tuple[Tuple[int]]
-    """
-    geom_labels = []
-    for shape in shapes:
-        assert isinstance(shape, DefinedShape)
-        new_geom_labels = []
-        for jordan in shape.jordans:
-            signal = 1 if float(jordan) > 0 else -1
-            if signal < 0:
-                jordan = ~jordan
-            curve = NurbsCurve.from_jordan(jordan)
-            new_geom_labels.append(signal * curve.label)
-        geom_labels.append(new_geom_labels)
-    return geom_labels
+    def winding(self, point: Tuple[float]) -> float:
+        # Verify if the point is at any vertex
+        vertices = self.eval(self.knots[:-1])
+        for i, vertex in enumerate(vertices):
+            if np.all(point == vertex):
+                vec_left = vertices[(i - 1) % len(vertices)] - point
+                vec_righ = vertices[(i + 1) % len(vertices)] - point
+                wind = 0.5 * np.arccos(np.inner(vec_left, vec_righ)) / np.pi
+                return wind
+
+        wind = 0
+        for vertexa, vertexb in zip(vertices, np.roll(vertices, -1, axis=0)):
+            sub_wind = integral.winding_number_linear(vertexa, vertexb, point)
+            if abs(sub_wind) == 0.5:
+                wind = 0.5
+                break
+            wind += sub_wind
+        return wind
diff --git a/src/compmec/section/dataio.py b/src/compmec/section/dataio.py
index 9d0245b..cc4c9a9 100644
--- a/src/compmec/section/dataio.py
+++ b/src/compmec/section/dataio.py
@@ -6,18 +6,19 @@
 
 import json
 import os
-from abc import ABC, abstractmethod
 from pathlib import Path
 from typing import Dict, Optional
 
 import jsonschema
 
+from .abcs import IFileIO
 from .curve import Curve, Node
+from .geometry import Geometry
 from .material import Material
 from .section import Section
 
 
-class FileIO(ABC):
+class FileIO(IFileIO):
     """
     Abstract Reader class that serves as basic interface to read/save
     informations from/to given file.
@@ -41,64 +42,6 @@ def filepath(self) -> str:
         """
         return self.__filepath
 
-    @abstractmethod
-    def is_open(self) -> bool:
-        """
-        Tells if the reader is open
-        """
-        raise NotImplementedError
-
-    @abstractmethod
-    def open(self, mode: str = "r"):
-        """
-        Opens the file
-        """
-        raise NotImplementedError
-
-    @abstractmethod
-    def close(self):
-        """
-        Closes the file
-        """
-        raise NotImplementedError
-
-    @abstractmethod
-    def read_nodes(self):
-        """
-        Saves all the nodes from file into Node class
-        """
-        raise NotImplementedError
-
-    @abstractmethod
-    def read_curves(self):
-        """
-        Creates all the curves instances from file
-        """
-        raise NotImplementedError
-
-    @abstractmethod
-    def read_materials(self):
-        """
-        Creates all the materials instances from file
-        """
-        raise NotImplementedError
-
-    @abstractmethod
-    def read_sections(self):
-        """
-        Creates all the sections instances from file
-        """
-        raise NotImplementedError
-
-    def read(self):
-        """
-        Read the file and create all instaces
-        """
-        self.read_nodes()
-        self.read_curves()
-        self.read_materials()
-        self.read_sections()
-
     def __enter__(self):
         self.open()
         return self
@@ -171,7 +114,7 @@ def save_json(self, dictionary: Dict):
         with open(self.filepath, "w", encoding="ascii") as file:
             file.write(json_object)
 
-    def read_nodes(self):
+    def load_nodes(self):
         """
         Reads a nodes json and returns the data inside it.
 
@@ -191,7 +134,7 @@ def read_nodes(self):
             Node.clear(labels)
         Node.insert_matrix(matrix)
 
-    def read_curves(self):
+    def load_curves(self):
         """
         Reads a curve json and returns the data inside it.
 
@@ -213,7 +156,7 @@ def read_curves(self):
             curve = Curve.new_instance("nurbs", infos)
             curve.label = label
 
-    def read_materials(self):
+    def load_materials(self):
         """
         Reads a material json and returns the data inside it.
 
@@ -233,7 +176,7 @@ def read_materials(self):
             material = Material.new_instance("isotropic", info)
             material.name = name
 
-    def read_sections(self):
+    def load_sections(self):
         """
         Reads a section json file and returns the data inside it.
 
@@ -253,5 +196,9 @@ def read_sections(self):
         if self.overwrite:
             Section.clear(sections.keys())
         for name, info in sections.items():
-            section = Section.from_dict(info)
+            geome_labels = info["geom_labels"]
+            geometries = tuple(map(Geometry, geome_labels))
+            geome_names = tuple(geome.name for geome in geometries)
+            mater_names = tuple(info["materials"])
+            section = Section(geome_names, mater_names)
             section.name = name
diff --git a/src/compmec/section/field.py b/src/compmec/section/field.py
new file mode 100644
index 0000000..c9e8874
--- /dev/null
+++ b/src/compmec/section/field.py
@@ -0,0 +1,203 @@
+"""
+This file contains classes responsibles to evaluate some
+specific values within the section, at any point.
+
+For example, ChargedField is responsible to compute the
+stress and strain of the section for every
+"""
+
+from typing import Tuple
+
+import numpy as np
+
+from .abcs import IField, ISection
+
+
+class Field(IField):
+    """
+    This is a base abstract class parent of others
+
+    It's responsible to decide if a given point is
+    inside/outside of given section.
+    """
+
+    def __call__(self, points: Tuple[Tuple[float]]) -> Tuple[Tuple[float]]:
+        try:
+            return self.eval(points)
+        except TypeError:
+            stress, strain = self.eval([points])
+            return stress[0], strain[0]
+
+
+class ChargedField(Field):
+    """
+    A stress and strain field evaluator when
+    given forces and momentums are applied
+
+    This evaluator returns the values
+    (S13, S23, S33, E33, E13, E23, E33, E11, E22)
+    """
+
+    def __init__(self, section: ISection):
+        self.section = section
+        self.__charges = np.zeros(6, dtype="float64")
+
+    @property
+    def ndata(self):
+        return 8
+
+    @property
+    def forces(self):
+        """
+        Forces used in the charged field.
+
+        * Fx: shear force
+        * Fy: shear force
+        * Fz: axial force
+
+        :getter: Returns the (Fx, Fy, Fz)
+        :setter: Sets the new forces
+        :type: Tuple[float]
+        """
+        return self.__charges[:3]
+
+    @property
+    def momentums(self):
+        """
+        Forces used in the charged field.
+
+        * Mx: bending momentum
+        * My: bending momentum
+        * Mz: torsion momentum
+
+        :getter: Returns the (Mx, My, Mz)
+        :setter: Sets the new momentums
+        :type: Tuple[float]
+        """
+        return self.__charges[3:]
+
+    @forces.setter
+    def forces(self, new_forces: Tuple[float]):
+        self.__charges[:3] = new_forces
+
+    @momentums.setter
+    def momentums(self, new_momentums: Tuple[float]):
+        self.__charges[3:] = new_momentums
+
+    def __stress_axial(self, winds):
+        sigmazz = self.forces[2] / self.section.area()
+        return np.any(winds, axis=1) * sigmazz
+
+    def __stress_bending(self, points, winds):
+        bend_center = self.section.bending_center()
+        ixx, ixy, iyy = self.section.second_moment(bend_center)
+        detii = ixx * iyy - ixy**2
+        matrix = np.array([[iyy, -ixy], [-ixy, ixx]])
+        momx, momy, _ = self.momentums
+        vector = np.dot(matrix, [-momy, momx]) / detii
+        return np.dot(points, vector) * np.any(winds, axis=1)
+
+    def __winding_numbers(
+        self, points: Tuple[Tuple[float]]
+    ) -> Tuple[Tuple[float]]:
+        """
+        Computes the winding number of every point,
+        for every geometry
+
+        """
+        geometries = tuple(self.section.geometries)
+        winds = np.zeros((len(points), len(geometries)), dtype="float64")
+        for i, point in enumerate(points):
+            for j, geome in enumerate(geometries):
+                winds[i, j] = geome.winding(point)
+        return winds
+
+    def __stress_eval(self, points, winds):
+        """
+        Evaluates the stress values
+
+        The stress tensor is given by
+            [  0    0  E13]
+        S = [  0    0  E23]
+            [E13  E23  E33]
+
+        Returned values are a matrix of shape (n, 3)
+        each line are the stress components: E13, E23, E33
+
+        :param points: The wanted points, a matrix of shape (n, 2)
+        :type points: Tuple[Tuple[float]]
+        :return: The strain matrix of shape (n, 3)
+        :rtype: Tuple[Tuple[float]]
+
+        """
+        results = np.zeros((len(points), 3), dtype="float64")
+        if self.forces[2]:  # Axial force
+            results[:, 2] += self.__stress_axial(winds)
+        if np.any(self.momentums[:2]):  # Bending moments
+            results[:, 2] += self.__stress_bending(points, winds)
+        if np.any(self.forces[:2]):  # Shear force
+            raise NotImplementedError
+        if self.momentums[2]:  # Torsion
+            raise NotImplementedError
+        return results
+
+    def __strain_eval(self, winds, stresses):
+        """
+        Evaluates the strain values from stress values by
+        using Hook's law for isotropic materials
+
+        The winds serves to know the position of the points,
+        to decide which material will be used
+
+        The strain tensor is given by
+            [E11    0  E13]
+        E = [  0  E22  E23]
+            [E13  E23  E33]
+        The values E22 and E11 are the same
+
+        Returned values are a matrix of shape (n, 4)
+        each line are the strain components: E11, E33, E13, E23
+
+        :param points: The wanted points, a matrix of shape (n, 2)
+        :type points: Tuple[Tuple[float]]
+        :return: The strain matrix of shape (n, 4)
+        :rtype: Tuple[Tuple[float]]
+
+        """
+        strain = np.zeros((len(winds), 4), dtype="float64")
+        for i, material in enumerate(self.section.materials):
+            subwinds = winds[:, i]
+            mask = subwinds != 0
+            young = material.young_modulus
+            poiss = material.poissons_ratio
+            fract11 = -poiss / young
+            fract13 = (1 + poiss) / young
+            strain[mask, 1] += subwinds[mask] * stresses[mask, 2] / young
+            strain[mask, 0] += fract11 * subwinds[mask] * stresses[mask, 2]
+            strain[mask, 2] += fract13 * subwinds[mask] * stresses[mask, 0]
+            strain[mask, 3] += fract13 * subwinds[mask] * stresses[mask, 1]
+        return strain
+
+    def eval(self, points):
+        """
+        Evaluate the stress and strain at given points
+
+        The inputs/outputs are matrices:
+        * points is a (n, 2) matrix
+        * stress is a (n, 3) matrix
+        * strain is a (n, 4) matrix
+
+        The stress components are S13, S23, S33, meaning
+        two shear stresses and one normal stress
+        The strain components are E11, E33, E13, E23
+
+        :param points: The wanted points, a matrix of shape (n, 2)
+        :type points: Tuple[Tuple[float]]
+        :return: The pair (stress, strain)
+        :rtype: Tuple[Tuple[Tuple[float]]]
+        """
+        points = np.array(points, dtype="float64")
+        winds = self.__winding_numbers(points)
+        stress = self.__stress_eval(points, winds)
+        strain = self.__strain_eval(winds, stress)
+        return stress, strain
diff --git a/src/compmec/section/geometry.py b/src/compmec/section/geometry.py
new file mode 100644
index 0000000..99a1d2b
--- /dev/null
+++ b/src/compmec/section/geometry.py
@@ -0,0 +1,80 @@
+"""
+File that contains Geometry class, which defines a planar region
+"""
+
+from __future__ import annotations
+
+from typing import Optional, Tuple
+
+import numpy as np
+from compmec.shape.shape import ConnectedShape, DefinedShape, SimpleShape
+
+from .abcs import NamedTracker
+from .curve import Curve, NurbsCurve
+
+
+class Geometry(NamedTracker):
+    """
+    Connected Geometry class that represents a bidimensional region
+    on the plane and has some functions such as to decide if points
+    are inside the region
+    """
+
+    def __init__(self, curve_labels: Tuple[int], name: Optional[str] = None):
+        curve_labels = tuple(map(int, curve_labels))
+        for label in curve_labels:
+            assert abs(label) in Curve.instances
+        self.__curve_labels = curve_labels
+        self.name = name
+
+    @property
+    def labels(self) -> Tuple[int]:
+        """
+        Gives the curve labels that defines the geometry
+        """
+        return self.__curve_labels
+
+    def winding(self, point: Tuple[float]) -> float:
+        """
+        Computes the winding number of the giving point
+
+        Normally winding number is one for each curve, but this method gives
+        the overall winding number
+        """
+        wind_tolerance = 1e-9
+        labels = np.array(sorted(self.labels), dtype="int32")
+        winds = []
+        for label in labels:
+            curve = Curve.instances[abs(label)]
+            wind = curve.winding(point)
+            if wind_tolerance < wind < 1 - wind_tolerance:
+                return wind
+            winds.append(wind)
+        winds = np.array(winds, dtype="int64")
+        mask = (winds < 0.5) & (labels < 0)
+        mask |= (winds > 0.5) & (labels > 0)
+        return float(np.all(mask))
+
+
+def shapes2geometries(shapes: Tuple[DefinedShape]) -> Tuple[Geometry]:
+    """
+    Transform shapes instances into geometry instances
+
+    :param shapes: The shape to be converted
+    :type shapes: Union[SimpleShape, ConnectedShape]
+    :return: The geometry used
+    :rtype: Tuple[]
+    """
+    geometries = []
+    for shape in shapes:
+        assert isinstance(shape, (SimpleShape, ConnectedShape))
+        curve_labels = []
+        for jordan in shape.jordans:
+            signal = 1 if float(jordan) > 0 else -1
+            if signal < 0:
+                jordan = ~jordan
+            curve = NurbsCurve.from_jordan(jordan)
+            curve_labels.append(signal * curve.label)
+        geometry = Geometry(curve_labels)
+        geometries.append(geometry)
+    return geometries
diff --git a/src/compmec/section/integral.py b/src/compmec/section/integral.py
index 7a7a094..5b6cb35 100644
--- a/src/compmec/section/integral.py
+++ b/src/compmec/section/integral.py
@@ -46,51 +46,29 @@ def comb(a, b):
     return prod
 
 
-# pylint: disable=invalid-name
-def integrate_polygon(
-    xverts: Tuple[float], yverts: Tuple[float], amax: int = 3, bmax: int = 3
-) -> Tuple[Tuple[float]]:
+def winding_number_linear(
+    pointa: Tuple[float], pointb: Tuple[float], center: Tuple[float]
+) -> float:
     """
-    Computes integrals, returning a matrix II such
-
-    II_{a, b} = int_D x^a * y^b dx dy
-
-    for a = [0, ..., amax] and b = [0, ..., bmax]
-
-    and D being a polygonal domain given by the vertices
-
-    Parameters
-    ----------
-    xverts: tuple[float]
-        Abssice values of the vertices
-    yverts: tuple[float]
-        Ordoned values of the vertices
-    amax: int, default 3
-        The maximum expoent multiplying 'x'
-    bmax: int, default 3
-        The maximum expoent multiplying 'y'
-    return: tuple[tuple[float]]
-        A matrix II of shape (amax+1, bmax+1) such
-        II_{ij} = int_D x^a * y^b dx dy
+    Computes the winding number for a straight line AB with respect to C
+
+    :param pointa: The pair (xa, ya)
+    :type pointa: Tuple[float]
+    :param pointb: The pair (xb, yb)
+    :type pointb: Tuple[float]
+    :param center: The pair (xc, yc)
+    :type center: Tuple[float]
+    :return: The winding number value in the interval [-0.5, 0.5]
+    :rtype: float
     """
-    xvan = np.vander(xverts, amax + 1, True)
-    yvan = np.vander(yverts, bmax + 1, True)
-
-    cross = xverts * np.roll(yverts, shift=-1)
-    cross -= yverts * np.roll(xverts, shift=-1)
-
-    geomprops = np.zeros((amax + 1, bmax + 1), dtype="float64")
-    for a in range(amax + 1):
-        for b in range(bmax + 1):
-            M = np.zeros((a + 1, b + 1), dtype="float64")
-            for i in range(a + 1):
-                for j in range(b + 1):
-                    M[i, j] = comb(i + j, i) * comb(a + b - i - j, b - j)
-            X = np.roll(xvan[:, : a + 1], shift=-1, axis=0) * xvan[:, a::-1]
-            Y = np.roll(yvan[:, : b + 1], shift=-1, axis=0) * yvan[:, b::-1]
-            geomprops[a, b] = np.einsum("k,ki,ij,kj", cross, X, M, Y)
-            geomprops[a, b] /= (a + b + 2) * (a + b + 1) * comb(a + b, a)
-    return geomprops
+    x1 = float(pointa[0] - center[0])
+    y1 = float(pointa[1] - center[1])
+    x2 = float(pointb[0] - center[0])
+    y2 = float(pointb[1] - center[1])
+    cross = x1 * y2 - x2 * y1
+    inner = x1 * x2 + y1 * y2
+    wind = np.arctan2(cross, inner) / math.tau
+    return wind
 
 
 class Integration:
@@ -278,6 +256,34 @@ def closed(npts: int) -> Tuple[Tuple[float]]:
         weights = np.array(weights[npts - 2], dtype="float64")
         return nodes, weights
 
+    @staticmethod
+    def opened(npts: int) -> Tuple[Tuple[float]]:
+        """
+        Open newton cotes formula
+
+        Parameters
+        ----------
+        npts: int
+            The number of points to use gauss integration.
+            Must be at least 1
+        return: tuple[tuple[float]]
+            The pair (nodes, weights), with each node in the interval [0, 1]
+        """
+        if not isinstance(npts, int) or npts < 1:
+            raise ValueError(f"npts invalid: {npts}")
+        if npts > 7:
+            raise NotImplementedError
+        nodes = tuple(num / (npts + 1) for num in range(1, npts + 1))
+        nodes = np.array(nodes, dtype="float64")
+        weights = (
+            (1.0,),
+            (0.5, 0.5),
+            (2 / 3, -1 / 3, 2 / 3),
+            (11 / 24, 1 / 24, 1 / 24, 11 / 24),
+        )
+        weights = np.array(weights[npts - 1], dtype="float64")
+        return nodes, weights
+
     @staticmethod
     def chebyshev(npts: int) -> Tuple[Tuple[float]]:
         """Chebyshev integration
@@ -329,3 +335,193 @@ def chebyshev(npts: int) -> Tuple[Tuple[float]]:
         nodes = np.array(nodes, dtype="float64")
         weights = np.array(all_weights[npts - 1], dtype="float64")
         return nodes, weights
+
+
+class Polynomial:
+    """
+    Class responsible to compute the integral
+
+    II(a, b) = int_D x^a * y^b dx dy
+
+    Which D is the domain defined by the interior of a curve
+    The (a, b) values are the expoents
+
+    """
+
+    expoents = [
+        (0, 0),
+        (0, 1),
+        (1, 0),
+        (0, 2),
+        (1, 1),
+        (2, 0),
+        (0, 3),
+        (2, 1),
+        (1, 2),
+        (3, 0),
+    ]
+
+    # pylint: disable=invalid-name
+    @staticmethod
+    def polygon(vertices: Tuple[Tuple[float]]) -> Tuple[float]:
+        """
+        Computes integrals, returning the values of the integral
+
+        II_{a, b} = int_D x^a * y^b dx dy
+
+        for (a, b) = [(0, 0), (0, 1), (1, 0), (0, 2), (1, 1), (2, 0),
+                      (0, 3), (2, 1), (1, 2), (3, 0)]
+
+        and D being a polygonal domain given by the vertices
+
+        Parameters
+        ----------
+
+        :param vertices: The vertices that defines the polygon
+        :type vertices: tuple[tuple[float]]
+        :return: A vector of lenght 10 containing the integral values
+        :rtype: tuple[float]
+        """
+        vertices = np.array(vertices, dtype="float64")
+        xvan = np.vander(vertices[:, 0], 4, True)
+        yvan = np.vander(vertices[:, 1], 4, True)
+
+        cross = vertices[:, 0] * np.roll(vertices[:, 1], shift=-1)
+        cross -= vertices[:, 1] * np.roll(vertices[:, 0], shift=-1)
+
+        expoents = Polynomial.expoents
+        geomprops = np.zeros(len(expoents), dtype="float64")
+        for k, (a, b) in enumerate(expoents):
+            M = np.zeros((a + 1, b + 1), dtype="float64")
+            for i in range(a + 1):
+                for j in range(b + 1):
+                    M[i, j] = comb(i + j, i) * comb(a + b - i - j, b - j)
+            X = np.roll(xvan[:, : a + 1], shift=-1, axis=0) * xvan[:, a::-1]
+            Y = np.roll(yvan[:, : b + 1], shift=-1, axis=0) * yvan[:, b::-1]
+            geomprops[k] = np.einsum("k,ki,ij,kj", cross, X, M, Y)
+            geomprops[k] /= (a + b + 2) * (a + b + 1) * comb(a + b, a)
+        return geomprops
+
+    @staticmethod
+    def adaptative(curve, tolerance: float = 1e-9) -> Tuple[float]:
+        """
+        Computes the polynomials integrals over the area defined by the curve
+
+        It uses an adaptative algorithm that allows computing the integrals
+        over smooth curves using milne's (open newton quadrature 3 points)
+
+        This function uses a recursive approach
+
+        :param curve: The boundary curve around the area
+        :type curve: Curve
+        :return: The integral values
+        :rtype: float
+
+        """
+        expoents = Polynomial.expoents
+        integrals = np.zeros(len(expoents), dtype="float64")
+        integrator = AdaptativePolynomialIntegrator(curve, tolerance=tolerance)
+        for k, (expx, expy) in enumerate(expoents):
+            value = integrator.integrate(expx, expy)
+            integrals[k] = value / (2 + expx + expy)
+        return integrals
+
+
+class AdaptativePolynomialIntegrator:
+    """
+    Adaptative Polynomial Integrator
+
+    Receives a curve and computes the polynomials integrals recursivelly
+    using open newton cotes quadrature with 3 points
+    """
+
+    integ_matrix = (
+        np.array([[0, 4, -1], [-4, 0, 4], [1, -4, 0]], dtype="float64") / 3
+    )
+
+    def __init__(self, curve, *, max_depth: int = 9, tolerance=1e-9):
+        self.curve = curve
+        self.max_depth = max_depth
+
+        mesh = set()
+        for t0, t1 in zip(curve.knots, curve.knots[1:]):
+            subtol = tolerance * (t1 - t0) / (curve.knots[-1] - curve.knots[0])
+            mesh |= self.find_mesh(t0, t1, tolerance=subtol)
+        self.mesh = tuple(sorted(mesh))
+
+    def find_mesh(
+        self, t0: float, t1: float, tolerance: float = 1e-9
+    ) -> Tuple[float]:
+        """
+        Finds the curve's subdivisions such the area computed
+        by the adaptative subdivision is lower than the tolerance
+        """
+        subknots = np.linspace(t0, t1, 5)
+        xvals, yvals = np.transpose(self.curve.eval(subknots))
+        area_mid = xvals[::2] @ self.integ_matrix @ yvals[::2]
+        area_lef = xvals[:3] @ self.integ_matrix @ yvals[:3]
+        area_rig = xvals[2:] @ self.integ_matrix @ yvals[2:]
+        diff = abs(area_lef + area_rig - area_mid)
+        mesh = {t0, t1}
+        if diff > tolerance:
+            mesh |= self.find_mesh(t0, subknots[2], tolerance / 2)
+            mesh |= self.find_mesh(subknots[2], t1, tolerance / 2)
+        return mesh
+
+    def milne_formula(self, t0: float, t1: float, a: int, b: int) -> float:
+        """
+        Computes the integral using Milne formula
+
+        int_{t0}^{t1} x^a * y^b * (x * dy - y * dx)
+
+        with (x(t), y(t)) being the curve
+
+        Parameters
+        ----------
+
+        :param t0: The lower interval's value
+        :type t0: float
+        :param t1: The upper interval's value
+        :type t1: float
+        :param a: The x expoent
+        :type a: int
+        :param b: The y expoent
+        :type b: int
+        :return: The interval's integral value
+        :rtype: float
+
+        """
+        ts = np.linspace(t0, t1, 5)
+        points = self.curve.eval(ts)
+        xs = points[:, 0]
+        ys = points[:, 1]
+        f1 = xs[1] ** a * ys[1] ** b
+        f2 = xs[2] ** a * ys[2] ** b
+        f3 = xs[3] ** a * ys[3] ** b
+        f1 *= xs[1] * (ys[2] - ys[0]) - ys[1] * (xs[2] - xs[0])
+        f2 *= xs[2] * (ys[3] - ys[1]) - ys[2] * (xs[3] - xs[1])
+        f3 *= xs[3] * (ys[4] - ys[2]) - ys[3] * (xs[4] - xs[2])
+        return (4 * f1 - 2 * f2 + 4 * f3) / 3
+
+    def integrate(self, expx: int, expy: int) -> float:
+        """
+        Computes the integral over the entire curve
+
+        int_{Gamma} x^expx * y^expy * (x * dy - y * dx)
+
+        with (x(t), y(t)) being the curve
+
+        Parameters
+        ----------
+
+        :param tolerance: The tolerance to how when stop adaptative quadrature
+        :type tolerance: float
+        :return: The interval's integral value
+        :rtype: float
+        """
+        value = 0
+        for t0, t1 in zip(self.mesh, self.mesh[1:]):
+            tm = 0.5 * (t0 + t1)
+            value += self.milne_formula(t0, tm, expx, expy)
+            value += self.milne_formula(tm, t1, expx, expy)
+        return value
diff --git a/src/compmec/section/material.py b/src/compmec/section/material.py
index 9184d69..8989877 100644
--- a/src/compmec/section/material.py
+++ b/src/compmec/section/material.py
@@ -4,30 +4,19 @@
 
 from __future__ import annotations
 
-from abc import ABC, abstractmethod
 from collections import OrderedDict
-from typing import Dict, Optional, Tuple
+from typing import Dict, Optional
 
+from .abcs import IMaterial, NamedTracker
 
-class Material(ABC):
+
+class Material(IMaterial, NamedTracker):
     """
     Material abstract class, the parent of other more specific materials
     """
 
     instances = OrderedDict()
 
-    @staticmethod
-    def clear(names: Optional[Tuple[str]] = None):
-        """
-        Removes all given instances of Curve
-        """
-        if names is None:
-            Material.instances.clear()
-            return
-        for name in names:
-            if name in Material.instances:
-                Material.instances.pop(name)
-
     @staticmethod
     def new_instance(tipo: str, dictionary: Dict) -> Material:
         """
@@ -44,64 +33,6 @@ def new_instance(tipo: str, dictionary: Dict) -> Material:
             raise NotImplementedError
         return tipo2class[tipo].from_dict(dictionary)
 
-    @staticmethod
-    def __next_available_name() -> str:
-        index = 1
-        while True:
-            name = f"custom-material-{index}"
-            if name not in Material.instances:
-                return name
-            index += 1
-
-    def __init__(self, name: Optional[str] = None):
-        if name is None:
-            name = Material.__next_available_name()
-        elif name in Material.instances:
-            msg = f"Cannot create material '{name}'! There's already one!"
-            raise ValueError(msg)
-        self.__name = name
-        self.__class__.instances[name] = self
-
-    @abstractmethod
-    def to_dict(self) -> Dict:
-        """
-        Converts the material to a dictionary
-        """
-        raise NotImplementedError
-
-    @classmethod
-    @abstractmethod
-    def from_dict(cls, dictionary: Dict) -> Material:
-        """
-        Converts the dictionary in a material instance
-        """
-        raise NotImplementedError
-
-    @property
-    def name(self) -> str:
-        """
-        Gives the material name
-
-        :getter: Returns the material's name
-        :setter: Attribuates a new name for material
-        :type: str
-
-        """
-        return self.__name
-
-    @name.setter
-    def name(self, new_name: str):
-        if self.name == new_name:
-            return
-        if new_name in self.instances:
-            msg = f"Cannot set the name '{new_name}' for the material "
-            msg += f"'{self.name}' cause there's already a material with "
-            msg += "the same name!\n"
-            msg += f"Cur: '{self.instances[new_name]}'\n"
-            msg += f"New: '{self}'"
-            raise ValueError(msg)
-        self.instances[new_name] = self.instances.pop(self.name)
-
 
 class Isotropic(Material):
     """
@@ -147,9 +78,9 @@ def __init__(
     ):
         Isotropic.__verify_young_modulus(young_modulus)
         Isotropic.__verify_poissons_ratio(poissons_ratio)
-        super().__init__(name)
         self.__young_modulus = young_modulus
         self.__poissons_ratio = poissons_ratio
+        self.name = name
 
     def __str__(self) -> str:
         values = [
diff --git a/src/compmec/section/section.py b/src/compmec/section/section.py
index 008d8cd..d59bb55 100644
--- a/src/compmec/section/section.py
+++ b/src/compmec/section/section.py
@@ -9,17 +9,20 @@
 from __future__ import annotations
 
 from collections import OrderedDict
-from typing import Dict, Optional, Tuple, Union
+from typing import Iterable, Optional, Tuple, Union
 
 import numpy as np
 from compmec.shape.shape import DefinedShape
 
-from .curve import Curve, shapes_to_curves
-from .integral import integrate_polygon
+from .abcs import ISection, NamedTracker
+from .curve import Curve
+from .field import ChargedField
+from .geometry import Geometry, shapes2geometries
+from .integral import Polynomial
 from .material import Material
 
 
-class BaseSection:
+class BaseSection(ISection, NamedTracker):
     """
     BaseSection class that is the base for others section classes
 
@@ -30,106 +33,47 @@ class BaseSection:
 
     instances = OrderedDict()
 
-    @staticmethod
-    def clear(names: Optional[Tuple[str]] = None):
-        """
-        Removes all given instances of Curve
-        """
-        if names is None:
-            BaseSection.instances.clear()
-            return
-        for name in names:
-            if name in BaseSection.instances:
-                BaseSection.instances.pop(name)
-
-    @staticmethod
-    def __next_available_name() -> str:
-        index = 1
-        while True:
-            name = f"custom-section-{index}"
-            if name not in Material.instances:
-                return name
-            index += 1
-
     def __init__(
         self,
-        geom_labels: Tuple[Tuple[int]],
+        geome_names: Tuple[str],
         mater_names: Tuple[str],
         name: Optional[str] = None,
     ):
-        for labels in geom_labels:
-            for label in labels:
-                assert abs(label) in Curve.instances
+        for geo_name in geome_names:
+            assert geo_name in Geometry.instances
         for mat_name in mater_names:
             assert mat_name in Material.instances
-        if name is None:
-            name = BaseSection.__next_available_name()
-        elif name in BaseSection.instances:
-            raise ValueError
-        self.__geom_labels = tuple(
-            tuple(map(int, labels)) for labels in geom_labels
-        )
+        self.name = name
+        self.__geome_names = tuple(geome_names)
         self.__mater_names = tuple(mater_names)
-        self.__name = name
-        self.instances[name] = self
 
     @property
-    def name(self) -> str:
-        """
-        Gives the material name
-
-        :getter: Returns the material's name
-        :setter: Attribuates a new name for material
-        :type: str
-
-        """
-        return self.__name
-
-    @name.setter
-    def name(self, new_name: str):
-        if self.name == new_name:
-            return
-        if new_name in self.instances:
-            msg = f"Section name '{new_name}' is already used"
-            raise ValueError(msg)
-        self.instances[new_name] = self.instances.pop(self.name)
-        self.__name = new_name
-
-    @property
-    def geom_labels(self) -> Tuple[Tuple[int]]:
+    def geometries(self) -> Iterable[Geometry]:
         """
         Geometric curves labels that defines shapes
 
         :getter: Returns the curve labels
         :type: Tuple[Tuple[int]]
         """
-        return self.__geom_labels
+        for name in self.__geome_names:
+            yield Geometry.instances[name]
 
     @property
-    def mater_names(self) -> Tuple[str]:
-        """
-        Material names in each shape
-
-        :getter: Returns the material names
-        :type: Tuple[str]
-        """
-        return self.__mater_names
-
-    @property
-    def materials(self) -> Tuple[Material]:
+    def materials(self) -> Iterable[Material]:
         """
         Used materials for every shape
 
         :getter: Returns the used materials, in the shapes' order
         :type: Tuple[Material]
         """
-        return tuple(Material.instances[name] for name in self.mater_names)
+        for name in self.__mater_names:
+            yield Material.instances[name]
 
     @classmethod
     def from_shapes(
         cls,
         shapes: Union[DefinedShape, Tuple[DefinedShape]],
-        materials=Union[Material, Tuple[Material]],
+        materials: Union[Material, Tuple[Material]],
     ) -> BaseSection:
         """
         Creates an Section instance based on given shapes and materials
@@ -144,23 +88,9 @@ def from_shapes(
             shapes = [shapes]
         if isinstance(materials, Material):
             materials = [materials]
-        geom_labels = shapes_to_curves(shapes)
+        geome_names = tuple(geom.name for geom in shapes2geometries(shapes))
         mater_names = tuple(material.name for material in materials)
-        return cls(geom_labels, mater_names)
-
-    @classmethod
-    def from_dict(cls, dictionary: Dict) -> BaseSection:
-        """
-        Transforms the dictionary, read from json, into a section instance
-
-        :param dictionary: The informations that defines the section
-        :type dictionary: Dict
-        :return: The created section instance
-        :rtype: BaseSection
-        """
-        geom_labels = dictionary["geom_labels"]
-        mater_names = dictionary["materials"]
-        return cls(geom_labels, mater_names)
+        return cls(geome_names, mater_names)
 
 
 class GeometricSection(BaseSection):
@@ -171,10 +101,10 @@ class GeometricSection(BaseSection):
 
     def __init__(
         self,
-        shapes: Tuple[Tuple[int]],
-        materials: Tuple[Material],
+        geome_names: Tuple[Geometry],
+        mater_names: Tuple[Material],
     ):
-        super().__init__(shapes, materials)
+        super().__init__(geome_names, mater_names)
         self.__geomintegs = None
 
     def __compute_geomintegs(self):
@@ -183,30 +113,18 @@ def __compute_geomintegs(self):
         creating the object __geomintegs
         """
         integrals = {}
-        all_labels = {val for labels in self.geom_labels for val in labels}
+        all_labels = set()
+        for geometry in self.geometries:
+            all_labels |= set(map(abs, geometry.labels))
         for label in all_labels:
             curve = Curve.instances[label]
-            vertices = curve.eval(curve.knots)
-            vertices = np.array(vertices, dtype="float64")
-            xverts, yverts = np.transpose(vertices)
-            integrals[label] = integrate_polygon(xverts, yverts)
-
-        integral = np.zeros((4, 4), dtype="float64")
-        for labels in self.geom_labels:
-            for label in labels:
+            integrals[label] = Polynomial.adaptative(curve)
+
+        geomintegs = np.zeros(10, dtype="float64")
+        for geometry in self.geometries:
+            for label in geometry.labels:
                 signal = 1 if label > 0 else -1
-                integral += signal * integrals[abs(label)]
-
-        geomintegs = [integral[0, 0]]
-        geomintegs += [integral[0, 1], integral[1, 0]]
-        geomintegs += [integral[0, 2], integral[1, 1], integral[2, 0]]
-        geomintegs += [
-            integral[0, 3],
-            integral[1, 2],
-            integral[2, 1],
-            integral[3, 0],
-        ]
-        geomintegs = tuple(map(float, geomintegs))
+                geomintegs += signal * integrals[abs(label)]
         self.__geomintegs = geomintegs
 
     def area(self) -> float:
@@ -368,6 +286,16 @@ def gyradius(self) -> Tuple[float]:
         ixx, _, iyy = self.second_moment()
         return np.sqrt(iyy / area), np.sqrt(ixx / area)
 
+    def charged_field(self) -> ChargedField:
+        """
+        Gives the charged field instance to evaluate stresses
+
+        :return: The field evaluator
+        :rtype: ChargedField
+
+        """
+        return ChargedField(self)
+
 
 class Section(GeometricSection):
     """
diff --git a/tests/json/steel_circle.json b/tests/json/steel_circle.json
index 86ad76f..da34a31 100644
--- a/tests/json/steel_circle.json
+++ b/tests/json/steel_circle.json
@@ -10,7 +10,8 @@
     "curves": {
         "2000": {"degree": 3,
                  "knotvector": [0, 0, 0, 0, 0.5, 0.5, 0.5, 1, 1, 1, 1],
-                 "ctrllabels": [1, 2, 3, 4, 5, 6, 1]}
+                 "ctrllabels": [1, 2, 3, 4, 5, 6, 1],
+                 "weights": [3, 1, 1, 3, 1, 1, 3]}
     },
     "materials": {
         "steel": {
diff --git a/tests/test_axial.py b/tests/test_axial.py
index baa114e..e354f9a 100644
--- a/tests/test_axial.py
+++ b/tests/test_axial.py
@@ -11,7 +11,6 @@
 
 
 @pytest.mark.order(4)
-@pytest.mark.skip(reason="Test only in the version 0.4.0")
 @pytest.mark.dependency(
     depends=[
         "tests/test_material.py::test_end",
@@ -36,11 +35,32 @@ def test_centered_square(self):
         side = 2
         geometry = Primitive.square(side)
         material = Isotropic(young_modulus=210e3, poissons_ratio=0.3)
-        section = Section(geometry, material)
+        section = Section.from_shapes(geometry, material)
         field = section.charged_field()
-        points = [(0, 0)]
-        values = field.eval(points)
-        assert np.all(np.abs(values) < 1e-9)
+
+        points = [(0, 0)]  # origin
+        points += [
+            (1, 0),
+            (1, 1),
+            (0, 1),
+            (-1, 1),
+            (-1, 0),
+            (-1, -1),
+            (0, -1),
+            (1, -1),
+        ]
+        field.forces = (0, 0, 0)
+        field.momentums = (0, 0, 0)
+        stress, strain = field.eval(points)
+        assert np.all(np.abs(stress) < 1e-9)  # no charge, all zero
+        assert np.all(np.abs(strain) < 1e-9)  # no charge, all zero
+
+        field.forces = (0, 0, 4)
+        field.momentums = (0, 0, 0)
+        stress, strain = field.eval(points)
+        assert np.all(np.abs(stress[:, :2]) < 1e-9)  # no shear
+        abs_diff = np.abs(stress[:, 2] - 1)
+        assert np.all(abs_diff < 1e-9)
 
     @pytest.mark.order(4)
     @pytest.mark.dependency(
@@ -66,11 +86,33 @@ def test_centered_square(self):
         geometry = Primitive.square(ext_side)
         geometry -= Primitive.square(int_side)
         material = Isotropic(young_modulus=210e3, poissons_ratio=0.3)
-        section = Section(geometry, material)
+        section = Section.from_shapes(geometry, material)
         field = section.charged_field()
+
         points = [(0, 0)]
-        values = field.eval(points)
-        assert np.all(np.abs(values) < 1e-9)
+        points += [
+            (1, 0),
+            (1, 1),
+            (0, 1),
+            (-1, 1),
+            (-1, 0),
+            (-1, -1),
+            (0, -1),
+            (1, -1),
+        ]
+        field.forces = (0, 0, 0)
+        field.momentums = (0, 0, 0)
+        stress, strain = field.eval(points)
+        assert np.all(np.abs(stress) < 1e-9)  # no charge, all zero
+        assert np.all(np.abs(strain) < 1e-9)  # no charge, all zero
+
+        field.forces = (0, 0, 3)
+        field.momentums = (0, 0, 0)
+        stress, strain = field.eval(points)
+        assert np.all(np.abs(stress[:, :2]) < 1e-9)  # no shear
+        good_normal_stress = (0, 1, 1, 1, 1, 1, 1, 1, 1)
+        abs_diff = np.abs(stress[:, 2] - good_normal_stress)
+        assert np.all(abs_diff < 1e-9)
 
     @pytest.mark.order(4)
     @pytest.mark.dependency(
diff --git a/tests/test_bending.py b/tests/test_bending.py
index da26746..53294bb 100644
--- a/tests/test_bending.py
+++ b/tests/test_bending.py
@@ -11,7 +11,6 @@
 
 
 @pytest.mark.order(4)
-@pytest.mark.skip(reason="Test only in the version 0.4.0")
 @pytest.mark.dependency(
     depends=[
         "tests/test_material.py::test_end",
@@ -36,11 +35,43 @@ def test_centered_square(self):
         side = 2
         geometry = Primitive.square(side)
         material = Isotropic(young_modulus=210e3, poissons_ratio=0.3)
-        section = Section(geometry, material)
+        section = Section.from_shapes(geometry, material)
         field = section.charged_field()
+
         points = [(0, 0)]
-        values = field.eval(points)
-        assert np.all(np.abs(values) < 1e-9)
+        points += [
+            (1, 0),
+            (1, 1),
+            (0, 1),
+            (-1, 1),
+            (-1, 0),
+            (-1, -1),
+            (0, -1),
+            (1, -1),
+        ]
+        field.forces = (0, 0, 0)
+        field.momentums = (0, 0, 0)
+        stress, strain = field.eval(points)
+        assert np.all(np.abs(stress) < 1e-9)  # no charge, all zero
+        assert np.all(np.abs(strain) < 1e-9)  # no charge, all zero
+
+        field.forces = (0, 0, 0)
+        field.momentums = (4 / 3, 0, 0)
+        stress, strain = field.eval(points)
+        assert np.all(np.abs(stress[:, :2]) < 1e-9)  # No shear
+        good_normal_stress = [0]
+        good_normal_stress += [0, 1, 1, 1, 0, -1, -1, -1]
+        abs_diff = np.abs(stress[:, 2] - good_normal_stress)
+        assert np.all(abs_diff < 1e-9)
+
+        field.forces = (0, 0, 0)
+        field.momentums = (0, 4 / 3, 0)
+        stress, strain = field.eval(points)
+        assert np.all(np.abs(stress[:, :2]) < 1e-9)  # No shear
+        good_normal_stress = [0]
+        good_normal_stress += [-1, -1, 0, 1, 1, 1, 0, -1]
+        abs_diff = np.abs(stress[:, 2] - good_normal_stress)
+        assert np.all(abs_diff < 1e-9)
 
     @pytest.mark.order(4)
     @pytest.mark.dependency(
@@ -66,11 +97,55 @@ def test_centered_square(self):
         geometry = Primitive.square(ext_side)
         geometry -= Primitive.square(int_side)
         material = Isotropic(young_modulus=210e3, poissons_ratio=0.3)
-        section = Section(geometry, material)
+        section = Section.from_shapes(geometry, material)
         field = section.charged_field()
-        points = [(0, 0)]
-        values = field.eval(points)
-        assert np.all(np.abs(values) < 1e-9)
+
+        points = [(0, 0)]  # origin
+        points += [
+            (0.5, 0),
+            (0.5, 0.5),
+            (0, 0.5),
+            (-0.5, 0.5),
+            (-0.5, 0),
+            (-0.5, -0.5),
+            (0, -0.5),
+            (0.5, -0.5),
+        ]
+        points += [
+            (1, 0),
+            (1, 1),
+            (0, 1),
+            (-1, 1),
+            (-1, 0),
+            (-1, -1),
+            (0, -1),
+            (1, -1),
+        ]
+        field.forces = (0, 0, 0)
+        field.momentums = (0, 0, 0)
+        stress, strain = field.eval(points)
+        assert np.all(np.abs(stress) < 1e-9)  # no charge, all zero
+        assert np.all(np.abs(strain) < 1e-9)  # no charge, all zero
+
+        field.forces = (0, 0, 0)
+        field.momentums = (1.25, 0, 0)
+        stress, strain = field.eval(points)
+        assert np.all(np.abs(stress[:, :2]) < 1e-9)  # No shear
+        good_normal_stress = [0]  # origin
+        good_normal_stress += [0, 0.5, 0.5, 0.5, 0, -0.5, -0.5, -0.5]
+        good_normal_stress += [0, 1, 1, 1, 0, -1, -1, -1]
+        abs_diff = np.abs(stress[:, 2] - good_normal_stress)
+        assert np.all(abs_diff < 1e-9)
+
+        field.forces = (0, 0, 0)
+        field.momentums = (0, 1.25, 0)
+        stress, strain = field.eval(points)
+        assert np.all(np.abs(stress[:, :2]) < 1e-9)  # No shear
+        good_normal_stress = [0]
+        good_normal_stress += [-0.5, -0.5, 0, 0.5, 0.5, 0.5, 0, -0.5]
+        good_normal_stress += [-1, -1, 0, 1, 1, 1, 0, -1]
+        abs_diff = np.abs(stress[:, 2] - good_normal_stress)
+        assert np.all(abs_diff < 1e-9)
 
     @pytest.mark.order(4)
     @pytest.mark.dependency(
diff --git a/tests/test_curve.py b/tests/test_curve.py
new file mode 100644
index 0000000..957e571
--- /dev/null
+++ b/tests/test_curve.py
@@ -0,0 +1,47 @@
+import numpy as np
+import pytest
+
+from compmec import nurbs
+from compmec.section.curve import NurbsCurve
+
+
+@pytest.mark.order(1)
+@pytest.mark.dependency()
+def test_begin():
+    pass
+
+
+@pytest.mark.order(1)
+@pytest.mark.timeout(10)
+@pytest.mark.dependency(depends=["test_begin"])
+def test_winding_square():
+    knotvector = nurbs.GeneratorKnotVector.uniform(1, 5)
+    vertices = [[1, 1], [-1, 1], [-1, -1], [1, -1], [1, 1]]
+    vertices = np.array(vertices, dtype="float64")
+    curve = nurbs.Curve(knotvector, vertices)
+    curve = NurbsCurve(curve)
+
+    assert curve.winding((0, 0)) == 1
+    assert curve.winding((1, 0)) == 0.5
+    assert curve.winding((1, 1)) == 0.25
+    assert curve.winding((0, 1)) == 0.5
+    assert curve.winding((-1, 1)) == 0.25
+    assert curve.winding((-1, 0)) == 0.5
+    assert curve.winding((-1, -1)) == 0.25
+    assert curve.winding((0, -1)) == 0.5
+    assert curve.winding((1, -1)) == 0.25
+    assert curve.winding((1, 0)) == 0.5
+
+
+@pytest.mark.order(1)
+@pytest.mark.dependency(
+    depends=[
+        "test_winding_square",
+    ]
+)
+def test_end():
+    pass
+
+
+if __name__ == "__main__":
+    test_winding_square()
diff --git a/tests/test_geometry.py b/tests/test_geometry.py
new file mode 100644
index 0000000..c120751
--- /dev/null
+++ b/tests/test_geometry.py
@@ -0,0 +1,100 @@
+"""
+File to test src/compmec/section/geometry.py module
+"""
+
+import numpy as np
+import pytest
+
+from compmec import nurbs
+from compmec.section.curve import NurbsCurve
+from compmec.section.geometry import Geometry
+
+
+@pytest.mark.order(1)
+@pytest.mark.dependency()
+def test_begin():
+    pass
+
+
+@pytest.mark.order(1)
+@pytest.mark.timeout(10)
+@pytest.mark.dependency(depends=["test_begin"])
+def test_full_square():
+    knotvector = nurbs.GeneratorKnotVector.uniform(1, 5)
+    vertices = [[1, 1], [-1, 1], [-1, -1], [1, -1], [1, 1]]
+    vertices = np.array(vertices, dtype="float64")
+    curve = nurbs.Curve(knotvector, vertices)
+    curve = NurbsCurve(curve)
+
+    geometry = Geometry([curve.label])
+
+    assert geometry.winding((0, 0)) == 1
+    assert geometry.winding((1, 0)) == 0.5
+    assert geometry.winding((1, 1)) == 0.25
+    assert geometry.winding((0, 1)) == 0.5
+    assert geometry.winding((-1, 1)) == 0.25
+    assert geometry.winding((-1, 0)) == 0.5
+    assert geometry.winding((-1, -1)) == 0.25
+    assert geometry.winding((0, -1)) == 0.5
+    assert geometry.winding((1, -1)) == 0.25
+    assert geometry.winding((1, 0)) == 0.5
+
+
+@pytest.mark.order(1)
+@pytest.mark.timeout(10)
+@pytest.mark.dependency(depends=["test_begin"])
+def test_hollow_square():
+    knotvector = nurbs.GeneratorKnotVector.uniform(1, 5)
+    vertices = [[3, 3], [-3, 3], [-3, -3], [3, -3], [3, 3]]
+    vertices = np.array(vertices, dtype="float64")
+    curve_ext = nurbs.Curve(knotvector, vertices)
+    curve_ext = NurbsCurve(curve_ext)
+
+    vertices = [[1, 1], [-1, 1], [-1, -1], [1, -1], [1, 1]]
+    vertices = np.array(vertices, dtype="float64")
+    curve_int = nurbs.Curve(knotvector, vertices)
+    curve_int = NurbsCurve(curve_int)
+
+    geometry = Geometry([curve_ext.label, -curve_int.label])
+
+    assert geometry.winding((0, 0)) == 0
+    assert geometry.winding((1, 0)) == 0.5
+    assert geometry.winding((1, 1)) == 0.25
+    assert geometry.winding((0, 1)) == 0.5
+    assert geometry.winding((-1, 1)) == 0.25
+    assert geometry.winding((-1, 0)) == 0.5
+    assert geometry.winding((-1, -1)) == 0.25
+    assert geometry.winding((0, -1)) == 0.5
+    assert geometry.winding((1, -1)) == 0.25
+    assert geometry.winding((1, 0)) == 0.5
+
+    assert geometry.winding((2, 0)) == 1
+    assert geometry.winding((2, 2)) == 1
+    assert geometry.winding((0, 2)) == 1
+    assert geometry.winding((-2, 2)) == 1
+    assert geometry.winding((-2, 0)) == 1
+    assert geometry.winding((-2, -2)) == 1
+    assert geometry.winding((0, -2)) == 1
+    assert geometry.winding((2, -2)) == 1
+    assert geometry.winding((2, 0)) == 1
+
+    assert geometry.winding((3, 0)) == 0.5
+    assert geometry.winding((3, 3)) == 0.25
+    assert geometry.winding((0, 3)) == 0.5
+    assert geometry.winding((-3, 3)) == 0.25
+    assert geometry.winding((-3, 0)) == 0.5
+    assert geometry.winding((-3, -3)) == 0.25
+    assert geometry.winding((0, -3)) == 0.5
+    assert geometry.winding((3, -3)) == 0.25
+    assert geometry.winding((3, 0)) == 0.5
+
+
+@pytest.mark.order(1)
+@pytest.mark.dependency(
+    depends=[
+        "test_full_square",
+        "test_hollow_square",
+    ]
+)
+def test_end():
+    pass
diff --git a/tests/test_geomprop.py b/tests/test_geomprop.py
index 28053f7..efc48e5 100644
--- a/tests/test_geomprop.py
+++ b/tests/test_geomprop.py
@@ -17,7 +17,11 @@
 
 @pytest.mark.order(3)
 @pytest.mark.dependency(
-    depends=["tests/test_material.py::test_end"],
+    depends=[
+        "tests/test_material.py::test_end",
+        "tests/test_curve.py::test_end",
+        "tests/test_geometry.py::test_end",
+    ],
     scope="session",
 )
 def test_begin():
@@ -111,9 +115,11 @@ def test_shifted_rectangle(self):
         assert area == width * height
         assert Qx == area * center[1]
         assert Qy == area * center[0]
-        assert Ixx == width * height**3 / 12 + area * center[1] * center[1]
+        good_Ixx = width * height**3 / 12 + area * center[1] * center[1]
+        good_Iyy = height * width**3 / 12 + area * center[0] * center[0]
+        assert abs(Ixx - good_Ixx) < 1e-6
         assert Ixy == area * center[0] * center[1]
-        assert Iyy == height * width**3 / 12 + area * center[0] * center[0]
+        assert abs(Iyy - good_Iyy) < 1e-6
 
     @pytest.mark.order(3)
     @pytest.mark.dependency(
@@ -142,10 +148,10 @@ def test_begin(self):
     def test_read_square(self):
         json_filepath = "tests/json/steel_square.json"
         with JsonIO(json_filepath) as reader:
-            reader.read_nodes()
-            reader.read_curves()
-            reader.read_materials()
-            reader.read_sections()
+            reader.load_nodes()
+            reader.load_curves()
+            reader.load_materials()
+            reader.load_sections()
         square = Section.instances["square"]
         area = square.area()
         Qx, Qy = square.first_moment()
@@ -153,31 +159,30 @@ def test_read_square(self):
         assert area == 4
         assert Qx == 0
         assert Qy == 0
-        assert Ixx == 4 / 3
+        assert abs(Ixx - 4 / 3) < 1e-9
         assert Ixy == 0
-        assert Iyy == 4 / 3
+        assert abs(Iyy - 4 / 3) < 1e-9
 
     @pytest.mark.order(3)
-    @pytest.mark.timeout(1)
-    @pytest.mark.skip(reason="Not yet implemented")
+    @pytest.mark.timeout(20)
     @pytest.mark.dependency(depends=["TestToFromJson::test_begin"])
     def test_read_circle(self):
         json_filepath = "tests/json/steel_circle.json"
         with JsonIO(json_filepath) as reader:
-            reader.read_nodes()
-            reader.read_curves()
-            reader.read_materials()
-            reader.read_sections()
+            reader.load_nodes()
+            reader.load_curves()
+            reader.load_materials()
+            reader.load_sections()
         square = Section.instances["square"]
         area = square.area()
         Qx, Qy = square.first_moment()
         Ixx, Ixy, Iyy = square.second_moment()
-        assert area == math.pi
-        assert Qx == 0
-        assert Qy == 0
-        assert Ixx == math.pi / 4
-        assert Ixy == 0
-        assert Iyy == math.pi / 4
+        assert abs(area - math.pi) < 2e-6
+        assert abs(Qx) < 1e-6
+        assert abs(Qy) < 1e-6
+        assert abs(Ixx - math.pi / 4) < 1e-6
+        assert abs(Ixy) < 1e-6
+        assert abs(Iyy - math.pi / 4) < 1e-6
 
     @pytest.mark.order(3)
     @pytest.mark.dependency(
diff --git a/tests/test_integral.py b/tests/test_integral.py
index 847895b..593fec6 100644
--- a/tests/test_integral.py
+++ b/tests/test_integral.py
@@ -1,7 +1,7 @@
 import numpy as np
 import pytest
 
-from compmec.section.integral import Integration
+from compmec.section.integral import Integration, Polynomial, comb
 
 
 @pytest.mark.order(1)
@@ -10,6 +10,47 @@ def test_begin():
     pass
 
 
+@pytest.mark.order(1)
+@pytest.mark.timeout(10)
+@pytest.mark.dependency(depends=["test_begin"])
+def test_comb():
+    """
+    Tests the binomial function
+    """
+    assert comb(1, 0) == 1
+    assert comb(1, 1) == 1
+
+    assert comb(2, 0) == 1
+    assert comb(2, 1) == 2
+    assert comb(2, 2) == 1
+
+    assert comb(3, 0) == 1
+    assert comb(3, 1) == 3
+    assert comb(3, 2) == 3
+    assert comb(3, 3) == 1
+
+    assert comb(4, 0) == 1
+    assert comb(4, 1) == 4
+    assert comb(4, 2) == 6
+    assert comb(4, 3) == 4
+    assert comb(4, 4) == 1
+
+    assert comb(5, 0) == 1
+    assert comb(5, 1) == 5
+    assert comb(5, 2) == 10
+    assert comb(5, 3) == 10
+    assert comb(5, 4) == 5
+    assert comb(5, 5) == 1
+
+    assert comb(6, 0) == 1
+    assert comb(6, 1) == 6
+    assert comb(6, 2) == 15
+    assert comb(6, 3) == 20
+    assert comb(6, 4) == 15
+    assert comb(6, 5) == 6
+    assert comb(6, 6) == 1
+
+
 @pytest.mark.order(1)
 @pytest.mark.timeout(10)
 @pytest.mark.dependency(depends=["test_begin"])
@@ -24,7 +65,6 @@ def test_closed_newton():
     """
     for nptsinteg in range(2, 8):
         nodes, weights = Integration.closed(nptsinteg)
-        print("npts = ", nptsinteg)
         for degree in range(nptsinteg):  # Integrate each basis
             good_integral = 1 / (1 + degree)
             fvalues = nodes**degree
@@ -48,6 +88,43 @@ def test_closed_newton():
             assert diff < 1e-15
 
 
+@pytest.mark.order(1)
+@pytest.mark.timeout(10)
+@pytest.mark.dependency(depends=["test_begin"])
+def test_opened_newton():
+    """
+    Tests if the given nodes and weights integrates exactly
+    a polynomial of degree p, with p+1 evaluation points
+
+    We start with P(x) = x^{p}
+    And after we evaluate P(x) = sum_i c_i * x^i
+    for random c_i
+    """
+    for nptsinteg in range(1, 4):
+        nodes, weights = Integration.opened(nptsinteg)
+        for degree in range(nptsinteg):  # Integrate each basis
+            good_integral = 1 / (1 + degree)
+            fvalues = nodes**degree
+            test_integral = np.inner(fvalues, weights)
+            diff = abs(test_integral - good_integral)
+            assert abs(diff) < 1e-5
+
+    for nptsinteg in range(1, 4):
+        nodes, weights = Integration.opened(nptsinteg)
+        fvalues = np.zeros(len(nodes))
+        for degree in range(nptsinteg):
+            coefs = np.random.uniform(-1, 1, degree + 1)
+            good_integral = sum(ci / (1 + i) for i, ci in enumerate(coefs))
+            fvalues.fill(0)
+            for i, node in enumerate(nodes):
+                for ck in coefs[::-1]:  # Horner's method
+                    fvalues[i] *= node
+                    fvalues[i] += ck
+            test_integral = np.inner(fvalues, weights)
+            diff = abs(test_integral - good_integral)
+            assert diff < 1e-15
+
+
 @pytest.mark.order(1)
 @pytest.mark.timeout(10)
 @pytest.mark.dependency(depends=["test_begin"])
@@ -161,13 +238,57 @@ def test_singular_logarithm():
             assert diff < 1e-15
 
 
+@pytest.mark.order(1)
+@pytest.mark.timeout(10)
+@pytest.mark.dependency(depends=["test_begin"])
+def test_integral_polygon():
+    """
+    Tests the polynomial integrals of polygons
+    """
+    vertices = ((1, 1), (-1, 1), (-1, -1), (1, -1))
+    geomprops = Polynomial.polygon(vertices)
+    assert geomprops[0] == 4  # area
+    assert geomprops[1] == 0  # Qx
+    assert geomprops[2] == 0  # Qy
+    assert geomprops[3] == 4 / 3  # Ixx
+    assert geomprops[4] == 0  # Ixy
+    assert geomprops[5] == 4 / 3  # Iyy
+    assert geomprops[6] == 0  # Ixxx
+    assert geomprops[7] == 0  # Ixxy
+    assert geomprops[8] == 0  # Ixyy
+    assert geomprops[9] == 0  # Iyyy
+
+
+@pytest.mark.order(1)
+@pytest.mark.timeout(10)
+@pytest.mark.dependency(depends=["test_begin"])
+def test_fail():
+    """
+    Tests the fail cases
+    """
+    with pytest.raises(ValueError):
+        Integration.opened(0)
+    with pytest.raises(ValueError):
+        Integration.closed(1)
+    with pytest.raises(ValueError):
+        Integration.chebyshev(0)
+    with pytest.raises(ValueError):
+        Integration.gauss(0)
+    with pytest.raises(ValueError):
+        Integration.log(0)
+
+
 @pytest.mark.order(1)
 @pytest.mark.dependency(
     depends=[
+        "test_comb",
         "test_closed_newton",
+        "test_opened_newton",
         "test_chebyshev",
         "test_gauss",
         "test_singular_logarithm",
+        "test_integral_polygon",
+        "test_fail",
     ]
 )
 def test_end():