36 #include <loki/Singleton.h>
44 namespace Collections {
73 int nPhys = m_stdExp->GetTotPoints();
74 int ntot = m_numElmt*nPhys;
75 int nmodes = m_stdExp->GetNcoeffs();
80 in[0] = entry0; in[1] = entry1;
83 output = (m_coordim == 3)? entry3: (m_coordim == 2)?
86 for(
int i = 0; i < m_dim; ++i)
88 tmp[i] = wsp + i*ntot;
92 for(
int i = 0; i < m_dim; ++i)
96 for(
int j = 1; j < m_coordim; ++j)
99 in[j],1, tmp[i], 1, tmp[i],1);
107 Blas::Dgemm(
'N',
'N', m_iProdWRTStdDBase[0]->GetRows(),
108 m_numElmt,m_iProdWRTStdDBase[0]->GetColumns(),
109 1.0, m_iProdWRTStdDBase[0]->GetRawPtr(),
110 m_iProdWRTStdDBase[0]->GetRows(),
111 tmp[0].
get(), nPhys, 0.0,
112 output.get(), nmodes);
115 for(
int i = 1; i < m_dim; ++i)
118 Blas::Dgemm(
'N',
'N', m_iProdWRTStdDBase[i]->GetRows(),
119 m_numElmt,m_iProdWRTStdDBase[i]->GetColumns(),
120 1.0, m_iProdWRTStdDBase[i]->GetRawPtr(),
121 m_iProdWRTStdDBase[i]->GetRows(),
122 tmp[i].
get(), nPhys, 1.0,
123 output.get(), nmodes);
133 ASSERTL0(
false,
"Not valid for this operator.");
145 vector<StdRegions::StdExpansionSharedPtr> pCollExp,
150 m_dim = PtsKey.size();
151 m_coordim = m_stdExp->GetCoordim();
153 int nqtot = m_stdExp->GetTotPoints();
154 int nmodes = m_stdExp->GetNcoeffs();
158 for(
int i = 0; i < m_dim; ++i)
163 for(
int j = 0; j < nqtot; ++j)
167 m_stdExp->IProductWRTDerivBase(i,tmp,tmp1);
169 &(m_iProdWRTStdDBase[i]->GetPtr())[0]
173 m_derivFac = pGeomData->GetDerivFactors(pCollExp);
174 m_jac = pGeomData->GetJac(pCollExp);
175 m_wspSize = m_dim*nqtot*m_numElmt;
180 OperatorKey IProductWRTDerivBase_StdMat::m_typeArr[] = {
183 IProductWRTDerivBase_StdMat::create,
184 "IProductWRTDerivBase_StdMat_Seg"),
187 IProductWRTDerivBase_StdMat::create,
188 "IProductWRTDerivBase_StdMat_Tri"),
191 IProductWRTDerivBase_StdMat::create,
192 "IProductWRTDerivBase_StdMat_NodalTri"),
195 IProductWRTDerivBase_StdMat::create,
196 "IProductWRTDerivBase_StdMat_Quad"),
199 IProductWRTDerivBase_StdMat::create,
200 "IProductWRTDerivBase_StdMat_Tet"),
203 IProductWRTDerivBase_StdMat::create,
204 "IProductWRTDerivBase_StdMat_NodalTet"),
207 IProductWRTDerivBase_StdMat::create,
208 "IProductWRTDerivBase_StdMat_Pyr"),
211 IProductWRTDerivBase_StdMat::create,
212 "IProductWRTDerivBase_StdMat_Prism"),
215 IProductWRTDerivBase_StdMat::create,
216 "IProductWRTDerivBase_StdMat_NodalPrism"),
219 IProductWRTDerivBase_StdMat::create,
220 "IProductWRTDerivBase_StdMat_Hex")
243 unsigned int nPhys = m_stdExp->GetTotPoints();
244 unsigned int ntot = m_numElmt*nPhys;
245 unsigned int nmodes = m_stdExp->GetNcoeffs();
246 unsigned int nmax = max(ntot,m_numElmt*nmodes);
251 in[0] = entry0; in[1] = entry1; in[2] = entry2;
253 output = (m_coordim == 3)? entry3: (m_coordim == 2)?
256 for(
int i = 0; i < m_dim; ++i)
258 tmp[i] = wsp + i*nmax;
262 for(
int i = 0; i < m_dim; ++i)
266 for(
int j = 1; j < m_coordim; ++j)
269 in[j],1, tmp[i], 1, tmp[i],1);
276 for(
int n = 0; n < m_numElmt; ++n)
278 m_stdExp->IProductWRTDerivBase(0,tmp[0]+n*nPhys,
279 tmp1 = output + n*nmodes);
283 for(
int i = 1; i < m_dim; ++i)
287 for(
int n = 0; n < m_numElmt; ++n)
289 m_stdExp->IProductWRTDerivBase(i,tmp[i]+n*nPhys,tmp[0]);
291 tmp1 = output+n*nmodes,1);
302 ASSERTL0(
false,
"Not valid for this operator.");
313 vector<StdRegions::StdExpansionSharedPtr> pCollExp,
318 m_dim = PtsKey.size();
319 m_coordim = m_stdExp->GetCoordim();
321 int nqtot = m_stdExp->GetTotPoints();
323 m_derivFac = pGeomData->GetDerivFactors(pCollExp);
324 m_jac = pGeomData->GetJac(pCollExp);
325 m_wspSize = m_dim*nqtot*m_numElmt;
330 OperatorKey IProductWRTDerivBase_IterPerExp::m_typeArr[] = {
333 IProductWRTDerivBase_IterPerExp::create,
334 "IProductWRTDerivBase_IterPerExp_Seg"),
337 IProductWRTDerivBase_IterPerExp::create,
338 "IProductWRTDerivBase_IterPerExp_Tri"),
341 IProductWRTDerivBase_IterPerExp::create,
342 "IProductWRTDerivBase_IterPerExp_NodalTri"),
345 IProductWRTDerivBase_IterPerExp::create,
346 "IProductWRTDerivBase_IterPerExp_Quad"),
349 IProductWRTDerivBase_IterPerExp::create,
350 "IProductWRTDerivBase_IterPerExp_Tet"),
353 IProductWRTDerivBase_IterPerExp::create,
354 "IProductWRTDerivBase_IterPerExp_NodalTet"),
357 IProductWRTDerivBase_IterPerExp::create,
358 "IProductWRTDerivBase_IterPerExp_Pyr"),
361 IProductWRTDerivBase_IterPerExp::create,
362 "IProductWRTDerivBase_IterPerExp_Prism"),
365 IProductWRTDerivBase_IterPerExp::create,
366 "IProductWRTDerivBase_IterPerExp_NodalPrism"),
369 IProductWRTDerivBase_IterPerExp::create,
370 "IProductWRTDerivBase_IterPerExp_Hex")
394 unsigned int nmodes = m_expList[0]->GetNcoeffs();
395 unsigned int nPhys = m_expList[0]->GetTotPoints();
400 in[0] = entry0; in[1] = entry1; in[2] = entry2;
402 output = (m_coordim == 3)? entry3: (m_coordim == 2)?
405 for(
int n = 0; n < m_numElmt; ++n)
407 m_expList[n]->IProductWRTDerivBase(0,
409 tmp1 = output + n * nmodes);
412 for(
int i = 1; i < m_dim; ++i)
414 for(
int n = 0; n < m_numElmt; ++n)
416 m_expList[n]->IProductWRTDerivBase(i,in[i]+n*nPhys,tmp);
419 tmp1 = output+n*nmodes,1);
430 ASSERTL0(
false,
"Not valid for this operator.");
440 vector<StdRegions::StdExpansionSharedPtr> pCollExp,
444 m_expList = pCollExp;
445 m_dim = pCollExp[0]->GetNumBases();
446 m_coordim = pCollExp[0]->GetCoordim();
451 OperatorKey IProductWRTDerivBase_NoCollection::m_typeArr[] = {
454 IProductWRTDerivBase_NoCollection::create,
455 "IProductWRTDerivBase_NoCollection_Seg"),
458 IProductWRTDerivBase_NoCollection::create,
459 "IProductWRTDerivBase_NoCollection_Tri"),
462 IProductWRTDerivBase_NoCollection::create,
463 "IProductWRTDerivBase_NoCollection_NodalTri"),
466 IProductWRTDerivBase_NoCollection::create,
467 "IProductWRTDerivBase_NoCollection_Quad"),
470 IProductWRTDerivBase_NoCollection::create,
471 "IProductWRTDerivBase_NoCollection_Tet"),
474 IProductWRTDerivBase_NoCollection::create,
475 "IProductWRTDerivBase_NoCollection_NodalTet"),
478 IProductWRTDerivBase_NoCollection::create,
479 "IProductWRTDerivBase_NoCollection_Pyr"),
482 IProductWRTDerivBase_NoCollection::create,
483 "IProductWRTDerivBase_NoCollection_Prism"),
486 IProductWRTDerivBase_NoCollection::create,
487 "IProductWRTDerivBase_NoCollection_NodalPrism"),
490 IProductWRTDerivBase_NoCollection::create,
491 "IProductWRTDerivBase_NoCollection_Hex")
516 Vmath::Vmul(m_numElmt*m_nquad0, m_jac, 1, input, 1, wsp, 1);
517 Vmath::Vmul(m_numElmt*m_nquad0, &m_derivFac[0][0], 1,
522 Blas::Dgemm(
'T',
'N', m_nmodes0, m_numElmt, m_nquad0,
523 1.0, m_derbase0.get(), m_nquad0,
524 &wsp[0], m_nquad0, 0.0,
525 &output[0], m_nmodes0);
534 ASSERTL0(
false,
"Not valid for this operator.");
546 vector<StdRegions::StdExpansionSharedPtr> pCollExp,
549 m_nquad0 (m_stdExp->GetNumPoints(0)),
550 m_nmodes0 (m_stdExp->GetBasisNumModes(0)),
551 m_derbase0(m_stdExp->GetBasis(0)->GetDbdata())
553 m_wspSize = m_numElmt*m_nquad0;
554 m_derivFac = pGeomData->GetDerivFactors(pCollExp);
555 m_jac = pGeomData->GetJacWithStdWeights(pCollExp);
561 RegisterCreatorFunction(
563 IProductWRTDerivBase_SumFac_Seg::create,
564 "IProductWRTDerivBase_SumFac_Seg");
586 unsigned int nPhys = m_stdExp->GetTotPoints();
587 unsigned int ntot = m_numElmt*nPhys;
588 unsigned int nmodes = m_stdExp->GetNcoeffs();
589 unsigned int nmax = max(ntot,m_numElmt*nmodes);
594 in[0] = entry0; in[1] = entry1; in[2] = entry2;
596 output = (m_coordim == 2)? entry2: entry3;
598 tmp[0] = wsp; tmp[1] = wsp + nmax;
602 for(
int i = 0; i < 2; ++i)
606 for(
int j = 1; j < m_coordim; ++j)
609 in[j],1, tmp[i], 1, tmp[i],1);
616 m_nmodes0, m_nmodes1,
618 m_jac, tmp[0], output, wsp1);
623 m_nmodes0, m_nmodes1,
625 m_jac, tmp[1], tmp[0], wsp1);
627 Vmath::Vadd(m_numElmt*nmodes,tmp[0],1,output,1,output,1);
636 ASSERTL0(
false,
"Not valid for this operator.");
656 vector<StdRegions::StdExpansionSharedPtr> pCollExp,
659 m_nquad0 (m_stdExp->GetNumPoints(0)),
660 m_nquad1 (m_stdExp->GetNumPoints(1)),
661 m_nmodes0 (m_stdExp->GetBasisNumModes(0)),
662 m_nmodes1 (m_stdExp->GetBasisNumModes(1)),
663 m_colldir0(m_stdExp->GetBasis(0)->Collocation()),
664 m_colldir1(m_stdExp->GetBasis(1)->Collocation()),
665 m_base0 (m_stdExp->GetBasis(0)->GetBdata()),
666 m_base1 (m_stdExp->GetBasis(1)->GetBdata()),
667 m_derbase0(m_stdExp->GetBasis(0)->GetDbdata()),
668 m_derbase1(m_stdExp->GetBasis(1)->GetDbdata())
671 m_coordim = m_stdExp->GetCoordim();
673 m_derivFac = pGeomData->GetDerivFactors(pCollExp);
674 m_jac = pGeomData->GetJacWithStdWeights(pCollExp);
675 m_wspSize = 4 * m_numElmt * (max(m_nquad0*m_nquad1,
676 m_nmodes0*m_nmodes1));
681 OperatorKey IProductWRTDerivBase_SumFac_Quad::m_type =
684 IProductWRTDerivBase_SumFac_Quad::create,
685 "IProductWRTDerivBase_IterPerExp_Quad");
707 unsigned int nPhys = m_stdExp->GetTotPoints();
708 unsigned int ntot = m_numElmt*nPhys;
709 unsigned int nmodes = m_stdExp->GetNcoeffs();
710 unsigned int nmax = max(ntot,m_numElmt*nmodes);
715 in[0] = entry0; in[1] = entry1; in[2] = entry2;
717 output = (m_coordim == 2)? entry2: entry3;
719 tmp[0] = wsp; tmp[1] = wsp + nmax;
742 for(
int i = 0; i < 2; ++i)
744 Vmath::Vmul (ntot,m_derivFac[i],1, in[0],1, tmp[i],1);
746 for(
int j = 1; j < m_coordim; ++j)
749 in[j],1, tmp[i], 1, tmp[i],1);
754 for (
int i = 0; i < m_numElmt; ++i)
757 Vmath::Vmul(nPhys,&m_fac0[0],1,tmp[0].
get()+i*nPhys,1,
758 tmp[0].
get()+i*nPhys,1);
762 tmp[0].
get()+i*nPhys,1,tmp[0].
get()+i*nPhys,1);
766 TriIProduct(m_sortTopVertex, m_numElmt, m_nquad0, m_nquad1,
767 m_nmodes0, m_nmodes1, m_derbase0, m_base1,
768 m_jac, tmp[0], output, wsp1);
771 TriIProduct(m_sortTopVertex, m_numElmt, m_nquad0, m_nquad1,
772 m_nmodes0, m_nmodes1, m_base0, m_derbase1,
773 m_jac, tmp[1], tmp[0], wsp1);
775 Vmath::Vadd(m_numElmt*nmodes,tmp[0],1,output,1,output,1);
784 ASSERTL0(
false,
"Not valid for this operator.");
807 vector<StdRegions::StdExpansionSharedPtr> pCollExp,
810 m_nquad0 (m_stdExp->GetNumPoints(0)),
811 m_nquad1 (m_stdExp->GetNumPoints(1)),
812 m_nmodes0 (m_stdExp->GetBasisNumModes(0)),
813 m_nmodes1 (m_stdExp->GetBasisNumModes(1)),
814 m_colldir0(m_stdExp->GetBasis(0)->Collocation()),
815 m_colldir1(m_stdExp->GetBasis(1)->Collocation()),
816 m_base0 (m_stdExp->GetBasis(0)->GetBdata()),
817 m_base1 (m_stdExp->GetBasis(1)->GetBdata()),
818 m_derbase0(m_stdExp->GetBasis(0)->GetDbdata()),
819 m_derbase1(m_stdExp->GetBasis(1)->GetDbdata())
822 m_coordim = m_stdExp->GetCoordim();
824 m_derivFac = pGeomData->GetDerivFactors(pCollExp);
825 m_jac = pGeomData->GetJacWithStdWeights(pCollExp);
826 m_wspSize = 4 * m_numElmt * (max(m_nquad0*m_nquad1,
827 m_nmodes0*m_nmodes1));
829 if(m_stdExp->GetBasis(0)->GetBasisType()
832 m_sortTopVertex =
true;
836 m_sortTopVertex =
false;
840 = m_stdExp->GetBasis(0)->GetZ();
842 = m_stdExp->GetBasis(1)->GetZ();
846 for (
int i = 0; i < m_nquad0; ++i)
848 for(
int j = 0; j < m_nquad1; ++j)
850 m_fac0[i+j*m_nquad0] = 2.0/(1-z1[j]);
856 for (
int i = 0; i < m_nquad0; ++i)
858 for(
int j = 0; j < m_nquad1; ++j)
860 m_fac1[i+j*m_nquad0] = (1+z0[i])/(1-z1[j]);
867 OperatorKey IProductWRTDerivBase_SumFac_Tri::m_type =
870 IProductWRTDerivBase_SumFac_Tri::create,
871 "IProductWRTDerivBase_IterPerExp_Tri");
893 unsigned int nPhys = m_stdExp->GetTotPoints();
894 unsigned int ntot = m_numElmt*nPhys;
895 unsigned int nmodes = m_stdExp->GetNcoeffs();
896 unsigned int nmax = max(ntot,m_numElmt*nmodes);
901 in[0] = entry0; in[1] = entry1;
906 for(
int i = 0; i < 3; ++i)
908 tmp[i] = wsp + i*nmax;
912 for(
int i = 0; i < 3; ++i)
916 for(
int j = 1; j < 3; ++j)
919 in[j],1, tmp[i], 1, tmp[i],1);
926 HexIProduct(
false,m_colldir1,m_colldir2, m_numElmt,
927 m_nquad0, m_nquad1, m_nquad2,
928 m_nmodes0, m_nmodes1, m_nmodes2,
929 m_derbase0, m_base1, m_base2,
930 m_jac,tmp[0],output,wsp1);
932 HexIProduct(m_colldir0,
false,m_colldir2, m_numElmt,
933 m_nquad0, m_nquad1, m_nquad2,
934 m_nmodes0, m_nmodes1, m_nmodes2,
935 m_base0, m_derbase1, m_base2,
936 m_jac,tmp[1],tmp[0],wsp1);
937 Vmath::Vadd(m_numElmt*nmodes,tmp[0],1,output,1,output,1);
939 HexIProduct(m_colldir0,m_colldir1,
false, m_numElmt,
940 m_nquad0, m_nquad1, m_nquad2,
941 m_nmodes0, m_nmodes1, m_nmodes2,
942 m_base0, m_base1, m_derbase2,
943 m_jac,tmp[2],tmp[0],wsp1);
944 Vmath::Vadd(m_numElmt*nmodes,tmp[0],1,output,1,output,1);
953 ASSERTL0(
false,
"Not valid for this operator.");
977 vector<StdRegions::StdExpansionSharedPtr> pCollExp,
980 m_nquad0 (m_stdExp->GetNumPoints(0)),
981 m_nquad1 (m_stdExp->GetNumPoints(1)),
982 m_nquad2 (m_stdExp->GetNumPoints(2)),
983 m_nmodes0 (m_stdExp->GetBasisNumModes(0)),
984 m_nmodes1 (m_stdExp->GetBasisNumModes(1)),
985 m_nmodes2 (m_stdExp->GetBasisNumModes(2)),
986 m_colldir0(m_stdExp->GetBasis(0)->Collocation()),
987 m_colldir1(m_stdExp->GetBasis(1)->Collocation()),
988 m_colldir2(m_stdExp->GetBasis(2)->Collocation()),
989 m_base0 (m_stdExp->GetBasis(0)->GetBdata()),
990 m_base1 (m_stdExp->GetBasis(1)->GetBdata()),
991 m_base2 (m_stdExp->GetBasis(2)->GetBdata()),
992 m_derbase0(m_stdExp->GetBasis(0)->GetDbdata()),
993 m_derbase1(m_stdExp->GetBasis(1)->GetDbdata()),
994 m_derbase2(m_stdExp->GetBasis(2)->GetDbdata())
997 m_jac = pGeomData->GetJacWithStdWeights(pCollExp);
998 m_wspSize = 6 * m_numElmt * (max(m_nquad0*m_nquad1*m_nquad2,
999 m_nmodes0*m_nmodes1*m_nmodes2));
1000 m_derivFac = pGeomData->GetDerivFactors(pCollExp);
1006 RegisterCreatorFunction(
1008 IProductWRTDerivBase_SumFac_Hex::create,
1009 "IProductWRTDerivBase_SumFac_Hex");
1020 virtual
void operator()(
1027 unsigned int nPhys = m_stdExp->GetTotPoints();
1028 unsigned int ntot = m_numElmt*nPhys;
1029 unsigned int nmodes = m_stdExp->GetNcoeffs();
1030 unsigned int nmax = max(ntot,m_numElmt*nmodes);
1035 in[0] = entry0; in[1] = entry1;
1040 for(
int i = 0; i < 3; ++i)
1042 tmp[i] = wsp + i*nmax;
1083 for(
int i = 0; i < 3; ++i)
1085 Vmath::Vmul (ntot,m_derivFac[i],1, in[0],1,tmp[i],1);
1086 for(
int j = 1; j < 3; ++j)
1089 in[j],1, tmp[i], 1, tmp[i],1);
1093 wsp1 = wsp + 3*nmax;
1096 for (
int i = 0; i < m_numElmt; ++i)
1100 tmp[2].
get() + i*nPhys, 1,
1105 tmp[0].get()+i*nPhys,1,tmp[0].get()+i*nPhys,1);
1108 Vmath::Vmul(nPhys,&m_fac0[0],1,tmp[0].
get()+i*nPhys,1,
1109 tmp[0].
get()+i*nPhys,1);
1112 Vmath::Vvtvp(nPhys,&m_fac3[0],1,tmp[2].
get()+i*nPhys,1,
1113 tmp[1].
get()+i*nPhys,1,tmp[1].
get()+i*nPhys,1);
1116 Vmath::Vmul(nPhys,&m_fac2[0],1,tmp[1].
get()+i*nPhys,1,
1117 tmp[1].
get()+i*nPhys,1);
1123 m_nquad0, m_nquad1, m_nquad2,
1124 m_nmodes0, m_nmodes1, m_nmodes2,
1125 m_derbase0, m_base1, m_base2,
1126 m_jac,tmp[0],output,wsp1);
1129 m_nquad0, m_nquad1, m_nquad2,
1130 m_nmodes0, m_nmodes1, m_nmodes2,
1131 m_base0, m_derbase1, m_base2,
1132 m_jac,tmp[1],tmp[0],wsp1);
1133 Vmath::Vadd(m_numElmt*nmodes,tmp[0],1,output,1,output,1);
1136 m_nquad0, m_nquad1, m_nquad2,
1137 m_nmodes0, m_nmodes1, m_nmodes2,
1138 m_base0, m_base1, m_derbase2,
1139 m_jac,tmp[2],tmp[0],wsp1);
1140 Vmath::Vadd(m_numElmt*nmodes,tmp[0],1,output,1,output,1);
1149 ASSERTL0(
false,
"Not valid for this operator.");
1175 vector<StdRegions::StdExpansionSharedPtr> pCollExp,
1178 m_nquad0 (m_stdExp->GetNumPoints(0)),
1179 m_nquad1 (m_stdExp->GetNumPoints(1)),
1180 m_nquad2 (m_stdExp->GetNumPoints(2)),
1181 m_nmodes0 (m_stdExp->GetBasisNumModes(0)),
1182 m_nmodes1 (m_stdExp->GetBasisNumModes(1)),
1183 m_nmodes2 (m_stdExp->GetBasisNumModes(2)),
1184 m_base0 (m_stdExp->GetBasis(0)->GetBdata()),
1185 m_base1 (m_stdExp->GetBasis(1)->GetBdata()),
1186 m_base2 (m_stdExp->GetBasis(2)->GetBdata()),
1187 m_derbase0(m_stdExp->GetBasis(0)->GetDbdata()),
1188 m_derbase1(m_stdExp->GetBasis(1)->GetDbdata()),
1189 m_derbase2(m_stdExp->GetBasis(2)->GetDbdata())
1192 m_jac = pGeomData->GetJacWithStdWeights(pCollExp);
1193 m_wspSize = 6 * m_numElmt * (max(m_nquad0*m_nquad1*m_nquad2,
1194 m_nmodes0*m_nmodes1*m_nmodes2));
1195 m_derivFac = pGeomData->GetDerivFactors(pCollExp);
1199 = m_stdExp->GetBasis(0)->GetZ();
1201 = m_stdExp->GetBasis(1)->GetZ();
1203 = m_stdExp->GetBasis(2)->GetZ();
1210 for (
int i = 0; i < m_nquad0; ++i)
1212 for(
int j = 0; j < m_nquad1; ++j)
1214 for(
int k = 0; k < m_nquad2; ++k)
1217 m_fac0[i + j*m_nquad0 + k*m_nquad0*m_nquad1]
1218 = 4.0/((1-z1[j])*(1-z2[k]));
1219 m_fac1[i + j*m_nquad0 + k*m_nquad0*m_nquad1]
1221 m_fac2[i + j*m_nquad0 + k*m_nquad0*m_nquad1]
1223 m_fac3[i + j*m_nquad0 + k*m_nquad0*m_nquad1]
1229 if(m_stdExp->GetBasis(0)->GetBasisType()
1232 m_sortTopEdge =
true;
1236 m_sortTopEdge =
false;
1244 RegisterCreatorFunction(
1246 IProductWRTDerivBase_SumFac_Tet::create,
1247 "IProductWRTDerivBase_SumFac_Tet");
1269 unsigned int nPhys = m_stdExp->GetTotPoints();
1270 unsigned int ntot = m_numElmt*nPhys;
1271 unsigned int nmodes = m_stdExp->GetNcoeffs();
1272 unsigned int nmax = max(ntot,m_numElmt*nmodes);
1277 in[0] = entry0; in[1] = entry1;
1282 for(
int i = 0; i < 3; ++i)
1284 tmp[i] = wsp + i*nmax;
1313 for(
int i = 0; i < 3; ++i)
1317 for(
int j = 1; j < 3; ++j)
1320 in[j],1, tmp[i], 1, tmp[i],1);
1323 wsp1 = wsp + 3*nmax;
1326 for (
int i = 0; i < m_numElmt; ++i)
1329 Vmath::Vmul(nPhys,&m_fac0[0],1,tmp[0].
get()+i*nPhys,1,
1330 tmp[0].
get()+i*nPhys,1);
1333 Vmath::Vvtvp(nPhys,&m_fac1[0],1,tmp[2].
get()+i*nPhys,1,
1334 tmp[0].
get()+i*nPhys,1,tmp[0].
get()+i*nPhys,1);
1339 m_nquad0, m_nquad1, m_nquad2,
1340 m_nmodes0, m_nmodes1, m_nmodes2,
1341 m_derbase0, m_base1, m_base2,
1342 m_jac,tmp[0],output,wsp1);
1345 m_nquad0, m_nquad1, m_nquad2,
1346 m_nmodes0, m_nmodes1, m_nmodes2,
1347 m_base0, m_derbase1, m_base2,
1348 m_jac,tmp[1],tmp[0],wsp1);
1349 Vmath::Vadd(m_numElmt*nmodes,tmp[0],1,output,1,output,1);
1352 m_nquad0, m_nquad1, m_nquad2,
1353 m_nmodes0, m_nmodes1, m_nmodes2,
1354 m_base0, m_base1, m_derbase2,
1355 m_jac,tmp[2],tmp[0],wsp1);
1356 Vmath::Vadd(m_numElmt*nmodes,tmp[0],1,output,1,output,1);
1365 ASSERTL0(
false,
"Not valid for this operator.");
1390 vector<StdRegions::StdExpansionSharedPtr> pCollExp,
1393 m_nquad0 (m_stdExp->GetNumPoints(0)),
1394 m_nquad1 (m_stdExp->GetNumPoints(1)),
1395 m_nquad2 (m_stdExp->GetNumPoints(2)),
1396 m_nmodes0 (m_stdExp->GetBasisNumModes(0)),
1397 m_nmodes1 (m_stdExp->GetBasisNumModes(1)),
1398 m_nmodes2 (m_stdExp->GetBasisNumModes(2)),
1399 m_base0 (m_stdExp->GetBasis(0)->GetBdata()),
1400 m_base1 (m_stdExp->GetBasis(1)->GetBdata()),
1401 m_base2 (m_stdExp->GetBasis(2)->GetBdata()),
1402 m_derbase0(m_stdExp->GetBasis(0)->GetDbdata()),
1403 m_derbase1(m_stdExp->GetBasis(1)->GetDbdata()),
1404 m_derbase2(m_stdExp->GetBasis(2)->GetDbdata())
1407 m_jac = pGeomData->GetJacWithStdWeights(pCollExp);
1408 m_wspSize = 6 * m_numElmt * (max(m_nquad0*m_nquad1*m_nquad2,
1409 m_nmodes0*m_nmodes1*m_nmodes2));
1410 m_derivFac = pGeomData->GetDerivFactors(pCollExp);
1412 if(m_stdExp->GetBasis(0)->GetBasisType()
1415 m_sortTopVertex =
true;
1419 m_sortTopVertex =
false;
1423 = m_stdExp->GetBasis(0)->GetZ();
1425 = m_stdExp->GetBasis(2)->GetZ();
1430 for (
int i = 0; i < m_nquad0; ++i)
1432 for(
int j = 0; j < m_nquad1; ++j)
1434 for(
int k = 0; k < m_nquad2; ++k)
1437 m_fac0[i + j*m_nquad0 + k*m_nquad0*m_nquad1]
1440 m_fac1[i + j*m_nquad0 + k*m_nquad0*m_nquad1]
1441 = (1+z0[i])/(1-z2[k]);
1451 RegisterCreatorFunction(
1453 IProductWRTDerivBase_SumFac_Prism::create,
1454 "IProductWRTDerivBase_SumFac_Prism");
virtual void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
Inner product WRT deriv base operator using sum-factorisation (Tri)
#define ASSERTL0(condition, msg)
Array< TwoD, const NekDouble > m_derivFac
Inner product WRT deriv base operator using sum-factorisation (Segment)
void HexIProduct(bool colldir0, bool colldir1, bool colldir2, int numElmt, int nquad0, int nquad1, int nquad2, int nmodes0, int nmodes1, int nmodes2, const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &base2, const Array< OneD, const NekDouble > &jac, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
IProductWRTDerivBase_StdMat(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
std::vector< PointsKey > PointsKeyVector
Array< OneD, NekDouble > m_fac0
boost::tuple< LibUtilities::ShapeType, OperatorType, ImplementationType, ExpansionIsNodal > OperatorKey
Key for describing an Operator.
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
Array< OneD, NekDouble > m_fac0
Array< TwoD, const NekDouble > m_derivFac
Array< OneD, const NekDouble > m_derbase1
Array< OneD, NekDouble > m_fac1
IProductWRTDerivBase_SumFac_Tri(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Principle Modified Functions .
Inner product WRT deriv base operator using LocalRegions implementation.
Inner product WRT deriv base operator using standard matrix approach.
Array< OneD, const NekDouble > m_base1
Array< TwoD, const NekDouble > m_derivFac
Array< OneD, const NekDouble > m_base0
Array< OneD, const NekDouble > m_derbase1
Array< OneD, const NekDouble > m_base1
Array< OneD, const NekDouble > m_jac
Array< OneD, const NekDouble > m_jac
IProductWRTDerivBase_SumFac_Quad(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
void QuadIProduct(bool colldir0, bool colldir1, int numElmt, int nquad0, int nquad1, int nmodes0, int nmodes1, const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &jac, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
Base class for operators on a collection of elements.
Array< OneD, const NekDouble > m_base2
virtual void operator()(const Array< OneD, const NekDouble > &entry0, Array< OneD, NekDouble > &entry1, Array< OneD, NekDouble > &entry2, Array< OneD, NekDouble > &entry3, Array< OneD, NekDouble > &wsp)
Perform operation.
virtual void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
virtual void operator()(const Array< OneD, const NekDouble > &entry0, Array< OneD, NekDouble > &entry1, Array< OneD, NekDouble > &entry2, Array< OneD, NekDouble > &entry3, Array< OneD, NekDouble > &wsp)
Perform operation.
virtual void operator()(const Array< OneD, const NekDouble > &entry0, Array< OneD, NekDouble > &entry1, Array< OneD, NekDouble > &entry2, Array< OneD, NekDouble > &entry3, Array< OneD, NekDouble > &wsp)
Perform operation.
virtual void operator()(const Array< OneD, const NekDouble > &entry0, Array< OneD, NekDouble > &entry1, Array< OneD, NekDouble > &entry2, Array< OneD, NekDouble > &entry3, Array< OneD, NekDouble > &wsp)
Perform operation.
void TetIProduct(bool sortTopEdge, int numElmt, int nquad0, int nquad1, int nquad2, int nmodes0, int nmodes1, int nmodes2, const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &base2, const Array< OneD, const NekDouble > &jac, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
virtual void operator()(const Array< OneD, const NekDouble > &entry0, Array< OneD, NekDouble > &entry1, Array< OneD, NekDouble > &entry2, Array< OneD, NekDouble > &entry3, Array< OneD, NekDouble > &wsp)
Perform operation.
Array< OneD, const NekDouble > m_jac
IProductWRTDerivBase_SumFac_Prism(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
Array< OneD, const NekDouble > m_jac
Array< OneD, NekDouble > m_fac2
void TriIProduct(bool sortTopVertex, int numElmt, int nquad0, int nquad1, int nmodes0, int nmodes1, const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &jac, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
Array< TwoD, const NekDouble > m_derivFac
virtual void operator()(const Array< OneD, const NekDouble > &entry0, Array< OneD, NekDouble > &entry1, Array< OneD, NekDouble > &entry2, Array< OneD, NekDouble > &entry3, Array< OneD, NekDouble > &wsp)
Perform operation.
Inner product WRT deriv base operator using sum-factorisation (Tet)
Array< OneD, const NekDouble > m_base2
virtual void operator()(const Array< OneD, const NekDouble > &entry0, Array< OneD, NekDouble > &entry1, Array< OneD, NekDouble > &entry2, Array< OneD, NekDouble > &entry3, Array< OneD, NekDouble > &wsp)
Perform operation.
virtual void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
Array< OneD, const NekDouble > m_derbase0
Array< OneD, const NekDouble > m_derbase0
virtual void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
virtual void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
IProductWRTDerivBase_SumFac_Seg(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
Array< OneD, const NekDouble > m_jac
virtual void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
Array< OneD, const NekDouble > m_base0
IProductWRTDerivBase_NoCollection(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
Inner product WRT deriv base operator using element-wise operation.
Inner product WRT deriv base operator using sum-factorisation (Quad)
IProductWRTDerivBase_IterPerExp(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
Inner product WRT deriv base operator using sum-factorisation (Hex)
OperatorFactory & GetOperatorFactory()
Returns the singleton Operator factory object.
virtual void operator()(const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &output1, Array< OneD, NekDouble > &output2, Array< OneD, NekDouble > &wsp)
Perform operation.
Array< OneD, const NekDouble > m_derbase0
Array< OneD, NekDouble > m_fac0
Array< OneD, const NekDouble > m_derbase2
void PrismIProduct(bool sortTopVertex, int numElmt, int nquad0, int nquad1, int nquad2, int nmodes0, int nmodes1, int nmodes2, const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &base2, const Array< OneD, const NekDouble > &jac, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
Array< OneD, const NekDouble > m_derbase1
Array< OneD, const NekDouble > m_base1
Array< OneD, NekDouble > m_fac1
Inner product WRT deriv base operator using sum-factorisation (Prism)
#define OPERATOR_CREATE(cname)
Array< OneD, NekDouble > m_fac1
IProductWRTDerivBase_SumFac_Hex(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
Array< OneD, const NekDouble > m_base2
Array< OneD, NekDouble > m_fac3
virtual void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
Array< OneD, const NekDouble > m_base0
Array< OneD, const NekDouble > m_derbase2
Array< OneD, DNekMatSharedPtr > m_iProdWRTStdDBase
Array< OneD, const NekDouble > m_derbase1
virtual void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
Array< TwoD, const NekDouble > m_derivFac
boost::shared_ptr< CoalescedGeomData > CoalescedGeomDataSharedPtr
virtual void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
Array< OneD, const NekDouble > m_derbase1
Array< OneD, const NekDouble > m_jac
Array< OneD, const NekDouble > m_derbase0
Array< TwoD, const NekDouble > m_derivFac
Array< OneD, const NekDouble > m_derbase0
Array< OneD, const NekDouble > m_base1
Array< TwoD, const NekDouble > m_derivFac
Array< OneD, const NekDouble > m_jac
Array< OneD, const NekDouble > m_base1
void Zero(int n, T *x, const int incx)
Zero vector.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Array< OneD, const NekDouble > m_jac
IProductWRTDerivBase_SumFac_Tet(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
Array< TwoD, const NekDouble > m_derivFac
vector< StdRegions::StdExpansionSharedPtr > m_expList
Array< OneD, const NekDouble > m_derbase2
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Array< OneD, const NekDouble > m_base0
Array< OneD, const NekDouble > m_base0
Array< OneD, const NekDouble > m_derbase0
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.