00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 #include "Epetra_config.h"
00019
00020 #ifdef HAVE_MPI
00021 #include "mpi.h"
00022 #include "Epetra_MpiComm.h"
00023 #else
00024 #include "Epetra_SerialComm.h"
00025 #endif
00026
00027 #include <Epetra_Operator.h>
00028 #include <Epetra_MultiVector.h>
00029 #include <Epetra_Vector.h>
00030 #include <Teuchos_ParameterList.hpp>
00031
00032 #include <algorithm>
00033 #include <iostream>
00034 #include <fstream>
00035 #include <sstream>
00036 #include <cstdlib>
00037 #include <cassert>
00038 #include <unistd.h>
00039
00040 #include "debugmeshbuilder.h"
00041 #include "tetmeshbuilder.h"
00042 #include "broadcastmeshbuilder.h"
00043 #include "netgenreader.h"
00044 #include "femaxmesh.h"
00045 #include "nedelecmesh.h"
00046 #include "lagrangemesh.h"
00047 #include "export_Epetra_CrsMatrix.h"
00048
00049 #include "eigsolv/AbstractEigsolvOperators.h"
00050 #include "eigsolv/LinearEigsolvOperators.h"
00051 #include "eigsolv/QuadraticEigsolvOperators.h"
00052 #include "eigsolv/BlockPCGSolver.h"
00053 #include "eigsolv/Projector.h"
00054 #include "eigsolv/JDBSYM.h"
00055 #include "eigsolv/CheckingTools.h"
00056
00057 #include <pbe.h>
00058
00059 using namespace std;
00060 using namespace mesh;
00061
00062 void usage(const Epetra_Comm& Comm, string err_msg) {
00063 if (Comm.MyPID() == 0) {
00064 cerr << err_msg << endl << endl;
00065 cerr << "Usage: matrix_generator_perlerd [-o order] mesh_file" << endl;
00066 }
00067 #ifdef HAVE_MPI
00068 MPI_Finalize();
00069 #endif
00070 exit(1);
00071 }
00072
00073
00074 int main(int argc, char *argv[])
00075 {
00076 #ifdef HAVE_MPI
00077 MPI_Init(&argc, &argv);
00078 Epetra_MpiComm Comm(MPI_COMM_WORLD);
00079 #else
00080 Epetra_SerialComm Comm;
00081 #endif
00082
00083 #ifdef __GNUC__
00084 std::set_terminate (__gnu_cxx::__verbose_terminate_handler);
00085 #endif
00086
00087 if (Comm.MyPID() == 0)
00088 cout << "running on " << Comm.NumProc() << " processor(s)..." << endl;
00089
00090
00091
00092
00093 const int nof_pbe_timers = 32;
00094 pbe_init(argv[0], nof_pbe_timers, 0, NULL);
00095 pbe_start(50, "whole_program");
00096
00097
00098
00099
00100
00101
00102 int element_order = 1;
00103
00104 extern char *optarg;
00105 extern int optind;
00106 int c;
00107 while ((c = getopt(argc, argv, "o:")) != -1) {
00108 switch (c) {
00109 case 'o':
00110 element_order = atoi(optarg);
00111 if (!(element_order == 1 || element_order == 2))
00112 usage(Comm, "Unsupported element order");
00113 break;
00114 case ':':
00115 case '?':
00116 usage(Comm, "Unsupported option");
00117 break;
00118 }
00119 }
00120
00121 if (argc - optind != 1) {
00122 usage(Comm, "");
00123 }
00124
00125 string file_name(argv[optind]);
00126
00127 if (access(file_name.c_str(), R_OK)) {
00128 usage(Comm, "Unable to access matrix file.");
00129 }
00130
00131 ifstream istr(file_name.c_str());
00132
00133
00134
00135
00136 TetMeshBuilder* builder = new TetMeshBuilder();
00137 #ifdef HAVE_MPI
00138 BroadcastMeshBuilder* broadcast = new BroadcastMeshBuilder(&Comm, builder);
00139 NetgenReader* reader(0);
00140 if (Comm.MyPID() == 0) {
00141 reader = new NetgenReader(&istr, broadcast);
00142 reader->read();
00143 } else
00144 broadcast->receiver();
00145 delete broadcast;
00146 if (Comm.MyPID() == 0)
00147 delete reader;
00148 #else
00149 NetgenReader* reader = new NetgenReader(&istr, builder);
00150 reader->read();
00151 delete reader;
00152 #endif
00153 TetMesh* tmesh = builder->get_mesh();
00154 delete builder;
00155
00156
00157 int sym_plane_config = 0x6;
00158 FemaxMesh mesh_mixed(*tmesh, element_order, sym_plane_config);
00159
00160
00161 AbstractEigsolvOperators* lop;
00162 Teuchos::ParameterList params;
00163 if (element_order == 1)
00164 lop = new LinearEigsolvOperators(mesh_mixed, params, Comm);
00165 else
00166 lop = new QuadraticEigsolvOperators(mesh_mixed, params, Comm);
00167
00168
00169
00170
00171
00172
00173 string::reverse_iterator it_dot = find(file_name.rbegin(), file_name.rend(), '.');
00174 string::reverse_iterator it_slash = find(file_name.rbegin(), file_name.rend(), '/');
00175 it_dot ++;
00176 string mesh_name(it_dot, it_slash);
00177 reverse(mesh_name.begin(), mesh_name.end());
00178
00179
00180 pbe_start(12, "matrix export");
00181 string m_file_name(mesh_name);
00182 m_file_name += "_M_";
00183 if (element_order == 1)
00184 m_file_name += "1";
00185 else
00186 m_file_name += "2";
00187 m_file_name += ".mtx";
00188
00189 export_Epetra_CrsMatrix(*(dynamic_cast<const Epetra_CrsMatrix*>(lop->getM())), 15, m_file_name);
00190 pbe_stop(12);
00191
00192
00193
00194
00195 pbe_stop(50);
00196 pbe_dump();
00197 pbe_finalize(1);
00198
00199
00200 delete lop;
00201 delete tmesh;
00202
00203 #ifdef HAVE_MPI
00204 MPI_Finalize();
00205 #endif
00206 return EXIT_SUCCESS;
00207 }