222 lines
6.6 KiB
C
222 lines
6.6 KiB
C
|
|
#ifdef newc
|
|
#include <algorithm>
|
|
#include <functional>
|
|
#include <vector>
|
|
#include <cstring>
|
|
#include <iostream>
|
|
#include <iomanip>
|
|
#include <fstream>
|
|
#include <cstdlib>
|
|
#include <cstdio>
|
|
#include <string>
|
|
#include <cmath>
|
|
#include <strstream>
|
|
using namespace std;
|
|
#else
|
|
#include <iostream.h>
|
|
#include <iomanip.h>
|
|
#include <fstream.h>
|
|
#include <stdlib.h>
|
|
#include <stdio.h>
|
|
#include <string.h>
|
|
#include <math.h>
|
|
#endif
|
|
|
|
#include "TwoPunctures.h"
|
|
|
|
inline string &lTrim(string &ss)
|
|
{
|
|
string::iterator p = find_if(ss.begin(), ss.end(), not1(ptr_fun<int, int>(isspace)));
|
|
ss.erase(ss.begin(), p);
|
|
return ss;
|
|
}
|
|
inline string &rTrim(string &ss)
|
|
{
|
|
string::reverse_iterator p = find_if(ss.rbegin(), ss.rend(), not1(ptr_fun<int, int>(isspace)));
|
|
ss.erase(p.base(), ss.end());
|
|
return ss;
|
|
}
|
|
inline string &Trim(string &st)
|
|
{
|
|
lTrim(rTrim(st));
|
|
return st;
|
|
}
|
|
|
|
int parse_parts(string str, string &sgrp, string &skey, string &sval, int &ind)
|
|
{
|
|
int pos1, pos2;
|
|
string s0;
|
|
|
|
ind = 0;
|
|
|
|
// remove comments
|
|
str = str.substr(0, str.find("#"));
|
|
if (rTrim(str).empty())
|
|
return 0; // continue;
|
|
|
|
// parse {group, key, val}
|
|
pos1 = str.find("::");
|
|
pos2 = str.find("=");
|
|
if (pos1 == string::npos || pos2 == string::npos)
|
|
return -1;
|
|
|
|
s0 = str.substr(0, pos1);
|
|
sgrp = lTrim(s0);
|
|
s0 = str.substr(pos1 + 2, pos2 - pos1 - 2);
|
|
skey = rTrim(s0);
|
|
s0 = str.substr(pos2 + 1);
|
|
sval = Trim(s0);
|
|
|
|
pos1 = sval.find("\"");
|
|
pos2 = sval.rfind("\"");
|
|
if (pos1 != string::npos)
|
|
{
|
|
sval = sval.substr(1, pos2 - 1);
|
|
}
|
|
|
|
pos1 = skey.find("[");
|
|
pos2 = skey.find("]");
|
|
if (pos1 != string::npos)
|
|
{
|
|
s0 = skey.substr(0, pos1);
|
|
ind = atoi(skey.substr(pos1 + 1, pos2 - pos1 - 1).c_str());
|
|
skey = s0;
|
|
}
|
|
|
|
return 1;
|
|
}
|
|
//=======================================
|
|
int main(int argc, char *argv[])
|
|
{
|
|
double mp, mm, b, Mp, Mm, admtol, Newtontol;
|
|
int nA, nB, nphi, Newtonmaxit;
|
|
double P_plusx, P_plusy, P_plusz;
|
|
double P_minusx, P_minusy, P_minusz;
|
|
double S_plusx, S_plusy, S_plusz;
|
|
double S_minusx, S_minusy, S_minusz;
|
|
// read parameter from file
|
|
{
|
|
const int LEN = 256;
|
|
char pline[LEN];
|
|
string str, sgrp, skey, sval;
|
|
int sind;
|
|
const char pname[] = "TwoPunctureinput.par";
|
|
ifstream inf(pname, ifstream::in);
|
|
if (!inf.good())
|
|
{
|
|
cout << "Can not open parameter file " << pname << endl;
|
|
exit(0);
|
|
}
|
|
|
|
for (int i = 1; inf.good(); i++)
|
|
{
|
|
inf.getline(pline, LEN);
|
|
str = pline;
|
|
|
|
int status = parse_parts(str, sgrp, skey, sval, sind);
|
|
if (status == -1)
|
|
{
|
|
cout << "error reading parameter file " << pname << " in line " << i << endl;
|
|
exit(0);
|
|
}
|
|
else if (status == 0)
|
|
continue;
|
|
// we assume input in Brugmann's convention
|
|
if (sgrp == "ABE")
|
|
{
|
|
if (skey == "mm")
|
|
mm = atof(sval.c_str());
|
|
else if (skey == "mp")
|
|
mp = atof(sval.c_str());
|
|
else if (skey == "b")
|
|
b = atof(sval.c_str());
|
|
else if (skey == "P_plusx")
|
|
P_plusy = -atof(sval.c_str());
|
|
else if (skey == "P_plusy")
|
|
P_plusx = atof(sval.c_str());
|
|
else if (skey == "P_plusz")
|
|
P_plusz = atof(sval.c_str());
|
|
else if (skey == "P_minusx")
|
|
P_minusy = -atof(sval.c_str());
|
|
else if (skey == "P_minusy")
|
|
P_minusx = atof(sval.c_str());
|
|
else if (skey == "P_minusz")
|
|
P_minusz = atof(sval.c_str());
|
|
else if (skey == "S_plusx")
|
|
S_plusy = -atof(sval.c_str());
|
|
else if (skey == "S_plusy")
|
|
S_plusx = atof(sval.c_str());
|
|
else if (skey == "S_plusz")
|
|
S_plusz = atof(sval.c_str());
|
|
else if (skey == "S_minusx")
|
|
S_minusy = -atof(sval.c_str());
|
|
else if (skey == "S_minusy")
|
|
S_minusx = atof(sval.c_str());
|
|
else if (skey == "S_minusz")
|
|
S_minusz = atof(sval.c_str());
|
|
else if (skey == "Mp")
|
|
Mp = atof(sval.c_str());
|
|
else if (skey == "Mm")
|
|
Mm = atof(sval.c_str());
|
|
else if (skey == "admtol")
|
|
admtol = atof(sval.c_str());
|
|
else if (skey == "Newtontol")
|
|
Newtontol = atof(sval.c_str());
|
|
else if (skey == "nA")
|
|
nA = atoi(sval.c_str());
|
|
else if (skey == "nB")
|
|
nB = atoi(sval.c_str());
|
|
else if (skey == "nphi")
|
|
nphi = atoi(sval.c_str());
|
|
else if (skey == "Newtonmaxit")
|
|
Newtonmaxit = atoi(sval.c_str());
|
|
}
|
|
}
|
|
inf.close();
|
|
}
|
|
// echo parameters
|
|
{
|
|
cout << "///////////////////////////////////////////////////////////////" << endl;
|
|
cout << " mp = " << mp << endl;
|
|
cout << " mm = " << mm << endl;
|
|
cout << " b = " << b << endl;
|
|
cout << " P_plusx = " << P_plusx << endl;
|
|
cout << " P_plusy = " << P_plusy << endl;
|
|
cout << " P_plusz = " << P_plusz << endl;
|
|
cout << " P_minusx = " << P_minusx << endl;
|
|
cout << " P_minusy = " << P_minusy << endl;
|
|
cout << " P_minusz = " << P_minusz << endl;
|
|
cout << " S_plusx = " << S_plusx << endl;
|
|
cout << " S_plusy = " << S_plusy << endl;
|
|
cout << " S_plusz = " << S_plusz << endl;
|
|
cout << " S_minusx = " << S_minusx << endl;
|
|
cout << " S_minusy = " << S_minusy << endl;
|
|
cout << " S_minusz = " << S_minusz << endl;
|
|
cout << " Mp = " << Mp << endl;
|
|
cout << " Mm = " << Mm << endl;
|
|
cout << " admtol = " << admtol << endl;
|
|
cout << " Newtontol = " << Newtontol << endl;
|
|
cout << " nA = " << nA << endl;
|
|
cout << " nB = " << nB << endl;
|
|
cout << " nphi = " << nphi << endl;
|
|
cout << "Newtonmaxit = " << Newtonmaxit << endl;
|
|
cout << "///////////////////////////////////////////////////////////////" << endl;
|
|
}
|
|
//===========================the computation body====================================================
|
|
TwoPunctures *ADM;
|
|
|
|
ADM = new TwoPunctures(mp, mm, b, P_plusx, P_plusy, P_plusz, S_plusx, S_plusy, S_plusz,
|
|
P_minusx, P_minusy, P_minusz, S_minusx, S_minusy, S_minusz,
|
|
nA, nB, nphi, Mp, Mm, admtol, Newtontol, Newtonmaxit);
|
|
ADM->Solve();
|
|
ADM->Save("Ansorg.psid");
|
|
|
|
delete ADM;
|
|
//=======================caculation done=============================================================
|
|
cout << "===============================================================" << endl;
|
|
cout << "Initial data is successfully producede!!" << endl;
|
|
|
|
exit(0);
|
|
}
|