Commit 99fa0425 authored by Gianmarco Mengaldo's avatar Gianmarco Mengaldo

Merge branch 'master' into feature/Isocontour3DH1D

parents 7caa333e 1c576935
......@@ -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);
......
......@@ -168,8 +168,8 @@ namespace Nektar
m_ncoeffs = ncoeffs_per_plane*nzplanes;
m_npoints = npoints_per_plane*nzplanes;
m_coeffs = Array<OneD, NekDouble> (m_ncoeffs);
m_phys = Array<OneD, NekDouble> (m_npoints);
m_coeffs = Array<OneD, NekDouble> (m_ncoeffs, 0.0);
m_phys = Array<OneD, NekDouble> (m_npoints, 0.0);
int nel = m_planes[0]->GetExpSize();
m_coeff_offset = Array<OneD,int>(nel*nzplanes);
......
......@@ -588,8 +588,24 @@ namespace Nektar
std::vector<NekDouble> HomoLen(2);
HomoLen[0] = m_lhom_y;
HomoLen[1] = m_lhom_z;
int nhom_modes_y = m_homogeneousBasis_y->GetNumModes();
int nhom_modes_z = m_homogeneousBasis_z->GetNumModes();
std::vector<unsigned int> yIDs;
std::vector<unsigned int> zIDs;
for(int n = 0; n < nhom_modes_z; ++n)
{
for(int m = 0; m < nhom_modes_y; ++m)
{
zIDs.push_back(n);
yIDs.push_back(m);
}
}
m_lines[0]->GeneralGetFieldDefinitions(returnval, 2, 1, HomoBasis, HomoLen);
m_lines[0]->GeneralGetFieldDefinitions(returnval, 2, 1, HomoBasis,
HomoLen, zIDs, yIDs);
return returnval;
}
......@@ -603,8 +619,24 @@ namespace Nektar
HomoLen[0] = m_lhom_y;
HomoLen[1] = m_lhom_z;
int nhom_modes_y = m_homogeneousBasis_y->GetNumModes();
int nhom_modes_z = m_homogeneousBasis_z->GetNumModes();
std::vector<unsigned int> yIDs;
std::vector<unsigned int> zIDs;
for(int n = 0; n < nhom_modes_z; ++n)
{
for(int m = 0; m < nhom_modes_y; ++m)
{
zIDs.push_back(n);
yIDs.push_back(m);
}
}
// enforce NumHomoDir == 1 by direct call
m_lines[0]->GeneralGetFieldDefinitions(fielddef, 2, 1, HomoBasis, HomoLen);
m_lines[0]->GeneralGetFieldDefinitions(fielddef, 2, 1, HomoBasis,
HomoLen, zIDs, yIDs);
}
void ExpListHomogeneous2D::v_AppendFieldData(LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector<NekDouble> &fielddata, Array<OneD, NekDouble> &coeffs)
......@@ -683,6 +715,29 @@ namespace Nektar
}
}
}
void ExpListHomogeneous2D::v_WriteVtkPieceData(std::ostream &outfile, int expansion,
std::string var)
{
int i;
int nq = (*m_exp)[expansion]->GetTotPoints();
int npoints_per_line = m_lines[0]->GetTotPoints();
// printing the fields of that zone
outfile << " <DataArray type=\"Float32\" Name=\""
<< var << "\">" << endl;
outfile << " ";
for (int n = 0; n < m_lines.num_elements(); ++n)
{
const Array<OneD, NekDouble> phys = m_phys + m_phys_offset[expansion] + n*npoints_per_line;
for(i = 0; i < nq; ++i)
{
outfile << (fabs(phys[i]) < NekConstants::kNekZeroTol ? 0 : phys[i]) << " ";
}
}
outfile << endl;
outfile << " </DataArray>" << endl;
}
void ExpListHomogeneous2D::v_PhysDeriv(const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &out_d0,
......
......@@ -197,6 +197,9 @@ namespace Nektar
virtual void v_ExtractDataToCoeffs(LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector<NekDouble> &fielddata, std::string &field, Array<OneD, NekDouble> &coeffs);
virtual void v_WriteVtkPieceData(std::ostream &outfile, int expansion,
std::string var);
virtual void v_HomogeneousFwdTrans(const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray,
CoeffState coeffstate = eLocal,
......
......@@ -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