@@ -546,7 +546,7 @@ namespace dftefe
546
546
std::vector<ValueTypeOperator> basisOverlapClassicalBlock (
547
547
dofsPerCellCFE * dofsPerCellCFE);
548
548
std::vector<ValueTypeOperator> JxWxNCellConj (
549
- dofsPerCellCFE * nQuadPointInCellClassicalBlock);
549
+ dofsPerCellCFE * nQuadPointInCellClassicalBlock, 0 );
550
550
551
551
size_type stride = 0 ;
552
552
size_type m = 1 , n = dofsPerCellCFE,
@@ -588,6 +588,98 @@ namespace dftefe
588
588
n,
589
589
linAlgOpContext);
590
590
591
+ std::vector<ValueTypeOperator> basisOverlapECBlockEnrich (
592
+ dofsPerCell * numEnrichmentIdsInCell);
593
+
594
+ std::vector<ValueTypeOperator> basisOverlapECBlockClass (
595
+ dofsPerCellCFE * numEnrichmentIdsInCell);
596
+
597
+ if (numEnrichmentIdsInCell > 0 )
598
+ {
599
+ JxWxNCellConj.resize (
600
+ dofsPerCell * nQuadPointInCellEnrichmentBlockEnrichment, 0 );
601
+
602
+ m = 1 , n = dofsPerCell, k = nQuadPointInCellEnrichmentBlockEnrichment;
603
+
604
+ linearAlgebra::blasLapack::scaleStridedVarBatched<ValueTypeOperator,
605
+ ValueTypeOperator,
606
+ memorySpace>(
607
+ 1 ,
608
+ linearAlgebra::blasLapack::Layout::ColMajor,
609
+ linearAlgebra::blasLapack::ScalarOp::Identity,
610
+ linearAlgebra::blasLapack::ScalarOp::Conj,
611
+ &stride,
612
+ &stride,
613
+ &stride,
614
+ &m,
615
+ &n,
616
+ &k,
617
+ cellJxWValuesEnrichmentBlockEnrichment.data (),
618
+ cumulativeEnrichmentBlockEnrichmentDofQuadPoints,
619
+ JxWxNCellConj.data (),
620
+ linAlgOpContext);
621
+
622
+ linearAlgebra::blasLapack::
623
+ gemm<ValueTypeOperand, ValueTypeOperand, memorySpace>(
624
+ linearAlgebra::blasLapack::Layout::ColMajor,
625
+ linearAlgebra::blasLapack::Op::NoTrans,
626
+ linearAlgebra::blasLapack::Op::NoTrans,
627
+ dofsPerCell,
628
+ numEnrichmentIdsInCell,
629
+ k,
630
+ (ValueTypeOperand)1.0 ,
631
+ JxWxNCellConj.data (),
632
+ n,
633
+ enrichmentValuesVec.data (),
634
+ k,
635
+ (ValueTypeOperand)0.0 ,
636
+ basisOverlapECBlockEnrich.data (),
637
+ dofsPerCell,
638
+ linAlgOpContext);
639
+
640
+ JxWxNCellConj.resize (
641
+ dofsPerCellCFE * nQuadPointInCellEnrichmentBlockClassical, 0 );
642
+
643
+ m = 1 , n = dofsPerCellCFE,
644
+ k = nQuadPointInCellEnrichmentBlockClassical;
645
+
646
+ linearAlgebra::blasLapack::scaleStridedVarBatched<ValueTypeOperator,
647
+ ValueTypeOperator,
648
+ memorySpace>(
649
+ 1 ,
650
+ linearAlgebra::blasLapack::Layout::ColMajor,
651
+ linearAlgebra::blasLapack::ScalarOp::Identity,
652
+ linearAlgebra::blasLapack::ScalarOp::Conj,
653
+ &stride,
654
+ &stride,
655
+ &stride,
656
+ &m,
657
+ &n,
658
+ &k,
659
+ cellJxWValuesEnrichmentBlockClassical.data (),
660
+ cumulativeEnrichmentBlockClassicalDofQuadPoints,
661
+ JxWxNCellConj.data (),
662
+ linAlgOpContext);
663
+
664
+ linearAlgebra::blasLapack::
665
+ gemm<ValueTypeOperand, ValueTypeOperator, memorySpace>(
666
+ linearAlgebra::blasLapack::Layout::ColMajor,
667
+ linearAlgebra::blasLapack::Op::NoTrans,
668
+ linearAlgebra::blasLapack::Op::Trans,
669
+ n,
670
+ numEnrichmentIdsInCell,
671
+ k,
672
+ (ValueTypeOperand)1.0 ,
673
+ JxWxNCellConj.data (),
674
+ n,
675
+ classicalComponentInQuadValuesEC.data (),
676
+ numEnrichmentIdsInCell,
677
+ (ValueTypeOperator)0.0 ,
678
+ basisOverlapECBlockClass.data (),
679
+ n,
680
+ linAlgOpContext);
681
+ }
682
+
591
683
for (unsigned int iNode = 0 ; iNode < dofsPerCell; iNode++)
592
684
{
593
685
for (unsigned int jNode = 0 ; jNode < dofsPerCell; jNode++)
@@ -618,85 +710,93 @@ namespace dftefe
618
710
else if (iNode >= dofsPerCellCFE &&
619
711
jNode < dofsPerCellCFE && calculateWings)
620
712
{
621
- ValueTypeOperator NpiNcj = (ValueTypeOperator)0 ,
622
- ciNciNcj = (ValueTypeOperator)0 ;
623
- // Ni_pristine*Ni_classical at quadpoints
624
- for (unsigned int qPoint = 0 ;
625
- qPoint < nQuadPointInCellEnrichmentBlockEnrichment;
626
- qPoint++)
627
- {
628
- NpiNcj +=
629
- *(enrichmentValuesVec.data () +
630
- (iNode - dofsPerCellCFE) *
631
- nQuadPointInCellEnrichmentBlockEnrichment +
632
- qPoint) *
633
- *(cumulativeEnrichmentBlockEnrichmentDofQuadPoints +
634
- dofsPerCell * qPoint + jNode
635
- /* nQuadPointInCellEnrichmentBlockEnrichment *
636
- jNode +
637
- qPoint*/ ) *
638
- cellJxWValuesEnrichmentBlockEnrichment[qPoint];
639
- }
713
+ // ValueTypeOperator NpiNcj = (ValueTypeOperator)0,
714
+ // ciNciNcj = (ValueTypeOperator)0;
715
+ // // Ni_pristine*Ni_classical at quadpoints
716
+ // for (unsigned int qPoint = 0;
717
+ // qPoint < nQuadPointInCellEnrichmentBlockEnrichment;
718
+ // qPoint++)
719
+ // {
720
+ // NpiNcj +=
721
+ // *(enrichmentValuesVec.data() +
722
+ // (iNode - dofsPerCellCFE) *
723
+ // nQuadPointInCellEnrichmentBlockEnrichment +
724
+ // qPoint) *
725
+ // *(cumulativeEnrichmentBlockEnrichmentDofQuadPoints +
726
+ // dofsPerCell * qPoint + jNode
727
+ // /*nQuadPointInCellEnrichmentBlockEnrichment *
728
+ // jNode +
729
+ // qPoint*/) *
730
+ // cellJxWValuesEnrichmentBlockEnrichment[qPoint];
731
+ // }
640
732
641
- // Ni_classical using Mc = d quadrature * interpolated
642
- // ci's in Ni_classicalQuadrature of Mc = d
643
- for (unsigned int qPoint = 0 ;
644
- qPoint < nQuadPointInCellEnrichmentBlockClassical;
645
- qPoint++)
646
- {
647
- ciNciNcj +=
648
- classicalComponentInQuadValuesEC
649
- [numEnrichmentIdsInCell * qPoint +
650
- (iNode - dofsPerCellCFE)] *
651
- *(cumulativeEnrichmentBlockClassicalDofQuadPoints +
652
- dofsPerCellCFE * qPoint + jNode
653
- /* nQuadPointInCellEnrichmentBlockClassical *
654
- jNode + qPoint*/ ) *
655
- cellJxWValuesEnrichmentBlockClassical[qPoint];
656
- }
657
- *basisOverlapTmpIter += NpiNcj - ciNciNcj;
733
+ // // Ni_classical using Mc = d quadrature * interpolated
734
+ // // ci's in Ni_classicalQuadrature of Mc = d
735
+ // for (unsigned int qPoint = 0;
736
+ // qPoint < nQuadPointInCellEnrichmentBlockClassical;
737
+ // qPoint++)
738
+ // {
739
+ // ciNciNcj +=
740
+ // classicalComponentInQuadValuesEC
741
+ // [numEnrichmentIdsInCell * qPoint +
742
+ // (iNode - dofsPerCellCFE)] *
743
+ // *(cumulativeEnrichmentBlockClassicalDofQuadPoints +
744
+ // dofsPerCellCFE * qPoint + jNode
745
+ // /*nQuadPointInCellEnrichmentBlockClassical *
746
+ // jNode + qPoint*/) *
747
+ // cellJxWValuesEnrichmentBlockClassical[qPoint];
748
+ // }
749
+ // *basisOverlapTmpIter += NpiNcj - ciNciNcj;
750
+
751
+ *basisOverlapTmpIter =
752
+ *(basisOverlapECBlockEnrich.data () + (iNode - dofsPerCellCFE) * dofsPerCell + jNode) -
753
+ *(basisOverlapECBlockClass.data () + (iNode - dofsPerCellCFE) * dofsPerCellCFE + jNode);
658
754
}
659
755
660
756
else if (iNode < dofsPerCellCFE &&
661
757
jNode >= dofsPerCellCFE && calculateWings)
662
758
{
663
- ValueTypeOperator NciNpj = (ValueTypeOperator)0 ,
664
- NcicjNcj = (ValueTypeOperator)0 ;
665
- // Ni_pristine*Ni_classical at quadpoints
666
- for (unsigned int qPoint = 0 ;
667
- qPoint < nQuadPointInCellEnrichmentBlockEnrichment;
668
- qPoint++)
669
- {
670
- NciNpj +=
671
- *(cumulativeEnrichmentBlockEnrichmentDofQuadPoints +
672
- dofsPerCell * qPoint + iNode
673
- /* nQuadPointInCellEnrichmentBlockEnrichment *
674
- iNode +
675
- qPoint*/ ) *
676
- *(enrichmentValuesVec.data () +
677
- (jNode - dofsPerCellCFE) *
678
- nQuadPointInCellEnrichmentBlockEnrichment +
679
- qPoint) *
680
- cellJxWValuesEnrichmentBlockEnrichment[qPoint];
681
- }
759
+ // ValueTypeOperator NciNpj = (ValueTypeOperator)0,
760
+ // NcicjNcj = (ValueTypeOperator)0;
761
+ // // Ni_pristine*Ni_classical at quadpoints
762
+ // for (unsigned int qPoint = 0;
763
+ // qPoint < nQuadPointInCellEnrichmentBlockEnrichment;
764
+ // qPoint++)
765
+ // {
766
+ // NciNpj +=
767
+ // *(cumulativeEnrichmentBlockEnrichmentDofQuadPoints +
768
+ // dofsPerCell * qPoint + iNode
769
+ // /*nQuadPointInCellEnrichmentBlockEnrichment *
770
+ // iNode +
771
+ // qPoint*/) *
772
+ // *(enrichmentValuesVec.data() +
773
+ // (jNode - dofsPerCellCFE) *
774
+ // nQuadPointInCellEnrichmentBlockEnrichment +
775
+ // qPoint) *
776
+ // cellJxWValuesEnrichmentBlockEnrichment[qPoint];
777
+ // }
682
778
683
- // Ni_classical using Mc = d quadrature * interpolated
684
- // ci's in Ni_classicalQuadrature of Mc = d
685
- for (unsigned int qPoint = 0 ;
686
- qPoint < nQuadPointInCellEnrichmentBlockClassical;
687
- qPoint++)
688
- {
689
- NcicjNcj +=
690
- *(cumulativeEnrichmentBlockClassicalDofQuadPoints +
691
- dofsPerCellCFE * qPoint + iNode
692
- /* nQuadPointInCellEnrichmentBlockClassical *
693
- iNode + qPoint*/ ) *
694
- classicalComponentInQuadValuesEC
695
- [numEnrichmentIdsInCell * qPoint +
696
- (jNode - dofsPerCellCFE)] *
697
- cellJxWValuesEnrichmentBlockClassical[qPoint];
698
- }
699
- *basisOverlapTmpIter += NciNpj - NcicjNcj;
779
+ // // Ni_classical using Mc = d quadrature * interpolated
780
+ // // ci's in Ni_classicalQuadrature of Mc = d
781
+ // for (unsigned int qPoint = 0;
782
+ // qPoint < nQuadPointInCellEnrichmentBlockClassical;
783
+ // qPoint++)
784
+ // {
785
+ // NcicjNcj +=
786
+ // *(cumulativeEnrichmentBlockClassicalDofQuadPoints +
787
+ // dofsPerCellCFE * qPoint + iNode
788
+ // /*nQuadPointInCellEnrichmentBlockClassical *
789
+ // iNode + qPoint*/) *
790
+ // classicalComponentInQuadValuesEC
791
+ // [numEnrichmentIdsInCell * qPoint +
792
+ // (jNode - dofsPerCellCFE)] *
793
+ // cellJxWValuesEnrichmentBlockClassical[qPoint];
794
+ // }
795
+ // *basisOverlapTmpIter += NciNpj - NcicjNcj;
796
+
797
+ *basisOverlapTmpIter =
798
+ *(basisOverlapECBlockEnrich.data () + (jNode - dofsPerCellCFE) * dofsPerCell + iNode) -
799
+ *(basisOverlapECBlockClass.data () + (jNode - dofsPerCellCFE) * dofsPerCellCFE + iNode);
700
800
}
701
801
702
802
else if (iNode >= dofsPerCellCFE && jNode >= dofsPerCellCFE)
0 commit comments