Skip to content

GitLab

Commit bbd76592 authored by Jose Manuel Morgado Chávez's avatar Jose Manuel Morgado Chávez
Browse files

Cleaning codes

parent 235a22d4
Hide whitespace changes
Inline Side-by-side
include/NumA/RadonTransform/RadonTransform.h
View file @ bbd76592
......@@ -82,14 +82,13 @@ public:
*/
Eigen::MatrixXd matrix_assembly( const std::vector<double> &x, const std::vector<double> &y, const std::vector<double> &xi);
/**
* Compute DD
*/
/* ---------------------------------------------------------------------------------------------------------------------------------------------------------------- */
// Compute DD
Eigen::VectorXd computeDD( const Eigen::VectorXd &GPD );
/**
* Compute GPD
*/
// Compute GPD
double computeGPD( const Eigen::VectorXd &DD, const double x, const double xi );
double computeGPD( const Eigen::VectorXd &DD, const double x, const double xi, const std::string& filename );
......@@ -98,9 +97,7 @@ std::vector<double> computeGPD( const Eigen::VectorXd &DD, const std::vector<dou
std::vector<double> computeGPD( const Eigen::VectorXd &DD, const std::vector<double> &xvec, const std::vector<double> &xivec, const std::string& filename );
/**
* Compute Dterm contribution.
*/
// Compute Dterm contribution.
double computeDterm( const Eigen::VectorXd &DD, const double x, const double xi );
};
......
include/NumA/triangulation/mesh.h
View file @ bbd76592
......@@ -49,7 +49,7 @@ class Mesh
// --------------------------------------------------
/* Set minimum area for elements */
void SetMinimalArea( float area );
void SetMaximalArea( float area );
/* Creation of the mesh */
void GenerateMesh();
......
src/NumA/RadonTransform/RadonTransform.cpp
View file @ bbd76592
......@@ -44,173 +44,23 @@ RadonTransform::~RadonTransform()
void RadonTransform::init()
{
// --------------------------------------------------------------------------------------------
// =============================================================================================
// =============================================================================================
// --------------------------------------------------------------------------------------------
// Set mesh (Proper implementation)
// --------------------------------------------------------------------------------------------
// mesh.SetMinimalArea ( 0.001 ); // Set minimal area for mesh elements.
// mesh.GenerateMesh(); // Generate triangulation for object mesh.
// mesh.Nodes( 1 ); // Nodes for a constant interpolating polynomial.
// mesh.BasisFunctions( 0 ); // Definition of interpolating basis functions.
// mesh.Nodes( 1 ); // Nodes for a constant interpolating polynomial.
// mesh.BasisFunctions( 0 ); // Definition of interpolating basis functions.
// =============================================================================================
// =============================================================================================
// --------------------------------------------------------------------------------------------
// Set mesh (Reading from file)
mesh.SetMaximalArea ( 0.001 ); // Set minimal area for mesh elements.
mesh.GenerateMesh(); // Generate triangulation for object mesh.
mesh.Nodes( 1 ); // Nodes for a constant interpolating polynomial.
mesh.BasisFunctions( 0 ); // Definition of interpolating basis functions.
mesh.Report( 1, 1, 0, "/home/jose/PhD_thesis/Seminars/Tutorial-RT/Example/Results/Mesh/" ); // Definition of interpolating basis functions.
// --------------------------------------------------------------------------------------------
// Set sampling lines
// --------------------------------------------------------------------------------------------
// Set elements
ifstream elem;
elem.open( "/usr/local/share/data/mesh/mesh3/elements.dat" ); // TODO: TURN TO A RELATIVE PATH.
if ( elem )
{
string lineelem, headelem;
int v1, v2, v3;
std::vector<int> ele1(3);
while ( getline(elem, lineelem) )
{
istringstream isse(lineelem);
if ( !(isse >> headelem >> v1 >> v2 >> v3) )
{
throw runtime_error("Input file 'elements.dat' has not the correct format: vertex0 \t vertex1 \t vertex2");
} else
{
ele1[0] = v1;
ele1[1] = v2;
ele1[2] = v3;
mesh.elements.push_back(ele1);
}
}
} else
{
throw runtime_error("Missing input file 'elements.dat'. It should be located in 'partons_root_directory/data/mesh/mesh1'");
}
elem.close();
// Set vertices
ifstream vert;
vert.open( "/usr/local/share/data/mesh/mesh3/vertices.dat" ); // TODO: CONVERT TO A RELATIVE PATH.
if ( vert )
{
string line, head;
double a, b;
NumA::point vertex;
while ( getline(vert, line) )
{
istringstream iss(line);
if ( !(iss >> head >> a >> b) )
{
throw runtime_error("Input file 'vertices.dat' has not the correct format: x-coord. \t y-coord.");
} else
{
vertex.x = a; vertex.y = b;
mesh.vertices.push_back(vertex);
}
}
} else
{
throw runtime_error("Missing input file 'vertices.dat'. It should be located in 'partons_root_directory/data/mesh/mesh3'");
}
vert.close();
// Set vertices' neighbors.
ifstream vneighbors;
vneighbors.open( "/usr/local/share/data/mesh/mesh3/neighborsv.dat" ); // TODO: CONVERT TO A RELATIVE PATH.
if ( vneighbors )
{
string linen;
int n1, n2, n3, n4, n5, n6, n7, n8, n9, n10;
std::vector<int> vn;
while ( getline(vneighbors, linen) )
{
vn.clear();
istringstream ss(linen);
ss >> n1;
vn.push_back(n1);
if ( ss >> n2 )
{
vn.push_back(n2);
if ( ss >> n3 )
{
vn.push_back(n3);
if ( ss >> n4 )
{
vn.push_back(n4);
if ( ss >> n5 )
{
vn.push_back(n5);
if ( ss >> n6 )
{
vn.push_back(n6);
if ( ss >> n7 )
{
vn.push_back(n7);
if ( ss >> n8 )
{
vn.push_back(n8);
if ( ss >> n9 )
{
vn.push_back(n9);
if ( ss >> n10 )
{
vn.push_back(n10);
}
}
}
}
}
}
}
}
}
mesh.vneighbors.push_back(vn);
}
} else
{
throw runtime_error("Missing input file 'neighborsv.dat'. It should be located in 'partons_root_directory/data/mesh/mesh3'");
}
vneighbors.close();
mesh.Nodes( 1 ); // Nodes for a constant interpolating polynomial.
mesh.BasisFunctions( 0 ); // Definition of interpolating basis functions.
// ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
// Sampling lines.
int m = 4*mesh.elements.size(); // No. sampling lines.
int m = 12*mesh.elements.size(); // No. sampling lines.
int nod = mesh.nodes.size(); // No. interpolation nodes.
x.resize(m); // Resize vectors containing sampling points accordingly.
......@@ -225,9 +75,25 @@ void RadonTransform::init()
x[i] = unif(re);
y[i] = unif(re);
xi[i] = x[i] * y[i]; // Randon sampling for \xi (DGLAP region).
// if ( xi[i] > 0.5 ) // Inconmplete RT - Uncomment this if-stament and comment the three above.
// {
// while ( xi[i] > 0.5 )
// {
// x[i] = unif(re);
// y[i] = unif(re);
// xi[i] = x[i] * y[i];
// }
// }
}
// Build Radon transform matrix - Set the corresponding attribute
// =============================================================================================
// =============================================================================================
// --------------------------------------------------------------------------------------------
// Build Radon transform matrix
// --------------------------------------------------------------------------------------------
RTMatrix = build_matrix( x, y, xi );
}
......@@ -236,7 +102,7 @@ void RadonTransform::init()
Eigen::MatrixXd RadonTransform::build_matrix( const std::vector<double> &x, const std::vector<double> &y, const std::vector<double> &xi )
{
const int m = x.size();
const int nod = mesh.nodes.size(); // Number of interpolation nodes // Number of sampling lines.
const int nod = mesh.nodes.size(); // Number of interpolation nodes // Number of sampling lines.
Eigen::MatrixXd MatrixPos(m, nod);
Eigen::MatrixXd MatrixNeg(m, nod);
......@@ -255,9 +121,9 @@ for ( int i = 0; i < mesh.vertices.size(); i++ )
mesh.vertices[i].y = - mesh.vertices[i].y;
}
Eigen::MatrixXd matrix = MatrixPos + MatrixNeg; // Compute Radon transform matrix and set the corresponding class attribute.
Eigen::MatrixXd matrix = MatrixPos + MatrixNeg; // Compute Radon transform matrix and set the corresponding class attribute.
for ( int i = 0; i < m; i++ ) // Handle P-scheme prescription.
for ( int i = 0; i < m; i++ ) // Handle P-scheme prescription.
{
matrix.row(i) *= (1-x[i]);
}
......@@ -400,12 +266,6 @@ Eigen::VectorXd RadonTransform::computeDD( const Eigen::VectorXd &GPD )
Eigen::VectorXd DD = cov*matrixT*GPD; // Solve system of equations as: Inverse(mTm)*mT*H
// ! Temporary step (print covariance matrix into an output file)
// ofstream mateval;
// mateval.open("results/Test_SaturatedModel/FullComputation/Cov-4noe.dat");
// mateval << cov << "\n";
// mateval.close();
return DD;
}
......@@ -458,13 +318,6 @@ std::vector<double> RadonTransform::computeGPD( const Eigen::VectorXd &DD, const
GPD.at(i) = GPD_mat(i);
}
// ! Temporary step (print covariance matrix into an output file)
// std::string pathtofile = "results/Test_SaturatedModel/FullComputation/MatEval-4noe.dat";
// ofstream mateval;
// mateval.open( pathtofile, std::ofstream::out | std::ofstream::app );
// mateval << matrix << "\n";
// mateval.close();
return GPD;
}
......@@ -487,16 +340,6 @@ double RadonTransform::computeGPD( const Eigen::VectorXd &DD, const double x, co
double GPD = GPD_mat(0);
// ! Temporary step (print covariance matrix into an output file)
// std::string pathtofile = "results/Test_SaturatedModel/FullComputation/";
// pathtofile.append(filename);
// pathtofile.append(".dat");
// ofstream mateval;
// mateval.open( pathtofile, std::ofstream::out | std::ofstream::app );
// mateval << matrix << "\n";
// mateval.close();
return GPD;
}
......@@ -525,16 +368,6 @@ std::vector<double> RadonTransform::computeGPD( const Eigen::VectorXd &DD, const
GPD.at(i) = GPD_mat(i);
}
// ! Temporary step (print covariance matrix into an output file)
// std::string pathtofile = "results/Test_SaturatedModel/FullComputation/";
// pathtofile.append(filename);
// pathtofile.append(".dat");
// ofstream mateval;
// mateval.open( pathtofile );
// mateval << matrix << "\n";
// mateval.close();
return GPD;
}
......
src/NumA/triangulation/mesh.cpp
View file @ bbd76592
......@@ -62,12 +62,12 @@ Mesh::~Mesh()
}
/**
* Mesh::SetMinimalArea( float m_area )
* Mesh::SetMaximalArea( float m_area )
*
* Setter for the minimal are of elements making up the mesh.
*
*/
void Mesh::SetMinimalArea( float area )
void Mesh::SetMaximalArea( float area )
{
if ( area > 0 && area < 1 )
{
......@@ -95,8 +95,6 @@ if ( area > 0 && area < 1 )
void Mesh::GenerateMesh()
{
cout << "GenerateMesh entered \n";
// ---------------------------------------------------------
// DEFINITION OF INPUT AND OPUTS FOR THE TRIANGLE FUNCTION
// ---------------------------------------------------------
......@@ -153,11 +151,9 @@ if ( area_constraint )
switches_string.append( "0.001"/*to_string( m_area )*/ ); // ! There is an issue with this implementation within PARTONS:
} else // ! to_string() converts float points into commas.
{
throw runtime_error("(mesh.cpp Line 172: m_area) Minimal area must be set. Otherwise no mesh would be generated.");
throw runtime_error("(mesh.cpp Line 151: m_area) Maximal area must be set. Otherwise no mesh would be generated.");
}
cout << "Switches to triangle: " << switches_string << "\n";
char* switches = const_cast<char*>(switches_string.c_str()); // triangulate() accepts 'char*' as
......
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