NOT WORKING

I need to fix something in another place and make some test, so this
commit does not comile.

CITATION.cff

@@ -1,50 +0,0 @@
-# This CITATION.cff file was generated with cffinit.
-# Visit https://bit.ly/cffinit to generate yours today!
-
-cff-version: 1.2.0
-title: Finite element Particle Kinetic Code
-message: >-
-  If you use this software, please cite it using the
-  metadata from this file.
-type: software
-authors:
-  - given-names: Jorge
-    family-names: Gonzalez
-    email: jorge.gonzalez@upm.es
-    affiliation: Universidad Politécnica de Madrid
-    orcid: 'https://orcid.org/0000-0001-7905-5001'
-repository-code: 'https://gitlab.com/JorgeGonz/fpakc'
-abstract: >-
-  Welcome to fpakc (Finite element PArticle Kinetic Code), a
-  modern object oriented Fortran open-source code for
-  particle simulations of plasma and gases. This code works
-  by simulating charged and neutral particles, following
-  their trajectories, collisions and boundary conditions
-  imposed by the user.
-
-  One of our aims is to make a code easy to maintain as well
-  as easy to use by a variety of reserchers and students.
-
-  This code is currenlty in very early steps of development.
-
-  The code aims to be easy to maintain and easy to use,
-  allowing its application from complex problems to easy
-  examples that can be used, for example, as teaching
-  exercises.
-
-  Parallelization techniques such as OpenMP, MPI will be
-  used to distribute the cpu load. We aim to make fpakc GPU
-  compatible in the future.
-
-  The codefpakc makes use of finite elements to generate
-  meshes in complex geometries. Particle properties are
-  deposited in the nodes and cells of the mesh. The
-  electromagnetic field, with the boundary conditions
-  imposed by the user, is solved also in this mesh.
-keywords:
-  - particle-in-cell
-  - plasma
-  - finite elements
-license: GPL-3.0
-version: beta
-date-released: '2025-10-01'

data/collisions/IO_e-Kr.dat

@@ -1,52 +0,0 @@
-# D108525 "refs": {"B56": {"note": "CLM-R294 (1989)"}}
-# Relative energy (eV) cross section (m^2)
-1.40E+01 0
-1.62E+01 7.249E-21
-1.88E+01 1.199E-20
-2.18E+01 1.644E-20
-2.53E+01 2.1E-20
-2.94E+01 2.542E-20
-3.41E+01 2.937E-20
-3.95E+01 3.26E-20
-4.58E+01 3.499E-20
-5.32E+01 3.653E-20
-6.17E+01 3.726E-20
-7.15E+01 3.728E-20
-8.29E+01 3.671E-20
-9.62E+01 3.566E-20
-1.12E+02 3.426E-20
-1.29E+02 3.259E-20
-1.50E+02 3.075E-20
-1.74E+02 2.881E-20
-2.02E+02 2.682E-20
-2.34E+02 2.484E-20
-2.72E+02 2.289E-20
-3.15E+02 2.101E-20
-3.65E+02 1.922E-20
-4.24E+02 1.751E-20
-4.91E+02 1.592E-20
-5.70E+02 1.443E-20
-6.61E+02 1.305E-20
-7.67E+02 1.177E-20
-8.89E+02 1.06E-20
-1.03E+03 9.526E-21
-1.20E+03 8.547E-21
-1.39E+03 7.658E-21
-1.61E+03 6.851E-21
-1.87E+03 6.121E-21
-2.16E+03 5.462E-21
-2.51E+03 4.868E-21
-2.91E+03 4.334E-21
-3.38E+03 3.855E-21
-3.92E+03 3.426E-21
-4.54E+03 3.041E-21
-5.27E+03 2.698E-21
-6.11E+03 2.391E-21
-7.09E+03 2.118E-21
-8.22E+03 1.875E-21
-9.53E+03 1.658E-21
-1.11E+04 1.466E-21
-1.28E+04 1.295E-21
-1.49E+04 1.143E-21
-1.72E+04 1.009E-21
-2.00E+04 8.898E-22

data/collisions/IO_e-Xe.dat

@@ -1,52 +0,0 @@
-# EL cross sections extracted from PROGRAM MAGBOLTZ, VERSION 7.1 JUNE 2004 www.lxcat.net/Biagi-v7.1
-# Relative energy (eV) cross section (m^2)
-1.21E+01 0
-1.41E+01 3.923E-21
-1.64E+01 1.194E-20
-1.91E+01 2.1E-20
-2.22E+01 2.946E-20
-2.58E+01 3.65E-20
-3.00E+01 4.185E-20
-3.49E+01 4.552E-20
-4.06E+01 4.766E-20
-4.72E+01 4.85E-20
-5.49E+01 4.828E-20
-6.39E+01 5.031E-20
-7.43E+01 5.1E-20
-8.64E+01 5.1E-20
-1.01E+02 5.032E-20
-1.17E+02 4.906E-20
-1.36E+02 4.732E-20
-1.58E+02 4.521E-20
-1.84E+02 4.283E-20
-2.14E+02 4.029E-20
-2.49E+02 3.764E-20
-2.90E+02 3.497E-20
-3.37E+02 3.233E-20
-3.92E+02 2.975E-20
-4.56E+02 2.726E-20
-5.31E+02 2.489E-20
-6.17E+02 2.266E-20
-7.18E+02 2.056E-20
-8.35E+02 1.861E-20
-9.72E+02 1.68E-20
-1.13E+03 1.514E-20
-1.32E+03 1.361E-20
-1.53E+03 1.221E-20
-1.78E+03 1.094E-20
-2.07E+03 9.781E-21
-2.41E+03 8.735E-21
-2.80E+03 7.789E-21
-3.26E+03 6.938E-21
-3.79E+03 6.171E-21
-4.41E+03 5.484E-21
-5.13E+03 4.868E-21
-5.97E+03 4.316E-21
-6.94E+03 3.824E-21
-8.07E+03 3.385E-21
-9.39E+03 2.994E-21
-1.09E+04 2.646E-21
-1.27E+04 2.336E-21
-1.48E+04 2.062E-21
-1.72E+04 1.818E-21
-2.00E+04 1.602E-21

data/see/constant_yield.dat

@@ -0,0 +1,4 @@
+#Relative energy (eV) yield ()
+0.000 1.000E-01
+1.000 1.000E-01
+

doc/user-manual/.gitignore

@@ -6,7 +6,6 @@
 *.aux
 *.ps
 bibliography.bib.bak
-bibliography.bib.sav
 *.bbl
 *.blg
 *.out
@@ -20,3 +19,4 @@ fpakc_UserManual-blx.bib
 *.gls
 *.ist
 fpakc_UserManual.run.xml
+*.bib.sav

doc/user-manual/figures/logos/PPartiC.eps

