Merge branch 'feature/doubleMesh' into 'development'

Improvements to mesh and double mesh

See merge request JorgeGonz/fpakc!13

doc/user-manual/bibliography.bib

@@ -1,4 +1,4 @@
-% Encoding: UTF-8
+% Encoding: UTF-8
 
 @InProceedings{boris1970relativistic,
   author    = {Boris, Jay P},
@@ -7,38 +7,48 @@
   year      = {1970},
   pages     = {3--67},
 }
-
+
-@Misc{gfortranURL,
+@article{higginson2020corrected,
-  author       = {GNU Project},
+  title={A corrected method for Coulomb scattering in arbitrarily weighted particle-in-cell plasma simulations},
-  title        = {gfortran - the GNU Fortran compiler},
+  author={Higginson, Drew Pitney and Holod, Ihor and Link, Anthony},
-  howpublished = {\url{https://gcc.gnu.org/wiki/GFortran}},
+  journal={Journal of Computational Physics},
-}
+  volume={413},
-
+  pages={109450},
-@Misc{ifortURL,
+  year={2020},
-  author       = {Intel\textsuperscript{\textregistered}},
+  publisher={Elsevier}
-  title        = {Intel\textsuperscript{\textregistered} Fortran Compiler},
+}
-  howpublished = {\url{https://software.intel.com/content/www/us/en/develop/tools/oneapi/components/fortran-compiler.html}},
+
-}
+@Misc{gfortranURL,
-
+  author       = {GNU Project},
-@Misc{openblasURL,
+  title        = {gfortran - the GNU Fortran compiler},
-  title        = {OpenBLAS, an optimized BLAS library},
+  howpublished = {\url{https://gcc.gnu.org/wiki/GFortran}},
-  howpublished = {\url{https://www.openblas.net/}},
+}
-}
+
-
+@Misc{ifortURL,
-@Misc{jsonfortranURL,
+  author       = {Intel\textsuperscript{\textregistered}},
-  title        = {JSON-Fortran},
+  title        = {Intel\textsuperscript{\textregistered} Fortran Compiler},
-  howpublished = {\url{https://github.com/jacobwilliams/json-fortran}},
+  howpublished = {\url{https://software.intel.com/content/www/us/en/develop/tools/oneapi/components/fortran-compiler.html}},
-}
+}
-
+
-@Misc{jsonURL,
+@Misc{openblasURL,
-  title        = {JSON, JavaScript Object Notation},
+  title        = {OpenBLAS, an optimized BLAS library},
-  howpublished = {\url{https://www.json.org/json-en.html}},
+  howpublished = {\url{https://www.openblas.net/}},
-}
+}
-
+
-@Misc{gmshURL,
+@Misc{jsonfortranURL,
-  author       = {Christophe Geuzaine and Jean-François Remacle},
+  title        = {JSON-Fortran},
-  title        = {Gmsh},
+  howpublished = {\url{https://github.com/jacobwilliams/json-fortran}},
-  howpublished = {\url{https://gmsh.info/}},
+}
-}
+
-
+@Misc{jsonURL,
-@Comment{jabref-meta: databaseType:bibtex;}
+  title        = {JSON, JavaScript Object Notation},
+  howpublished = {\url{https://www.json.org/json-en.html}},
+}
+
+@Misc{gmshURL,
+  author       = {Christophe Geuzaine and Jean-François Remacle},
+  title        = {Gmsh},
+  howpublished = {\url{https://gmsh.info/}},
+}
+
+@Comment{jabref-meta: databaseType:bibtex;}

doc/user-manual/fpakc_UserManual.tex

@@ -21,6 +21,10 @@
     Author   = {Jorge Gonzalez}
   }
 }
+% Allows breaking of URL in bibliography.
+\setcounter{biburllcpenalty}{7000}
+\setcounter{biburlucpenalty}{8000}
+
 \makeglossaries
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -31,6 +35,8 @@
 \newacronym{cpu}{CPU}{Central Processing Unit}
 \newacronym{json}{JSON}{JavaScript Object Notation}
 \newacronym{io}{I/O}{input/output}
+\newacronym{mcc}{MCC}{Monte-Carlo Collisions}
+\newacronym{cs}{CS}{Coulomb Scattering}
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \newglossaryentry{openmp}{name={OpenMP},description={Shared-memory parallelization}}
@@ -60,7 +66,9 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   \chapter{Introduction}
     \section{About fpakc}
-      The \Gls{fpakc} is a simulation tool that models species in plasmas (ions, electrons and neutrals) following the trajectories of macro-particles as they move and interact between them and the boundaries of the domain.
+      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 very early steps of development and further improvements are expected very soon.
     
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -78,11 +86,11 @@
               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 to be ease to use.
-              Input files are required to be in a \textit{human} format, meaning that the different options can be easily understander without constant reference to the user guide.
+              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 a variety of scientist: from very established ones to newcomers to the field and 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 the code and code development and should always be followed, at least for the releases in the official repository.
+      These are foundation stones of \acrshort{fpakc} and its development and should always be followed, at least for the releases in the official repository.
     
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     \section{How to collaborate}
@@ -92,22 +100,32 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   \chapter{Operation Scheme}
     \section{The Particle Method}
-      \acrshort{fpakc} uses macro-particles to simulate the dynamics of different plasma species (mainly ions, electrons and neutrals).
+      \Gls{fpakc} uses macro-particles to simulate the dynamics of different plasma species (mainly ions, electrons and neutrals).
-      These macro-particles represent a large amount of real particles.
+      These macro-particles could represent a large amount of real particles.
       For now own, macro-particles will be referred as just particles by abusing of language.
-      In the evolution of these particles, external forces (as the electromagnetic field), interaction between particles (as collisions) and interaction with the boundaries of the domain are included.
+      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.
+              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 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.
+              Monte-Carlo collisions (elastic, ionization, charge-exchange$\ldots$) can be carried out in a specific mesh, to better adjust to the cell size required, similar to the mean-free path.
+              Although not yet implemented, Coulomb scattering will always be performed in the mesh used for scattering, which cell size should be in the order of the Debye length.
+        \item A new array containing all particles inside the numerical domain is obtained.
+        \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 is completed.
+      \end{enumerate}
     
-      At each time step, particles are first pushed accounting for possible acceleration by external forces.
+      \Gls{fpakc} has the capability to configure all the behavior of the simulation via the input file.
-      Then, the cell in which the particle ends up is located.
+      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}.
-      If a boundary is encountered, the interaction between the particle and the boundary is calculated.
-      Next, collisions for the particles in each cell are carried on.
-      This may include different collision processes for each particle.
-      Finally, the particles properties are scattered into the mesh nodes.
-      These properties are density, momentum and the stress tensor.
-      Non-dimensional units are used for this, but output files are converted into dimensional units.
-      If requested, the electromagnetic field is computed.
-    
-      More in depth explanation of the different steps are given in the following sections.
     
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     \section{Injection of new particles}
@@ -122,85 +140,93 @@
       Particles are pushed in the selected domain.
       Velocity and position are updated according to the old particle values and the external forces.
       All the push routines for the different geometries can be found in \lstinline|moduleSolver|.
+      The pushers included in \Gls{fpakc} are:
       
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+      \begin{itemize}
-      \subsection{3D Cartesian pusher}
+      \item 3D Cartesian pusher.
-        Moving particles in a simple 3D Cartesian space.
+            Moving particles in a simple 3D Cartesian space.
         
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+      \item 2D Cylindrical.
-      \subsection{2D Cylindrical}
+            When a 2D cylindrical geometry is used ($z$, $r$), a Boris solver\cite{boris1970relativistic} is used to move particles accounting for the effect of the symmetry axis.
-        When a 2D cylindrical geometry is used ($z$, $r$), a Boris solver\cite{boris1970relativistic} is used to move particles accounting for the effect of the symmetry axis.
+            This pusher removes the issue with particles going to infinite velocity when $r \rightarrow 0$ by pushing the particles in Cartesian space and then converting it to $r-z$ geometry.
-        This pusher removes the issue with particles going to infinite velocity when $r \rightarrow 0$ by pushing the particles in Cartesian space and then converting it to $r-z$ geometry.
+            Velocity in the $\theta$ direction is updated for collision processes, although the dynamic in the angular direction is assumed as symmetric.
-        Velocity in the $\theta$ direction is updated for collision processes, although the dynamic in the angular direction is assumed as symmetric.
         
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+      \item 2D Cartesian pusher.
-      \subsection{2D Cartesian pusher}
+            Moving particles in a simple 2D Cartesian space.
-        Moving particles in a simple 2D Cartesian space.
         
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+      \item 1D Radial pusher.
-      \subsection{1D Radial pusher}
+            Same implementation as 2D cylindrical pusher but direction $z$ is ignored.
-        Same implementation as 2D cylindrical pusher but direction $z$ is ignored.
 
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+      \item 1D Cartesian pusher.
-      \subsection{1D Cartesian pusher}
+            Moving particles in a simple 1D Cartesian space. Same implementation as in 2D Cartesian but $y$ direction is ignored.
-        Moving particles in a simple 1D Cartesian space. Same implementation as in 2D Cartesian but $y$ direction is ignored.
       
+      \end{itemize}
+    
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     \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 neighbor search, starting from the previous cell containing the particle.
       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 (reflected) or might be eliminated from it (absorbed).
+      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.
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     \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}.
-      This schemes detect when a particle change cells and modify its weight accordingly.
+      These schemes detect when a particle change cells and modify its weight accordingly.
       To avoid particles having a larger weight than the rest, particle can be split in multiple particles if weight become too large.
-                  
+      The use of a Variable Weighting Scheme is defined by the user in the input file.
+    
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    \section{Interaction between species}\label{ssec:collisions}
+      For each cell, interaction among the particles in it are carried out.
+      \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 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 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 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}:
+        \begin{itemize}
+          \item Elastic.
+                In this type of collision, particles exchange energy due to hard-sphere model.
+                Total energy is conserved.
+                Resulting velocity directions are chosen from Maxwellian distribution functions.
+                This interaction is useful for short-range collisions as neutral-neutral and charged-neutral elastic collisions.
+      
+          \item Charge Exchange.
+                When an ion interacts with a neutral particle, an electron is exchanged between the two particles with no exchange of energy.
+                This is called a resonant charge-exchange.
+        
+          \item Electron Impact Ionization.
+                When the relative energy between a neutral and an electron is above the ionization threshold, there is a probability that the neutral particle will become ionized.
+                This ionization emits and additional electron
+
+          \item Recombination.
+                When an electron and an ion interact, there is a possibility for them to be recombined into a neutral particle.
+                The photon emitted by this process is not modeled yet.
+        \end{itemize}
+      
+      \subsection{\acrlong{cs}}
+        Although not yet implement, a first approach will be soon implemented using Ref.~\cite{higginson2020corrected} as a guideline.
+        
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     \section{Reset of particle array}
       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 parallelisation and is very CPU consuming.
-    
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-    \section{Interaction between species}\label{ssec:collisions}
-      For each cell, interaction among the particles in it are carried out.
-      The type of interaction between the different particles is defined by the user.
-      In general, the maximum number of interaction in a cell is computed.
-      For each collision, a pair of particles is selected.
-      A loop over all possible collisions for the pair of particles is performed.
-      If a random number generated is above the probability of collision for the type divided by the maximum one, the collision take place.
-      
-      Collisions can change the velocity of the particles involved (elastic), create new particles (ionization-recombination) or change the type of particle (charge-exchange).
-      
-      Below are described the type of collision process implemented between species.
-      
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-      \subsection{Elastic collision}
-        In this type of collision, particles exchange energy due to hard-sphere model.
-        Total energy is conserved.
-        Resulting velocity directions are chosen from Maxwellian distribution functions.
-        This interaction is useful for short-range collisions as neutral-neutral and charged-neutral elastic collisions.
-      
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-      \subsection{Charge Exchange}
-        When an ion interacts with a neutral particle, an electron is exchanged between the two particles with no exchange of energy.
-        This is called a resonant charge-exchange.
-        
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-      \subsection{Electron Impact Ionization}
-        When the relative energy between a neutral and an electron is above the ionization threshold, there is a probability that the neutral particle will become ionized.
-        This ionization emits and additional electron
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-      \subsection{Recombination}
-        When an electron and an ion interact, there is a possibility for them to be recombined into a neutral particle.
-        The photon emitted by this process is not modeled yet.
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     \section{Scattering}
@@ -358,7 +384,10 @@ make
                                \end{itemize}
           \item \textbf{meshType}: Character.
                                    Format of mesh file.
-                                   Currently, only the value \textbf{gmsh} is accepted, which makes reference to \Gls{gmsh} v2.0 output format.
+                                   Accepted formats are:
+                                   \begin{itemize}
+                                     \item \textbf{gmsh2}: \Gls{gmsh} file format in version 2.0.
+                                   \end{itemize}
           \item \textbf{meshFile}: Character.
                                    Mesh filename.
                                    This file is searched in the path \textbf{output.path} and must contain the file extension.
@@ -617,9 +646,14 @@ make
         \begin{itemize}
           \item \textbf{folderCollisions}: Character.
                                            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}.
+                                         If this value is not present, the mesh defined in \textbf{geometry.meshFile} is used for \acrshort{mcc}.
+                                         The format of this mesh needs to be the same as the one defined in \textbf{geometry.meshType}.
           \item \textbf{collisions}: Object.
                                      Array.
-                                     Contains the different binary collisions.
+                                     Contains the different short range interactions (\acrshort{mcc}).
                                      Multiple collision types can be defined for each pair of species.
                                      Each object in the array is defined by:
                                      \begin{itemize}
@@ -662,21 +696,22 @@ make
         \end{itemize}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   \chapter{Example runs\label{ch:exampleRuns}}
-    \section{1D Cathode}
+    \section{1D Emissive Cathode (1D\_Cathode)}
       Emission from a 1D cathond in both, cartesian and radial coordinates.
       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 ilustrate hoy \acrshort{fpakc} can deal with different geometries by just modifiying some parameters in the input file.
       The same mesh file (\lstinline|mesh.msh|) is used for both cases.
     
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-    \section{ALPHIE Grid system}
+    \section{ALPHIE Grid system (ALPHIE\_Grid)}
       Two-dimensional axialsymmetry case to study the counterflow of electrons and Argon ions going trhough the ALPHIE grid system.
       A \lstinline|mesh.geo| file is provided to easily modify the parameters of the grid system and generate a new mesh with \Gls{gmsh}.
       
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-    \section{Flow around cylinder}
+    \section{Flow around cylinder (cylFlow)}
       Simple case of neutral Argon flow around a cylinder in a 2D axialsymmetry geometry.
-      Elastic collisions between argon particles are included as default.
+      Elastic collisions between argon particles are included.
+      Two cases are presented here: one in which the same mesh (meshSingle.msh) for scattering and collisions is used (input.json) and a second one (inputDualMesh.json) in which a mesh is used for scattering (mesh.msh) and a second one is used only for collisions (meshColl.msh).
     
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   \printglossaries

runs/1D_Cathode/inputCart.json

@@ -15,8 +15,8 @@
   },
   "geometry": {
     "type":     "1DCart",
-    "meshType": "gmsh",
+    "meshType": "gmsh2",
-      "meshFile": "mesh.msh"
+    "meshFile": "mesh.msh"
   },
   "species": [
     {"name": "Electron", "type": "charged", "mass": 9.109e-31, "charge":-1.0, "weight": 1.0e1},

runs/1D_Cathode/inputRadEmission.json

@@ -15,8 +15,8 @@
   },
   "geometry": {
     "type":     "1DRad",
-    "meshType": "gmsh",
+    "meshType": "gmsh2",
-      "meshFile": "mesh.msh"
+    "meshFile": "mesh.msh"
   },
   "species": [
     {"name": "Electron", "type": "charged", "mass": 9.109e-31, "charge":-1.0, "weight": 1.0e1},

runs/ALPHIE_Grid/inputBaseCase.json

@@ -0,0 +1,76 @@
+{
+  "output": {
+    "path": "./runs/ALPHIE_Grid/",
+    "triggerOutput": 500,
+    "cpuTime": false,
+    "numColl": false,
+    "EMField": true,
+	"folder": "base_case"
+  },
+  "geometry": {
+    "type":     "2DCyl",
+    "meshType": "gmsh2",
+    "meshFile": "mesh.msh"
+  },
+  "species": [
+    {"name": "Argon+",   "type": "charged", "mass": 6.633e-26, "charge": 1.0, "weight": 1.0e1},
+    {"name": "Electron", "type": "charged", "mass": 9.109e-31, "charge":-1.0, "weight": 1.0e2}
+  ],
+  "boundary": [ 
+    {"name": "Ionization Chanber", "physicalSurface": 1, "bTypes": [
+                                                           {"type": "transparent"},
+                                                           {"type": "transparent"}
+                                                          ]},
+    {"name": "Vacuum Chamber",     "physicalSurface": 2, "bTypes": [
+                                                           {"type": "transparent"},
+                                                           {"type": "transparent"}
+                                                          ]},
+    {"name": "Exterior",           "physicalSurface": 3, "bTypes": [
+                                                           {"type": "reflection"},
+                                                           {"type": "reflection"}
+                                                          ]},
+    {"name": "Grid Extraction",    "physicalSurface": 4, "bTypes": [
+                                                           {"type": "absorption"},
+                                                           {"type": "absorption"}
+                                                          ]},
+    {"name": "Grid Acceleration",  "physicalSurface": 5, "bTypes": [
+                                                           {"type": "absorption"},
+                                                           {"type": "absorption"}
+                                                          ]},
+    {"name": "Axis",               "physicalSurface": 6, "bTypes": [
+                                                           {"type": "axis"},
+                                                           {"type": "axis"}
+                                                          ]}
+  ],
+  "boundaryEM": [
+    {"name": "Extraction Grid",    "type": "dirichlet", "potential": -150.0, "physicalSurface": 4},
+    {"name": "Acceleration Grid",  "type": "dirichlet", "potential": -600.0, "physicalSurface": 5},
+    {"name": "Ionization Chamber", "type": "dirichlet", "potential":    0.0, "physicalSurface": 1}
+  ],
+  "inject": [
+    {"name": "Ionization Argon+",   "species": "Argon+",   "flow": 27.0e-6, "units": "A", "v":    322.0, "T": [ 500.0,  500.0,  500.0],
+                                    "velDist": ["Maxwellian", "Maxwellian", "Maxwellian"], "n": [ 1, 0, 0], "physicalSurface": 1},
+    {"name": "Ionization Electron", "species": "Electron", "flow": 27.0e-6, "units": "A", "v":  87000.0, "T": [ 500.0,  500.0,  500.0],
+                                    "velDist": ["Maxwellian", "Maxwellian", "Maxwellian"], "n": [ 1, 0, 0], "physicalSurface": 1},
+    {"name": "Cathode Electron",    "species": "Electron", "flow":  9.0e-5, "units": "A", "v":  87000.0, "T": [2500.0, 2500.0, 2500.0],
+                                    "velDist": ["Maxwellian", "Maxwellian", "Maxwellian"], "n": [-1, 0, 0], "physicalSurface": 2}
+  ],
+  "reference": {
+    "density": 1.0e16,
+    "mass": 9.109e-31,
+    "temperature": 2500.0,
+    "radius": 1.88e-10
+  },
+  "case": {
+    "tau":  [1.0e-9, 1.0e-11],
+    "time": 1.0e-6,
+    "pusher": ["2DCylCharged", "2DCylCharged"],
+    "WeightingScheme": "Volume",
+    "EMSolver": "Electrostatic"
+  },
+  "parallel": {
+    "OpenMP":{
+      "nThreads": 24
+    }
+  }
+}

runs/ALPHIE_Grid/inputIonization_0.10mA.json

@@ -0,0 +1,73 @@
+{
+  "output": {
+    "path": "./runs/ALPHIE_Grid/",
+    "triggerOutput": 500,
+    "cpuTime": false,
+    "numColl": false,
+    "EMField": true,
+	"folder": "ionization_0.10mA"
+  },
+  "geometry": {
+    "type":     "2DCyl",
+    "meshType": "gmsh2",
+    "meshFile": "mesh.msh"
+  },
+  "species": [
+    {"name": "Argon+",   "type": "charged", "mass": 6.633e-26, "charge": 1.0, "weight": 1.0e1},
+    {"name": "Electron", "type": "charged", "mass": 9.109e-31, "charge":-1.0, "weight": 1.0e2}
+  ],
+  "boundary": [ 
+    {"name": "Ionization Chanber", "physicalSurface": 1, "bTypes": [
+                                                           {"type": "transparent"},
+                                                           {"type": "ionization", "neutral": {"ion": "Argon+", "mass": 6.633e-26, "density": 1.0e17, "velocity": [323, 0, 0], "temperature": 300},
+                                                                                  "effectiveTime": 5.0e-6,"energyThreshold": 15.76, "crossSection": "./data/collisions/IO_e-Ar.dat"}
+                                                          ]},
+    {"name": "Vacuum Chamber",     "physicalSurface": 2, "bTypes": [
+                                                           {"type": "transparent"},
+                                                           {"type": "transparent"}
+                                                          ]},
+    {"name": "Exterior",           "physicalSurface": 3, "bTypes": [
+                                                           {"type": "reflection"},
+                                                           {"type": "reflection"}
+                                                          ]},
+    {"name": "Grid Extraction",    "physicalSurface": 4, "bTypes": [
+                                                           {"type": "absorption"},
+                                                           {"type": "absorption"}
+                                                          ]},
+    {"name": "Grid Acceleration",  "physicalSurface": 5, "bTypes": [
+                                                           {"type": "absorption"},
+                                                           {"type": "absorption"}
+                                                          ]},
+    {"name": "Axis",               "physicalSurface": 6, "bTypes": [
+                                                           {"type": "axis"},
+                                                           {"type": "axis"}
+                                                          ]}
+  ],
+  "boundaryEM": [
+    {"name": "Extraction Grid",    "type": "dirichlet", "potential": -150.0, "physicalSurface": 4},
+    {"name": "Acceleration Grid",  "type": "dirichlet", "potential": -600.0, "physicalSurface": 5},
+    {"name": "Ionization Chamber", "type": "dirichlet", "potential":    0.0, "physicalSurface": 1}
+  ],
+  "inject": [
+    {"name": "Cathode Electron",    "species": "Electron", "flow":  1.0e-4, "units": "A", "v":  87000.0, "T": [2500.0, 2500.0, 2500.0],
+                                    "velDist": ["Maxwellian", "Maxwellian", "Maxwellian"], "n": [-1, 0, 0], "physicalSurface": 2}
+  ],
+  "reference": {
+    "density": 1.0e16,
+    "mass": 9.109e-31,
+    "temperature": 2500.0,
+    "radius": 1.88e-10
+  },
+  "case": {
+    "tau":  [1.0e-9, 1.0e-11],
+    "time": 1.0e-6,
+    "pusher": ["2DCylCharged", "2DCylCharged"],
+    "WeightingScheme": "Volume",
+    "EMSolver": "Electrostatic"
+  },
+  "parallel": {
+    "OpenMP":{
+      "nThreads": 24
+    }
+  }
+}

runs/Argon_Expansion/CX_case.json

@@ -3,12 +3,12 @@
     "path": "./runs/Argon_Expansion/",
 	"triggerOutput": 10,
 	"cpuTime": false,
-	"numColl": false,
+	"numColl": true,
 	"folder": "CX_case"
   },
   "geometry": {
     "type":     "2DCyl",
-    "meshType": "gmsh",
+    "meshType": "gmsh2",
     "meshFile": "mesh.msh"
   },
   "species": [

runs/Argon_Expansion/elastic_case.json

@@ -8,7 +8,7 @@
   },
   "geometry": {
     "type":     "2DCyl",
-    "meshType": "gmsh",
+    "meshType": "gmsh2",
     "meshFile": "mesh.msh"
   },
   "species": [

runs/Argon_Expansion/nocoll_case.json

@@ -0,0 +1,55 @@
+{
+  "output": {
+    "path": "./runs/Argon_Expansion/",
+	"triggerOutput": 10,
+	"cpuTime": false,
+	"numColl": false,
+	"folder": "Nocoll_case"
+  },
+  "geometry": {
+    "type":     "2DCyl",
+    "meshType": "gmsh2",
+    "meshFile": "mesh.msh"
+  },
+  "species": [
+    {"name": "Argon",  "type": "neutral", "mass": 6.633e-26, "weight": 1.0e8,                "ion": "Argon+"},
+    {"name": "Argon+", "type": "charged", "mass": 6.633e-26, "weight": 1.0e8, "charge": 1.0, "neutral": "Argon"}
+  ],
+  "boundary": [ 
+    {"name": "Injection", "physicalSurface": 1, "bTypes": [
+                                                  {"type": "transparent"},
+                                                  {"type": "transparent"}
+                                                 ]},
+    {"name": "Exterior",  "physicalSurface": 2, "bTypes": [
+                                                  {"type": "transparent"},
+                                                  {"type": "transparent"}
+                                                 ]},
+    {"name": "Axis",      "physicalSurface": 3, "bTypes": [
+                                                  {"type": "axis"},
+                                                  {"type": "axis"}
+                                                 ]}
+  ],
+  "inject": [
+    {"name": "Exhausts Ar", "species": "Argon", "flow": 0.7, "units": "sccm",     "v": 300.0, "T": [300.0, 300.0, 300.0],
+                            "velDist": ["Maxwellian", "Maxwellian", "Maxwellian"], "n": [1, 0, 0], "physicalSurface": 1},
+    {"name": "Exhausts Ar+", "species": "Argon+", "flow": 0.3, "units": "sccm",     "v": 40000.0, "T": [300.0, 300.0, 300.0],
+                             "velDist": ["Maxwellian", "Maxwellian", "Maxwellian"], "n": [1, 0, 0], "physicalSurface": 1}
+  ],
+  "reference": {
+    "density": 1.0e19,
+    "mass": 6.633e-26,
+    "temperature": 300.0,
+    "radius": 1.88e-10
+  },
+  "case": {
+    "tau":  [1.0e-6, 1.0e-6],
+    "time": 4.0e-3,
+    "pusher": ["2DCylNeutral", "2DCylNeutral"],
+    "WeightingScheme": "Volume"
+  },
+  "parallel": {
+    "OpenMP":{
+      "nThreads": 24
+    }
+  }
+}

runs/cylFlow/input.json

@@ -3,12 +3,12 @@
     "path": "./runs/cylFlow/",
 	"triggerOutput": 10,
 	"cpuTime": true,
-	"numColl": false
+	"numColl": true
   },
   "geometry": {
     "type":     "2DCyl",
-    "meshType": "gmsh",
+    "meshType": "gmsh2",
-    "meshFile": "mesh.msh"
+    "meshFile": "meshSingle.msh"
   },
   "species": [
     {"name": "Argon",  "type": "neutral", "mass": 6.633e-26, "weight": 5.0e8}
@@ -57,7 +57,7 @@
   },
   "parallel": {
     "OpenMP":{
-      "nThreads": 8
+      "nThreads": 24
     }
   }
 }

runs/cylFlow/inputDualMesh.json

@@ -0,0 +1,64 @@
+{
+  "output": {
+    "path": "./runs/cylFlow/",
+	"triggerOutput": 10,
+	"cpuTime": true,
+	"numColl": true
+  },
+  "geometry": {
+    "type":     "2DCyl",
+    "meshType": "gmsh2",
+    "meshFile": "mesh.msh"
+  },
+  "species": [
+    {"name": "Argon",  "type": "neutral", "mass": 6.633e-26, "weight": 5.0e8}
+  ],
+  "boundary": [ 
+    {"name": "Injection",      "physicalSurface": 1, "bTypes": [
+                                                       {"type": "transparent"}
+                                                      ]},
+    {"name": "Chamber Walls",  "physicalSurface": 2, "bTypes": [
+                                                       {"type": "reflection"}
+                                                      ]},
+    {"name": "Exterior",       "physicalSurface": 3, "bTypes": [
+                                                       {"type": "transparent"}
+                                                      ]},
+    {"name": "Cylinder Walls", "physicalSurface": 4, "bTypes": [
+                                                       {"type": "reflection"}
+                                                      ]},
+    {"name": "Axis",           "physicalSurface": 5, "bTypes": [
+                                                       {"type": "axis"}
+                                                      ]}
+  ],
+  "inject": [
+    {"name": "Nozzle", "species": "Argon", "flow": 10.0, "units": "sccm",     "v": 300.0, "T": [300.0, 300.0, 300.0],
+                       "velDist": ["Maxwellian", "Maxwellian", "Maxwellian"], "n": [1, 0, 0], "physicalSurface": 1}
+  ],
+  "reference": {
+    "density": 1.0e20,
+    "mass": 6.633e-26,
+    "temperature": 300.0,
+    "radius": 1.88e-10
+  },
+  "case": {
+    "tau":  [5.0e-7],
+    "time":  1.0e-3,
+    "pusher": ["2DCylNeutral"],
+    "WeightingScheme": "Volume"
+  },
+  "interactions": {
+    "folderCollisions": "./data/collisions/",
+    "meshCollisions": "meshColl.msh",
+    "collisions": [
+      {"species_i": "Argon", "species_j": "Argon", 
+       "cTypes": [
+         {"type": "elastic", "crossSection": "EL_Ar-Ar.dat"}
+       ]}
+    ]
+  },
+  "parallel": {
+    "OpenMP":{
+      "nThreads": 24
+    }
+  }
+}

runs/cylFlow/mesh.geo

@@ -62,10 +62,13 @@ Physical Surface(4) = {4};
 Physical Surface(5) = {5};
 
 Transfinite Line {12, 2, 4, 6} = cyl_h/Lcell + 1 Using Progression 1;
-Transfinite Line {1, 13, 10} = cyl_s/Lcell + 1 Using Progression 1;
+Transfinite Line {10} = cyl_s/Lcell + 1 Using Progression 1;
-Transfinite Line {11, 16, 15, 7} = (dom_h - cyl_h)/Lcell + 1 Using Progression 1;
+Transfinite Line {1, 13} = cyl_s/Lcell + 1 Using Progression 0.95;
+Transfinite Line {11, 7} = (dom_h - cyl_h)/Lcell + 1 Using Progression 1;
+Transfinite Line {16,-15} = (dom_h - cyl_h)/Lcell + 1 Using Progression 0.95;
 Transfinite Line {3, 9} = cyl_l/Lcell + 1 Using Progression 1;
-Transfinite Line {5, 14, 8} = (dom_l - cyl_e)/Lcell + 1 Using Progression 1;
+Transfinite Line {8} = (dom_l - cyl_e)/Lcell + 1 Using Progression 1;
+Transfinite Line {-5, -14} = (dom_l - cyl_e)/Lcell + 1 Using Progression 0.95;
 
 Transfinite Surface{1};
 Recombine Surface {1};

runs/cylFlow/meshColl.geo

@@ -0,0 +1,73 @@
+cl__1 = 1;
+cyl_h = 0.005;
+cyl_l = 0.02;
+cyl_s = 0.03;
+cyl_e = cyl_s + cyl_l;
+dom_h = 0.03;
+dom_l = 0.07;
+
+Lcell = 0.001;
+
+Point(1)  = {0,     0,     0, cl__1};
+Point(2)  = {cyl_s, 0,     0, cl__1};
+Point(3)  = {cyl_s, cyl_h, 0, cl__1};
+Point(4)  = {cyl_e, cyl_h, 0, cl__1};
+Point(5)  = {cyl_e, 0,     0, cl__1};
+Point(6)  = {dom_l, 0,     0, cl__1};
+Point(7)  = {dom_l, cyl_h, 0, cl__1};
+Point(8)  = {dom_l, dom_h, 0, cl__1};
+Point(9)  = {cyl_e, dom_h, 0, cl__1};
+Point(10) = {cyl_s, dom_h, 0, cl__1};
+Point(11) = {0,     dom_h, 0, cl__1};
+Point(12) = {0,     cyl_h, 0, cl__1};
+
+Line(1)  = {1,  2};
+Line(2)  = {2,  3};
+Line(3)  = {3,  4};
+Line(4)  = {4,  5};
+Line(5)  = {5,  6};
+Line(6)  = {6,  7};
+Line(7)  = {7,  8};
+Line(8)  = {8,  9};
+Line(9)  = {9,  10};
+Line(10) = {10, 11};
+Line(11) = {11, 12};
+Line(12) = {12, 1};
+Line(13) = {12, 3};
+Line(14) = {4,  7};
+Line(15) = {4,  9};
+Line(16) = {10, 3};
+
+Line Loop(1) = {1, 2, -13, 12};
+Plane Surface(1) = {1};
+Line Loop(2) = {13, -16, 10, 11};
+Plane Surface(2) = {2};
+Line Loop(3) = {3, 15, 9, 16};
+Plane Surface(3) = {3};
+Line Loop(4) = {5, 6, -14, 4};
+Plane Surface(4) = {4};
+Line Loop(5) = {14, 7, 8, -15};
+Plane Surface(5) = {5};
+
+Physical Surface(1) = {1};
+Physical Surface(2) = {2};
+Physical Surface(3) = {3};
+Physical Surface(4) = {4};
+Physical Surface(5) = {5};
+
+Transfinite Line {12, 2, 4, 6} = cyl_h/Lcell + 1 Using Progression 1;
+Transfinite Line {1, 13, 10} = cyl_s/Lcell + 1 Using Progression 1;
+Transfinite Line {11, 16, 15, 7} = (dom_h - cyl_h)/Lcell + 1 Using Progression 1;
+Transfinite Line {3, 9} = cyl_l/Lcell + 1 Using Progression 1;
+Transfinite Line {5, 14, 8} = (dom_l - cyl_e)/Lcell + 1 Using Progression 1;
+
+Transfinite Surface{1};
+Recombine Surface {1};
+Transfinite Surface{2};
+Recombine Surface {2};
+Transfinite Surface{3};
+Recombine Surface {3};
+Transfinite Surface{4};
+Recombine Surface {4};
+Transfinite Surface{5};
+Recombine Surface {5};

runs/cylFlow/meshSingle.geo

@@ -0,0 +1,79 @@
+cl__1 = 1;
+cyl_h = 0.005;
+cyl_l = 0.02;
+cyl_s = 0.03;
+cyl_e = cyl_s + cyl_l;
+dom_h = 0.03;
+dom_l = 0.07;
+
+Lcell = 0.001;
+
+Point(1)  = {0,     0,     0, cl__1};
+Point(2)  = {cyl_s, 0,     0, cl__1};
+Point(3)  = {cyl_s, cyl_h, 0, cl__1};
+Point(4)  = {cyl_e, cyl_h, 0, cl__1};
+Point(5)  = {cyl_e, 0,     0, cl__1};
+Point(6)  = {dom_l, 0,     0, cl__1};
+Point(7)  = {dom_l, cyl_h, 0, cl__1};
+Point(8)  = {dom_l, dom_h, 0, cl__1};
+Point(9)  = {cyl_e, dom_h, 0, cl__1};
+Point(10) = {cyl_s, dom_h, 0, cl__1};
+Point(11) = {0,     dom_h, 0, cl__1};
+Point(12) = {0,     cyl_h, 0, cl__1};
+
+Line(1)  = {1,  2};
+Line(2)  = {2,  3};
+Line(3)  = {3,  4};
+Line(4)  = {4,  5};
+Line(5)  = {5,  6};
+Line(6)  = {6,  7};
+Line(7)  = {7,  8};
+Line(8)  = {8,  9};
+Line(9)  = {9,  10};
+Line(10) = {10, 11};
+Line(11) = {11, 12};
+Line(12) = {12, 1};
+Line(13) = {12, 3};
+Line(14) = {4,  7};
+Line(15) = {4,  9};
+Line(16) = {10, 3};
+
+Line Loop(1) = {1, 2, -13, 12};
+Plane Surface(1) = {1};
+Line Loop(2) = {13, -16, 10, 11};
+Plane Surface(2) = {2};
+Line Loop(3) = {3, 15, 9, 16};
+Plane Surface(3) = {3};
+Line Loop(4) = {5, 6, -14, 4};
+Plane Surface(4) = {4};
+Line Loop(5) = {14, 7, 8, -15};
+Plane Surface(5) = {5};
+
+Physical Line(1) = {12, 11};
+Physical Line(2) = {10, 9, 8};
+Physical Line(3) = {7, 6};
+Physical Line(4) = {2, 3, 4};
+Physical Line(5) = {1, 5};
+
+Physical Surface(1) = {1};
+Physical Surface(2) = {2};
+Physical Surface(3) = {3};
+Physical Surface(4) = {4};
+Physical Surface(5) = {5};
+
+Transfinite Line {12, 2, 4, 6} = cyl_h/Lcell + 1 Using Progression 1;
+Transfinite Line {1, 13, 10} = cyl_s/Lcell + 1 Using Progression 1;
+Transfinite Line {11, 16, 15, 7} = (dom_h - cyl_h)/Lcell + 1 Using Progression 1;
+Transfinite Line {3, 9} = cyl_l/Lcell + 1 Using Progression 1;
+Transfinite Line {5, 14, 8} = (dom_l - cyl_e)/Lcell + 1 Using Progression 1;
+
+Transfinite Surface{1};
+Recombine Surface {1};
+Transfinite Surface{2};
+Recombine Surface {2};
+Transfinite Surface{3};
+Recombine Surface {3};
+Transfinite Surface{4};
+Recombine Surface {4};
+Transfinite Surface{5};
+Recombine Surface {5};

src/fpakc.f90

@@ -4,6 +4,7 @@ PROGRAM fpakc
   USE moduleErrors
   USE moduleInject
   USE moduleSolver
+  USE moduleMesh
   USE moduleCompTime
   USE moduleCaseParam
   USE OMP_LIB
@@ -67,10 +68,19 @@ PROGRAM fpakc
     tColl = omp_get_wtime()
     !$OMP END SINGLE
 
-    CALL doCollisions()
+    IF (ASSOCIATED(meshForMCC)) CALL meshForMCC%doCollisions()
 
     !$OMP SINGLE
     tColl = omp_get_wtime() - tColl
+    
+    !Coulomb scattering
+    tCoul = omp_get_wTime()
+    !$OMP END SINGLE
+
+    IF (ASSOCIATED(mesh%doCoulomb)) CALL mesh%doCoulomb()
+
+    !$OMP SINGLE
+    tCoul = omp_get_wTime() - tCoul
 
     !Reset particles
     tReset = omp_get_wtime()

src/makefile

@@ -3,12 +3,13 @@ OBJECTS = $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleMeshBoundary.o $(OBJDIR)/module
 			  	$(OBJDIR)/moduleErrors.o $(OBJDIR)/moduleList.o $(OBJDIR)/moduleOutput.o \
 			  	$(OBJDIR)/moduleBoundary.o $(OBJDIR)/moduleCaseParam.o $(OBJDIR)/moduleRefParam.o \
 			  	$(OBJDIR)/moduleCollisions.o $(OBJDIR)/moduleTable.o $(OBJDIR)/moduleParallel.o \
-			  	$(OBJDIR)/moduleEM.o $(OBJDIR)/moduleRandom.o  $(OBJDIR)/moduleMath.o \
+			  	$(OBJDIR)/moduleEM.o $(OBJDIR)/moduleRandom.o $(OBJDIR)/moduleMath.o \
-			  	$(OBJDIR)/moduleMesh3DCart.o $(OBJDIR)/moduleMesh3DCartRead.o	\
+			  	$(OBJDIR)/moduleMeshInputGmsh2.o $(OBJDIR)/moduleMeshOutputGmsh2.o \
-			  	$(OBJDIR)/moduleMesh2DCyl.o  $(OBJDIR)/moduleMesh2DCylRead.o \
+			  	$(OBJDIR)/moduleMesh3DCart.o \
-			  	$(OBJDIR)/moduleMesh2DCart.o $(OBJDIR)/moduleMesh2DCartRead.o	\
+			  	$(OBJDIR)/moduleMesh2DCyl.o  \
-          $(OBJDIR)/moduleMesh1DCart.o $(OBJDIR)/moduleMesh1DCartRead.o \
+			  	$(OBJDIR)/moduleMesh2DCart.o \
-          $(OBJDIR)/moduleMesh1DRad.o  $(OBJDIR)/moduleMesh1DRadRead.o 
+          $(OBJDIR)/moduleMesh1DRad.o  \
+          $(OBJDIR)/moduleMesh1DCart.o
 
 all: $(OUTPUT)
 

src/modules/mesh/1DCart/makefile

@@ -1,8 +1,5 @@
-all: moduleMesh1DCart.o moduleMesh1DCartRead.o
+all: moduleMesh1DCart.o
 
 moduleMesh1DCart.o: moduleMesh1DCart.f90
 	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
 
-moduleMesh1DCartRead.o: moduleMesh1DCart.o moduleMesh1DCartRead.f90
-	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
-

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

@@ -71,7 +71,7 @@ MODULE moduleMesh1DCart
     !Connectivity to nodes
     CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL()
     !Connectivity to adjacent elements
-    CLASS(*), POINTER:: e1 => NULL(), e2 => NULL()
+    CLASS(meshElement), POINTER:: e1 => NULL(), e2 => NULL()
     REAL(8):: arNodes(1:2)
     CONTAINS
       PROCEDURE, PASS::   init        => initVol1DCartSegm
@@ -198,18 +198,19 @@ MODULE moduleMesh1DCart
     !VOLUME FUNCTIONS
     !SEGMENT FUNCTIONS
     !Init segment element
-    SUBROUTINE initVol1DCartSegm(self, n, p)
+    SUBROUTINE initVol1DCartSegm(self, n, p, nodes)
       USE moduleRefParam
       IMPLICIT NONE
 
       CLASS(meshVol1DCartSegm), INTENT(out):: self
       INTEGER, INTENT(in):: n
       INTEGER, INTENT(in):: p(:)
+      TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:)
       REAL(8), DIMENSION(1:3):: r1, r2
 
       self%n = n
-      self%n1 => mesh%nodes(p(1))%obj
+      self%n1 => nodes(p(1))%obj
-      self%n2 => mesh%nodes(p(2))%obj
+      self%n2 => nodes(p(2))%obj
       !Get element coordinates
       r1 = self%n1%getCoordinates()
       r2 = self%n2%getCoordinates()
@@ -308,27 +309,27 @@ MODULE moduleMesh1DCart
     END SUBROUTINE partialDerSegm
 
     !Computes local stiffness matrix
-    FUNCTION elemKSegm(self) RESULT(ke)
+    PURE FUNCTION elemKSegm(self) RESULT(localK)
       IMPLICIT NONE
 
       CLASS(meshVol1DCartSegm), INTENT(in):: self
-      REAL(8):: ke(1:2,1:2)
+      REAL(8), ALLOCATABLE:: localK(:,:)
       REAL(8):: Xii(1:3)
       REAL(8):: dPsi(1:1, 1:2)
       REAL(8):: invJ(1), detJ
       INTEGER:: l
 
-      ke      = 0.D0
+      ALLOCATE(localK(1:2,1:2))
+      localK      = 0.D0
       Xii     = 0.D0
-
       DO l = 1, 3
         xii(1) = corSeg(l)
         dPsi    = self%dPsi(Xii)
         detJ    = self%detJac(Xii, dPsi)
         invJ    = self%invJac(Xii, dPsi)
-        ke      = ke + MATMUL(RESHAPE(MATMUL(invJ,dPsi), (/ 2, 1/)), &
+        localK  = localK + MATMUL(RESHAPE(MATMUL(invJ,dPsi), (/ 2, 1/)), &
-                              RESHAPE(MATMUL(invJ,dPsi), (/ 1, 2/)))* &
+                                  RESHAPE(MATMUL(invJ,dPsi), (/ 1, 2/)))* &
-                       wSeg(l)/detJ
+                           wSeg(l)/detJ
 
       END DO
 
@@ -395,12 +396,12 @@ MODULE moduleMesh1DCart
       w_p = self%weight(part%xi)
       tensorS = outerProduct(part%v, part%v)
 
-      vertex => self%n1%output(part%sp)
+      vertex => self%n1%output(part%species%n)
       vertex%den          = vertex%den          + part%weight*w_p(1)
       vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(1)*part%v(:)
       vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(1)*tensorS
 
-      vertex => self%n2%output(part%sp)
+      vertex => self%n2%output(part%species%n)
       vertex%den          = vertex%den          + part%weight*w_p(2)
       vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(2)*part%v(:)
       vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(2)*tensorS
@@ -459,7 +460,7 @@ MODULE moduleMesh1DCart
 
       CLASS(meshVol1DCartSegm), INTENT(in):: self
       REAL(8), INTENT(in):: xi(1:3)
-      CLASS(*), POINTER, INTENT(out):: nextElement
+      CLASS(meshElement), POINTER, INTENT(out):: nextElement
       
       NULLIFY(nextElement)
       IF     (xi(1) < -1.D0) THEN
@@ -522,6 +523,121 @@ MODULE moduleMesh1DCart
 
     END FUNCTION invJ1DCart
 
+    SUBROUTINE connectMesh1DCart(self)
+      IMPLICIT NONE
+
+      CLASS(meshGeneric), INTENT(inout):: self
+      INTEGER:: e, et
+
+      DO e = 1, self%numVols
+        !Connect Vol-Vol
+        DO et = 1, self%numVols
+          IF (e /= et) THEN
+            CALL connectVolVol(self%vols(e)%obj, self%vols(et)%obj)
+
+          END IF
+
+        END DO
+
+        SELECT TYPE(self)
+        TYPE IS(meshParticles)
+          !Connect Vol-Edge
+          DO et = 1, self%numEdges
+            CALL connectVolEdge(self%vols(e)%obj, self%edges(et)%obj)
+
+          END DO
+
+      END SELECT
+
+      END DO
+
+    END SUBROUTINE connectMesh1DCart
+
+    SUBROUTINE connectVolVol(elemA, elemB)
+      IMPLICIT NONE
+
+      CLASS(meshVol), INTENT(inout):: elemA
+      CLASS(meshVol), INTENT(inout):: elemB
+
+        SELECT TYPE(elemA)
+        TYPE IS(meshVol1DCartSegm)
+          SELECT TYPE(elemB)
+          TYPE IS(meshVol1DCartSegm)
+            CALL connectSegmSegm(elemA, elemB)
+
+          END SELECT
+
+        END SELECT
+
+    END SUBROUTINE connectVolVol
+
+    SUBROUTINE connectSegmSegm(elemA, elemB)
+      IMPLICIT NONE
+
+      CLASS(meshVol1DCartSegm), INTENT(inout), TARGET:: elemA
+      CLASS(meshVol1DCartSegm), INTENT(inout), TARGET:: elemB
+
+      IF (.NOT. ASSOCIATED(elemA%e1) .AND. &
+          elemA%n2%n == elemB%n1%n) THEN
+        elemA%e1 => elemB
+        elemB%e2 => elemA
+        
+      END IF
+
+      IF (.NOT. ASSOCIATED(elemA%e2) .AND. &
+          elemA%n1%n == elemB%n2%n) THEN
+
+        elemA%e2 => elemB
+        elemB%e1 => elemA
+
+      END IF
+
+    END SUBROUTINE connectSegmSegm
+
+    SUBROUTINE connectVolEdge(elemA, elemB)
+      IMPLICIT NONE
+
+      CLASS(meshVol),  INTENT(inout):: elemA
+      CLASS(meshEdge), INTENT(inout):: elemB
+
+      SELECT TYPE(elemA)
+      TYPE IS (meshVol1DCartSegm)
+        SELECT TYPE(elemB)
+        CLASS IS(meshEdge1DCart)
+          CALL connectSegmEdge(elemA, elemB)
+
+        END SELECT
+
+      END SELECT
+
+    END SUBROUTINE connectVolEdge
+
+    SUBROUTINE connectSegmEdge(elemA, elemB)
+      IMPLICIT NONE
+
+      CLASS(meshVol1DCartSegm), INTENT(inout), TARGET:: elemA
+      CLASS(meshEdge1DCart),    INTENT(inout), TARGET:: elemB
+
+      IF (.NOT. ASSOCIATED(elemA%e1) .AND. &
+          elemA%n2%n == elemB%n1%n) THEN
+
+        elemA%e1 => elemB
+        elemB%e2 => elemA
+ 
+        !Revers the normal to point inside the domain
+        elemB%normal = - elemB%normal
+
+      END IF
+
+      IF (.NOT. ASSOCIATED(elemA%e2) .AND. &
+          elemA%n1%n == elemB%n1%n) THEN
+
+        elemA%e2 => elemB
+        elemB%e1 => elemA
+
+      END IF
+
+    END SUBROUTINE connectSegmEdge
 
 END MODULE moduleMesh1DCart
 

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

@@ -1,264 +0,0 @@
-MODULE moduleMesh1DCartRead
-  USE moduleMesh
-  USE moduleMesh1DCart
-  
-  !TODO: make this abstract to allow different mesh formats
-  TYPE, EXTENDS(meshGeneric):: mesh1DCartGeneric
-    CONTAINS
-      PROCEDURE, PASS:: init     => init1DCartMesh
-      PROCEDURE, PASS:: readMesh => readMesh1DCart
-
-  END TYPE
-
-  INTERFACE connected
-    MODULE PROCEDURE connectedVolVol, connectedVolEdge
-
-  END INTERFACE connected
-
-  CONTAINS
-    !Init 1D mesh
-    SUBROUTINE init1DCartMesh(self, meshFormat)
-      USE moduleMesh
-      USE moduleErrors
-      IMPLICIT NONE
-
-      CLASS(mesh1DCartGeneric), INTENT(out):: self
-      CHARACTER(:), ALLOCATABLE, INTENT(in):: meshFormat
-
-      SELECT CASE(meshFormat)
-      CASE ("gmsh")
-        self%printOutput => printOutputGmsh
-        self%printColl   => printCollGmsh
-        self%printEM     => printEMGmsh
-
-      CASE DEFAULT
-        CALL criticalError("Mesh type " // meshFormat // " not supported.", "init1DCart")
-        
-      END SELECT
-
-    END SUBROUTINE init1DCartMesh
-
-    !Reads 1D mesh
-    SUBROUTINE readMesh1DCart(self, filename)
-      USE  moduleBoundary
-      IMPLICIT NONE
-
-      CLASS(mesh1DCartGeneric), INTENT(inout):: self
-      CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
-      REAL(8):: x
-      INTEGER:: p(1:2)
-      INTEGER:: e, et, n, eTemp, elemType, bt
-      INTEGER:: totalNumElem
-      INTEGER:: boundaryType
-
-      !Open file mesh
-      OPEN(10, FILE=TRIM(filename))
-      !Skip header
-      READ(10, *)
-      READ(10, *)
-      READ(10, *)
-      READ(10, *)
-      !Read number of nodes
-      READ(10, *) self%numNodes
-      !Allocate required matrices and vectors
-      ALLOCATE(self%nodes(1:self%numNodes))
-      ALLOCATE(self%K(1:self%numNodes, 1:self%numNodes))
-      ALLOCATE(self%IPIV(1:self%numNodes, 1:self%numNodes))
-      self%K = 0.D0
-      self%IPIV = 0
-      !Read nodes coordinates. Only relevant for x
-      DO e = 1, self%numNodes
-        READ(10, *) n, x
-        ALLOCATE(meshNode1DCart:: self%nodes(n)%obj)
-        CALL self%nodes(n)%obj%init(n, (/ x, 0.D0, 0.D0 /))
-
-      END DO
-      !Skips comments
-      READ(10, *)
-      READ(10, *)
-      !Reads the total number of elements (edges+vol)
-      READ(10, *) totalNumElem
-      self%numEdges = 0
-      DO e = 1, totalNumElem
-        READ(10, *) eTemp, elemType
-        IF (elemType == 15) THEN !15 is physical node in GMSH
-          self%numEdges = e
-
-        END IF
-
-      END DO
-
-      !Substract the number of edges to the total number of elements
-      !to obtain the number of volume elements
-      self%numVols = totalNumelem - self%numEdges
-      !Allocates arrays
-      ALLOCATE(self%edges(1:self%numEdges))
-      ALLOCATE(self%vols(1:self%numVols))
-
-      !Go back to the beginning of reading elements
-      DO e = 1, totalNumelem
-        BACKSPACE(10)
-
-      END DO
-
-      !Reads edges
-      DO e = 1, self%numEdges
-        READ(10, *) n, elemType, eTemp, boundaryType, eTemp, p(1)
-        !Associate boundary condition 
-        bt = getBoundaryId(boundaryType)
-
-        ALLOCATE(meshEdge1DCart:: self%edges(e)%obj)
-
-        CALL self%edges(e)%obj%init(n, p(1:1), bt, boundaryType)
-
-      END DO
-
-      !Read and initialize volumes
-      DO e = 1, self%numVols
-        READ(10, *) n, elemType, eTemp, eTemp, eTemp, p(1:2)
-        ALLOCATE(meshVol1DCartSegm:: self%vols(e)%obj)
-        CALL self%vols(e)%obj%init(n - self%numEdges, p(1:2))
-
-      END DO
-
-      CLOSE(10)
-
-      !Build connectivity between elements
-      DO e = 1, self%numVols
-        !Connectivity between volumes
-        DO et = 1, self%numVols
-          IF (e /= et) THEN
-            CALL connected(self%vols(e)%obj, self%vols(et)%obj)
-
-          END IF
-
-        END DO
-
-        !Connectivity betwen vols and edges
-        DO et = 1, self%numEdges
-          CALL connected(self%vols(e)%obj, self%edges(et)%obj)
-
-        END DO
-
-        !Constructs the global K matrix
-        CALL constructGlobalK(self%K, self%vols(e)%obj)
-
-      END DO
-
-    END SUBROUTINE readMesh1DCart
-
-    SUBROUTINE connectedVolVol(elemA, elemB)
-      IMPLICIT NONE
-
-      CLASS(meshVol), INTENT(inout):: elemA
-      CLASS(meshVol), INTENT(inout):: elemB
-
-        SELECT TYPE(elemA)
-        TYPE IS(meshVol1DCartSegm)
-          SELECT TYPE(elemB)
-          TYPE IS(meshVol1DCartSegm)
-            CALL connectedSegmSegm(elemA, elemB)
-
-          END SELECT
-
-        END SELECT
-
-    END SUBROUTINE connectedVolVol
-
-    SUBROUTINE connectedSegmSegm(elemA, elemB)
-      IMPLICIT NONE
-
-      CLASS(meshVol1DCartSegm), INTENT(inout), TARGET:: elemA
-      CLASS(meshVol1DCartSegm), INTENT(inout), TARGET:: elemB
-
-      IF (.NOT. ASSOCIATED(elemA%e1) .AND. &
-          elemA%n2%n == elemB%n1%n) THEN
-        elemA%e1 => elemB
-        elemB%e2 => elemA
-        
-      END IF
-
-      IF (.NOT. ASSOCIATED(elemA%e2) .AND. &
-          elemA%n1%n == elemB%n2%n) THEN
-
-        elemA%e2 => elemB
-        elemB%e1 => elemA
-
-      END IF
-
-    END SUBROUTINE connectedSegmSegm
-
-    SUBROUTINE connectedVolEdge(elemA, elemB)
-      IMPLICIT NONE
-
-      CLASS(meshVol),  INTENT(inout):: elemA
-      CLASS(meshEdge), INTENT(inout):: elemB
-
-      SELECT TYPE(elemA)
-      TYPE IS (meshVol1DCartSegm)
-        SELECT TYPE(elemB)
-        CLASS IS(meshEdge1DCart)
-          CALL connectedSegmEdge(elemA, elemB)
-
-        END SELECT
-
-      END SELECT
-
-    END SUBROUTINE connectedVolEdge
-
-    SUBROUTINE connectedSegmEdge(elemA, elemB)
-      IMPLICIT NONE
-
-      CLASS(meshVol1DCartSegm), INTENT(inout), TARGET:: elemA
-      CLASS(meshEdge1DCart),    INTENT(inout), TARGET:: elemB
-
-      IF (.NOT. ASSOCIATED(elemA%e1) .AND. &
-          elemA%n2%n == elemB%n1%n) THEN
-
-        elemA%e1 => elemB
-        elemB%e2 => elemA
- 
-        !Revers the normal to point inside the domain
-        elemB%normal = - elemB%normal
-
-      END IF
-
-      IF (.NOT. ASSOCIATED(elemA%e2) .AND. &
-          elemA%n1%n == elemB%n1%n) THEN
-
-        elemA%e2 => elemB
-        elemB%e1 => elemA
-
-      END IF
-
-    END SUBROUTINE connectedSegmEdge
-
-    SUBROUTINE constructGlobalK(K, elem)
-      IMPLICIT NONE
-
-      REAL(8), INTENT(inout):: K(1:,1:)
-      CLASS(meshVol), INTENT(in):: elem
-      REAL(8):: localK(1:2,1:2)
-      INTEGER:: i, j
-      INTEGER:: n(1:2)
-
-      SELECT TYPE(elem)
-      TYPE IS(meshVol1DCartSegm)
-        localK = elem%elemK()
-        n = (/ elem%n1%n, elem%n2%n /)
-
-      CLASS DEFAULT
-        n = 0
-        localK = 0.D0
-
-      END SELECT
-      
-      DO i = 1, 2
-        DO j = 1, 2
-          K(n(i), n(j)) = K(n(i), n(j)) + localK(i, j)
-        END DO
-      END DO
-
-    END SUBROUTINE constructGlobalK
-
-END MODULE moduleMesh1DCartRead

src/modules/mesh/1DRad/makefile

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

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

@@ -72,7 +72,7 @@ MODULE moduleMesh1DRad
     !Connectivity to nodes
     CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL()
     !Connectivity to adjacent elements
-    CLASS(*), POINTER:: e1 => NULL(), e2 => NULL()
+    CLASS(meshElement), POINTER:: e1 => NULL(), e2 => NULL()
     REAL(8):: arNodes(1:2)
     CONTAINS
       PROCEDURE, PASS::   init        => initVol1DRadSegm
@@ -200,18 +200,19 @@ MODULE moduleMesh1DRad
     !VOLUME FUNCTIONS
     !SEGMENT FUNCTIONS
     !Init segment element
-    SUBROUTINE initVol1DRadSegm(self, n, p)
+    SUBROUTINE initVol1DRadSegm(self, n, p, nodes)
       USE moduleRefParam
       IMPLICIT NONE
 
       CLASS(meshVol1DRadSegm), INTENT(out):: self
       INTEGER, INTENT(in):: n
       INTEGER, INTENT(in):: p(:)
+      TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:)
       REAL(8), DIMENSION(1:3):: r1, r2
 
       self%n = n
-      self%n1 => mesh%nodes(p(1))%obj
+      self%n1 => nodes(p(1))%obj
-      self%n2 => mesh%nodes(p(2))%obj
+      self%n2 => nodes(p(2))%obj
       !Get element coordinates
       r1 = self%n1%getCoordinates()
       r2 = self%n2%getCoordinates()
@@ -312,19 +313,20 @@ MODULE moduleMesh1DRad
     END SUBROUTINE partialDerRad
 
     !Computes local stiffness matrix
-    PURE FUNCTION elemKRad(self) RESULT(ke)
+    PURE FUNCTION elemKRad(self) RESULT(localK)
       USE moduleConstParam, ONLY: PI2
       IMPLICIT NONE
 
       CLASS(meshVol1DRadSegm), INTENT(in):: self
-      REAL(8):: ke(1:2,1:2)
+      REAL(8), ALLOCATABLE:: localK(:,:)
       REAL(8):: Xii(1:3)
       REAL(8):: dPsi(1:1, 1:2)
       REAL(8):: invJ(1), detJ
       REAL(8):: r, fPsi(1:2)
       INTEGER:: l
 
-      ke      = 0.D0
+      ALLOCATE(localK(1:2, 1:2))
+      localK  = 0.D0
       Xii     = 0.D0
       DO l = 1, 3
         xii(1) = corSeg(l)
@@ -333,13 +335,13 @@ MODULE moduleMesh1DRad
         invJ    = self%invJac(Xii, dPsi)
         fPsi    = self%fPsi(Xii)
         r       = DOT_PRODUCT(fPsi, self%r)
-        ke      = ke + MATMUL(RESHAPE(MATMUL(invJ,dPsi), (/ 2, 1/)), &
+        localK  = localK + MATMUL(RESHAPE(MATMUL(invJ,dPsi), (/ 2, 1/)), &
-                              RESHAPE(MATMUL(invJ,dPsi), (/ 1, 2/)))* &
+                                  RESHAPE(MATMUL(invJ,dPsi), (/ 1, 2/)))* &
-                       r*wSeg(l)/detJ
+                           r*wSeg(l)/detJ
 
       END DO
 
-      ke      = ke*PI2
+      localK = localK*PI2
 
     END FUNCTION elemKRad
 
@@ -406,12 +408,12 @@ MODULE moduleMesh1DRad
       w_p = self%weight(part%xi)
       tensorS = outerProduct(part%v, part%v)
 
-      vertex => self%n1%output(part%sp)
+      vertex => self%n1%output(part%species%n)
       vertex%den          = vertex%den          + part%weight*w_p(1)
       vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(1)*part%v(:)
       vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(1)*tensorS
 
-      vertex => self%n2%output(part%sp)
+      vertex => self%n2%output(part%species%n)
       vertex%den          = vertex%den          + part%weight*w_p(2)
       vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(2)*part%v(:)
       vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(2)*tensorS
@@ -470,7 +472,7 @@ MODULE moduleMesh1DRad
 
       CLASS(meshVol1DRadSegm), INTENT(in):: self
       REAL(8), INTENT(in):: xi(1:3)
-      CLASS(*), POINTER, INTENT(out):: nextElement
+      CLASS(meshElement), POINTER, INTENT(out):: nextElement
       
       NULLIFY(nextElement)
       IF     (xi(1) < -1.D0) THEN
@@ -532,5 +534,121 @@ MODULE moduleMesh1DRad
 
     END FUNCTION invJ1DRad
 
+    SUBROUTINE connectMesh1DRad(self)
+      IMPLICIT NONE
+
+      CLASS(meshGeneric), INTENT(inout):: self
+      INTEGER:: e, et
+
+      DO e = 1, self%numVols
+        !Connect Vol-Vol
+        DO et = 1, self%numVols
+          IF (e /= et) THEN
+            CALL connectVolVol(self%vols(e)%obj, self%vols(et)%obj)
+
+          END IF
+
+        END DO
+
+        SELECT TYPE(self)
+        TYPE IS(meshParticles)
+          !Connect Vol-Edge
+          DO et = 1, self%numEdges
+            CALL connectVolEdge(self%vols(e)%obj, self%edges(et)%obj)
+
+          END DO
+
+        END SELECT
+
+      END DO
+
+    END SUBROUTINE connectMesh1DRad
+
+    SUBROUTINE connectVolVol(elemA, elemB)
+      IMPLICIT NONE
+
+      CLASS(meshVol), INTENT(inout):: elemA
+      CLASS(meshVol), INTENT(inout):: elemB
+
+        SELECT TYPE(elemA)
+        TYPE IS(meshVol1DRadSegm)
+          SELECT TYPE(elemB)
+          TYPE IS(meshVol1DRadSegm)
+            CALL connectSegmSegm(elemA, elemB)
+
+          END SELECT
+
+        END SELECT
+
+    END SUBROUTINE connectVolVol
+
+    SUBROUTINE connectSegmSegm(elemA, elemB)
+      IMPLICIT NONE
+
+      CLASS(meshVol1DRadSegm), INTENT(inout), TARGET:: elemA
+      CLASS(meshVol1DRadSegm), INTENT(inout), TARGET:: elemB
+
+      IF (.NOT. ASSOCIATED(elemA%e1) .AND. &
+          elemA%n2%n == elemB%n1%n) THEN
+        elemA%e1 => elemB
+        elemB%e2 => elemA
+        
+      END IF
+
+      IF (.NOT. ASSOCIATED(elemA%e2) .AND. &
+          elemA%n1%n == elemB%n2%n) THEN
+
+        elemA%e2 => elemB
+        elemB%e1 => elemA
+
+      END IF
+
+    END SUBROUTINE connectSegmSegm
+
+    SUBROUTINE connectVolEdge(elemA, elemB)
+      IMPLICIT NONE
+
+      CLASS(meshVol),  INTENT(inout):: elemA
+      CLASS(meshEdge), INTENT(inout):: elemB
+
+      SELECT TYPE(elemA)
+      TYPE IS (meshVol1DRadSegm)
+        SELECT TYPE(elemB)
+        CLASS IS(meshEdge1DRad)
+          CALL connectSegmEdge(elemA, elemB)
+
+        END SELECT
+
+      END SELECT
+
+    END SUBROUTINE connectVolEdge
+
+    SUBROUTINE connectSegmEdge(elemA, elemB)
+      IMPLICIT NONE
+
+      CLASS(meshVol1DRadSegm), INTENT(inout), TARGET:: elemA
+      CLASS(meshEdge1DRad),    INTENT(inout), TARGET:: elemB
+
+      IF (.NOT. ASSOCIATED(elemA%e1) .AND. &
+          elemA%n2%n == elemB%n1%n) THEN
+
+        elemA%e1 => elemB
+        elemB%e2 => elemA
+
+        !Revers the normal to point inside the domain
+        elemB%normal = - elemB%normal
+
+      END IF
+
+      IF (.NOT. ASSOCIATED(elemA%e2) .AND. &
+          elemA%n1%n == elemB%n1%n) THEN
+
+        elemA%e2 => elemB
+        elemB%e1 => elemA
+
+      END IF
+
+    END SUBROUTINE connectSegmEdge
+
 END MODULE moduleMesh1DRad
 

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

@@ -1,264 +0,0 @@
-MODULE moduleMesh1DRadRead
-  USE moduleMesh
-  USE moduleMesh1DRad
-  
-  !TODO: make this abstract to allow different mesh formats
-  TYPE, EXTENDS(meshGeneric):: mesh1DRadGeneric
-    CONTAINS
-      PROCEDURE, PASS:: init     => init1DRadMesh
-      PROCEDURE, PASS:: readMesh => readMesh1DRad
-
-  END TYPE
-
-  INTERFACE connected
-    MODULE PROCEDURE connectedVolVol, connectedVolEdge
-
-  END INTERFACE connected
-
-  CONTAINS
-    !Init 1D mesh
-    SUBROUTINE init1DRadMesh(self, meshFormat)
-      USE moduleMesh
-      USE moduleErrors
-      IMPLICIT NONE
-
-      CLASS(mesh1DRadGeneric), INTENT(out):: self
-      CHARACTER(:), ALLOCATABLE, INTENT(in):: meshFormat
-
-      SELECT CASE(meshFormat)
-      CASE ("gmsh")
-        self%printOutput => printOutputGmsh
-        self%printColl   => printCollGmsh
-        self%printEM     => printEMGmsh
-
-      CASE DEFAULT
-        CALL criticalError("Mesh type " // meshFormat // " not supported.", "init1DRad")
-        
-      END SELECT
-
-    END SUBROUTINE init1DRadMesh
-
-    !Reads 1D mesh
-    SUBROUTINE readMesh1DRad(self, filename)
-      USE  moduleBoundary
-      IMPLICIT NONE
-
-      CLASS(mesh1DRadGeneric), INTENT(inout):: self
-      CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
-      REAL(8):: x
-      INTEGER:: p(1:2)
-      INTEGER:: e, et, n, eTemp, elemType, bt
-      INTEGER:: totalNumElem
-      INTEGER:: boundaryType
-
-      !Open file mesh
-      OPEN(10, FILE=TRIM(filename))
-      !Skip header
-      READ(10, *)
-      READ(10, *)
-      READ(10, *)
-      READ(10, *)
-      !Read number of nodes
-      READ(10, *) self%numNodes
-      !Allocate required matrices and vectors
-      ALLOCATE(self%nodes(1:self%numNodes))
-      ALLOCATE(self%K(1:self%numNodes, 1:self%numNodes))
-      ALLOCATE(self%IPIV(1:self%numNodes, 1:self%numNodes))
-      self%K = 0.D0
-      self%IPIV = 0
-      !Read nodes coordinates. Only relevant for x
-      DO e = 1, self%numNodes
-        READ(10, *) n, x
-        ALLOCATE(meshNode1DRad:: self%nodes(n)%obj)
-        CALL self%nodes(n)%obj%init(n, (/ x, 0.D0, 0.D0 /))
-
-      END DO
-      !Skips comments
-      READ(10, *)
-      READ(10, *)
-      !Reads the total number of elements (edges+vol)
-      READ(10, *) totalNumElem
-      self%numEdges = 0
-      DO e = 1, totalNumElem
-        READ(10, *) eTemp, elemType
-        IF (elemType == 15) THEN !15 is physical node in GMSH
-          self%numEdges = e
-
-        END IF
-
-      END DO
-
-      !Substract the number of edges to the total number of elements
-      !to obtain the number of volume elements
-      self%numVols = totalNumelem - self%numEdges
-      !Allocates arrays
-      ALLOCATE(self%edges(1:self%numEdges))
-      ALLOCATE(self%vols(1:self%numVols))
-
-      !Go back to the beginning of reading elements
-      DO e = 1, totalNumelem
-        BACKSPACE(10)
-
-      END DO
-
-      !Reads edges
-      DO e = 1, self%numEdges
-        READ(10, *) n, elemType, eTemp, boundaryType, eTemp, p(1)
-        !Associate boundary condition 
-        bt = getBoundaryId(boundaryType)
-
-        ALLOCATE(meshEdge1DRad:: self%edges(e)%obj)
-
-        CALL self%edges(e)%obj%init(n, p(1:1), bt, boundaryType)
-
-      END DO
-
-      !Read and initialize volumes
-      DO e = 1, self%numVols
-        READ(10, *) n, elemType, eTemp, eTemp, eTemp, p(1:2)
-        ALLOCATE(meshVol1DRadSegm:: self%vols(e)%obj)
-        CALL self%vols(e)%obj%init(n - self%numEdges, p(1:2))
-
-      END DO
-
-      CLOSE(10)
-
-      !Build connectivity between elements
-      DO e = 1, self%numVols
-        !Connectivity between volumes
-        DO et = 1, self%numVols
-          IF (e /= et) THEN
-            CALL connected(self%vols(e)%obj, self%vols(et)%obj)
-
-          END IF
-
-        END DO
-
-        !Connectivity betwen vols and edges
-        DO et = 1, self%numEdges
-          CALL connected(self%vols(e)%obj, self%edges(et)%obj)
-
-        END DO
-
-        !Constructs the global K matrix
-        CALL constructGlobalK(self%K, self%vols(e)%obj)
-
-      END DO
-
-    END SUBROUTINE readMesh1DRad
-
-    SUBROUTINE connectedVolVol(elemA, elemB)
-      IMPLICIT NONE
-
-      CLASS(meshVol), INTENT(inout):: elemA
-      CLASS(meshVol), INTENT(inout):: elemB
-
-        SELECT TYPE(elemA)
-        TYPE IS(meshVol1DRadSegm)
-          SELECT TYPE(elemB)
-          TYPE IS(meshVol1DRadSegm)
-            CALL connectedSegmSegm(elemA, elemB)
-
-          END SELECT
-
-        END SELECT
-
-    END SUBROUTINE connectedVolVol
-
-    SUBROUTINE connectedSegmSegm(elemA, elemB)
-      IMPLICIT NONE
-
-      CLASS(meshVol1DRadSegm), INTENT(inout), TARGET:: elemA
-      CLASS(meshVol1DRadSegm), INTENT(inout), TARGET:: elemB
-
-      IF (.NOT. ASSOCIATED(elemA%e1) .AND. &
-          elemA%n2%n == elemB%n1%n) THEN
-        elemA%e1 => elemB
-        elemB%e2 => elemA
-        
-      END IF
-
-      IF (.NOT. ASSOCIATED(elemA%e2) .AND. &
-          elemA%n1%n == elemB%n2%n) THEN
-
-        elemA%e2 => elemB
-        elemB%e1 => elemA
-
-      END IF
-
-    END SUBROUTINE connectedSegmSegm
-
-    SUBROUTINE connectedVolEdge(elemA, elemB)
-      IMPLICIT NONE
-
-      CLASS(meshVol),  INTENT(inout):: elemA
-      CLASS(meshEdge), INTENT(inout):: elemB
-
-      SELECT TYPE(elemA)
-      TYPE IS (meshVol1DRadSegm)
-        SELECT TYPE(elemB)
-        CLASS IS(meshEdge1DRad)
-          CALL connectedSegmEdge(elemA, elemB)
-
-        END SELECT
-
-      END SELECT
-
-    END SUBROUTINE connectedVolEdge
-
-    SUBROUTINE connectedSegmEdge(elemA, elemB)
-      IMPLICIT NONE
-
-      CLASS(meshVol1DRadSegm), INTENT(inout), TARGET:: elemA
-      CLASS(meshEdge1DRad),    INTENT(inout), TARGET:: elemB
-
-      IF (.NOT. ASSOCIATED(elemA%e1) .AND. &
-          elemA%n2%n == elemB%n1%n) THEN
-
-        elemA%e1 => elemB
-        elemB%e2 => elemA
-
-        !Revers the normal to point inside the domain
-        elemB%normal = - elemB%normal
-
-      END IF
-
-      IF (.NOT. ASSOCIATED(elemA%e2) .AND. &
-          elemA%n1%n == elemB%n1%n) THEN
-
-        elemA%e2 => elemB
-        elemB%e1 => elemA
-
-      END IF
-
-    END SUBROUTINE connectedSegmEdge
-
-    SUBROUTINE constructGlobalK(K, elem)
-      IMPLICIT NONE
-
-      REAL(8), INTENT(inout):: K(1:,1:)
-      CLASS(meshVol), INTENT(in):: elem
-      REAL(8):: localK(1:2,1:2)
-      INTEGER:: i, j
-      INTEGER:: n(1:2)
-
-      SELECT TYPE(elem)
-      TYPE IS(meshVol1DRadSegm)
-        localK = elem%elemK()
-        n = (/ elem%n1%n, elem%n2%n /)
-
-      CLASS DEFAULT
-        n = 0
-        localK = 0.D0
-
-      END SELECT
-      
-      DO i = 1, 2
-        DO j = 1, 2
-          K(n(i), n(j)) = K(n(i), n(j)) + localK(i, j)
-        END DO
-      END DO
-
-    END SUBROUTINE constructGlobalK
-
-END MODULE moduleMesh1DRadRead

src/modules/mesh/2DCart/makefile

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

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

@@ -77,7 +77,7 @@ MODULE moduleMesh2DCart
     !Connectivity to nodes
     CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL(), n3 => NULL(), n4 => NULL()
     !Connectivity to adjacent elements
-    CLASS(*), POINTER:: e1 => NULL(), e2 => NULL(), e3 => NULL(), e4 => NULL()
+    CLASS(meshElement), POINTER:: e1 => NULL(), e2 => NULL(), e3 => NULL(), e4 => NULL()
     REAL(8):: arNodes(1:4) = 0.D0
     CONTAINS
       PROCEDURE, PASS::   init        => initVolQuad2DCart
@@ -107,7 +107,7 @@ MODULE moduleMesh2DCart
     !Connectivity to nodes
     CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL(), n3 => NULL()
     !Connectivity to adjacent elements
-    CLASS(*), POINTER:: e1 => NULL(), e2 => NULL(), e3 => NULL()
+    CLASS(meshElement), POINTER:: e1 => NULL(), e2 => NULL(), e3 => NULL()
     REAL(8):: arNodes(1:3) = 0.D0
 
     CONTAINS
@@ -283,20 +283,21 @@ MODULE moduleMesh2DCart
     !VOLUME FUNCTIONS
     !QUAD FUNCTIONS
     !Inits quadrilateral element
-    SUBROUTINE initVolQuad2DCart(self, n, p)
+    SUBROUTINE initVolQuad2DCart(self, n, p, nodes)
       USE moduleRefParam
       IMPLICIT NONE
 
       CLASS(meshVol2DCartQuad), INTENT(out):: self
       INTEGER, INTENT(in):: n
       INTEGER, INTENT(in):: p(:)
+      TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:)
       REAL(8), DIMENSION(1:3):: r1, r2, r3, r4
 
       self%n  = n
-      self%n1 => mesh%nodes(p(1))%obj
+      self%n1 => nodes(p(1))%obj
-      self%n2 => mesh%nodes(p(2))%obj
+      self%n2 => nodes(p(2))%obj
-      self%n3 => mesh%nodes(p(3))%obj
+      self%n3 => nodes(p(3))%obj
-      self%n4 => mesh%nodes(p(4))%obj
+      self%n4 => nodes(p(4))%obj
       !Get element coordinates
       r1 = self%n1%getCoordinates()
       r2 = self%n2%getCoordinates()
@@ -414,17 +415,18 @@ MODULE moduleMesh2DCart
     END SUBROUTINE partialDerQuad
 
     !Computes element local stiffness matrix
-    PURE FUNCTION elemKQuad(self) RESULT(ke)
+    PURE FUNCTION elemKQuad(self) RESULT(localK)
       IMPLICIT NONE
 
       CLASS(meshVol2DCartQuad), INTENT(in):: self
+      REAL(8), ALLOCATABLE:: localK(:,:)
       REAL(8):: xi(1:3)
       REAL(8):: fPsi(1:4), dPsi(1:2,1:4)
-      REAL(8):: ke(1:4,1:4)
       REAL(8):: invJ(1:2,1:2), detJ
       INTEGER:: l, m
 
-      ke=0.D0
+      ALLOCATE(localK(1:4, 1:4))
+      localK=0.D0
       xi=0.D0
       !Start 2D Gauss Quad Integral
       DO l=1, 3
@@ -436,7 +438,7 @@ MODULE moduleMesh2DCart
           fPsi      = self%fPsi(xi)
           detJ      = self%detJac(xi,dPsi)
           invJ      = self%invJac(xi,dPsi)
-          ke        = ke + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*wQuad(l)*wQuad(m)/detJ
+          localK    = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*wQuad(l)*wQuad(m)/detJ
 
         END DO
       END DO
@@ -510,22 +512,22 @@ MODULE moduleMesh2DCart
       w_p = self%weight(part%xi)
       tensorS = outerProduct(part%v, part%v)
 
-      vertex => self%n1%output(part%sp)
+      vertex => self%n1%output(part%species%n)
       vertex%den          = vertex%den          + part%weight*w_p(1)
       vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(1)*part%v(:)
       vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(1)*tensorS
 
-      vertex => self%n2%output(part%sp)
+      vertex => self%n2%output(part%species%n)
       vertex%den          = vertex%den          + part%weight*w_p(2)
       vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(2)*part%v(:)
       vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(2)*tensorS
 
-      vertex => self%n3%output(part%sp)
+      vertex => self%n3%output(part%species%n)
       vertex%den          = vertex%den          + part%weight*w_p(3)
       vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(3)*part%v(:)
       vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(3)*tensorS
 
-      vertex => self%n4%output(part%sp)
+      vertex => self%n4%output(part%species%n)
       vertex%den          = vertex%den          + part%weight*w_p(4)
       vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(4)*part%v(:)
       vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(4)*tensorS
@@ -607,7 +609,7 @@ MODULE moduleMesh2DCart
 
       CLASS(meshVol2DCartQuad), INTENT(in):: self
       REAL(8), INTENT(in):: xi(1:3)
-      CLASS(*), POINTER, INTENT(out):: nextElement
+      CLASS(meshElement), POINTER, INTENT(out):: nextElement
       REAL(8):: xiArray(1:4)
       INTEGER:: nextInt
 
@@ -630,23 +632,23 @@ MODULE moduleMesh2DCart
 
     !TRIA ELEMENT
     !Init tria element
-    SUBROUTINE initVolTria2DCart(self, n, p)
+    SUBROUTINE initVolTria2DCart(self, n, p, nodes)
       USE moduleRefParam
       IMPLICIT NONE
 
       CLASS(meshVol2DCartTria), INTENT(out):: self
       INTEGER, INTENT(in):: n
       INTEGER, INTENT(in):: p(:)
+      TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:)
       REAL(8), DIMENSION(1:3):: r1, r2, r3
-      REAL(8):: A
 
       !Assign node index
       self%n = n
 
       !Assign nodes to element
-      self%n1 => mesh%nodes(p(1))%obj
+      self%n1 => nodes(p(1))%obj
-      self%n2 => mesh%nodes(p(2))%obj
+      self%n2 => nodes(p(2))%obj
-      self%n3 => mesh%nodes(p(3))%obj
+      self%n3 => nodes(p(3))%obj
       !Get element coordinates
       r1 = self%n1%getCoordinates()
       r2 = self%n2%getCoordinates()
@@ -777,27 +779,28 @@ MODULE moduleMesh2DCart
     END SUBROUTINE partialDerTria
 
     !Computes element local stiffness matrix
-    PURE FUNCTION elemKTria(self) RESULT(ke)
+    PURE FUNCTION elemKTria(self) RESULT(localK)
       IMPLICIT NONE
 
       CLASS(meshVol2DCartTria), INTENT(in):: self
+      REAL(8), ALLOCATABLE:: localK(:,:)
       REAL(8):: xi(1:3)
       REAL(8):: fPsi(1:3), dPsi(1:2,1:3)
-      REAL(8):: ke(1:3,1:3)
       REAL(8):: invJ(1:2,1:2), detJ
       INTEGER:: l
 
-      ke=0.D0
+      ALLOCATE(localK(1:4, 1:4))
+      localK=0.D0
       xi=0.D0
       !Start 2D Gauss Quad Integral
       DO l=1, 4
-        xi(1) = xi1Tria(l)
+        xi(1)  = xi1Tria(l)
-        xi(2) = xi2Tria(l)
+        xi(2)  = xi2Tria(l)
-        dPsi  = self%dPsi(xi)
+        dPsi   = self%dPsi(xi)
-        detJ  = self%detJac(xi,dPsi)
+        detJ   = self%detJac(xi,dPsi)
-        invJ  = self%invJac(xi,dPsi)
+        invJ   = self%invJac(xi,dPsi)
-        fPsi  = self%fPsi(xi)
+        fPsi   = self%fPsi(xi)
-        ke    = ke + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*wTria(l)/detJ
+        localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*wTria(l)/detJ
 
       END DO
 
@@ -869,17 +872,17 @@ MODULE moduleMesh2DCart
       w_p = self%weight(part%xi)
       tensorS = outerProduct(part%v, part%v)
 
-      vertex => self%n1%output(part%sp)
+      vertex => self%n1%output(part%species%n)
       vertex%den          = vertex%den          + part%weight*w_p(1)
       vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(1)*part%v(:)
       vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(1)*tensorS
 
-      vertex => self%n2%output(part%sp)
+      vertex => self%n2%output(part%species%n)
       vertex%den          = vertex%den          + part%weight*w_p(2)
       vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(2)*part%v(:)
       vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(2)*tensorS
 
-      vertex => self%n3%output(part%sp)
+      vertex => self%n3%output(part%species%n)
       vertex%den          = vertex%den          + part%weight*w_p(3)
       vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(3)*part%v(:)
       vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(3)*tensorS
@@ -950,7 +953,7 @@ MODULE moduleMesh2DCart
 
       CLASS(meshVol2DCartTria), INTENT(in):: self
       REAL(8), INTENT(in):: xi(1:3)
-      CLASS(*), POINTER, INTENT(out):: nextElement
+      CLASS(meshElement), POINTER, INTENT(out):: nextElement
       REAL(8):: xiArray(1:3)
       INTEGER:: nextInt
 
@@ -1018,4 +1021,453 @@ MODULE moduleMesh2DCart
 
     END FUNCTION invJ2DCart
 
+    SUBROUTINE connectMesh2DCart(self)
+      IMPLICIT NONE
+
+      CLASS(meshGeneric), INTENT(inout):: self
+      INTEGER:: e, et
+
+      DO e = 1, self%numVols
+        !Connect Vol-Vol
+        DO et = 1, self%numVols
+          IF (e /= et) THEN
+            CALL connectVolVol(self%vols(e)%obj, self%vols(et)%obj)
+
+          END IF
+
+        END DO
+
+        SELECT TYPE(self)
+        TYPE IS(meshParticles)
+          !Connect Vol-Edge
+          DO et = 1, self%numEdges
+            CALL connectVolEdge(self%vols(e)%obj, self%edges(et)%obj)
+
+          END DO
+
+        END SELECT
+
+      END DO
+
+    END SUBROUTINE connectMesh2DCart
+
+    !Selects type of elements to build connection
+    SUBROUTINE connectVolVol(elemA, elemB)
+      IMPLICIT NONE
+
+      CLASS(meshVol), INTENT(inout):: elemA
+      CLASS(meshVol), INTENT(inout):: elemB
+
+      SELECT TYPE(elemA)
+      TYPE IS(meshVol2DCartQuad)
+        !Element A is a quadrilateral
+        SELECT TYPE(elemB)
+        TYPE IS(meshVol2DCartQuad)
+          !Element B is a quadrilateral
+          CALL connectQuadQuad(elemA, elemB)
+
+        TYPE IS(meshVol2DCartTria)
+          !Element B is a triangle
+          CALL connectQuadTria(elemA, elemB)
+
+        END SELECT
+
+      TYPE IS(meshVol2DCartTria)
+        !Element A is a Triangle
+        SELECT TYPE(elemB)
+        TYPE IS(meshVol2DCartQuad)
+          !Element B is a quadrilateral
+          CALL connectQuadTria(elemB, elemA)
+
+        TYPE IS(meshVol2DCartTria)
+          !Element B is a triangle
+          CALL connectTriaTria(elemA, elemB)
+
+        END SELECT
+
+      END SELECT
+
+    END SUBROUTINE connectVolVol
+
+    SUBROUTINE connectVolEdge(elemA, elemB)
+      IMPLICIT NONE
+
+      CLASS(meshVol),  INTENT(inout):: elemA
+      CLASS(meshEdge), INTENT(inout):: elemB
+
+      SELECT TYPE(elemB)
+      CLASS IS(meshEdge2DCart)
+        SELECT TYPE(elemA)
+        TYPE IS(meshVol2DCartQuad)
+          !Element A is a quadrilateral
+          CALL connectQuadEdge(elemA, elemB)
+
+        TYPE IS(meshVol2DCartTria)
+          !Element A is a triangle
+          CALL connectTriaEdge(elemA, elemB)
+
+        END SELECT
+
+      END SELECT
+
+    END SUBROUTINE connectVolEdge
+
+    SUBROUTINE connectQuadQuad(elemA, elemB)
+      IMPLICIT NONE
+
+      CLASS(meshVol2DCartQuad), INTENT(inout), TARGET:: elemA
+      CLASS(meshVol2DCartQuad), INTENT(inout), TARGET:: elemB
+
+      !Check direction 1
+      IF (.NOT. ASSOCIATED(elemA%e1) .AND. &
+          elemA%n1%n == elemB%n4%n   .AND. &
+          elemA%n2%n == elemB%n3%n) THEN
+        elemA%e1 => elemB
+        elemB%e3 => elemA
+
+      END IF
+
+      !Check direction 2
+      IF (.NOT. ASSOCIATED(elemA%e2) .AND. &
+          elemA%n2%n == elemB%n1%n   .AND. &
+          elemA%n3%n == elemB%n4%n) THEN
+        elemA%e2 => elemB
+        elemB%e4 => elemA
+
+      END IF
+
+      !Check direction 3
+      IF (.NOT. ASSOCIATED(elemA%e3) .AND. &
+          elemA%n3%n == elemB%n2%n   .AND. &
+          elemA%n4%n == elemB%n1%n) THEN
+        elemA%e3 => elemB
+        elemB%e1 => elemA
+
+      END IF
+
+      !Check direction 4
+      IF (.NOT. ASSOCIATED(elemA%e4) .AND. &
+          elemA%n4%n == elemB%n3%n   .AND. &
+          elemA%n1%n == elemB%n2%n) THEN
+        elemA%e4 => elemB
+        elemB%e2 => elemA
+
+      END IF
+
+    END SUBROUTINE connectQuadQuad
+
+    SUBROUTINE connectQuadTria(elemA, elemB)
+      IMPLICIT NONE
+
+      CLASS(meshVol2DCartQuad), INTENT(inout), TARGET:: elemA
+      CLASS(meshVol2DCartTria), INTENT(inout), TARGET:: elemB
+
+      !Check direction 1
+      IF (.NOT. ASSOCIATED(elemA%e1)) THEN
+        IF     (elemA%n1%n == elemB%n1%n .AND. &
+                elemA%n2%n == elemB%n3%n) THEN
+          elemA%e1 => elemB
+          elemB%e3 => elemA
+
+        ELSEIF (elemA%n1%n == elemB%n3%n .AND. &
+                elemA%n2%n == elemB%n2%n) THEN
+          elemA%e1 => elemB
+          elemB%e2 => elemA
+
+        ELSEIF (elemA%n1%n == elemB%n2%n .AND. &
+                elemA%n2%n == elemB%n1%n) THEN
+          elemA%e1 => elemB
+          elemB%e1 => elemA
+
+        END IF
+
+      END IF
+
+      !Check direction 2
+      IF (.NOT. ASSOCIATED(elemA%e2)) THEN
+        IF     (elemA%n2%n == elemB%n1%n .AND. &
+                elemA%n3%n == elemB%n3%n) THEN
+          elemA%e2 => elemB
+          elemB%e3 => elemA
+
+        ELSEIF (elemA%n2%n == elemB%n3%n .AND. &
+                elemA%n3%n == elemB%n2%n) THEN
+          elemA%e2 => elemB
+          elemB%e2 => elemA
+
+        ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
+                elemA%n3%n == elemB%n1%n) THEN
+          elemA%e2 => elemB
+          elemB%e1 => elemA
+
+        END IF
+
+      END IF
+
+      !Check direction 3
+      IF (.NOT. ASSOCIATED(elemA%e3)) THEN
+        IF     (elemA%n3%n == elemB%n1%n .AND. &
+                elemA%n4%n == elemB%n3%n) THEN
+          elemA%e3 => elemB
+          elemB%e3 => elemA
+
+        ELSEIF (elemA%n3%n == elemB%n3%n .AND. &
+                elemA%n4%n == elemB%n2%n) THEN
+          elemA%e3 => elemB
+          elemB%e2 => elemA
+
+        ELSEIF (elemA%n3%n == elemB%n2%n .AND. &
+                elemA%n4%n == elemB%n1%n) THEN
+          elemA%e3 => elemB
+          elemB%e1 => elemA
+
+        END IF
+
+      END IF
+
+      !Check direction 4
+      IF (.NOT. ASSOCIATED(elemA%e4)) THEN
+        IF     (elemA%n4%n == elemB%n1%n .AND. &
+                elemA%n1%n == elemB%n3%n) THEN
+          elemA%e4 => elemB
+          elemB%e3 => elemA
+
+        ELSEIF (elemA%n4%n == elemB%n3%n .AND. &
+                elemA%n1%n == elemB%n2%n) THEN
+          elemA%e4 => elemB
+          elemB%e2 => elemA
+
+        ELSEIF (elemA%n4%n == elemB%n2%n .AND. &
+                elemA%n1%n == elemB%n1%n) THEN
+          elemA%e4 => elemB
+          elemB%e1 => elemA
+
+        END IF
+
+      END IF
+
+    END SUBROUTINE connectQuadTria
+
+    SUBROUTINE connectTriaTria(elemA, elemB)
+      IMPLICIT NONE
+
+      CLASS(meshVol2DCartTria), INTENT(inout), TARGET:: elemA
+      CLASS(meshVol2DCartTria), INTENT(inout), TARGET:: elemB
+
+      !Check direction 1
+      IF (.NOT. ASSOCIATED(elemA%e1)) THEN
+        IF     (elemA%n1%n == elemB%n1%n .AND. &
+                elemA%n2%n == elemB%n3%n) THEN
+          elemA%e1 => elemB
+          elemB%e3 => elemA
+
+        ELSEIF (elemA%n1%n == elemB%n2%n .AND. &
+                elemA%n2%n == elemB%n1%n) THEN
+          elemA%e1 => elemB
+          elemB%e1 => elemA
+
+        ELSEIF (elemA%n1%n == elemB%n3%n .AND. &
+                elemA%n2%n == elemB%n2%n) THEN
+          elemA%e1 => elemB
+          elemB%e2 => elemA
+
+        END IF
+
+      END IF
+
+      !Check direction 2
+      IF (.NOT. ASSOCIATED(elemA%e2)) THEN
+        IF     (elemA%n2%n == elemB%n1%n .AND. &
+                elemA%n3%n == elemB%n3%n) THEN
+          elemA%e2 => elemB
+          elemB%e3 => elemA
+
+        ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
+                elemA%n3%n == elemB%n1%n) THEN
+          elemA%e2 => elemB
+          elemB%e1 => elemA
+
+        ELSEIF (elemA%n2%n == elemB%n3%n .AND. &
+                elemA%n3%n == elemB%n2%n) THEN
+          elemA%e2 => elemB
+          elemB%e2 => elemA
+
+        END IF
+
+
+      END IF
+
+      !Check direction 3
+      IF (.NOT. ASSOCIATED(elemA%e3)) THEN
+        IF     (elemA%n3%n == elemB%n1%n .AND. &
+                elemA%n1%n == elemB%n3%n) THEN
+          elemA%e3 => elemB
+          elemB%e3 => elemA
+
+        ELSEIF (elemA%n3%n == elemB%n2%n .AND. &
+                elemA%n1%n == elemB%n1%n) THEN
+          elemA%e3 => elemB
+          elemB%e1 => elemA
+
+        ELSEIF (elemA%n3%n == elemB%n3%n .AND. &
+                elemA%n1%n == elemB%n2%n) THEN
+          elemA%e3 => elemB
+          elemB%e2 => elemA
+
+        END IF
+
+
+      END IF
+
+    END SUBROUTINE connectTriaTria
+
+    SUBROUTINE connectQuadEdge(elemA, elemB)
+      IMPLICIT NONE
+
+      CLASS(meshVol2DCartQuad), INTENT(inout), TARGET:: elemA
+      CLASS(meshEdge2DCart),    INTENT(inout), TARGET:: elemB
+
+      !Check direction 1
+      IF (.NOT. ASSOCIATED(elemA%e1)) THEN
+        IF     (elemA%n1%n == elemB%n1%n .AND. &
+                elemA%n2%n == elemB%n2%n) THEN
+          elemA%e1 => elemB
+          elemB%e1 => elemA
+
+        ELSEIF (elemA%n1%n == elemB%n2%n .AND. &
+                elemA%n2%n == elemB%n1%n) THEN
+          elemA%e1 => elemB
+          elemB%e2 => elemA
+
+          !Revers the normal to point inside the domain
+          elemB%normal = - elemB%normal
+
+        END IF
+
+      END IF
+
+      !Check direction 2
+      IF (.NOT. ASSOCIATED(elemA%e2)) THEN
+        IF     (elemA%n2%n == elemB%n1%n .AND. &
+                elemA%n3%n == elemB%n2%n) THEN
+          elemA%e2 => elemB
+          elemB%e1 => elemA
+
+        ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
+                elemA%n3%n == elemB%n1%n) THEN
+          elemA%e2 => elemB
+          elemB%e2 => elemA
+
+          !Revers the normal to point inside the domain
+          elemB%normal = - elemB%normal
+
+        END IF
+
+      END IF
+
+      !Check direction 3
+      IF (.NOT. ASSOCIATED(elemA%e3)) THEN
+        IF     (elemA%n3%n == elemB%n1%n .AND. &
+                elemA%n4%n == elemB%n2%n) THEN
+          elemA%e3 => elemB
+          elemB%e1 => elemA
+
+        ELSEIF (elemA%n3%n == elemB%n2%n .AND. &
+                elemA%n4%n == elemB%n1%n) THEN
+          elemA%e3 => elemB
+          elemB%e2 => elemA
+
+          !Revers the normal to point inside the domain
+          elemB%normal = - elemB%normal
+
+        END IF
+
+      END IF
+
+      !Check direction 4
+      IF (.NOT. ASSOCIATED(elemA%e4)) THEN
+        IF     (elemA%n4%n == elemB%n1%n .AND. &
+                elemA%n1%n == elemB%n2%n) THEN
+          elemA%e4 => elemB
+          elemB%e1 => elemA
+
+        ELSEIF (elemA%n4%n == elemB%n2%n .AND. &
+                elemA%n1%n == elemB%n1%n) THEN
+          elemA%e4 => elemB
+          elemB%e2 => elemA
+
+          !Revers the normal to point inside the domain
+          elemB%normal = - elemB%normal
+
+        END IF
+
+      END IF
+
+    END SUBROUTINE connectQuadEdge
+
+    SUBROUTINE connectTriaEdge(elemA, elemB)
+      IMPLICIT NONE
+
+      CLASS(meshVol2DCartTria), INTENT(inout), TARGET:: elemA
+      CLASS(meshEdge2DCart),    INTENT(inout), TARGET:: elemB
+
+      !Check direction 1
+      IF (.NOT. ASSOCIATED(elemA%e1)) THEN
+        IF     (elemA%n1%n == elemB%n1%n .AND. &
+                elemA%n2%n == elemB%n2%n) THEN
+          elemA%e1 => elemB
+          elemB%e1 => elemA
+
+        ELSEIF (elemA%n1%n == elemB%n2%n .AND. &
+                elemA%n2%n == elemB%n1%n) THEN
+          elemA%e1 => elemB
+          elemB%e2 => elemA
+
+          !Revers the normal to point inside the domain
+          elemB%normal = - elemB%normal
+          
+        END IF
+
+      END IF
+
+      !Check direction 2
+      IF (.NOT. ASSOCIATED(elemA%e2)) THEN
+        IF     (elemA%n2%n == elemB%n1%n .AND. &
+                elemA%n3%n == elemB%n2%n) THEN
+          elemA%e2 => elemB
+          elemB%e1 => elemA
+
+        ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
+                elemA%n3%n == elemB%n1%n) THEN
+          elemA%e2 => elemB
+          elemB%e2 => elemA
+
+          !Revers the normal to point inside the domain
+          elemB%normal = - elemB%normal
+
+        END IF
+
+      END IF
+
+      !Check direction 3
+      IF (.NOT. ASSOCIATED(elemA%e3)) THEN
+        IF     (elemA%n3%n == elemB%n1%n .AND. &
+                elemA%n1%n == elemB%n2%n) THEN
+          elemA%e3 => elemB
+          elemB%e1 => elemA
+
+        ELSEIF (elemA%n3%n == elemB%n2%n .AND. &
+                elemA%n1%n == elemB%n1%n) THEN
+          elemA%e3 => elemB
+          elemB%e2 => elemA
+
+          !Revers the normal to point inside the domain
+          elemB%normal = - elemB%normal
+
+        END IF
+
+      END IF
+
+    END SUBROUTINE connectTriaEdge
+
 END MODULE moduleMesh2DCart

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

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

src/modules/mesh/2DCyl/makefile

@@ -1,8 +1,5 @@
-all : moduleMesh2DCyl.o moduleMesh2DCylRead.o
+all : moduleMesh2DCyl.o
 
 moduleMesh2DCyl.o: moduleMesh2DCyl.f90
 	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
 
-moduleMesh2DCylRead.o: moduleMesh2DCyl.o moduleMesh2DCylRead.f90
-	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
-

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

@@ -77,7 +77,7 @@ MODULE moduleMesh2DCyl
     !Connectivity to nodes
     CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL(), n3 => NULL(), n4 => NULL()
     !Connectivity to adjacent elements
-    CLASS(*), POINTER:: e1 => NULL(), e2 => NULL(), e3 => NULL(), e4 => NULL()
+    CLASS(meshElement), POINTER:: e1 => NULL(), e2 => NULL(), e3 => NULL(), e4 => NULL()
     REAL(8):: arNodes(1:4) = 0.D0
 
     CONTAINS
@@ -108,7 +108,7 @@ MODULE moduleMesh2DCyl
     !Connectivity to nodes
     CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL(), n3 => NULL()
     !Connectivity to adjacent elements
-    CLASS(*), POINTER:: e1 => NULL(), e2 => NULL(), e3 => NULL()
+    CLASS(meshElement), POINTER:: e1 => NULL(), e2 => NULL(), e3 => NULL()
     REAL(8):: arNodes(1:3) = 0.D0
 
     CONTAINS
@@ -271,20 +271,21 @@ MODULE moduleMesh2DCyl
     !VOLUME FUNCTIONS
     !QUAD FUNCTIONS
     !Inits quadrilateral element
-    SUBROUTINE initVolQuad2DCyl(self, n, p)
+    SUBROUTINE initVolQuad2DCyl(self, n, p, nodes)
       USE moduleRefParam
       IMPLICIT NONE
 
       CLASS(meshVol2DCylQuad), INTENT(out):: self
       INTEGER, INTENT(in):: n
       INTEGER, INTENT(in):: p(:)
+      TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:)
       REAL(8), DIMENSION(1:3):: r1, r2, r3, r4
 
       self%n  = n
-      self%n1 => mesh%nodes(p(1))%obj
+      self%n1 => nodes(p(1))%obj
-      self%n2 => mesh%nodes(p(2))%obj
+      self%n2 => nodes(p(2))%obj
-      self%n3 => mesh%nodes(p(3))%obj
+      self%n3 => nodes(p(3))%obj
-      self%n4 => mesh%nodes(p(4))%obj
+      self%n4 => nodes(p(4))%obj
       !Get element coordinates
       r1 = self%n1%getCoordinates()
       r2 = self%n2%getCoordinates()
@@ -427,18 +428,19 @@ MODULE moduleMesh2DCyl
     END FUNCTION randposVolQuad
 
     !Computes element local stiffness matrix
-    PURE FUNCTION elemKQuad(self) RESULT(ke)
+    PURE FUNCTION elemKQuad(self) RESULT(localK)
       USE moduleConstParam, ONLY: PI2
       IMPLICIT NONE
 
       CLASS(meshVol2DCylQuad), INTENT(in):: self
+      REAL(8), ALLOCATABLE:: localK(:,:)
       REAL(8):: r, xi(1:3)
       REAL(8):: fPsi(1:4), dPsi(1:2,1:4)
-      REAL(8):: ke(1:4,1:4)
       REAL(8):: invJ(1:2,1:2), detJ
       INTEGER:: l, m
 
-      ke=0.D0
+      ALLOCATE(localK(1:4, 1:4))
+      localK=0.D0
       xi=0.D0
       !Start 2D Gauss Quad Integral
       DO l=1, 3
@@ -451,13 +453,13 @@ MODULE moduleMesh2DCyl
           detJ      = self%detJac(xi,dPsi)
           invJ      = self%invJac(xi,dPsi)
           r         = DOT_PRODUCT(fPsi,self%r)
-          ke        = ke + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)), &
+          localK    = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)), &
                                             MATMUL(invJ,dPsi))* &
                            r*wQuad(l)*wQuad(m)/detJ
 
         END DO
       END DO
-      ke = ke*PI2
+      localK = localK*PI2
 
     END FUNCTION elemKQuad
 
@@ -531,22 +533,22 @@ MODULE moduleMesh2DCyl
       w_p = self%weight(part%xi)
       tensorS = outerProduct(part%v, part%v)
 
-      vertex => self%n1%output(part%sp)
+      vertex => self%n1%output(part%species%n)
       vertex%den          = vertex%den          + part%weight*w_p(1)
       vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(1)*part%v(:)
       vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(1)*tensorS
 
-      vertex => self%n2%output(part%sp)
+      vertex => self%n2%output(part%species%n)
       vertex%den          = vertex%den          + part%weight*w_p(2)
       vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(2)*part%v(:)
       vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(2)*tensorS
 
-      vertex => self%n3%output(part%sp)
+      vertex => self%n3%output(part%species%n)
       vertex%den          = vertex%den          + part%weight*w_p(3)
       vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(3)*part%v(:)
       vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(3)*tensorS
 
-      vertex => self%n4%output(part%sp)
+      vertex => self%n4%output(part%species%n)
       vertex%den          = vertex%den          + part%weight*w_p(4)
       vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(4)*part%v(:)
       vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(4)*tensorS
@@ -628,7 +630,7 @@ MODULE moduleMesh2DCyl
 
       CLASS(meshVol2DCylQuad), INTENT(in):: self
       REAL(8), INTENT(in):: xi(1:3)
-      CLASS(*), POINTER, INTENT(out):: nextElement
+      CLASS(meshElement), POINTER, INTENT(out):: nextElement
       REAL(8):: xiArray(1:4)
       INTEGER:: nextInt
 
@@ -651,22 +653,23 @@ MODULE moduleMesh2DCyl
 
     !TRIA ELEMENT
     !Init tria element
-    SUBROUTINE initVolTria2DCyl(self, n, p)
+    SUBROUTINE initVolTria2DCyl(self, n, p, nodes)
       USE moduleRefParam
       IMPLICIT NONE
 
       CLASS(meshVol2DCylTria), INTENT(out):: self
       INTEGER, INTENT(in):: n
       INTEGER, INTENT(in):: p(:)
+      TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:)
       REAL(8), DIMENSION(1:3):: r1, r2, r3
 
       !Assign node index
       self%n = n
 
       !Assign nodes to element
-      self%n1 => mesh%nodes(p(1))%obj
+      self%n1 => nodes(p(1))%obj
-      self%n2 => mesh%nodes(p(2))%obj
+      self%n2 => nodes(p(2))%obj
-      self%n3 => mesh%nodes(p(3))%obj
+      self%n3 => nodes(p(3))%obj
       !Get element coordinates
       r1 = self%n1%getCoordinates()
       r2 = self%n2%getCoordinates()
@@ -800,32 +803,33 @@ MODULE moduleMesh2DCyl
     END SUBROUTINE partialDerTria
 
     !Computes element local stiffness matrix
-    PURE FUNCTION elemKTria(self) RESULT(ke)
+    PURE FUNCTION elemKTria(self) RESULT(localK)
       USE moduleConstParam, ONLY: PI2
       IMPLICIT NONE
 
       CLASS(meshVol2DCylTria), INTENT(in):: self
+      REAL(8), ALLOCATABLE:: localK(:,:)
       REAL(8):: r, xi(1:3)
       REAL(8):: fPsi(1:3), dPsi(1:2,1:3)
-      REAL(8):: ke(1:3,1:3)
       REAL(8):: invJ(1:2,1:2), detJ
       INTEGER:: l
 
-      ke=0.D0
+      ALLOCATE(localK(1:4, 1:4))
+      localK=0.D0
       xi=0.D0
       !Start 2D Gauss Quad Integral
       DO l=1, 4
-        xi(1) = xi1Tria(l)
+        xi(1)  = xi1Tria(l)
-        xi(2) = xi2Tria(l)
+        xi(2)  = xi2Tria(l)
-        dPsi  = self%dPsi(xi)
+        dPsi   = self%dPsi(xi)
-        detJ  = self%detJac(xi,dPsi)
+        detJ   = self%detJac(xi,dPsi)
-        invJ  = self%invJac(xi,dPsi)
+        invJ   = self%invJac(xi,dPsi)
-        fPsi  = self%fPsi(xi)
+        fPsi   = self%fPsi(xi)
-        r     = DOT_PRODUCT(fPsi,self%r)
+        r      = DOT_PRODUCT(fPsi,self%r)
-        ke    = ke + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*r*wTria(l)/detJ
+        localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*r*wTria(l)/detJ
 
       END DO
-      ke = ke*PI2
+      localK = localK*PI2
 
     END FUNCTION elemKTria
 
@@ -898,17 +902,17 @@ MODULE moduleMesh2DCyl
       w_p = self%weight(part%xi)
       tensorS = outerProduct(part%v, part%v)
 
-      vertex => self%n1%output(part%sp)
+      vertex => self%n1%output(part%species%n)
       vertex%den          = vertex%den          + part%weight*w_p(1)
       vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(1)*part%v(:)
       vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(1)*tensorS
 
-      vertex => self%n2%output(part%sp)
+      vertex => self%n2%output(part%species%n)
       vertex%den          = vertex%den          + part%weight*w_p(2)
       vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(2)*part%v(:)
       vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(2)*tensorS
 
-      vertex => self%n3%output(part%sp)
+      vertex => self%n3%output(part%species%n)
       vertex%den          = vertex%den          + part%weight*w_p(3)
       vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(3)*part%v(:)
       vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(3)*tensorS
@@ -979,7 +983,7 @@ MODULE moduleMesh2DCyl
 
       CLASS(meshVol2DCylTria), INTENT(in):: self
       REAL(8), INTENT(in):: xi(1:3)
-      CLASS(*), POINTER, INTENT(out):: nextElement
+      CLASS(meshElement), POINTER, INTENT(out):: nextElement
       REAL(8):: xiArray(1:3)
       INTEGER:: nextInt
 
@@ -1047,4 +1051,453 @@ MODULE moduleMesh2DCyl
 
     END FUNCTION invJ2DCyl
 
+    SUBROUTINE connectMesh2DCyl(self)
+      IMPLICIT NONE
+
+      CLASS(meshGeneric), INTENT(inout):: self
+      INTEGER:: e, et
+
+      DO e = 1, self%numVols
+        !Connect Vol-Vol
+        DO et = 1, self%numVols
+          IF (e /= et) THEN
+            CALL connectVolVol(self%vols(e)%obj, self%vols(et)%obj)
+
+          END IF
+
+        END DO
+
+        SELECT TYPE(self)
+        TYPE IS(meshParticles)
+          !Connect Vol-Edge
+          DO et = 1, self%numEdges
+            CALL connectVolEdge(self%vols(e)%obj, self%edges(et)%obj)
+
+          END DO
+
+        END SELECT
+
+      END DO
+
+    END SUBROUTINE connectMesh2DCyl
+
+    !Selects type of elements to build connection
+    SUBROUTINE connectVolVol(elemA, elemB)
+      IMPLICIT NONE
+
+      CLASS(meshVol), INTENT(inout):: elemA
+      CLASS(meshVol), INTENT(inout):: elemB
+
+      SELECT TYPE(elemA)
+      TYPE IS(meshVol2DCylQuad)
+        !Element A is a quadrilateral
+        SELECT TYPE(elemB)
+        TYPE IS(meshVol2DCylQuad)
+          !Element B is a quadrilateral
+          CALL connectQuadQuad(elemA, elemB)
+
+        TYPE IS(meshVol2DCylTria)
+          !Element B is a triangle
+          CALL connectQuadTria(elemA, elemB)
+
+        END SELECT
+
+      TYPE IS(meshVol2DCylTria)
+        !Element A is a Triangle
+        SELECT TYPE(elemB)
+        TYPE IS(meshVol2DCylQuad)
+          !Element B is a quadrilateral
+          CALL connectQuadTria(elemB, elemA)
+
+        TYPE IS(meshVol2DCylTria)
+          !Element B is a triangle
+          CALL connectTriaTria(elemA, elemB)
+
+        END SELECT
+
+      END SELECT
+
+    END SUBROUTINE connectVolVol
+
+    SUBROUTINE connectVolEdge(elemA, elemB)
+      IMPLICIT NONE
+
+      CLASS(meshVol),  INTENT(inout):: elemA
+      CLASS(meshEdge), INTENT(inout):: elemB
+
+      SELECT TYPE(elemB)
+      CLASS IS(meshEdge2DCyl)
+        SELECT TYPE(elemA)
+        TYPE IS(meshVol2DCylQuad)
+          !Element A is a quadrilateral
+          CALL connectQuadEdge(elemA, elemB)
+
+        TYPE IS(meshVol2DCylTria)
+          !Element A is a triangle
+          CALL connectTriaEdge(elemA, elemB)
+
+        END SELECT
+
+      END SELECT
+
+    END SUBROUTINE connectVolEdge
+
+    SUBROUTINE connectQuadQuad(elemA, elemB)
+      IMPLICIT NONE
+
+      CLASS(meshVol2DCylQuad), INTENT(inout), TARGET:: elemA
+      CLASS(meshVol2DCylQuad), INTENT(inout), TARGET:: elemB
+
+      !Check direction 1
+      IF (.NOT. ASSOCIATED(elemA%e1) .AND. &
+          elemA%n1%n == elemB%n4%n   .AND. &
+          elemA%n2%n == elemB%n3%n) THEN
+        elemA%e1 => elemB
+        elemB%e3 => elemA
+
+      END IF
+
+      !Check direction 2
+      IF (.NOT. ASSOCIATED(elemA%e2) .AND. &
+          elemA%n2%n == elemB%n1%n   .AND. &
+          elemA%n3%n == elemB%n4%n) THEN
+        elemA%e2 => elemB
+        elemB%e4 => elemA
+
+      END IF
+
+      !Check direction 3
+      IF (.NOT. ASSOCIATED(elemA%e3) .AND. &
+          elemA%n3%n == elemB%n2%n   .AND. &
+          elemA%n4%n == elemB%n1%n) THEN
+        elemA%e3 => elemB
+        elemB%e1 => elemA
+
+      END IF
+
+      !Check direction 4
+      IF (.NOT. ASSOCIATED(elemA%e4) .AND. &
+          elemA%n4%n == elemB%n3%n   .AND. &
+          elemA%n1%n == elemB%n2%n) THEN
+        elemA%e4 => elemB
+        elemB%e2 => elemA
+
+      END IF
+
+    END SUBROUTINE connectQuadQuad
+
+    SUBROUTINE connectQuadTria(elemA, elemB)
+      IMPLICIT NONE
+
+      CLASS(meshVol2DCylQuad), INTENT(inout), TARGET:: elemA
+      CLASS(meshVol2DCylTria), INTENT(inout), TARGET:: elemB
+
+      !Check direction 1
+      IF (.NOT. ASSOCIATED(elemA%e1)) THEN
+        IF     (elemA%n1%n == elemB%n1%n .AND. &
+                elemA%n2%n == elemB%n3%n) THEN
+          elemA%e1 => elemB
+          elemB%e3 => elemA
+
+        ELSEIF (elemA%n1%n == elemB%n3%n .AND. &
+                elemA%n2%n == elemB%n2%n) THEN
+          elemA%e1 => elemB
+          elemB%e2 => elemA
+
+        ELSEIF (elemA%n1%n == elemB%n2%n .AND. &
+                elemA%n2%n == elemB%n1%n) THEN
+          elemA%e1 => elemB
+          elemB%e1 => elemA
+
+        END IF
+
+      END IF
+
+      !Check direction 2
+      IF (.NOT. ASSOCIATED(elemA%e2)) THEN
+        IF     (elemA%n2%n == elemB%n1%n .AND. &
+                elemA%n3%n == elemB%n3%n) THEN
+          elemA%e2 => elemB
+          elemB%e3 => elemA
+
+        ELSEIF (elemA%n2%n == elemB%n3%n .AND. &
+                elemA%n3%n == elemB%n2%n) THEN
+          elemA%e2 => elemB
+          elemB%e2 => elemA
+
+        ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
+                elemA%n3%n == elemB%n1%n) THEN
+          elemA%e2 => elemB
+          elemB%e1 => elemA
+
+        END IF
+
+      END IF
+
+      !Check direction 3
+      IF (.NOT. ASSOCIATED(elemA%e3)) THEN
+        IF     (elemA%n3%n == elemB%n1%n .AND. &
+                elemA%n4%n == elemB%n3%n) THEN
+          elemA%e3 => elemB
+          elemB%e3 => elemA
+
+        ELSEIF (elemA%n3%n == elemB%n3%n .AND. &
+                elemA%n4%n == elemB%n2%n) THEN
+          elemA%e3 => elemB
+          elemB%e2 => elemA
+
+        ELSEIF (elemA%n3%n == elemB%n2%n .AND. &
+                elemA%n4%n == elemB%n1%n) THEN
+          elemA%e3 => elemB
+          elemB%e1 => elemA
+
+        END IF
+
+      END IF
+
+      !Check direction 4
+      IF (.NOT. ASSOCIATED(elemA%e4)) THEN
+        IF     (elemA%n4%n == elemB%n1%n .AND. &
+                elemA%n1%n == elemB%n3%n) THEN
+          elemA%e4 => elemB
+          elemB%e3 => elemA
+
+        ELSEIF (elemA%n4%n == elemB%n3%n .AND. &
+                elemA%n1%n == elemB%n2%n) THEN
+          elemA%e4 => elemB
+          elemB%e2 => elemA
+
+        ELSEIF (elemA%n4%n == elemB%n2%n .AND. &
+                elemA%n1%n == elemB%n1%n) THEN
+          elemA%e4 => elemB
+          elemB%e1 => elemA
+
+        END IF
+
+      END IF
+
+    END SUBROUTINE connectQuadTria
+
+    SUBROUTINE connectTriaTria(elemA, elemB)
+      IMPLICIT NONE
+
+      CLASS(meshVol2DCylTria), INTENT(inout), TARGET:: elemA
+      CLASS(meshVol2DCylTria), INTENT(inout), TARGET:: elemB
+
+      !Check direction 1
+      IF (.NOT. ASSOCIATED(elemA%e1)) THEN
+        IF     (elemA%n1%n == elemB%n1%n .AND. &
+                elemA%n2%n == elemB%n3%n) THEN
+          elemA%e1 => elemB
+          elemB%e3 => elemA
+
+        ELSEIF (elemA%n1%n == elemB%n2%n .AND. &
+                elemA%n2%n == elemB%n1%n) THEN
+          elemA%e1 => elemB
+          elemB%e1 => elemA
+
+        ELSEIF (elemA%n1%n == elemB%n3%n .AND. &
+                elemA%n2%n == elemB%n2%n) THEN
+          elemA%e1 => elemB
+          elemB%e2 => elemA
+
+        END IF
+
+      END IF
+
+      !Check direction 2
+      IF (.NOT. ASSOCIATED(elemA%e2)) THEN
+        IF     (elemA%n2%n == elemB%n1%n .AND. &
+                elemA%n3%n == elemB%n3%n) THEN
+          elemA%e2 => elemB
+          elemB%e3 => elemA
+
+        ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
+                elemA%n3%n == elemB%n1%n) THEN
+          elemA%e2 => elemB
+          elemB%e1 => elemA
+
+        ELSEIF (elemA%n2%n == elemB%n3%n .AND. &
+                elemA%n3%n == elemB%n2%n) THEN
+          elemA%e2 => elemB
+          elemB%e2 => elemA
+
+        END IF
+
+
+      END IF
+
+      !Check direction 3
+      IF (.NOT. ASSOCIATED(elemA%e3)) THEN
+        IF     (elemA%n3%n == elemB%n1%n .AND. &
+                elemA%n1%n == elemB%n3%n) THEN
+          elemA%e3 => elemB
+          elemB%e3 => elemA
+
+        ELSEIF (elemA%n3%n == elemB%n2%n .AND. &
+                elemA%n1%n == elemB%n1%n) THEN
+          elemA%e3 => elemB
+          elemB%e1 => elemA
+
+        ELSEIF (elemA%n3%n == elemB%n3%n .AND. &
+                elemA%n1%n == elemB%n2%n) THEN
+          elemA%e3 => elemB
+          elemB%e2 => elemA
+
+        END IF
+
+
+      END IF
+
+    END SUBROUTINE connectTriaTria
+
+    SUBROUTINE connectQuadEdge(elemA, elemB)
+      IMPLICIT NONE
+
+      CLASS(meshVol2DCylQuad), INTENT(inout), TARGET:: elemA
+      CLASS(meshEdge2DCyl),    INTENT(inout), TARGET:: elemB
+
+      !Check direction 1
+      IF (.NOT. ASSOCIATED(elemA%e1)) THEN
+        IF     (elemA%n1%n == elemB%n1%n .AND. &
+                elemA%n2%n == elemB%n2%n) THEN
+          elemA%e1 => elemB
+          elemB%e1 => elemA
+
+        ELSEIF (elemA%n1%n == elemB%n2%n .AND. &
+                elemA%n2%n == elemB%n1%n) THEN
+          elemA%e1 => elemB
+          elemB%e2 => elemA
+
+          !Revers the normal to point inside the domain
+          elemB%normal = - elemB%normal
+
+        END IF
+
+      END IF
+
+      !Check direction 2
+      IF (.NOT. ASSOCIATED(elemA%e2)) THEN
+        IF     (elemA%n2%n == elemB%n1%n .AND. &
+                elemA%n3%n == elemB%n2%n) THEN
+          elemA%e2 => elemB
+          elemB%e1 => elemA
+
+        ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
+                elemA%n3%n == elemB%n1%n) THEN
+          elemA%e2 => elemB
+          elemB%e2 => elemA
+
+          !Revers the normal to point inside the domain
+          elemB%normal = - elemB%normal
+
+        END IF
+
+      END IF
+
+      !Check direction 3
+      IF (.NOT. ASSOCIATED(elemA%e3)) THEN
+        IF     (elemA%n3%n == elemB%n1%n .AND. &
+                elemA%n4%n == elemB%n2%n) THEN
+          elemA%e3 => elemB
+          elemB%e1 => elemA
+
+        ELSEIF (elemA%n3%n == elemB%n2%n .AND. &
+                elemA%n4%n == elemB%n1%n) THEN
+          elemA%e3 => elemB
+          elemB%e2 => elemA
+
+          !Revers the normal to point inside the domain
+          elemB%normal = - elemB%normal
+
+        END IF
+
+      END IF
+
+      !Check direction 4
+      IF (.NOT. ASSOCIATED(elemA%e4)) THEN
+        IF     (elemA%n4%n == elemB%n1%n .AND. &
+                elemA%n1%n == elemB%n2%n) THEN
+          elemA%e4 => elemB
+          elemB%e1 => elemA
+
+        ELSEIF (elemA%n4%n == elemB%n2%n .AND. &
+                elemA%n1%n == elemB%n1%n) THEN
+          elemA%e4 => elemB
+          elemB%e2 => elemA
+
+          !Revers the normal to point inside the domain
+          elemB%normal = - elemB%normal
+
+        END IF
+
+      END IF
+
+    END SUBROUTINE connectQuadEdge
+
+    SUBROUTINE connectTriaEdge(elemA, elemB)
+      IMPLICIT NONE
+
+      CLASS(meshVol2DCylTria), INTENT(inout), TARGET:: elemA
+      CLASS(meshEdge2DCyl),    INTENT(inout), TARGET:: elemB
+
+      !Check direction 1
+      IF (.NOT. ASSOCIATED(elemA%e1)) THEN
+        IF     (elemA%n1%n == elemB%n1%n .AND. &
+                elemA%n2%n == elemB%n2%n) THEN
+          elemA%e1 => elemB
+          elemB%e1 => elemA
+
+        ELSEIF (elemA%n1%n == elemB%n2%n .AND. &
+                elemA%n2%n == elemB%n1%n) THEN
+          elemA%e1 => elemB
+          elemB%e2 => elemA
+
+          !Revers the normal to point inside the domain
+          elemB%normal = - elemB%normal
+
+        END IF
+
+      END IF
+
+      !Check direction 2
+      IF (.NOT. ASSOCIATED(elemA%e2)) THEN
+        IF     (elemA%n2%n == elemB%n1%n .AND. &
+                elemA%n3%n == elemB%n2%n) THEN
+          elemA%e2 => elemB
+          elemB%e1 => elemA
+
+        ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
+                elemA%n3%n == elemB%n1%n) THEN
+          elemA%e2 => elemB
+          elemB%e2 => elemA
+
+          !Revers the normal to point inside the domain
+          elemB%normal = - elemB%normal
+
+        END IF
+
+      END IF
+
+      !Check direction 3
+      IF (.NOT. ASSOCIATED(elemA%e3)) THEN
+        IF     (elemA%n3%n == elemB%n1%n .AND. &
+                elemA%n1%n == elemB%n2%n) THEN
+          elemA%e3 => elemB
+          elemB%e1 => elemA
+
+        ELSEIF (elemA%n3%n == elemB%n2%n .AND. &
+                elemA%n1%n == elemB%n1%n) THEN
+          elemA%e3 => elemB
+          elemB%e2 => elemA
+
+          !Revers the normal to point inside the domain
+          elemB%normal = - elemB%normal
+
+        END IF
+
+      END IF
+
+    END SUBROUTINE connectTriaEdge
+
 END MODULE moduleMesh2DCyl

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

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

