src/fem/lagrangemesh.cpp

Go to the documentation of this file.

/***************************************************************************
                          lagrangemesh.cpp  -  description
                             -------------------
    begin                : Fri Feb 13 2004
    copyright            : (C) 2004 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 <cassert>
#include "lagrangemesh.h"
#include "tetmesh/point.h"
#include "tetmesh/edge.h"

LagrangeMesh::LagrangeMesh(mesh::TetMesh* mesh, int order, int sym_plane_config)
    : _mesh(mesh), _order(order) {
    
    _pmesh = dynamic_cast<mesh::ParallelTetMesh*>(_mesh);
    
    // Order should be equal to 1 or 2.    
    assert(1 <= order && order <= 2);
    
    _nof_orig_dofs.resize(order + 1);
    _nof_orig_dofs[1] = _mesh->get_nof_points();
    switch (order) {
    case 1:
        _nof_orig_dofs[0] = _nof_orig_dofs[1];
        break;
    case 2:
        _nof_orig_dofs[2] = _mesh->get_nof_edges();
        _nof_orig_dofs[0] = _nof_orig_dofs[1] + _nof_orig_dofs[2];
        break;
    }
    
    //
    // Global dofs
    //
    _num_global_dofs.resize(order + 1);
    _num_global_dofs[1] = _pmesh->get_num_global_points();
    switch (order) {
    case 1:
        _num_global_dofs[0] = _num_global_dofs[1];
        break;
    case 2: 
        _num_global_dofs[2] = _pmesh->get_num_global_edges();
        _num_global_dofs[0] = _num_global_dofs[1] + _num_global_dofs[2];
    }
    
    set_map(sym_plane_config);
}
  
LagrangeMesh::~LagrangeMesh()
{ }

mesh::id_t LagrangeMesh::get_num_my_dofs(int order) const {
    return _nof_orig_dofs.at(order);   
}

mesh::id_t LagrangeMesh::get_num_global_mdofs(int order) const {
    return _nof_mapped_dofs.at(order);
}

mesh::id_t LagrangeMesh::get_num_global_dofs(int order) const {
    return _num_global_dofs.at(order);
}

void LagrangeMesh::set_map(unsigned int sym_plane_config) {
    rDebugAll("num_my_dofs in lagrange mesh = %d", get_num_my_dofs(0));

    // Create global elimination array
    std::vector<bool> elim(get_num_my_dofs(0));
    for (mesh::id_t i = 0; i < _mesh->get_nof_points(); i ++) {
        mesh::Point* point = _mesh->get_point(i);
        if (point->on_boundary()) {
            elim[i] = true;
        } else if ((sym_plane_config ^ 0x7) & point->on_symmetry_planes()) {
            elim[i] = true;
        } else {
            elim[i] = false;
        }
    }

    if (_order > 1) {
        for (mesh::id_t i = 0; i < _mesh->get_nof_edges(); i ++) {
            mesh::Edge* edge = _mesh->get_edge(i);
            if (edge->on_boundary()) {
                elim[get_num_my_dofs(1) + i] = true;
            } else if ((sym_plane_config ^ 0x7) & edge->on_symmetry_planes()) {
                elim[get_num_my_dofs(1) + i] = true;
            } else {
                elim[get_num_my_dofs(1) + i] = false;
            }
        }
    } 

    _nof_mapped_dofs.resize(_order + 1);

    //
    // 1st order gdofs only
    //
    std::vector<mesh::id_t> eliminated_dofs;
    for (mesh::id_t ldof = 0; ldof < get_num_my_dofs(1); ++ ldof) {
        if(elim[ldof])
            eliminated_dofs.push_back(gdof(ldof));
    }
    parallel::compute_mgid_map(get_num_global_dofs(1), _pmesh->get_node_gids(), eliminated_dofs, 
                               &_nof_mapped_dofs[1], _mgid_map, _pmesh->get_comm());
                               
    if (_order == 1) {
        _nof_mapped_dofs[0] = _nof_mapped_dofs[1];
    } else {
        //
        // Recompute mgids using all gdofs
        //
        rAssert(_order == 2);
        for (mesh::id_t ldof = get_num_my_dofs(1); ldof < get_num_my_dofs(0); ++ ldof) {
            if(elim[ldof])
                eliminated_dofs.push_back(gdof(ldof));
        }
        std::vector<mesh::id_t> all_my_gdofs;
        for (mesh::id_t ldof = 0; ldof < get_num_my_dofs(0); ++ ldof)
            all_my_gdofs.push_back(gdof(ldof));
        parallel::compute_mgid_map(get_num_global_dofs(0), all_my_gdofs, eliminated_dofs, 
                                   &_nof_mapped_dofs[0], _mgid_map, _pmesh->get_comm());
        _nof_mapped_dofs[2] = _nof_mapped_dofs[0] - _nof_mapped_dofs[1];
    }
    
    rDebugAll("lagrange: num_my_dofs=%d, num_my_eliminated_dofs=%d", 
              get_num_my_dofs(0) , static_cast<int>(eliminated_dofs.size()));
    rDebugAll("lagrange mesh: #gdofs=%d, #mgdofs=%d", get_num_global_dofs(0), _nof_mapped_dofs[0]);
    
    // debug
#if 0
    for (std::map<id_t, id_t>::iterator it=_mgid_map.begin(); it != _mgid_map.end(); ++ it) {
        rDebugAll("lagrange mesh _mgid_map[%d] = %d", (*it).first, (*it).second);
    }
    rDebugAll("lagrange mesh num_global_mgids = %d", num_global_mgids);
    for (int i = 0; i < num_my_dofs; ++ i) {
        rDebugAll("lagrange mesh lid:%d, gid:%d, mgid:%d", i, _pmesh->point_gid(i), map_dof(i));
    }
#endif
}

mesh::id_t LagrangeMesh::gdof(mesh::id_t ldof) const {
    if (ldof < _mesh->get_nof_points())
        return _pmesh->point_gid(ldof);
    else
        return _pmesh->get_num_global_points() \
            + _pmesh->edge_gid(ldof - _mesh->get_nof_points());
}

bool LagrangeMesh::is_owned_gdof(mesh::id_t gdof) const {
    rAssert(0 <= gdof && gdof < get_num_global_dofs(0));
    
    if (gdof < get_num_global_dofs(1)) {
        return _pmesh->is_owned_point(gdof);    
    } else {
        return _pmesh->is_owned_edge(gdof - get_num_global_dofs(1));
    }
}

mesh::id_t LagrangeMesh::map_dof(mesh::id_t ldof) const {
    std::map<mesh::id_t, mesh::id_t>::const_iterator found = _mgid_map.find(gdof(ldof));
    if (found == _mgid_map.end())
        return mesh::ID_NONE;
    else
        return (*found).second;
}

---