At this stage, the input file should always start as follows:

import sys, os, os.path
filePath = os.path.abspath(os.path.dirname(sys.argv[0])) 
fileName = os.path.splitext(os.path.basename(__file__))[0]
toolsPath = ('..') + os.sep + 'tools' 
sys.path.append(toolsPath)

This bunch of commands basically defines “where I am”, “what's my name” and “where I can find the tools I need”. This information will be used by Python in order to find the files and tools needed to run the simulation.
For Python pros, you should understand their meaning at a glance. For all the others, don't worry: there is no need to go any further!

Before proceeding to the proper definition of the problem, some tools employed by POL(I)FLOWS have to be imported, as follows:

import pfemutils 
pfemutils.findbins()  
import pfem as w  
import pfemtools as wt

If, later on, we want to use the graphic interface, written in python and called viewer, we have to import it as well. It can be done as follows:

import viewer as v

All the problem properties and simulation options have to be defined inside the 'main' body of the Pyhton input file. The definition of the 'main' function can be performed as follows:

def main():

NB: remember that in the Python synthax the body of a function (everything that is described in the following sections, in this case) is defined by indentation, i.e. no parentheses required, but every line has to be indented with respect to `def main():`.

As already mentioned, for the moment the geometry and the mesh are built in Gmsh. They can be imported by POL(I)FLOWS as follows:

mshFile = filePath+os.sep+'mesh_file_name.msh'
msh = w.MshData()
msh.load(mshFile) 

If one wants the main information about the imported mesh to be printed on the screen, the following command can be added:

print msh

NB: MshData is the class managing all the data structure and the user should not have any need to interact with it.
An exception exists, nonetheless:

msh.useVoronoiCells = True/False

If useVoronoiCells is set to True, the Voronoi diagram of the domain is created after each remeshing, or just updated, if no remeshing took place.
The default value is set to False and, even if there are many occasions in which this option can come in handy, its use should be reserved to “expert” users.

After importing the mesh, the following command should be used in order to create a workspace directory for the specific simulation:

pfemutils.load(fileName)

The creation of the workspace directory is very important since all the results and outputs (see Results gathering and post-processing and Results archiving sections) will be saved in this folder.
Besides a cleaner storage of the results, by creating it the user is allowed to run multiple simulations at the same time, without overwriting the results coming from other simulations.

The workspace folder will be called as the input file (without '.py') and placed in the global 'workspace' folder.

NB: If a folder with the same name already exists, i.e. the simulation has already been run in the past, all the data contained in there will be deleted and replaced by the new ones!

The problem can be defined as follows:

pbl = w.Problem()
		   

Then, values can be assigned to problem variables in this way:

pbl.variable_name = variable_value 

The list of variables that can be modified by the user is reported in the table below:

mu
Global problem viscosity, usually used for scaling (in general set it equal to the one of the main medium).

rho0
Global problem density, usually used for scaling (in general set it equal to the one of the main medium).

alpha
Alpha parameter used in the alpha-shape technique, i.e.:

• = 0. : no elements are conserved;
• = INF : all the elements are conserved.

remeshing
Remeshing option:

• = True : a remeshing is performed at each time step (default);
• = False: remeshing is never performed and the code works as a classical Finite Element code.

beta
Distortion factor used to decide whether remesh or not, i.e.:

• = 0. : remesh at each time step;
• = INF : never remesh.

gravity
Gravity acceleration.
NB: always applied in the y-direction!

bodyForceX
Body force along x-direction.

extP
External pressure value.

nonLinAlgorithm
Non-linear algorithm:

• = 1 : Picard algorithm (default);
• = 2 : Newton-Raphson algorithm.

solScheme
Solution scheme (for the moment limited to incompressible flows):

• = 1 : monolithic scheme with pressure-stabilizing Petrov-Galerking technique (default);
• = 2 : fractional step scheme with full pressure segregation for stabilization.

scalingU
Global scaling velocity used in the definition of the stabilizing parameter of the Petrov-Galerking pressure stabilization technique.

scalingP
Global scaling pressure used in the numerical definition of the tangent stiffness matrix of the Newton-Raphson scheme.

