00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
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
00064 assert(1 <= order && order <= 2);
00065
00066
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
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
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
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
00116
00117
00118
00119
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
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
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
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
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
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
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
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
00262
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
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
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
00298
00299
00300
00301
00302 std::sort(all_dofs.begin(), all_dofs.end(), less_node_lex);
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
00308 }
00309
00310 std::sort(all_dofs.begin(), all_dofs.end(), less_id);
00311
00312
00313
00314
00315
00316
00317
00318
00319 rDebug("Generate sorted A and M matrices");
00320 sort_matrix(A, A_s, all_dofs);
00321
00322
00323
00324 sort_matrix(M, M_s, all_dofs);
00325
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
00337 int length;
00338 double* values;
00339 int* indices;
00340
00341
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
00353 assert(info >= 0);
00354 }
00355
00356
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
00418 std::vector<mesh::Tet *> tets;
00419 _mesh->find_tets_by_point(x, tets);
00420
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
00427 for (unsigned int j = 0; j < tets.size(); j ++){
00428 mesh::Tet* tet = tets[j];
00429
00430 NedelecElement* elem = get_element(tet->get_id());
00431
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
00442 curField = elem->eval(tet, x, dofs);
00443 sumField += curField;
00444 }
00445
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
00456 std::vector<mesh::Tet *> tets;
00457 _mesh->find_tets_by_point(x, tets);
00458
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
00465 for (unsigned int j = 0; j < tets.size(); j ++){
00466 mesh::Tet* tet = tets[j];
00467
00468 NedelecElement* elem = get_element(tet->get_id());
00469
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
00480 curField = elem->eval_curl(tet, x, dofs);
00481 sumField += curField;
00482 }
00483
00484 return sumField / static_cast<double>(tets.size());
00485 }
00486
00487 double NedelecMesh::q_factor(double lambda, const double* q) const {
00488
00489 const double pi = ::acos(-1.0);
00490
00491 const double cLight = 2.99792458e8;
00492
00493 const double mu0 = 4.0*pi / 1e7;
00494
00495 const double eps0 = 8.8542e-12;
00496
00497 const double sigmaCu = 6.26e7;
00498
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
00507 mesh::id_t tet = face->get_tet_id();
00508 NedelecElement* elem = get_element(tet);
00509
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
00520 integral += elem->surface_integral_curl(_mesh->get_tet(tet), face, dofs);
00521 }
00522
00523 const int nofSym = 0;
00524 integral *= (1 << nofSym);
00525 double omega = ::sqrt(lambda) * cLight;
00526
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;
00531 double U = (1 << nofSym) * eps0 / 2;
00532 double Q = omega * U / W;
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 }