src/fem/nedelecmesh.cpp

Go to the documentation of this file.

/***************************************************************************
                          nedelecmesh.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 <cassert>
#include <cmath>
#include <new>
#include <vector>
#include <iostream>
#include <iomanip>
#include <algorithm>

#include "Epetra_Map.h"
#include "Epetra_Vector.h"

#include "rlog/rlog.h"

#include "nedelecelement1.h"
#include "nedelecelement2.h"
#include "nedelecmesh.h"

#include "matrix.h"

#include "fmxxphysicomath.h"
#include "fmxxtopology.h"

using namespace std;
using namespace colarray;
using namespace mesh;

typedef std::pair<mesh::id_t, mesh::id_t> nodepair;
typedef utility::triple<nodepair, mesh::id_t, mesh::id_t> dof_with_ids;

static bool less_node_lex(const dof_with_ids& a, const dof_with_ids& b) {
    if (a.first.first != b.first.first)                
        return a.first.first < b.first.first;                                    
    else
        return a.first.second < b.first.second;    
}

static bool less_id(const dof_with_ids& a, const dof_with_ids& b) {
    return a.second < b.second;
}     

NedelecMesh::NedelecMesh(mesh::TetMesh* mesh, int order, unsigned 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);
    
    // Locally stored dofs
    _nof_orig_dofs.resize(order + 1);
    _nof_orig_dofs[1] = _mesh->get_nof_edges();
    switch (order) {
    case 1:
        _nof_orig_dofs[0] = _nof_orig_dofs[1];
        break;
    case 2:
        _nof_orig_dofs[2] = _mesh->get_nof_edges() + 3*_mesh->get_nof_faces();
        _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_edges();
    switch (order) {
    case 1:
        _num_global_dofs[0] = _num_global_dofs[1];
        break;
    case 2: 
        _num_global_dofs[2] = _pmesh->get_num_global_edges() + 3*_pmesh->get_num_global_faces();
        _num_global_dofs[0] = _num_global_dofs[1] + _num_global_dofs[2];
    }    

    if (order == 1)
        _element = new NedelecElement1();
    else // order == 2
        _element = new NedelecElement2();
  
    this->set_map(sym_plane_config);
}

NedelecMesh::~NedelecMesh() {
    delete _element;
}

void NedelecMesh::assembleAM(Epetra_FECrsMatrix* A,
                             Epetra_FECrsMatrix* M) {
    int info;
    Matrix<double> Ae, Me;

    // loop over locally stored elements
    TetMesh::tet_iterator it;
    for (it = _mesh->tet_begin(); it != _mesh->tet_end(); ++ it) {
        Tet& tet = *it;
        NedelecElement* elem = get_element(tet.get_id());
        int nof_local_dofs = elem->get_nof_dofs();
        const int *dof_ids = elem->get_dof_ids(_mesh, &tet);
//        rDebugAll("dof_lids: %u %u %u %u %u %u",
//                  dof_ids[0], dof_ids[1], dof_ids[2], dof_ids[3], dof_ids[4], dof_ids[5]);
//        rDebugAll("dof_mgids: %d %d %d %d %d %d",
//                  map_dof(dof_ids[0]), map_dof(dof_ids[1]), map_dof(dof_ids[2]), 
//                  map_dof(dof_ids[3]), map_dof(dof_ids[4]), map_dof(dof_ids[5]));

        elem->get_element_matrix(&tet, Ae, Me);
       
        for (int i_element = 0; i_element < nof_local_dofs; i_element ++) {
            int i_mapped = this->map_dof(dof_ids[i_element]);
            if (i_mapped != -1) {
    
                for (int j_element = 0; j_element < nof_local_dofs; j_element ++) {
                    int j_mapped = this->map_dof(dof_ids[j_element]);
                    if (j_mapped != -1) {
                        double val;
            
                        val = Ae(i_element, j_element);     
                        if (val != 0.0) {
                            info = A->SumIntoGlobalValues(1, &i_mapped, 1, &j_mapped, &val);
                            assert(info >= 0);
                            if (info > 0) {
                                info = A->InsertGlobalValues(1, &i_mapped, 1, &j_mapped, &val);
                                assert(info >= 0);
                            }                            
                        }

                        val = Me(i_element, j_element);
                        if (val != 0.0) {            
                            info = M->SumIntoGlobalValues(1, &i_mapped, 1, &j_mapped, &val);
                            assert(info >= 0);                            
                            if (info > 0) {
                                info = M->InsertGlobalValues(1, &i_mapped, 1, &j_mapped, &val);
                                assert(info >= 0);
                            }                            

                        }
                    }
                }
            }
        }    
    }
    A->GlobalAssemble(false);
    M->GlobalAssemble(false);
}

void NedelecMesh::set_map(unsigned int sym_plane_config) {
    // Create elimination array
    vector<bool> elim(get_num_my_dofs(0));
    for (mesh::id_t e = 0; e < _mesh->get_nof_tets(); e ++) {
        NedelecElement* element = get_element(e);
        mesh::Tet* tet = _mesh->get_tet(e);
        const int* dof_ids = element->get_dof_ids(_mesh, tet);
        const bool* elim_dofs = element->get_eliminated_dofs(_mesh, tet, sym_plane_config);
        for (int k = 0; k < element->get_nof_dofs(); k ++)
            elim[dof_ids[k]] = elim_dofs[k];
    }

    if (_order > 1) {
        // Eliminate additional dofs due to linear dependece of face dofs
        for (mesh::id_t i = 0; i < _mesh->get_nof_faces(); i ++)
            elim[2*_mesh->get_nof_edges() + 3*i] = true;
    }

    _nof_mapped_dofs.resize(_order + 1);
    
    //
    // 1st order gdofs only
    //
    vector<mesh::id_t> eliminated_dofs;
    for (mesh::id_t i = 0; i < get_num_my_dofs(1); ++ i) {
        if(elim[i])
            eliminated_dofs.push_back(_pmesh->edge_gid(i));
    }
    parallel::compute_mgid_map(get_num_global_dofs(1), _pmesh->get_edge_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));
        }
        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("nedelec: num_my_dofs=%d, num_my_eliminated_dofs=%d", 
              get_num_my_dofs(0) , static_cast<int>(eliminated_dofs.size()));
    rDebugAll("nedelec mesh: #gdofs=%d, #mgdofs=%d", get_num_global_dofs(0), _nof_mapped_dofs[0]);
    
    // debug
#if 0
    for (std::map<mesh::id_t, mesh::id_t>::iterator it=_mgid_map.begin();
         it != _mgid_map.end();
         ++ it)
    {
        rDebugAll("_mgid_map[%d] = %d", (*it).first, (*it).second);
    }
    {
        MPI_Comm comm = MPI_COMM_WORLD;
        int rank, size;
        MPI_Comm_size(comm, &size);
        MPI_Comm_rank(comm, &rank);
        for (int p = 0; p < size; ++ p) {
            if (p == rank) {
                for (int i = 0; i < num_my_dofs; ++ i) {
                    rDebugAll("--lid:%d, gid:%d, mgid:%d", i, _pmesh->edge_gid(i), map_dof(i));
                }
                fsync(1);
            }
            MPI_Barrier(comm);
        }
    }
    rDebugAll("edge mgid to point gid mapping...");
    for (mesh::id_t mgid = 0; mgid < num_global_mgids; ++ mgid) {
        // Find gid
        std::map<mesh::id_t, mesh::id_t>::iterator found;
        found = std::lower_bound(_mgid_map.begin(), _mgid_map.end(), mgid, less_than_value);
        if (found != _mgid_map.end() && (*found).second == mgid) {
            mesh::id_t edge_gid = (*found).first;
            // Find edge lid
            mesh::id_t edge_lid;
            for (edge_lid = 0; _pmesh->edge_gid(edge_lid) != edge_gid; ++ edge_lid);
            assert(_pmesh->edge_gid(edge_lid) == edge_gid);
          
            rDebugAll("edge=mgid=%u, edge_gid=%u: point_gid0=%u, point_gid1=%u", 
                      mgid, 
                      edge_gid, 
                      _pmesh->point_gid(_pmesh->get_edge(edge_lid)->get_start_id()),
                      _pmesh->point_gid(_pmesh->get_edge(edge_lid)->get_end_id()));
        }
    }
#endif
}

void NedelecMesh::generate_sorted_AM(const Epetra_FECrsMatrix* A, const Epetra_FECrsMatrix* M, Epetra_FECrsMatrix* A_s, Epetra_FECrsMatrix* M_s) {
    // Generate vector of point pairs
    // Equals owned edges not removed in 1st order case.
    vector<dof_with_ids> my_dofs;        
    
    mesh::id_t my_nof_edges = _pmesh->get_nof_edges();
    for (mesh::id_t edge_lid = 0; edge_lid < my_nof_edges; edge_lid++) {
        mesh::id_t edge_gid = _pmesh->edge_gid(edge_lid);        
        mesh::id_t dof_id = map_dof(edge_lid);
        if (dof_id != -1) {
            if (_pmesh->is_owned_edge(edge_gid)) {                
                Edge* edge = _mesh->get_edge(edge_lid);                           
                nodepair dof(_pmesh->point_gid(edge->get_start_id()), _pmesh->point_gid(edge->get_end_id()));
                my_dofs.push_back(dof_with_ids(dof, dof_id, ID_NONE));                
            }        
        }
    }
    
    for (mesh::id_t i = 0; i < my_dofs.size(); i++) {        
//         rDebug("my_dofs[%2d] = ((%2d, %2d), %2d, ?)", i, my_dofs[i].first.first, my_dofs[i].first.second, my_dofs[i].second);
    }
    
    
    
    vector<int> recvcnts, displs;    
    int buf_size = parallel::set_counts_displs(my_dofs.size(), recvcnts, displs, _pmesh->get_comm());
    
    std::vector<dof_with_ids> all_dofs(buf_size);
    //Do the byte stuff.
    vector<int> recvcnts_bytes, displs_bytes;
    for (mesh::id_t i = 0; i < recvcnts.size(); i++) {
        recvcnts_bytes.push_back(recvcnts[i] * sizeof(dof_with_ids));
        displs_bytes.push_back(displs[i] * sizeof(dof_with_ids));
    }
    MPI_Allgatherv(&my_dofs[0], my_dofs.size() * sizeof(dof_with_ids), MPI_BYTE, 
                   &all_dofs[0], &recvcnts_bytes[0], &displs_bytes[0], MPI_BYTE, _pmesh->get_comm());
    
    //Wrong!                                    
  /*  for (mesh::id_t i = 0; i < all_dofs.size(); i++) {            
        rDebug("Created all_dofs[%2d] = ((%2d, %2d), %2d, ?)", i, all_dofs[i].first.first, all_dofs[i].first.second, all_dofs[i].second);
    }
  */      
    std::sort(all_dofs.begin(), all_dofs.end(), less_node_lex); //Sort by point ID's
    
    for (mesh::id_t i = 0; i < all_dofs.size(); i++) {
        dof_with_ids d = all_dofs[i];
        all_dofs[i] = dof_with_ids(d.first, d.second, i);
//         rDebug("Node lex-sorted all_dofs[%2d] = ((%2d, %2d), %2d, %2d)", i, all_dofs[i].first.first, all_dofs[i].first.second, all_dofs[i].second, all_dofs[i].third);
    }
    
    std::sort(all_dofs.begin(), all_dofs.end(), less_id); //Sort by 2nd column.
    
