diff --git a/PyGear/examples.tessel/Readme b/PyGear/examples.tessel/Readme
new file mode 100644
index 0000000..7ed6e6d
--- /dev/null
+++ b/PyGear/examples.tessel/Readme
@@ -0,0 +1,12 @@
+
+./test_mpi.py
+
+
+./plot_triangles.py triangles.dat
+./plot_voronoi.py voronoi.dat
+
+
+
+
+#
+test_evol.py
diff --git a/PyGear/examples.tessel/TODO b/PyGear/examples.tessel/TODO
index 6c44987..f8956c8 100644
--- a/PyGear/examples.tessel/TODO
+++ b/PyGear/examples.tessel/TODO
@@ -1,34 +1,98 @@
+
+
+- dump triangles
+
+[ {"id":idx, "coord", [[x,y,z],[x,y,z],[x,y,z]] } ]
+
+
+3x
+n*(3*3)
+
+
+
+
+
+
+########################################
+#
+
+
+########################
+# draw the Voronoi
+
+
+for i in xrange(nbtot.nbody):
+ vpos = gadget.GetvPointsForOnePoint(i)
+ #pt.scatter(nb.pos[i,0],nb.pos[i,1],lw=0,s=50)
+
+ plot.draw_cell(vpos,color=[rho[i],rho[i],rho[i]],alpha=0.8)
+
+
+
+########################
+# draw the triangles
+
+TriangleList = gadget.GetAllDelaunayTriangles()
+
+i = 0
+for Triangle in TriangleList:
+ P1 = Triangle['coord'][0]
+ P2 = Triangle['coord'][1]
+ P3 = Triangle['coord'][2]
+ plot.draw_triangle(P1,P2,P3,c='r')
+
+ #cm = 1/3.*(P1+P2+P3)
+ #pt.text(cm[0],cm[1],Triangle['id'],fontsize=12,horizontalalignment='center',verticalalignment='center')
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
# EnergyInt EnergyPot EnergyKin
- pas de baisse de pression...
- viscosit� artificielle (Hess)
python test_evol.py
- snapshot des cellules
- graphs des cellules -> film
- conservation d'�nergie
- rapatrier EnergyInt (Entropy)
-> initializer correctement Entropy
-
- pas de temps -> vsig
si on construit 2x delaunay, �a plante....
-- parallelisme
diff --git a/PyGear/examples.tessel/plot_triangles.py b/PyGear/examples.tessel/plot_triangles.py
new file mode 100755
index 0000000..2927b94
--- /dev/null
+++ b/PyGear/examples.tessel/plot_triangles.py
@@ -0,0 +1,65 @@
+#!/usr/bin/env python
+
+
+from numpy import *
+from pNbody import io
+import sys
+
+
+file = sys.argv[1]
+
+
+
+def read_triangles(filename):
+
+
+ f = open(file,'r')
+
+ header = io.ReadBlock(f,(int32))
+ n = header[0]
+
+ triangles = io.ReadDataBlock(f,float,shape=(n,3,3))
+ num = io.ReadDataBlock(f,int32,shape=(n,))
+
+ return triangles,num
+
+
+
+triangles,num = read_triangles(file)
+
+
+######################################################################
+#
+# PLOT
+#
+######################################################################
+#sys.exit()
+
+import Ptools as pt
+import plot
+
+
+########################
+# draw a box
+plot.draw_box(x=array([0,1,1,0]),y=array([0,0,1,1]))
+
+########################
+# draw the triangles
+
+
+for Triangle in triangles:
+ P1 = Triangle[0]
+ P2 = Triangle[1]
+ P3 = Triangle[2]
+ plot.draw_triangle(P1,P2,P3,c='r')
+
+
+
+
+
+
+
+
+pt.axis([-0.1,1.1,-0.1,1.1])
+pt.show()
+
diff --git a/PyGear/examples.tessel/plot_voronoi.py b/PyGear/examples.tessel/plot_voronoi.py
new file mode 100755
index 0000000..6ecf3f7
--- /dev/null
+++ b/PyGear/examples.tessel/plot_voronoi.py
@@ -0,0 +1,66 @@
+#!/usr/bin/env python
+
+
+from numpy import *
+from pNbody import io
+import sys
+
+
+file = sys.argv[1]
+
+
+
+def read_voronoi(filename):
+
+
+ f = open(file,'r')
+
+ header = io.ReadBlock(f,(int32))
+ n = header[0]
+
+ cells = n*[""]
+
+ for i in xrange(n):
+
+ header = io.ReadBlock(f,(int32))
+ np = header[0]
+
+ points = io.ReadDataBlock(f,float,shape=(np,3))
+
+ cells[i]=points
+
+
+ return cells
+
+
+
+cells = read_voronoi(file)
+
+
+
+######################################################################
+#
+# PLOT
+#
+######################################################################
+#sys.exit()
+
+import Ptools as pt
+import plot
+
+
+########################
+# draw a box
+plot.draw_box(x=array([0,1,1,0]),y=array([0,0,1,1]))
+
+
+
+for vpos in cells:
+ plot.draw_cell(vpos,alpha=0.8)
+
+
+
+
+pt.axis([-0.1,1.1,-0.1,1.1])
+pt.show()
+
diff --git a/PyGear/examples.tessel/test_evol.py b/PyGear/examples.tessel/test_evol.py
index 2c8443d..6018417 100755
--- a/PyGear/examples.tessel/test_evol.py
+++ b/PyGear/examples.tessel/test_evol.py
@@ -1,335 +1,339 @@
#!/usr/bin/env python
from pNbody import *
from pNbody import ic
from pNbody import libutil
from numpy import *
import sys
import time
import copy
from PyGear import gadget
def SaveFile(gadget,i):
nb.pos = gadget.GetAllPositions()
nb.vel = gadget.GetAllVelocities()
nb.mass = gadget.GetAllMasses()
nb.num = gadget.GetAllID()
nb.tpe = None
#nb.u = gadget.GetAllEnergySpec() # not defined
nb.rho = gadget.GetAllvDensities() # note the v
nb.rsp = gadget.GetAllvVolumes()
nb.u = gadget.GetAllvEnergySpec()
nb.init()
# save
nb.rename('snap/snap%05d.dat'%i)
nb.set_pio("yes")
nb.write()
def SaveEnergy(gadget,time):
nb.pos = gadget.GetAllPositions()
nb.vel = gadget.GetAllVelocities()
nb.mass = gadget.GetAllMasses()
nb.num = gadget.GetAllID()
nb.tpe = None
#nb.u = gadget.GetAllEnergySpec() # not defined
nb.rho = gadget.GetAllvDensities() # note the v
nb.rsp = gadget.GetAllvVolumes()
nb.u = gadget.GetAllvEnergySpec()
nb.init()
Ekin = nb.Ekin()
Eint = sum(nb.u)
Epot = 0
efile = 'snap/energy.txt'
if os.path.isfile(efile):
f = open(efile,'a')
else:
f = open(efile,'w')
f.write("# Time EnergyInt EnergyPot EnergyKin \n")
f.write( "%g %g %g %g\n"%(time,Eint,Epot,Ekin) )
f.close()
nb = Nbody("snap.dat",ftype="gadget")
###################################
# init PyGear
gadget.InitMPI() # init MPI
gadget.InitDefaultParameters() # init default parameters
gadget.Info()
params = {}
params['ErrTolTheta'] = 0.7
params['DesNumNgb'] = 15
params['MaxNumNgbDeviation'] = 3
params['UnitLength_in_cm'] = 3.085e+21
params['UnitMass_in_g'] = 4.435693e+44
params['UnitVelocity_in_cm_per_s'] = 97824708.2699
params['BoxSize'] = 1.0
params['PartAllocFactor'] = 10
gadget.SetParameters(params)
###################################
# load particles
gadget.LoadParticles(array(nb.npart),nb.pos,nb.vel,nb.mass,nb.num,nb.tpe,nb.u)
###################################
# compute Delaunay
t1 = time.time()
gadget.ConstructDelaunay()
t2 = time.time()
print "Delaunay","in",(t2-t1),"s"
###################################
# compute Voronoi
t1 = time.time()
gadget.ComputeVoronoi()
t2 = time.time()
print "Voronoi","in",(t2-t1),"s"
###################################
# now, we have vDensity, we can compute
# the entropy
gadget.InitvEntropy()
###################################
# retrieve ghost points
gpos = gadget.GetAllGhostPositions()
gnb = Nbody(pos=gpos,ftype='gadget')
gnb.u = zeros(gnb.nbody)
gnb.rho = gadget.GetAllGhostvDensities()
gnb.rsp = gadget.GetAllGhostvVolumes()
gnb.init()
# save
gnb.rename('ghost.dat')
gnb.set_pio("yes")
gnb.write()
###################################
# retrieve points
nb.pos = gadget.GetAllPositions()
nb.vel = gadget.GetAllVelocities()
nb.mass = gadget.GetAllMasses()
nb.num = gadget.GetAllID()
nb.tpe = None
#nb.u = gadget.GetAllEnergySpec() # not defined
nb.rho = gadget.GetAllvDensities() # note the v
nb.rsp = gadget.GetAllvVolumes() # note the v
nb.u = gadget.GetAllvEnergySpec() # not the v
nb.init()
# save
nb.rename('snap2.dat')
nb.set_pio("yes")
nb.write()
###################################
# first accel
gadget.ComputeTesselAccelerations()
acc = gadget.GetAllvAccelerations()
acc = sqrt(acc[:,0]**2+acc[:,1]**2+acc[:,2]**2)
dt = gadget.tessel_get_timestep()
# first half drift step
gadget.tessel_drift(dt/2)
###################################################
#
# now loop
#
###################################################
for step in xrange(10000):
print ""
print "###########################"
print "Step %d dt=%g"%(step,dt )
print ""
gadget.domain_Decomposition()
gadget.gravity_tree()
###################################
# compute Delaunay
t1 = time.time()
gadget.ConstructDelaunay()
t2 = time.time()
print "Delaunay","in",(t2-t1),"s"
###################################
# compute Voronoi
t1 = time.time()
gadget.ComputeVoronoi()
t2 = time.time()
print "Voronoi","in",(t2-t1),"s"
###################################
# compute acceleration
gadget.ComputeTesselAccelerations()
acc = gadget.GetAllvAccelerations()
acc = sqrt(acc[:,0]**2+acc[:,1]**2+acc[:,2]**2)
# leap frog
gadget.tessel_kick(dt/2) # the syst. is inchronized here
# new timestep
dt = gadget.tessel_get_timestep()
gadget.tessel_kick(dt/2)
gadget.tessel_drift(dt)
# save output
if fmod(step,10)==0:
SaveFile(gadget,step)
# save energy
SaveEnergy(gadget,step)
-
+
+
+ gadget.tessel_dump_triangles('snap/triangles%04d.dat'%step)
+ gadget.tessel_dump_voronoi('snap/voronoi%04d.dat'%step)
+
######################################################################
# PLOT
#
######################################################################
# save all
#nbtot = nb + gnb
nbtot = copy.deepcopy(nb) # this is important to avoid to mix particles
nbtot.append(gnb,do_not_sort=True) # then they are incompatible with vpos below
nbtot.rename('snap+ghost.dat')
nbtot.write()
###################################
# get some values
TriangleList = gadget.GetAllDelaunayTriangles()
rho,mn,mx,cd = libutil.set_ranges(nbtot.rho,scale='log')
rho = rho/255.
#rho = clip(nbtot.rho,0,1)
import Ptools as pt
import plot
########################
# draw a box
plot.draw_box(x=array([0,1,1,0]),y=array([0,0,1,1]))
########################
# draw the Voronoi
cmap = pt.GetColormap('rainbow4',revesed=False)
# single point with vpoints
for i in xrange(nbtot.nbody):
vpos = gadget.GetvPointsForOnePoint(i)
#pt.scatter(nb.pos[i,0],nb.pos[i,1],lw=0,s=50)
plot.draw_cell(vpos,color=[rho[i],rho[i],rho[i]],alpha=0.8)
pt.axis([-0.1,1.1,-0.1,1.1])
pt.show()
diff --git a/PyGear/examples.tessel/test_mpi.py b/PyGear/examples.tessel/test_mpi.py
index a36ef68..6b0afe7 100755
--- a/PyGear/examples.tessel/test_mpi.py
+++ b/PyGear/examples.tessel/test_mpi.py
@@ -1,191 +1,194 @@
#!/usr/bin/env python
from pNbody import *
from pNbody import ic
from pNbody import libutil
from numpy import *
import sys
import time
import copy
from PyGear import gadget
nb = Nbody("snap.dat",ftype="gadget")
###################################
# init PyGear
gadget.InitMPI() # init MPI
gadget.InitDefaultParameters() # init default parameters
gadget.Info()
params = {}
params['ErrTolTheta'] = 0.7
params['DesNumNgb'] = 15
params['MaxNumNgbDeviation'] = 3
params['UnitLength_in_cm'] = 3.085e+21
params['UnitMass_in_g'] = 4.435693e+44
params['UnitVelocity_in_cm_per_s'] = 97824708.2699
params['BoxSize'] = 1.0
params['PartAllocFactor'] = 10
gadget.SetParameters(params)
###################################
# load particles
gadget.LoadParticles(array(nb.npart),nb.pos,nb.vel,nb.mass,nb.num,nb.tpe)
###################################
# compute Delaunay
t1 = time.time()
gadget.ConstructDelaunay()
t2 = time.time()
print "Delaunay","in",(t2-t1),"s"
###################################
# compute Voronoi
t1 = time.time()
gadget.ComputeVoronoi()
t2 = time.time()
print "Voronoi","in",(t2-t1),"s"
#sys.exit()
###################################
# retrieve ghost points
gpos = gadget.GetAllGhostPositions()
gnb = Nbody(pos=gpos,ftype='gadget')
gnb.u = zeros(gnb.nbody)
gnb.rho = gadget.GetAllGhostvDensities()
gnb.rsp = gadget.GetAllGhostvVolumes()
gnb.init()
# save
gnb.rename('ghost.dat')
gnb.set_pio("yes")
gnb.write()
###################################
# retrieve points
nb.pos = gadget.GetAllPositions()
nb.vel = gadget.GetAllVelocities()
nb.mass = gadget.GetAllMasses()
nb.num = gadget.GetAllID()
nb.tpe = None
nb.u = gadget.GetAllEnergySpec() # not defined
nb.rho = gadget.GetAllvDensities() # note the v
nb.rsp = gadget.GetAllvVolumes()
nb.rsp = gadget.GetAllDensities() # put sph densities in rsp
nb.init()
# save
nb.rename('snap2.dat')
nb.set_pio("yes")
nb.write()
# save all
#nbtot = nb + gnb
nbtot = copy.deepcopy(nb) # this is important to avoid to mix particles
nbtot.append(gnb,do_not_sort=True) # then they are incompatible with vpos below
nbtot.rename('snap+ghost.dat')
nbtot.write()
###################################
# get some values
TriangleList = gadget.GetAllDelaunayTriangles()
rho,mn,mx,cd = libutil.set_ranges(nbtot.rho,scale='log')
rho = rho/255.
#rho = clip(nbtot.rho,0,1)
+gadget.tessel_dump_triangles('triangles.dat')
+gadget.tessel_dump_voronoi('voronoi.dat')
+
######################################################################
#
# PLOT
#
######################################################################
#sys.exit()
import Ptools as pt
import plot
########################
# draw a box
plot.draw_box(x=array([0,1,1,0]),y=array([0,0,1,1]))
########################
# draw the triangles
'''
i = 0
for Triangle in TriangleList:
P1 = Triangle['coord'][0]
P2 = Triangle['coord'][1]
P3 = Triangle['coord'][2]
plot.draw_triangle(P1,P2,P3,c='r')
#cm = 1/3.*(P1+P2+P3)
#pt.text(cm[0],cm[1],Triangle['id'],fontsize=12,horizontalalignment='center',verticalalignment='center')
'''
########################
# draw the points
#pt.scatter( nb.pos[:,0], nb.pos[:,1],lw=0,s=50,color='r')
#pt.scatter(gnb.pos[:,0],gnb.pos[:,1],lw=0,s=50)
########################
# draw the Voronoi
cmap = pt.GetColormap('rainbow4',revesed=False)
# single point with vpoints
for i in xrange(nbtot.nbody):
vpos = gadget.GetvPointsForOnePoint(i)
#pt.scatter(nb.pos[i,0],nb.pos[i,1],lw=0,s=50)
plot.draw_cell(vpos,color=[rho[i],rho[i],rho[i]],alpha=0.8)
pt.axis([-0.1,1.1,-0.1,1.1])
pt.show()
diff --git a/src/proto.h b/src/proto.h
index cc1f65f..055b5a0 100644
--- a/src/proto.h
+++ b/src/proto.h
@@ -1,527 +1,532 @@
/*! \file proto.h
* \brief this file contains all function prototypes of the code
*/
#ifndef ALLVARS_H
#include "allvars.h"
#endif
#ifdef HAVE_HDF5
#include <hdf5.h>
#endif
void advance_and_find_timesteps(void);
void allocate_commbuffers(void);
void allocate_memory(void);
void begrun(void);
int blockpresent(enum iofields blocknr);
#ifdef BLOCK_SKIPPING
int blockabsent(enum iofields blocknr);
#endif
void catch_abort(int sig);
void catch_fatal(int sig);
void check_omega(void);
void close_outputfiles(void);
int compare_key(const void *a, const void *b);
void compute_accelerations(int mode);
void compute_global_quantities_of_system(void);
void compute_potential(void);
int dens_compare_key(const void *a, const void *b);
void density(int mode);
void density_decouple(void);
void density_evaluate(int i, int mode);
#ifdef CHIMIE
int stars_dens_compare_key(const void *a, const void *b);
void stars_density(void);
void stars_density_evaluate(int i, int mode);
#endif
void distribute_file(int nfiles, int firstfile, int firsttask, int lasttask, int *filenr, int *master, int *last);
double dmax(double, double);
double dmin(double, double);
void do_box_wrapping(void);
void domain_Decomposition(void);
int domain_compare_key(const void *a, const void *b);
int domain_compare_key(const void *a, const void *b);
int domain_compare_toplist(const void *a, const void *b);
void domain_countToGo(void);
void domain_decompose(void);
void domain_determineTopTree(void);
void domain_exchangeParticles(int partner, int sphflag, int send_count, int recv_count);
void domain_findExchangeNumbers(int task, int partner, int sphflag, int *send, int *recv);
void domain_findExtent(void);
int domain_findSplit(int cpustart, int ncpu, int first, int last);
int domain_findSplityr(int cpustart, int ncpu, int first, int last);
void domain_shiftSplit(void);
void domain_shiftSplityr(void);
void domain_sumCost(void);
void domain_topsplit(int node, peanokey startkey);
void domain_topsplit_local(int node, peanokey startkey);
double drift_integ(double a, void *param);
void dump_particles(void);
void empty_read_buffer(enum iofields blocknr, int offset, int pc, int type);
void endrun(int);
void energy_statistics(void);
#ifdef ADVANCEDSTATISTICS
void advanced_energy_statistics(void);
#endif
void every_timestep_stuff(void);
void ewald_corr(double dx, double dy, double dz, double *fper);
void ewald_force(int ii, int jj, int kk, double x[3], double force[3]);
void ewald_init(void);
double ewald_pot_corr(double dx, double dy, double dz);
double ewald_psi(double x[3]);
void fill_Tab_IO_Labels(void);
void fill_write_buffer(enum iofields blocknr, int *pindex, int pc, int type);
void find_dt_displacement_constraint(double hfac);
int find_files(char *fname);
int find_next_outputtime(int time);
void find_next_sync_point_and_drift(void);
void force_create_empty_nodes(int no, int topnode, int bits, int x, int y, int z, int *nodecount, int *nextfree);
void force_exchange_pseudodata(void);
void force_flag_localnodes(void);
void force_insert_pseudo_particles(void);
void force_setupnonrecursive(int no);
void force_treeallocate(int maxnodes, int maxpart);
int force_treebuild(int npart);
int force_treebuild_single(int npart);
int force_treeevaluate(int target, int mode, double *ewaldcountsum);
int force_treeevaluate_direct(int target, int mode);
int force_treeevaluate_ewald_correction(int target, int mode, double pos_x, double pos_y, double pos_z, double aold);
void force_treeevaluate_potential(int target, int type);
void force_treeevaluate_potential_shortrange(int target, int mode);
int force_treeevaluate_shortrange(int target, int mode);
void force_treefree(void);
void force_treeupdate_pseudos(void);
void force_update_hmax(void);
void force_update_len(void);
void force_update_node(int no, int flag);
void force_update_node_hmax_local(void);
void force_update_node_hmax_toptree(void);
void force_update_node_len_local(void);
void force_update_node_len_toptree(void);
void force_update_node_recursive(int no, int sib, int father);
void force_update_pseudoparticles(void);
void force_update_size_of_parent_node(int no);
void free_memory(void);
int get_bytes_per_blockelement(enum iofields blocknr);
void get_dataset_name(enum iofields blocknr, char *buf);
int get_datatype_in_block(enum iofields blocknr);
double get_drift_factor(int time0, int time1);
double get_gravkick_factor(int time0, int time1);
double get_hydrokick_factor(int time0, int time1);
int get_particles_in_block(enum iofields blocknr, int *typelist);
double get_random_number(int id);
#ifdef SFR
double get_StarFormation_random_number(int id);
#endif
#ifdef FEEDBACK_WIND
double get_FeedbackWind_random_number(int id);
#endif
#ifdef CHIMIE
double get_Chimie_random_number(int id);
#endif
#ifdef CHIMIE_KINETIC_FEEDBACK
double get_ChimieKineticFeedback_random_number(int id);
#endif
int get_timestep(int p, double *a, int flag);
int get_values_per_blockelement(enum iofields blocknr);
int grav_tree_compare_key(const void *a, const void *b);
void gravity_forcetest(void);
void gravity_tree(void);
void gravity_tree_shortrange(void);
double gravkick_integ(double a, void *param);
int hydro_compare_key(const void *a, const void *b);
void hydro_evaluate(int target, int mode);
void hydro_force(void);
double hydrokick_integ(double a, void *param);
int imax(int, int);
int imin(int, int);
void init(void);
void init_drift_table(void);
void init_peano_map(void);
#ifdef COSMICTIME
void init_cosmictime_table(void);
double get_cosmictime_difference(int time0, int time1);
#endif
void long_range_force(void);
void long_range_init(void);
void long_range_init_regionsize(void);
void move_particles(int time0, int time1);
size_t my_fread(void *ptr, size_t size, size_t nmemb, FILE * stream);
size_t my_fwrite(void *ptr, size_t size, size_t nmemb, FILE * stream);
int ngb_clear_buf(FLOAT searchcenter[3], FLOAT hguess, int numngb);
void ngb_treeallocate(int npart);
void ngb_treebuild(void);
int ngb_treefind_pairs(FLOAT searchcenter[3], FLOAT hsml, int phase, int *startnode);
#ifdef MULTIPHASE
int ngb_treefind_phase_pairs(FLOAT searchcenter[3], FLOAT hsml, int phase, int *startnode);
int ngb_treefind_sticky_collisions(FLOAT searchcenter[3], FLOAT hguess, int phase, int *startnode);
#endif
int ngb_treefind_variable(FLOAT searchcenter[3], FLOAT hguess, int phase, int *startnode);
#ifdef CHIMIE
int ngb_treefind_variable_for_chimie(FLOAT searchcenter[3], FLOAT hguess, int *startnode);
#endif
void ngb_treefree(void);
void ngb_treesearch(int);
void ngb_treesearch_pairs(int);
void ngb_update_nodes(void);
void open_outputfiles(void);
peanokey peano_hilbert_key(int x, int y, int z, int bits);
void peano_hilbert_order(void);
void pm_init_nonperiodic(void);
void pm_init_nonperiodic_allocate(int dimprod);
void pm_init_nonperiodic_free(void);
void pm_init_periodic(void);
void pm_init_periodic_allocate(int dimprod);
void pm_init_periodic_free(void);
void pm_init_regionsize(void);
void pm_setup_nonperiodic_kernel(void);
int pmforce_nonperiodic(int grnr);
void pmforce_periodic(void);
int pmpotential_nonperiodic(int grnr);
void pmpotential_periodic(void);
double pow(double, double); /* on some old DEC Alphas, the correct prototype for pow() is missing, even when math.h is included */
void read_file(char *fname, int readTask, int lastTask);
void read_header_attributes_in_hdf5(char *fname);
void read_ic(char *fname);
int read_outputlist(char *fname);
void read_parameter_file(char *fname);
void readjust_timebase(double TimeMax_old, double TimeMax_new);
void reorder_gas(void);
void reorder_particles(void);
#ifdef STELLAR_PROP
void reorder_stars(void);
void reorder_st(void);
#endif
void restart(int mod);
void run(void);
void savepositions(int num);
double second(void);
void seed_glass(void);
void set_random_numbers(void);
void set_softenings(void);
void set_units(void);
void init_local_sys_state(void);
void setup_smoothinglengths(void);
#ifdef CHIMIE
void stars_setup_smoothinglengths(void);
#endif
void statistics(void);
void terminate_processes(void);
double timediff(double t0, double t1);
#ifdef HAVE_HDF5
void write_header_attributes_in_hdf5(hid_t handle);
#endif
void write_file(char *fname, int readTask, int lastTask);
void write_pid_file(void);
#ifdef COOLING
int init_cooling(FLOAT metallicity);
int init_cooling_with_metals();
double cooling_function(double temperature);
double cooling_function_with_metals(double temperature,double metal);
void init_from_new_redshift(double Redshift);
double J_0();
double J_nu(double e);
double sigma_rad_HI(double e);
double sigma_rad_HeI(double e);
double sigma_rad_HeII(double e);
double cooling_bremstrahlung_HI(double T);
double cooling_bremstrahlung_HeI(double T);
double cooling_bremstrahlung_HeII(double T);
double cooling_ionization_HI(double T);
double cooling_ionization_HeI(double T);
double cooling_ionization_HeII(double T);
double cooling_recombination_HI(double T);
double cooling_recombination_HeI(double T);
double cooling_recombination_HeII(double T);
double cooling_dielectric_recombination(double T);
double cooling_excitation_HI(double T);
double cooling_excitation_HII(double T);
double cooling_compton(double T);
double A_HII(double T);
double A_HeIId(double T);
double A_HeII(double T);
double A_HeIII(double T);
double G_HI(double T);
double G_HeI(double T);
double G_HeII(double T);
double G_gHI();
double G_gHeI();
double G_gHeII();
double G_gHI_t(double J0);
double G_gHeI_t(double J0);
double G_gHeII_t(double J0);
double G_gHI_w();
double G_gHeI_w();
double G_gHeII_w();
double heating_radiative_HI();
double heating_radiative_HeI();
double heating_radiative_HeII();
double heating_radiative_HI_t(double J0);
double heating_radiative_HeI_t(double J0);
double heating_radiative_HeII_t(double J0);
double heating_radiative_HI_w();
double heating_radiative_HeI_w();
double heating_radiative_HeII_w();
double heating_compton();
void print_cooling(double T,double c1,double c2,double c3,double c4,double c5,double c6,double c7,double c8,double c9,double c10,double c11,double c12,double c13,double h1, double h2, double h3, double h4);
void compute_densities(double T,double X,double* n_H, double* n_HI,double* n_HII,double* n_HEI,double* n_HEII,double* n_HEIII,double* n_E,double* mu);
void compute_cooling_from_T_and_Nh(double T,double X,double n_H,double *c1,double *c2,double *c3,double *c4,double *c5,double *c6,double *c7,double *c8,double *c9,double *c10,double *c11,double *c12,double *c13,double *h1, double *h2, double *h3, double *h4);
double compute_cooling_from_Egyspec_and_Density(double Egyspec,double Density, double *MeanWeight);
double DoCooling(FLOAT Density,FLOAT Entropy,int Phase,int i,FLOAT DtEntropyVisc, double dt, double hubble_a);
void CoolingForOne(int i,int t0,int t1,double a3inv,double hubble_a);
void cooling();
double lambda(FLOAT density,FLOAT egyspec, int phase, int i);
#endif
#ifdef HEATING
void heating();
double gamma_fct(FLOAT Density,FLOAT Entropy,int i);
#endif
#ifdef AGN_HEATING
void agn_heating();
double gamma_fct(FLOAT density,double r, double SpecPower);
double HeatingRadialDependency(double r);
#endif
#ifdef MULTIPHASE
void update_phase(void);
void init_sticky(void);
void sticky(void);
void sticky_compute_energy_kin(int mode);
void sticky_collisions(void);
void sticky_collisions2(int loop);
void sticky_evaluate(int target, int mode, int loop);
int sticky_compare_key(const void *a, const void *b);
#endif
#ifdef FEEDBACK_WIND
void feedbackwind_compute_energy_kin(int mode);
#endif
#ifdef CHIMIE
void init_chimie(void);
void check_chimie(void);
void chimie(void);
void do_chimie(void);
void chimie_evaluate(int target, int mode);
int chimie_compare_key(const void *a, const void *b);
int get_nelts();
char* get_Element(i);
float get_SolarAbundance(i);
#if defined(CHIMIE_THERMAL_FEEDBACK) && defined(CHIMIE_COMPUTE_THERMAL_FEEDBACK_ENERGY)
void chimie_compute_energy_int(int mode);
#endif
#if defined(CHIMIE_KINETIC_FEEDBACK) && defined(CHIMIE_COMPUTE_KINETIC_FEEDBACK_ENERGY)
void chimie_compute_energy_kin(int mode);
#endif
#ifdef CHIMIE_KINETIC_FEEDBACK
void chimie_apply_wind(void);
#endif
#endif
#ifdef OUTERPOTENTIAL
void init_outer_potential(void);
void outer_forces(void);
void outer_potential(void);
#ifdef NFW
void init_outer_potential_nfw(void);
void outer_forces_nfw(void);
void outer_potential_nfw(void);
#endif
#ifdef PLUMMER
void init_outer_potential_plummer(void);
void outer_forces_plummer(void);
void outer_potential_plummer(void);
#endif
#ifdef PISOTHERM
void init_outer_potential_pisotherm(void);
void outer_forces_pisotherm(void);
void outer_potential_pisotherm(void);
double potential_f(double r, void * params);
double get_potential(double r);
#endif
#ifdef CORIOLIS
void init_outer_potential_coriolis(void);
void set_outer_potential_coriolis(void);
void outer_forces_coriolis(void);
void outer_potential_coriolis(void);
#endif
#endif
#ifdef SFR
void star_formation(void);
void rearrange_particle_sequence(void);
void sfr_compute_energy_int(int mode);
void sfr_check_number_of_stars(int mode);
#endif
#ifdef AGN_ACCRETION
void compute_agn_accretion(void);
#endif
#ifdef BUBBLES
void init_bubble(void);
void make_bubble(void);
void create_bubble(int sign);
#endif
#ifdef BONDI_ACCRETION
void bondi_accretion(void);
#endif
#ifdef PNBODY
void init_pnbody();
void finalize_pnbody();
void compute_pnbody();
#endif
#ifdef AB_TURB
void init_turb();
#endif
#if defined(ART_VISCO_MM)|| defined(ART_VISCO_RO) || defined(ART_VISCO_CD)
void move_art_visc(int i,double dt_drift);
#ifdef ART_VISCO_CD
void art_visc_allocate();
void art_visc_free();
void compute_art_visc(int i);
#endif
#endif
#ifdef TESSEL
void ConstructDelaunay();
void ComputeVoronoi();
void setup_searching_radius();
int ngb_treefind_variable_for_tessel(FLOAT searchcenter[3], FLOAT hsml, int phase, int *startnode);
void ghost();
void tessel_compute_accelerations();
void tessel_convert_energy_to_entropy();
void tessel_kick(float dt_kick);
void tessel_drift(float dt_drift);
double tessel_get_timestep();
int CheckCompletenessForThisPoint(int i);
int ghost_compare_key(const void *a, const void *b);
void CheckTriangles();
void AddGhostPoints(int istart,int nadd);
+
+void dump_triangles(char *filename);
+void dump_voronoi(char *filename);
+
+
#ifdef PY_INTERFACE
#include <Python.h>
PyObject *gadget_GetAllDelaunayTriangles(self, args);
PyObject *gadget_GetAllvPoints(self, args);
PyObject *gadget_GetAllvDensities(PyObject* self);
PyObject *gadget_GetAllvVolumes(PyObject* self);
PyObject *gadget_GetAllvPressures(PyObject* self);
PyObject *gadget_GetAllvEnergySpec(PyObject* self);
PyObject *gadget_GetAllvAccelerations(PyObject* self);
PyObject *gadget_GetvPointsForOnePoint(self, args);
PyObject *gadget_GetNgbPointsForOnePoint(self, args);
PyObject *gadget_GetNgbPointsAndFacesForOnePoint(self, args);
PyObject *gadget_GetAllGhostPositions(PyObject* self);
PyObject *gadget_GetAllGhostvDensities(PyObject* self);
PyObject *gadget_GetAllGhostvVolumes(PyObject* self);
#endif
#endif
#ifdef PY_INTERFACE
void allocate_commbuffersQ(void);
void density_sub(void);
void density_evaluate_sub(int i, int mode);
void do_box_wrappingQ(void);
void domain_DecompositionQ(void);
void domain_decomposeQ(void);
int domain_findSplitQ(int cpustart, int ncpu, int first, int last);
void domain_shiftSplitQ(void);
void domain_findExchangeNumbersQ(int task, int partner, int sphflag, int *send, int *recv);
void domain_exchangeParticlesQ(int partner, int sphflag, int send_count, int recv_count);
void domain_countToGoQ(void);
void domain_walktoptreeQ(int no);
void domain_sumCostQ(void);
void domain_findExtentQ(void);
void domain_determineTopTreeQ(void);
void domain_topsplit_localQ(int node, peanokey startkey);
void domain_topsplitQ(int node, peanokey startkey);
int force_treeevaluate_sub(int target, int mode, double *ewaldcountsum);
void force_treeevaluate_potential_sub(int target, int type);
void force_treeevaluate_potential_shortrange_sub(int target, int mode);
int force_treeevaluate_shortrange_sub(int target, int mode);
void gravity_tree_sub(void);
void sph(void);
void sph_evaluate(int target, int mode);
void sph_sub(void);
void sph_evaluate_sub(int target, int mode);
void sph_thermal_conductivity(void);
void sph_evaluate_thermal_conductivity(int target, int mode);
int sph_compare_key(const void *a, const void *b);
void peano_hilbert_orderQ(void);
void reorder_gasQ(void);
void reorder_particlesQ(void);
void setup_smoothinglengths_sub(void);
#endif
diff --git a/src/python_interface.c b/src/python_interface.c
index d25f0a9..dd9e4f0 100644
--- a/src/python_interface.c
+++ b/src/python_interface.c
@@ -1,3336 +1,3368 @@
#ifdef PY_INTERFACE
#include <Python.h>
#include <math.h>
#include <string.h>
#include <stdio.h>
#include <numpy/arrayobject.h>
#include <mpi.h>
#include "allvars.h"
#include "proto.h"
#define TO_INT(a) ( (PyArrayObject*) PyArray_CastToType(a, PyArray_DescrFromType(NPY_INT) ,0) )
#define TO_DOUBLE(a) ( (PyArrayObject*) PyArray_CastToType(a, PyArray_DescrFromType(NPY_DOUBLE) ,0) )
#define TO_FLOAT(a) ( (PyArrayObject*) PyArray_CastToType(a, PyArray_DescrFromType(NPY_FLOAT) ,0) )
static int Init()
{
/* main.c */
RestartFlag = 0;
All.CPU_TreeConstruction = All.CPU_TreeWalk = All.CPU_Gravity = All.CPU_Potential = All.CPU_Domain =
All.CPU_Snapshot = All.CPU_Total = All.CPU_CommSum = All.CPU_Imbalance = All.CPU_Hydro =
All.CPU_HydCompWalk = All.CPU_HydCommSumm = All.CPU_HydImbalance =
All.CPU_EnsureNgb = All.CPU_Predict = All.CPU_TimeLine = All.CPU_PM = All.CPU_Peano = 0;
CPUThisRun = 0;
/* from init.c, after read ic */
int i, j;
double a3;
All.Time = All.TimeBegin;
All.Ti_Current = 0;
if(All.ComovingIntegrationOn)
{
All.Timebase_interval = (log(All.TimeMax) - log(All.TimeBegin)) / TIMEBASE;
a3 = All.Time * All.Time * All.Time;
}
else
{
All.Timebase_interval = (All.TimeMax - All.TimeBegin) / TIMEBASE;
a3 = 1;
}
set_softenings();
All.NumCurrentTiStep = 0; /* setup some counters */
All.SnapshotFileCount = 0;
if(RestartFlag == 2)
All.SnapshotFileCount = atoi(All.InitCondFile + strlen(All.InitCondFile) - 3) + 1;
All.TotNumOfForces = 0;
All.NumForcesSinceLastDomainDecomp = 0;
if(All.ComovingIntegrationOn)
if(All.PeriodicBoundariesOn == 1)
check_omega();
All.TimeLastStatistics = All.TimeBegin - All.TimeBetStatistics;
if(All.ComovingIntegrationOn) /* change to new velocity variable */
{
for(i = 0; i < NumPart; i++)
for(j = 0; j < 3; j++)
P[i].Vel[j] *= sqrt(All.Time) * All.Time;
}
for(i = 0; i < NumPart; i++) /* start-up initialization */
{
for(j = 0; j < 3; j++)
P[i].GravAccel[j] = 0;
#ifdef PMGRID
for(j = 0; j < 3; j++)
P[i].GravPM[j] = 0;
#endif
P[i].Ti_endstep = 0;
P[i].Ti_begstep = 0;
P[i].OldAcc = 0;
P[i].GravCost = 1;
P[i].Potential = 0;
}
#ifdef PMGRID
All.PM_Ti_endstep = All.PM_Ti_begstep = 0;
#endif
#ifdef FLEXSTEPS
All.PresentMinStep = TIMEBASE;
for(i = 0; i < NumPart; i++) /* start-up initialization */
{
P[i].FlexStepGrp = (int) (TIMEBASE * get_random_number(P[i].ID));
}
#endif
for(i = 0; i < N_gas; i++) /* initialize sph_properties */
{
for(j = 0; j < 3; j++)
{
SphP[i].VelPred[j] = P[i].Vel[j];
SphP[i].HydroAccel[j] = 0;
}
SphP[i].DtEntropy = 0;
if(RestartFlag == 0)
{
SphP[i].Hsml = 0;
SphP[i].Density = -1;
}
}
ngb_treeallocate(MAX_NGB);
force_treeallocate(All.TreeAllocFactor * All.MaxPart, All.MaxPart);
All.NumForcesSinceLastDomainDecomp = 1 + All.TotNumPart * All.TreeDomainUpdateFrequency;
Flag_FullStep = 1; /* to ensure that Peano-Hilber order is done */
domain_Decomposition(); /* do initial domain decomposition (gives equal numbers of particles) */
ngb_treebuild(); /* will build tree */
setup_smoothinglengths();
#ifdef CHIMIE
#ifndef CHIMIE_INPUT_ALL
stars_setup_smoothinglengths();
#endif
#endif
#ifdef TESSEL
setup_searching_radius();
#endif
TreeReconstructFlag = 1;
#ifndef TESSEL /*do not convert if tessel is used. In tessel, the conversion is done with tessel_convert_energy_to_entropy() */
/* at this point, the entropy variable normally contains the
* internal energy, read in from the initial conditions file, unless the file
* explicitly signals that the initial conditions contain the entropy directly.
* Once the density has been computed, we can convert thermal energy to entropy.
*/
#ifndef ISOTHERM_EQS
if(header.flag_entropy_instead_u == 0)
for(i = 0; i < N_gas; i++)
SphP[i].Entropy = GAMMA_MINUS1 * SphP[i].Entropy / pow(SphP[i].Density / a3, GAMMA_MINUS1);
#endif
#endif /* TESSEL */
return 1;
}
static void Begrun1()
{
struct global_data_all_processes all;
if(ThisTask == 0)
{
printf("\nThis is pyGadget, version `%s'.\n", GADGETVERSION);
printf("\nRunning on %d processors.\n", NTask);
}
//read_parameter_file(ParameterFile); /* ... read in parameters for this run */
allocate_commbuffers(); /* ... allocate buffer-memory for particle
exchange during force computation */
set_units();
#if defined(PERIODIC) && (!defined(PMGRID) || defined(FORCETEST))
ewald_init();
#endif
//open_outputfiles();
random_generator = gsl_rng_alloc(gsl_rng_ranlxd1);
gsl_rng_set(random_generator, 42); /* start-up seed */
#ifdef PMGRID
long_range_init();
#endif
All.TimeLastRestartFile = CPUThisRun;
if(RestartFlag == 0 || RestartFlag == 2)
{
set_random_numbers();
}
else
{
all = All; /* save global variables. (will be read from restart file) */
restart(RestartFlag); /* ... read restart file. Note: This also resets
all variables in the struct `All'.
However, during the run, some variables in the parameter
file are allowed to be changed, if desired. These need to
copied in the way below.
Note: All.PartAllocFactor is treated in restart() separately.
*/
All.MinSizeTimestep = all.MinSizeTimestep;
All.MaxSizeTimestep = all.MaxSizeTimestep;
All.BufferSize = all.BufferSize;
All.BunchSizeForce = all.BunchSizeForce;
All.BunchSizeDensity = all.BunchSizeDensity;
All.BunchSizeHydro = all.BunchSizeHydro;
All.BunchSizeDomain = all.BunchSizeDomain;
All.TimeLimitCPU = all.TimeLimitCPU;
All.ResubmitOn = all.ResubmitOn;
All.TimeBetSnapshot = all.TimeBetSnapshot;
All.TimeBetStatistics = all.TimeBetStatistics;
All.CpuTimeBetRestartFile = all.CpuTimeBetRestartFile;
All.ErrTolIntAccuracy = all.ErrTolIntAccuracy;
All.MaxRMSDisplacementFac = all.MaxRMSDisplacementFac;
All.ErrTolForceAcc = all.ErrTolForceAcc;
All.TypeOfTimestepCriterion = all.TypeOfTimestepCriterion;
All.TypeOfOpeningCriterion = all.TypeOfOpeningCriterion;
All.NumFilesWrittenInParallel = all.NumFilesWrittenInParallel;
All.TreeDomainUpdateFrequency = all.TreeDomainUpdateFrequency;
All.SnapFormat = all.SnapFormat;
All.NumFilesPerSnapshot = all.NumFilesPerSnapshot;
All.MaxNumNgbDeviation = all.MaxNumNgbDeviation;
All.ArtBulkViscConst = all.ArtBulkViscConst;
All.OutputListOn = all.OutputListOn;
All.CourantFac = all.CourantFac;
All.OutputListLength = all.OutputListLength;
memcpy(All.OutputListTimes, all.OutputListTimes, sizeof(double) * All.OutputListLength);
strcpy(All.ResubmitCommand, all.ResubmitCommand);
strcpy(All.OutputListFilename, all.OutputListFilename);
strcpy(All.OutputDir, all.OutputDir);
strcpy(All.RestartFile, all.RestartFile);
strcpy(All.EnergyFile, all.EnergyFile);
strcpy(All.InfoFile, all.InfoFile);
strcpy(All.CpuFile, all.CpuFile);
strcpy(All.TimingsFile, all.TimingsFile);
strcpy(All.SnapshotFileBase, all.SnapshotFileBase);
if(All.TimeMax != all.TimeMax)
readjust_timebase(All.TimeMax, all.TimeMax);
}
}
static void Begrun2()
{
if(RestartFlag == 0 || RestartFlag == 2)
Init(); /* ... read in initial model */
#ifdef PMGRID
long_range_init_regionsize();
#endif
if(All.ComovingIntegrationOn)
init_drift_table();
//if(RestartFlag == 2)
// All.Ti_nextoutput = find_next_outputtime(All.Ti_Current + 1);
//else
// All.Ti_nextoutput = find_next_outputtime(All.Ti_Current);
All.TimeLastRestartFile = CPUThisRun;
}
/************************************************************/
/* PYTHON INTERFACE */
/************************************************************/
static PyObject *gadget_Info(PyObject *self, PyObject *args, PyObject *kwds)
{
printf("I am proc %d among %d procs.\n",ThisTask,NTask);
return Py_BuildValue("i",1);
}
static PyObject *gadget_InitMPI(PyObject *self, PyObject *args, PyObject *kwds)
{
//MPI_Init(0, 0); /* this is done in mpi4py */
MPI_Comm_rank(MPI_COMM_WORLD, &ThisTask);
MPI_Comm_size(MPI_COMM_WORLD, &NTask);
for(PTask = 0; NTask > (1 << PTask); PTask++);
return Py_BuildValue("i",1);
}
static PyObject * gadget_InitDefaultParameters(PyObject* self)
{
/* list of Gadget parameters */
/*
All.InitCondFile ="ICs/cluster_littleendian.dat";
All.OutputDir ="cluster/";
All.EnergyFile ="energy.txt";
All.InfoFile ="info.txt";
All.TimingsFile ="timings.txt";
All.CpuFile ="cpu.txt";
All.RestartFile ="restart";
All.SnapshotFileBase ="snapshot";
All.OutputListFilename ="parameterfiles/outputs_lcdm_gas.txt";
*/
/* CPU time -limit */
All.TimeLimitCPU = 36000; /* = 10 hours */
All.ResubmitOn = 0;
//All.ResubmitCommand = "my-scriptfile";
All.ICFormat = 1;
All.SnapFormat = 1;
All.ComovingIntegrationOn = 0;
All.TypeOfTimestepCriterion = 0;
All.OutputListOn = 0;
All.PeriodicBoundariesOn = 0;
/* Caracteristics of run */
All.TimeBegin = 0.0; /*% Begin of the simulation (z=23)*/
All.TimeMax = 1.0;
All.Omega0 = 0;
All.OmegaLambda = 0;
All.OmegaBaryon = 0;
All.HubbleParam = 0;
All.BoxSize = 0;
/* Output frequency */
All.TimeBetSnapshot = 0.1;
All.TimeOfFirstSnapshot = 0.0; /*% 5 constant steps in log(a) */
All.CpuTimeBetRestartFile = 36000.0; /* here in seconds */
All.TimeBetStatistics = 0.05;
All.NumFilesPerSnapshot = 1;
All.NumFilesWrittenInParallel = 1;
/* Accuracy of time integration */
All.ErrTolIntAccuracy = 0.025;
All.MaxRMSDisplacementFac = 0.2;
All.CourantFac = 0.15;
All.MaxSizeTimestep = 0.03;
All.MinSizeTimestep = 0.0;
/* Tree algorithm, force accuracy, domain update frequency */
All.ErrTolTheta = 0.7;
All.TypeOfOpeningCriterion = 0;
All.ErrTolForceAcc = 0.005;
All.TreeDomainUpdateFrequency = 0.1;
/* Further parameters of SPH */
All.DesNumNgb = 50;
All.MaxNumNgbDeviation = 2;
All.ArtBulkViscConst = 0.8;
All.InitGasTemp = 0;
All.MinGasTemp = 0;
/* Memory allocation */
All.PartAllocFactor = 2.0;
All.TreeAllocFactor = 2.0;
All.BufferSize = 30;
/* System of units */
All.UnitLength_in_cm = 3.085678e21; /* 1.0 kpc */
All.UnitMass_in_g = 1.989e43; /* 1.0e10 solar masses */
All.UnitVelocity_in_cm_per_s = 1e5; /* 1 km/sec */
All.GravityConstantInternal = 0;
/* Softening lengths */
All.MinGasHsmlFractional = 0.25;
All.SofteningGas = 0.5;
All.SofteningHalo = 0.5;
All.SofteningDisk = 0.5;
All.SofteningBulge = 0.5;
All.SofteningStars = 0.5;
All.SofteningBndry = 0.5;
All.SofteningGasMaxPhys = 0.5;
All.SofteningHaloMaxPhys = 0.5;
All.SofteningDiskMaxPhys = 0.5;
All.SofteningBulgeMaxPhys = 0.5;
All.SofteningStarsMaxPhys = 0.5;
All.SofteningBndryMaxPhys = 0.5;
return Py_BuildValue("i",1);
}
static PyObject * gadget_GetParameters()
{
PyObject *dict;
PyObject *key;
PyObject *value;
dict = PyDict_New();
/* CPU time -limit */
key = PyString_FromString("TimeLimitCPU");
value = PyFloat_FromDouble(All.TimeLimitCPU);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("ResubmitOn");
value = PyFloat_FromDouble(All.ResubmitOn);
PyDict_SetItem(dict,key,value);
//All.ResubmitCommand
key = PyString_FromString("ICFormat");
value = PyInt_FromLong(All.ICFormat);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("SnapFormat");
value = PyInt_FromLong(All.SnapFormat);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("ComovingIntegrationOn");
value = PyInt_FromLong(All.ComovingIntegrationOn);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("TypeOfTimestepCriterion");
value = PyInt_FromLong(All.TypeOfTimestepCriterion);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("OutputListOn");
value = PyInt_FromLong(All.OutputListOn);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("PeriodicBoundariesOn");
value = PyInt_FromLong(All.PeriodicBoundariesOn);
PyDict_SetItem(dict,key,value);
/* Caracteristics of run */
key = PyString_FromString("TimeBegin");
value = PyFloat_FromDouble(All.TimeBegin);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("TimeMax");
value = PyFloat_FromDouble(All.TimeMax);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("Omega0");
value = PyFloat_FromDouble(All.Omega0);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("OmegaLambda");
value = PyFloat_FromDouble(All.OmegaLambda);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("OmegaBaryon");
value = PyFloat_FromDouble(All.OmegaBaryon);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("HubbleParam");
value = PyFloat_FromDouble(All.HubbleParam);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("BoxSize");
value = PyFloat_FromDouble(All.BoxSize);
PyDict_SetItem(dict,key,value);
/* Output frequency */
key = PyString_FromString("TimeBetSnapshot");
value = PyFloat_FromDouble(All.TimeBetSnapshot);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("TimeOfFirstSnapshot");
value = PyFloat_FromDouble(All.TimeOfFirstSnapshot);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("CpuTimeBetRestartFile");
value = PyFloat_FromDouble(All.CpuTimeBetRestartFile);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("TimeBetStatistics");
value = PyFloat_FromDouble(All.TimeBetStatistics);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("NumFilesPerSnapshot");
value = PyInt_FromLong(All.NumFilesPerSnapshot);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("NumFilesWrittenInParallel");
value = PyInt_FromLong(All.NumFilesWrittenInParallel);
PyDict_SetItem(dict,key,value);
/* Accuracy of time integration */
key = PyString_FromString("ErrTolIntAccuracy");
value = PyFloat_FromDouble(All.ErrTolIntAccuracy);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("MaxRMSDisplacementFac");
value = PyFloat_FromDouble(All.MaxRMSDisplacementFac);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("CourantFac");
value = PyFloat_FromDouble(All.CourantFac);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("MaxSizeTimestep");
value = PyFloat_FromDouble(All.MaxSizeTimestep);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("MinSizeTimestep");
value = PyFloat_FromDouble(All.MinSizeTimestep);
PyDict_SetItem(dict,key,value);
/* Tree algorithm, force accuracy, domain update frequency */
key = PyString_FromString("ErrTolTheta");
value = PyFloat_FromDouble(All.ErrTolTheta);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("TypeOfOpeningCriterion");
value = PyInt_FromLong(All.TypeOfOpeningCriterion);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("ErrTolForceAcc");
value = PyFloat_FromDouble(All.ErrTolForceAcc);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("TreeDomainUpdateFrequency");
value = PyFloat_FromDouble(All.TreeDomainUpdateFrequency);
PyDict_SetItem(dict,key,value);
/* Further parameters of SPH */
key = PyString_FromString("DesNumNgb");
value = PyInt_FromLong(All.DesNumNgb);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("MaxNumNgbDeviation");
value = PyInt_FromLong(All.MaxNumNgbDeviation);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("ArtBulkViscConst");
value = PyInt_FromLong(All.ArtBulkViscConst);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("InitGasTemp");
value = PyInt_FromLong(All.InitGasTemp);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("MinGasTemp");
value = PyInt_FromLong(All.MinGasTemp);
PyDict_SetItem(dict,key,value);
/* Memory allocation */
key = PyString_FromString("PartAllocFactor");
value = PyFloat_FromDouble(All.PartAllocFactor);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("TreeAllocFactor");
value = PyFloat_FromDouble(All.TreeAllocFactor);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("BufferSize");
value = PyInt_FromLong(All.BufferSize);
PyDict_SetItem(dict,key,value);
/* System of units */
key = PyString_FromString("UnitLength_in_cm");
value = PyFloat_FromDouble(All.UnitLength_in_cm);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("UnitMass_in_g");
value = PyFloat_FromDouble(All.UnitMass_in_g);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("UnitVelocity_in_cm_per_s");
value = PyFloat_FromDouble(All.UnitVelocity_in_cm_per_s);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("GravityConstantInternal");
value = PyFloat_FromDouble(All.GravityConstantInternal);
PyDict_SetItem(dict,key,value);
/* Softening lengths */
key = PyString_FromString("MinGasHsmlFractional");
value = PyFloat_FromDouble(All.MinGasHsmlFractional);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("SofteningGas");
value = PyFloat_FromDouble(All.SofteningGas);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("SofteningHalo");
value = PyFloat_FromDouble(All.SofteningHalo);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("SofteningDisk");
value = PyFloat_FromDouble(All.SofteningDisk);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("SofteningBulge");
value = PyFloat_FromDouble(All.SofteningBulge);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("SofteningStars");
value = PyFloat_FromDouble(All.SofteningStars);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("SofteningBndry");
value = PyFloat_FromDouble(All.SofteningBndry);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("SofteningGasMaxPhys");
value = PyFloat_FromDouble(All.SofteningGasMaxPhys);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("SofteningHaloMaxPhys");
value = PyFloat_FromDouble(All.SofteningHaloMaxPhys);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("SofteningDiskMaxPhys");
value = PyFloat_FromDouble(All.SofteningDiskMaxPhys);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("SofteningBulgeMaxPhys");
value = PyFloat_FromDouble(All.SofteningBulgeMaxPhys);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("SofteningStarsMaxPhys");
value = PyFloat_FromDouble(All.SofteningStarsMaxPhys);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("SofteningBndryMaxPhys");
value = PyFloat_FromDouble(All.SofteningBndryMaxPhys);
PyDict_SetItem(dict,key,value);
/*
key = PyString_FromString("OutputInfo");
value = PyFloat_FromDouble(All.OutputInfo);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("PeanoHilbertOrder");
value = PyFloat_FromDouble(All.PeanoHilbertOrder);
PyDict_SetItem(dict,key,value);
*/
return Py_BuildValue("O",dict);
}
static PyObject * gadget_SetParameters(PyObject *self, PyObject *args)
{
PyObject *dict;
PyObject *key;
PyObject *value;
int ivalue;
float fvalue;
double dvalue;
/* here, we can have either arguments or dict directly */
if(PyDict_Check(args))
{
dict = args;
}
else
{
if (! PyArg_ParseTuple(args, "O",&dict))
return NULL;
}
/* check that it is a PyDictObject */
if(!PyDict_Check(dict))
{
PyErr_SetString(PyExc_AttributeError, "argument is not a dictionary.");
return NULL;
}
Py_ssize_t pos=0;
while(PyDict_Next(dict,&pos,&key,&value))
{
if(PyString_Check(key))
{
/* CPU time -limit */
if(strcmp(PyString_AsString(key), "TimeLimitCPU")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.TimeLimitCPU = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "ResubmitOn")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.ResubmitOn = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "ICFormat")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.ICFormat = PyInt_AsLong(value);
}
if(strcmp(PyString_AsString(key), "SnapFormat")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.SnapFormat = PyInt_AsLong(value);
}
if(strcmp(PyString_AsString(key), "ComovingIntegrationOn")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.ComovingIntegrationOn = PyInt_AsLong(value);
}
if(strcmp(PyString_AsString(key), "TypeOfTimestepCriterion")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.TypeOfTimestepCriterion = PyInt_AsLong(value);
}
if(strcmp(PyString_AsString(key), "OutputListOn")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.OutputListOn = PyInt_AsLong(value);
}
if(strcmp(PyString_AsString(key), "PeriodicBoundariesOn")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.PeriodicBoundariesOn = PyInt_AsLong(value);
}
/* Caracteristics of run */
if(strcmp(PyString_AsString(key), "TimeBegin")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.TimeBegin = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "TimeMax")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.TimeMax = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "Omega0")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.Omega0 = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "OmegaLambda")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.OmegaLambda = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "OmegaBaryon")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.OmegaBaryon = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "HubbleParam")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.HubbleParam = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "BoxSize")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.BoxSize = PyFloat_AsDouble(value);
}
/* Output frequency */
if(strcmp(PyString_AsString(key), "TimeBetSnapshot")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.TimeBetSnapshot = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "TimeOfFirstSnapshot")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.TimeOfFirstSnapshot = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "CpuTimeBetRestartFile")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.CpuTimeBetRestartFile = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "TimeBetStatistics")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.TimeBetStatistics = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "NumFilesPerSnapshot")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.NumFilesPerSnapshot = PyInt_AsLong(value);
}
if(strcmp(PyString_AsString(key), "NumFilesWrittenInParallel")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.NumFilesWrittenInParallel = PyInt_AsLong(value);
}
/* Accuracy of time integration */
if(strcmp(PyString_AsString(key), "ErrTolIntAccuracy")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.ErrTolIntAccuracy = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "MaxRMSDisplacementFac")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.MaxRMSDisplacementFac = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "CourantFac")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.CourantFac = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "MaxSizeTimestep")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.MaxSizeTimestep = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "MinSizeTimestep")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.MinSizeTimestep = PyFloat_AsDouble(value);
}
/* Tree algorithm, force accuracy, domain update frequency */
if(strcmp(PyString_AsString(key), "ErrTolTheta")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.ErrTolTheta = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "TypeOfOpeningCriterion")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.TypeOfOpeningCriterion = PyInt_AsLong(value);
}
if(strcmp(PyString_AsString(key), "ErrTolForceAcc")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.ErrTolForceAcc = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "TreeDomainUpdateFrequency")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.TreeDomainUpdateFrequency = PyFloat_AsDouble(value);
}
/* Further parameters of SPH */
if(strcmp(PyString_AsString(key), "DesNumNgb")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.DesNumNgb = PyInt_AsLong(value);
}
if(strcmp(PyString_AsString(key), "MaxNumNgbDeviation")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.MaxNumNgbDeviation = PyInt_AsLong(value);
}
if(strcmp(PyString_AsString(key), "ArtBulkViscConst")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.ArtBulkViscConst = PyInt_AsLong(value);
}
if(strcmp(PyString_AsString(key), "InitGasTemp")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.InitGasTemp = PyInt_AsLong(value);
}
if(strcmp(PyString_AsString(key), "MinGasTemp")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.MinGasTemp = PyInt_AsLong(value);
}
/* Memory allocation */
if(strcmp(PyString_AsString(key), "PartAllocFactor")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.PartAllocFactor = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "TreeAllocFactor")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.TreeAllocFactor = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "BufferSize")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.BufferSize = PyInt_AsLong(value);
}
/* System of units */
if(strcmp(PyString_AsString(key), "UnitLength_in_cm")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.UnitLength_in_cm = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "UnitMass_in_g")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.UnitMass_in_g = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "UnitVelocity_in_cm_per_s")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.UnitVelocity_in_cm_per_s = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "GravityConstantInternal")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.GravityConstantInternal = PyFloat_AsDouble(value);
}
/* Softening lengths */
if(strcmp(PyString_AsString(key), "MinGasHsmlFractional")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.MinGasHsmlFractional = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "SofteningGas")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.SofteningGas = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "SofteningHalo")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.SofteningHalo = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "SofteningDisk")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.SofteningDisk = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "SofteningBulge")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.SofteningBulge = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "SofteningStars")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.SofteningStars = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "SofteningBndry")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.SofteningBndry = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "SofteningGasMaxPhys")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.SofteningGasMaxPhys = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "SofteningHaloMaxPhys")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.SofteningHaloMaxPhys = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "SofteningDiskMaxPhys")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.SofteningDiskMaxPhys = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "SofteningBulgeMaxPhys")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.SofteningBulgeMaxPhys = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "SofteningStarsMaxPhys")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.SofteningStarsMaxPhys = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "SofteningBndryMaxPhys")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.SofteningBndryMaxPhys = PyFloat_AsDouble(value);
}
}
}
return Py_BuildValue("i",1);
}
static PyObject *
gadget_check_parser(PyObject *self, PyObject *args, PyObject *keywds)
{
int voltage;
char *state = "a stiff";
char *action = "voom";
char *type = "Norwegian Blue";
static char *kwlist[] = {"voltage", "state", "action", "type", NULL};
if (!PyArg_ParseTupleAndKeywords(args, keywds, "i|sss", kwlist,
&voltage, &state, &action, &type))
return NULL;
printf("-- This parrot wouldn't %s if you put %i Volts through it.\n",
action, voltage);
printf("-- Lovely plumage, the %s -- It's %s!\n", type, state);
Py_INCREF(Py_None);
return Py_None;
}
static PyObject *gadget_Free(PyObject *self, PyObject *args, PyObject *kwds)
{
free_memory();
ngb_treefree();
force_treefree();
free(Exportflag );
free(DomainStartList);
free(DomainEndList);
free(TopNodes);
free(DomainWork);
free(DomainCount);
free(DomainCountSph);
free(DomainTask);
free(DomainNodeIndex);
free(DomainTreeNodeLen);
free(DomainHmax);
free(DomainMoment);
free(CommBuffer);
gsl_rng_free(random_generator);
return Py_BuildValue("i",1);
}
static PyObject *gadget_LoadParticles(PyObject *self, PyObject *args, PyObject *kwds)
{
int i,j;
size_t bytes;
PyArrayObject *ntype=Py_None;
PyArrayObject *pos=Py_None;
PyArrayObject *vel=Py_None;
PyArrayObject *mass=Py_None;
PyArrayObject *num=Py_None;
PyArrayObject *tpe=Py_None;
PyArrayObject *u=Py_None;
PyArrayObject *rho=Py_None;
static char *kwlist[] = {"npart", "pos","vel","mass","num","tpe","u","rho", NULL};
if (! PyArg_ParseTupleAndKeywords(args,kwds, "OOOOOO|OO",kwlist,&ntype,&pos,&vel,&mass,&num,&tpe,&u,&rho ))
return NULL;
/* check type */
if (!(PyArray_Check(pos)))
{
PyErr_SetString(PyExc_ValueError,"aruments 2 must be array.");
return NULL;
}
/* check type */
if (!(PyArray_Check(mass)))
{
PyErr_SetString(PyExc_ValueError,"aruments 3 must be array.");
return NULL;
}
/* check dimension */
if ( (pos->nd!=2))
{
PyErr_SetString(PyExc_ValueError,"Dimension of argument 2 must be 2.");
return NULL;
}
/* check dimension */
if ( (mass->nd!=1))
{
PyErr_SetString(PyExc_ValueError,"Dimension of argument 3 must be 1.");
return NULL;
}
/* check size */
if ( (pos->dimensions[1]!=3))
{
PyErr_SetString(PyExc_ValueError,"First size of argument must be 3.");
return NULL;
}
/* check size */
if ( (pos->dimensions[0]!=mass->dimensions[0]))
{
PyErr_SetString(PyExc_ValueError,"Size of argument 2 must be similar to argument 3.");
return NULL;
}
/* ensure double */
ntype = TO_INT(ntype);
pos = TO_FLOAT(pos);
vel = TO_FLOAT(vel);
mass = TO_FLOAT(mass);
num = TO_INT(num);
tpe = TO_FLOAT(tpe);
/* optional variables */
if (u!=Py_None)
u = TO_FLOAT(u);
if (rho!=Py_None)
rho = TO_FLOAT(rho);
/***************************************
* some inits *
/***************************************/
RestartFlag = 0;
Begrun1();
/***************************************
* LOAD PARTILES *
/***************************************/
NumPart = 0;
N_gas = *(int*) (ntype->data + 0*(ntype->strides[0]));
for (i = 0; i < 6; i++)
NumPart += *(int*) (ntype->data + i*(ntype->strides[0]));
if (NumPart!=pos->dimensions[0])
{
PyErr_SetString(PyExc_ValueError,"Numpart != pos->dimensions[0].");
return NULL;
}
MPI_Allreduce(&NumPart, &All.TotNumPart, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
MPI_Allreduce(&N_gas, &All.TotN_gas, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
All.MaxPart = All.PartAllocFactor * (All.TotNumPart / NTask);
All.MaxPartSph = All.PartAllocFactor * (All.TotN_gas / NTask);
/*********************/
/* allocate P */
/*********************/
if(!(P = malloc(bytes = All.MaxPart * sizeof(struct particle_data))))
{
printf("failed to allocate memory for `P' (%g MB).\n", bytes / (1024.0 * 1024.0));
endrun(1);
}
if(!(SphP = malloc(bytes = All.MaxPartSph * sizeof(struct sph_particle_data))))
{
printf("failed to allocate memory for `SphP' (%g MB) %d.\n", bytes / (1024.0 * 1024.0), sizeof(struct sph_particle_data));
endrun(1);
}
/*********************/
/* init P */
/*********************/
for (i = 0; i < NumPart; i++)
{
P[i].Pos[0] = *(float *) (pos->data + i*(pos->strides[0]) + 0*pos->strides[1]);
P[i].Pos[1] = *(float *) (pos->data + i*(pos->strides[0]) + 1*pos->strides[1]);
P[i].Pos[2] = *(float *) (pos->data + i*(pos->strides[0]) + 2*pos->strides[1]);
P[i].Vel[0] = *(float *) (vel->data + i*(vel->strides[0]) + 0*vel->strides[1]);
P[i].Vel[1] = *(float *) (vel->data + i*(vel->strides[0]) + 1*vel->strides[1]);
P[i].Vel[2] = *(float *) (vel->data + i*(vel->strides[0]) + 2*vel->strides[1]);
P[i].Mass = *(float *) (mass->data + i*(mass->strides[0]));
P[i].ID = *(unsigned int *) (num->data + i*(num->strides[0]));
P[i].Type = *(int *) (tpe->data + i*(tpe->strides[0]));
//P[i].Active = 1;
//printf(" P.ID=%d%09d\n", (int) (P[i].ID / 1000000000), (int) (P[i].ID % 1000000000));
#ifdef TESSEL
P[i].iPref = -1;
#endif
}
/*********************/
/* init SphP */
/*********************/
if (u!=Py_None)
for (i = 0; i < NumPart; i++)
{
SphP[i].Entropy = *(float *) (u->data + i*(u->strides[0]));
//#ifndef ISOTHERM_EQS
// SphP[i].Entropy = GAMMA_MINUS1 * SphP[i].Entropy / pow(SphP[i].Density / a3, GAMMA_MINUS1);
//#endif
}
if (rho!=Py_None)
for (i = 0; i < NumPart; i++)
{
SphP[i].Density = *(float *) (rho->data + i*(rho->strides[0]));
}
/***************************************
* END LOAD PARTICLES *
/***************************************/
/***************************************
* finish inits *
/***************************************/
Begrun2();
/***************************************
* free memory *
/***************************************/
/* free the memory allocated,
because these vectors where casted
and their memory is now handeled by the C part */
Py_DECREF(ntype);
Py_DECREF(pos);
Py_DECREF(vel);
Py_DECREF(mass);
Py_DECREF(num);
Py_DECREF(tpe);
/* optional variables */
if (u!=Py_None)
Py_DECREF(u);
if (rho!=Py_None)
Py_DECREF(rho);
return Py_BuildValue("i",1);
}
static PyObject *gadget_LoadParticlesQ(PyObject *self, PyObject *args, PyObject *kwds)
{
int i,j;
size_t bytes;
PyArrayObject *ntype,*pos,*vel,*mass,*num,*tpe;
static char *kwlist[] = {"npart", "pos","vel","mass","num","tpe", NULL};
if (! PyArg_ParseTupleAndKeywords(args, kwds, "|OOOOOO",kwlist,&ntype,&pos,&vel,&mass,&num,&tpe))
return Py_BuildValue("i",1);
/* check type */
if (!(PyArray_Check(pos)))
{
PyErr_SetString(PyExc_ValueError,"aruments 1 must be array.");
return NULL;
}
/* check type */
if (!(PyArray_Check(mass)))
{
PyErr_SetString(PyExc_ValueError,"aruments 2 must be array.");
return NULL;
}
/* check dimension */
if ( (pos->nd!=2))
{
PyErr_SetString(PyExc_ValueError,"Dimension of argument 1 must be 2.");
return NULL;
}
/* check dimension */
if ( (mass->nd!=1))
{
PyErr_SetString(PyExc_ValueError,"Dimension of argument 2 must be 1.");
return NULL;
}
/* check size */
if ( (pos->dimensions[1]!=3))
{
PyErr_SetString(PyExc_ValueError,"First size of argument must be 3.");
return NULL;
}
/* check size */
if ( (pos->dimensions[0]!=mass->dimensions[0]))
{
PyErr_SetString(PyExc_ValueError,"Size of argument 1 must be similar to argument 2.");
return NULL;
}
/* ensure double */
ntype = TO_INT(ntype);
pos = TO_FLOAT(pos);
vel = TO_FLOAT(vel);
mass = TO_FLOAT(mass);
num = TO_FLOAT(num);
tpe = TO_FLOAT(tpe);
/***************************************
* some inits *
/***************************************/
allocate_commbuffersQ();
/***************************************
* LOAD PARTILES *
/***************************************/
NumPartQ = 0;
N_gasQ = *(int*) (ntype->data + 0*(ntype->strides[0]));
for (i = 0; i < 6; i++)
NumPartQ += *(int*) (ntype->data + i*(ntype->strides[0]));
if (NumPartQ!=pos->dimensions[0])
{
PyErr_SetString(PyExc_ValueError,"NumpartQ != pos->dimensions[0].");
return NULL;
}
MPI_Allreduce(&NumPartQ, &All.TotNumPartQ, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
MPI_Allreduce(&N_gasQ, &All.TotN_gasQ, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
All.MaxPartQ = All.PartAllocFactor * (All.TotNumPartQ / NTask);
All.MaxPartSphQ = All.PartAllocFactor * (All.TotN_gasQ / NTask);
/*********************/
/* allocate Q */
/*********************/
if(!(Q = malloc(bytes = All.MaxPartQ * sizeof(struct particle_data))))
{
printf("failed to allocate memory for `Q' (%g MB).\n", bytes / (1024.0 * 1024.0));
endrun(1);
}
if(!(SphQ = malloc(bytes = All.MaxPartSphQ * sizeof(struct sph_particle_data))))
{
printf("failed to allocate memory for `SphQ' (%g MB) %d.\n", bytes / (1024.0 * 1024.0), sizeof(struct sph_particle_data));
endrun(1);
}
/*********************/
/* init P */
/*********************/
for (i = 0; i < NumPartQ; i++)
{
Q[i].Pos[0] = *(float *) (pos->data + i*(pos->strides[0]) + 0*pos->strides[1]);
Q[i].Pos[1] = *(float *) (pos->data + i*(pos->strides[0]) + 1*pos->strides[1]);
Q[i].Pos[2] = *(float *) (pos->data + i*(pos->strides[0]) + 2*pos->strides[1]);
Q[i].Vel[0] = *(float *) (vel->data + i*(vel->strides[0]) + 0*vel->strides[1]);
Q[i].Vel[1] = *(float *) (vel->data + i*(vel->strides[0]) + 1*vel->strides[1]);
Q[i].Vel[2] = *(float *) (vel->data + i*(vel->strides[0]) + 2*vel->strides[1]);
Q[i].Mass = *(float *) (mass->data + i*(mass->strides[0]));
Q[i].ID = *(unsigned int *) (num->data + i*(num->strides[0]));
Q[i].Type = *(int *) (tpe->data + i*(tpe->strides[0]));
//Q[i].Active = 1;
}
/***************************************
* END LOAD PARTILES *
/***************************************/
domain_DecompositionQ();
/***************************************
* finish inits *
/***************************************/
return Py_BuildValue("i",1);
}
static PyObject *gadget_AllPotential(PyObject *self)
{
compute_potential();
return Py_BuildValue("i",1);
}
static PyObject * gadget_GetAllPotential(PyObject* self)
{
PyArrayObject *pot;
npy_intp ld[1];
int i;
ld[0] = NumPart;
pot = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_FLOAT);
for (i = 0; i < pot->dimensions[0]; i++)
{
*(float *) (pot->data + i*(pot->strides[0])) = P[i].Potential;
}
return PyArray_Return(pot);
}
static PyObject *gadget_AllAcceleration(PyObject *self)
{
NumForceUpdate = NumPart;
gravity_tree();
return Py_BuildValue("i",1);
}
static PyObject * gadget_GetAllAcceleration(PyObject* self)
{
PyArrayObject *acc;
npy_intp ld[2];
int i;
ld[0] = NumPart;
ld[1] = 3;
acc = (PyArrayObject *) PyArray_SimpleNew(2,ld,PyArray_FLOAT);
for (i = 0; i < acc->dimensions[0]; i++)
{
*(float *) (acc->data + i*(acc->strides[0]) + 0*acc->strides[1]) = P[i].GravAccel[0];
*(float *) (acc->data + i*(acc->strides[0]) + 1*acc->strides[1]) = P[i].GravAccel[1];
*(float *) (acc->data + i*(acc->strides[0]) + 2*acc->strides[1]) = P[i].GravAccel[2];
}
return PyArray_Return(acc);
}
static PyObject *gadget_GetAllDensities(PyObject* self)
{
PyArrayObject *rho;
npy_intp ld[1];
int i;
ld[0] = N_gas;
rho = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_FLOAT);
for (i = 0; i < rho->dimensions[0]; i++)
{
*(float *) (rho->data + i*(rho->strides[0])) = SphP[i].Density;
}
return PyArray_Return(rho);
}
static PyObject *gadget_GetAllHsml(PyObject* self)
{
PyArrayObject *hsml;
npy_intp ld[1];
int i;
ld[0] = N_gas;
hsml = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_FLOAT);
for (i = 0; i < hsml->dimensions[0]; i++)
{
*(float *) (hsml->data + i*(hsml->strides[0])) = SphP[i].Hsml;
}
return PyArray_Return(hsml);
}
static PyObject *gadget_GetAllPositions(PyObject* self)
{
PyArrayObject *pos;
npy_intp ld[2];
int i;
ld[0] = NumPart;
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[i].Pos[0];
*(float *) (pos->data + i*(pos->strides[0]) + 1*pos->strides[1]) = P[i].Pos[1];
*(float *) (pos->data + i*(pos->strides[0]) + 2*pos->strides[1]) = P[i].Pos[2];
}
return PyArray_Return(pos);
}
static PyObject *gadget_GetAllVelocities(PyObject* self)
{
PyArrayObject *vel;
npy_intp ld[2];
int i;
ld[0] = NumPart;
ld[1] = 3;
vel = (PyArrayObject *) PyArray_SimpleNew(2,ld,PyArray_FLOAT);
for (i = 0; i < vel->dimensions[0]; i++)
{
*(float *) (vel->data + i*(vel->strides[0]) + 0*vel->strides[1]) = P[i].Vel[0];
*(float *) (vel->data + i*(vel->strides[0]) + 1*vel->strides[1]) = P[i].Vel[1];
*(float *) (vel->data + i*(vel->strides[0]) + 2*vel->strides[1]) = P[i].Vel[2];
}
return PyArray_Return(vel);
}
static PyObject *gadget_GetAllMasses(PyObject* self)
{
PyArrayObject *mass;
npy_intp ld[1];
int i;
ld[0] = NumPart;
mass = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_FLOAT);
for (i = 0; i < mass->dimensions[0]; i++)
{
*(float *) (mass->data + i*(mass->strides[0])) = P[i].Mass;
}
return PyArray_Return(mass);
}
static PyObject *gadget_GetAllEntropy(PyObject* self)
{
PyArrayObject *entropy;
npy_intp ld[1];
int i;
ld[0] = NumPart;
entropy = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_FLOAT);
for (i = 0; i < entropy->dimensions[0]; i++)
{
*(float *) (entropy->data + i*(entropy->strides[0])) = SphP[i].Entropy;
}
return PyArray_Return(entropy);
}
static PyObject *gadget_GetAllEnergySpec(PyObject* self)
{
PyArrayObject *energy;
npy_intp ld[1];
int i;
double a3inv;
ld[0] = NumPart;
energy = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_FLOAT);
if(All.ComovingIntegrationOn)
{
a3inv = 1 / (All.Time * All.Time * All.Time);
}
else
a3inv = 1;
for (i = 0; i < energy->dimensions[0]; i++)
{
#ifdef ISOTHERM_EQS
*(float *) (energy->data + i*(energy->strides[0])) = SphP[i].Entropy;
#else
a3inv = 1.;
*(float *) (energy->data + i*(energy->strides[0])) = dmax(All.MinEgySpec,SphP[i].Entropy / GAMMA_MINUS1 * pow(SphP[i].Density * a3inv, GAMMA_MINUS1));
#endif
}
return PyArray_Return(energy);
}
static PyObject *gadget_GetAllID(PyObject* self)
{
PyArrayObject *id;
npy_intp ld[1];
int i;
ld[0] = NumPart;
id = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_INT);
for (i = 0; i < id->dimensions[0]; i++)
{
*(float *) (id->data + i*(id->strides[0])) = P[i].ID;
}
return PyArray_Return(id);
}
static PyObject *gadget_GetAllTypes(PyObject* self)
{
PyArrayObject *type;
npy_intp ld[1];
int i;
ld[0] = NumPart;
type = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_INT);
for (i = 0; i < type->dimensions[0]; i++)
{
*(int *) (type->data + i*(type->strides[0])) = P[i].Type;
}
return PyArray_Return(type);
}
static PyObject *gadget_GetAllPositionsQ(PyObject* self)
{
PyArrayObject *pos;
npy_intp ld[2];
int i;
ld[0] = NumPartQ;
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]) = Q[i].Pos[0];
*(float *) (pos->data + i*(pos->strides[0]) + 1*pos->strides[1]) = Q[i].Pos[1];
*(float *) (pos->data + i*(pos->strides[0]) + 2*pos->strides[1]) = Q[i].Pos[2];
}
return PyArray_Return(pos);
}
static PyObject *gadget_GetAllVelocitiesQ(PyObject* self)
{
PyArrayObject *vel;
npy_intp ld[2];
int i;
ld[0] = NumPartQ;
ld[1] = 3;
vel = (PyArrayObject *) PyArray_SimpleNew(2,ld,PyArray_FLOAT);
for (i = 0; i < vel->dimensions[0]; i++)
{
*(float *) (vel->data + i*(vel->strides[0]) + 0*vel->strides[1]) = Q[i].Vel[0];
*(float *) (vel->data + i*(vel->strides[0]) + 1*vel->strides[1]) = Q[i].Vel[1];
*(float *) (vel->data + i*(vel->strides[0]) + 2*vel->strides[1]) = Q[i].Vel[2];
}
return PyArray_Return(vel);
}
static PyObject *gadget_GetAllMassesQ(PyObject* self)
{
PyArrayObject *mass;
npy_intp ld[1];
int i;
ld[0] = NumPartQ;
mass = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_FLOAT);
for (i = 0; i < mass->dimensions[0]; i++)
{
*(float *) (mass->data + i*(mass->strides[0])) = Q[i].Mass;
}
return PyArray_Return(mass);
}
static PyObject *gadget_GetAllIDQ(PyObject* self)
{
PyArrayObject *id;
npy_intp ld[1];
int i;
ld[0] = NumPartQ;
id = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_INT);
for (i = 0; i < id->dimensions[0]; i++)
{
*(float *) (id->data + i*(id->strides[0])) = Q[i].ID;
}
return PyArray_Return(id);
}
static PyObject *gadget_GetAllTypesQ(PyObject* self)
{
PyArrayObject *type;
npy_intp ld[1];
int i;
ld[0] = NumPartQ;
type = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_INT);
for (i = 0; i < type->dimensions[0]; i++)
{
*(int *) (type->data + i*(type->strides[0])) = Q[i].Type;
}
return PyArray_Return(type);
}
static PyObject *gadget_GetPos(PyObject *self, PyObject *args, PyObject *kwds)
{
int i,j;
size_t bytes;
PyArrayObject *pos;
if (! PyArg_ParseTuple(args, "O",&pos))
return PyString_FromString("error : GetPos");
for (i = 0; i < pos->dimensions[0]; i++)
{
*(float *) (pos->data + i*(pos->strides[0]) + 0*pos->strides[1]) = P[i].Pos[0];
*(float *) (pos->data + i*(pos->strides[0]) + 1*pos->strides[1]) = P[i].Pos[1];
*(float *) (pos->data + i*(pos->strides[0]) + 2*pos->strides[1]) = P[i].Pos[2];
}
//return PyArray_Return(Py_None);
return Py_BuildValue("i",1);
}
static PyObject * gadget_Potential(PyObject* self, PyObject *args)
{
PyArrayObject *pos;
float eps;
if (! PyArg_ParseTuple(args, "Of",&pos,&eps))
return PyString_FromString("error");
PyArrayObject *pot;
int i;
npy_intp ld[1];
int input_dimension;
size_t bytes;
input_dimension =pos->nd;
if (input_dimension != 2)
PyErr_SetString(PyExc_ValueError,"dimension of first argument must be 2");
pos = TO_FLOAT(pos);
/* create a NumPy object */
ld[0]=pos->dimensions[0];
pot = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_FLOAT);
NumPartQ = pos->dimensions[0];
All.ForceSofteningQ = eps;
if(!(Q = malloc(bytes = NumPartQ * sizeof(struct particle_data))))
{
printf("failed to allocate memory for `Q' (%g MB).\n", bytes / (1024.0 * 1024.0));
endrun(1);
}
if(!(SphQ = malloc(bytes = NumPartQ * sizeof(struct sph_particle_data))))
{
printf("failed to allocate memory for `SphQ' (%g MB).\n", bytes / (1024.0 * 1024.0));
endrun(1);
}
for (i = 0; i < pos->dimensions[0]; i++)
{
Q[i].Pos[0] = *(float *) (pos->data + i*(pos->strides[0]) + 0*pos->strides[1]);
Q[i].Pos[1] = *(float *) (pos->data + i*(pos->strides[0]) + 1*pos->strides[1]);
Q[i].Pos[2] = *(float *) (pos->data + i*(pos->strides[0]) + 2*pos->strides[1]);
Q[i].Type = 0;
Q[i].Mass = 0;
Q[i].Potential = 0;
}
compute_potential_sub();
for (i = 0; i < pos->dimensions[0]; i++)
{
*(float *)(pot->data + i*(pot->strides[0])) = Q[i].Potential;
}
free(Q);
free(SphQ);
return PyArray_Return(pot);
}
static PyObject * gadget_Acceleration(PyObject* self, PyObject *args)
{
PyArrayObject *pos;
float eps;
if (! PyArg_ParseTuple(args, "Of",&pos,&eps))
return PyString_FromString("error");
PyArrayObject *acc;
int i;
int input_dimension;
size_t bytes;
input_dimension =pos->nd;
if (input_dimension != 2)
PyErr_SetString(PyExc_ValueError,"dimension of first argument must be 2");
pos = TO_FLOAT(pos);
/* create a NumPy object */
acc = (PyArrayObject *) PyArray_SimpleNew(pos->nd,pos->dimensions,PyArray_FLOAT);
NumPartQ = pos->dimensions[0];
All.ForceSofteningQ = eps;
if(!(Q = malloc(bytes = NumPartQ * sizeof(struct particle_data))))
{
printf("failed to allocate memory for `Q' (%g MB).\n", bytes / (1024.0 * 1024.0));
endrun(1);
}
if(!(SphQ = malloc(bytes = NumPartQ * sizeof(struct sph_particle_data))))
{
printf("failed to allocate memory for `SphQ' (%g MB).\n", bytes / (1024.0 * 1024.0));
endrun(1);
}
for (i = 0; i < pos->dimensions[0]; i++)
{
Q[i].Pos[0] = *(float *) (pos->data + i*(pos->strides[0]) + 0*pos->strides[1]);
Q[i].Pos[1] = *(float *) (pos->data + i*(pos->strides[0]) + 1*pos->strides[1]);
Q[i].Pos[2] = *(float *) (pos->data + i*(pos->strides[0]) + 2*pos->strides[1]);
Q[i].Type = 0;
Q[i].Mass = 0;
Q[i].GravAccel[0] = 0;
Q[i].GravAccel[1] = 0;
Q[i].GravAccel[2] = 0;
}
gravity_tree_sub();
for (i = 0; i < pos->dimensions[0]; i++)
{
*(float *)(acc->data + i*(acc->strides[0]) + 0*acc->strides[1]) = Q[i].GravAccel[0];
*(float *)(acc->data + i*(acc->strides[0]) + 1*acc->strides[1]) = Q[i].GravAccel[1];
*(float *)(acc->data + i*(acc->strides[0]) + 2*acc->strides[1]) = Q[i].GravAccel[2];
}
free(Q);
free(SphQ);
return PyArray_Return(acc);
}
static PyObject * gadget_InitHsml(PyObject* self, PyObject *args)
{
PyArrayObject *pos,*hsml;
if (! PyArg_ParseTuple(args, "OO",&pos,&hsml))
return PyString_FromString("error");
int i;
int input_dimension;
size_t bytes;
int ld[1];
PyArrayObject *vden,*vhsml;
input_dimension =pos->nd;
if (input_dimension != 2)
PyErr_SetString(PyExc_ValueError,"dimension of first argument must be 2");
if (pos->dimensions[0] != hsml->dimensions[0])
PyErr_SetString(PyExc_ValueError,"pos and hsml must have the same dimension.");
pos = TO_FLOAT(pos);
hsml = TO_FLOAT(hsml);
/* create a NumPy object */
ld[0]=pos->dimensions[0];
vden = (PyArrayObject *) PyArray_SimpleNew(1,pos->dimensions,pos->descr->type_num);
vhsml = (PyArrayObject *) PyArray_SimpleNew(1,pos->dimensions,pos->descr->type_num);
NumPartQ = pos->dimensions[0];
N_gasQ = NumPartQ;
All.Ti_Current=1; /* need to flag active particles */
if(!(Q = malloc(bytes = NumPartQ * sizeof(struct particle_data))))
{
printf("failed to allocate memory for `Q' (%g MB).\n", bytes / (1024.0 * 1024.0));
endrun(1);
}
if(!(SphQ = malloc(bytes = NumPartQ * sizeof(struct sph_particle_data))))
{
printf("failed to allocate memory for `SphQ' (%g MB).\n", bytes / (1024.0 * 1024.0));
endrun(1);
}
for (i = 0; i < pos->dimensions[0]; i++)
{
Q[i].Pos[0] = *(float *) (pos->data + i*(pos->strides[0]) + 0*pos->strides[1]);
Q[i].Pos[1] = *(float *) (pos->data + i*(pos->strides[0]) + 1*pos->strides[1]);
Q[i].Pos[2] = *(float *) (pos->data + i*(pos->strides[0]) + 2*pos->strides[1]);
SphQ[i].Hsml = *(float *) (hsml->data + i*(hsml->strides[0]));
}
setup_smoothinglengths_sub();
for (i = 0; i < pos->dimensions[0]; i++)
{
*(float *)(vhsml->data + i*(vhsml->strides[0])) = SphQ[i].Hsml;
*(float *)(vden->data + i*(vden->strides[0])) = SphQ[i].Density;
}
free(Q);
free(SphQ);
return Py_BuildValue("OO",vden,vhsml);
}
static PyObject * gadget_Density(PyObject* self, PyObject *args)
{
PyArrayObject *pos,*hsml;
if (! PyArg_ParseTuple(args, "OO",&pos,&hsml))
return PyString_FromString("error");
int i;
int input_dimension;
size_t bytes;
int ld[1];
PyArrayObject *vden,*vhsml;
input_dimension =pos->nd;
if (input_dimension != 2)
PyErr_SetString(PyExc_ValueError,"dimension of first argument must be 2");
if (pos->dimensions[0] != hsml->dimensions[0])
PyErr_SetString(PyExc_ValueError,"pos and hsml must have the same dimension.");
pos = TO_FLOAT(pos);
hsml = TO_FLOAT(hsml);
/* create a NumPy object */
ld[0]=pos->dimensions[0];
vden = (PyArrayObject *) PyArray_SimpleNew(1,pos->dimensions,pos->descr->type_num);
vhsml = (PyArrayObject *) PyArray_SimpleNew(1,pos->dimensions,pos->descr->type_num);
NumPartQ = pos->dimensions[0];
N_gasQ = NumPartQ;
All.Ti_Current=1; /* need to flag active particles */
if(!(Q = malloc(bytes = NumPartQ * sizeof(struct particle_data))))
{
printf("failed to allocate memory for `Q' (%g MB).\n", bytes / (1024.0 * 1024.0));
endrun(1);
}
if(!(SphQ = malloc(bytes = NumPartQ * sizeof(struct sph_particle_data))))
{
printf("failed to allocate memory for `SphQ' (%g MB).\n", bytes / (1024.0 * 1024.0));
endrun(1);
}
for (i = 0; i < pos->dimensions[0]; i++)
{
Q[i].Pos[0] = *(float *) (pos->data + i*(pos->strides[0]) + 0*pos->strides[1]);
Q[i].Pos[1] = *(float *) (pos->data + i*(pos->strides[0]) + 1*pos->strides[1]);
Q[i].Pos[2] = *(float *) (pos->data + i*(pos->strides[0]) + 2*pos->strides[1]);
SphQ[i].Hsml = *(float *) (hsml->data + i*(hsml->strides[0]));
}
density_sub();
for (i = 0; i < pos->dimensions[0]; i++)
{
*(float *)(vhsml->data + i*(vhsml->strides[0])) = SphQ[i].Hsml;
*(float *)(vden->data + i*(vden->strides[0])) = SphQ[i].Density;
}
free(Q);
free(SphQ);
return Py_BuildValue("OO",vden,vhsml);
}
static PyObject * gadget_SphEvaluate(PyObject* self, PyObject *args)
{
PyArrayObject *pos,*hsml,*obs;
if (! PyArg_ParseTuple(args, "OOO",&pos,&hsml,&obs))
return PyString_FromString("error");
int i;
int input_dimension;
size_t bytes;
int ld[1];
PyArrayObject *vobs;
input_dimension =pos->nd;
if (input_dimension != 2)
PyErr_SetString(PyExc_ValueError,"dimension of first argument must be 2");
if (pos->dimensions[0] != hsml->dimensions[0])
PyErr_SetString(PyExc_ValueError,"pos and hsml must have the same dimension.");
if (obs->nd != 1)
PyErr_SetString(PyExc_ValueError,"dimension of obs must be 1.");
if (obs->dimensions[0] != NumPart)
PyErr_SetString(PyExc_ValueError,"The size of obs must be NumPart.");
pos = TO_FLOAT(pos);
hsml = TO_FLOAT(hsml);
obs = TO_FLOAT(obs);
/* create a NumPy object */
ld[0]=pos->dimensions[0];
vobs = (PyArrayObject *) PyArray_SimpleNew(1,pos->dimensions,pos->descr->type_num);
NumPartQ = pos->dimensions[0];
N_gasQ = NumPartQ;
All.Ti_Current=1; /* need to flag active particles */
if(!(Q = malloc(bytes = NumPartQ * sizeof(struct particle_data))))
{
printf("failed to allocate memory for `Q' (%g MB).\n", bytes / (1024.0 * 1024.0));
endrun(1);
}
if(!(SphQ = malloc(bytes = NumPartQ * sizeof(struct sph_particle_data))))
{
printf("failed to allocate memory for `SphQ' (%g MB).\n", bytes / (1024.0 * 1024.0));
endrun(1);
}
for (i = 0; i < pos->dimensions[0]; i++)
{
Q[i].Pos[0] = *(float *) (pos->data + i*(pos->strides[0]) + 0*pos->strides[1]);
Q[i].Pos[1] = *(float *) (pos->data + i*(pos->strides[0]) + 1*pos->strides[1]);
Q[i].Pos[2] = *(float *) (pos->data + i*(pos->strides[0]) + 2*pos->strides[1]);
SphQ[i].Hsml = *(float *) (hsml->data + i*(hsml->strides[0]));
}
/* now, give observable value for P */
for (i = 0; i < NumPart; i++)
{
SphP[i].Observable = *(float *) (obs->data + i*(obs->strides[0]));
}
sph_sub();
for (i = 0; i < pos->dimensions[0]; i++)
{
*(float *)(vobs->data + i*(vobs->strides[0])) = SphQ[i].Observable;
}
free(Q);
free(SphQ);
return Py_BuildValue("O",vobs);
}
static PyObject * gadget_Ngbs(PyObject* self, PyObject *args)
{
PyArrayObject *pos;
float eps;
if (! PyArg_ParseTuple(args, "Of",&pos,&eps))
return PyString_FromString("error");
PyArrayObject *poss;
int i,j,n,nn;
int input_dimension;
size_t bytes;
int startnode,numngb;
int phase=0;
FLOAT searchcenter[3];
double dx,dy,dz,r2,eps2;
input_dimension =pos->nd;
if (input_dimension != 1)
PyErr_SetString(PyExc_ValueError,"dimension of first argument must be 1");
pos = TO_FLOAT(pos);
eps2 = eps*eps;
searchcenter[0] = (FLOAT)*(float *) (pos->data + 0*(pos->strides[0]));
searchcenter[1] = (FLOAT)*(float *) (pos->data + 1*(pos->strides[0]));
searchcenter[2] = (FLOAT)*(float *) (pos->data + 2*(pos->strides[0]));
startnode = All.MaxPart;
/* ici, il faut faire une fct qui fonctionne en //, cf hydra --> Exportflag */
numngb = ngb_treefind_pairs(&searchcenter[0], (FLOAT)eps, phase, &startnode);
nn=0;
for(n = 0;n < numngb; n++)
{
j = Ngblist[n];
dx = searchcenter[0] - P[j].Pos[0];
dy = searchcenter[1] - P[j].Pos[1];
dz = searchcenter[2] - P[j].Pos[2];
r2 = dx * dx + dy * dy + dz * dz;
if (r2<=eps2)
{
printf("%d r=%g\n",nn,sqrt(r2));
nn++;
}
}
return PyArray_Return(pos);
}
static PyObject * gadget_SphEvaluateOrigAll(PyObject* self, PyObject *args)
{
PyArrayObject *obs;
if (! PyArg_ParseTuple(args, "O",&obs))
return PyString_FromString("error");
int i;
size_t bytes;
int ld[1];
PyArrayObject *vobs;
if (obs->nd != 1)
PyErr_SetString(PyExc_ValueError,"dimension of obs must be 1.");
if (obs->dimensions[0] != NumPart)
PyErr_SetString(PyExc_ValueError,"The size of obs must be NumPart.");
obs = TO_FLOAT(obs);
/* create a NumPy object */
ld[0]=NumPart;
vobs = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_FLOAT);
/* flag all particles as active */
for (i = 0; i < NumPart; i++)
P[i].Ti_endstep == All.Ti_Current;
/* now, give observable value for P */
for (i = 0; i < NumPart; i++)
{
SphP[i].Observable = *(float *) (obs->data + i*(obs->strides[0]));
}
sph_orig();
for (i = 0; i < NumPart; i++)
{
*(float *)(vobs->data + i*(vobs->strides[0])) = SphP[i].Observable;
}
return Py_BuildValue("O",vobs);
}
static PyObject * gadget_SphEvaluateAll(PyObject* self, PyObject *args)
{
PyArrayObject *obs;
if (! PyArg_ParseTuple(args, "O",&obs))
return PyString_FromString("error");
int i;
size_t bytes;
int ld[1];
PyArrayObject *vobs;
if (obs->nd != 1)
PyErr_SetString(PyExc_ValueError,"dimension of obs must be 1.");
if (obs->dimensions[0] != NumPart)
PyErr_SetString(PyExc_ValueError,"The size of obs must be NumPart.");
obs = TO_FLOAT(obs);
/* create a NumPy object */
ld[0]=NumPart;
vobs = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_FLOAT);
/* flag all particles as active */
for (i = 0; i < NumPart; i++)
P[i].Ti_endstep == All.Ti_Current;
/* now, give observable value for P */
for (i = 0; i < NumPart; i++)
{
SphP[i].Observable = *(float *) (obs->data + i*(obs->strides[0]));
}
sph();
for (i = 0; i < NumPart; i++)
{
*(float *)(vobs->data + i*(vobs->strides[0])) = SphP[i].Observable;
}
return Py_BuildValue("O",vobs);
}
static PyObject * gadget_SphEvaluateGradientAll(PyObject* self, PyObject *args)
{
PyArrayObject *obs;
if (! PyArg_ParseTuple(args, "O",&obs))
return PyString_FromString("error");
int i;
size_t bytes;
npy_intp ld[2];
PyArrayObject *grad;
if (obs->nd != 1)
PyErr_SetString(PyExc_ValueError,"dimension of obs must be 1.");
if (obs->dimensions[0] != NumPart)
PyErr_SetString(PyExc_ValueError,"The size of obs must be NumPart.");
obs = TO_FLOAT(obs);
All.Ti_Current=1; /* need to flag active particles */
/* now, give observable value for P */
for (i = 0; i < NumPart; i++)
{
SphP[i].Observable = *(float *) (obs->data + i*(obs->strides[0]));
}
sph();
/* create a NumPy object */
ld[0] = NumPart;
ld[1] = 3;
grad = (PyArrayObject *) PyArray_SimpleNew(2,ld,PyArray_FLOAT);
for (i = 0; i < NumPart; i++)
{
*(float *) (grad->data + i*(grad->strides[0]) + 0*grad->strides[1]) = SphP[i].GradObservable[0];
*(float *) (grad->data + i*(grad->strides[0]) + 1*grad->strides[1]) = SphP[i].GradObservable[1];
*(float *) (grad->data + i*(grad->strides[0]) + 2*grad->strides[1]) = SphP[i].GradObservable[2];
}
return PyArray_Return(grad);
}
static PyObject * gadget_DensityEvaluateAll(PyObject* self, PyObject *args)
{
int i;
/* flag all particles as active */
//All.Ti_Current=1; /* need to flag active particles */
for (i = 0; i < NumPart; i++)
P[i].Ti_endstep == All.Ti_Current;
density(0);
return Py_BuildValue("i",1);
}
static PyObject * gadget_DensityEvaluateGradientAll(PyObject* self, PyObject *args)
{
PyArrayObject *obs;
if (! PyArg_ParseTuple(args, "O",&obs))
return PyString_FromString("error");
int i;
size_t bytes;
npy_intp ld[2];
PyArrayObject *grad;
if (obs->nd != 1)
PyErr_SetString(PyExc_ValueError,"dimension of obs must be 1.");
if (obs->dimensions[0] != NumPart)
PyErr_SetString(PyExc_ValueError,"The size of obs must be NumPart.");
obs = TO_FLOAT(obs);
//All.Ti_Current=1; /* need to flag active particles */
for (i = 0; i < NumPart; i++)
P[i].Ti_endstep == All.Ti_Current;
/* now, give observable value for P */
for (i = 0; i < NumPart; i++)
{
SphP[i].Observable = *(float *) (obs->data + i*(obs->strides[0]));
}
density_sph_gradient();
/* create a NumPy object */
ld[0] = NumPart;
ld[1] = 3;
grad = (PyArrayObject *) PyArray_SimpleNew(2,ld,PyArray_FLOAT);
for (i = 0; i < NumPart; i++)
{
*(float *) (grad->data + i*(grad->strides[0]) + 0*grad->strides[1]) = SphP[i].GradObservable[0];
*(float *) (grad->data + i*(grad->strides[0]) + 1*grad->strides[1]) = SphP[i].GradObservable[1];
*(float *) (grad->data + i*(grad->strides[0]) + 2*grad->strides[1]) = SphP[i].GradObservable[2];
}
return PyArray_Return(grad);
}
static PyObject * gadget_EvaluateThermalConductivity(PyObject* self, PyObject *args)
{
PyArrayObject *obs;
if (! PyArg_ParseTuple(args, "O",&obs))
return PyString_FromString("error");
int i;
size_t bytes;
npy_intp ld[2];
PyArrayObject *grad;
if (obs->nd != 1)
PyErr_SetString(PyExc_ValueError,"dimension of obs must be 1.");
if (obs->dimensions[0] != NumPart)
PyErr_SetString(PyExc_ValueError,"The size of obs must be NumPart.");
obs = TO_FLOAT(obs);
//All.Ti_Current=1; /* need to flag active particles */
for (i = 0; i < NumPart; i++)
P[i].Ti_endstep == All.Ti_Current;
/* now, give observable value for P */
for (i = 0; i < NumPart; i++)
{
SphP[i].Observable = *(float *) (obs->data + i*(obs->strides[0]));
}
density_sph_gradient();
// equivalent of SphP[i].GradEnergyInt[0] is now in
//SphP[i].GradObservable[0]
//SphP[i].GradObservable[1]
//SphP[i].GradObservable[2]
/* now, compute thermal conductivity */
sph_thermal_conductivity();
/* create a NumPy object */
ld[0] = NumPart;
ld[1] = 3;
grad = (PyArrayObject *) PyArray_SimpleNew(2,ld,PyArray_FLOAT);
for (i = 0; i < NumPart; i++)
{
*(float *) (grad->data + i*(grad->strides[0]) + 0*grad->strides[1]) = SphP[i].GradObservable[0];
*(float *) (grad->data + i*(grad->strides[0]) + 1*grad->strides[1]) = SphP[i].GradObservable[1];
*(float *) (grad->data + i*(grad->strides[0]) + 2*grad->strides[1]) = SphP[i].GradObservable[2];
}
return PyArray_Return(grad);
}
static PyObject *gadget_LoadParticles2(PyObject *self, PyObject *args, PyObject *kwds)
{
int i,j;
size_t bytes;
PyArrayObject *ntype,*pos,*vel,*mass,*num,*tpe;
static char *kwlist[] = {"npart", "pos","vel","mass","num","tpe", NULL};
if (! PyArg_ParseTupleAndKeywords(args, kwds, "|OOOOOO",kwlist,&ntype,&pos,&vel,&mass,&num,&tpe))
return Py_BuildValue("i",1);
/* check type */
if (!(PyArray_Check(pos)))
{
PyErr_SetString(PyExc_ValueError,"aruments 1 must be array.");
return NULL;
}
/* check type */
if (!(PyArray_Check(mass)))
{
PyErr_SetString(PyExc_ValueError,"aruments 2 must be array.");
return NULL;
}
/* check dimension */
if ( (pos->nd!=2))
{
PyErr_SetString(PyExc_ValueError,"Dimension of argument 1 must be 2.");
return NULL;
}
/* check dimension */
if ( (mass->nd!=1))
{
PyErr_SetString(PyExc_ValueError,"Dimension of argument 2 must be 1.");
return NULL;
}
/* check size */
if ( (pos->dimensions[1]!=3))
{
PyErr_SetString(PyExc_ValueError,"First size of argument must be 3.");
return NULL;
}
/* check size */
if ( (pos->dimensions[0]!=mass->dimensions[0]))
{
PyErr_SetString(PyExc_ValueError,"Size of argument 1 must be similar to argument 2.");
return NULL;
}
/* ensure double */
// ntype = TO_INT(ntype);
// pos = TO_FLOAT(pos);
// vel = TO_FLOAT(vel);
// mass = TO_FLOAT(mass);
// num = TO_FLOAT(num);
// tpe = TO_FLOAT(tpe);
/***************************************
* some inits *
/***************************************/
RestartFlag = 0;
Begrun1();
/***************************************
* LOAD PARTILES *
/***************************************/
NumPart = 0;
N_gas = *(int*) (ntype->data + 0*(ntype->strides[0]));
for (i = 0; i < 6; i++)
NumPart += *(int*) (ntype->data + i*(ntype->strides[0]));
if (NumPart!=pos->dimensions[0])
{
PyErr_SetString(PyExc_ValueError,"Numpart != pos->dimensions[0].");
return NULL;
}
MPI_Allreduce(&NumPart, &All.TotNumPart, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
MPI_Allreduce(&N_gas, &All.TotN_gas, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
All.MaxPart = All.PartAllocFactor * (All.TotNumPart / NTask);
All.MaxPartSph = All.PartAllocFactor * (All.TotN_gas / NTask);
/*********************/
/* allocate P */
/*********************/
if(!(P = malloc(bytes = All.MaxPart * sizeof(struct particle_data))))
{
printf("failed to allocate memory for `P' (%g MB).\n", bytes / (1024.0 * 1024.0));
endrun(1);
}
if(!(SphP = malloc(bytes = All.MaxPartSph * sizeof(struct sph_particle_data))))
{
printf("failed to allocate memory for `SphP' (%g MB) %d.\n", bytes / (1024.0 * 1024.0), sizeof(struct sph_particle_data));
endrun(1);
}
/*********************/
/* init P */
/*********************/
float * fpt;
for (i = 0; i < NumPart; i++)
{
//P[i].Pos[0] = *(float *) (pos->data + i*(pos->strides[0]) + 0*pos->strides[1]);
//P[i].Pos[1] = *(float *) (pos->data + i*(pos->strides[0]) + 1*pos->strides[1]);
//P[i].Pos[2] = *(float *) (pos->data + i*(pos->strides[0]) + 2*pos->strides[1]);
//&P[i].Pos[0] = (float *) (pos->data + i*(pos->strides[0]) + 0*pos->strides[1]);
//&P[i].Pos[1] = (float *) (pos->data + i*(pos->strides[0]) + 1*pos->strides[1]);
//&P[i].Pos[2] = (float *) (pos->data + i*(pos->strides[0]) + 2*pos->strides[1]);
fpt = (float *) (pos->data + i*(pos->strides[0]) + 0*pos->strides[1]);
P[i].Vel[0] = *(float *) (vel->data + i*(vel->strides[0]) + 0*vel->strides[1]);
P[i].Vel[1] = *(float *) (vel->data + i*(vel->strides[0]) + 1*vel->strides[1]);
P[i].Vel[2] = *(float *) (vel->data + i*(vel->strides[0]) + 2*vel->strides[1]);
P[i].Mass = *(float *) (mass->data + i*(mass->strides[0]));
P[i].ID = *(unsigned int *) (num->data + i*(num->strides[0]));
P[i].Type = *(int *) (tpe->data + i*(tpe->strides[0]));
//P[i].Active = 1;
}
/***************************************
* END LOAD PARTILES *
/***************************************/
/***************************************
* finish inits *
/***************************************/
Begrun2();
return Py_BuildValue("i",1);
}
static PyObject *gadget_domain_Decomposition(PyObject *self, PyObject *args, PyObject *kwds)
{
All.NumForcesSinceLastDomainDecomp = 1 + All.TotNumPart * All.TreeDomainUpdateFrequency;
Flag_FullStep = 1; /* to ensure that Peano-Hilber order is done */
domain_Decomposition();
return Py_BuildValue("i",1);
}
static PyObject *gadget_gravity_tree(PyObject *self, PyObject *args, PyObject *kwds)
{
gravity_tree();
return Py_BuildValue("i",1);
}
#ifdef TESSEL
static PyObject *gadget_ConstructDelaunay(PyObject *self, PyObject *args, PyObject *kwds)
{
ConstructDelaunay();
return Py_BuildValue("i",1);
}
static PyObject *gadget_ComputeVoronoi(PyObject *self, PyObject *args, PyObject *kwds)
{
ComputeVoronoi();
return Py_BuildValue("i",1);
}
static PyObject *gadget_ComputeTesselAccelerations(PyObject *self, PyObject *args, PyObject *kwds)
{
tessel_compute_accelerations();
return Py_BuildValue("i",1);
}
static PyObject *gadget_InitvEntropy(PyObject *self, PyObject *args, PyObject *kwds)
{
tessel_convert_energy_to_entropy();
return Py_BuildValue("i",1);
}
static PyObject *gadget_tessel_drift(PyObject *self, PyObject *args, PyObject *kwds)
{
float dt;
if (!PyArg_ParseTuple(args,"f",&dt))
return NULL;
tessel_drift(dt);
return Py_BuildValue("i",1);
}
static PyObject *gadget_tessel_kick(PyObject *self, PyObject *args, PyObject *kwds)
{
float dt;
if (!PyArg_ParseTuple(args,"f",&dt))
return NULL;
tessel_kick(dt);
return Py_BuildValue("i",1);
}
static PyObject *gadget_tessel_get_timestep(PyObject *self, PyObject *args, PyObject *kwds)
{
double dt;
dt = tessel_get_timestep();
return Py_BuildValue("f",(float)dt);
}
+static PyObject *gadget_tessel_dump_triangles(PyObject *self, PyObject *args, PyObject *kwds)
+ {
+
+ char *filename;
+ if (!PyArg_ParseTuple(args, "s", &filename))
+ return NULL;
+
+ dump_triangles(filename);
+ return Py_BuildValue("i",1);
+ }
+
+static PyObject *gadget_tessel_dump_voronoi(PyObject *self, PyObject *args, PyObject *kwds)
+ {
+
+ char *filename;
+ if (!PyArg_ParseTuple(args, "s", &filename))
+ return NULL;
+
+ dump_voronoi(filename);
+ return Py_BuildValue("i",1);
+ }
#endif /*TESSEL*/
/* definition of the method table */
static PyMethodDef gadgetMethods[] = {
{"Info", (PyCFunction)gadget_Info, METH_VARARGS,
"give some info"},
{"InitMPI", (PyCFunction)gadget_InitMPI, METH_VARARGS,
"Init MPI"},
{"InitDefaultParameters", (PyCFunction)gadget_InitDefaultParameters, METH_VARARGS,
"Init default parameters"},
{"GetParameters", (PyCFunction)gadget_GetParameters, METH_VARARGS,
"get gadget parameters"},
{"SetParameters", (PyCFunction)gadget_SetParameters, METH_VARARGS,
"Set gadget parameters"},
{"check_parser", (PyCFunction)gadget_check_parser, METH_VARARGS|METH_KEYWORDS,
"check the parser"},
{"LoadParticles", (PyCFunction)gadget_LoadParticles, METH_VARARGS|METH_KEYWORDS,
"LoadParticles partilces"},
{"LoadParticlesQ", (PyCFunction)gadget_LoadParticlesQ, METH_VARARGS,
"LoadParticles partilces Q"},
{"LoadParticles2", (PyCFunction)gadget_LoadParticles2, METH_VARARGS,
"LoadParticles partilces"},
{"AllPotential", (PyCFunction)gadget_AllPotential, METH_VARARGS,
"Computes the potential for each particle"},
{"AllAcceleration", (PyCFunction)gadget_AllAcceleration, METH_VARARGS,
"Computes the gravitational acceleration for each particle"},
{"GetAllAcceleration", (PyCFunction)gadget_GetAllAcceleration, METH_VARARGS,
"get the gravitational acceleration for each particle"},
{"GetAllPotential", (PyCFunction)gadget_GetAllPotential, METH_VARARGS,
"get the potential for each particle"},
{"GetAllDensities", (PyCFunction)gadget_GetAllDensities, METH_VARARGS,
"get the densities for each particle"},
{"GetAllHsml", (PyCFunction)gadget_GetAllHsml, METH_VARARGS,
"get the sph smoothing length for each particle"},
{"GetAllPositions", (PyCFunction)gadget_GetAllPositions, METH_VARARGS,
"get the position for each particle"},
{"GetAllVelocities", (PyCFunction)gadget_GetAllVelocities, METH_VARARGS,
"get the velocities for each particle"},
{"GetAllMasses", (PyCFunction)gadget_GetAllMasses, METH_VARARGS,
"get the mass for each particle"},
{"GetAllEnergySpec", (PyCFunction)gadget_GetAllEnergySpec, METH_VARARGS,
"get the specific energy for each particle"},
{"GetAllEntropy", (PyCFunction)gadget_GetAllEntropy, METH_VARARGS,
"get the entropy for each particle"},
{"GetAllID", (PyCFunction)gadget_GetAllID, METH_VARARGS,
"get the ID for each particle"},
{"GetAllTypes", (PyCFunction)gadget_GetAllTypes, METH_VARARGS,
"get the type for each particle"},
{"GetPos", (PyCFunction)gadget_GetPos, METH_VARARGS,
"get the position for each particle (no memory overhead)"},
{"Potential", (PyCFunction)gadget_Potential, METH_VARARGS,
"get the potential for a givent sample of points"},
{"Acceleration", (PyCFunction)gadget_Acceleration, METH_VARARGS,
"get the acceleration for a givent sample of points"},
{"SphEvaluateOrigAll", (PyCFunction)gadget_SphEvaluateOrigAll, METH_VARARGS,
"run the original sph routine."},
{"SphEvaluateAll", (PyCFunction)gadget_SphEvaluateAll, METH_VARARGS,
"compute mean value of a given field based on the sph convolution for all points."},
{"SphEvaluateGradientAll", (PyCFunction)gadget_SphEvaluateGradientAll, METH_VARARGS,
"compute the gradient of a given field based on the sph convolution for all points."},
{"DensityEvaluateAll", (PyCFunction)gadget_DensityEvaluateAll, METH_VARARGS,
"simply run the density function"},
{"DensityEvaluateGradientAll", (PyCFunction)gadget_DensityEvaluateGradientAll, METH_VARARGS,
"run the modified density function and compute a gradient from the given value"},
{"EvaluateThermalConductivity", (PyCFunction)gadget_EvaluateThermalConductivity, METH_VARARGS,
"evaluate the thermal coductivity for each particle"},
{"SphEvaluate", (PyCFunction)gadget_SphEvaluate, METH_VARARGS,
"compute mean value based on the sph convolution for a given number of points"},
{"InitHsml", (PyCFunction)gadget_InitHsml, METH_VARARGS,
"Init hsml based on the three for a given number of points"},
{"Density", (PyCFunction)gadget_Density, METH_VARARGS,
"compute Density based on the three for a given number of points"},
{"Ngbs", (PyCFunction)gadget_Ngbs, METH_VARARGS,
"return the position of the neighbors for a given point"},
{"GetAllPositionsQ", (PyCFunction)gadget_GetAllPositionsQ, METH_VARARGS,
"get the position for each particle Q"},
{"GetAllVelocitiesQ", (PyCFunction)gadget_GetAllVelocitiesQ, METH_VARARGS,
"get the velocities for each particle Q"},
{"GetAllMassesQ", (PyCFunction)gadget_GetAllMassesQ, METH_VARARGS,
"get the mass for each particle Q"},
{"GetAllIDQ", (PyCFunction)gadget_GetAllIDQ, METH_VARARGS,
"get the ID for each particle Q"},
{"GetAllTypesQ", (PyCFunction)gadget_GetAllTypesQ, METH_VARARGS,
"get the type for each particle Q"},
{"Free", (PyCFunction)gadget_Free, METH_VARARGS,
"release memory"},
{"domain_Decomposition", (PyCFunction)gadget_domain_Decomposition, METH_VARARGS,
"call domain_Decomposition"},
{"gravity_tree", (PyCFunction)gadget_gravity_tree, METH_VARARGS,
"call gravity_tree"},
#ifdef TESSEL
{"ConstructDelaunay", (PyCFunction)gadget_ConstructDelaunay, METH_VARARGS,
"Construct the Delaunay tesselation"},
{"ComputeVoronoi", (PyCFunction)gadget_ComputeVoronoi, METH_VARARGS,
"Compute the Voronoi tesselation"},
{"ComputeTesselAccelerations", (PyCFunction)gadget_ComputeTesselAccelerations, METH_VARARGS,
"Compute the acceleration using the tesselation"},
{"InitvEntropy", (PyCFunction)gadget_InitvEntropy, METH_VARARGS,
"Initialize vEntropy using internal energy."},
{"tessel_drift", (PyCFunction)gadget_tessel_drift, METH_VARARGS,
"drift particles"},
{"tessel_kick", (PyCFunction)gadget_tessel_kick, METH_VARARGS,
"kick particles"},
{"tessel_get_timestep", (PyCFunction)gadget_tessel_get_timestep, METH_VARARGS,
"get timestep for tessel"},
+ {"tessel_dump_triangles", (PyCFunction)gadget_tessel_dump_triangles, METH_VARARGS,
+ "dump triangles"},
+
+ {"tessel_dump_voronoi", (PyCFunction)gadget_tessel_dump_voronoi, METH_VARARGS,
+ "dump voronoi"},
+
+
+
{"GetAllDelaunayTriangles", (PyCFunction)gadget_GetAllDelaunayTriangles, METH_VARARGS,
"Get all the Delaunay Triangles"},
{"GetAllvPoints", gadget_GetAllvPoints, METH_VARARGS,
"Get voronoi points"},
{"GetAllvDensities", gadget_GetAllvDensities, METH_VARARGS,
"Get voronoi density of all points"},
{"GetAllvVolumes", gadget_GetAllvVolumes, METH_VARARGS,
"Get voronoi volumes of all points"},
{"GetAllvPressures", gadget_GetAllvPressures, METH_VARARGS,
"Get voronoi presures of all points"},
{"GetAllvEnergySpec", gadget_GetAllvEnergySpec, METH_VARARGS,
"Get voronoi energyspec of all points"},
{"GetAllvAccelerations", gadget_GetAllvAccelerations, METH_VARARGS,
"Get voronoi accelerations of all points"},
{"GetvPointsForOnePoint", gadget_GetvPointsForOnePoint, METH_VARARGS,
"Get voronoi points for a given point"},
{"GetNgbPointsForOnePoint", gadget_GetNgbPointsForOnePoint, METH_VARARGS,
"Get neighbors points for a given point"},
{"GetNgbPointsAndFacesForOnePoint", gadget_GetNgbPointsAndFacesForOnePoint, METH_VARARGS,
"Get neighbors points and faces for a given point"},
{"GetAllGhostPositions", gadget_GetAllGhostPositions, METH_VARARGS,
"Get all positions of the ghosts points"},
{"GetAllGhostvDensities", gadget_GetAllGhostvDensities, METH_VARARGS,
"Get all densities of the ghosts points"},
{"GetAllGhostvVolumes", gadget_GetAllGhostvVolumes, METH_VARARGS,
- "Get all volumes of the ghosts points"},
+ "Get all volumes of the ghosts points"},
+
+
+
#endif
{NULL, NULL, 0, NULL} /* Sentinel */
};
void initgadget(void)
{
(void) Py_InitModule("gadget", gadgetMethods);
import_array();
}
#endif /*PY_INTERFACE*/
diff --git a/src/tessel.c b/src/tessel.c
index 632ecb6..aebbcf6 100644
--- a/src/tessel.c
+++ b/src/tessel.c
@@ -1,2819 +1,2985 @@
#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");
printf("Pt x=%g y=%g z=%g\n",Pt->Pos[0],Pt->Pos[1],Pt->Pos[2]);
endrun(-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);
//printf("Find triangle...\n");
iT= NewFindTriangle(Pt);
/* create the new triangles */
//printf("Split triangle...\n");
SplitTriangle(&Triangles[iT],Pt);
/* test the new triangles and divide and modify if necessary */
//printf("Do triangles in stack...\n");
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]);
//printf("check after %d\n",i);
CheckTriangles();
}
if (ThisTask==0)
printf("%d triangles\n",nT);
CheckTriangles();
printf("check ok\n");
}
/*! 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_Median=-1; /* index towards next median */
int number_of_vPoints=0;
int number_of_Medians=0;
int next;
int initialvPoint;
T0 = Tstart;
iP0 = iPstart;
initialPoint = T0->P[iP0];
initialPoint->ivPoint=nvPoints;
initialvPoint=nvPoints;
//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].Pos[2] = 0.5;
vPoints[nvPoints].next = nvPoints+1;
/* 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 */
nvPoints++;
number_of_vPoints++;
if (T1==Tstart) /* end of loop */
{
//printf(" end of loop\n");
initialPoint->nvPoints = number_of_vPoints;
if (nvPoints-2<initialvPoint)
{
printf("this is bad... please check.\n");
endrun(-72354);
}
initialPoint->ivPoint = nvPoints-2; /* the current point points towards the previous-1 one */
vPoints[nvPoints-1].next = initialvPoint; /* close the loop : the last points towards the first */
/* 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;
/* 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;
P[i].Pressure = (SphP[i].Entropy) * pow(P[i].Density, GAMMA); /* SphP[i].Entropy this is bad, as we use SphP */
}
/* do the ghost points now */
for(i=NumPart;i<NumPart+NumgPart;i++)
{
P[i].Volume = P[P[i].iPref].Volume;
P[i].Density= P[P[i].iPref].Density;
P[i].Pressure = P[P[i].iPref].Pressure;
}
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;
nT=0;
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)*2;
Pe[j].Pos[1] = tesselDomainCenter[1] + tesselDomainRadius * sin(2./3.*PI*j)*2;
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]);
if (ThisTask==0)
printf("Start local AddPoint in Delaunay.\n");
/* loop over all points */
for (i=0;i<NumPart;i++)
{
AddPoint(&P[i]);
//CheckPoints(i-1);
}
if (ThisTask==0)
printf("%d triangles\n",nT);
/* add ghost points */
CheckTriangles();
printf("check after AddPoint local is ok\n");
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");
}
/*! Compute accelerations
*/
void tessel_convert_energy_to_entropy()
{
#ifndef ISOTHERM_EQS
int i;
double a3=1;
if(All.ComovingIntegrationOn)
a3 = (All.Time * All.Time * All.Time);
else
a3 = 1;
if(header.flag_entropy_instead_u == 0)
for(i = 0; i < N_gas; i++)
{
P[i].Entropy = GAMMA_MINUS1 * SphP[i].Entropy / pow(P[i].Density / a3, GAMMA_MINUS1);
}
#endif
}
/*! Compute accelerations
*/
void tessel_compute_accelerations()
{
int i;
int np;
/* first triangle */
int iT,iTstart;
int p,iP0,iP1,iPnext0;
struct TriangleInList *T,*Tnext,*Tstart;
double acc[3];
double dx, dy, dz;
double Rij,Aij;
double eij[3],mpij[3],cmij[3],cij[3];
struct particle_data *Pi,*Pj;
struct vPoint *vP1,*vP2;
int nextvp1,nextvp2;
double PjpPi,PjmPi;
for(i=0;i<N_gas;i++)
{
/* some init */
acc[0] = acc[1] = acc[2] = 0;
Pi=&P[i];
np=0;
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(-712354);
iT = iTstart; /* index of triangle */
iP0 = p; /* index of the ref point "0" in the triangle */
T = Tstart;
np = 0;
nextvp1 = Pi->ivPoint;
/* visit all triangles around the point */
while (1)
{
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(-71235445);
return -1; /* need to increase the size for this point */
}
/* we have a ngb point */
np++;
Pj = T->P[iP1];
nextvp2 = vPoints[nextvp1].next;
vP1 = &vPoints[nextvp1];
vP2 = &vPoints[nextvp2];
dx = Pi->Pos[0] - Pj->Pos[0];
dy = Pi->Pos[1] - Pj->Pos[1];
dz = Pi->Pos[2] - Pj->Pos[2];
Rij = sqrt(dx * dx + dy * dy + dz * dz);
eij[0] = (Pj->Pos[0]-Pi->Pos[0])/Rij ;
eij[1] = (Pj->Pos[1]-Pi->Pos[1])/Rij ;
eij[2] = (Pj->Pos[2]-Pi->Pos[2])/Rij ;
/* mid point between i and j */
mpij[0] = 0.5*(Pj->Pos[0]+Pi->Pos[0]) ;
mpij[1] = 0.5*(Pj->Pos[1]+Pi->Pos[1]) ;
mpij[2] = 0.5*(Pj->Pos[2]+Pi->Pos[2]) ;
#ifdef TWODIMS
/* center of mass of the face */
cmij[0] = 0.5*(vP1->Pos[0]+vP2->Pos[0]);
cmij[1] = 0.5*(vP1->Pos[1]+vP2->Pos[1]);
cmij[2] = 0.5*(vP1->Pos[2]+vP2->Pos[2]);
/* area of the face */
dx = vP2->Pos[0]-vP1->Pos[0];
dy = vP2->Pos[1]-vP1->Pos[1];
dz = vP2->Pos[2]-vP1->Pos[2];
Aij = sqrt(dx * dx + dy * dy + dz * dz);
#else
need to compute mass center here
#endif
/* cij (mid point between ij, towards mass center of the face) cmij - mpij */
cij[0] = cmij[0]-mpij[0];
cij[1] = cmij[1]-mpij[1];
cij[2] = cmij[2]-mpij[2];
/* compute pressure */
PjpPi = Pj->Pressure + Pi->Pressure;
PjmPi = Pj->Pressure - Pi->Pressure;
/* now, the accelection */
acc[0] += -Aij*( 0.5*PjpPi*eij[0] + PjmPi*cij[0]/Rij);
acc[1] += -Aij*( 0.5*PjpPi*eij[1] + PjmPi*cij[1]/Rij);
acc[2] += -Aij*( 0.5*PjpPi*eij[2] + PjmPi*cij[2]/Rij);
/* 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(-7123544);
}
/* end of loop */
if (Tnext==Tstart)
break;
iP0 = iPnext0;
T = Tnext;
nextvp1 = nextvp2;
}
P[i].tesselAccel[0] = acc[0];
P[i].tesselAccel[1] = acc[1];
P[i].tesselAccel[2] = acc[2];
}
}
void tessel_kick(float dt_kick)
{
int i,j;
for(i = 0; i < NumPart; i++)
{
for(j = 0; j < 3; j++)
P[i].Vel[j] += P[i].tesselAccel[j] * dt_kick;
}
}
void tessel_drift(float dt_drift)
{
int i,j;
for(i = 0; i < NumPart; i++)
{
for(j = 0; j < 3; j++)
P[i].Pos[j] += P[i].Vel[j] * dt_drift;
}
}
double tessel_get_timestep()
{
double ax, ay, az, ac, csnd, hsml;
int i;
double dt,dtmin;
dtmin = 1e32;
for(i = 0; i < NumPart; i++)
{
//ax = P[i].tesselAccel[0];
//ay = P[i].tesselAccel[1];
//az = P[i].tesselAccel[2];
//ac = sqrt(ax * ax + ay * ay + az * az);
csnd = sqrt(GAMMA * P[i].Pressure / P[i].Density);
#ifdef TWODIMS
hsml = sqrt(P[i].Volume/PI);
#else
endrun("-787899");
#endif
dt = 2 * All.CourantFac * hsml / csnd;
if (dt<dtmin)
dtmin = dt;
}
return dtmin;
}
+/************************************************************/
+/* */
+/* IO ROUTINES */
+/* */
+/************************************************************/
+
+
+#define SKIP {my_fwrite(&blksize,sizeof(int),1,fd);}
+
+
+
+
+
+
+void dump_triangles(char *fname)
+ {
+
+ int iT;
+ FILE *fd = 0;
+ int blksize;
+ double f;
+ int idx;
+
+ //sprintf(fname,"%s%s_%04d","./","tria",1);
+
+
+ if(!(fd = fopen(fname, "w")))
+ {
+ printf("can't open file `%s' for writing triangles.\n", fname);
+ }
+ else
+ {
+ /* number of triangles */
+ blksize=sizeof(int);
+ SKIP;
+ my_fwrite(&nT, sizeof(nT), 1, fd);
+ SKIP;
+
+
+ /* triangles */
+ blksize=sizeof(f)*nT*3*3;
+ SKIP;
+
+ for (iT=0;iT<nT;iT++)
+ {
+ f = Triangles[iT].P[0]->Pos[0];
+ my_fwrite(&f, sizeof(f), 1, fd);
+ f = Triangles[iT].P[0]->Pos[1];
+ my_fwrite(&f, sizeof(f), 1, fd);
+ f = Triangles[iT].P[0]->Pos[2];
+ my_fwrite(&f, sizeof(f), 1, fd);
+
+ f = Triangles[iT].P[1]->Pos[0];
+ my_fwrite(&f, sizeof(f), 1, fd);
+ f = Triangles[iT].P[1]->Pos[1];
+ my_fwrite(&f, sizeof(f), 1, fd);
+ f = Triangles[iT].P[1]->Pos[2];
+ my_fwrite(&f, sizeof(f), 1, fd);
+
+ f = Triangles[iT].P[2]->Pos[0];
+ my_fwrite(&f, sizeof(f), 1, fd);
+ f = Triangles[iT].P[2]->Pos[1];
+ my_fwrite(&f, sizeof(f), 1, fd);
+ f = Triangles[iT].P[2]->Pos[2];
+ my_fwrite(&f, sizeof(f), 1, fd);
+ }
+
+ SKIP;
+
+ /* triangles */
+ blksize=sizeof(idx)*nT;
+ SKIP;
+
+ for (iT=0;iT<nT;iT++)
+ {
+ idx = Triangles[iT].idx;
+ my_fwrite(&idx, sizeof(idx), 1, fd);
+ }
+
+ SKIP;
+
+ close(fd);
+ }
+ }
+
+
+
+
+void dump_voronoi(char *fname)
+ {
+
+ int i,j;
+ FILE *fd = 0;
+ int blksize;
+ double f;
+ int np;
+ int n;
+
+ int next;
+
+
+ //sprintf(fname,"%s%s_%04d","./","tria",1);
+
+
+ if(!(fd = fopen(fname, "w")))
+ {
+ printf("can't open file `%s' for writing voronoi.\n", fname);
+ }
+ else
+ {
+
+
+ n = NumPart+NumgPart;
+
+ /* number of points */
+ blksize=sizeof(int);
+ SKIP;
+ my_fwrite(&n, sizeof(n), 1, fd);
+ SKIP;
+
+
+
+ for (i=0;i<n;i++)
+ {
+
+
+ next = P[i].ivPoint;
+ np = P[i].nvPoints;
+
+ /* number of points in the cell */
+ blksize=sizeof(int);
+ SKIP;
+ my_fwrite(&np, sizeof(np), 1, fd);
+ SKIP;
+
+
+ blksize=sizeof(f)*np*3;
+ SKIP;
+
+ for (j = 0; j < np; j++)
+ {
+
+ f = vPoints[next].Pos[0];
+ my_fwrite(&f, sizeof(f), 1, fd);
+ f = vPoints[next].Pos[1];
+ my_fwrite(&f, sizeof(f), 1, fd);
+ f = 0.5;
+ my_fwrite(&f, sizeof(f), 1, fd);
+
+ next = vPoints[next].next;
+
+ }
+
+ SKIP;
+
+ }
+
+
+
+ close(fd);
+ }
+ }
+
+
/************************************************************/
/* */
/* 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.5;
}
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_GetAllvPressures(PyObject* self)
{
import_array(); /* needed but strange no ? */
PyArrayObject *pressure;
npy_intp ld[1];
int i;
ld[0] = NumPart;
pressure = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE);
for (i = 0; i < pressure->dimensions[0]; i++)
{
*(double *) (pressure->data + i*(pressure->strides[0])) = P[i].Pressure;
}
return PyArray_Return(pressure);
}
PyObject *gadget_GetAllvEnergySpec(PyObject* self)
{
import_array(); /* needed but strange no ? */
PyArrayObject *energy;
npy_intp ld[1];
int i;
double a3inv;
double EnergySpec;
ld[0] = NumPart;
energy = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_DOUBLE);
if(All.ComovingIntegrationOn)
{
a3inv = 1 / (All.Time * All.Time * All.Time);
}
else
a3inv = 1;
for (i = 0; i < energy->dimensions[0]; i++)
{
#ifdef ISOTHERM_EQS
*(double *) (energy->data + i*(energy->strides[0])) = P[i].Entropy;
#else
EnergySpec=P[i].Entropy / GAMMA_MINUS1 * pow(P[i].Density * a3inv, GAMMA_MINUS1);
if (EnergySpec<All.MinEgySpec)
*(double *) (energy->data + i*(energy->strides[0])) = All.MinEgySpec;
else
*(double *) (energy->data + i*(energy->strides[0])) = P[i].Entropy / GAMMA_MINUS1 * pow(P[i].Density * a3inv, GAMMA_MINUS1);
#endif
}
return PyArray_Return(energy);
}
PyObject *gadget_GetAllvAccelerations(PyObject* self)
{
import_array(); /* needed but strange no ? */
PyArrayObject *acc;
npy_intp ld[2];
int i;
ld[0] = NumPart;
ld[1] = 3;
acc = (PyArrayObject *) PyArray_SimpleNew(2,ld,PyArray_FLOAT);
for (i = 0; i < acc->dimensions[0]; i++)
{
*(float *) (acc->data + i*(acc->strides[0]) + 0*acc->strides[1]) = P[i].tesselAccel[0];
*(float *) (acc->data + i*(acc->strides[0]) + 1*acc->strides[1]) = P[i].tesselAccel[1];
*(float *) (acc->data + i*(acc->strides[0]) + 2*acc->strides[1]) = P[i].tesselAccel[2];
}
return PyArray_Return(acc);
}
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.5;
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_GetNgbPointsForOnePoint(self, args)
PyObject *self;
PyObject *args;
{
import_array(); /* needed but strange no ? */
PyArrayObject *pos;
npy_intp ld[2];
int i,k;
int np=0;
int d=1; /* 1=counterclockwise 2=clockwise */
/* first triangle */
int iT,iTstart;
int p,iP0,iP1,iPnext0;
struct TriangleInList *T,*Tnext,*Tstart;
if (!PyArg_ParseTuple(args,"i",&i))
return NULL;
/*
first loop to count np (number of ngbs)
*/
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(-712354);
iT = iTstart; /* index of triangle */
iP0 = p; /* index of the ref point "0" in the triangle */
T = Tstart;
np = 0;
/* visit all triangles around the point */
while (1)
{
iP1 = (iP0+d) % 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(-71235445);
return -1; /* need to increase the size for this point */
}
/* we have a ngb point */
np++;
/* now, find the index of the ref. point in the next triangle */
iPnext0 = (T->idxe[iP1] + d) % 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(-7123544);
}
/* end of loop */
if (Tnext==Tstart)
break;
iP0 = iPnext0;
T = Tnext;
}
ld[0] = np;
ld[1] = 3;
pos = (PyArrayObject *) PyArray_SimpleNew(2,ld,PyArray_FLOAT);
/*
second loop to record positions
*/
iTstart = P[i].iT;
Tstart=&Triangles[iTstart];
k=0;
/* 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(-712354);
iT = iTstart; /* index of triangle */
iP0 = p; /* index of the ref point "0" in the triangle */
T = Tstart;
np = 0;
/* visit all triangles around the point */
while (1)
{
iP1 = (iP0+d) % 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(-71235445);
return -1; /* need to increase the size for this point */
}
/* we have a ngb point */
if (T->P[iP1]->iPref==-1)
{
*(float *) (pos->data + k*(pos->strides[0]) + 0*pos->strides[1]) = T->P[iP1]->Pos[0];
*(float *) (pos->data + k*(pos->strides[0]) + 1*pos->strides[1]) = T->P[iP1]->Pos[1];
*(float *) (pos->data + k*(pos->strides[0]) + 2*pos->strides[1]) = 0;
}
else /* a ghost point */
{
*(float *) (pos->data + k*(pos->strides[0]) + 0*pos->strides[1]) = P[T->P[iP1]->iPref].Pos[0];
*(float *) (pos->data + k*(pos->strides[0]) + 1*pos->strides[1]) = P[T->P[iP1]->iPref].Pos[1];
*(float *) (pos->data + k*(pos->strides[0]) + 2*pos->strides[1]) = 0;
}
k++;
/* now, find the index of the ref. point in the next triangle */
iPnext0 = (T->idxe[iP1] + d) % 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(-7123544);
}
/* end of loop */
if (Tnext==Tstart)
break;
iP0 = iPnext0;
T = Tnext;
}
return PyArray_Return(pos);
}
PyObject *gadget_GetNgbPointsAndFacesForOnePoint(self, args)
PyObject *self;
PyObject *args;
{
import_array(); /* needed but strange no ? */
PyArrayObject *list;
PyArrayObject *elt;
PyArrayObject *pos1,*pos2,*pos3;
npy_intp ld[2];
int i;
int np=0;
int d=1; /* 1=counterclockwise 2=clockwise */
/* first triangle */
int iT,iTstart;
int p,iP0,iP1,iPnext0;
struct TriangleInList *T,*Tnext,*Tstart;
int nextv,nextv2;
if (!PyArg_ParseTuple(args,"i",&i))
return NULL;
list = PyList_New(0);
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(-712354);
iT = iTstart; /* index of triangle */
iP0 = p; /* index of the ref point "0" in the triangle */
T = Tstart;
np = 0;
nextv = P[i].ivPoint;
/* visit all triangles around the point */
while (1)
{
iP1 = (iP0+d) % 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(-71235445);
return -1; /* need to increase the size for this point */
}
/* we have a ngb point */
elt = PyTuple_New(3);
pos1 = PyTuple_New(3);
pos2 = PyTuple_New(3);
pos3 = PyTuple_New(3);
PyTuple_SetItem(pos1, (Py_ssize_t)0, PyFloat_FromDouble((double)T->P[iP1]->Pos[0]));
PyTuple_SetItem(pos1, (Py_ssize_t)1, PyFloat_FromDouble((double)T->P[iP1]->Pos[1]));
PyTuple_SetItem(pos1, (Py_ssize_t)2, PyFloat_FromDouble((double)0.5));
PyTuple_SetItem(pos2, (Py_ssize_t)0, PyFloat_FromDouble((double)vPoints[nextv].Pos[0]));
PyTuple_SetItem(pos2, (Py_ssize_t)1, PyFloat_FromDouble((double)vPoints[nextv].Pos[1]));
PyTuple_SetItem(pos2, (Py_ssize_t)2, PyFloat_FromDouble((double)0.5));
nextv2 = vPoints[nextv].next;
PyTuple_SetItem(pos3, (Py_ssize_t)0, PyFloat_FromDouble((double)vPoints[nextv2].Pos[0]));
PyTuple_SetItem(pos3, (Py_ssize_t)1, PyFloat_FromDouble((double)vPoints[nextv2].Pos[1]));
PyTuple_SetItem(pos3, (Py_ssize_t)2, PyFloat_FromDouble((double)0.5));
PyTuple_SetItem(elt,(Py_ssize_t)0,pos1);
PyTuple_SetItem(elt,(Py_ssize_t)1,pos2);
PyTuple_SetItem(elt,(Py_ssize_t)2,pos3);
PyList_Append(list,elt);
nextv = vPoints[nextv].next;
/* now, find the index of the ref. point in the next triangle */
iPnext0 = (T->idxe[iP1] + d) % 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(-7123544);
}
/* end of loop */
if (Tnext==Tstart)
break;
iP0 = iPnext0;
T = Tnext;
}
return Py_BuildValue("O",list);
}
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