| 534 | |
| 535 | template <class SCAL> |
| 536 | void H1AMG_Preconditioner<SCAL> :: AddElementMatrix (FlatArray<int> dnums, |
| 537 | FlatMatrix<SCAL> elmat, |
| 538 | ElementId id, |
| 539 | LocalHeap & lh) |
| 540 | { |
| 541 | HeapReset hr(lh); |
| 542 | // vertex weights |
| 543 | // static Timer t("h1amg - addelmat"); |
| 544 | // static Timer t1("h1amg - addelmat calc v-schur"); |
| 545 | // static Timer t3("h1amg - addelmat calc e-schur"); |
| 546 | // static Timer t5("h1amg - addelmat invert"); |
| 547 | |
| 548 | // RegionTimer reg (t); |
| 549 | |
| 550 | size_t ndof = dnums.Size(); |
| 551 | BitArray used(ndof, lh); |
| 552 | |
| 553 | FlatMatrix<SCAL> ext_elmat(ndof+1, ndof+1, lh); |
| 554 | |
| 555 | { |
| 556 | // RegionTimer reg (t5); |
| 557 | ext_elmat.Rows(0,ndof).Cols(0,ndof) = elmat; |
| 558 | ext_elmat.Row(ndof) = 1; |
| 559 | ext_elmat.Col(ndof) = 1; |
| 560 | ext_elmat(ndof, ndof) = 0; |
| 561 | CalcInverse (ext_elmat); |
| 562 | } |
| 563 | |
| 564 | { |
| 565 | // RegionTimer reg (t1); |
| 566 | for (size_t i = 0; i < dnums.Size(); i++) |
| 567 | { |
| 568 | Mat<2,2,SCAL> ai; |
| 569 | ai(0,0) = ext_elmat(i,i); |
| 570 | ai(0,1) = ai(1,0) = ext_elmat(i, ndof); |
| 571 | ai(1,1) = ext_elmat(ndof, ndof); |
| 572 | ai = Inv(ai); |
| 573 | double weight = fabs(ai(0,0)); |
| 574 | vertex_weights_ht.Do(IVec<1>(dnums[i]), [weight] (auto & v) { v += weight; }); |
| 575 | } |
| 576 | } |
| 577 | { |
| 578 | // RegionTimer reg (t3); |
| 579 | for (size_t i = 0; i < dnums.Size(); i++) |
| 580 | for (size_t j = 0; j < i; j++) |
| 581 | { |
| 582 | Mat<3,3,SCAL> ai; |
| 583 | ai(0,0) = ext_elmat(i,i); |
| 584 | ai(1,1) = ext_elmat(j,j); |
| 585 | ai(0,1) = ai(1,0) = ext_elmat(i,j); |
| 586 | ai(2,2) = ext_elmat(ndof,ndof); |
| 587 | ai(0,2) = ai(2,0) = ext_elmat(i,ndof); |
| 588 | ai(1,2) = ai(2,1) = ext_elmat(j,ndof); |
| 589 | ai = Inv(ai); |
| 590 | double weight = fabs(ai(0,0)); |
| 591 | edge_weights_ht.Do(IVec<2>(dnums[j], dnums[i]).Sort(), [weight] (auto & v) { v += weight; }); |
| 592 | } |
| 593 | } |