@@ -0,0 +1,177 @@
+%!PS-Adobe-3.0 EPSF-3.0
+%%Creator: cairo 1.16.0 (https://cairographics.org)
+%%CreationDate: Wed Oct 07 12:01:13 2020
+%%Pages: 1
+%%DocumentData: Clean7Bit
+%%LanguageLevel: 2
+%%BoundingBox: 6 7 205 57
+%%EndComments
+%%BeginProlog
+50 dict begin
+/q { gsave } bind def
+/Q { grestore } bind def
+/cm { 6 array astore concat } bind def
+/w { setlinewidth } bind def
+/J { setlinecap } bind def
+/j { setlinejoin } bind def
+/M { setmiterlimit } bind def
+/d { setdash } bind def
+/m { moveto } bind def
+/l { lineto } bind def
+/c { curveto } bind def
+/h { closepath } bind def
+/re { exch dup neg 3 1 roll 5 3 roll moveto 0 rlineto
+      0 exch rlineto 0 rlineto closepath } bind def
+/S { stroke } bind def
+/f { fill } bind def
+/f* { eofill } bind def
+/n { newpath } bind def
+/W { clip } bind def
+/W* { eoclip } bind def
+/BT { } bind def
+/ET { } bind def
+/BDC { mark 3 1 roll /BDC pdfmark } bind def
+/EMC { mark /EMC pdfmark } bind def
+/cairo_store_point { /cairo_point_y exch def /cairo_point_x exch def } def
+/Tj { show currentpoint cairo_store_point } bind def
+/TJ {
+  {
+    dup
+    type /stringtype eq
+    { show } { -0.001 mul 0 cairo_font_matrix dtransform rmoveto } ifelse
+  } forall
+  currentpoint cairo_store_point
+} bind def
+/cairo_selectfont { cairo_font_matrix aload pop pop pop 0 0 6 array astore
+    cairo_font exch selectfont cairo_point_x cairo_point_y moveto } bind def
+/Tf { pop /cairo_font exch def /cairo_font_matrix where
+      { pop cairo_selectfont } if } bind def
+/Td { matrix translate cairo_font_matrix matrix concatmatrix dup
+      /cairo_font_matrix exch def dup 4 get exch 5 get cairo_store_point
+      /cairo_font where { pop cairo_selectfont } if } bind def
+/Tm { 2 copy 8 2 roll 6 array astore /cairo_font_matrix exch def
+      cairo_store_point /cairo_font where { pop cairo_selectfont } if } bind def
+/g { setgray } bind def
+/rg { setrgbcolor } bind def
+/d1 { setcachedevice } bind def
+/cairo_data_source {
+  CairoDataIndex CairoData length lt
+    { CairoData CairoDataIndex get /CairoDataIndex CairoDataIndex 1 add def }
+    { () } ifelse
+} def
+/cairo_flush_ascii85_file { cairo_ascii85_file status { cairo_ascii85_file flushfile } if } def
+/cairo_image { image cairo_flush_ascii85_file } def
+/cairo_imagemask { imagemask cairo_flush_ascii85_file } def
+%%EndProlog
+%%BeginSetup
+%%EndSetup
+%%Page: 1 1
+%%BeginPageSetup
+%%PageBoundingBox: 6 7 205 57
+%%EndPageSetup
+q 6 7 199 50 rectclip
+1 0 0 -1 0 71 cm q
+0 g
+49.766 27.305 m 49.766 29.898 49.078 31.813 47.703 33.039 c 46.328 34.27
+ 44.297 34.883 41.609 34.883 c 36.625 34.883 l 36.625 47.945 l 32.438 47.945
+ l 32.438 14.352 l 41.516 14.352 l 44.234 14.352 46.285 14.953 47.672 16.148
+ c 49.066 17.348 49.766 19.242 49.766 21.836 c h
+45.594 21.398 m 45.594 20.18 45.313 19.301 44.75 18.758 c 44.195 18.219
+ 43.313 17.945 42.094 17.945 c 36.625 17.945 l 36.625 31.289 l 42.094 31.289
+ l 43.313 31.289 44.195 31.008 44.75 30.445 c 45.313 29.883 45.594 29 45.594
+ 27.789 c h
+45.594 21.398 m f
+75.641 27.305 m 75.641 29.898 74.953 31.813 73.578 33.039 c 72.203 34.27
+ 70.172 34.883 67.484 34.883 c 62.5 34.883 l 62.5 47.945 l 58.313 47.945
+ l 58.313 14.352 l 67.391 14.352 l 70.109 14.352 72.16 14.953 73.547 16.148
+ c 74.941 17.348 75.641 19.242 75.641 21.836 c h
+71.469 21.398 m 71.469 20.18 71.188 19.301 70.625 18.758 c 70.07 18.219
+ 69.188 17.945 67.969 17.945 c 62.5 17.945 l 62.5 31.289 l 67.969 31.289
+ l 69.188 31.289 70.07 31.008 70.625 30.445 c 71.188 29.883 71.469 29 71.469
+ 27.789 c h
+71.469 21.398 m f
+99.891 47.945 m 99.473 47.945 99.016 47.82 98.516 47.57 c 98.023 47.313
+ 97.602 46.957 97.25 46.508 c 96.508 47.469 95.535 47.945 94.328 47.945 
+c 90.047 47.945 l 87.617 47.945 85.852 47.461 84.75 46.492 c 83.645 45.516
+ 83.094 43.836 83.094 41.461 c 83.094 40.508 l 83.094 35.551 85.426 33.07
+ 90.094 33.07 c 96.047 33.07 l 96.047 30.664 l 96.047 29.676 95.742 28.91
+ 95.141 28.367 c 94.535 27.816 93.641 27.539 92.453 27.539 c 85.016 27.539
+ l 85.016 23.945 l 91.969 23.945 l 94.75 23.945 96.801 24.555 98.125 25.773
+ c 99.457 26.984 100.125 28.938 100.125 31.633 c 100.125 43.148 l 100.125
+ 43.629 100.305 43.988 100.672 44.227 c 101.047 44.469 101.664 44.586 102.531
+ 44.586 c 102.531 47.945 l h
+93.844 44.352 m 94.801 44.352 95.406 44.129 95.656 43.68 c 95.914 43.234
+ 96.047 42.656 96.047 41.945 c 96.047 36.664 l 90.047 36.664 l 89.211 36.664
+ 88.523 36.941 87.984 37.492 c 87.441 38.035 87.172 38.719 87.172 39.539
+ c 87.172 41.945 l 87.172 42.781 87.363 43.391 87.75 43.773 c 88.133 44.16
+ 88.742 44.352 89.578 44.352 c h
+93.844 44.352 m f
+114.047 44.352 m 114.047 27.539 l 110.547 27.539 l 110.547 23.945 l 117.891
+ 23.945 l 117.891 25.242 l 118.953 24.379 120.234 23.945 121.734 23.945 
+c 126.672 23.945 l 126.672 27.789 l 121.25 27.789 l 120.227 27.789 119.461
+ 28.117 118.953 28.773 c 118.441 29.43 118.172 30.094 118.141 30.758 c 118.141
+ 44.352 l 123.359 44.352 l 123.359 47.945 l 110.547 47.945 l 110.547 44.352
+ l h
+114.047 44.352 m f
+146.359 47.945 m 144.723 47.945 143.473 47.422 142.609 46.367 c 141.754
+ 45.305 141.328 43.91 141.328 42.18 c 141.328 27.539 l 136.766 27.539 l 
+136.766 23.945 l 141.328 23.945 l 141.328 18.18 l 145.406 18.18 l 145.406
+ 23.945 l 152.406 23.945 l 152.406 27.539 l 145.406 27.539 l 145.406 42.43
+ l 145.406 43.129 145.578 43.625 145.922 43.914 c 146.273 44.207 146.82 
+44.352 147.563 44.352 c 152.406 44.352 l 152.406 47.945 l h
+146.359 47.945 m f
+172.719 14.352 m 172.719 18.664 l 168.156 18.664 l 168.156 14.352 l h
+161.344 27.539 m 161.344 23.945 l 172.469 23.945 l 172.469 44.352 l 179.578
+ 44.352 l 179.578 47.945 l 168.391 47.945 l 168.391 27.539 l h
+161.344 27.539 m f
+197.25 47.945 m 194.469 47.945 192.332 47.27 190.844 45.914 c 189.352 44.551
+ 188.609 42.508 188.609 39.789 c 188.609 22.508 l 188.609 19.789 189.352
+ 17.754 190.844 16.398 c 192.332 15.035 194.469 14.352 197.25 14.352 c 204.063
+ 14.352 l 204.063 17.992 l 196.766 17.992 l 195.578 17.992 194.613 18.348
+ 193.875 19.055 c 193.145 19.754 192.781 20.598 192.781 21.586 c 192.781
+ 40.695 l 192.781 41.688 193.145 42.535 193.875 43.242 c 194.613 43.953 
+195.578 44.305 196.766 44.305 c 204.063 44.305 l 204.063 47.945 l h
+197.25 47.945 m f
+0.831373 0 0 rg
+17.277 53.707 m 17.277 56.629 14.91 58.996 11.988 58.996 c 9.066 58.996
+ 6.695 56.629 6.695 53.707 c 6.695 50.785 9.066 48.414 11.988 48.414 c 14.91
+ 48.414 17.277 50.785 17.277 53.707 c f
+0 0 0.501961 rg
+23.906 61.84 m 23.906 62.969 22.992 63.883 21.863 63.883 c 20.734 63.883
+ 19.82 62.969 19.82 61.84 c 19.82 60.715 20.734 59.801 21.863 59.801 c 22.992
+ 59.801 23.906 60.715 23.906 61.84 c f
+0.596078 g
+0.751181 w
+0 J
+0 j
+[] 0.0 d
+4 M q 1 0 0 1 0 0 cm
+16.941 50.758 m 29.465 47.813 l S Q
+27.848 47.258 m 30.453 47.578 l 28.262 49.023 l 28.523 48.41 28.355 47.699
+ 27.848 47.258 c h
+27.848 47.258 m f*
+0.137102 w
+1 j
+q -1 0.235294 -0.235294 -1 0 0 cm
+-15.851 -50.987 m -18.248 -51.872 l -15.849 -52.753 l -16.234 -52.23 -16.233
+ -51.519 -15.851 -50.987 c h
+-15.851 -50.987 m S Q
+0.6 g
+0.751178 w
+0 j
+q 1 0 0 1 0 0 cm
+23.234 59.613 m 30.398 49.301 l S Q
+28.824 49.969 m 30.973 48.461 l 30.313 51 l 30.098 50.375 29.496 49.957
+ 28.824 49.969 c h
+28.824 49.969 m f*
+0.115667 w
+1 j
+q -0.694805 1 -1 -0.694805 0 0 cm
+20.193 -42.855 m 18.17 -43.597 l 20.191 -44.342 l 19.87 -43.904 19.87 -43.302
+ 20.193 -42.855 c h
+20.193 -42.855 m S Q
+Q Q
+showpage
+%%Trailer
+end
+%%EOF

doc/user-manual/figures/logos/PPartiC.svg

@@ -0,0 +1,159 @@
+<?xml version="1.0" encoding="UTF-8" standalone="no"?>
+<svg
+   xmlns:dc="http://purl.org/dc/elements/1.1/"
+   xmlns:cc="http://creativecommons.org/ns#"
+   xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#"
+   xmlns:svg="http://www.w3.org/2000/svg"
+   xmlns="http://www.w3.org/2000/svg"
+   xmlns:sodipodi="http://sodipodi.sourceforge.net/DTD/sodipodi-0.dtd"
+   xmlns:inkscape="http://www.inkscape.org/namespaces/inkscape"
+   width="80mm"
+   height="25mm"
+   viewBox="0 0 80 25"
+   version="1.1"
+   id="svg8"
+   inkscape:version="1.0 (4035a4fb49, 2020-05-01)"
+   sodipodi:docname="PPartiC.svg">
+  <defs
+     id="defs2">
+    <marker
+       inkscape:stockid="Arrow2Mend"
+       orient="auto"
+       refY="0.0"
+       refX="0.0"
+       id="Arrow2Mend"
+       style="overflow:visible;"
+       inkscape:isstock="true">
+      <path
+         id="path881"
+         style="fill-rule:evenodd;stroke-width:0.625;stroke-linejoin:round;stroke:#000000;stroke-opacity:1;fill:#000000;fill-opacity:1"
+         d="M 8.7185878,4.0337352 L -2.2072895,0.016013256 L 8.7185884,-4.0017078 C 6.9730900,-1.6296469 6.9831476,1.6157441 8.7185878,4.0337352 z "
+         transform="scale(0.6) rotate(180) translate(0,0)" />
+    </marker>
+    <marker
+       inkscape:stockid="Arrow2Send"
+       orient="auto"
+       refY="0.0"
+       refX="0.0"
+       id="Arrow2Send"
+       style="overflow:visible;"
+       inkscape:isstock="true">
+      <path
+         id="path887"
+         style="fill-rule:evenodd;stroke-width:0.625;stroke-linejoin:round;stroke:#989898;stroke-opacity:1;fill:#989898;fill-opacity:1"
+         d="M 8.7185878,4.0337352 L -2.2072895,0.016013256 L 8.7185884,-4.0017078 C 6.9730900,-1.6296469 6.9831476,1.6157441 8.7185878,4.0337352 z "
+         transform="scale(0.3) rotate(180) translate(-2.3,0)" />
+    </marker>
+    <marker
+       inkscape:stockid="Arrow1Mend"
+       orient="auto"
+       refY="0.0"
+       refX="0.0"
+       id="Arrow1Mend"
+       style="overflow:visible;"
+       inkscape:isstock="true">
+      <path
+         id="path863"
+         d="M 0.0,0.0 L 5.0,-5.0 L -12.5,0.0 L 5.0,5.0 L 0.0,0.0 z "
+         style="fill-rule:evenodd;stroke:#000000;stroke-width:1pt;stroke-opacity:1;fill:#000000;fill-opacity:1"
+         transform="scale(0.4) rotate(180) translate(10,0)" />
+    </marker>
+    <marker
+       inkscape:stockid="Arrow1Lend"
+       orient="auto"
+       refY="0.0"
+       refX="0.0"
+       id="Arrow1Lend"
+       style="overflow:visible;"
+       inkscape:isstock="true">
+      <path
+         id="path857"
+         d="M 0.0,0.0 L 5.0,-5.0 L -12.5,0.0 L 5.0,5.0 L 0.0,0.0 z "
+         style="fill-rule:evenodd;stroke:#000000;stroke-width:1pt;stroke-opacity:1;fill:#000000;fill-opacity:1"
+         transform="scale(0.8) rotate(180) translate(12.5,0)" />
+    </marker>
+    <marker
+       inkscape:isstock="true"
+       style="overflow:visible"
+       id="Arrow2Send-8"
+       refX="0"
+       refY="0"
+       orient="auto"
+       inkscape:stockid="Arrow2Send">
+      <path
+         transform="matrix(-0.3,0,0,-0.3,0.69,0)"
+         d="M 8.7185878,4.0337352 -2.2072895,0.01601326 8.7185884,-4.0017078 c -1.7454984,2.3720609 -1.7354408,5.6174519 -6e-7,8.035443 z"
+         style="fill:#999999;fill-opacity:1;fill-rule:evenodd;stroke:#999999;stroke-width:0.625;stroke-linejoin:round;stroke-opacity:1"
+         id="path887-3" />
+    </marker>
+  </defs>
+  <sodipodi:namedview
+     id="base"
+     pagecolor="#ffffff"
+     bordercolor="#666666"
+     borderopacity="1.0"
+     inkscape:pageopacity="0.0"
+     inkscape:pageshadow="2"
+     inkscape:zoom="2.8"
+     inkscape:cx="-108.81501"
+     inkscape:cy="51.644646"
+     inkscape:document-units="mm"
+     inkscape:current-layer="layer1"
+     inkscape:document-rotation="0"
+     showgrid="false"
+     inkscape:window-width="2560"
+     inkscape:window-height="1378"
+     inkscape:window-x="-9"
+     inkscape:window-y="26"
+     inkscape:window-maximized="1" />
+  <metadata
+     id="metadata5">
+    <rdf:RDF>
+      <cc:Work
+         rdf:about="">
+        <dc:format>image/svg+xml</dc:format>
+        <dc:type
+           rdf:resource="http://purl.org/dc/dcmitype/StillImage" />
+        <dc:title></dc:title>
+      </cc:Work>
+    </rdf:RDF>
+  </metadata>
+  <g
+     inkscape:label="Layer 1"
+     inkscape:groupmode="layer"
+     id="layer1">
+    <text
+       xml:space="preserve"
+       style="font-style:normal;font-weight:normal;font-size:16.9333px;line-height:1.25;font-family:sans-serif;fill:#000000;fill-opacity:1;stroke:none;stroke-width:0.264583"
+       x="9.921875"
+       y="16.914433"
+       id="text12"><tspan
+         sodipodi:role="line"
+         id="tspan10"
+         x="9.921875"
+         y="16.914433"
+         style="font-style:italic;font-variant:normal;font-weight:normal;font-stretch:normal;font-size:16.9333px;font-family:'Share Tech Mono';-inkscape-font-specification:'Share Tech Mono Italic';stroke-width:0.264583">PPartiC</tspan></text>
+    <circle
+       style="fill:#d40000;stroke-width:0.264583"
+       id="path14"
+       cx="4.2286086"
+       cy="18.946056"
+       r="1.8662573" />
+    <circle
+       style="fill:#000080;stroke-width:0.576413"
+       id="path16"
+       cx="7.7130766"
+       cy="21.816313"
+       r="0.7205171" />
+    <path
+       style="fill:none;fill-opacity:1;stroke:#989898;stroke-width:0.265;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1;marker-end:url(#Arrow2Send)"
+       d="M 5.9767484,17.906622 10.394345,16.867188"
+       id="path18"
+       sodipodi:nodetypes="cc" />
+    <path
+       sodipodi:nodetypes="cc"
+       id="path18-6"
+       d="m 8.195956,21.029626 2.527716,-3.638021"
+       style="fill:none;stroke:#999999;stroke-width:0.264999;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1;marker-end:url(#Arrow2Send-8)" />
+  </g>
+</svg>

doc/user-manual/figures/logos/fpakc.eps

@@ -0,0 +1,294 @@
+%!PS-Adobe-3.0 EPSF-3.0
+%%Creator: cairo 1.16.0 (https://cairographics.org)
+%%CreationDate: Sat Nov 21 21:23:19 2020
+%%Pages: 1
+%%DocumentData: Clean7Bit
+%%LanguageLevel: 2
+%%BoundingBox: 8 15 191 78
+%%EndComments
+%%BeginProlog
+50 dict begin
+/q { gsave } bind def
+/Q { grestore } bind def
+/cm { 6 array astore concat } bind def
+/w { setlinewidth } bind def
+/J { setlinecap } bind def
+/j { setlinejoin } bind def
+/M { setmiterlimit } bind def
+/d { setdash } bind def
+/m { moveto } bind def
+/l { lineto } bind def
+/c { curveto } bind def
+/h { closepath } bind def
+/re { exch dup neg 3 1 roll 5 3 roll moveto 0 rlineto
+      0 exch rlineto 0 rlineto closepath } bind def
+/S { stroke } bind def
+/f { fill } bind def
+/f* { eofill } bind def
+/n { newpath } bind def
+/W { clip } bind def
+/W* { eoclip } bind def
+/BT { } bind def
+/ET { } bind def
+/BDC { mark 3 1 roll /BDC pdfmark } bind def
+/EMC { mark /EMC pdfmark } bind def
+/cairo_store_point { /cairo_point_y exch def /cairo_point_x exch def } def
+/Tj { show currentpoint cairo_store_point } bind def
+/TJ {
+  {
+    dup
+    type /stringtype eq
+    { show } { -0.001 mul 0 cairo_font_matrix dtransform rmoveto } ifelse
+  } forall
+  currentpoint cairo_store_point
+} bind def
+/cairo_selectfont { cairo_font_matrix aload pop pop pop 0 0 6 array astore
+    cairo_font exch selectfont cairo_point_x cairo_point_y moveto } bind def
+/Tf { pop /cairo_font exch def /cairo_font_matrix where
+      { pop cairo_selectfont } if } bind def
+/Td { matrix translate cairo_font_matrix matrix concatmatrix dup
+      /cairo_font_matrix exch def dup 4 get exch 5 get cairo_store_point
+      /cairo_font where { pop cairo_selectfont } if } bind def
+/Tm { 2 copy 8 2 roll 6 array astore /cairo_font_matrix exch def
+      cairo_store_point /cairo_font where { pop cairo_selectfont } if } bind def
+/g { setgray } bind def
+/rg { setrgbcolor } bind def
+/d1 { setcachedevice } bind def
+/cairo_data_source {
+  CairoDataIndex CairoData length lt
+    { CairoData CairoDataIndex get /CairoDataIndex CairoDataIndex 1 add def }
+    { () } ifelse
+} def
+/cairo_flush_ascii85_file { cairo_ascii85_file status { cairo_ascii85_file flushfile } if } def
+/cairo_image { image cairo_flush_ascii85_file } def
+/cairo_imagemask { imagemask cairo_flush_ascii85_file } def
+%%EndProlog
+%%BeginSetup
+%%BeginResource: font f-0-0
+%!FontType1-1.1 f-0-0 1.0
+11 dict begin
+/FontName /f-0-0 def
+/PaintType 0 def
+/FontType 1 def
+/FontMatrix [0.001 0 0 0.001 0 0] readonly def
+/FontBBox {30 -240 640 735                                   } readonly def
+/Encoding 256 array
+0 1 255 {1 index exch /.notdef put} for
+dup 97 /a put
+dup 99 /c put
+dup 102 /f put
+dup 107 /k put
+dup 112 /p put
+readonly def
+currentdict end
+currentfile eexec
+f983ef0097ece636fb4a96c74d26ab84185f6dfa4a16a7a1c27bbe3f1156aea698df336d20b467
+b10e7f33846656653c5ac6962759d3056cbdb3190bac614b984bf5a132dc418192443014ba63de
+800d392b6fea026574bb2535fd7bb5338f35bf15a88ea328fdaa49670c7852e3d060f3c5d6b07f
+2ef6d0f22646c5d18e19a2ae3ee120390f6dd96f76dcf1e127de5e9299077a00c17c0d71e36e5b
+9d5ec58fceda57739a6a4214d4b79d6c48d2784b60c320323c7acddddf34db833cac0cf109f799
+69d114a330d372e5c978a66acc84e3fe5557f6240856a013ffaa0199444e5c5036f775eba4a5c5
+8cde66cf604b9aca2178431127b8a1ff7ed633a65c04600af5f573483112251caad907dcd8c61e
+23500065b1568be79a17d2379b63d656e8f7715cdfdf357f0e30d9ab91d113c88231d875d60897
+1dee27eb5da34a09d2f1b1a2e7ab8be1036393cd6ae53eff0d77f452e9bf45eccc9e4836ae77c5
+b74262fa43bfccfd9a147a18dcdae6e50cc518129a64f669a5fae69c8082dec571d7d01a4e9f05
+6d3e52de0ea004acdbd1b3bf9b19aa17f839a75365a77b7275442a967093ffdf1694a72f9978a2
+304ac331d401b48e3e62a3a92cd39516dd480f088980d1ad8993a1f6fefb07e0d86b6f0293bb41
+68ac465726267cacb7516a0e910fe67c2dbfef06d8b64a9811506650d32fa182a0adcab8e2e21e
+ca6d0dc81959c25ea2d3f7ccec13e0cb4a7ef88e97c36e74fa13010220d6835ebdcbabdb507d84
+239e5483e8a8b7a52d6e1ea4ea1f5e6bef4534710c4055265aaa86fb445f3b2fc62cfdd9e283d8
+8bd083d09f0971cde00f2031b58b304d5f647f02aabf7ba9062c33979cd391f692c72ee179b7a8
+16f9c9e668d20021bd2b6a0f0114898729c6228be2895a696aaef0ebbcc842e64d5e72cc1d9b75
+44314028987a238f8fc4c18a0db3546c9ea42194b6bbdc45587e36d605fe2b7608d9292ddc0c9b
+be3e420b36fd52f9aef97f13533e101f34d4f882848f4845a7a824da815a710abe11a1ec8363c1
+06daa18dffb5a451af7bf3b20d79a63c94e305050ea893b7a61cc8cfdb2e8a593a073c2021b298
+40863c70742ab1734e1d6811cf1927832da10f562f6895575b50044179588ddd9ec4c413f68c3e
+3063f9594dc94115af4d9d4d6259c2ebb5afa796131772de3d297a8cc04f7f10398acc9142b1aa
+2da9741ad314918ff1553dafe4751b4c0efdc9d6a549acbf1b3d209f6ebe8f6561d627f37bbce1
+7213b92bf332c27718ca9f868f1724cda0774ea4c3a5a2ba99509eb9128c456e5526f234dc3adc
+37ac61ead9dafba1b5d58a9443ceb92474535cd3515e9ce357420b230fe927e81f06b2363c70aa
+b6e00858a44972ad3f8759069235bba0b8ae2c65a59fe3ee5642f88a8550a765907eb4f9432ac4
+9e896114d0bc969bc2c14acd9a50c31e2095133b6b4fc11a1136dbba4b515eaabf0cba23ffe795
+3532a1fca89780a841f3a5fe2514d31fd6d41fbcb5b8caadd53c1fa7b06506963f37006269d0ae
+fc1d5d6bd7f6788544e01a77bdf35aaacdd41dd4fc16237759516c60ee9b57e7b56606e0fb5a25
+9e6cb3b2e22d3ac4db73c228fda1b327cace7cbe25594ea2a9445efeda7826604e3daf18ca977a
+9a2788cddcea95b5460b648c12bebc3c39302a07481fe18ceed4c9e12ac7d51f8104bd589cf9f6
+68371d3aa8510b5ff08c972b84fd4e48c356f9098b9d130085f59bc7e748ac59e9715060d29b38
+d828e479d8a85e0aaaa7b4ab34b547643b9c3417305b29be240ba9514cd4656079f26f378a6af9
+f7c694b1332a14baa203faa5c6745eea7c9cc0b3fc13f722ff941857d141db5665b32ada8b7831
+06f5c01361bd4b161450546ad8812d364b394658c4d3e1e2cf1daba58315787fbe299a27a387f2
+1545b5d4bc3022e6a6803e971698504ebc44c9f7c02256a5ae01c9fb3194baf05fb7c9cd862427
+95ca9cfd57ba76d6697fc842ee19ff35ce8b9d3299a185a4d28537058b5ea471a818f521183ffa
+4b171f078577707cfbb3b2c6de4659b6dbacb952a0213b0586bd5879096ec53540ad07463e5b25
+1a38d715de6a0bb6f92d3e20e1b25be654fe7aad2e057428b68e2d6830f2ee0ca384fb862813ba
+d48832ec89991fa5e6bf9d3135378f0000000000000000000000000000000000000000000000000000000000000000
+0000000000000000000000000000000000000000000000000000000000000000
+0000000000000000000000000000000000000000000000000000000000000000
+0000000000000000000000000000000000000000000000000000000000000000
+0000000000000000000000000000000000000000000000000000000000000000
+0000000000000000000000000000000000000000000000000000000000000000
+0000000000000000000000000000000000000000000000000000000000000000
+0000000000000000000000000000000000000000000000000000000000000000
+cleartomark
+%%EndResource
+%%EndSetup
+%%Page: 1 1
+%%BeginPageSetup
+%%PageBoundingBox: 8 15 191 78
+%%EndPageSetup
+q 8 15 183 63 rectclip
+1 0 0 -1 0 86 cm q
+0 g
+BT
+56.000126 0 0 -56.000126 52.419201 52.791732 Tm
+/f-0-0 1 Tf
+[(fpak)10(c)]TJ
+ET
+0.749999 w
+0 J
+0 j
+[] 0.0 d
+4 M q 1 0 0 1 0 0 cm
+8.996 27.238 m 20.621 9.359 l 34.656 9.363 l 34.621 14.746 l 25.941 18.195
+ l 17.332 33.629 l 17.41 55.852 l 9.316 55.785 l h
+8.996 27.238 m S Q
+q 1 0 0 1 0 0 cm
+21.156 69.262 m 21.156 38.387 l 40.305 38.387 l 40.305 53.027 l 27.457 
+53.027 l 27.457 69.52 l h
+21.156 69.262 m S Q
+q 1 0 0 1 0 0 cm
+27.898 43.988 m 27.898 47.73 l 33.426 47.73 l 33.426 43.898 l h
+27.898 43.988 m S Q
+0.283465 w
+q 1 0 0 1 0 0 cm
+9.094 50.535 m 17.414 50.297 l S Q
+q 1 0 0 1 0 0 cm
+9.34 43.07 m 17.094 44.93 l S Q
+q 1 0 0 1 0 0 cm
+9.43 34.727 m 17.348 41.363 l S Q
+q 1 0 0 1 0 0 cm
+17.332 33.629 m 8.996 27.238 l S Q
+q 1 0 0 1 0 0 cm
+13.297 31.055 m 11.32 36.062 l 12.719 43.309 l 14.023 50.113 l 12.02 55.516
+ l S Q
+q 1 0 0 1 0 0 cm
+17.332 33.629 m 12.203 21.801 l S Q
+q 1 0 0 1 0 0 cm
+21.246 26.223 m 12.203 21.801 l S Q
+q 1 0 0 1 0 0 cm
+15.379 16.828 m 21.246 26.223 l S Q
+q 1 0 0 1 0 0 cm
+23.562 22.152 m 15.379 16.828 l S Q
+q 1 0 0 1 0 0 cm
+25.941 18.195 m 15.379 16.828 l S Q
+q 1 0 0 1 0 0 cm
+13.297 31.055 m 14.234 26.23 l S Q
+q 1 0 0 1 0 0 cm
+14.234 26.23 m 21.246 26.223 l S Q
+q 1 0 0 1 0 0 cm
+25.941 18.195 m 17.918 13.906 l S Q
+q 1 0 0 1 0 0 cm
+20.621 9.359 m 25.941 18.195 l S Q
+q 1 0 0 1 0 0 cm
+34.656 9.363 m 25.941 18.195 l S Q
+q 1 0 0 1 0 0 cm
+25.941 18.195 m 26.527 9.582 l S Q
+1 0 0 rg
+23.523 67.027 m 23.523 67.027 23.523 67.031 23.523 67.031 c 23.52 67.031
+ 23.52 67.027 23.52 67.027 c 23.52 67.027 23.52 67.027 23.523 67.027 c h
+23.523 67.027 m f
+0 g
+0.749999 w
+q 1 0 0 1 0 0 cm
+23.523 67.027 m 23.523 67.027 23.523 67.031 23.523 67.031 c 23.52 67.031
+ 23.52 67.027 23.52 67.027 c 23.52 67.027 23.52 67.027 23.523 67.027 c h
+23.523 67.027 m S Q
+1 0 0 rg
+24.66 63.297 m 24.66 63.297 24.656 63.297 24.656 63.297 c 24.656 63.293
+ 24.66 63.293 24.66 63.293 c 24.66 63.293 24.66 63.293 24.66 63.297 c h
+24.66 63.297 m f
+0 g
+q 1 0 0 1 0 0 cm
+24.66 63.297 m 24.66 63.297 24.656 63.297 24.656 63.297 c 24.656 63.293
+ 24.66 63.293 24.66 63.293 c 24.66 63.293 24.66 63.293 24.66 63.297 c h
+24.66 63.297 m S Q
+0 0 1 rg
+23.797 64.039 m 23.797 64.273 23.609 64.465 23.371 64.465 c 23.137 64.465
+ 22.949 64.273 22.949 64.039 c 22.949 63.805 23.137 63.613 23.371 63.613
+ c 23.609 63.613 23.797 63.805 23.797 64.039 c h
+23.797 64.039 m f
+1 0 0 rg
+37.844 42.566 m 37.844 42.945 37.535 43.254 37.156 43.254 c 36.773 43.254
+ 36.465 42.945 36.465 42.566 c 36.465 42.184 36.773 41.875 37.156 41.875
+ c 37.535 41.875 37.844 42.184 37.844 42.566 c h
+37.844 42.566 m f
+0 g
+0.141732 w
+q 1 0 0 1 0 0 cm
+35.977 42.188 m 31.02 39.777 l S Q
+31.531 40.027 m 31.91 39.895 l 30.895 39.715 l 31.66 40.406 l h
+31.531 40.027 m f*
+0.0679747 w
+q 1 0.486446 -0.486446 1 0 0 cm
+41.243 19.965 m 41.497 19.708 l 40.605 19.963 l 41.496 20.22 l h
+41.243 19.965 m S Q
+0.141732 w
+q 1 0 0 1 0 0 cm
+29.828 39.73 m 23.582 45.012 l S Q
+24.016 44.645 m 24.047 44.246 l 23.473 45.102 l 24.414 44.68 l h
+24.016 44.645 m f*
+0.0577315 w
+q 1 -0.845211 0.845211 1 0 0 cm
+-8.002 37.881 m -7.787 37.664 l -8.544 37.88 l -7.787 38.098 l h
+-8.002 37.881 m S Q
+0.141732 w
+q 1 0 0 1 0 0 cm
+23.82 46.355 m 25.027 51.145 l S Q
+24.887 50.594 m 24.543 50.387 l 25.059 51.281 l 25.094 50.25 l h
+24.887 50.594 m f*
+0.0732958 w
+q -0.252177 -1 1 -0.252177 0 0 cm
+-53.469 11.403 m -53.193 11.129 l -54.156 11.402 l -53.195 11.679 l h
+-53.469 11.403 m S Q
+0.141732 w
+q 1 0 0 1 0 0 cm
+23.453 63.156 m 26.16 58.465 l S Q
+25.875 58.953 m 25.98 59.34 l 26.23 58.34 l 25.488 59.059 l h
+25.875 58.953 m f*
+0.0654792 w
+q -0.576787 1 -1 -0.576787 0 0 cm
+33.038 -44.931 m 33.282 -45.177 l 32.424 -44.932 l 33.284 -44.686 l h
+33.038 -44.931 m S Q
+0.141732 w
+q 1 0 0 1 0 0 cm
+26.215 57.066 m 26.281 52.215 l S Q
+26.273 52.785 m 26.555 53.07 l 26.285 52.074 l 25.984 53.062 l h
+26.273 52.785 m f*
+0.0755832 w
+q -0.0137977 1 -1 -0.0137977 0 0 cm
+52.413 -26.997 m 52.694 -27.282 l 51.702 -26.999 l 52.694 -26.711 l h
+52.413 -26.997 m S Q
+0.141732 w
+q 1 0 0 1 0 0 cm
+26.5 50.707 m 34.664 48.957 l S Q
+34.109 49.074 m 33.895 49.41 l 34.805 48.926 l 33.773 48.859 l h
+34.109 49.074 m f*
+0.0739077 w
+q -1 0.214603 -0.214603 -1 0 0 cm
+-22.54 -53.911 m -22.266 -54.188 l -23.235 -53.912 l -22.263 -53.637 l 
+h
+-22.54 -53.911 m S Q
+0.141732 w
+q 1 0 0 1 0 0 cm
+36.809 49.66 m 38.371 43.828 l S Q
+38.223 44.379 m 38.426 44.727 l 38.406 43.691 l 37.879 44.578 l h
+38.223 44.379 m f*
+0.0730145 w
+q -0.26796 1 -1 -0.26796 0 0 cm
+31.85 -46.757 m 32.123 -47.034 l 31.163 -46.757 l 32.122 -46.486 l h
+31.85 -46.757 m S Q
+Q Q
+showpage
+%%Trailer
+end
+%%EOF

doc/user-manual/figures/scatteringQuad.eps

@@ -0,0 +1,303 @@
+%!PS-Adobe-3.0 EPSF-3.0
+%%Creator: cairo 1.16.0 (https://cairographics.org)
+%%CreationDate: Thu Oct 08 15:41:58 2020
+%%Pages: 1
+%%DocumentData: Clean7Bit
+%%LanguageLevel: 2
+%%BoundingBox: 13 30 125 153
+%%EndComments
+%%BeginProlog
+50 dict begin
+/q { gsave } bind def
+/Q { grestore } bind def
+/cm { 6 array astore concat } bind def
+/w { setlinewidth } bind def
+/J { setlinecap } bind def
+/j { setlinejoin } bind def
+/M { setmiterlimit } bind def
+/d { setdash } bind def
+/m { moveto } bind def
+/l { lineto } bind def
+/c { curveto } bind def
+/h { closepath } bind def
+/re { exch dup neg 3 1 roll 5 3 roll moveto 0 rlineto
+      0 exch rlineto 0 rlineto closepath } bind def
+/S { stroke } bind def
+/f { fill } bind def
+/f* { eofill } bind def
+/n { newpath } bind def
+/W { clip } bind def
+/W* { eoclip } bind def
+/BT { } bind def
+/ET { } bind def
+/BDC { mark 3 1 roll /BDC pdfmark } bind def
+/EMC { mark /EMC pdfmark } bind def
+/cairo_store_point { /cairo_point_y exch def /cairo_point_x exch def } def
+/Tj { show currentpoint cairo_store_point } bind def
+/TJ {
+  {
+    dup
+    type /stringtype eq
+    { show } { -0.001 mul 0 cairo_font_matrix dtransform rmoveto } ifelse
+  } forall
+  currentpoint cairo_store_point
+} bind def
+/cairo_selectfont { cairo_font_matrix aload pop pop pop 0 0 6 array astore
+    cairo_font exch selectfont cairo_point_x cairo_point_y moveto } bind def
+/Tf { pop /cairo_font exch def /cairo_font_matrix where
+      { pop cairo_selectfont } if } bind def
+/Td { matrix translate cairo_font_matrix matrix concatmatrix dup
+      /cairo_font_matrix exch def dup 4 get exch 5 get cairo_store_point
+      /cairo_font where { pop cairo_selectfont } if } bind def
+/Tm { 2 copy 8 2 roll 6 array astore /cairo_font_matrix exch def
+      cairo_store_point /cairo_font where { pop cairo_selectfont } if } bind def
+/g { setgray } bind def
+/rg { setrgbcolor } bind def
+/d1 { setcachedevice } bind def
+/cairo_data_source {
+  CairoDataIndex CairoData length lt
+    { CairoData CairoDataIndex get /CairoDataIndex CairoDataIndex 1 add def }
+    { () } ifelse
+} def
+/cairo_flush_ascii85_file { cairo_ascii85_file status { cairo_ascii85_file flushfile } if } def
+/cairo_image { image cairo_flush_ascii85_file } def
+/cairo_imagemask { imagemask cairo_flush_ascii85_file } def
+%%EndProlog
+%%BeginSetup
+%%EndSetup
+%%Page: 1 1
+%%BeginPageSetup
+%%PageBoundingBox: 13 30 125 153
+%%EndPageSetup
+q 13 30 112 123 rectclip
+1 0 0 -1 0 171 cm q
+0 0.784314 0 rg
+20.836 84.285 52.273 43.371 re f
+0 g
+0.751178 w
+0 J
+0 j
+[] 0.0 d
+4 M q 1 0 0 1 0 0 cm
+20.836 84.285 52.273 43.371 re S Q
+1 0.627451 0 rg
+73.109 84.285 43.184 43.371 re f
+0 g
+q 1 0 0 1 0 0 cm
+73.109 84.285 43.184 43.371 re S Q
+0.0431373 0 1 rg
+73.109 32.199 43.184 52.086 re f
+0 g
+q 1 0 0 1 0 0 cm
+73.109 32.199 43.184 52.086 re S Q
+1 0 0 rg
+20.836 32.199 52.273 52.086 re f
+0 g
+q 1 0 0 1 0 0 cm
+20.836 32.199 52.273 52.086 re S Q
+q 1 0 0 1 0 0 cm
+20.836 32.199 95.457 95.457 re S Q
+0 0 1 rg
+22.902 127.77 m 22.902 128.914 21.973 129.844 20.824 129.844 c 19.68 129.844
+ 18.75 128.914 18.75 127.77 c 18.75 126.621 19.68 125.691 20.824 125.691
+ c 21.973 125.691 22.902 126.621 22.902 127.77 c h
+22.902 127.77 m f
+q 1 0 0 1 0 0 cm
+22.902 127.77 m 22.902 128.914 21.973 129.844 20.824 129.844 c 19.68 129.844
+ 18.75 128.914 18.75 127.77 c 18.75 126.621 19.68 125.691 20.824 125.691
+ c 21.973 125.691 22.902 126.621 22.902 127.77 c h
+22.902 127.77 m S Q
+1 0 0 rg
+118.059 127.23 m 118.059 128.379 117.129 129.309 115.98 129.309 c 114.836
+ 129.309 113.906 128.379 113.906 127.23 c 113.906 126.086 114.836 125.156
+ 115.98 125.156 c 117.129 125.156 118.059 126.086 118.059 127.23 c f
+q 1 0 0 1 0 0 cm
+118.059 127.23 m 118.059 128.379 117.129 129.309 115.98 129.309 c 114.836
+ 129.309 113.906 128.379 113.906 127.23 c 113.906 126.086 114.836 125.156
+ 115.98 125.156 c 117.129 125.156 118.059 126.086 118.059 127.23 c S Q
+0 0.784314 0 rg
+117.789 32.41 m 117.789 33.559 116.859 34.488 115.715 34.488 c 114.566 
+34.488 113.637 33.559 113.637 32.41 c 113.637 31.266 114.566 30.336 115.715
+ 30.336 c 116.859 30.336 117.789 31.266 117.789 32.41 c f
+q 1 0 0 1 0 0 cm
+117.789 32.41 m 117.789 33.559 116.859 34.488 115.715 34.488 c 114.566 
+34.488 113.637 33.559 113.637 32.41 c 113.637 31.266 114.566 30.336 115.715
+ 30.336 c 116.859 30.336 117.789 31.266 117.789 32.41 c S Q
+1 0.627451 0 rg
+23.238 32.41 m 23.238 33.559 22.309 34.488 21.16 34.488 c 20.016 34.488
+ 19.086 33.559 19.086 32.41 c 19.086 31.266 20.016 30.336 21.16 30.336 c
+ 22.309 30.336 23.238 31.266 23.238 32.41 c f
+q 1 0 0 1 0 0 cm
+23.238 32.41 m 23.238 33.559 22.309 34.488 21.16 34.488 c 20.016 34.488
+ 19.086 33.559 19.086 32.41 c 19.086 31.266 20.016 30.336 21.16 30.336 c
+ 22.309 30.336 23.238 31.266 23.238 32.41 c S Q
+0 g
+76.465 84.191 m 76.465 86.262 74.785 87.941 72.715 87.941 c 70.645 87.941
+ 68.965 86.262 68.965 84.191 c 68.965 82.121 70.645 80.441 72.715 80.441
+ c 74.785 80.441 76.465 82.121 76.465 84.191 c h
+76.465 84.191 m f
+q 1 0 0 1 0 0 cm
+76.465 84.191 m 76.465 86.262 74.785 87.941 72.715 87.941 c 70.645 87.941
+ 68.965 86.262 68.965 84.191 c 68.965 82.121 70.645 80.441 72.715 80.441
+ c 74.785 80.441 76.465 82.121 76.465 84.191 c h
+76.465 84.191 m S Q
+0.749999 w
+q 1 0 0 1 0 0 cm
+20.355 84.105 m 115.98 84.105 l S Q
+q 1 0 0 1 0 0 cm
+72.922 127.277 m 72.922 32.008 l S Q
+14.742 139.539 m 16.68 139.539 l 16.68 132.867 l 14.57 133.289 l 14.57 
+132.211 l 16.664 131.789 l 17.852 131.789 l 17.852 139.539 l 19.789 139.539
+ l 19.789 140.539 l 14.742 140.539 l h
+14.742 139.539 m f
+120.676 136.695 m 124.816 136.695 l 124.816 137.695 l 119.254 137.695 l
+ 119.254 136.695 l 119.699 136.238 120.309 135.617 121.082 134.836 c 121.863
+ 134.047 122.352 133.535 122.551 133.305 c 122.934 132.891 123.199 132.535
+ 123.348 132.242 c 123.504 131.941 123.582 131.648 123.582 131.367 c 123.582
+ 130.898 123.414 130.52 123.082 130.227 c 122.758 129.938 122.336 129.789
+ 121.816 129.789 c 121.441 129.789 121.043 129.852 120.629 129.977 c 120.223
+ 130.102 119.785 130.301 119.316 130.57 c 119.316 129.367 l 119.793 129.18
+ 120.238 129.039 120.645 128.945 c 121.059 128.844 121.441 128.789 121.785
+ 128.789 c 122.691 128.789 123.414 129.02 123.957 129.477 c 124.496 129.926
+ 124.77 130.531 124.77 131.289 c 124.77 131.645 124.699 131.984 124.566 
+132.305 c 124.43 132.629 124.184 133.008 123.832 133.445 c 123.727 133.563
+ 123.414 133.891 122.895 134.43 c 122.371 134.973 121.633 135.727 120.676
+ 136.695 c h
+120.676 136.695 m f
+122.875 24.07 m 123.438 24.195 123.875 24.453 124.188 24.836 c 124.508 
+25.211 124.672 25.68 124.672 26.242 c 124.672 27.109 124.375 27.781 123.781
+ 28.258 c 123.188 28.727 122.344 28.961 121.25 28.961 c 120.883 28.961 120.504
+ 28.922 120.109 28.852 c 119.723 28.781 119.328 28.672 118.922 28.523 c 
+118.922 27.383 l 119.242 27.57 119.598 27.719 119.984 27.82 c 120.379 27.914
+ 120.789 27.961 121.219 27.961 c 121.957 27.961 122.52 27.816 122.906 27.523
+ c 123.301 27.234 123.5 26.805 123.5 26.242 c 123.5 25.734 123.316 25.332
+ 122.953 25.039 c 122.586 24.75 122.086 24.602 121.453 24.602 c 120.422 
+24.602 l 120.422 23.633 l 121.5 23.633 l 122.07 23.633 122.516 23.52 122.828
+ 23.289 c 123.141 23.051 123.297 22.711 123.297 22.273 c 123.297 21.828 
+123.133 21.484 122.813 21.242 c 122.5 21.004 122.047 20.883 121.453 20.883
+ c 121.117 20.883 120.766 20.922 120.391 20.992 c 120.023 21.055 119.617
+ 21.16 119.172 21.305 c 119.172 20.258 l 119.629 20.133 120.051 20.039 120.438
+ 19.977 c 120.832 19.914 121.203 19.883 121.547 19.883 c 122.453 19.883 
+123.164 20.086 123.688 20.492 c 124.207 20.898 124.469 21.453 124.469 22.148
+ c 124.469 22.641 124.328 23.051 124.047 23.383 c 123.773 23.719 123.383
+ 23.945 122.875 24.07 c h
+122.875 24.07 m f
+17.223 19.746 m 14.238 24.418 l 17.223 24.418 l h
+16.91 18.715 m 18.41 18.715 l 18.41 24.418 l 19.66 24.418 l 19.66 25.402
+ l 18.41 25.402 l 18.41 27.465 l 17.223 27.465 l 17.223 25.402 l 13.285 
+25.402 l 13.285 24.262 l h
+16.91 18.715 m f
+39.906 47.727 m 38.297 52.07 l 41.516 52.07 l h
+39.234 46.555 m 40.578 46.555 l 43.906 55.305 l 42.672 55.305 l 41.875 
+53.055 l 37.938 53.055 l 37.141 55.305 l 35.891 55.305 l h
+39.234 46.555 m f
+45.5 57.066 m 48.188 57.066 l 48.188 57.707 l 44.578 57.707 l 44.578 57.066
+ l 44.867 56.766 45.266 56.359 45.766 55.848 c 46.266 55.34 46.582 55.012
+ 46.719 54.863 c 46.957 54.582 47.125 54.348 47.219 54.16 c 47.32 53.965
+ 47.375 53.777 47.375 53.598 c 47.375 53.297 47.266 53.051 47.047 52.863
+ c 46.836 52.668 46.566 52.566 46.234 52.566 c 45.992 52.566 45.738 52.609
+ 45.469 52.691 c 45.195 52.777 44.91 52.902 44.609 53.066 c 44.609 52.301
+ l 44.922 52.176 45.207 52.082 45.469 52.02 c 45.738 51.957 45.988 51.926
+ 46.219 51.926 c 46.801 51.926 47.27 52.074 47.625 52.363 c 47.977 52.656
+ 48.156 53.047 48.156 53.535 c 48.156 53.777 48.109 54 48.016 54.207 c 47.93
+ 54.418 47.773 54.66 47.547 54.941 c 47.484 55.016 47.281 55.23 46.938 55.582
+ c 46.594 55.938 46.113 56.434 45.5 57.066 c h
+45.5 57.066 m f
+93.52 47.684 m 91.91 52.027 l 95.129 52.027 l h
+92.848 46.512 m 94.191 46.512 l 97.52 55.262 l 96.285 55.262 l 95.488 53.012
+ l 91.551 53.012 l 90.754 55.262 l 89.504 55.262 l h
+92.848 46.512 m f
+98.582 57.02 m 99.832 57.02 l 99.832 52.676 l 98.473 52.957 l 98.473 52.254
+ l 99.832 51.973 l 100.598 51.973 l 100.598 57.02 l 101.848 57.02 l 101.848
+ 57.66 l 98.582 57.66 l h
+98.582 57.02 m f
+45.98 98.063 m 44.371 102.406 l 47.59 102.406 l h
+45.309 96.891 m 46.652 96.891 l 49.98 105.641 l 48.746 105.641 l 47.949
+ 103.391 l 44.012 103.391 l 43.215 105.641 l 41.965 105.641 l h
+45.309 96.891 m f
+53.23 104.98 m 53.605 105.055 53.895 105.215 54.105 105.465 c 54.313 105.715
+ 54.418 106.023 54.418 106.387 c 54.418 106.949 54.219 107.387 53.824 107.699
+ c 53.438 108.004 52.891 108.152 52.184 108.152 c 51.941 108.152 51.699 
+108.125 51.449 108.074 c 51.199 108.035 50.938 107.965 50.668 107.871 c 
+50.668 107.137 l 50.875 107.254 51.105 107.348 51.355 107.418 c 51.613 107.48
+ 51.887 107.512 52.168 107.512 c 52.645 107.512 53.012 107.418 53.262 107.23
+ c 53.52 107.035 53.652 106.754 53.652 106.387 c 53.652 106.055 53.531 105.793
+ 53.293 105.605 c 53.051 105.418 52.723 105.324 52.309 105.324 c 51.652 
+105.324 l 51.652 104.684 l 52.34 104.684 l 52.723 104.684 53.016 104.613
+ 53.215 104.465 c 53.41 104.309 53.512 104.09 53.512 103.809 c 53.512 103.52
+ 53.406 103.293 53.199 103.137 c 52.988 102.98 52.691 102.902 52.309 102.902
+ c 52.098 102.902 51.875 102.93 51.637 102.98 c 51.395 103.023 51.129 103.09
+ 50.84 103.184 c 50.84 102.496 l 51.129 102.414 51.402 102.355 51.652 102.324
+ c 51.91 102.285 52.156 102.262 52.387 102.262 c 52.969 102.262 53.426 102.398
+ 53.762 102.668 c 54.105 102.93 54.277 103.285 54.277 103.73 c 54.277 104.043
+ 54.184 104.309 53.996 104.527 c 53.816 104.746 53.563 104.898 53.23 104.98
+ c h
+53.23 104.98 m f
+93.707 97.875 m 92.098 102.219 l 95.316 102.219 l h
+93.035 96.703 m 94.379 96.703 l 97.707 105.453 l 96.473 105.453 l 95.676
+ 103.203 l 91.738 103.203 l 90.941 105.453 l 89.691 105.453 l h
+93.035 96.703 m f
+100.754 102.836 m 98.801 105.867 l 100.754 105.867 l h
+100.551 102.164 m 101.52 102.164 l 101.52 105.867 l 102.316 105.867 l 102.316
+ 106.508 l 101.52 106.508 l 101.52 107.852 l 100.754 107.852 l 100.754 106.508
+ l 98.176 106.508 l 98.176 105.773 l h
+100.551 102.164 m f
+27.949 77.395 m 27.949 80.301 l 27.043 80.301 l 27.043 72.754 l 27.949 
+72.754 l 27.949 73.582 l 28.137 73.262 28.371 73.02 28.652 72.863 c 28.941
+ 72.707 29.293 72.629 29.699 72.629 c 30.363 72.629 30.902 72.895 31.309
+ 73.426 c 31.723 73.949 31.934 74.637 31.934 75.488 c 31.934 76.355 31.723
+ 77.051 31.309 77.582 c 30.902 78.105 30.363 78.363 29.699 78.363 c 29.293
+ 78.363 28.941 78.285 28.652 78.129 c 28.371 77.973 28.137 77.73 27.949 
+77.395 c h
+31.012 75.488 m 31.012 74.832 30.871 74.316 30.59 73.941 c 30.316 73.566
+ 29.949 73.379 29.48 73.379 c 29 73.379 28.625 73.566 28.355 73.941 c 28.082
+ 74.316 27.949 74.832 27.949 75.488 c 27.949 76.156 28.082 76.676 28.355
+ 77.051 c 28.625 77.426 29 77.613 29.48 77.613 c 29.949 77.613 30.316 77.426
+ 30.59 77.051 c 30.871 76.676 31.012 76.156 31.012 75.488 c h
+31.012 75.488 m f
+35.578 70.629 m 35.148 71.379 34.828 72.125 34.609 72.863 c 34.398 73.594
+ 34.297 74.332 34.297 75.082 c 34.297 75.832 34.398 76.578 34.609 77.316
+ c 34.828 78.059 35.148 78.793 35.578 79.535 c 34.797 79.535 l 34.316 78.773
+ 33.953 78.027 33.703 77.285 c 33.461 76.547 33.344 75.813 33.344 75.082
+ c 33.344 74.355 33.461 73.625 33.703 72.895 c 33.941 72.156 34.305 71.402
+ 34.797 70.629 c h
+35.578 70.629 m f
+41.875 72.754 m 39.906 75.41 l 41.984 78.223 l 40.922 78.223 l 39.328 76.066
+ l 37.734 78.223 l 36.672 78.223 l 38.797 75.363 l 36.859 72.754 l 37.922
+ 72.754 l 39.375 74.707 l 40.813 72.754 l h
+41.875 72.754 m f
+43.121 79.676 m 44.168 79.676 l 44.168 76.066 l 43.027 76.301 l 43.027 
+75.707 l 44.152 75.488 l 44.793 75.488 l 44.793 79.676 l 45.84 79.676 l 
+45.84 80.223 l 43.121 80.223 l h
+43.121 79.676 m f
+47.621 76.988 m 48.652 76.988 l 48.652 77.816 l 47.855 79.379 l 47.215 
+79.379 l 47.621 77.816 l h
+47.621 76.988 m f
+55.117 72.754 m 53.148 75.41 l 55.227 78.223 l 54.164 78.223 l 52.57 76.066
+ l 50.977 78.223 l 49.914 78.223 l 52.039 75.363 l 50.102 72.754 l 51.164
+ 72.754 l 52.617 74.707 l 54.055 72.754 l h
+55.117 72.754 m f
+56.801 79.676 m 59.035 79.676 l 59.035 80.223 l 56.02 80.223 l 56.02 79.676
+ l 56.27 79.426 56.602 79.094 57.02 78.676 c 57.434 78.25 57.699 77.973 
+57.816 77.848 c 58.023 77.621 58.164 77.426 58.238 77.27 c 58.32 77.105 
+58.363 76.941 58.363 76.785 c 58.363 76.535 58.273 76.332 58.098 76.176 
+c 57.918 76.02 57.691 75.941 57.41 75.941 c 57.211 75.941 56.996 75.98 56.77
+ 76.051 c 56.551 76.113 56.316 76.219 56.066 76.363 c 56.066 75.707 l 56.316
+ 75.605 56.551 75.527 56.77 75.473 c 56.996 75.422 57.207 75.395 57.395 
+75.395 c 57.883 75.395 58.273 75.52 58.566 75.77 c 58.855 76.012 59.004 
+76.34 59.004 76.754 c 59.004 76.941 58.965 77.125 58.895 77.301 c 58.82 
+77.48 58.691 77.688 58.504 77.926 c 58.449 77.988 58.277 78.168 57.988 78.457
+ c 57.707 78.75 57.309 79.156 56.801 79.676 c h
+56.801 79.676 m f
+60.484 70.629 m 61.266 70.629 l 61.754 71.402 62.117 72.156 62.359 72.895
+ c 62.609 73.625 62.734 74.355 62.734 75.082 c 62.734 75.813 62.609 76.547
+ 62.359 77.285 c 62.117 78.027 61.754 78.773 61.266 79.535 c 60.484 79.535
+ l 60.922 78.793 61.242 78.059 61.453 77.316 c 61.672 76.578 61.781 75.832
+ 61.781 75.082 c 61.781 74.332 61.672 73.594 61.453 72.863 c 61.242 72.125
+ 60.922 71.379 60.484 70.629 c h
+60.484 70.629 m f
+Q Q
+showpage
+%%Trailer
+end
+%%EOF

doc/user-manual/fpakc_UserManual.tex

@@ -1,5 +1,5 @@
 \documentclass[10pt,a4paper,twoside]{book}
-%\usepackage[latin1]{inputenc}
+\usepackage[latin1]{inputenc}
 \usepackage{amsmath}
 \usepackage{amsfonts}
 \usepackage{amssymb}
@@ -72,7 +72,7 @@
       The \Gls{fpakc} is a simulation tool that models species in plasma (ions, electrons and neutrals) following the trajectories of macro-particles as they move and interact between them and the boundaries of the domain.
       Particles properties are scattered into a finite element mesh in 1, 2 or three dimensions, with the possibility to choose different geometries.
       The official repository can be found at: \url{https://gitlab.com/JorgeGonz/fpakc.git}.
-      The code is currently in the very early steps of development and further refinements are expected very soon.
+      The code is currently in very early steps of development and further improvements are expected very soon.
     
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     \section{Main Guidelines}
@@ -86,11 +86,11 @@
         \item \acrshort{fpakc} is coded in a \textit{understandable} way.
               This means that the code is required to be written in a clear way that is easy to understand and maintain.
               Variables and procedure names need to be self-understanding.
-              This eases the process of fixing bugs and improving the codes by a large team of developers.
+              This ease the process of fixing bugs and improving the codes by a large team of developers.
               For more information, please refer to the \acrshort{fpakc} Coding Style document.
-        \item \acrshort{fpakc} requires being ease to use.
+        \item \acrshort{fpakc} requires to be ease to use.
               Input files are required to be in a \textit{human} format, meaning that the different options can be easily understood without constant reference to the user guide.
-              \acrshort{fpakc} is aimed to be used in a wide range of applications and by various scientists: from well-established ones to newcomers to the field and also students.
+              \acrshort{fpakc} is aimed to be used in a wide range of applications and by a variety of scientist: from very established ones to newcomers to the field and also students.
       \end{enumerate}
     
       These are foundation stones of \acrshort{fpakc} and its development and should always be followed, at least for the releases in the official repository.
@@ -105,16 +105,16 @@
     \section{The Particle Method}
       \Gls{fpakc} uses macro-particles to simulate the dynamics of different plasma species (mainly ions, electrons and neutrals).
       These macro-particles could represent a large amount of real particles.
-      For now own, macro-particles will be referred as just particles by abuse of language.
+      For now own, macro-particles will be referred as just particles by abusing of language.
       During the initiation phase, the input and mesh file(s) are reading.
       If an initial distribution for a species is specified in the input file, particles to match that distribution are loaded into the cells.
       
       The general steps performed in each iteration are:
       \begin{enumerate}
         \item Firstly, new particles are introduced into the domain as specified in the input file.
-        \item Particles are then pushed, accounting for possible acceleration by external forces.
+        \item Particles are then pushed accounting for possible acceleration by external forces.
-              During this process, if a particle changes cell, it is found using the connectivity between elements.
+              During this process, if a particle changes cell it is found using the connectivity between elements.
-              If a particle encounters a boundary instead a new cell, the interaction between the boundary and the wall are computed.
+              If a particle encounters a boundary instead a new cell, the interaction between the boundary and the wall is computed.
               A particle may abandon the computational domain and is no longer accounted for.
         \item Next, collisions for the particles inside each cell are carried out.
               This may include different collision processes for each particle.
@@ -124,10 +124,10 @@
         \item Finally, particle properties are scattered among the mesh nodes.
               These properties are density, momentum and the stress tensor.
         \item If requested, the electromagnetic field is computed.
-        \item If the number of iteration requires writing output files, it is done after all steps for the particles are completed.
+        \item If the number of iteration requires writing output files, it is done after all steps for the particles is completed.
       \end{enumerate}
     
-      \Gls{fpakc} has the capability to configure all the behaviour of the simulation via the input file.
+      \Gls{fpakc} has the capability to configure all the behavior of the simulation via the input file.
       Parameters as injection, the kind of pusher used for each species, boundary conditions or collisions are user-input parameters and will be described in Chap.~\ref{ch:input_file}.
     
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -168,8 +168,8 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     \section{Find new cell}      
       Once the position and velocity of the particle are updated, the new cell that contains the particle is searched.
-      This is done by a neighbour search, starting from the previous cell containing the particle.
+      This is done by a neighbor search, starting from the previous cell containing the particle.
-      In the process of finding the new cell, a particle might encounter a boundary.
+      In the process of finding the new cell, it is possible that a particle encounters a boundary.
       When the particle interacts with the boundary, the particle may continue its life in the simulation or might be eliminated from it.
       Once that the new cell is found or that the particle life has been terminated, the pushing is complete.
       If a secondary mesh is used for the Monte-Carlo Collision method, the new cell in that mesh in which the particle reside is also found by the same method, although no interaction with the boundaries is accounted for this step.
@@ -178,7 +178,7 @@
     \section{Variable Weighting Scheme\label{sec:weightingScheme}}
       One of the issues in particle simulations, specially for axial-symmetrical cases, is that due to the disparate volume of cells, specially close to the axis, the statistics in some cells is usually poor.
       To try to fix that, the possibility to include a Variable Weighting Scheme in the simulations is available in \Gls{fpakc}.
-      These schemes detect when a particle changes cells and split it if necessary to improve statistics.
+      These schemes detect when a particle change cells and split it if necessary to improve statistics.
       The use of a Variable Weighting Scheme is defined by the user in the input file.
       
       Beware that this can increase the number of particles in the simulation and increase computational time.
@@ -189,16 +189,16 @@
       \Gls{fpakc} distinguish between two types of interactions: \acrfull{mcc} and \acrfull{cs}.
       \acrshort{mcc} refers to the process in which two particles interact in short range.
       These processes include, but are not limited to: elastic collisions, ionization/recombination, charge-exchange, excitation/de-excitation\ldots
-      A secondary mesh, with cell sizes in the range of the mean-free path, can be used for this type of collision.
+      A secondary mesh, with cell sizes in the range of the mean-free path, can be used for this type of collisions.
       \acrshort{cs} refers to the large range interaction that a charged species suffer do to the charge of other particles.
       The interactions between the different species is defined by the user.
       
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       \subsection{\acrlong{mcc}}
-        For each cell, the maximum number of collisions between particle is computed.
+        For each cell the maximum number of collisions between particle is computed.
         For each collision, a random pair of particles is chosen.
         A loop over all possible collisions for the pair of particles chosen is performed.
-        If a random number is above the probability of collision for that specific type, the collision takes place.
+        If a random number is above the probability of collision for that specific type, the collision take place.
         If not, the next type for the particle pair is checked.
               
         Below are described the type of collision process implemented in \acrshort{fpakc}:
@@ -219,7 +219,7 @@
 
           \item Recombination.
                 When an electron and an ion interact, there is a possibility for them to be recombined into a neutral particle.
-                The photons emitted by this process are not modelled yet.
+                The photon emitted by this process is not modelled yet.
         \end{itemize}
       
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -238,7 +238,7 @@
       Once that the pushing is complete, the array of particles that remain inside the domain is copied to a new array.
       The new array containing only the particles inside the domain will be the one used in the next steps.
       In this section, particles are assigned to the list of particles inside each individual cell.
-      Unfortunately, this is done right now without parallelization and is very CPU consuming.
+      Unfortunately, this is done right now without parallelisation and is very CPU consuming.
       
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     \section{Probing}\label{sec:probing}
@@ -250,21 +250,21 @@
       The user can decide the grid width and the number of points in each direction.
       The distribution function will be calculated and wrote with a time step decided by the user.
       
-      If a particle velocity resides outside the velocity grid (in any direction), it will not be added to the tally of the distribution function.
+      If a particle velocity resides outside of the velocity grid (in any direction), it wont be added to the tally of the distribution function.
-      Due to the limitation of only considering particles in the cell, and not neighbour particles, two probes for the same species at different positions but in the same cell will output the same results.
+      Due to the limitation of only taking into account particles in the cell, and not neighbour particles, two probes for the same species at different positions but in the same cell will output the same results.
-      A more advance method considering distance between the particles and the probe position as well as particles in neighbour cells could be implemented to improve the statistics of the distribution function.
+      A more advance method taking into account distance between the particles and the probe position as well as particles in neighbour cells could be implemented to improve the statistics of the distribution function.
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     \section{Scattering}
       The properties of each particle are deposited in the nodes from the containing cell.
-      This process depends on the cell type, but in general, each node receives a proportional part of the particle properties as a function of the particle position inside the cell.
+      This process depend on the cell type, but in general, each node receive a proportional part of the particle properties as a function of the particle position inside the cell.
-      The figure \ref{fig:scatteringQuad} shows how a particle at a generic position $p(x_1, x_2)$ inside the cell is scattered to the four nodes.
+      Figure \ref{fig:scatteringQuad} shows how a particle at a generic position $p(x_1, x_2)$ inside the cell is scattered to the four nodes.
       \begin{wrapfigure}{l}{0.4\textwidth}
         \centering
         \includegraphics{figures/scatteringQuad}
         \caption{\label{fig:scatteringQuad}Example of how a particle is weighted in a quadrilateral cell.}
       \end{wrapfigure}
-      Each node receives a proportional part of the area formed by dividing the cell in for rectangles, using as an additional vertex the particle position.
+      Each node receives a proportional part of the area formed by dividing the cell in for rectangles using as an additional vertex the particle position.
       These properties are dimensionless, but they are converted to the correct units once the output is printed.
     
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -273,11 +273,11 @@
     
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     \section{Average scheme}
-      Particle-in-cell codes have an intrinsic statistical noise associated with them.
+      Particle-in-cell codes has an intrinsic statistical noise associated with them.
       Although this can be reduced by increasing the number of particles, this also increases the CPU requirements of the case.
       
       It is quite common that most cases reach a quasi-steady state after a number of iterations and time-average results can be obtained after to improve analysis, plotting and restarting the case using these time-average results as new species backgrounds.
-      Although this is possible to do once the simulation is finished with post-processing tools, this is limited to the number of iterations printed.
+      Although this is possible to do once the simulation is finished with post-processing tools, this is limited to the amount of iterations printed.
       
       \Gls{fpakc} implements a simple average scheme that, after a start time provided by the user, scores a mean and standard deviation of all the main species properties, and the electromagnetic field.
       This scheme is based on the Welford's online algorithm~\cite{welford1962note}.
@@ -286,7 +286,7 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
   \chapter{Installation}
     \section{Required Packages}
-      To properly compile \gls{fpakc}, the following packages are required.
+      In order to properly compile \gls{fpakc}, the following packages are required.
       \subsection{Gfortran}
         The \Gls{opensource} free compiler \Gls{gfortran}\cite{gfortranURL} from GCC is the basic way to compile \acrshort{fpakc}.
         It is distributed with all GNU/Linux distributions.
@@ -369,7 +369,7 @@ make
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     \section{Case file}
       The required format for the case file is \Gls{json}.
-      \Gls{json} is a case-sensitive format, so input must be written with the correct capitalization.
+      \Gls{json} is a case-sensitive format, so input must be written with the correct capitalisation.
       The basic structure and options available for the case file are explained below.
       The order of the objects and variables is irrelevant, but the structure needs to be maintained.
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -380,9 +380,9 @@ make
           \item \textbf{path}: Character.
                                Path for the output files. This path is also used to locate the mesh input file.
           \item \textbf{folder}: Character.
-                                 Base name of the folder in which output files are placed.
+                                 Base name of the folder in wich output files are placed.
                                  The date and time is appended to this name.
-                                 If none is provided, only the date and time is written as the folder name.
+                                 If none is provided, only the date and time is writted as the folder name.
           \item \textbf{triggerOutput}: Integer.
                                         Determines the number of iterations between writing output files for macroscopic quantities.
           \item \textbf{cpuTime}: Logical. 
@@ -460,10 +460,6 @@ make
                                    \begin{itemize}
                                      \item \textbf{gmsh2}: \Gls{gmsh} file format in version 2.0. This has to be in ASCII format.
                                      \item \textbf{vtu}: \Gls{vtu} file format. This has to be in ASCII format.
-                                     \item \textbf{text}: Plain text file format only intended for 1D cases.
-                                                          This has to be in ASCII format and comma separated.
-                                                          The first column represents the position and the second column the physical ID of the node.
-                                                          Values have to be $1$ (left boundary), $2$ (right boundary), or $0$ (no boundary.)
                                    \end{itemize}
           \item \textbf{meshFile}: Character.
                                    Mesh filename.
@@ -519,7 +515,7 @@ make
                                   \item \textbf{absorption}: Particle is eliminated from the domain.
                                                              The particle is first moved into the edge and its properties are scattered among the edge nodes.
                                   \item \textbf{transparent}: Particle abandon the numerical domain.
-                                  \item \textbf{wallTemperature}: Reflective wall with constant temperature that exchange heat with particles.
+                                  \item \textbf{wallTemperature}: Reflective wall with cosntant temperature that exchange heat with particles.
                                                                   Required parameters are:
                                                                   \begin{itemize}
                                                                     \item \textbf{temperature}: Real.
@@ -530,8 +526,8 @@ make
                                                                                                  Specific heat capacity of the material.
                                                                   \end{itemize}
                                   \item \textbf{ionization}: Per each particle crossing the surface with this type of boundary, a number of ionization events are calculated.
-                                                             A pair of ion-electron is generated for each ionization event, taking as a reference a neutral background.
+                                                             A pair of ion-electron is generated for each ionization event taking as a reference a neutral background.
-                                                             The secondary electron is taken as the same type as the incident particle.
+                                                             Secondary electron is taken as same type as incident particle.
                                                              The available input is:
                                                              \begin{itemize}
                                                                \item \textbf{neutral}: Object.
@@ -544,7 +540,7 @@ make
                                                                                          \item \textbf{mass}: Real.
                                                                                                               Units in $\unit{kg}$.
                                                                                                               Mass of neutral species.
-                                                                                                              If missing, the mass of the ion is used
+                                                                                                              If missing, the mass of the ion is ussed
                                                                                          \item \textbf{density}: Real.
                                                                                                                  Units in $\unit{m^{-3}}$.
                                                                                                                  Density of neutral background.
@@ -562,18 +558,18 @@ make
                                                                                        \end{itemize}
                                                                \item \textbf{effectiveTime}: Real.
                                                                                              Units in $\unit{s}$.
-                                                                                             As the particle is no longer simulated once it crossed the boundary, this time represents the effective time in which the particle produces ionization processes in the neutral background.
+                                                                                             As the particle is no longer simulated once it crossed the boundary, this time represent the effective time in which the particle produces ionization processes in the neutral background.
                                                                                              Required parameter.
                                                                \item \textbf{energyThreashold}: Real.
                                                                                                 Units in $\unit{eV}$.
                                                                                                 Ionization energy threshold for the simulated process.
                                                                                                 Required parameter.
                                                                \item \textbf{crossSection}: Character.
-                                                                                            Complete path to the cross-section data for the ionization process.
+                                                                                            Complete path to the cross section data for the ionization process.
                                                                                              
                                                              \end{itemize}
                                   \item \textbf{axis}: Identifies the symmetry axis for 2D cylindrical simulations.
-                                                       If , for some reason, a particle interacts with this axis, it is reflected.
+                                                       If for some reason a particle interact with this axis, it is reflected.
                                 \end{itemize}
         \end{itemize}
       
@@ -589,26 +585,18 @@ make
                                Type of boundary.
                                 Accepted values are:
                                \begin{itemize}
-                                 \item \textbf{dirichlet}: Constant value of electric potential on the surface.
+                                 \item \textbf{dirichlet}: Elastic reflection of particles.
-                                 \item \textbf{dirichletTime}: Constant value of the electric potential with a time variable profile.
-                                                               The value of \textbf{boundaryEM.potential} will be multiplied for the corresponding value in the file \textbf{boundaryEM.temporalProfile}.
                                 \end{itemize}
-          \item \textbf{potential}: Real.
+           \item \textbf{potential}: Real.
-                                    Fixed potential for Dirichlet boundary condition.
+                                     Fixed potential for Dirichlet boundary condition.
           \item \textbf{physicalSurface}: Integer.
                                           Identification of the edge in the mesh file.
-          \item \textbf{temporalProfile}: Character.
-                                          Filename of the 2 column file containing the time variable profile.
-                                          File must be located in \textbf{output.path}.
-                                          The first column is the time in $\unit{s}$.
-                                          The second column is the factor that will multiply the value of the boundary.
-        
         \end{itemize}
       
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       \subsection{inject}
         The array \textbf{inject} specifies the injection of particles from different surfaces.
-        The injection of particles needs to be associated to a physicalSurface in the mesh file.
+        The injection of particles need to be associated to a physicalSurface in the mesh file.
         Multiple injections can be associated to the same surface.
         \begin{itemize}
           \item \textbf{name}: Character.
@@ -622,9 +610,7 @@ make
                                 Available values are:
                                 \begin{itemize}
                                   \item \textbf{A}: Ampere.
-                                  \item \textbf{Am2}: Ampere per square meter.
+                                  \item \textbf{sccm}: Standard cubic centimeter.
-                                                      This value will be multiplied by the area of injection.
-                                  \item \textbf{sccm}: Standard cubic centimetre.
                                   \item \textbf{part/s}: Particles (real) per second.
                                 \end{itemize}
           \item \textbf{v}: Real.
@@ -641,7 +627,7 @@ make
                                   \begin{itemize}
                                     \item \textbf{Maxwellian}: Maxwellian distribution of temperature \textbf{T} and mean \textbf{v} times the value of \textbf{n} in the specified direction.
                                     \item \textbf{Half-Maxwellian}: Half-Maxwellian distribution of temperature \textbf{T} and mean \textbf{v} times the value of \textbf{n} in the specified direction.
-                                                                    Only considers the positive part of the half-Maxwellian.
+                                                                    Only takes into account the positive part of the half-Maxwellian.
                                     \item \textbf{Delta}: Dirac's delta distribution function. All particles are injected with velocity \textbf{v} times the value of \textbf{n} in the specified direction.
                                   \end{itemize}
           \item \textbf{T}: Real.
@@ -650,11 +636,6 @@ make
                             Temperature in each direction.
           \item \textbf{physicalSurface}: Integer.
                                           Identification of the edge in the mesh file.
-          \item \textbf{particlesPerEdge}: Integer.
-                                           Optional.
-                                           Number of particles to be injected by each edge in the numerical domain.
-                                           The weight of the particles for each edge will modified by the surface of the edge to ensure the right flux is injected.
-                                           If no value is provided, the number of particles to inject per edge will be calculated with the species weight and the surface of the edge respect to the total one.
         \end{itemize}        
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       \subsection{reference}
@@ -670,7 +651,7 @@ make
           \item \textbf{radius}: Real.
                                  Reference atomic radius in $\unit{m}$.
           \item \textbf{crossSection}: Real.
-                                       Reference cross-section in $\unit{m^2}$.
+                                       Reference cross section in $\unit{m^2}$.
                                        If this value is present, radius is ignored.
         \end{itemize}
       
@@ -696,8 +677,8 @@ make
                                  Indicates the type of pusher used for each species:
                                  \begin{itemize}
                                    \item \textbf{Neutral}: Pushes a particle without any external force.
-                                   \item \textbf{Electrostatic}: Pushes a particle, including the effect of the electrostatic field.
+                                   \item \textbf{Electrostatic}: Pushes a particle including the effect of the electrostatic field.
-                                   \item \textbf{Electromagnetic}: Pushes a particle, accounting for the electromagnetic field.
+                                   \item \textbf{Electromagnetic}: Pushes particles accounting for the electromagnetic field.
                                  \end{itemize}
           \item \textbf{WeightingScheme}: Character.
                                           Indicates the variable weighting scheme to be used in the simulation.
@@ -725,15 +706,11 @@ make
                                   \begin{itemize}
                                     \item \textbf{species}: Character.
                                                             Name of species as defined in the object \textbf{species}.
-                                    \item \textbf{file}: Character.
+                                    \item \textbf{file}:    Character.
-                                                         Output file from previous run used as an initial state for the species.
+                                                            Output file from previous run used as an initial state for the species.
-                                                         The file format must be the same as in \textbf{geometry.meshType}
+                                                            The file format must be the same as in \textbf{geometry.meshType}
-                                                         Initial particles are assumed to have a Maxwellian distribution.
+                                                            Initial particles are assumed to have a Maxwellian distribution.
-                                                         File must be located in \textbf{output.path}.
+                                                            File must be located at \textbf{output.path}.
-                                    \item \textbf{particlesPerCell}: Integer.
-                                                                     Optional.
-                                                                     Initial number of particles per cell.
-                                                                     If not, the number of particles per cell will be assigned based on the species weight and the cell volume.
                                   \end{itemize}
         \end{itemize}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -749,11 +726,11 @@ make
         \end{itemize}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       \subsection{interactions}\label{ssec:input_interactions}
-        This object determines the different interactions among species.
+        This object determine the different interactions among species.
         Acceptable values are:
         \begin{itemize}
           \item \textbf{folderCollisions}: Character.
-                                           Indicates the path to in which the cross-section tables are allocated.
+                                           Indicates the path to in which the cross section tables are allocated.
           \item \textbf{meshCollisions}: Character.
                                          Determines a specific mesh for \acrshort{mcc} processes.
                                          The file needs to be located in the folder \textbf{output.folder}.
@@ -780,18 +757,13 @@ make
                                                                                      Accepted values are \textbf{elastic}, \textbf{chargeExchange}, \textbf{ionization} and \textbf{recombination}.
                                                                                      Please refer to Sec.~\ref{ssec:collisions} for a description of the different collision types.
                                                                 \item \textbf{crossSection}: Character.
-                                                                                             File in \textbf{interactions.folderCollisions} that contains the cross-section data as a 1D table of relative energy (in $\unit{eV}$) and cross-section (in $\unit{m^-2}$).
+                                                                                             File in \textbf{interactions.folderCollisions} that contains the cross section data as a 1D table of relative energy (in $\unit{eV}$) and cross section (in $\unit{m^-2}$).
                                                                 \item \textbf{energyThreshold}: Real.
                                                                                                 Energy threshold of the collisional process in $\unit{eV}$.
                                                                                                 Only valid for \textbf{ionization} and \textbf{recombination} processes.
                                                                 \item \textbf{electron}: Character.
                                                                                          Name of species designed as electrons.
                                                                                          Only valid for \textbf{ionization} and \textbf{recombination} processes.
-                                                                \item \textbf{electronSecondary}: Character.
-                                                                                                  Optional.
-                                                                                                  Name of species designed as secondary electrons.
-                                                                                                  If none provided, \textbf{electron} is used.
-                                                                                                  Only valid for \textbf{ionization}.
                                                               \end{itemize}
                                      \end{itemize}
           \item \textbf{Coulomb}: Array of objects.
@@ -801,7 +773,7 @@ make
                                   \begin{itemize}
                                     \item \textbf{species\_i}, \textbf{species\_j}: Character.
                                                                                     Define the two species involved in the collision processes.
-                                                                                    Order is indifferent.
+                                                                                    Order is indiferent.
                                   
                                   \end{itemize}
         \end{itemize}
@@ -827,9 +799,9 @@ make
     
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     \section{1D Emissive Cathode (1D\_Cathode)}
-      Emission from a 1D cathode in both, Cartesian and radial coordinates.
+      Emission from a 1D cathode in both, cartesian and radial coordinates.
-      Both cases insert the same number of electrons from the minimum coordinate and have the same boundary conditions for particles and the electrostatic field.
+      Both cases insert the same amount of electrons from the minimum coordinate and have the same boundary conditions for particles and the electrostatic field.
-      This case is useful to illustrate how \acrshort{fpakc} can deal with different geometries by just modifying some parameters in the input file.
+      This case is useful to ilustrate hoy \acrshort{fpakc} can deal with different geometries by just modifying some parameters in the input file.
       The same mesh file (\lstinline|mesh.msh|) is used for both cases.
       
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

src/fpakc.f90

@@ -11,12 +11,12 @@ PROGRAM fpakc
   USE OMP_LIB
   IMPLICIT NONE
 
+  ! t = time step
+  INTEGER:: t
   ! arg1 = Input argument 1 (input file)
   CHARACTER(200):: arg1
   ! inputFile = path+name of input file
   CHARACTER(:), ALLOCATABLE:: inputFile
-  ! generic integer for time step
-  INTEGER:: t
 
   tStep = omp_get_wtime()
   !Gets the input file
@@ -27,18 +27,11 @@ PROGRAM fpakc
 
   !Reads the json configuration file
   CALL readConfig(inputFile)
-
-  !Create output folder and initial files
-  CALL initOutput(inputFile)
-
   !Do '0' iteration
-  timeStep = tInitial
+  t = tInitial
 
   !$OMP PARALLEL DEFAULT(SHARED)
   !$OMP SINGLE
-  ! Initial reset of probes
-  CALL resetProbes()
-
   CALL verboseError("Initial scatter of particles...")
   !$OMP END SINGLE
   CALL doScatter()
@@ -52,21 +45,19 @@ PROGRAM fpakc
   tStep = omp_get_wtime() - tStep
 
   !Output initial state
-  CALL doOutput()
+  CALL doOutput(t)
   CALL verboseError('Starting main loop...')
   !$OMP PARALLEL DEFAULT(SHARED) 
   DO t = tInitial + 1, tFinal
+    !Insert new particles and push them
     !$OMP SINGLE
     tStep = omp_get_wtime()
 
-    ! Update global time step index
-    timeStep = t
-
     !Checks if a species needs to me moved in this iteration
-    CALL solver%updatePushSpecies()
+    CALL solver%updatePushSpecies(t)
 
     !Checks if probes need to be calculated this iteration
-    CALL resetProbes()
+    CALL resetProbes(t)
     tPush = omp_get_wtime()
     !$OMP END SINGLE
 
@@ -84,7 +75,7 @@ PROGRAM fpakc
     !$OMP END SINGLE
 
     IF (doMCCollisions) THEN
-      CALL meshForMCC%doCollisions()
+      CALL meshForMCC%doCollisions(t)
 
     END IF
 
@@ -129,12 +120,12 @@ PROGRAM fpakc
     !$OMP SINGLE
     tEMField = omp_get_wtime() - tEMField
 
-    CALL doAverage()
+    CALL doAverage(t)
 
     tStep = omp_get_wtime() - tStep
 
     !Output data
-    CALL doOutput()
+    CALL doOutput(t)
     !$OMP END SINGLE
 
   END DO

src/makefile

@@ -9,7 +9,6 @@ OBJECTS = $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleMeshBoundary.o $(OBJDIR)/module
 			  	$(OBJDIR)/moduleMeshInputVTU.o $(OBJDIR)/moduleMeshOutputVTU.o \
 			  	$(OBJDIR)/moduleMeshInputGmsh2.o $(OBJDIR)/moduleMeshOutputGmsh2.o \
 			  	$(OBJDIR)/moduleMeshInput0D.o $(OBJDIR)/moduleMeshOutput0D.o \
-				$(OBJDIR)/moduleMeshInputText.o $(OBJDIR)/moduleMeshOutputText.o \
 			  	$(OBJDIR)/moduleMesh3DCart.o \
 			  	$(OBJDIR)/moduleMesh2DCyl.o  \
 			  	$(OBJDIR)/moduleMesh2DCart.o \

src/modules/common/moduleCaseParam.f90

@@ -2,13 +2,8 @@
 MODULE moduleCaseParam
   !Final and initial iterations
   INTEGER:: tFinal, tInitial = 0
-  ! Global index of current iteration
-  INTEGER:: timeStep
-  ! Time step for all species
   REAL(8), ALLOCATABLE:: tau(:)
-  ! Minimum time step
   REAL(8):: tauMin
-  ! Time step for Monte-Carlo Collisions
   REAL(8):: tauColl
 
 END MODULE moduleCaseParam

src/modules/common/moduleRandom.f90

@@ -40,46 +40,47 @@ MODULE moduleRandom
       INTEGER:: rnd
       REAL(8):: rnd01
 
-      rnd = 0
+      rnd = 0.D0
       CALL RANDOM_NUMBER(rnd01)
 
-      rnd = a + FLOOR((b+1-a)*rnd01)
+      rnd = INT(REAL(b - a) * rnd01) + 1
 
     END FUNCTION randomIntAB
 
-    !Returns a random number in a Maxwellian distribution of mean 0 and width 1 with the Box-Muller Method
-    function randomMaxwellian() result(rnd)
-      USE moduleConstParam, only: pi
-      implicit none
-
-      real(8):: rnd
-      real(8):: v1, v2, Rsquare
-
-      v1 = 0.d0
-      do while (v1 <= 0.d0)
-        v1 = random()
-
-      end do
-      v2 = random()
-
-      rnd = sqrt(-2.d0*log(v1))*cos(2*pi*v2)
-
-    end function randomMaxwellian
-
     !Returns a random number in a Maxwellian distribution of mean 0 and width 1
-    FUNCTION randomHalfMaxwellian() RESULT(rnd)
+    FUNCTION randomMaxwellian() RESULT(rnd)
+      USE moduleConstParam, ONLY: PI
       IMPLICIT NONE
 
       REAL(8):: rnd
-      REAL(8):: x
+      REAL(8):: x, y
 
       rnd = 0.D0
       x = 0.D0
       DO WHILE (x == 0.D0)
         CALL RANDOM_NUMBER(x)
       END DO
+      CALL RANDOM_NUMBER(y)
 
-      rnd = DSQRT(-DLOG(x))
+      rnd = DSQRT(-2.D0*DLOG(x))*DCOS(2.D0*PI*y)
+
+    END FUNCTION randomMaxwellian
+
+    FUNCTION randomHalfMaxwellian() RESULT(rnd)
+      USE moduleConstParam, ONLY: PI
+      IMPLICIT NONE
+
+      REAL(8):: rnd
+      REAL(8):: x, y
+
+      rnd = 0.D0
+      x = 0.D0
+      DO WHILE (x == 0.D0)
+        CALL RANDOM_NUMBER(x)
+      END DO
+      y = random(-PI, PI)
+
+      rnd = DSQRT(-2.D0*DLOG(x))*DCOS(y)
 
     END FUNCTION randomHalfMaxwellian
 
@@ -90,21 +91,10 @@ MODULE moduleRandom
       REAL(8), INTENT(in):: cumWeight(1:)
       REAL(8), INTENT(in):: sumWeight
       REAL(8):: rnd0b
-      INTEGER:: rnd, i
+      INTEGER:: rnd
 
-      rnd0b = random()
+      rnd0b = random(0.D0, sumWeight)
-      i = 1
+      rnd = MINLOC(DABS(rnd0b - cumWeight), 1)
-      DO 
-        IF (rnd0b <= cumWeight(i)/sumWeight) THEN
-          rnd = i
-          EXIT
-
-        ELSE
-          i = i +1
-
-        END IF
-      END DO
-      ! rnd = MINLOC(DABS(rnd0b - cumWeight), 1)
 
     END FUNCTION randomWeighted
 

src/modules/common/moduleTable.f90

@@ -30,6 +30,8 @@ MODULE moduleTable
         READ(id, '(A)', iostat = stat) dummy
         !If EOF or error, exit file
         IF (stat /= 0) EXIT
+        !If empty line, exit file
+        IF (dummy == "") EXIT
         !Skip comment
         IF (INDEX(dummy,'#') /= 0) CYCLE
         !Add row
@@ -55,6 +57,7 @@ MODULE moduleTable
         !TODO: Make this a function
         IF (INDEX(dummy,'#') /= 0) CYCLE
         IF (stat /= 0) EXIT
+        IF (dummy == "") EXIT
         !Add data
         !TODO: substitute with extracting information from dummy
         BACKSPACE(id)
@@ -93,7 +96,7 @@ MODULE moduleTable
         f = self%fMax
 
       ELSE
-        i = MINLOC(ABS(x - self%x), 1)
+        i = MINLOC(x - self%x, 1)
         deltaX = x - self%x(i)
         IF (deltaX < 0 ) THEN 
           i = i - 1

src/modules/init/moduleInput.f90

@@ -84,6 +84,20 @@ MODULE moduleInput
       CALL readParallel(config)
       CALL checkStatus(config, "readParallel")
 
+      !If everything is correct, creates the output folder
+      CALL EXECUTE_COMMAND_LINE('mkdir ' // path // folder )
+      !Copies input file to output folder
+      CALL EXECUTE_COMMAND_LINE('cp ' // inputFile // ' ' // path // folder)
+      !Copies particle mesh
+      IF (mesh%dimen > 0) THEN
+        CALL EXECUTE_COMMAND_LINE('cp ' // pathMeshParticle // ' ' // path // folder)
+        IF (doubleMesh) THEN
+          CALL EXECUTE_COMMAND_LINE('cp ' // pathMeshColl // ' ' // path // folder)
+
+        END IF
+
+      END IF
+
     END SUBROUTINE readConfig
     
     !Checks the status of the JSON case file and, if failed, exits the execution.
@@ -259,17 +273,13 @@ MODULE moduleInput
           !Read BC
           CALL readEMBoundary(config)
 
-        CASE("ElectrostaticBoltzmann")
-          !Read BC
-          CALL readEMBoundary(config)
-
         CASE("ConstantB")
           !Read BC
           CALL readEMBoundary(config)
           !Read constant magnetic field
           DO i = 1, 3
-            WRITE(iString, '(i2)') i
+            WRITE(istring, '(i2)') i
-            CALL config%get(object // '.B(' // iString // ')', B(i), found)
+            CALL config%get(object // '.B(' // istring // ')', B(i), found)
             IF (.NOT. found) THEN
               CALL criticalError('Constant magnetic field not provided in direction ' // iString, 'readSolver')
 
@@ -312,7 +322,7 @@ MODULE moduleInput
       LOGICAL:: found
       CHARACTER(:), ALLOCATABLE:: object
       INTEGER:: nInitial
-      INTEGER:: i, p, e
+      INTEGER:: i, j, p, e
       CHARACTER(LEN=2):: iString
       CHARACTER(:), ALLOCATABLE:: spName
       INTEGER:: sp
@@ -328,12 +338,10 @@ MODULE moduleInput
       REAL(8):: densityCen
       !Mean velocity and temperature at particle position
       REAL(8):: velocityXi(1:3), temperatureXi
-      INTEGER:: nNewPart = 0
+      INTEGER:: nNewPart = 0.D0
-      REAL(8):: weight = 0.D0
-      CLASS(meshCell), POINTER:: cell
       TYPE(particle), POINTER:: partNew
+      CLASS(meshCell), POINTER:: cell
       REAL(8):: vTh
-      TYPE(lNode), POINTER:: partCurr, partNext
 
       CALL config%info('solver.initial', found, n_children = nInitial)
 
@@ -348,9 +356,6 @@ MODULE moduleInput
           !Reads node values at the nodes
           filename = path // spFile
           CALL mesh%readInitial(filename, density, velocity, temperature)
-          !Check if initial number of particles is given
-          CALL config%get(object // '.particlesPerCell', nNewPart, found)
-
           !For each volume in the node, create corresponding particles
           DO e = 1, mesh%numCells
             !Scale variables
@@ -363,15 +368,14 @@ MODULE moduleInput
             densityCen = mesh%cells(e)%obj%gatherF((/ 0.D0, 0.D0, 0.D0 /), nNodes, sourceScalar)
 
             !Calculate number of particles
-            IF (.NOT. found) THEN
+            nNewPart = INT(densityCen * (mesh%cells(e)%obj%volume*Vol_ref) / species(sp)%obj%weight)
-              nNewPart = FLOOR(densityCen * (mesh%cells(e)%obj%volume*Vol_ref) / species(sp)%obj%weight)
 
-            END IF
+            !Allocate array of particles for this species
-            weight = densityCen * (mesh%cells(e)%obj%volume*Vol_ref) / REAL(nNewPart)
+            ALLOCATE(partOld(sp)%p(1:nNewPart))
 
-            !Allocate new particles
+            !Create new particles
             DO p = 1, nNewPart
-              ALLOCATE(partNew)
+              partNew = partOld(sp)%p(p)
               partNew%species => species(sp)%obj
               partNew%r       =  mesh%cells(e)%obj%randPos()
               partNew%Xi      =  mesh%cells(e)%obj%phy2log(partNew%r)
@@ -404,10 +408,7 @@ MODULE moduleInput
 
               partNew%n_in   = .TRUE.
 
-              partNew%weight = weight
+              partNew%weight = species(sp)%obj%weight
-
-              !Assign particle to temporal list of particles
-              CALL partInitial%add(partNew)
 
               !Assign particle to list in volume
               IF (listInCells) THEN
@@ -428,30 +429,10 @@ MODULE moduleInput
 
           END DO
 
+          nPartOld(sp) = SIZE(partOld(sp)%p)
+
         END DO
 
-        !Convert temporal list of particles into initial partOld array
-        !Deallocate the list of initial particles
-        nNewPart = partInitial%amount
-        IF (nNewPart > 0) THEN
-          ALLOCATE(partOld(1:nNewPart))
-          partCurr => partInitial%head
-          DO p = 1, nNewPart
-            partNext => partCurr%next
-            partOld(p) = partCurr%part
-            DEALLOCATE(partCurr)
-            partCurr => partNext
-
-          END DO
-
-          IF (ASSOCIATED(partInitial%head)) NULLIFY(partInitial%head)
-          IF (ASSOCIATED(partInitial%tail)) NULLIFY(partInitial%tail)
-          partInitial%amount = 0
-
-        END IF
-
-        nPartOld = SIZE(partOld)
-
       END IF
 
     END SUBROUTINE readInitial
@@ -598,7 +579,10 @@ MODULE moduleInput
 
       END DO
 
+      !Allocate the wrapper array that contains particles
+      ALLOCATE(partOld(1:nSpecies))
       !Set number of particles to 0 for init state
+      ALLOCATE(nPartOld(1:nSpecies))
       nPartOld = 0
 
       !Initialize the lock for the non-analogue (NA) list of particles
@@ -632,7 +616,7 @@ MODULE moduleInput
       INTEGER:: i, k, ij
       INTEGER:: pt_i, pt_j
       REAL(8):: energyThreshold, energyBinding
-      CHARACTER(:), ALLOCATABLE:: electron, electronSecondary
+      CHARACTER(:), ALLOCATABLE:: electron
       INTEGER:: e
       CLASS(meshCell), POINTER:: cell
 
@@ -709,16 +693,8 @@ MODULE moduleInput
                 IF (.NOT. found) CALL criticalError('energyThreshold not found for collision' // object, 'readInteractions')
                 CALL config%get(object // '.electron', electron, found)
                 IF (.NOT. found) CALL criticalError('electron not found for collision' // object, 'readInteractions')
-                CALL config%get(object // '.electronSecondary', electronSecondary, found)
+                CALL initBinaryIonization(interactionMatrix(ij)%collisions(k)%obj, &
-                IF (found) THEN
+                                          crossSecFilePath, energyThreshold, electron)
-                  CALL initBinaryIonization(interactionMatrix(ij)%collisions(k)%obj, &
-                                            crossSecFilePath, energyThreshold, electron, electronSecondary)
-
-                ELSE
-                  CALL initBinaryIonization(interactionMatrix(ij)%collisions(k)%obj, &
-                                            crossSecFilePath, energyThreshold, electron)
-
-                END IF
 
               CASE ('recombination')
                 !Electorn impact ionization
@@ -803,7 +779,7 @@ MODULE moduleInput
 
       TYPE(json_file), INTENT(inout):: config
       INTEGER:: i, s
-      CHARACTER(2):: iString, sString
+      CHARACTER(2):: istring, sString
       CHARACTER(:), ALLOCATABLE:: object, bType
       REAL(8):: Tw, cw !Wall temperature and specific heat
       !Neutral Properties
@@ -811,16 +787,16 @@ MODULE moduleInput
       REAL(8), DIMENSION(:), ALLOCATABLE:: v0
       REAL(8):: effTime
       REAL(8):: eThreshold !Energy threshold
-      INTEGER:: speciesID, electronSecondaryID
+      INTEGER:: speciesID
-      CHARACTER(:), ALLOCATABLE:: speciesName, crossSection, electronSecondary
+      CHARACTER(:), ALLOCATABLE:: speciesName, crossSection, yield
       LOGICAL:: found
       INTEGER:: nTypes
 
       CALL config%info('boundary', found, n_children = nBoundary)
       ALLOCATE(boundary(1:nBoundary))
       DO i = 1, nBoundary
-        WRITE(iString, '(i2)') i
+        WRITE(istring, '(i2)') i
-        object = 'boundary(' // TRIM(iString) // ')'
+        object = 'boundary(' // TRIM(istring) // ')'
 
         boundary(i)%n = i
         CALL config%get(object // '.name',            boundary(i)%name,            found)
@@ -829,77 +805,71 @@ MODULE moduleInput
         IF (nTypes /= nSpecies) CALL criticalError('Not enough boundary types defined in ' // object, 'readBoundary')
         ALLOCATE(boundary(i)%bTypes(1:nSpecies))
         DO s = 1, nSpecies
-          associate(bound => boundary(i)%bTypes(s)%obj)
+          WRITE(sString,'(i2)') s
-            WRITE(sString,'(i2)') s
+          object = 'boundary(' // TRIM(iString) // ').bTypes(' // TRIM(sString) // ')'
-            object = 'boundary(' // TRIM(iString) // ').bTypes(' // TRIM(sString) // ')'
+          CALL config%get(object // '.type', bType, found)
-            CALL config%get(object // '.type', bType, found)
+          SELECT CASE(bType)
-            SELECT CASE(bType)
+          CASE('reflection')
-            CASE('reflection')
+            ALLOCATE(boundaryReflection:: boundary(i)%bTypes(s)%obj)
-              ALLOCATE(boundaryReflection:: bound)
 
-            CASE('absorption')
+          CASE('absorption')
-              ALLOCATE(boundaryAbsorption:: bound)
+            ALLOCATE(boundaryAbsorption:: boundary(i)%bTypes(s)%obj)
 
-            CASE('transparent')
+          CASE('transparent')
-              ALLOCATE(boundaryTransparent:: bound)
+            ALLOCATE(boundaryTransparent:: boundary(i)%bTypes(s)%obj)
 
-            CASE('axis')
+          CASE('ionization')
-              ALLOCATE(boundaryAxis:: bound)
+            !Neutral parameters
+            CALL config%get(object // '.neutral.ion', speciesName, found)
+            IF (.NOT. found) CALL criticalError("missing parameter 'ion' for neutrals in ionization", 'readBoundary')
+            speciesID = speciesName2Index(speciesName)
+            CALL config%get(object // '.neutral.mass', m0, found)
+            IF (.NOT. found) THEN
+              m0 = species(s)%obj%m*m_ref
+            END IF
+            CALL config%get(object // '.neutral.density', n0, found)
+            IF (.NOT. found) CALL criticalError("missing parameter 'density' for neutrals in ionization", 'readBoundary')
+            CALL config%get(object // '.neutral.velocity', v0, found)
+            IF (.NOT. found) CALL criticalError("missing parameter 'velocity' for neutrals in ionization", 'readBoundary')
+            CALL config%get(object // '.neutral.temperature', T0, found)
+            IF (.NOT. found) CALL criticalError("missing parameter 'temperature' for neutrals in ionization", 'readBoundary')
 
-            CASE('wallTemperature')
+            CALL config%get(object // '.effectiveTime', effTime, found)
-              CALL config%get(object // '.temperature', Tw, found)
+            IF (.NOT. found) CALL criticalError("missing parameter 'effectiveTime' for ionization", 'readBoundary')
-              IF (.NOT. found) CALL criticalError("temperature not found for wallTemperature boundary type", 'readBoundary')
-              CALL config%get(object // '.specificHeat', cw, found)
-              IF (.NOT. found) CALL criticalError("specificHeat not found for wallTemperature boundary type", 'readBoundary')
 
-              CALL initWallTemperature(bound, Tw, cw)
+            CALL config%get(object // '.energyThreshold', eThreshold, found)
+            IF (.NOT. found) CALL criticalError("missing parameter 'eThreshold' in ionization", 'readBoundary')
 
-            CASE('ionization')
+            CALL config%get(object // '.crossSection', crossSection, found)
-              !Neutral parameters
+            IF (.NOT. found) CALL criticalError("missing parameter 'crossSection' for neutrals in ionization", 'readBoundary')
-              CALL config%get(object // '.neutral.ion', speciesName, found)
-              IF (.NOT. found) CALL criticalError("missing parameter 'ion' for neutrals in ionization", 'readBoundary')
-              speciesID = speciesName2Index(speciesName)
-              CALL config%get(object // '.neutral.mass', m0, found)
-              IF (.NOT. found) THEN
-                m0 = species(s)%obj%m*m_ref
-              END IF
-              CALL config%get(object // '.neutral.density', n0, found)
-              IF (.NOT. found) CALL criticalError("missing parameter 'density' for neutrals in ionization", 'readBoundary')
-              CALL config%get(object // '.neutral.velocity', v0, found)
-              IF (.NOT. found) CALL criticalError("missing parameter 'velocity' for neutrals in ionization", 'readBoundary')
-              CALL config%get(object // '.neutral.temperature', T0, found)
-              IF (.NOT. found) CALL criticalError("missing parameter 'temperature' for neutrals in ionization", 'readBoundary')
 
-              CALL config%get(object // '.effectiveTime', effTime, found)
+            CALL initIonization(boundary(i)%bTypes(s)%obj, species(s)%obj%m, m0, n0, v0, T0, &
-              IF (.NOT. found) CALL criticalError("missing parameter 'effectiveTime' for ionization", 'readBoundary')
+                                                           speciesID, effTime, crossSection, eThreshold)
 
-              CALL config%get(object // '.energyThreshold', eThreshold, found)
+          CASE('wallTemperature')
-              IF (.NOT. found) CALL criticalError("missing parameter 'eThreshold' in ionization", 'readBoundary')
+            CALL config%get(object // '.temperature', Tw, found)
+            IF (.NOT. found) CALL criticalError("temperature not found for wallTemperature boundary type", 'readBoundary')
+            CALL config%get(object // '.specificHeat', cw, found)
+            IF (.NOT. found) CALL criticalError("specificHeat not found for wallTemperature boundary type", 'readBoundary')
 
-              CALL config%get(object // '.crossSection', crossSection, found)
+            CALL initWallTemperature(boundary(i)%bTypes(s)%obj, Tw, cw)
-              IF (.NOT. found) CALL criticalError("missing parameter 'crossSection' for neutrals in ionization", 'readBoundary')
 
-              CALL config%get(object // '.electronSecondary', electronSecondary, found)
+          CASE('secondaryEmission')
-              electronSecondaryID = speciesName2Index(electronSecondary)
+            CALL config%get(object // '.yield', yield, found)
-              IF (found) THEN
+            IF (.NOT. found) CALL criticalError("missing parameter 'yield' for secondary emission", 'readBoundary')
-                CALL initIonization(bound, species(s)%obj%m, m0, n0, v0, T0, &
+            CALL config%get(object // '.electron', speciesName, found)
-                                                               speciesID, effTime, crossSection, eThreshold,electronSecondaryID)
+            IF (.NOT. found) CALL criticalError("missing parameter 'electron' for secondary emission", 'readBoundary')
+            speciesID = speciesName2Index(speciesName)
 
-              ELSE
+            CALL initSEE(boundary(i)%bTypes(s)%obj, yield, speciesID)
-                CALL initIonization(bound, species(s)%obj%m, m0, n0, v0, T0, &
-                                                               speciesID, effTime, crossSection, eThreshold)
 
-              END IF
+          CASE('axis')
+            ALLOCATE(boundaryAxis:: boundary(i)%bTypes(s)%obj)
 
-            case('outflowAdaptive')
+          CASE DEFAULT
-              allocate(boundaryOutflowAdaptive:: bound)
+            CALL criticalError('Boundary type ' // bType // ' undefined', 'readBoundary')
 
-            CASE DEFAULT
+          END SELECT
-              CALL criticalError('Boundary type ' // bType // ' undefined', 'readBoundary')
-
-            END SELECT
-
-          end associate
 
         END DO
 
@@ -916,7 +886,6 @@ MODULE moduleInput
       USE moduleMeshInputGmsh2, ONLY: initGmsh2
       USE moduleMeshInputVTU,   ONLY: initVTU
       USE moduleMeshInput0D,    ONLY: init0D
-      USE moduleMeshInputText,  ONLY: initText
       USE moduleMesh3DCart
       USE moduleMesh2DCyl
       USE moduleMesh2DCart
@@ -934,6 +903,7 @@ MODULE moduleInput
       LOGICAL:: found
       CHARACTER(:), ALLOCATABLE:: meshFormat, meshFile
       REAL(8):: volume
+      CHARACTER(:), ALLOCATABLE:: meshFileVTU !Temporary to test VTU OUTPUT
 
       object = 'geometry'
 
@@ -975,9 +945,9 @@ MODULE moduleInput
         !Read the 0D mesh
         CALL mesh%readMesh(pathMeshParticle)
 
-        !Get the volume
+        !Get the volumne
         CALL config%get(object // '.volume', volume, found)
-        !Rescale the volume
+        !Rescale the volumne
         IF (found) THEN
           mesh%cells(1)%obj%volume = mesh%cells(1)%obj%volume*volume / Vol_ref
           mesh%nodes(1)%obj%v     = mesh%cells(1)%obj%volume
@@ -1065,20 +1035,6 @@ MODULE moduleInput
 
           END IF
 
-        case ("text")
-          !Check if the geometry is right.
-          if (mesh%dimen /= 1) then
-            call criticalError("Text mesh is only allowed for 1D geometries", 'readGeometry')
-          
-          end if
-
-          !Read the mesh
-          call initText(mesh)
-          if (doubleMesh) then
-            call initText(meshColl)
-
-          end if
-
         CASE DEFAULT
           CALL criticalError('Mesh format ' // meshFormat // ' not defined.', 'readGeometry')
 
@@ -1125,13 +1081,13 @@ MODULE moduleInput
       TYPE(json_file), INTENT(inout):: config
       CHARACTER(:), ALLOCATABLE:: object
       LOGICAL:: found
-      CHARACTER(2):: iString
+      CHARACTER(2):: istring
       INTEGER:: i
       CHARACTER(:), ALLOCATABLE:: speciesName
       REAL(8), ALLOCATABLE, DIMENSION(:):: r
       REAL(8), ALLOCATABLE, DIMENSION(:):: v1, v2, v3
       INTEGER, ALLOCATABLE, DIMENSION(:):: points
-      REAL(8):: everyTimeStep
+      REAL(8):: timeStep
 
       CALL config%info('output.probes', found, n_children = nProbes)
 
@@ -1139,7 +1095,7 @@ MODULE moduleInput
 
       DO i = 1, nProbes
         WRITE(iString, '(I2)') i
-        object = 'output.probes(' // trim(iString) // ')'
+        object = 'output.probes(' // trim(istring) // ')'
 
         CALL config%get(object // '.species', speciesName, found)
         CALL config%get(object // '.position', r, found)
@@ -1147,14 +1103,16 @@ MODULE moduleInput
         CALL config%get(object // '.velocity_2', v2, found)
         CALL config%get(object // '.velocity_3', v3, found)
         CALL config%get(object // '.points', points, found)
-        CALL config%get(object // '.timeStep', everyTimeStep, found)
+        CALL config%get(object // '.timeStep', timeStep, found)
 
         IF (ANY(points < 2)) CALL criticalError("Number of points in probe " // iString // " incorrect", 'readProbes')
 
-        CALL probe(i)%init(i, speciesName, r, v1, v2, v3, points, everyTimeStep)
+        CALL probe(i)%init(i, speciesName, r, v1, v2, v3, points, timeStep)
 
       END DO
 
+      CALL resetProbes(tInitial)
+
     END SUBROUTINE readProbes
 
     SUBROUTINE readEMBoundary(config)
@@ -1162,6 +1120,7 @@ MODULE moduleInput
       USE moduleOutput
       USE moduleErrors
       USE moduleEM
+      USE moduleRefParam
       USE moduleSpecies
       USE json_module
       IMPLICIT NONE
@@ -1169,72 +1128,34 @@ MODULE moduleInput
       TYPE(json_file), INTENT(inout):: config
       CHARACTER(:), ALLOCATABLE:: object
       LOGICAL:: found
-      CHARACTER(:), ALLOCATABLE:: typeEM
+      CHARACTER(2):: istring
-      REAL(8):: potential
+      INTEGER:: i, e, s
-      INTEGER:: physicalSurface
-      CHARACTER(:), ALLOCATABLE:: temporalProfile, temporalProfilePath
-      INTEGER:: b, s, n, ni
-      CHARACTER(2):: bString
       INTEGER:: info
       EXTERNAL:: dgetrf
 
       CALL config%info('boundaryEM', found, n_children = nBoundaryEM)
       
-      IF (found) THEN
+      IF (found) ALLOCATE(boundEM(1:nBoundaryEM))
-        ALLOCATE(boundaryEM(1:nBoundaryEM))
-
-      END IF
       
-      DO b = 1, nBoundaryEM
+      DO i = 1, nBoundaryEM
-        WRITE(bString, '(I2)') b
+        WRITE(istring, '(I2)') i
-        object = 'boundaryEM(' // TRIM(bString) // ')'
+        object = 'boundaryEM(' // trim(istring) // ')'
 
-        CALL config%get(object // '.type', typeEM, found)
+        CALL config%get(object // '.type', boundEM(i)%typeEM, found)
 
-        SELECT CASE(typeEM)
+        SELECT CASE(boundEM(i)%typeEM)
         CASE ("dirichlet")
-          CALL config%get(object // '.potential', potential, found)
+          CALL config%get(object // '.potential', boundEM(i)%potential, found)
-          IF (.NOT. found) THEN
+          IF (.NOT. found) &
             CALL criticalError('Required parameter "potential" for Dirichlet boundary condition not found', 'readEMBoundary')
+          boundEM(i)%potential = boundEM(i)%potential/Volt_ref
 
-          END IF
+          CALL config%get(object // '.physicalSurface', boundEM(i)%physicalSurface, found)
-
+          IF (.NOT. found) &
-          CALL config%get(object // '.physicalSurface', physicalSurface, found)
+            CALL criticalError('Required parameter "physicalSurface" for Dirichlet boundary condition not found', 'readEMBoundary')
-          IF (.NOT. found) THEN
-            CALL criticalError('Required parameter "physicalSurface" for Dirichlet boundary condition not found', &
-                               'readEMBoundary')
-
-          END IF
-
-          CALL initDirichlet(boundaryEM(b)%obj, physicalSurface, potential)
-
-        CASE ("dirichletTime")
-          CALL config%get(object // '.potential', potential, found)
-          IF (.NOT. found) THEN
-            CALL criticalError('Required parameter "potential" for Dirichlet Time boundary condition not found', &
-                               'readEMBoundary')
-
-          END IF
-
-          CALL config%get(object // '.temporalProfile', temporalProfile, found)
-          IF (.NOT. found) THEN
-            CALL criticalError('Required parameter "temporalProfile" for Dirichlet Time boundary condition not found', &
-                               'readEMBoundary')
-
-          END IF
-          temporalProfilePath = path // temporalProfile
-
-          CALL config%get(object // '.physicalSurface', physicalSurface, found)
-          IF (.NOT. found) THEN
-            CALL criticalError('Required parameter "physicalSurface" for Dirichlet Time boundary condition not found', &
-                               'readEMBoundary')
-
-          END IF
-
-          CALL initDirichletTime(boundaryEM(b)%obj, physicalSurface, potential, temporalProfilePath)
 
         CASE DEFAULT
-          CALL criticalError('Boundary type ' // typeEM // ' not yet supported', 'readEMBoundary')
+          CALL criticalError('Boundary type ' // boundEM(i)%typeEM // ' not yet supported', 'readEMBoundary')
 
         END SELECT
 
@@ -1253,28 +1174,18 @@ MODULE moduleInput
 
       END DO
 
-      ! Modify K matrix due to boundary conditions
+      IF (ALLOCATED(boundEM)) THEN
-      DO b = 1, nBoundaryEM
+        DO e = 1, mesh%numEdges
-        SELECT TYPE(boundary => boundaryEM(b)%obj)
+          IF (ANY(mesh%edges(e)%obj%physicalSurface == boundEM(:)%physicalSurface)) THEN
-        TYPE IS(boundaryEMDirichlet)
+            DO i = 1, nBoundaryEM
-          DO n = 1, boundary%nNodes
+              IF (mesh%edges(e)%obj%physicalSurface == boundEM(i)%physicalSurface) THEN
-            ni = boundary%nodes(n)%obj%n
+                CALL boundEM(i)%apply(mesh%edges(e)%obj)
-            mesh%K(ni, :)  = 0.D0
-            mesh%K(ni, ni) = 1.D0
 
-          END DO
+              END IF
-
+            END DO
-        TYPE IS(boundaryEMDirichletTime)
+          END IF
-          DO n = 1, boundary%nNodes
+        END DO
-            ni = boundary%nodes(n)%obj%n
+      END IF
-            mesh%K(ni, :)  = 0.D0
-            mesh%K(ni, ni) = 1.D0
-
-          END DO
-
-        END SELECT
-
-      END DO
 
       !Compute the PLU factorization of K once boundary conditions have been read
       CALL dgetrf(mesh%numNodes, mesh%numNodes, mesh%K, mesh%numNodes, mesh%IPIV, info)
@@ -1295,25 +1206,24 @@ MODULE moduleInput
 
       TYPE(json_file), INTENT(inout):: config
       INTEGER:: i
-      CHARACTER(2):: iString
+      CHARACTER(2):: istring
       CHARACTER(:), ALLOCATABLE:: object
       LOGICAL:: found
       CHARACTER(:), ALLOCATABLE:: speciesName
       CHARACTER(:), ALLOCATABLE:: name
       REAL(8):: v
-      REAL(8), ALLOCATABLE:: temperature(:), normal(:)
+      REAL(8), ALLOCATABLE:: T(:), normal(:)
       REAL(8):: flow
       CHARACTER(:), ALLOCATABLE:: units
       INTEGER:: physicalSurface
-      INTEGER:: particlesPerEdge
       INTEGER:: sp
 
       CALL config%info('inject', found, n_children = nInject)
       ALLOCATE(inject(1:nInject))
       nPartInj = 0
       DO i = 1, nInject
-        WRITE(iString, '(i2)') i
+        WRITE(istring, '(i2)') i
-        object = 'inject(' // trim(iString) // ')'
+        object = 'inject(' // trim(istring) // ')'
 
         !Find species
         CALL config%get(object // '.species', speciesName, found)
@@ -1321,7 +1231,7 @@ MODULE moduleInput
 
         CALL config%get(object // '.name',            name,            found)
         CALL config%get(object // '.v',               v,               found)
-        CALL config%get(object // '.T',               temperature,     found)
+        CALL config%get(object // '.T',               T,               found)
         CALL config%get(object // '.n',               normal,          found)
         IF (.NOT. found) THEN
           ALLOCATE(normal(1:3))
@@ -1330,10 +1240,8 @@ MODULE moduleInput
         CALL config%get(object // '.flow',            flow,            found)
         CALL config%get(object // '.units',           units,           found)
         CALL config%get(object // '.physicalSurface', physicalSurface, found)
-        particlesPerEdge = 0
-        CALL config%get(object // '.particlesPerEdge', particlesPerEdge, found)
 
-        CALL inject(i)%init(i, v, normal, temperature, flow, units, sp, physicalSurface, particlesPerEdge)
+        CALL inject(i)%init(i, v, normal, T, flow, units, sp, physicalSurface)
 
         CALL readVelDistr(config, inject(i), object)
         
@@ -1352,7 +1260,6 @@ MODULE moduleInput
       USE moduleCaseParam, ONLY: tauMin
       USE moduleMesh, ONLY: mesh
       USE moduleSpecies, ONLY: nSpecies
-      USE moduleRefParam, ONLY: ti_ref
       IMPLICIT NONE
 
       TYPE(json_file), INTENT(inout):: config
@@ -1366,11 +1273,8 @@ MODULE moduleInput
         CALL config%get('average.startTime', tStart, found)
 
         IF (found) THEN
-          tAverageStart = INT(tStart / ti_ref / tauMin)
+          tAverageStart = INT(tStart / tauMin)
 
-        ELSE
-          tAverageStart = 0
-          
         END IF
 
         ALLOCATE(averageScheme(1:mesh%numNodes))
@@ -1397,28 +1301,28 @@ MODULE moduleInput
       TYPE(injectGeneric), INTENT(inout):: inj
       CHARACTER(:), ALLOCATABLE, INTENT(in):: object
       INTEGER:: i
-      CHARACTER(2):: iString
+      CHARACTER(2):: istring
       CHARACTER(:), ALLOCATABLE:: fvType
       LOGICAL:: found
-      REAL(8):: v, temperature, m
+      REAL(8):: v, T, m
 
       !Reads species mass
       m = inj%species%m
       !Reads distribution functions for velocity
       DO i = 1, 3
-        WRITE(iString, '(i2)') i
+        WRITE(istring, '(i2)') i
-        CALL config%get(object // '.velDist('// TRIM(iString) //')', fvType, found)
+        CALL config%get(object // '.velDist('// TRIM(istring) //')', fvType, found)
-        IF (.NOT. found) CALL criticalError("No velocity distribution in direction " // iString // &
+        IF (.NOT. found) CALL criticalError("No velocity distribution in direction " // istring // &
                                             " found for " // object, 'readVelDistr')
 
         SELECT CASE(fvType)
         CASE ("Maxwellian")
-          temperature = inj%temperature(i)
+          T = inj%T(i)
-          CALL initVelDistMaxwellian(inj%v(i)%obj, temperature, m)
+          CALL initVelDistMaxwellian(inj%v(i)%obj, t, m)
 
         CASE ("Half-Maxwellian")
-          temperature = inj%temperature(i)
+          T = inj%T(i)
-          CALL initVelDistHalfMaxwellian(inj%v(i)%obj, temperature, m)
+          CALL initVelDistHalfMaxwellian(inj%v(i)%obj, t, m)
 
         CASE ("Delta")
           v = inj%vMod*inj%n(i)
@@ -1459,37 +1363,5 @@ MODULE moduleInput
 
     END SUBROUTINE readParallel
 
-    SUBROUTINE initOutput(inputFile)
-      USE moduleRefParam
-      USE moduleMesh,   ONLY: mesh, doubleMesh, pathMeshParticle, pathMeshColl
-      USE moduleOutput, ONLY: path, folder
-      IMPLICIT NONE
-
-      CHARACTER(:), ALLOCATABLE, INTENT(in):: inputFile
-      INTEGER:: fileReference = 30
-      !If everything is correct, creates the output folder
-      CALL EXECUTE_COMMAND_LINE('mkdir ' // path // folder )
-      !Copies input file to output folder
-      CALL EXECUTE_COMMAND_LINE('cp ' // inputFile // ' ' // path // folder)
-      !Copies particle mesh
-      IF (mesh%dimen > 0) THEN
-        CALL EXECUTE_COMMAND_LINE('cp ' // pathMeshParticle // ' ' // path // folder)
-        IF (doubleMesh) THEN
-          CALL EXECUTE_COMMAND_LINE('cp ' // pathMeshColl // ' ' // path // folder)
-
-        END IF
-
-      END IF
-
-      ! Write commit of fpakc
-      CALL SYSTEM('git rev-parse HEAD > ' // path // folder // '/' // 'fpakc_commit.txt')
-
-      ! Write file with reference values
-      OPEN (fileReference, file=path // folder // '/' // 'reference.txt')
-      WRITE(fileReference, "(7(1X,A20))")      'L_ref', 'v_ref', 'ti_ref', 'Vol_ref', 'EF_ref', 'Volt_ref', 'B_ref'
-      WRITE(fileReference, "(7(1X,ES20.6E3))")  L_ref,   v_ref,   ti_ref,   Vol_ref,   EF_ref,   Volt_ref,   B_ref
-      CLOSE(fileReference)
-
-    END SUBROUTINE initOutput
 
 END MODULE moduleInput

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

@@ -122,8 +122,6 @@ MODULE moduleMesh1DCart
 
       self%x = r1(1)
 
-      self%surface = 1.D0
-
       self%normal = (/ 1.D0, 0.D0, 0.D0 /)
 
       !Boundary index

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

@@ -122,8 +122,6 @@ MODULE moduleMesh1DRad
 
       self%r = r1(1)
 
-      self%surface = 1.D0
-
       self%normal = (/ 1.D0, 0.D0, 0.D0 /)
 
       !Boundary index

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

@@ -163,7 +163,7 @@ MODULE moduleMesh2DCart
       r2 = self%n2%getCoordinates()
       self%x  = (/r1(1), r2(1)/)
       self%y  = (/r1(2), r2(2)/)
-      self%surface = SQRT((self%x(2) - self%x(1))**2 + (self%y(2) - self%y(1))**2)
+      self%weight = 1.D0
       !Normal vector
       self%normal = (/ -(self%y(2)-self%y(1)), &
                          self%x(2)-self%x(1) , &
@@ -318,8 +318,6 @@ MODULE moduleMesh2DCart
       INTEGER, INTENT(in):: nNodes
       REAL(8):: fPsi(1:nNodes)
 
-      fPsi = 0.D0
-
       fPsi = (/ (1.D0 - Xi(1)) * (1.D0 - Xi(2)), &
                 (1.D0 + Xi(1)) * (1.D0 - Xi(2)), &
                 (1.D0 + Xi(1)) * (1.D0 + Xi(2)), &
@@ -494,36 +492,34 @@ MODULE moduleMesh2DCart
 
     END FUNCTION insideQuad
 
-    !Transform physical coordinates to element coordinates with a Taylor series
+    !Transform physical coordinates to element coordinates
     PURE FUNCTION phy2logQuad(self,r) RESULT(Xi) 
       IMPLICIT NONE
 
       CLASS(meshCell2DCartQuad), INTENT(in):: self
       REAL(8), INTENT(in):: r(1:3)
       REAL(8):: Xi(1:3)
-      REAL(8):: Xi0(1:3), detJ, pDerInv(1:2,1:2), deltaR(1:2), x0(1:2)
+      REAL(8):: XiO(1:3), detJ, invJ(1:3,1:3), f(1:3)
       REAL(8):: dPsi(1:3,1:4), fPsi(1:4)
       REAL(8):: pDer(1:3, 1:3)
       REAL(8):: conv
 
       !Iterative newton method to transform coordinates
-      conv  = 1.D0
+      conv = 1.D0
-      Xi0   = 0.D0
+      XiO  = 0.D0
-      Xi(3) = 0.D0
 
       DO WHILE(conv > 1.D-4)
-        fPsi           = self%fPsi(Xi0, 4)
+        dPsi = self%dPsi(XiO, 4)
-        x0(1)          = dot_product(fPsi, self%x)
+        pDer = self%partialDer(4, dPsi)
-        x0(2)          = dot_product(fPsi, self%y)
+        detJ = self%detJac(pDer)
-        deltaR         = r(1:2) - x0
+        invJ = self%invJac(pDer)
-        dPsi           = self%dPsi(Xi0, 4)
+        fPsi = self%fPsi(XiO, 4)
-        pDer           = self%partialDer(4, dPsi)
+        f = (/ DOT_PRODUCT(fPsi,self%x), &
-        detJ           = self%detJac(pDer)
+               DOT_PRODUCT(fPsi,self%y), &
-        pDerInv(1,1:2) = (/  pDer(2,2), -pDer(1,2) /)
+               0.D0 /) - r
-        pDerInv(2,1:2) = (/ -pDer(2,1),  pDer(1,1) /)
+        Xi = XiO - MATMUL(invJ, f)/detJ
-        Xi(1:2)        = Xi0(1:2) + MATMUL(pDerInv, deltaR)/detJ
+        conv = MAXVAL(DABS(Xi-XiO),1)
-        conv           = MAXVAL(DABS(Xi(1:2)-Xi0(1:2)),1)
+        XiO = Xi
-        Xi0(1:2)       = Xi(1:2)
 
       END DO
 
@@ -573,7 +569,6 @@ MODULE moduleMesh2DCart
       pDer = self%partialDer(4, dPsi)
       detJ = self%detJac(pDer)
       fPsi = self%fPsi(Xi, 4)
-
       !Compute total volume of the cell
       self%volume  = detJ*4.D0
       !Compute volume per node
@@ -680,8 +675,8 @@ MODULE moduleMesh2DCart
 
       dPsi = 0.D0
 
-      dPsi(1,1:3) = (/ -1.D0, 1.D0, 0.D0 /)
+      dPsi(1,:) = (/ -1.D0, 1.D0, 0.D0 /)
-      dPsi(2,1:3) = (/ -1.D0, 0.D0, 1.D0 /)
+      dPsi(2,:) = (/ -1.D0, 0.D0, 1.D0 /)
 
     END FUNCTION dPsiTria
 
@@ -767,7 +762,6 @@ MODULE moduleMesh2DCart
         pDer   = self%partialDer(3, dPsi)
         detJ   = self%detJac(pDer)
         invJ   = self%invJac(pDer)
-
         localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*wTria(l)/detJ
 
       END DO
@@ -826,19 +820,19 @@ MODULE moduleMesh2DCart
       CLASS(meshCell2DCartTria), INTENT(in):: self
       REAL(8), INTENT(in):: r(1:3)
       REAL(8):: Xi(1:3)
-      REAL(8):: detJ, pDerInv(1:2,1:2), deltaR(1:2)
+      REAL(8):: deltaR(1:3)
-      REAL(8):: dPsi(1:3,1:4)
+      REAL(8):: dPsi(1:3, 1:3)
       REAL(8):: pDer(1:3, 1:3)
+      REAL(8):: invJ(1:3,  1:3), detJ
 
       !Direct method to convert coordinates
-      Xi(3)          = 0.D0
+      Xi     = 0.D0
-      deltaR         = (/ r(1) - self%x(1), r(2) - self%y(1) /)
+      deltaR = (/ r(1) - self%x(1), r(2) - self%y(1), 0.D0 /)
-      dPsi           = self%dPsi(Xi, 3)
+      dPsi   = self%dPsi(Xi, 3)
-      pDer           = self%partialDer(3, dPsi)
+      pDer   = self%partialDer(3, dPsi)
-      detJ           = self%detJac(pDer)
+      invJ   = self%invJac(pDer)
-      pDerInv(1,1:2) = (/  pDer(2,2), -pDer(1,2) /)
+      detJ   = self%detJac(pDer)
-      pDerInv(2,1:2) = (/ -pDer(2,1),  pDer(1,1) /)
+      Xi     = MATMUL(invJ,deltaR)/detJ
-      Xi(1:2)        = MATMUL(pDerInv,deltaR)/detJ
 
     END FUNCTION phy2logTria
 
@@ -913,8 +907,8 @@ MODULE moduleMesh2DCart
 
       invJ = 0.D0
 
-      invJ(1, 1:2) = (/  pDer(2,2), -pDer(2,1) /)
+      invJ(1, 1:2) = (/  pDer(2,2), -pDer(1,2) /)
-      invJ(2, 1:2) = (/ -pDer(1,2),  pDer(1,1) /)
+      invJ(2, 1:2) = (/ -pDer(2,1),  pDer(1,1) /)
       invJ(3, 3)   = 1.D0
 
     END FUNCTION invJ2DCart

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

@@ -144,7 +144,6 @@ MODULE moduleMesh2DCyl
       USE moduleSpecies
       USE moduleBoundary
       USE moduleErrors
-      USE moduleConstParam, ONLY: PI
       IMPLICIT NONE
 
       CLASS(meshEdge2DCyl), INTENT(out):: self
@@ -164,15 +163,7 @@ MODULE moduleMesh2DCyl
       r2 = self%n2%getCoordinates()
       self%z  = (/r1(1), r2(1)/)
       self%r  = (/r1(2), r2(2)/)
-      !Edge surface
+      self%weight = r2(2)**2 - r1(2)**2
-      IF (self%z(2) /= self%z(1)) THEN
-        self%surface = ABS(self%r(2) + self%r(1))*ABS(self%z(2) - self%z(1))
-
-      ELSE
-        self%surface = ABS(self%r(2)**2 - self%r(1)**2)
-
-      END IF
-      self%surface = self%surface * PI
       !Normal vector
       self%normal = (/ -(self%r(2)-self%r(1)), &
                          self%z(2)-self%z(1) , &
@@ -232,13 +223,21 @@ MODULE moduleMesh2DCyl
       CLASS(meshEdge2DCyl), INTENT(in):: self
       REAL(8):: rnd
       REAL(8):: r(1:3)
-      REAL(8):: p1(1:2), p2(1:2)
+      REAL(8):: dr, dz
 
       rnd = random()
+      dr  = self%r(2) - self%r(1)
+      dz  = self%z(2) - self%z(1)
+      IF (dr /= 0.D0) THEN
+        r(2) = dr * DSQRT(rnd) + self%r(1)
+        r(1) = dz * (r(2) - self%r(1))/dr + self%z(1)
+
+      ELSE
+        r(2) = self%r(1)
+        r(1) = dz * rnd + self%z(1)
+
+      END IF
 
-      p1 = (/self%z(1), self%r(1) /)
-      p2 = (/self%z(2), self%r(2) /)
-      r(1:2) = (1.D0 - rnd)*p1 + rnd*p2
       r(3)   =  0.D0
 
     END FUNCTION randPosEdge
@@ -247,6 +246,7 @@ MODULE moduleMesh2DCyl
     !QUAD FUNCTIONS
     !Init element
     SUBROUTINE initCellQuad2DCyl(self, n, p, nodes)
+      USE moduleRefParam
       IMPLICIT NONE
 
       CLASS(meshCell2DCylQuad), INTENT(out):: self
@@ -326,8 +326,6 @@ MODULE moduleMesh2DCyl
       INTEGER, INTENT(in):: nNodes
       REAL(8):: fPsi(1:nNodes)
 
-      fPsi = 0.D0
-
       fPsi = (/ (1.D0 - Xi(1)) * (1.D0 - Xi(2)), &
                 (1.D0 + Xi(1)) * (1.D0 - Xi(2)), &
                 (1.D0 + Xi(1)) * (1.D0 + Xi(2)), &
@@ -498,7 +496,7 @@ MODULE moduleMesh2DCyl
 
     END FUNCTION elemFQuad
 
-    !Check if Xi is inside the element
+    !Checks if Xi is inside the element
     PURE FUNCTION insideQuad(Xi) RESULT(ins)
       IMPLICIT NONE
 
@@ -510,36 +508,34 @@ MODULE moduleMesh2DCyl
 
     END FUNCTION insideQuad
 
-    !Transform physical coordinates to element coordinates with a Taylor series
+    !Transform physical coordinates to element coordinates
     PURE FUNCTION phy2logQuad(self,r) RESULT(Xi) 
       IMPLICIT NONE
 
       CLASS(meshCell2DCylQuad), INTENT(in):: self
       REAL(8), INTENT(in):: r(1:3)
       REAL(8):: Xi(1:3)
-      REAL(8):: Xi0(1:3), detJ, pDerInv(1:2,1:2), deltaR(1:2), x0(1:2)
+      REAL(8):: XiO(1:3), detJ, invJ(1:3,1:3), f(1:3)
       REAL(8):: dPsi(1:3,1:4), fPsi(1:4)
       REAL(8):: pDer(1:3, 1:3)
       REAL(8):: conv
 
       !Iterative newton method to transform coordinates
-      conv  = 1.D0
+      conv = 1.D0
-      Xi0   = 0.D0
+      XiO  = 0.D0
-      Xi(3) = 0.D0
 
       DO WHILE(conv > 1.D-4)
-        fPsi           = self%fPsi(Xi0, 4)
+        dPsi = self%dPsi(XiO, 4)
-        x0(1)          = dot_product(fPsi, self%z)
+        pDer = self%partialDer(4, dPsi)
-        x0(2)          = dot_product(fPsi, self%r)
+        detJ = self%detJac(pDer)
-        deltaR         = r(1:2) - x0
+        invJ = self%invJac(pDer)
-        dPsi           = self%dPsi(Xi0, 4)
+        fPsi = self%fPsi(XiO, 4)
-        pDer           = self%partialDer(4, dPsi)
+        f = (/ DOT_PRODUCT(fPsi,self%z), &
-        detJ           = self%detJac(pDer)
+               DOT_PRODUCT(fPsi,self%r), &
-        pDerInv(1,1:2) = (/  pDer(2,2), -pDer(1,2) /)
+               0.D0 /) - r
-        pDerInv(2,1:2) = (/ -pDer(2,1),  pDer(1,1) /)
+        Xi = XiO - MATMUL(invJ, f)/detJ
-        Xi(1:2)        = Xi0(1:2) + MATMUL(pDerInv, deltaR)/detJ
+        conv = MAXVAL(DABS(Xi-XiO),1)
-        conv           = MAXVAL(DABS(Xi(1:2)-Xi0(1:2)),1)
+        XiO = Xi
-        Xi0(1:2)       = Xi(1:2)
 
       END DO
 
@@ -557,7 +553,7 @@ MODULE moduleMesh2DCyl
 
       XiArray = (/ -Xi(2), Xi(1), Xi(2), -Xi(1) /)
       nextInt = MAXLOC(XiArray,1)
-      !Select the higher value of directions and searches in that direction
+      !Selects the higher value of directions and searches in that direction
       NULLIFY(neighbourElement)
       SELECT CASE (nextInt)
       CASE (1)
@@ -585,7 +581,6 @@ MODULE moduleMesh2DCyl
       REAL(8):: dPsi(1:3, 1:4), pDer(1:3, 1:3)
 
       self%volume  = 0.D0
-
       !2D 1 point Gauss Quad Integral
       Xi = 0.D0
       dPsi = self%dPsi(Xi, 4)
@@ -594,18 +589,18 @@ MODULE moduleMesh2DCyl
       fPsi = self%fPsi(Xi, 4)
       r = DOT_PRODUCT(fPsi,self%r)
       !Computes total volume of the cell
-      self%volume  = r*detJ*PI8 !2*pi * 4 (weight of 1 point 2D-Gaussian integral)
+      self%volume  = r*detJ*PI8 !4*2*pi
-      !Computes volume per node. Change the radius point to calculate the area to improve accuracy near the axis.
+      !Computes volume per node
-      Xi = (/-5.D-1, -0.25D0, 0.D0/)
+      Xi = (/-5.D-1, -5.D-1, 0.D0/)
       r = self%gatherF(Xi, 4, self%r)
       self%n1%v = self%n1%v + fPsi(1)*r*detJ*PI8
-      Xi = (/ 5.D-1, -0.25D0, 0.D0/)
+      Xi = (/ 5.D-1, -5.D-1, 0.D0/)
       r = self%gatherF(Xi, 4, self%r)
       self%n2%v = self%n2%v + fPsi(2)*r*detJ*PI8
-      Xi = (/ 5.D-1,  0.75D0, 0.D0/)
+      Xi = (/ 5.D-1,  5.D-1, 0.D0/)
       r = self%gatherF(Xi, 4, self%r)
       self%n3%v = self%n3%v + fPsi(3)*r*detJ*PI8
-      Xi = (/-5.D-1,  0.75D0, 0.D0/)
+      Xi = (/-5.D-1,  5.D-1, 0.D0/)
       r = self%gatherF(Xi, 4, self%r)
       self%n4%v = self%n4%v + fPsi(4)*r*detJ*PI8
 
@@ -707,8 +702,8 @@ MODULE moduleMesh2DCyl
 
       dPsi = 0.D0
 
-      dPsi(1,1:3) = (/ -1.D0, 1.D0, 0.D0 /)
+      dPsi(1,:) = (/ -1.D0, 1.D0, 0.D0 /)
-      dPsi(2,1:3) = (/ -1.D0, 0.D0, 1.D0 /)
+      dPsi(2,:) = (/ -1.D0, 0.D0, 1.D0 /)
 
     END FUNCTION dPsiTria
 
@@ -860,19 +855,19 @@ MODULE moduleMesh2DCyl
       CLASS(meshCell2DCylTria), INTENT(in):: self
       REAL(8), INTENT(in):: r(1:3)
       REAL(8):: Xi(1:3)
-      REAL(8):: detJ, pDerInv(1:2,1:2), deltaR(1:2)
+      REAL(8):: deltaR(1:3)
-      REAL(8):: dPsi(1:3,1:4)
+      REAL(8):: dPsi(1:3, 1:3)
       REAL(8):: pDer(1:3, 1:3)
+      REAL(8):: invJ(1:3,  1:3), detJ
 
       !Direct method to convert coordinates
-      Xi(3)          = 0.D0
+      Xi     = 0.D0
-      deltaR         = (/ r(1) - self%z(1), r(2) - self%r(1) /)
+      deltaR = (/ r(1) - self%z(1), r(2) - self%r(1), 0.D0 /)
-      dPsi           = self%dPsi(Xi, 3)
+      dPsi   = self%dPsi(Xi, 3)
-      pDer           = self%partialDer(3, dPsi)
+      pDer   = self%partialDer(3, dPsi)
-      detJ           = self%detJac(pDer)
+      invJ   = self%invJac(pDer)
-      pDerInv(1,1:2) = (/  pDer(2,2), -pDer(1,2) /)
+      detJ   = self%detJac(pDer)
-      pDerInv(2,1:2) = (/ -pDer(2,1),  pDer(1,1) /)
+      Xi     = MATMUL(invJ,deltaR)/detJ
-      Xi(1:2)        = MATMUL(pDerInv,deltaR)/detJ
 
     END FUNCTION phy2logTria
 
@@ -950,8 +945,8 @@ MODULE moduleMesh2DCyl
 
       invJ = 0.D0
 
-      invJ(1, 1:2) = (/  pDer(2,2), -pDer(2,1) /)
+      invJ(1, 1:2) = (/  pDer(2,2), -pDer(1,2) /)
-      invJ(2, 1:2) = (/ -pDer(1,2),  pDer(1,1) /)
+      invJ(2, 1:2) = (/ -pDer(2,1),  pDer(1,1) /)
       invJ(3, 3)   = 1.D0
 
     END FUNCTION invJ2DCyl

src/modules/mesh/3DCart/moduleMesh3DCart.f90

@@ -109,7 +109,6 @@ MODULE moduleMesh3DCart
       USE moduleBoundary
       USE moduleErrors
       USE moduleMath
-      USE moduleRefParam, ONLY: L_ref
       IMPLICIT NONE
 
       CLASS(meshEdge3DCartTria), INTENT(out):: self
@@ -143,8 +142,6 @@ MODULE moduleMesh3DCart
       self%normal = crossProduct(vec1, vec2)
       self%normal = normalize(self%normal)
 
-      self%surface = 1.D0/L_ref**2 !TODO: FIX THIS WHEN MOVING TO 3D
-
       !Boundary index
       self%boundary => boundary(bt)
       ALLOCATE(self%fBoundary(1:nSpecies))

src/modules/mesh/inout/0D/moduleMeshOutput0D.f90

@@ -1,22 +1,22 @@
 MODULE moduleMeshOutput0D
 
   CONTAINS
-    SUBROUTINE printOutput0D(self)
+    SUBROUTINE printOutput0D(self, t)
       USE moduleMesh
       USE moduleRefParam
       USE moduleSpecies
       USE moduleOutput
-      USE moduleCaseParam, ONLY: timeStep
       IMPLICIT NONE
 
       CLASS(meshParticles), INTENT(in):: self
+      INTEGER, INTENT(in):: t
       INTEGER:: i
       TYPE(outputFormat):: output
       CHARACTER(:), ALLOCATABLE:: fileName
 
       DO i = 1, nSpecies
         fileName='OUTPUT_' // species(i)%obj%name // '.dat'
-        IF (timeStep == 0) THEN
+        IF (t == 0) THEN
           OPEN(20, file = path // folder // '/' // fileName, action = 'write')
           WRITE(20, "(A1, 14X, A5, A20, 40X, A20, 2(A20))") "#","t (s)","density (m^-3)", "velocity (m/s)", &
                                                                         "pressure (Pa)", "temperature (K)"
@@ -27,17 +27,14 @@ MODULE moduleMeshOutput0D
         
         OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write')
         CALL calculateOutput(self%nodes(1)%obj%output(i), output, self%nodes(1)%obj%v, species(i)%obj)
-        WRITE(20, "(7(ES20.6E3))") REAL(timeStep)*tauMin*ti_ref, output%density,    &
+        WRITE(20, "(7(ES20.6E3))") REAL(t)*tauMin*ti_ref, output%density, output%velocity, output%pressure, output%temperature
-                                                                 output%velocity,   &
-                                                                 output%pressure,   &
-                                                                 output%temperature
         CLOSE(20)
 
       END DO
 
     END SUBROUTINE printOutput0D
 
-    SUBROUTINE printColl0D(self)
+    SUBROUTINE printColl0D(self, t)
       USE moduleMesh
       USE moduleRefParam
       USE moduleCaseParam
@@ -46,11 +43,12 @@ MODULE moduleMeshOutput0D
       IMPLICIT NONE
 
       CLASS(meshGeneric), INTENT(in):: self
+      INTEGER, INTENT(in):: t
       CHARACTER(:), ALLOCATABLE:: fileName
       INTEGER:: k
 
       fileName='OUTPUT_Collisions.dat'
-      IF (timeStep == tInitial) THEN
+      IF (t == tInitial) THEN
         OPEN(20, file = path // folder // '/' // fileName, action = 'write')
         WRITE(20, "(A1, 14X, A5, A20)") "#","t (s)","collisions"
         WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
@@ -59,12 +57,12 @@ MODULE moduleMeshOutput0D
       END IF
       
       OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write')
-      WRITE(20, "(ES20.6E3, 10I20)") REAL(timeStep)*tauMin*ti_ref, (self%cells(1)%obj%tallyColl(k)%tally, k=1,nCollPairs)
+      WRITE(20, "(ES20.6E3, 10I20)") REAL(t)*tauMin*ti_ref, (self%cells(1)%obj%tallyColl(k)%tally, k=1,nCollPairs)
       CLOSE(20)
 
     END SUBROUTINE printColl0D
 
-    SUBROUTINE printEM0D(self)
+    SUBROUTINE printEM0D(self, t)
       USE moduleMesh
       USE moduleRefParam
       USE moduleCaseParam
@@ -72,6 +70,7 @@ MODULE moduleMeshOutput0D
       IMPLICIT NONE
 
       CLASS(meshParticles), INTENT(in):: self
+      INTEGER, INTENT(in):: t
 
     END SUBROUTINE printEM0D
 

src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90

@@ -108,7 +108,6 @@ MODULE moduleMeshInputGmsh2
       READ(10, *) totalNumElem
 
       !Count edges and volume elements
-      numEdges = 0
       SELECT TYPE(self)
       TYPE IS(meshParticles)
         self%numEdges = 0
@@ -329,7 +328,7 @@ MODULE moduleMeshInputGmsh2
 
       DO i = 1, numNodes
         !Reads the density
-        READ(10, *) e, density(i)
+        READ(10, *), e, density(i)
 
       END DO
 
@@ -340,7 +339,7 @@ MODULE moduleMeshInputGmsh2
 
       DO i = 1, numNodes
         !Reads the velocity
-        READ(10, *) e, velocity(i, 1:3)
+        READ(10, *), e, velocity(i, 1:3)
 
       END DO
 

src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90

@@ -80,50 +80,50 @@ MODULE moduleMeshOutputGmsh2
     END SUBROUTINE writeGmsh2FooterElementData
 
     !Prints the scattered properties of particles into the nodes
-    SUBROUTINE printOutputGmsh2(self)
+    SUBROUTINE printOutputGmsh2(self, t)
       USE moduleMesh
       USE moduleRefParam
       USE moduleSpecies
       USE moduleOutput
       USE moduleMeshInoutCommon
-      USE moduleCaseParam, ONLY: timeStep
       IMPLICIT NONE
 
       CLASS(meshParticles), INTENT(in):: self
+      INTEGER, INTENT(in):: t
       INTEGER:: n, i
       TYPE(outputFormat):: output(1:self%numNodes)
       REAL(8):: time
       CHARACTER(:), ALLOCATABLE:: fileName
 
-      time = DBLE(timeStep)*tauMin*ti_ref
+      time = DBLE(t)*tauMin*ti_ref
 
       DO i = 1, nSpecies
-        fileName = formatFileName(prefix, species(i)%obj%name, 'msh', timeStep)
+        fileName = formatFileName(prefix, species(i)%obj%name, 'msh', t)
         WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
         OPEN (60, file = path // folder // '/' //  fileName)
 
         CALL writeGmsh2HeaderMesh(60)
 
-        CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' density (m^-3)', timeStep, time, 1, self%numNodes)
+        CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' density (m^-3)', t, time, 1, self%numNodes)
         DO n=1, self%numNodes
           CALL calculateOutput(self%nodes(n)%obj%output(i), output(n), self%nodes(n)%obj%v, species(i)%obj)
           WRITE(60, "(I6,ES20.6E3)") n, output(n)%density
         END DO
         CALL writeGmsh2FooterNodeData(60)
 
-        CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' velocity (m s^-1)', timeStep, time, 3, self%numNodes)
+        CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' velocity (m s^-1)', t, time, 3, self%numNodes)
         DO n=1, self%numNodes
           WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%velocity
         END DO
         CALL writeGmsh2FooterNodeData(60)
 
-        CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' Pressure (Pa)', timeStep, time, 1, self%numNodes)
+        CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' Pressure (Pa)', t, time, 1, self%numNodes)
         DO n=1, self%numNodes
           WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%pressure
         END DO
         CALL writeGmsh2FooterNodeData(60)
 
-        CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' Temperature (K)', timeStep, time, 1, self%numNodes)
+        CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' Temperature (K)', t, time, 1, self%numNodes)
         DO n=1, self%numNodes
           WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%temperature
         END DO
@@ -135,7 +135,7 @@ MODULE moduleMeshOutputGmsh2
     END SUBROUTINE printOutputGmsh2
 
     !Prints the number of collisions into the volumes
-    SUBROUTINE printCollGmsh2(self)
+    SUBROUTINE printCollGmsh2(self, t)
       USE moduleMesh
       USE moduleRefParam
       USE moduleCaseParam
@@ -145,6 +145,7 @@ MODULE moduleMeshOutputGmsh2
       IMPLICIT NONE
 
       CLASS(meshGeneric), INTENT(in):: self
+      INTEGER, INTENT(in):: t
       INTEGER:: numEdges
       INTEGER:: k, c
       INTEGER:: n
@@ -166,9 +167,9 @@ MODULE moduleMeshOutputGmsh2
       END SELECT
 
       IF (collOutput) THEN
-        time = DBLE(timeStep)*tauMin*ti_ref
+        time = DBLE(t)*tauMin*ti_ref
 
-        fileName = formatFileName(prefix, 'Collisions', 'msh', timeStep)
+        fileName = formatFileName(prefix, 'Collisions', 'msh', t)
         WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
         OPEN (60, file = path // folder // '/' //  fileName)
 
@@ -178,7 +179,7 @@ MODULE moduleMeshOutputGmsh2
           DO c = 1, interactionMatrix(k)%amount
             WRITE(cString, "(I2)") c
             title = '"Pair ' // interactionMatrix(k)%sp_i%name // '-' // interactionMatrix(k)%sp_j%name // ' collision ' // cString
-            CALL writeGmsh2HeaderElementData(60, title, timeStep, time, 1, self%numCells)
+            CALL writeGmsh2HeaderElementData(60, title, t, time, 1, self%numCells)
             DO n=1, self%numCells
               WRITE(60, "(I6,I10)") n + numEdges, self%cells(n)%obj%tallyColl(k)%tally(c)
             END DO
@@ -195,7 +196,7 @@ MODULE moduleMeshOutputGmsh2
     END SUBROUTINE printCollGmsh2
 
     !Prints the electrostatic EM properties into the nodes and volumes
-    SUBROUTINE printEMGmsh2(self)
+    SUBROUTINE printEMGmsh2(self, t)
       USE moduleMesh
       USE moduleRefParam
       USE moduleCaseParam
@@ -204,6 +205,7 @@ MODULE moduleMeshOutputGmsh2
       IMPLICIT NONE
 
       CLASS(meshParticles), INTENT(in):: self
+      INTEGER, INTENT(in):: t
       INTEGER:: n, e
       REAL(8):: time
       CHARACTER(:), ALLOCATABLE:: fileName
@@ -212,27 +214,27 @@ MODULE moduleMeshOutputGmsh2
       Xi = (/ 0.D0, 0.D0, 0.D0 /)
 
       IF (emOutput) THEN
-        time = DBLE(timeStep)*tauMin*ti_ref
+        time = DBLE(t)*tauMin*ti_ref
 
-        fileName = formatFileName(prefix, 'EMField', 'msh', timeStep)
+        fileName = formatFileName(prefix, 'EMField', 'msh', t)
         WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
         OPEN (20, file = path // folder // '/' //  fileName)
 
         CALL writeGmsh2HeaderMesh(20)
 
-        CALL writeGmsh2HeaderNodeData(20, 'Potential (V)', timeStep, time, 1, self%numNodes)
+        CALL writeGmsh2HeaderNodeData(20, 'Potential (V)', t, time, 1, self%numNodes)
         DO n=1, self%numNodes
           WRITE(20, *) n, self%nodes(n)%obj%emData%phi*Volt_ref
         END DO
         CALL writeGmsh2FooterNodeData(20)
         
-        CALL writeGmsh2HeaderElementData(20, 'Electric Field (V m^-1)', timeStep, time, 3, self%numCells)
+        CALL writeGmsh2HeaderElementData(20, 'Electric Field (V m^-1)', t, time, 3, self%numCells)
         DO e=1, self%numCells
           WRITE(20, *) e+self%numEdges, self%cells(e)%obj%gatherElectricField(Xi)*EF_ref
         END DO
         CALL writeGmsh2FooterElementData(20)
 
-        CALL writeGmsh2HeaderNodeData(20, 'Magnetic Field (T)', timeStep, time, 3, self%numNodes)
+        CALL writeGmsh2HeaderNodeData(20, 'Magnetic Field (T)', t, time, 3, self%numNodes)
         DO n=1, self%numNodes
           WRITE(20, *) n, self%nodes(n)%obj%emData%B * B_ref
         END DO

src/modules/mesh/inout/makefile

@@ -1,4 +1,4 @@
-all: vtu.o gmsh2.o 0D.o text.o
+all: vtu.o gmsh2.o 0D.o
 
 vtu.o: moduleMeshInoutCommon.o
 	$(MAKE) -C vtu all
@@ -9,8 +9,5 @@ gmsh2.o:
 0D.o:
 	$(MAKE) -C 0D all
 
-text.o:
-	$(MAKE) -C text all
-
 %.o: %.f90
 	$(FC) $(FCFLAGS) -c $< -o $(OBJDIR)/$@

src/modules/mesh/inout/moduleMeshInoutCommon.f90

@@ -3,17 +3,17 @@ MODULE moduleMeshInoutCommon
   CHARACTER(LEN=4):: prefix = 'Step'
 
   CONTAINS
-    PURE FUNCTION formatFileName(prefix, suffix, extension, timeStep) RESULT(fileName)
+    PURE FUNCTION formatFileName(prefix, suffix, extension, t) RESULT(fileName)
       USE moduleOutput
       IMPLICIT NONE
 
       CHARACTER(*), INTENT(in):: prefix, suffix, extension
-      INTEGER, INTENT(in), OPTIONAL:: timeStep
+      INTEGER, INTENT(in), OPTIONAL:: t
       CHARACTER (LEN=iterationDigits):: tString
       CHARACTER(:), ALLOCATABLE:: fileName
 
-      IF (PRESENT(timeStep)) THEN
+      IF (PRESENT(t)) THEN
-        WRITE(tString, iterationFormat) timeStep
+        WRITE(tString, iterationFormat) t
         fileName = prefix //  '_' // tString // '_' // suffix // '.' // extension
 
       ELSE

src/modules/mesh/inout/text/makefile

@@ -1,7 +0,0 @@
-all: moduleMeshInputText.o moduleMeshOutputText.o
-
-moduleMeshInputText.o: moduleMeshOutputText.o moduleMeshInputText.f90
-	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
-
-%.o: %.f90
-	$(FC) $(FCFLAGS) -c $< -o $(OBJDIR)/$@

src/modules/mesh/inout/text/moduleMeshInputText.f90

@@ -1,232 +0,0 @@
-module moduleMeshInputText
-  !The mesh is stored as a column-wise text file.
-  !Aimed for simple geometries in 1D
-
-  contains
-    !Inits the text mesh
-    subroutine initText(self)
-      use moduleMesh
-      use moduleMeshOutputText
-      implicit none
-
-      class(meshGeneric), intent(inout), target:: self 
-
-      if (associated(meshForMCC,self)) then
-        self%printColl => printCollText
-
-      end if
-
-      select type(self)
-      type is (meshParticles)
-        self%printOutput  => printOutputText
-        self%printEM      => printEMText
-        self%printAverage => printAverageText
-
-        self%readInitial  => readInitialText
-
-      end select
-
-      self%readMesh => readText
-
-    end subroutine initText
-
-    !Reads the text mesh
-    subroutine readText(self, filename)
-      use moduleMesh
-      use moduleMesh1DCart
-      use moduleMesh1DRad
-      use moduleErrors
-      implicit none
-
-      class(meshGeneric), intent(inout):: self
-      character(:), allocatable, intent(in):: filename !Dummy file, not used
-      integer:: fileID, reason
-      character(len=256):: line
-      integer:: nNodes
-      real(8):: r(1:3) !dummy 3D coordinate
-      integer:: physicalID
-      integer:: n, c
-      integer, allocatable:: p(:)
-      integer:: bt
-
-      fileID = 10
-
-      open(fileID, file=trim(filename))
-
-      !Skip header
-      read(fileID, *)
-
-      !Get number of nodes
-      nNodes = 0
-      do
-        read(fileID, *, iostat=reason) line
-
-        if (reason > 0) then
-          call criticalError('Error reading mesh file', 'readText')
-
-        else if (reason < 0) then
-          exit
-
-        else if (len(line) > 0) then
-          nNodes = nNodes + 1
-
-        end if
-
-      end do
-
-      if (nNodes == 0) then
-        call criticalError('No nodes read in mesh file', 'readText')
-
-      end if
-
-      self%numNodes = nNodes
-      allocate(self%nodes(1:self%numNodes))
-
-      SELECT TYPE(self)
-      TYPE IS(meshParticles)
-        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
-
-      END SELECT
-
-      self%numCells = nNodes - 1
-      allocate(self%cells(1:self%numCells))
-
-      select type(self)
-      type is (meshParticles)
-        self%numEdges = 2
-
-        allocate(self%edges(1:self%numEdges))
-
-      end select
-
-      !Read the mesh now
-      rewind(fileID)
-      
-      !Skip header
-      read(fileID, *)
-
-      !Allocate nodes and edges
-      do n = 1, self%numNodes
-        r = 0.D0
-
-        read(fileID, *) r(1), physicalID
-
-        select case(self%geometry)
-        case("Cart")
-          allocate(meshNode1DCart:: self%nodes(n)%obj)
-
-        case("Rad")
-          allocate(meshNode1DRad:: self%nodes(n)%obj)
-
-        end select
-
-        !Init nodes
-        call self%nodes(n)%obj%init(n, r)
-
-        !Allocate edges if required)
-        select type(self)
-        type is (meshParticles)
-          if ((physicalID == 1) .or. (physicalID == 2)) then
-            select case(self%geometry)
-            case("Cart")
-              allocate(meshEdge1DCart:: self%edges(physicalID)%obj)
-
-            case("Rad")
-              allocate(meshEdge1DRad:: self%edges(physicalID)%obj)
-            
-            end select
-
-            allocate(p(1))
-            p(1) = n
-            bt = getBoundaryId(physicalID)
-            call self%edges(physicalID)%obj%init(physicalID, p, physicalID, physicalID)
-            deallocate(p)
-
-          end if
-
-        end select
-
-      end do
-
-      !Allocate cells
-      n = 1
-      allocate(p(1:2))
-      do c = 1, self%numCells
-        p(1) = n
-        n = n + 1
-        p(2) = n
-        
-        select case(self%geometry)
-        case("Cart")
-          allocate(meshCell1DCartSegm:: self%cells(c)%obj)
-
-        case("Rad")
-          allocate(meshCell1DRadSegm:: self%cells(c)%obj)
-
-        end select
-
-        call self%cells(c)%obj%init(c, p, self%nodes)
-
-
-      end do
-      deallocate(p)
-
-      close(fileID)
-
-      !Call mesh connectivity
-      CALL self%connectMesh
-
-    end subroutine readText
-
-    subroutine readInitialText(filename, density, velocity, temperature)
-      use moduleErrors
-      implicit none
-
-      character(:), allocatable, intent(in):: filename
-      real(8), allocatable, intent(out), dimension(:):: density
-      real(8), allocatable, intent(out), dimension(:,:):: velocity
-      real(8), allocatable, intent(out), dimension(:):: temperature
-      integer:: fileID, reason
-      character(len=256):: line
-      integer:: nNodes
-      integer:: n
-
-      fileID = 10
-
-      open(fileID, file=trim(filename))
-
-      do 
-        read(fileID, *, iostat=reason) line
-      
-        if (reason > 0) then
-          call criticalError('Error reading mesh file', 'readText')
-
-        else if (reason < 0) then
-          exit
-
-        else if (len(line) > 0) then
-          nNodes = nNodes + 1
-
-        end if
-
-      end do
-
-      allocate(density(1:nNodes))
-      allocate(velocity(1:nNodes, 1:3))
-      allocate(temperature(1:nNodes))
-
-      rewind(fileID)
-
-      do n = 1, nNodes
-        read(fileID, *) density(n), velocity(n, 1:3), temperature(n)
-
-      end do
-
-      close(fileID)
-
-    end subroutine readInitialText
-
-end module moduleMeshInputText

src/modules/mesh/inout/text/moduleMeshOutputText.f90

@@ -1,265 +0,0 @@
-module moduleMeshOutputText
-  contains
-
-    subroutine writeSpeciesOutput(self, fileID, speciesIndex)
-      use moduleMesh
-      use moduleOutput
-      use moduleRefParam, only: L_ref
-      implicit none
-
-      class(meshParticles), INTENT(in):: self
-      integer, intent(in):: fileID
-      integer, intent(in):: speciesIndex
-      real(8):: r(1:3)
-      type(outputFormat):: output
-      integer:: n
-
-      do n = 1, self%numNodes
-        r = self%nodes(n)%obj%getCoordinates()
-        call calculateOutput(self%nodes(n)%obj%output(speciesIndex), output, self%nodes(n)%obj%v, species(speciesIndex)%obj)
-
-        write(fileID, '(5(ES0.6E3,","),ES0.6E3)') r(1)*L_ref, output%density, output%velocity, output%temperature
-
-      end do
-
-    end subroutine writeSpeciesOutput
-
-    subroutine writeCollOutput(self, fileID)
-      use moduleMesh
-      use moduleCollisions
-      use moduleRefParam, only: L_ref
-      implicit none
-
-      class(meshGeneric), intent(in):: self
-      integer, intent(in):: fileID
-      integer:: n, k, c
-      
-      do n = 1, self%numCells
-        write(fileID, '(I0)', advance='no') n
-
-        do k = 1, nCollPairs
-          do c = 1, interactionMatrix(k)%amount
-            write(fileID, '(",",I0)', advance='no') self%cells(n)%obj%tallyColl(k)%tally(c)
-
-          end do
-
-        end do
-        write(fileID, *)
-
-      end do
-
-    end subroutine writeCollOutput
-
-    subroutine writeEMOutput(self, fileID)
-      use moduleMesh
-      use moduleRefParam, only: L_ref, Volt_ref, B_ref, EF_ref
-      implicit none
-
-      class(meshParticles), intent(in):: self
-      integer, intent(in):: fileID
-      integer:: n, c
-      real(8):: r(1:3), Xi(1:3)
-
-      do n = 1, self%numNodes
-        r = self%nodes(n)%obj%getCoordinates()
-        if (n == self%numNodes) then
-          Xi = (/ 1.D0, 0.D0, 0.D0 /)
-          c = self%numNodes - 1
-
-        else
-          Xi = (/ 0.D0, 0.D0, 0.D0 /)
-          c = n
-
-        end if
-
-        associate(output => self%nodes(n)%obj%emData)
-          write(fileID, '(7(ES0.6E3,","),ES0.6E3)') r(1)*L_ref, &
-                                                    output%phi*Volt_ref, &
-                                                    self%cells(c)%obj%gatherElectricField(Xi)*EF_ref, &
-                                                    output%B*B_ref
-
-        end associate
-
-      end do
-
-    end subroutine writeEMOutput
-
-    subroutine writeAverage(self, fileIDMean, &
-                                  fileIDDeviation, &
-                                  speciesIndex)
-
-      use moduleMesh
-      use moduleOutput
-      use moduleAverage
-      use moduleRefParam, only: L_ref
-      implicit none
-
-      class(meshParticles), intent(in):: self
-      integer, intent(in):: fileIDMean, fileIDDeviation
-      INTEGER, intent(in):: speciesIndex
-      real(8):: r(1:3)
-      type(outputFormat):: outputMean
-      type(outputFormat):: outputDeviation
-      integer:: n
-
-      do n = 1, self%numNodes
-        r = self%nodes(n)%obj%getCoordinates()
-
-        call calculateOutput(averageScheme(n)%mean%output(speciesIndex),      outputMean,      &
-                             self%nodes(n)%obj%v, species(speciesIndex)%obj)
-
-        write(fileIDMean, '(5(ES0.6E3,","),ES0.6E3)') r(1)*L_ref, outputMean%density, outputMean%velocity, outputMean%temperature
-
-        call calculateOutput(averageScheme(n)%deviation%output(speciesIndex), outputDeviation, &
-                             self%nodes(n)%obj%v, species(speciesIndex)%obj)
-
-        write(fileIDDeviation, '(5(ES0.6E3,","),ES0.6E3)') r(1)*L_ref, outputDeviation%density, outputDeviation%velocity, outputDeviation%temperature
-
-      end do
-
-    end subroutine writeAverage
-
-    subroutine printOutputText(self)
-      use moduleMesh
-      use moduleSpecies
-      use moduleMeshInoutCommon
-      use moduleCaseParam, ONLY: timeStep
-      implicit none
-      
-      class(meshParticles), intent(in):: self
-
-      INTEGER:: s, fileID
-      character(:), allocatable:: fileName
-
-      fileID = 60
-
-      do s = 1, nSpecies
-        fileName = formatFileName(prefix, species(s)%obj%name, 'csv', timeStep)
-        write(*, "(6X,A15,A)") "Creating file: ", fileName
-        open (fileID, file = path // folder // '/' //  fileName)
-
-        write(fileID, '(5(A,","),A)') 'Position (m)', &
-                                      'Density (m^-3)', &
-                                      'Velocity (m s^-1):0', 'Velocity (m s^-1):1', 'Velocity (m s^-1):2', &
-                                      'Temperature (K)'
-
-        call writeSpeciesOutput(self, fileID, s)
-
-        close(fileID)
-
-      end do
-
-    end subroutine printOutputText
-
-    subroutine printCollText(self)
-      use moduleMesh
-      use moduleOutput
-      use moduleMeshInoutCommon
-      use moduleCaseParam, only: timeStep
-      implicit none
-      
-      class(meshGeneric), intent(in):: self
-      integer:: fileID
-      character(:), allocatable:: fileName
-      integer:: k, c
-      character (len=2):: cString
-
-      fileID = 62
-
-      if (collOutput) then
-        fileName = formatFileName(prefix, 'Collisions', 'csv', timeStep)
-        write(*, "(6X,A15,A)") "Creating file: ", fileName
-        open (fileID, file = path // folder // '/' //  fileName)
-
-        write(fileID, '(A)', advance='no') "Cell"
-        do k = 1, nCollPairs
-          do c = 1, interactionMatrix(k)%amount
-            write(cString, "(I2)") c
-            write(fileID, '(",",A)', advance='no') 'Pair ' // interactionMatrix(k)%sp_i%name // '-' // interactionMatrix(k)%sp_j%name // ' collision ' // cString
-
-          end do
-        end do
-
-        write(fileID, *)
-
-        call writeCollOutput(self, fileID)
-
-        close(fileID)
-
-      end if
-
-    end subroutine printCollText
-
-    subroutine printEMText(self)
-      use moduleMesh
-      use moduleMeshInoutCommon
-      use moduleCaseParam, only: timeStep
-      implicit none
-      
-      class(meshParticles), intent(in):: self
-      integer:: fileID
-      character(:), allocatable:: fileName
-
-      fileID = 64
-
-      if (emOutput) then
-        fileName = formatFileName(prefix, 'EMField', 'csv', timeStep)
-        write(*, "(6X,A15,A)") "Creating file: ", fileName
-        open (fileID, file = path // folder // '/' //  fileName)
-
-        write(fileID, '(8(A,","),A)') 'Position (m)', &
-                                      'Potential (V)', &
-                                      'Electric Field (V m^-1):0', 'Electric Field (V m^-1):1', 'Electric Field (V m^-1):2', &
-                                      'Magnetic Field (T):0',      'Magnetic Field (T):1',      'Magnetic Field (T):2'
-
-        call writeEMOutput(self, fileID)
-
-        close(fileID)
-
-      end if
-
-    end subroutine printEMText
-
-    subroutine printAverageText(self)
-      use moduleMesh
-      use moduleSpecies
-      use moduleMeshInoutCommon
-      implicit none
-      
-      class(meshParticles), intent(in):: self
-      integer:: s
-      integer:: fileIDMean, fileIDDeviation
-      character(:), allocatable:: fileNameMean, fileNameDeviation
-
-      fileIDMean = 66
-      fileIDDeviation = 67
-
-      do s = 1, nSpecies
-        fileNameMean = formatFileName('Average_mean', species(s)%obj%name, 'csv', timeStep)
-        write(*, "(6X,A15,A)") "Creating file: ", fileNameMean
-        open (fileIDMean, file = path // folder // '/' //  fileNameMean)
-
-        write(fileIDMean, '(5(A,","),A)') 'Position (m)', &
-                                          'Density, mean (m^-3)', &
-                                          'Velocity, mean (m s^-1):0', 'Velocity (m s^-1):1', 'Velocity (m s^-1):2', &
-                                          'Temperature, mean (K)'
-
-        fileNameDeviation = formatFileName('Average_deviation', species(s)%obj%name, 'csv', timeStep)
-        write(*, "(6X,A15,A)") "Creating file: ", fileNameDeviation
-        open (fileIDDeviation, file = path // folder // '/' //  fileNameDeviation)
-
-        write(fileIDDeviation, '(5(A,","),A)') 'Position (m)', &
-                                               'Density, deviation (m^-3)', &
-                                               'Velocity, deviation (m s^-1):0', 'Velocity (m s^-1):1', 'Velocity (m s^-1):2', &
-                                               'Temperature, deviation (K)'
-
-        call writeAverage(self, fileIDMean, fileIDDeviation, s)
-
-        close(fileIDMean)
-        close(fileIDDeviation)
-
-      end do
-
-    end subroutine printAverageText
-
-end module moduleMeshOutputText

src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90

@@ -167,7 +167,7 @@ MODULE moduleMeshInputVTU
       CLASS(meshGeneric), INTENT(inout):: self
       CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
       REAL(8):: r(1:3) !3 generic coordinates
-      INTEGER:: fileID
+      INTEGER:: fileID, error, found
       CHARACTER(LEN=256):: line
       INTEGER:: numNodes, numElements, numEdges
       INTEGER, ALLOCATABLE, DIMENSION(:):: entitiesID, offsets, connectivity, types
@@ -275,7 +275,6 @@ MODULE moduleMeshInputVTU
       END DO
 
       !Count the number of edges
-      numEdges = 0
       SELECT CASE(self%dimen)
       CASE(3)
         !Edges are triangles, type 5 in VTK
@@ -496,7 +495,7 @@ MODULE moduleMeshInputVTU
         END SELECT
 
       END DO
-
+        
       !Call mesh connectivity
       CALL self%connectMesh
 
@@ -548,8 +547,6 @@ MODULE moduleMeshInputVTU
       CALL readDataBlock(fileID, numNodes, temperature)
       REWIND(fileID)
 
-      close(fileID)
-
     END SUBROUTINE readInitialVTU
 
 END MODULE moduleMeshInputVTU

src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90

@@ -11,7 +11,7 @@ MODULE moduleMeshOutputVTU
 
       WRITE(fileID,"(A)") '<?xml version="1.0"?>'
       WRITE(fileID,"(2X, A)") '<VTKFile type="UnstructuredGrid">'
-      WRITE(fileID,"(4X, A)") '<UnstructuredGrid>'
+      WRITE(fileID,"(4X, A,ES20.6E3,A)") '<UnstructuredGrid>'
       WRITE(fileID,"(6X, A, I10, A, I10, A)") '<Piece NumberOfPoints="', nNodes, '" NumberOfCells="', nCells, '">'
 
     END SUBROUTINE writeHeader
@@ -209,22 +209,23 @@ MODULE moduleMeshOutputVTU
       WRITE(fileID,"(8X,A)") '<CellData>'
       !Electric field
       WRITE(fileID,"(10X,A, A, A)") '<DataArray type="Float64" Name="Electric Field (V m^-1)" NumberOfComponents="3">'
-      WRITE(fileID,"(6(ES20.6E3))") (self%cells(n)%obj%gatherElectricField(Xi)*EF_ref, n = 1, self%numCells)
+      WRITE(fileID, "(6(ES20.6E3))") (self%cells(n)%obj%gatherElectricField(Xi)*EF_ref, n = 1, self%numCells)
       WRITE(fileID,"(10X,A)") '</DataArray>'
       WRITE(fileID,"(8X,A)") '</CellData>'
 
     END SUBROUTINE writeEM
 
-    SUBROUTINE writeCollection(fileID, fileNameStep, fileNameCollection)
+    SUBROUTINE writeCollection(fileID, t, fileNameStep, fileNameCollection)
       USE moduleCaseParam
       USE moduleOutput
       USE moduleRefParam
       IMPLICIT NONE
 
       INTEGER:: fileID
+      INTEGER, INTENT(in):: t
       CHARACTER(*):: fileNameStep, fileNameCollection
 
-      IF (timeStep == tInitial) THEN
+      IF (t == tInitial) THEN
         !Create collection file
         WRITE(*, "(6X,A15,A)") "Creating file: ", fileNameCollection
         OPEN (fileID + 1, file = path // folder  // '/' // fileNameCollection)
@@ -236,11 +237,10 @@ MODULE moduleMeshOutputVTU
 
       !Write iteration file in collection
       OPEN (fileID + 1, file = path // folder  // '/' // fileNameCollection, ACCESS='APPEND')
-      WRITE(fileID + 1, "(4X, A, ES20.6E3, A, A, A)") &
+      WRITE(fileID + 1, "(4X, A, ES20.6E3, A, A, A)") '<DataSet timestep="', DBLE(t)*tauMin*ti_ref,'" file="', fileNameStep,'"/>'
-            '<DataSet timestep="', DBLE(timeStep)*tauMin*ti_ref,'" file="', fileNameStep,'"/>'
 
       !Close collection file
-      IF (timeStep == tFinal) THEN
+      IF (t == tFinal) THEN
         WRITE (fileID + 1, "(2X, A)") '</Collection>'
         WRITE (fileID + 1, "(A)") '</VTKFile>'
 
@@ -307,21 +307,22 @@ MODULE moduleMeshOutputVTU
 
     END SUBROUTINE writeAverage
 
-    SUBROUTINE printOutputVTU(self)
+    SUBROUTINE printOutputVTU(self,t)
       USE moduleMesh
       USE moduleSpecies
       USE moduleMeshInoutCommon
-      USE moduleCaseParam, ONLY: timeStep
       IMPLICIT NONE
 
       CLASS(meshParticles), INTENT(in):: self
-      INTEGER:: i, fileID
+      INTEGER, INTENT(in):: t
+      INTEGER:: n, i, fileID
       CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection
+      TYPE(outputFormat):: output(1:self%numNodes)
 
       fileID = 60
 
       DO i = 1, nSpecies
-        fileName = formatFileName(prefix, species(i)%obj%name, 'vtu', timeStep)
+        fileName = formatFileName(prefix, species(i)%obj%name, 'vtu', t)
         WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
         OPEN (fileID, file = path // folder // '/' //  fileName)
 
@@ -337,27 +338,29 @@ MODULE moduleMeshOutputVTU
 
         !Write collection file for time plotting
         fileNameCollection = formatFileName('Collection', species(i)%obj%name, 'pvd')
-        CALL writeCollection(fileID, fileName, filenameCollection)
+        CALL writeCollection(fileID, t, fileName, filenameCollection)
 
       END DO
 
     END SUBROUTINE printOutputVTU
 
-    SUBROUTINE printCollVTU(self)
+    SUBROUTINE printCollVTU(self,t)
       USE moduleMesh
       USE moduleOutput
       USE moduleMeshInoutCommon
-      USE moduleCaseParam, ONLY: timeStep
       IMPLICIT NONE
 
       CLASS(meshGeneric), INTENT(in):: self
-      INTEGER:: fileID
+      INTEGER, INTENT(in):: t
+      INTEGER:: n, i, fileID
       CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection
+      CHARACTER (LEN=iterationDigits):: tstring
+      TYPE(outputFormat):: output(1:self%numNodes)
 
       fileID = 62
 
       IF (collOutput) THEN
-        fileName = formatFileName(prefix, 'Collisions', 'vtu', timeStep)
+        fileName = formatFileName(prefix, 'Collisions', 'vtu', t)
         WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
         OPEN (fileID, file = path // folder // '/' //  fileName)
 
@@ -373,26 +376,26 @@ MODULE moduleMeshOutputVTU
 
         !Write collection file for time plotting
         fileNameCollection = formatFileName('Collection', 'Collisions', 'pvd')
-        CALL writeCollection(fileID, fileName, filenameCollection)
+        CALL writeCollection(fileID, t, fileName, filenameCollection)
 
       END IF
 
     END SUBROUTINE printCollVTU
 
-    SUBROUTINE printEMVTU(self)
+    SUBROUTINE printEMVTU(self, t)
       USE moduleMesh
       USE moduleMeshInoutCommon
-      USE moduleCaseParam, ONLY: timeStep
       IMPLICIT NONE
 
       CLASS(meshParticles), INTENT(in):: self
+      INTEGER, INTENT(in):: t
       INTEGER:: fileID
       CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection
 
       fileID = 64
 
       IF (emOutput) THEN
-        fileName = formatFileName(prefix, 'EMField', 'vtu', timeStep)
+        fileName = formatFileName(prefix, 'EMField', 'vtu', t)
         WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
         OPEN (fileID, file = path // folder // '/' //  fileName)
 
@@ -408,7 +411,7 @@ MODULE moduleMeshOutputVTU
 
         !Write collection file for time plotting
         fileNameCollection = formatFileName('Collection', 'EMField', 'pvd')
-        CALL writeCollection(fileID, fileName, filenameCollection)
+        CALL writeCollection(fileID, t, fileName, filenameCollection)
 
       END IF
 
@@ -421,8 +424,9 @@ MODULE moduleMeshOutputVTU
       IMPLICIT NONE
 
       CLASS(meshParticles), INTENT(in):: self
-      INTEGER:: i, fileIDMean, fileIDDeviation
+      INTEGER:: n, i, fileIDMean, fileIDDeviation
       CHARACTER(:), ALLOCATABLE:: fileNameMean, fileNameDeviation
+      TYPE(outputFormat):: output(1:self%numNodes)
 
       fileIDMean = 66
       fileIDDeviation = 67

src/modules/mesh/moduleMesh.f90

@@ -59,13 +59,6 @@ MODULE moduleMesh
 
   END TYPE meshNodeCont
 
-  ! Array of pointers to nodes.
-  TYPE:: meshNodePointer
-    CLASS(meshNode), POINTER:: obj
-    CONTAINS
-
-  END TYPE meshNodePointer
-
   !Type for array of boundary functions (one per species)
   TYPE, PUBLIC:: fBoundaryGeneric
     PROCEDURE(boundary_interface), POINTER, NOPASS:: apply => NULL()
@@ -83,8 +76,8 @@ MODULE moduleMesh
     CLASS(meshCell), POINTER:: eColl => NULL()
     !Normal vector
     REAL(8):: normal(1:3)
-    ! Surface of edge
+    !Weight for random injection of particles
-    REAL(8):: surface = 0.D0
+    REAL(8):: weight = 1.D0
     !Pointer to boundary type
     TYPE(boundaryCont), POINTER:: boundary
     !Array of functions for boundary conditions
@@ -344,10 +337,10 @@ MODULE moduleMesh
     !Array of cell elements
     TYPE(meshCellCont),  ALLOCATABLE:: cells(:)
     !PROCEDURES SPECIFIC OF FILE TYPE
-    PROCEDURE(readMesh_interface),    POINTER, PASS::   readMesh     => NULL()
+    PROCEDURE(readMesh_interface),     POINTER, PASS::   readMesh     => NULL()
-    PROCEDURE(readInitial_interface), POINTER, NOPASS:: readInitial  => NULL()
+    PROCEDURE(readInitial_interface),  POINTER, NOPASS:: readInitial  => NULL()
-    PROCEDURE(connectMesh_interface), POINTER, PASS::   connectMesh  => NULL()
+    PROCEDURE(connectMesh_interface),  POINTER, PASS::   connectMesh  => NULL()
-    PROCEDURE(printColl_interface),   POINTER, PASS::   printColl    => NULL()
+    PROCEDURE(printColl_interface),    POINTER, PASS::   printColl    => NULL()
     CONTAINS
       !GENERIC PROCEDURES
       PROCEDURE, PASS:: doCollisions
@@ -379,9 +372,10 @@ MODULE moduleMesh
     END SUBROUTINE connectMesh_interface
 
     !Prints number of collisions in each cell
-    SUBROUTINE printColl_interface(self)
+    SUBROUTINE printColl_interface(self, t)
       IMPORT meshGeneric
       CLASS(meshGeneric), INTENT(in):: self
+      INTEGER, INTENT(in):: t
 
     END SUBROUTINE printColl_interface
 
@@ -409,16 +403,18 @@ MODULE moduleMesh
 
   ABSTRACT INTERFACE
     !Prints Species data
-    SUBROUTINE printOutput_interface(self)
+    SUBROUTINE printOutput_interface(self, t)
       IMPORT meshParticles
       CLASS(meshParticles), INTENT(in):: self
+      INTEGER, INTENT(in):: t
 
     END SUBROUTINE printOutput_interface
 
     !Prints EM info
-    SUBROUTINE printEM_interface(self)
+    SUBROUTINE printEM_interface(self, t)
       IMPORT meshParticles
       CLASS(meshParticles), INTENT(in):: self
+      INTEGER, INTENT(in):: t
 
     END SUBROUTINE printEM_interface
 
@@ -499,29 +495,28 @@ MODULE moduleMesh
       IMPLICIT NONE
 
       CLASS(meshParticles), INTENT(inout):: self
-      INTEGER:: c
+      INTEGER:: e
+      INTEGER:: nNodes
       INTEGER, ALLOCATABLE:: n(:)
       REAL(8), ALLOCATABLE:: localK(:,:)
       INTEGER:: i, j
 
-      DO c = 1, self%numCells
+      DO e = 1, self%numCells
-        associate(nNodes => self%cells(c)%obj%nNodes)
+        nNodes = self%cells(e)%obj%nNodes
-          ALLOCATE(n(1:nNodes))
+        ALLOCATE(n(1:nNodes))
-          ALLOCATE(localK(1:nNodes, 1:nNodes))
+        ALLOCATE(localK(1:nNodes, 1:nNodes))
-          n      = self%cells(c)%obj%getNodes(nNodes)
+        n      = self%cells(e)%obj%getNodes(nNodes)
-          localK = self%cells(c)%obj%elemK(nNodes)
+        localK = self%cells(e)%obj%elemK(nNodes)
 
-          DO i = 1, nNodes
+        DO i = 1, nNodes
-            DO j = 1, nNodes
+          DO j = 1, nNodes
-              self%K(n(i), n(j)) = self%K(n(i), n(j)) + localK(i, j)
+            self%K(n(i), n(j)) = self%K(n(i), n(j)) + localK(i, j)
-
-            END DO
 
           END DO
 
-          DEALLOCATE(n, localK)
+        END DO
 
-        end associate
+        DEALLOCATE(n, localK)
 
       END DO
 
@@ -618,7 +613,6 @@ MODULE moduleMesh
       INTEGER:: sp
       INTEGER:: i
       CLASS(meshNode), POINTER:: node
-      REAL(8):: pFraction !Particle fraction
 
       cellNodes = self%getNodes(nNodes)
       fPsi = self%fPsi(part%Xi, nNodes)
@@ -629,11 +623,10 @@ MODULE moduleMesh
 
       DO i = 1, nNodes
         node => mesh%nodes(cellNodes(i))%obj
-        pFraction = fPsi(i)*part%weight
         CALL OMP_SET_LOCK(node%lock)
-        node%output(sp)%den          = node%output(sp)%den          + pFraction
+        node%output(sp)%den          = node%output(sp)%den          + part%weight*fPsi(i)
-        node%output(sp)%mom(:)       = node%output(sp)%mom(:)       + pFraction*part%v(:)
+        node%output(sp)%mom(:)       = node%output(sp)%mom(:)       + part%weight*fPsi(i)*part%v(:)
-        node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + pFraction*tensorS
+        node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*fPsi(i)*tensorS
         CALL OMP_UNSET_LOCK(node%lock)
 
       END DO
@@ -794,7 +787,7 @@ MODULE moduleMesh
     END FUNCTION findCellBrute
 
     !Computes collisions in element
-    SUBROUTINE doCollisions(self)
+    SUBROUTINE doCollisions(self, t)
       USE moduleCollisions
       USE moduleSpecies
       USE moduleList
@@ -802,10 +795,10 @@ MODULE moduleMesh
       USE moduleRandom
       USE moduleOutput
       USE moduleMath
-      USE moduleCaseParam, ONLY: timeStep
       IMPLICIT NONE
 
       CLASS(meshGeneric), INTENT(inout), TARGET:: self
+      INTEGER, INTENT(in):: t
       INTEGER:: e
       CLASS(meshCell), POINTER:: cell
       INTEGER:: k, i, j
@@ -821,11 +814,10 @@ MODULE moduleMesh
       REAL(8):: rnd_real !Random number for collision
       INTEGER:: rnd_int  !Random number for collision
 
-      IF (MOD(timeStep, everyColl) == 0) THEN
+      IF (MOD(t, everyColl) == 0) THEN
         !Collisions need to be performed in this iteration
-        !$OMP DO SCHEDULE(DYNAMIC) PRIVATE(part_i, part_j, partTemp_i, partTemp_j)
+        !$OMP DO SCHEDULE(DYNAMIC) PRIVATE(part_i, part_j, partTemp_i, partTemp_j, cell)
         DO e=1, self%numCells
-
           cell => self%cells(e)%obj
 
           !TODO: Simplify this, to many sublevels
@@ -892,7 +884,7 @@ MODULE moduleMesh
                   !Obtain the cross sections for the different processes
                   !TODO: From here it might be a procedure in interactionMatrix
                   vRel = NORM2(part_i%v-part_j%v)
-                  rMass = reducedMass(part_i%weight*part_i%species%m, part_j%weight*part_j%species%m)
+                  rMass = reducedMass(part_i%weight*part_i%species%mass, part_j%weight*part_j%species%mass)
                   eRel = rMass*vRel**2
                   CALL interactionMatrix(k)%getSigmaVrel(vRel, eRel, sigmaVrelTotal, sigmaVrel)
 
@@ -918,9 +910,7 @@ MODULE moduleMesh
                     !Loop over collisions
                     DO c = 1, interactionMatrix(k)%amount
                       IF (rnd_real <= probabilityColl(c)) THEN
-                        !$OMP CRITICAL
                         CALL interactionMatrix(k)%collisions(c)%obj%collide(part_i, part_j, vRel)
-                        !$OMP END CRITICAL
 
                         !If collisions are gonna be output, count the collision
                         IF (collOutput) THEN
@@ -1030,9 +1020,6 @@ MODULE moduleMesh
           ALLOCATE(deltaV_ij(1:cell%listPart_in(i)%amount, 1:3))
           ALLOCATE(p_ij(1:cell%listPart_in(i)%amount, 1:3))
           ALLOCATE(mass_ij(1:cell%listPart_in(i)%amount))
-          deltaV_ij = 0.D0
-          p_ij      = 0.D0
-          mass_ij   = 0.D0
           !Loop over particles of species_i
           partTemp => cell%listPart_in(i)%head
           p = 1
@@ -1092,7 +1079,7 @@ MODULE moduleMesh
 
             !Compute changes in velocity for each particle
             deltaV_ij(p,1:3) = MATMUL(rotation, W) + velocity - partTemp%part%v
-            mass_ij(p)       = pair%sp_i%m*partTemp%part%weight
+            mass_ij(p)       = pair%sp_i%mass*partTemp%part%weight
             p_ij(p,1:3)      = mass_ij(p)*partTemp%part%v
 
             !Move to the next particle in the list
@@ -1117,9 +1104,6 @@ MODULE moduleMesh
             ALLOCATE(deltaV_ji(1:cell%listPart_in(j)%amount, 1:3))
             ALLOCATE(p_ji(1:cell%listPart_in(j)%amount, 1:3))
             ALLOCATE(mass_ji(1:cell%listPart_in(j)%amount))
-            deltaV_ji = 0.D0
-            p_ji      = 0.D0
-            mass_ji   = 0.D0
             !Loop over particles of species_j
             partTemp => cell%listPart_in(j)%head
             p = 1
@@ -1178,7 +1162,7 @@ MODULE moduleMesh
 
               !Compute changes in velocity for each particle
               deltaV_ji(p,1:3) = MATMUL(rotation, W) + velocity - partTemp%part%v
-              mass_ji(p)       = pair%sp_j%m*partTemp%part%weight
+              mass_ji(p)       = pair%sp_j%mass*partTemp%part%weight
               p_ji(p,1:3)      = mass_ji(p)*partTemp%part%v
 
               !Move to the next particle in the list

src/modules/mesh/moduleMeshBoundary.f90

@@ -77,20 +77,6 @@ MODULE moduleMeshBoundary
 
     END SUBROUTINE transparent
 
-    !Symmetry axis. Reflects particles.
-    !Although this function should never be called, it is set as a reflective boundary
-    !to properly deal with possible particles reaching a corner and selecting this boundary.
-    SUBROUTINE symmetryAxis(edge, part)
-      USE moduleSpecies
-      IMPLICIT NONE
-
-      CLASS(meshEdge), INTENT(inout):: edge
-      CLASS(particle), INTENT(inout):: part
-
-      CALL reflection(edge, part)
-
-    END SUBROUTINE symmetryAxis
-
     !Wall with temperature
     SUBROUTINE wallTemperature(edge, part)
       USE moduleSpecies
@@ -139,7 +125,7 @@ MODULE moduleMeshBoundary
 
       SELECT TYPE(bound => edge%boundary%bTypes(part%species%n)%obj)
       TYPE IS(boundaryIonization)
-        mRel = reducedMass(bound%m0, part%species%m)
+        mRel = reducedMass(bound%m0, part%species%mass)
         vRel = SUM(DABS(part%v-bound%v0))
         eRel = mRel*vRel**2*5.D-1
 
@@ -161,13 +147,7 @@ MODULE moduleMeshBoundary
             ALLOCATE(newElectron)
             ALLOCATE(newIon)
 
-            IF (ASSOCIATED(bound%electronSecondary)) THEN
+            newElectron%species => part%species
-              newElectron%species => bound%electronSecondary
-
-            ELSE
-              newElectron%species => part%species
-
-            END IF
             newIon%species      => bound%species
 
             newElectron%v = v0 + (1.D0 + bound%deltaV*v0/NORM2(v0))
@@ -213,32 +193,134 @@ MODULE moduleMeshBoundary
 
       END SELECT
 
-      !Removes ionizing electron regardless the number of pair created
+      !Removes ionizing electron regardless the number of pairs created
       part%n_in = .FALSE.
 
     END SUBROUTINE ionization
 
-    subroutine outflowAdaptive(edge, part)
+    !Symmetry axis. Reflects particles.
-      use moduleRandom
+    !Although this function should never be called, it is set as a reflective boundary
-      implicit none
+    !to properly deal with possible particles reaching a corner and selecting this boundary.
+    SUBROUTINE symmetryAxis(edge, part)
+      USE moduleSpecies
+      IMPLICIT NONE
 
-      class(meshEdge), intent(inout):: edge
+      CLASS(meshEdge), INTENT(inout):: edge
-      class(particle), intent(inout):: part
+      CLASS(particle), INTENT(inout):: part
 
-      select type(bound => edge%boundary%bTypes(part%species%n)%obj)
+      CALL reflection(edge, part)
-      type is(boundaryOutflowAdaptive)
 
-        if (random() < 0.844d0) then
+    END SUBROUTINE symmetryAxis
-          call reflection(edge, part)
 
-        else
+    !Secondary emission of electrons by impacting particle.
-          call transparent(edge, part)
+    SUBROUTINE secondaryEmission(edge, part)
+      USE moduleSpecies
+      USE moduleRandom
+      USE moduleConstParam
+      IMPLICIT NONE
 
-        end if
+      CLASS(meshEdge), INTENT(inout):: edge
+      CLASS(particle), INTENT(inout):: part
+      REAL(8):: vRel, eRel
+      REAL(8), DIMENSION(1:3):: rElectron, XiElectron!Position of new electrons (impacting particle)
+      REAL(8), DIMENSION(1:3):: rIncident !Vector from imapcting particle position to particle position
+      REAL(8):: theta !incident angle
+      REAL(8):: yield
+      REAL(8):: energy !incident energy corrected by threshold and 
+      INTEGER:: nNewElectrons
+      REAL(8), ALLOCATABLE:: weight(:)
+      INTEGER:: p
+      INTEGER:: cell
+      TYPE(particle), POINTER:: secondaryElectron
 
-      end select
+      SELECT TYPE(bound => edge%boundary%bTypes(part%species%n)%obj)
+      TYPE IS(boundarySEE)
+        !Get relative velocity
+        vRel = NORM2(part%v)
+        !Convert to relative energy
+        eRel = part%species%mass*vRel**2*5.D-1
 
-    end subroutine outflowAdaptive
+        !If energy is abound threshold calculate secondary electrons
+        IF (eRel >= bound%energyThreshold) THEN
+
+          !position of impacting ion
+          rElectron  = edge%intersection(part%r)
+          XiElectron = mesh%cells(part%cell)%obj%phy2log(rElectron)
+
+          !Calculate incident angle
+          rIncident = part%r - rElectron
+          theta = ACOS(DOT_PRODUCT(rIncident, edge%normal) / (NORM2(rIncident) * NORM2(edge%normal)))
+
+          !Calculate incident energy
+          energy = (eRel - bound%energyThreshold) * PI2 / (PI2 + theta**2) + bound%energyThreshold
+
+          !Get number of secondary electrons particles
+          yield = part%weight*bound%yield%get(eRel) * (1.D0 * theta**2 / PI2) !Check equation!
+          
+          !Convert the number to macro-particles
+          nNewElectrons = FLOOR(yield / bound%electron%weight)
+
+          !If the weight of new macro-particles is below the yield, correct adding a particle
+          IF (REAL(nNewElectrons) * bound%electron%weight < yield) THEN
+            nNewElectrons = nNewElectrons + 1
+            ALLOCATE(weight(1:nNewElectrons))
+            weight(1:nNewElectrons-1) = bound%electron%weight
+            weight(nNewElectrons)     = yield - SUM(weight(1:nNewElectrons-1))
+
+          ELSE
+            ALLOCATE(weight(1:nNewElectrons))
+            weight(1:nNewElectrons) = bound%electron%weight
+          
+          END IF
+
+          !New cell of origin
+          IF (ASSOCIATED(edge%e1)) THEN
+            cell = edge%e1%n
+
+          ELSEIF (ASSOCIATED(edge%e2)) THEN
+            cell = edge%e2%n
+
+          END IF
+
+          !Create the new electron macro-particles
+          DO p = 1, nNewElectrons
+            !Create new macro-particle
+            ALLOCATE(secondaryElectron)
+
+            !Assign species to electron
+            secondaryElectron%species => bound%electron
+
+            !Assign position to particle
+            secondaryElectron%r  = rElectron
+            secondaryElectron%Xi = XiElectron
+
+            !Assign cell to electron
+            secondaryElectron%cell = cell
+
+            !Assign weight
+            secondaryElectron%weight = weight(p)
+
+            !Assume particle is inside the numerical domain
+            secondaryElectron%n_in = .TRUE.
+
+            !Assign velocity
+            secondaryElectron%v = 2.D0 * edge%normal + 1.D-1 * (/ randomMaxwellian(), randomMaxwellian(), randomMaxwellian() /)
+
+            !Add particle to list
+            CALL partSurfaces%setLock()
+            CALL partSurfaces%add(secondaryElectron)
+            CALL partSurfaces%unsetLock()
+
+          END DO
+
+        END IF
+
+        !Absorb impacting particle
+        CALL absorption(edge, part)
+
+      END SELECT
+
+    END SUBROUTINE secondaryEmission
 
     !Points the boundary function to specific type
     SUBROUTINE pointBoundaryFunction(edge, s)
@@ -258,17 +340,17 @@ MODULE moduleMeshBoundary
       TYPE IS(boundaryTransparent)
         edge%fBoundary(s)%apply => transparent
 
-      TYPE IS(boundaryAxis)
-        edge%fBoundary(s)%apply => symmetryAxis
-
       TYPE IS(boundaryWallTemperature)
         edge%fBoundary(s)%apply => wallTemperature
 
       TYPE IS(boundaryIonization)
         edge%fBoundary(s)%apply => ionization
 
-      type is(boundaryOutflowAdaptive)
+      TYPE IS(boundaryAxis)
-        edge%fBoundary(s)%apply => outflowAdaptive
+        edge%fBoundary(s)%apply => symmetryAxis
+
+      TYPE IS(boundarySEE)
+        edge%fBoundary(s)%apply => secondaryEmission
 
       CLASS DEFAULT
         CALL criticalError("Boundary type not defined in this geometry", 'pointBoundaryFunction')

src/modules/moduleBoundary.f90

@@ -26,12 +26,6 @@ MODULE moduleBoundary
 
   END TYPE boundaryTransparent
 
-  !Symmetry axis
-  TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryAxis
-    CONTAINS
-
-  END TYPE boundaryAxis
-
   !Wall Temperature boundary
   TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryWallTemperature
     !Thermal velocity of the wall: square root(Wall temperature X specific heat)
@@ -44,7 +38,6 @@ MODULE moduleBoundary
   TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryIonization
     REAL(8):: m0, n0, v0(1:3), vTh !Properties of background neutrals.
     CLASS(speciesGeneric), POINTER:: species !Ion species
-    CLASS(speciesCharged), POINTER:: electronSecondary !Pointer to species considerer as secondary electron
     TYPE(table1D):: crossSection
     REAL(8):: effectiveTime
     REAL(8):: eThreshold
@@ -53,13 +46,21 @@ MODULE moduleBoundary
 
   END TYPE boundaryIonization
 
-  !Boundary for quasi-neutral outflow adjusting reflection coefficient
+  !Secondary electron emission (by ion impact)
-  type, public, extends(boundaryGeneric):: boundaryOutflowAdaptive
+  TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundarySEE
-    real(8):: outflowCurrent
+    !Yield as a function of ion energy
-    real(8):: reflectionFraction
+    TYPE(table1D):: yield
-    contains
+    CLASS(speciesGeneric), POINTER:: electron !Electron species for secondary emission
+    REAL(8):: energyThreshold
+    CONTAINS
 
-  end type boundaryOutflowAdaptive
+  END TYPE boundarySEE
+
+  !Symmetry axis
+  TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryAxis
+    CONTAINS
+
+  END TYPE boundaryAxis
 
   !Wrapper for boundary types (one per species)
   TYPE:: bTypesCont
@@ -112,19 +113,17 @@ MODULE moduleBoundary
 
     END SUBROUTINE initWallTemperature
 
-    SUBROUTINE initIonization(boundary, me, m0, n0, v0, T0, ion, effTime, crossSection, eThreshold, electronSecondary)
+    SUBROUTINE initIonization(boundary, me, m0, n0, v0, T0, speciesID, effTime, crossSection, eThreshold)
       USE moduleRefParam
       USE moduleSpecies
       USE moduleCaseParam
       USE moduleConstParam
-      USE moduleErrors
       IMPLICIT NONE
 
       CLASS(boundaryGeneric), ALLOCATABLE, INTENT(out):: boundary
       REAL(8), INTENT(in):: me !Electron mass
       REAL(8), INTENT(in):: m0, n0, v0(1:3), T0 !Neutral properties
-      INTEGER, INTENT(in):: ion
+      INTEGER:: speciesID
-      INTEGER, OPTIONAL, INTENT(in):: electronSecondary
       REAL(8):: effTime
       CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSection
       REAL(8), INTENT(in):: eThreshold
@@ -137,22 +136,7 @@ MODULE moduleBoundary
         boundary%n0 = n0 * Vol_ref
         boundary%v0 = v0 / v_ref
         boundary%vTh = DSQRT(kb*T0/m0)/v_ref
-        boundary%species => species(ion)%obj
+        boundary%species => species(speciesID)%obj
-        IF (PRESENT(electronSecondary)) THEN
-          SELECT TYPE(sp => species(electronSecondary)%obj)
-          TYPE IS(speciesCharged)
-            boundary%electronSecondary => sp
-
-          CLASS DEFAULT
-            CALL criticalError("Species " // sp%name // " chosen for " // &
-                               "secondary electron is not a charged species", 'initIonization')
-
-          END SELECT
-
-        ELSE
-          boundary%electronSecondary => NULL()
-
-        END IF
         boundary%effectiveTime = effTime / ti_ref
         CALL boundary%crossSection%init(crossSection)
         CALL boundary%crossSection%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2)
@@ -163,4 +147,36 @@ MODULE moduleBoundary
 
     END SUBROUTINE initIonization
 
+    SUBROUTINE initSEE(boundary, tableFile, speciesID)
+      USE moduleRefParam
+      USE moduleConstParam
+      USE moduleSpecies
+      IMPLICIT NONE
+
+      CLASS(boundaryGeneric), ALLOCATABLE, INTENT(out):: boundary
+      CHARACTER(:), ALLOCATABLE, INTENT(in):: tableFile
+      INTEGER:: speciesID
+      INTEGER:: i
+
+      ALLOCATE(boundarySEE:: boundary)
+      SELECT TYPE(boundary)
+      TYPE IS(boundarySEE)
+        CALL boundary%yield%init(tableFile)
+        CALL boundary%yield%convert(eV2J/(m_ref*v_ref**2), 1.D0)
+        boundary%electron => species(speciesID)%obj
+        !Search for the threshold energy in the table
+        DO i = 1, SIZE(boundary%yield%f)
+          IF (boundary%yield%f(i) > 0.D0) THEN
+            boundary%energyThreshold = boundary%yield%x(i)
+
+            EXIT
+
+          END IF
+
+        END DO
+
+      END SELECT
+
+    END SUBROUTINE initSEE
+
 END MODULE moduleBoundary

src/modules/moduleCollisions.f90

@@ -43,8 +43,7 @@ MODULE moduleCollisions
   TYPE, EXTENDS(collisionBinary):: collisionBinaryIonization
     REAL(8):: eThreshold !Minimum energy (non-dimensional units) required for ionization
     REAL(8):: deltaV !Change in velocity due to exchange of eThreshold
-    CLASS(speciesCharged), POINTER:: electron          !Pointer to species considerer as electrons
+    CLASS(speciesCharged), POINTER:: electron !Pointer to species considerer as electrons
-    CLASS(speciesCharged), POINTER:: electronSecondary !Pointer to species considerer as secondary electron
     CONTAINS
       PROCEDURE, PASS:: collide => collideBinaryIonization
 
@@ -165,8 +164,8 @@ MODULE moduleCollisions
 
       self%amount = amount
 
-      mass_i = species(i)%obj%m
+      mass_i = species(i)%obj%mass
-      mass_j = species(j)%obj%m
+      mass_j = species(j)%obj%mass
 
       ALLOCATE(self%collisions(1:self%amount))
 
@@ -228,8 +227,8 @@ MODULE moduleCollisions
       REAL(8):: m_i, m_j
       REAL(8), DIMENSION(1:3):: vCM, vp
 
-      m_i = part_i%species%m
+      m_i = part_i%species%mass
-      m_j = part_j%species%m
+      m_j = part_j%species%mass
       !Applies the collision
       vCM = velocityCM(part_i%weight*m_i, part_i%v, part_j%weight*m_j, part_j%v)
       vp  = vRel*randomDirectionVHS()
@@ -242,7 +241,7 @@ MODULE moduleCollisions
 
     !ELECTRON IMPACT IONIZATION
     !Inits electron impact ionization
-    SUBROUTINE initBinaryIonization(collision, crossSectionFilename, energyThreshold, electron, electronSecondary)
+    SUBROUTINE initBinaryIonization(collision, crossSectionFilename, energyThreshold, electron)
       USE moduleTable
       USE moduleRefParam
       USE moduleConstParam
@@ -254,8 +253,7 @@ MODULE moduleCollisions
       CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSectionFilename
       REAL(8), INTENT(in):: energyThreshold
       CHARACTER(:), ALLOCATABLE, INTENT(in):: electron
-      CHARACTER(:), ALLOCATABLE, OPTIONAL, INTENT(in):: electronSecondary
+      INTEGER:: electronIndex
-      INTEGER:: electronIndex, electronSecondaryIndex
 
       ALLOCATE(collisionBinaryIonization:: collision)
 
@@ -280,29 +278,12 @@ MODULE moduleCollisions
 
         CLASS DEFAULT
           CALL criticalError("Species " // sp%name // " chosen for " // &
-                             "impacting electron is not a charged species", 'initBinaryIonization')
+                             "secondary electron is not a charged species", 'initBinaryIonization')
 
         END SELECT
 
-        IF (PRESENT(electronSecondary)) THEN
-          electronSecondaryIndex = speciesName2Index(electronSecondary)
-          SELECT TYPE(sp => species(electronSecondaryIndex)%obj)
-          TYPE IS(speciesCharged)
-            collision%electronSecondary => sp
-
-          CLASS DEFAULT
-            CALL criticalError("Species " // sp%name // " chosen for " // &
-                               "secondary electron is not a charged species", 'initBinaryIonization')
-
-          END SELECT
-
-        ELSE
-          collision%electronSecondary => NULL()
-
-        END IF
-
         !momentum change per ionization process
-        collision%deltaV = sqrt(collision%eThreshold / collision%electron%m)
+        collision%deltaV = sqrt(collision%eThreshold / collision%electron%mass)
 
       END SELECT
 
@@ -326,7 +307,7 @@ MODULE moduleCollisions
       REAL(8), DIMENSION(1:3):: vChange
       TYPE(particle), POINTER:: newElectron => NULL(), remainingNeutral => NULL()
 
-      rMass = reducedMass(part_i%weight*part_i%species%m, part_j%weight*part_j%species%m)
+      rMass = reducedMass(part_i%weight*part_i%species%mass, part_j%weight*part_j%species%mass)
       eRel = rMass*vRel**2
       !Relative energy must be higher than threshold
       IF (eRel > self%eThreshold) THEN
@@ -355,12 +336,6 @@ MODULE moduleCollisions
         !Copy basic information from primary electron
         newElectron = electron
 
-        !If secondary electron species indicates, convert
-        IF (ASSOCIATED(self%electronSecondary)) THEN
-          newElectron%species => self%electronSecondary
-
-        END IF
-
         !Secondary electorn gains energy from ionization
         newElectron%v       =  vChange
 
@@ -387,7 +362,7 @@ MODULE moduleCollisions
           CALL sp%ionize(neutral)
 
         CLASS DEFAULT
-          CALL criticalError(sp%name // " is not a neutral", 'collideBinaryIonization')
+          ! CALL criticalError(sp%name // " is not a neutral", 'collideBinaryIonization')
           RETURN
 
         END SELECT

src/modules/moduleCoulomb.f90

@@ -61,7 +61,7 @@ MODULE moduleCoulomb
       self%sp_i => species(i)%obj
       self%sp_j => species(j)%obj
 
-      self%one_plus_massRatio_ij = 1.D0 + self%sp_i%m/self%sp_j%m
+      self%one_plus_massRatio_ij = 1.D0 + self%sp_i%mass/self%sp_j%mass
 
       Z_i = 0.D0
       Z_j = 0.D0
@@ -87,11 +87,11 @@ MODULE moduleCoulomb
 
       scaleFactor = (n_ref * qe**4 * ti_ref) / (eps_0**2 * m_ref**2 * v_ref**3)
 
-      self%A_i = Z_i**2*Z_j**2*self%lnCoulomb / (2.D0 * PI**2 * self%sp_i%m**2) * scaleFactor !Missing density because it's cell dependent
+      self%A_i = Z_i**2*Z_j**2*self%lnCoulomb / (2.D0 * PI**2 * self%sp_i%mass**2) * scaleFactor !Missing density because it's cell dependent
-      self%A_j = Z_j**2*Z_i**2*self%lnCoulomb / (2.D0 * PI**2 * self%sp_j%m**2) * scaleFactor !Missing density because it's cell dependent
+      self%A_j = Z_j**2*Z_i**2*self%lnCoulomb / (2.D0 * PI**2 * self%sp_j%mass**2) * scaleFactor !Missing density because it's cell dependent
 
-      self%l2_j = self%sp_j%m / 2.D0 !Missing temperature because it's cell dependent
+      self%l2_j = self%sp_j%mass / 2.D0 !Missing temperature because it's cell dependent
-      self%l2_i = self%sp_i%m / 2.D0 !Missing temperature because it's cell dependent
+      self%l2_i = self%sp_i%mass / 2.D0 !Missing temperature because it's cell dependent
 
     END SUBROUTINE initInteractionCoulomb
 

src/modules/moduleInject.f90

@@ -54,16 +54,15 @@ MODULE moduleInject
     INTEGER:: id
     CHARACTER(:), ALLOCATABLE:: name
     REAL(8):: vMod !Velocity (module)
-    REAL(8):: temperature(1:3) !Temperature
+    REAL(8):: T(1:3) !Temperature
     REAL(8):: n(1:3) !Direction of injection
     LOGICAL:: fixDirection !The injection of particles has a fix direction defined by n
     INTEGER:: nParticles !Number of particles to introduce each time step
     CLASS(speciesGeneric), POINTER:: species !Species of injection
     INTEGER:: nEdges
     INTEGER, ALLOCATABLE:: edges(:) !Array with edges
-    INTEGER, ALLOCATABLE:: particlesPerEdge(:) ! Particles per edge
+    REAL(8), ALLOCATABLE:: cumWeight(:) !Array of cummulative probability
-    REAL(8), ALLOCATABLE:: weightPerEdge(:) ! Weight per edge
+    REAL(8):: sumWeight
-    REAL(8):: surface ! Total surface of injection
     TYPE(velDistCont):: v(1:3) !Velocity distribution function in each direction
     CONTAINS
       PROCEDURE, PASS:: init         => initInject
@@ -76,7 +75,7 @@ MODULE moduleInject
 
   CONTAINS
     !Initialize an injection of particles
-    SUBROUTINE initInject(self, i, v, n, temperature, flow, units, sp, physicalSurface, particlesPerEdge)
+    SUBROUTINE initInject(self, i, v, n, T, flow, units, sp, physicalSurface)
       USE moduleMesh
       USE moduleRefParam
       USE moduleConstParam
@@ -87,29 +86,49 @@ MODULE moduleInject
 
       CLASS(injectGeneric), INTENT(inout):: self
       INTEGER, INTENT(in):: i
-      REAL(8), INTENT(in):: v, n(1:3), temperature(1:3)
+      REAL(8), INTENT(in):: v, n(1:3), T(1:3)
-      INTEGER, INTENT(in):: sp, physicalSurface, particlesPerEdge
+      INTEGER, INTENT(in):: sp, physicalSurface
       REAL(8):: tauInject
       REAL(8), INTENT(in):: flow
       CHARACTER(:), ALLOCATABLE, INTENT(in):: units
       INTEGER:: e, et
       INTEGER:: phSurface(1:mesh%numEdges)
       INTEGER:: nVolColl
-      REAL(8):: fluxPerStep = 0.D0
 
-      self%id          =  i
+      self%id         =  i
-      self%vMod        =  v / v_ref
+      self%vMod       =  v / v_ref
-      self%n           =  n / NORM2(n)
+      self%n          =  n / NORM2(n)
-      self%temperature =  temperature / T_ref
+      self%T          =  T / T_ref
+      self%species    => species(sp)%obj
+      tauInject = tau(self%species%n)
+      SELECT CASE(units)
+      CASE ("sccm")
+        !Standard cubic centimeter per minute
+        self%nParticles = INT(flow*sccm2atomPerS*tauInject*ti_ref/species(sp)%obj%weight)
+
+      CASE ("A")
+        !Input current in Ampers
+        self%nParticles = INT(flow*tauInject*ti_ref/(qe*species(sp)%obj%weight))
+
+      CASE ("part/s")
+        !Input current in Ampers
+        self%nParticles = INT(flow*tauInject*ti_ref/species(sp)%obj%weight)
+
+      CASE DEFAULT
+        CALL criticalError("No support for units: " // units, 'initInject')
+        
+      END SELECT
+      !Scale particles for different species steps
+      IF (self%nParticles == 0) CALL criticalError("The number of particles for inject is 0.", 'initInject')
+
       !Gets the edge elements from which particles are injected
       DO e = 1, mesh%numEdges
         phSurface(e) = mesh%edges(e)%obj%physicalSurface
 
       END DO
+
       self%nEdges = COUNT(phSurface == physicalSurface)
-      ALLOCATE(self%edges(1:self%nEdges))
+      ALLOCATE(inject(i)%edges(1:self%nEdges))
-      ALLOCATE(self%particlesPerEdge(1:self%nEdges))
-      ALLOCATE(self%weightPerEdge(1:self%nEdges))
       et = 0
       DO e=1, mesh%numEdges
         IF (mesh%edges(e)%obj%physicalSurface == physicalSurface) THEN
@@ -141,82 +160,15 @@ MODULE moduleInject
 
       END DO
 
-      !Calculates total area
+      !Calculates cumulative probability
-      self%surface = 0.D0
+      ALLOCATE(self%cumWeight(1:self%nEdges))
-      DO et = 1, self%nEdges
+      et = 1
-        self%surface = self%surface + mesh%edges(self%edges(et))%obj%surface
+      self%cumWeight(1) = mesh%edges(self%edges(et))%obj%weight
+      DO et = 2, self%nEdges
+        self%cumWeight(et) = mesh%edges(self%edges(et))%obj%weight + self%cumWeight(et-1)
 
       END DO
-
+      self%sumWeight = self%cumWeight(self%nEdges)
-      ! Information about species and flux
-      self%species => species(sp)%obj
-      tauInject    = tau(self%species%n)
-      ! Convert units
-      SELECT CASE(units)
-      CASE ("sccm")
-        !Standard cubic centimeter per minute
-        fluxPerStep = flow*sccm2atomPerS
-
-      CASE ("A")
-        !Current in Ampers
-        SELECT TYPE(sp => self%species)
-        CLASS IS(speciesCharged)
-          fluxPerStep = flow/(qe*abs(sp%q))
-
-        CLASS DEFAULT
-          call criticalError('Attempted to assign a flux in "A" to a species without charge.', 'initInject')
-
-        END SELECT
-
-      CASE ("Am2")
-        !Input current in Ampers per square meter
-        SELECT TYPE(sp => self%species)
-        CLASS IS(speciesCharged)
-          fluxPerStep = flow*self%surface*L_ref**2/(qe*abs(sp%q))
-
-        CLASS DEFAULT
-          call criticalError('Attempted to assign a flux in "Am2" to a species without charge.', 'initInject')
-
-        END SELECT
-
-
-      CASE ("part/s")
-        !Input current in Ampers
-        fluxPerStep = flow
-
-      CASE DEFAULT
-        CALL criticalError("No support for units: " // units, 'initInject')
-
-      END SELECT
-      fluxPerStep = fluxPerStep * tauInject * ti_ref / self%surface
-
-      !Assign particles per edge
-      IF (particlesPerEdge > 0) THEN
-        ! Particles per edge defined by the user
-        self%particlesPerEdge = particlesPerEdge
-        DO et = 1, self%nEdges
-          self%weightPerEdge(et) = fluxPerStep*mesh%edges(self%edges(et))%obj%surface / REAL(particlesPerEdge)
-
-        END DO
-
-        self%nParticles = SUM(self%particlesPerEdge)
-
-      ELSE
-        ! No particles assigned per edge, use the species weight
-        self%weightPerEdge = self%species%weight
-        DO et = 1, self%nEdges
-          self%particlesPerEdge(et) = max(1,FLOOR(fluxPerStep*mesh%edges(self%edges(et))%obj%surface / self%species%weight))
-        END DO
-
-        self%nParticles = SUM(self%particlesPerEdge)
-
-        !Rescale weight to match flux
-        self%weightPerEdge = fluxPerStep * self%surface / (real(self%nParticles))
-
-      END IF
-
-      !Scale particles for different species steps
-      IF (self%nParticles == 0) CALL criticalError("The number of particles for inject is 0.", 'initInject')
 
     END SUBROUTINE initInject
 
@@ -227,16 +179,21 @@ MODULE moduleInject
       IMPLICIT NONE
 
       INTEGER:: i
+      INTEGER, DIMENSION(1:nInject):: nMin, nMax
       
       !$OMP SINGLE
       nPartInj = 0
       DO i = 1, nInject
         IF (solver%pusher(inject(i)%species%n)%pushSpecies) THEN 
+          nMin(i)  = nPartInj + 1 
           nPartInj = nPartInj + inject(i)%nParticles
+          nMax(i)  = nPartInj
 
         END IF
 
       END DO
+      PRINT *, nMin
+      PRINT *, nMax
 
       IF (ALLOCATED(partInj)) DEALLOCATE(partInj)
       ALLOCATE(partInj(1:nPartInj))
@@ -244,30 +201,30 @@ MODULE moduleInject
 
       DO i=1, nInject
         IF (solver%pusher(inject(i)%species%n)%pushSpecies) THEN 
-          CALL inject(i)%addParticles()
+          CALL inject(i)%addParticles(nMin(i), nMax(i))
 
         END IF
       END DO
 
     END SUBROUTINE doInjects
 
-    SUBROUTINE initVelDistMaxwellian(velDist, temperature, m)
+    SUBROUTINE initVelDistMaxwellian(velDist, T, m)
       IMPLICIT NONE
 
       CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist
-      REAL(8), INTENT(in):: temperature, m
+      REAL(8), INTENT(in):: T, m
 
-      velDist = velDistMaxwellian(vTh = DSQRT(2.d0*temperature/m))
+      velDist = velDistMaxwellian(vTh = DSQRT(T/m))
 
     END SUBROUTINE initVelDistMaxwellian
 
-    SUBROUTINE initVelDistHalfMaxwellian(velDist, temperature, m)
+    SUBROUTINE initVelDistHalfMaxwellian(velDist, T, m)
       IMPLICIT NONE
 
       CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist
-      REAL(8), INTENT(in):: temperature, m
+      REAL(8), INTENT(in):: T, m
 
-      velDist = velDistHalfMaxwellian(vTh = DSQRT(2.d0*temperature/m))
+      velDist = velDistHalfMaxwellian(vTh = DSQRT(T/m))
 
     END SUBROUTINE initVelDistHalfMaxwellian
 
@@ -289,7 +246,7 @@ MODULE moduleInject
       REAL(8):: v
       v = 0.D0
 
-      v = self%vTh*randomMaxwellian()/sqrt(2.d0)
+      v = self%vTh*randomMaxwellian()
 
     END FUNCTION randomVelMaxwellian
 
@@ -302,7 +259,7 @@ MODULE moduleInject
       REAL(8):: v
       v = 0.D0
 
-      v = self%vTh*randomHalfMaxwellian()/sqrt(2.d0)
+      v = self%vTh*randomHalfMaxwellian()
 
     END FUNCTION randomVelHalfMaxwellian
 
@@ -327,82 +284,68 @@ MODULE moduleInject
       IMPLICIT NONE
 
       CLASS(injectGeneric), INTENT(in):: self
-      INTEGER, SAVE:: nMin
+      INTEGER, INTENT(in):: nMin, nMax !Min and Max index in partInj array
-      INTEGER:: i, e
+      INTEGER:: randomX
+      INTEGER:: i
       INTEGER:: n, sp
       CLASS(meshEdge), POINTER:: randomEdge
       REAL(8):: direction(1:3)
 
       !Insert particles
       !$OMP SINGLE
-      nMin = 0
+      !Assign weight to particle.
-      DO i = 1, self%id -1
+      partInj(nMin:nMax)%weight = self%species%weight
-        IF (solver%pusher(inject(i)%species%n)%pushSpecies) THEN 
+      !Particle is considered to be outside the domain
-          nMin = nMin + inject(i)%nParticles
+      partInj(nMin:nMax)%n_in = .FALSE.
-
-        END IF
-
-      END DO
-      nMin = nMin + 1
       !$OMP END SINGLE
 
       !$OMP DO
-      DO e = 1, self%nEdges
+      DO n = nMin, nMax
-        ! Select edge for injection
+        randomX = randomWeighted(self%cumWeight, self%sumWeight)
-        randomEdge => mesh%edges(self%edges(e))%obj
-        ! Inject particles in edge
-        DO i = 1, self%particlesPerEdge(e)
-          ! Index in the global partInj array
-          n = nMin - 1 + SUM(self%particlesPerEdge(1:e-1)) + i
-          !Particle is considered to be outside the domain
-          partInj(n)%n_in = .FALSE.
-          !Random position in edge
-          partInj(n)%r = randomEdge%randPos()
-          !Assign weight to particle.
-          partInj(n)%weight = self%weightPerEdge(e)
-          !Volume associated to the edge:
-          IF (ASSOCIATED(randomEdge%e1)) THEN
-            partInj(n)%cell = randomEdge%e1%n
 
-          ELSEIF (ASSOCIATED(randomEdge%e2)) THEN
+        randomEdge => mesh%edges(self%edges(randomX))%obj
-            partInj(n)%cell = randomEdge%e2%n
+        !Random position in edge
+        partInj(n)%r = randomEdge%randPos()
+        !Volume associated to the edge:
+        IF (ASSOCIATED(randomEdge%e1)) THEN
+          partInj(n)%cell = randomEdge%e1%n
 
-          ELSE
+        ELSEIF (ASSOCIATED(randomEdge%e2)) THEN
-            CALL criticalError("No Volume associated to edge", 'addParticles')
+          partInj(n)%cell = randomEdge%e2%n
 
-          END IF
+        ELSE
-          partInj(n)%cellColl = randomEdge%eColl%n
+          CALL criticalError("No Volume associated to edge", 'addParticles') 
-          sp = self%species%n
 
-          !Assign particle type
+        END IF
-          partInj(n)%species => self%species
+        partInj(n)%cellColl = randomEdge%eColl%n
+        sp = self%species%n
 
-          if (all(self%n == 0.D0)) then
+        !Assign particle type
-            direction =  randomEdge%normal
+        partInj(n)%species => self%species
 
-          else
+        direction = self%n
-            direction = self%n
 
-          end if
+        !Sample initial velocity
+        partInj(n)%v = self%vMod*direction +  (/ self%v(1)%obj%randomVel(), &
+                                                 self%v(2)%obj%randomVel(), &
+                                                 self%v(3)%obj%randomVel() /)
 
-          partInj(n)%v = 0.D0
+        !For each direction, velocities have to agree with the direction of injection
+        DO i = 1, 3
+          DO WHILE (partInj(n)%v(i)*direction(i) < 0)
+            partInj(n)%v(i) = self%vMod*direction(i) + self%v(i)%obj%randomVel()
 
-          do while(dot_product(partInj(n)%v, direction) <= 0.d0)
+          END DO
-            partInj(n)%v = self%vMod*direction +  (/ self%v(1)%obj%randomVel(), &
-                                                     self%v(2)%obj%randomVel(), &
-                                                     self%v(3)%obj%randomVel() /)
-          end do
-
-
-          !Obtain natural coordinates of particle in cell
-          partInj(n)%Xi = mesh%cells(partInj(n)%cell)%obj%phy2log(partInj(n)%r)
-          !Push new particle with the minimum time step
-          CALL solver%pusher(sp)%pushParticle(partInj(n), tau(sp))
-          !Assign cell to new particle
-          CALL solver%updateParticleCell(partInj(n))
 
         END DO
 
+        !Obtain natural coordinates of particle in cell
+        partInj(n)%Xi = mesh%cells(partInj(n)%cell)%obj%phy2log(partInj(n)%r)
+        !Push new particle with the minimum time step
+        CALL solver%pusher(sp)%pushParticle(partInj(n), tau(sp))
+        !Assign cell to new particle
+        CALL solver%updateParticleCell(partInj(n))
+
       END DO
       !$OMP END DO
 

src/modules/moduleProbe.f90

@@ -27,7 +27,7 @@ MODULE moduleProbe
 
   CONTAINS
     !Functions for probeDistFunc type
-    SUBROUTINE init(self, id, speciesName, r, v1, v2, v3, points, everyTimeStep)
+    SUBROUTINE init(self, id, speciesName, r, v1, v2, v3, points, timeStep)
       USE moduleCaseParam
       USE moduleRefParam
       USE moduleSpecies
@@ -41,7 +41,7 @@ MODULE moduleProbe
       REAL(8), INTENT(in):: r(1:3)
       REAL(8), INTENT(in):: v1(1:2), v2(1:2), v3(1:2)
       INTEGER, INTENT(in):: points(1:3)
-      REAL(8), INTENT(in):: everyTimeStep
+      REAL(8), INTENT(in):: timeStep
       INTEGER:: sp, i
       REAL(8):: dv(1:3)
 
@@ -91,17 +91,17 @@ MODULE moduleProbe
                       1:self%nv(3)))
 
       !Number of iterations between output
-      IF (everyTimeStep == 0.D0) THEN
+      IF (timeStep == 0.D0) THEN
         self%every = 1
 
       ELSE
-        self%every = NINT(everyTimeStep/ tauMin / ti_ref)
+        self%every = NINT(timeStep/ tauMin / ti_ref)
 
       END IF
       
       !Maximum radius
       !TODO: Make this an input parameter
-      self%maxR = 1.D-2/L_ref
+      self%maxR = 1.D0
 
       !Init the probe lock
       CALL OMP_INIT_LOCK(self%lock)
@@ -148,7 +148,7 @@ MODULE moduleProbe
         deltaR = NORM2(self%r - part%r)
 
         !Only include particle if it is inside the maximum radius
-        ! IF (deltaR < self%maxR) THEN
+        IF (deltaR < self%maxR) THEN
           !find lower index for all dimensions
           CALL self%findLowerIndex(part%v, i, j, k, inside)
 
@@ -162,40 +162,40 @@ MODULE moduleProbe
             fk  = self%vk(k+1) - part%v(3)
             fk1 = part%v(3)    - self%vk(k)
 
-            weight = part%weight * DEXP(-deltaR/self%maxR)
+            ! weight = part%weight * DEXP(deltaR/self%maxR)
-            ! weight = part%weight
+            weight = part%weight
 
             !Lock the probe
             CALL OMP_SET_LOCK(self%lock)
             
             !Assign particle weight to distribution function
-            self%f(i  , j  , k  ) = self%f(i  , j  , k  ) + fi  * fj  * fk  * weight
+            self%f(i  , j  , k  ) = fi  * fj  * fk  * weight
-            self%f(i+1, j  , k  ) = self%f(i+1, j  , k  ) + fi1 * fj  * fk  * weight
+            self%f(i+1, j  , k  ) = fi1 * fj  * fk  * weight
-            self%f(i  , j+1, k  ) = self%f(i  , j+1, k  ) + fi  * fj1 * fk  * weight
+            self%f(i  , j+1, k  ) = fi  * fj1 * fk  * weight
-            self%f(i+1, j+1, k  ) = self%f(i+1, j+1, k  ) + fi1 * fj1 * fk  * weight
+            self%f(i+1, j+1, k  ) = fi1 * fj1 * fk  * weight
-            self%f(i  , j  , k+1) = self%f(i  , j  , k+1) + fi  * fj  * fk1 * weight
+            self%f(i  , j  , k+1) = fi  * fj  * fk1 * weight
-            self%f(i+1, j  , k+1) = self%f(i+1, j  , k+1) + fi1 * fj  * fk1 * weight
+            self%f(i+1, j  , k+1) = fi1 * fj  * fk1 * weight
-            self%f(i  , j+1, k+1) = self%f(i  , j+1, k+1) + fi  * fj1 * fk1 * weight
+            self%f(i  , j+1, k+1) = fi  * fj1 * fk1 * weight
-            self%f(i+1, j+1, k+1) = self%f(i+1, j+1, k+1) + fi1 * fj1 * fk1 * weight
+            self%f(i+1, j+1, k+1) = fi1 * fj1 * fk1 * weight
 
             !Unlock the probe
             CALL OMP_UNSET_LOCK(self%lock)
 
           END IF
 
-        ! END IF
+        END IF
 
       END IF
           
     END SUBROUTINE calculate
 
-    SUBROUTINE output(self)
+    SUBROUTINE output(self, t)
       USE moduleOutput
       USE moduleRefParam
-      USE moduleCaseParam, ONLY: timeStep
       IMPLICIT NONE
 
       CLASS(probeDistFunc), INTENT(inout):: self
+      INTEGER, INTENT(in):: t
       CHARACTER (LEN=iterationDigits):: tstring
       CHARACTER (LEN=3):: pstring
       CHARACTER(:), ALLOCATABLE:: filename
@@ -204,14 +204,14 @@ MODULE moduleProbe
       !Divide by the velocity cube volume
       self%f = self%f * self%dvInv
 
-      WRITE(tstring, iterationFormat) timeStep
+      WRITE(tstring, iterationFormat) t
       WRITE(pstring, "(I3.3)") self%id
       fileName='Probe_' // tstring// '_f_' // pstring // '.dat'
       WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
       OPEN (10, file = path // folder // '/' //  fileName)
       WRITE(10, "(A1, 1X, A)")               "# ", self%species%name
       WRITE(10, "(A6, 3(ES15.6E3), A2)") "# r = ", self%r(:)*L_ref, " m"
-      WRITE(10, "(A6, ES15.6E3, A2)")    "# t = ", REAL(timeStep)*tauMin*ti_ref, " s"
+      WRITE(10, "(A6, ES15.6E3, A2)")    "# t = ", REAL(t)*tauMin*ti_ref, " s"
       WRITE(10, "(A1, A19, 3(A20))")     "#", "v1 (m s^-1)", "v2 (m s^-1)", "v3 (m s^-1)", "f"
       DO i = 1, self%nv(1)
         DO j = 1, self%nv(2)
@@ -252,15 +252,15 @@ MODULE moduleProbe
 
     END SUBROUTINE doProbes
 
-    SUBROUTINE outputProbes()
+    SUBROUTINE outputProbes(t)
-      USE moduleCaseParam, ONLY: timeStep
       IMPLICIT NONE
 
+      INTEGER, INTENT(in):: t
       INTEGER:: i
 
       DO i = 1, nProbes
         IF (probe(i)%update) THEN
-          CALL probe(i)%output()
+          CALL probe(i)%output(t)
 
         END IF
 
@@ -268,15 +268,15 @@ MODULE moduleProbe
 
     END SUBROUTINE outputProbes
 
-    SUBROUTINE resetProbes()
+    SUBROUTINE resetProbes(t)
-      USE moduleCaseParam, ONLY: timeStep
       IMPLICIT NONE
 
+      INTEGER, INTENT(in):: t
       INTEGER:: i
 
       DO i = 1, nProbes
         probe(i)%f = 0.D0
-        probe(i)%update = timeStep == tFinal .OR. timeStep == tInitial .OR. MOD(timeStep, probe(i)%every) == 0
+        probe(i)%update = t == tFinal .OR. t == tInitial .OR. MOD(t, probe(i)%every) == 0
 
       END DO
 

src/modules/moduleSpecies.f90

@@ -4,12 +4,54 @@ MODULE moduleSpecies
   USE OMP_LIB
   IMPLICIT NONE
   
+  !Basic type that defines a macro-particle
+  TYPE:: particle
+    REAL(8):: r(1:3)                         !Position
+    REAL(8):: v(1:3)                         !Velocity
+    CLASS(speciesGeneric), POINTER:: species !Pointer to particle's species
+    INTEGER:: cell                           !Index of element in which the particle is located TODO: Make these pointers
+    INTEGER:: cellColl                       !Index of element in which the particle is located in the Collision Mesh
+    REAL(8):: Xi(1:3)                        !Logical coordinates of particle in cell
+    LOGICAL:: n_in                           !Flag that indicates if a particle is in the domain
+    REAL(8):: weight=0.D0                    !weight of particle
+
+  END TYPE particle
+
+  !Wrapper to store the particles per species
+  TYPE:: particleArray
+    TYPE(particle), ALLOCATABLE, DIMENSION(:):: p
+
+  END TYPE particleArray
+
+  !Array of pointers for the particles to be pushed
+  TYPE:: particlePointer
+    TYPE(particle), POINTER:: p
+
+  END TYPE particlePointer
+
+  !Arrays that contain the particles
+  TYPE(particleArray),   ALLOCATABLE, TARGET, DIMENSION(:):: particles           !array of particles in the domain
+  TYPE(particle),        ALLOCATABLE, TARGET, DIMENSION(:):: particlesInjection  !array of inject particles
+  TYPE(particlePointer), ALLOCATABLE,         DIMENSION(:):: particlesToPush     !particles pushed in each iteration
+  !Integers to track number of particles in domain
+  INTEGER, ALLOCATABLE, DIMENSION(:):: nParticles
+  INTEGER:: nParticlesTotal
+  !Number of injected particles
+  INTEGER, ALLOCATABLE, DIMENSION(:):: nParticlesInject
+  INTEGER:: nPariclesInjectTotal
+
+  !Generic species type
   TYPE, ABSTRACT:: speciesGeneric
-    CHARACTER(:), ALLOCATABLE:: name
+    INTEGER:: n=0 !Index of species
-    REAL(8):: m=0.D0, weight=0.D0, qm=0.D0
+    CHARACTER(:), ALLOCATABLE:: name !Name of species
-    INTEGER:: n=0
+    !Mass, default weight of species and charge over mass
+    REAL(8):: mass=0.D0, weight=0.D0, qm=0.D0
+    INTEGER:: every !How many interations between advancing the species
+    LOGICAL:: pushSpecies !Boolean to indicate if the species is moved in the iteration
+
   END TYPE speciesGeneric
 
+  !Neutral species
   TYPE, EXTENDS(speciesGeneric):: speciesNeutral
     CLASS(speciesGeneric), POINTER:: ion => NULL()
     CONTAINS
@@ -17,6 +59,7 @@ MODULE moduleSpecies
 
   END TYPE speciesNeutral
 
+  !Charged species
   TYPE, EXTENDS(speciesGeneric):: speciesCharged
     REAL(8):: q=0.D0
     CLASS(speciesGeneric), POINTER:: ion => NULL(), neutral => NULL()
@@ -26,34 +69,17 @@ MODULE moduleSpecies
 
   END TYPE speciesCharged
 
+  !Wrapper for species
   TYPE:: speciesCont
     CLASS(speciesGeneric), ALLOCATABLE:: obj
 
   END TYPE
 
+  !Number of species
   INTEGER:: nSpecies
+  !Array for species
   TYPE(speciesCont), ALLOCATABLE, TARGET:: species(:)
 
-  TYPE particle
-    REAL(8):: r(1:3) !Position
-    REAL(8):: v(1:3) !Velocity
-    CLASS(speciesGeneric), POINTER:: species !Pointer to species associated with this particle
-    INTEGER:: cell !Index of element in which the particle is located
-    INTEGER:: cellColl !Index of element in which the particle is located in the Collision Mesh
-    REAL(8):: Xi(1:3) !Logical coordinates of particle in element e_p.
-    LOGICAL:: n_in !Flag that indicates if a particle is in the domain
-    REAL(8):: weight=0.D0 !weight of particle
-
-  END TYPE particle
-
-  !Number of old particles
-  INTEGER:: nPartOld
-  !Number of injected particles
-  INTEGER:: nPartInj
-  !Arrays that contain the particles
-  TYPE(particle), ALLOCATABLE, DIMENSION(:), TARGET:: partOld !array of particles from previous iteration
-  TYPE(particle), ALLOCATABLE, DIMENSION(:), TARGET:: partInj !array of inject particles
-
   CONTAINS
     FUNCTION speciesName2Index(speciesName) RESULT(sp)
       USE moduleErrors

src/modules/output/moduleOutput.f90

@@ -22,6 +22,7 @@ MODULE moduleOutput
 
   !Type for EM data in node
   TYPE emNode
+    CHARACTER(:), ALLOCATABLE:: type
     REAL(8):: phi
     REAL(8):: B(1:3)
 
@@ -147,8 +148,8 @@ MODULE moduleOutput
         formatValues%density     = rawValues%den*tempVol
         formatValues%velocity(:) = tempVel
         IF (tensorTrace(tensorTemp) > 0.D0) THEN
-         formatValues%pressure    = speciesIn%m*tensorTrace(tensorTemp)*tempVol/3.D0
+          formatValues%pressure    = speciesIn%mass*tensorTrace(tensorTemp)*tempVol/3.D0
-         formatValues%temperature = formatValues%pressure/(formatValues%density*kb)
+          formatValues%temperature = formatValues%pressure/(formatValues%density*kb)
 
         END IF
       END IF
@@ -159,12 +160,12 @@ MODULE moduleOutput
 
     END SUBROUTINE calculateOutput
 
-    SUBROUTINE printTime(first)
+    SUBROUTINE printTime(t, first)
       USE moduleSpecies
       USE moduleCompTime
-      USE moduleCaseParam, ONLY: timeStep
       IMPLICIT NONE
 
+      INTEGER, INTENT(in):: t
       LOGICAL, INTENT(in), OPTIONAL:: first
       CHARACTER(:), ALLOCATABLE:: fileName
 
@@ -186,7 +187,7 @@ MODULE moduleOutput
 
         OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write')
 
-        WRITE (20, "(I10, I10, 7(ES20.6E3))") timeStep, nPartOld, tStep, tPush, tReset, tColl, tCoul, tWeight, tEMField
+        WRITE (20, "(I10, I10, 7(ES20.6E3))") t, nParticlesTotal, tStep, tPush, tReset, tColl, tCoul, tWeight, tEMField
 
         CLOSE(20)
 

src/modules/solver/electromagnetic/moduleEM.f90

@@ -1,202 +1,55 @@
 !Module to solve the electromagnetic field
 MODULE moduleEM
-  USE moduleMesh
-  USE moduleTable
   IMPLICIT NONE
 
-  ! Generic type for electromagnetic boundary conditions
+  TYPE:: boundaryEM
-  TYPE, PUBLIC, ABSTRACT:: boundaryEMGeneric
+    CHARACTER(:), ALLOCATABLE:: typeEM
-    INTEGER:: nNodes
+    INTEGER:: physicalSurface
-    TYPE(meshNodePointer), ALLOCATABLE:: nodes(:)
-
-    CONTAINS
-      PROCEDURE(applyEM_interface), DEFERRED, PASS:: apply
-
-  END TYPE boundaryEMGeneric
-
-  ABSTRACT INTERFACE
-    ! Apply boundary condition to the load vector for the Poission equation
-    SUBROUTINE applyEM_interface(self, vectorF)
-      IMPORT boundaryEMGeneric
-      CLASS(boundaryEMGeneric), INTENT(in):: self
-      REAL(8), INTENT(inout):: vectorF(:)
-
-    END SUBROUTINE applyEM_interface
-
-  END INTERFACE
-
-  TYPE, EXTENDS(boundaryEMGeneric):: boundaryEMDirichlet
     REAL(8):: potential
 
     CONTAINS
-      ! boundaryEMGeneric DEFERRED PROCEDURES
+      PROCEDURE, PASS:: apply
-      PROCEDURE, PASS:: apply => applyDirichlet
 
-  END TYPE boundaryEMDirichlet
+  END TYPE boundaryEM
-
-  TYPE, EXTENDS(boundaryEMGeneric):: boundaryEMDirichletTime
-    REAL(8):: potential
-    TYPE(table1D):: temporalProfile
-
-    CONTAINS
-      ! boundaryEMGeneric DEFERRED PROCEDURES
-      PROCEDURE, PASS:: apply => applyDirichletTime
-
-  END TYPE boundaryEMDirichletTime
-
-  ! Container for boundary conditions
-  TYPE:: boundaryEMCont
-    CLASS(boundaryEMGeneric), ALLOCATABLE:: obj
-
-  END TYPE boundaryEMCont
 
   INTEGER:: nBoundaryEM
-  TYPE(boundaryEMCont), ALLOCATABLE:: boundaryEM(:)
+  TYPE(boundaryEM), ALLOCATABLE:: boundEM(:)
 
   !Information of charge and reference parameters for rho vector
   REAL(8), ALLOCATABLE:: qSpecies(:)
 
   CONTAINS
-    SUBROUTINE findNodes(self, physicalSurface)
+    !Apply boundary conditions to the K matrix for Poisson's equation
+    SUBROUTINE apply(self, edge)
       USE moduleMesh
       IMPLICIT NONE
 
-      CLASS(boundaryEMGeneric), INTENT(inout):: self
+      CLASS(boundaryEM), INTENT(in):: self
-      INTEGER, INTENT(in):: physicalSurface
+      CLASS(meshEdge):: edge
-      CLASS(meshEdge), POINTER:: edge
+      INTEGER:: nNodes
-      INTEGER, ALLOCATABLE:: nodes(:), nodesEdge(:)
+      INTEGER, ALLOCATABLE:: nodes(:)
-      INTEGER:: nNodes, nodesNew
+      INTEGER:: n
-      INTEGER:: e, n
 
-      !Temporal array to hold nodes
+      nNodes = edge%nNodes
-      ALLOCATE(nodes(0))
+      nodes = edge%getNodes(nNodes)
 
-      ! Loop thorugh the edges and identify those that are part of the boundary
-      DO e = 1, mesh%numEdges
-        edge => mesh%edges(e)%obj
-        IF (edge%physicalSurface == physicalSurface) THEN
-          ! Edge is of the right boundary index
-          ! Get nodes in the edge
-          nNodes = edge%nNodes
-          nodesEdge = edge%getNodes(nNodes)
-          ! Collect all nodes that are not already in the temporal array
-          DO n = 1, nNodes
-            IF (ANY(nodes == nodesEdge(n))) THEN
-              ! Node already in array, skip
-              CYCLE
-
-            ELSE
-              ! If not, add element to array of nodes
-              nodes = [nodes, nodesEdge(n)]
-
-            END IF
-
-          END DO
-
-        END IF
-
-      END DO
-
-      ! Point boundary to nodes
-      nNodes = SIZE(nodes)
-      ALLOCATE(self%nodes(nNodes))
-      self%nNodes = nNodes
       DO n = 1, nNodes
-        self%nodes(n)%obj => mesh%nodes(nodes(n))%obj
+        SELECT CASE(self%typeEM)
+        CASE ("dirichlet")
+          mesh%K(nodes(n), :)        = 0.D0
+          mesh%K(nodes(n), nodes(n)) = 1.D0
+        
+          mesh%nodes(nodes(n))%obj%emData%type = self%typeEM
+          mesh%nodes(nodes(n))%obj%emData%phi  = self%potential
+
+        END SELECT
 
       END DO
 
-    END SUBROUTINE findNodes
+    END SUBROUTINE
-
-    ! Initialize Dirichlet boundary condition
-    SUBROUTINE initDirichlet(self, physicalSurface, potential)
-      USE moduleRefParam, ONLY: Volt_ref
-      IMPLICIT NONE
-
-      CLASS(boundaryEMGeneric), ALLOCATABLE, INTENT(out):: self
-      INTEGER, INTENT(in):: physicalSurface
-      REAL(8), INTENT(in):: potential
-
-      ! Allocate boundary edge
-      ALLOCATE(boundaryEMDirichlet:: self)
-
-      SELECT TYPE(self)
-      TYPE IS(boundaryEMDirichlet)
-        self%potential = potential / Volt_ref
-
-        CALL findNodes(self, physicalSurface)
-
-      END SELECT
-
-    END SUBROUTINE initDirichlet
-
-    ! Initialize Dirichlet boundary condition
-    SUBROUTINE initDirichletTime(self, physicalSurface, potential, temporalProfile)
-      USE moduleRefParam, ONLY: Volt_ref, ti_ref
-      IMPLICIT NONE
-
-      CLASS(boundaryEMGeneric), ALLOCATABLE, INTENT(out):: self
-      INTEGER, INTENT(in):: physicalSurface
-      REAL(8), INTENT(in):: potential
-      CHARACTER(:), ALLOCATABLE, INTENT(in):: temporalProfile
-
-      ! Allocate boundary edge
-      ALLOCATE(boundaryEMDirichletTime:: self)
-
-      SELECT TYPE(self)
-      TYPE IS(boundaryEMDirichletTime)
-        self%potential = potential / Volt_ref
-
-        CALL findNodes(self, physicalSurface)
-
-        CALL self%temporalProfile%init(temporalProfile)
-
-        CALL self%temporalProfile%convert(1.D0/ti_ref, 1.D0)
-
-      END SELECT
-
-    END SUBROUTINE initDirichletTime
-
-    !Apply Dirichlet boundary condition to the poisson equation
-    SUBROUTINE applyDirichlet(self, vectorF)
-      USE moduleMesh
-      IMPLICIT NONE
-
-      CLASS(boundaryEMDirichlet), INTENT(in):: self
-      REAL(8), INTENT(inout):: vectorF(:)
-      INTEGER:: n, ni
-
-      DO n = 1, self%nNodes
-        self%nodes(n)%obj%emData%phi  = self%potential
-        vectorF(self%nodes(n)%obj%n)  = self%nodes(n)%obj%emData%phi
-
-      END DO
-
-    END SUBROUTINE applyDirichlet
-
-    !Apply Dirichlet boundary condition with time temporal profile
-    SUBROUTINE applyDirichletTime(self, vectorF)
-      USE moduleMesh
-      USE moduleCaseParam, ONLY: timeStep, tauMin
-      IMPLICIT NONE
-
-      CLASS(boundaryEMDirichletTime), INTENT(in):: self
-      REAL(8), INTENT(inout):: vectorF(:)
-      REAL(8):: timeFactor
-      INTEGER:: n, ni
-
-      timeFactor = self%temporalProfile%get(DBLE(timeStep)*tauMin)
-
-      DO n = 1, self%nNodes
-        self%nodes(n)%obj%emData%phi  = self%potential * timeFactor
-        vectorF(self%nodes(n)%obj%n)  = self%nodes(n)%obj%emData%phi
-
-      END DO
-
-    END SUBROUTINE applyDirichletTime
 
     !Assemble the source vector based on the charge density to solve Poisson's equation
-    SUBROUTINE assembleSourceVector(vectorF, n_e)
+    SUBROUTINE assembleSourceVector(vectorF)
       USE moduleMesh
       USE moduleRefParam
       IMPLICIT NONE
@@ -205,9 +58,8 @@ MODULE moduleEM
       REAL(8), ALLOCATABLE:: localF(:)
       INTEGER, ALLOCATABLE:: nodes(:)
       REAL(8), ALLOCATABLE:: rho(:)
-      REAL(8), INTENT(in), OPTIONAL:: n_e(1:mesh%numNodes)
       INTEGER:: nNodes
-      INTEGER:: e, i, ni, b
+      INTEGER:: e, i, ni
       CLASS(meshNode), POINTER:: node
 
       !$OMP SINGLE
@@ -225,10 +77,6 @@ MODULE moduleEM
           ni = nodes(i)
           node => mesh%nodes(ni)%obj
           rho(i) = DOT_PRODUCT(qSpecies(:), node%output(:)%den/(vol_ref*node%v*n_ref))
-          IF (PRESENT(n_e)) THEN
-            rho(i) = rho(i) - n_e(i)
-
-          END IF
 
         END DO
 
@@ -249,12 +97,18 @@ MODULE moduleEM
       !$OMP END DO
 
       !Apply boundary conditions
-      !$OMP SINGLE
+      !$OMP DO
-      do b = 1, nBoundaryEM
+      DO i = 1, mesh%numNodes
-        call boundaryEM(b)%obj%apply(vectorF)
+        node => mesh%nodes(i)%obj
 
-      end do
+        SELECT CASE(node%emData%type)
-      !$OMP END SINGLE
+        CASE ("dirichlet")
+          vectorF(i) = node%emData%phi
+
+        END SELECT
+
+      END DO
+      !$OMP END DO
 
     END SUBROUTINE assembleSourceVector
 
@@ -302,86 +156,4 @@ MODULE moduleEM
 
     END SUBROUTINE solveElecField
 
-    FUNCTION BoltzmannElectron(phi, n) RESULT(n_e)
-      USE moduleRefParam
-      USE moduleConstParam
-      IMPLICIT NONE
-
-      INTEGER, INTENT(in):: n
-      REAL(8), INTENT(in):: phi(1:n)
-      REAL(8):: n_e(1:n)
-      REAL(8):: n_e0 = 1.0D16, phi_0 = -500.0D0, T_e = 11604.0
-      INTEGER:: i
-
-      n_e = n_e0 / n_ref * exp(qe * (phi*Volt_ref - phi_0) / (kb * T_e))
-
-      RETURN
-
-    END FUNCTION BoltzmannElectron
-
-    SUBROUTINE solveElecFieldBoltzmann()
-      USE moduleMesh
-      USE moduleErrors
-      IMPLICIT NONE
-
-      INTEGER, SAVE:: INFO
-      INTEGER:: n
-      REAL(8), ALLOCATABLE, SAVE:: tempF(:)
-      REAL(8), ALLOCATABLE, SAVE:: n_e(:), phi_old(:), phi(:)
-      INTEGER:: k
-      EXTERNAL:: dgetrs
-
-      !$OMP SINGLE
-      ALLOCATE(tempF(1:mesh%numNodes))
-      ALLOCATE(n_e(1:mesh%numNodes))
-      ALLOCATE(phi_old(1:mesh%numNodes))
-      ALLOCATE(phi(1:mesh%numNodes))
-      !$OMP END SINGLE
-
-      !$OMP DO
-      DO n = 1, mesh%numNodes
-        phi_old(n) = mesh%nodes(n)%obj%emData%phi
-
-      END DO
-      !$OMP END DO
-
-      !$OMP SINGLE
-      DO k = 1, 100
-        n_e = BoltzmannElectron(phi_old, mesh%numNodes)
-        CALL assembleSourceVector(tempF, n_e)
-
-        CALL dgetrs('N', mesh%numNodes, 1, mesh%K, mesh%numNodes, &
-                                mesh%IPIV, tempF,  mesh%numNodes, info)
-        phi = tempF
-
-        PRINT *, MAXVAL(n_e), MINVAL(n_e)
-        PRINT *, MAXVAL(phi), MINVAL(phi)
-        PRINT*, k, "diff = ", MAXVAL(ABS(phi - phi_old))
-        phi_old = phi
-
-      END DO
-      !$OMP END SINGLE
-
-      IF (info == 0) THEN
-        !Suscessful resolution of Poission equation
-        !$OMP DO
-        DO n = 1, mesh%numNodes
-          mesh%nodes(n)%obj%emData%phi = phi_old(n)
-
-        END DO
-        !$OMP END DO
-
-      ELSE
-        !$OMP SINGLE
-        CALL criticalError('Poisson equation failed', 'solveElecFieldBoltzmann')
-        !$OMP END SINGLE
-
-      END IF
-
-      !$OMP SINGLE
-      DEALLOCATE(tempF, n_e, phi_old, phi)
-      !$OMP END SINGLE
-
-    END SUBROUTINE solveElecFieldBoltzmann
-
 END MODULE moduleEM

src/modules/solver/moduleSolver.f90

@@ -3,7 +3,7 @@ MODULE moduleSolver
 
   !Generic type for pusher of particles
   TYPE, PUBLIC:: pusherGeneric
-    PROCEDURE(push_interafece), POINTER, NOPASS:: pushParticle => NULL()
+    PROCEDURE(push_interface), POINTER, NOPASS:: pushParticle => NULL()
     !Boolean to indicate if the species is moved in the iteration
     LOGICAL:: pushSpecies
     !How many interations between advancing the species
@@ -29,13 +29,13 @@ MODULE moduleSolver
 
   INTERFACE
     !Push a particle
-    PURE SUBROUTINE push_interafece(part, tauIn)
+    PURE SUBROUTINE push_interface(part, tauIn)
       USE moduleSpecies
 
       TYPE(particle), INTENT(inout):: part
       REAL(8), INTENT(in):: tauIn
 
-    END SUBROUTINE push_interafece
+    END SUBROUTINE push_interface
 
     !Solve the electromagnetic field
     SUBROUTINE solveEM_interface()
@@ -123,7 +123,12 @@ MODULE moduleSolver
       END SELECT
 
       self%pushSpecies = .FALSE.
-      self%every = INT(tauSp/tau)
+      self%every = FLOOR(tauSp/tau)
+      !Correct value if not fulfilled
+      IF (tau*REAL(self%every) < tauSp) THEN
+        self%every = self%every + 1
+
+      END IF
 
     END SUBROUTINE initPusher
 
@@ -138,9 +143,6 @@ MODULE moduleSolver
       CASE('Electrostatic','ConstantB')
         self%solveEM => solveElecField
 
-      CASE('ElectrostaticBoltzmann')
-        self%solveEM => solveElecFieldBoltzmann
-
       END SELECT
 
     END SUBROUTINE initEM
@@ -167,24 +169,75 @@ MODULE moduleSolver
       USE moduleMesh
       IMPLICIT NONE
 
-      INTEGER:: n
+      INTEGER:: p
-      INTEGER:: sp
+      INTEGER:: s, sp
+      INTEGER:: e
 
-      !$OMP DO
+      !$OMP SINGLE
-      DO n = 1, nPartOld
+      !Allocate the array of particles to push
-        !Select species type
+      nSpeciesToPush = COUNT(solver%pusher(:)%pushSpecies)
-        sp = partOld(n)%species%n
+      ALLOCATE(particlesToPush(1:nSpeciesToPush))
+      ALLOCATE(nPartOldToPush(1:nSpeciesToPush))
+      !Point the arrays to the location of the particles
+      sp = 0
+      DO s = 1, nSpecies
         !Checks if the species sp is update this iteration
         IF (solver%pusher(sp)%pushSpecies) THEN
-          !Push particle
+          sp = sp + 1
-          CALL solver%pusher(sp)%pushParticle(partOld(n), tau(sp))
+          particlesToPush(sp)%p = partOld(s)%p
-          !Find cell in wich particle reside
+          nPartOldToPush(sp)    = nPartOld(s)
-          CALL solver%updateParticleCell(partOld(n))
 
         END IF
 
       END DO
-      !$OMP END DO
+
+      !Delete list of particles in cells for collisions if the particle is pushed
+      IF (listInCells) THEN
+        DO e = 1, mesh%numCells
+          DO s = 1, nSpecies
+            IF (solver%pusher(s)%pushSpecies) THEN 
+              CALL mesh%cells(e)%obj%listPart_in(s)%erase()
+              mesh%cells(e)%obj%totalWeight(s) = 0.D0
+
+            END IF
+
+          END DO
+
+        END DO
+
+      END IF
+
+      !Erase the list of particles inside the cell in coll mesh if the particle is pushed
+      IF (doubleMesh) THEN
+        DO e = 1, meshColl%numCells
+          DO s = 1, nSpecies
+            IF (solver%pusher(s)%pushSpecies) THEN 
+              CALL meshColl%cells(e)%obj%listPart_in(s)%erase()
+              meshColl%cells(e)%obj%totalWeight(s) = 0.D0
+
+            END IF
+
+          END DO
+
+        END DO
+    
+      END IF
+      !$OMP END SINGLE
+
+      !Now, push particles
+      !$OMP DO
+      DO sp = 1, nSpeciesToPush
+        DO p = 1, nPartOldToPush(sp)
+          !Push particle
+          CALL solver%pusher(sp)%pushParticle(particlesToPush(sp)%p(p), tau(sp))
+          !Find cell in which particle reside
+          CALL solver%updateParticleCell(particlesToPush(sp)%p(p))
+
+        END DO
+
+      END DO
+      !$END OMP DO
+
 
     END SUBROUTINE doPushes
 
@@ -245,7 +298,10 @@ MODULE moduleSolver
       !$OMP SECTION
       nOldIn = 0
       IF (ALLOCATED(partOld)) THEN
-        nOldIn = COUNT(partOld%n_in)
+        DO s = 1, nSpecies
+          nOldIn = nOldin + COUNT(partOld(s)%p(:)%n_in)
+
+        END DO
 
       END IF
       !$OMP SECTION
@@ -322,40 +378,6 @@ MODULE moduleSolver
 
       END DO
 
-      !$OMP SECTION
-      !Erase the list of particles inside the cell if particles have been pushed
-      IF (listInCells) THEN
-        DO s = 1, nSpecies
-          DO e = 1, mesh%numCells
-            IF (solver%pusher(s)%pushSpecies) THEN 
-              CALL mesh%cells(e)%obj%listPart_in(s)%erase()
-              mesh%cells(e)%obj%totalWeight(s) = 0.D0
-
-            END IF
-
-          END DO
-
-        END DO
-
-      END IF
-
-      !$OMP SECTION
-      !Erase the list of particles inside the cell in coll mesh
-      IF (doubleMesh) THEN
-        DO s = 1, nSpecies
-          DO e = 1, meshColl%numCells
-            IF (solver%pusher(s)%pushSpecies) THEN 
-              CALL meshColl%cells(e)%obj%listPart_in(s)%erase()
-              meshColl%cells(e)%obj%totalWeight(s) = 0.D0
-
-            END IF
-
-          END DO
-
-        END DO
-      
-      END IF
-
       !$OMP END SECTIONS
 
       !$OMP SINGLE
@@ -370,14 +392,17 @@ MODULE moduleSolver
       USE moduleMesh
       IMPLICIT NONE
 
-      INTEGER:: n
+      INTEGER:: n, sp
       CLASS(meshCell), POINTER:: cell
 
       !Loops over the particles to scatter them
       !$OMP DO
-      DO n = 1, nPartOld
+      DO sp = 1, nSpeciesToPush
-        cell => mesh%cells(partOld(n)%cell)%obj
+        DO n = 1, nPartOldToPush(sp)
-        CALL cell%scatter(cell%nNodes, partOld(n))
+          cell => mesh%cells(particlesToPush(sp)%p(n)%cell)%obj
+          CALL cell%scatter(cell%nNodes, particlesToPush(sp)%p(n))
+
+        END DO
 
       END DO
       !$OMP END DO
@@ -494,46 +519,52 @@ MODULE moduleSolver
     END SUBROUTINE updateParticleCell
 
     !Update the information about if a species needs to be moved this iteration
-    SUBROUTINE updatePushSpecies(self)
+    SUBROUTINE updatePushSpecies(self, t)
       USE moduleSpecies
-      USE moduleCaseparam, ONLY: timeStep
       IMPLICIT NONE
 
       CLASS(solverGeneric), INTENT(inout):: self
+      INTEGER, INTENT(in):: t
       INTEGER:: s
 
       DO s=1, nSpecies
-        self%pusher(s)%pushSpecies = MOD(timeStep, self%pusher(s)%every) == 0
+        self%pusher(s)%pushSpecies = MOD(t, self%pusher(s)%every) == 0
 
       END DO
 
     END SUBROUTINE updatePushSpecies
 
     !Output the different data and information
-    SUBROUTINE doOutput()
+    SUBROUTINE doOutput(t)
       USE moduleMesh
       USE moduleOutput
       USE moduleSpecies
       USE moduleCompTime
       USE moduleProbe
-      USE moduleCaseParam, ONLY: timeStep
       IMPLICIT NONE
 
-      CALL outputProbes()
+      INTEGER, INTENT(in):: t
+
+      IF (t == tInitial) THEN
+        CALL SYSTEM('git rev-parse HEAD > ' // path // folder // '/' // 'fpack_commit.txt')
+
+      END IF
+
+      CALL outputProbes(t)
 
       counterOutput = counterOutput + 1
       IF (counterOutput >= triggerOutput .OR. &
-          timeStep == tFinal .OR. timeStep == tInitial) THEN
+          t == tFinal .OR. t == tInitial) THEN
 
         !Resets output counter
         counterOutput=0
 
-        CALL mesh%printOutput()
+        CALL mesh%printOutput(t)
-        IF (ASSOCIATED(meshForMCC)) CALL meshForMCC%printColl()
+        IF (ASSOCIATED(meshForMCC)) CALL meshForMCC%printColl(t)
-        CALL mesh%printEM()
+        CALL mesh%printEM(t)
-        WRITE(*, "(5X,A21,I10,A1,I10)") "t/tFinal: ", timeStep, "/", tFinal
+        WRITE(*, "(5X,A21,I10,A1,I10)") "t/tFinal: ", t, "/", tFinal
         WRITE(*, "(5X,A21,I10)")        "Particles: ", nPartOld
-        IF (timeStep == 0) THEN
+        IF (t == 0) THEN
           WRITE(*, "(5X,A21,F8.1,A2)")   "     init time: ",    1.D3*tStep, "ms"
 
         ELSE
@@ -541,8 +572,8 @@ MODULE moduleSolver
 
         END IF
 
-        IF (nPartOld > 0) THEN
+        IF (nPartOldTotal > 0) THEN
-          WRITE(*, "(5X,A21,F8.1,A2)") "avg t/particle: ", 1.D9*tStep/DBLE(nPartOld), "ns"
+          WRITE(*, "(5X,A21,F8.1,A2)") "avg t/particle: ", 1.D9*tStep/DBLE(nPartOldTotal), "ns"
 
         END IF
         WRITE(*,*)
@@ -551,32 +582,34 @@ MODULE moduleSolver
 
       counterCPUTime = counterCPUTime + 1
       IF (counterCPUTime >= triggerCPUTime .OR. &
-          timeStep == tFinal .OR. timeStep == tInitial) THEN
+          t == tFinal .OR. t == tInitial) THEN
 
         !Reset CPU Time counter
         counterCPUTime = 0
 
-        CALL printTime(timeStep == 0)
+        CALL printTime(t, t == 0)
 
       END IF
 
       !Output average values
-      IF (useAverage .AND. timeStep == tFinal) THEN
+      IF (useAverage .AND. t == tFinal) THEN
         CALL mesh%printAverage()
 
       END IF
 
     END SUBROUTINE doOutput
 
-    SUBROUTINE doAverage()
+    SUBROUTINE doAverage(t)
       USE moduleAverage
       USE moduleMesh
       IMPLICIT NONE
 
+      INTEGER, INTENT(in):: t
       INTEGER:: tAverage, n
 
+
       IF (useAverage) THEN
-        tAverage = timeStep - tAverageStart
+        tAverage = t - tAverageStart
 
         IF (tAverage == 1) THEN
           !First iteration in which average scheme is used

src/modules/solver/pusher/modulePusher.f90

@@ -14,7 +14,7 @@ MODULE modulePusher
     END SUBROUTINE pushCartNeutral
 
     PURE SUBROUTINE pushCartElectrostatic(part, tauIn)
-      USE moduleSPecies
+      USE moduleSpecies
       USE moduleMesh
       IMPLICIT NONE