scalingU_NR
Global scaling velocity used in the numerical definition of the tangent stiffness matrix of the Newton-Raphson scheme.

matrixScaling
Matrix scaling option:

• = True : matrices rows (i.e. nodal equations) are scaled with respect to some specific nodal quantities;
• = False : no specific matrix scaling is performed (only global, through mu and rho0) (default).

matrixScalingParameter
If matrixScaling = True, this option sets the type of scaling to be performed:

• = 1 : scaling based on the sum of node's neighbouring elements jacobians;
• = 2 : scaling based on the number of node's neighbouring elements;
• = 3 : scaling based on the area of nodal Voronoi cell.

imposedPressureAtFreeSurface
Imposed pressure at free surface option:

• `= True` : pressure is strongly imposed at the free surface;
• `= False` : pressure at free surface is introduced in a consistent way as a boundary term in the weak form (default).

NB: The strong imposition of the pressure at free surfaces leads to a violation of the principle of objectivity. The user is strongly discouraged to use this option, unless he/she is fully aware of what he/she is doing!
NB2: when imposing the value of the pressure at the free surface, only a zero-pressure condition can be imposed. This means that, consistently, the extP parameter should be put to zero as well, when imposedPressureAtFreeSurface = True!

surfTensionAtFreeSurface
Surface tension value at “free” surface (e.g. at water-air interface).
For the moment, if surfTensionAtFreeSurface = 0., surface tension is not introduced in the formulation at all.

The time integration scheme can be defined as follows:

scheme = w.TimeIntegrationSchemeName(msh, pbl) 

The list of available schemes is reported in the table below:

Then, values can be assigned to time integration scheme parameters in this way:

scheme.parameterName = parameterValue

The available time integration scheme parameters are collected in the table below:

ttot
Total time of the simulation.

dt
Time step (fixed, for the moment).

t
Current time.
NB:it does not make any sense for the user to interact directly with it, except very special occasions!

savefreq
Results archiving frequency, i.e. results are archived every savefreq time steps.

nthreads
Number of threads, for parallel computing.
NB: for the moment only the most important routines (matrix building and assembly) are parallelized, using Intel® Threading Building Blocks (TBB) library.

addRemoveNodesOption
Adding and removing nodes during the simulation, in order to keep a good spatial discretization of the domain:

• = True : at each remeshing, nodes can be added and/or removed according to geometric criteria defined through gamma and omega parameters (see below);
• = False : nodes are never added nor removed during the simulation (default).

gamma
If addRemoveNodesOption = True, a node can be removed if too close to another, if their relative distance h is lower than a reference distance h0 multiplied by gamma. So:

• gamma = 0. : nodes are never removed;
• gamma = INF : all the nodes are removed.

omega
If addRemoveNodesOption = True, a node can be added at the center of a triangle, if its area s is greater than a reference area s0 multiplied by omega. So:

• omega = 0. : a node is added at the center of every triangle;
• omega = INF : no nodes are added.

updatePressureAfterRemesh
Update pressure values at nodes with the ones obtained from the new equilibrium solution or not, i.e.:

• = True : pressure at nodes is regularly updated, even after a remeshing (default);
• = False : pressure at nodes is regularly updated, except when a remeshing takes place; in this case the pressure at nodes is kept the same as the last equilibrium configuration before the remeshing.

NB : only the “output” pressure values are affected by this parameter! Behind, the code always computes the “true” pressures, but these are inaccessible to the user or to any external software (e.g. a coupled solid solver)!
NB2 : this is a very artificial feature, that can be used in order to get rid of some annoying pressure oscillations appearing after remeshing. Be aware, if you decide to use it!
NB3 : of course, the use of this option only makes sense if the `pbl.remeshing` option is activated and the `pbl.beta` parameter is both different from 0 (remeshing at each time step) and from any value high enough to make the code never perform a remeshing!
NB4 : if two or more subsequent remeshings take place, pressure will never be updated again, until there are at least two subsequent time steps performed without remeshing. In this case results can quickly become meaningless from a physical point of view!

In a general problem, there can be many bodies interacting with each other, therefore there can be different materials.

In POL(I)FLOWS a material is called a `Medium`. Materials can be created and assigned to the different domains in the problem as follows:

