First implementation of 2D Cartesian space.

Files and types with 'Cyl' have been changed to '2DCyl' to better
differentiate between the two types of 2D geometry.

Solvers for charged and neutral particles in 2D Cartesian space.

Added solveds for 1D neutral particles (this branch is not the place to
    do it, but it was a minor change).

User Manual updated with the new accepted options.

src/modules/mesh/2DCart/makefile

@@ -0,0 +1,11 @@
+all : moduleMesh2DCart.o moduleMesh2DCartBoundary.o moduleMesh2DCartRead.o
+
+moduleMesh2DCart.o: moduleMesh2DCart.f90
+	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
+
+moduleMesh2DCartBoundary.o: moduleMesh2DCart.o moduleMesh2DCartBoundary.f90
+	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
+
+moduleMesh2DCartRead.o: moduleMesh2DCart.o moduleMesh2DCartBoundary.o moduleMesh2DCartRead.f90
+	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
+

src/modules/mesh/2DCart/moduleMesh2DCartBoundary.f90

@@ -0,0 +1,154 @@
+!moduleMesh2DCartBoundary: Boundary functions for cylindrical coordinates
+SUBMODULE (moduleMesh2DCart) moduleMesh2DCartBoundary
+  USE moduleMesh2DCart
+
+  CONTAINS
+    SUBROUTINE reflection(edge, part)
+      USE moduleSpecies
+      IMPLICIT NONE
+
+      CLASS(meshEdge), INTENT(inout):: edge
+      CLASS(particle), INTENT(inout):: part
+      REAL(8):: edgeNorm, cosT, sinT, rp(1:2), rpp(1:2), vpp(1:2)
+
+      !TODO: Try to do this without select
+      SELECT TYPE(edge)
+      TYPE IS(meshEdge2DCart)
+        edgeNorm = DSQRT((edge%y(2)-edge%y(1))**2 + (edge%x(2)-edge%x(1))**2)
+        cosT = (edge%x(2)-edge%x(1))/edgeNorm
+        sinT = DSQRT(1-cosT**2)
+
+        rp(1) = part%r(1) - edge%x(1);
+        rp(2) = part%r(2) - edge%y(1);
+
+        rpp(1) = cosT*rp(1) - sinT*rp(2)
+        rpp(2) = sinT*rp(1) + cosT*rp(2)
+        rpp(2) = -rpp(2)
+
+        vpp(1) = cosT*part%v(1) - sinT*part%v(2)
+        vpp(2) = sinT*part%v(1) + cosT*part%v(2)
+        vpp(2) = -vpp(2)
+
+        part%r(1) =  cosT*rpp(1) + sinT*rpp(2) + edge%x(1);
+        part%r(2) = -sinT*rpp(1) + cosT*rpp(2) + edge%y(1);
+        part%v(1) =  cosT*vpp(1) + sinT*vpp(2)
+        part%v(2) = -sinT*vpp(1) + cosT*vpp(2)
+
+      END SELECT
+
+      part%n_in = .TRUE.
+
+    END SUBROUTINE reflection
+
+    !Absoption in a surface
+    SUBROUTINE absorption(edge, part)
+      USE moduleSpecies
+      IMPLICIT NONE
+
+      CLASS(meshEdge), INTENT(inout):: edge
+      CLASS(particle), INTENT(inout):: part
+      REAL(8):: rEdge(1:2) !Position of particle projected to the edge
+      REAL(8):: a, b, c
+      REAL(8):: a2b2
+      REAL(8):: d !Distance from particle to edge
+
+      SELECT TYPE(edge)
+      TYPE IS(meshEdge2DCart)
+        a = (edge%x(1) - edge%x(2))
+        b = (edge%y(1) - edge%y(2))
+        c =  edge%x(1)*edge%y(2) - edge%x(2)*edge%y(1)
+
+        a2b2 = a**2 + b**2
+
+        rEdge(1) = (b*( b*part%r(1) - a*part%r(2)) - a*c)/a2b2
+        rEdge(2) = (a*(-b*part%r(1) + a*part%r(2)) - b*c)/a2b2
+
+        d = NORM2(rEdge - part%r(1:2))
+        !Reduce weight of particle by the distance to the edge and move it to the edge
+        IF (d > 0.D0) THEN
+          part%weight = part%weight / d
+          part%r(1:2) = rEdge
+
+        END IF
+
+        !Scatter particle in associated volume
+        IF (ASSOCIATED(edge%e1)) THEN
+          CALL edge%e1%scatter(part)
+
+        ELSE
+          CALL edge%e2%scatter(part)
+
+        END IF
+
+      END SELECT
+
+      !Remove particle from the domain
+      part%n_in = .FALSE.
+
+    END SUBROUTINE absorption
+
+    !Transparent boundary condition
+    SUBROUTINE transparent(edge, part)
+      USE moduleSpecies
+      IMPLICIT NONE
+
+      CLASS(meshEdge), INTENT(inout):: edge
+      CLASS(particle), INTENT(inout):: part
+
+      !Removes particle from domain
+      part%n_in = .FALSE.
+
+    END SUBROUTINE transparent
+
+    !Wall with temperature
+    SUBROUTINE wallTemperature(edge, part)
+      USE moduleSpecies
+      USE moduleBoundary
+      USE moduleRandom
+      IMPLICIT NONE
+
+      CLASS(meshEdge), INTENT(inout):: edge
+      CLASS(particle), INTENT(inout):: part
+      REAL(8):: edgeNorm, cosT, sinT, rp(1:2), rpp(1:2), vpp(1:2)
+      INTEGER:: i
+
+      !Modifies particle velocity according to wall temperature
+      SELECT TYPE(bound => edge%boundary%bTypes(part%sp)%obj)
+      TYPE IS(boundaryWallTemperature)
+        DO i = 1, 3
+          part%v(i) = part%v(i) + bound%vTh*randomMaxwellian()
+
+        END DO
+          
+      END SELECT
+
+      !Reflects particle in the edge
+      SELECT TYPE(edge)
+      TYPE IS(meshEdge2DCart)
+        edgeNorm = DSQRT((edge%y(2)-edge%y(1))**2 + (edge%x(2)-edge%x(1))**2)
+        cosT = (edge%x(2)-edge%x(1))/edgeNorm
+        sinT = DSQRT(1-cosT**2)
+
+        rp(1) = part%r(1) - edge%x(1);
+        rp(2) = part%r(2) - edge%y(1);
+
+        rpp(1) = cosT*rp(1) - sinT*rp(2)
+        rpp(2) = sinT*rp(1) + cosT*rp(2)
+        rpp(2) = -rpp(2)
+
+        vpp(1) = cosT*part%v(1) - sinT*part%v(2)
+        vpp(2) = sinT*part%v(1) + cosT*part%v(2)
+        vpp(2) = -vpp(2)
+
+        part%r(1) =  cosT*rpp(1) + sinT*rpp(2) + edge%x(1);
+        part%r(2) = -sinT*rpp(1) + cosT*rpp(2) + edge%y(1);
+        part%v(1) =  cosT*vpp(1) + sinT*vpp(2)
+        part%v(2) = -sinT*vpp(1) + cosT*vpp(2)
+
+      END SELECT
+
+      part%n_in = .TRUE.
+
+    END SUBROUTINE wallTemperature
+
+END SUBMODULE moduleMesh2DCartBoundary

