• Main Page
  • Related Pages
  • Namespaces
  • Classes
  • Files
  • File List
  • File Members

MinSurfaceMPS.cc

Go to the documentation of this file.
00001 //*****************************************************************//
00002 //    Albany 2.0:  Copyright 2012 Sandia Corporation               //
00003 //    This Software is released under the BSD license detailed     //
00004 //    in the file "license.txt" in the top-level Albany directory  //
00005 //*****************************************************************//
00006 
00007 #include "topology/Topology.h"
00008 #include <iostream>
00009 #include <fstream>
00010 #include <stdio.h>
00011 #include <string>
00012 #include <sstream>
00013 
00014 typedef stk::mesh::Entity Entity;
00015 
00016 #include <boost/numeric/ublas/matrix.hpp>
00017 #include <boost/numeric/ublas/io.hpp>
00018 #include <boost/numeric/ublas/matrix_proxy.hpp>
00019 #include <boost/numeric/ublas/vector.hpp>
00020 #include <boost/numeric/ublas/matrix_sparse.hpp>
00021 
00022 
00023 // \brief Function that helps creating text files with different names
00024 std::string itoa(int num);
00025 
00026 
00027 int main (int ac, char* av[]){
00028 
00029   std::string input_file;
00030   std::string normals_file;
00031   typedef stk::mesh::Entity Entity;
00032   std::string output_file;// = "output.e";
00033 
00034   using namespace std;
00035 
00036 
00037   if(ac == 2) {
00038     std::cout << "Generating MPS header file with mesh " << av[1] << std::endl;
00039     stringstream ss;
00040     ss << av[1]; //insert the char
00041     ss >> input_file; //extract -convert char into string
00042 
00043   } else if (ac == 3){
00044     std::cout <<  "Generating MPS rhs files with mesh " << av[1] << " and normals from file: " << av[2] << std::endl;
00045     stringstream ss1;
00046     stringstream ss2;
00047     ss1 << av[1]; //insert the char
00048     ss1 >> input_file; //extract -convert char into string
00049     ss2 << av[2];
00050     ss2 >> normals_file;
00051 
00052     //Check if third parameter from command line is a .txt file
00053     string s2 = ".txt";
00054     if (normals_file.find(s2) == std::string::npos) {
00055       std::cout << "Missing extension .txt for the third input" << '\n';
00056       throw(0);
00057     }
00058   } else {
00059     std::cout << "Usage: MinsurfaceMPS Mesh.exo or MinsurfaceMPS Mesh.exo Normals.txt" << std::endl; //USAGE OF THIS FILE FROM THE COMMAND LINE
00060     return 1;
00061   }
00062 
00063   //
00064   // Read the mesh
00065   //
00066   // Copied from Partition.cc
00067   Teuchos::GlobalMPISession mpiSession(&ac, &av);
00068   LCM::Topology topology(input_file, output_file);
00069 
00070 
00071   //--------------------------------------------------------------------------------------------------------
00072   //                        mps file OUTPUT for the MinsurfaceSolver.cpp
00073   //--------------------------------------------------------------------------------------------------------
00074   //NOTE: When analyzing a single mesh, the coefficientsObjFunction, and the BoundaryOperator
00075   //remain the same. This is set in the first part of the MPS file.The only one parameter that varies is the
00076   //BoundaryVector, which is defined in the second part of the MPS file
00077   //--------------------------------------------------------------------------------------------------------
00078 
00079   //--------------------------------------------------------------------------------------------------------
00080   //CREATE THE MPS FILE
00081   //--------------------------------------------------------------------------------------------------------
00082   stringstream ss;
00083   std::string file_name_;
00084   ss << av[1]; //insert the char
00085   ss >> file_name_; //extract -convert char into string
00086   string exoFile = ".exo"; //Define the string with the extension file of input mesh
00087   std::size_t Position = file_name_.find(exoFile); //Find the position where .exo starts
00088 
00089   //First part of MPS file
00090   //Create the name of the text file that contains the header of the MPS file
00091   std::string file_name1 = "MPS_"+file_name_.substr(0, Position)+"_header.txt";
00092   std::string file_name2 = file_name_.substr(0, Position); //name of the problem (important for the LP solver that reads the mps file)
00093 
00094   //Compute the header of the mps file if 2 command line inputs are given
00095   if (ac == 2)
00096   {
00097     ofstream mpsFile(file_name1.c_str(),std::ofstream::out);
00098     if (mpsFile.is_open())
00099     {
00100 
00101       mpsFile << "NAME" <<" "<< "MPS_"<<file_name2<<"_" <<endl; //The linear problem name comes from the second argument in the command line
00102       mpsFile << "ROWS" <<endl;
00103       mpsFile <<" "<<"N"<<" "<<"Z"<<endl; //Define the variable that holds the objective function. "N" means that it is unbounded
00104 
00105       //Define the rows and their format
00106       std::vector<std::vector<int> > edgesDirec =  topology.edgesDirections();
00107       for (unsigned int i = 0; i< edgesDirec.size();i++)
00108       {
00109         mpsFile<<" "<<"E"<<" "<<"R"<<i+1<<endl;
00110       }
00111       //Define the number of columns
00112       mpsFile << "COLUMNS" <<endl;
00113       //Obtain the boundary operator and the areas of all faces
00114       std::vector<std::vector<double> > outputForMpsFile = topology.outputForMpsFile();
00115       for (unsigned int i = 0; i < outputForMpsFile.size(); i+=4)
00116       {
00117         //The columns of the boundary operator need to be saved in order X1, X2, X3,...etc. Also their corresponding areas
00118         mpsFile <<" "<<"X"<<outputForMpsFile[i+1][0]<<" "<<"Z"<<" "<<setprecision(numeric_limits<long double>::digits10)<<outputForMpsFile[i][0]<<endl;
00119         mpsFile <<" "<<"X"<<outputForMpsFile[i+1][0]<<" "<<"R"<<outputForMpsFile[i+1][1]<<" "<<outputForMpsFile[i+1][2]<<endl;
00120         mpsFile <<" "<<"X"<<outputForMpsFile[i+2][0]<<" "<<"R"<<outputForMpsFile[i+2][1]<<" "<<outputForMpsFile[i+2][2]<<endl;
00121         mpsFile <<" "<<"X"<<outputForMpsFile[i+3][0]<<" "<<"R"<<outputForMpsFile[i+3][1]<<" "<<outputForMpsFile[i+3][2]<<endl;
00122         mpsFile <<" "<<"X"<<outputForMpsFile[i+1][3]<<" "<<"Z"<<" "<<setprecision(numeric_limits<long double>::digits10)<<outputForMpsFile[i][0]<<endl;
00123         mpsFile <<" "<<"X"<<outputForMpsFile[i+1][3]<<" "<<"R"<<outputForMpsFile[i+1][4]<<" "<<outputForMpsFile[i+1][5]<<endl;
00124         mpsFile <<" "<<"X"<<outputForMpsFile[i+2][3]<<" "<<"R"<<outputForMpsFile[i+2][4]<<" "<<outputForMpsFile[i+2][5]<<endl;
00125         mpsFile <<" "<<"X"<<outputForMpsFile[i+3][3]<<" "<<"R"<<outputForMpsFile[i+3][4]<<" "<<outputForMpsFile[i+3][5]<<endl;
00126       }
00127     }
00128   }
00129 
00130   //Second part of the MPS file: Define the r vector into the mps file
00131   //Compute the normals below when the number of inputs is 3
00132   else
00133   {
00134     ifstream normals_(normals_file.c_str()); //Read the file that contains all normals
00135     int numberNormals_; //Variable to define the number of rows of the matrix mNormals
00136     if (normals_.is_open())
00137     {
00138       normals_ >> numberNormals_; //Obtain the total number of input normals
00139     }
00140 
00141     //mNormal is a matrix containing all the normal vectors
00142     std::vector<std::vector<double> > mNormals(numberNormals_,std::vector<double>(3,0));//numberNormals_
00143 
00144     if (normals_.is_open())
00145     {
00146       for (int i = 0; i < numberNormals_ ; i++)
00147       {
00148         normals_ >> mNormals[i][0] >> mNormals[i][1]>> mNormals[i][2];
00149       }
00150     }
00151 
00152     //Extract the numbers from the matrix that contains the normals
00153     std::vector<std::vector<int> >BoundaryVector;
00154     for(int i = 0; i < numberNormals_; i++)
00155     {
00156       std::string file_name = "MPS_"+file_name2+"_rhs_" + itoa(i+1) + ".txt"; //file_name2 comes from the second input in the command line
00157 
00158       ofstream mpsFile(file_name.c_str(),std::ofstream::out);
00159       if (mpsFile.is_open())
00160       {
00161         mpsFile <<"RHS"<<endl;
00162 
00163         //Obtain the vector that defines the 1D boundary
00164         std::vector<std::vector<int> >BoundaryVector;
00165 
00166         std::vector<double> normalToPlane;
00167         normalToPlane.push_back(mNormals.at(i).at(0));
00168         normalToPlane.push_back(mNormals.at(i).at(1));
00169         normalToPlane.push_back(mNormals.at(i).at(2));
00170 
00171         //Finds coordinates of 3 vectors normal to the normal vector.
00172         std::vector<std::vector<double> > pointsOnPlane = topology.getCoordinatesOfTriangle(normalToPlane);
00173 
00174 
00175         //Finds the Node (Entity of rank 0) that is closest to each of the points on the plane
00176         std::vector<Entity*> closestNodes = topology.getClosestNodesOnSurface(pointsOnPlane);
00177 
00178         //Finds the identifiers of the nodes (entity rank 0) along the shortest
00179         //path connecting the three points
00180         std::vector<std::vector<int> > ShortestPathFinal = topology.shortestpath(closestNodes);
00181         BoundaryVector = topology.boundaryVector(ShortestPathFinal);//.boundaryVectorOuterSurface
00182         for (unsigned int i =0; i < BoundaryVector.size();i++)
00183         {
00184           mpsFile <<" "<<"RHS1"<<" "<<"R"<<BoundaryVector[i][0]<<" "<<BoundaryVector[i][1]<<endl;
00185         }
00186 
00187         //Define the bounds of the variables X1, X2, ...etc
00188         mpsFile << "BOUNDS"<<endl;
00189         mpsFile << "ENDATA"<<endl;
00190       }//End of second part of mps file
00191     }//End of for loop
00192   }
00193   return 0;
00194 }
00195 
00196 // \brief Function that helps creating text files with different names
00197 std::string itoa(int num) {
00198   std::stringstream converter;
00199   converter << num;
00200   return converter.str();
00201 }
00202 

Generated on Wed Mar 26 2014 18:36:40 for Albany: a Trilinos-based PDE code by  doxygen 1.7.1