/*    for (mesh::id_t i = 0; i < all_dofs.size(); i++) {
        rDebug("ID-sorted all_dofs[%2d] = ((%2d, %2d), %2d, %2d)", i, all_dofs[i].first.first, all_dofs[i].first.second, all_dofs[i].second, all_dofs[i].third);
    }*/
    
    //Create sorted A and M.
    
    
    rDebug("Generate sorted A and M matrices");  
    sort_matrix(A, A_s, all_dofs);      
//     cout << A_s;
    
//     rDebug("Generate sorted M");
    sort_matrix(M, M_s, all_dofs);
//     cout << M_s;
    
}
    
void NedelecMesh::sort_matrix(const Epetra_FECrsMatrix* M, Epetra_FECrsMatrix* M_s, const vector<dof_with_ids>& all_dofs) {
    int info;    
    for (int lrid = 0; lrid < M->NumMyRows(); lrid++) {        
        int grid = M->GRID(lrid);        
        assert (grid < all_dofs.size() && grid > -1);        
        int sgrid = all_dofs[grid].third;
        assert (sgrid < all_dofs.size() && sgrid > -1);        
//         int nof_entries = M->NumMyEntries(lrid);        
        int length;
        double* values;// = new double[nof_entries];
        int* indices;// = new int[nof_entries];
//         rDebugAll("lrid = %2d, grid = %2d, sgrid = %2d, length = %2d", lrid, grid, sgrid, length);        
//         info = M->ExtractMyRowCopy(lrid, nof_entries, length, values, indices);
        info = M->ExtractMyRowView(lrid, length, values, indices);
        assert(info >= 0);        
        
        for (int j = 0; j < length; j++) {            
            int gcid = M->GCID(indices[j]);
            assert (gcid < all_dofs.size() && gcid > -1);
            int sgcid = all_dofs[gcid].third;
            assert (sgcid < all_dofs.size() && sgcid > -1);
            double val = values[j];            
            info = M_s->InsertGlobalValues(1, &sgrid, 1, &sgcid, &val);
//             M_s->InsertGlobalValues(sgrid, 1, &val, &sgcid);
            assert(info >= 0);
        }
//         delete[] values;
//         delete[] indices;                            
    }
    M_s->GlobalAssemble(false);    
}


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

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

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

