Commit da2e0b0e authored by Douglas Serson's avatar Douglas Serson

Output extra plane in Homogeneous1D to fill the domain

parent f6157317
......@@ -367,9 +367,24 @@ namespace Nektar
int expansion,
int istrip)
{
// If there is only one plane (e.g. HalfMode), we write a 2D plane.
if (m_planes.num_elements() == 1)
{
m_planes[0]->WriteVtkPieceHeader(outfile, expansion);
return;
}
// If we are using Fourier points, output extra plane to fill domain
int outputExtraPlane = 0;
if ( m_homogeneousBasis->GetBasisType() == LibUtilities::eFourier
&& m_homogeneousBasis->GetPointsType() ==
LibUtilities::eFourierEvenlySpaced)
{
outputExtraPlane = 1;
}
int i, j;
int nquad0 = (*m_exp)[expansion]->GetNumPoints(0);
int nquad1 = m_planes.num_elements();
int nquad1 = m_planes.num_elements() + outputExtraPlane;
int ntot = nquad0*nquad1;
int ntotminus = (nquad0-1)*(nquad1-1);
......@@ -379,6 +394,20 @@ namespace Nektar
coords[2] = Array<OneD,NekDouble>(ntot);
GetCoords(expansion,coords[0],coords[1],coords[2]);
if (outputExtraPlane)
{
// Copy coords[0] and coords[1] to extra plane
Array<OneD,NekDouble> tmp;
Vmath::Vcopy (nquad0, coords[0], 1,
tmp = coords[0] + (nquad1-1)*nquad0, 1);
Vmath::Vcopy (nquad0, coords[1], 1,
tmp = coords[1] + (nquad1-1)*nquad0, 1);
// Fill coords[2] for extra plane
NekDouble z = coords[2][nquad0*m_planes.num_elements()-1] +
(coords[2][nquad0] - coords[2][0]);
Vmath::Fill(nquad0, z, tmp = coords[2] + (nquad1-1)*nquad0, 1);
}
outfile << " <Piece NumberOfPoints=\""
<< ntot << "\" NumberOfCells=\""
<< ntotminus << "\">" << endl;
......
......@@ -371,10 +371,19 @@ namespace Nektar
return;
}
// If we are using Fourier points, output extra plane to fill domain
int outputExtraPlane = 0;
if ( m_homogeneousBasis->GetBasisType() == LibUtilities::eFourier
&& m_homogeneousBasis->GetPointsType() ==
LibUtilities::eFourierEvenlySpaced)
{
outputExtraPlane = 1;
}
int i,j,k;
int nq0 = (*m_exp)[expansion]->GetNumPoints(0);
int nq1 = (*m_exp)[expansion]->GetNumPoints(1);
int nq2 = m_planes.num_elements();
int nq2 = m_planes.num_elements() + outputExtraPlane;
int ntot = nq0*nq1*nq2;
int ntotminus = (nq0-1)*(nq1-1)*(nq2-1);
......@@ -384,6 +393,20 @@ namespace Nektar
coords[2] = Array<OneD,NekDouble>(ntot);
GetCoords(expansion,coords[0],coords[1],coords[2]);
if (outputExtraPlane)
{
// Copy coords[0] and coords[1] to extra plane
Array<OneD,NekDouble> tmp;
Vmath::Vcopy (nq0*nq1, coords[0], 1,
tmp = coords[0] + (nq2-1)*nq0*nq1, 1);
Vmath::Vcopy (nq0*nq1, coords[1], 1,
tmp = coords[1] + (nq2-1)*nq0*nq1, 1);
// Fill coords[2] for extra plane
NekDouble z = coords[2][nq0*nq1*m_planes.num_elements()-1] +
(coords[2][nq0*nq1] - coords[2][0]);
Vmath::Fill(nq0*nq1, z, tmp = coords[2] + (nq2-1)*nq0*nq1, 1);
}
NekDouble DistStrip;
m_session->LoadParameter("DistStrip", DistStrip, 0);
// Reset the z-coords for homostrips
......
......@@ -880,6 +880,35 @@ namespace Nektar
int nq = (*m_exp)[expansion]->GetTotPoints();
int npoints_per_plane = m_planes[0]->GetTotPoints();
// If we are using Fourier points, output extra plane to fill domain
int outputExtraPlane = 0;
Array<OneD, NekDouble> extraPlane;
if ( m_homogeneousBasis->GetBasisType() == LibUtilities::eFourier
&& m_homogeneousBasis->GetPointsType() ==
LibUtilities::eFourierEvenlySpaced)
{
outputExtraPlane = 1;
// Get extra plane data
if (m_StripZcomm->GetSize() == 1)
{
extraPlane = m_phys + m_phys_offset[expansion];
}
else
{
// Determine to and from rank for communication
int size = m_StripZcomm->GetSize();
int rank = m_StripZcomm->GetRank();
int fromRank = (rank+1) % size;
int toRank = (rank == 0) ? size-1 : rank-1;
// Communicate using SendRecv
extraPlane = Array<OneD, NekDouble>(nq);
Array<OneD, NekDouble> send (nq,
m_phys + m_phys_offset[expansion]);
m_StripZcomm->SendRecv(toRank, send,
fromRank, extraPlane);
}
}
// printing the fields of that zone
outfile << " <DataArray type=\"Float64\" Name=\""
<< var << "\">" << endl;
......@@ -892,6 +921,14 @@ namespace Nektar
outfile << (fabs(phys[i]) < NekConstants::kNekZeroTol ? 0 : phys[i]) << " ";
}
}
if (outputExtraPlane)
{
for(i = 0; i < nq; ++i)
{
outfile << (fabs(extraPlane[i]) < NekConstants::kNekZeroTol ?
0 : extraPlane[i]) << " ";
}
}
outfile << endl;
outfile << " </DataArray>" << endl;
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment