First implementation of 1D radial case. Only charged particles taked

into account (as in 1D Cartesian case).

The 1D Cathode example case has been modified, having now 2 input files:
  - inputCart.json: Used for Cartesian coordinates
  - inputRad.json:  Used for Radial coordinates

Pusher is a Boris pusher but without z direction.

runs/1D_Cathode/inputCart.json

@@ -25,12 +25,11 @@
     {"name": "Infinite", "type": "absorption", "physicalSurface": 2}
   ],
   "boundaryEM": [
-    {"name": "Cathode",  "type": "dirichlet", "potential":    0.0, "physicalSurface": 1},
+    {"name": "Cathode",  "type": "dirichlet", "potential":  -10.0, "physicalSurface": 1}
-    {"name": "Infinite", "type": "dirichlet", "potential":  100.0, "physicalSurface": 2}
   ],
   "inject": [
     {"name": "Cathode Electron", "species": "Electron", "flow":  9.0e-5, "units": "A", "v": 27500.0, "T": [2500.0, 2500.0, 2500.0],
-                                 "velDist": ["Delta", "Maxwellian", "Maxwellian"],     "n": [ 1, 0, 0], "physicalSurface": 1}
+                                 "velDist": ["Maxwellian", "Maxwellian", "Maxwellian"],     "n": [ 1, 0, 0], "physicalSurface": 1}
   ],
   "case": {
     "tau":  [1.0e-11],

runs/1D_Cathode/inputRad.json

@@ -0,0 +1,45 @@
+
+{
+  "output": {
+    "path": "./runs/1D_Cathode/",
+    "triggerOutput": 100,
+    "cpuTime": false,
+    "numColl": false,
+    "EMField": true
+  },
+  "reference": {
+    "density": 1.0e16,
+    "mass": 9.109e-31,
+    "temperature": 2500.0
+  },
+  "geometry": {
+    "type":     "1DRad",
+    "meshType": "gmsh",
+      "meshFile": "mesh.msh"
+  },
+  "species": [
+    {"name": "Electron", "type": "charged", "mass": 9.109e-31, "charge":-1.0, "weight": 1.0e1}
+  ],
+  "boundary": [ 
+    {"name": "Cathode",  "type": "absorption", "physicalSurface": 1},
+    {"name": "Infinite", "type": "absorption", "physicalSurface": 2}
+  ],
+  "boundaryEM": [
+    {"name": "Cathode",  "type": "dirichlet", "potential":  -10.0, "physicalSurface": 1}
+  ],
+  "inject": [
+    {"name": "Cathode Electron", "species": "Electron", "flow":  9.0e-5, "units": "A", "v": 27500.0, "T": [2500.0, 2500.0, 2500.0],
+                                 "velDist": ["Maxwellian", "Maxwellian", "Maxwellian"],     "n": [ 1, 0, 0], "physicalSurface": 1}
+  ],
+  "case": {
+    "tau":  [1.0e-11],
+    "time": 1.0e-6,
+    "pusher": ["1DRadCharged"],
+    "EMSolver": "Electrostatic"
+  },
+  "parallel": {
+    "OpenMP":{
+      "nThreads": 1
+    }
+  }
+}

src/makefile

@@ -1,11 +1,12 @@
 OBJECTS = $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleCompTime.o $(OBJDIR)/moduleSolver.o \
 			  	$(OBJDIR)/moduleSpecies.o $(OBJDIR)/moduleInject.o $(OBJDIR)/moduleInput.o \
 			  	$(OBJDIR)/moduleErrors.o $(OBJDIR)/moduleList.o $(OBJDIR)/moduleOutput.o \
-			  	$(OBJDIR)/moduleMeshCyl.o	$(OBJDIR)/moduleMeshCylRead.o	$(OBJDIR)/moduleMeshCylBoundary.o \
 			  	$(OBJDIR)/moduleBoundary.o $(OBJDIR)/moduleCaseParam.o $(OBJDIR)/moduleRefParam.o \
 			  	$(OBJDIR)/moduleCollisions.o $(OBJDIR)/moduleTable.o $(OBJDIR)/moduleParallel.o \
-			  	$(OBJDIR)/moduleEM.o $(OBJDIR)/moduleMesh1D.o $(OBJDIR)/moduleMesh1DRead.o \
+			  	$(OBJDIR)/moduleEM.o \
-					$(OBJDIR)/moduleMesh1DBoundary.o
+			  	$(OBJDIR)/moduleMeshCyl.o	$(OBJDIR)/moduleMeshCylRead.o	$(OBJDIR)/moduleMeshCylBoundary.o \
+          $(OBJDIR)/moduleMesh1DCart.o $(OBJDIR)/moduleMesh1DCartRead.o $(OBJDIR)/moduleMesh1DCartBoundary.o \
+          $(OBJDIR)/moduleMesh1DRad.o  $(OBJDIR)/moduleMesh1DRadRead.o  $(OBJDIR)/moduleMesh1DRadBoundary.o
 
 all: $(OUTPUT)
 

src/modules/mesh/1DCart/makefile

@@ -1,11 +1,11 @@
-all: moduleMesh1D.o moduleMesh1DBoundary.o moduleMesh1DRead.o
+all: moduleMesh1DCart.o moduleMesh1DCartBoundary.o moduleMesh1DCartRead.o
 
-moduleMesh1D.o: moduleMesh1D.f90
+moduleMesh1DCart.o: moduleMesh1DCart.f90
 	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
 
-moduleMesh1DBoundary.o: moduleMesh1D.o moduleMesh1DBoundary.f90
+moduleMesh1DCartBoundary.o: moduleMesh1DCart.o moduleMesh1DCartBoundary.f90
 	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
 
-moduleMesh1DRead.o: moduleMesh1D.o moduleMesh1DBoundary.o moduleMesh1DRead.f90
+moduleMesh1DCartRead.o: moduleMesh1DCart.o moduleMesh1DCartBoundary.o moduleMesh1DCartRead.f90
 	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
 

src/modules/mesh/1DCart/moduleMesh1DCart.f90

@@ -2,40 +2,40 @@
 !              x == x
 !              y == unused
 !              z == unused
-MODULE moduleMesh1D
+MODULE moduleMesh1DCart
   USE moduleMesh
   IMPLICIT NONE
 
-  TYPE, PUBLIC, EXTENDS(meshNode):: meshNode1D
+  TYPE, PUBLIC, EXTENDS(meshNode):: meshNode1DCart
     !Element coordinates
     REAL(8):: x = 0.D0
     CONTAINS
-      PROCEDURE, PASS:: init => initNode1D
+      PROCEDURE, PASS:: init => initNode1DCart
-      PROCEDURE, PASS:: getCoordinates => getCoord1D
+      PROCEDURE, PASS:: getCoordinates => getCoord1DCart
 
-  END TYPE meshNode1D
+  END TYPE meshNode1DCart
 
-  TYPE, PUBLIC, ABSTRACT, EXTENDS(meshEdge):: meshEdge1D
+  TYPE, PUBLIC, ABSTRACT, EXTENDS(meshEdge):: meshEdge1DCart
     !Element coordinates
     REAL(8):: x = 0.D0
     !Connectivity to nodes
     CLASS(meshNode), POINTER:: n1 => NULL()
     CONTAINS
-      PROCEDURE, PASS:: init     => initEdge1D
+      PROCEDURE, PASS:: init     => initEdge1DCart
-      PROCEDURE, PASS:: getNodes => getNodes1D
+      PROCEDURE, PASS:: getNodes => getNodes1DCart
-      PROCEDURE, PASS:: randPos  => randPos1D
+      PROCEDURE, PASS:: randPos  => randPos1DCart
 
-  END TYPE meshEdge1D
+  END TYPE meshEdge1DCart
 
-  TYPE, PUBLIC, ABSTRACT, EXTENDS(meshVol):: meshVol1D
+  TYPE, PUBLIC, ABSTRACT, EXTENDS(meshVol):: meshVol1DCart
     CONTAINS
-      PROCEDURE, PASS:: detJac => detJ1D
+      PROCEDURE, PASS:: detJac => detJ1DCart
-      PROCEDURE, PASS:: invJac => invJ1D
+      PROCEDURE, PASS:: invJac => invJ1DCart
       PROCEDURE(fPsi_interface), DEFERRED, NOPASS:: fPsi
       PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi
       PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer
   
-  END TYPE meshVol1D
+  END TYPE meshVol1DCart
 
   ABSTRACT INTERFACE
     PURE FUNCTION fPsi_interface(xi) RESULT(fPsi)
@@ -51,8 +51,8 @@ MODULE moduleMesh1D
     END FUNCTION dPsi_interface
 
     PURE SUBROUTINE partialDer_interface(self, dPsi, dx)
-      IMPORT meshVol1D
+      IMPORT meshVol1DCart
-      CLASS(meshVol1D), INTENT(in):: self
+      CLASS(meshVol1DCart), INTENT(in):: self
       REAL(8), INTENT(in):: dPsi(1:,1:)
       REAL(8), INTENT(out), DIMENSION(1):: dx
 
@@ -60,7 +60,7 @@ MODULE moduleMesh1D
 
   END INTERFACE
 
-  TYPE, PUBLIC, EXTENDS(meshVol1D):: meshVol1DSegm
+  TYPE, PUBLIC, EXTENDS(meshVol1DCart):: meshVol1DCartSegm
     !Element coordinates
     REAL(8):: x(1:2)
     !Connectivity to nodes
@@ -69,7 +69,7 @@ MODULE moduleMesh1D
     CLASS(*), POINTER:: e1 => NULL(), e2 => NULL()
     REAL(8):: arNodes(1:2)
     CONTAINS
-      PROCEDURE, PASS::   init        => initVol1DSegm
+      PROCEDURE, PASS::   init        => initVol1DCartSegm
       PROCEDURE, PASS::   area        => areaSegm
       PROCEDURE, NOPASS:: fPsi        => fPsiSegm
       PROCEDURE, NOPASS:: dPsi        => dPsiSegm
@@ -85,17 +85,17 @@ MODULE moduleMesh1D
       PROCEDURE, PASS::   nextElement => nextElementSegm
       PROCEDURE, PASS::   resetOutput => resetOutputSegm
 
-  END TYPE meshVol1DSegm
+  END TYPE meshVol1DCartSegm
 
   CONTAINS
     !NODE FUNCTIONS
     !Init node element
-    SUBROUTINE initNode1D(self, n, r)
+    SUBROUTINE initNode1DCart(self, n, r)
       USE moduleSpecies
       USE moduleRefParam
       IMPLICIT NONE
 
-      CLASS(meshNode1D), INTENT(out):: self
+      CLASS(meshNode1DCart), INTENT(out):: self
       INTEGER, INTENT(in):: n
       REAL(8), INTENT(in):: r(1:3)
 
@@ -107,24 +107,24 @@ MODULE moduleMesh1D
       !Allocates output
       ALLOCATE(self%output(1:nSpecies))
 
-    END SUBROUTINE initNode1D
+    END SUBROUTINE initNode1DCart
 
-    PURE FUNCTION getCoord1D(self) RESULT(r)
+    PURE FUNCTION getCoord1DCart(self) RESULT(r)
       IMPLICIT NONE
 
-      CLASS(meshNode1D), INTENT(in):: self
+      CLASS(meshNode1DCart), INTENT(in):: self
       REAL(8):: r(1:3)
 
       r = (/ self%x, 0.D0, 0.D0 /)
 
-    END FUNCTION getCoord1D
+    END FUNCTION getCoord1DCart
 
     !EDGE FUNCTIONS
     !Inits edge element
-    SUBROUTINE initEdge1D(self, n, p, bt, physicalSurface)
+    SUBROUTINE initEdge1DCart(self, n, p, bt, physicalSurface)
       IMPLICIT NONE
 
-      CLASS(meshEdge1D), INTENT(out):: self
+      CLASS(meshEdge1DCart), INTENT(out):: self
       INTEGER, INTENT(in):: n
       INTEGER, INTENT(in):: p(:)
       INTEGER, INTENT(in):: bt
@@ -145,37 +145,37 @@ MODULE moduleMesh1D
       !Physical Surface
       self%physicalSurface = physicalSurface
 
-    END SUBROUTINE initEdge1D
+    END SUBROUTINE initEdge1DCart
 
     !Get nodes from edge
-    PURE FUNCTION getNodes1D(self) RESULT(n)
+    PURE FUNCTION getNodes1DCart(self) RESULT(n)
       IMPLICIT NONE
 
-      CLASS(meshEdge1D), INTENT(in):: self
+      CLASS(meshEdge1DCart), INTENT(in):: self
       INTEGER, ALLOCATABLE:: n(:)
 
       ALLOCATE(n(1))
       n = (/ self%n1%n /)
 
-    END FUNCTION getNodes1D
+    END FUNCTION getNodes1DCart
 
     !Calculates a 'random' position in edge
-    FUNCTION randPos1D(self) RESULT(r)
+    FUNCTION randPos1DCart(self) RESULT(r)
-      CLASS(meshEdge1D), INTENT(in):: self
+      CLASS(meshEdge1DCart), INTENT(in):: self
       REAL(8):: r(1:3)
 
       r = (/ self%x, 0.D0, 0.D0 /)
 
-    END FUNCTION randPos1D
+    END FUNCTION randPos1DCart
 
     !VOLUME FUNCTIONS
     !SEGMENT FUNCTIONS
     !Init segment element
-    SUBROUTINE initVol1DSegm(self, n, p)
+    SUBROUTINE initVol1DCartSegm(self, n, p)
       USE moduleRefParam
       IMPLICIT NONE
 
-      CLASS(meshVol1DSegm), INTENT(out):: self
+      CLASS(meshVol1DCartSegm), INTENT(out):: self
       INTEGER, INTENT(in):: n
       INTEGER, INTENT(in):: p(:)
       REAL(8), DIMENSION(1:3):: r1, r2
@@ -197,13 +197,13 @@ MODULE moduleMesh1D
 
       CALL OMP_INIT_LOCK(self%lock)
 
-    END SUBROUTINE initVol1DSegm
+    END SUBROUTINE initVol1DCartSegm
 
     !Computes element area
     PURE SUBROUTINE areaSegm(self)
       IMPLICIT NONE
 
-      CLASS(meshVol1DSegm), INTENT(inout):: self
+      CLASS(meshVol1DCartSegm), INTENT(inout):: self
       REAL(8):: l !element length
       REAL(8):: fPsi(1:2)
       REAL(8):: detJ
@@ -254,7 +254,7 @@ MODULE moduleMesh1D
     PURE SUBROUTINE partialDerSegm(self, dPsi, dx)
       IMPLICIT NONE
 
-      CLASS(meshVol1DSegm), INTENT(in):: self
+      CLASS(meshVol1DCartSegm), INTENT(in):: self
       REAL(8), INTENT(in):: dPsi(1:,1:)
       REAL(8), INTENT(out), DIMENSION(1):: dx
 
@@ -266,7 +266,7 @@ MODULE moduleMesh1D
     PURE FUNCTION elemKSegm(self) RESULT(ke)
       IMPLICIT NONE
 
-      CLASS(meshVol1DSegm), INTENT(in):: self
+      CLASS(meshVol1DCartSegm), INTENT(in):: self
       REAL(8):: ke(1:2,1:2)
       REAL(8):: Xii(1:3)
       REAL(8):: dPsi(1:1, 1:2)
@@ -285,7 +285,7 @@ MODULE moduleMesh1D
     PURE FUNCTION elemFSegm(self, source) RESULT(localF)
       IMPLICIT NONE
 
-      CLASS(meshVol1DSegm), INTENT(in):: self
+      CLASS(meshVol1DCartSegm), INTENT(in):: self
       REAL(8), INTENT(in):: source(1:)
       REAL(8), ALLOCATABLE:: localF(:)
       REAL(8):: fPsi(1:2)
@@ -326,7 +326,7 @@ MODULE moduleMesh1D
       USE moduleSpecies
       IMPLICIT NONE
 
-      CLASS(meshVol1DSegm), INTENT(in):: self
+      CLASS(meshVol1DCartSegm), INTENT(in):: self
       CLASS(particle), INTENT(in):: part
       TYPE(outputNode), POINTER:: vertex
       REAL(8):: w_p(1:2)
@@ -351,7 +351,7 @@ MODULE moduleMesh1D
     PURE FUNCTION gatherEFSegm(self, xi) RESULT(EF)
       IMPLICIT NONE
 
-      CLASS(meshVol1DSegm), INTENT(in):: self
+      CLASS(meshVol1DCartSegm), INTENT(in):: self
       REAL(8), INTENT(in):: xi(1:3)
       REAL(8):: dPsi(1, 1:2)
       REAL(8):: phi(1:2)
@@ -373,7 +373,7 @@ MODULE moduleMesh1D
     PURE FUNCTION getNodesSegm(self) RESULT(n)
       IMPLICIT NONE
 
-      CLASS(meshVol1DSegm), INTENT(in):: self
+      CLASS(meshVol1DCartSegm), INTENT(in):: self
       INTEGER, ALLOCATABLE:: n(:)
 
       ALLOCATE(n(1:2))
@@ -384,7 +384,7 @@ MODULE moduleMesh1D
     PURE FUNCTION phy2logSegm(self, r) RESULT(xN)
       IMPLICIT NONE
 
-      CLASS(meshVol1DSegm), INTENT(in):: self
+      CLASS(meshVol1DCartSegm), INTENT(in):: self
       REAL(8), INTENT(in):: r(1:3)
       REAL(8):: xN(1:3)
 
@@ -397,7 +397,7 @@ MODULE moduleMesh1D
     SUBROUTINE nextElementSegm(self, xi, nextElement)
       IMPLICIT NONE
 
-      CLASS(meshVol1DSegm), INTENT(in):: self
+      CLASS(meshVol1DCartSegm), INTENT(in):: self
       REAL(8), INTENT(in):: xi(1:3)
       CLASS(*), POINTER, INTENT(out):: nextElement
       
@@ -418,7 +418,7 @@ MODULE moduleMesh1D
       USE moduleOutput
       IMPLICIT NONE
 
-      CLASS(meshVol1DSegm), INTENT(inout):: self
+      CLASS(meshVol1DCartSegm), INTENT(inout):: self
       INTEGER:: k
 
       DO k = 1, nSpecies
@@ -437,10 +437,10 @@ MODULE moduleMesh1D
 
     !COMMON FUNCTIONS FOR 1D VOLUME ELEMENTS
     !Computes the element Jacobian determinant
-    PURE FUNCTION detJ1D(self, xi, dPsi_in) RESULT(dJ)
+    PURE FUNCTION detJ1DCart(self, xi, dPsi_in) RESULT(dJ)
       IMPLICIT NONE
 
-      CLASS(meshVol1D), INTENT(in):: self
+      CLASS(meshVol1DCart), INTENT(in):: self
       REAL(8), INTENT(in):: xi(1:3)
       REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:,1:)
       REAL(8), ALLOCATABLE:: dPsi(:,:)