src/modules/mesh/3DCart/makefile

@@ -1,8 +1,5 @@
-all : moduleMesh3DCart.o moduleMesh3DCartRead.o
+all : moduleMesh3DCart.o
 
 moduleMesh3DCart.o: moduleMesh3DCart.f90
 	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
 
-moduleMesh3DCartRead.o: moduleMesh3DCart.o moduleMesh3DCartRead.f90
-	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
-

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

@@ -71,7 +71,7 @@ MODULE moduleMesh3DCart
     !Connectivity to nodes
     CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL(), n3 => NULL(), n4 => NULL()
     !Connectivity to adjacent elements
-    CLASS(*), POINTER:: e1 => NULL(), e2 => NULL(), e3 => NULL(), e4 => NULL()
+    CLASS(meshElement), POINTER:: e1 => NULL(), e2 => NULL(), e3 => NULL(), e4 => NULL()
     CONTAINS
       PROCEDURE, PASS::   init        => initVolTetra3DCart
       PROCEDURE, PASS::   randPos     => randPosVolTetra
@@ -248,21 +248,22 @@ MODULE moduleMesh3DCart
     !VOLUME FUNCTIONS
     !TETRA FUNCTIONS
     !Inits tetrahedron element
-    SUBROUTINE initVolTetra3DCart(self, n, p)
+    SUBROUTINE initVolTetra3DCart(self, n, p, nodes)
       USE moduleRefParam
       IMPLICIT NONE
 
       CLASS(meshVol3DCartTetra), INTENT(out):: self
       INTEGER, INTENT(in):: n
       INTEGER, INTENT(in):: p(:)