mesh::id_t NedelecMesh::gdof(mesh::id_t ldof) const {
    if (ldof < _mesh->get_nof_edges()) {
        return _pmesh->edge_gid(ldof);
    } else if (ldof < 2*_mesh->get_nof_edges()) {
        return _pmesh->edge_gid(ldof - _mesh->get_nof_edges()) + get_num_global_dofs(1);
    } else {
        mesh::id_t face_lid = (ldof - 2*_mesh->get_nof_edges()) / 3;
        return 2*get_num_global_dofs(1) \
            + 3*_pmesh->face_gid(face_lid) \
            + (ldof - 2*_mesh->get_nof_edges()) % 3;
    }
}

bool NedelecMesh::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_edge(gdof);    
    } else if (gdof < 2*get_num_global_dofs(1)) {
        return _pmesh->is_owned_edge(gdof - get_num_global_dofs(1));
    } else {
        return _pmesh->is_owned_face((gdof - 2*get_num_global_dofs(1))/3);
    }
}

mesh::id_t NedelecMesh::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;
}

int NedelecMesh::map_dof_rev(int mapped_dof) const {
    assert(0);
    return _map_rev[mapped_dof];
}

mesh::Vector3 NedelecMesh::eval(const mesh::Vector3& x, const colarray::Vector<double>& q) const {
  
    Vector3 curField;
    Vector3 sumField(0.0, 0.0, 0.0);
    
    // identify tetrahedr(on|a) containing x using octree data structure
    std::vector<mesh::Tet *> tets;
    _mesh->find_tets_by_point(x, tets);
    // warn if no x is not in mesh
    if (tets.size() == 0) {
        rWarning("eval: No tet found at location (%g,%g,%g).", x.x, x.y, x.z);
        return sumField;
    }

    // get the field for each tetrahedron
    for (unsigned int j = 0; j < tets.size(); j ++){
        mesh::Tet* tet = tets[j];
        // get corresponding NedelecElement
        NedelecElement* elem = get_element(tet->get_id());
        // build vector containing all local DOFs
        const int* dof_ids = elem->get_dof_ids(_mesh, tet);
        colarray::Vector<double> dofs(elem->get_nof_dofs());
        for (int i = 0; i < elem->get_nof_dofs(); i ++) {
            int mapped_id = map_dof(dof_ids[i]);
            if (mapped_id == -1)
                dofs(i) = 0.0;
            else
                dofs(i) = q(mapped_id);
        }
        // evaluate field
        curField = elem->eval(tet, x, dofs);
        sumField += curField;
    }
    // return the field-vector with average field value
    return sumField / static_cast<double>(tets.size());
}