@@ -458,13 +458,13 @@ MODULE moduleMesh1D
       CALL self%partialDer(dPsi, dx)
       dJ = dx(1)
 
-    END FUNCTION detJ1D
+    END FUNCTION detJ1DCart
 
     !Computes the invers Jacobian
-    PURE FUNCTION invJ1D(self, xi, dPsi_in) RESULT(invJ)
+    PURE FUNCTION invJ1DCart(self, xi, dPsi_in) RESULT(invJ)
       IMPLICIT NONE
 
-      CLASS(meshVol1D), INTENT(in):: self
+      CLASS(meshVol1DCart), INTENT(in):: self
       REAL(8), INTENT(in):: xi(1:3)
       REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:,1:)
       REAL(8), ALLOCATABLE:: dPsi(:,:)
@@ -482,8 +482,8 @@ MODULE moduleMesh1D
       CALL self%partialDer(dPsi, dx)
       invJ = 1.D0/dx(1)
 
-    END FUNCTION invJ1D
+    END FUNCTION invJ1DCart
 
 
-END MODULE moduleMesh1D
+END MODULE moduleMesh1DCart
 

src/modules/mesh/1DCart/moduleMesh1DCartBoundary.f90

@@ -1,24 +1,24 @@
-MODULE moduleMesh1DBoundary
+MODULE moduleMesh1DCartBoundary
-  USE moduleMesh1D
+  USE moduleMesh1DCart
 