+      TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:)
       REAL(8), DIMENSION(1:3):: r1, r2, r3, r4 !Positions of each node
       REAL(8):: volNodes(1:4) !Volume of each node
 
       self%n  = n
-      self%n1 => mesh%nodes(p(1))%obj
+      self%n1 => nodes(p(1))%obj
-      self%n2 => mesh%nodes(p(2))%obj
+      self%n2 => nodes(p(2))%obj
-      self%n3 => mesh%nodes(p(3))%obj
+      self%n3 => nodes(p(3))%obj
-      self%n4 => mesh%nodes(p(4))%obj
+      self%n4 => nodes(p(4))%obj
       !Get element coordinates
       r1 = self%n1%getCoordinates()
       r2 = self%n2%getCoordinates()
@@ -417,23 +418,25 @@ MODULE moduleMesh3DCart
 
     END SUBROUTINE partialDerTetra
 
-    PURE FUNCTION elemKTetra(self) RESULT(ke)
+    PURE FUNCTION elemKTetra(self) RESULT(localK)
       IMPLICIT NONE
 
       CLASS(meshVol3DCartTetra), INTENT(in):: self
+      REAL(8), ALLOCATABLE:: localK(:,:)
       REAL(8):: xii(1:3)
       REAL(8):: fPsi(1:4), dPsi(1:3, 1:4)
