Commit 526ba95d authored by Dave Moxey's avatar Dave Moxey
Browse files

Merge branch 'fix/GetFacePhys' into 'master'

Fix/get face phys

Fix to use correct memory when interpolating faces in tets, prisms and pyramids.

See merge request !393
parents e6dd7e22 2aa9b85e
This diff is collapsed.
......@@ -580,7 +580,7 @@ namespace Nektar
int nquad1 = m_base[1]->GetNumPoints();
int nquad2 = m_base[2]->GetNumPoints();
Array<OneD,NekDouble> o_tmp(nquad0*nquad1*nquad2);
Array<OneD,NekDouble> o_tmp(GetFaceNumPoints(face));
if (orient == StdRegions::eNoOrientation)
{
......@@ -593,14 +593,14 @@ namespace Nektar
if(orient == StdRegions::eDir1FwdDir1_Dir2FwdDir2)
{
//Directions A and B positive
Vmath::Vcopy(nquad0*nquad1,&(inarray[0]),1,&(outarray[0]),1);
Vmath::Vcopy(nquad0*nquad1,&(inarray[0]),1,&(o_tmp[0]),1);
}
else if(orient == StdRegions::eDir1BwdDir1_Dir2FwdDir2)
{
//Direction A negative and B positive
for (int j=0; j<nquad1; j++)
{
Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0-1)+j*nquad0,-1,&(outarray[0])+(j*nquad0),1);
Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0-1)+j*nquad0,-1,&(o_tmp[0])+(j*nquad0),1);
}
}
else if(orient == StdRegions::eDir1FwdDir1_Dir2BwdDir2)
......@@ -608,7 +608,7 @@ namespace Nektar
//Direction A positive and B negative
for (int j=0; j<nquad1; j++)
{
Vmath::Vcopy(nquad0,&(inarray[0])+nquad0*(nquad1-1-j),1,&(outarray[0])+(j*nquad0),1);
Vmath::Vcopy(nquad0,&(inarray[0])+nquad0*(nquad1-1-j),1,&(o_tmp[0])+(j*nquad0),1);
}
}
else if(orient == StdRegions::eDir1BwdDir1_Dir2BwdDir2)
......@@ -616,7 +616,7 @@ namespace Nektar
//Direction A negative and B negative
for(int j=0; j<nquad1; j++)
{
Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0*nquad1-1-j*nquad0),-1,&(outarray[0])+(j*nquad0),1);
Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0*nquad1-1-j*nquad0),-1,&(o_tmp[0])+(j*nquad0),1);
}
}
else if(orient == StdRegions::eDir1FwdDir2_Dir2FwdDir1)
......@@ -624,7 +624,7 @@ namespace Nektar
//Transposed, Direction A and B positive
for (int i=0; i<nquad0; i++)
{
Vmath::Vcopy(nquad1,&(inarray[0])+i,nquad0,&(outarray[0])+(i*nquad1),1);
Vmath::Vcopy(nquad1,&(inarray[0])+i,nquad0,&(o_tmp[0])+(i*nquad1),1);
}
}
else if(orient == StdRegions::eDir1FwdDir2_Dir2BwdDir1)
......@@ -632,7 +632,7 @@ namespace Nektar
//Transposed, Direction A positive and B negative
for (int i=0; i<nquad0; i++)
{
Vmath::Vcopy(nquad1,&(inarray[0])+(nquad0-1-i),nquad0,&(outarray[0])+(i*nquad1),1);
Vmath::Vcopy(nquad1,&(inarray[0])+(nquad0-1-i),nquad0,&(o_tmp[0])+(i*nquad1),1);
}
}
else if(orient == StdRegions::eDir1BwdDir2_Dir2FwdDir1)
......@@ -640,7 +640,7 @@ namespace Nektar
//Transposed, Direction A negative and B positive
for (int i=0; i<nquad0; i++)
{
Vmath::Vcopy(nquad1,&(inarray[0])+i+nquad0*(nquad1-1),-nquad0,&(outarray[0])+(i*nquad1),1);
Vmath::Vcopy(nquad1,&(inarray[0])+i+nquad0*(nquad1-1),-nquad0,&(o_tmp[0])+(i*nquad1),1);
}
}
else if(orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1)
......@@ -648,10 +648,9 @@ namespace Nektar
//Transposed, Direction A and B negative
for (int i=0; i<nquad0; i++)
{
Vmath::Vcopy(nquad1,&(inarray[0])+(nquad0*nquad1-1-i),-nquad0,&(outarray[0])+(i*nquad1),1);
Vmath::Vcopy(nquad1,&(inarray[0])+(nquad0*nquad1-1-i),-nquad0,&(o_tmp[0])+(i*nquad1),1);
}
}
o_tmp=outarray;
//interpolate
LibUtilities::Interp2D(m_base[0]->GetPointsKey(), m_base[1]->GetPointsKey(), o_tmp,
FaceExp->GetBasis(0)->GetPointsKey(),FaceExp->GetBasis(1)->GetPointsKey(),outarray);
......@@ -662,18 +661,16 @@ namespace Nektar
//Direction A and B positive
for (int k=0; k<nquad2; k++)
{
Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0*nquad1*k),1,&(outarray[0])+(k*nquad0),1);
Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0*nquad1*k),1,&(o_tmp[0])+(k*nquad0),1);
}
o_tmp=outarray;
}
else
{
//Direction A negative and B positive
for (int k=0; k<nquad2; k++)
{
Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0-1)+(nquad0*nquad1*k),-1,&(outarray[0])+(k*nquad0),1);
Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0-1)+(nquad0*nquad1*k),-1,&(o_tmp[0])+(k*nquad0),1);
}
o_tmp=outarray;
}
//interpolate
......@@ -685,7 +682,7 @@ namespace Nektar
{
//Directions A and B positive
Vmath::Vcopy(nquad0*nquad2,&(inarray[0])+(nquad0-1),
nquad0,&(outarray[0]),1);
nquad0,&(o_tmp[0]),1);
}
else if(orient == StdRegions::eDir1BwdDir1_Dir2FwdDir2)
{
......@@ -693,7 +690,7 @@ namespace Nektar
for (int k=0; k<nquad2; k++)
{
Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0*nquad1-1)+(k*nquad0*nquad1),
-nquad0,&(outarray[0])+(k*nquad0),1);
-nquad0,&(o_tmp[0])+(k*nquad0),1);
}
}
else if(orient == StdRegions::eDir1FwdDir1_Dir2BwdDir2)
......@@ -702,7 +699,7 @@ namespace Nektar
for (int k=0; k<nquad2; k++)
{
Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0-1)+(nquad0*nquad1*(nquad2-1-k)),
nquad0,&(outarray[0])+(k*nquad0),1);
nquad0,&(o_tmp[0])+(k*nquad0),1);
}
}
else if(orient == StdRegions::eDir1BwdDir1_Dir2BwdDir2)
......@@ -711,7 +708,7 @@ namespace Nektar
for (int k=0; k<nquad2; k++)
{
Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0*nquad1-1)+(nquad0*nquad1*(nquad2-1-k)),
-nquad0,&(outarray[0])+(k*nquad0),1);
-nquad0,&(o_tmp[0])+(k*nquad0),1);
}
}
else if(orient == StdRegions::eDir1FwdDir2_Dir2FwdDir1)
......@@ -720,7 +717,7 @@ namespace Nektar
for (int j=0; j<nquad1; j++)
{
Vmath::Vcopy(nquad2,&(inarray[0])+(nquad0-1)+(j*nquad0),
nquad0*nquad1,&(outarray[0])+(j*nquad2),1);
nquad0*nquad1,&(o_tmp[0])+(j*nquad2),1);
}
}
else if(orient == StdRegions::eDir1FwdDir2_Dir2BwdDir1)
......@@ -729,7 +726,7 @@ namespace Nektar
for (int j=0; j<nquad0; j++)
{
Vmath::Vcopy(nquad2,&(inarray[0])+(nquad0*nquad1-1-j*nquad0),
nquad0*nquad1,&(outarray[0])+(j*nquad2),1);
nquad0*nquad1,&(o_tmp[0])+(j*nquad2),1);
}
}
else if(orient == StdRegions::eDir1BwdDir2_Dir2FwdDir1)
......@@ -738,7 +735,7 @@ namespace Nektar
for (int j=0; j<nquad0; j++)
{
Vmath::Vcopy(nquad2,&(inarray[0])+nquad0*nquad1*(nquad2-1)+nquad0+j*nquad0,
-nquad0*nquad1,&(outarray[0])+(j*nquad2),1);
-nquad0*nquad1,&(o_tmp[0])+(j*nquad2),1);
}
}
else if(orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1)
......@@ -747,10 +744,9 @@ namespace Nektar
for (int j=0; j<nquad0; j++)
{
Vmath::Vcopy(nquad2,&(inarray[0])+(nquad0*nquad1*nquad2-1-j*nquad0),
-nquad0*nquad1,&(outarray[0])+(j*nquad2),1);
-nquad0*nquad1,&(o_tmp[0])+(j*nquad2),1);
}
}
o_tmp=outarray;
//interpolate
LibUtilities::Interp2D(m_base[1]->GetPointsKey(), m_base[2]->GetPointsKey(), o_tmp,
FaceExp->GetBasis(0)->GetPointsKey(),FaceExp->GetBasis(1)->GetPointsKey(),outarray);
......@@ -762,7 +758,7 @@ namespace Nektar
for (int k=0; k<nquad2; k++)
{
Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0*(nquad1-1))+(k*nquad0*nquad1),
1,&(outarray[0])+(k*nquad0),1);
1,&(o_tmp[0])+(k*nquad0),1);
}
}
else
......@@ -771,10 +767,9 @@ namespace Nektar
for (int k=0; k<nquad2; k++)
{
Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0*nquad1-1)+(k*nquad0*nquad1),
-1,&(outarray[0])+(k*nquad0),1);
-1,&(o_tmp[0])+(k*nquad0),1);
}
}
o_tmp=outarray;
//interpolate
LibUtilities::Interp2D(m_base[0]->GetPointsKey(), m_base[2]->GetPointsKey(), o_tmp,
FaceExp->GetBasis(0)->GetPointsKey(),FaceExp->GetBasis(1)->GetPointsKey(),outarray);
......@@ -784,7 +779,7 @@ namespace Nektar
if(orient == StdRegions::eDir1FwdDir1_Dir2FwdDir2)
{
//Directions A and B positive
Vmath::Vcopy(nquad1*nquad2,&(inarray[0]),nquad0,&(outarray[0]),1);
Vmath::Vcopy(nquad1*nquad2,&(inarray[0]),nquad0,&(o_tmp[0]),1);
}
else if(orient == StdRegions::eDir1BwdDir1_Dir2FwdDir2)
{
......@@ -792,7 +787,7 @@ namespace Nektar
for (int k=0; k<nquad2; k++)
{
Vmath::Vcopy(nquad1,&(inarray[0])+nquad0*(nquad1-1)+(k*nquad0*nquad1),
-nquad0,&(outarray[0])+(k*nquad1),1);
-nquad0,&(o_tmp[0])+(k*nquad1),1);
}
}
else if(orient == StdRegions::eDir1FwdDir1_Dir2BwdDir2)
......@@ -801,7 +796,7 @@ namespace Nektar
for (int k=0; k<nquad2; k++)
{
Vmath::Vcopy(nquad1,&(inarray[0])+(nquad0*nquad1*(nquad2-1-k)),
nquad0,&(outarray[0])+(k*nquad1),1);
nquad0,&(o_tmp[0])+(k*nquad1),1);
}
}
else if(orient == StdRegions::eDir1BwdDir1_Dir2BwdDir2)
......@@ -810,7 +805,7 @@ namespace Nektar
for (int k=0; k<nquad2; k++)
{
Vmath::Vcopy(nquad1,&(inarray[0])+nquad0*(nquad1-1)+(nquad0*nquad1*(nquad2-1-k)),
-nquad0,&(outarray[0])+(k*nquad1),1);
-nquad0,&(o_tmp[0])+(k*nquad1),1);
}
}
else if(orient == StdRegions::eDir1FwdDir2_Dir2FwdDir1)
......@@ -819,7 +814,7 @@ namespace Nektar
for (int j=0; j<nquad1; j++)
{
Vmath::Vcopy(nquad2,&(inarray[0])+j*nquad0,nquad0*nquad1,
&(outarray[0])+(j*nquad2),1);
&(o_tmp[0])+(j*nquad2),1);
}
}
else if(orient == StdRegions::eDir1FwdDir2_Dir2BwdDir1)
......@@ -828,7 +823,7 @@ namespace Nektar
for (int j=0; j<nquad1; j++)
{
Vmath::Vcopy(nquad2,&(inarray[0])+(nquad0*(nquad1-1)-j*nquad0),
nquad0*nquad1,&(outarray[0])+(j*nquad2),1);
nquad0*nquad1,&(o_tmp[0])+(j*nquad2),1);
}
}
else if(orient == StdRegions::eDir1BwdDir2_Dir2FwdDir1)
......@@ -837,7 +832,7 @@ namespace Nektar
for (int j=0; j<nquad1; j++)
{
Vmath::Vcopy(nquad2,&(inarray[0])+nquad0*nquad1*(nquad2-1)+j*nquad0,
-nquad0*nquad1,&(outarray[0])+(j*nquad2),1);
-nquad0*nquad1,&(o_tmp[0])+(j*nquad2),1);
}
}
else if(orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1)
......@@ -846,10 +841,9 @@ namespace Nektar
for (int j=0; j<nquad1; j++)
{
Vmath::Vcopy(nquad2,&(inarray[0])+(nquad0*(nquad1*nquad2-1)-j*nquad0),
-nquad0*nquad1,&(outarray[0])+(j*nquad2),1);
-nquad0*nquad1,&(o_tmp[0])+(j*nquad2),1);
}
}
o_tmp=outarray;
//interpolate
LibUtilities::Interp2D(m_base[1]->GetPointsKey(), m_base[2]->GetPointsKey(), o_tmp,
FaceExp->GetBasis(0)->GetPointsKey(),FaceExp->GetBasis(1)->GetPointsKey(),outarray);
......@@ -869,8 +863,19 @@ namespace Nektar
const Array<TwoD, const NekDouble> &df = geomFactors->GetDerivFactors(ptsKeys);
const Array<OneD, const NekDouble> &jac = geomFactors->GetJac(ptsKeys);
int nq0 = ptsKeys[0].GetNumPoints();
int nq1 = ptsKeys[1].GetNumPoints();
int nq2 = ptsKeys[2].GetNumPoints();
int nq01 = nq0*nq1;
int nqtot;
LibUtilities::BasisKey tobasis0 = DetFaceBasisKey(face,0);
LibUtilities::BasisKey tobasis1 = DetFaceBasisKey(face,1);
// Number of quadrature points in face expansion.
int nq = m_base[0]->GetNumPoints()*m_base[0]->GetNumPoints();
int nq_face = tobasis0.GetNumPoints()*tobasis1.GetNumPoints();
int vCoordDim = GetCoordim();
int i;
......@@ -878,7 +883,7 @@ namespace Nektar
Array<OneD, Array<OneD, NekDouble> > &normal = m_faceNormals[face];
for (i = 0; i < vCoordDim; ++i)
{
normal[i] = Array<OneD, NekDouble>(nq);
normal[i] = Array<OneD, NekDouble>(nq_face);
}
// Regular geometry case
......@@ -893,7 +898,7 @@ namespace Nektar
{
for(i = 0; i < vCoordDim; ++i)
{
Vmath::Fill(nq,-df[3*i+2][0],normal[i],1);
normal[i][0] = -df[3*i+2][0];;
}
break;
}
......@@ -901,7 +906,7 @@ namespace Nektar
{
for(i = 0; i < vCoordDim; ++i)
{
Vmath::Fill(nq,-df[3*i+1][0],normal[i],1);
normal[i][0] = -df[3*i+1][0];
}
break;
}
......@@ -909,7 +914,7 @@ namespace Nektar
{
for(i = 0; i < vCoordDim; ++i)
{
Vmath::Fill(nq,df[3*i][0]+df[3*i+2][0],normal[i],1);
normal[i][0] = df[3*i][0]+df[3*i+2][0];
}
break;
}
......@@ -917,7 +922,7 @@ namespace Nektar
{
for(i = 0; i < vCoordDim; ++i)
{
Vmath::Fill(nq,df[3*i+1][0],normal[i],1);
normal[i][0] = df[3*i+1][0];
}
break;
}
......@@ -925,7 +930,7 @@ namespace Nektar
{
for(i = 0; i < vCoordDim; ++i)
{
Vmath::Fill(nq,-df[3*i][0],normal[i],1);
normal[i][0] = -df[3*i][0];
}
break;
}
......@@ -942,21 +947,15 @@ namespace Nektar
fac = 1.0/sqrt(fac);
for (i = 0; i < vCoordDim; ++i)
{
Vmath::Smul(nq,fac,normal[i],1,normal[i],1);
Vmath::Fill(nq_face,fac*normal[i][0],normal[i],1);
}
}
else
{
// Set up deformed normals.
int j, k;
int nq0 = ptsKeys[0].GetNumPoints();
int nq1 = ptsKeys[1].GetNumPoints();
int nq2 = ptsKeys[2].GetNumPoints();
int nq01 = nq0*nq1;
int nqtot;
// Determine number of quadrature points on the face.
// Determine number of quadrature points on the face of 3D elmt
if (face == 0)
{
nqtot = nq0*nq1;
......@@ -973,7 +972,6 @@ namespace Nektar
LibUtilities::PointsKey points0;
LibUtilities::PointsKey points1;
Array<OneD, NekDouble> work (nq, 0.0);
Array<OneD, NekDouble> normals(vCoordDim*nqtot,0.0);
// Extract Jacobian along face and recover local derivatives
......@@ -1083,37 +1081,39 @@ namespace Nektar
ASSERTL0(false,"face is out of range (face < 4)");
}
Array<OneD, NekDouble> work (nq_face, 0.0);
// Interpolate Jacobian and invert
LibUtilities::Interp2D(points0, points1, jac,
m_base[0]->GetPointsKey(),
m_base[0]->GetPointsKey(),
tobasis0.GetPointsKey(),
tobasis1.GetPointsKey(),
work);
Vmath::Sdiv(nq, 1.0, &work[0], 1, &work[0], 1);
Vmath::Sdiv(nq_face, 1.0, &work[0], 1, &work[0], 1);
// Interpolate normal and multiply by inverse Jacobian.
for(i = 0; i < vCoordDim; ++i)
{
LibUtilities::Interp2D(points0, points1,
&normals[i*nqtot],
m_base[0]->GetPointsKey(),
m_base[0]->GetPointsKey(),
tobasis0.GetPointsKey(),
tobasis1.GetPointsKey(),
&normal[i][0]);
Vmath::Vmul(nq,work,1,normal[i],1,normal[i],1);
Vmath::Vmul(nq_face,work,1,normal[i],1,normal[i],1);
}
// Normalise to obtain unit normals.
Vmath::Zero(nq,work,1);
Vmath::Zero(nq_face,work,1);
for(i = 0; i < GetCoordim(); ++i)
{
Vmath::Vvtvp(nq,normal[i],1,normal[i],1,work,1,work,1);
Vmath::Vvtvp(nq_face,normal[i],1,normal[i],1,work,1,work,1);
}
Vmath::Vsqrt(nq,work,1,work,1);
Vmath::Sdiv (nq,1.0,work,1,work,1);
Vmath::Vsqrt(nq_face,work,1,work,1);
Vmath::Sdiv (nq_face,1.0,work,1,work,1);
for(i = 0; i < GetCoordim(); ++i)
{
Vmath::Vmul(nq,normal[i],1,work,1,normal[i],1);
Vmath::Vmul(nq_face,normal[i],1,work,1,normal[i],1);
}
}
}
......
......@@ -365,7 +365,7 @@ namespace Nektar
int nq1 = m_base[1]->GetNumPoints();
int nq2 = m_base[2]->GetNumPoints();
Array<OneD,NekDouble> o_tmp(nq0*nq1*nq2);
Array<OneD,NekDouble> o_tmp(GetFaceNumPoints(face));
if (orient == StdRegions::eNoOrientation)
{
......@@ -510,8 +510,12 @@ namespace Nektar
const Array<TwoD, const NekDouble> &df = geomFactors->GetDerivFactors(ptsKeys);
const Array<OneD, const NekDouble> &jac = geomFactors->GetJac(ptsKeys);
LibUtilities::BasisKey tobasis0 = DetFaceBasisKey(face,0);
LibUtilities::BasisKey tobasis1 = DetFaceBasisKey(face,1);
// Number of quadrature points in face expansion.
int nq = m_base[0]->GetNumPoints()*m_base[0]->GetNumPoints();
int nq_face = tobasis0.GetNumPoints()*tobasis1.GetNumPoints();
int vCoordDim = GetCoordim();
int i;
......@@ -519,7 +523,7 @@ namespace Nektar
Array<OneD, Array<OneD, NekDouble> > &normal = m_faceNormals[face];
for (i = 0; i < vCoordDim; ++i)
{
normal[i] = Array<OneD, NekDouble>(nq);
normal[i] = Array<OneD, NekDouble>(nq_face);
}
// Regular geometry case
......@@ -534,7 +538,7 @@ namespace Nektar
{
for(i = 0; i < vCoordDim; ++i)
{
Vmath::Fill(nq,-df[3*i+2][0],normal[i],1);
normal[i][0] = -df[3*i+2][0];
}
break;
}
......@@ -542,7 +546,7 @@ namespace Nektar
{
for(i = 0; i < vCoordDim; ++i)
{
Vmath::Fill(nq,-df[3*i+1][0],normal[i],1);
normal[i][0] = -df[3*i+1][0];
}
break;
}
......@@ -550,7 +554,7 @@ namespace Nektar
{
for(i = 0; i < vCoordDim; ++i)
{
Vmath::Fill(nq,df[3*i][0]+df[3*i+2][0],normal[i],1);
normal[i][0] = df[3*i][0]+df[3*i+2][0];
}
break;
}
......@@ -558,7 +562,7 @@ namespace Nektar
{
for(i = 0; i < vCoordDim; ++i)
{
Vmath::Fill(nq,df[3*i+1][0]+df[3*i+2][0],normal[i],1);
normal[i][0] = df[3*i+1][0]+df[3*i+2][0];
}
break;
}
......@@ -566,7 +570,7 @@ namespace Nektar
{
for(i = 0; i < vCoordDim; ++i)
{
Vmath::Fill(nq,-df[3*i][0],normal[i],1);
normal[i][0] = -df[3*i][0];
}
break;
}
......@@ -583,7 +587,7 @@ namespace Nektar
fac = 1.0/sqrt(fac);
for (i = 0; i < vCoordDim; ++i)
{
Vmath::Smul(nq,fac,normal[i],1,normal[i],1);
Vmath::Fill(nq_face,fac*normal[i][0],normal[i],1);
}
}
else
......@@ -614,7 +618,6 @@ namespace Nektar
LibUtilities::PointsKey points0;
LibUtilities::PointsKey points1;
Array<OneD, NekDouble> work (nq, 0.0);
Array<OneD, NekDouble> normals(vCoordDim*nqtot,0.0);
// Extract Jacobian along face and recover local derivatives
......@@ -724,37 +727,38 @@ namespace Nektar
ASSERTL0(false,"face is out of range (face < 4)");
}
Array<OneD, NekDouble> work (nq_face, 0.0);
// Interpolate Jacobian and invert
LibUtilities::Interp2D(points0, points1, jac,
m_base[0]->GetPointsKey(),
m_base[0]->GetPointsKey(),
tobasis0.GetPointsKey(),
tobasis1.GetPointsKey(),
work);
Vmath::Sdiv(nq, 1.0, &work[0], 1, &work[0], 1);
Vmath::Sdiv(nq_face, 1.0, &work[0], 1, &work[0], 1);
// Interpolate normal and multiply by inverse Jacobian.
for(i = 0; i < vCoordDim; ++i)
{
LibUtilities::Interp2D(points0, points1,
&normals[i*nqtot],
m_base[0]->GetPointsKey(),
m_base[0]->GetPointsKey(),
tobasis0.GetPointsKey(),
tobasis1.GetPointsKey(),
&normal[i][0]);
Vmath::Vmul(nq,work,1,normal[i],1,normal[i],1);
Vmath::Vmul(nq_face,work,1,normal[i],1,normal[i],1);
}
// Normalise to obtain unit normals.
Vmath::Zero(nq,work,1);
Vmath::Zero(nq_face,work,1);
for(i = 0; i < GetCoordim(); ++i)
{
Vmath::Vvtvp(nq,normal[i],1,normal[i],1,work,1,work,1);
Vmath::Vvtvp(nq_face,normal[i],1,normal[i],1,work,1,work,1);
}
Vmath::Vsqrt(nq,work,1,work,1);
Vmath::Sdiv (nq,1.0,work,1,work,1);
Vmath::Vsqrt(nq_face,work,1,work,1);
Vmath::Sdiv (nq_face,1.0,work,1,work,1);
for(i = 0; i < GetCoordim(); ++i)
{
Vmath::Vmul(nq,normal[i],1,work,1,normal[i],1);
Vmath::Vmul(nq_face,normal[i],1,work,1,normal[i],1);
}
}
}
......
......@@ -641,8 +641,12 @@ namespace Nektar
//Directions A and B positive
Vmath::Vcopy(nquad0*nquad1,inarray.get(),1,o_tmp.get(),1);
//interpolate
LibUtilities::Interp2D(m_base[0]->GetPointsKey(), m_base[1]->GetPointsKey(), o_tmp.get(),
FaceExp->GetBasis(0)->GetPointsKey(),FaceExp->GetBasis(1)->GetPointsKey(),o_tmp2.get());
LibUtilities::Interp2D(m_base[0]->GetPointsKey(),
m_base[1]->GetPointsKey(),
o_tmp.get(),
FaceExp->GetBasis(0)->GetPointsKey(),
FaceExp->GetBasis(1)->GetPointsKey(),
o_tmp2.get());