mesh::Vector3 NedelecMesh::eval_curl(const mesh::Vector3& x, const colarray::Vector<double>& q) const {
  
    Vector3 curField;
    Vector3 sumField(0.0, 0.0, 0.0);
    
    // identify tetrahedr(on|a) cotaining using octree data structure
    std::vector<mesh::Tet *> tets;
    _mesh->find_tets_by_point(x, tets);
    // warn if no x is not in mesh
    if (tets.size() == 0) {
        rWarning("eval_curl: No tet found at location (%g,%g,%g).", x.x, x.y, x.z);
        return sumField;
    }

    // get the curl for each tetrahedron
    for (unsigned int j = 0; j < tets.size(); j ++){
        mesh::Tet* tet = tets[j];
        // get corresponding NedelecElement
        NedelecElement* elem = get_element(tet->get_id());
        // build vector containing all local DOFs
        const int* dof_ids = elem->get_dof_ids(_mesh, tet);
        colarray::Vector<double> dofs(elem->get_nof_dofs());
        for (int i = 0; i < elem->get_nof_dofs(); i ++) {
            int mapped_id = map_dof(dof_ids[i]);
            if (mapped_id == -1)
                dofs(i) = 0.0;
            else
                dofs(i) = q(mapped_id);
        }
        // evaluate curl
        curField = elem->eval_curl(tet, x, dofs);
        sumField += curField;
    }
    // return the curl-vector with average curl value
    return sumField / static_cast<double>(tets.size());
}