-      REAL(8):: ke(1:4,1:4)
       REAL(8):: invJ(1:3,1:3), detJ
 
+      ALLOCATE(localK(1:4,1:4))
+      localK = 0.D0
+      xii    = 0.D0
       !TODO: One point Gauss integral. Upgrade when possible
-      ke = 0.D0
+      xii    = (/ 0.25D0, 0.25D0, 0.25D0 /)
-      xii = (/ 0.25D0, 0.25D0, 0.25D0 /)
+      dPsi   = self%dPsi(xii)
-      dPsi = self%dPsi(xii)
+      detJ   = self%detJac(xii, dPsi)
-      detJ = self%detJac(xii, dPsi)
+      invJ   = self%invJac(xii, dPsi)
-      invJ = self%invJac(xii, dPsi)
+      fPsi   = self%fPsi(xii)
-      fPsi = self%fPsi(xii)
+      localK = MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*1.D0/detJ
-      ke   = ke + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*1.D0/detJ
 
     END FUNCTION elemKTetra
 
@@ -456,7 +459,7 @@ MODULE moduleMesh3DCart
       detJ = self%detJac(xii, dPsi)
       fPsi = self%fPsi(xii)
       f    = DOT_PRODUCT(fPsi, source)
-      localF = localF + f*fPsi*1.D0*detJ
+      localF = f*fPsi*1.D0*detJ
 
     END FUNCTION elemFTetra
 
