examples/matrix_generator_perlerd.cpp

Go to the documentation of this file.
00001 /***************************************************************************
00002                           main.cpp  -  description
00003                              -------------------
00004     begin                : Fri Dec 12 13:07:39 CET 2003
00005     copyright            : (C) 2003 by Roman Geus
00006     email                : roman.geus@psi.ch
00007 ***************************************************************************/
00008 
00009 /***************************************************************************
00010  *                                                                         *
00011  *   This program is free software; you can redistribute it and/or modify  *
00012  *   it under the terms of the GNU General Public License as published by  *
00013  *   the Free Software Foundation; either version 2 of the License, or     *
00014  *   (at your option) any later version.                                   *
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     // Init PBE
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     // parse command-line options
00099     //
00100 
00101     // default
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     // Generate mesh data from file input on all processors
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     // create FE mesh data structures
00157     int sym_plane_config = 0x6;
00158     FemaxMesh mesh_mixed(*tmesh, element_order, sym_plane_config);
00159 
00160     // create operators
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     // matrix output
00170     //
00171 
00172     // get mesh name
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     // M
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     // Finalize PBE
00194     //
00195     pbe_stop(50);
00196     pbe_dump();
00197     pbe_finalize(1);
00198 
00199     // free memory
00200     delete lop;
00201     delete tmesh;
00202 
00203 #ifdef HAVE_MPI
00204     MPI_Finalize();
00205 #endif
00206     return EXIT_SUCCESS;
00207 }

Generated on Fri Oct 26 13:35:11 2007 for FEMAXX (Finite Element Maxwell Eigensolver) by  doxygen 1.4.7