-  TYPE, PUBLIC, EXTENDS(meshEdge1D):: meshEdge1DRef
+  TYPE, PUBLIC, EXTENDS(meshEdge1DCart):: meshEdge1DCartRef
     CONTAINS
       PROCEDURE, PASS:: fBoundary => reflection
 
-  END TYPE meshEdge1DRef
+  END TYPE meshEdge1DCartRef
 
-  TYPE, PUBLIC, EXTENDS(meshEdge1D):: meshEdge1DAbs
+  TYPE, PUBLIC, EXTENDS(meshEdge1DCart):: meshEdge1DCartAbs
     CONTAINS
       PROCEDURE, PASS:: fBoundary => absorption
 
-  END TYPE meshEdge1DAbs
+  END TYPE meshEdge1DCartAbs
 
   CONTAINS
     SUBROUTINE reflection(self, part)
       USE moduleSpecies
       IMPLICIT NONE
 
-      CLASS(meshEdge1DRef), INTENT(inout):: self
+      CLASS(meshEdge1DCartRef), INTENT(inout):: self
       CLASS(particle), INTENT(inout):: part
 
       part%v(1) = -part%v(1)
@@ -30,11 +30,11 @@ MODULE moduleMesh1DBoundary
       USE moduleSpecies
       IMPLICIT NONE
 
-      CLASS(meshEdge1DAbs), INTENT(inout):: self
+      CLASS(meshEdge1DCartAbs), INTENT(inout):: self
       CLASS(particle), INTENT(inout):: part
 
       part%n_in = .FALSE.
 
     END SUBROUTINE absorption
 
-END MODULE moduleMesh1DBoundary
+END MODULE moduleMesh1DCartBoundary

src/modules/mesh/1DCart/moduleMesh1DCartRead.f90

