| 642 | } |
| 643 | |
| 644 | void Transfer2TPMesh(const CoefficientFunction * cfstd, GridFunction* gfutp,LocalHeap & lh) |
| 645 | { |
| 646 | static Timer tall("TPHighOrderFESpace::Transfer2TPMesh"); RegionTimer rall(tall); |
| 647 | const shared_ptr<TPHighOrderFESpace> tpfes = dynamic_pointer_cast<TPHighOrderFESpace>(gfutp->GetFESpace()); |
| 648 | const shared_ptr<FESpace> & fesx = tpfes->Space(-1); |
| 649 | const shared_ptr<FESpace> & fesy = tpfes->Space(0); |
| 650 | IterateElementsTP(*tpfes,VOL,lh, |
| 651 | [&] (ElementId ei0,ElementId ei1,LocalHeap & lh) |
| 652 | { |
| 653 | HeapReset hr(lh); |
| 654 | ArrayMem<int,2> ind(2); |
| 655 | ind[0] = ei0.Nr(); ind[1] = ei1.Nr(); |
| 656 | int elnr = tpfes->GetIndex(ind); |
| 657 | TPHighOrderFE & tpfel = dynamic_cast<TPHighOrderFE&>(tpfes->GetFE(ElementId(elnr),lh)); |
| 658 | auto & felx = fesx->GetFE(ei0,lh); |
| 659 | auto & fely = fesy->GetFE(ei1,lh); |
| 660 | int ndofx = felx.GetNDof(); |
| 661 | int ndofy = fely.GetNDof(); |
| 662 | |
| 663 | ArrayMem<const IntegrationRule *,2> irs(tpfel.elements.Size()); |
| 664 | for(int s=0;s<irs.Size();s++) |
| 665 | irs[s] = &SelectIntegrationRule(tpfel.elements[s]->ElementType(),2*tpfel.elements[s]->Order()); |
| 666 | TPIntegrationRule ir(irs); |
| 667 | const ElementTransformation & tptrafo = tpfes->GetTrafo(ElementId(elnr),lh); |
| 668 | const ElementTransformation & trafox = fesx->GetMeshAccess()->GetTrafo(ei0,lh); |
| 669 | const ElementTransformation & trafoy = fesy->GetMeshAccess()->GetTrafo(ei1,lh);; |
| 670 | TPMappedIntegrationRule & tpmir = dynamic_cast<TPMappedIntegrationRule & >(tptrafo(ir, lh)); |
| 671 | BaseMappedIntegrationRule & mirx = trafox(*irs[0],lh); |
| 672 | BaseMappedIntegrationRule & miry = trafoy(*irs[1],lh); |
| 673 | int tpnip = irs[0]->Size()*irs[1]->Size(); |
| 674 | |
| 675 | FlatMatrix<> result(tpnip,fesy->GetEvaluator(VOL)->Dim()*fesx->GetEvaluator(VOL)->Dim(),lh); |
| 676 | cfstd->Evaluate(tpmir,result); |
| 677 | FlatMatrix<double,ColMajor> shapesx(mirx.Size()*fesx->GetEvaluator(VOL)->Dim(),ndofx,lh); |
| 678 | FlatMatrix<double,ColMajor> shapesy(miry.Size()*fesy->GetEvaluator(VOL)->Dim(),ndofy,lh); |
| 679 | FlatMatrix<double,ColMajor> shapesx2(mirx.Size()*fesx->GetEvaluator(VOL)->Dim(),ndofx,lh); |
| 680 | FlatMatrix<double,ColMajor> shapesy2(miry.Size()*fesy->GetEvaluator(VOL)->Dim(),ndofy,lh); |
| 681 | |
| 682 | fesx->GetEvaluator(VOL)->CalcMatrix(felx,mirx,shapesx,lh); |
| 683 | fesy->GetEvaluator(VOL)->CalcMatrix(fely,miry,shapesy,lh); |
| 684 | shapesx2 = shapesx; |
| 685 | shapesy2 = shapesy; |
| 686 | for(int i=0;i<mirx.Size();i++) |
| 687 | shapesx.Rows(i*fesx->GetEvaluator(VOL)->Dim(),(i+1)*fesx->GetEvaluator(VOL)->Dim())*=mirx[i].GetWeight(); |
| 688 | for(int i=0;i<miry.Size();i++) |
| 689 | shapesy.Rows(i*fesy->GetEvaluator(VOL)->Dim(),(i+1)*fesy->GetEvaluator(VOL)->Dim())*=miry[i].GetWeight(); |
| 690 | |
| 691 | FlatMatrix<> resultasmat(mirx.Size(),miry.Size()*fesy->GetEvaluator(VOL)->Dim(),&result(0)); |
| 692 | FlatMatrix<> elmatrhs(ndofx,ndofy,lh); |
| 693 | elmatrhs = Trans(shapesx)*resultasmat*shapesy; |
| 694 | |
| 695 | FlatMatrix<> elmatx(ndofx,lh),elmaty(ndofy,lh); |
| 696 | elmatx = Trans(shapesx2)*shapesx; |
| 697 | elmaty = Trans(shapesy2)*shapesy; |
| 698 | CalcInverse(elmatx); |
| 699 | CalcInverse(elmaty); |
| 700 | FlatMatrix<> elmatout(ndofx,ndofy,lh); |
| 701 | elmatout = elmatx*elmatrhs*elmaty; |
no test coverage detected