In outflow

Makefile

@@ -28,7 +28,7 @@ DEBUGS := -O3
 # -I is a preprocessor flag, not a compiler flag
 
 #CXXFLAGS := $(DEBUGS) -std=c++11 -Wall -Wextra -MP -MMD
-#CXXFLAGS := $(DEBUGS) -std=c++11 -Wextra -MP -MMD -fopenmp -lm
+#CXXFLAGS := $(DEBUGS) -std=c++11 -Wall -Wextra -MP -MMD -fopenmp -lm -fsanitize=leak -fanalyzer
 CXXFLAGS := $(DEBUGS) -std=c++11 -MP -MMD -fopenmp -lm -D_GNU_SOURCE
 CPPFLAGS := -I $(INCLUDE_DIR) -Ieigen -Ijson/single_include/nlohmann -Ilibigl/include
 

include/MpsBoundaryCondition.h

@@ -38,9 +38,8 @@ class MpsBoundaryCondition {
 	 * @details    None or Periodic.
 	 * @param      PSystem    The particle system
 	 * @param      Particles  The particles data
-	 * @param      Buckets    The buckets data
 	 */
-	void setBucketBC(MpsParticleSystem *PSystem, MpsParticle *Particles, MpsBucket *Buckets);
+	void setBucketBC(MpsParticleSystem *PSystem, MpsParticle *Particles);
 
 	/**
 	 * @brief      Updates type of particle

include/MpsInflowOutflow.h

@@ -0,0 +1,125 @@
+/**
+ * @defgroup   MPSINFLOWOUTFLOW Mps Inflow Outflow
+ *
+ * @brief      This file maneges Inflow/Outflow boundary conditions.
+ *
+ * @author     Rubens Amaro
+ * @date       2022
+ */
+
+#ifndef MPS_INCLUDE_INFLOWOUTFLOW_H_
+#define MPS_INCLUDE_INFLOWOUTFLOW_H_
+
+#include "MpsParticleSystem.h"
+#include "MpsParticle.h"
+#include "MpsBucket.h"
+
+/**
+ * @brief      This class describes the Inflow/Outflow engine.
+ * @details    Creates or removes particles, updates array of particles, and 
+ * manages Inflow/Outflow.
+ */
+class MpsInflowOutflow {
+public:
+	/**
+	 * @brief      Constructs a new instance.
+	 */
+	MpsInflowOutflow();
+
+	/**
+	 * @brief      Destroys the object.
+	 */
+	virtual ~MpsInflowOutflow();
+
+	/**
+	 * @brief      This class describes a plane.
+	 * @details    The normalized plane equation is written component wise as ax + by + cz − d = 0, 
+	 * where normal = (a, b, c) and d = dot(normal,p) for a given point p on the plane
+	 * @se
+	 */
+	class Plane {
+	public:
+		double a, b, c, d;	///< Normalized plane equation components ax + by + cz − d = 0
+		double normal[3];	///< Normalized normal vector
+		int ID;				///< Plane ID
+		double press;		///< Plane inOutflow cte press
+		double velocity[3];	///< Plane inOutflow cte velocity vector
+	};
+	/**
+	 * Plane of the iNFLOW/oUTFLOW plane (interface)
+	 */
+	Plane Pio;
+
+	/**
+	 * @brief      Constructs a plane given a point and normal
+	 * @details    Constructs a plane in 3D space specifying a single point on the plane and 
+	 * the surface normal pointing toward the fluid domain
+	 * @param      PSystem  The p system
+	 * @param[in]  ii       The new value
+	 */
+	void setPlaneFromPointNormal(MpsParticleSystem *PSystem, const int ii);
+
+	/**
+	 * @brief      Calculates the euclidean signed distance between point P = (Px, Py, Pz) and the plane.
+	 * @param[in]  Px    X component of Point P
+	 * @param[in]  Py    Y component of Point P
+	 * @param[in]  Pz    Z component of Point P
+	 *
+	 * @return     The signed distance.
+	 */
+	double calcSignedDistance(const double Px, const double Py, const double Pz);
+
+
+	// Impose motion to particles
+	void imposeMotionParticles(MpsParticleSystem *PSystem, MpsParticle *Particles);
+
+	// Set initial values to inOutflow variables
+	void setInOutflowVariables(MpsParticleSystem *PSystem, MpsParticle *Particles);
+
+	// Check particles in the Inflow/Outflow region, create real and IO particles, or delete real particles
+	void checkCreateDeleteParticlesInOutflow(MpsParticleSystem *PSystem, MpsParticle *Particles);
+
+	// Swap the data between Real particles in the array Particles->numRealAndIOParticles and the array Particles->numParticles
+	void swapIdRealAndIOParticlesInOutflow(MpsParticleSystem *PSystem, MpsParticle *Particles);
+
+	/**
+	 * @brief      Verify overlaped IO-IO or IO-Real particles, assign overlaped 
+	 * IO particles as Ghost and updates the current number of ghost particles
+	 * @details    Used only in the simulations with Inflow/Ouftlow Boundary Condition
+	 * @param      PSystem  The p system
+	 */
+
+	/**
+	 * @brief      Verify overlaped IO-IO or IO-Real particles, assign overlaped 
+	 * IO particles as Ghost and updates the current number of ghost particles
+	 * @details    Used only in the simulations with Inflow/Ouftlow Boundary Condition
+	 * @param      PSystem    The physical system
+	 * @param      Particles  The particles data
+	 * @param      Buckets    The buckets data
+	 */
+	void checkOverlapedIOParticles(MpsParticleSystem *PSystem, MpsParticle *Particles, MpsBucket *Buckets);
+
+	// Check IO particles in the Inflow/Outflow region
+	// void checkIOParticlesInOutflow(MpsParticleSystem *PSystem, MpsParticle *Particles);
+
+	// Clear variables only of the inOutflow (IO) particles
+	// void clearVariablesInOutflowParticles(MpsParticleSystem *PSystem, MpsParticle *Particles);
+
+	/**
+	 * @brief      Computes the correction matrix
+	 * @details    Contribution of the neighboring particles.
+	 * @param      PSystem    The physical system
+	 * @param      Particles  The particles data
+	 * @param      Buckets    The buckets data
+	 * @see        <a href="https://doi.org/10.1016/j.cma.2010.12.001" target="_blank">Step-by-step improvement of MPS method in simulating violent free-surface motions and impact-loads</a>
+	 */
+	//void correctionMatrix(MpsParticleSystem *PSystem, MpsParticle *Particles, MpsBucket *Buckets);
+
+
+protected:
+
+private:
+
+};
+
+#endif // MPS_INCLUDE_INFLOWOUTFLOW_H_

include/MpsInputOutput.h

@@ -18,6 +18,7 @@
 #include "MpsParticleSystem.h"
 #include "MpsParticle.h"
 #include "MpsPressure.h"
+#include "MpsInflowOutflow.h"
 
 /**
  * @brief      This class manage input/output files.
@@ -46,6 +47,21 @@ class MpsInputOutput {
 	 */
 	void readInputFile(MpsParticleSystem *PSystem, MpsParticle *Particles);
 
+	/**
+	 * @brief      Check if the answer is ok
+	 * @param      answer  The answer
+	 * @param      result  The result
+	 * @return     True if ok
+	 */
+	bool checkAnswerOK(std::string& answer, bool& result);
+
+	/**
+	 * @brief      Ask yes or no
+	 * @param      message  The message
+	 * @return     Return the user answer
+	 */
+	bool question_yesno(std::string const& message);
+
 	/**
 	 * @brief      Reads data from input file .grid to class MpsParticle and allocate memory to its attributes.
 	 * @param      PSystem    The particle system
@@ -113,9 +129,8 @@ class MpsInputOutput {
 	/**
 	 * @brief      Write header for vtu files.
 	 * @param      PSystem    The particle system
-	 * @param      Particles  The particles data
 	 */
-	void writePvd(MpsParticleSystem *PSystem, MpsParticle *Particles);
+	void writePvd(MpsParticleSystem *PSystem);
 
 	/**
 	 * @brief      Writes buckets in .vtk file (Paraview).
@@ -124,11 +139,12 @@ class MpsInputOutput {
 	 * data assigned to it as a linked-list.
 	 * @see        <a href="https://doi.org/10.1007/s40571-015-0059-2" target="_blank">Performance improvements of differential operators code for MPS method on GPU</a>
 	 * @param      PSystem    The particle system
-	 * @param      Particles  The particles data
 	 */
-	void writeBuckets(MpsParticleSystem *PSystem, MpsParticle *Particles);
+	void writeBuckets(MpsParticleSystem *PSystem);
+
+
+	void writeInOutFlowPlan(MpsParticleSystem *PSystem, MpsInflowOutflow *inOutFlow);
 
-	// Pressure sensors
 
 	/**
 	 * @brief      Writes pressure at sensors in .txt file..
@@ -148,11 +164,9 @@ class MpsInputOutput {
 
 	/**
 	 * @brief      Writes header in .txt files.
-	 * @param      PSystem    The particle system
-	 * @param      Particles  The particles data
 	 * @warning    Not working.
 	 */
-	void writeHeaderTxtFiles(MpsParticleSystem *PSystem, MpsParticle *Particles);
+	void writeHeaderTxtFiles();
 
 	/**
 	 * @brief      Convert String to Char Array
@@ -164,6 +178,12 @@ class MpsInputOutput {
 	 * @brief      Delete all files inside the output folder.
 	 */
 	void deleteDirectoryFiles();
+
+	/**
+	 * @brief      Allocation of memory for Inflow/Outflow planes (interfaces)
+	 * @param      PSystem    The particle system
+	 */
+	void allocateMemoryInOutflow(MpsParticleSystem *PSystem);
 
 	std::string meshRigidFilename;
 	std::string meshDeformableFilename;

include/MpsParticle.h

@@ -106,6 +106,36 @@ class MpsParticle {
 	 * @param      PSystem    The physical system
 	 */
 	void checkParticleOutDomain(MpsParticleSystem *PSystem);
+
+	/**
+	 * @brief      Sets the particle data to ghost.
+	 * @param[in]  i        The ghost particle
+	 * @param      PSystem  The physical system
+	 */
+	void setParticleDataToGhost(const int i, MpsParticleSystem *PSystem);
+
+	/**
+	 * @brief      Move ghost particles to the last positions of the array
+	 * @param      PSystem  The physical system
+	 */
+	void moveGhostToLastPosArray(MpsParticleSystem *PSystem);
+
+	/**
+	 * @brief      Swap the data between "Last Real Particle" and i-th Ghost particle
+	 * @param[in]  i          ID of the i-th Ghost particle
+	 * @param[in]  iLastReal  ID of the Last Real particle
+	 * @param      PSystem    The physical system
+	 */
+	void swapDataLastRealPartAndGhostPart(const int i, const int iLastReal, MpsParticleSystem *PSystem);
+
+	/**
+	 * @brief      Swap the data between "Last Real+IO Particle" and "Last Real Particle"
+	 * @details    Used only in the simulations with Inflow/Ouftlow Boundary Condition
+	 * @param[in]  iLastReal    ID of the Last Real Particle previously assigned as ghost
+	 * @param[in]  iLastRealIO  The last Real+IO Particle
+	 * @param      PSystem      The p system
+	 */
+	void swapDataLastRealIOPartAndLastRealPart(const int iLastReal, const int iLastRealIO, MpsParticleSystem *PSystem);
 
 	/**
 	 * @brief      Sets the parameters.
@@ -137,13 +167,25 @@ class MpsParticle {
 	///////////// INPUTS /////////////
 	// Scalars
 	int numParticlesZero;	///< Number of particles at the initial instant of simulation
-	int numParticles;		///< Number of particles during simulation
+	int numParticles;		///< Number of real (effective) particles during simulation
+
+	int numParticlesMemory;	///< Number of particles used to allocate memory (Important for inOutflow)
+	int numIOParticles;		///< Number of inOutflow (IO) particles
+	int numCreatedParticles;///< Number of created particles in the current step
+	int numDeletedParticles;///< Number of deleted particles in the current step
+	int numRealAndIOParticles;	///< Number of particles including inOutflow (IO) particles
+	int numGhostParticles;	///< Number of ghost particles (out of domain) in the current step
+	double memoryFactor;	///< Value to multiply the initial number of particles. Used to increase the allocated memory
+
 
 	///////////// INPUTS - ARRAYS /////////////
 	// Scalars
 	int    *particleType;	///< Particle type
+	int    *particleID;		///< Global and static particle ID. The particle carries its ID value over all simulation
 	double *press;			///< Particle pressure
 	double *pressAverage;	///< Time averaged particle pressure
+	double *pressMin;		///< Particle minimum pressure
+	double *pressMax;		///< Particle maximum pressure
 	// Vectors
 	double *pos;				///< Particle position
 	double *vel;				///< Particle velocity
@@ -155,7 +197,19 @@ class MpsParticle {
 	double *correcMatrixRow2;	///< Correction matrix - Row 2
 	double *correcMatrixRow3;	///< Correction matrix - Row 3
 	double *normal;				///< Particle normal
-	double *dvelCollision;		///< Variation of velocity due collision
+	double *dvelCollision;		///< Variation of velocity due to Particle Collision (PC)
+	double *dVelShift;			///< Variation of velocity due to Particle Shifting (PS)
+	double *dPosShift;			///< Variation of position due to Particle Shifting (PS)
+
+	double *signDist;			///< Signed distance between particle and inOutflow plane
+	bool *isInIORegion;			///< True if the particle is in the inOutflow region
+	bool realParticleCreated;	///< Flag that indicates that at least one real particle was created
+	int *motherID;				///< ID of the real (effective) particle that generates the IO particle
+	int *ppeID;					///< ID of the inner Particles PPE solver
+	int *ioflowID;				///< ID of the inOutflow Plane near the IO particle
+	double *velGradientRow1;	///< Gradient of the shifted velocity x shifted position - Row 1
+	double *velGradientRow2;	///< Gradient of the shifted velocity x shifted position - Row 2
+	double *velGradientRow3;	///< Gradient of the shifted velocity x shifted position - Row 3
 
 	// Polygons
 	// Scalars
@@ -212,6 +266,7 @@ class MpsParticle {
 	double *npcdDeviation2;	///< NPCD deviation modulus
 	double *concentration;	///< Concentration
 	double *velDivergence;	///< Divergence of velocity
+	double *velDivergenceki;	///< Divergence of velocity, step k
 	double *diffusiveTerm;	///< Diffusive term
 	double *Dns, *invDns;	///< Density and its inverse