Commit f88f7967 authored by Gianmarco Mengaldo's avatar Gianmarco Mengaldo

Merge branch 'fix/ExtractSurface3DH1D' into 'master'

Fix/extractsurface3DH1D

This MR extends the extract surface module of FieldConvert to work with 3DH1D expansions. The addnormals option is also extended to 2D and 3DH1D expansions.

I had to force FieldConvert to declare all expansions using NewField=true, since otherwise it would not work. This seems to be a problem with the copy constructor from ContField3DH1D (maybe with the boundary expansions?).

See merge request !521
parents 56652230 cdb3b162
......@@ -157,6 +157,14 @@ namespace Nektar
}
}
void ContField3DHomogeneous1D::v_FillBndCondFromField(void)
{
for(int n = 0; n < m_planes.num_elements(); ++n)
{
m_planes[n]->FillBndCondFromField();
}
}
/**
*
*/
......
......@@ -73,6 +73,8 @@ namespace Nektar
private:
virtual void v_ImposeDirichletConditions(Array<OneD,NekDouble>& outarray);
virtual void v_FillBndCondFromField();
/// Template method virtual forwarded for LocalToGlobal()
virtual void v_LocalToGlobal(void);
......
......@@ -68,6 +68,7 @@ struct Field {
Field() : m_verbose(false),
m_declareExpansionAsContField(false),
m_declareExpansionAsDisContField(false),
m_declareAsNewField(false),
m_writeBndFld(false),
m_fldToBnd(false),
m_addNormals(false),
......@@ -88,6 +89,7 @@ struct Field {
bool m_declareExpansionAsContField;
bool m_declareExpansionAsDisContField;
bool m_declareAsNewField;
bool m_useFFT;
......@@ -391,6 +393,11 @@ struct Field {
string var = "DefaultVar",
bool NewField = false)
{
if (m_declareAsNewField)
{
NewField = true;
var = m_session->GetVariables()[0];
}
MultiRegions::ExpListSharedPtr tmp;
switch (m_graph->GetMeshDimension())
{
......@@ -398,8 +405,8 @@ struct Field {
{
if (NumHomogeneousDir == 1)
{
ASSERTL0(m_declareExpansionAsContField ||
m_declareExpansionAsDisContField,
ASSERTL0( !(m_declareExpansionAsContField ||
m_declareExpansionAsDisContField),
"ContField2DHomogeneous1D or "
"DisContField2DHomogenenous1D has not been "
"implemented");
......
......@@ -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)
{
......
......@@ -62,6 +62,7 @@ ProcessBoundaryExtract::ProcessBoundaryExtract(FieldSharedPtr f) : ProcessModule
f->m_writeBndFld = true;
f->m_declareExpansionAsContField = true;
f->m_declareAsNewField = true;
// check for correct input files
if((f->m_inputfiles.count("xml") == 0)&&(f->m_inputfiles.count("xml.gz") == 0))
......
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