src/modules/mesh/2DCart/moduleMesh2DCartRead.f90

@@ -0,0 +1,599 @@
+MODULE moduleMesh2DCartRead
+  USE moduleMesh
+  USE moduleMesh2DCart
+
+  TYPE, EXTENDS(meshGeneric):: mesh2DCartGeneric
+    CONTAINS
+      PROCEDURE, PASS:: init     => init2DCartMesh
+      PROCEDURE, PASS:: readMesh => readMesh2DCartGmsh
+
+  END TYPE
+
+  INTERFACE connected
+    MODULE PROCEDURE connectedVolVol, connectedVolEdge
+
+  END INTERFACE connected
+
+  CONTAINS
+    !Init mesh
+    SUBROUTINE init2DCartMesh(self, meshFormat)
+      USE moduleMesh
+      USE moduleErrors
+      IMPLICIT NONE
+
+      CLASS(mesh2DCartGeneric), INTENT(out):: self
+      CHARACTER(:), ALLOCATABLE, INTENT(in):: meshFormat
+
+      SELECT CASE(meshFormat)
+      CASE ("gmsh")
+        self%printOutput => printOutputGmsh
+        self%printColl   => printCollGmsh
+        self%printEM     => printEMGmsh
+
+      CASE DEFAULT
+        CALL criticalError("Mesh type " // meshFormat // " not supported.", "init2DCartMesh")
+        
+      END SELECT
+
+    END SUBROUTINE init2DCartMesh
+
+    !Read mesh from gmsh file
+    SUBROUTINE readMesh2DCartGmsh(self, filename)
+      USE moduleBoundary
+      IMPLICIT NONE
+
+      CLASS(mesh2DCartGeneric), INTENT(inout):: self
+      CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
+      REAL(8):: x, y
+      INTEGER:: p(1:4)
+      INTEGER:: e=0, et=0, n=0, eTemp=0, elemType=0, bt = 0
+      INTEGER:: totalNumElem
+      INTEGER:: boundaryType
+
+      !Read msh
+      OPEN(10, FILE=TRIM(filename))
+      !Skip header
+      READ(10, *)
+      READ(10, *)
+      READ(10, *)
+      READ(10, *)
+      !Read number of nodes
+      READ(10, *) self%numNodes
+      !Allocate required matrices and vectors
+      ALLOCATE(self%nodes(1:self%numNodes))
+      ALLOCATE(self%K(1:self%numNodes,1:self%numNodes))
+      ALLOCATE(self%IPIV(1:self%numNodes,1:self%numNodes))
+      self%K = 0.D0
+      self%IPIV = 0
+      !Read nodes cartesian coordinates (x=x, y=y, z=null)
+      DO e=1, self%numNodes
+        READ(10, *) n, x, y
+        ALLOCATE(meshNode2DCart:: self%nodes(n)%obj)
+        CALL self%nodes(n)%obj%init(n, (/x, y, 0.D0 /))
+
+      END DO
+      !Skips comments
+      READ(10, *)
+      READ(10, *)
+      !Reads Totalnumber of elements
+      READ(10, *) TotalnumElem
+      !counts edges and volume elements
+      self%numEdges = 0
+      DO e=1, TotalnumElem
+        READ(10,*) eTemp, elemType
+        IF (elemType==1) THEN
+          self%numEdges=e
+        END IF
+      END DO
+      !Substract the number of edges to the total number of elements
+      !to obtain the number of volume elements
+      self%numVols = TotalnumElem - self%numEdges
+      !Allocates arrays
+      ALLOCATE(self%edges(1:self%numEdges))
+      ALLOCATE(self%vols(1:self%numVols))
+
+      !Go back to the beggining to read elements
+      DO e=1, totalNumElem
+        BACKSPACE(10)
+      END DO
+
+      !Reads edges
+      DO e=1, self%numEdges
+        READ(10,*) n, elemType, eTemp, boundaryType, eTemp, p(1:2)
+        !Associate boundary condition procedure.
+        bt = getBoundaryId(boundaryType)
+
+        ALLOCATE(meshEdge2DCart:: self%edges(e)%obj)
+
+        CALL self%edges(e)%obj%init(n, p(1:2), bt, boundaryType)
+
+      END DO
+
+      !Read and initialize volumes
+      DO e=1, self%numVols
+        READ(10,*) n, elemType
+        BACKSPACE(10)
+
+        SELECT CASE(elemType)
+        CASE (2)
+          !Triangular element
+          READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:3)
+          ALLOCATE(meshVol2DCartTria:: self%vols(e)%obj)
+          CALL self%vols(e)%obj%init(n - self%numEdges, p(1:3))
+
+        CASE (3)
+          !Quadrilateral element
+          READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:4)
+          ALLOCATE(meshVol2DCartQuad:: self%vols(e)%obj)
+          CALL self%vols(e)%obj%init(n - self%numEdges, p(1:4))
+
+        END SELECT
+
+      END DO
+
+      CLOSE(10)
+
+      !Build connectivity between elements
+      DO e = 1, self%numVols
+        !Connectivity between volumes
+        DO et = 1, self%numVols
+          IF (e /= et) THEN
+            CALL connected(self%vols(e)%obj, self%vols(et)%obj)
+
+          END IF
+        END DO
+
+        !Connectivity between vols and edges
+        DO et = 1, self%numEdges
+          CALL connected(self%vols(e)%obj, self%edges(et)%obj)
+
+        END DO
+
+        !Constructs the global K matrix
+        CALL constructGlobalK(self%K, self%vols(e)%obj)
+
+      END DO
+
+    END SUBROUTINE readMesh2DCartGmsh
+
+    !Selects type of elements to build connection
+    SUBROUTINE connectedVolVol(elemA, elemB)
+      IMPLICIT NONE
+
+      CLASS(meshVol), INTENT(inout):: elemA
+      CLASS(meshVol), INTENT(inout):: elemB
+
+      SELECT TYPE(elemA)
+      TYPE IS(meshVol2DCartQuad)
+        !Element A is a quadrilateral
+        SELECT TYPE(elemB)
+        TYPE IS(meshVol2DCartQuad)
+          !Element B is a quadrilateral
+          CALL connectedQuadQuad(elemA, elemB)
+
+        TYPE IS(meshVol2DCartTria)
+          !Element B is a triangle
+          CALL connectedQuadTria(elemA, elemB)
+
+        END SELECT
+
+      TYPE IS(meshVol2DCartTria)
+        !Element A is a Triangle
+        SELECT TYPE(elemB)
+        TYPE IS(meshVol2DCartQuad)
+          !Element B is a quadrilateral
+          CALL connectedQuadTria(elemB, elemA)
+
+        TYPE IS(meshVol2DCartTria)
+          !Element B is a triangle
+          CALL connectedTriaTria(elemA, elemB)
+
+        END SELECT
+
+      END SELECT
+
+    END SUBROUTINE connectedVolVol
+
+    SUBROUTINE connectedVolEdge(elemA, elemB)
+      IMPLICIT NONE
+
+      CLASS(meshVol),  INTENT(inout):: elemA
+      CLASS(meshEdge), INTENT(inout):: elemB
+
+      SELECT TYPE(elemB)
+      CLASS IS(meshEdge2DCart)
+        SELECT TYPE(elemA)
+        TYPE IS(meshVol2DCartQuad)
+          !Element A is a quadrilateral
+          CALL connectedQuadEdge(elemA, elemB)
+
+        TYPE IS(meshVol2DCartTria)
+          !Element A is a triangle
+          CALL connectedTriaEdge(elemA, elemB)
+
+        END SELECT
+
+      END SELECT
+
+    END SUBROUTINE connectedVolEdge
+
+    SUBROUTINE connectedQuadQuad(elemA, elemB)
+      IMPLICIT NONE
+
+      CLASS(meshVol2DCartQuad), INTENT(inout), TARGET:: elemA
+      CLASS(meshVol2DCartQuad), INTENT(inout), TARGET:: elemB
+
+      !Check direction 1
+      IF (.NOT. ASSOCIATED(elemA%e1) .AND. &
+          elemA%n1%n == elemB%n4%n   .AND. &
+          elemA%n2%n == elemB%n3%n) THEN
+        elemA%e1 => elemB
+        elemB%e3 => elemA
+
+      END IF
+
+      !Check direction 2
+      IF (.NOT. ASSOCIATED(elemA%e2) .AND. &
+          elemA%n2%n == elemB%n1%n   .AND. &
+          elemA%n3%n == elemB%n4%n) THEN
+        elemA%e2 => elemB
+        elemB%e4 => elemA
+
+      END IF
+
+      !Check direction 3
+      IF (.NOT. ASSOCIATED(elemA%e3) .AND. &
+          elemA%n3%n == elemB%n2%n   .AND. &
+          elemA%n4%n == elemB%n1%n) THEN
+        elemA%e3 => elemB
+        elemB%e1 => elemA
+
+      END IF
+
+      !Check direction 4
+      IF (.NOT. ASSOCIATED(elemA%e4) .AND. &
+          elemA%n4%n == elemB%n3%n   .AND. &
+          elemA%n1%n == elemB%n2%n) THEN
+        elemA%e4 => elemB
+        elemB%e2 => elemA
+
+      END IF
+
+    END SUBROUTINE connectedQuadQuad
+
+    SUBROUTINE connectedQuadTria(elemA, elemB)
+      IMPLICIT NONE
+
+      CLASS(meshVol2DCartQuad), INTENT(inout), TARGET:: elemA
+      CLASS(meshVol2DCartTria), INTENT(inout), TARGET:: elemB
+
+      !Check direction 1
+      IF (.NOT. ASSOCIATED(elemA%e1)) THEN
+        IF     (elemA%n1%n == elemB%n1%n .AND. &
+                elemA%n2%n == elemB%n3%n) THEN
+          elemA%e1 => elemB
+          elemB%e3 => elemA
+
+        ELSEIF (elemA%n1%n == elemB%n3%n .AND. &
+                elemA%n2%n == elemB%n2%n) THEN
+          elemA%e1 => elemB
+          elemB%e2 => elemA
+
+        ELSEIF (elemA%n1%n == elemB%n2%n .AND. &
+                elemA%n2%n == elemB%n1%n) THEN
+          elemA%e1 => elemB
+          elemB%e1 => elemA
+
+        END IF
+
+      END IF
+
+      !Check direction 2
+      IF (.NOT. ASSOCIATED(elemA%e2)) THEN
+        IF     (elemA%n2%n == elemB%n1%n .AND. &
+                elemA%n3%n == elemB%n3%n) THEN
+          elemA%e2 => elemB
+          elemB%e3 => elemA
+
+        ELSEIF (elemA%n2%n == elemB%n3%n .AND. &
+                elemA%n3%n == elemB%n2%n) THEN
+          elemA%e2 => elemB
+          elemB%e2 => elemA
+
+        ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
+                elemA%n3%n == elemB%n1%n) THEN
+          elemA%e2 => elemB
+          elemB%e1 => elemA
+
+        END IF
+
+      END IF
+
+      !Check direction 3
+      IF (.NOT. ASSOCIATED(elemA%e3)) THEN
+        IF     (elemA%n3%n == elemB%n1%n .AND. &
+                elemA%n4%n == elemB%n3%n) THEN
+          elemA%e3 => elemB
+          elemB%e3 => elemA
+
+        ELSEIF (elemA%n3%n == elemB%n3%n .AND. &
+                elemA%n4%n == elemB%n2%n) THEN
+          elemA%e3 => elemB
+          elemB%e2 => elemA
+
+        ELSEIF (elemA%n3%n == elemB%n2%n .AND. &
+                elemA%n4%n == elemB%n1%n) THEN
+          elemA%e3 => elemB
+          elemB%e1 => elemA
+
+        END IF
+
+      END IF
+
+      !Check direction 4
+      IF (.NOT. ASSOCIATED(elemA%e4)) THEN
+        IF     (elemA%n4%n == elemB%n1%n .AND. &
+                elemA%n1%n == elemB%n3%n) THEN
+          elemA%e4 => elemB
+          elemB%e3 => elemA
+
+        ELSEIF (elemA%n4%n == elemB%n3%n .AND. &
+                elemA%n1%n == elemB%n2%n) THEN
+          elemA%e4 => elemB
+          elemB%e2 => elemA
+
+        ELSEIF (elemA%n4%n == elemB%n2%n .AND. &
+                elemA%n1%n == elemB%n1%n) THEN
+          elemA%e4 => elemB
+          elemB%e1 => elemA
+
+        END IF
+
+      END IF
+
+    END SUBROUTINE connectedQuadTria
+
+    SUBROUTINE connectedTriaTria(elemA, elemB)
+      IMPLICIT NONE
+
+      CLASS(meshVol2DCartTria), INTENT(inout), TARGET:: elemA
+      CLASS(meshVol2DCartTria), INTENT(inout), TARGET:: elemB
+
+      !Check direction 1
+      IF (.NOT. ASSOCIATED(elemA%e1)) THEN
+        IF     (elemA%n1%n == elemB%n1%n .AND. &
+                elemA%n2%n == elemB%n3%n) THEN
+          elemA%e1 => elemB
+          elemB%e3 => elemA
+
+        ELSEIF (elemA%n1%n == elemB%n2%n .AND. &
+                elemA%n2%n == elemB%n1%n) THEN
+          elemA%e1 => elemB
+          elemB%e1 => elemA
+
+        ELSEIF (elemA%n1%n == elemB%n3%n .AND. &
+                elemA%n2%n == elemB%n2%n) THEN
+          elemA%e1 => elemB
+          elemB%e2 => elemA
+
+        END IF
+
+      END IF
+
+      !Check direction 2
+      IF (.NOT. ASSOCIATED(elemA%e2)) THEN
+        IF     (elemA%n2%n == elemB%n1%n .AND. &
+                elemA%n3%n == elemB%n3%n) THEN
+          elemA%e2 => elemB
+          elemB%e3 => elemA
+
+        ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
+                elemA%n3%n == elemB%n1%n) THEN
+          elemA%e2 => elemB
+          elemB%e1 => elemA
+
+        ELSEIF (elemA%n2%n == elemB%n3%n .AND. &
+                elemA%n3%n == elemB%n2%n) THEN
+          elemA%e2 => elemB
+          elemB%e2 => elemA
+
+        END IF
+
+
+      END IF
+
+      !Check direction 3
+      IF (.NOT. ASSOCIATED(elemA%e3)) THEN
+        IF     (elemA%n3%n == elemB%n1%n .AND. &
+                elemA%n1%n == elemB%n3%n) THEN
+          elemA%e3 => elemB
+          elemB%e3 => elemA
+
+        ELSEIF (elemA%n3%n == elemB%n2%n .AND. &
+                elemA%n1%n == elemB%n1%n) THEN
+          elemA%e3 => elemB
+          elemB%e1 => elemA
+
+        ELSEIF (elemA%n3%n == elemB%n3%n .AND. &
+                elemA%n1%n == elemB%n2%n) THEN
+          elemA%e3 => elemB
+          elemB%e2 => elemA
+
+        END IF
+
+
+      END IF
+
+    END SUBROUTINE connectedTriaTria
+
+    SUBROUTINE connectedQuadEdge(elemA, elemB)
+      IMPLICIT NONE
+
+      CLASS(meshVol2DCartQuad), INTENT(inout), TARGET:: elemA
+      CLASS(meshEdge2DCart),    INTENT(inout), TARGET:: elemB
+
+      !Check direction 1
+      IF (.NOT. ASSOCIATED(elemA%e1)) THEN
+        IF     (elemA%n1%n == elemB%n1%n .AND. &
+                elemA%n2%n == elemB%n2%n) THEN
+          elemA%e1 => elemB
+          elemB%e2 => elemA
+
+        ELSEIF (elemA%n1%n == elemB%n2%n .AND. &
+                elemA%n2%n == elemB%n1%n) THEN
+          elemA%e1 => elemB
+          elemB%e1 => elemA
+
+        END IF
+
+      END IF
+
+      !Check direction 2
+      IF (.NOT. ASSOCIATED(elemA%e2)) THEN
+        IF     (elemA%n2%n == elemB%n1%n .AND. &
+                elemA%n3%n == elemB%n2%n) THEN
+          elemA%e2 => elemB
+          elemB%e2 => elemA
+
+        ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
+                elemA%n3%n == elemB%n1%n) THEN
+          elemA%e2 => elemB
+          elemB%e1 => elemA
+
+        END IF
+
+      END IF
+
+      !Check direction 3
+      IF (.NOT. ASSOCIATED(elemA%e3)) THEN
+        IF     (elemA%n3%n == elemB%n1%n .AND. &
+                elemA%n4%n == elemB%n2%n) THEN
+          elemA%e3 => elemB
+          elemB%e2 => elemA
+
+        ELSEIF (elemA%n3%n == elemB%n2%n .AND. &
+                elemA%n4%n == elemB%n1%n) THEN
+          elemA%e3 => elemB
+          elemB%e1 => elemA
+
+        END IF
+
+      END IF
+
+      !Check direction 4
+      IF (.NOT. ASSOCIATED(elemA%e4)) THEN
+        IF     (elemA%n4%n == elemB%n1%n .AND. &
+                elemA%n1%n == elemB%n2%n) THEN
+          elemA%e4 => elemB
+          elemB%e2 => elemA
+
+        ELSEIF (elemA%n4%n == elemB%n2%n .AND. &
+                elemA%n1%n == elemB%n1%n) THEN
+          elemA%e4 => elemB
+          elemB%e1 => elemA
+
+        END IF
+
+      END IF
+
+    END SUBROUTINE connectedQuadEdge
+
+    SUBROUTINE connectedTriaEdge(elemA, elemB)
+      IMPLICIT NONE
+
+      CLASS(meshVol2DCartTria), INTENT(inout), TARGET:: elemA
+      CLASS(meshEdge2DCart),    INTENT(inout), TARGET:: elemB
+
+      !Check direction 1
+      IF (.NOT. ASSOCIATED(elemA%e1)) THEN
+        IF     (elemA%n1%n == elemB%n1%n .AND. &
+                elemA%n2%n == elemB%n2%n) THEN
+          elemA%e1 => elemB
+          elemB%e2 => elemA
+
+        ELSEIF (elemA%n1%n == elemB%n2%n .AND. &
+                elemA%n2%n == elemB%n1%n) THEN
+          elemA%e1 => elemB
+          elemB%e1 => elemA
+
+        END IF
+
+      END IF
+
+      !Check direction 2
+      IF (.NOT. ASSOCIATED(elemA%e2)) THEN
+        IF     (elemA%n2%n == elemB%n1%n .AND. &
+                elemA%n3%n == elemB%n2%n) THEN
+          elemA%e2 => elemB
+          elemB%e2 => elemA
+
+        ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
+                elemA%n3%n == elemB%n1%n) THEN
+          elemA%e2 => elemB
+          elemB%e1 => elemA
+
+        END IF
+
+      END IF
+
+      !Check direction 3
+      IF (.NOT. ASSOCIATED(elemA%e3)) THEN
+        IF     (elemA%n3%n == elemB%n1%n .AND. &
+                elemA%n1%n == elemB%n2%n) THEN
+          elemA%e3 => elemB
+          elemB%e2 => elemA
+
+        ELSEIF (elemA%n3%n == elemB%n2%n .AND. &
+                elemA%n1%n == elemB%n1%n) THEN
+          elemA%e3 => elemB
+          elemB%e1 => elemA
+
+        END IF
+
+      END IF
+
+    END SUBROUTINE connectedTriaEdge
+
+    SUBROUTINE constructGlobalK(K, elem)
+      IMPLICIT NONE
+
+      REAL(8), INTENT(inout):: K(1:,1:)
+      CLASS(meshVol), INTENT(in):: elem
+      REAL(8), ALLOCATABLE:: localK(:,:)
+      INTEGER:: nNodes, i, j
+      INTEGER, ALLOCATABLE:: n(:)
+
+      SELECT TYPE(elem)
+      TYPE IS(meshVol2DCartQuad)
+        nNodes = 4
+        ALLOCATE(localK(1:nNodes,1:nNodes))
+        localK = elem%elemK()
+        ALLOCATE(n(1:nNodes))
+        n = (/ elem%n1%n, elem%n2%n, &
+               elem%n3%n, elem%n4%n /)
+
+      TYPE IS(meshVol2DCartTria)
+        nNodes = 3
+        ALLOCATE(localK(1:nNodes,1:nNodes))
+        localK = elem%elemK()
+        ALLOCATE(n(1:nNodes))
+        n = (/ elem%n1%n, elem%n2%n, elem%n3%n /)
+
+      CLASS DEFAULT
+        nNodes = 0
+        ALLOCATE(localK(1:1, 1:1))
+        localK = 0.D0
+        ALLOCATE(n(1:1))
+        n = 0
+
+      END SELECT
+      
+      DO i = 1, nNodes
+        DO j = 1, nNodes
+          K(n(i), n(j)) = K(n(i), n(j)) + localK(i, j)
+        END DO
+      END DO
+
+    END SUBROUTINE constructGlobalK
+
+END MODULE moduleMesh2DCartRead