@@ -0,0 +1,268 @@
+MODULE moduleMesh1DCartRead
+  USE moduleMesh
+  USE moduleMesh1DCart
+  USE moduleMesh1DCartBoundary
+  
+  !TODO: make this abstract to allow different mesh formats
+  TYPE, EXTENDS(meshGeneric):: mesh1DCartGeneric
+    CONTAINS
+      PROCEDURE, PASS:: init     => init1DCartMesh
+      PROCEDURE, PASS:: readMesh => readMesh1DCart
+
+  END TYPE
+
+  INTERFACE connected
+    MODULE PROCEDURE connectedVolVol, connectedVolEdge
+
+  END INTERFACE connected
+
+  CONTAINS
+    !Init 1D mesh
+    SUBROUTINE init1DCartMesh(self, meshFormat)
+      USE moduleMesh
+      USE moduleErrors
+      IMPLICIT NONE
+
+      CLASS(mesh1DCartGeneric), 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.", "init1DCart")
+        
+      END SELECT
+
+    END SUBROUTINE init1DCartMesh
+
+    !Reads 1D mesh
+    SUBROUTINE readMesh1DCart(self, filename)
+      USE  moduleBoundary
+      IMPLICIT NONE
+
+      CLASS(mesh1DCartGeneric), INTENT(inout):: self
+      CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
+      REAL(8):: x
+      INTEGER:: p(1:2)
+      INTEGER:: e, et, n, eTemp, elemType, bt
+      INTEGER:: totalNumElem
+      INTEGER:: boundaryType
+
+      !Open file mesh
+      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 coordinates. Only relevant for x
+      DO e = 1, self%numNodes
+        READ(10, *) n, x
+        ALLOCATE(meshNode1DCart:: self%nodes(n)%obj)
+        CALL self%nodes(n)%obj%init(n, (/ x, 0.D0, 0.D0 /))
+
+      END DO
+      !Skips comments
+      READ(10, *)
+      READ(10, *)
+      !Reads the total number of elements (edges+vol)
+      READ(10, *) totalNumElem
+      self%numEdges = 0
+      DO e = 1, totalNumElem
+        READ(10, *) eTemp, elemType
+        IF (elemType == 15) THEN !15 is physical node in GMSH
+          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 beginning of reading 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)
+        !Associate boundary condition 
+        bt = getBoundaryId(boundaryType)
+        SELECT CASE(boundary(bt)%obj%boundaryType)
+        CASE ('reflection')
+          ALLOCATE(meshEdge1DCartRef:: self%edges(e)%obj)
+
+        CASE ('absorption')
+          ALLOCATE(meshEdge1DCartAbs:: self%edges(e)%obj)
+
+        END SELECT
+
+        CALL self%edges(e)%obj%init(n, p(1:1), bt, boundaryType)
+
+      END DO
+
+      !Read and initialize volumes
+      DO e = 1, self%numVols
+        READ(10, *) n, elemType, eTemp, eTemp, eTemp, p(1:2)
+        ALLOCATE(meshVol1DCartSegm:: self%vols(e)%obj)
+        CALL self%vols(e)%obj%init(n - self%numEdges, p(1:2))
+
+      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 betwen 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 readMesh1DCart
+
+    SUBROUTINE connectedVolVol(elemA, elemB)
+      IMPLICIT NONE
+
+      CLASS(meshVol), INTENT(inout):: elemA
+      CLASS(meshVol), INTENT(inout):: elemB
+
+        SELECT TYPE(elemA)
+        TYPE IS(meshVol1DCartSegm)
+          SELECT TYPE(elemB)
+          TYPE IS(meshVol1DCartSegm)
+            CALL connectedSegmSegm(elemA, elemB)
+
+          END SELECT
+
+        END SELECT
+
+    END SUBROUTINE connectedVolVol
+
+    SUBROUTINE connectedSegmSegm(elemA, elemB)
+      IMPLICIT NONE
+
+      CLASS(meshVol1DCartSegm), INTENT(inout), TARGET:: elemA
+      CLASS(meshVol1DCartSegm), INTENT(inout), TARGET:: elemB
+
+      IF (.NOT. ASSOCIATED(elemA%e1) .AND. &
+          elemA%n2%n == elemB%n1%n) THEN
+        elemA%e1 => elemB
+        elemB%e2 => elemA
+        
+      END IF
+
+      IF (.NOT. ASSOCIATED(elemA%e2) .AND. &
+          elemA%n1%n == elemB%n2%n) THEN
+
+        elemA%e2 => elemB
+        elemB%e1 => elemA
+
+      END IF
+
+    END SUBROUTINE connectedSegmSegm
+
+    SUBROUTINE connectedVolEdge(elemA, elemB)
+      IMPLICIT NONE
+
+      CLASS(meshVol),  INTENT(inout):: elemA
+      CLASS(meshEdge), INTENT(inout):: elemB
+
+      SELECT TYPE(elemA)
+      TYPE IS (meshVol1DCartSegm)
+        SELECT TYPE(elemB)
+        CLASS IS(meshEdge1DCart)
+          CALL connectedSegmEdge(elemA, elemB)
+
+        END SELECT
+
+      END SELECT
+
+    END SUBROUTINE connectedVolEdge
+
+    SUBROUTINE connectedSegmEdge(elemA, elemB)
+      IMPLICIT NONE
+
+      CLASS(meshVol1DCartSegm), INTENT(inout), TARGET:: elemA
+      CLASS(meshEdge1DCart),    INTENT(inout), TARGET:: elemB
+
+      IF (.NOT. ASSOCIATED(elemA%e1) .AND. &
+          elemA%n2%n == elemB%n1%n) THEN
+
+        elemA%e1 => elemB
+        elemB%e2 => elemA
+
+      END IF
+
+      IF (.NOT. ASSOCIATED(elemA%e2) .AND. &
+          elemA%n1%n == elemB%n1%n) THEN
+
+        elemA%e2 => elemB
+        elemB%e1 => elemA
+
+      END IF
+
+    END SUBROUTINE connectedSegmEdge
+
+    SUBROUTINE constructGlobalK(K, elem)
+      IMPLICIT NONE
+
+      REAL(8), INTENT(inout):: K(1:,1:)
+      CLASS(meshVol), INTENT(in):: elem
+      REAL(8):: localK(1:2,1:2)
+      INTEGER:: i, j
+      INTEGER:: n(1:2)
+
+      SELECT TYPE(elem)
+      TYPE IS(meshVol1DCartSegm)
+        localK = elem%elemK()
+        n = (/ elem%n1%n, elem%n2%n /)
+
+      CLASS DEFAULT
+        n = 0
+        localK = 0.D0
+
+      END SELECT
+      
+      DO i = 1, 2
+        DO j = 1, 2
+          K(n(i), n(j)) = K(n(i), n(j)) + localK(i, j)
+        END DO
+      END DO
+
+    END SUBROUTINE constructGlobalK
+
+END MODULE moduleMesh1DCartRead

src/modules/mesh/1DRad/makefile

@@ -0,0 +1,11 @@
+all: moduleMesh1DRad.o moduleMesh1DRadBoundary.o moduleMesh1DRadRead.o
+
+moduleMesh1DRad.o: moduleMesh1DRad.f90
+	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
+
+moduleMesh1DRadBoundary.o: moduleMesh1DRad.o moduleMesh1DRadBoundary.f90
+	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
+
+moduleMesh1DRadRead.o: moduleMesh1DRad.o moduleMesh1DRadBoundary.o moduleMesh1DRadRead.f90
+	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
+

src/modules/mesh/1DRad/moduleMesh1DRad.f90

