Commit 99596333 authored by jclamber's avatar jclamber
Browse files

daily commit

git-svn-id: http://svn.oamp.fr/repos/uns_projects/trunk@98 f264a43e-d52d-4b82-913b-c2bd5215a18a
parent 06576e2d
......@@ -25,26 +25,23 @@ extern "C" {
using namespace uns_proj;
using namespace jclut;
using namespace std;
using namespace uns;
//------------------------------------------------------------------------------
// PUBLIC CRectify
CRectify::CRectify(uns::CunsIn * _uns)
{
unsin = _uns;
init();
}
//------------------------------------------------------------------------------
//
CRectify::CRectify()
{
init();
}
//------------------------------------------------------------------------------
//
bool CRectify::rectify(const int _nbody,const float _time,
float * _pos, float * _vel,float * _mass, float * _rho,
const bool _use_rho, const bool _rectify,
std::string _cod_file,std::string _rect_file, const float _threshold)
{
init();
nbody= _nbody;
pos = _pos;
......@@ -60,97 +57,83 @@ bool CRectify::rectify(const int _nbody,const float _time,
rect = _rectify;
use_rho = _use_rho;
process2();
}
bool CRectify::rectify(int _nbody, float * _pos, float * _vel,float * _mass, float * _rho,
float _threshold, std::string _cod_file)
{
nbody= _nbody;
pos = _pos;
vel = _vel;
mass = _mass;
rho = _rho;
threshold = _threshold;
cod_file = _cod_file;
process();
}
bool CRectify::rectify(std::string component,float threshold, std::string cod_file)
{
}
// PUBLIC CDataIndex
//------------------------------------------------------------------------------
//
bool CDataIndex::sortData(const CDataIndex &a, const CDataIndex &b)
{
return a.data < b.data;
}
// PRIVATE
//------------------------------------------------------------------------------
//
void CRectify::init()
{
density=NULL;
nbody=0;
initOldFrame();
}
void CRectify::process2()
//------------------------------------------------------------------------------
//
void CRectify::process()
{
if (cod_file.length()!=0 &&!CSnaptools::isFileExist(cod_file)) {
std::cerr << "Requested codfile["<<codf<<"] does not exist, abort\n";
std::cerr << "\nRequested codfile["<<cod_file<<"] does not exist, abort\n\n";
std::exit(1);
}
if (density) delete density;
if (use_rho) {
processRho();
}
findCenter();
findMoment();
snapTransform();
}
void CRectify::process()
{
if (density) delete density;
findCenter();
findMoment();
snapTransform();
computeVectors();
if (rect) {
snapTransform();
}
}
//------------------------------------------------------------------------------
//
void CRectify::findCenter()
{
bool is_codfile=false;
if (cod_file.length()!=0) { // read cod from file
bool ok=CSnaptools::getTimeDataFile(cod_file,time,6,cod);
if (! ok ) {
is_codfile=CSnaptools::getTimeDataFile(cod_file,time,6,cod);
if (! is_codfile ) {
std::cerr << "From file ["<<cod_file<<"] unable to find corresponding time ["<<time<<"]"
<< "aborting program....\n";
std::exit(1);
}
} else {
if (use_rho) {
}
if (use_rho) { // compute density, and COD from rho
// compute density
processRho();
// compute cod if ! codfile
if (! is_codfile)
CSnaptools::moveToCod<float>(vmass.size(),&vpos[0],&vvel[0],&vmass[0],&vrho[0],cod,false); // compute cod don't move
}
}
else {
std::cerr << "WE should be here.....FILE" << __FILE__<<" Line : "
<<__LINE__<< "\n";
std::exit(1);
}
if (!is_codfile && !use_rho) {
std::cerr << "Program aborted, because you MUST provide a valid COD file or/and you MUST enable \"rho=t\" variable\n";
std::exit(1);
}
// cod double to float
for (int i=0; i<6; i++) {
codf[i] = (float) cod[i];
}
}
//------------------------------------------------------------------------------
//
void CRectify::processRho()
{
......@@ -212,36 +195,51 @@ void CRectify::processRho()
}
}
//------------------------------------------------------------------------------
//
void CRectify::findMoment()
{
CLRM(w_qpole);
float w_sum=0.0;
for (unsigned int i=0; i<vmass.size(); i++) {
w_sum += vrho[i] * vmass[i];
}
if ( use_rho) { // data are stored in Vectors
for (unsigned int i=0; i<vmass.size(); i++) {
w_sum += vrho[i] * vmass[i];
}
float tmpv[3], pos_b[3];
matrix tmpm;
for (unsigned int i=0; i<vmass.size(); i++) {
float w_b = vrho[i] * vmass[i];
vectutils::subv(pos_b, &vpos[i*3], codf); // shift to cod
MULVS(tmpv, pos_b, w_b);
OUTVP(tmpm, tmpv, pos_b);
ADDM(w_qpole, w_qpole, tmpm);
float tmpv[3], pos_b[3];
matrix tmpm;
for (unsigned int i=0; i<vmass.size(); i++) {
float w_b = vrho[i] * vmass[i];
vectutils::subv(pos_b, &vpos[i*3], codf); // shift to cod
MULVS(tmpv, pos_b, w_b);
OUTVP(tmpm, tmpv, pos_b);
ADDM(w_qpole, w_qpole, tmpm);
}
} else { // data are stored in arrays
for (unsigned int i=0; i<nbody; i++) {
w_sum += mass[i];
}
float tmpv[3], pos_b[3];
matrix tmpm;
for (unsigned int i=0; i<nbody; i++) {
float w_b = mass[i];
vectutils::subv(pos_b, &pos[i*3], codf); // shift to cod
MULVS(tmpv, pos_b, w_b);
OUTVP(tmpm, tmpv, pos_b);
ADDM(w_qpole, w_qpole, tmpm);
}
}
DIVMS(w_qpole, w_qpole, w_sum);
std::cerr << " find moment, done\n";
}
void CRectify::snapTransform()
//------------------------------------------------------------------------------
//
void CRectify::computeVectors()
{
float oldframe[][3] = {
{ 1.0, 0.0, 0.0, },
{ 0.0, 1.0, 0.0, },
{ 0.0, 0.0, 1.0, }};
float pos_b[3], vel_b[3];
//initOldFrame(); // !!!!!
// computing Eigen vectors according to moment w_qpole
eigenFrame(frame, w_qpole);
float tmp;
vectutils::dotvp(tmp,oldframe[0], frame[0]);
......@@ -254,9 +252,16 @@ void CRectify::snapTransform()
printvec("e_x:", frame[0]);
printvec("e_y:", frame[1]);
printvec("e_z:", frame[2]);
//for (b = btab; b < btab+nbody; b++) {
// !!!!
// update old frame vectors
for (int i = 0; i < NDIM; i++)
vectutils::setv(oldframe[i], frame[i]);
}
//------------------------------------------------------------------------------
//
void CRectify::snapTransform()
{
float pos_b[3], vel_b[3];
// LOOP on all particles not only those from threshold !!!
for (unsigned int i=0; i<nbody; i++) {
vectutils::subv(&pos[i*3], &pos[i*3], codf);
......@@ -270,10 +275,10 @@ void CRectify::snapTransform()
vectutils::setv(&vel[i*3], vel_b);
//SETV(Acc(b), acc_b);
}
for (int i = 0; i < NDIM; i++)
vectutils::setv(oldframe[i], frame[i]);
}
}
//------------------------------------------------------------------------------
//
void CRectify::eigenFrame(float frame[3][3], matrix mat)
{
float **q, *d, **v;
......@@ -291,6 +296,8 @@ void CRectify::eigenFrame(float frame[3][3], matrix mat)
for (j = 1; j <= 3; j++)
frame[i-1][j-1] = v[j][i];
}
//------------------------------------------------------------------------------
//
void CRectify::printvec(std::string name, float vec[3])
{
double PI=acos(-1.);
......@@ -300,7 +307,8 @@ void CRectify::printvec(std::string name, float vec[3])
name.c_str(), rtp[0], vec[0], vec[1], vec[2],
rtp[1]*180.0/PI, rtp[2]*180.0/PI);
}
//------------------------------------------------------------------------------
//
void CRectify::xyz2rtp(float xyz[3], float rtp[3])
{
double PI=acos(-1.);
......
......@@ -17,7 +17,6 @@
#ifndef CRECTIFY_H
#define CRECTIFY_H
#include <string>
#include "uns.h"
#include "cfalcon.h"
#define SINGLEPREC
#define real float
......@@ -49,25 +48,27 @@ public:
class CRectify {
public:
CRectify(uns::CunsIn * uns);
CRectify();
bool rectify(std::string component,float theshold=70.0, std::string cod_file="");
bool rectify(int nbody,float * pos, float * vel,float * mass, float * rho,
float threshold=70.0, std::string cod_file="");
bool rectify(const int nbody,const float time,
float * pos, float * vel,float * mass, float * weight,
const bool _use_rho, const bool rectify=false,
std::string cod_file="",std::string rect_file="", const float threshold=0.0);
void process();
void process2();
void initOldFrame() {
float init_frame[][3] = {
{ 1.0, 0.0, 0.0, },
{ 0.0, 1.0, 0.0, },
{ 0.0, 0.0, 1.0, }};
memcpy(oldframe,init_frame,9*sizeof(float));
}
private:
uns::CunsIn * unsin;
jclut::CDensity * density;
float time, * pos, * vel, * mass, * rho;
float threshold;
float oldframe[3][3];
std::string cod_file, rect_file;
bool rect, use_rho;
int nbody;
......@@ -77,6 +78,7 @@ private:
void findCenter();
void findMoment();
void snapTransform();
void computeVectors();
void eigenFrame(float frame[3][3], matrix mat);
void xyz2rtp(float xyz[3], float rtp[3]);
void printvec(std::string name, float vec[3]);
......
......@@ -28,7 +28,7 @@ namespace jclut {
template <class T> static void moveToCod(const int nbody,T * pos, T * Vel, T * mass, T * rho, double cod[6], bool move, bool verbose=false);
template <class T> static void moveToCom(const int nbody,T * pos, T * mass, bool verbose=false);
static bool getTimeDataFile(std::string _file,const float time,const int n,double data[],const double offset=0.0000001,const bool verbose=false);
static bool getTimeDataFile(std::string _file,const float time,const int n,double data[],const double offset=0.0001,const bool verbose=false);
static bool isFileExist(std::string);
static std::string basename(const std::string);
static std::string dirname(const std::string);
......
......@@ -37,6 +37,7 @@ const char * defv[] = {
"codfile=\n use cod from file to recenter",
"threshold=70.0\n percentage of density kept >= threshold",
"times=all\n selected time",
"first=f\n add a trailing numbering to the first output file",
"verbose=f\n verbose on/off",
"VERSION=1.0\n compiled on <"__DATE__"> JCL ",
NULL,
......@@ -44,7 +45,7 @@ const char * defv[] = {
const char * usage="Rectify an UNS snapshot";
using namespace uns_proj;
using namespace std;
//------------------------------------------------------------------------------
// M A I N
//------------------------------------------------------------------------------
......@@ -63,24 +64,33 @@ int main(int argc, char ** argv )
std::string select_c (getparam ((char *) "select" ));
std::string select_t (getparam ((char *) "times" ));
float threshold(getdparam ((char *) "threshold"));
bool first (getbparam((char *) "first" ));
bool verbose = getbparam((char *) "verbose" );
std::cerr << "codfile["<<codf<<"] length="<< codf.length()<<"\n";
double cod[6];
bool okcodf=CSnaptools::getTimeDataFile(codf,5.,6,cod);
std::cerr << "ok ocdf=" << okcodf <<"\n";
for (int i=0;i<6;i++)
std::cerr << " " << cod[i];
std::cerr << "\n";
std::exit(1);
if (codf.length()==0 && !use_rho) {
std::cerr << "\nProgram aborted, because you MUST provide a valid COD file\n"
<<" \"codfile=....\", or/and you MUST enable \"rho=t\" variable\n\n";
std::exit(1);
}
bool one_file=false;
bool stop=false;
bool special_nemo=false;
if (outname=="-" || outname==".") special_nemo=true;
uns::CunsOut * unsout=NULL; // out object
bool first_out=true;
// Rectify object
CRectify * crectify = new CRectify();
// -----------------------------------------------
// instantiate a new UNS input object (for reading)
uns::CunsIn * uns = new uns::CunsIn(simname,select_c,select_t,verbose);
if (uns->isValid()) { // input file is known by UNS lib
while(uns->snapshot->nextFrame()) { // there is a new frame
int cpt=0;
while(uns->snapshot->nextFrame()&&!stop) { // there is a new frame
std::cerr << "Input file is of type :"<<uns->snapshot->getInterfaceType()<<"\n";
bool ok;
int cnbody,nbody;
......@@ -106,21 +116,49 @@ int main(int argc, char ** argv )
ok=uns->snapshot->getData("hsml",&cnbody,&hsml);
std::cerr << "nbody=" << nbody << " time="<<time <<"\n";
CRectify * crectify = new CRectify();
crectify->rectify(nbody,pos,vel,mass,rho, threshold);
uns::CunsOut * unsout = new uns::CunsOut(outname,"nemo",verbose);
unsout->snapshot->setData("time",time);
unsout->snapshot->setData("pos" ,nbody,pos,false);
unsout->snapshot->setData("mass" ,nbody,mass,false);
if (rho)
unsout->snapshot->setData("rho" ,nbody,rho ,false);
if (hsml)
unsout->snapshot->setData("hsml",nbody,hsml,false);
unsout->snapshot->save();
delete unsout;
std::cerr << "Snapshot ="<< uns->snapshot->getFileName() << "\n";
crectify->rectify(nbody,time,pos,vel,mass,rho,use_rho,true,codf,rectf, threshold);
if (outname.length()!=0) {
stringstream number;
number << cpt++;
std::string out_name=std::string(outname);;
if (! special_nemo) { // ! standard output && ! "."
if (one_file || (cpt==1 && !first)) {
out_name=std::string(outname);
if (one_file) stop = true; // do not continue
} else {
stringstream ss;
ss << std::string(outname) << "." << setw(5) << setfill('0') << number.str();
//out_name=std::string(outname)+"."+number.str();
out_name=ss.str();
}
// create a new UNS out object
unsout = new uns::CunsOut(out_name,"nemo",verbose);
} else {
if (first_out) {
first_out = false;
// instantiate only once unsout, because outname="-"
unsout = new uns::CunsOut(out_name,"nemo",verbose);
}
}
std::cerr << "output filename=["<<out_name<<"]\n";
unsout->snapshot->setData("time",time);
unsout->snapshot->setData("pos" ,nbody,pos,false);
unsout->snapshot->setData("mass" ,nbody,mass,false);
if (rho)
unsout->snapshot->setData("rho" ,nbody,rho ,false);
if (hsml)
unsout->snapshot->setData("hsml",nbody,hsml,false);
unsout->snapshot->save();
if (!special_nemo) {
delete unsout;
}
}
}
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment