Commit 0f998d37 authored by Dave Moxey's avatar Dave Moxey
Browse files

Merge remote-tracking branch 'upstream/master' into fix/cfs-dgopt

parents e8357f01 480c850d
......@@ -3,22 +3,35 @@ Changelog
v4.4.0
------
**Library:**
- Add support for variable polynomial order for 3D simulations with continuous
Galerkin discretisation (!604)
**IncNavierStokesSolver:**
- Add ability to simulate additional scalar fields (!624)
**NekMesh:**
- Modify curve module to allow for spline input (!628)
**FieldConvert:**
- Add module to stretch homogeneous direction (!609)
v4.3.3
------
**Library**:
- Fix filters when using adaptive driver to avoid output being overwritten after
each adaptive update (!588)
- Minor fix to suppress Xxt output unless `--verbose` is specified (!642)
- Fix of DirectFull solver in case where only Neumann boundary conditions
are imposed. (!655)
**FieldConvert**:
- Fix to avoid repeated import of field file (!649)
- Fix issue with C^0 projection (!644)
**CompressibleFlowSolver**:
- Fix issue with residual output (!647)
- Issues with 1D Euler solver fixed (!565)
- Fix deadlocking issue with boundary conditions (!657)
**Packaging**:
......
......@@ -102,6 +102,7 @@ possibly also Reynolds stresses) into single file;
\item \inltt{equispacedoutput}: Write data as equi-spaced output using simplices to represent the data for connecting points;
\item \inltt{extract}: Extract a boundary field;
\item \inltt{homplane}: Extract a plane from 3DH1D expansions;
\item \inltt{homstretch}: Stretch a 3DH1D expansion by an integer factor;
\item \inltt{innerproduct}: take the inner product between one or a series of fields with another field (or series of fields).
\item \inltt{interpfield}: Interpolates one field to another, requires fromxml, fromfld to be defined;
\item \inltt{interppointdatatofld}: Interpolates given discrete data using a finite difference approximation to a fld file given an xml file;
......@@ -327,6 +328,21 @@ The output file \inltt{file-plane.fld} can be processed in a similar
way as described in section \ref{s:utilities:fieldconvert:sub:convert}
to visualise it either in Tecplot or in Paraview.
\subsection{Stretch a 3DH1D expansion: \textit{homstretch} module}
To stretch a 3DH1D expansion in the z-direction, use the command:
\begin{lstlisting}[style=BashInputStyle]
FieldConvert -m homstretch:factor=value file.xml file.fld file-stretch.fld
\end{lstlisting}
The number of modes in the resulting field can be chosen using the command-line
parameter \inltt{output-points-hom-z}. Note that the output for
this module should always be a \inltt{.fld} file and this should not
be used in combination with other modules using a single command.
The output file \inltt{file-stretch.fld} can be processed in a similar
way as described in section \ref{s:utilities:fieldconvert:sub:convert}
to visualise it either in Tecplot or in Paraview.
\subsection{Inner Product of a single or series of fields with respect to a single or series of fields: \textit{innerproduct} module}
You can take the inner product of one field with another field using
......
......@@ -38,6 +38,7 @@
#include <LocalRegions/Expansion.h>
#include <LocalRegions/Expansion2D.h>
#include <LocalRegions/Expansion3D.h>
#include <LibUtilities/BasicUtils/ShapeType.hpp>
#include <boost/config.hpp>
......@@ -832,8 +833,8 @@ namespace Nektar
for(j = 0; j < nEdges; ++j)
{
nEdgeIntCoeffs = exp->GetEdgeNcoeffs(j) - 2;
meshEdgeId = exp->GetGeom()->GetEid(j);
nEdgeIntCoeffs = EdgeSize[meshEdgeId];
if(graph[1].count(meshEdgeId) == 0)
{
if(tempGraph[1].count(meshEdgeId) == 0)
......@@ -1301,12 +1302,12 @@ namespace Nektar
: AssemblyMap(pSession, variable)
{
int i, j, k, l;
int p, q, numModes0, numModes1;
int cnt = 0;
int intDofCnt;
int meshVertId, meshEdgeId, meshFaceId;
int globalId;
int nEdgeInteriorCoeffs;
int nFaceInteriorCoeffs;
int firstNonDirGraphVertId;
LibUtilities::CommSharedPtr vComm = m_comm->GetRowComm();
LocalRegions::ExpansionSharedPtr exp, bndExp;
......@@ -1326,6 +1327,8 @@ namespace Nektar
// Stores vertex, edge and face reordered vertices.
DofGraph graph(3);
DofGraph dofs(3);
vector<map<int, int> > faceModes(2);
map<int, LibUtilities::ShapeType> faceType;
set<int> extraDirVerts, extraDirEdges;
BottomUpSubStructuredGraphSharedPtr bottomUpGraph;
......@@ -1336,25 +1339,26 @@ namespace Nektar
{
exp = locExpVector[i];
for(j = 0; j < locExpVector[i]->GetNverts(); ++j)
for(j = 0; j < exp->GetNverts(); ++j)
{
dofs[0][exp->GetGeom()->GetVid(j)] = 1;
}
for(j = 0; j < locExpVector[i]->GetNedges(); ++j)
for(j = 0; j < exp->GetNedges(); ++j)
{
if (dofs[1].count(exp->GetGeom()->GetEid(j)) > 0)
{
if (dofs[1][exp->GetGeom()->GetEid(j)] !=
locExpVector[i]->GetEdgeNcoeffs(j)-2)
exp->GetEdgeNcoeffs(j)-2)
{
ASSERTL0( (exp->GetEdgeBasisType(j) == LibUtilities::eModified_A) ||
(exp->GetEdgeBasisType(j) == LibUtilities::eModified_B),
(exp->GetEdgeBasisType(j) == LibUtilities::eModified_B) ||
(exp->GetEdgeBasisType(j) == LibUtilities::eModified_C),
"CG with variable order only available with modal expansion");
}
dofs[1][exp->GetGeom()->GetEid(j)] =
min(dofs[1][exp->GetGeom()->GetEid(j)],
locExpVector[i]->GetEdgeNcoeffs(j)-2);
exp->GetEdgeNcoeffs(j)-2);
}
else
{
......@@ -1363,30 +1367,37 @@ namespace Nektar
}
}
for(j = 0; j < locExpVector[i]->GetNfaces(); ++j)
for(j = 0; j < exp->GetNfaces(); ++j)
{
if (dofs[2].count(exp->GetGeom()->GetFid(j)) > 0)
faceOrient = exp->GetGeom()->GetForient(j);
meshFaceId = exp->GetGeom()->GetFid(j);
exp->GetFaceNumModes(j, faceOrient, numModes0, numModes1);
if (faceModes[0].count(meshFaceId) > 0)
{
if (dofs[2][exp->GetGeom()->GetFid(j)] !=
exp->GetFaceIntNcoeffs(j))
{
ASSERTL0( false,
"CG with variable order not available in 3D");
}
dofs[2][exp->GetGeom()->GetFid(j)] =
min(dofs[2][exp->GetGeom()->GetFid(j)],
exp->GetFaceIntNcoeffs(j));
faceModes[0][meshFaceId] =
min(faceModes[0][meshFaceId], numModes0);
faceModes[1][meshFaceId] =
min(faceModes[1][meshFaceId], numModes1);
}
else
{
dofs[2][exp->GetGeom()->GetFid(j)] =
exp->GetFaceIntNcoeffs(j);
faceModes[0][meshFaceId] = numModes0;
faceModes[1][meshFaceId] = numModes1;
// Get shape of this face
SpatialDomains::Geometry3DSharedPtr geom;
geom = boost::dynamic_pointer_cast<SpatialDomains::
Geometry3D> (exp->GetGeom());
faceType[meshFaceId] =
geom->GetFace(j)->GetShapeType();
}
}
}
// Now use information from all partitions to determine
// the correct size
map<int, int>::iterator dofIt;
map<int, int>::iterator dofIt, dofIt2;
// edges
Array<OneD, long> edgeId (dofs[1].size());
Array<OneD, NekDouble> edgeDof (dofs[1].size());
......@@ -1403,19 +1414,41 @@ namespace Nektar
dofs[1][edgeId[i]] = (int) (edgeDof[i]+0.5);
}
// faces
Array<OneD, long> faceId (dofs[2].size());
Array<OneD, NekDouble> faceDof (dofs[2].size());
for(dofIt = dofs[2].begin(), i=0; dofIt != dofs[2].end(); dofIt++, i++)
Array<OneD, long> faceId (faceModes[0].size());
Array<OneD, NekDouble> faceP (faceModes[0].size());
Array<OneD, NekDouble> faceQ (faceModes[0].size());
for(dofIt = faceModes[0].begin(), dofIt2 = faceModes[1].begin(),i=0;
dofIt != faceModes[0].end(); dofIt++, dofIt2++, i++)
{
faceId[i] = dofIt->first;
faceDof[i] = (NekDouble) dofIt->second;
faceP[i] = (NekDouble) dofIt->second;
faceQ[i] = (NekDouble) dofIt2->second;
}
Gs::gs_data *tmp2 = Gs::Init(faceId, vComm);
Gs::Gather(faceDof, Gs::gs_min, tmp2);
Gs::Gather(faceP, Gs::gs_min, tmp2);
Gs::Gather(faceQ, Gs::gs_min, tmp2);
Gs::Finalise(tmp2);
for (i=0; i < dofs[2].size(); i++)
int P, Q;
for (i=0; i < faceModes[0].size(); i++)
{
dofs[2][faceId[i]] = (int) (faceDof[i]+0.5);
faceModes[0][faceId[i]] = (int) (faceP[i]+0.5);
faceModes[1][faceId[i]] = (int) (faceQ[i]+0.5);
P = faceModes[0][faceId[i]];
Q = faceModes[1][faceId[i]];
if (faceType[faceId[i]] == LibUtilities::eQuadrilateral)
{
// Quad face
dofs[2][faceId[i]] =
LibUtilities::StdQuadData::getNumberOfCoefficients(P,Q) -
LibUtilities::StdQuadData::getNumberOfBndCoefficients(P,Q);
}
else
{
// Tri face
dofs[2][faceId[i]] =
LibUtilities::StdTriData::getNumberOfCoefficients(P,Q) -
LibUtilities::StdTriData::getNumberOfBndCoefficients(P,Q);
}
}
Array<OneD, const BndCond> bndCondVec(1, bndConditions);
......@@ -1478,7 +1511,6 @@ namespace Nektar
for(j = 0; j < exp->GetNfaces(); ++j)
{
nFaceInteriorCoeffs = exp->GetFaceIntNcoeffs(j);
meshFaceId = exp->GetGeom()->GetFid(j);
graphVertOffset[graph[2][meshFaceId]+1] = dofs[2][meshFaceId];
}
......@@ -1590,7 +1622,6 @@ namespace Nektar
for(j = 0; j < exp->GetNfaces(); ++j)
{
nFaceInteriorCoeffs = exp->GetFaceIntNcoeffs(j);
faceOrient = exp->GetGeom()->GetForient(j);
meshFaceId = exp->GetGeom()->GetFid(j);
......@@ -1605,29 +1636,81 @@ namespace Nektar
exp->GetFaceInteriorMap(j,faceOrient,faceInteriorMap,faceInteriorSign);
// Set the global DOF's for the interior modes of face j
for(k = 0; k < dofs[2][exp->GetGeom()->GetFid(j)]; ++k)
{
m_localToGlobalMap[cnt+faceInteriorMap[k]] =
graphVertOffset[graph[2][meshFaceId]]+k;
}
for(k = dofs[2][exp->GetGeom()->GetFid(j)]; k < nFaceInteriorCoeffs; ++k)
exp->GetFaceNumModes(j, faceOrient, numModes0, numModes1);
switch(faceType[meshFaceId])
{
m_localToGlobalMap[cnt+faceInteriorMap[k]] =
graphVertOffset[graph[2][meshFaceId]];
}
if(m_signChange)
case LibUtilities::eQuadrilateral:
{
for(k = 0; k < dofs[2][exp->GetGeom()->GetFid(j)]; ++k)
int kLoc=0;
k = 0;
for( q = 2; q < numModes1; q++)
{
m_localToGlobalSign[cnt+faceInteriorMap[k]] = (NekDouble) faceInteriorSign[k];
for( p = 2; p < numModes0; p++)
{
if( (p < faceModes[0][meshFaceId]) &&
(q < faceModes[1][meshFaceId]))
{
m_localToGlobalMap[cnt+faceInteriorMap[kLoc]] =
graphVertOffset[graph[2][meshFaceId]]+k;
if(m_signChange)
{
m_localToGlobalSign[cnt+faceInteriorMap[kLoc]] =
(NekDouble) faceInteriorSign[kLoc];
}
k++;
}
else
{
m_localToGlobalMap[cnt+faceInteriorMap[kLoc]] =
graphVertOffset[graph[2][meshFaceId]];
if(m_signChange)
{
m_localToGlobalSign[cnt+faceInteriorMap[kLoc]] = 0.0;
}
}
kLoc++;
}
}
for(k = dofs[2][exp->GetGeom()->GetFid(j)]; k < nFaceInteriorCoeffs; ++k)
}
break;
case LibUtilities::eTriangle:
{
int kLoc=0;
k = 0;
for( p = 2; p < numModes0; p++)
{
m_localToGlobalSign[cnt+faceInteriorMap[k]] = 0;
for( q = 1; q < numModes1-p; q++)
{
if( (p < faceModes[0][meshFaceId]) &&
(p+q < faceModes[1][meshFaceId]))
{
m_localToGlobalMap[cnt+faceInteriorMap[kLoc]] =
graphVertOffset[graph[2][meshFaceId]]+k;
if(m_signChange)
{
m_localToGlobalSign[cnt+faceInteriorMap[kLoc]] =
(NekDouble) faceInteriorSign[kLoc];
}
k++;
}
else
{
m_localToGlobalMap[cnt+faceInteriorMap[kLoc]] =
graphVertOffset[graph[2][meshFaceId]];
if(m_signChange)
{
m_localToGlobalSign[cnt+faceInteriorMap[kLoc]] = 0.0;
}
}
kLoc++;
}
}
}
break;
default:
ASSERTL0(false,"Shape not recognised");
break;
}
}
}
......@@ -1688,19 +1771,28 @@ namespace Nektar
bndExp->GetEdgeInteriorMap(
k,edgeOrient,edgeInteriorMap,edgeInteriorSign);
for(l = 0; l < nEdgeInteriorCoeffs; ++l)
for(l = 0; l < dofs[1][meshEdgeId]; ++l)
{
m_bndCondCoeffsToGlobalCoeffsMap[cnt+edgeInteriorMap[l]] =
graphVertOffset[graph[1][meshEdgeId]]+l;
}
for(l = dofs[1][meshEdgeId]; l < nEdgeInteriorCoeffs; ++l)
{
m_bndCondCoeffsToGlobalCoeffsMap[cnt+edgeInteriorMap[l]] =
graphVertOffset[graph[1][meshEdgeId]];
}
// Fill the sign vector if required
if(m_signChange)
{
for(l = 0; l < nEdgeInteriorCoeffs; ++l)
for(l = 0; l < dofs[1][meshEdgeId]; ++l)
{
m_bndCondCoeffsToGlobalCoeffsSign[cnt+edgeInteriorMap[l]] = (NekDouble) edgeInteriorSign[l];
}
for(l = dofs[1][meshEdgeId]; l < nEdgeInteriorCoeffs; ++l)
{
m_bndCondCoeffsToGlobalCoeffsSign[cnt+edgeInteriorMap[l]] = 0.0;
}
}
if (bndConditions[i]->GetBoundaryConditionType() !=
......@@ -1714,7 +1806,7 @@ namespace Nektar
foundExtraEdges.count(meshEdgeId) == 0 &&
nEdgeInteriorCoeffs > 0)
{
for(l = 0; l < nEdgeInteriorCoeffs; ++l)
for(l = 0; l < dofs[1][meshEdgeId]; ++l)
{
int loc = bndCondExp[i]->GetCoeff_Offset(j) +
edgeInteriorMap[l];
......@@ -1723,6 +1815,15 @@ namespace Nektar
ExtraDirDof t(loc, gid, edgeInteriorSign[l]);
m_extraDirDofs[i].push_back(t);
}
for(l = dofs[1][meshEdgeId]; l < nEdgeInteriorCoeffs; ++l)
{
int loc = bndCondExp[i]->GetCoeff_Offset(j) +
edgeInteriorMap[l];
int gid = graphVertOffset[
graph[1][meshEdgeId]];
ExtraDirDof t(loc, gid, 0.0);
m_extraDirDofs[i].push_back(t);
}
foundExtraEdges.insert(meshEdgeId);
}
}
......@@ -2007,7 +2108,7 @@ namespace Nektar
int maxIntDof = 0;
int dof = 0;
int cnt;
int i,j,k;
int i,j,k,l;
int meshVertId;
int meshEdgeId;
int meshFaceId;
......@@ -2108,14 +2209,22 @@ namespace Nektar
dof = exp->GetEdgeNcoeffs(j)-2;
// Set the global DOF's for the interior modes of edge j
// run backwards because of variable P case "ghost" modes
for(k = dof-1; k >= 0; --k)
// for varP, ignore modes with sign == 0
for(k = 0, l = 0; k < dof; ++k)
{
if (m_signChange)
{
if (m_localToGlobalSign[cnt+edgeInteriorMap[k]]==0)
{
continue;
}
}
vGlobalId = m_localToGlobalMap[cnt+edgeInteriorMap[k]];
m_globalToUniversalMap[vGlobalId]
= nVert + meshEdgeId * maxEdgeDof + k + 1;
= nVert + meshEdgeId * maxEdgeDof + l + 1;
m_globalToUniversalBndMap[vGlobalId]=m_globalToUniversalMap[vGlobalId];
maxBndGlobalId = (vGlobalId > maxBndGlobalId ? vGlobalId : maxBndGlobalId);
l++;
}
}
......@@ -2140,16 +2249,23 @@ namespace Nektar
exp->GetFaceInteriorMap(j,faceOrient,faceInteriorMap,faceInteriorSign);
dof = exp->GetFaceIntNcoeffs(j);
for(k = dof-1; k >= 0; --k)
for(k = 0, l = 0; k < dof; ++k)
{
if (m_signChange)
{
if (m_localToGlobalSign[cnt+faceInteriorMap[k]]==0)
{
continue;
}
}
vGlobalId = m_localToGlobalMap[cnt+faceInteriorMap[k]];
m_globalToUniversalMap[vGlobalId]
= nVert + nEdge*maxEdgeDof + meshFaceId * maxFaceDof
+ k + 1;
+ l + 1;
m_globalToUniversalBndMap[vGlobalId]=m_globalToUniversalMap[vGlobalId];
maxBndGlobalId = (vGlobalId > maxBndGlobalId ? vGlobalId : maxBndGlobalId);
l++;
}
}
......
......@@ -339,7 +339,14 @@ namespace Nektar
for(int j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); ++j)
{
sign = m_locToGloMap->GetBndCondCoeffsToGlobalCoeffsSign(bndcnt);
inout[map[bndcnt++]] = sign * coeffs[j];
if(sign)
{
inout[map[bndcnt++]] = sign * coeffs[j];
}
else
{
bndcnt++;
}
}
}
else
......@@ -465,7 +472,14 @@ namespace Nektar
{
sign = m_locToGloMap->GetBndCondCoeffsToGlobalCoeffsSign(
bndcnt);
tmp[bndMap[bndcnt++]] = sign * coeffs[j];
if (sign)
{
tmp[bndMap[bndcnt++]] = sign * coeffs[j];
}
else
{
bndcnt++;
}
}
}
else
......
......@@ -895,11 +895,6 @@ namespace Nektar
Array<OneD, NekDouble> &Fwd,
Array<OneD, NekDouble> &Bwd)
{
// Expansion casts
LocalRegions::Expansion1DSharedPtr exp1D;
LocalRegions::Expansion1DSharedPtr exp1DFirst;
LocalRegions::Expansion1DSharedPtr exp1DLast;
// Counter variables
int n, v;
......@@ -911,10 +906,7 @@ namespace Nektar
Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
&elmtToTrace = m_traceMap->GetElmtToTrace();
// Basis shared pointer
LibUtilities::BasisSharedPtr Basis;
// Set forward and backard state to zero
Vmath::Zero(Fwd.num_elements(), Fwd, 1);
Vmath::Zero(Bwd.num_elements(), Bwd, 1);
......@@ -924,12 +916,8 @@ namespace Nektar
// Loop on the elements
for (cnt = n = 0; n < nElements; ++n)
{
exp1D = (*m_exp)[n]->as<LocalRegions::Expansion1D>();
// Set the offset of each element
phys_offset = GetPhys_Offset(n);
Basis = (*m_exp)[n]->GetBasis(0);
for(v = 0; v < 2; ++v, ++cnt)
{
......@@ -1052,24 +1040,19 @@ namespace Nektar
Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
&elmtToTrace = m_traceMap->GetElmtToTrace();
// Basis shared pointer
LibUtilities::BasisSharedPtr Basis;
vector<bool> negatedFluxNormal = GetNegatedFluxNormal();
for (n = 0; n < GetExpSize(); ++n)
{
// Basis definition on each element
Basis = (*m_exp)[n]->GetBasis(0);
// Number of coefficients on each element
int e_ncoeffs = (*m_exp)[n]->GetNcoeffs();
offset = GetCoeff_Offset(n);
// Implementation for every points except Gauss points
if (Basis->GetBasisType() != LibUtilities::eGauss_Lagrange)
if ((*m_exp)[n]->GetBasis(0)->GetBasisType() !=
LibUtilities::eGauss_Lagrange)
{
t_offset = GetTrace()->GetCoeff_Offset(elmtToTrace[n][0]->GetElmtId());
if(negatedFluxNormal[2*n])
......
......@@ -65,7 +65,7 @@ namespace Nektar
const Array<OneD,const NekDouble> &pInput,
Array<OneD, NekDouble> &pOutput,
const AssemblyMapSharedPtr &locToGloMap,
const int pNumDir = 0);
const int pNumDir);
};
}
}
......
......@@ -146,7 +146,7 @@ namespace Nektar
}
else
{
SolveLinearSystem(nDirDofs, pInput, pOutput, pLocToGloMap);
SolveLinearSystem(nGlobDofs, pInput, pOutput, pLocToGloMap, nDirDofs);
}
}
......