examples/matrix_generator_perlerd.cpp

Go to the documentation of this file.

/***************************************************************************
                          main.cpp  -  description
                             -------------------
    begin                : Fri Dec 12 13:07:39 CET 2003
    copyright            : (C) 2003 by Roman Geus
    email                : roman.geus@psi.ch
***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/

#include "Epetra_config.h"

#ifdef HAVE_MPI
#include "mpi.h"
#include "Epetra_MpiComm.h"
#else
#include "Epetra_SerialComm.h"
#endif

#include <Epetra_Operator.h>
#include <Epetra_MultiVector.h>
#include <Epetra_Vector.h>
#include <Teuchos_ParameterList.hpp>

#include <algorithm>
#include <iostream>
#include <fstream>
#include <sstream>
#include <cstdlib>
#include <cassert>
#include <unistd.h>

#include "debugmeshbuilder.h"
#include "tetmeshbuilder.h"
#include "broadcastmeshbuilder.h"
#include "netgenreader.h"
#include "femaxmesh.h"
#include "nedelecmesh.h"
#include "lagrangemesh.h"
#include "export_Epetra_CrsMatrix.h"

#include "eigsolv/AbstractEigsolvOperators.h"
#include "eigsolv/LinearEigsolvOperators.h"
#include "eigsolv/QuadraticEigsolvOperators.h"
#include "eigsolv/BlockPCGSolver.h"
#include "eigsolv/Projector.h"
#include "eigsolv/JDBSYM.h"
#include "eigsolv/CheckingTools.h"

#include <pbe.h>

using namespace std;
using namespace mesh;

void usage(const Epetra_Comm& Comm, string err_msg) {
    if (Comm.MyPID() == 0) {
        cerr << err_msg << endl << endl;
        cerr << "Usage: matrix_generator_perlerd [-o order] mesh_file" << endl;
    }
#ifdef HAVE_MPI
    MPI_Finalize();
#endif
    exit(1);
}


int main(int argc, char *argv[])
{
#ifdef HAVE_MPI
    MPI_Init(&argc, &argv);
    Epetra_MpiComm Comm(MPI_COMM_WORLD);
#else
    Epetra_SerialComm Comm;
#endif

#ifdef __GNUC__
    std::set_terminate (__gnu_cxx::__verbose_terminate_handler);
#endif

    if (Comm.MyPID() == 0)
        cout << "running on " << Comm.NumProc() << " processor(s)..." << endl;

    //
    // Init PBE
    //
    const int nof_pbe_timers = 32;
    pbe_init(argv[0], nof_pbe_timers, 0, NULL);
    pbe_start(50, "whole_program");

    //
    // parse command-line options
    //

    // default
    int element_order = 1;

    extern char *optarg;
    extern int optind;
    int c;
    while ((c = getopt(argc, argv, "o:")) != -1) {
        switch (c) {
        case 'o':
            element_order = atoi(optarg);
            if (!(element_order == 1 || element_order == 2))
                usage(Comm, "Unsupported element order");
            break;
        case ':':
        case '?':
            usage(Comm, "Unsupported option");
            break;
        }
    }

    if (argc - optind != 1) {
        usage(Comm, "");
    }

    string file_name(argv[optind]);

    if (access(file_name.c_str(), R_OK)) {
        usage(Comm, "Unable to access matrix file.");
    }

    ifstream istr(file_name.c_str());

    //
    // Generate mesh data from file input on all processors
    //
    TetMeshBuilder* builder = new TetMeshBuilder();
#ifdef HAVE_MPI
    BroadcastMeshBuilder* broadcast = new BroadcastMeshBuilder(&Comm, builder);
    NetgenReader* reader(0);
    if (Comm.MyPID() == 0) {
        reader = new NetgenReader(&istr, broadcast);
        reader->read();
    } else
        broadcast->receiver();
    delete broadcast;
    if (Comm.MyPID() == 0)
        delete reader;
#else
    NetgenReader* reader = new NetgenReader(&istr, builder);
    reader->read();
    delete reader;
#endif
    TetMesh* tmesh = builder->get_mesh();
    delete builder;

    // create FE mesh data structures
    int sym_plane_config = 0x6;
    FemaxMesh mesh_mixed(*tmesh, element_order, sym_plane_config);

    // create operators
    AbstractEigsolvOperators* lop;
    Teuchos::ParameterList params;
    if (element_order == 1)
        lop = new LinearEigsolvOperators(mesh_mixed, params, Comm);
    else
        lop = new QuadraticEigsolvOperators(mesh_mixed, params, Comm);

    //
    // matrix output
    //

    // get mesh name
    string::reverse_iterator it_dot = find(file_name.rbegin(), file_name.rend(), '.');
    string::reverse_iterator it_slash = find(file_name.rbegin(), file_name.rend(), '/');
    it_dot ++;
    string mesh_name(it_dot, it_slash);
    reverse(mesh_name.begin(), mesh_name.end());

    // M
    pbe_start(12, "matrix export");
    string m_file_name(mesh_name);
    m_file_name += "_M_";
    if (element_order == 1)
        m_file_name += "1";
    else
        m_file_name += "2";
    m_file_name += ".mtx";

    export_Epetra_CrsMatrix(*(dynamic_cast<const Epetra_CrsMatrix*>(lop->getM())), 15, m_file_name);
    pbe_stop(12);

    //
    // Finalize PBE
    //
    pbe_stop(50);
    pbe_dump();
    pbe_finalize(1);

    // free memory
    delete lop;
    delete tmesh;

#ifdef HAVE_MPI
    MPI_Finalize();
#endif
    return EXIT_SUCCESS;
}

---