src/fem/nedelecelement2.cpp

Go to the documentation of this file.
00001 /***************************************************************************
00002                           nedelecelement2.cpp  -  description
00003                              -------------------
00004     begin                : Tue Dec 23 2003
00005     copyright            : (C) 2003 by Roman Geus
00006     email                : roman.geus@psi.ch
00007 ***************************************************************************/
00008 
00009 /***************************************************************************
00010  *                                                                         *
00011  *   This program is free software; you can redistribute it and/or modify  *
00012  *   it under the terms of the GNU General Public License as published by  *
00013  *   the Free Software Foundation; either version 2 of the License, or     *
00014  *   (at your option) any later version.                                   *
00015  *                                                                         *
00016  ***************************************************************************/
00017 
00018 #include "nedelecelement2.h"
00019 
00020 using namespace colarray;
00021 using namespace mesh;
00022 
00023 NedelecElement2::NedelecElement2() {
00024 }
00025 
00026 NedelecElement2::~NedelecElement2() {
00027 }
00028 
00029 const int* NedelecElement2::get_dof_ids(mesh::TetMesh* m, mesh::Tet* tet) const {
00030     static int dof_ids[24];
00031     using mesh::id_t;
00032 
00033     // loop over local edges
00034     for (int i = 0; i < 6; i ++) {
00035         dof_ids[i] = tet->get_edge_id(i);
00036         dof_ids[i+6] = dof_ids[i] + m->get_nof_edges();
00037     }
00038     // loop over local faces
00039     for (int i = 0; i < 4; i ++) {
00040         id_t face_id = tet->get_face_id(i);
00041         int base_face_dof = 3*face_id + 2*m->get_nof_edges();
00042         id_t corners[3];
00043         if (i == 0) {
00044             corners[0] = tet->get_corner_id(0);
00045             corners[1] = tet->get_corner_id(1);
00046             corners[2] = tet->get_corner_id(2);
00047         } else if (i == 1) {
00048             corners[0] = tet->get_corner_id(1);
00049             corners[1] = tet->get_corner_id(2);
00050             corners[2] = tet->get_corner_id(3);
00051         } else if (i == 2) {
00052             corners[0] = tet->get_corner_id(2);
00053             corners[1] = tet->get_corner_id(3);
00054             corners[2] = tet->get_corner_id(0);
00055         } else if (i == 3) {
00056             corners[0] = tet->get_corner_id(3);
00057             corners[1] = tet->get_corner_id(0);
00058             corners[2] = tet->get_corner_id(1);
00059         }
00060         // loop over all three face corners
00061         for (int k = 0; k < 3; k ++) {
00062             id_t p_id = m->get_face(face_id)->get_corner_id(k);
00063             int l;
00064             for (l = 0; l < 3; l ++)
00065                 if (corners[l] == p_id)
00066                     break;
00067             assert(l < 3);
00068             dof_ids[i + 12+ 4*l] = base_face_dof + k;
00069         }
00070     }
00071 #if 0
00072     for (int k = 0; k < 24; k ++)
00073         std::cout << dof_ids[k] << " ";
00074     std::cout << "\n";
00075 #endif
00076   
00077     return dof_ids;
00078 }
00079 
00080 const bool* NedelecElement2::get_eliminated_dofs(mesh::TetMesh* m,
00081                                                  mesh::Tet* tet,
00082                                                  unsigned int sym_plane_config) const {
00083     static bool eliminated_dofs[24];
00084     // loop over local edges
00085     for (int i = 0; i < 6; i ++) {
00086         if (tet->get_edge(i)->on_boundary()) {
00087             eliminated_dofs[i+6] = eliminated_dofs[i] = true;
00088         } else if ((sym_plane_config ^ 0x7) & tet->get_edge(i)->on_symmetry_planes()) {
00089             eliminated_dofs[i+6] = eliminated_dofs[i] = true;
00090         } else {
00091             eliminated_dofs[i+6] = eliminated_dofs[i] = false;
00092         }
00093     }
00094     // loop over local faces
00095     for (int i = 0; i < 4; i ++) {
00096         if (tet->get_face(i)->on_boundary()) {
00097             // eliminate all three dofs if face if on boundary
00098             eliminated_dofs[i + 12] = true;
00099             eliminated_dofs[i + 16] = true;
00100             eliminated_dofs[i + 20] = true;      
00101         } else if ((sym_plane_config ^ 0x7) & tet->get_face(i)->on_symmetry_planes()) {
00102             // eliminate all three dofs if face if on symmetry plane with E x n = 0
00103             eliminated_dofs[i + 12] = true;
00104             eliminated_dofs[i + 16] = true;
00105             eliminated_dofs[i + 20] = true;
00106         } else {
00107             // otherwise arbitrarily eliminate one of the three dofs, due to linear dependece
00108             eliminated_dofs[i + 12] = false;
00109             eliminated_dofs[i + 16] = false;
00110             eliminated_dofs[i + 20] = false;
00111         }
00112     }
00113     return eliminated_dofs;
00114 }
00115 
00116 void NedelecElement2::get_element_matrix(mesh::Tet* tet,
00117                                          Matrix<double>& Ae,
00118                                          Matrix<double>& Me) {
00119 
00120     if (Ae.get_m() != 24 || Ae.get_n() != 24) {
00121         // inappropriate Matrix shape: reassign
00122         Ae = Matrix<double>(24, 24);
00123     }
00124     if (Me.get_m() != 24 || Me.get_n() != 24) {
00125         // inappropriate Matrix shape: reassign
00126         Me = Matrix<double>(24, 24);
00127     }
00128 
00129 #ifndef GAUSSIANQUADRATURE  
00130     Tet::GradientMatrix J(*tet);
00131     get_Ae(tet, 2, J, Ae);
00132     get_Me(tet, 2, J, Me);
00133 #else
00134     get_AeMe(tet, 2, Ae, Me);
00135 #endif
00136 }
00137 
00138 mesh::Vector3 NedelecElement2::eval(Tet* tet, const Vector3& x, const Vector<double>& dof) {
00139     mesh::Tet::GradientMatrix J(*tet);
00140     return eval_cartesian(tet, 2, dof, J, x);
00141 }
00142 mesh::Vector3 NedelecElement2::eval_curl(Tet* tet, const Vector3& x, const Vector<double>& dof) {
00143     mesh::Tet::GradientMatrix J(*tet);
00144     return eval_curl_cartesian(tet, 2, dof, J, x);
00145 }

Generated on Fri Oct 26 13:35:12 2007 for FEMAXX (Finite Element Maxwell Eigensolver) by  doxygen 1.4.7