@@ -0,0 +1,500 @@
+!moduleMesh1DRad: 1D radial module
+!              x == r
+!              y == theta (unused)
+!              z == unused
+MODULE moduleMesh1DRad
+  USE moduleMesh
+  IMPLICIT NONE
+
+  TYPE, PUBLIC, EXTENDS(meshNode):: meshNode1DRad
+    !Element coordinates
+    REAL(8):: r = 0.D0
+    CONTAINS
+      PROCEDURE, PASS:: init => initNode1DRad
+      PROCEDURE, PASS:: getCoordinates => getCoord1DRad
+
+  END TYPE meshNode1DRad
+
+  TYPE, PUBLIC, ABSTRACT, EXTENDS(meshEdge):: meshEdge1DRad
+    !Element coordinates
+    REAL(8):: r = 0.D0
+    !Connectivity to nodes
+    CLASS(meshNode), POINTER:: n1 => NULL()
+    CONTAINS
+      PROCEDURE, PASS:: init     => initEdge1DRad
+      PROCEDURE, PASS:: getNodes => getNodes1DRad
+      PROCEDURE, PASS:: randPos  => randPos1DRad
+
+  END TYPE meshEdge1DRad
+
+  TYPE, PUBLIC, ABSTRACT, EXTENDS(meshVol):: meshVol1DRad
+    CONTAINS
+      PROCEDURE, PASS:: detJac => detJ1DRad
+      PROCEDURE, PASS:: invJac => invJ1DRad
+      PROCEDURE(fPsi_interface), DEFERRED, NOPASS:: fPsi
+      PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi
+      PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer
+  
+  END TYPE meshVol1Drad
+
+  ABSTRACT INTERFACE
+    PURE FUNCTION fPsi_interface(xi) RESULT(fPsi)
+      REAL(8), INTENT(in):: xi(1:3)
+      REAL(8), ALLOCATABLE:: fPsi(:)
+
+    END FUNCTION fPsi_interface
+
+    PURE FUNCTION dPsi_interface(xi) RESULT(dPsi)
+      REAL(8), INTENT(in):: xi(1:3)
+      REAL(8), ALLOCATABLE:: dPsi(:,:)
+
+    END FUNCTION dPsi_interface
+
+    PURE SUBROUTINE partialDer_interface(self, dPsi, dx)
+      IMPORT meshVol1DRad
+
+      CLASS(meshVol1DRad), INTENT(in):: self
+      REAL(8), INTENT(in):: dPsi(1:,1:)
+      REAL(8), INTENT(out), DIMENSION(1):: dx
+
+    END SUBROUTINE partialDer_interface
+
+  END INTERFACE
+
+  TYPE, PUBLIC, EXTENDS(meshVol1DRad):: meshVol1DRadSegm
+    !Element coordinates
+    REAL(8):: r(1:2)
+    !Connectivity to nodes
+    CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL()
+    !Connectivity to adjacent elements
+    CLASS(*), POINTER:: e1 => NULL(), e2 => NULL()
+    REAL(8):: arNodes(1:2)
+    CONTAINS
+      PROCEDURE, PASS::   init        => initVol1DRadSegm
+      PROCEDURE, PASS::   area        => areaRad
+      PROCEDURE, NOPASS:: fPsi        => fPsiRad
+      PROCEDURE, NOPASS:: dPsi        => dPsiRad
+      PROCEDURE, PASS::   partialDer  => partialDerRad
+      PROCEDURE, PASS::   elemK       => elemKRad
+      PROCEDURE, PASS::   elemF       => elemFRad
+      PROCEDURE, NOPASS:: weight      => weightRad
+      PROCEDURE, NOPASS:: inside      => insideRad
+      PROCEDURE, PASS::   scatter     => scatterRad
+      PROCEDURE, PASS::   gatherEF    => gatherEFRad
+      PROCEDURE, PASS::   getNodes    => getNodesRad
+      PROCEDURE, PASS::   phy2log     => phy2logRad
+      PROCEDURE, PASS::   nextElement => nextElementRad
+      PROCEDURE, PASS::   resetOutput => resetOutputRad
+
+  END TYPE meshVol1DRadSegm
+
+  CONTAINS
+    !NODE FUNCTIONS
+    !Init node element
+    SUBROUTINE initNode1DRad(self, n, r)
+      USE moduleSpecies
+      USE moduleRefParam
+      IMPLICIT NONE
+
+      CLASS(meshNode1DRad), INTENT(out):: self
+      INTEGER, INTENT(in):: n
+      REAL(8), INTENT(in):: r(1:3)
+
+      self%n = n
+      self%r = r(1)/L_ref
+      !Node volume, to be determined in mesh
+      self%v = 0.D0
+
+      !Allocates output
+      ALLOCATE(self%output(1:nSpecies))
+
+    END SUBROUTINE initNode1DRad
+
+    PURE FUNCTION getCoord1DRad(self) RESULT(r)
+      IMPLICIT NONE
+
+      CLASS(meshNode1DRad), INTENT(in):: self
+      REAL(8):: r(1:3)
+
+      r = (/ self%r, 0.D0, 0.D0 /)
+
+    END FUNCTION getCoord1DRad
+
+    !EDGE FUNCTIONS
+    !Inits edge element
+    SUBROUTINE initEdge1DRad(self, n, p, bt, physicalSurface)
+      IMPLICIT NONE
+
+      CLASS(meshEdge1DRad), INTENT(out):: self
+      INTEGER, INTENT(in):: n
+      INTEGER, INTENT(in):: p(:)
+      INTEGER, INTENT(in):: bt
+      INTEGER, INTENT(in):: physicalSurface
+      REAL(8), DIMENSION(1:3):: r1
+
+      self%n = n
+      self%n1 => mesh%nodes(p(1))%obj
+      !Get element coordinates
+      r1 = self%n1%getCoordinates()
+
+      self%r = r1(1)
+
+      self%normal = (/ 1.D0, 0.D0, 0.D0 /)
+
+      !Boundary index
+      self%bt = bt
+      !Physical Surface
+      self%physicalSurface = physicalSurface
+
+    END SUBROUTINE initEdge1DRad
+
+    !Get nodes from edge
+    PURE FUNCTION getNodes1DRad(self) RESULT(n)
+      IMPLICIT NONE
+
+      CLASS(meshEdge1DRad), INTENT(in):: self
+      INTEGER, ALLOCATABLE:: n(:)
+
+      ALLOCATE(n(1))
+      n = (/ self%n1%n /)
+
+    END FUNCTION getNodes1DRad
+
+    !Calculates a 'random' position in edge
+    FUNCTION randPos1DRad(self) RESULT(r)
+      CLASS(meshEdge1DRad), INTENT(in):: self
+      REAL(8):: r(1:3)
+
+      r = (/ self%r, 0.D0, 0.D0 /)
+
+    END FUNCTION randPos1DRad
+
+    !VOLUME FUNCTIONS
+    !SEGMENT FUNCTIONS
+    !Init segment element
+    SUBROUTINE initVol1DRadSegm(self, n, p)
+      USE moduleRefParam
+      IMPLICIT NONE
+
+      CLASS(meshVol1DRadSegm), INTENT(out):: self
+      INTEGER, INTENT(in):: n
+      INTEGER, INTENT(in):: p(:)
+      REAL(8), DIMENSION(1:3):: r1, r2
+
+      self%n = n
+      self%n1 => mesh%nodes(p(1))%obj
+      self%n2 => mesh%nodes(p(2))%obj
+      !Get element coordinates
+      r1 = self%n1%getCoordinates()
+      r2 = self%n2%getCoordinates()
+      self%r = (/ r1(1), r2(1) /)
+
+      !Assign node volume
+      CALL self%area()
+      self%n1%v = self%n1%v + self%arNodes(1)
+      self%n2%v = self%n2%v + self%arNodes(2)
+
+      self%sigmaVrelMax = sigma_ref/L_ref**2
+
+      CALL OMP_INIT_LOCK(self%lock)
+
+    END SUBROUTINE initVol1DRadSegm
+
+    !Computes element area
+    PURE SUBROUTINE areaRad(self)
+      IMPLICIT NONE
+
+      CLASS(meshVol1DRadSegm), INTENT(inout):: self
+      REAL(8):: l !element length
+      REAL(8):: fPsi(1:2)
+      REAL(8):: r
+      REAL(8):: detJ
+      REAL(8):: Xii(1:3)
+
+      self%volume  = 0.D0
+      self%arNodes = 0.D0
+      !1 point Gauss integral
+      Xii = 0.D0
+      fPsi = self%fPsi(Xii)
+      detJ = self%detJac(Xii)
+      r    = DOT_PRODUCT(fPsi, self%r)
+      l = 2.D0*detJ
+      self%volume  =      r*l
+      self%arNodes = fPsi*r*l
+
+    END SUBROUTINE areaRad
+
+    !Computes element functions at point xii
+    PURE FUNCTION fPsiRad(xi) RESULT(fPsi)
+      IMPLICIT NONE
+
+      REAL(8), INTENT(in):: xi(1:3)
+      REAL(8), ALLOCATABLE:: fPsi(:)
+
+      ALLOCATE(fPsi(1:2))
+
+      fPsi(1) = 1.D0 - xi(1)
+      fPsi(2) = 1.D0 + xi(1)
+      fPsi    = fPsi * 5.D-1
+
+    END FUNCTION fPsiRad
+
+    !Computes element derivative shape function at Xii
+    PURE FUNCTION dPsiRad(xi) RESULT(dPsi)
+      IMPLICIT NONE
+
+      REAL(8), INTENT(in):: xi(1:3)
+      REAL(8), ALLOCATABLE:: dPsi(:,:)
+
+      ALLOCATE(dPsi(1:1, 1:2))
+
+      dPsi(1, 1) = -5.D-1
+      dPsi(1, 2) =  5.D-1
+
+    END FUNCTION dPsiRad
+
+    !Computes partial derivatives of coordinates
+    PURE SUBROUTINE partialDerRad(self, dPsi, dx)
+      IMPLICIT NONE
+
+      CLASS(meshVol1DRadSegm), INTENT(in):: self
+      REAL(8), INTENT(in):: dPsi(1:,1:)
+      REAL(8), INTENT(out), DIMENSION(1):: dx
+
+      dx(1) = DOT_PRODUCT(dPsi(1,:), self%r)
+
+    END SUBROUTINE partialDerRad
+
+    !Computes local stiffness matrix
+    PURE FUNCTION elemKRad(self) RESULT(ke)
+      USE moduleConstParam, ONLY: PI
+      IMPLICIT NONE
+
+      CLASS(meshVol1DRadSegm), INTENT(in):: self
+      REAL(8):: ke(1:2,1:2)
+      REAL(8):: Xii(1:3)
+      REAL(8):: dPsi(1:1, 1:2)
+      REAL(8):: invJ
+      REAL(8):: r, fPsi(1:2)
+
+      ke      = 0.D0
+      Xii     = 0.D0
+      dPsi    = self%dPsi(Xii)
+      invJ    = self%invJac(Xii, dPsi)
+      fPsi    = self%fPsi(Xii)
+      r       = DOT_PRODUCT(fPsi, self%r)
+      ke(1,:) = (/ dPsi(1,1)*dPsi(1,1), dPsi(1,1)*dPsi(1,2) /) 
+      ke(2,:) = (/ dPsi(1,2)*dPsi(1,1), dPsi(1,2)*dPsi(1,2) /) 
+      ke      = 2.D0*ke*invJ
+      ke      = ke*r*2.D0*PI
+
+    END FUNCTION elemKRad
+
+    PURE FUNCTION elemFRad(self, source) RESULT(localF)
+      USE moduleConstParam, ONLY: PI
+      IMPLICIT NONE
+
+      CLASS(meshVol1DRadSegm), INTENT(in):: self
+      REAL(8), INTENT(in):: source(1:)
+      REAL(8), ALLOCATABLE:: localF(:)
+      REAL(8):: fPsi(1:2)
+      REAL(8):: r
+      REAL(8):: detJ
+      REAL(8):: Xii(1:3)
+
+      Xii    = 0.D0
+      fPsi   = self%fPsi(Xii)
+      detJ   = self%detJac(Xii)
+      r      = DOT_PRODUCT(fPsi,self%r)
+      ALLOCATE(localF(1:2))
+      localF = 2.D0*DOT_PRODUCT(fPsi, source)*detJ
+      localF = localF*r*2.D0*PI
+
+    END FUNCTION elemFRad
+
+    PURE FUNCTION weightRad(xi) RESULT(w)
+      IMPLICIT NONE
+
+      REAL(8), INTENT(in):: xi(1:3)
+      REAL(8):: w(1:2)
+
+      w = fPsiRad(xi)
+
+    END FUNCTION weightRad
+
+    PURE FUNCTION insideRad(xi) RESULT(ins)
+      IMPLICIT NONE
+
+      REAL(8), INTENT(in):: xi(1:3)
+      LOGICAL:: ins
+
+      ins = xi(1) >=-1.D0 .AND. &
+            xi(1) <= 1.D0
+
+    END FUNCTION insideRad
+
+    SUBROUTINE scatterRad(self, part)
+      USE moduleOutput
+      USE moduleSpecies
+      IMPLICIT NONE
+
+      CLASS(meshVol1DRadSegm), INTENT(in):: self
+      CLASS(particle), INTENT(in):: part
+      TYPE(outputNode), POINTER:: vertex
+      REAL(8):: w_p(1:2)
+      REAL(8):: tensorS(1:3,1:3)
+
+      w_p = self%weight(part%xi)
+      tensorS = outerProduct(part%v, part%v)
+
+      vertex => self%n1%output(part%sp)
+      vertex%den          = vertex%den          + part%weight*w_p(1)
+      vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(1)*part%v(:)
+      vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(1)*tensorS
+
+      vertex => self%n2%output(part%sp)
+      vertex%den          = vertex%den          + part%weight*w_p(2)
+      vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(2)*part%v(:)
+      vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(2)*tensorS
+
+    END SUBROUTINE scatterRad
+
+    !Gathers EF at position Xii
+    PURE FUNCTION gatherEFRad(self, xi) RESULT(EF)
+      IMPLICIT NONE
+
+      CLASS(meshVol1DRadSegm), INTENT(in):: self
+      REAL(8), INTENT(in):: xi(1:3)
+      REAL(8):: dPsi(1, 1:2)
+      REAL(8):: phi(1:2)
+      REAL(8):: EF(1:3)
+      REAL(8):: invJ
+
+      phi = (/ self%n1%emData%phi, &
+               self%n2%emData%phi /)
+
+      dPsi  =  self%dPsi(xi)
+      invJ  =  self%invJac(xi, dPsi)
+      EF(1) = -DOT_PRODUCT(dPsi(1, :), phi)*invJ
+      EF(2) =  0.D0
+      EF(3) =  0.D0
+
+    END FUNCTION gatherEFRad
+
+    !Get nodes from 1D volume
+    PURE FUNCTION getNodesRad(self) RESULT(n)
+      IMPLICIT NONE
+
+      CLASS(meshVol1DRadSegm), INTENT(in):: self
+      INTEGER, ALLOCATABLE:: n(:)
+
+      ALLOCATE(n(1:2))
+      n = (/ self%n1%n, self%n2%n /)
+
+    END FUNCTION getNodesRad
+
+    PURE FUNCTION phy2logRad(self, r) RESULT(xN)
+      IMPLICIT NONE
+
+      CLASS(meshVol1DRadSegm), INTENT(in):: self
+      REAL(8), INTENT(in):: r(1:3)
+      REAL(8):: xN(1:3)
+
+      xN = 0.D0
+      xN(1) = 2.D0*(r(1) - self%r(1))/(self%r(2) - self%r(1)) - 1.D0
+
+    END FUNCTION phy2logRad
+
+    !Get next element for a logical position xi
+    SUBROUTINE nextElementRad(self, xi, nextElement)
+      IMPLICIT NONE
+
+      CLASS(meshVol1DRadSegm), INTENT(in):: self
+      REAL(8), INTENT(in):: xi(1:3)
+      CLASS(*), POINTER, INTENT(out):: nextElement
+      
+      NULLIFY(nextElement)
+      IF     (xi(1) < -1.D0) THEN
+        nextElement => self%e2
+
+      ELSEIF (xi(1) >  1.D0) THEN
+        nextElement => self%e1
+
+      END IF
+
+    END SUBROUTINE nextElementRad
+
+    !Reset the output of nodes in element
+    PURE SUBROUTINE resetOutputRad(self)
+      USE moduleSpecies
+      USE moduleOutput
+      IMPLICIT NONE
+
+      CLASS(meshVol1DRadSegm), INTENT(inout):: self
+      INTEGER:: k
+
+      DO k = 1, nSpecies
+        self%n1%output(k)%den     = 0.D0
+        self%n1%output(k)%mom     = 0.D0
+        self%n1%output(k)%tensorS = 0.D0
+
+        self%n2%output(k)%den     = 0.D0
+        self%n2%output(k)%mom     = 0.D0
+        self%n2%output(k)%tensorS = 0.D0
+
+      END DO
+
+    END SUBROUTINE resetOutputRad
+
+
+    !COMMON FUNCTIONS FOR 1D VOLUME ELEMENTS
+    !Computes the element Jacobian determinant
+    PURE FUNCTION detJ1DRad(self, xi, dPsi_in) RESULT(dJ)
+      IMPLICIT NONE
+
+      CLASS(meshVol1DRad), INTENT(in):: self
+      REAL(8), INTENT(in):: xi(1:3)
+      REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:,1:)
+      REAL(8), ALLOCATABLE:: dPsi(:,:)
+      REAL(8):: dJ
+      REAL(8):: dx(1)
+
+      IF (PRESENT(dPsi_in)) THEN
+        dPsi = dPsi_in
+
+      ELSE
+        dPsi = self%dPsi(xi)
+
+      END IF
+
+      CALL self%partialDer(dPsi, dx)
+      dJ = dx(1)
+
+    END FUNCTION detJ1DRad
+
+    !Computes the invers Jacobian
+    PURE FUNCTION invJ1DRad(self, xi, dPsi_in) RESULT(invJ)
+      IMPLICIT NONE
+
+      CLASS(meshVol1DRad), INTENT(in):: self
+      REAL(8), INTENT(in):: xi(1:3)
+      REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:,1:)
+      REAL(8), ALLOCATABLE:: dPsi(:,:)
+      REAL(8):: dx(1)
+      REAL(8):: invJ
+
+      IF (PRESENT(dPsi_in)) THEN
+        dPsi = dPsi_in
+
+      ELSE
+        dPsi = self%dPsi(xi)
+
+      END IF
+
+      CALL self%partialDer(dPsi, dx)
+      invJ = 1.D0/dx(1)
+
+    END FUNCTION invJ1DRad
+
+END MODULE moduleMesh1DRad
+

