Commit c049aa42 authored by Gianmarco Mengaldo's avatar Gianmarco Mengaldo
Browse files

Added interpolation routines for FR using Gauss points. Needs to be improved.


git-svn-id: https://gforge.sci.utah.edu/svn/nektar/trunk@4042 305cdda6-5ce1-45b3-a98d-dfc68c8b3305
parent abceaded
......@@ -728,45 +728,60 @@ namespace Nektar
* @param Bwd Backward state.
*/
void DisContField1D::v_GetFwdBwdTracePhys(const Array<OneD, const NekDouble> &field,
Array<OneD, NekDouble> &Fwd,
Array<OneD, NekDouble> &Bwd)
Array<OneD, NekDouble> &Fwd,
Array<OneD, NekDouble> &Bwd)
{
/// Counter variables
// Counter variables
int n, p, i;
/// Number of elements
// Number of elements
int nElements = GetExpSize();
/// Number of quadrature points of each element
int nLocalQuadraturePts;
// Number of solution points of each element
int nLocalSolutionPts;
/// Initial index of each element
// Initial index of each element
int phys_offset;
/// Index of each interface between adjacent elements
// Index of each interface between adjacent elements
int interface_offset;
/// Index in case of subdomains
// Index in case of subdomains
int subdomain_offset;
/// Coordinate of each standard element (-1, 1)
// Coordinate of each standard element (-1, 1)
NekDouble vertex_coord;
/// Set forward and backard state to zero
// 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);
/// Loop on the elements
for(n = 0; n < nElements; ++n)
// Loop on the elements
for (n = 0; n < nElements; ++n)
{
/// Set the offset of each element
// Set the offset of each element
phys_offset = GetPhys_Offset(n);
/// Set the number of quadrature points of each element
nLocalQuadraturePts = (*m_exp)[n]->GetNumPoints(0);
// Set the number of solution points of each element
nLocalSolutionPts = (*m_exp)[n]->GetNumPoints(0);
for(p = 0; p < 2; ++p)
// Temporary vector for interpolation routine
Array<OneD, NekDouble> tmp(nLocalSolutionPts, 0.0);
// Partition the field vector in local vectors
tmp = field + (n * nLocalSolutionPts);
// Basis definition on each element
Basis = (*m_exp)[n]->GetBasis(0);
for (p = 0; p < 2; ++p)
{
// Coordinate vector for interpolation routine
Array<OneD, NekDouble> interface_coord(3, 0.0);
vertex_coord = 0.0;
interface_offset = (*m_exp)[n]->GetGeom1D()->GetVid(p);
subdomain_offset = (*m_exp)[0]->GetGeom1D()->GetVid(0);
......@@ -777,25 +792,51 @@ namespace Nektar
vertex_coord += ((*m_exp)[n]->GetVertexNormal(p))[i][0];
}
if(vertex_coord >= 0.0)
// Set the x-coordinate of the standard interface point
interface_coord[0] = vertex_coord;
// Implementation for every points except Gauss points
if (Basis->GetPointsType() != LibUtilities::eGaussGaussLegendre)
{
Fwd[interface_offset] = field[phys_offset+nLocalQuadraturePts-1];
}
if(vertex_coord < 0.0)
if(vertex_coord >= 0.0)
{
Fwd[interface_offset] = field[phys_offset+nLocalSolutionPts-1];
}
if(vertex_coord < 0.0)
{
Bwd[interface_offset] = field[phys_offset];
}
}
// Implementation for Gauss points
// (it doesn't work for WeakDG)
else
{
Bwd[interface_offset] = field[phys_offset];
StdRegions::StdSegExp StdSeg(Basis->GetBasisKey());
if(vertex_coord >= 0.0)
{
Fwd[interface_offset] =
//(*m_exp)[n]->PhysEvaluate(interface_coord, tmp);
StdSeg.PhysEvaluate(interface_coord, tmp);
}
if(vertex_coord < 0.0)
{
Bwd[interface_offset] =
//(*m_exp)[n]->PhysEvaluate(interface_coord, tmp);
StdSeg.PhysEvaluate(interface_coord, tmp);
}
}
}
}
/// Fill boundary conditions into the missing elements
// Fill boundary conditions into the missing elements
int id1 = 0;
int id2 = 0;
int firstVertex = (*m_exp)[0]->GetGeom1D()->GetVid(0);
int lastVertex = (*m_exp)[nElements-1]->GetGeom1D()->GetVid(1);
Array<OneD, NekDouble> processed(m_bndCondExpansions.num_elements()+1, -1.0);
/// Bug temporary fixed for Periodic boundary conditions
// Bug temporary fixed for Periodic boundary conditions
if(SpatialDomains::ePeriodic)
{
int nTracePts = Fwd.num_elements();
......@@ -805,7 +846,7 @@ namespace Nektar
for(n = 0; n < m_bndCondExpansions.num_elements(); ++n)
{
/// Check if the current boundary condition belongs to the current subdomain
// Check if the current boundary condition belongs to the current subdomain
if((m_bndCondExpansions[n]->GetVertex()->GetVid() >= firstVertex)
&& (m_bndCondExpansions[n]->GetVertex()->GetVid() <= lastVertex)
&& (processed[n] != m_bndCondExpansions[n]->GetVertex()->GetVid()))
......
......@@ -63,5 +63,13 @@ namespace Nektar
{
v_Advect(nConvectiveFields, fields, advVel, inarray, outarray);
}
void Advection::InterpToInterface(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
Array<OneD, NekDouble> &total,
Array<OneD, NekDouble> &InterfaceValue)
{
v_InterpToInterface(fields, total, InterfaceValue);
}
}
}
......@@ -68,6 +68,11 @@ namespace Nektar
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray);
SOLVER_UTILS_EXPORT void InterpToInterface(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
Array<OneD, NekDouble> &total,
Array<OneD, NekDouble> &InterfaceValue);
template<typename FuncPointerT, typename ObjectPointerT>
void SetFluxVector(FuncPointerT func, ObjectPointerT obj)
{
......@@ -94,6 +99,14 @@ namespace Nektar
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray) = 0;
virtual void v_InterpToInterface(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
Array<OneD, NekDouble> &total,
Array<OneD, NekDouble> &InterfaceValue)
{
};
AdvectionFluxVecCB m_fluxVector;
RiemannSolverSharedPtr m_riemann;
};
......
......@@ -162,6 +162,9 @@ namespace Nektar
// Vector to store the solution in physical space
Array<OneD, Array<OneD, NekDouble> > physfield (nConvectiveFields);
// Vector to store the flux interpolated at the interfaces
Array<OneD, Array<OneD, NekDouble> > interpolatedFlux (nConvectiveFields);
// Resize each column of the flux vector to the number of
// solution points
......@@ -173,7 +176,8 @@ namespace Nektar
for(i = 0; i < nConvectiveFields; ++i)
{
physfield[i] = inarray[i];
physfield[i] = inarray[i];
interpolatedFlux[i] = Array<OneD, NekDouble>(2*nTracePts-2);
}
// Get the discontinuous flux FD
......@@ -200,6 +204,16 @@ namespace Nektar
}
}
// Interpolation routine for Gauss points
if (Basis->GetPointsType() == LibUtilities::eGaussGaussLegendre)
{
for (j = 0; j < nConvectiveFields; j++)
{
v_InterpToInterface(fields, fluxvector[j], interpolatedFlux[j]);
}
}
// Store forwards/backwards space along trace space.
Array<OneD, Array<OneD, NekDouble> > Fwd (nConvectiveFields);
Array<OneD, Array<OneD, NekDouble> > Bwd (nConvectiveFields);
......@@ -213,6 +227,7 @@ namespace Nektar
fields[i]->GetFwdBwdTracePhys(inarray[i], Fwd[i], Bwd[i]);
}
// Computing the Riemann flux at each flux (interface) point
m_riemann->Solve(Fwd, Bwd, numflux);
......@@ -252,6 +267,22 @@ namespace Nektar
numfluxjumpsRight[0][i] = numflux[0][i+1] - fluxvector[0][offsetEnd];
}
// Loop to compute the left and the right jumps of the flux
// with the new interpolation order used for Gauss points
if (Basis->GetPointsType() == LibUtilities::eGaussGaussLegendre)
{
j = 0;
for (i = 0; i < nElements; i++)
{
numfluxjumpsLeft[0][i] = numflux[0][i] - interpolatedFlux[0][j];
j++;
numfluxjumpsRight[0][i] = numflux[0][i+1] - interpolatedFlux[0][j];
j++;
}
}
for (i = 0; i < nElements; i++)
{
Vmath::Smul(nSolutionPts/nElements,
......@@ -311,5 +342,68 @@ namespace Nektar
}
}
}
void AdvectionFR::v_InterpToInterface(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
Array<OneD, NekDouble> &total,
Array<OneD, NekDouble> &InterfaceValue)
{
// Counter variables
int n, i, j, cnt, numero;
// Number of elements
int nElements = fields[0]->GetExpSize();
// Initial index of each element
int phys_offset;
// Number of quadrature points of each element
int nLocalSolutionPts;
// Number of vertices of each element
int nLocalVertices;
// Index of each interface between adjacent elements
int interface_offset;
cnt = 0;
for (n = 0; n < nElements; n++)
{
// Getting the basis
LibUtilities::BasisSharedPtr Basis;
Basis = fields[0]->GetExp(n)->GetBasis(0);
StdRegions::StdSegExp StdSeg(Basis->GetBasisKey());
// Set the offset of each element
phys_offset = fields[0]->GetPhys_Offset(n);
// Set the number of solution points of each element
nLocalSolutionPts = fields[0]->GetExp(n)->GetNumPoints(0);
// Set the number of vertices of each element
nLocalVertices = fields[0]->GetExp(n)->GetNverts();
Array<OneD, NekDouble> tmp(nLocalSolutionPts, 0.0);
tmp = total + n*nLocalSolutionPts;
for (i = 0; i < nLocalVertices; i++)
{
Array<OneD, NekDouble> interface_coord(3, 0.0);
interface_offset = fields[0]->GetExp(n)->GetGeom1D()->GetVid(i);
for (j = 0; j < (fields[0]->GetExp(n)->GetVertexNormal(i)).num_elements(); j++)
{
interface_coord[0] += (fields[0]->GetExp(n)->GetVertexNormal(i))[j][0];
//InterfaceValue[cnt] = fields[0]->GetExp(n)->PhysEvaluate(interface_coord, total);
InterfaceValue[cnt] = StdSeg.PhysEvaluate(interface_coord, tmp);
cnt++;
}
}
}
}// End v_InterpToInterface
}
}
......@@ -78,6 +78,11 @@ namespace Nektar
const Array<OneD, Array<OneD, NekDouble> > &advVel,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray);
virtual void v_InterpToInterface(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
Array<OneD, NekDouble> &total,
Array<OneD, NekDouble> &InterfaceValue);
};
}
}
......
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