Commit 7ada241e authored by Douglas Serson's avatar Douglas Serson Committed by Chris Cantwell

Fix wss module with partitioned domains

(cherry picked from commit 43ecd9b8)
parent fdce2f8d
......@@ -126,7 +126,7 @@ void ProcessWSS::Process(po::variables_map &vm)
Array<OneD, Array<OneD, NekDouble> > fshear(nshear);
Array<OneD, MultiRegions::ExpListSharedPtr> BndExp(newfields);
Array<OneD, MultiRegions::ExpListSharedPtr> BndElmtExp(nfields);
Array<OneD, MultiRegions::ExpListSharedPtr> BndElmtExp(spacedim);
// Extract original fields to boundary (for output)
for (int i = 0; i < m_f->m_exp.size(); ++i)
......@@ -155,144 +155,160 @@ void ProcessWSS::Process(po::variables_map &vm)
m_f->m_fielddef[0]->m_fields.push_back("Shear_mag");
}
// Create map of boundary ids for partitioned domains
SpatialDomains::BoundaryConditions bcs(m_f->m_session,
m_f->m_exp[0]->GetGraph());
const SpatialDomains::BoundaryRegionCollection bregions =
bcs.GetBoundaryRegions();
SpatialDomains::BoundaryRegionCollection::const_iterator breg_it;
map<int,int> BndRegionMap;
int cnt =0;
for(breg_it = bregions.begin(); breg_it != bregions.end();
++breg_it, ++cnt)
{
BndRegionMap[breg_it->first] = cnt;
}
// Loop over boundaries to Write
for(int b = 0; b < m_f->m_bndRegionsToWrite.size(); ++b)
{
int bnd = m_f->m_bndRegionsToWrite[b];
// Get expansion list for boundary and for elements containing this bnd
for(i = 0; i < newfields; i++)
{
BndExp[i] = m_f->m_exp[nfields + i]->UpdateBndCondExpansion(bnd);
}
for(i = 0; i < spacedim; i++)
if(BndRegionMap.count(m_f->m_bndRegionsToWrite[i]) == 1)
{
m_f->m_exp[i]->GetBndElmtExpansion(bnd, BndElmtExp[i]);
}
// Get number of points in expansions
int nqb = BndExp[0]->GetTotPoints();
int nqe = BndElmtExp[0]->GetTotPoints();
// Initialise local arrays for the velocity gradients, and stress components
// size of total number of quadrature points for elements in this bnd
for(i = 0; i < ngrad; ++i)
{
grad[i] = Array<OneD, NekDouble>(nqe);
}
int bnd = BndRegionMap[m_f->m_bndRegionsToWrite[i]];
// Get expansion list for boundary and for elements containing this bnd
for(i = 0; i < newfields; i++)
{
BndExp[i] = m_f->m_exp[nfields + i]->UpdateBndCondExpansion(bnd);
}
for(i = 0; i < spacedim; i++)
{
m_f->m_exp[i]->GetBndElmtExpansion(bnd, BndElmtExp[i]);
}
for(i = 0; i < nstress; ++i)
{
stress[i] = Array<OneD, NekDouble>(nqe);
}
// Get number of points in expansions
int nqb = BndExp[0]->GetTotPoints();
int nqe = BndElmtExp[0]->GetTotPoints();
// initialise arrays in the boundary
for(i = 0; i < nstress; ++i)
{
fstress[i] = Array<OneD, NekDouble>(nqb);
}
// Initialise local arrays for the velocity gradients, and stress components
// size of total number of quadrature points for elements in this bnd
for(i = 0; i < ngrad; ++i)
{
grad[i] = Array<OneD, NekDouble>(nqe);
}
for(i = 0; i < ngrad; ++i)
{
fgrad[i] = Array<OneD, NekDouble>(nqb);
}
for(i = 0; i < nstress; ++i)
{
stress[i] = Array<OneD, NekDouble>(nqe);
}
for(i = 0; i < nshear; ++i)
{
fshear[i] = Array<OneD, NekDouble>(nqb, 0.0);
}
// initialise arrays in the boundary
for(i = 0; i < nstress; ++i)
{
fstress[i] = Array<OneD, NekDouble>(nqb);
}
//Extract Velocities
for(i = 0; i < spacedim; ++i)
{
velocity[i] = BndElmtExp[i]->GetPhys();
}
for(i = 0; i < ngrad; ++i)
{
fgrad[i] = Array<OneD, NekDouble>(nqb);
}
//Compute gradients (velocity correction scheme method)
for(i = 0; i < spacedim; ++i)
{
if (spacedim == 2)
for(i = 0; i < nshear; ++i)
{
BndElmtExp[i]->PhysDeriv(velocity[i],grad[i*spacedim+0],
grad[i*spacedim+1]);
fshear[i] = Array<OneD, NekDouble>(nqb, 0.0);
}
else
//Extract Velocities
for(i = 0; i < spacedim; ++i)
{
BndElmtExp[i]->PhysDeriv(velocity[i],grad[i*spacedim+0],
grad[i*spacedim+1],
grad[i*spacedim+2]);
velocity[i] = BndElmtExp[i]->GetPhys();
}
}
//Compute stress component terms tau_ij = mu*(u_i,j + u_j,i)
for(i = 0; i < spacedim; ++i)
{
for(j = 0; j < spacedim; ++j)
//Compute gradients (velocity correction scheme method)
for(i = 0; i < spacedim; ++i)
{
Vmath::Vadd(nqe, grad[i*spacedim+j], 1,
grad[j*spacedim+i], 1,
stress[i*spacedim+j], 1);
Vmath::Smul(nqe, kinvis, stress[i*spacedim+j], 1,
stress[i*spacedim+j], 1);
if (spacedim == 2)
{
BndElmtExp[i]->PhysDeriv(velocity[i],grad[i*spacedim+0],
grad[i*spacedim+1]);
}
else
{
BndElmtExp[i]->PhysDeriv(velocity[i],grad[i*spacedim+0],
grad[i*spacedim+1],
grad[i*spacedim+2]);
}
}
}
// Get boundary stress values.
for(j = 0; j < nstress; ++j)
{
m_f->m_exp[0]->ExtractElmtToBndPhys(bnd, stress[j],fstress[j]);
}
//Get normals
Array<OneD, Array<OneD, NekDouble> > normals;
m_f->m_exp[0]->GetBoundaryNormals(bnd, normals);
// Reverse normals, to get correct orientation for the body
for(i = 0; i < spacedim; ++i)
{
Vmath::Neg(nqb, normals[i], 1);
}
//Compute stress component terms tau_ij = mu*(u_i,j + u_j,i)
for(i = 0; i < spacedim; ++i)
{
for(j = 0; j < spacedim; ++j)
{
Vmath::Vadd(nqe, grad[i*spacedim+j], 1,
grad[j*spacedim+i], 1,
stress[i*spacedim+j], 1);
Vmath::Smul(nqe, kinvis, stress[i*spacedim+j], 1,
stress[i*spacedim+j], 1);
}
}
//calculate wss, and update coeffs in the boundary expansion
// S = tau_ij * n_j
for(i = 0; i < spacedim; ++i)
{
for(j = 0; j < spacedim; ++j)
// Get boundary stress values.
for(j = 0; j < nstress; ++j)
{
Vmath::Vvtvp(nqb,normals[j],1,fstress[i*spacedim+j],1,
fshear[i],1,
fshear[i],1);
m_f->m_exp[0]->ExtractElmtToBndPhys(bnd, stress[j],fstress[j]);
}
}
// T = S - (S.n)n
for(i = 0; i < spacedim; ++i)
{
Vmath::Vvtvp(nqb,normals[i],1,fshear[i],1,
fshear[nshear-1],1,
fshear[nshear-1],1);
}
Vmath::Smul(nqb, -1.0, fshear[nshear-1], 1, fshear[nshear-1], 1);
//Get normals
Array<OneD, Array<OneD, NekDouble> > normals;
m_f->m_exp[0]->GetBoundaryNormals(bnd, normals);
// Reverse normals, to get correct orientation for the body
for(i = 0; i < spacedim; ++i)
{
Vmath::Neg(nqb, normals[i], 1);
}
for (i = 0; i < spacedim; i++)
{
Vmath::Vvtvp(nqb,normals[i], 1, fshear[nshear-1], 1,
fshear[i], 1,
fshear[i], 1);
BndExp[i]->FwdTrans(fshear[i],
BndExp[i]->UpdateCoeffs());
}
//calculate wss, and update coeffs in the boundary expansion
// S = tau_ij * n_j
for(i = 0; i < spacedim; ++i)
{
for(j = 0; j < spacedim; ++j)
{
Vmath::Vvtvp(nqb,normals[j],1,fstress[i*spacedim+j],1,
fshear[i],1,
fshear[i],1);
}
}
// Tw
Vmath::Zero(nqb, fshear[nshear-1], 1);
for(i = 0; i < spacedim; ++i)
{
Vmath::Vvtvp(nqb,fshear[i],1,fshear[i],1,
fshear[nshear-1],1,
fshear[nshear-1],1);
// T = S - (S.n)n
for(i = 0; i < spacedim; ++i)
{
Vmath::Vvtvp(nqb,normals[i],1,fshear[i],1,
fshear[nshear-1],1,
fshear[nshear-1],1);
}
Vmath::Smul(nqb, -1.0, fshear[nshear-1], 1, fshear[nshear-1], 1);
for (i = 0; i < spacedim; i++)
{
Vmath::Vvtvp(nqb,normals[i], 1, fshear[nshear-1], 1,
fshear[i], 1,
fshear[i], 1);
BndExp[i]->FwdTrans(fshear[i],
BndExp[i]->UpdateCoeffs());
}
// Tw
Vmath::Zero(nqb, fshear[nshear-1], 1);
for(i = 0; i < spacedim; ++i)
{
Vmath::Vvtvp(nqb,fshear[i],1,fshear[i],1,
fshear[nshear-1],1,
fshear[nshear-1],1);
}
Vmath::Vsqrt(nqb, fshear[nshear-1], 1, fshear[nshear-1], 1);
BndExp[nshear-1]->FwdTrans(fshear[nshear-1],
BndExp[nshear-1]->UpdateCoeffs());
}
Vmath::Vsqrt(nqb, fshear[nshear-1], 1, fshear[nshear-1], 1);
BndExp[nshear-1]->FwdTrans(fshear[nshear-1],
BndExp[nshear-1]->UpdateCoeffs());
}
}
......
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