src/modules/mesh/1DRad/moduleMesh1DRadBoundary.f90

@@ -0,0 +1,40 @@
+MODULE moduleMesh1DRadBoundary
+  USE moduleMesh1DRad
+
+  TYPE, PUBLIC, EXTENDS(meshEdge1DRad):: meshEdge1DRadRef
+    CONTAINS
+      PROCEDURE, PASS:: fBoundary => reflection
+
+  END TYPE meshEdge1DRadRef
+
+  TYPE, PUBLIC, EXTENDS(meshEdge1DRad):: meshEdge1DRadAbs
+    CONTAINS
+      PROCEDURE, PASS:: fBoundary => absorption
+
+  END TYPE meshEdge1DRadAbs
+
+  CONTAINS
+    SUBROUTINE reflection(self, part)
+      USE moduleSpecies
+      IMPLICIT NONE
+
+      CLASS(meshEdge1DRadRef), INTENT(inout):: self
+      CLASS(particle), INTENT(inout):: part
+
+      part%v(1) = -part%v(1)
+      part%r(1) =  2.D0*self%r - part%r(1)
+
+    END SUBROUTINE reflection
+
+    SUBROUTINE absorption(self, part)
+      USE moduleSpecies
+      IMPLICIT NONE
+
+      CLASS(meshEdge1DRadAbs), INTENT(inout):: self
+      CLASS(particle), INTENT(inout):: part
+
+      part%n_in = .FALSE.
+
+    END SUBROUTINE absorption
+
+END MODULE moduleMesh1DRadBoundary

