src/fem/femaxmesh.cpp

Go to the documentation of this file.

/***************************************************************************
                          femaxmesh.cpp  -  description
                             -------------------
    begin                : Fri Dec 19 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_Map.h"
#include "femaxmesh.h"


FemaxMesh::FemaxMesh(mesh::TetMesh& mesh, int order, int sym_plane_config)
    : mesh_(mesh),
      order_(order),
      mesh_nedelec_(new NedelecMesh(&mesh, order, sym_plane_config)),
      mesh_lagrange_(new LagrangeMesh(&mesh, order, sym_plane_config))
{ 
    _pmesh = dynamic_cast<mesh::ParallelTetMesh*>(&mesh_);
}

FemaxMesh::~FemaxMesh() {
    delete mesh_nedelec_;
    delete mesh_lagrange_;
}

void FemaxMesh::constructY(Epetra_FECrsMatrix& Y, Epetra_Map& domain_map, Epetra_Map& range_map) {
    double one = 1.0;
    double minus_one = -1.0;

    // row and column offset for (2,2)-block
    int i_offset22 = mesh_nedelec_->get_num_global_mdofs(1);
    int j_offset22 = mesh_lagrange_->get_num_global_mdofs(1);
    
    mesh::TetMesh::edge_iterator iter = mesh_.edge_begin();
    mesh::TetMesh::edge_iterator iter_end = mesh_.edge_end();

    // loop over all edges
    for (; iter != iter_end; ++ iter) {
        int edge_id_mapped = mesh_nedelec_->map_dof((*iter).get_id());
        int node1_id_mapped = mesh_lagrange_->map_dof((*iter).get_start_id());
        int node2_id_mapped = mesh_lagrange_->map_dof((*iter).get_end_id());

        if (edge_id_mapped != -1) {
            // start node
            if (node1_id_mapped != -1) {
                Y.InsertGlobalValues(1, &edge_id_mapped, 1, &node1_id_mapped, &minus_one);
            }
            // end node            
            if (node2_id_mapped != -1) {
                Y.InsertGlobalValues(1, &edge_id_mapped, 1, &node2_id_mapped, &one);
            }
            
            // (2,2)-block: identity matrix
            if (order_ == 2) {
                int row = i_offset22 + edge_id_mapped;
                int col = j_offset22 + edge_id_mapped;
                Y.InsertGlobalValues(1, &row, 1, &col, &one);
            }
        }
    }

    Y.GlobalAssemble(false);
    Y.FillComplete(domain_map, range_map);
    
    // Patch entries, which have been inserted multiple times
    for (int i_row = 0; i_row < Y.NumMyRows(); i_row++) {
        int num_entries;
        double* values;
        Y.ExtractMyRowView(i_row, num_entries, values);        
        for (int j_entry = 0; j_entry < num_entries; j_entry++) {
            // assert (values[j_entry] >= 1 || values[j_entry] <= -1);
            if (!(values[j_entry] >= 1 || values[j_entry] <= -1))
                rWarningAll("Y[%d,%d] = %g", i_row, j_entry, values[j_entry]);
            if (values[j_entry] > 0.0)
                values[j_entry] = 1.0;
            else if (values[j_entry] < 0.0)
                values[j_entry] = -1.0;
        }
    }
    // cerr <<endl<<endl<< "!!!!! Y Matrix !!!!!!"<<endl<<Y;
}

void FemaxMesh::constructY_transp(Epetra_FECrsMatrix& Y_transp, Epetra_Map& domain_map, Epetra_Map& range_map) {
    double one = 1.0;
    double minus_one = -1.0;

    // row and column offset for (2,2)-block
    int i_offset22 = mesh_lagrange_->get_num_global_mdofs(1);
    int j_offset22 = mesh_nedelec_->get_num_global_mdofs(1);

    mesh::TetMesh::edge_iterator iter = mesh_.edge_begin();
    mesh::TetMesh::edge_iterator iter_end = mesh_.edge_end();

    // loop over all edges
    for (; iter != iter_end; ++ iter) {
        int edge_id_mapped = mesh_nedelec_->map_dof((*iter).get_id());
        int node1_id_mapped = mesh_lagrange_->map_dof((*iter).get_start_id());
        int node2_id_mapped = mesh_lagrange_->map_dof((*iter).get_end_id());
        
        if (edge_id_mapped != -1) {
            // start node
            if (node1_id_mapped != -1)
                Y_transp.InsertGlobalValues(1, &node1_id_mapped, 1, &edge_id_mapped, &minus_one);
            // end node
            if (node2_id_mapped != -1)
                Y_transp.InsertGlobalValues(1, &node2_id_mapped, 1, &edge_id_mapped, &one);
                
            // (2,2)-block: identity matrix
            if (order_ == 2) {
                int row = i_offset22 + edge_id_mapped;
                int col = j_offset22 + edge_id_mapped;
                Y_transp.InsertGlobalValues(1, &row, 1, &col, &one);
            }
        }
    }

    Y_transp.GlobalAssemble(false);
    Y_transp.FillComplete(domain_map, range_map);
    
    // Patch entries, which have been inserted multiple times
    for (int i_row = 0; i_row < Y_transp.NumMyRows(); i_row++) {
        int num_entries;
        double* values;
        Y_transp.ExtractMyRowView(i_row, num_entries, values);        
        for (int j_entry = 0; j_entry < num_entries; j_entry++) {
            // assert (values[j_entry] >= 1 || values[j_entry] <= -1);
            if (!(values[j_entry] >= 1 || values[j_entry] <= -1))
                rWarningAll("Y'[%d,%d] = %g", i_row, j_entry, values[j_entry]);
            if (values[j_entry] > 0.0)
                values[j_entry] = 1.0;
            else if (values[j_entry] < 0.0)
                values[j_entry] = -1.0;
        }
    }
    // cerr <<endl<<endl<< "!!!!! Y_transp Matrix !!!!!!"<<endl<<Y_transp;
}

---