src/fem/nedelecelement2.cpp

Go to the documentation of this file.

/***************************************************************************
                          nedelecelement2.cpp  -  description
                             -------------------
    begin                : Tue Dec 23 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 "nedelecelement2.h"

using namespace colarray;
using namespace mesh;

NedelecElement2::NedelecElement2() {
}

NedelecElement2::~NedelecElement2() {
}

const int* NedelecElement2::get_dof_ids(mesh::TetMesh* m, mesh::Tet* tet) const {
    static int dof_ids[24];
    using mesh::id_t;

    // loop over local edges
    for (int i = 0; i < 6; i ++) {
        dof_ids[i] = tet->get_edge_id(i);
        dof_ids[i+6] = dof_ids[i] + m->get_nof_edges();
    }
    // loop over local faces
    for (int i = 0; i < 4; i ++) {
        id_t face_id = tet->get_face_id(i);
        int base_face_dof = 3*face_id + 2*m->get_nof_edges();
        id_t corners[3];
        if (i == 0) {
            corners[0] = tet->get_corner_id(0);
            corners[1] = tet->get_corner_id(1);
            corners[2] = tet->get_corner_id(2);
        } else if (i == 1) {
            corners[0] = tet->get_corner_id(1);
            corners[1] = tet->get_corner_id(2);
            corners[2] = tet->get_corner_id(3);
        } else if (i == 2) {
            corners[0] = tet->get_corner_id(2);
            corners[1] = tet->get_corner_id(3);
            corners[2] = tet->get_corner_id(0);
        } else if (i == 3) {
            corners[0] = tet->get_corner_id(3);
            corners[1] = tet->get_corner_id(0);
            corners[2] = tet->get_corner_id(1);
        }
        // loop over all three face corners
        for (int k = 0; k < 3; k ++) {
            id_t p_id = m->get_face(face_id)->get_corner_id(k);
            int l;
            for (l = 0; l < 3; l ++)
                if (corners[l] == p_id)
                    break;
            assert(l < 3);
            dof_ids[i + 12+ 4*l] = base_face_dof + k;
        }
    }
#if 0
    for (int k = 0; k < 24; k ++)
        std::cout << dof_ids[k] << " ";
    std::cout << "\n";
#endif
  
    return dof_ids;
}

const bool* NedelecElement2::get_eliminated_dofs(mesh::TetMesh* m,
                                                 mesh::Tet* tet,
                                                 unsigned int sym_plane_config) const {
    static bool eliminated_dofs[24];
    // loop over local edges
    for (int i = 0; i < 6; i ++) {
        if (tet->get_edge(i)->on_boundary()) {
            eliminated_dofs[i+6] = eliminated_dofs[i] = true;
        } else if ((sym_plane_config ^ 0x7) & tet->get_edge(i)->on_symmetry_planes()) {
            eliminated_dofs[i+6] = eliminated_dofs[i] = true;
        } else {
            eliminated_dofs[i+6] = eliminated_dofs[i] = false;
        }
    }
    // loop over local faces
    for (int i = 0; i < 4; i ++) {
        if (tet->get_face(i)->on_boundary()) {
            // eliminate all three dofs if face if on boundary
            eliminated_dofs[i + 12] = true;
            eliminated_dofs[i + 16] = true;
            eliminated_dofs[i + 20] = true;      
        } else if ((sym_plane_config ^ 0x7) & tet->get_face(i)->on_symmetry_planes()) {
            // eliminate all three dofs if face if on symmetry plane with E x n = 0
            eliminated_dofs[i + 12] = true;
            eliminated_dofs[i + 16] = true;
            eliminated_dofs[i + 20] = true;
        } else {
            // otherwise arbitrarily eliminate one of the three dofs, due to linear dependece
            eliminated_dofs[i + 12] = false;
            eliminated_dofs[i + 16] = false;
            eliminated_dofs[i + 20] = false;
        }
    }
    return eliminated_dofs;
}

void NedelecElement2::get_element_matrix(mesh::Tet* tet,
                                         Matrix<double>& Ae,
                                         Matrix<double>& Me) {

    if (Ae.get_m() != 24 || Ae.get_n() != 24) {
        // inappropriate Matrix shape: reassign
        Ae = Matrix<double>(24, 24);
    }
    if (Me.get_m() != 24 || Me.get_n() != 24) {
        // inappropriate Matrix shape: reassign
        Me = Matrix<double>(24, 24);
    }

#ifndef GAUSSIANQUADRATURE  
    Tet::GradientMatrix J(*tet);
    get_Ae(tet, 2, J, Ae);
    get_Me(tet, 2, J, Me);
#else
    get_AeMe(tet, 2, Ae, Me);
#endif
}

mesh::Vector3 NedelecElement2::eval(Tet* tet, const Vector3& x, const Vector<double>& dof) {
    mesh::Tet::GradientMatrix J(*tet);
    return eval_cartesian(tet, 2, dof, J, x);
}
mesh::Vector3 NedelecElement2::eval_curl(Tet* tet, const Vector3& x, const Vector<double>& dof) {
    mesh::Tet::GradientMatrix J(*tet);
    return eval_curl_cartesian(tet, 2, dof, J, x);
}

---