src/modules/mesh/1DRad/moduleMesh1DRadRead.f90

@@ -1,13 +1,13 @@
-MODULE moduleMesh1DRead
+MODULE moduleMesh1DRadRead
   USE moduleMesh
-  USE moduleMesh1D
+  USE moduleMesh1DRad
-  USE moduleMesh1DBoundary
+  USE moduleMesh1DRadBoundary
   
   !TODO: make this abstract to allow different mesh formats
-  TYPE, EXTENDS(meshGeneric):: mesh1DGeneric
+  TYPE, EXTENDS(meshGeneric):: mesh1DRadGeneric
     CONTAINS
-      PROCEDURE, PASS:: init     => init1DMesh
+      PROCEDURE, PASS:: init     => init1DRadMesh
-      PROCEDURE, PASS:: readMesh => readMesh1D
+      PROCEDURE, PASS:: readMesh => readMesh1DRad
 
   END TYPE
 
@@ -18,12 +18,12 @@ MODULE moduleMesh1DRead
 
   CONTAINS
     !Init 1D mesh
-    SUBROUTINE init1DMesh(self, meshFormat)
+    SUBROUTINE init1DRadMesh(self, meshFormat)
       USE moduleMesh
       USE moduleErrors
       IMPLICIT NONE
 
-      CLASS(mesh1DGeneric), INTENT(out):: self
+      CLASS(mesh1DRadGeneric), INTENT(out):: self
       CHARACTER(:), ALLOCATABLE, INTENT(in):: meshFormat
 
       SELECT CASE(meshFormat)
@@ -33,18 +33,18 @@ MODULE moduleMesh1DRead
         self%printEM     => printEMGmsh
 
       CASE DEFAULT
