src/fem/nedelecmesh.cpp

Go to the documentation of this file.
00001 /***************************************************************************
00002                           nedelecmesh.cpp  -  description
00003                              -------------------
00004     begin                : Fri Dec 19 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 <cassert>
00019 #include <cmath>
00020 #include <new>
00021 #include <vector>
00022 #include <iostream>
00023 #include <iomanip>
00024 #include <algorithm>
00025 
00026 #include "Epetra_Map.h"
00027 #include "Epetra_Vector.h"
00028 
00029 #include "rlog/rlog.h"
00030 
00031 #include "nedelecelement1.h"
00032 #include "nedelecelement2.h"
00033 #include "nedelecmesh.h"
00034 
00035 #include "matrix.h"
00036 
00037 #include "fmxxphysicomath.h"
00038 #include "fmxxtopology.h"
00039 
00040 using namespace std;
00041 using namespace colarray;
00042 using namespace mesh;
00043 
00044 typedef std::pair<mesh::id_t, mesh::id_t> nodepair;
00045 typedef utility::triple<nodepair, mesh::id_t, mesh::id_t> dof_with_ids;
00046 
00047 static bool less_node_lex(const dof_with_ids& a, const dof_with_ids& b) {
00048     if (a.first.first != b.first.first)                
00049         return a.first.first < b.first.first;                                    
00050     else
00051         return a.first.second < b.first.second;    
00052 }
00053 
00054 static bool less_id(const dof_with_ids& a, const dof_with_ids& b) {
00055     return a.second < b.second;
00056 }     
00057 
00058 NedelecMesh::NedelecMesh(mesh::TetMesh* mesh, int order, unsigned int sym_plane_config)
00059   : _mesh(mesh), _order(order){
00060     
00061     _pmesh = dynamic_cast<mesh::ParallelTetMesh*>(_mesh);
00062     
00063     // Order should be equal to 1 or 2.    
00064     assert(1 <= order && order <= 2);
00065     
00066     // Locally stored dofs
00067     _nof_orig_dofs.resize(order + 1);
00068     _nof_orig_dofs[1] = _mesh->get_nof_edges();
00069     switch (order) {
00070     case 1:
00071         _nof_orig_dofs[0] = _nof_orig_dofs[1];
00072         break;
00073     case 2:
00074         _nof_orig_dofs[2] = _mesh->get_nof_edges() + 3*_mesh->get_nof_faces();
00075         _nof_orig_dofs[0] = _nof_orig_dofs[1] + _nof_orig_dofs[2];
00076         break;
00077     }
00078     
00079     // Global dofs
00080     _num_global_dofs.resize(order + 1);
00081     _num_global_dofs[1] = _pmesh->get_num_global_edges();
00082     switch (order) {
00083     case 1:
00084         _num_global_dofs[0] = _num_global_dofs[1];
00085         break;
00086     case 2: 
00087         _num_global_dofs[2] = _pmesh->get_num_global_edges() + 3*_pmesh->get_num_global_faces();
00088         _num_global_dofs[0] = _num_global_dofs[1] + _num_global_dofs[2];
00089     }    
00090 
00091     if (order == 1)
00092         _element = new NedelecElement1();
00093     else // order == 2
00094         _element = new NedelecElement2();
00095   
00096     this->set_map(sym_plane_config);
00097 }
00098 
00099 NedelecMesh::~NedelecMesh() {
00100     delete _element;
00101 }
00102 
00103 void NedelecMesh::assembleAM(Epetra_FECrsMatrix* A,
00104                              Epetra_FECrsMatrix* M) {
00105     int info;
00106     Matrix<double> Ae, Me;
00107 
00108     // loop over locally stored elements
00109     TetMesh::tet_iterator it;
00110     for (it = _mesh->tet_begin(); it != _mesh->tet_end(); ++ it) {
00111         Tet& tet = *it;
00112         NedelecElement* elem = get_element(tet.get_id());
00113         int nof_local_dofs = elem->get_nof_dofs();
00114         const int *dof_ids = elem->get_dof_ids(_mesh, &tet);
00115 //        rDebugAll("dof_lids: %u %u %u %u %u %u",
00116 //                  dof_ids[0], dof_ids[1], dof_ids[2], dof_ids[3], dof_ids[4], dof_ids[5]);
00117 //        rDebugAll("dof_mgids: %d %d %d %d %d %d",
00118 //                  map_dof(dof_ids[0]), map_dof(dof_ids[1]), map_dof(dof_ids[2]), 
00119 //                  map_dof(dof_ids[3]), map_dof(dof_ids[4]), map_dof(dof_ids[5]));
00120 
00121         elem->get_element_matrix(&tet, Ae, Me);
00122        
00123         for (int i_element = 0; i_element < nof_local_dofs; i_element ++) {
00124             int i_mapped = this->map_dof(dof_ids[i_element]);
00125             if (i_mapped != -1) {
00126     
00127                 for (int j_element = 0; j_element < nof_local_dofs; j_element ++) {
00128                     int j_mapped = this->map_dof(dof_ids[j_element]);
00129                     if (j_mapped != -1) {
00130                         double val;
00131             
00132                         val = Ae(i_element, j_element);     
00133                         if (val != 0.0) {
00134                             info = A->SumIntoGlobalValues(1, &i_mapped, 1, &j_mapped, &val);
00135                             assert(info >= 0);
00136                             if (info > 0) {
00137                                 info = A->InsertGlobalValues(1, &i_mapped, 1, &j_mapped, &val);
00138                                 assert(info >= 0);
00139                             }                            
00140                         }
00141 
00142                         val = Me(i_element, j_element);
00143                         if (val != 0.0) {            
00144                             info = M->SumIntoGlobalValues(1, &i_mapped, 1, &j_mapped, &val);
00145                             assert(info >= 0);                            
00146                             if (info > 0) {
00147                                 info = M->InsertGlobalValues(1, &i_mapped, 1, &j_mapped, &val);
00148                                 assert(info >= 0);
00149                             }                            
00150 
00151                         }
00152                     }
00153                 }
00154             }
00155         }    
00156     }
00157     A->GlobalAssemble(false);
00158     M->GlobalAssemble(false);
00159 }
00160 
00161 void NedelecMesh::set_map(unsigned int sym_plane_config) {
00162     // Create elimination array
00163     vector<bool> elim(get_num_my_dofs(0));
00164     for (mesh::id_t e = 0; e < _mesh->get_nof_tets(); e ++) {
00165         NedelecElement* element = get_element(e);
00166         mesh::Tet* tet = _mesh->get_tet(e);
00167         const int* dof_ids = element->get_dof_ids(_mesh, tet);
00168         const bool* elim_dofs = element->get_eliminated_dofs(_mesh, tet, sym_plane_config);
00169         for (int k = 0; k < element->get_nof_dofs(); k ++)
00170             elim[dof_ids[k]] = elim_dofs[k];
00171     }
00172 
00173     if (_order > 1) {
00174         // Eliminate additional dofs due to linear dependece of face dofs
00175         for (mesh::id_t i = 0; i < _mesh->get_nof_faces(); i ++)
00176             elim[2*_mesh->get_nof_edges() + 3*i] = true;
00177     }
00178 
00179     _nof_mapped_dofs.resize(_order + 1);
00180     
00181     //
00182     // 1st order gdofs only
00183     //
00184     vector<mesh::id_t> eliminated_dofs;
00185     for (mesh::id_t i = 0; i < get_num_my_dofs(1); ++ i) {
00186         if(elim[i])
00187             eliminated_dofs.push_back(_pmesh->edge_gid(i));
00188     }
00189     parallel::compute_mgid_map(get_num_global_dofs(1), _pmesh->get_edge_gids(), eliminated_dofs, 
00190                                &_nof_mapped_dofs[1], _mgid_map, _pmesh->get_comm());
00191                                
00192     if (_order == 1) {
00193         _nof_mapped_dofs[0] = _nof_mapped_dofs[1];
00194     } else {
00195         //
00196         // Recompute mgids using all gdofs
00197         //
00198         rAssert(_order == 2);
00199         for (mesh::id_t ldof = get_num_my_dofs(1); ldof < get_num_my_dofs(0); ++ ldof) {
00200             if(elim[ldof])
00201                 eliminated_dofs.push_back(gdof(ldof));
00202         }
00203         vector<mesh::id_t> all_my_gdofs;
00204         for (mesh::id_t ldof = 0; ldof < get_num_my_dofs(0); ++ ldof)
00205             all_my_gdofs.push_back(gdof(ldof));
00206         parallel::compute_mgid_map(get_num_global_dofs(0), all_my_gdofs, eliminated_dofs, 
00207                                    &_nof_mapped_dofs[0], _mgid_map, _pmesh->get_comm());
00208         _nof_mapped_dofs[2] = _nof_mapped_dofs[0] - _nof_mapped_dofs[1];
00209     }
00210     
00211     rDebugAll("nedelec: num_my_dofs=%d, num_my_eliminated_dofs=%d", 
00212               get_num_my_dofs(0) , static_cast<int>(eliminated_dofs.size()));
00213     rDebugAll("nedelec mesh: #gdofs=%d, #mgdofs=%d", get_num_global_dofs(0), _nof_mapped_dofs[0]);
00214     
00215     // debug
00216 #if 0
00217     for (std::map<mesh::id_t, mesh::id_t>::iterator it=_mgid_map.begin();
00218          it != _mgid_map.end();
00219          ++ it)
00220     {
00221         rDebugAll("_mgid_map[%d] = %d", (*it).first, (*it).second);
00222     }
00223     {
00224         MPI_Comm comm = MPI_COMM_WORLD;
00225         int rank, size;
00226         MPI_Comm_size(comm, &size);
00227         MPI_Comm_rank(comm, &rank);
00228         for (int p = 0; p < size; ++ p) {
00229             if (p == rank) {
00230                 for (int i = 0; i < num_my_dofs; ++ i) {
00231                     rDebugAll("--lid:%d, gid:%d, mgid:%d", i, _pmesh->edge_gid(i), map_dof(i));
00232                 }
00233                 fsync(1);
00234             }
00235             MPI_Barrier(comm);
00236         }
00237     }
00238     rDebugAll("edge mgid to point gid mapping...");
00239     for (mesh::id_t mgid = 0; mgid < num_global_mgids; ++ mgid) {
00240         // Find gid
00241         std::map<mesh::id_t, mesh::id_t>::iterator found;
00242         found = std::lower_bound(_mgid_map.begin(), _mgid_map.end(), mgid, less_than_value);
00243         if (found != _mgid_map.end() && (*found).second == mgid) {
00244             mesh::id_t edge_gid = (*found).first;
00245             // Find edge lid
00246             mesh::id_t edge_lid;
00247             for (edge_lid = 0; _pmesh->edge_gid(edge_lid) != edge_gid; ++ edge_lid);
00248             assert(_pmesh->edge_gid(edge_lid) == edge_gid);
00249           
00250             rDebugAll("edge=mgid=%u, edge_gid=%u: point_gid0=%u, point_gid1=%u", 
00251                       mgid, 
00252                       edge_gid, 
00253                       _pmesh->point_gid(_pmesh->get_edge(edge_lid)->get_start_id()),
00254                       _pmesh->point_gid(_pmesh->get_edge(edge_lid)->get_end_id()));
00255         }
00256     }
00257 #endif
00258 }
00259 
00260 void NedelecMesh::generate_sorted_AM(const Epetra_FECrsMatrix* A, const Epetra_FECrsMatrix* M, Epetra_FECrsMatrix* A_s, Epetra_FECrsMatrix* M_s) {
00261     // Generate vector of point pairs
00262     // Equals owned edges not removed in 1st order case.
00263     vector<dof_with_ids> my_dofs;        
00264     
00265     mesh::id_t my_nof_edges = _pmesh->get_nof_edges();
00266     for (mesh::id_t edge_lid = 0; edge_lid < my_nof_edges; edge_lid++) {
00267         mesh::id_t edge_gid = _pmesh->edge_gid(edge_lid);        
00268         mesh::id_t dof_id = map_dof(edge_lid);
00269         if (dof_id != -1) {
00270             if (_pmesh->is_owned_edge(edge_gid)) {                
00271                 Edge* edge = _mesh->get_edge(edge_lid);                           
00272                 nodepair dof(_pmesh->point_gid(edge->get_start_id()), _pmesh->point_gid(edge->get_end_id()));
00273                 my_dofs.push_back(dof_with_ids(dof, dof_id, ID_NONE));                
00274             }        
00275         }
00276     }
00277     
00278     for (mesh::id_t i = 0; i < my_dofs.size(); i++) {        
00279 //         rDebug("my_dofs[%2d] = ((%2d, %2d), %2d, ?)", i, my_dofs[i].first.first, my_dofs[i].first.second, my_dofs[i].second);
00280     }
00281     
00282     
00283     
00284     vector<int> recvcnts, displs;    
00285     int buf_size = parallel::set_counts_displs(my_dofs.size(), recvcnts, displs, _pmesh->get_comm());
00286     
00287     std::vector<dof_with_ids> all_dofs(buf_size);
00288     //Do the byte stuff.
00289     vector<int> recvcnts_bytes, displs_bytes;
00290     for (mesh::id_t i = 0; i < recvcnts.size(); i++) {
00291         recvcnts_bytes.push_back(recvcnts[i] * sizeof(dof_with_ids));
00292         displs_bytes.push_back(displs[i] * sizeof(dof_with_ids));
00293     }
00294     MPI_Allgatherv(&my_dofs[0], my_dofs.size() * sizeof(dof_with_ids), MPI_BYTE, 
00295                    &all_dofs[0], &recvcnts_bytes[0], &displs_bytes[0], MPI_BYTE, _pmesh->get_comm());
00296     
00297     //Wrong!                                    
00298   /*  for (mesh::id_t i = 0; i < all_dofs.size(); i++) {            
00299         rDebug("Created all_dofs[%2d] = ((%2d, %2d), %2d, ?)", i, all_dofs[i].first.first, all_dofs[i].first.second, all_dofs[i].second);
00300     }
00301   */      
00302     std::sort(all_dofs.begin(), all_dofs.end(), less_node_lex); //Sort by point ID's
00303     
00304     for (mesh::id_t i = 0; i < all_dofs.size(); i++) {
00305         dof_with_ids d = all_dofs[i];
00306         all_dofs[i] = dof_with_ids(d.first, d.second, i);
00307 //         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);
00308     }
00309     
00310     std::sort(all_dofs.begin(), all_dofs.end(), less_id); //Sort by 2nd column.
00311     
00312 /*    for (mesh::id_t i = 0; i < all_dofs.size(); i++) {
00313         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);
00314     }*/
00315     
00316     //Create sorted A and M.
00317     
00318     
00319     rDebug("Generate sorted A and M matrices");  
00320     sort_matrix(A, A_s, all_dofs);      
00321 //     cout << A_s;
00322     
00323 //     rDebug("Generate sorted M");
00324     sort_matrix(M, M_s, all_dofs);
00325 //     cout << M_s;
00326     
00327 }
00328     
00329 void NedelecMesh::sort_matrix(const Epetra_FECrsMatrix* M, Epetra_FECrsMatrix* M_s, const vector<dof_with_ids>& all_dofs) {
00330     int info;    
00331     for (int lrid = 0; lrid < M->NumMyRows(); lrid++) {        
00332         int grid = M->GRID(lrid);        
00333         assert (grid < all_dofs.size() && grid > -1);        
00334         int sgrid = all_dofs[grid].third;
00335         assert (sgrid < all_dofs.size() && sgrid > -1);        
00336 //         int nof_entries = M->NumMyEntries(lrid);        
00337         int length;
00338         double* values;// = new double[nof_entries];
00339         int* indices;// = new int[nof_entries];
00340 //         rDebugAll("lrid = %2d, grid = %2d, sgrid = %2d, length = %2d", lrid, grid, sgrid, length);        
00341 //         info = M->ExtractMyRowCopy(lrid, nof_entries, length, values, indices);
00342         info = M->ExtractMyRowView(lrid, length, values, indices);
00343         assert(info >= 0);        
00344         
00345         for (int j = 0; j < length; j++) {            
00346             int gcid = M->GCID(indices[j]);
00347             assert (gcid < all_dofs.size() && gcid > -1);
00348             int sgcid = all_dofs[gcid].third;
00349             assert (sgcid < all_dofs.size() && sgcid > -1);
00350             double val = values[j];            
00351             info = M_s->InsertGlobalValues(1, &sgrid, 1, &sgcid, &val);
00352 //             M_s->InsertGlobalValues(sgrid, 1, &val, &sgcid);
00353             assert(info >= 0);
00354         }
00355 //         delete[] values;
00356 //         delete[] indices;                            
00357     }
00358     M_s->GlobalAssemble(false);    
00359 }
00360 
00361 
00362 mesh::id_t NedelecMesh::get_num_my_dofs(int order) const {
00363     return _nof_orig_dofs.at(order);
00364 }
00365 
00366 mesh::id_t NedelecMesh::get_num_global_mdofs(int order) const {
00367     return _nof_mapped_dofs.at(order);
00368 }
00369 
00370 mesh::id_t NedelecMesh::get_num_global_dofs(int order) const {
00371     return _num_global_dofs.at(order);
00372 }
00373 
00374 mesh::id_t NedelecMesh::gdof(mesh::id_t ldof) const {
00375     if (ldof < _mesh->get_nof_edges()) {
00376         return _pmesh->edge_gid(ldof);
00377     } else if (ldof < 2*_mesh->get_nof_edges()) {
00378         return _pmesh->edge_gid(ldof - _mesh->get_nof_edges()) + get_num_global_dofs(1);
00379     } else {
00380         mesh::id_t face_lid = (ldof - 2*_mesh->get_nof_edges()) / 3;
00381         return 2*get_num_global_dofs(1) \
00382             + 3*_pmesh->face_gid(face_lid) \
00383             + (ldof - 2*_mesh->get_nof_edges()) % 3;
00384     }
00385 }
00386 
00387 bool NedelecMesh::is_owned_gdof(mesh::id_t gdof) const {
00388     rAssert(0 <= gdof && gdof < get_num_global_dofs(0));
00389     
00390     if (gdof < get_num_global_dofs(1)) {
00391         return _pmesh->is_owned_edge(gdof);    
00392     } else if (gdof < 2*get_num_global_dofs(1)) {
00393         return _pmesh->is_owned_edge(gdof - get_num_global_dofs(1));
00394     } else {
00395         return _pmesh->is_owned_face((gdof - 2*get_num_global_dofs(1))/3);
00396     }
00397 }
00398 
00399 mesh::id_t NedelecMesh::map_dof(mesh::id_t ldof) const {
00400     std::map<mesh::id_t, mesh::id_t>::const_iterator found = _mgid_map.find(gdof(ldof));
00401     if (found == _mgid_map.end())
00402         return mesh::ID_NONE;
00403     else
00404         return (*found).second;
00405 }
00406 
00407 int NedelecMesh::map_dof_rev(int mapped_dof) const {
00408     assert(0);
00409     return _map_rev[mapped_dof];
00410 }
00411 
00412 mesh::Vector3 NedelecMesh::eval(const mesh::Vector3& x, const colarray::Vector<double>& q) const {
00413   
00414     Vector3 curField;
00415     Vector3 sumField(0.0, 0.0, 0.0);
00416     
00417     // identify tetrahedr(on|a) containing x using octree data structure
00418     std::vector<mesh::Tet *> tets;
00419     _mesh->find_tets_by_point(x, tets);
00420     // warn if no x is not in mesh
00421     if (tets.size() == 0) {
00422         rWarning("eval: No tet found at location (%g,%g,%g).", x.x, x.y, x.z);
00423         return sumField;
00424     }
00425 
00426     // get the field for each tetrahedron
00427     for (unsigned int j = 0; j < tets.size(); j ++){
00428         mesh::Tet* tet = tets[j];
00429         // get corresponding NedelecElement
00430         NedelecElement* elem = get_element(tet->get_id());
00431         // build vector containing all local DOFs
00432         const int* dof_ids = elem->get_dof_ids(_mesh, tet);
00433         colarray::Vector<double> dofs(elem->get_nof_dofs());
00434         for (int i = 0; i < elem->get_nof_dofs(); i ++) {
00435             int mapped_id = map_dof(dof_ids[i]);
00436             if (mapped_id == -1)
00437                 dofs(i) = 0.0;
00438             else
00439                 dofs(i) = q(mapped_id);
00440         }
00441         // evaluate field
00442         curField = elem->eval(tet, x, dofs);
00443         sumField += curField;
00444     }
00445     // return the field-vector with average field value
00446     return sumField / static_cast<double>(tets.size());
00447 }
00448 
00449 
00450 mesh::Vector3 NedelecMesh::eval_curl(const mesh::Vector3& x, const colarray::Vector<double>& q) const {
00451   
00452     Vector3 curField;
00453     Vector3 sumField(0.0, 0.0, 0.0);
00454     
00455     // identify tetrahedr(on|a) cotaining using octree data structure
00456     std::vector<mesh::Tet *> tets;
00457     _mesh->find_tets_by_point(x, tets);
00458     // warn if no x is not in mesh
00459     if (tets.size() == 0) {
00460         rWarning("eval_curl: No tet found at location (%g,%g,%g).", x.x, x.y, x.z);
00461         return sumField;
00462     }
00463 
00464     // get the curl for each tetrahedron
00465     for (unsigned int j = 0; j < tets.size(); j ++){
00466         mesh::Tet* tet = tets[j];
00467         // get corresponding NedelecElement
00468         NedelecElement* elem = get_element(tet->get_id());
00469         // build vector containing all local DOFs
00470         const int* dof_ids = elem->get_dof_ids(_mesh, tet);
00471         colarray::Vector<double> dofs(elem->get_nof_dofs());
00472         for (int i = 0; i < elem->get_nof_dofs(); i ++) {
00473             int mapped_id = map_dof(dof_ids[i]);
00474             if (mapped_id == -1)
00475                 dofs(i) = 0.0;
00476             else
00477                 dofs(i) = q(mapped_id);
00478         }
00479         // evaluate curl
00480         curField = elem->eval_curl(tet, x, dofs);
00481         sumField += curField;
00482     }
00483     // return the curl-vector with average curl value
00484     return sumField / static_cast<double>(tets.size());
00485 }
00486 
00487 double NedelecMesh::q_factor(double lambda, const double* q) const {
00488     // pi
00489     const double pi = ::acos(-1.0);
00490     // speed of light [m / s]
00491     const double cLight = 2.99792458e8;
00492     // magnetic field constant [V * s / A / m]
00493     const double mu0 = 4.0*pi / 1e7;
00494     // electric field constant [A * s / V / m]
00495     const double eps0 = 8.8542e-12;
00496     // conductivity [1 / Ohm / m]
00497     const double sigmaCu = 6.26e7;    // value published in CERN summer school notes
00498     // const double sigmaCu = 5.8e7;  // value used at PSI for copper cavity
00499 
00500     double integral = 0.0;
00501 
00502     mesh::TetMesh::surface_iterator it;
00503     for (it = _mesh->surface_begin(); it != _mesh->surface_end(); ++ it) {
00504         Face* face = *it;
00505 
00506         // find tet and elem belonging to boundary face
00507         mesh::id_t tet = face->get_tet_id();
00508         NedelecElement* elem = get_element(tet);
00509         // dofs for elem
00510         colarray::Vector<double> dofs(elem->get_nof_dofs());
00511         for (int i = 0; i < elem->get_nof_dofs(); ++ i) {
00512             const int* dof_ids = elem->get_dof_ids(_mesh, _mesh->get_tet(tet));
00513             int mapped_id = map_dof(dof_ids[i]);
00514             if (mapped_id == -1)
00515                 dofs(i) = 0.0;
00516             else
00517                 dofs(i) = q[mapped_id];
00518         }
00519         // calculate int_(face)(curl N_i * curl N_j, ds)
00520         integral += elem->surface_integral_curl(_mesh->get_tet(tet), face, dofs);
00521     }
00522   
00523     const int nofSym = 0;           // nofSym is not necessary since it cancels out
00524     integral *= (1 << nofSym);
00525     double omega = ::sqrt(lambda) * cLight;
00526     // double f = omega / 2.0 / pi;
00527     double intJpower2 = integral / mu0 / mu0 / omega / omega;
00528     double delta = sqrt(2.0 / (mu0 * omega * sigmaCu));
00529     double reZs = 1.0 / sigmaCu / delta;
00530     double W = reZs / 2 * intJpower2;          // Wall power loss
00531     double U = (1 << nofSym) * eps0 / 2;       // stored total energy
00532     double Q = omega * U / W;                  // unloaded quality factor
00533 #if 0
00534     cout << "lambda=" << lambda << endl;
00535     cout << "integral=" << integral << endl;
00536     cout << "omega=" << omega << endl;
00537     cout << "intJpower2=" << intJpower2 << endl;
00538     cout << "delta=" << delta << endl;
00539     cout << "reZs=" << reZs << endl;
00540     cout << "W=" << W << endl;
00541     cout << "U=" << U << endl;
00542 #endif
00543     return Q;
00544 }

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