@@ -496,22 +499,22 @@ MODULE moduleMesh3DCart
       w_p = self%weight(part%xi)
       tensorS = outerProduct(part%v, part%v)
 
-      vertex => self%n1%output(part%sp)
+      vertex => self%n1%output(part%species%n)
       vertex%den          = vertex%den          + part%weight*w_p(1)
       vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(1)*part%v(:)
       vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(1)*tensorS
 
-      vertex => self%n2%output(part%sp)
+      vertex => self%n2%output(part%species%n)
       vertex%den          = vertex%den          + part%weight*w_p(2)
       vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(2)*part%v(:)
       vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(2)*tensorS
 
-      vertex => self%n3%output(part%sp)
+      vertex => self%n3%output(part%species%n)
       vertex%den          = vertex%den          + part%weight*w_p(3)
       vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(3)*part%v(:)
       vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(3)*tensorS
 
-      vertex => self%n4%output(part%sp)
+      vertex => self%n4%output(part%species%n)
       vertex%den          = vertex%den          + part%weight*w_p(4)
       vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(4)*part%v(:)
       vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(4)*tensorS
@@ -579,7 +582,7 @@ MODULE moduleMesh3DCart
 
       CLASS(meshVol3DCartTetra), INTENT(in):: self
       REAL(8), INTENT(in):: xi(1:3)
-      CLASS(*), POINTER, INTENT(out):: nextElement
+      CLASS(meshElement), POINTER, INTENT(out):: nextElement
       REAL(8):: xiArray(1:4)
       INTEGER:: nextInt
 
@@ -662,5 +665,373 @@ MODULE moduleMesh3DCart
 
     END FUNCTION invJ3DCart
 
+    !Selects type of elements to build connection
+    SUBROUTINE connectVolVol(elemA, elemB)
+      IMPLICIT NONE
+
+      CLASS(meshVol), INTENT(inout):: elemA
+      CLASS(meshVol), INTENT(inout):: elemB
+
+      SELECT TYPE(elemA)
+      TYPE IS(meshVol3DCartTetra)
+        !Element A is a tetrahedron
+        SELECT TYPE(elemB)
+        TYPE IS(meshVol3DCartTetra)
+          !Element B is a tetrahedron
+          CALL connectTetraTetra(elemA, elemB)
+
+        END SELECT
+
+      END SELECT
+
+    END SUBROUTINE connectVolVol
+
+    SUBROUTINE connectVolEdge(elemA, elemB)
+      IMPLICIT NONE
+
+      CLASS(meshVol),  INTENT(inout):: elemA
+      CLASS(meshEdge), INTENT(inout):: elemB
+
+      SELECT TYPE(elemB)
+      CLASS IS(meshEdge3DCartTria)
+        SELECT TYPE(elemA)
+        TYPE IS(meshVol3DCartTetra)
+          !Element A is a tetrahedron
+          CALL connectTetraEdge(elemA, elemB)
+
+        END SELECT
+
+      END SELECT
+
+    END SUBROUTINE connectVolEdge
+
+    SUBROUTINE connectMesh3DCart(self)
+      IMPLICIT NONE
+
+      CLASS(meshGeneric), INTENT(inout):: self
+      INTEGER:: e, et
+
+      DO e = 1, self%numVols
+        !Connect Vol-Vol
+        DO et = 1, self%numVols
+          IF (e /= et) THEN
+            CALL connectVolVol(self%vols(e)%obj, self%vols(et)%obj)
+
+          END IF
+
+        END DO
+
+        SELECT TYPE(self)
+        TYPE IS(meshParticles)
+          !Connect Vol-Edge
+          DO et = 1, self%numEdges
+            CALL connectVolEdge(self%vols(e)%obj, self%edges(et)%obj)
+
+          END DO
+
+        END SELECT
+
+      END DO
+
+    END SUBROUTINE connectMesh3DCart
+
+    !Checks if two sets of nodes are coincidend in any order
+    PURE FUNCTION coincidentNodes(nodesA, nodesB) RESULT(coincident)
+      IMPLICIT NONE
+
+      INTEGER, DIMENSION(1:3), INTENT(in):: nodesA, nodesB
+      LOGICAL:: coincident
+      INTEGER:: i
+
+      coincident = .FALSE.
+      DO i = 1, 3
+        IF (ANY(nodesA(i) == nodesB)) THEN
+          coincident = .TRUE.
+
+        ELSE
+          coincident = .FALSE.
+          EXIT
+
+        END IF
+
+      END DO
+
+    END FUNCTION coincidentNodes
+
+    SUBROUTINE connectTetraTetra(elemA, elemB)
+      IMPLICIT NONE
+
+      CLASS(meshVol3DCartTetra), INTENT(inout), TARGET:: elemA
+      CLASS(meshVol3DCartTetra), INTENT(inout), TARGET:: elemB
+
+      !Check surface 1
+      IF (.NOT. ASSOCIATED(elemA%e1)) THEN
+        IF     (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n3%n/), &
+                                (/elemB%n1%n, elemB%n2%n, elemB%n3%n/))) THEN
+
+          elemA%e1 => elemB
+          elemB%e1 => elemA
+
+        ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n3%n/), &
+                                (/elemB%n2%n, elemB%n3%n, elemB%n4%n/))) THEN
+
+          elemA%e1 => elemB
+          elemB%e2 => elemA
+
+        ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n3%n/), &
+                                (/elemB%n1%n, elemB%n2%n, elemB%n4%n/))) THEN
+
+          elemA%e1 => elemB
+          elemB%e3 => elemA
+
+        ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n3%n/), &
+                                (/elemB%n1%n, elemB%n3%n, elemB%n4%n/))) THEN
+
+          elemA%e1 => elemB
+          elemB%e4 => elemA
+
+        END IF
+
+      END IF
+
+      !Check surface 2
+      IF (.NOT. ASSOCIATED(elemA%e2)) THEN
+        IF     (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), &
+                                (/elemB%n1%n, elemB%n2%n, elemB%n3%n/))) THEN
+
+          elemA%e2 => elemB
+          elemB%e1 => elemA
+
+        ELSEIF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), &
+                                (/elemB%n2%n, elemB%n3%n, elemB%n4%n/))) THEN
+
+          elemA%e2 => elemB
+          elemB%e2 => elemA
+
+        ELSEIF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), &
+                                (/elemB%n1%n, elemB%n2%n, elemB%n4%n/))) THEN
+
+          elemA%e2 => elemB
+          elemB%e3 => elemA
+
+        ELSEIF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), &
+                                (/elemB%n1%n, elemB%n3%n, elemB%n4%n/))) THEN
+
+          elemA%e2 => elemB
+          elemB%e4 => elemA
+
+        END IF
+
+      END IF
+
+      !Check surface 3
+      IF (.NOT. ASSOCIATED(elemA%e3)) THEN
+        IF     (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), &
+                                (/elemB%n1%n, elemB%n2%n, elemB%n3%n/))) THEN
+
+          elemA%e3 => elemB
+          elemB%e1 => elemA
+
+        ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), &
+                                (/elemB%n2%n, elemB%n3%n, elemB%n4%n/))) THEN
+
+          elemA%e3 => elemB
+          elemB%e2 => elemA
+
+        ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), &
+                                (/elemB%n1%n, elemB%n2%n, elemB%n4%n/))) THEN
+
+          elemA%e3 => elemB
+          elemB%e3 => elemA
+
+        ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), &
+                                (/elemB%n1%n, elemB%n3%n, elemB%n4%n/))) THEN
+
+          elemA%e3 => elemB
+          elemB%e4 => elemA
+
+        END IF
+
+      END IF
+
+      !Check surface 4
+      IF (.NOT. ASSOCIATED(elemA%e4)) THEN
+        IF     (coincidentNodes((/elemA%n1%n, elemA%n3%n, elemA%n4%n/), &
+                                (/elemB%n1%n, elemB%n2%n, elemB%n3%n/))) THEN
+
+          elemA%e4 => elemB
+          elemB%e1 => elemA
+
+        ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n3%n, elemA%n4%n/), &
+                                (/elemB%n2%n, elemB%n3%n, elemB%n4%n/))) THEN
+
+          elemA%e4 => elemB
+          elemB%e2 => elemA
+
+        ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n3%n, elemA%n4%n/), &
+                                (/elemB%n1%n, elemB%n2%n, elemB%n4%n/))) THEN
+
+          elemA%e4 => elemB
+          elemB%e3 => elemA
+
+        ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n3%n, elemA%n4%n/), &
+                                (/elemB%n1%n, elemB%n3%n, elemB%n4%n/))) THEN
+
+          elemA%e4 => elemB
+          elemB%e4 => elemA
+
+        END IF
+
+      END IF
+
+    END SUBROUTINE connectTetraTetra
+
+    SUBROUTINE connectTetraEdge(elemA, elemB)
+      USE moduleMath
+      IMPLICIT NONE
+
+      CLASS(meshVol3DCartTetra), INTENT(inout), TARGET:: elemA
+      CLASS(meshEdge3DCartTria), INTENT(inout), TARGET:: elemB
+      INTEGER:: nodesEdge(1:3)
+      REAL(8), DIMENSION(1:3):: vec1, vec2
+      REAL(8):: normVol(1:3)
+
+      nodesEdge = (/ elemB%n1%n, elemB%n2%n, elemB%n3%n /)
+
+      !Check surface 1
+      IF (.NOT. ASSOCIATED(elemA%e1)) THEN
+        IF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elema%n3%n/), &
+                              nodesEdge)) THEN
+
+          vec1 = (/ elemA%x(2) - elemA%x(1), &
+                    elemA%y(2) - elemA%y(1), &
+                    elemA%z(2) - elemA%z(1) /)
+          vec2 = (/ elemA%x(3) - elemA%x(1), &
+                    elemA%y(3) - elemA%y(1), &
+                    elemA%z(3) - elemA%z(1) /)
+          normVol = crossProduct(vec1, vec2)
+          normVol = normalize(normVol)
+
+          IF (DOT_PRODUCT(elemB%normal, normVol) == -1.D0) THEN
+
+            elemA%e1 => elemB
+            elemB%e1 => elemA
+
+          ELSE
+
+            elemA%e1 => elemB
+            elemB%e2 => elemA
+
+            !Revers the normal to point inside the domain
+            elemB%normal = -elemB%normal
+
+          END IF
+
+        END IF
+
+      END IF
+
+      !Check surface 2
+      IF (.NOT. ASSOCIATED(elemA%e2)) THEN
+        IF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), &
+                              nodesEdge)) THEN
+
+          vec1 = (/ elemA%x(3) - elemA%x(2), &
+                    elemA%y(3) - elemA%y(2), &
+                    elemA%z(3) - elemA%z(2) /)
+          vec2 = (/ elemA%x(4) - elemA%x(2), &
+                    elemA%y(4) - elemA%y(2), &
+                    elemA%z(4) - elemA%z(2) /)
+          normVol = crossProduct(vec1, vec2)
+          normVol = normalize(normVol)
+
+          IF (DOT_PRODUCT(elemB%normal, normVol) == -1.D0) THEN
+
+            elemA%e2 => elemB
+            elemB%e1 => elemA
+
+          ELSE
+
+            elemA%e2 => elemB
+            elemB%e2 => elemA
+
+            !Revers the normal to point inside the domain
+            elemB%normal = -elemB%normal
+
+          END IF
+
+        END IF
+
+      END IF
+
+      !Check surface 3
+      IF (.NOT. ASSOCIATED(elemA%e3)) THEN
+        IF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elema%n4%n/), &
+                              nodesEdge)) THEN
+
+          vec1 = (/ elemA%x(2) - elemA%x(1), &
+                    elemA%y(2) - elemA%y(1), &
+                    elemA%z(2) - elemA%z(1) /)
+          vec2 = (/ elemA%x(4) - elemA%x(1), &
+                    elemA%y(4) - elemA%y(1), &
+                    elemA%z(4) - elemA%z(1) /)
+          normVol = crossProduct(vec1, vec2)
+          normVol = normalize(normVol)
+
+          IF (DOT_PRODUCT(elemB%normal, normVol) == -1.D0) THEN
+
+            elemA%e3 => elemB
+            elemB%e1 => elemA
+
+          ELSE
+
+            elemA%e3 => elemB
+            elemB%e2 => elemA
+
+            !Revers the normal to point inside the domain
+            elemB%normal = -elemB%normal
+
+          END IF
+
+        END IF
+
+      END IF
+
+      !Check surface 4
+      IF (.NOT. ASSOCIATED(elemA%e4)) THEN
+        IF (coincidentNodes((/elemA%n1%n, elemA%n3%n, elema%n4%n/), &
+                              nodesEdge)) THEN
+
+          vec1 = (/ elemA%x(3) - elemA%x(1), &
+                    elemA%y(3) - elemA%y(1), &
+                    elemA%z(3) - elemA%z(1) /)
+          vec2 = (/ elemA%x(4) - elemA%x(1), &
+                    elemA%y(4) - elemA%y(1), &
+                    elemA%z(4) - elemA%z(1) /)
+          normVol = crossProduct(vec1, vec2)
+          normVol = normalize(normVol)
+
+          IF (DOT_PRODUCT(elemB%normal, normVol) == -1.D0) THEN
+
+            elemA%e4 => elemB
+            elemB%e1 => elemA
+
+
+          ELSE
+
+            elemA%e4 => elemB
+            elemB%e2 => elemA
+
+            !Revers the normal to point inside the domain
+            elemB%normal = -elemB%normal
+
+          END IF
+
+        END IF
+
+      END IF
+
+    END SUBROUTINE connectTetraEdge
+
 END MODULE moduleMesh3DCart
 

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

@@ -1,544 +0,0 @@
-MODULE moduleMesh3DCartRead
-  USE moduleMesh
-  USE moduleMesh3DCart
-
-  TYPE, EXTENDS(meshGeneric):: mesh3DCartGeneric
-    CONTAINS
-      PROCEDURE, PASS:: init     => init3DCartMesh
-      PROCEDURE, PASS:: readMesh => readMesh3DCartGmsh
-
-  END TYPE
-
-  INTERFACE connected
-    MODULE PROCEDURE connectedVolVol, connectedVolEdge
-
-  END INTERFACE connected
-
-  CONTAINS
-    !Init mesh
-    SUBROUTINE init3DCartMesh(self, meshFormat)
-      USE moduleMesh
-      USE moduleErrors
-      IMPLICIT NONE
-
-      CLASS(mesh3DCartGeneric), INTENT(out):: self
-      CHARACTER(:), ALLOCATABLE, INTENT(in):: meshFormat
-
-      SELECT CASE(meshFormat)
-      CASE ("gmsh")
-        self%printOutput => printOutputGmsh
-        self%printColl   => printCollGmsh
-        self%printEM     => printEMGmsh
-
-      CASE DEFAULT
-        CALL criticalError("Mesh type " // meshFormat // " not supported.", "init3DCartMesh")
-
-      END SELECT
-
-    END SUBROUTINE init3DCartMesh
-
-    !Read mesh from gmsh file
-    SUBROUTINE readMesh3DCartGmsh(self, filename)
-      USE moduleBoundary
-      IMPLICIT NONE
-
-      CLASS(mesh3DCartGeneric), INTENT(inout):: self
-      CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
-      REAL(8):: x, y, z
-      INTEGER:: p(1:4)
-      INTEGER:: e = 0, et = 0, n = 0, eTemp = 0, elemType = 0, bt = 0
-      INTEGER:: totalNumElem
-      INTEGER:: boundaryType
-
-      !Read mesh
-      OPEN(10, FILE=TRIM(filename))
-      !Skip header
-      READ(10, *)
-      READ(10, *)
-      READ(10, *)
-      READ(10, *)
-      !Read number of nodes
-      READ(10, *) self%numNodes
-      !Allocate required matrices and vectors
-      ALLOCATE(self%nodes(1:self%numNodes))
-      ALLOCATE(self%K(1:self%numNodes, 1:self%numNodes))
-      ALLOCATE(self%IPIV(1:self%numNodes, 1:self%numNodes))
-      self%K = 0.D0
-      self%IPIV = 0
-
-      !Read node cartesian coordinates (x = x, y = y, z = z)
-      DO e = 1, self%numNodes
-        READ(10, *) n, x, y, z
-        ALLOCATE(meshNode3Dcart::self%nodes(n)%obj)
-        CALL self%nodes(n)%obj%init(n, (/x, y, z /))
-
-      END DO
-
-      !Skip comments
-      READ(10, *)
-      READ(10, *)
-
-      !Reads total number of elements
-      READ(10, *) totalNumElem
-      !conts edges and volume elements
-      self%numEdges = 0
-      DO e = 1, totalNumElem
-        READ(10, *) eTemp, elemType
-        IF (elemType == 2) THEN
-          self%numEdges = e
-
-        END IF
-
-      END DO
-
-      !Substract the number of edges to the total number of elements to obtain the number
-      !of volume elements
-      self%numVols = totalNumElem - self%numEdges
-
-      !Allocate required arrays
-      ALLOCATE(self%edges(1:self%numEdges))
-      ALLOCATE(self%vols(1:self%numVols))
-
-      !Go back to the beggining to read each specific element
-      DO e = 1, totalNumElem
-        BACKSPACE(10)
-
-      END DO
-
-      !Reads surfaces
-      DO e = 1, self%numEdges
-        READ(10, *) n, elemType
-        BACKSPACE(10)
-
-        SELECT CASE(elemType)
-        CASE(2)
-          !Triangular surface
-          READ(10, *) n, elemType, eTemp, boundaryType, eTemp, p(1:3)
-          bt = getBoundaryID(boundaryType)
-
-          ALLOCATE(meshEdge3DCartTria:: self%edges(e)%obj)
-
-          CALL self%edges(e)%obj%init(n, p(1:3), bt, boundaryType)
-
-        END SELECT
-
-      END DO
-
-      !Read and initialize volumes
-      DO e = 1, self%numVols
-        READ(10, *) n, elemType
-        BACKSPACE(10)
-
-        SELECT CASE(elemType)
-        CASE(4)
-          !Tetrahedron element
-          READ(10, *) n, elemType, eTemp, eTemp, eTemp, p(1:4)
-          ALLOCATE(meshVol3DCartTetra::  self%vols(e)%obj)
-          CALL self%vols(e)%obj%init(n - self%numEdges, p(1:4))
-
-        END SELECT
-
-      END DO
-
-      CLOSE(10)
-
-      !Build connectivy between elements
-      DO e = 1, self%numVols
-        !Connectivity between volumes
-        DO et = 1, self%numVols
-          IF (e /= et) THEN
-            CALL connected(self%vols(e)%obj, self%vols(et)%obj)
-
-          END IF
-
-        END DO
-
-        !Connectivity between vols and surfaces
-        DO et = 1, self%numEdges
-          CALL connected(self%vols(e)%obj, self%edges(et)%obj)
-
-        END DO
-
-        !Constructs the global K matrix
-        CALL constructGlobalK(self%K, self%vols(e)%obj)
-        
-      END DO
-
-    END SUBROUTINE readMesh3DCartGmsh
-
-    !Selects type of elements to build connection
-    SUBROUTINE connectedVolVol(elemA, elemB)
-      IMPLICIT NONE
-
-      CLASS(meshVol), INTENT(inout):: elemA
-      CLASS(meshVol), INTENT(inout):: elemB
-
-      SELECT TYPE(elemA)
-      TYPE IS(meshVol3DCartTetra)
-        !Element A is a tetrahedron
-        SELECT TYPE(elemB)
-        TYPE IS(meshVol3DCartTetra)
-          !Element B is a tetrahedron
-          CALL connectedTetraTetra(elemA, elemB)
-
-        END SELECT
-
-      END SELECT
-
-    END SUBROUTINE connectedVolVol
-
-
-    SUBROUTINE connectedVolEdge(elemA, elemB)
-      IMPLICIT NONE
-
-      CLASS(meshVol),  INTENT(inout):: elemA
-      CLASS(meshEdge), INTENT(inout):: elemB
-
-      SELECT TYPE(elemB)
-      CLASS IS(meshEdge3DCartTria)
-        SELECT TYPE(elemA)
-        TYPE IS(meshVol3DCartTetra)
-          !Element A is a tetrahedron
-          CALL connectedTetraEdge(elemA, elemB)
-
-        END SELECT
-
-      END SELECT
-
-    END SUBROUTINE connectedVolEdge
-
-    PURE FUNCTION coincidentNodes(nodesA, nodesB) RESULT(coincident)
-      IMPLICIT NONE
-
-      INTEGER, DIMENSION(1:3), INTENT(in):: nodesA, nodesB
-      LOGICAL:: coincident
-      INTEGER:: i
-
-      coincident = .FALSE.
-      DO i = 1, 3
-        IF (ANY(nodesA(i) == nodesB)) THEN
-          coincident = .TRUE.
-
-        ELSE
-          coincident = .FALSE.
-          EXIT
-
-        END IF
-
-      END DO
-
-    END FUNCTION coincidentNodes
-
-    SUBROUTINE connectedTetraTetra(elemA, elemB)
-      IMPLICIT NONE
-
-      CLASS(meshVol3DCartTetra), INTENT(inout), TARGET:: elemA
-      CLASS(meshVol3DCartTetra), INTENT(inout), TARGET:: elemB
-
-      !TODO: Try to find a much clear way to do this
-
-      !Check surface 1
-      IF (.NOT. ASSOCIATED(elemA%e1)) THEN
-        IF     (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n3%n/), &
-                                (/elemB%n1%n, elemB%n2%n, elemB%n3%n/))) THEN
-
-          elemA%e1 => elemB
-          elemB%e1 => elemA
-
-        ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n3%n/), &
-                                (/elemB%n2%n, elemB%n3%n, elemB%n4%n/))) THEN
-
-          elemA%e1 => elemB
-          elemB%e2 => elemA
-
-        ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n3%n/), &
-                                (/elemB%n1%n, elemB%n2%n, elemB%n4%n/))) THEN
-
-          elemA%e1 => elemB
-          elemB%e3 => elemA
-
-        ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n3%n/), &
-                                (/elemB%n1%n, elemB%n3%n, elemB%n4%n/))) THEN
-
-          elemA%e1 => elemB
-          elemB%e4 => elemA
-
-        END IF
-
-      END IF
-
-      !Check surface 2
-      IF (.NOT. ASSOCIATED(elemA%e2)) THEN
-        IF     (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), &
-                                (/elemB%n1%n, elemB%n2%n, elemB%n3%n/))) THEN
-
-          elemA%e2 => elemB
-          elemB%e1 => elemA
-
-        ELSEIF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), &
-                                (/elemB%n2%n, elemB%n3%n, elemB%n4%n/))) THEN
-
-          elemA%e2 => elemB
-          elemB%e2 => elemA
-
-        ELSEIF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), &
-                                (/elemB%n1%n, elemB%n2%n, elemB%n4%n/))) THEN
-
-          elemA%e2 => elemB
-          elemB%e3 => elemA
-
-        ELSEIF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), &
-                                (/elemB%n1%n, elemB%n3%n, elemB%n4%n/))) THEN
-
-          elemA%e2 => elemB
-          elemB%e4 => elemA
-
-        END IF
-
-      END IF
-
-      !Check surface 3
-      IF (.NOT. ASSOCIATED(elemA%e3)) THEN
-        IF     (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), &
-                                (/elemB%n1%n, elemB%n2%n, elemB%n3%n/))) THEN
-
-          elemA%e3 => elemB
-          elemB%e1 => elemA
-
-        ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), &
-                                (/elemB%n2%n, elemB%n3%n, elemB%n4%n/))) THEN
-
-          elemA%e3 => elemB
-          elemB%e2 => elemA
-
-        ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), &
-                                (/elemB%n1%n, elemB%n2%n, elemB%n4%n/))) THEN
-
-          elemA%e3 => elemB
-          elemB%e3 => elemA
-
-        ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), &
-                                (/elemB%n1%n, elemB%n3%n, elemB%n4%n/))) THEN
-
-          elemA%e3 => elemB
-          elemB%e4 => elemA
-
-        END IF
-
-      END IF
-
-      !Check surface 4
-      IF (.NOT. ASSOCIATED(elemA%e4)) THEN
-        IF     (coincidentNodes((/elemA%n1%n, elemA%n3%n, elemA%n4%n/), &
-                                (/elemB%n1%n, elemB%n2%n, elemB%n3%n/))) THEN
-
-          elemA%e4 => elemB
-          elemB%e1 => elemA
-
-        ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n3%n, elemA%n4%n/), &
-                                (/elemB%n2%n, elemB%n3%n, elemB%n4%n/))) THEN
-
-          elemA%e4 => elemB
-          elemB%e2 => elemA
-
-        ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n3%n, elemA%n4%n/), &
-                                (/elemB%n1%n, elemB%n2%n, elemB%n4%n/))) THEN
-
-          elemA%e4 => elemB
-          elemB%e3 => elemA
-
-        ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n3%n, elemA%n4%n/), &
-                                (/elemB%n1%n, elemB%n3%n, elemB%n4%n/))) THEN
-
-          elemA%e4 => elemB
-          elemB%e4 => elemA
-
-        END IF
-
-      END IF
-
-    END SUBROUTINE connectedTetraTetra
-
-    SUBROUTINE connectedTetraEdge(elemA, elemB)
-      USE moduleMath
-      IMPLICIT NONE
-
-      CLASS(meshVol3DCartTetra), INTENT(inout), TARGET:: elemA
-      CLASS(meshEdge3DCartTria), INTENT(inout), TARGET:: elemB
-      INTEGER:: nodesEdge(1:3)
-      REAL(8), DIMENSION(1:3):: vec1, vec2
-      REAL(8):: normVol(1:3)
-
-      nodesEdge = (/ elemB%n1%n, elemB%n2%n, elemB%n3%n /)
-
-      !Check surface 1
-      IF (.NOT. ASSOCIATED(elemA%e1)) THEN
-        IF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elema%n3%n/), &
-                              nodesEdge)) THEN
-
-          vec1 = (/ elemA%x(2) - elemA%x(1), &
-                    elemA%y(2) - elemA%y(1), &
-                    elemA%z(2) - elemA%z(1) /)
-          vec2 = (/ elemA%x(3) - elemA%x(1), &
-                    elemA%y(3) - elemA%y(1), &
-                    elemA%z(3) - elemA%z(1) /)
-          normVol = crossProduct(vec1, vec2)
-          normVol = normalize(normVol)
-
-          IF (DOT_PRODUCT(elemB%normal, normVol) == -1.D0) THEN
-
-            elemA%e1 => elemB
-            elemB%e1 => elemA
-
-          ELSE
-
-            elemA%e1 => elemB
-            elemB%e2 => elemA
-
-            !Revers the normal to point inside the domain
-            elemB%normal = -elemB%normal
-
-          END IF
-
-        END IF
-
-      END IF
-
-      !Check surface 2
-      IF (.NOT. ASSOCIATED(elemA%e2)) THEN
-        IF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), &
-                              nodesEdge)) THEN
-
-          vec1 = (/ elemA%x(3) - elemA%x(2), &
-                    elemA%y(3) - elemA%y(2), &
-                    elemA%z(3) - elemA%z(2) /)
-          vec2 = (/ elemA%x(4) - elemA%x(2), &
-                    elemA%y(4) - elemA%y(2), &
-                    elemA%z(4) - elemA%z(2) /)
-          normVol = crossProduct(vec1, vec2)
-          normVol = normalize(normVol)
-
-          IF (DOT_PRODUCT(elemB%normal, normVol) == -1.D0) THEN
-
-            elemA%e2 => elemB
-            elemB%e1 => elemA
-
-          ELSE
-
-            elemA%e2 => elemB
-            elemB%e2 => elemA
-
-            !Revers the normal to point inside the domain
-            elemB%normal = -elemB%normal
-
-          END IF
-
-        END IF
-
-      END IF
-
-      !Check surface 3
-      IF (.NOT. ASSOCIATED(elemA%e3)) THEN
-        IF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elema%n4%n/), &
-                              nodesEdge)) THEN
-
-          vec1 = (/ elemA%x(2) - elemA%x(1), &
-                    elemA%y(2) - elemA%y(1), &
-                    elemA%z(2) - elemA%z(1) /)
-          vec2 = (/ elemA%x(4) - elemA%x(1), &
-                    elemA%y(4) - elemA%y(1), &
-                    elemA%z(4) - elemA%z(1) /)
-          normVol = crossProduct(vec1, vec2)
-          normVol = normalize(normVol)
-
-          IF (DOT_PRODUCT(elemB%normal, normVol) == -1.D0) THEN
-
-            elemA%e3 => elemB
-            elemB%e1 => elemA
-
-          ELSE
-
-            elemA%e3 => elemB
-            elemB%e2 => elemA
-
-            !Revers the normal to point inside the domain
-            elemB%normal = -elemB%normal
-
-          END IF
-
-        END IF
-
-      END IF
-
-      !Check surface 4
-      IF (.NOT. ASSOCIATED(elemA%e4)) THEN
-        IF (coincidentNodes((/elemA%n1%n, elemA%n3%n, elema%n4%n/), &
-                              nodesEdge)) THEN
-
-          vec1 = (/ elemA%x(3) - elemA%x(1), &
-                    elemA%y(3) - elemA%y(1), &
-                    elemA%z(3) - elemA%z(1) /)
-          vec2 = (/ elemA%x(4) - elemA%x(1), &
-                    elemA%y(4) - elemA%y(1), &
-                    elemA%z(4) - elemA%z(1) /)
-          normVol = crossProduct(vec1, vec2)
-          normVol = normalize(normVol)
-
-          IF (DOT_PRODUCT(elemB%normal, normVol) == -1.D0) THEN
-
-            elemA%e4 => elemB
-            elemB%e1 => elemA
-
-
-          ELSE
-
-            elemA%e4 => elemB
-            elemB%e2 => elemA
-
-            !Revers the normal to point inside the domain
-            elemB%normal = -elemB%normal
-
-          END IF
-
-        END IF
-
-      END IF
-
-    END SUBROUTINE connectedTetraEdge
-
-    SUBROUTINE constructGlobalK(K, elem)
-      IMPLICIT NONE
-
-      REAL(8), INTENT(inout):: K(1:, 1:)
-      CLASS(meshVol), INTENT(in):: elem
-      REAL(8), ALLOCATABLE:: localK(:,:)
-      INTEGER:: nNodes, i, j
-      INTEGER, ALLOCATABLE:: n(:)
-
-      SELECT TYPE(elem)
-      TYPE IS(meshVol3DCartTetra)
-        nNodes = 4
-        ALLOCATE(localK(1:nNodes,1:nNodes))
-        localK = elem%elemK()
-        ALLOCATE(n(1:nNodes))
-        n = (/ elem%n1%n, elem%n2%n, &
-               elem%n3%n, elem%n4%n /)
-
-      CLASS DEFAULT
-        nNodes = 0
-        ALLOCATE(localK(1:1, 1:1))
-        localK = 0.D0
-        ALLOCATE(n(1:1))
-        n = 0
-
-      END SELECT
-      
-      DO i = 1, nNodes
-        DO j = 1, nNodes
-          K(n(i), n(j)) = K(n(i), n(j)) + localK(i, j)
-        END DO
-      END DO
-
-    END SUBROUTINE constructGlobalK
-
-END MODULE moduleMesh3DCartRead

