Commit c7443b35 authored by Douglas Serson's avatar Douglas Serson

Extend addnormals to 2D and 3DH1D

parent 4d6b9c53
......@@ -161,14 +161,11 @@ void OutputFld::Process(po::variables_map &vm)
if(m_f->m_addNormals)
{
ASSERTL0(m_f->m_exp[0]->GetCoordim(0) == 3,
"Add normals to extracted boundaries only set up in 3 dimensions");
int normdim = 3; // currently assuming 3D normals;
int normdim = m_f->m_graph->GetMeshDimension();
string normstr[3] = {"Norm_x","Norm_y","Norm_z"};
// Add normal information
StdRegions::StdExpansionSharedPtr elmt;
StdRegions::StdExpansion2DSharedPtr bc;
Array<OneD, int> BoundarytoElmtID, BoundarytoTraceID;
m_f->m_exp[0]->GetBoundaryToElmtMap(BoundarytoElmtID,
......@@ -178,46 +175,105 @@ void OutputFld::Process(po::variables_map &vm)
int cnt = 0;
for(int n = 0; n < Border; ++n)
{
cnt += BndExp[0][n]->GetExpSize();
if (m_f->m_fielddef[0]->m_numHomogeneousDir == 1)
{
cnt += BndExp[0][n]->GetPlane(0)->GetExpSize();
}
else
{
cnt += BndExp[0][n]->GetExpSize();
}
}
Array<OneD, NekDouble> tmp_array;
Array<OneD, Array<OneD, NekDouble> > NormCoeff(normdim);
Array<OneD, Array<OneD, NekDouble> > NormPhys(normdim);
for(int j = 0; j < normdim; ++j)
{
NormCoeff[j] = Array<OneD, NekDouble>(BndExp[0][Border]->GetNcoeffs(),0.0);
NormPhys[j] = Array<OneD, NekDouble>(BndExp[0][Border]->GetTotPoints(),0.0);
}
// setup coeff arrays of normals.
// setup phys arrays of normals
for(int j = 0; j < BndExp[0][Border]->GetExpSize();
++j, cnt++)
++j)
{
// find element and face of this expansion.
int elmtid = BoundarytoElmtID[cnt];
elmt = m_f->m_exp[0]->GetExp(elmtid);
int cnt2;
if (m_f->m_fielddef[0]->m_numHomogeneousDir == 1)
{
int exp_per_plane = BndExp[0][Border]->GetPlane(0)->
GetExpSize();
cnt2 = cnt + (j%exp_per_plane);
}
else
{
cnt2 = cnt+j;
}
int elmtid = BoundarytoElmtID[cnt2];
// Get face 2D expansion from element expansion
bc = boost::dynamic_pointer_cast<StdRegions::StdExpansion2D> (BndExp[0][Border]->GetExp(j));
elmt = m_f->m_exp[0]->GetExp(elmtid);
//identify boundary of element looking at.
int boundary = BoundarytoTraceID[cnt];
//Get face normals
const Array<OneD, const Array<OneD, NekDouble> > normals = elmt->GetFaceNormal(boundary);
int boundary = BoundarytoTraceID[cnt2];
for(int k = 0; k < normdim; ++k)
// Dimension specific part
switch(normdim)
{
bc->FwdTrans(normals[k],tmp_array = NormCoeff[k]+BndExp[0][Border]->GetCoeff_Offset(j));
case 2:
{
// Get edge 1D expansion from element expansion
StdRegions::StdExpansion1DSharedPtr bc;
bc = boost::dynamic_pointer_cast
<StdRegions::StdExpansion1D>
(BndExp[0][Border]->GetExp(j));
// Get edge normals
const Array<OneD, const Array<OneD, NekDouble> >
normals = elmt->GetEdgeNormal(boundary);
for(int k = 0; k < normdim; ++k)
{
Vmath::Vcopy(bc->GetTotPoints(),
normals[k], 1,
tmp_array = NormPhys[k]+
BndExp[0][Border]->
GetPhys_Offset(j), 1);
}
}
break;
case 3:
{
// Get face 2D expansion from element expansion
StdRegions::StdExpansion2DSharedPtr bc;
bc = boost::dynamic_pointer_cast
<StdRegions::StdExpansion2D>
(BndExp[0][Border]->GetExp(j));
//Get face normals
const Array<OneD, const Array<OneD, NekDouble> >
normals = elmt->GetFaceNormal(boundary);
for(int k = 0; k < normdim; ++k)
{
Vmath::Vcopy(bc->GetTotPoints(),
normals[k], 1,
tmp_array = NormPhys[k]+
BndExp[0][Border]->
GetPhys_Offset(j), 1);
}
}
break;
default:
ASSERTL0(false, "Addnormals requires expdim >=2.");
break;
}
}
// add normal coefficients to list to be dumped
for (int j = 0; j < normdim; ++j)
{
Vmath::Vcopy(BndExp[0][Border]->GetNcoeffs(),
NormCoeff[j],1,
BndExp[0][Border]->UpdateCoeffs(),1);
BndExp[0][Border]->FwdTrans( NormPhys[j],
BndExp[0][Border]->UpdateCoeffs());
for (int k = 0; k < FieldDef.size(); ++k)
{
......
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