-        CALL criticalError("Mesh type " // meshFormat // " not supported.", "init1D")
+        CALL criticalError("Mesh type " // meshFormat // " not supported.", "init1DRad")
         
       END SELECT
 
-    END SUBROUTINE init1DMesh
+    END SUBROUTINE init1DRadMesh
 
     !Reads 1D mesh
-    SUBROUTINE readMesh1D(self, filename)
+    SUBROUTINE readMesh1DRad(self, filename)
       USE  moduleBoundary
       IMPLICIT NONE
 
-      CLASS(mesh1DGeneric), INTENT(inout):: self
+      CLASS(mesh1DRadGeneric), INTENT(inout):: self
       CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
       REAL(8):: x
       INTEGER:: p(1:2)
@@ -70,7 +70,7 @@ MODULE moduleMesh1DRead
       !Read nodes coordinates. Only relevant for x
       DO e = 1, self%numNodes
         READ(10, *) n, x
-        ALLOCATE(meshNode1D:: self%nodes(n)%obj)
+        ALLOCATE(meshNode1DRad:: self%nodes(n)%obj)
         CALL self%nodes(n)%obj%init(n, (/ x, 0.D0, 0.D0 /))
 
       END DO
@@ -109,10 +109,10 @@ MODULE moduleMesh1DRead
         bt = getBoundaryId(boundaryType)
         SELECT CASE(boundary(bt)%obj%boundaryType)
         CASE ('reflection')
-          ALLOCATE(meshEdge1DRef:: self%edges(e)%obj)
+          ALLOCATE(meshEdge1DRadRef:: self%edges(e)%obj)
 
         CASE ('absorption')
-          ALLOCATE(meshEdge1DAbs:: self%edges(e)%obj)
+          ALLOCATE(meshEdge1DRadAbs:: self%edges(e)%obj)
 
         END SELECT
 
@@ -123,7 +123,7 @@ MODULE moduleMesh1DRead
       !Read and initialize volumes
       DO e = 1, self%numVols
         READ(10, *) n, elemType, eTemp, eTemp, eTemp, p(1:2)
-        ALLOCATE(meshVol1DSegm:: self%vols(e)%obj)
+        ALLOCATE(meshVol1DRadSegm:: self%vols(e)%obj)
         CALL self%vols(e)%obj%init(n - self%numEdges, p(1:2))
 
       END DO
@@ -152,7 +152,7 @@ MODULE moduleMesh1DRead
 
       END DO
 
-    END SUBROUTINE readMesh1D
+    END SUBROUTINE readMesh1DRad
 
     SUBROUTINE connectedVolVol(elemA, elemB)
       IMPLICIT NONE
@@ -161,9 +161,9 @@ MODULE moduleMesh1DRead
       CLASS(meshVol), INTENT(inout):: elemB
 
         SELECT TYPE(elemA)
-        TYPE IS(meshVol1DSegm)
+        TYPE IS(meshVol1DRadSegm)
           SELECT TYPE(elemB)
-          TYPE IS(meshVol1DSegm)
+          TYPE IS(meshVol1DRadSegm)
             CALL connectedSegmSegm(elemA, elemB)
 
           END SELECT
@@ -175,8 +175,8 @@ MODULE moduleMesh1DRead
     SUBROUTINE connectedSegmSegm(elemA, elemB)
       IMPLICIT NONE
 
-      CLASS(meshVol1DSegm), INTENT(inout), TARGET:: elemA
+      CLASS(meshVol1DRadSegm), INTENT(inout), TARGET:: elemA
-      CLASS(meshVol1DSegm), INTENT(inout), TARGET:: elemB
+      CLASS(meshVol1DRadSegm), INTENT(inout), TARGET:: elemB
 
       IF (.NOT. ASSOCIATED(elemA%e1) .AND. &
           elemA%n2%n == elemB%n1%n) THEN
@@ -202,9 +202,9 @@ MODULE moduleMesh1DRead
       CLASS(meshEdge), INTENT(inout):: elemB
 
       SELECT TYPE(elemA)
-      TYPE IS (meshVol1DSegm)
+      TYPE IS (meshVol1DRadSegm)
         SELECT TYPE(elemB)
-        CLASS IS(meshEdge1D)
+        CLASS IS(meshEdge1DRad)
           CALL connectedSegmEdge(elemA, elemB)
 
         END SELECT
@@ -216,8 +216,8 @@ MODULE moduleMesh1DRead
     SUBROUTINE connectedSegmEdge(elemA, elemB)
       IMPLICIT NONE
 
-      CLASS(meshVol1DSegm), INTENT(inout), TARGET:: elemA
+      CLASS(meshVol1DRadSegm), INTENT(inout), TARGET:: elemA
-      CLASS(meshEdge1D),    INTENT(inout), TARGET:: elemB
+      CLASS(meshEdge1DRad),    INTENT(inout), TARGET:: elemB
 
       IF (.NOT. ASSOCIATED(elemA%e1) .AND. &
           elemA%n2%n == elemB%n1%n) THEN
@@ -247,7 +247,7 @@ MODULE moduleMesh1DRead
       INTEGER:: n(1:2)
 
       SELECT TYPE(elem)
-      TYPE IS(meshVol1DSegm)
+      TYPE IS(meshVol1DRadSegm)
         localK = elem%elemK()
         n = (/ elem%n1%n, elem%n2%n /)
 
@@ -265,4 +265,4 @@ MODULE moduleMesh1DRead
 
     END SUBROUTINE constructGlobalK
 
-END MODULE moduleMesh1DRead
+END MODULE moduleMesh1DRadRead

src/modules/mesh/2DCyl/moduleMeshCyl.f90

@@ -360,7 +360,7 @@ MODULE moduleMeshCyl
 
     !Computes element local stiffness matrix
     PURE FUNCTION elemKQuad(self) RESULT(ke)
-      USE moduleConstParam
+      USE moduleConstParam, ONLY: PI
       IMPLICIT NONE
 
       CLASS(meshVolCylQuad), INTENT(in):: self

src/modules/mesh/makefile

@@ -1,4 +1,4 @@
-all: moduleMesh.o 2DCyl.o 1DCart.o
+all: moduleMesh.o 2DCyl.o 1DRad.o 1DCart.o
 
 2DCyl.o:
 	$(MAKE) -C 2DCyl all
@@ -6,6 +6,9 @@ all: moduleMesh.o 2DCyl.o 1DCart.o
 1DCart.o:
 	$(MAKE) -C 1DCart all
 
+1DRad.o:
+	$(MAKE) -C 1DRad all
+
 moduleMesh.o: moduleMesh.f90
 	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
 

src/modules/moduleInput.f90

@@ -361,7 +361,8 @@ MODULE moduleInput
     SUBROUTINE readGeometry(config)
       USE moduleMesh
       USE moduleMeshCylRead,    ONLY: meshCylGeneric
-      USE moduleMesh1DRead, ONLY: mesh1DGeneric
+      USE moduleMesh1DCartRead, ONLY: mesh1DCartGeneric
+      USE moduleMesh1DRadRead,  ONLY: mesh1DRadGeneric
       USE moduleErrors
       USE moduleOutput
       USE json_module
@@ -381,7 +382,11 @@ MODULE moduleInput
 
       CASE ("1DCart")
         !Creates a 1D cartesian mesh
-        ALLOCATE(mesh1DGeneric:: mesh)
+        ALLOCATE(mesh1DCartGeneric:: mesh)
+
+      CASE ("1DRad")
+        !Creates a 1D cartesian mesh
+        ALLOCATE(mesh1DRadGeneric:: mesh)
 
       CASE DEFAULT
         CALL criticalError("Geometry type " // geometryType // " not supported.", "readGeometry")

src/modules/moduleSolver.f90

@@ -73,10 +73,13 @@ MODULE moduleSolver
         self%pushParticle => pushCylCharged
 
       CASE('1DCartCharged')
-        self%pushParticle => push1DCharged
+        self%pushParticle => push1DCartCharged
+
+      CASE('1DRadCharged')
+        self%pushParticle => push1DRadCharged
 
       CASE DEFAULT
-        CALL criticalError('Solver ' // pusherType // ' not found','readCase')
+        CALL criticalError('Pusher ' // pusherType // ' not found','initPusher')
 
       END SELECT
 
@@ -225,7 +228,7 @@ MODULE moduleSolver
     END SUBROUTINE pushCylCharged
 
     !Push charged particles in 1D cartesian coordinates
-    PURE SUBROUTINE push1DCharged(part)
+    PURE SUBROUTINE push1DCartCharged(part)
       USE moduleSPecies
       USE moduleEM
       IMPLICIT NONE
@@ -249,7 +252,47 @@ MODULE moduleSolver
 
       part = part_temp
 
-    END SUBROUTINE push1DCharged
+    END SUBROUTINE push1DCartCharged
+
+    !Push one particle. Boris pusher for 1D Radial Charged particle
+    PURE SUBROUTINE push1DRadCharged(part)
+      USE moduleSpecies
+      USE moduleEM
+      IMPLICIT NONE
+
+      TYPE(particle), INTENT(inout):: part
+      REAL(8):: v_p_oh_star(1:2)
+      TYPE(particle):: part_temp
+      REAL(8):: x_new, y_new, r, sin_alpha, cos_alpha
+      REAL(8):: tauSp
+      REAL(8):: qmEFt(1:3)!charge*tauSp*EF/mass
+
+      part_temp = part
+      !Time step for the species
+      tauSp = tau(part_temp%sp)
+      !Get electric field at particle position
+      qmEFt = part_temp%qm*gatherElecField(part_temp)*tauSp
+      !r,theta
+      v_p_oh_star(1) = part%v(1) + qmEFt(1)
+      x_new = part%r(1) + v_p_oh_star(1)*tauSp
+      v_p_oh_star(2) = part%v(2) + qmEFt(2)
+      y_new =             v_p_oh_star(2)*tauSp
+      r = DSQRT(x_new**2+y_new**2)
+      part_temp%r(1) = r
+      IF (r > 0.D0) THEN
+        sin_alpha = y_new/r
+        cos_alpha = x_new/r
+      ELSE
+        sin_alpha = 0.D0
+        cos_alpha = 1.D0
+      END IF
+      part_temp%v(1) =  cos_alpha*v_p_oh_star(1)+sin_alpha*v_p_oh_star(2)
+      part_temp%v(2) = -sin_alpha*v_p_oh_star(1)+cos_alpha*v_p_oh_star(2)
+      part_temp%n_in = .FALSE. !Assume particle is outside until cell is found
+      !Copy temporal particle to particle
+      part=part_temp
+
+    END SUBROUTINE push1DRadCharged
 
     !Do the collisions in all the cells
     SUBROUTINE doCollisions()