48 namespace MultiRegions
53 string PreconditionerLowEnergy::className
56 PreconditionerLowEnergy::create,
57 "LowEnergy Preconditioning");
66 PreconditionerLowEnergy::PreconditionerLowEnergy(
67 const boost::shared_ptr<GlobalLinSys> &plinsys,
71 m_locToGloMap(pLocToGloMap)
81 boost::shared_ptr<MultiRegions::ExpList>
82 expList=((
m_linsys.lock())->GetLocMat()).lock();
86 locExpansion = expList->GetExp(0);
88 int nDim = locExpansion->GetShapeDimension();
91 "Preconditioner type only valid in 3D");
123 boost::shared_ptr<MultiRegions::ExpList>
124 expList=((
m_linsys.lock())->GetLocMat()).lock();
129 int nVerts, nEdges,nFaces;
130 int eid, fid, n, cnt, nedgemodes, nfacemodes;
133 int vMap1, vMap2, sign1, sign2;
134 int m, v, eMap1, eMap2, fMap1, fMap2;
135 int offset, globalrow, globalcol, nCoeffs;
141 expList->GetPeriodicEntities(periodicVerts,periodicEdges,periodicFaces);
172 map<LibUtilities::ShapeType,DNekScalMatSharedPtr> transmatrixmap;
173 map<LibUtilities::ShapeType,DNekScalMatSharedPtr> transposedtransmatrixmap;
185 int n_exp = expList->GetNumElmts();
190 int numBlks = 1+nNonDirEdgeIDs+nNonDirFaceIDs;
192 for(i = 0; i < numBlks; ++i)
196 n_blks[0]=nNonDirVerts;
200 map<int,int> uniqueEdgeMap;
201 map<int,int> uniqueFaceMap;
220 for(i=0; i<extradiredges.num_elements(); ++i)
222 meshEdgeId=extradiredges[i];
223 edgeDirMap.insert(meshEdgeId);
227 for(i = 0; i < bndCondExp.num_elements(); i++)
230 for(j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
232 bndCondFaceExp = boost::dynamic_pointer_cast<
234 if (bndConditions[i]->GetBoundaryConditionType() ==
237 for(k = 0; k < bndCondFaceExp->GetNedges(); k++)
240 if(edgeDirMap.count(meshEdgeId) == 0)
242 edgeDirMap.insert(meshEdgeId);
246 faceDirMap.insert(meshFaceId);
254 int nlocalNonDirEdges=0;
255 int nlocalNonDirFaces=0;
257 int edgematrixlocation=0;
258 int ntotaledgeentries=0;
260 map<int,int> EdgeSize;
261 map<int,int> FaceSize;
264 for(n = 0; n < n_exp; ++n)
266 eid = expList->GetOffset_Elmt_Id(n);
267 locExpansion = expList->GetExp(eid);
269 nEdges = locExpansion->GetNedges();
270 for(j = 0; j < nEdges; ++j)
272 int nEdgeInteriorCoeffs = locExpansion->GetEdgeNcoeffs(j) - 2;
274 EdgeSize[meshEdgeId] = nEdgeInteriorCoeffs;
277 nFaces = locExpansion->GetNfaces();
278 for(j = 0; j < nFaces; ++j)
280 int nFaceInteriorCoeffs = locExpansion->GetFaceIntNcoeffs(j);
282 FaceSize[meshFaceId] = nFaceInteriorCoeffs;
286 m_comm = expList->GetComm();
292 PeriodicMap::const_iterator pIt;
293 for (pIt = periodicEdges.begin(); pIt != periodicEdges.end(); ++pIt)
295 meshEdgeId = pIt->first;
297 if(edgeDirMap.count(meshEdgeId)==0)
299 dof = EdgeSize[meshEdgeId];
300 if(uniqueEdgeMap.count(meshEdgeId)==0 && dof > 0)
302 bool SetUpNewEdge =
true;
305 for (i = 0; i < pIt->second.size(); ++i)
307 if (!pIt->second[i].isLocal)
312 int meshEdgeId2 = pIt->second[i].id;
314 if(edgeDirMap.count(meshEdgeId2)==0)
316 if(uniqueEdgeMap.count(meshEdgeId2)!=0)
319 uniqueEdgeMap[meshEdgeId] =
320 uniqueEdgeMap[meshEdgeId2];
321 SetUpNewEdge =
false;
326 edgeDirMap.insert(meshEdgeId);
327 SetUpNewEdge =
false;
333 uniqueEdgeMap[meshEdgeId]=edgematrixlocation;
335 edgeglobaloffset[edgematrixlocation]+=ntotaledgeentries;
337 edgemodeoffset[edgematrixlocation]=dof*dof;
339 ntotaledgeentries+=dof*dof;
341 n_blks[1+edgematrixlocation++]=dof;
349 for(cnt=n=0; n < n_exp; ++n)
351 eid = expList->GetOffset_Elmt_Id(n);
352 locExpansion = expList->GetExp(eid);
354 for (j = 0; j < locExpansion->GetNedges(); ++j)
357 dof = EdgeSize[meshEdgeId];
358 maxEdgeDof = (dof > maxEdgeDof ? dof : maxEdgeDof);
360 if(edgeDirMap.count(meshEdgeId)==0)
362 if(uniqueEdgeMap.count(meshEdgeId)==0 && dof > 0)
365 uniqueEdgeMap[meshEdgeId]=edgematrixlocation;
367 edgeglobaloffset[edgematrixlocation]+=ntotaledgeentries;
369 edgemodeoffset[edgematrixlocation]=dof*dof;
371 ntotaledgeentries+=dof*dof;
373 n_blks[1+edgematrixlocation++]=dof;
376 nlocalNonDirEdges+=dof*dof;
381 int facematrixlocation=0;
382 int ntotalfaceentries=0;
387 for (pIt = periodicFaces.begin(); pIt != periodicFaces.end(); ++pIt)
389 meshFaceId = pIt->first;
391 if(faceDirMap.count(meshFaceId)==0)
393 dof = FaceSize[meshFaceId];
395 if(uniqueFaceMap.count(meshFaceId) == 0 && dof > 0)
397 bool SetUpNewFace =
true;
399 if(pIt->second[0].isLocal)
401 int meshFaceId2 = pIt->second[0].id;
403 if(faceDirMap.count(meshFaceId2)==0)
405 if(uniqueFaceMap.count(meshFaceId2)!=0)
408 uniqueFaceMap[meshFaceId] =
409 uniqueFaceMap[meshFaceId2];
410 SetUpNewFace =
false;
415 faceDirMap.insert(meshFaceId);
416 SetUpNewFace =
false;
422 uniqueFaceMap[meshFaceId]=facematrixlocation;
424 facemodeoffset[facematrixlocation]=dof*dof;
426 faceglobaloffset[facematrixlocation]+=ntotalfaceentries;
428 ntotalfaceentries+=dof*dof;
430 n_blks[1+nNonDirEdgeIDs+facematrixlocation++]=dof;
437 for(cnt=n=0; n < n_exp; ++n)
439 eid = expList->GetOffset_Elmt_Id(n);
441 locExpansion = expList->GetExp(eid);
443 for (j = 0; j < locExpansion->GetNfaces(); ++j)
447 dof = FaceSize[meshFaceId];
448 maxFaceDof = (dof > maxFaceDof ? dof : maxFaceDof);
451 if(faceDirMap.count(meshFaceId)==0)
453 if(uniqueFaceMap.count(meshFaceId)==0 && dof > 0)
455 uniqueFaceMap[meshFaceId]=facematrixlocation;
457 facemodeoffset[facematrixlocation]=dof*dof;
459 faceglobaloffset[facematrixlocation]+=ntotalfaceentries;
461 ntotalfaceentries+=dof*dof;
463 n_blks[1+nNonDirEdgeIDs+facematrixlocation++]=dof;
466 nlocalNonDirFaces+=dof*dof;
485 int edgematrixoffset=0;
486 int facematrixoffset=0;
489 for(n=0; n < n_exp; ++n)
491 eid = expList->GetOffset_Elmt_Id(n);
493 locExpansion = expList->GetExp(eid);
496 for(j = 0; j < locExpansion->GetNedges(); ++j)
503 if(edgeDirMap.count(meshEdgeId)==0)
506 int edgeOffset = edgeglobaloffset[uniqueEdgeMap[meshEdgeId]];
509 int uniOffset = meshEdgeId;
510 pIt = periodicEdges.find(meshEdgeId);
511 if (pIt != periodicEdges.end())
513 for (
int l = 0; l < pIt->second.size(); ++l)
515 uniOffset = min(uniOffset, pIt->second[l].id);
518 uniOffset = uniOffset *maxEdgeDof*maxEdgeDof;
520 for(k=0; k<nedgemodes*nedgemodes; ++k)
522 vGlobal=edgeOffset+k;
523 localEdgeToGlobalMatrixMap[edgematrixoffset+k]=vGlobal;
524 EdgeBlockToUniversalMap[vGlobal] = uniOffset + k + 1;
526 edgematrixoffset+=nedgemodes*nedgemodes;
533 for(j = 0; j < locExpansion->GetNfaces(); ++j)
541 if(faceDirMap.count(meshFaceId)==0)
544 int faceOffset = faceglobaloffset[uniqueFaceMap[meshFaceId]];
547 int uniOffset = meshFaceId;
549 pIt = periodicFaces.find(meshFaceId);
550 if (pIt != periodicFaces.end())
552 uniOffset = min(uniOffset, pIt->second[0].id);
554 uniOffset = uniOffset * maxFaceDof * maxFaceDof;
556 for(k=0; k<nfacemodes*nfacemodes; ++k)
558 vGlobal=faceOffset+k;
560 localFaceToGlobalMatrixMap[facematrixoffset+k]
563 FaceBlockToUniversalMap[vGlobal] = uniOffset + k + 1;
565 facematrixoffset+=nfacemodes*nfacemodes;
588 for(cnt=n=0; n < n_exp; ++n)
590 eid = expList->GetOffset_Elmt_Id(n);
592 locExpansion = expList->GetExp(eid);
593 nCoeffs=locExpansion->NumBndryCoeffs();
597 R=(*(transmatrixmap[eType]));
598 RT=(*(transposedtransmatrixmap[eType]));
601 (nCoeffs, nCoeffs, zero, storage);
605 (nCoeffs, nCoeffs, zero, storage);
608 nVerts=locExpansion->GetGeom()->GetNumVerts();
609 nEdges=locExpansion->GetGeom()->GetNumEdges();
610 nFaces=locExpansion->GetGeom()->GetNumFaces();
613 loc_mat = (
m_linsys.lock())->GetStaticCondBlock(n);
616 bnd_mat=loc_mat->GetBlock(0,0);
619 offset = bnd_mat->GetRows();
631 for (v=0; v<nVerts; ++v)
633 vMap1=locExpansion->GetVertexMap(v);
637 GetLocalToGlobalBndMap(cnt+vMap1)-nDirBnd;
641 for (m=0; m<nVerts; ++m)
643 vMap2=locExpansion->GetVertexMap(m);
648 GetLocalToGlobalBndMap(cnt+vMap2)-nDirBnd;
651 if (globalcol == globalrow)
655 GetLocalToGlobalBndSign(cnt + vMap1);
657 GetLocalToGlobalBndSign(cnt + vMap2);
660 += sign1*sign2*RSRT(vMap1,vMap2);
665 pIt = periodicVerts.find(meshVertId);
666 if (pIt != periodicVerts.end())
668 for (k = 0; k < pIt->second.size(); ++k)
670 meshVertId = min(meshVertId, pIt->second[k].id);
674 VertBlockToUniversalMap[globalrow]
682 for (eid=0; eid<nEdges; ++eid)
684 nedgemodes=locExpansion->GetEdgeNcoeffs(eid)-2;
688 (nedgemodes,nedgemodes,zero,storage);
693 if(edgeDirMap.count(meshEdgeId)==0)
695 for (v=0; v<nedgemodes; ++v)
697 eMap1=edgemodearray[v];
699 for (m=0; m<nedgemodes; ++m)
701 eMap2=edgemodearray[m];
705 GetLocalToGlobalBndSign(cnt + eMap1);
707 GetLocalToGlobalBndSign(cnt + eMap2);
709 NekDouble globalEdgeValue = sign1*sign2*RSRT(eMap1,eMap2);
711 EdgeBlockArray[edgematrixoffset+v*nedgemodes+m]=globalEdgeValue;
714 edgematrixoffset+=nedgemodes*nedgemodes;
719 for (fid=0; fid<nFaces; ++fid)
721 nfacemodes = locExpansion->GetFaceIntNcoeffs(fid);
725 (nfacemodes,nfacemodes,zero,storage);
729 if(faceDirMap.count(meshFaceId)==0)
734 pIt = periodicFaces.find(meshFaceId);
735 if (pIt != periodicFaces.end())
737 if(meshFaceId == min(meshFaceId, pIt->second[0].id))
739 facemodearray = locExpansion->GetFaceInverseBoundaryMap(fid,faceOrient);
744 facemodearray = locExpansion->GetFaceInverseBoundaryMap(fid,faceOrient);
747 for (v=0; v<nfacemodes; ++v)
749 fMap1=facemodearray[v];
751 for (m=0; m<nfacemodes; ++m)
753 fMap2=facemodearray[m];
757 GetLocalToGlobalBndSign(cnt + fMap1);
759 GetLocalToGlobalBndSign(cnt + fMap2);
762 NekDouble globalFaceValue = sign1*sign2*RSRT(fMap1,fMap2);
765 FaceBlockArray[facematrixoffset+v*nfacemodes+m]=globalFaceValue;
768 facematrixoffset+=nfacemodes*nfacemodes;
776 m_RBlk->SetBlock(n,n, transmatrixmap[eType]);
777 m_RTBlk->SetBlock(n,n, transposedtransmatrixmap[eType]);
789 if(ntotaledgeentries)
794 localEdgeToGlobalMatrixMap,
803 if(ntotalfaceentries)
808 localFaceToGlobalMatrixMap,
817 for (
int i = 0; i < nNonDirVerts; ++i)
819 VertBlk->SetValue(i,i,1.0/vertArray[i]);
828 for(
int loc=0; loc<nNonDirEdgeIDs; ++loc)
830 nedgemodes = n_blks[1+loc];
832 (nedgemodes,nedgemodes,zero,storage);
834 for (v=0; v<nedgemodes; ++v)
836 for (m=0; m<nedgemodes; ++m)
838 NekDouble EdgeValue = GlobalEdgeBlock[offset+v*nedgemodes+m];
839 gmat->SetValue(v,m,EdgeValue);
844 m_BlkMat->SetBlock(1+loc,1+loc, gmat);
845 offset+=edgemodeoffset[loc];
852 for(
int loc=0; loc<nNonDirFaceIDs; ++loc)
854 nfacemodes=n_blks[1+nNonDirEdgeIDs+loc];
856 (nfacemodes,nfacemodes,zero,storage);
858 for (v=0; v<nfacemodes; ++v)
860 for (m=0; m<nfacemodes; ++m)
862 NekDouble FaceValue = GlobalFaceBlock[offset+v*nfacemodes+m];
863 gmat->SetValue(v,m,FaceValue);
867 m_BlkMat->SetBlock(1+nNonDirEdgeIDs+loc,1+nNonDirEdgeIDs+loc, gmat);
868 offset+=facemodeoffset[loc];
871 int totblks=
m_BlkMat->GetNumberOfBlockRows();
872 for (i=1; i< totblks; ++i)
874 unsigned int nmodes=
m_BlkMat->GetNumberOfRowsInBlockRow(i);
879 (nmodes,nmodes,zero,storage);
900 int nNonDir = nGlobal-nDir;
917 boost::shared_ptr<MultiRegions::ExpList>
918 expList=((
m_linsys.lock())->GetLocMat()).lock();
926 int n_exp=expList->GetNumElmts();
929 map<LibUtilities::ShapeType,DNekScalMatSharedPtr> transmatrixmap;
930 map<LibUtilities::ShapeType,DNekScalMatSharedPtr> transposedtransmatrixmap;
931 map<LibUtilities::ShapeType,DNekScalMatSharedPtr> invtransmatrixmap;
932 map<LibUtilities::ShapeType,DNekScalMatSharedPtr> invtransposedtransmatrixmap;
966 for(n=0; n < n_exp; ++n)
968 nel = expList->GetOffset_Elmt_Id(n);
970 locExpansion = expList->GetExp(nel);
974 m_RBlk->SetBlock(n,n, transmatrixmap[eType]);
977 m_RTBlk->SetBlock(n,n, transposedtransmatrixmap[eType]);
980 m_InvRBlk->SetBlock(n,n, invtransmatrixmap[eType]);
983 m_InvRTBlk->SetBlock(n,n, invtransposedtransmatrixmap[eType]);
1001 int nDirBndDofs =
m_locToGloMap->GetNumGlobalDirBndCoeffs();
1002 int nGlobHomBndDofs = nGlobBndDofs - nDirBndDofs;
1018 Vmath::Vcopy(nGlobBndDofs, pInOut.get(), 1, tmp.get(), 1);
1024 F_LocBnd=R*F_LocBnd;
1042 int nDirBndDofs =
m_locToGloMap->GetNumGlobalDirBndCoeffs();
1043 int nGlobHomBndDofs = nGlobBndDofs - nDirBndDofs;
1047 ASSERTL1(pInput.num_elements() >= nGlobHomBndDofs,
1048 "Input array is greater than the nGlobHomBndDofs");
1049 ASSERTL1(pOutput.num_elements() >= nGlobHomBndDofs,
1050 "Output array is greater than the nGlobHomBndDofs");
1067 Vmath::Vcopy(nGlobHomBndDofs, pInput.get(), 1, tmp.get() + nDirBndDofs, 1);
1073 F_LocBnd=R*F_LocBnd;
1092 int nDirBndDofs =
m_locToGloMap->GetNumGlobalDirBndCoeffs();
1093 int nGlobHomBndDofs = nGlobBndDofs - nDirBndDofs;
1096 ASSERTL1(pInOut.num_elements() >= nGlobBndDofs,
1097 "Output array is greater than the nGlobBndDofs");
1111 m_locToGloMap->GlobalToLocalBnd(V_GlobHomBnd,V_LocBnd, nDirBndDofs);
1114 V_LocBnd=RT*V_LocBnd;
1124 Vmath::Vcopy(nGlobBndDofs-nDirBndDofs, tmp.get() + nDirBndDofs, 1, pInOut.get() + nDirBndDofs, 1);
1135 int nDirBndDofs =
m_locToGloMap->GetNumGlobalDirBndCoeffs();
1136 int nGlobHomBndDofs = nGlobBndDofs - nDirBndDofs;
1139 ASSERTL1(pInput.num_elements() >= nGlobHomBndDofs,
1140 "Input array is greater than the nGlobHomBndDofs");
1141 ASSERTL1(pOutput.num_elements() >= nGlobHomBndDofs,
1142 "Output array is greater than the nGlobHomBndDofs");
1159 Vmath::Vcopy(nGlobHomBndDofs, pInput.get(), 1, tmp.get() + nDirBndDofs, 1);
1165 F_LocBnd=invR*F_LocBnd;
1180 int nDirBndDofs =
m_locToGloMap->GetNumGlobalDirBndCoeffs();
1181 int nGlobHomBndDofs = nGlobBndDofs - nDirBndDofs;
1184 ASSERTL1(pInput.num_elements() >= nGlobHomBndDofs,
1185 "Input array is greater than the nGlobHomBndDofs");
1186 ASSERTL1(pOutput.num_elements() >= nGlobHomBndDofs,
1187 "Output array is greater than the nGlobHomBndDofs");
1200 m_locToGloMap->GlobalToLocalBnd(pInput,pLocal, nDirBndDofs);
1203 F_LocBnd=invRT*F_LocBnd;
1221 const boost::shared_ptr<DNekScalMat > &loc_mat)
1223 boost::shared_ptr<MultiRegions::ExpList>
1224 expList=((
m_linsys.lock())->GetLocMat()).lock();
1227 locExpansion = expList->GetExp(offset);
1228 unsigned int nbnd=locExpansion->NumBndryCoeffs();
1234 map<LibUtilities::ShapeType,DNekScalMatSharedPtr> transmatrixmap;
1235 map<LibUtilities::ShapeType,DNekScalMatSharedPtr> transposedtransmatrixmap;
1251 (expList->GetExp(offset))->DetShapeType();
1255 DNekScalMat &RT = (*(transposedtransmatrixmap[eType]));
1276 unsigned int nGlobalBnd =
m_locToGloMap->GetNumGlobalBndCoeffs();
1277 unsigned int nEntries =
m_locToGloMap->GetNumLocalBndCoeffs();
1290 for (i = 0; i < nEntries; ++i)
1292 vCounts[vMap[i]] += 1.0;
1300 for (i = 0; i < nEntries; ++i)
1313 int nGlobHomBnd = nGlobalBnd - nDirBnd;
1338 const int nVerts = 6;
1339 const double point[][3] = {
1340 {-1,-1,0}, {1,-1,0}, {1,1,0},
1341 {-1,1,0}, {0,-1,sqrt(
double(3))}, {0,1,sqrt(
double(3))},
1346 for(
int i=0; i < nVerts; ++i)
1349 ( three, i, point[i][0], point[i][1], point[i][2] );
1351 const int nEdges = 9;
1352 const int vertexConnectivity[][2] = {
1353 {0,1}, {1,2}, {3,2}, {0,3}, {0,4},
1354 {1,4}, {2,5}, {3,5}, {4,5}
1359 for(
int i=0; i < nEdges; ++i){
1361 for(
int j=0; j<2; ++j)
1363 vertsArray[j] = verts[vertexConnectivity[i][j]];
1372 const int nFaces = 5;
1374 const int quadEdgeConnectivity[][4] = { {0,1,2,3}, {1,6,8,5}, {3,7,8,4} };
1375 const bool isQuadEdgeFlipped[][4] = { {0,0,1,1}, {0,0,1,1}, {0,0,1,1} };
1377 const int quadId[] = { 0,-1,1,-1,2 };
1380 const int triEdgeConnectivity[][3] = { {0,5,4}, {2,6,7} };
1381 const bool isTriEdgeFlipped[][3] = { {0,0,1}, {0,0,1} };
1383 const int triId[] = { -1,0,-1,1,-1 };
1387 for(
int f = 0; f < nFaces; ++f){
1388 if(f == 1 || f == 3) {
1392 for(
int j = 0; j < 3; ++j){
1393 edgeArray[j] = edges[triEdgeConnectivity[i][j]];
1402 for(
int j=0; j < 4; ++j){
1403 edgeArray[j] = edges[quadEdgeConnectivity[i][j]];
1429 const int nVerts = 4;
1430 const double point[][3] = {
1431 {-1,-1/sqrt(
double(3)),-1/sqrt(
double(6))},
1432 {1,-1/sqrt(
double(3)),-1/sqrt(
double(6))},
1433 {0,2/sqrt(
double(3)),-1/sqrt(
double(6))},
1434 {0,0,3/sqrt(
double(6))}};
1436 boost::shared_ptr<SpatialDomains::PointGeom> verts[4];
1437 for(i=0; i < nVerts; ++i)
1442 ( three, i, point[i][0], point[i][1], point[i][2] );
1450 const int nEdges = 6;
1451 const int vertexConnectivity[][2] = {
1452 {0,1},{1,2},{0,2},{0,3},{1,3},{2,3}
1457 for(i=0; i < nEdges; ++i)
1459 boost::shared_ptr<SpatialDomains::PointGeom>
1463 vertsArray[j] = verts[vertexConnectivity[i][j]];
1474 const int nFaces = 4;
1475 const int edgeConnectivity[][3] = {
1476 {0,1,2}, {0,4,3}, {1,5,4}, {2,5,3}
1478 const bool isEdgeFlipped[][3] = {
1479 {0,0,1}, {0,0,1}, {0,0,1}, {0,0,1}
1484 for(i=0; i < nFaces; ++i)
1488 for(j=0; j < 3; ++j)
1490 edgeArray[j] = edges[edgeConnectivity[i][j]];
1491 eorientArray[j] = isEdgeFlipped[i][j] ?
1521 const int nVerts = 8;
1522 const double point[][3] = {
1523 {0,0,0}, {1,0,0}, {1,1,0}, {0,1,0},
1524 {0,0,1}, {1,0,1}, {1,1,1}, {0,1,1}
1529 for(
int i = 0; i < nVerts; ++i ) {
1532 point[i][1], point[i][2]);
1540 const int nEdges = 12;
1541 const int vertexConnectivity[][2] = {
1542 {0,1}, {1,2}, {2,3}, {0,3}, {0,4}, {1,5},
1543 {2,6}, {3,7}, {4,5}, {5,6}, {6,7}, {4,7}
1548 for(
int i = 0; i < nEdges; ++i ) {
1550 for(
int j = 0; j < 2; ++j ) {
1551 vertsArray[j] = verts[vertexConnectivity[i][j]];
1561 const int nFaces = 6;
1562 const int edgeConnectivity[][4] = {
1563 {0,1,2,3}, {0,5,8,4}, {1,6,9,5},
1564 {2,7,10,6}, {3,7,11,4}, {8,9,10,11}
1566 const bool isEdgeFlipped[][4] = {
1567 {0,0,0,1}, {0,0,1,1}, {0,0,1,1},
1568 {0,0,1,1}, {0,0,1,1}, {0,0,0,1}
1573 for(
int i = 0; i < nFaces; ++i ) {
1576 for(
int j = 0; j < 4; ++j ) {
1577 edgeArray[j] = edges[edgeConnectivity[i][j]];
1578 eorientArray[j] = isEdgeFlipped[i][j] ?
1605 boost::shared_ptr<MultiRegions::ExpList>
1606 expList=((
m_linsys.lock())->GetLocMat()).lock();
1610 locExpansion = expList->GetExp(0);
1617 DNekMatSharedPtr Rtettmp, RTtettmp, Rhextmp, RThextmp, Rprismtmp, RTprismtmp ;
1633 int nummodes=locExpansion->GetBasisNumModes(0);
1693 StdRegions::VarCoeffMap::const_iterator x;
1694 cnt = expList->GetPhys_Offset(0);
1698 vVarCoeffMap[x->first] = x->second + cnt;
1765 m_Rtet = TetExp->GetLocMatrix(TetR);
1768 m_RTtet = TetExp->GetLocMatrix(TetRT);
1772 Rtettmp=TetExp->BuildInverseTransformationMatrix(m_Rtet);
1775 RTtettmp=TetExp->BuildInverseTransformationMatrix(m_Rtet);
1776 RTtettmp->Transpose();
1788 m_Rhex = HexExp->GetLocMatrix(HexR);
1790 m_RThex = HexExp->GetLocMatrix(HexRT);
1794 Rhextmp=HexExp->BuildInverseTransformationMatrix(
m_Rhex);
1796 RThextmp=HexExp->BuildInverseTransformationMatrix(
m_Rhex);
1797 RThextmp->Transpose();
1809 Rprismoriginal = PrismExp->GetLocMatrix(PrismR);
1811 RTprismoriginal = PrismExp->GetLocMatrix(PrismRT);
1813 unsigned int nRows=Rprismoriginal->GetRows();
1822 for(i=0; i<nRows; ++i)
1824 for(j=0; j<nRows; ++j)
1826 Rvalue=(*Rprismoriginal)(i,j);
1827 RTvalue=(*RTprismoriginal)(i,j);
1828 Rtmpprism->SetValue(i,j,Rvalue);
1829 RTtmpprism->SetValue(i,j,RTvalue);
1845 Rprismtmp=PrismExp->BuildInverseTransformationMatrix(
m_Rprism);
1848 RTprismtmp=PrismExp->BuildInverseTransformationMatrix(
m_Rprism);
1849 RTprismtmp->Transpose();
1886 int TetVertex0=TetExp->GetVertexMap(0);
1887 int TetVertex1=TetExp->GetVertexMap(1);
1888 int TetVertex2=TetExp->GetVertexMap(2);
1889 int TetVertex3=TetExp->GetVertexMap(3);
1906 int PrismVertex0=PrismExp->GetVertexMap(0);
1907 int PrismVertex1=PrismExp->GetVertexMap(1);
1908 int PrismVertex2=PrismExp->GetVertexMap(2);
1909 int PrismVertex3=PrismExp->GetVertexMap(3);
1910 int PrismVertex4=PrismExp->GetVertexMap(4);
1911 int PrismVertex5=PrismExp->GetVertexMap(5);
1915 PrismExp->GetEdgeInverseBoundaryMap(0);
1917 PrismExp->GetEdgeInverseBoundaryMap(1);
1919 PrismExp->GetEdgeInverseBoundaryMap(2);
1921 PrismExp->GetEdgeInverseBoundaryMap(3);
1923 PrismExp->GetEdgeInverseBoundaryMap(4);
1925 PrismExp->GetEdgeInverseBoundaryMap(5);
1927 PrismExp->GetEdgeInverseBoundaryMap(6);
1929 PrismExp->GetEdgeInverseBoundaryMap(7);
1931 PrismExp->GetEdgeInverseBoundaryMap(8);
1935 PrismExp->GetFaceInverseBoundaryMap(1);
1937 PrismExp->GetFaceInverseBoundaryMap(3);
1939 PrismExp->GetFaceInverseBoundaryMap(0);
1941 PrismExp->GetFaceInverseBoundaryMap(2);
1943 PrismExp->GetFaceInverseBoundaryMap(4);
1946 for(i=0; i< PrismEdge0.num_elements(); ++i)
1948 Rvalue=(*m_Rtet)(TetVertex0,TetEdge0[i]);
1949 Rmodprism->SetValue(PrismVertex0,PrismEdge0[i],Rvalue);
1950 Rvalue=(*m_Rtet)(TetVertex0,TetEdge2[i]);
1951 Rmodprism->SetValue(PrismVertex0,PrismEdge3[i],Rvalue);
1952 Rvalue=(*m_Rtet)(TetVertex0,TetEdge3[i]);
1953 Rmodprism->SetValue(PrismVertex0,PrismEdge4[i],Rvalue);
1956 RTvalue=(*m_RTtet)(TetEdge0[i],TetVertex0);
1957 RTmodprism->SetValue(PrismEdge0[i],PrismVertex0,RTvalue);
1958 RTvalue=(*m_RTtet)(TetEdge2[i],TetVertex0);
1959 RTmodprism->SetValue(PrismEdge3[i],PrismVertex0,RTvalue);
1960 RTvalue=(*m_RTtet)(TetEdge3[i],TetVertex0);
1961 RTmodprism->SetValue(PrismEdge4[i],PrismVertex0,RTvalue);
1965 for(i=0; i< PrismEdge1.num_elements(); ++i)
1967 Rvalue=(*m_Rtet)(TetVertex1,TetEdge0[i]);
1968 Rmodprism->SetValue(PrismVertex1,PrismEdge0[i],Rvalue);
1969 Rvalue=(*m_Rtet)(TetVertex1,TetEdge1[i]);
1970 Rmodprism->SetValue(PrismVertex1,PrismEdge1[i],Rvalue);
1971 Rvalue=(*m_Rtet)(TetVertex1,TetEdge4[i]);
1972 Rmodprism->SetValue(PrismVertex1,PrismEdge5[i],Rvalue);
1975 RTvalue=(*m_RTtet)(TetEdge0[i],TetVertex1);
1976 RTmodprism->SetValue(PrismEdge0[i],PrismVertex1,RTvalue);
1977 RTvalue=(*m_RTtet)(TetEdge1[i],TetVertex1);
1978 RTmodprism->SetValue(PrismEdge1[i],PrismVertex1,RTvalue);
1979 RTvalue=(*m_RTtet)(TetEdge4[i],TetVertex1);
1980 RTmodprism->SetValue(PrismEdge5[i],PrismVertex1,RTvalue);
1984 for(i=0; i< PrismEdge2.num_elements(); ++i)
1986 Rvalue=(*m_Rtet)(TetVertex2,TetEdge1[i]);
1987 Rmodprism->SetValue(PrismVertex2,PrismEdge1[i],Rvalue);
1988 Rvalue=(*m_Rtet)(TetVertex1,TetEdge0[i]);
1989 Rmodprism->SetValue(PrismVertex2,PrismEdge2[i],Rvalue);
1990 Rvalue=(*m_Rtet)(TetVertex2,TetEdge5[i]);
1991 Rmodprism->SetValue(PrismVertex2,PrismEdge6[i],Rvalue);
1994 RTvalue=(*m_RTtet)(TetEdge1[i],TetVertex2);
1995 RTmodprism->SetValue(PrismEdge1[i],PrismVertex2,RTvalue);
1996 RTvalue=(*m_RTtet)(TetEdge0[i],TetVertex1);
1997 RTmodprism->SetValue(PrismEdge2[i],PrismVertex2,RTvalue);
1998 RTvalue=(*m_RTtet)(TetEdge5[i],TetVertex2);
1999 RTmodprism->SetValue(PrismEdge6[i],PrismVertex2,RTvalue);
2003 for(i=0; i< PrismEdge3.num_elements(); ++i)
2005 Rvalue=(*m_Rtet)(TetVertex2,TetEdge2[i]);
2006 Rmodprism->SetValue(PrismVertex3,PrismEdge3[i],Rvalue);
2007 Rvalue=(*m_Rtet)(TetVertex0,TetEdge0[i]);
2008 Rmodprism->SetValue(PrismVertex3,PrismEdge2[i],Rvalue);
2009 Rvalue=(*m_Rtet)(TetVertex2,TetEdge5[i]);
2010 Rmodprism->SetValue(PrismVertex3,PrismEdge7[i],Rvalue);
2013 RTvalue=(*m_RTtet)(TetEdge2[i],TetVertex2);
2014 RTmodprism->SetValue(PrismEdge3[i],PrismVertex3,RTvalue);
2015 RTvalue=(*m_RTtet)(TetEdge0[i],TetVertex0);
2016 RTmodprism->SetValue(PrismEdge2[i],PrismVertex3,RTvalue);
2017 RTvalue=(*m_RTtet)(TetEdge5[i],TetVertex2);
2018 RTmodprism->SetValue(PrismEdge7[i],PrismVertex3,RTvalue);
2022 for(i=0; i< PrismEdge4.num_elements(); ++i)
2024 Rvalue=(*m_Rtet)(TetVertex3,TetEdge3[i]);
2025 Rmodprism->SetValue(PrismVertex4,PrismEdge4[i],Rvalue);
2026 Rvalue=(*m_Rtet)(TetVertex3,TetEdge4[i]);
2027 Rmodprism->SetValue(PrismVertex4,PrismEdge5[i],Rvalue);
2028 Rvalue=(*m_Rtet)(TetVertex0,TetEdge2[i]);
2029 Rmodprism->SetValue(PrismVertex4,PrismEdge8[i],Rvalue);
2032 RTvalue=(*m_RTtet)(TetEdge3[i],TetVertex3);
2033 RTmodprism->SetValue(PrismEdge4[i],PrismVertex4,RTvalue);
2034 RTvalue=(*m_RTtet)(TetEdge4[i],TetVertex3);
2035 RTmodprism->SetValue(PrismEdge5[i],PrismVertex4,RTvalue);
2036 RTvalue=(*m_RTtet)(TetEdge2[i],TetVertex0);
2037 RTmodprism->SetValue(PrismEdge8[i],PrismVertex4,RTvalue);
2041 for(i=0; i< PrismEdge5.num_elements(); ++i)
2043 Rvalue=(*m_Rtet)(TetVertex3,TetEdge3[i]);
2044 Rmodprism->SetValue(PrismVertex5,PrismEdge6[i],Rvalue);
2045 Rvalue=(*m_Rtet)(TetVertex3,TetEdge4[i]);
2046 Rmodprism->SetValue(PrismVertex5,PrismEdge7[i],Rvalue);
2047 Rvalue=(*m_Rtet)(TetVertex2,TetEdge2[i]);
2048 Rmodprism->SetValue(PrismVertex5,PrismEdge8[i],Rvalue);
2051 RTvalue=(*m_RTtet)(TetEdge3[i],TetVertex3);
2052 RTmodprism->SetValue(PrismEdge6[i],PrismVertex5,RTvalue);
2053 RTvalue=(*m_RTtet)(TetEdge4[i],TetVertex3);
2054 RTmodprism->SetValue(PrismEdge7[i],PrismVertex5,RTvalue);
2055 RTvalue=(*m_RTtet)(TetEdge2[i],TetVertex2);
2056 RTmodprism->SetValue(PrismEdge8[i],PrismVertex5,RTvalue);
2060 for(i=0; i< PrismFace1.num_elements(); ++i)
2062 Rvalue=(*m_Rtet)(TetVertex0,TetFace[i]);
2063 Rmodprism->SetValue(PrismVertex0,PrismFace1[i],Rvalue);
2064 Rvalue=(*m_Rtet)(TetVertex1,TetFace[i]);
2065 Rmodprism->SetValue(PrismVertex1,PrismFace1[i],Rvalue);
2066 Rvalue=(*m_Rtet)(TetVertex3,TetFace[i]);
2067 Rmodprism->SetValue(PrismVertex4,PrismFace1[i],Rvalue);
2070 RTvalue=(*m_RTtet)(TetFace[i],TetVertex0);
2071 RTmodprism->SetValue(PrismFace1[i],PrismVertex0,RTvalue);
2072 RTvalue=(*m_RTtet)(TetFace[i],TetVertex1);
2073 RTmodprism->SetValue(PrismFace1[i],PrismVertex1,RTvalue);
2074 RTvalue=(*m_RTtet)(TetFace[i],TetVertex3);
2075 RTmodprism->SetValue(PrismFace1[i],PrismVertex4,RTvalue);
2079 for(i=0; i< PrismFace3.num_elements(); ++i)
2081 Rvalue=(*m_Rtet)(TetVertex1,TetFace[i]);
2082 Rmodprism->SetValue(PrismVertex2,PrismFace3[i],Rvalue);
2083 Rvalue=(*m_Rtet)(TetVertex0,TetFace[i]);
2084 Rmodprism->SetValue(PrismVertex3,PrismFace3[i],Rvalue);
2085 Rvalue=(*m_Rtet)(TetVertex3,TetFace[i]);
2086 Rmodprism->SetValue(PrismVertex5,PrismFace3[i],Rvalue);
2089 RTvalue=(*m_RTtet)(TetFace[i],TetVertex1);
2090 RTmodprism->SetValue(PrismFace3[i],PrismVertex2,RTvalue);
2091 RTvalue=(*m_RTtet)(TetFace[i],TetVertex0);
2092 RTmodprism->SetValue(PrismFace3[i],PrismVertex3,RTvalue);
2093 RTvalue=(*m_RTtet)(TetFace[i],TetVertex3);
2094 RTmodprism->SetValue(PrismFace3[i],PrismVertex5,RTvalue);
2098 for(i=0; i< PrismFace1.num_elements(); ++i)
2100 for(j=0; j<PrismEdge0.num_elements(); ++j)
2102 Rvalue=(*m_Rtet)(TetEdge0[j],TetFace[i]);
2103 Rmodprism->SetValue(PrismEdge0[j],PrismFace1[i],Rvalue);
2104 Rvalue=(*m_Rtet)(TetEdge3[j],TetFace[i]);
2105 Rmodprism->SetValue(PrismEdge4[j],PrismFace1[i],Rvalue);
2106 Rvalue=(*m_Rtet)(TetEdge4[j],TetFace[i]);
2107 Rmodprism->SetValue(PrismEdge5[j],PrismFace1[i],Rvalue);
2110 RTvalue=(*m_RTtet)(TetFace[i],TetEdge0[j]);
2111 RTmodprism->SetValue(PrismFace1[i],PrismEdge0[j],RTvalue);
2112 RTvalue=(*m_RTtet)(TetFace[i],TetEdge3[j]);
2113 RTmodprism->SetValue(PrismFace1[i],PrismEdge4[j],RTvalue);
2114 RTvalue=(*m_RTtet)(TetFace[i],TetEdge4[j]);
2115 RTmodprism->SetValue(PrismFace1[i],PrismEdge5[j],RTvalue);
2120 for(i=0; i< PrismFace3.num_elements(); ++i)
2122 for(j=0; j<PrismEdge2.num_elements(); ++j)
2124 Rvalue=(*m_Rtet)(TetEdge0[j],TetFace[i]);
2125 Rmodprism->SetValue(PrismEdge2[j],PrismFace3[i],Rvalue);
2126 Rvalue=(*m_Rtet)(TetEdge4[j],TetFace[i]);
2127 Rmodprism->SetValue(PrismEdge6[j],PrismFace3[i],Rvalue);
2128 Rvalue=(*m_Rtet)(TetEdge3[j],TetFace[i]);
2129 Rmodprism->SetValue(PrismEdge7[j],PrismFace3[i],Rvalue);
2131 RTvalue=(*m_RTtet)(TetFace[i],TetEdge0[j]);
2132 RTmodprism->SetValue(PrismFace3[i],PrismEdge2[j],RTvalue);
2133 RTvalue=(*m_RTtet)(TetFace[i],TetEdge4[j]);
2134 RTmodprism->SetValue(PrismFace3[i],PrismEdge6[j],RTvalue);
2135 RTvalue=(*m_RTtet)(TetFace[i],TetEdge3[j]);
2136 RTmodprism->SetValue(PrismFace3[i],PrismEdge7[j],RTvalue);
const StdRegions::VarCoeffMap & GetVarCoeffs() const
LibUtilities::CommSharedPtr m_comm
virtual void v_DoPreconditioner(const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
#define ASSERTL0(condition, msg)
Principle Modified Functions .
static void Gather(Nektar::Array< OneD, NekDouble > pU, gs_op pOp, gs_data *pGsh, Nektar::Array< OneD, NekDouble > pBuffer=NullNekDouble1DArray)
Performs a gather-scatter operation of the provided values.
virtual void v_DoMultiplybyInverseTransposedTransformationMatrix(const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
Multiply by the block tranposed inverse transformation matrix.
void Gathr(int n, const T *x, const int *y, T *z)
Gather vector z[i] = x[y[i]].
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
#define sign(a, b)
return the sign(b)*a
boost::shared_ptr< AssemblyMap > AssemblyMapSharedPtr
DNekScalMatSharedPtr m_Rprism
Array< OneD, NekDouble > m_multiplicity
Array< OneD, NekDouble > m_locToGloSignMult
DNekScalBlkMatSharedPtr m_RTBlk
DNekScalMatSharedPtr m_RThex
boost::shared_ptr< QuadGeom > QuadGeomSharedPtr
virtual void v_InitObject()
static gs_data * Init(const Nektar::Array< OneD, long > pId, const LibUtilities::CommSharedPtr &pComm)
Initialise Gather-Scatter map.
DNekScalBlkMatSharedPtr m_RBlk
Principle Modified Functions .
PreconFactory & GetPreconFactory()
void Sdiv(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha/y.
boost::shared_ptr< HexExp > HexExpSharedPtr
DNekBlkMatSharedPtr m_BlkMat
const StdRegions::ConstFactorMap & GetConstFactors() const
Returns all the constants.
boost::shared_ptr< HexGeom > HexGeomSharedPtr
int GetEdgeNcoeffs(const int i) const
This function returns the number of expansion coefficients belonging to the i-th edge.
boost::shared_ptr< DNekMat > DNekMatSharedPtr
DNekScalMatSharedPtr m_RTinvtet
boost::shared_ptr< TetExp > TetExpSharedPtr
DNekScalMatSharedPtr m_Rinvprism
boost::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
DNekScalBlkMatSharedPtr m_InvRTBlk
Gauss Radau pinned at x=-1, .
SpatialDomains::HexGeomSharedPtr CreateRefHexGeom(void)
Sets up the reference hexahedral element needed to construct a low energy basis.
DNekScalMatSharedPtr m_Rhex
DNekScalMatSharedPtr m_RTprism
boost::shared_ptr< StdExpansion2D > StdExpansion2DSharedPtr
virtual DNekScalMatSharedPtr v_TransformedSchurCompl(int offset, const boost::shared_ptr< DNekScalMat > &loc_mat)
Set up the transformed block matrix system.
virtual void v_BuildPreconditioner()
Construct the low energy preconditioner from .
DNekScalMatSharedPtr m_RTinvhex
int GetNVarCoeffs() const
int GetFaceIntNcoeffs(const int i) const
void SetUpReferenceElements(void)
Sets up the reference elements needed by the preconditioner.
std::map< StdRegions::VarCoeffType, Array< OneD, NekDouble > > VarCoeffMap
boost::shared_ptr< SegGeom > SegGeomSharedPtr
DNekScalBlkMatSharedPtr m_InvRBlk
Principle Modified Functions .
StdRegions::Orientation DeterminePeriodicFaceOrient(StdRegions::Orientation faceOrient, StdRegions::Orientation perFaceOrient)
Determine relative orientation between two faces.
boost::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
virtual void v_DoTransformToLowEnergy(Array< OneD, NekDouble > &pInOut, int offset)
Transform the solution vector vector to low energy.
Defines a specification for a set of points.
void Assmb(int n, const T *x, const int *y, T *z)
Assemble z[y[i]] += x[i]; z should be zero'd first.
boost::shared_ptr< Expansion > ExpansionSharedPtr
boost::shared_ptr< PrismExp > PrismExpSharedPtr
boost::shared_ptr< T > as()
std::map< int, std::vector< PeriodicEntity > > PeriodicMap
DNekScalMatSharedPtr m_Rinvtet
virtual void v_DoMultiplybyInverseTransformationMatrix(const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
Multiply by the block inverse transformation matrix.
Describe a linear system.
virtual void v_DoTransformFromLowEnergy(Array< OneD, NekDouble > &pInOut)
transform the solution vector from low energy back to the original basis.
boost::shared_ptr< AssemblyMap > m_locToGloMap
StdRegions::MatrixType GetMatrixType() const
Return the matrix type.
boost::shared_ptr< Geometry2D > Geometry2DSharedPtr
SpatialDomains::TetGeomSharedPtr CreateRefTetGeom(void)
Sets up the reference tretrahedral element needed to construct a low energy basis.
DNekScalMatSharedPtr m_Rinvhex
void ModifyPrismTransformationMatrix(LocalRegions::TetExpSharedPtr TetExp, LocalRegions::PrismExpSharedPtr PrismExp, DNekMatSharedPtr Rmodprism, DNekMatSharedPtr RTmodprism)
Modify the prism transformation matrix to align with the tetrahedral modes.
DNekScalMatSharedPtr m_Rtet
boost::shared_ptr< PrismGeom > PrismGeomSharedPtr
Gauss Radau pinned at x=-1, .
DNekScalMatSharedPtr m_RTtet
void SetupBlockTransformationMatrix(void)
boost::shared_ptr< TetGeom > TetGeomSharedPtr
boost::shared_ptr< TriGeom > TriGeomSharedPtr
SpatialDomains::PrismGeomSharedPtr CreateRefPrismGeom(void)
Sets up the reference prismatic element needed to construct a low energy basis.
const boost::weak_ptr< GlobalLinSys > m_linsys
void CreateMultiplicityMap(void)
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Describes the specification for a Basis.
boost::shared_ptr< PointGeom > PointGeomSharedPtr
1D Gauss-Lobatto-Legendre quadrature points
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.
DNekScalMatSharedPtr m_RTinvprism
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.