w.Medium(msh, "physicalEntity", materialViscosity, materialDensity, materialTypeValue)
w.Medium(msh, physicalEntity2, materialViscosity2, materialDensity2, materialTypeValue2)
...

• physicalEntityName can be a string corresponding to the name given to a physical entity in the Gmsh file, representing a given body, domain or portion of a domain, or an `int` representing the number of a node or the reference number of a group of nodes previously defined. NB : this should be a 2D entity (a surface), in 2D simulations!

• materialViscosity is a double corresponding to the material viscosity;

• materialDensity is a double corresponding to the material initial density;

• materialTypeValue is an int corresponding to the material type, to be assigned according to the following table:

NB: since the domain is continuously rebuilt, a given element material is defined according to the materials composing the nodes belonging to the element. Hence, every node has to be assigned a material! This is done automatically at the very beginning of the simulation, but in order to do it correctly the order in which materials are defined is very important: all the SOLID materials have to be defined first, then SLAVE_FLUID materials and in the end the MASTER_FLUID material.

Since solving a differential equation without imposing the right boundary conditions (BCs) is meaningless, boundaries have to be defined on which the proper BCs can be applied.

In POL(I)FLOWS boundaries can be defined and assigned simple BCs as follows:

w.Boundary(msh, "physicalEntity", boundaryType, boundaryConditionValue)
w.Boundary(msh, physicalEntity2, boundaryType2, boundaryConditionValue2)
...

• physicalEntityName can be a string corresponding to the name given to a physical entity in the Gmsh file, representing a given boundary, or an int representing the number of a node or the reference number of a group of nodes previously defined. NB : in this case, this should be a 1D entity (a line), in 2D simulations!

• boundaryType is an int corresponding to the boundary condition type, to be assigned according to the following table:

• boundaryConditionValue is a double corresponding to the imposed velocity or pressure.
  

NB: in this way only space and time constant boundary conditions can be imposed. More complex BCs can be imposed through what in POL(I)FLOWS are called Loadings.

In POL(I)FLOWS complex boundary or initial conditions are dealt with through what we call `Loadings`.
A Loading is everything that can be applied as a boundary or intial condition, from a constant translation, to a arbitrarily-time-dependent rotation, to an initial velocity, and so on…

In POL(I)FLOWS, loadings are managed by a LoadingSet, to which different loadings can be added, allowing for a combination of many of them.

A LoadingSet can be defined as follows:

loadingset = w.LoadingSet(msh)

Then, different loadings can be added as:

loadingset.add(nb1,w.LoadingType(msh, "physicalEntity", options))
loadingset.add(nb2,wt.LoadingType(msh, physicalEntity, options))
...

• nb is an int defining the reference number of the loading. This number is arbitrary BUT each loading has to have a unique reference number! Otherwise, the last loading to be defined will overwrite all the others with the same reference number.

• LoadingType is the type of loading. This can come from one of the base classes defined in the core C++ implementation (that is `w.LoadingType(…)`) or from some derivations of the base classes in Python, usually performed in the `pfemtools.py` file (that is `wt.LoadingType(…)`). See loadings C++ base classes or Python derived classes for more explanations.

• physicalEntity can be a `string` corresponding to the name given to a physical entity in the Gmsh file, representing a given boundary or body, or an int representing the number of a node or the reference number of a group of nodes previously defined. NB : depending on the LoadingType, this can be either a 1D (line) or 2D (surface) entity.

• options depend on the loading (see loadings C++ base classes or Python derived classes)

As already mentioned, loadings can be defined either through some base classes defined in the core C++ implementation, or through some Python derived classes usually contained in the pfemtools.py file.
Below, you will find a list of the already available loadings both as base classes and derived classes.

Remember that base classes are called by POL(I)FLOWS directly, i.e. using a 'w.', while derived classes are called by pfemtools, i.e. using a 'wt.'.

Note that one of the amazing advantages in using a Python interface connected to a C++ implementation is that classes can be derived, or even created, directly in Python, without modifying the C++ source code! This means that if you need a special loading which is not in the list, you can built it your self, as explained in the Creating your own loading case section.

C++ base classes for Loadings

Python derived classes for Loadings