src/modules/mesh/inout/gmsh2/makefile

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

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

@@ -0,0 +1,291 @@
+MODULE moduleMeshInputGmsh2
+
+  CONTAINS
+    !Inits a mesh to use Gmsh2 format
+    SUBROUTINE initGmsh2(self)
+      USE moduleMesh
+      USE moduleMeshOutputGmsh2
+      IMPLICIT NONE
+
+      CLASS(meshGeneric), INTENT(inout), TARGET:: self
+
+      IF (ASSOCIATED(meshForMCC, self)) self%printColl => printCollGmsh2
+      SELECT TYPE(self)
+      TYPE IS(meshParticles)
+        self%printOutput => printOutputGmsh2
+        self%printEM     => printEMGmsh2
+
+      END SELECT
+      self%readMesh    => readGmsh2
+
+    END SUBROUTINE initGmsh2
+
+    !Reads a Gmsh 2 format
+    SUBROUTINE readGmsh2(self, filename)
+      USE moduleMesh3DCart
+      USE moduleMesh2DCyl
+      USE moduleMesh2DCart
+      USE moduleMesh1DRad
+      USE moduleMesh1DCart
+      USE moduleBoundary
+      IMPLICIT NONE
+
+      CLASS(meshGeneric), INTENT(inout):: self
+      CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
+      REAL(8):: r(1:3) !3 generic coordinates
+      INTEGER, ALLOCATABLE:: p(:) !Array for nodes
+      INTEGER:: e = 0, n = 0, eTemp = 0, elemType = 0, bt = 0
+      INTEGER:: totalNumElem
+      INTEGER:: numEdges
+      INTEGER:: boundaryType
+
+      !Read mesh
+      OPEN(10, FILE=TRIM(filename))
+
+      !Skip header
+      READ(10, *)
+      READ(10, *)
+      READ(10, *)
+      READ(10, *)
+
+      !Read number of nodes
+      READ(10, *) self%numNodes
+
+      !Allocate required matrices and vectors
+      ALLOCATE(self%nodes(1:self%numNodes))
+      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
+
+      !Read the nodes information
+      DO e = 1, self%numNodes
+        READ(10, *) n, r(1), r(2), r(3)
+        SELECT CASE(self%geometry)
+        CASE("3DCart")
+          ALLOCATE(meshNode3Dcart::self%nodes(n)%obj)
+
+        CASE("2DCyl")
+          ALLOCATE(meshNode2DCyl:: self%nodes(n)%obj)
+          r(3) = 0.D0
+
+        CASE("2DCart")
+          ALLOCATE(meshNode2DCart:: self%nodes(n)%obj)
+          r(3) = 0.D0
+
+        CASE("1DRad")
+          ALLOCATE(meshNode1DRad:: self%nodes(n)%obj)
+          r(2:3) = 0.D0
+
+        CASE("1DCart")
+          ALLOCATE(meshNode1DCart:: self%nodes(n)%obj)
+          r(2:3) = 0.D0
+
+        END SELECT
+        CALL self%nodes(n)%obj%init(n, r)
+
+      END DO
+
+
+      !Skip comments
+      READ(10, *)
+      READ(10, *)
+
+      !Reads total number of elements (no nodes)
+      READ(10, *) totalNumElem
+
+      !conts edges and volume elements
+      SELECT TYPE(self)
+      TYPE IS(meshParticles)
+        self%numEdges = 0
+        DO e = 1, totalNumElem
+          READ(10, *) eTemp, elemType
+          SELECT CASE(self%geometry)
+          CASE("3DCart")
+            !Element type 2 is triangle in gmsh
+            IF (elemType == 2)  self%numEdges = e
+
+          CASE("2DCyl","2DCart")
+            !Element type 1 is segment in Gmsh
+            IF (elemType == 1)  self%numEdges = e
+
+          CASE("1DRad","1DCart")
+            !Element type 15 is physical point in Gmsh
+            IF (elemType == 15) self%numEdges = e
+
+          END SELECT
+
+        END DO
+
+        !Substract the number of edges to the total number of elements
+        !to obtain the number of volume elements
+        self%numVols = TotalnumElem - self%numEdges
+        ALLOCATE(self%edges(1:self%numEdges))
+        numEdges = self%numEdges
+
+        !Go back to the beggining to read elements
+        DO e=1, totalNumElem
+          BACKSPACE(10)
+        END DO
+
+      TYPE IS(meshCollisions)
+        self%numVols = TotalnumElem
+        numEdges = 0
+
+      END SELECT
+
+      !Allocates arrays
+      ALLOCATE(self%vols(1:self%numVols))
+
+      SELECT TYPE(self)
+      TYPE IS(meshParticles)
+        !Reads edges
+        DO e=1, self%numEdges
+          !Reads the edge according to the geometry
+          SELECT CASE(self%geometry)
+          CASE("3DCart")
+            READ(10, *) n, elemType, eTemp, boundaryType
+            BACKSPACE(10)
+
+            !Associate boundary condition procedure.
+            bt = getBoundaryID(boundaryType)
+
+            SELECT CASE(elemType)
+            CASE(2)
+              !Triangular surface
+              ALLOCATE(p(1:3))
+
+              READ(10, *) n, elemType, eTemp, boundaryType, eTemp, p(1:3)
+
+              ALLOCATE(meshEdge3DCartTria:: self%edges(e)%obj)
+
+            END SELECT
+
+          CASE("2DCyl")
+            ALLOCATE(p(1:2))
+
+            READ(10,*) n, elemType, eTemp, boundaryType, eTemp, p(1:2)
+            !Associate boundary condition procedure.
+            bt = getBoundaryId(boundaryType)
+
+            ALLOCATE(meshEdge2DCyl:: self%edges(e)%obj)
+
+          CASE("2DCart")
+            ALLOCATE(p(1:2))
+
+            READ(10,*) n, elemType, eTemp, boundaryType, eTemp, p(1:2)
+            !Associate boundary condition procedure.
+            bt = getBoundaryId(boundaryType)
+
+            ALLOCATE(meshEdge2DCart:: self%edges(e)%obj)
+
+          CASE("1DRad")
+            ALLOCATE(p(1:1))
+
+            READ(10, *) n, elemType, eTemp, boundaryType, eTemp, p(1)
+            !Associate boundary condition 
+            bt = getBoundaryId(boundaryType)
+
+            ALLOCATE(meshEdge1DRad:: self%edges(e)%obj)
+
+          CASE("1DCart")
+            ALLOCATE(p(1:1))
+
+            READ(10, *) n, elemType, eTemp, boundaryType, eTemp, p(1)
+            !Associate boundary condition 
+            bt = getBoundaryId(boundaryType)
+
+            ALLOCATE(meshEdge1DCart:: self%edges(e)%obj)
+
+          END SELECT
+
+          CALL self%edges(e)%obj%init(n, p, bt, boundaryType)
+          DEALLOCATE(p)
+
+        END DO
+
+      END SELECT
+
+      !Read and initialize volumes
+      DO e = 1, self%numVols
+        !Reads the volume according to the geometry
+        SELECT CASE(self%geometry)
+        CASE("3DCart")
+          READ(10, *) n, elemType
+          BACKSPACE(10)
+
+          SELECT CASE(elemType)
+          CASE(4)
+            !Tetrahedron element
+            ALLOCATE(p(1:4))
+            READ(10, *) n, elemType, eTemp, eTemp, eTemp, p(1:4)
+            ALLOCATE(meshVol3DCartTetra::  self%vols(e)%obj)
+
+          END SELECT
+
+        CASE("2DCyl")
+          READ(10,*) n, elemType
+          BACKSPACE(10)
+
+          SELECT CASE(elemType)
+          CASE (2)
+            !Triangular element
+            ALLOCATE(p(1:3))
+            READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:3)
+            ALLOCATE(meshVol2DCylTria:: self%vols(e)%obj)
+
+          CASE (3)
+            !Quadrilateral element
+            ALLOCATE(p(1:4))
+            READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:4)
+            ALLOCATE(meshVol2DCylQuad:: self%vols(e)%obj)
+
+          END SELECT
+
+        CASE("2DCart")
+          READ(10,*) n, elemType
+          BACKSPACE(10)
+
+          SELECT CASE(elemType)
+          CASE (2)
+            !Triangular element
+            ALLOCATE(p(1:3))
+            READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:3)
+            ALLOCATE(meshVol2DCartTria:: self%vols(e)%obj)
+
+          CASE (3)
+            !Quadrilateral element
+            ALLOCATE(p(1:4))
+            READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:4)
+            ALLOCATE(meshVol2DCartQuad:: self%vols(e)%obj)
+
+          END SELECT
+
+        CASE("1DRad")
+          ALLOCATE(p(1:2))
+
+          READ(10, *) n, elemType, eTemp, eTemp, eTemp, p(1:2)
+          ALLOCATE(meshVol1DRadSegm:: self%vols(e)%obj)
+
+        CASE("1DCart")
+          ALLOCATE(p(1:2))
+
+          READ(10, *) n, elemType, eTemp, eTemp, eTemp, p(1:2)
+          ALLOCATE(meshVol1DCartSegm:: self%vols(e)%obj)
+
+        END SELECT
+
+        CALL self%vols(e)%obj%init(n - numEdges, p, self%nodes)
+        DEALLOCATE(p)
+
+      END DO
+
+      CLOSE(10)
+
+    END SUBROUTINE readGmsh2
+
+END MODULE moduleMeshInputGmsh2

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

