Page MenuHomec4science

No OneTemporary

File Metadata

Created
Sun, Nov 24, 21:20
diff --git a/PyGear/examples.tessel/mksnap.py b/PyGear/examples.tessel/mksnap.py
index 406c8b2..85db02b 100755
--- a/PyGear/examples.tessel/mksnap.py
+++ b/PyGear/examples.tessel/mksnap.py
@@ -1,29 +1,29 @@
#!/usr/bin/env python
from pNbody import ic
from pNbody import libutil
from numpy import *
import sys
import time
import copy
from PyGear import gadget
###################################
# generate model
-n = 100
+n = 50
nb = ic.box(n,0.5,0.5,0.5,ftype='gadget',irand=0)
#nb = ic.plummer(n,1,1,1,eps=0.1,rmax=1.,ftype='gadget')
nb.pos[:,2] = nb.pos[:,2]*0
nb.pos = nb.pos + 0.5
nb.rename('snap.dat')
nb.write()
diff --git a/src/allvars.h b/src/allvars.h
index 06ebb65..d52a171 100644
--- a/src/allvars.h
+++ b/src/allvars.h
@@ -1,1922 +1,1923 @@
/*! \file allvars.h
* \brief declares global variables.
*
* This file declares all global variables. Further variables should be added here, and declared as
* 'extern'. The actual existence of these variables is provided by the file 'allvars.c'. To produce
* 'allvars.c' from 'allvars.h', do the following:
*
* - Erase all #define's, typedef's, and enum's
* - add #include "allvars.h", delete the #ifndef ALLVARS_H conditional
* - delete all keywords 'extern'
* - delete all struct definitions enclosed in {...}, e.g.
* "extern struct global_data_all_processes {....} All;"
* becomes "struct global_data_all_processes All;"
*/
#ifndef ALLVARS_H
#define ALLVARS_H
#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_integration.h>
#include "tags.h"
#define GADGETVERSION "2.0" /*!< code version string */
#define TIMEBASE (1<<28) /*!< The simulated timespan is mapped onto the integer interval [0,TIMESPAN],
* where TIMESPAN needs to be a power of 2. Note that (1<<28) corresponds to 2^29
*/
#define MAXTOPNODES 200000 /*!< Maximum number of nodes in the top-level tree used for domain decomposition */
typedef long long peanokey; /*!< defines the variable type used for Peano-Hilbert keys */
#define BITS_PER_DIMENSION 18 /*!< Bits per dimension available for Peano-Hilbert order.
Note: If peanokey is defined as type int, the allowed maximum is 10.
If 64-bit integers are used, the maximum is 21 */
#define PEANOCELLS (((peanokey)1)<<(3*BITS_PER_DIMENSION)) /*!< The number of different Peano-Hilbert cells */
#define RNDTABLE 3000 /*!< gives the length of a table with random numbers, refreshed at every timestep.
This is used to allow application of random numbers to a specific particle
in a way that is independent of the number of processors used. */
#define MAX_REAL_NUMBER 1e37
#define MIN_REAL_NUMBER 1e-37
#define MAXLEN_FILENAME 100 /*!< Maximum number of characters for filenames (including the full path) */
#ifdef ISOTHERM_EQS
#define GAMMA (1.0) /*!< index for isothermal gas */
#else
#define GAMMA (5.0/3) /*!< adiabatic index of simulated gas */
#endif
#define GAMMA_MINUS1 (GAMMA-1)
#define HYDROGEN_MASSFRAC 0.76 /*!< mass fraction of hydrogen, relevant only for radiative cooling */
/* Some physical constants in cgs units */
#define GRAVITY 6.672e-8 /*!< Gravitational constant (in cgs units) */
#define SOLAR_MASS 1.989e33
#define SOLAR_LUM 3.826e33
#define RAD_CONST 7.565e-15
#define AVOGADRO 6.0222e23
#define BOLTZMANN 1.3806e-16
#define GAS_CONST 8.31425e7
#define C 2.9979e10
#define PLANCK 6.6262e-27
#define CM_PER_MPC 3.085678e24
#define PROTONMASS 1.6726e-24
#define ELECTRONMASS 9.10953e-28
#define THOMPSON 6.65245e-25
#define ELECTRONCHARGE 4.8032e-10
#define HUBBLE 3.2407789e-18 /* in h/sec */
#define YEAR_IN_SECOND 31536000.0 /* year in sec */
#define FEH_SOLAR 0.00181 /* used only if cooling with metal is on and chimie is off */
#define PI 3.1415926535897931
#define TWOPI 6.2831853071795862
/* Some conversion factors */
#define SEC_PER_MEGAYEAR 3.155e13
#define SEC_PER_YEAR 3.155e7
#ifndef ASMTH
#define ASMTH 1.25 /*!< ASMTH gives the scale of the short-range/long-range force split in units of FFT-mesh cells */
#endif
#ifndef RCUT
#define RCUT 4.5 /*!< RCUT gives the maximum distance (in units of the scale used for the force split) out to
which short-range forces are evaluated in the short-range tree walk. */
#endif
#define MAX_NGB 20000 /*!< defines maximum length of neighbour list */
#define MAXLEN_OUTPUTLIST 500 /*!< maxmimum number of entries in list of snapshot output times */
#define DRIFT_TABLE_LENGTH 1000 /*!< length of the lookup table used to hold the drift and kick factors */
#ifdef COSMICTIME
#define COSMICTIME_TABLE_LENGTH 1000 /*!< length of the lookup table used for the cosmic time computation */
#endif
#define MAXITER 1000 /*!< maxmimum number of steps for SPH neighbour iteration */
#ifdef DOUBLEPRECISION /*!< If defined, the variable type FLOAT is set to "double", otherwise to FLOAT */
#define FLOAT double
#else
#define FLOAT float
#endif
#ifndef TWODIMS
#define NUMDIMS 3 /*!< For 3D-normalized kernel */
#define KERNEL_COEFF_1 2.546479089470 /*!< Coefficients for SPH spline kernel and its derivative */
#define KERNEL_COEFF_2 15.278874536822
#define KERNEL_COEFF_3 45.836623610466
#define KERNEL_COEFF_4 30.557749073644
#define KERNEL_COEFF_5 5.092958178941
#define KERNEL_COEFF_6 (-15.278874536822)
#define NORM_COEFF 4.188790204786 /*!< Coefficient for kernel normalization. Note: 4.0/3 * PI = 4.188790204786 */
#else
#define NUMDIMS 2 /*!< For 2D-normalized kernel */
#define KERNEL_COEFF_1 (5.0/7*2.546479089470) /*!< Coefficients for SPH spline kernel and its derivative */
#define KERNEL_COEFF_2 (5.0/7*15.278874536822)
#define KERNEL_COEFF_3 (5.0/7*45.836623610466)
#define KERNEL_COEFF_4 (5.0/7*30.557749073644)
#define KERNEL_COEFF_5 (5.0/7*5.092958178941)
#define KERNEL_COEFF_6 (5.0/7*(-15.278874536822))
#define NORM_COEFF M_PI /*!< Coefficient for kernel normalization. */
#endif
#ifdef MULTIPHASE
#define GAS_SPH 0
#define GAS_STICKY 1
#define GAS_DARK 2
#endif
#if defined(SFR) || defined(STELLAR_PROP)
#define ST 1
#endif
#ifdef CHIMIE
#define NELEMENTS 5
#define MAXNELEMENTS 64
#define FIRST_ELEMENT "Fe"
#define FE 0
#endif
#ifdef COOLING
#define COOLING_NMETALICITIES 9
#define COOLING_NTEMPERATURES 171
#endif
#ifdef COMPUTE_VELOCITY_DISPERSION
#define VELOCITY_DISPERSION_SIZE 3
#endif
extern int SetMinTimeStepForActives;
extern int ThisTask; /*!< the rank of the local processor */
extern int NTask; /*!< number of processors */
extern int PTask; /*!< smallest integer such that NTask <= 2^PTask */
extern int NumPart; /*!< number of particles on the LOCAL processor */
extern int N_gas; /*!< number of gas particles on the LOCAL processor */
#if defined(SFR) || defined(STELLAR_PROP)
extern int N_stars; /*!< number of stars particle on the LOCAL processor */
#endif
#ifdef MULTIPHASE
extern int N_sph;
extern int N_sticky;
extern int N_stickyflaged;
extern int N_dark;
extern int NumColPotLocal; /*!< local number of potentially collisional particles */
extern int NumColPot; /*!< total number of potentially collisional particles */
extern int NumColLocal; /*!< local number of collisions */
extern int NumCol; /*!< total number of collisions */
extern int NumNoColLocal;
extern int NumNoCol;
#endif
extern long long Ntype[6]; /*!< total number of particles of each type */
extern int NtypeLocal[6]; /*!< local number of particles of each type */
extern int NumForceUpdate; /*!< number of active particles on local processor in current timestep */
extern int NumSphUpdate; /*!< number of active SPH particles on local processor in current timestep */
#ifdef CHIMIE
extern int NumStUpdate;
#endif
#ifdef TESSEL
extern int NumPTUpdate;
#endif
extern double CPUThisRun; /*!< Sums the CPU time for the process (current submission only) */
#ifdef SPLIT_DOMAIN_USING_TIME
extern double CPU_Gravity;
#endif
extern int RestartFlag; /*!< taken from command line used to start code. 0 is normal start-up from
initial conditions, 1 is resuming a run from a set of restart files, while 2
marks a restart from a snapshot file. */
extern char *Exportflag; /*!< Buffer used for flagging whether a particle needs to be exported to another process */
extern int *Ngblist; /*!< Buffer to hold indices of neighbours retrieved by the neighbour search routines */
extern int TreeReconstructFlag; /*!< Signals that a new tree needs to be constructed */
#ifdef SFR
extern int RearrangeParticlesFlag;/*!< Signals that particles must be rearanged */
#endif
extern int Flag_FullStep; /*!< This flag signals that the current step involves all particles */
extern gsl_rng *random_generator; /*!< the employed random number generator of the GSL library */
extern double RndTable[RNDTABLE]; /*!< Hold a table with random numbers, refreshed every timestep */
#ifdef SFR
extern double StarFormationRndTable[RNDTABLE]; /*!< Hold a table with random numbers, refreshed every timestep */
#endif
#ifdef FEEDBACK_WIND
extern double FeedbackWindRndTable[RNDTABLE]; /*!< Hold a table with random numbers, refreshed every timestep */
#endif
#ifdef CHIMIE
extern double ChimieRndTable[RNDTABLE]; /*!< Hold a table with random numbers, refreshed every timestep */
#endif
#ifdef CHIMIE_KINETIC_FEEDBACK
extern double ChimieKineticFeedbackRndTable[RNDTABLE]; /*!< Hold a table with random numbers, refreshed every timestep */
#endif
#ifdef AB_TURB
//Ornstein-Uhlenbeck variables
extern double StOUVar;
extern double* StOUPhases;
extern gsl_rng* StRng;
//forcing field in fourie space
extern double* StAmpl;
extern double* StAka; //phases (real part)
extern double* StAkb; //phases (imag part)
extern double* StMode;
extern int StNModes;
//integertime StTPrev; (yr : ask ?)
extern int StTPrev;
extern double StSolWeightNorm;
#endif
#ifdef PY_INTERFACE
extern int NumPartQ;
extern int N_gasQ;
extern long long NtypeQ[6]; /*!< total number of particles of each type */
extern int NtypeLocalQ[6]; /*!< local number of particles of each type */
extern double DomainCornerQ[3]; /*!< gives the lower left corner of simulation volume */
extern double DomainCenterQ[3]; /*!< gives the center of simulation volume */
extern double DomainLenQ; /*!< gives the (maximum) side-length of simulation volume */
extern double DomainFacQ; /*!< factor used for converting particle coordinates to a Peano-Hilbert mesh covering the simulation volume */
extern int DomainMyStartQ; /*!< first domain mesh cell that resides on the local processor */
extern int DomainMyLastQ; /*!< last domain mesh cell that resides on the local processor */
extern int *DomainStartListQ; /*!< a table that lists the first domain mesh cell for all processors */
extern int *DomainEndListQ; /*!< a table that lists the last domain mesh cell for all processors */
extern double *DomainWorkQ; /*!< a table that gives the total "work" due to the particles stored by each processor */
extern int *DomainCountQ; /*!< a table that gives the total number of particles held by each processor */
extern int *DomainCountSphQ; /*!< a table that gives the total number of SPH particles held by each processor */
extern int *DomainTaskQ; /*!< this table gives for each leaf of the top-level tree the processor it was assigned to */
extern peanokey *DomainKeyBufQ; /*!< this points to a buffer used during the exchange of particle data */
extern int NTopnodesQ; /*!< total number of nodes in top-level tree */
extern int NTopleavesQ; /*!< number of leaves in top-level tree. Each leaf can be assigned to a different processor */
extern void *CommBufferQ; /*!< points to communication buffer, which is used in the domain decomposition, the
parallel tree-force computation, the SPH routines, etc. */
#endif
extern double DomainCorner[3]; /*!< gives the lower left corner of simulation volume */
extern double DomainCenter[3]; /*!< gives the center of simulation volume */
extern double DomainLen; /*!< gives the (maximum) side-length of simulation volume */
extern double DomainFac; /*!< factor used for converting particle coordinates to a Peano-Hilbert mesh covering the simulation volume */
extern int DomainMyStart; /*!< first domain mesh cell that resides on the local processor */
extern int DomainMyLast; /*!< last domain mesh cell that resides on the local processor */
extern int *DomainStartList; /*!< a table that lists the first domain mesh cell for all processors */
extern int *DomainEndList; /*!< a table that lists the last domain mesh cell for all processors */
extern double *DomainWork; /*!< a table that gives the total "work" due to the particles stored by each processor */
extern int *DomainCount; /*!< a table that gives the total number of particles held by each processor */
extern int *DomainCountSph; /*!< a table that gives the total number of SPH particles held by each processor */
extern int *DomainTask; /*!< this table gives for each leaf of the top-level tree the processor it was assigned to */
extern int *DomainNodeIndex; /*!< this table gives for each leaf of the top-level tree the corresponding node of the gravitational tree */
extern FLOAT *DomainTreeNodeLen; /*!< this table gives for each leaf of the top-level tree the side-length of the corresponding node of the gravitational tree */
extern FLOAT *DomainHmax; /*!< this table gives for each leaf of the top-level tree the maximum SPH smoothing length among the particles of the corresponding node of the gravitational tree */
extern struct DomainNODE
{
FLOAT s[3]; /*!< center-of-mass coordinates */
FLOAT vs[3]; /*!< center-of-mass velocities */
FLOAT mass; /*!< mass of node */
#ifdef STELLAR_FLUX
FLOAT starlum; /*!< star luminosity of node */
#endif
#ifdef UNEQUALSOFTENINGS
#ifndef ADAPTIVE_GRAVSOFT_FORGAS
int bitflags; /*!< this bit-field encodes the particle type with the largest softening among the particles of the nodes, and whether there are particles with different softening in the node */
#else
FLOAT maxsoft; /*!< hold the maximum gravitational softening of particles in the
node if the ADAPTIVE_GRAVSOFT_FORGAS option is selected */
#endif
#endif
}
*DomainMoment; /*!< this table stores for each node of the top-level tree corresponding node data from the gravitational tree */
extern peanokey *DomainKeyBuf; /*!< this points to a buffer used during the exchange of particle data */
extern peanokey *Key; /*!< a table used for storing Peano-Hilbert keys for particles */
extern peanokey *KeySorted; /*!< holds a sorted table of Peano-Hilbert keys for all particles, used to construct top-level tree */
extern int NTopnodes; /*!< total number of nodes in top-level tree */
extern int NTopleaves; /*!< number of leaves in top-level tree. Each leaf can be assigned to a different processor */
extern struct topnode_data
{
int Daughter; /*!< index of first daughter cell (out of 8) of top-level node */
int Pstart; /*!< for the present top-level node, this gives the index of the first node in the concatenated list of topnodes collected from all processors */
int Blocks; /*!< for the present top-level node, this gives the number of corresponding nodes in the concatenated list of topnodes collected from all processors */
int Leaf; /*!< if the node is a leaf, this gives its number when all leaves are traversed in Peano-Hilbert order */
peanokey Size; /*!< number of Peano-Hilbert mesh-cells represented by top-level node */
peanokey StartKey; /*!< first Peano-Hilbert key in top-level node */
long long Count; /*!< counts the number of particles in this top-level node */
}
#ifdef PY_INTERFACE
*TopNodesQ,
#endif
*TopNodes; /*!< points to the root node of the top-level tree */
extern double TimeOfLastTreeConstruction; /*!< holds what it says, only used in connection with FORCETEST */
/* variables for input/output, usually only used on process 0 */
extern char ParameterFile[MAXLEN_FILENAME]; /*!< file name of parameterfile used for starting the simulation */
extern FILE *FdInfo; /*!< file handle for info.txt log-file. */
extern FILE *FdLog; /*!< file handle for log.txt log-file. */
extern FILE *FdEnergy; /*!< file handle for energy.txt log-file. */
#ifdef SYSTEMSTATISTICS
extern FILE *FdSystem;
#endif
extern FILE *FdTimings; /*!< file handle for timings.txt log-file. */
extern FILE *FdCPU; /*!< file handle for cpu.txt log-file. */
#ifdef FORCETEST
extern FILE *FdForceTest; /*!< file handle for forcetest.txt log-file. */
#endif
#ifdef SFR
extern FILE *FdSfr; /*!< file handle for sfr.txt log-file. */
#endif
#ifdef CHIMIE
extern FILE *FdChimie; /*!< file handle for chimie log-file. */
#endif
#ifdef MULTIPHASE
extern FILE *FdPhase; /*!< file handle for pase.txt log-file. */
extern FILE *FdSticky; /*!< file handle for sticky.txt log-file. */
#endif
#ifdef AGN_ACCRETION
extern FILE *FdAccretion; /*!< file handle for accretion.txt log-file. */
#endif
#ifdef BONDI_ACCRETION
extern FILE *FdBondi; /*!< file handle for bondi.txt log-file. */
#endif
#ifdef BUBBLES
extern FILE *FdBubble; /*!< file handle for bubble.txt log-file. */
#endif
extern double DriftTable[DRIFT_TABLE_LENGTH]; /*!< table for the cosmological drift factors */
extern double GravKickTable[DRIFT_TABLE_LENGTH]; /*!< table for the cosmological kick factor for gravitational forces */
extern double HydroKickTable[DRIFT_TABLE_LENGTH]; /*!< table for the cosmological kick factor for hydrodynmical forces */
#ifdef COSMICTIME
extern double CosmicTimeTable[COSMICTIME_TABLE_LENGTH]; /*!< table for the computation of cosmic time */
#endif
extern void *CommBuffer; /*!< points to communication buffer, which is used in the domain decomposition, the
parallel tree-force computation, the SPH routines, etc. */
/*! This structure contains data which is the SAME for all tasks (mostly code parameters read from the
* parameter file). Holding this data in a structure is convenient for writing/reading the restart file, and
* it allows the introduction of new global variables in a simple way. The only thing to do is to introduce
* them into this structure.
*/
extern struct global_data_all_processes
{
long long TotNumPart; /*!< total particle numbers (global value) */
long long TotN_gas; /*!< total gas particle number (global value) */
#ifdef PY_INTERFACE
long long TotNumPartQ; /*!< total particle numbers (global value) */
long long TotN_gasQ; /*!< total gas particle number (global value) */
int MaxPartQ; /*!< This gives the maxmimum number of particles that can be stored on one processor. */
int MaxPartSphQ; /*!< This gives the maxmimum number of SPH particles that can be stored on one processor. */
int BunchSizeSph;
int BunchSizeDensitySph;
double ForceSofteningQ;
#endif
#if defined(SFR) || defined(STELLAR_PROP)
long long TotN_stars; /*!< total stars particle number (global value) */
#endif
#ifdef MULTIPHASE
long long TotN_sph; /*!< total sph particle number (global value) */
long long TotN_sticky; /*!< total sticky particle number (global value) */
long long TotN_stickyflaged; /*!< total sticky flaged particle number (global value) */
long long TotN_stickyactive; /*!< total sticky active particle number (global value) */
long long TotN_dark; /*!< total dark particle number (global value) */
#endif
int MaxPart; /*!< This gives the maxmimum number of particles that can be stored on one processor. */
int MaxPartSph; /*!< This gives the maxmimum number of SPH particles that can be stored on one processor. */
#ifdef TESSEL
int MaxgPart;
#endif
#ifdef STELLAR_PROP
int MaxPartStars; /*!< This gives the maxmimum number of Star particles that can be stored on one processor. */
#endif
double BoxSize; /*!< Boxsize in case periodic boundary conditions are used */
int ICFormat; /*!< selects different versions of IC file-format */
int SnapFormat; /*!< selects different versions of snapshot file-formats */
int NumFilesPerSnapshot; /*!< number of files in multi-file snapshot dumps */
int NumFilesWrittenInParallel;/*!< maximum number of files that may be written simultaneously when
writing/reading restart-files, or when writing snapshot files */
int BufferSize; /*!< size of communication buffer in MB */
int BunchSizeForce; /*!< number of particles fitting into the buffer in the parallel tree-force algorithm */
int BunchSizeDensity; /*!< number of particles fitting into the communication buffer in the density computation */
int BunchSizeHydro; /*!< number of particles fitting into the communication buffer in the SPH hydrodynamical force computation */
int BunchSizeDomain; /*!< number of particles fitting into the communication buffer in the domain decomposition */
#ifdef MULTIPHASE
int BunchSizeSticky; /*!< number of particles fitting into the communication buffer in the Chimie computation */
#endif
#ifdef CHIMIE
int BunchSizeChimie; /*!< number of particles fitting into the communication buffer in the Chimie computation */
int BunchSizeStarsDensity; /*!< number of particles fitting into the communication buffer in the star density computation */
#endif
#ifdef TESSEL
int BunchSizeGhost;
#endif
double PartAllocFactor; /*!< in order to maintain work-load balance, the particle load will usually
NOT be balanced. Each processor allocates memory for PartAllocFactor times
the average number of particles to allow for that */
double TreeAllocFactor; /*!< Each processor allocates a number of nodes which is TreeAllocFactor times
the maximum(!) number of particles. Note: A typical local tree for N
particles needs usually about ~0.65*N nodes. */
#ifdef SFR
double StarsAllocFactor; /*!< Estimated fraction of gas particles that will form stars during the simulation
This allow to reduce the memory stored for stellar particles */
#endif
/* some SPH parameters */
double DesNumNgb; /*!< Desired number of SPH neighbours */
double MaxNumNgbDeviation; /*!< Maximum allowed deviation neighbour number */
double ArtBulkViscConst; /*!< Sets the parameter \f$\alpha\f$ of the artificial viscosity */
#ifdef ART_CONDUCTIVITY
double ArtCondConst; /*!< Sets the parameter \f$\alpha\f$ of the artificial conductivity */
double ArtCondThreshold;
#endif
double InitGasTemp; /*!< may be used to set the temperature in the IC's */
double MinGasTemp; /*!< may be used to set a floor for the gas temperature */
double MinEgySpec; /*!< the minimum allowed temperature expressed as energy per unit mass */
/* Usefull constants */
double Boltzmann;
double ProtonMass;
double mumh;
#ifdef COOLING
/* Cooling parameters */
double *logT;
double *logL;
gsl_interp_accel *acc_cooling_spline;
gsl_spline *cooling_spline;
double CoolingType;
char CoolingFile[MAXLEN_FILENAME]; /*!< cooling file */
double CutofCoolingTemperature;
/*
new metal dependent cooling
*/
double CoolingParameters_zmin;
double CoolingParameters_zmax;
double CoolingParameters_slz;
double CoolingParameters_tmin;
double CoolingParameters_tmax;
double CoolingParameters_slt;
double CoolingParameters_FeHSolar;
double CoolingParameters_cooling_data_max;
double CoolingParameters_cooling_data[COOLING_NMETALICITIES][COOLING_NTEMPERATURES];
int CoolingParameters_p;
int CoolingParameters_q;
#endif
#ifdef CHIMIE
int ChimieNumberOfParameterFiles;
char ChimieParameterFile[MAXLEN_FILENAME]; /*!< chimie parameter file */
double ChimieSupernovaEnergy;
double ChimieKineticFeedbackFraction;
double ChimieWindSpeed;
double ChimieWindTime;
double ChimieSNIaThermalTime;
double ChimieSNIIThermalTime;
double ChimieMaxSizeTimestep;
#endif
#if defined (CHIMIE) || defined (COOLING)
double InitGasMetallicity;
#endif
#if !defined (HEATING_PE)
double HeatingPeElectronFraction;
#endif
#if !defined (HEATING_PE) || defined (STELLAR_FLUX) || defined (EXTERNAL_FLUX)
double HeatingPeSolarEnergyDensity;
#endif
#if !defined (HEATING_PE) || defined (STELLAR_FLUX)
double HeatingPeLMRatioGas;
double HeatingPeLMRatioHalo;
double HeatingPeLMRatioDisk;
double HeatingPeLMRatioBulge;
double HeatingPeLMRatioStars;
double HeatingPeLMRatioBndry;
double HeatingPeLMRatio[6];
#endif
#ifdef EXTERNAL_FLUX
double HeatingExternalFLuxEnergyDensity;
#endif
#ifdef MULTIPHASE
double CriticalTemperature;
double CriticalEgySpec;
double CriticalNonCollisionalTemperature;
double CriticalNonCollisionalEgySpec;
#ifdef COLDGAS_CYCLE
double ColdGasCycleTransitionTime;
double ColdGasCycleTransitionParameter;
#endif
#endif
#ifdef MULTIPHASE
/* some STICKY parameters */
int StickyUseGridForCollisions;
double StickyTime; /*!< Cooling time of sticky particle collision */
double StickyCollisionTime;
double StickyLastCollisionTime;
double StickyIdleTime;
double StickyMinVelocity;
double StickyMaxVelocity;
int StickyGridNx;
int StickyGridNy;
int StickyGridNz;
double StickyGridXmin;
double StickyGridXmax;
double StickyGridYmin;
double StickyGridYmax;
double StickyGridZmin;
double StickyGridZmax;
double StickyLambda;
double StickyDensity;
double StickyDensityPower;
double StickyBetaR;
double StickyBetaT;
double StickyRsphFact; /*!< Fraction of the sph radius used in sticky particle */
#endif
#ifdef OUTERPOTENTIAL
#ifdef NFW
double HaloConcentration;
double HaloMass;
double GasMassFraction;
double NFWPotentialCte;
double Rs;
#endif
#ifdef PLUMMER
double PlummerMass;
double PlummerSoftenning;
double PlummerPotentialCte;
#endif
#ifdef MIYAMOTONAGAI
double MiyamotoNagaiMass;
double MiyamotoNagaiHr;
double MiyamotoNagaiHz;
double MiyamotoNagaiPotentialCte;
#endif
#ifdef PISOTHERM
double Rho0;
double Rc;
double PisothermPotentialCte;
double GasMassFraction;
double PotentialInf;
gsl_function PotentialF;
gsl_integration_workspace *Potentialw;
#endif
#ifdef CORIOLIS
double CoriolisOmegaX;
double CoriolisOmegaY;
double CoriolisOmegaZ;
double CoriolisOmegaX0;
double CoriolisOmegaY0;
double CoriolisOmegaZ0;
#endif
#endif
#ifdef SFR
int StarFormationNStarsFromGas;
double StarFormationStarMass;
double StarFormationMgMsFraction;
int StarFormationType;
double StarFormationCstar;
double StarFormationTime;
double StarFormationDensity;
double StarFormationTemperature;
double ThresholdDensity;
#endif
#ifdef FEEDBACK
double SupernovaTime;
#endif
#ifdef FEEDBACK_WIND
double SupernovaWindEgySpecPerMassUnit;
double SupernovaWindFractionInEgyKin;
double SupernovaWindParameter;
double SupernovaWindSpeed;
double SupernovaWindIntAccuracy;
#endif
#ifdef AGN_ACCRETION
double TimeBetAccretion;
double AccretionRadius;
double AGNFactor;
double MinMTotInRa;
double TimeLastAccretion;
double LastMTotInRa;
double MTotInRa;
double dMTotInRa;
#endif
#ifdef BUBBLES
char BubblesInitFile[MAXLEN_FILENAME]; /*!< bubble file */
double *BubblesTime;
double *BubblesD;
double *BubblesR;
double *BubblesE;
double *BubblesA;
double *BubblesB;
int BubblesIndex;
double BubblesAlpha;
double BubblesBeta;
double BubblesDelta;
double BubblesRadiusFactor;
double EnergyBubbles;
#endif
#ifdef AGN_HEATING
double AGNHeatingPower;
double AGNHeatingRmax;
#endif
#ifdef BONDI_ACCRETION
double BondiEfficiency;
double BondiBlackHoleMass;
double BondiHsmlFactor;
double BondiPower;
double BondiTimeBet;
double BondiTimeLast;
#endif
#if defined (AGN_ACCRETION) || defined (BONDI_ACCRETION)
double LightSpeed;
#endif
#if defined(ART_VISCO_MM)|| defined(ART_VISCO_RO) || defined(ART_VISCO_CD)
double ArtBulkViscConstMin;
double ArtBulkViscConstMax;
double ArtBulkViscConstL;
#endif
#ifdef AB_TURB
double StDecay;
double StEnergy;
double StDtFreq;
double StKmin;
double StKmax;
double StSolWeight;
double StAmplFac;
int StSpectForm;
int StSeed;
#endif
/* some force counters */
long long TotNumOfForces; /*!< counts total number of force computations */
long long NumForcesSinceLastDomainDecomp; /*!< count particle updates since last domain decomposition */
/* system of units */
double G; /*!< Gravity-constant in internal units */
double UnitTime_in_s; /*!< factor to convert internal time unit to seconds/h */
double UnitMass_in_g; /*!< factor to convert internal mass unit to grams/h */
double UnitVelocity_in_cm_per_s; /*!< factor to convert intqernal velocity unit to cm/sec */
double UnitLength_in_cm; /*!< factor to convert internal length unit to cm/h */
double UnitPressure_in_cgs; /*!< factor to convert internal pressure unit to cgs units (little 'h' still around!) */
double UnitDensity_in_cgs; /*!< factor to convert internal length unit to g/cm^3*h^2 */
double UnitCoolingRate_in_cgs; /*!< factor to convert internal cooling rate to cgs units */
double UnitEnergy_in_cgs; /*!< factor to convert internal energy to cgs units */
double UnitTime_in_Megayears; /*!< factor to convert internal time to megayears/h */
double GravityConstantInternal; /*!< If set to zero in the parameterfile, the internal value of the
gravitational constant is set to the Newtonian value based on the system of
units specified. Otherwise the value provided is taken as internal gravity constant G. */
/* Cosmological parameters */
double Hubble; /*!< Hubble-constant in internal units */
double Omega0; /*!< matter density in units of the critical density (at z=0)*/
double OmegaLambda; /*!< vaccum energy density relative to crictical density (at z=0) */
double OmegaBaryon; /*!< baryon density in units of the critical density (at z=0)*/
double HubbleParam; /*!< little `h', i.e. Hubble constant in units of 100 km/s/Mpc. Only needed to get absolute physical values for cooling physics */
/* Code options */
int ComovingIntegrationOn; /*!< flags that comoving integration is enabled */
int PeriodicBoundariesOn; /*!< flags that periodic boundaries are enabled */
int ResubmitOn; /*!< flags that automatic resubmission of job to queue system is enabled */
int TypeOfOpeningCriterion; /*!< determines tree cell-opening criterion: 0 for Barnes-Hut, 1 for relative criterion */
int TypeOfTimestepCriterion; /*!< gives type of timestep criterion (only 0 supported right now - unlike gadget-1.1) */
int OutputListOn; /*!< flags that output times are listed in a specified file */
/* Parameters determining output frequency */
int SnapshotFileCount; /*!< number of snapshot that is written next */
double TimeBetSnapshot; /*!< simulation time interval between snapshot files */
double TimeOfFirstSnapshot; /*!< simulation time of first snapshot files */
double CpuTimeBetRestartFile; /*!< cpu-time between regularly generated restart files */
double TimeLastRestartFile; /*!< cpu-time when last restart-file was written */
double TimeBetStatistics; /*!< simulation time interval between computations of energy statistics */
double TimeLastStatistics; /*!< simulation time when the energy statistics was computed the last time */
int NumCurrentTiStep; /*!< counts the number of system steps taken up to this point */
/* Current time of the simulation, global step, and end of simulation */
double Time; /*!< current time of the simulation */
double TimeBegin; /*!< time of initial conditions of the simulation */
double TimeStep; /*!< difference between current times of previous and current timestep */
double TimeMax; /*!< marks the point of time until the simulation is to be evolved */
/* variables for organizing discrete timeline */
double Timebase_interval; /*!< factor to convert from floating point time interval to integer timeline */
int Ti_Current; /*!< current time on integer timeline */
int Ti_nextoutput; /*!< next output time on integer timeline */
#ifdef FLEXSTEPS
int PresentMinStep; /*!< If FLEXSTEPS is used, particle timesteps are chosen as multiples of the present minimum timestep. */
int PresentMaxStep; /*!< If FLEXSTEPS is used, this is the maximum timestep in timeline units, rounded down to the next power 2 division */
#endif
#ifdef PMGRID
int PM_Ti_endstep; /*!< begin of present long-range timestep */
int PM_Ti_begstep; /*!< end of present long-range timestep */
#endif
/* Placement of PM grids */
#ifdef PMGRID
double Asmth[2]; /*!< Gives the scale of the long-range/short-range split (in mesh-cells), both for the coarse and the high-res mesh */
double Rcut[2]; /*!< Gives the maximum radius for which the short-range force is evaluated with the tree (in mesh-cells), both for the coarse and the high-res mesh */
double Corner[2][3]; /*!< lower left corner of coarse and high-res PM-mesh */
double UpperCorner[2][3]; /*!< upper right corner of coarse and high-res PM-mesh */
double Xmintot[2][3]; /*!< minimum particle coordinates both for coarse and high-res PM-mesh */
double Xmaxtot[2][3]; /*!< maximum particle coordinates both for coarse and high-res PM-mesh */
double TotalMeshSize[2]; /*!< total extension of coarse and high-res PM-mesh */
#endif
/* Variables that keep track of cumulative CPU consumption */
double TimeLimitCPU; /*!< CPU time limit as defined in parameterfile */
double CPU_TreeConstruction; /*!< time spent for constructing the gravitational tree */
double CPU_TreeWalk; /*!< actual time spent for pure tree-walks */
double CPU_Gravity; /*!< cumulative time used for gravity computation (tree-algorithm only) */
double CPU_Potential; /*!< time used for computing gravitational potentials */
double CPU_Domain; /*!< cumulative time spent for domain decomposition */
double CPU_Snapshot; /*!< time used for writing snapshot files */
double CPU_Total; /*!< cumulative time spent for domain decomposition */
double CPU_CommSum; /*!< accumulated time used for communication, and for collecting partial results, in tree-gravity */
double CPU_Imbalance; /*!< cumulative time lost accross all processors as work-load imbalance in gravitational tree */
double CPU_HydCompWalk; /*!< time used for actual SPH computations, including neighbour search */
double CPU_HydCommSumm; /*!< cumulative time used for communication in SPH, and for collecting partial results */
double CPU_HydImbalance; /*!< cumulative time lost due to work-load imbalance in SPH */
double CPU_Hydro; /*!< cumulative time spent for SPH related computations */
#ifdef SFR
double CPU_StarFormation; /*!< cumulative time spent for star formation computations */
#endif
#ifdef CHIMIE
double CPU_Chimie; /*!< cumulative time spent for chimie computations */
double CPU_ChimieDensCompWalk;
double CPU_ChimieDensCommSumm;
double CPU_ChimieDensImbalance;
double CPU_ChimieDensEnsureNgb;
double CPU_ChimieCompWalk;
double CPU_ChimieCommSumm;
double CPU_ChimieImbalance;
#endif
#ifdef MULTIPHASE
double CPU_Sticky; /*!< cumulative time spent for sticky computations */
#endif
double CPU_EnsureNgb; /*!< time needed to iterate on correct neighbour numbers */
double CPU_Predict; /*!< cumulative time to drift the system forward in time, including dynamic tree updates */
double CPU_TimeLine; /*!< time used for determining new timesteps, and for organizing the timestepping, including kicks of active particles */
double CPU_PM; /*!< time used for long-range gravitational force */
double CPU_Peano; /*!< time required to establish Peano-Hilbert order */
#ifdef DETAILED_CPU_DOMAIN
double CPU_Domain_findExtend;
double CPU_Domain_determineTopTree;
double CPU_Domain_sumCost;
double CPU_Domain_findSplit;
double CPU_Domain_shiftSplit;
double CPU_Domain_countToGo;
double CPU_Domain_exchange;
#endif
#ifdef DETAILED_CPU_GRAVITY
double CPU_Gravity_TreeWalk1;
double CPU_Gravity_TreeWalk2;
double CPU_Gravity_CommSum1;
double CPU_Gravity_CommSum2;
double CPU_Gravity_Imbalance1;
double CPU_Gravity_Imbalance2;
#endif
#ifdef COOLING
double CPU_Cooling;
#endif
#ifdef DETAILED_CPU
double CPU_Leapfrog;
double CPU_Physics;
double CPU_Residual;
double CPU_Accel;
double CPU_Begrun;
#endif
/* tree code opening criterion */
double ErrTolTheta; /*!< BH tree opening angle */
double ErrTolForceAcc; /*!< parameter for relative opening criterion in tree walk */
/* adjusts accuracy of time-integration */
double ErrTolIntAccuracy; /*!< accuracy tolerance parameter \f$ \eta \f$ for timestep criterion. The
timestep is \f$ \Delta t = \sqrt{\frac{2 \eta eps}{a}} \f$ */
double MinSizeTimestep; /*!< minimum allowed timestep. Normally, the simulation terminates if the
timestep determined by the timestep criteria falls below this limit. */
double MaxSizeTimestep; /*!< maximum allowed timestep */
double MaxRMSDisplacementFac; /*!< this determines a global timestep criterion for cosmological simulations
in comoving coordinates. To this end, the code computes the rms velocity
of all particles, and limits the timestep such that the rms displacement
is a fraction of the mean particle separation (determined from the
particle mass and the cosmological parameters). This parameter specifies
this fraction. */
double CourantFac; /*!< SPH-Courant factor */
/* frequency of tree reconstruction/domain decomposition */
double TreeDomainUpdateFrequency; /*!< controls frequency of domain decompositions */
/* Gravitational and hydrodynamical softening lengths (given in terms of an `equivalent' Plummer softening length).
* Five groups of particles are supported 0="gas", 1="halo", 2="disk", 3="bulge", 4="stars", 5="bndry"
*/
double MinGasHsmlFractional; /*!< minimum allowed SPH smoothing length in units of SPH gravitational softening length */
double MinGasHsml; /*!< minimum allowed SPH smoothing length */
double SofteningGas; /*!< comoving gravitational softening lengths for type 0 */
double SofteningHalo; /*!< comoving gravitational softening lengths for type 1 */
double SofteningDisk; /*!< comoving gravitational softening lengths for type 2 */
double SofteningBulge; /*!< comoving gravitational softening lengths for type 3 */
double SofteningStars; /*!< comoving gravitational softening lengths for type 4 */
double SofteningBndry; /*!< comoving gravitational softening lengths for type 5 */
double SofteningGasMaxPhys; /*!< maximum physical softening length for type 0 */
double SofteningHaloMaxPhys; /*!< maximum physical softening length for type 1 */
double SofteningDiskMaxPhys; /*!< maximum physical softening length for type 2 */
double SofteningBulgeMaxPhys; /*!< maximum physical softening length for type 3 */
double SofteningStarsMaxPhys; /*!< maximum physical softening length for type 4 */
double SofteningBndryMaxPhys; /*!< maximum physical softening length for type 5 */
double SofteningTable[6]; /*!< current (comoving) gravitational softening lengths for each particle type */
double ForceSoftening[6]; /*!< the same, but multiplied by a factor 2.8 - at that scale the force is Newtonian */
double MassTable[6]; /*!< Table with particle masses for particle types with equal mass.
If particle masses are all equal for one type, the corresponding entry in MassTable
is set to this value, allowing the size of the snapshot files to be reduced. */
/* some filenames */
char InitCondFile[MAXLEN_FILENAME]; /*!< filename of initial conditions */
char OutputDir[MAXLEN_FILENAME]; /*!< output directory of the code */
char SnapshotFileBase[MAXLEN_FILENAME]; /*!< basename to construct the names of snapshotf files */
char EnergyFile[MAXLEN_FILENAME]; /*!< name of file with energy statistics */
#ifdef SYSTEMSTATISTICS
char SystemFile[MAXLEN_FILENAME];
#endif
char CpuFile[MAXLEN_FILENAME]; /*!< name of file with cpu-time statistics */
char InfoFile[MAXLEN_FILENAME]; /*!< name of log-file with a list of the timesteps taken */
char LogFile[MAXLEN_FILENAME]; /*!< name of log-file with varied info */
#ifdef SFR
char SfrFile[MAXLEN_FILENAME]; /*!< name of file with sfr records */
#endif
#ifdef CHIMIE
char ChimieFile[MAXLEN_FILENAME]; /*!< name of file with chimie records */
#endif
#ifdef MULTIPHASE
char PhaseFile[MAXLEN_FILENAME]; /*!< name of file with phase records */
char StickyFile[MAXLEN_FILENAME]; /*!< name of file with sticky records */
#endif
#ifdef AGN_ACCRETION
char AccretionFile[MAXLEN_FILENAME]; /*!< name of file with accretion records */
#endif
#ifdef BONDI_ACCRETION
char BondiFile[MAXLEN_FILENAME]; /*!< name of file with bondi records */
#endif
#ifdef BUBBLES
char BubbleFile[MAXLEN_FILENAME]; /*!< name of file with bubble records */
#endif
char TimingsFile[MAXLEN_FILENAME]; /*!< name of file with performance metrics of gravitational tree algorithm */
char RestartFile[MAXLEN_FILENAME]; /*!< basename of restart-files */
char ResubmitCommand[MAXLEN_FILENAME]; /*!< name of script-file that will be executed for automatic restart */
char OutputListFilename[MAXLEN_FILENAME]; /*!< name of file with list of desired output times */
double OutputListTimes[MAXLEN_OUTPUTLIST]; /*!< table with desired output times */
int OutputListLength; /*!< number of output times stored in the table of desired output times */
#ifdef RANDOMSEED_AS_PARAMETER
int RandomSeed; /*!< initial random seed >*/
#endif
}
All; /*!< a container variable for global variables that are equal on all processors */
/*! This structure holds all the information that is
* stored for each particle of the simulation.
*/
extern struct particle_data
{
FLOAT Pos[3]; /*!< particle position at its current time */
FLOAT Mass; /*!< particle mass */
FLOAT Vel[3]; /*!< particle velocity at its current time */
FLOAT GravAccel[3]; /*!< particle acceleration due to gravity */
#ifdef PMGRID
FLOAT GravPM[3]; /*!< particle acceleration due to long-range PM gravity force*/
#endif
#ifdef FORCETEST
FLOAT GravAccelDirect[3]; /*!< particle acceleration when computed with direct summation */
#endif
FLOAT Potential; /*!< gravitational potential */
FLOAT OldAcc; /*!< magnitude of old gravitational force. Used in relative opening criterion */
#ifndef LONGIDS
unsigned int ID; /*!< particle identifier */
#else
unsigned long long ID; /*!< particle identifier */
#endif
int Type; /*!< flags particle type. 0=gas, 1=halo, 2=disk, 3=bulge, 4=stars, 5=bndry */
int Ti_endstep; /*!< marks start of current timestep of particle on integer timeline */
int Ti_begstep; /*!< marks end of current timestep of particle on integer timeline */
#ifdef FLEXSTEPS
int FlexStepGrp; /*!< a random 'offset' on the timeline to create a smooth groouping of particles */
#endif
float GravCost; /*!< weight factor used for balancing the work-load */
#ifdef PSEUDOSYMMETRIC
float AphysOld; /*!< magnitude of acceleration in last timestep. Used to make a first order
prediction of the change of acceleration expected in the future, thereby
allowing to guess whether a decrease/increase of the timestep should occur
in the timestep that is started. */
#endif
#ifdef PARTICLE_FLAG
float Flag;
#endif
#ifdef STELLAR_PROP
unsigned int StPIdx; /*!< index to the corresponding StP particle */
#endif
#ifdef TESSEL
int iT; /*!< index of a triangle to which the point belong to */
int IsDone;
int IsAdded; /*!< if the point has already be added in the tesselation */
int ivPoint; /*!< index of first voronoi point */
int nvPoints; /*!< number of voronoi points */
int iMedian;
int nMedians;
double Volume;
double Density;
- double rSearch; /* radius in which particles must search for ngbs */
+ double rSearch; /*!< radius in which particles must search for ngbs */
+ struct particle_data* Pref; /*!< for a ghost point, points towards the reference point */
#endif
}
*P, /*!< holds particle data on local processor */
#ifdef PY_INTERFACE
*Q,
*DomainPartBufQ, /*!< buffer for particle data used in domain decomposition */
#endif
*DomainPartBuf; /*!< buffer for particle data used in domain decomposition */
/* the following struture holds data that is stored for each SPH particle in addition to the collisionless
* variables.
*/
extern struct sph_particle_data
{
FLOAT Entropy; /*!< current value of entropy (actually entropic function) of particle */
FLOAT Density; /*!< current baryonic mass density of particle */
FLOAT Hsml; /*!< current smoothing length */
FLOAT Left; /*!< lower bound in iterative smoothing length search */
FLOAT Right; /*!< upper bound in iterative smoothing length search */
FLOAT NumNgb; /*!< weighted number of neighbours found */
#ifdef AVOIDNUMNGBPROBLEM
FLOAT OldNumNgb;
#endif
FLOAT Pressure; /*!< current pressure */
FLOAT DtEntropy; /*!< rate of change of entropy */
#ifdef STELLAR_FLUX
FLOAT EnergyFlux; /*!< current value of local energy flux - Sph particles */
#endif
#ifdef AGN_HEATING
FLOAT EgySpecAGNHeat; /*!< current value of specific energy radiated of particle - Sph particles */
FLOAT DtEgySpecAGNHeat; /*!< rate of change of specific radiated energy - Sph particles */
FLOAT DtEntropyAGNHeat;
#endif
#ifdef MULTIPHASE
FLOAT StickyTime;
int StickyFlag;
#ifdef COUNT_COLLISIONS
float StickyCollisionNumber;
#endif
#endif
#ifdef FEEDBACK
FLOAT EgySpecFeedback;
FLOAT DtEgySpecFeedback;
FLOAT EnergySN;
FLOAT EnergySNrem;
FLOAT TimeSN;
FLOAT FeedbackVel[3]; /*!< kick due to feedback force */
#endif
#ifdef FEEDBACK_WIND
FLOAT FeedbackWindVel[3]; /*!< kick due to feedback force */
#endif
FLOAT HydroAccel[3]; /*!< acceleration due to hydrodynamical force */
FLOAT VelPred[3]; /*!< predicted SPH particle velocity at the current time */
FLOAT DivVel; /*!< local velocity divergence */
FLOAT CurlVel; /*!< local velocity curl */
FLOAT Rot[3]; /*!< local velocity curl */
FLOAT DhsmlDensityFactor; /*!< correction factor needed in the equation of motion of the conservative entropy formulation of SPH */
FLOAT MaxSignalVel; /*!< maximum "signal velocity" occuring for this particle */
#ifdef MULTIPHASE
int Phase;
int StickyIndex;
int StickyNgb;
int StickyMaxID;
float StickyMaxFs;
FLOAT StickyNewVel[3];
#endif
#ifdef OUTPUTOPTVAR1
FLOAT OptVar1; /*!< optional variable 1 */
#endif
#ifdef OUTPUTOPTVAR2
FLOAT OptVar2; /*!< optional variable 2 */
#endif
#ifdef COMPUTE_VELOCITY_DISPERSION
FLOAT VelocityDispersion[VELOCITY_DISPERSION_SIZE]; /*!< velocity dispersion */
#endif
#ifdef CHIMIE
FLOAT Metal[NELEMENTS];
FLOAT dMass; /*!< mass variation due to mass transfere */
#ifdef CHIMIE_THERMAL_FEEDBACK
FLOAT DeltaEgySpec;
FLOAT SNIaThermalTime; /*!< flag particles that got energy from SNIa */
FLOAT SNIIThermalTime; /*!< flag particles that got energy from SNII */
double NumberOfSNIa;
double NumberOfSNII;
#endif
#ifdef CHIMIE_KINETIC_FEEDBACK
FLOAT WindTime; /*!< flag particles that belongs to the wind */
unsigned int WindFlag; /*!< flag particles that will be part of the wind */
#endif
#endif /*CHIMIE*/
#ifdef ENTROPYPRED
FLOAT EntropyPred; /*!< predicted entropy at the current time */
#endif
#ifdef ART_CONDUCTIVITY
FLOAT EnergyIntPred;
FLOAT GradEnergyInt[3];
#endif
#ifdef AB_TURB
FLOAT TurbAccel[3];
#endif
#if defined(ART_VISCO_MM)|| defined(ART_VISCO_RO) || defined(ART_VISCO_CD)
double ArtBulkViscConst;
#ifdef ART_VISCO_CD
double DmatCD[3][3];
double TmatCD[3][3];
double DiVelAccurate;
double DiVelTemp;
double ArtBulkViscConstOld;
double R_CD;
FLOAT MaxSignalVelCD;
#endif
#endif
#if PY_INTERFACE
FLOAT Observable;
FLOAT ObsMoment0;
FLOAT ObsMoment1;
FLOAT GradObservable[3];
#endif
}
*SphP, /*!< holds SPH particle data on local processor */
#ifdef PY_INTERFACE
*SphQ,
*DomainSphBufQ, /*!< buffer for SPH particle data in domain decomposition */
#endif
*DomainSphBuf; /*!< buffer for SPH particle data in domain decomposition */
#ifdef STELLAR_PROP
/* the following struture holds data that is stored for each SPH particle in addition to the collisionless
* variables.
*/
extern struct st_particle_data
{
#ifdef CHECK_ID_CORRESPONDENCE
unsigned int ID; /*!< particle identifier (must be the same as P[].ID) only used to check ID correspondance */
#endif
FLOAT FormationTime; /*!< star formation time of particle */
FLOAT InitialMass; /*!< initial stellar mass */
#ifndef LONGIDS
unsigned int IDProj; /*!< id of projenitor particle */
#else
unsigned long long IDProj; /*!< id of projenitor particle */
#endif
FLOAT Metal[NELEMENTS];
FLOAT Density; /*!< current baryonic mass density of particle */
FLOAT Volume; /*!< current volume of particle */
FLOAT Hsml; /*!< current smoothing length */
FLOAT Left; /*!< lower bound in iterative smoothing length search */
FLOAT Right; /*!< upper bound in iterative smoothing length search */
FLOAT NumNgb; /*!< weighted number of neighbours found */
unsigned int PIdx; /*!< index to the corresponding particle */
#ifdef AVOIDNUMNGBPROBLEM
FLOAT OldNumNgb;
#endif
FLOAT DhsmlDensityFactor; /*!< correction factor needed in the equation of motion of the conservative entropy formulation of SPH */
double TotalEjectedGasMass;
double TotalEjectedEltMass[NELEMENTS];
double TotalEjectedEgySpec;
double NumberOfSNIa;
double NumberOfSNII;
#ifdef CHIMIE_KINETIC_FEEDBACK
double NgbMass; /*!< mass of neighbours */
#endif
#ifdef CHIMIE
unsigned int Flag;
#endif
}
*StP, /*!< holds ST particle data on local processor */
*DomainStBuf; /*!< buffer for ST particle data in domain decomposition */
#endif
/* Variables for Tree
*/
extern int MaxNodes; /*!< maximum allowed number of internal nodes */
extern int Numnodestree; /*!< number of (internal) nodes in each tree */
extern struct NODE
{
FLOAT len; /*!< sidelength of treenode */
FLOAT center[3]; /*!< geometrical center of node */
#ifdef ADAPTIVE_GRAVSOFT_FORGAS
FLOAT maxsoft; /*!< hold the maximum gravitational softening of particles in the
node if the ADAPTIVE_GRAVSOFT_FORGAS option is selected */
#endif
#ifdef STELLAR_FLUX
FLOAT starlum ; /*!< star luminosity of node */
#endif
union
{
int suns[8]; /*!< temporary pointers to daughter nodes */
struct
{
FLOAT s[3]; /*!< center of mass of node */
FLOAT mass; /*!< mass of node */
int bitflags; /*!< a bit-field with various information on the node */
int sibling; /*!< this gives the next node in the walk in case the current node can be used */
int nextnode; /*!< this gives the next node in case the current node needs to be opened */
int father; /*!< this gives the parent node of each node (or -1 if we have the root node) */
}
d;
}
u;
}
*Nodes_base, /*!< points to the actual memory allocted for the nodes */
*Nodes; /*!< this is a pointer used to access the nodes which is shifted such that Nodes[All.MaxPart]
gives the first allocated node */
extern int *Nextnode; /*!< gives next node in tree walk */
extern int *Father; /*!< gives parent node in tree */
extern struct extNODE /*!< this structure holds additional tree-node information which is not needed in the actual gravity computation */
{
FLOAT hmax; /*!< maximum SPH smoothing length in node. Only used for gas particles */
FLOAT vs[3]; /*!< center-of-mass velocity */
}
*Extnodes_base, /*!< points to the actual memory allocted for the extended node information */
*Extnodes; /*!< provides shifted access to extended node information, parallel to Nodes/Nodes_base */
/*! Header for the standard file format.
*/
extern struct io_header
{
int npart[6]; /*!< number of particles of each type in this file */
double mass[6]; /*!< mass of particles of each type. If 0, then the masses are explicitly
stored in the mass-block of the snapshot file, otherwise they are omitted */
double time; /*!< time of snapshot file */
double redshift; /*!< redshift of snapshot file */
int flag_sfr; /*!< flags whether the simulation was including star formation */
int flag_feedback; /*!< flags whether feedback was included (obsolete) */
unsigned int npartTotal[6]; /*!< total number of particles of each type in this snapshot. This can be
different from npart if one is dealing with a multi-file snapshot. */
int flag_cooling; /*!< flags whether cooling was included */
int num_files; /*!< number of files in multi-file snapshot */
double BoxSize; /*!< box-size of simulation in case periodic boundaries were used */
double Omega0; /*!< matter density in units of critical density */
double OmegaLambda; /*!< cosmological constant parameter */
double HubbleParam; /*!< Hubble parameter in units of 100 km/sec/Mpc */
int flag_stellarage; /*!< flags whether the file contains formation times of star particles */
int flag_metals; /*!< flags whether the file contains metallicity values for gas and star particles */
unsigned int npartTotalHighWord[6]; /*!< High word of the total number of particles of each type */
int flag_entropy_instead_u; /*!< flags that IC-file contains entropy instead of u */
int flag_chimie_extraheader; /*!< flags that IC-file contains extra-header for chimie */
#ifdef MULTIPHASE
double critical_energy_spec;
#ifdef MESOMACHINE
char fill[38];
#else
char fill[48]; /* use 42 with regor... */
#endif
#else
char fill[56]; /*!< fills to 256 Bytes */
#endif
}
header; /*!< holds header for snapshot files */
#ifdef CHIMIE_EXTRAHEADER
/*! Header for the chimie part.
*/
extern struct io_chimie_extraheader
{
int nelts; /*!< number of chemical element followed */
float SolarAbundances[NELEMENTS];
char labels[256-4-4*(NELEMENTS)];
}
chimie_extraheader;
#endif
#define IO_NBLOCKS 24 /*!< total number of defined information blocks for snapshot files.
Must be equal to the number of entries in "enum iofields" */
enum iofields /*!< this enumeration lists the defined output blocks in snapshot files. Not all of them need to be present. */
{
IO_POS,
IO_VEL,
IO_ID,
IO_MASS,
IO_U,
IO_RHO,
IO_HSML,
IO_POT,
IO_ACCEL,
IO_DTENTR,
IO_TSTP,
IO_ERADSPH,
IO_ERADSTICKY,
IO_ERADFEEDBACK,
IO_ENERGYFLUX,
IO_METALS,
IO_STAR_FORMATIONTIME,
IO_INITIAL_MASS,
IO_STAR_IDPROJ,
IO_STAR_RHO,
IO_STAR_HSML,
IO_STAR_METALS,
IO_OPTVAR1,
IO_OPTVAR2
};
extern char Tab_IO_Labels[IO_NBLOCKS][4]; /*<! This table holds four-byte character tags used for fileformat 2 */
/* global state of system, used for global statistics
*/
extern struct state_of_system
{
double Mass;
double EnergyKin;
double EnergyPot;
double EnergyInt;
#ifdef COOLING
double EnergyRadSph;
#endif
#ifdef AGN_HEATING
double EnergyAGNHeat;
#endif
#ifdef MULTIPHASE
double EnergyRadSticky;
#endif
#ifdef FEEDBACK_WIND
double EnergyFeedbackWind;
#endif
#ifdef BUBBLES
double EnergyBubbles;
#endif
#ifdef CHIMIE_THERMAL_FEEDBACK
double EnergyThermalFeedback;
#endif
#ifdef CHIMIE_KINETIC_FEEDBACK
double EnergyKineticFeedback;
#endif
double EnergyTot;
double Momentum[4];
double AngMomentum[4];
double CenterOfMass[4];
double MassComp[6];
double EnergyKinComp[6];
double EnergyPotComp[6];
double EnergyIntComp[6];
#ifdef COOLING
double EnergyRadSphComp[6];
#endif
#ifdef AGN_HEATING
double EnergyAGNHeatComp[6];
#endif
#ifdef MULTIPHASE
double EnergyRadStickyComp[6];
#endif
#ifdef FEEDBACK_WIND
double EnergyFeedbackWindComp[6];
#endif
#ifdef BUBBLES
double EnergyBubblesComp[6];
#endif
#ifdef CHIMIE_THERMAL_FEEDBACK
double EnergyThermalFeedbackComp[6];
#endif
#ifdef CHIMIE_KINETIC_FEEDBACK
double EnergyKineticFeedbackComp[6];
#endif
double EnergyTotComp[6];
double MomentumComp[6][4];
double AngMomentumComp[6][4];
double CenterOfMassComp[6][4];
}
SysState; /*<! Structure for storing some global statistics about the simulation. */
/*! This structure contains data related to the energy budget.
These values are different for each task. It need to be stored
in the restart flag.
*/
extern struct local_state_of_system
{
double EnergyTest;
double EnergyInt1;
double EnergyInt2;
double EnergyKin1;
double EnergyKin2;
#ifdef COOLING
double RadiatedEnergy;
#endif
#ifdef SFR
double StarEnergyInt;
#ifdef FEEDBACK
double StarEnergyFeedback;
#endif
#endif
#ifdef CHIMIE_THERMAL_FEEDBACK
double EnergyThermalFeedback;
#endif
#ifdef CHIMIE_KINETIC_FEEDBACK
double EnergyKineticFeedback;
#endif
#ifdef MULTIPHASE
double EnergyRadSticky;
#endif
#ifdef FEEDBACK_WIND
double EnergyFeedbackWind;
#endif
}
LocalSysState; /*<! Structure for storing some local statistics about the simulation. */
/* Various structures for communication
*/
extern struct gravdata_in
{
union
{
FLOAT Pos[3];
FLOAT Acc[3];
FLOAT Potential;
}
u;
#if defined(UNEQUALSOFTENINGS) || defined(STELLAR_FLUX)
int Type;
#ifdef ADAPTIVE_GRAVSOFT_FORGAS
FLOAT Soft;
#endif
#endif
#ifdef STELLAR_FLUX
FLOAT EnergyFlux;
#endif
union
{
FLOAT OldAcc;
int Ninteractions;
}
w;
}
*GravDataIn, /*!< holds particle data to be exported to other processors */
*GravDataGet, /*!< holds particle data imported from other processors */
*GravDataResult, /*!< holds the partial results computed for imported particles. Note: We use GravDataResult = GravDataGet, such that the result replaces the imported data */
*GravDataOut; /*!< holds partial results received from other processors. This will overwrite the GravDataIn array */
extern struct gravdata_index
{
int Task;
int Index;
int SortIndex;
}
*GravDataIndexTable; /*!< the particles to be exported are grouped by task-number. This table allows the results to be disentangled again and to be assigned to the correct particle */
extern struct densdata_in
{
FLOAT Pos[3];
FLOAT Vel[3];
FLOAT Hsml;
#ifdef MULTIPHASE
int Phase;
#endif
int Index;
int Task;
#ifdef ART_CONDUCTIVITY
FLOAT EnergyIntPred;
#endif
}
*DensDataIn, /*!< holds particle data for SPH density computation to be exported to other processors */
*DensDataGet; /*!< holds imported particle data for SPH density computation */
extern struct densdata_out
{
FLOAT Rho;
FLOAT Div, Rot[3];
FLOAT DhsmlDensity;
FLOAT Ngb;
#ifdef ART_CONDUCTIVITY
FLOAT GradEnergyInt[3];
#endif
}
*DensDataResult, /*!< stores the locally computed SPH density results for imported particles */
*DensDataPartialResult; /*!< imported partial SPH density results from other processors */
extern struct hydrodata_in
{
FLOAT Pos[3];
FLOAT Vel[3];
FLOAT Hsml;
#ifdef FEEDBACK
FLOAT EnergySN;
#endif
#ifdef MULTIPHASE
int Phase;
FLOAT Entropy;
int StickyFlag;
#endif
FLOAT Mass;
FLOAT Density;
FLOAT Pressure;
FLOAT F1;
FLOAT DhsmlDensityFactor;
int Timestep;
int Task;
int Index;
#ifdef WITH_ID_IN_HYDRA
int ID;
#endif
#ifdef ART_CONDUCTIVITY
FLOAT NormGradEnergyInt;
#endif
#if defined(ART_VISCO_MM)|| defined(ART_VISCO_RO) || defined(ART_VISCO_CD)
double ArtBulkViscConst;
#endif
}
*HydroDataIn, /*!< holds particle data for SPH hydro-force computation to be exported to other processors */
*HydroDataGet; /*!< holds imported particle data for SPH hydro-force computation */
extern struct hydrodata_out
{
FLOAT Acc[3];
FLOAT DtEntropy;
#ifdef FEEDBACK
FLOAT DtEgySpecFeedback;
FLOAT FeedbackAccel[3]; /*!< acceleration due to feedback force */
#endif
FLOAT MaxSignalVel;
#ifdef COMPUTE_VELOCITY_DISPERSION
FLOAT VelocityDispersion[VELOCITY_DISPERSION_SIZE];
#endif
#ifdef MULTIPHASE
FLOAT StickyDVel[3]; /*!< differences in velocities induced by sticky collisions */
#endif
#ifdef OUTPUT_CONDUCTIVITY
FLOAT OptVar2;
#endif
#ifdef ART_VISCO_CD
double DmatCD[3][3];
double TmatCD[3][3];
double R_CD;
FLOAT MaxSignalVelCD;
#endif
}
*HydroDataResult, /*!< stores the locally computed SPH hydro results for imported particles */
*HydroDataPartialResult; /*!< imported partial SPH hydro-force results from other processors */
#ifdef MULTIPHASE
extern struct stickydata_in
{
FLOAT Pos[3];
FLOAT Vel[3];
FLOAT Mass;
FLOAT Hsml;
int ID;
int StickyMaxID;
int StickyNgb;
float StickyMaxFs;
int Task;
int Index;
}
*StickyDataIn, /*!< holds particle data for sticky computation to be exported to other processors */
*StickyDataGet; /*!< holds imported particle data for sticky computation */
extern struct stickydata_out
{
int StickyMaxID;
int StickyNgb;
float StickyMaxFs;
FLOAT StickyNewVel[3];
}
*StickyDataResult, /*!< stores the locally computed sticky results for imported particles */
*StickyDataPartialResult; /*!< imported partial sticky results from other processors */
extern struct Sticky_index
{
int Index;
int CellIndex;
int Flag;
}
*StickyIndex;
#endif
#ifdef CHIMIE
extern struct chimiedata_in
{
FLOAT Pos[3];
FLOAT Vel[3];
#ifndef LONGIDS
unsigned int ID; /*!< particle identifier */
#else
unsigned long long ID; /*!< particle identifier */
#endif
FLOAT Hsml;
#ifdef FEEDBACK
FLOAT EnergySN;
#endif
#ifdef MULTIPHASE
int Phase;
FLOAT Entropy;
int StickyFlag;
#endif
FLOAT Density;
FLOAT Volume;
FLOAT Pressure;
FLOAT F1;
FLOAT DhsmlDensityFactor;
int Timestep;
int Task;
int Index;
#ifdef WITH_ID_IN_HYDRA
int ID;
#endif
double TotalEjectedGasMass;
double TotalEjectedEltMass[NELEMENTS];
double TotalEjectedEgySpec;
double NumberOfSNIa;
double NumberOfSNII;
#ifdef CHIMIE_KINETIC_FEEDBACK
FLOAT NgbMass;
#endif
}
*ChimieDataIn, /*!< holds particle data for Chimie computation to be exported to other processors */
*ChimieDataGet; /*!< holds imported particle data for Chimie computation */
extern struct chimiedata_out
{
FLOAT Acc[3];
FLOAT DtEntropy;
#ifdef FEEDBACK
FLOAT DtEgySpecFeedback;
FLOAT FeedbackAccel[3]; /*!< acceleration due to feedback force */
#endif
FLOAT MaxSignalVel;
#ifdef COMPUTE_VELOCITY_DISPERSION
FLOAT VelocityDispersion[VELOCITY_DISPERSION_SIZE];
#endif
#ifdef MULTIPHASE
FLOAT StickyDVel[3]; /*!< differences in velocities induced by sticky collisions */
#endif
}
*ChimieDataResult, /*!< stores the locally computed Chimie results for imported particles */
*ChimieDataPartialResult; /*!< imported partial Chimie results from other processors */
extern struct starsdensdata_in
{
FLOAT Pos[3];
FLOAT Hsml;
int Index;
int Task;
}
*StarsDensDataIn, /*!< holds particle data for SPH density computation to be exported to other processors */
*StarsDensDataGet; /*!< holds imported particle data for SPH density computation */
extern struct starsdensdata_out
{
FLOAT Rho;
FLOAT Volume;
FLOAT DhsmlDensity;
FLOAT Ngb;
#ifdef CHIMIE_KINETIC_FEEDBACK
FLOAT NgbMass;
#endif
}
*StarsDensDataResult, /*!< stores the locally computed SPH density results for imported particles */
*StarsDensDataPartialResult; /*!< imported partial SPH density results from other processors */
#endif /*CHIMIE*/
#ifdef TESSEL
extern struct ghostdata_in
{
FLOAT Pos[3];
FLOAT rSearch;
int Index;
int Task;
}
*GhostDataIn, /*!< holds particle data for SPH density computation to be exported to other processors */
*GhostDataGet; /*!< holds imported particle data for SPH density computation */
extern struct ghostdata_out
{
FLOAT Value;
}
*GhostDataResult, /*!< stores the locally computed SPH density results for imported particles */
*GhostDataPartialResult; /*!< imported partial SPH density results from other processors */
/* ghost particles */
//extern struct ghost_particle_data
//{
// FLOAT Pos[3]; /*!< particle position at its current time */
// FLOAT Mass; /*!< particle mass */
//}
// *gP;
extern int NumgPart;
#endif /* TESSEL */
#ifdef PY_INTERFACE
extern struct denssphdata_in
{
FLOAT Pos[3];
FLOAT Vel[3];
FLOAT Hsml;
int Index;
int Task;
FLOAT Observable;
}
*DensSphDataIn, /*!< holds particle data for SPH density computation to be exported to other processors */
*DensSphDataGet; /*!< holds imported particle data for SPH density computation */
extern struct denssphdata_out
{
FLOAT Rho;
FLOAT Div, Rot[3];
FLOAT DhsmlDensity;
FLOAT Ngb;
FLOAT GradObservable[3];
}
*DensSphDataResult, /*!< stores the locally computed SPH density results for imported particles */
*DensSphDataPartialResult; /*!< imported partial SPH density results from other processors */
extern struct sphdata_in
{
FLOAT Pos[3];
FLOAT Vel[3];
FLOAT Hsml;
FLOAT ObsMoment0;
FLOAT ObsMoment1;
FLOAT Observable;
int Task;
int Index;
}
*SphDataIn, /*!< holds particle data for SPH hydro-force computation to be exported to other processors */
*SphDataGet; /*!< holds imported particle data for SPH hydro-force computation */
extern struct sphdata_out
{
FLOAT ObsMoment0;
FLOAT ObsMoment1;
FLOAT GradObservable[3];
}
*SphDataResult, /*!< stores the locally computed SPH hydro results for imported particles */
*SphDataPartialResult; /*!< imported partial SPH hydro-force results from other processors */
#endif /*PY_INTERFACE*/
#endif
diff --git a/src/ghost.c b/src/ghost.c
index 387ed8c..1bb4c9b 100644
--- a/src/ghost.c
+++ b/src/ghost.c
@@ -1,559 +1,560 @@
#ifdef TESSEL
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <mpi.h>
#include "allvars.h"
#include "proto.h"
/*! \file ghost.c
* \brief
*
*/
#ifdef PERIODIC
static double boxSize, boxHalf;
#ifdef LONG_X
static double boxSize_X, boxHalf_X;
#else
#define boxSize_X boxSize
#define boxHalf_X boxHalf
#endif
#ifdef LONG_Y
static double boxSize_Y, boxHalf_Y;
#else
#define boxSize_Y boxSize
#define boxHalf_Y boxHalf
#endif
#ifdef LONG_Z
static double boxSize_Z, boxHalf_Z;
#else
#define boxSize_Z boxSize
#define boxHalf_Z boxHalf
#endif
#endif
static int ifirstGhost,nGhostToAdd;
/*! This function computes the local density for each active SPH particle,
* the number of neighbours in the current smoothing radius, and the
* divergence and curl of the velocity field. The pressure is updated as
* well. If a particle with its smoothing region is fully inside the
* local domain, it is not exported to the other processors. The function
* also detects particles that have a number of neighbours outside the
* allowed tolerance range. For these particles, the smoothing length is
* adjusted accordingly, and the density computation is executed again.
* Note that the smoothing length is not allowed to fall below the lower
* bound set by MinGasHsml.
*/
void ghost(int mode)
{
long long ntot, ntotleft;
int *noffset, *nbuffer, *nsend, *nsend_local, *numlist, *ndonelist;
int i, j, n, ndone, npleft, maxfill, source, iter = 0;
int level, ngrp, sendTask, recvTask, place, nexport;
double dt_entr, tstart, tend, tstart_ngb = 0, tend_ngb = 0;
double sumt, sumcomm, timengb, sumtimengb;
double timecomp = 0, timeimbalance = 0, timecommsumm = 0, sumimbalance;
MPI_Status status;
int completeness;
#ifdef DETAILED_CPU_OUTPUT_IN_DENSITY
double *timengblist;
double *timecomplist;
double *timecommsummlist;
double *timeimbalancelist;
#endif
#ifdef DETAILED_CPU
double t0=0,t1=0;
#endif
#ifdef PERIODIC
boxSize = All.BoxSize;
boxHalf = 0.5 * All.BoxSize;
#ifdef LONG_X
boxHalf_X = boxHalf * LONG_X;
boxSize_X = boxSize * LONG_X;
#endif
#ifdef LONG_Y
boxHalf_Y = boxHalf * LONG_Y;
boxSize_Y = boxSize * LONG_Y;
#endif
#ifdef LONG_Z
boxHalf_Z = boxHalf * LONG_Z;
boxSize_Z = boxSize * LONG_Z;
#endif
#endif
#ifdef DETAILED_CPU
if (mode==1)
t0 = second();
#endif
noffset = malloc(sizeof(int) * NTask); /* offsets of bunches in common list */
nbuffer = malloc(sizeof(int) * NTask);
nsend_local = malloc(sizeof(int) * NTask);
nsend = malloc(sizeof(int) * NTask * NTask);
ndonelist = malloc(sizeof(int) * NTask);
for(n = 0, NumPTUpdate = 0; n < N_gas; n++)
{
P[n].IsDone = 0;
P[n].IsAdded = 0;
NumPTUpdate++;
}
numlist = malloc(NTask * sizeof(int) * NTask);
MPI_Allgather(&NumPTUpdate, 1, MPI_INT, numlist, 1, MPI_INT, MPI_COMM_WORLD);
for(i = 0, ntot = 0; i < NTask; i++)
ntot += numlist[i];
free(numlist);
/* we will repeat the whole thing for those particles where we didn't
* find enough neighbours
*/
do
{
i = 0; /* beginn with this index */
ntotleft = ntot; /* particles left for all tasks together */
nGhostToAdd=0;
ifirstGhost=NumPart+NumgPart;
while(ntotleft > 0)
{
for(j = 0; j < NTask; j++)
nsend_local[j] = 0;
/* do local particles and prepare export list */
tstart = second();
for(nexport = 0, ndone = 0; i < N_gas && nexport < All.BunchSizeGhost - NTask; i++)
if(P[i].IsDone == 0)
{
ndone++;
for(j = 0; j < NTask; j++)
Exportflag[j] = 0;
ghost_evaluate(i, 0);
for(j = 0; j < NTask; j++)
{
if(Exportflag[j])
{
GhostDataIn[nexport].Pos[0] = P[i].Pos[0];
GhostDataIn[nexport].Pos[1] = P[i].Pos[1];
GhostDataIn[nexport].Pos[2] = P[i].Pos[2];
GhostDataIn[nexport].Index = i;
GhostDataIn[nexport].rSearch = P[i].rSearch;
nexport++;
nsend_local[j]++;
}
}
}
tend = second();
timecomp += timediff(tstart, tend);
qsort(GhostDataIn, nexport, sizeof(struct ghostdata_in), ghost_compare_key);
for(j = 1, noffset[0] = 0; j < NTask; j++)
noffset[j] = noffset[j - 1] + nsend_local[j - 1];
tstart = second();
MPI_Allgather(nsend_local, NTask, MPI_INT, nsend, NTask, MPI_INT, MPI_COMM_WORLD);
tend = second();
timeimbalance += timediff(tstart, tend);
/* now do the particles that need to be exported */
for(level = 1; level < (1 << PTask); level++)
{
tstart = second();
for(j = 0; j < NTask; j++)
nbuffer[j] = 0;
for(ngrp = level; ngrp < (1 << PTask); ngrp++)
{
maxfill = 0;
for(j = 0; j < NTask; j++)
{
if((j ^ ngrp) < NTask)
if(maxfill < nbuffer[j] + nsend[(j ^ ngrp) * NTask + j])
maxfill = nbuffer[j] + nsend[(j ^ ngrp) * NTask + j];
}
if(maxfill >= All.BunchSizeGhost)
break;
sendTask = ThisTask;
recvTask = ThisTask ^ ngrp;
if(recvTask < NTask)
{
if(nsend[ThisTask * NTask + recvTask] > 0 || nsend[recvTask * NTask + ThisTask] > 0)
{
/* get the particles */
MPI_Sendrecv(&GhostDataIn[noffset[recvTask]],
nsend_local[recvTask] * sizeof(struct ghostdata_in), MPI_BYTE,
recvTask, TAG_DENS_A,
&GhostDataGet[nbuffer[ThisTask]],
nsend[recvTask * NTask + ThisTask] * sizeof(struct ghostdata_in),
MPI_BYTE, recvTask, TAG_DENS_A, MPI_COMM_WORLD, &status);
}
}
for(j = 0; j < NTask; j++)
if((j ^ ngrp) < NTask)
nbuffer[j] += nsend[(j ^ ngrp) * NTask + j];
}
tend = second();
timecommsumm += timediff(tstart, tend);
tstart = second();
for(j = 0; j < nbuffer[ThisTask]; j++)
ghost_evaluate(j, 1);
tend = second();
timecomp += timediff(tstart, tend);
/* do a block to explicitly measure imbalance */
tstart = second();
MPI_Barrier(MPI_COMM_WORLD);
tend = second();
timeimbalance += timediff(tstart, tend);
/* get the result */
tstart = second();
for(j = 0; j < NTask; j++)
nbuffer[j] = 0;
for(ngrp = level; ngrp < (1 << PTask); ngrp++)
{
maxfill = 0;
for(j = 0; j < NTask; j++)
{
if((j ^ ngrp) < NTask)
if(maxfill < nbuffer[j] + nsend[(j ^ ngrp) * NTask + j])
maxfill = nbuffer[j] + nsend[(j ^ ngrp) * NTask + j];
}
if(maxfill >= All.BunchSizeGhost)
break;
sendTask = ThisTask;
recvTask = ThisTask ^ ngrp;
if(recvTask < NTask)
{
if(nsend[ThisTask * NTask + recvTask] > 0 || nsend[recvTask * NTask + ThisTask] > 0)
{
/* send the results */
MPI_Sendrecv(&GhostDataResult[nbuffer[ThisTask]],
nsend[recvTask * NTask + ThisTask] * sizeof(struct ghostdata_out),
MPI_BYTE, recvTask, TAG_DENS_B,
&GhostDataPartialResult[noffset[recvTask]],
nsend_local[recvTask] * sizeof(struct ghostdata_out),
MPI_BYTE, recvTask, TAG_DENS_B, MPI_COMM_WORLD, &status);
/* add the result to the particles */
for(j = 0; j < nsend_local[recvTask]; j++)
{
source = j + noffset[recvTask];
place = GhostDataIn[source].Index;
// this is needed if we want to get some values
//SphP[place].Value = GhostDataPartialResult[source].Value;
}
}
}
for(j = 0; j < NTask; j++)
if((j ^ ngrp) < NTask)
nbuffer[j] += nsend[(j ^ ngrp) * NTask + j];
}
tend = second();
timecommsumm += timediff(tstart, tend);
level = ngrp - 1;
}
MPI_Allgather(&ndone, 1, MPI_INT, ndonelist, 1, MPI_INT, MPI_COMM_WORLD);
for(j = 0; j < NTask; j++)
ntotleft -= ndonelist[j];
}
/* do final operations on results */
tstart = second();
/* add all ghost particles to the tesselation */
AddGhostPoints(ifirstGhost,nGhostToAdd);
//CheckTriangles();
for(i = 0, npleft = 0; i < N_gas; i++)
{
if(!P[i].IsDone)
{
completeness = CheckCompletenessForThisPoint(i);
/* now check whether we had enough neighbours */
if (completeness!=0)
{
/* need to redo this particle */
npleft++;
}
else
P[i].IsDone = 1; /* Mark as inactive */
}
}
tend = second();
timecomp += timediff(tstart, tend);
numlist = malloc(NTask * sizeof(int) * NTask);
MPI_Allgather(&npleft, 1, MPI_INT, numlist, 1, MPI_INT, MPI_COMM_WORLD);
for(i = 0, ntot = 0; i < NTask; i++)
ntot += numlist[i];
free(numlist);
if(ntot > 0)
{
if(iter == 0)
tstart_ngb = second();
iter++;
if(iter > 0 && ThisTask == 0)
{
printf("ngb tessel iteration %d: need to repeat for %d%09d particles.\n", iter,
(int) (ntot / 1000000000), (int) (ntot % 1000000000));
fflush(stdout);
}
if(iter > MAXITER)
{
printf("failed to converge in neighbour iteration in density()\n");
fflush(stdout);
endrun(1155);
}
}
else
tend_ngb = second();
}
while(ntot > 0);
/* mark as active again */
//for(i = 0; i < NumPart; i++)
// P[i].IsDone = 0;
free(ndonelist);
free(nsend);
free(nsend_local);
free(nbuffer);
free(noffset);
/* collect some timing information */
if(iter > 0)
timengb = timediff(tstart_ngb, tend_ngb);
else
timengb = 0;
MPI_Reduce(&timengb, &sumtimengb, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
MPI_Reduce(&timecomp, &sumt, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
MPI_Reduce(&timecommsumm, &sumcomm, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
MPI_Reduce(&timeimbalance, &sumimbalance, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
}
#define NGB_PERIODIC_SIGN_X(x) (xtmp=(x),(xtmp>boxHalf_X)?(-1.0):((xtmp<-boxHalf_X)?(1.0):0))
#define NGB_PERIODIC_SIGN_Y(x) (xtmp=(x),(xtmp>boxHalf_Y)?(-1.0):((xtmp<-boxHalf_Y)?(1.0):0))
#define NGB_PERIODIC_SIGN_Z(x) (xtmp=(x),(xtmp>boxHalf_Z)?(-1.0):((xtmp<-boxHalf_Z)?(1.0):0))
#define ADD_TO_QUADRANT(i,qflag)(1<<i | qflag)
#define QUADRANT_IS_EMPTY(i,qflag)(((1<<i & qflag)==0)?1:0)
#define QUADRANT(x,y,z)( x + 3*y + 9*z )
/*! This function represents the core of the SPH density computation. The
* target particle may either be local, or reside in the communication
* buffer.
*/
void ghost_evaluate(int target, int mode)
{
int j, n, startnode, numngb, numngb_inbox;
double h, h2, hinv, hinv3, hinv4;
double dx, dy, dz, r, r2;
FLOAT *pos;
FLOAT posj[3];
int phase=0;
#ifdef PERIODIC
double xtmp;
double sgnx,sgny,sgnz;
int quadrant;
#endif
if(mode == 0)
{
pos = P[target].Pos;
h = P[target].rSearch;
}
else
{
pos = GhostDataGet[target].Pos;
h = GhostDataGet[target].rSearch;
}
h2 = h * h;
hinv = 1.0 / h;
#ifndef TWODIMS
hinv3 = hinv * hinv * hinv;
#else
hinv3 = hinv * hinv / boxSize_Z;
#endif
hinv4 = hinv3 * hinv;
startnode = All.MaxPart;
numngb = 0;
do
{
numngb_inbox = ngb_treefind_variable_for_tessel(&pos[0], h, phase, &startnode);
for(n = 0; n < numngb_inbox; n++)
{
j = Ngblist[n];
#ifdef PERIODIC
sgnx = NGB_PERIODIC_SIGN_X(P[j].Pos[0]-pos[0]);
sgny = NGB_PERIODIC_SIGN_Y(P[j].Pos[1]-pos[1]);
sgnz = NGB_PERIODIC_SIGN_Z(P[j].Pos[2]-pos[2]);
quadrant = QUADRANT((int)(sgnx+1),(int)(sgny+1),(int)(sgnz+1));
#endif
/* check if a particle is already in a quadrant */
//printf("%g %g %g\n",sgnx+1,sgny+1,sgnz+1);
//printf("quadrant=%d %d\n",quadrant,QUADRANT_IS_EMPTY(quadrant,P[j].IsAdded));
if( QUADRANT_IS_EMPTY(quadrant,P[j].IsAdded) )
{
/* add offset */
posj[0] = P[j].Pos[0] + boxSize_X*sgnx;
posj[1] = P[j].Pos[1] + boxSize_Y*sgny;
posj[2] = P[j].Pos[2] + boxSize_Z*sgnz;
dx = pos[0] - posj[0];
dy = pos[1] - posj[1];
dz = pos[2] - posj[2];
r2 = dx * dx + dy * dy + dz * dz;
if(r2 < h2)
{
/* we can add a point j to the tesselation */
if (NumgPart+NumgPart>=All.MaxPart)
{
printf("NumgPart+NumgPart %d >= MaxPart %d !!!\n",NumgPart+NumgPart,All.MaxPart);
endrun(-123);
}
else
{
//printf("this point is added %g %g %g (%g %g %g)\n",posj[0],posj[1],posj[2],P[j].Pos[0],P[j].Pos[1],P[j].Pos[2]);
P[j].IsAdded = ADD_TO_QUADRANT(quadrant,P[j].IsAdded);
P[NumPart+NumgPart].Pos[0] = posj[0];
P[NumPart+NumgPart].Pos[1] = posj[1];
P[NumPart+NumgPart].Pos[2] = posj[2];
P[NumPart+NumgPart].Mass = P[j].Mass;
P[NumPart+NumgPart].ID = NumPart+NumgPart; /* this is no longer needed */
+ P[NumPart+NumgPart].Pref = &P[j];
//printf("this point is added %d %g %g %g (%g %g %g)\n",NumPart+NumgPart,P[NumPart+NumgPart].Pos[0],P[NumPart+NumgPart].Pos[1],P[NumPart+NumgPart].Pos[2],posj[0],posj[1],posj[2]);
NumgPart++;
nGhostToAdd++;
}
}
}
}
}
while(startnode >= 0);
if(mode == 0)
{
//SphP[target].Value = value;
}
else
{
//GhostDataResult[target].Value = value;
}
}
/*! This routine is a comparison kernel used in a sort routine to group
* particles that are exported to the same processor.
*/
int ghost_compare_key(const void *a, const void *b)
{
if(((struct ghostdata_in *) a)->Task < (((struct ghostdata_in *) b)->Task))
return -1;
if(((struct ghostdata_in *) a)->Task > (((struct ghostdata_in *) b)->Task))
return +1;
return 0;
}
#endif
diff --git a/src/tessel.c b/src/tessel.c
index 1749710..3312602 100644
--- a/src/tessel.c
+++ b/src/tessel.c
@@ -1,2057 +1,2065 @@
#ifdef TESSEL
#include <math.h>
#include <string.h>
#include <stdio.h>
#include <allvars.h>
struct Face /* a list of edges (only for 3d tesselation) */
{
};
struct Triangle
{
struct particle_data Pt1[3];
struct particle_data Pt2[3];
struct particle_data Pt3[3];
};
struct TriangleInList
{
int idx; /* index of current triangle (used for checks) */
struct particle_data* P[3]; /* pointers towards the 3 point */
struct TriangleInList* T[3]; /* pointers towards the 3 triangles */
int idxe[3]; /* index of point in the first triangle, opposite to the common edge */
struct Median* Med[3]; /* pointers towards 3 medians */
};
struct Median
{
double a; /* params for the equation of the segemnt */
double b; /* params for the equation of the segemnt */
double c; /* params for the equation of the segemnt */
struct vPoint *vPs; /* pointer towards starting vPoint of the segment */
struct vPoint *vPe; /* pointer towards stoping vPoint of the segment */
struct particle_data *Pa; /* pointer towards point A*/
struct particle_data *Pb; /* pointer towards point B*/
};
/* a voronoi point */
struct vPoint
{
double Pos[3];
int next;
};
/*
LOCAL VARIABLES
*/
#define MAXNUMTRIANGLES 1000000
int nT=0,numTinStack=0; /* number of triangles in the list */
struct TriangleInList Triangles[MAXNUMTRIANGLES]; /* list of triangles */
struct TriangleInList *TStack[MAXNUMTRIANGLES]; /* index of triangles to check */
struct Median MediansList[3*MAXNUMTRIANGLES][3];
int nvPoints=0; /* number of Voronoi Points */
int nMedians=0; /* number of Medians */
struct vPoint vPoints[5*MAXNUMTRIANGLES];
struct Median Medians[5*MAXNUMTRIANGLES];
double tesselDomainRadius,tesselDomainCenter[3];
struct particle_data Pe[3]; /* edges */
void setup_searching_radius()
{
int i;
for(i=0;i<NumPart;i++)
P[i].rSearch=0.1/10000;
}
void lines_intersections(double a0, double b0, double c0, double a1, double b1, double c1, double *x, double *y)
{
*x = (c1*b0 - c0*b1)/(a0*b1 - a1*b0);
*y = (c1*a0 - c0*a1)/(a1*b0 - a0*b1);
}
/*!
*/
struct Triangle TriangleInList2Triangle(struct TriangleInList Tl)
{
struct Triangle T;
T.Pt1->Pos[0] = Tl.P[0]->Pos[0];
T.Pt1->Pos[1] = Tl.P[0]->Pos[1];
T.Pt2->Pos[0] = Tl.P[1]->Pos[0];
T.Pt2->Pos[1] = Tl.P[1]->Pos[1];
T.Pt3->Pos[0] = Tl.P[2]->Pos[0];
T.Pt3->Pos[1] = Tl.P[2]->Pos[1];
return T;
}
/*! For a set of three points, construct a triangle
*/
struct Triangle MakeTriangleFromPoints(struct particle_data Pt1,struct particle_data Pt2,struct particle_data Pt3)
{
struct Triangle T;
T.Pt1->Pos[0] = Pt1.Pos[0];
T.Pt1->Pos[1] = Pt1.Pos[1];
T.Pt2->Pos[0] = Pt2.Pos[0];
T.Pt2->Pos[1] = Pt2.Pos[1];
T.Pt3->Pos[0] = Pt3.Pos[0];
T.Pt3->Pos[1] = Pt3.Pos[1];
return T;
}
/*! For a set of three points, this function computes the 3 medians.
*/
void TriangleMedians(struct particle_data Pt1,struct particle_data Pt2,struct particle_data Pt3,struct particle_data *Pmm1,struct particle_data *Pmm2,struct particle_data *Pmm3,struct particle_data *Pme1,struct particle_data *Pme2,struct particle_data *Pme3)
{
double ma1,mb1,mc1;
double ma2,mb2,mc2;
double ma3,mb3,mc3;
/* median 0-1 */
ma1 = 2*(Pt2.Pos[0] - Pt1.Pos[0]);
mb1 = 2*(Pt2.Pos[1] - Pt1.Pos[1]);
mc1 = (Pt1.Pos[0]*Pt1.Pos[0]) - (Pt2.Pos[0]*Pt2.Pos[0]) + (Pt1.Pos[1]*Pt1.Pos[1]) - (Pt2.Pos[1]*Pt2.Pos[1]);
/* median 1-2 */
ma2 = 2*(Pt3.Pos[0] - Pt2.Pos[0]);
mb2 = 2*(Pt3.Pos[1] - Pt2.Pos[1]);
mc2 = (Pt2.Pos[0]*Pt2.Pos[0]) - (Pt3.Pos[0]*Pt3.Pos[0]) + (Pt2.Pos[1]*Pt2.Pos[1]) - (Pt3.Pos[1]*Pt3.Pos[1]);
/* median 2-0 */
ma3 = 2*(Pt1.Pos[0] - Pt3.Pos[0]);
mb3 = 2*(Pt1.Pos[1] - Pt3.Pos[1]);
mc3 = (Pt3.Pos[0]*Pt3.Pos[0]) - (Pt1.Pos[0]*Pt1.Pos[0]) + (Pt3.Pos[1]*Pt3.Pos[1]) - (Pt1.Pos[1]*Pt1.Pos[1]);
/* intersection m0-1 -- m1-2 */
Pmm1->Pos[0] = (mc2*mb1 - mc1*mb2)/(ma1*mb2 - ma2*mb1);
Pmm1->Pos[1] = (mc2*ma1 - mc1*ma2)/(ma2*mb1 - ma1*mb2);
/* intersection m1-2 -- m2-0 */
Pmm2->Pos[0] = (mc2*mb1 - mc1*mb2)/(ma1*mb2 - ma2*mb1);
Pmm2->Pos[1] = (mc2*ma1 - mc1*ma2)/(ma2*mb1 - ma1*mb2);
/* intersection m2-0 -- m0-1 */
Pmm3->Pos[0] = (mc2*mb1 - mc1*mb2)/(ma1*mb2 - ma2*mb1);
Pmm3->Pos[1] = (mc2*ma1 - mc1*ma2)/(ma2*mb1 - ma1*mb2);
/* intersection m1-2 -- e1-2 */
Pme1->Pos[0] = 0.5*(Pt1.Pos[0] + Pt2.Pos[0]);
Pme1->Pos[1] = 0.5*(Pt1.Pos[1] + Pt2.Pos[1]);
/* intersection m2-3 -- e3-1 */
Pme2->Pos[0] = 0.5*(Pt2.Pos[0] + Pt3.Pos[0]);
Pme2->Pos[1] = 0.5*(Pt2.Pos[1] + Pt3.Pos[1]);
/* intersection m3-1 -- e1-2 */
Pme3->Pos[0] = 0.5*(Pt3.Pos[0] + Pt1.Pos[0]);
Pme3->Pos[1] = 0.5*(Pt3.Pos[1] + Pt1.Pos[1]);
}
/*! For a set of three points, this function computes their cirum-circle.
* Its radius is return, while the center is return using pointers.
*/
double CircumCircleProperties(struct particle_data Pt1,struct particle_data Pt2,struct particle_data Pt3, double *xc, double *yc)
{
double r;
double x21,x32,y21,y32;
double x12mx22,y12my22,x22mx32,y22my32;
double c1,c2;
x21 = Pt2.Pos[0]-Pt1.Pos[0];
x32 = Pt3.Pos[0]-Pt2.Pos[0];
y21 = Pt2.Pos[1]-Pt1.Pos[1];
y32 = Pt3.Pos[1]-Pt2.Pos[1];
x12mx22 = (Pt1.Pos[0]*Pt1.Pos[0])-(Pt2.Pos[0]*Pt2.Pos[0]);
y12my22 = (Pt1.Pos[1]*Pt1.Pos[1])-(Pt2.Pos[1]*Pt2.Pos[1]);
x22mx32 = (Pt2.Pos[0]*Pt2.Pos[0])-(Pt3.Pos[0]*Pt3.Pos[0]);
y22my32 = (Pt2.Pos[1]*Pt2.Pos[1])-(Pt3.Pos[1]*Pt3.Pos[1]);
c1 = x12mx22 + y12my22;
c2 = x22mx32 + y22my32;
*xc = (y32*c1 - y21*c2)/2.0/( x32*y21 - x21*y32 );
*yc = (x32*c1 - x21*c2)/2.0/( x21*y32 - x32*y21 );
r = sqrt( (Pt1.Pos[0]-*xc)*(Pt1.Pos[0]-*xc) + (Pt1.Pos[1]-*yc)*(Pt1.Pos[1]-*yc) ) ;
return r;
}
/*! For a given triangle T, the routine tells if the point P4
is in the circum circle of the triangle or not.
*/
int InCircumCircle(struct Triangle T,struct particle_data Pt4)
{
double a,b,c;
double d,e,f;
double g,h,i;
double det;
/*
a = T.Pt1->Pos[0] - Pt4.Pos[0];
b = T.Pt1->Pos[1] - Pt4.Pos[1];
c = (T.Pt1->Pos[0]*T.Pt1->Pos[0] - Pt4.Pos[0]*Pt4.Pos[0]) + (T.Pt1->Pos[1]*T.Pt1->Pos[1] - Pt4.Pos[1]*Pt4.Pos[1]);
d = T.Pt2->Pos[0] - Pt4.Pos[0];
e = T.Pt2->Pos[1] - Pt4.Pos[1];
f = (T.Pt2->Pos[0]*T.Pt2->Pos[0] - Pt4.Pos[0]*Pt4.Pos[0]) + (T.Pt2->Pos[1]*T.Pt2->Pos[1] - Pt4.Pos[1]*Pt4.Pos[1]);
g = T.Pt3->Pos[0] - Pt4.Pos[0];
h = T.Pt3->Pos[1] - Pt4.Pos[1];
i = (T.Pt3->Pos[0]*T.Pt3->Pos[0] - Pt4.Pos[0]*Pt4.Pos[0]) + (T.Pt3->Pos[1]*T.Pt3->Pos[1] - Pt4.Pos[1]*Pt4.Pos[1]);
*/
/*
Volker Formula
*/
a = T.Pt2->Pos[0] - T.Pt1->Pos[0];
b = T.Pt2->Pos[1] - T.Pt1->Pos[1];
c = a*a + b*b;
d = T.Pt3->Pos[0] - T.Pt1->Pos[0];
e = T.Pt3->Pos[1] - T.Pt1->Pos[1];
f = d*d + e*e;
g = Pt4.Pos[0] - T.Pt1->Pos[0];
h = Pt4.Pos[1] - T.Pt1->Pos[1];
i = g*g + h*h;
det = a*e*i - a*f*h - b*d*i + b*f*g + c*d*h - c*e*g;
if (det<0)
return 1; /* inside */
else
return 0; /* outside */
}
/*! For a given triangle T, the routine tells if the point P4
lie inside the triangle or not.
*/
int InTriangle(struct Triangle T,struct particle_data Pt4)
{
double c1,c2,c3;
/* here, we use the cross product */
c1 = (T.Pt2->Pos[0]-T.Pt1->Pos[0])*(Pt4.Pos[1]-T.Pt1->Pos[1]) - (T.Pt2->Pos[1]-T.Pt1->Pos[1])*(Pt4.Pos[0]-T.Pt1->Pos[0]);
c2 = (T.Pt3->Pos[0]-T.Pt2->Pos[0])*(Pt4.Pos[1]-T.Pt2->Pos[1]) - (T.Pt3->Pos[1]-T.Pt2->Pos[1])*(Pt4.Pos[0]-T.Pt2->Pos[0]);
c3 = (T.Pt1->Pos[0]-T.Pt3->Pos[0])*(Pt4.Pos[1]-T.Pt3->Pos[1]) - (T.Pt1->Pos[1]-T.Pt3->Pos[1])*(Pt4.Pos[0]-T.Pt3->Pos[0]);
if ( (c1>0) && (c2>0) && (c3>0) ) /* inside */
return 1;
else
return 0;
}
int InTriangleOrOutside(struct Triangle T,struct particle_data Pt4)
{
double c1,c2,c3;
c1 = (T.Pt2->Pos[0]-T.Pt1->Pos[0])*(Pt4.Pos[1]-T.Pt1->Pos[1]) - (T.Pt2->Pos[1]-T.Pt1->Pos[1])*(Pt4.Pos[0]-T.Pt1->Pos[0]);
if (c1<0)
return 2; /* to triangle T[2] */
c2 = (T.Pt3->Pos[0]-T.Pt2->Pos[0])*(Pt4.Pos[1]-T.Pt2->Pos[1]) - (T.Pt3->Pos[1]-T.Pt2->Pos[1])*(Pt4.Pos[0]-T.Pt2->Pos[0]);
if (c2<0)
return 0; /* to triangle T[1] */
c3 = (T.Pt1->Pos[0]-T.Pt3->Pos[0])*(Pt4.Pos[1]-T.Pt3->Pos[1]) - (T.Pt1->Pos[1]-T.Pt3->Pos[1])*(Pt4.Pos[0]-T.Pt3->Pos[0]);
if (c3<0)
return 1; /* to triangle T[0] */
return -1; /* the point is inside */
}
/*! For a given triangle, orient it positively.
*/
struct Triangle OrientTriangle(struct Triangle T)
{
double a,b,c,d;
double det;
struct particle_data Ptsto;
a = T.Pt2->Pos[0] - T.Pt1->Pos[0];
b = T.Pt2->Pos[1] - T.Pt1->Pos[1];
c = T.Pt3->Pos[0] - T.Pt1->Pos[0];
d = T.Pt3->Pos[1] - T.Pt1->Pos[1];
det = (a*d) - (b*c);
if (det<0)
{
Ptsto.Pos[0] = T.Pt1->Pos[0];
Ptsto.Pos[1] = T.Pt1->Pos[1];
T.Pt1->Pos[0] = T.Pt3->Pos[0];
T.Pt1->Pos[1] = T.Pt3->Pos[1];
T.Pt3->Pos[0] = Ptsto.Pos[0];
T.Pt3->Pos[1] = Ptsto.Pos[1];
T = OrientTriangle(T);
}
return T;
}
/*! For a given triangle, orient it positively.
*/
struct TriangleInList OrientTriangleInList(struct TriangleInList T)
{
double a,b,c,d;
double det;
struct particle_data Ptsto;
a = T.P[1]->Pos[0] - T.P[0]->Pos[0];
b = T.P[1]->Pos[1] - T.P[0]->Pos[1];
c = T.P[2]->Pos[0] - T.P[0]->Pos[0];
d = T.P[2]->Pos[1] - T.P[0]->Pos[1];
det = (a*d) - (b*c);
if (det<0)
{
Ptsto.Pos[0] = T.P[0]->Pos[0];
Ptsto.Pos[1] = T.P[0]->Pos[1];
T.P[0]->Pos[0] = T.P[2]->Pos[0];
T.P[0]->Pos[1] = T.P[2]->Pos[1];
T.P[2]->Pos[0] = Ptsto.Pos[0];
T.P[2]->Pos[1] = Ptsto.Pos[1];
T = OrientTriangleInList(T);
}
return T;
}
/*! This function computes the extension of the domain.
* It computes:
* len : max-min
* tesselDomainCenter : min + 0.5*len
* tesselDomainRadius = len*1.5;
*/
void FindTesselExtent()
{
int i,j;
double xmin[3], xmax[3],len;
/* determine local extension */
for(j = 0; j < 3; j++)
{
xmin[j] = MAX_REAL_NUMBER;
xmax[j] = -MAX_REAL_NUMBER;
}
for(i = 0; i < NumPart; i++)
{
for(j = 0; j < 3; j++)
{
if(xmin[j] > P[i].Pos[j])
xmin[j] = P[i].Pos[j];
if(xmax[j] < P[i].Pos[j])
xmax[j] = P[i].Pos[j];
}
}
len = 0;
for(j = 0; j < 3; j++)
{
if(xmax[j] - xmin[j] > len)
len = xmax[j] - xmin[j];
}
for(j = 0; j < 3; j++)
tesselDomainCenter[j] = xmin[j] + len/2.;
tesselDomainRadius = len*2;
printf("tesselDomainRadius = %g\n",tesselDomainRadius);
printf("tesselDomainCenter = (%g %g)\n",tesselDomainCenter[0],tesselDomainCenter[1]);
}
int FindSegmentInTriangle(struct TriangleInList *T,double v,struct particle_data P[3])
{
double v0,v1,v2;
double x0,x1,x2;
double y0,y1,y2;
double f;
double x,y;
int iP;
/* if the triangle as an edge point, do nothing */
if ( (T->P[0]==&Pe[0]) || (T->P[1]==&Pe[0]) || (T->P[2]==&Pe[0]) )
return 0;
/* if the triangle as an edge point, do nothing */
if ( (T->P[0]==&Pe[1]) || (T->P[1]==&Pe[1]) || (T->P[2]==&Pe[1]) )
return 0;
/* if the triangle as an edge point, do nothing */
if ( (T->P[0]==&Pe[2]) || (T->P[1]==&Pe[2]) || (T->P[2]==&Pe[2]) )
return 0;
iP = 0;
v0 = T->P[0]->Mass;
v1 = T->P[1]->Mass;
v2 = T->P[2]->Mass;
//printf("Triangle %d : %g %g %g\n",T->idx,v0,v1,v2);
/* we could also use the sign v-v0 * v-v1 ??? */
if (( ((v>v0)&&(v<v1)) || ((v>v1)&&(v<v0)) )&& (v0 != v1)) /* in 0-1 */
{
x0 = T->P[0]->Pos[0];
y0 = T->P[0]->Pos[1];
x1 = T->P[1]->Pos[0];
y1 = T->P[1]->Pos[1];
f = (v-v0)/(v1-v0);
P[iP].Pos[0] = f*(x1-x0) + x0;
P[iP].Pos[1] = f*(y1-y0) + y0;
iP++;
}
if (( ((v>v1)&&(v<v2)) || ((v>v2)&&(v<v1)) )&& (v1 != v2)) /* in 1-2 */
{
x0 = T->P[1]->Pos[0];
y0 = T->P[1]->Pos[1];
x1 = T->P[2]->Pos[0];
y1 = T->P[2]->Pos[1];
f = (v-v1)/(v2-v1);
P[iP].Pos[0] = f*(x1-x0) + x0;
P[iP].Pos[1] = f*(y1-y0) + y0;
iP++;
}
if (( ((v>v2)&&(v<v0)) || ((v>v0)&&(v<v2)) )&& (v2 != v0)) /* in 2-0 */
{
x0 = T->P[2]->Pos[0];
y0 = T->P[2]->Pos[1];
x1 = T->P[0]->Pos[0];
y1 = T->P[0]->Pos[1];
f = (v-v2)/(v0-v2);
P[iP].Pos[0] = f*(x1-x0) + x0;
P[iP].Pos[1] = f*(y1-y0) + y0;
iP++;
}
return iP;
}
void CheckTriangles(void)
{
int iT;
struct TriangleInList *T,*Te;
for (iT=0;iT<nT;iT++)
{
T = &Triangles[iT];
Te = T->T[0];
if (Te!=NULL)
{
if ((Te->T[0]!=NULL)&&(Te->T[0] == T))
{
}
else
if ((Te->T[1]!=NULL)&&(Te->T[1] == T))
{
}
else
if ((Te->T[2]!=NULL)&&(Te->T[2] == T))
{
}
else
{
printf("Triangle %d does not point towards %d, while T->T2=%d\n",Te->idx,T->idx,T->T[0]->idx);
exit(-1);
}
}
Te = T->T[1];
if (Te!=NULL)
{
if ((Te->T[0]!=NULL)&&(Te->T[0] == T))
{
}
else
if ((Te->T[1]!=NULL)&&(Te->T[1] == T))
{
}
else
if ((Te->T[2]!=NULL)&&(Te->T[2] == T))
{
}
else
{
printf("Triangle %d does not point towards %d, while T->T2=%d\n",Te->idx,T->idx,T->T[1]->idx);
exit(-1);
}
}
Te = T->T[2];
if (Te!=NULL)
{
if ((Te->T[0]!=NULL)&&(Te->T[0] == T))
{
}
else
if ((Te->T[1]!=NULL)&&(Te->T[1] == T))
{
}
else
if ((Te->T[2]!=NULL)&&(Te->T[2] == T))
{
}
else
{
printf("Triangle %d does not point towards %d, while T->T2=%d\n",Te->idx,T->idx,T->T[2]->idx);
exit(-1);
}
}
}
}
void CheckPoints(int n)
{
int i,p,iT;
for (i=0;i<n;i++)
{
/* find p */
iT = P[i].iT;
p=-1;
if (Triangles[iT].P[0]==&P[i])
p=0;
else
if (Triangles[iT].P[1]==&P[i])
p=1;
else
if (Triangles[iT].P[2]==&P[i])
p=2;
// printf(" ---> i=%d (id=%d) iT=%d \n",i,P[i].ID,iT);
// printf(" P.ID=%d%09d\n", (int) (Triangles[iT].P[0]->ID / 1000000000), (int) (Triangles[iT].P[0]->ID % 1000000000));
// printf(" P.ID=%d%09d\n", (int) (Triangles[iT].P[1]->ID / 1000000000), (int) (Triangles[iT].P[1]->ID % 1000000000));
// printf(" P.ID=%d%09d\n", (int) (Triangles[iT].P[2]->ID / 1000000000), (int) (Triangles[iT].P[2]->ID % 1000000000));
if (p==-1)
endrun(-123555);
//printf("%d ok\n",i);
}
}
/*! Flip two triangles.
Te = T.T[i]
*/
void FlipTriangle(int i,struct TriangleInList *T,struct TriangleInList *Te,struct TriangleInList *T1,struct TriangleInList *T2)
{
struct TriangleInList Ts1,Ts2;
int i0,i1,i2;
int j0,j1,j2;
int j;
Ts1 = *T; /* save the content of the pointed triangle */
Ts2 = *Te; /* save the content of the pointed triangle */
j = T->idxe[i]; /* index of point opposite to i */
i0= i;
i1= (i+1) % 3;
i2= (i+2) % 3;
j0= j;
j1= (j+1) % 3;
j2= (j+2) % 3;
/* triangle 1 */
T1->P[0] = Ts1.P[i0];
T1->P[1] = Ts1.P[i1];
T1->P[2] = Ts2.P[j0];
T1->T[0] = Ts2.T[j1];
T1->T[1] = T2;
T1->T[2] = Ts1.T[i2];
T1->idxe[0] = Ts2.idxe[j1];
T1->idxe[1] = 1;
T1->idxe[2] = Ts1.idxe[i2];
/* triangle 2 */
T2->P[0] = Ts2.P[j0];
T2->P[1] = Ts2.P[j1];
T2->P[2] = Ts1.P[i0];
/* restore links with points */
// if (Ts1.P[i0]->iT==Ts1.idx)
// Ts1.P[i0]->iT=T1->idx;
// if (Ts1.P[i1]->iT==Ts1.idx)
// Ts1.P[i1]->iT=T1->idx;
if (Ts1.P[i2]->iT==Ts1.idx)
Ts1.P[i2]->iT=T2->idx;
/* if (Ts2.P[j0]->iT==Ts2.idx)
Ts2.P[j0]->iT=T2->idx;
if (Ts2.P[j1]->iT==Ts2.idx)
Ts2.P[j1]->iT=T2->idx; */
if (Ts2.P[j2]->iT==Ts2.idx)
Ts2.P[j2]->iT=T1->idx;
T2->T[0] = Ts1.T[i1];
T2->T[1] = T1;
T2->T[2] = Ts2.T[j2];
T2->idxe[0] = Ts1.idxe[i1];
T2->idxe[1] = 1;
T2->idxe[2] = Ts2.idxe[j2];
/* restore links with adjacents triangles */
if (Ts1.T[i1]!=NULL)
{
Ts1.T[i1]->T[ Ts1.idxe[i1] ] = T2;
Ts1.T[i1]->idxe[ Ts1.idxe[i1] ] = 0;
}
if (Ts1.T[i2] !=NULL)
{
Ts1.T[i2]->T[ Ts1.idxe[i2] ] = T1;
Ts1.T[i2]->idxe[ Ts1.idxe[i2] ] = 2;
}
if (Ts2.T[j1] !=NULL)
{
Ts2.T[j1]->T[ Ts2.idxe[j1] ] = T1;
Ts2.T[j1]->idxe[ Ts2.idxe[j1] ] = 0;
}
if (Ts2.T[j2] !=NULL)
{
Ts2.T[j2]->T[ Ts2.idxe[j2] ] = T2;
Ts2.T[j2]->idxe[ Ts2.idxe[j2] ] = 2;
}
}
void DoTrianglesInStack(void)
{
struct TriangleInList *T,*Te,*T1,*T2,*Tee;
struct TriangleInList Ts1,Ts2;
struct particle_data Pt;
int istack;
int idx1,idx2;
int i;
istack=0;
while(numTinStack>0)
{
int insphere=0;
T = TStack[istack];
//printf(" DoInStack T=%d (istack=%d, numTinStack=%d)\n",T->idx,istack,numTinStack);
/* find the opposite point of the 3 adjacent triangles */
/*******************/
/* triangle 1 */
/*******************/
i = 0;
Te = T->T[i];
if (Te!=NULL)
{
/* index of opposite point */
Pt = *Te->P[T->idxe[i]];
insphere = InCircumCircle(TriangleInList2Triangle(*T),Pt);
if (insphere)
{
//printf("insphere (1)... %g %g %g in T=%d\n",Pt.Pos[0],Pt.Pos[1],Pt.Pos[2],T->idx);
/* index of the new triangles */
idx1 = T->idx;
idx2 = Te->idx;
T1 = &Triangles[idx1];
T2 = &Triangles[idx2];
FlipTriangle(i,T,Te,T1,T2);
/* add triangles in stack */
if (numTinStack+1>MAXNUMTRIANGLES)
{
printf("\nNo more memory !\n");
printf("numTinStack+1=%d > MAXNUMTRIANGLES=%d\n",numTinStack+1,MAXNUMTRIANGLES);
printf("You should increase MAXNUMTRIANGLES\n\n");
exit(-1);
}
TStack[istack ] = T1;
TStack[istack+numTinStack] = T2;
numTinStack++;
continue;
}
}
/*******************/
/* triangle 2 */
/*******************/
i = 1;
Te = T->T[i];
if (Te!=NULL)
{
/* index of opposite point */
Pt = *Te->P[T->idxe[i]];
insphere = InCircumCircle(TriangleInList2Triangle(*T),Pt);
if (insphere)
{
//printf("insphere (2)... %g %g %g in T=%d\n",Pt.Pos[0],Pt.Pos[1],Pt.Pos[2],T->idx);
/* index of the new triangles */
idx1 = T->idx;
idx2 = Te->idx;
T1 = &Triangles[idx1];
T2 = &Triangles[idx2];
FlipTriangle(i,T,Te,T1,T2);
/* add triangles in stack */
if (numTinStack+1>MAXNUMTRIANGLES)
{
printf("\nNo more memory !\n");
printf("numTinStack+1=%d > MAXNUMTRIANGLES=%d\n",numTinStack+1,MAXNUMTRIANGLES);
printf("You should increase MAXNUMTRIANGLES\n\n");
exit(-1);
}
TStack[istack ] = T1;
TStack[istack+numTinStack] = T2;
numTinStack++;
continue;
}
}
/*******************/
/* triangle 3 */
/*******************/
i = 2;
Te = T->T[i];
if (Te!=NULL)
{
/* index of opposite point */
Pt = *Te->P[T->idxe[i]];
insphere = InCircumCircle(TriangleInList2Triangle(*T),Pt);
if (insphere)
{
//printf("insphere (3)... %g %g %g in T=%d\n",Pt.Pos[0],Pt.Pos[1],Pt.Pos[2],T->idx);
/* index of the new triangles */
idx1 = T->idx;
idx2 = Te->idx;
T1 = &Triangles[idx1];
T2 = &Triangles[idx2];
FlipTriangle(i,T,Te,T1,T2);
/* add triangles in stack */
if (numTinStack+1>MAXNUMTRIANGLES)
{
printf("\nNo more memory !\n");
printf("numTinStack+1=%d > MAXNUMTRIANGLES=%d\n",numTinStack+1,MAXNUMTRIANGLES);
printf("You should increase MAXNUMTRIANGLES\n\n");
exit(-1);
}
TStack[istack ] = T1;
TStack[istack+numTinStack] = T2;
numTinStack++;
continue;
}
}
numTinStack--;
istack++;
//printf("one triangle less...(istack=%d numTinStack=%d)\n",istack,numTinStack);
}
}
void Check(void)
{
int iT;
printf("===========================\n");
for(iT=0;iT<nT;iT++)
{
printf("* T %d\n",Triangles[iT].idx);
printf("pt1 %g %g %g\n",Triangles[iT].P[0]->Pos[0],Triangles[iT].P[0]->Pos[1],Triangles[iT].P[0]->Pos[2]);
printf("pt2 %g %g %g\n",Triangles[iT].P[1]->Pos[0],Triangles[iT].P[1]->Pos[1],Triangles[iT].P[1]->Pos[2]);
printf("pt3 %g %g %g\n",Triangles[iT].P[2]->Pos[0],Triangles[iT].P[2]->Pos[1],Triangles[iT].P[2]->Pos[2]);
if (Triangles[iT].T[0]!=NULL)
printf("T1 %d\n",Triangles[iT].T[0]->idx);
else
printf("T1 x\n");
if (Triangles[iT].T[1]!=NULL)
printf("T2 %d\n",Triangles[iT].T[1]->idx);
else
printf("T2 x\n");
if (Triangles[iT].T[2]!=NULL)
printf("T3 %d\n",Triangles[iT].T[2]->idx);
else
printf("T3 x\n");
}
printf("===========================\n");
}
/*! Split a triangle in 3, using the point P inside it.
Update the global list.
*/
void SplitTriangle(struct TriangleInList *pT,struct particle_data *Pt)
{
struct TriangleInList T,*T0,*T1,*T2,*Te;
int idx,idx0,idx1,idx2;
T = *pT; /* save the content of the pointed triangle */
idx = T.idx;
/* index of the new triangles */
idx0 = idx;
idx1 = nT;
idx2 = nT+1;
/* increment counter */
nT=nT+2;
/* check memory */
if (nT>MAXNUMTRIANGLES)
{
printf("\nNo more memory !\n");
printf("nT=%d > MAXNUMTRIANGLES=%d\n",nT,MAXNUMTRIANGLES);
printf("You should increase MAXNUMTRIANGLES\n\n");
exit(-1);
}
/* create pointers towards the triangles */
T0 = &Triangles[idx0];
T1 = &Triangles[idx1];
T2 = &Triangles[idx2];
/* first */
T0->idx = idx0;
T0->P[0] = T.P[0];
T0->P[1] = T.P[1];
T0->P[2] = Pt;
/* second */
T1->idx = idx1;
T1->P[0] = T.P[1];
T1->P[1] = T.P[2];
T1->P[2] = Pt;
/* third */
T2->idx = idx2;
T2->P[0] = T.P[2];
T2->P[1] = T.P[0];
T2->P[2] = Pt;
/* link points with triangles */
/* the new point Pt belongs to the
3 new triangle, link it to the first */
Pt->iT = idx0;
/* if P[0] points towards the old Triangle, link it with T0 */
if (T.P[0]->iT==T.idx)
{
T.P[0]->iT=T0->idx;
}
/* if P[1] points towards the old Triangle, link it with T1 */
if (T.P[1]->iT==T.idx)
{
T.P[1]->iT=T1->idx;
}
/* if P[2] points towards the old Triangle, link it with T2 */
if (T.P[2]->iT==T.idx)
{
T.P[2]->iT=T2->idx;
}
/* add adjacents */
T0->T[0] = T1;
T0->T[1] = T2;
T0->T[2] = T.T[2];
T1->T[0] = T2;
T1->T[1] = T0;
T1->T[2] = T.T[0];
T2->T[0] = T0;
T2->T[1] = T1;
T2->T[2] = T.T[1];
/* add ext point */
T0->idxe[0] = 1;
T0->idxe[1] = 0;
T0->idxe[2] = T.idxe[2];
T1->idxe[0] = 1;
T1->idxe[1] = 0;
T1->idxe[2] = T.idxe[0];
T2->idxe[0] = 1;
T2->idxe[1] = 0;
T2->idxe[2] = T.idxe[1];
/* restore links with adgacents triangles */
Te = T0->T[2];
if (Te!=NULL)
{
Te->T[ T0->idxe[2]] = T0;
Te->idxe[T0->idxe[2]] = 2;
}
Te = T1->T[2];
if (Te!=NULL)
{
Te->T[ T1->idxe[2]] = T1;
Te->idxe[T1->idxe[2]] = 2;
}
Te = T2->T[2];
if (Te!=NULL)
{
Te->T[ T2->idxe[2]] = T2;
Te->idxe[T2->idxe[2]] = 2;
}
/* add the new triangles in the stack */
TStack[numTinStack] = T0;
numTinStack++;
TStack[numTinStack] = T1;
numTinStack++;
TStack[numTinStack] = T2;
numTinStack++;
//printf("--> add in stack %d %d %d\n",T0->idx,T1->idx,T2->idx);
}
int FindTriangle(struct particle_data *Pt)
{
int iT;
/* find triangle containing the point */
for(iT=0;iT<nT;iT++) /* loop over all triangles */
{
if (InTriangle(TriangleInList2Triangle( Triangles[iT] ),*Pt))
break;
}
return iT;
}
int NewFindTriangle(struct particle_data *Pt)
{
int iT;
struct TriangleInList *T;
int e;
iT = 0; /* star with first triangle */
T = &Triangles[iT];
while (1)
{
/* test position of the point relative to the triangle */
e = InTriangleOrOutside(TriangleInList2Triangle( *T ),*Pt);
//printf("T=%d e=%d Te=%d\n",T->idx,e,T->T[e]->idx);
if (e==-1) /* the point is inside */
break;
T = T->T[e];
if (T==NULL)
{
printf("point lie outside the limits.\n");
exit(-1);
}
}
//printf("done with find triangle (T=%d)\n",T->idx);
return T->idx;
}
/*! Add a new point in the tesselation
*/
void AddPoint(struct particle_data *Pt)
{
int iT;
/* find the triangle that contains the point P */
//iT= FindTriangle(Pt);
iT= NewFindTriangle(Pt);
/* create the new triangles */
SplitTriangle(&Triangles[iT],Pt);
/* test the new triangles and divide and modify if necessary */
DoTrianglesInStack();
/* check */
//CheckTriangles();
}
void AddGhostPoints(int istart,int nadd)
{
int i;
if (ThisTask==0)
printf("%d new ghost points\n",nadd);
/* loop over all ghostpoints */
for (i=istart;i<istart+nadd;i++)
{
AddPoint(&P[i]);
}
}
/*! Add gost points to the tesselation
*/
void ExtendWithGhostPoints()
{
size_t bytes;
/* all particles try to find ngbs particles that lie in another proc */
if (ThisTask==0)
printf("ExtendWithGhostPoints...\n");
All.MaxgPart=All.MaxPart;
NumgPart=0;
/* if(!(gP = malloc(bytes = All.MaxgPart * sizeof(struct ghost_particle_data))))
{
printf("failed to allocate memory for `gP' (%g MB).\n", bytes / (1024.0 * 1024.0));
endrun(1);
} */
ghost();
/* free(gP);*/
if (ThisTask==0)
printf("ExtendWithGhostPoints done.\n");
}
/*! Compute all medians properties (a,b,c)
* For each triangle, for each edge, the function computes the
* median properties which is stored in MediansList
*/
void ComputeMediansProperties()
{
int iT;
/* loop over all triangles */
for(iT=0;iT<nT;iT++)
{
struct particle_data Pt0,Pt1,Pt2;
Pt0.Pos[0] = Triangles[iT].P[0]->Pos[0];
Pt0.Pos[1] = Triangles[iT].P[0]->Pos[1];
Pt1.Pos[0] = Triangles[iT].P[1]->Pos[0];
Pt1.Pos[1] = Triangles[iT].P[1]->Pos[1];
Pt2.Pos[0] = Triangles[iT].P[2]->Pos[0];
Pt2.Pos[1] = Triangles[iT].P[2]->Pos[1];
/* median 0-1 */
MediansList[iT][2].a = 2*(Pt1.Pos[0] - Pt0.Pos[0]);
MediansList[iT][2].b = 2*(Pt1.Pos[1] - Pt0.Pos[1]);
MediansList[iT][2].c = (Pt0.Pos[0]*Pt0.Pos[0]) - (Pt1.Pos[0]*Pt1.Pos[0]) + (Pt0.Pos[1]*Pt0.Pos[1]) - (Pt1.Pos[1]*Pt1.Pos[1]);
/* median 1-2 */
MediansList[iT][0].a = 2*(Pt2.Pos[0] - Pt1.Pos[0]);
MediansList[iT][0].b = 2*(Pt2.Pos[1] - Pt1.Pos[1]);
MediansList[iT][0].c = (Pt1.Pos[0]*Pt1.Pos[0]) - (Pt2.Pos[0]*Pt2.Pos[0]) + (Pt1.Pos[1]*Pt1.Pos[1]) - (Pt2.Pos[1]*Pt2.Pos[1]);
/* median 2-0 */
MediansList[iT][1].a = 2*(Pt0.Pos[0] - Pt2.Pos[0]);
MediansList[iT][1].b = 2*(Pt0.Pos[1] - Pt2.Pos[1]);
MediansList[iT][1].c = (Pt2.Pos[0]*Pt2.Pos[0]) - (Pt0.Pos[0]*Pt0.Pos[0]) + (Pt2.Pos[1]*Pt2.Pos[1]) - (Pt0.Pos[1]*Pt0.Pos[1]);
/* link The triangle with the MediansList */
Triangles[iT].Med[0] = &MediansList[iT][0]; /* median 1-2 */
Triangles[iT].Med[1] = &MediansList[iT][1]; /* median 2-0 */
Triangles[iT].Med[2] = &MediansList[iT][2]; /* median 0-1 */
}
}
/*! Compute the intersetions of medians around a point of index p (index of the point in the triangle T)
*
*/
void ComputeMediansAroundPoint(struct TriangleInList *Tstart,int iPstart)
{
/*
Tstart : pointer to first triangle
iPstart : index of master point relative to triangle Tstart
if p = 0:
T1 = T0->T[iTn]; pn=1
if p = 1:
T1 = T0->T[iTn]; pn=2
if p = 0:
T1 = T0->T[iTn]; pn=3
iTn = (p+1) % 3;
*/
double x,y;
struct TriangleInList *T0,*T1;
int iP0,iP1;
int iT1;
struct particle_data *initialPoint;
int iM0,iM1;
int next_vPoint=-1; /* index towards next voronoi point */
int next_Median=-1; /* index towards next median */
int number_of_vPoints=0;
int number_of_Medians=0;
int next;
T0 = Tstart;
iP0 = iPstart;
initialPoint = T0->P[iP0];
//printf("\n--> rotating around T=%d p=%d\n",T0->idx,iP0);
/* rotate around the point */
while (1)
{
/* next triangle */
iT1= (iP0+1) % 3;
T1 = T0->T[iT1];
if (T1==NULL)
{
//printf("reach an edge\n");
T0->P[iP0]->IsDone=2;
//printf("%g %g\n",T0->P[iP0]->Pos[0],T0->P[iP0]->Pos[1]);
return;
}
//printf(" next triangle = %d\n",T1->idx);
/* index of point in the triangle */
iP1 = T0->idxe[iT1]; /* index of point opposite to iTn */
iP1 = (iP1+1) % 3; /* next index of point opposite to iTn */
//printf(" initial point=%g %g current point =%g %g iP1=%d\n",initialPoint->Pos[0],initialPoint->Pos[1],T1->P[iP1]->Pos[0],T1->P[iP1]->Pos[1],iP1);
/* check */
if (initialPoint!=T1->P[iP1])
{
printf(" problem : initial point=%g %g current point =%g %g iP1=%d\n",initialPoint->Pos[0],initialPoint->Pos[1],T1->P[iP1]->Pos[0],T1->P[iP1]->Pos[1],iP1);
exit(-1);
}
/* compute the intersection of the two medians */
iM0 = (iP0+1) % 3;
iM1 = (iP1+1) % 3;
lines_intersections(T0->Med[iM0]->a,T0->Med[iM0]->b,T0->Med[iM0]->c,T1->Med[iM1]->a,T1->Med[iM1]->b,T1->Med[iM1]->c,&x,&y);
/* create a new vPoint and put it to the vPoints list */
vPoints[nvPoints].Pos[0] = x;
vPoints[nvPoints].Pos[1] = y;
vPoints[nvPoints].next = next_vPoint;
/* end point for T0 */
T0->Med[iM0]->vPe = &vPoints[nvPoints];
/* here, we could add the point to T0->Med[(iM0+2) % 3] as Ps */
/* start point for T0 */
T1->Med[iM1]->vPs = &vPoints[nvPoints];
/* here, we could add the point to T1->Med[(iM1+2) % 3] as Ps */
/* increment vPoints */
next_vPoint=nvPoints;
nvPoints++;
number_of_vPoints++;
if (T1==Tstart) /* end of loop */
{
//printf(" end of loop\n");
initialPoint->ivPoint = next_vPoint;
initialPoint->nvPoints = number_of_vPoints;
/* create the median list */
/* first vPoint */
// next = initialPoint->ivPoint;
//
// for (j = 0; j < initialPoint->nvPoints; j++)
// {
// Medians[nMedians].vPe = vPoints[prev];
// Medians[nMedians].vPs = vPoints[next];
// next = vPoints[next].next;
// }
break;
}
T0 = T1;
iP0 = iP1;
}
}
/*! Compute all medians intersections and define Ps and Pe
* For each triangle, compute the medians
*/
// void oldComputeMediansIntersections()
// {
// int i,p,iT;
//
// for (i=0;i<NumPart+NumgPart;i++)
// P[i].IsDone = 0;
//
// /* loop over all triangles */
// for(iT=0;iT<nT;iT++)
// {
// /* loop over points in triangle */
// for(p=0;p<3;p++)
// {
// if (!(Triangles[iT].P[p]->IsDone))
// {
// //printf("in Triangle T %d do point %d\n",iT,p);
// Triangles[iT].P[p]->IsDone = 1;
// ComputeMediansAroundPoint(&Triangles[iT],p);
// }
// }
// }
// }
void ComputeMediansIntersections()
{
int i,p,iT;
for (i=0;i<NumPart+NumgPart;i++)
{
/* find p */
iT = P[i].iT;
p=-1;
if (Triangles[iT].P[0]==&P[i])
p=0;
else
if (Triangles[iT].P[1]==&P[i])
p=1;
else
if (Triangles[iT].P[2]==&P[i])
p=2;
if (p==-1)
endrun(-12354);
ComputeMediansAroundPoint(&Triangles[P[i].iT],p);
}
// /* do the same for edges */
// for (i=0;i<3;i++)
// {
// /* find p */
// iT = Pe[i].iT;
// p=-1;
// if (Triangles[iT].P[0]->ID==Pe[i].ID)
// p=0;
// else
// if (Triangles[iT].P[1]->ID==Pe[i].ID)
// p=1;
// else
// if (Triangles[iT].P[2]->ID==Pe[i].ID)
// p=2;
//
//
// if (p==-1)
// endrun(-12354);
//
// ComputeMediansAroundPoint(&Triangles[Pe[i].iT],p);
// }
}
/*! Compute the density for all particles
*/
int ComputeDensity()
{
int i,j;
int next; /* next voronoi point */
int np;
double x0,y0,x1,y1;
double area;
- for(i=0;i<NumPart+NumgPart;i++)
+ /* do the normal points first */
+ for(i=0;i<NumPart;i++)
{
next = P[i].ivPoint;
np = P[i].nvPoints;
x0 = 0; /* this ensure the first loop to give 0 */
y0 = 0; /* this ensure the first loop to give 0 */
area = 0;
for (j = 0; j < np; j++)
{
x1 = vPoints[next].Pos[0];
y1 = vPoints[next].Pos[1];
area = area + (x0*y1 - x1*y0);
x0 = x1;
y0 = y1;
next = vPoints[next].next;
}
/* connect the last with the first */
next = P[i].ivPoint;
x1 = vPoints[next].Pos[0];
y1 = vPoints[next].Pos[1];
area = area + (x0*y1 - x1*y0);
/* */
area = 0.5*fabs(area);
P[i].Volume=area;
P[i].Density=P[i].Mass/area;
}
+ /* do the ghost points now */
+ for(i=NumPart;i<NumPart+NumgPart;i++)
+ {
+ P[i].Volume = P[i].Pref->Volume;
+ P[i].Density= P[i].Pref->Density;
+ }
+
return 0;
}
/*! Compute the density for all particles
*/
int CheckCompletenessForThisPoint(int i)
{
/* loop over the triangles which contains the point */
/* first triangle */
int iT,iTstart;
int p,iP0,iP1,iPnext0;
struct TriangleInList *T,*Tnext,*Tstart;
double rc,xc,yc,rcmax;
iTstart = P[i].iT;
Tstart=&Triangles[iTstart];
/* find index in the triangle */
p=-1;
if (Tstart->P[0]==&P[i])
p=0;
else
if (Tstart->P[1]==&P[i])
p=1;
else
if (Tstart->P[2]==&P[i])
p=2;
if (p==-1)
endrun(-312354);
iT = iTstart; /* index of triangle */
iP0 = p; /* index of the ref point "0" in the triangle */
T = Tstart;
/* visit all triangles around the point */
rcmax=0;
while (1)
{
/* compute circumcircle properties of the current tirangle T */
rc = CircumCircleProperties(*T->P[0],*T->P[1],*T->P[2], &xc, &yc);
if (rc>rcmax)
rcmax=rc;
iP1 = (iP0+1) % 3 ; /* point 1 in the current triangle */
Tnext = T->T[iP1]; /* triangle opposit to iP1 */
if (Tnext==NULL) /* we reach an edge */
{
printf("mmmh... why are we here ?");
endrun(-31235445);
return -1; /* need to increase the size for this point */
}
/* now, find the index of the ref. point in the next triangle */
iPnext0 = (T->idxe[iP1] + 1) % 3 ;
/* check */
if (&P[i]!=Tnext->P[iPnext0])
//if (P[i].ID!=Tnext->P[iPnext0]->ID)
{
printf(" problem : initial point %d (id=%d x=%g y=%g z=%g) next-> (id=%d x=%g y=%g z=%g)\n",i,P[i].ID,P[i].Pos[0],P[i].Pos[1],P[i].Pos[2],Tnext->P[iPnext0]->ID,Tnext->P[i]->Pos[0],Tnext->P[i]->Pos[1],Tnext->P[i]->Pos[2]);
endrun(-3123544);
}
/* end of loop */
if (Tnext==Tstart)
break;
iP0 = iPnext0;
T = Tnext;
}
if (P[i].rSearch>2*rcmax)
return 0; /* ok, no need to */
else
{
/* need to increase hi */
//printf("need to increase hi for i=%d %g %g \n",i,P[i].rSearch,2*rcmax);
P[i].rSearch *= 1.5;
return 1;
}
}
/*! Construct the Delaunay tesselation
*/
int ConstructDelaunay()
{
int i,j;
if (ThisTask==0)
printf("Start ConstructDelaunay...\n");
/* find domain extent */
FindTesselExtent();
/* set edges Pe, the 3 points are in an equilateral triangle around all particles */
for (j=0;j<3;j++)
{
Pe[j].Pos[0] = tesselDomainCenter[0] + tesselDomainRadius * cos(2./3.*PI*j);
Pe[j].Pos[1] = tesselDomainCenter[1] + tesselDomainRadius * sin(2./3.*PI*j);
Pe[j].Pos[2] = 0;
Pe[j].Mass = 0;
}
/* Triangle list */
Triangles[0].idx = 0;
Triangles[0].P[0] = &Pe[0];
Triangles[0].P[1] = &Pe[1];
Triangles[0].P[2] = &Pe[2];
Triangles[0].T[0] = NULL;
Triangles[0].T[1] = NULL;
Triangles[0].T[2] = NULL;
Triangles[0].idxe[0] = -1;
Triangles[0].idxe[1] = 1;
Triangles[0].idxe[2] = -1;
nT++;
OrientTriangleInList(Triangles[0]);
/* loop over all points */
for (i=0;i<NumPart;i++)
{
AddPoint(&P[i]);
//CheckPoints(i-1);
}
/* add ghost points */
ExtendWithGhostPoints();
/* check */
CheckTriangles();
CheckPoints(NumPart+NumgPart);
if (ThisTask==0)
printf("ConstructDelaunay done.\n");
}
/*! Compute the Voronoi tesselation from the Delonay one
*/
int ComputeVoronoi()
{
if (ThisTask==0)
printf("Start ComputeVoronoi...\n");
ComputeMediansProperties();
ComputeMediansIntersections();
ComputeDensity();
if (ThisTask==0)
printf("ComputeVoronoi done.\n");
}
/************************************************************/
/* */
/* PYTHON INTERFACE */
/* */
/************************************************************/
#ifdef PY_INTERFACE
#include <Python.h>
#include <numpy/arrayobject.h>
PyObject *gadget_GetAllDelaunayTriangles(self, args)
PyObject *self;
PyObject *args;
{
import_array(); /* needed but strange no ? */
PyObject *OutputList;
PyObject *OutputDict;
PyArrayObject *tri = NULL;
npy_intp dim[2];
int iT;
/* loop over all triangles */
OutputList = PyList_New(0);
for (iT=0;iT<nT;iT++)
{
/* 3x3 vector */
dim[0]=3;
dim[1]=3;
tri = (PyArrayObject *) PyArray_SimpleNew(2,dim,PyArray_DOUBLE);
*(double *) (tri->data + 0*(tri->strides[0]) + 0*tri->strides[1]) = Triangles[iT].P[0]->Pos[0];
*(double *) (tri->data + 0*(tri->strides[0]) + 1*tri->strides[1]) = Triangles[iT].P[0]->Pos[1];
*(double *) (tri->data + 0*(tri->strides[0]) + 2*tri->strides[1]) = 0;
*(double *) (tri->data + 1*(tri->strides[0]) + 0*tri->strides[1]) = Triangles[iT].P[1]->Pos[0];
*(double *) (tri->data + 1*(tri->strides[0]) + 1*tri->strides[1]) = Triangles[iT].P[1]->Pos[1];
*(double *) (tri->data + 1*(tri->strides[0]) + 2*tri->strides[1]) = 0;
*(double *) (tri->data + 2*(tri->strides[0]) + 0*tri->strides[1]) = Triangles[iT].P[2]->Pos[0];
*(double *) (tri->data + 2*(tri->strides[0]) + 1*tri->strides[1]) = Triangles[iT].P[2]->Pos[1];
*(double *) (tri->data + 2*(tri->strides[0]) + 2*tri->strides[1]) = 0;
OutputDict = PyDict_New();
PyDict_SetItem(OutputDict,PyString_FromString("id"),PyInt_FromLong(Triangles[iT].idx) );
PyDict_SetItem(OutputDict,PyString_FromString("coord"),(PyObject*)tri);
//(PyObject*)tri
PyList_Append(OutputList, OutputDict );
}
return Py_BuildValue("O",OutputList);
}
PyObject *gadget_GetAllvPoints(self, args)
PyObject *self;
PyObject *args;
{
import_array(); /* needed but strange no ? */
PyArrayObject *pos;
npy_intp ld[2];
int i;
ld[0] = nvPoints;
ld[1] = 3;
pos = (PyArrayObject *) PyArray_SimpleNew(2,ld,PyArray_FLOAT);
for (i = 0; i < pos->dimensions[0]; i++)
{
*(float *) (pos->data + i*(pos->strides[0]) + 0*pos->strides[1]) = vPoints[i].Pos[0];
*(float *) (pos->data + i*(pos->strides[0]) + 1*pos->strides[1]) = vPoints[i].Pos[1];
*(float *) (pos->data + i*(pos->strides[0]) + 2*pos->strides[1]) = 0;
}
return PyArray_Return(pos);
}
PyObject *gadget_GetAllvDensities(PyObject* self)
{
import_array(); /* needed but strange no ? */
PyArrayObject *rho;
npy_intp ld[1];
int i;
ld[0] = NumPart;
rho = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE);
for (i = 0; i < rho->dimensions[0]; i++)
{
*(double *) (rho->data + i*(rho->strides[0])) = P[i].Density;
}
return PyArray_Return(rho);
}
PyObject *gadget_GetAllvVolumes(PyObject* self)
{
import_array(); /* needed but strange no ? */
PyArrayObject *volume;
npy_intp ld[1];
int i;
ld[0] = NumPart;
volume = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE);
for (i = 0; i < volume->dimensions[0]; i++)
{
*(double *) (volume->data + i*(volume->strides[0])) = P[i].Volume;
}
return PyArray_Return(volume);
}
PyObject *gadget_GetvPointsForOnePoint(self, args)
PyObject *self;
PyObject *args;
{
import_array(); /* needed but strange no ? */
PyArrayObject *pos;
npy_intp ld[2];
int i;
int np=0;
if (!PyArg_ParseTuple(args,"i",&i))
return NULL;
int next;
next = P[i].ivPoint;
np = P[i].nvPoints;
ld[0] = np;
ld[1] = 3;
pos = (PyArrayObject *) PyArray_SimpleNew(2,ld,PyArray_FLOAT);
for (i = 0; i < pos->dimensions[0]; i++)
{
*(float *) (pos->data + i*(pos->strides[0]) + 0*pos->strides[1]) = vPoints[next].Pos[0];
*(float *) (pos->data + i*(pos->strides[0]) + 1*pos->strides[1]) = vPoints[next].Pos[1];
*(float *) (pos->data + i*(pos->strides[0]) + 2*pos->strides[1]) = 0;
next = vPoints[next].next;
}
next = vPoints[next].next;
if (next!=0)
{
printf("error in tessel_get_vPointsForOnePoint\n");
return NULL;
}
return PyArray_Return(pos);
}
PyObject *gadget_GetAllGhostPositions(PyObject* self)
{
PyArrayObject *pos;
npy_intp ld[2];
int i;
import_array(); /* needed but strange no ? */
ld[0] = NumgPart;
ld[1] = 3;
pos = (PyArrayObject *) PyArray_SimpleNew(2,ld,PyArray_FLOAT);
for (i = 0; i < pos->dimensions[0]; i++)
{
*(float *) (pos->data + i*(pos->strides[0]) + 0*pos->strides[1]) = P[NumPart+i].Pos[0];
*(float *) (pos->data + i*(pos->strides[0]) + 1*pos->strides[1]) = P[NumPart+i].Pos[1];
*(float *) (pos->data + i*(pos->strides[0]) + 2*pos->strides[1]) = P[NumPart+i].Pos[2];
}
return PyArray_Return(pos);
}
PyObject *gadget_GetAllGhostvDensities(PyObject* self)
{
import_array(); /* needed but strange no ? */
PyArrayObject *rho;
npy_intp ld[1];
int i;
ld[0] = NumgPart;
rho = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE);
for (i = 0; i < rho->dimensions[0]; i++)
{
*(double *) (rho->data + i*(rho->strides[0])) = P[NumPart+i].Density;
}
return PyArray_Return(rho);
}
PyObject *gadget_GetAllGhostvVolumes(PyObject* self)
{
import_array(); /* needed but strange no ? */
PyArrayObject *volume;
npy_intp ld[1];
int i;
ld[0] = NumgPart;
volume = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE);
for (i = 0; i < volume->dimensions[0]; i++)
{
*(double *) (volume->data + i*(volume->strides[0])) = P[NumPart+i].Volume;
}
return PyArray_Return(volume);
}
#endif /* PY_INTERFACE */
#endif

Event Timeline