second derivative of \sum_ipt wprime * B(u)
| 737 | |
| 738 | // second derivative of \sum_ipt wprime * B(u) |
| 739 | void CalcHessianAdd (const FiniteElement & inner_fel, |
| 740 | const BaseMappedIntegrationRule & mir, |
| 741 | SliceMatrix<> wprime, |
| 742 | BareSliceVector<> elvecu, |
| 743 | SliceMatrix<> hessian, |
| 744 | LocalHeap & lh) const override |
| 745 | { |
| 746 | // a first simple implementation by numerical differentiation .... |
| 747 | static Timer t("interpolateDiffOp, Hessian"); |
| 748 | RegionTracer reg(TaskManager::GetThreadId(), t); |
| 749 | HeapReset hr(lh); |
| 750 | |
| 751 | /* |
| 752 | size_t ndof = inner_fel.GetNDof(); |
| 753 | double eps = 1e-6; |
| 754 | FlatVector<> wprimevec(wprime.Height()*wprime.Width(), lh); |
| 755 | for (int i = 0, ii = 0; i < wprime.Height(); i++) |
| 756 | for (int j = 0; j < wprime.Width(); j++, ii++) |
| 757 | wprimevec(ii) = wprime(i,j) * mir[i].GetWeight(); |
| 758 | |
| 759 | FlatMatrix<double,ColMajor> bmatl(mir.Size()*diffop->Dim(), ndof, lh); |
| 760 | FlatMatrix<double,ColMajor> bmatr(mir.Size()*diffop->Dim(), ndof, lh); |
| 761 | FlatMatrix<double,ColMajor> dbmat(mir.Size()*diffop->Dim(), ndof, lh); |
| 762 | |
| 763 | for (size_t i = 0; i < ndof; i++) |
| 764 | { |
| 765 | FlatVector<> elvecur(ndof, lh), elvecul(ndof, lh); |
| 766 | elvecur = elvecu; |
| 767 | elvecul = elvecu; |
| 768 | elvecur(i) += eps; |
| 769 | elvecul(i) -= eps; |
| 770 | CalcLinearizedMatrix(inner_fel, mir, elvecul, bmatl, lh); |
| 771 | CalcLinearizedMatrix(inner_fel, mir, elvecur, bmatr, lh); |
| 772 | dbmat = 1/(2*eps) * (bmatr-bmatl); |
| 773 | hessian.Row(i) += Trans(dbmat) * wprimevec; |
| 774 | } |
| 775 | // cout << "hessian 1 = " << endl << hessian << endl; |
| 776 | */ |
| 777 | |
| 778 | const ElementTransformation & trafo = mir.GetTransformation(); |
| 779 | ElementId ei = trafo.GetElementId(); |
| 780 | auto & interpol_fel = fes->GetFE(ei, lh); |
| 781 | |
| 782 | FlatMatrix<double> elmat(interpol_fel.GetNDof(), lh); |
| 783 | elmat = 0.0; |
| 784 | bool symmetric_so_far = false; |
| 785 | try |
| 786 | { |
| 787 | for (auto & sbfi : single_bli) |
| 788 | sbfi->CalcElementMatrixAdd (interpol_fel, trafo, elmat, symmetric_so_far, lh); |
| 789 | } |
| 790 | catch (const ExceptionNOSIMD& e) |
| 791 | { |
| 792 | cout << IM(6) << e.What() << endl |
| 793 | << "switching to scalar evaluation" << endl; |
| 794 | for (auto & sbfi : single_bli) |
| 795 | sbfi->SetSimdEvaluate(false); |
| 796 | for (auto & sbfi : m3_bli) |
no test coverage detected