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)
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)
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)
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)
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)
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]);
781 expList->GetSession()->DefinesCmdLineArgument(
"verbose");
792 if(ntotaledgeentries)
797 localEdgeToGlobalMatrixMap,
806 if(ntotalfaceentries)
811 localFaceToGlobalMatrixMap,
820 for (
int i = 0; i < nNonDirVerts; ++i)
822 VertBlk->SetValue(i,i,1.0/vertArray[i]);
831 for(
int loc=0; loc<nNonDirEdgeIDs; ++loc)
833 nedgemodes = n_blks[1+loc];
835 (nedgemodes,nedgemodes,zero,storage);
837 for (v=0; v<nedgemodes; ++v)
839 for (m=0; m<nedgemodes; ++m)
841 NekDouble EdgeValue = GlobalEdgeBlock[offset+v*nedgemodes+m];
842 gmat->SetValue(v,m,EdgeValue);
847 m_BlkMat->SetBlock(1+loc,1+loc, gmat);
848 offset+=edgemodeoffset[loc];
855 for(
int loc=0; loc<nNonDirFaceIDs; ++loc)
857 nfacemodes=n_blks[1+nNonDirEdgeIDs+loc];
859 (nfacemodes,nfacemodes,zero,storage);
861 for (v=0; v<nfacemodes; ++v)
863 for (m=0; m<nfacemodes; ++m)
865 NekDouble FaceValue = GlobalFaceBlock[offset+v*nfacemodes+m];
866 gmat->SetValue(v,m,FaceValue);
870 m_BlkMat->SetBlock(1+nNonDirEdgeIDs+loc,1+nNonDirEdgeIDs+loc, gmat);
871 offset+=facemodeoffset[loc];
874 int totblks=
m_BlkMat->GetNumberOfBlockRows();
875 for (i=1; i< totblks; ++i)
877 unsigned int nmodes=
m_BlkMat->GetNumberOfRowsInBlockRow(i);
882 (nmodes,nmodes,zero,storage);
903 int nNonDir = nGlobal-nDir;
920 boost::shared_ptr<MultiRegions::ExpList>
921 expList=((
m_linsys.lock())->GetLocMat()).lock();
929 int n_exp=expList->GetNumElmts();
932 map<LibUtilities::ShapeType,DNekScalMatSharedPtr> transmatrixmap;
933 map<LibUtilities::ShapeType,DNekScalMatSharedPtr> transposedtransmatrixmap;
934 map<LibUtilities::ShapeType,DNekScalMatSharedPtr> invtransmatrixmap;
935 map<LibUtilities::ShapeType,DNekScalMatSharedPtr> invtransposedtransmatrixmap;
969 for(n=0; n < n_exp; ++n)
973 locExpansion = expList->GetExp(nel);
977 m_RBlk->SetBlock(n,n, transmatrixmap[eType]);
980 m_RTBlk->SetBlock(n,n, transposedtransmatrixmap[eType]);
983 m_InvRBlk->SetBlock(n,n, invtransmatrixmap[eType]);
986 m_InvRTBlk->SetBlock(n,n, invtransposedtransmatrixmap[eType]);
1004 int nDirBndDofs =
m_locToGloMap->GetNumGlobalDirBndCoeffs();
1005 int nGlobHomBndDofs = nGlobBndDofs - nDirBndDofs;
1021 Vmath::Vcopy(nGlobBndDofs, pInOut.get(), 1, tmp.get(), 1);
1027 F_LocBnd=R*F_LocBnd;
1045 int nDirBndDofs =
m_locToGloMap->GetNumGlobalDirBndCoeffs();
1046 int nGlobHomBndDofs = nGlobBndDofs - nDirBndDofs;
1050 ASSERTL1(pInput.num_elements() >= nGlobHomBndDofs,
1051 "Input array is greater than the nGlobHomBndDofs");
1052 ASSERTL1(pOutput.num_elements() >= nGlobHomBndDofs,
1053 "Output array is greater than the nGlobHomBndDofs");
1070 Vmath::Vcopy(nGlobHomBndDofs, pInput.get(), 1, tmp.get() + nDirBndDofs, 1);
1076 F_LocBnd=R*F_LocBnd;
1095 int nDirBndDofs =
m_locToGloMap->GetNumGlobalDirBndCoeffs();
1096 int nGlobHomBndDofs = nGlobBndDofs - nDirBndDofs;
1099 ASSERTL1(pInOut.num_elements() >= nGlobBndDofs,
1100 "Output array is greater than the nGlobBndDofs");
1114 m_locToGloMap->GlobalToLocalBnd(V_GlobHomBnd,V_LocBnd, nDirBndDofs);
1117 V_LocBnd=RT*V_LocBnd;
1127 Vmath::Vcopy(nGlobBndDofs-nDirBndDofs, tmp.get() + nDirBndDofs, 1, pInOut.get() + nDirBndDofs, 1);
1138 int nDirBndDofs =
m_locToGloMap->GetNumGlobalDirBndCoeffs();
1139 int nGlobHomBndDofs = nGlobBndDofs - nDirBndDofs;
1142 ASSERTL1(pInput.num_elements() >= nGlobHomBndDofs,
1143 "Input array is greater than the nGlobHomBndDofs");
1144 ASSERTL1(pOutput.num_elements() >= nGlobHomBndDofs,
1145 "Output array is greater than the nGlobHomBndDofs");
1162 Vmath::Vcopy(nGlobHomBndDofs, pInput.get(), 1, tmp.get() + nDirBndDofs, 1);
1168 F_LocBnd=invR*F_LocBnd;
1183 int nDirBndDofs =
m_locToGloMap->GetNumGlobalDirBndCoeffs();
1184 int nGlobHomBndDofs = nGlobBndDofs - nDirBndDofs;
1187 ASSERTL1(pInput.num_elements() >= nGlobHomBndDofs,
1188 "Input array is greater than the nGlobHomBndDofs");
1189 ASSERTL1(pOutput.num_elements() >= nGlobHomBndDofs,
1190 "Output array is greater than the nGlobHomBndDofs");
1203 m_locToGloMap->GlobalToLocalBnd(pInput,pLocal, nDirBndDofs);
1206 F_LocBnd=invRT*F_LocBnd;
1224 const boost::shared_ptr<DNekScalMat > &loc_mat)
1226 boost::shared_ptr<MultiRegions::ExpList>
1227 expList=((
m_linsys.lock())->GetLocMat()).lock();
1230 locExpansion = expList->GetExp(offset);
1231 unsigned int nbnd=locExpansion->NumBndryCoeffs();
1237 map<LibUtilities::ShapeType,DNekScalMatSharedPtr> transmatrixmap;
1238 map<LibUtilities::ShapeType,DNekScalMatSharedPtr> transposedtransmatrixmap;
1254 (expList->GetExp(offset))->DetShapeType();
1258 DNekScalMat &RT = (*(transposedtransmatrixmap[eType]));
1279 unsigned int nGlobalBnd =
m_locToGloMap->GetNumGlobalBndCoeffs();
1280 unsigned int nEntries =
m_locToGloMap->GetNumLocalBndCoeffs();
1293 for (i = 0; i < nEntries; ++i)
1295 vCounts[vMap[i]] += 1.0;
1303 for (i = 0; i < nEntries; ++i)
1316 int nGlobHomBnd = nGlobalBnd - nDirBnd;
1341 const int nVerts = 6;
1342 const double point[][3] = {
1343 {-1,-1,0}, {1,-1,0}, {1,1,0},
1344 {-1,1,0}, {0,-1,sqrt(
double(3))}, {0,1,sqrt(
double(3))},
1349 for(
int i=0; i < nVerts; ++i)
1352 ( three, i, point[i][0], point[i][1], point[i][2] );
1354 const int nEdges = 9;
1355 const int vertexConnectivity[][2] = {
1356 {0,1}, {1,2}, {3,2}, {0,3}, {0,4},
1357 {1,4}, {2,5}, {3,5}, {4,5}
1362 for(
int i=0; i < nEdges; ++i){
1364 for(
int j=0; j<2; ++j)
1366 vertsArray[j] = verts[vertexConnectivity[i][j]];
1375 const int nFaces = 5;
1377 const int quadEdgeConnectivity[][4] = { {0,1,2,3}, {1,6,8,5}, {3,7,8,4} };
1378 const bool isQuadEdgeFlipped[][4] = { {0,0,1,1}, {0,0,1,1}, {0,0,1,1} };
1380 const int quadId[] = { 0,-1,1,-1,2 };
1383 const int triEdgeConnectivity[][3] = { {0,5,4}, {2,6,7} };
1384 const bool isTriEdgeFlipped[][3] = { {0,0,1}, {0,0,1} };
1386 const int triId[] = { -1,0,-1,1,-1 };
1390 for(
int f = 0; f < nFaces; ++f){
1391 if(f == 1 || f == 3) {
1395 for(
int j = 0; j < 3; ++j){
1396 edgeArray[j] = edges[triEdgeConnectivity[i][j]];
1405 for(
int j=0; j < 4; ++j){
1406 edgeArray[j] = edges[quadEdgeConnectivity[i][j]];
1432 const int nVerts = 4;
1433 const double point[][3] = {
1434 {-1,-1/sqrt(
double(3)),-1/sqrt(
double(6))},
1435 {1,-1/sqrt(
double(3)),-1/sqrt(
double(6))},
1436 {0,2/sqrt(
double(3)),-1/sqrt(
double(6))},
1437 {0,0,3/sqrt(
double(6))}};
1439 boost::shared_ptr<SpatialDomains::PointGeom> verts[4];
1440 for(i=0; i < nVerts; ++i)
1445 ( three, i, point[i][0], point[i][1], point[i][2] );
1453 const int nEdges = 6;
1454 const int vertexConnectivity[][2] = {
1455 {0,1},{1,2},{0,2},{0,3},{1,3},{2,3}
1460 for(i=0; i < nEdges; ++i)
1462 boost::shared_ptr<SpatialDomains::PointGeom>
1466 vertsArray[j] = verts[vertexConnectivity[i][j]];
1477 const int nFaces = 4;
1478 const int edgeConnectivity[][3] = {
1479 {0,1,2}, {0,4,3}, {1,5,4}, {2,5,3}
1481 const bool isEdgeFlipped[][3] = {
1482 {0,0,1}, {0,0,1}, {0,0,1}, {0,0,1}
1487 for(i=0; i < nFaces; ++i)
1491 for(j=0; j < 3; ++j)
1493 edgeArray[j] = edges[edgeConnectivity[i][j]];
1494 eorientArray[j] = isEdgeFlipped[i][j] ?
1524 const int nVerts = 8;
1525 const double point[][3] = {
1526 {0,0,0}, {1,0,0}, {1,1,0}, {0,1,0},
1527 {0,0,1}, {1,0,1}, {1,1,1}, {0,1,1}
1532 for(
int i = 0; i < nVerts; ++i ) {
1535 point[i][1], point[i][2]);
1543 const int nEdges = 12;
1544 const int vertexConnectivity[][2] = {
1545 {0,1}, {1,2}, {2,3}, {0,3}, {0,4}, {1,5},
1546 {2,6}, {3,7}, {4,5}, {5,6}, {6,7}, {4,7}
1551 for(
int i = 0; i < nEdges; ++i ) {
1553 for(
int j = 0; j < 2; ++j ) {
1554 vertsArray[j] = verts[vertexConnectivity[i][j]];
1564 const int nFaces = 6;
1565 const int edgeConnectivity[][4] = {
1566 {0,1,2,3}, {0,5,8,4}, {1,6,9,5},
1567 {2,7,10,6}, {3,7,11,4}, {8,9,10,11}
1569 const bool isEdgeFlipped[][4] = {
1570 {0,0,0,1}, {0,0,1,1}, {0,0,1,1},
1571 {0,0,1,1}, {0,0,1,1}, {0,0,0,1}
1576 for(
int i = 0; i < nFaces; ++i ) {
1579 for(
int j = 0; j < 4; ++j ) {
1580 edgeArray[j] = edges[edgeConnectivity[i][j]];
1581 eorientArray[j] = isEdgeFlipped[i][j] ?
1608 boost::shared_ptr<MultiRegions::ExpList>
1609 expList=((
m_linsys.lock())->GetLocMat()).lock();
1613 locExpansion = expList->GetExp(0);
1620 DNekMatSharedPtr Rtettmp, RTtettmp, Rhextmp, RThextmp, Rprismtmp, RTprismtmp ;
1636 int nummodes=locExpansion->GetBasisNumModes(0);
1696 StdRegions::VarCoeffMap::const_iterator x;
1697 cnt = expList->GetPhys_Offset(0);
1701 vVarCoeffMap[x->first] = x->second + cnt;
1768 m_Rtet = TetExp->GetLocMatrix(TetR);
1771 m_RTtet = TetExp->GetLocMatrix(TetRT);
1775 Rtettmp=TetExp->BuildInverseTransformationMatrix(m_Rtet);
1778 RTtettmp=TetExp->BuildInverseTransformationMatrix(m_Rtet);
1779 RTtettmp->Transpose();
1791 m_Rhex = HexExp->GetLocMatrix(HexR);
1793 m_RThex = HexExp->GetLocMatrix(HexRT);
1797 Rhextmp=HexExp->BuildInverseTransformationMatrix(
m_Rhex);
1799 RThextmp=HexExp->BuildInverseTransformationMatrix(
m_Rhex);
1800 RThextmp->Transpose();
1812 Rprismoriginal = PrismExp->GetLocMatrix(PrismR);
1814 RTprismoriginal = PrismExp->GetLocMatrix(PrismRT);
1816 unsigned int nRows=Rprismoriginal->GetRows();
1825 for(i=0; i<nRows; ++i)
1827 for(j=0; j<nRows; ++j)
1829 Rvalue=(*Rprismoriginal)(i,j);
1830 RTvalue=(*RTprismoriginal)(i,j);
1831 Rtmpprism->SetValue(i,j,Rvalue);
1832 RTtmpprism->SetValue(i,j,RTvalue);
1848 Rprismtmp=PrismExp->BuildInverseTransformationMatrix(
m_Rprism);
1851 RTprismtmp=PrismExp->BuildInverseTransformationMatrix(
m_Rprism);
1852 RTprismtmp->Transpose();
1889 int TetVertex0=TetExp->GetVertexMap(0);
1890 int TetVertex1=TetExp->GetVertexMap(1);
1891 int TetVertex2=TetExp->GetVertexMap(2);
1892 int TetVertex3=TetExp->GetVertexMap(3);
1909 int PrismVertex0=PrismExp->GetVertexMap(0);
1910 int PrismVertex1=PrismExp->GetVertexMap(1);
1911 int PrismVertex2=PrismExp->GetVertexMap(2);
1912 int PrismVertex3=PrismExp->GetVertexMap(3);
1913 int PrismVertex4=PrismExp->GetVertexMap(4);
1914 int PrismVertex5=PrismExp->GetVertexMap(5);
1918 PrismExp->GetEdgeInverseBoundaryMap(0);
1920 PrismExp->GetEdgeInverseBoundaryMap(1);
1922 PrismExp->GetEdgeInverseBoundaryMap(2);
1924 PrismExp->GetEdgeInverseBoundaryMap(3);
1926 PrismExp->GetEdgeInverseBoundaryMap(4);
1928 PrismExp->GetEdgeInverseBoundaryMap(5);
1930 PrismExp->GetEdgeInverseBoundaryMap(6);
1932 PrismExp->GetEdgeInverseBoundaryMap(7);
1934 PrismExp->GetEdgeInverseBoundaryMap(8);
1938 PrismExp->GetFaceInverseBoundaryMap(1);
1940 PrismExp->GetFaceInverseBoundaryMap(3);
1942 PrismExp->GetFaceInverseBoundaryMap(0);
1944 PrismExp->GetFaceInverseBoundaryMap(2);
1946 PrismExp->GetFaceInverseBoundaryMap(4);
1949 for(i=0; i< PrismEdge0.num_elements(); ++i)
1951 Rvalue=(*m_Rtet)(TetVertex0,TetEdge0[i]);
1952 Rmodprism->SetValue(PrismVertex0,PrismEdge0[i],Rvalue);
1953 Rvalue=(*m_Rtet)(TetVertex0,TetEdge2[i]);
1954 Rmodprism->SetValue(PrismVertex0,PrismEdge3[i],Rvalue);
1955 Rvalue=(*m_Rtet)(TetVertex0,TetEdge3[i]);
1956 Rmodprism->SetValue(PrismVertex0,PrismEdge4[i],Rvalue);
1959 RTvalue=(*m_RTtet)(TetEdge0[i],TetVertex0);
1960 RTmodprism->SetValue(PrismEdge0[i],PrismVertex0,RTvalue);
1961 RTvalue=(*m_RTtet)(TetEdge2[i],TetVertex0);
1962 RTmodprism->SetValue(PrismEdge3[i],PrismVertex0,RTvalue);
1963 RTvalue=(*m_RTtet)(TetEdge3[i],TetVertex0);
1964 RTmodprism->SetValue(PrismEdge4[i],PrismVertex0,RTvalue);
1968 for(i=0; i< PrismEdge1.num_elements(); ++i)
1970 Rvalue=(*m_Rtet)(TetVertex1,TetEdge0[i]);
1971 Rmodprism->SetValue(PrismVertex1,PrismEdge0[i],Rvalue);
1972 Rvalue=(*m_Rtet)(TetVertex1,TetEdge1[i]);
1973 Rmodprism->SetValue(PrismVertex1,PrismEdge1[i],Rvalue);
1974 Rvalue=(*m_Rtet)(TetVertex1,TetEdge4[i]);
1975 Rmodprism->SetValue(PrismVertex1,PrismEdge5[i],Rvalue);
1978 RTvalue=(*m_RTtet)(TetEdge0[i],TetVertex1);
1979 RTmodprism->SetValue(PrismEdge0[i],PrismVertex1,RTvalue);
1980 RTvalue=(*m_RTtet)(TetEdge1[i],TetVertex1);
1981 RTmodprism->SetValue(PrismEdge1[i],PrismVertex1,RTvalue);
1982 RTvalue=(*m_RTtet)(TetEdge4[i],TetVertex1);
1983 RTmodprism->SetValue(PrismEdge5[i],PrismVertex1,RTvalue);
1987 for(i=0; i< PrismEdge2.num_elements(); ++i)
1989 Rvalue=(*m_Rtet)(TetVertex2,TetEdge1[i]);
1990 Rmodprism->SetValue(PrismVertex2,PrismEdge1[i],Rvalue);
1991 Rvalue=(*m_Rtet)(TetVertex1,TetEdge0[i]);
1992 Rmodprism->SetValue(PrismVertex2,PrismEdge2[i],Rvalue);
1993 Rvalue=(*m_Rtet)(TetVertex2,TetEdge5[i]);
1994 Rmodprism->SetValue(PrismVertex2,PrismEdge6[i],Rvalue);
1997 RTvalue=(*m_RTtet)(TetEdge1[i],TetVertex2);
1998 RTmodprism->SetValue(PrismEdge1[i],PrismVertex2,RTvalue);
1999 RTvalue=(*m_RTtet)(TetEdge0[i],TetVertex1);
2000 RTmodprism->SetValue(PrismEdge2[i],PrismVertex2,RTvalue);
2001 RTvalue=(*m_RTtet)(TetEdge5[i],TetVertex2);
2002 RTmodprism->SetValue(PrismEdge6[i],PrismVertex2,RTvalue);
2006 for(i=0; i< PrismEdge3.num_elements(); ++i)
2008 Rvalue=(*m_Rtet)(TetVertex2,TetEdge2[i]);
2009 Rmodprism->SetValue(PrismVertex3,PrismEdge3[i],Rvalue);
2010 Rvalue=(*m_Rtet)(TetVertex0,TetEdge0[i]);
2011 Rmodprism->SetValue(PrismVertex3,PrismEdge2[i],Rvalue);
2012 Rvalue=(*m_Rtet)(TetVertex2,TetEdge5[i]);
2013 Rmodprism->SetValue(PrismVertex3,PrismEdge7[i],Rvalue);
2016 RTvalue=(*m_RTtet)(TetEdge2[i],TetVertex2);
2017 RTmodprism->SetValue(PrismEdge3[i],PrismVertex3,RTvalue);
2018 RTvalue=(*m_RTtet)(TetEdge0[i],TetVertex0);
2019 RTmodprism->SetValue(PrismEdge2[i],PrismVertex3,RTvalue);
2020 RTvalue=(*m_RTtet)(TetEdge5[i],TetVertex2);
2021 RTmodprism->SetValue(PrismEdge7[i],PrismVertex3,RTvalue);
2025 for(i=0; i< PrismEdge4.num_elements(); ++i)
2027 Rvalue=(*m_Rtet)(TetVertex3,TetEdge3[i]);
2028 Rmodprism->SetValue(PrismVertex4,PrismEdge4[i],Rvalue);
2029 Rvalue=(*m_Rtet)(TetVertex3,TetEdge4[i]);
2030 Rmodprism->SetValue(PrismVertex4,PrismEdge5[i],Rvalue);
2031 Rvalue=(*m_Rtet)(TetVertex0,TetEdge2[i]);
2032 Rmodprism->SetValue(PrismVertex4,PrismEdge8[i],Rvalue);
2035 RTvalue=(*m_RTtet)(TetEdge3[i],TetVertex3);
2036 RTmodprism->SetValue(PrismEdge4[i],PrismVertex4,RTvalue);
2037 RTvalue=(*m_RTtet)(TetEdge4[i],TetVertex3);
2038 RTmodprism->SetValue(PrismEdge5[i],PrismVertex4,RTvalue);
2039 RTvalue=(*m_RTtet)(TetEdge2[i],TetVertex0);
2040 RTmodprism->SetValue(PrismEdge8[i],PrismVertex4,RTvalue);
2044 for(i=0; i< PrismEdge5.num_elements(); ++i)
2046 Rvalue=(*m_Rtet)(TetVertex3,TetEdge3[i]);
2047 Rmodprism->SetValue(PrismVertex5,PrismEdge6[i],Rvalue);
2048 Rvalue=(*m_Rtet)(TetVertex3,TetEdge4[i]);
2049 Rmodprism->SetValue(PrismVertex5,PrismEdge7[i],Rvalue);
2050 Rvalue=(*m_Rtet)(TetVertex2,TetEdge2[i]);
2051 Rmodprism->SetValue(PrismVertex5,PrismEdge8[i],Rvalue);
2054 RTvalue=(*m_RTtet)(TetEdge3[i],TetVertex3);
2055 RTmodprism->SetValue(PrismEdge6[i],PrismVertex5,RTvalue);
2056 RTvalue=(*m_RTtet)(TetEdge4[i],TetVertex3);
2057 RTmodprism->SetValue(PrismEdge7[i],PrismVertex5,RTvalue);
2058 RTvalue=(*m_RTtet)(TetEdge2[i],TetVertex2);
2059 RTmodprism->SetValue(PrismEdge8[i],PrismVertex5,RTvalue);
2063 for(i=0; i< PrismFace1.num_elements(); ++i)
2065 Rvalue=(*m_Rtet)(TetVertex0,TetFace[i]);
2066 Rmodprism->SetValue(PrismVertex0,PrismFace1[i],Rvalue);
2067 Rvalue=(*m_Rtet)(TetVertex1,TetFace[i]);
2068 Rmodprism->SetValue(PrismVertex1,PrismFace1[i],Rvalue);
2069 Rvalue=(*m_Rtet)(TetVertex3,TetFace[i]);
2070 Rmodprism->SetValue(PrismVertex4,PrismFace1[i],Rvalue);
2073 RTvalue=(*m_RTtet)(TetFace[i],TetVertex0);
2074 RTmodprism->SetValue(PrismFace1[i],PrismVertex0,RTvalue);
2075 RTvalue=(*m_RTtet)(TetFace[i],TetVertex1);
2076 RTmodprism->SetValue(PrismFace1[i],PrismVertex1,RTvalue);
2077 RTvalue=(*m_RTtet)(TetFace[i],TetVertex3);
2078 RTmodprism->SetValue(PrismFace1[i],PrismVertex4,RTvalue);
2082 for(i=0; i< PrismFace3.num_elements(); ++i)
2084 Rvalue=(*m_Rtet)(TetVertex1,TetFace[i]);
2085 Rmodprism->SetValue(PrismVertex2,PrismFace3[i],Rvalue);
2086 Rvalue=(*m_Rtet)(TetVertex0,TetFace[i]);
2087 Rmodprism->SetValue(PrismVertex3,PrismFace3[i],Rvalue);
2088 Rvalue=(*m_Rtet)(TetVertex3,TetFace[i]);
2089 Rmodprism->SetValue(PrismVertex5,PrismFace3[i],Rvalue);
2092 RTvalue=(*m_RTtet)(TetFace[i],TetVertex1);
2093 RTmodprism->SetValue(PrismFace3[i],PrismVertex2,RTvalue);
2094 RTvalue=(*m_RTtet)(TetFace[i],TetVertex0);
2095 RTmodprism->SetValue(PrismFace3[i],PrismVertex3,RTvalue);
2096 RTvalue=(*m_RTtet)(TetFace[i],TetVertex3);
2097 RTmodprism->SetValue(PrismFace3[i],PrismVertex5,RTvalue);
2101 for(i=0; i< PrismFace1.num_elements(); ++i)
2103 for(j=0; j<PrismEdge0.num_elements(); ++j)
2105 Rvalue=(*m_Rtet)(TetEdge0[j],TetFace[i]);
2106 Rmodprism->SetValue(PrismEdge0[j],PrismFace1[i],Rvalue);
2107 Rvalue=(*m_Rtet)(TetEdge3[j],TetFace[i]);
2108 Rmodprism->SetValue(PrismEdge4[j],PrismFace1[i],Rvalue);
2109 Rvalue=(*m_Rtet)(TetEdge4[j],TetFace[i]);
2110 Rmodprism->SetValue(PrismEdge5[j],PrismFace1[i],Rvalue);
2113 RTvalue=(*m_RTtet)(TetFace[i],TetEdge0[j]);
2114 RTmodprism->SetValue(PrismFace1[i],PrismEdge0[j],RTvalue);
2115 RTvalue=(*m_RTtet)(TetFace[i],TetEdge3[j]);
2116 RTmodprism->SetValue(PrismFace1[i],PrismEdge4[j],RTvalue);
2117 RTvalue=(*m_RTtet)(TetFace[i],TetEdge4[j]);
2118 RTmodprism->SetValue(PrismFace1[i],PrismEdge5[j],RTvalue);
2123 for(i=0; i< PrismFace3.num_elements(); ++i)
2125 for(j=0; j<PrismEdge2.num_elements(); ++j)
2127 Rvalue=(*m_Rtet)(TetEdge0[j],TetFace[i]);
2128 Rmodprism->SetValue(PrismEdge2[j],PrismFace3[i],Rvalue);
2129 Rvalue=(*m_Rtet)(TetEdge4[j],TetFace[i]);
2130 Rmodprism->SetValue(PrismEdge6[j],PrismFace3[i],Rvalue);
2131 Rvalue=(*m_Rtet)(TetEdge3[j],TetFace[i]);
2132 Rmodprism->SetValue(PrismEdge7[j],PrismFace3[i],Rvalue);
2134 RTvalue=(*m_RTtet)(TetFace[i],TetEdge0[j]);
2135 RTmodprism->SetValue(PrismFace3[i],PrismEdge2[j],RTvalue);
2136 RTvalue=(*m_RTtet)(TetFace[i],TetEdge4[j]);
2137 RTmodprism->SetValue(PrismFace3[i],PrismEdge6[j],RTvalue);
2138 RTvalue=(*m_RTtet)(TetFace[i],TetEdge3[j]);
2139 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()
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
static gs_data * Init(const Nektar::Array< OneD, long > pId, const LibUtilities::CommSharedPtr &pComm, bool verbose=true)
Initialise Gather-Scatter map.
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.
bg::model::point< double, 3, bg::cs::cartesian > point
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.