NB: for the moment, only time dependent loadings are implemented, which means that each of them is constant in space over the physicalEntity on which it is applied: this still allows for the generation of space-variable loadings but with a lot of work and care for the user (e.g. by creating plenty of different physical entities and applying a different loading on each of them)! The automatic managing of some more complex kinds of loading will be added in the future.

NB2: the introduction of functions in the options could look quite weird and freightening/hard to understand at the beginning, but it is not at all! On the contrary, it is a very simple and powerful tool! To make sure of this, take a look at the Using functions section!

In POL(I)FLOWS, it is possible to define different criteria to trigger the moment when a remeshing should be performed.

In POL(I)FLOWS, criteria are managed by a CriterionSet, to which different criteria can be added, allowing for a combination of many of them.

A CriterionSet can be defined as follows:

criterionset = w.CriterionSet(scheme)

Then, different criteria can be added as:

criterionset.add(nb1,w.CriterionName(msh, options))
criterionset.add(nb2,wt.CriterionName(msh, options))
...

• nb is an int defining the reference number of the criterion. This number is arbitrary BUT each criterion has to have a unique reference number! Otherwise, the last criterion to be defined will overwrite all the others with the same reference number.

• CriterionName is the type of criterion. This can come from one of the base classes defined in the core C++ implementation (that is `w.CriterionName(…)`) or from some derivations of the base classes in Python, usually performed in the `pfemtools.py` file (that is `wt.CriterionName(…)`). See criteria C++ base classes or Python derived classes for more explanations.

• options depend on the loading (see loadings C++ base classes or Python derived classes)

As for loadings, criteria can be defined either through some base classes defined in the core C++ implementation, or through some Python derived classes usually contained in the pfemtools.py file.
Below, you will find a list of the already available criteria both as base classes and derived classes.

Remember that base classes are called by POL(I)FLOWS directly, i.e. using a 'w.', while derived classes are called by pfemtools, i.e. using a 'wt.'.

C++ base classes for Remeshing criteria

Python derived classes for Remeshing criteria

As already explained, in POL(I)FLOWS different kinds of results can be stored (see Results gathering and post-processing section).
.res and .vtk are global results files and they are always archived, according to the save frequency defined in the time integration scheme (see Time integration scheme options).

In a different way, more specific data are archived in .ascii files, and stored only if the user says so. Moreover it is completely up the user to choose the kind and location of the data to be stored. This section illustrates how to tell the code which quantities have to be archived.

In POL(I)FLOWS, in order to obtain a quantity of interest a so-called extractor has to be called.

Extractors are managed by an ExtractorsManager, that can be defined as follows:

extManager = w.ExtractorsManager(msh)

Then, different extractors can be added to the manager as:

extManager.add(nb1,w.ExtractorName(msh,"physicalEntity"))
extManager.add(nb2,w.ExtractorName2(msh,physicalEntity))
...

• nb is an int defining the reference number of the extractor. This number is arbitrary BUT each extractor has to have a unique reference number! Otherwise, the last extractor to be defined will overwrite all the others with the same reference number.

• ExtractorName is the type of extractor (force extractor, velocity extractor, …). The different types are given in the Extractors list.

• physicalEntity can be a string corresponding to the name given to a physical entity in the Gmsh file, representing a given boundary or body, or an int representing the number of a node or the reference number of a group of nodes previously defined. NB : this can be either a 1D (line) or 2D (surface) entity.

Extractors list

* NB : These extractors are completely defined in Python! Actually, as for loadings, also extractors can be defined either through some base classes defined in the core C++ implementation, or through some user-defined Python derived classes usually contained in the pfemtools.py file.

Remember that base classes are called by POL(I)FLOWS directly, i.e. using a 'w.', while derived classes are called by pfemtools, i.e. using a 'wt.'.

In order to make the simulation start, the input file should always end with:

scheme.start()

NB: this instruction can be slightly modified if we want to launch the graphical interface before running the simulation (see the Activate and deactivate the POL(I)FLOWS-GUI section).

Finally, run a simulation is good, but some times you have to wait until the end before analizing the results and making sure that everything's fine.
Learn how to take advantage of the POL(I)FLOWS Graphic User Interface to enjoy a real-time visualization of your results!