double NedelecMesh::q_factor(double lambda, const double* q) const {
    // pi
    const double pi = ::acos(-1.0);
    // speed of light [m / s]
    const double cLight = 2.99792458e8;
    // magnetic field constant [V * s / A / m]
    const double mu0 = 4.0*pi / 1e7;
    // electric field constant [A * s / V / m]
    const double eps0 = 8.8542e-12;
    // conductivity [1 / Ohm / m]
    const double sigmaCu = 6.26e7;    // value published in CERN summer school notes
    // const double sigmaCu = 5.8e7;  // value used at PSI for copper cavity

    double integral = 0.0;

    mesh::TetMesh::surface_iterator it;
    for (it = _mesh->surface_begin(); it != _mesh->surface_end(); ++ it) {
        Face* face = *it;

        // find tet and elem belonging to boundary face
        mesh::id_t tet = face->get_tet_id();
        NedelecElement* elem = get_element(tet);
        // dofs for elem
        colarray::Vector<double> dofs(elem->get_nof_dofs());
        for (int i = 0; i < elem->get_nof_dofs(); ++ i) {
            const int* dof_ids = elem->get_dof_ids(_mesh, _mesh->get_tet(tet));
            int mapped_id = map_dof(dof_ids[i]);
            if (mapped_id == -1)
                dofs(i) = 0.0;
            else
                dofs(i) = q[mapped_id];
        }
        // calculate int_(face)(curl N_i * curl N_j, ds)
        integral += elem->surface_integral_curl(_mesh->get_tet(tet), face, dofs);
    }
  
    const int nofSym = 0;           // nofSym is not necessary since it cancels out
    integral *= (1 << nofSym);
    double omega = ::sqrt(lambda) * cLight;
    // double f = omega / 2.0 / pi;
    double intJpower2 = integral / mu0 / mu0 / omega / omega;
    double delta = sqrt(2.0 / (mu0 * omega * sigmaCu));
    double reZs = 1.0 / sigmaCu / delta;
    double W = reZs / 2 * intJpower2;          // Wall power loss
    double U = (1 << nofSym) * eps0 / 2;       // stored total energy
    double Q = omega * U / W;                  // unloaded quality factor
#if 0
    cout << "lambda=" << lambda << endl;
    cout << "integral=" << integral << endl;
    cout << "omega=" << omega << endl;
    cout << "intJpower2=" << intJpower2 << endl;
    cout << "delta=" << delta << endl;
    cout << "reZs=" << reZs << endl;
    cout << "W=" << W << endl;
    cout << "U=" << U << endl;
#endif
    return Q;
}

---