@@ -0,0 +1,209 @@
+MODULE moduleMeshOutputGmsh2
+
+  CONTAINS
+    !Prints the scattered properties of particles into the nodes
+    SUBROUTINE printOutputGmsh2(self, t)
+      USE moduleMesh
+      USE moduleRefParam
+      USE moduleSpecies
+      USE moduleOutput
+      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
+      CHARACTER (LEN=iterationDigits):: tstring
+
+      time = DBLE(t)*tauMin*ti_ref
+
+      DO i = 1, nSpecies
+        WRITE(tstring, iterationFormat) t
+        fileName='OUTPUT_' // tstring// '_' // species(i)%obj%name // '.msh'
+        WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
+        OPEN (60, file = path // folder // '/' //  fileName)
+        WRITE(60, "(A)") '$MeshFormat'
+        WRITE(60, "(A)") '2.2 0 8'
+        WRITE(60, "(A)") '$EndMeshFormat'
+        WRITE(60, "(A)") '$NodeData'
+        WRITE(60, "(A)") '1'
+        WRITE(60, "(A)") '"Density ' // species(i)%obj%name // ' (m^-3)"'
+        WRITE(60, *) 1
+        WRITE(60, *) time
+        WRITE(60, *) 3
+        WRITE(60, *) t
+        WRITE(60, *) 1
+        WRITE(60, *) 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
+        WRITE(60, "(A)") '$EndNodeData'
+        WRITE(60, "(A)") '$NodeData'
+        WRITE(60, "(A)") '1'
+        WRITE(60, "(A)") '"Velocity ' // species(i)%obj%name // ' (m/s)"'
+        WRITE(60, *) 1
+        WRITE(60, *) time
+        WRITE(60, *) 3
+        WRITE(60, *) t
+        WRITE(60, *) 3
+        WRITE(60, *) self%numNodes
+        DO n=1, self%numNodes
+          WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%velocity
+        END DO
+        WRITE(60, "(A)") '$EndNodeData'
+        WRITE(60, "(A)") '$NodeData'
+        WRITE(60, "(A)") '1'
+        WRITE(60, "(A)") '"Pressure ' // species(i)%obj%name // ' (Pa)"'
+        WRITE(60, *) 1
+        WRITE(60, *) time
+        WRITE(60, *) 3
+        WRITE(60, *) t
+        WRITE(60, *) 1
+        WRITE(60, *) self%numNodes
+        DO n=1, self%numNodes
+          WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%pressure
+        END DO
+        WRITE(60, "(A)") '$EndNodeData'
+        WRITE(60, "(A)") '$NodeData'
+        WRITE(60, "(A)") '1'
+        WRITE(60, "(A)") '"Temperature ' // species(i)%obj%name // ' (K)"'
+        WRITE(60, *) 1
+        WRITE(60, *) time
+        WRITE(60, *) 3
+        WRITE(60, *) t
+        WRITE(60, *) 1
+        WRITE(60, *) self%numNodes
+        DO n=1, self%numNodes
+          WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%temperature
+        END DO
+        WRITE(60, "(A)") '$EndNodeData'
+        CLOSE (60)
+
+      END DO
+
+    END SUBROUTINE printOutputGmsh2
+
+    !Prints the number of collisions into the volumes
+    SUBROUTINE printCollGmsh2(self, t)
+      USE moduleMesh
+      USE moduleRefParam
+      USE moduleCaseParam
+      USE moduleCollisions
+      USE moduleOutput
+      IMPLICIT NONE
+
+      CLASS(meshGeneric), INTENT(in):: self
+      INTEGER:: numEdges
+      INTEGER, INTENT(in):: t
+      INTEGER:: n
+      REAL(8):: time
+      CHARACTER(:), ALLOCATABLE:: fileName
+      CHARACTER (LEN=iterationDigits):: tstring
+
+      SELECT TYPE(self)
+      TYPE IS(meshParticles)
+        numEdges = self%numEdges
+
+      TYPE IS(meshCollisions)
+        numEdges = 0
+
+      CLASS DEFAULT
+        numEdges = 0
+
+      END SELECT
+
+      IF (collOutput) THEN
+        time = DBLE(t)*tauMin*ti_ref
+        WRITE(tstring, iterationFormat) t
+
+        fileName='OUTPUT_' // tstring// '_Collisions.msh'
+        WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
+        OPEN (60, file = path // folder // '/' //  fileName)
+        WRITE(60, "(A)") '$MeshFormat'
+        WRITE(60, "(A)") '2.2 0 8'
+        WRITE(60, "(A)") '$EndMeshFormat'
+        WRITE(60, "(A)") '$ElementData'
+        WRITE(60, "(A)") '1'
+        WRITE(60, "(A)") '"Collisions"'
+        WRITE(60, *) 1
+        WRITE(60, *) time
+        WRITE(60, *) 3
+        WRITE(60, *) t
+        WRITE(60, *) 1
+        WRITE(60, *) self%numVols
+        DO n=1, self%numVols
+          WRITE(60, "(I6,I10)") n + numEdges, self%vols(n)%obj%nColl
+        END DO
+        WRITE(60, "(A)") '$EndElementData'
+
+        CLOSE(60)
+
+      END IF
+
+    END SUBROUTINE printCollGmsh2
+
+    !Prints the electrostatic EM properties into the nodes and volumes
+    SUBROUTINE printEMGmsh2(self, t)
+      USE moduleMesh
+      USE moduleRefParam
+      USE moduleCaseParam
+      USE moduleOutput
+      IMPLICIT NONE
+
+      CLASS(meshParticles), INTENT(in):: self
+      INTEGER, INTENT(in):: t
+      INTEGER:: n, e
+      REAL(8):: time
+      CHARACTER(:), ALLOCATABLE:: fileName
+      CHARACTER (LEN=iterationDigits):: tstring
+      REAL(8):: xi(1:3)
+
+      xi = (/ 0.D0, 0.D0, 0.D0 /)
+
+      IF (emOutput) THEN
+        time = DBLE(t)*tauMin*ti_ref
+        WRITE(tstring, iterationFormat) t
+
+        fileName='OUTPUT_' // tstring// '_EMField.msh'
+        WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
+        OPEN (20, file = path // folder // '/' //  fileName)
+        WRITE(20, "(A)") '$MeshFormat'
+        WRITE(20, "(A)") '2.2 0 8'
+        WRITE(20, "(A)") '$EndMeshFormat'
+        WRITE(20, "(A)") '$NodeData'
+        WRITE(20, "(A)") '1'
+        WRITE(20, "(A)") '"Potential (V)"'
+        WRITE(20, *) 1
+        WRITE(20, *) time
+        WRITE(20, *) 3
+        WRITE(20, *) t
+        WRITE(20, *) 1
+        WRITE(20, *) self%numNodes
+        DO n=1, self%numNodes
+          WRITE(20, *) n, self%nodes(n)%obj%emData%phi*Volt_ref
+        END DO
+        WRITE(20, "(A)") '$EndNodeData'
+        
+        WRITE(20, "(A)") '$ElementData'
+        WRITE(20, "(A)") '1'
+        WRITE(20, "(A)") '"Electric Field (V/m)"'
+        WRITE(20, *) 1
+        WRITE(20, *) time
+        WRITE(20, *) 3
+        WRITE(20, *) t
+        WRITE(20, *) 3
+        WRITE(20, *) self%numVols
+        DO e=1, self%numVols
+          WRITE(20, *) e+self%numEdges, self%vols(e)%obj%gatherEF(xi)*EF_ref
+        END DO
+        WRITE(20, "(A)") '$EndElementData'
+        CLOSE(20)
+
+      END IF
+
+    END SUBROUTINE printEMGmsh2
+
+END MODULE moduleMeshOutputGmsh2

src/modules/mesh/inout/makefile

@@ -0,0 +1,4 @@
+all: gmsh2.o
+
+gmsh2.o:
+	$(MAKE) -C gmsh2 all

src/modules/mesh/makefile

@@ -1,18 +1,18 @@
-all: moduleMesh.o moduleMeshBoundary.o 3DCart.o 2DCyl.o 2DCart.o 1DRad.o 1DCart.o
+all: moduleMesh.o moduleMeshBoundary.o inout.o 3DCart.o 2DCyl.o 2DCart.o 1DRad.o 1DCart.o
 
-3DCart.o:
+3DCart.o: moduleMesh.o
 	$(MAKE) -C 3DCart all
 
-2DCyl.o:
+2DCyl.o: moduleMesh.o
 	$(MAKE) -C 2DCyl all
 
-2DCart.o:
+2DCart.o: moduleMesh.o
 	$(MAKE) -C 2DCart all
 
-1DCart.o:
+1DCart.o: moduleMesh.o
 	$(MAKE) -C 1DCart all
 
-1DRad.o:
+1DRad.o: moduleMesh.o
 	$(MAKE) -C 1DRad all
 
 moduleMesh.o: moduleMesh.f90
@@ -20,3 +20,6 @@ moduleMesh.o: moduleMesh.f90
 
 moduleMeshBoundary.o: moduleMesh.o moduleMeshBoundary.f90
 	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
+
+inout.o: 3DCart.o 2DCyl.o 2DCart.o 1DRad.o 1DCart.o
+	$(MAKE) -C inout all

src/modules/mesh/moduleMesh.f90

@@ -5,10 +5,16 @@ MODULE moduleMesh
   USE moduleBoundary
   IMPLICIT NONE
 
-  !Parent of Node element
+  !Generic mesh element
-  TYPE, PUBLIC, ABSTRACT:: meshNode
+  TYPE, PUBLIC, ABSTRACT:: meshElement
-    !Node index
+    !Index
     INTEGER:: n = 0
+    CONTAINS
+
+  END TYPE meshElement
+
+  !Parent of Node element
+  TYPE, PUBLIC, ABSTRACT, EXTENDS(meshElement):: meshNode
     !Node volume
     REAL(8):: v = 0.D0
     !Output values
@@ -44,26 +50,28 @@ MODULE moduleMesh
   !Containers for nodes in the mesh
   TYPE:: meshNodeCont
     CLASS(meshNode), ALLOCATABLE:: obj
+    CONTAINS
 
   END TYPE meshNodeCont
 
   !Type for array of boundary functions (one per species)
   TYPE, PUBLIC:: fBoundaryGeneric
     PROCEDURE(boundary_interface), POINTER, NOPASS:: apply => NULL()
+    CONTAINS
 
   END TYPE
 
   !Parent of Edge element
-  TYPE, PUBLIC, ABSTRACT:: meshEdge
+  TYPE, PUBLIC, ABSTRACT, EXTENDS(meshElement):: meshEdge
-    !Element index
-    INTEGER:: n = 0
     !Connectivity to vols
     CLASS(meshVol), POINTER:: e1 => NULL(), e2 => NULL()
+    !Connectivity to vols in meshColl
+    CLASS(meshVol), POINTER:: eColl => NULL()
     !Normal vector
     REAL(8):: normal(1:3)
     !Weight for random injection of particles
     REAL(8):: weight = 1.D0
-    !Pointer to boundary element
+    !Pointer to boundary type
     TYPE(boundaryCont), POINTER:: boundary
     !Array of functions for boundary conditions
     TYPE(fBoundaryGeneric), ALLOCATABLE:: fBoundary(:)
@@ -98,7 +106,7 @@ MODULE moduleMesh
 
     END FUNCTION getNodesEdge_interface
 
-    !Returns the intersecction between an edge and a line defined by point r0 and vector v0
+    !Returns the intersecction between an edge and a line defined by point r0
     PURE FUNCTION intersectionEdge_interface(self, r0) RESULT(r)
       IMPORT:: meshEdge
       CLASS(meshEdge), INTENT(in):: self
@@ -136,9 +144,7 @@ MODULE moduleMesh
   END TYPE meshEdgeCont
 
   !Parent of Volume element
-  TYPE, PUBLIC, ABSTRACT:: meshVol
+  TYPE, PUBLIC, ABSTRACT, EXTENDS(meshElement):: meshVol
-    !Volume index
-    INTEGER:: n = 0
     !Maximum collision rate
     REAL(8):: sigmaVrelMax = 0.D0
     !Volume
@@ -149,8 +155,6 @@ MODULE moduleMesh
     INTEGER(KIND=OMP_LOCK_KIND):: lock
     !Number of collisions per volume
     INTEGER:: nColl = 0
-    !Collisional fraction
-    REAL(8):: collFrac = 0.D0
     !Total weight of particles inside cell
     REAL(8):: totalWeight = 0.D0
     CONTAINS
@@ -159,21 +163,23 @@ MODULE moduleMesh
       PROCEDURE(randPosVol_interface),  DEFERRED, PASS::   randPos
       PROCEDURE(scatter_interface),     DEFERRED, PASS::   scatter
       PROCEDURE(gatherEF_interface),    DEFERRED, PASS::   gatherEF
+      PROCEDURE(elemK_interface),       DEFERRED, PASS::   elemK
       PROCEDURE(elemF_interface),       DEFERRED, PASS::   elemF
       PROCEDURE,                                  PASS::   findCell
       PROCEDURE(phy2log_interface),     DEFERRED, PASS::   phy2log
       PROCEDURE(inside_interface),      DEFERRED, NOPASS:: inside
       PROCEDURE(nextElement_interface), DEFERRED, PASS::   nextElement
-      PROCEDURE,                                  PASS::   collision
 
   END TYPE meshVol
 
   ABSTRACT INTERFACE
-    SUBROUTINE initVol_interface(self, n, p)
+    SUBROUTINE initVol_interface(self, n, p, nodes)
       IMPORT:: meshVol
+      IMPORT meshNodeCont
       CLASS(meshVol), INTENT(out):: self
       INTEGER, INTENT(in):: n
       INTEGER, INTENT(in):: p(:)
+      TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:)
 
     END SUBROUTINE initVol_interface
 
@@ -202,6 +208,13 @@ MODULE moduleMesh
 
     END FUNCTION getNodesVol_interface
 
+    PURE FUNCTION elemK_interface(self) RESULT(localK)
+      IMPORT:: meshVol
+      CLASS(meshVol), INTENT(in):: self
+      REAL(8), ALLOCATABLE:: localK(:,:)
+
+    END FUNCTION elemK_interface
+
     PURE FUNCTION elemF_interface(self, source) RESULT(localF)
       IMPORT:: meshVol
       CLASS(meshVol), INTENT(in):: self
@@ -211,10 +224,10 @@ MODULE moduleMesh
     END FUNCTION elemF_interface
 
     SUBROUTINE nextElement_interface(self, xi, nextElement)
-      IMPORT:: meshVol
+      IMPORT:: meshVol, meshElement
       CLASS(meshVol), INTENT(in):: self
       REAL(8), INTENT(in):: xi(1:3)
-      CLASS(*), POINTER, INTENT(out):: nextElement
+      CLASS(meshElement), POINTER, INTENT(out):: nextElement
 
     END SUBROUTINE nextElement_interface
 
@@ -248,38 +261,25 @@ MODULE moduleMesh
 
   END TYPE meshVolCont
 
-  !Abstract type of mesh
+  !Generic mesh type
-  TYPE, PUBLIC, ABSTRACT:: meshGeneric
+  TYPE, ABSTRACT:: meshGeneric
-    INTEGER:: numEdges, numNodes, numVols
+    !Geometry of the mesh
+    CHARACTER(:), ALLOCATABLE:: geometry
+    !Number of elements
+    INTEGER:: numNodes, numVols
     !Array of nodes
     TYPE(meshNodeCont), ALLOCATABLE:: nodes(:)
-    !Array of boundary elements
-    TYPE(meshEdgeCont), ALLOCATABLE:: edges(:)
     !Array of volume elements
     TYPE(meshVolCont),  ALLOCATABLE:: vols(:)
-    !Global stiffness matrix
+    PROCEDURE(readMesh_interface),     POINTER, PASS:: readMesh     => NULL()
-    REAL(8), ALLOCATABLE, DIMENSION(:,:):: K 
+    PROCEDURE(connectMesh_interface),  POINTER, PASS:: connectMesh  => NULL()
-    !Permutation matrix for P L U factorization
+    PROCEDURE(printColl_interface),    POINTER, PASS:: printColl    => NULL()
-    INTEGER, ALLOCATABLE, DIMENSION(:,:):: IPIV
-    PROCEDURE(printOutput_interface), POINTER,  PASS:: printOutput => NULL()
-    PROCEDURE(printColl_interface),   POINTER,  PASS:: printColl   => NULL()
-    PROCEDURE(printEM_interface),     POINTER,  PASS:: printEM     => NULL()
     CONTAINS
-      PROCEDURE(initMesh_interface),    DEFERRED, PASS:: init
+      PROCEDURE, PASS:: doCollisions
-      PROCEDURE(readMesh_interface),    DEFERRED, PASS:: readMesh
 
-  END TYPE meshGeneric
+  END TYPE
 
   ABSTRACT INTERFACE
-    !Inits the mesh
-    SUBROUTINE initMesh_interface(self, meshFormat)
-      IMPORT meshGeneric
-
-      CLASS(meshGeneric), INTENT(out):: self
-      CHARACTER(:), ALLOCATABLE, INTENT(in):: meshFormat
-
-    END SUBROUTINE initMesh_interface
-
     !Reads the mesh from a file
     SUBROUTINE readMesh_interface(self, filename)
       IMPORT meshGeneric
@@ -289,16 +289,15 @@ MODULE moduleMesh
 
     END SUBROUTINE readMesh_interface
 
-    !Prints Species data
+    !Connects volume and edges to the mesh
-    SUBROUTINE printOutput_interface(self, t)
+    SUBROUTINE connectMesh_interface(self)
       IMPORT meshGeneric
 
-      CLASS(meshGeneric), INTENT(in):: self
+      CLASS(meshGeneric), INTENT(inout):: self
-      INTEGER, INTENT(in):: t
 
-    END SUBROUTINE printOutput_interface
+    END SUBROUTINE connectMesh_interface
 
-    !Prints number of collisions
+    !Prints number of collisions in each volume
     SUBROUTINE printColl_interface(self, t)
       IMPORT meshGeneric
 
@@ -307,21 +306,127 @@ MODULE moduleMesh
 
     END SUBROUTINE printColl_interface
 
+  END INTERFACE
+
+  !Particle mesh
+  TYPE, EXTENDS(meshGeneric), PUBLIC:: meshParticles
+    INTEGER:: numEdges
+    !Array of boundary elements
+    TYPE(meshEdgeCont), ALLOCATABLE:: edges(:)
+    !Global stiffness matrix
+    REAL(8), ALLOCATABLE, DIMENSION(:,:):: K 
+    !Permutation matrix for P L U factorization
+    INTEGER, ALLOCATABLE, DIMENSION(:,:):: IPIV
+    PROCEDURE(printOutput_interface), POINTER, PASS:: printOutput => NULL()
+    PROCEDURE(printEM_interface),     POINTER, PASS:: printEM     => NULL()
+    PROCEDURE(doCoulomb_interface),   POINTER, PASS:: doCoulomb   => NULL()
+    CONTAINS
+      PROCEDURE, PASS:: constructGlobalK
+
+  END TYPE meshParticles
+
+  ABSTRACT INTERFACE
+    !Perform Coulomb Scattering
+    SUBROUTINE doCoulomb_interface(self)
+      IMPORT meshParticles
+
+      CLASS(meshParticles), INTENT(inout):: self
+
+    END SUBROUTINE doCoulomb_interface
+    
+    !Prints Species data
+    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, t)
-      IMPORT meshGeneric
+      IMPORT meshParticles
 
-      CLASS(meshGeneric), INTENT(in):: self
+      CLASS(meshParticles), INTENT(in):: self
       INTEGER, INTENT(in):: t
 
     END SUBROUTINE printEM_interface
 
+
   END INTERFACE
 
-  !Generic mesh
+  TYPE(meshParticles), TARGET:: mesh
-  CLASS(meshGeneric), ALLOCATABLE, TARGET:: mesh
+
+  !Collision (MCC) mesh
+  TYPE, EXTENDS(meshGeneric):: meshCollisions
+    CONTAINS
+
+  END TYPE meshCollisions
+
+  TYPE(meshCollisions), TARGET:: meshColl
+
+  ABSTRACT INTERFACE
+    SUBROUTINE readMeshColl_interface(self, filename)
+      IMPORT meshCollisions
+
+      CLASS(meshCollisions), INTENT(inout):: self
+      CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
+
+    END SUBROUTINE readMeshColl_interface
+
+    SUBROUTINE connectMeshColl_interface(self)
+      IMPORT meshParticles
+
+      CLASS(meshParticles), INTENT(inout):: self
+
+    END SUBROUTINE connectMeshColl_interface
+
+  END INTERFACE
+
+  !Pointer to mesh used for MC collisions
+  CLASS(meshGeneric), POINTER:: meshForMCC => NULL()
+
+  !Procedure to find a volume for a particle in meshColl
+  PROCEDURE(findCellColl_interface), POINTER:: findCellColl => NULL()
+
+  ABSTRACT INTERFACE
+    SUBROUTINE findCellColl_interface(part)
+      USE moduleSpecies
+
+      TYPE(particle), INTENT(inout):: part
+
+    END SUBROUTINE findCellColl_interface
+
+  END INTERFACE
 
   CONTAINS
+    !Constructs the global K matrix
+    SUBROUTINE constructGlobalK(self)
+      IMPLICIT NONE
+
+      CLASS(meshParticles), INTENT(inout):: self
+      INTEGER:: e
+      INTEGER, ALLOCATABLE:: n(:)
+      REAL(8), ALLOCATABLE:: localK(:,:)
+      INTEGER:: nNodes, i, j
+
+      DO e = 1, self%numVols
+        n      = self%vols(e)%obj%getNodes()
+        localK = self%vols(e)%obj%elemK()
+        nNodes = SIZE(n)
+
+        DO i = 1, nNodes
+          DO j = 1, nNodes
+            self%K(n(i), n(j)) = self%K(n(i), n(j)) + localK(i, j)
+
+          END DO
+
+        END DO
+
+      END DO
+
+    END SUBROUTINE constructGlobalK
+
     !Reset the output of node
     PURE SUBROUTINE resetOutput(self)
       USE moduleSpecies
@@ -343,14 +448,15 @@ MODULE moduleMesh
     !Find next cell for particle
     RECURSIVE SUBROUTINE findCell(self, part, oldCell)
       USE moduleSpecies
+      USE moduleErrors
       USE OMP_LIB
       IMPLICIT NONE
 
       CLASS(meshVol), INTENT(inout):: self
-      CLASS(meshVol), OPTIONAL, INTENT(in):: oldCell
       CLASS(particle), INTENT(inout), TARGET:: part
+      CLASS(meshVol), OPTIONAL, INTENT(in):: oldCell
       REAL(8):: xi(1:3)
-      CLASS(*), POINTER:: nextElement
+      CLASS(meshElement), POINTER:: nextElement
 
       xi = self%phy2log(part%r)
       !Checks if particle is inside 'self' cell
@@ -375,7 +481,7 @@ MODULE moduleMesh
 
         CLASS IS (meshEdge)
           !Particle encountered a surface, apply boundary
-          CALL nextElement%fBoundary(part%sp)%apply(nextElement,part)
+          CALL nextElement%fBoundary(part%species%n)%apply(nextElement,part)
 
           !If particle is still inside the domain, call findCell
           IF (part%n_in) THEN 
@@ -389,15 +495,100 @@ MODULE moduleMesh
           END IF
 
         CLASS DEFAULT
-          WRITE(*,*) "ERROR, CHECK findCell"
+          WRITE (*, "(A, I6)") "Element = ", self%n
+          CALL criticalError("No connectivity found for element", "findCell")
 
         END SELECT
+
       END IF
 
     END SUBROUTINE findCell
 
+    !If Coll and Particle are the same, simply copy the part%vol into part%volColl
+    SUBROUTINE findCellSameMesh(part)
+      USE moduleSpecies
+      IMPLICIT NONE
+
+      TYPE(particle), INTENT(inout):: part
+
+      part%volColl = part%vol
+
+    END SUBROUTINE findCellSameMesh
+
+    !TODO: try to combine this with the findCell for a regular mesh
+    !Find the volume in which particle reside in the mesh for collisions
+    SUBROUTINE findCellCollMesh(part)
+      USE moduleSpecies
+      IMPLICIT NONE
+
+      TYPE(particle), INTENT(inout):: part
+      LOGICAL:: found
+      CLASS(meshVol), POINTER:: vol
+      REAL(8), DIMENSION(1:3):: xii
+      CLASS(meshElement), POINTER:: nextElement
+
+      found = .FALSE.
+      
+      vol => meshColl%vols(part%volColl)%obj
+      DO WHILE(.NOT. found)
+        xii = vol%phy2log(part%r)
+        IF (vol%inside(xii)) THEN
+          part%volColl = vol%n
+          CALL OMP_SET_LOCK(vol%lock)
+          CALL vol%listPart_in%add(part)
+          vol%totalWeight = vol%totalWeight + part%weight
+          CALL OMP_UNSET_LOCK(vol%lock)
+          found = .TRUE.
+
+        ELSE
+          CALL vol%nextElement(xii, nextElement)
+          SELECT TYPE(nextElement)
+          CLASS IS(meshVol)
+            !Try next element
+            vol => nextElement
+
+          CLASS DEFAULT
+            !Should never happend, but just in case, stops loops
+            found = .TRUE.
+
+          END SELECT
+
+        END IF
+
+      END DO
+
+    END SUBROUTINE findCellCollMesh
+
+    !Returns index of volume associated to a position (if any)
+    !If no voulme is found, returns 0
+    !WARNING: This function is slow and should only be used in initialization phase
+    FUNCTION findCellBrute(self, r) RESULT(nVol)
+      USE moduleSpecies
+      IMPLICIT NONE
+
+      CLASS(meshGeneric), INTENT(in):: self
+      REAL(8), DIMENSION(1:3), INTENT(in):: r
+      INTEGER:: nVol
+      INTEGER:: e
+      REAL(8), DIMENSION(1:3):: xii
+
+      !Inits RESULT
+      nVol = 0
+
+      DO e = 1, self%numVols
+        xii = self%vols(e)%obj%phy2log(r)
+        IF(self%vols(e)%obj%inside(xii)) THEN
+          nVol = self%vols(e)%obj%n
+          EXIT
+
+        END IF
+
+      END DO
+
+    END FUNCTION findCellBrute
+
     !Computes collisions in element
-    SUBROUTINE collision(self)
+    SUBROUTINE doCollisions(self)
       USE moduleCollisions
       USE moduleSpecies
       USE moduleList
@@ -405,7 +596,9 @@ MODULE moduleMesh
       USE moduleRandom
       IMPLICIT NONE
 
-      CLASS(meshVol), INTENT(inout):: self
+      CLASS(meshGeneric), INTENT(inout), TARGET:: self
+      INTEGER:: e
+      CLASS(meshVol), POINTER:: vol
       INTEGER:: nPart !Number of particles inside the cell
       REAL(8):: pMax !Maximum probability of collision
       INTEGER:: rnd !random index
@@ -415,237 +608,57 @@ MODULE moduleMesh
       REAL(8):: sigmaVrelMaxNew
       TYPE(pointerArray), ALLOCATABLE:: partTemp(:)
 
-      nPart = self%listPart_in%amount
+      !$OMP DO SCHEDULE(DYNAMIC)
-      !Computes iterations if there is more than one particle in the cell
+      DO e=1, self%numVols
-      IF (nPart > 1) THEN
+        vol => self%vols(e)%obj
-        !Probability of collision
+        nPart = vol%listPart_in%amount
-        pMax = self%totalWeight*self%sigmaVrelMax*tauMin/self%volume
+        !Computes iterations if there is more than one particle in the cell
+        IF (nPart > 1) THEN
+          !Probability of collision
+          pMax = vol%totalWeight*vol%sigmaVrelMax*tauMin/vol%volume
 
-        !Increases the collisional fraction of the cell
+          !Number of collisions in the cell
-        self%collFrac = self%collFrac + REAL(nPart)*pMax*0.5D0
+          vol%nColl = NINT(REAL(nPart)*pMax*0.5D0)
-        
-        !Number of collisions in the cell
-        self%nColl = FLOOR(self%collFrac)
 
-        IF (self%nColl > 0) THEN
+          IF (vol%nColl > 0) THEN
-          !Converts the list of particles to an array for easy access
+            !Converts the list of particles to an array for easy access
-          partTemp = self%listPart_in%convert2Array()
+            partTemp = vol%listPart_in%convert2Array()
-
-        END IF
-
-        DO n = 1, self%nColl
-          !Select random numbers
-          rnd = random(1, nPart)
-          part_i => partTemp(rnd)%part
-          rnd = random(1, nPart)
-          part_j => partTemp(rnd)%part
-          ij = interactionIndex(part_i%sp, part_j%sp)
-          sigmaVrelMaxNew = 0.D0
-          DO k = 1, interactionMatrix(ij)%amount
-            CALL interactionMatrix(ij)%collisions(k)%obj%collide(self%sigmaVrelMax, sigmaVrelMaxNew, part_i, part_j)
-
-          END DO
-
-          !Update maximum cross section*v_rel per each collision
-          IF (sigmaVrelMaxNew > self%sigmaVrelMax) THEN
-            self%sigmaVrelMax = sigmaVrelMaxNew
 
           END IF
 
-          !Removes one collision from the collisional fraction
+          DO n = 1, vol%nColl
-          self%collFrac = self%collFrac - 1.D0
+            !Select random numbers
-            
+            rnd = random(1, nPart)
-        END DO
+            part_i => partTemp(rnd)%part
+            rnd = random(1, nPart)
+            part_j => partTemp(rnd)%part
+            ij = interactionIndex(part_i%species%n, part_j%species%n)
+            sigmaVrelMaxNew = 0.D0
+            DO k = 1, interactionMatrix(ij)%amount
+              CALL interactionMatrix(ij)%collisions(k)%obj%collide(vol%sigmaVrelMax, sigmaVrelMaxNew, part_i, part_j)
 
-      END IF
+            END DO
 
-    END SUBROUTINE collision
+            !Update maximum cross section*v_rel per each collision
+            IF (sigmaVrelMaxNew > vol%sigmaVrelMax) THEN
+              vol%sigmaVrelMax = sigmaVrelMaxNew
 
-    SUBROUTINE printOutputGmsh(self, t)
+            END IF
-      USE moduleRefParam
-      USE moduleSpecies
-      USE moduleOutput
-      IMPLICIT NONE
 
-      CLASS(meshGeneric), INTENT(in):: self
+          END DO
-      INTEGER, INTENT(in):: t
-      INTEGER:: n, i
-      TYPE(outputFormat):: output(1:self%numNodes)
-      REAL(8):: time
-      CHARACTER(:), ALLOCATABLE:: fileName
-      CHARACTER (LEN=iterationDigits):: tstring
 
-      time = DBLE(t)*tauMin*ti_ref
+        END IF
-
-      DO i = 1, nSpecies
-        WRITE(tstring, iterationFormat) t
-        fileName='OUTPUT_' // tstring// '_' // species(i)%obj%name // '.msh'
-        WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
-        OPEN (60, file = path // folder // '/' //  fileName)
-        WRITE(60, "(A)") '$MeshFormat'
-        WRITE(60, "(A)") '2.2 0 8'
-        WRITE(60, "(A)") '$EndMeshFormat'
-        WRITE(60, "(A)") '$NodeData'
-        WRITE(60, "(A)") '1'
-        WRITE(60, "(A)") '"Density ' // species(i)%obj%name // ' (m^-3)"'
-        WRITE(60, *) 1
-        WRITE(60, *) time
-        WRITE(60, *) 3
-        WRITE(60, *) t
-        WRITE(60, *) 1
-        WRITE(60, *) 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
-        WRITE(60, "(A)") '$EndNodeData'
-        WRITE(60, "(A)") '$NodeData'
-        WRITE(60, "(A)") '1'
-        WRITE(60, "(A)") '"Velocity ' // species(i)%obj%name // ' (m/s)"'
-        WRITE(60, *) 1
-        WRITE(60, *) time
-        WRITE(60, *) 3
-        WRITE(60, *) t
-        WRITE(60, *) 3
-        WRITE(60, *) self%numNodes
-        DO n=1, self%numNodes
-          WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%velocity
-        END DO
-        WRITE(60, "(A)") '$EndNodeData'
-        WRITE(60, "(A)") '$NodeData'
-        WRITE(60, "(A)") '1'
-        WRITE(60, "(A)") '"Pressure ' // species(i)%obj%name // ' (Pa)"'
-        WRITE(60, *) 1
-        WRITE(60, *) time
-        WRITE(60, *) 3
-        WRITE(60, *) t
-        WRITE(60, *) 1
-        WRITE(60, *) self%numNodes
-        DO n=1, self%numNodes
-          WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%pressure
-        END DO
-        WRITE(60, "(A)") '$EndNodeData'
-        WRITE(60, "(A)") '$NodeData'
-        WRITE(60, "(A)") '1'
-        WRITE(60, "(A)") '"Temperature ' // species(i)%obj%name // ' (K)"'
-        WRITE(60, *) 1
-        WRITE(60, *) time
-        WRITE(60, *) 3
-        WRITE(60, *) t
-        WRITE(60, *) 1
-        WRITE(60, *) self%numNodes
-        DO n=1, self%numNodes
-          WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%temperature
-        END DO
-        WRITE(60, "(A)") '$EndNodeData'
-        CLOSE (60)
 
       END DO
+      !$OMP END DO
 
-    END SUBROUTINE printOutputGmsh
+    END SUBROUTINE doCollisions
 
-    SUBROUTINE printCollGmsh(self, t)
+    SUBROUTINE doCoulomb(self)
-      USE moduleRefParam
+      IMPORT meshParticles
-      USE moduleCaseParam
-      USE moduleCollisions
-      USE moduleOutput
-      IMPLICIT NONE
 
-      CLASS(meshGeneric), INTENT(in):: self
+      CLASS(meshParticles), INTENT(inout):: self
-      INTEGER, INTENT(in):: t
-      INTEGER:: n
-      REAL(8):: time
-      CHARACTER(:), ALLOCATABLE:: fileName
-      CHARACTER (LEN=iterationDigits):: tstring
 
+    END SUBROUTINE doCoulomb
 
-      IF (collOutput) THEN
-        time = DBLE(t)*tauMin*ti_ref
-        WRITE(tstring, iterationFormat) t
-
-        fileName='OUTPUT_' // tstring// '_Collisions.msh'
-        WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
-        OPEN (60, file = path // folder // '/' //  fileName)
-        WRITE(60, "(A)") '$MeshFormat'
-        WRITE(60, "(A)") '2.2 0 8'
-        WRITE(60, "(A)") '$EndMeshFormat'
-        WRITE(60, "(A)") '$ElementData'
-        WRITE(60, "(A)") '1'
-        WRITE(60, "(A)") '"Collisions"'
-        WRITE(60, *) 1
-        WRITE(60, *) time
-        WRITE(60, *) 3
-        WRITE(60, *) t
-        WRITE(60, *) 1
-        WRITE(60, *) self%numVols
-        DO n=1, self%numVols
-          WRITE(60, "(I6,I10)") n + self%numEdges, self%vols(n)%obj%nColl
-        END DO
-        WRITE(60, "(A)") '$EndElementData'
-
-        CLOSE(60)
-
-      END IF
-
-    END SUBROUTINE printCollGmsh
-
-    SUBROUTINE printEMGmsh(self, t)
-      USE moduleRefParam
-      USE moduleCaseParam
-      USE moduleOutput
-      IMPLICIT NONE
-
-      CLASS(meshGeneric), INTENT(in):: self
-      INTEGER, INTENT(in):: t
-      INTEGER:: n, e
-      REAL(8):: time
-      CHARACTER(:), ALLOCATABLE:: fileName
-      CHARACTER (LEN=iterationDigits):: tstring
-      REAL(8):: xi(1:3)
-
-      xi = (/ 0.D0, 0.D0, 0.D0 /)
-
-      IF (emOutput) THEN
-        time = DBLE(t)*tauMin*ti_ref
-        WRITE(tstring, iterationFormat) t
-
-        fileName='OUTPUT_' // tstring// '_EMField.msh'
-        WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
-        OPEN (20, file = path // folder // '/' //  fileName)
-        WRITE(20, "(A)") '$MeshFormat'
-        WRITE(20, "(A)") '2.2 0 8'
-        WRITE(20, "(A)") '$EndMeshFormat'
-        WRITE(20, "(A)") '$NodeData'
-        WRITE(20, "(A)") '1'
-        WRITE(20, "(A)") '"Potential (V)"'
-        WRITE(20, *) 1
-        WRITE(20, *) time
-        WRITE(20, *) 3
-        WRITE(20, *) t
-        WRITE(20, *) 1
-        WRITE(20, *) self%numNodes
-        DO n=1, self%numNodes
-          WRITE(20, *) n, self%nodes(n)%obj%emData%phi*Volt_ref
-        END DO
-        WRITE(20, "(A)") '$EndNodeData'
-        
-        WRITE(20, "(A)") '$ElementData'
-        WRITE(20, "(A)") '1'
-        WRITE(20, "(A)") '"Electric Field (V/m)"'
-        WRITE(20, *) 1
-        WRITE(20, *) time
-        WRITE(20, *) 3
-        WRITE(20, *) t
-        WRITE(20, *) 3
-        WRITE(20, *) self%numVols
-        DO e=1, self%numVols
-          WRITE(20, *) e+self%numEdges, self%vols(e)%obj%gatherEF(xi)*EF_ref
-        END DO
-        WRITE(20, "(A)") '$EndElementData'
-        CLOSE(20)
-
-      END IF
-
-    END SUBROUTINE printEMGmsh
-  
 END MODULE moduleMesh

src/modules/mesh/moduleMeshBoundary.f90

@@ -14,8 +14,6 @@ MODULE moduleMeshBoundary
       !rpp = final position of particle
       !vpp = final velocity of particle
       REAL(8), DIMENSION(1:3):: rp, vpp
-      REAL(8):: tI
-      REAL(8):: taup !time step for reflecting process
 
       !Reflect particle velocity
       vpp = part%v - 2.D0*DOT_PRODUCT(part%v, edge%normal)*edge%normal
@@ -91,7 +89,7 @@ MODULE moduleMeshBoundary
       INTEGER:: i
 
       !Modifies particle velocity according to wall temperature
-      SELECT TYPE(bound => edge%boundary%bTypes(part%sp)%obj)
+      SELECT TYPE(bound => edge%boundary%bTypes(part%species%n)%obj)
       TYPE IS(boundaryWallTemperature)
         DO i = 1, 3
           part%v(i) = part%v(i) + bound%vTh*randomMaxwellian()
@@ -123,9 +121,9 @@ MODULE moduleMeshBoundary
       TYPE(particle), POINTER:: newElectron
       TYPE(particle), POINTER:: newIon
 
-      SELECT TYPE(bound => edge%boundary%bTypes(part%sp)%obj)
+      SELECT TYPE(bound => edge%boundary%bTypes(part%species%n)%obj)
       TYPE IS(boundaryIonization)
-        mRel = (bound%m0*species(part%sp)%obj%m)*(bound%m0+species(part%sp)%obj%m)
+        mRel = (bound%m0*part%species%m)*(bound%m0+part%species%m)
         vRel = SUM(DABS(part%v-bound%v0))
         eRel = mRel*vRel**2*5.D-1
 
@@ -133,7 +131,7 @@ MODULE moduleMeshBoundary
           ionizationRate = part%weight*bound%n0*bound%crossSection%get(eRel)*vRel
 
           !Rounds the number of particles up
-          ionizationPair = NINT(ionizationRate*bound%effectiveTime/species(bound%sp)%obj%weight)
+          ionizationPair = NINT(ionizationRate*bound%effectiveTime/bound%species%weight)
 
           !Create the new pair of particles
           DO p = 1, ionizationPair
@@ -146,8 +144,8 @@ MODULE moduleMeshBoundary
             ALLOCATE(newElectron)
             ALLOCATE(newIon)
 
-            newElectron%sp = part%sp
+            newElectron%species => part%species
-            newIon%sp      = bound%sp
+            newIon%species      => bound%species
 
             newElectron%v = v0 + (1.D0 + bound%deltaV*v0/NORM2(v0))
             newIon%v      = v0
@@ -162,13 +160,13 @@ MODULE moduleMeshBoundary
             newIon%xi      = newElectron%xi
 
             newElectron%qm = part%qm
-            SELECT TYPE(spe => species(bound%sp)%obj)
+            SELECT TYPE(spe => bound%species)
             TYPE IS(speciesCharged)
               newIon%qm = spe%qm
 
             END SELECT
 
-            newElectron%weight = species(bound%sp)%obj%weight
+            newElectron%weight = bound%species%weight
             newIon%weight      = newElectron%weight
 
             newElectron%n_in = .TRUE.

src/modules/moduleBoundary.f90

@@ -1,5 +1,6 @@
 MODULE moduleBoundary
   USE moduleTable
+  USE moduleSpecies
 
   !Generic type for boundaries
   TYPE, PUBLIC:: boundaryGeneric
@@ -36,12 +37,11 @@ MODULE moduleBoundary
   !Ionization boundary
   TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryIonization
     REAL(8):: m0, n0, v0(1:3), vTh !Properties of background neutrals.
-    INTEGER:: sp !Ion species
+    CLASS(speciesGeneric), POINTER:: species !Ion species
     TYPE(table1D):: crossSection
     REAL(8):: effectiveTime
     REAL(8):: eThreshold
     REAL(8):: deltaV
-
     CONTAINS
 
   END TYPE boundaryIonization
@@ -52,23 +52,26 @@ MODULE moduleBoundary
 
   END TYPE boundaryAxis
 
+  !Wrapper for boundary types (one per species)
   TYPE:: bTypesCont
     CLASS(boundaryGeneric), ALLOCATABLE:: obj
+    CONTAINS
 
   END TYPE bTypesCont
 
+  !Wrapper for boundary conditions
   TYPE:: boundaryCont
-    INTEGER:: id = 0
+    INTEGER:: n = 0
-    CHARACTER(:), ALLOCATABLE:: name
+    CHARACTER(:), ALLOCATABLE:: name 
-    INTEGER:: physicalSurface = 0
+    INTEGER:: physicalSurface = 0 !Physical surface as defined in the mesh file
-    CLASS(bTypesCont), ALLOCATABLE:: bTypes(:)
+    CLASS(bTypesCont), ALLOCATABLE:: bTypes(:) !Array for boundary per species
     CONTAINS
 
   END TYPE boundaryCont
 
   !Number of boundaries
   INTEGER:: nBoundary = 0
-  !Array for boundary information
+  !Array for boundaries
   TYPE(boundaryCont), ALLOCATABLE, TARGET:: boundary(:)
 
   CONTAINS
@@ -81,7 +84,7 @@ MODULE moduleBoundary
 
       id = 0
       DO i = 1, nBoundary
-        IF (physicalSurface == boundary(i)%physicalSurface) id = boundary(i)%id
+        IF (physicalSurface == boundary(i)%physicalSurface) id = boundary(i)%n
 
       END DO
 
@@ -123,7 +126,7 @@ MODULE moduleBoundary
         boundary%n0 = n0 * Vol_ref
         boundary%v0 = v0 / v_ref
         boundary%vTh = DSQRT(kb*T0/m0)/v_ref
-        boundary%sp = speciesID
+        boundary%species => species(speciesID)%obj
         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)

src/modules/moduleCollisions.f90

@@ -199,8 +199,8 @@ MODULE moduleCollisions
       sigmaVrel = self%crossSec%get(eRel)*vRel
       sigmaVrelMaxNew = sigmaVrelMaxNew + sigmaVrel
       IF (sigmaVrelMaxNew/sigmaVrelMax > random()) THEN
-        m_i = species(part_i%sp)%obj%m
+        m_i = part_i%species%m
-        m_j = species(part_j%sp)%obj%m
+        m_j = part_j%species%m
         !Applies the collision
         vCM = velocityCM(m_i, part_i%v, m_j, part_j%v)
         vp  = vRel*randomDirectionVHS()
@@ -292,11 +292,11 @@ MODULE moduleCollisions
         sigmaVrelMaxNew = sigmaVrelMaxNew + sigmaVrel
         IF (sigmaVrelMaxNew/sigmaVrelMax > random()) THEN
           !Find which particle is the ionizing electron
-          IF    (part_i%sp == self%electron%sp) THEN
+          IF    (ASSOCIATED(part_i%species, self%electron)) THEN
             electron => part_i
             neutral  => part_j
 
-          ELSEIF(part_j%sp == self%electron%sp) THEN
+          ELSEIF(ASSOCIATED(part_j%species, self%electron)) THEN
             electron => part_j
             neutral  => part_i
 
@@ -314,17 +314,18 @@ MODULE moduleCollisions
 
           !Creates a new electron from ionization
           ALLOCATE(newElectron)
-          newElectron%sp = electron%sp
+          newElectron%species => electron%species
-          newElectron%v  = vp_n
+          newElectron%v       =  vp_n
-          newElectron%r  = neutral%r
+          newElectron%r       =  neutral%r
-          newElectron%xi = neutral%xi
+          newElectron%xi      = neutral%xi
-          newElectron%n_in = .TRUE.
+          newElectron%n_in    = .TRUE.
-          newElectron%vol = neutral%vol
+          newElectron%vol     = neutral%vol
-          newElectron%weight = neutral%weight
+          newElectron%volColl = neutral%volColl
-          newElectron%qm = electron%qm
+          newElectron%weight  = neutral%weight
+          newElectron%qm      = electron%qm
 
           !Ionize neutral particle
-          SELECT TYPE(sp => species(neutral%sp)%obj)
+          SELECT TYPE(sp => neutral%species)
           TYPE IS(speciesNeutral)
             CALL sp%ionize(neutral)
 
@@ -421,11 +422,11 @@ MODULE moduleCollisions
       sigmaVrelMaxNew = sigmaVrelMaxNew + sigmaVrel
       IF (sigmaVrelMaxNew/sigmaVrelMax > random()) THEN
         !Find which particle is the ionizing electron
-        IF    (part_i%sp == self%electron%sp) THEN
+        IF    (ASSOCIATED(part_i%species, self%electron)) THEN
           electron => part_i
           ion      => part_j
 
-        ELSEIF(part_j%sp == self%electron%sp) THEN
+        ELSEIF(ASSOCIATED(part_j%species, self%electron)) THEN
           electron => part_j
           ion      => part_i
 
@@ -442,7 +443,7 @@ MODULE moduleCollisions
         electron%n_in = .FALSE.
 
         !Neutralize ion particle
-        SELECT TYPE(sp => species(ion%sp)%obj)
+        SELECT TYPE(sp => ion%species)
         TYPE IS(speciesCharged)
           CALL sp%neutralize(ion)
 
@@ -501,7 +502,7 @@ MODULE moduleCollisions
       sigmaVrel = self%crossSec%get(eRel)*vRel
       sigmaVrelMaxNew = sigmaVrelMaxNew + sigmaVrel
       IF (sigmaVrelMaxNew/sigmaVrelMax > random()) THEN
-        SELECT TYPE(sp => species(part_i%sp)%obj)
+        SELECT TYPE(sp => part_i%species)
         TYPE IS (speciesNeutral)
           !Species i is neutral, ionize particle i
           CALL sp%ionize(part_i)
@@ -512,7 +513,7 @@ MODULE moduleCollisions
 
         END SELECT
 
-        SELECT TYPE(sp => species(part_j%sp)%obj)
+        SELECT TYPE(sp => part_j%species)
         TYPE IS (speciesNeutral)
           !Species j is neutral, ionize particle j
           CALL sp%ionize(part_j)

src/modules/moduleCompTime.f90

@@ -1,11 +1,12 @@
 !Information to calculate computation time
 MODULE moduleCompTime
-  REAL(8):: tStep=0.D0
+  REAL(8):: tStep    = 0.D0
-  REAL(8):: tPush=0.D0
+  REAL(8):: tPush    = 0.D0
-  REAL(8):: tReset=0.D0
+  REAL(8):: tReset   = 0.D0
-  REAL(8):: tColl=0.D0
+  REAL(8):: tColl    = 0.D0
-  REAL(8):: tWeight=0.D0
+  REAL(8):: tCoul    = 0.D0
-  REAL(8):: tEMField=0.D0
+  REAL(8):: tWeight  = 0.D0
+  REAL(8):: tEMField = 0.D0
 
 END MODULE moduleCompTime
 

src/modules/moduleInject.f90

@@ -1,5 +1,6 @@
 !injection of particles
 MODULE moduleInject
+  USE moduleSpecies
 
   !Generic type for velocity distribution function
   TYPE, ABSTRACT:: velDistGeneric
@@ -51,7 +52,7 @@ MODULE moduleInject
     REAL(8):: T(1:3) !Temperature
     REAL(8):: n(1:3) !Direction of injection
     INTEGER:: nParticles !Number of particles to introduce each time step
-    INTEGER:: sp !Species of injection
+    CLASS(speciesGeneric), POINTER:: species !Species of injection
     INTEGER:: nEdges
     INTEGER, ALLOCATABLE:: edges(:) !Array with edges
     REAL(8), ALLOCATABLE:: cumWeight(:) !Array of cummulative probability
@@ -85,12 +86,13 @@ MODULE moduleInject
       CHARACTER(:), ALLOCATABLE, INTENT(in):: units
       INTEGER:: e, et
       INTEGER:: phSurface(1:mesh%numEdges)
+      INTEGER:: nVolColl
 
-      self%id         = i
+      self%id         =  i
-      self%vMod       = v/v_ref
+      self%vMod       =  v/v_ref
-      self%n          = n
+      self%n          =  n
-      self%T          = T/T_ref
+      self%T          =  T/T_ref
-      self%sp  = sp
+      self%species    => species(sp)%obj
       SELECT CASE(units)
       CASE ("sccm")
         !Standard cubic centimeter per minute
@@ -124,6 +126,27 @@ MODULE moduleInject
         IF (mesh%edges(e)%obj%physicalSurface == physicalSurface) THEN
           et = et + 1
           self%edges(et) = mesh%edges(e)%obj%n
+          !Assign connectivity between injection edge and meshColl volume
+          IF (ASSOCIATED(meshForMCC, meshColl)) THEN
+            nVolColl = findCellBrute(meshColl, mesh%edges(e)%obj%randPos())
+            IF (nVolColl > 0) THEN
+              mesh%edges(e)%obj%eColl => meshColl%vols(nVolColl)%obj
+
+            ELSE
+              CALL criticalError("No connection between edge and meshColl", "initInject")
+
+            END IF
+
+          ELSE
+            IF (ASSOCIATED(mesh%edges(e)%obj%e1)) THEN
+              mesh%edges(e)%obj%eColl => mesh%edges(e)%obj%e1
+
+            ELSE
+              mesh%edges(e)%obj%eColl => mesh%edges(e)%obj%e2
+
+            END IF
+
+          END IF
 
         END IF
 
@@ -224,14 +247,12 @@ MODULE moduleInject
       !$OMP SINGLE
       nMin = SUM(inject(1:(self%id-1))%nParticles) + 1
       nMax = nMin + self%nParticles - 1
-      !Assign particle type
-      partInj(nMin:nMax)%sp = self%sp
       !Assign weight to particle.
-      partInj(nMin:nMax)%weight = species(self%sp)%obj%weight
+      partInj(nMin:nMax)%weight = self%species%weight
       !Particle is considered to be outside the domain
       partInj(nMin:nMax)%n_in = .FALSE.
       !Assign charge/mass to charged particle.
-      SELECT TYPE(sp => species(self%sp)%obj)
+      SELECT TYPE(sp => self%species)
       TYPE IS(speciesCharged)
         partInj(nMin:nMax)%qm = sp%qm
 
@@ -256,6 +277,10 @@ MODULE moduleInject
           CALL criticalError("No Volume associated to edge", 'addParticles') 
 
         END IF
+        partInj(n)%volColl = randomEdge%eColl%n
+
+        !Assign particle type
+        partInj(n)%species => self%species
 
         partInj(n)%v = (/ self%v(1)%obj%randomVel(), &
                           self%v(2)%obj%randomVel(), &
@@ -264,7 +289,7 @@ MODULE moduleInject
         !Obtain natural coordinates of particle in cell
         partInj(n)%xi = mesh%vols(partInj(n)%vol)%obj%phy2log(partInj(n)%r)
         !Push new particle with the minimum time step
-        CALL solver%pusher(self%sp)%pushParticle(partInj(n), tauMin)
+        CALL solver%pusher(self%species%n)%pushParticle(partInj(n), tauMin)
         !Assign cell to new particle
         CALL solver%updateParticleCell(partInj(n))
 

src/modules/moduleInput.f90

@@ -59,7 +59,7 @@ MODULE moduleInput
       CALL checkStatus(config, "readCase")
 
       !Read injection of particles
-      CALL verboseError('Reading Interactions between species...')
+      CALL verboseError('Reading injection of particles from boundaries...')
       CALL readInject(config)
       CALL checkStatus(config, "readInject")
 
@@ -283,15 +283,15 @@ MODULE moduleInput
             !Allocate new particles
             DO p = 1, nNewPart
               ALLOCATE(partNew)
-              partNew%sp     = sp
+              partNew%species => species(sp)%obj
-              partNew%v(1)   = velocity(1) + vTh*randomMaxwellian()
+              partNew%v(1)    =  velocity(1) + vTh*randomMaxwellian()
-              partNew%v(2)   = velocity(2) + vTh*randomMaxwellian()
+              partNew%v(2)    =  velocity(2) + vTh*randomMaxwellian()
-              partNew%v(3)   = velocity(3) + vTh*randomMaxwellian()
+              partNew%v(3)    =  velocity(3) + vTh*randomMaxwellian()
-              partNew%vol    = e
+              partNew%vol     =  e
-              partNew%r      = mesh%vols(e)%obj%randPos()
+              partNew%r       =  mesh%vols(e)%obj%randPos()
-              partNew%xi     = mesh%vols(e)%obj%phy2log(partNew%r)
+              partNew%xi      =  mesh%vols(e)%obj%phy2log(partNew%r)
-              partNew%n_in   = .TRUE.
+              partNew%n_in    =  .TRUE.
-              partNew%weight = species(sp)%obj%weight
+              partNew%weight  =  species(sp)%obj%weight
               !If charged species, add qm to particle
               SELECT TYPE(sp => species(sp)%obj)
               TYPE IS (speciesCharged)
@@ -442,8 +442,8 @@ MODULE moduleInput
         !Assign shared parameters for all species
         CALL config%get(object // '.name',   species(i)%obj%name,   found)
         CALL config%get(object // '.weight', species(i)%obj%weight, found)
-        species(i)%obj%sp = i
+        species(i)%obj%n = i
-        species(i)%obj%m  = mass
+        species(i)%obj%m = mass
 
       END DO
 
@@ -483,7 +483,6 @@ MODULE moduleInput
       END DO
 
       !Set number of particles to 0 for init state
-      !TODO: In a future, this should include the particles from init states
       nPartOld = 0
 
       !Initialize the lock for the non-analogue (NA) list of particles
@@ -497,6 +496,7 @@ MODULE moduleInput
       USE moduleList
       USE moduleCollisions
       USE moduleErrors
+      USE moduleMesh
       USE OMP_LIB
       USE json_module
       IMPLICIT NONE
@@ -515,74 +515,98 @@ MODULE moduleInput
       REAL(8):: energyThreshold, energyBinding
       CHARACTER(:), ALLOCATABLE:: electron
 
-      CALL initInteractionMatrix(interactionMatrix)
+      !Firstly, checks if the object 'interactions' exists
+      CALL config%info('interactions', found)
+      IF (found) THEN
+        !Checks if MC collisions have been defined
+        CALL config%info('interactions.collisions', found)
+        IF (found) THEN
+          !Checks if a mesh for collisions has been defined
+          !The mesh will be initialized and reader in readGeometry
+          CALL config%info('interactions.meshCollisions', found)
+          IF (found) THEN
+            !Points meshForMCC to the specific mesh defined
+            meshForMCC => meshColl
 
-      !Path for collision cross-section data files
+          ELSE
-      CALL config%get('interactions.folderCollisions', pathCollisions, found)
+            !Points the meshForMCC pointer to the Particles Mesh
+            meshForMCC => mesh
 
-      !Inits lock for list of particles
+          END IF
-      CALL OMP_INIT_LOCK(lockCollisions)
 
-      CALL config%info('interactions.collisions', found, n_children = nInteractions)
+          !Inits the MCC matrix
-      DO i = 1, nInteractions
+          CALL initInteractionMatrix(interactionMatrix)
-        WRITE(iString, '(I2)') i
-        object = 'interactions.collisions(' // TRIM(iString) // ')'
-        CALL config%get(object // '.species_i', species_i, found)
-        pt_i = speciesName2Index(species_i)
-        CALL config%get(object // '.species_j', species_j, found)
-        pt_j = speciesName2Index(species_j)
-        CALL config%info(object // '.cTypes', found, n_children = nCollisions)
-        ij = interactionIndex(pt_i,pt_j)
-        !Allocates the required number of collisions per each pair of species ij
-        CALL interactionMatrix(ij)%init(nCollisions)
 
-        DO k = 1, nCollisions
+          !Path for collision cross-section data files
-          WRITE (kString, '(I2)') k
+          CALL config%get('interactions.folderCollisions', pathCollisions, found)
-          object = 'interactions.collisions(' // TRIM(iString) // ').cTypes(' // TRIM(kString) // ')'
-          !Reads the cross section file
-          CALL config%get(object // '.crossSection', crossSecFile, found)
-          crossSecFilePath = pathCollisions // crossSecFile
-          IF (.NOT. found) CALL criticalError('crossSection not found for ' // object, 'readInteractions')
-          !Reads the type of collision
-          CALL config%get(object // '.type', cType, found)
-          !Initialize collision type and reads required additional data
-          SELECT CASE(cType)
-          CASE ('elastic')
-            !Elastic collision
-            CALL initBinaryElastic(interactionMatrix(ij)%collisions(k)%obj, &
-                                   crossSecFilePath, species(pt_i)%obj%m, species(pt_j)%obj%m)
 
-          CASE ('chargeExchange')
+          !Inits lock for list of particles
-            !Resonant charge exchange
+          CALL OMP_INIT_LOCK(lockCollisions)
-            CALL initBinaryChargeExchange(interactionMatrix(ij)%collisions(k)%obj, &
-                                          crossSecFilePath, species(pt_i)%obj%m, species(pt_j)%obj%m)
 
-          CASE ('ionization')
+          CALL config%info('interactions.collisions', found, n_children = nInteractions)
-            !Electorn impact ionization
+          DO i = 1, nInteractions
-            CALL config%get(object // '.energyThreshold', energyThreshold, found)
+            WRITE(iString, '(I2)') i
-            IF (.NOT. found) CALL criticalError('energyThreshold not found for collision' // object, 'readInteractions')
+            object = 'interactions.collisions(' // TRIM(iString) // ')'
-            CALL config%get(object // '.electron', electron, found)
+            CALL config%get(object // '.species_i', species_i, found)
-            IF (.NOT. found) CALL criticalError('electron not found for collision' // object, 'readInteractions')
+            pt_i = speciesName2Index(species_i)
-            CALL initBinaryIonization(interactionMatrix(ij)%collisions(k)%obj, &
+            CALL config%get(object // '.species_j', species_j, found)
-                                      crossSecFilePath, energyThreshold, species(pt_i)%obj%m, species(pt_j)%obj%m, electron)
+            pt_j = speciesName2Index(species_j)
+            CALL config%info(object // '.cTypes', found, n_children = nCollisions)
+            ij = interactionIndex(pt_i,pt_j)
+            !Allocates the required number of collisions per each pair of species ij
+            CALL interactionMatrix(ij)%init(nCollisions)
 
-          CASE ('recombination')
+            DO k = 1, nCollisions
-            !Electorn impact ionization
+              WRITE (kString, '(I2)') k
-            CALL config%get(object // '.energyBinding', energyBinding, found)
+              object = 'interactions.collisions(' // TRIM(iString) // ').cTypes(' // TRIM(kString) // ')'
-            IF (.NOT. found) CALL criticalError('energyThreshold not found for collision' // object, 'readInteractions')
+              !Reads the cross section file
-            CALL config%get(object // '.electron', electron, found)
+              CALL config%get(object // '.crossSection', crossSecFile, found)
-            IF (.NOT. found) CALL criticalError('electron not found for collision' // object, 'readInteractions')
+              crossSecFilePath = pathCollisions // crossSecFile
-            CALL initBinaryRecombination(interactionMatrix(ij)%collisions(k)%obj, &
+              IF (.NOT. found) CALL criticalError('crossSection not found for ' // object, 'readInteractions')
-                                      crossSecFilePath, energyBinding, species(pt_i)%obj%m, species(pt_j)%obj%m, electron)
+              !Reads the type of collision
+              CALL config%get(object // '.type', cType, found)
+              !Initialize collision type and reads required additional data
+              SELECT CASE(cType)
+              CASE ('elastic')
+                !Elastic collision
+                CALL initBinaryElastic(interactionMatrix(ij)%collisions(k)%obj, &
+                                       crossSecFilePath, species(pt_i)%obj%m, species(pt_j)%obj%m)
 
-          CASE DEFAULT
+              CASE ('chargeExchange')
-            CALL criticalError('Collision type' // cType // 'not defined yet', 'readInteractions')
+                !Resonant charge exchange
+                CALL initBinaryChargeExchange(interactionMatrix(ij)%collisions(k)%obj, &
+                                              crossSecFilePath, species(pt_i)%obj%m, species(pt_j)%obj%m)
 
-          END SELECT
+              CASE ('ionization')
+                !Electorn impact ionization
+                CALL config%get(object // '.energyThreshold', energyThreshold, found)
+                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 initBinaryIonization(interactionMatrix(ij)%collisions(k)%obj, &
+                                          crossSecFilePath, energyThreshold, species(pt_i)%obj%m, species(pt_j)%obj%m, electron)
 
-        END DO
+              CASE ('recombination')
+                !Electorn impact ionization
+                CALL config%get(object // '.energyBinding', energyBinding, found)
+                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 initBinaryRecombination(interactionMatrix(ij)%collisions(k)%obj, &
+                                          crossSecFilePath, energyBinding, species(pt_i)%obj%m, species(pt_j)%obj%m, electron)
 
-      END DO
+              CASE DEFAULT
+                CALL criticalError('Collision type' // cType // 'not defined yet', 'readInteractions')
+
+              END SELECT
+
+            END DO
+
+          END DO
+
+        END IF
+
+      END IF
 
     END SUBROUTINE readInteractions
 
@@ -616,7 +640,7 @@ MODULE moduleInput
         WRITE(istring, '(i2)') i
         object = 'boundary(' // TRIM(istring) // ')'
 
-        boundary(i)%id = i
+        boundary(i)%n = i
         CALL config%get(object // '.name',            boundary(i)%name,            found)
         CALL config%get(object // '.physicalSurface', boundary(i)%physicalSurface, found)
         CALL config%info(object // '.bTypes', found, n_children = nTypes)
@@ -689,11 +713,12 @@ MODULE moduleInput
     !Read the geometry (mesh) for the case
     SUBROUTINE readGeometry(config)
       USE moduleMesh
-      USE moduleMesh3DCartRead, ONLY: mesh3DCartGeneric
+      USE moduleMeshInputGmsh2, ONLY: initGmsh2
-      USE moduleMesh2DCylRead,  ONLY: mesh2DCylGeneric
+      USE moduleMesh3DCart, ONLY: connectMesh3DCart
-      USE moduleMesh2DCartRead, ONLY: mesh2DCartGeneric
+      USE moduleMesh2DCyl,  ONLY: connectMesh2DCyl
-      USE moduleMesh1DCartRead, ONLY: mesh1DCartGeneric
+      USE moduleMesh2DCart, ONLY: connectMesh2DCart
-      USE moduleMesh1DRadRead,  ONLY: mesh1DRadGeneric
+      USE moduleMesh1DRad,  ONLY: connectMesh1DRad
+      USE moduleMesh1DCart, ONLY: connectMesh1DCart
       USE moduleErrors
       USE moduleOutput
       USE json_module
@@ -701,44 +726,79 @@ MODULE moduleInput
 
       TYPE(json_file), INTENT(inout):: config
       LOGICAL:: found
-      CHARACTER(:), ALLOCATABLE:: geometryType, meshFormat, meshFile
+      LOGICAL:: doubleMesh
+      CHARACTER(:), ALLOCATABLE:: meshFormat, meshFile
       CHARACTER(:), ALLOCATABLE:: fullPath
 
+      !Firstly, indicates if a specific mesh for MC collisions is being use
+      doubleMesh = ASSOCIATED(meshForMCC, meshColl)
+
       !Selects the type of geometry.
-      CALL config%get('geometry.type', geometryType, found)
+      CALL config%get('geometry.type', mesh%geometry, found)
-      SELECT CASE(geometryType)
+      IF (doubleMesh) meshColl%geometry = mesh%geometry
-      CASE ("3DCart")
-        !Creates a 3D cylindrical mesh
-        ALLOCATE(mesh3DCartGeneric:: mesh)
-
-      CASE ("2DCyl")
-        !Creates a 2D cylindrical mesh
-        ALLOCATE(mesh2DCylGeneric:: mesh)
-
-      CASE ("2DCart")
-        !Creates a 2D cylindrical mesh
-        ALLOCATE(mesh2DCartGeneric:: mesh)
-
-      CASE ("1DCart")
-        !Creates a 1D cartesian mesh
-        ALLOCATE(mesh1DCartGeneric:: mesh)
-
-      CASE ("1DRad")
-        !Creates a 1D cartesian mesh
-        ALLOCATE(mesh1DRadGeneric:: mesh)
-
-      CASE DEFAULT
-        CALL criticalError("Geometry type " // geometryType // " not supported.", "readGeometry")
-
-      END SELECT
 
       !Gets the type of mesh
       CALL config%get('geometry.meshType', meshFormat, found)
-      CALL mesh%init(meshFormat)
+      SELECT CASE(meshFormat)
-      !Reads the mesh
+      CASE ("gmsh2")
+        CALL initGmsh2(mesh)
+        IF (doubleMesh) CALL initGmsh2(meshColl)
+
+      CASE DEFAULT
+        CALL criticalError("Mesh format " // meshFormat // " not recogniced", "readGeometry")
+
+      END SELECT
+
+      !Reads the mesh file
       CALL config%get('geometry.meshFile', meshFile, found)
       fullpath = path // meshFile
       CALL mesh%readMesh(fullPath)
+      DEALLOCATE(fullPath, meshFile)
+      IF (doubleMesh) THEN
+        !Reads the mesh file for collisions
+        CALL config%get('interactions.meshCollisions', meshFile, found)
+        fullpath = path // meshFile
+        CALL meshColl%readMesh(fullPath)
+
+      END IF
+
+      !Creates the connectivity between elements
+      SELECT CASE(mesh%geometry)
+      CASE("3DCart")
+        mesh%connectMesh => connectMesh3DCart
+
+      CASE("2DCyl")
+        mesh%connectMesh => connectMesh2DCyl
+
+      CASE("2DCart")
+        mesh%connectMesh => connectMesh2DCart
+
+      CASE("1DRad")
+        mesh%connectMesh => connectMesh1DRad
+
+      CASE("1DCart")
+        mesh%connectMesh => connectMesh1DCart
+
+      END SELECT
+      CALL mesh%connectMesh
+
+      IF (doubleMesh) THEN
+        meshColl%connectMesh => mesh%connectMesh
+        CALL meshColl%connectMesh
+
+      END IF
+
+      !Builds the K matrix for the Particles mesh
+      CALL mesh%constructGlobalK()
+
+      !Assign the procedure to find a volume for meshColl
+      IF (doubleMesh) THEN
+        findCellColl => findCellCollMesh
+
+      ELSE
+        findCellColl => findCellSameMesh
+
+      END IF
 
     END SUBROUTINE readGeometry
 
@@ -894,7 +954,7 @@ MODULE moduleInput
       REAL(8):: v, T, m
 
       !Reads species mass
-      m = species(inj%sp)%obj%m
+      m = inj%species%m
       !Reads distribution functions for velocity
       DO i = 1, 3
         WRITE(istring, '(i2)') i

src/modules/moduleList.f90

@@ -10,7 +10,7 @@ MODULE moduleList
   END TYPE lNode
 
   TYPE listNode
-    INTEGER:: amount = 0!TODO: Make private
+    INTEGER:: amount = 0
     TYPE(lNode),POINTER:: head => NULL()
     TYPE(lNode),POINTER:: tail => NULL()
     CONTAINS

src/modules/moduleMath.f90

@@ -39,4 +39,14 @@ MODULE moduleMath
 
     END FUNCTION normalize
 
+    !Norm 1 of vector
+    PURE FUNCTION norm1(a) RESULT(b)
+      IMPLICIT NONE
+      REAL(8), DIMENSION(:), INTENT(in):: a
+      REAL(8):: b
+
+      b = SUM(DABS(a))
+
+    END FUNCTION norm1
+
 END MODULE moduleMath

src/modules/moduleOutput.f90

@@ -109,7 +109,7 @@ MODULE moduleOutput
 
         END IF
 
-        WRITE (20, "(I10, I10, 6(ES20.6E3))") t, nPartOld, tStep, tPush, tReset, tColl, tWeight, tEMField
+        WRITE (20, "(I10, I10, 7(ES20.6E3))") t, nPartOld, tStep, tPush, tReset, tColl, tCoul, tWeight, tEMField
 
         CLOSE(20)
 

src/modules/moduleSolver.f90

@@ -154,7 +154,7 @@ MODULE moduleSolver
       !$OMP DO
       DO n=1, nPartOld
         !Select species type
-        sp = partOld(n)%sp
+        sp = partOld(n)%species%n
         !Checks if the species sp is update this iteration
         IF (solver%pusher(sp)%pushSpecies) THEN
           !Push particle
@@ -483,21 +483,6 @@ MODULE moduleSolver
 
     END SUBROUTINE push1DRadCharged
 
-    !Do the collisions in all the cells
-    SUBROUTINE doCollisions()
-      USE moduleMesh
-      IMPLICIT NONE
-
-      INTEGER:: e
-
-      !$OMP DO SCHEDULE(DYNAMIC)
-      DO e=1, mesh%numVols
-        CALL mesh%vols(e)%obj%collision()
-      END DO
-      !$OMP END DO
-
-    END SUBROUTINE doCollisions
-    
     SUBROUTINE doReset()
       USE moduleSpecies
       USE moduleMesh
@@ -623,6 +608,14 @@ MODULE moduleSolver
 
       END DO
 
+      !$OMP SECTION
+      !Erase the list of particles inside the cell in coll mesh
+      DO e = 1, meshColl%numVols
+        meshColl%vols(e)%obj%totalWeight = 0.D0
+        CALL meshColl%vols(e)%obj%listPart_in%erase()
+
+      END DO
+
       !$OMP END SECTIONS
 
       !$OMP SINGLE
@@ -722,7 +715,7 @@ MODULE moduleSolver
 
         part%weight = part%weight * fractionVolume
 
-        fractionWeight = part%weight / species(part%sp)%obj%weight
+        fractionWeight = part%weight / part%species%weight
 
         IF (fractionWeight >= 2.D0) THEN
           nSplit = FLOOR(fractionWeight)
@@ -788,6 +781,7 @@ MODULE moduleSolver
 
       volOld => mesh%vols(part%vol)%obj
       CALL volOld%findCell(part)
+      CALL findCellColl(part)
       volNew => mesh%vols(part%vol)%obj
       !Call the NA shcme
       IF (ASSOCIATED(self%weightingScheme)) THEN
@@ -831,7 +825,7 @@ MODULE moduleSolver
         counterOutput=0
 
         CALL mesh%printOutput(t)
-        CALL mesh%printColl(t)
+        IF (ASSOCIATED(meshForMCC)) CALL meshForMCC%printColl(t)
         CALL mesh%printEM(t)
         WRITE(*, "(5X,A21,I10,A1,I10)") "t/tmax: ",    t,       "/", tmax
         WRITE(*, "(5X,A21,I10)")        "Particles: ", nPartOld

src/modules/moduleSpecies.f90

@@ -7,7 +7,7 @@ MODULE moduleSpecies
   TYPE, ABSTRACT:: speciesGeneric
     CHARACTER(:), ALLOCATABLE:: name
     REAL(8):: m=0.D0, weight=0.D0
-    INTEGER:: sp=0
+    INTEGER:: n=0
   END TYPE speciesGeneric
 
   TYPE, EXTENDS(speciesGeneric):: speciesNeutral
@@ -37,8 +37,9 @@ MODULE moduleSpecies
   TYPE particle
     REAL(8):: r(1:3) !Position
     REAL(8):: v(1:3) !Velocity
-    INTEGER:: sp !Particle species id
+    CLASS(speciesGeneric), POINTER:: species !Pointer to species associated with this particle
     INTEGER:: vol !Index of element in which the particle is located
+    INTEGER:: volColl !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
@@ -48,6 +49,7 @@ MODULE moduleSpecies
 
   !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
@@ -60,15 +62,18 @@ MODULE moduleSpecies
 
       CHARACTER(:), ALLOCATABLE:: speciesName
       INTEGER:: sp 
-      INTEGER:: n
+      INTEGER:: s
 
       sp = 0
-      DO n = 1, nSpecies
+      DO s = 1, nSpecies
-        IF (speciesName == species(n)%obj%name) THEN
+        IF (speciesName == species(s)%obj%name) THEN
-          sp = species(n)%obj%sp
+          sp = species(s)%obj%n
-          EXIT
+          EXIT !If a species is found, exit the loop
+
         END IF
+
       END DO
+
       !If no species is found, call a critical error
       IF (sp == 0) CALL criticalError('Species ' // speciesName // ' not found.', 'speciesName2Index')
 
@@ -83,7 +88,7 @@ MODULE moduleSpecies
       TYPE(particle), INTENT(inout):: part
 
       IF (ASSOCIATED(self%ion)) THEN
-        part%sp = self%ion%sp
+        part%species => self%ion
 
       ELSE
         CALL criticalError('No ion defined for species' // self%name, 'ionizeNeutral')
@@ -101,7 +106,7 @@ MODULE moduleSpecies
       TYPE(particle), INTENT(inout):: part
 
       IF (ASSOCIATED(self%ion)) THEN
-        part%sp = self%ion%sp
+        part%species => self%ion
 
       ELSE
         CALL criticalError('No ion defined for species' // self%name, 'ionizeCharged')
@@ -119,7 +124,7 @@ MODULE moduleSpecies
       TYPE(particle), INTENT(inout):: part
 
       IF (ASSOCIATED(self%neutral)) THEN
-        part%sp = self%neutral%sp
+        part%species => self%neutral
 
       ELSE
         CALL criticalError('No neutral defined for species' // self%name, 'neutralizeCharged')