Commit 50349b86 authored by Gianmarco Mengaldo's avatar Gianmarco Mengaldo

Added support for BCs from file. At least 1 bug (i.e. var and varName always...

Added support for BCs from file. At least 1 bug (i.e. var and varName always as default values, don't know why).
parent 36763032
......@@ -102,9 +102,9 @@ namespace Nektar
m_trace = boost::dynamic_pointer_cast<ExpList>(trace);
m_traceMap = MemoryManager<AssemblyMapDG>::
AllocateSharedPtr(pSession,graph1D,trace,*this,
m_bndCondExpansions,m_bndConditions,periodicVertices,
variable);
AllocateSharedPtr(pSession, graph1D, trace, *this,
m_bndCondExpansions, m_bndConditions,
periodicVertices, variable);
tmpBndSol = Array<OneD,NekDouble>
(m_traceMap->GetNumLocalBndCoeffs());
......@@ -1230,9 +1230,12 @@ namespace Nektar
* @param bndCondExpansions List of boundary expansions.
* @param bndConditions Information about the boundary conditions.
*/
void DisContField1D::v_EvaluateBoundaryConditions(const NekDouble time,
const NekDouble x2_in,
const NekDouble x3_in)
void DisContField1D::v_EvaluateBoundaryConditions(
const NekDouble time,
int var,
std::string varName,
const NekDouble x2_in,
const NekDouble x3_in)
{
int i;
......@@ -1295,7 +1298,7 @@ namespace Nektar
}
else
{
ASSERTL0(false,"This type of BC not implemented yet");
ASSERTL0(false, "This type of BC not implemented yet");
}
}
}
......
......@@ -224,6 +224,8 @@ namespace Nektar
/// Evaluate all boundary conditions at a given time..
virtual void v_EvaluateBoundaryConditions(
const NekDouble time = 0.0,
int var = 0,
std::string varName = "",
const NekDouble x2_in = NekConstants::kNekUnsetDouble,
const NekDouble x3_in = NekConstants::kNekUnsetDouble);
......
......@@ -2116,9 +2116,12 @@ namespace Nektar
* boundary conditions unless time == 0.0 which is the
* case when the method is called from the constructor.
*/
void DisContField2D::v_EvaluateBoundaryConditions(const NekDouble time,
const NekDouble x2_in,
const NekDouble x3_in)
void DisContField2D::v_EvaluateBoundaryConditions(
const NekDouble time,
int var,
std::string varName,
const NekDouble x2_in,
const NekDouble x3_in)
{
int i;
int npoints;
......@@ -2141,7 +2144,7 @@ namespace Nektar
// Homogeneous input case for x2.
if (x2_in == NekConstants::kNekUnsetDouble)
{
locExpList->GetCoords(x0,x1,x2);
locExpList->GetCoords(x0, x1, x2);
}
else
{
......@@ -2159,32 +2162,50 @@ namespace Nektar
if (filebcs != "")
{
string varString = filebcs.substr(
0, filebcs.find_last_of("."));
int len = varString.length();
varString = varString.substr(len-1, len);
cout << "Boundary condition from file:"
<< filebcs << endl;
std::vector<LibUtilities::
FieldDefinitionsSharedPtr> FieldDef;
std::vector<std::vector<NekDouble> > FieldData;
LibUtilities::FieldIO f(m_session->GetComm());
f.Import(filebcs, FieldDef, FieldData);
// copy FieldData into locExpList
locExpList->ExtractDataToCoeffs(
FieldDef[0], FieldData[0],
FieldDef[0]->m_fields[0],
locExpList->UpdateCoeffs());
locExpList->BwdTrans_IterPerExp(
locExpList->GetCoeffs(),
locExpList->UpdatePhys());
string varString = filebcs.substr(
0, filebcs.find_last_of("."));
int len = varString.length();
varString = varString.substr(len-1, len);
cout << "Boundary condition from file:"
<< filebcs << endl;
cout << "var = " << var << endl;
cout << "varName = " << varName << endl;
std::vector<LibUtilities::
FieldDefinitionsSharedPtr> FieldDef;
std::vector<std::vector<NekDouble> > FieldData;
LibUtilities::FieldIO f(m_session->GetComm());
f.Import(filebcs, FieldDef, FieldData);
locExpList->FwdTrans_BndConstrained(
locExpList->GetPhys(),
locExpList->UpdateCoeffs());
const std::string varNameFld = FieldDef[0]->
m_fields[var];
cout << "varNameFld = " << varNameFld << endl;
if (varNameFld == varName)
{
// Copy FieldData into locExpList
locExpList->ExtractDataToCoeffs(
FieldDef[0], FieldData[0],
FieldDef[0]->m_fields[var],
locExpList->UpdateCoeffs());
locExpList->BwdTrans_IterPerExp(
locExpList->GetCoeffs(),
locExpList->UpdatePhys());
locExpList->FwdTrans_BndConstrained(
locExpList->GetPhys(),
locExpList->UpdateCoeffs());
}
else
{
ASSERTL0(
false,
"BCs fields order in session file do not "
"match variable order in the fld file");
}
}
else
{
......@@ -2211,44 +2232,50 @@ namespace Nektar
if (filebcs != "")
{
string var = filebcs.substr(
0, filebcs.find_last_of("."));
int len=var.length();
var = var.substr(len-1,len);
cout << "Boundary condition from file: "
<< filebcs << endl;
std::vector<LibUtilities::
FieldDefinitionsSharedPtr> FieldDef;
std::vector<std::vector<NekDouble> > FieldData;
LibUtilities::FieldIO f(m_session->GetComm());
f.Import(filebcs, FieldDef, FieldData);
// copy FieldData into locExpList
locExpList->ExtractDataToCoeffs(
FieldDef[0], FieldData[0],
FieldDef[0]->m_fields[0],
locExpList->UpdateCoeffs());
locExpList->BwdTrans_IterPerExp(
locExpList->GetCoeffs(),
locExpList->UpdatePhys());
/*
Array<OneD, NekDouble> x(locExpList->GetTotPoints(),0.0);
Array<OneD, NekDouble> y(locExpList->GetTotPoints(),0.0);
locExpList->GetCoords(x,y);
for(int i=0; i< locExpList->GetTotPoints(); i++)
{
cout<<i<<" "<<x[i]<<" "<<y[i]<<" "
<<locExpList->GetPhys()[i]<<endl;
}
*/
locExpList->IProductWRTBase(
locExpList->GetPhys(),
locExpList->UpdateCoeffs());
string varString = filebcs.substr(
0, filebcs.find_last_of("."));
int len = varString.length();
varString = varString.substr(len-1, len);
cout << "Boundary condition from file:"
<< filebcs << endl;
cout << "var = " << var << endl;
cout << "varName = " << varName << endl;
std::vector<LibUtilities::
FieldDefinitionsSharedPtr> FieldDef;
std::vector<std::vector<NekDouble> > FieldData;
LibUtilities::FieldIO f(m_session->GetComm());
f.Import(filebcs, FieldDef, FieldData);
const std::string varNameFld = FieldDef[0]->
m_fields[var];
cout << "varNameFld = " << varNameFld << endl;
if (varNameFld == varName)
{
// Copy FieldData into locExpList
locExpList->ExtractDataToCoeffs(
FieldDef[0], FieldData[0],
FieldDef[0]->m_fields[var],
locExpList->UpdateCoeffs());
locExpList->BwdTrans_IterPerExp(
locExpList->GetCoeffs(),
locExpList->UpdatePhys());
locExpList->FwdTrans_BndConstrained(
locExpList->GetPhys(),
locExpList->UpdateCoeffs());
}
else
{
ASSERTL0(
false,
"BCs fields order in session file do not "
"match variable order in the fld file");
}
}
else
{
......@@ -2274,11 +2301,11 @@ namespace Nektar
if (filebcs != "")
{
//Never tested!!!
string var = filebcs.substr(
// Never tested
string varString = filebcs.substr(
0, filebcs.find_last_of("."));
int len = var.length();
var = var.substr(len-1,len);
int len = varString.length();
varString = varString.substr(len-1, len);
std::vector<LibUtilities::
FieldDefinitionsSharedPtr> FieldDef;
......@@ -2286,28 +2313,43 @@ namespace Nektar
LibUtilities::FieldIO f(m_session->GetComm());
f.Import(filebcs, FieldDef, FieldData);
// copy FieldData into locExpList
locExpList->ExtractDataToCoeffs(
FieldDef[0], FieldData[0],
FieldDef[0]->m_fields[0],
locExpList->UpdateCoeffs());
locExpList->BwdTrans_IterPerExp(
locExpList->GetCoeffs(),
locExpList->UpdatePhys());
locExpList->IProductWRTBase(
locExpList->GetPhys(),
locExpList->UpdateCoeffs());
LibUtilities::Equation coeff =
boost::static_pointer_cast<
SpatialDomains::RobinBoundaryCondition
>(m_bndConditions[i])->m_robinPrimitiveCoeff;
// Array<OneD,NekDouble> timeArray(npoints, time);
// put primitive coefficient into the physical space
// storage
coeff.Evaluate(x0,x1,x2,time,
locExpList->UpdatePhys());
const std::string varNameFld = FieldDef[0]->
m_fields[var];
if (varNameFld == varName)
{
// Copy FieldData into locExpList
locExpList->ExtractDataToCoeffs(
FieldDef[0], FieldData[0],
FieldDef[0]->m_fields[var],
locExpList->UpdateCoeffs());
locExpList->BwdTrans_IterPerExp(
locExpList->GetCoeffs(),
locExpList->UpdatePhys());
locExpList->IProductWRTBase(
locExpList->GetPhys(),
locExpList->UpdateCoeffs());
LibUtilities::Equation coeff =
boost::static_pointer_cast<SpatialDomains::
RobinBoundaryCondition>(
m_bndConditions[i])->
m_robinPrimitiveCoeff;
// Array<OneD,NekDouble> timeArray(npoints, time);
// Put primitive coefficient into the physical
// space storage
coeff.Evaluate(x0, x1, x2, time,
locExpList->UpdatePhys());
}
else
{
ASSERTL0(
false,
"BCs fields order in session file do not "
"match variable order in the fld file");
}
}
else
{
......@@ -2327,8 +2369,8 @@ namespace Nektar
locExpList->GetPhys(),
locExpList->UpdateCoeffs());
// put primitive coefficient into the physical space
// storage
// put primitive coefficient into the physical
// space storage
coeff.Evaluate(x0, x1, x2, time,
locExpList->UpdatePhys());
}
......
......@@ -250,6 +250,8 @@ namespace Nektar
virtual void v_EvaluateBoundaryConditions(
const NekDouble time = 0.0,
int var = 0,
std::string varName = "",
const NekDouble x2_in = NekConstants::kNekUnsetDouble,
const NekDouble x3_in = NekConstants::kNekUnsetDouble);
......
......@@ -2392,110 +2392,130 @@
* @param bndCondExpansions List of boundary conditions.
* @param bndConditions Information about the boundary conditions.
*/
void DisContField3D::v_EvaluateBoundaryConditions(const NekDouble time,
const NekDouble x2_in,
const NekDouble x3_in)
void DisContField3D::v_EvaluateBoundaryConditions(
const NekDouble time,
int var,
std::string varName,
const NekDouble x2_in,
const NekDouble x3_in)
{
int i;
int npoints;
int nbnd = m_bndCondExpansions.num_elements();
MultiRegions::ExpListSharedPtr locExpList;
for(i = 0; i < nbnd; ++i)
for (i = 0; i < nbnd; ++i)
{
if(time == 0.0 || m_bndConditions[i]->GetUserDefined() ==
SpatialDomains::eTimeDependent)
if (time == 0.0 || m_bndConditions[i]->GetUserDefined() ==
SpatialDomains::eTimeDependent)
{
locExpList = m_bndCondExpansions[i];
npoints = locExpList->GetNpoints();
npoints = locExpList->GetNpoints();
Array<OneD,NekDouble> x0(npoints,0.0);
Array<OneD,NekDouble> x1(npoints,0.0);
Array<OneD,NekDouble> x2(npoints,0.0);
Array<OneD, NekDouble> x0(npoints, 0.0);
Array<OneD, NekDouble> x1(npoints, 0.0);
Array<OneD, NekDouble> x2(npoints, 0.0);
locExpList->GetCoords(x0,x1,x2);
locExpList->GetCoords(x0, x1, x2);
if(m_bndConditions[i]->GetBoundaryConditionType()
== SpatialDomains::eDirichlet)
if (m_bndConditions[i]->GetBoundaryConditionType()
== SpatialDomains::eDirichlet)
{
string filebcs = boost::static_pointer_cast<
SpatialDomains::DirichletBoundaryCondition>(
m_bndConditions[i])->m_filename;
if(filebcs != "")
if (filebcs != "")
{
string var = filebcs.substr(
0, filebcs.find_last_of("."));
int len = var.length();
var = var.substr(len-1,len);
cout << "Boundary condition from file:"
<< filebcs << endl;
std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
std::vector<std::vector<NekDouble> > FieldData;
Import(filebcs,FieldDef, FieldData);
// copy FieldData into locExpList
locExpList->ExtractDataToCoeffs(
FieldDef[0], FieldData[0],
FieldDef[0]->m_fields[0], locExpList->UpdateCoeffs());
locExpList->BwdTrans_IterPerExp(
locExpList->GetCoeffs(),
locExpList->UpdatePhys());
locExpList->FwdTrans_BndConstrained(
locExpList->GetPhys(),
locExpList->UpdateCoeffs());
std::vector<LibUtilities::
FieldDefinitionsSharedPtr> FieldDef;
std::vector<std::vector<NekDouble> > FieldData;
Import(filebcs, FieldDef, FieldData);
const std::string varNameFld = FieldDef[0]->
m_fields[var];
if (varNameFld == varName)
{
// Copy FieldData into locExpList
locExpList->ExtractDataToCoeffs(
FieldDef[0], FieldData[0],
FieldDef[0]->m_fields[var],
locExpList->UpdateCoeffs());
locExpList->BwdTrans_IterPerExp(
locExpList->GetCoeffs(),
locExpList->UpdatePhys());
locExpList->FwdTrans_BndConstrained(
locExpList->GetPhys(),
locExpList->UpdateCoeffs());
}
else
{
ASSERTL0(
false,
"BCs fields order in session file do not "
"match variable order in the fld file");
}
}
else
{
LibUtilities::Equation condition = boost::static_pointer_cast<
SpatialDomains::DirichletBoundaryCondition >(m_bndConditions[i])->m_dirichletCondition;
LibUtilities::Equation condition = boost::static_pointer_cast<SpatialDomains::
DirichletBoundaryCondition >(
m_bndConditions[i])->m_dirichletCondition;
condition.Evaluate(x0,x1,x2,time,locExpList->UpdatePhys());
condition.Evaluate(x0, x1, x2, time,
locExpList->UpdatePhys());
locExpList->FwdTrans_BndConstrained(locExpList->GetPhys(),
locExpList->UpdateCoeffs());
locExpList->FwdTrans_BndConstrained(
locExpList->GetPhys(),
locExpList->UpdateCoeffs());
}
}
else if(m_bndConditions[i]->GetBoundaryConditionType()
== SpatialDomains::eNeumann)
else if (m_bndConditions[i]->GetBoundaryConditionType()
== SpatialDomains::eNeumann)
{
LibUtilities::Equation condition = boost::static_pointer_cast<
SpatialDomains::NeumannBoundaryCondition
>(m_bndConditions[i])->m_neumannCondition;
LibUtilities::Equation condition = boost::
static_pointer_cast<SpatialDomains::
NeumannBoundaryCondition>(
m_bndConditions[i])->m_neumannCondition;
condition.Evaluate(x0,x1,x2,time,locExpList->UpdatePhys());
condition.Evaluate(x0, x1, x2, time,
locExpList->UpdatePhys());
locExpList->IProductWRTBase(locExpList->GetPhys(),
locExpList->UpdateCoeffs());
}
else if(m_bndConditions[i]->GetBoundaryConditionType()
== SpatialDomains::eRobin)
else if (m_bndConditions[i]->GetBoundaryConditionType()
== SpatialDomains::eRobin)
{
LibUtilities::Equation condition = boost::static_pointer_cast<
SpatialDomains::RobinBoundaryCondition
>(m_bndConditions[i])->m_robinFunction;
LibUtilities::Equation condition = boost::
static_pointer_cast<SpatialDomains::
RobinBoundaryCondition>(
m_bndConditions[i])->m_robinFunction;
LibUtilities::Equation coeff =
boost::static_pointer_cast<
SpatialDomains::RobinBoundaryCondition
>(m_bndConditions[i])->m_robinPrimitiveCoeff;
LibUtilities::Equation coeff = boost::
static_pointer_cast<SpatialDomains::
RobinBoundaryCondition>(
m_bndConditions[i])->m_robinPrimitiveCoeff;
condition.Evaluate(x0,x1,x2,time,locExpList->UpdatePhys());
condition.Evaluate(x0, x1, x2, time,
locExpList->UpdatePhys());
locExpList->IProductWRTBase(locExpList->GetPhys(),
locExpList->UpdateCoeffs());
// put primitive coefficient into the physical space
// storage
coeff.Evaluate(x0,x1,x2,time,
// Put primitive coefficient into the physical
// space storage
coeff.Evaluate(x0, x1, x2, time,
locExpList->UpdatePhys());
}
else
{
ASSERTL0(false,"This type of BC not implemented yet");
ASSERTL0(false, "This type of BC not implemented yet");
}
}
}
......
......@@ -238,6 +238,8 @@ namespace Nektar
virtual void v_EvaluateBoundaryConditions(
const NekDouble time = 0.0,
int var = 0,
std::string varName = "",
const NekDouble x2_in = NekConstants::kNekUnsetDouble,
const NekDouble x3_in = NekConstants::kNekUnsetDouble);
......
......@@ -220,7 +220,9 @@ namespace Nektar
}
void DisContField3DHomogeneous1D::EvaluateBoundaryConditions(
const NekDouble time)
const NekDouble time,
int var,
std::string varName)
{
int n;
const Array<OneD, const NekDouble> z = m_homogeneousBasis->GetZ();
......@@ -234,7 +236,7 @@ namespace Nektar
for (n = 0; n < m_planes.num_elements(); ++n)
{
m_planes[n]->EvaluateBoundaryConditions(
time,0.5*m_lhom*(1.0+local_z[n]));
time, var, varName, 0.5*m_lhom*(1.0+local_z[n]));
}
// Fourier transform coefficient space boundary values
......@@ -302,10 +304,12 @@ namespace Nektar
void DisContField3DHomogeneous1D::v_EvaluateBoundaryConditions(
const NekDouble time,
int var,
std::string varName,
const NekDouble x2_in,
const NekDouble x3_in)
{
EvaluateBoundaryConditions(time);
EvaluateBoundaryConditions(time, var, varName);
}
boost::shared_ptr<ExpList> &DisContField3DHomogeneous1D::
......
......@@ -118,7 +118,9 @@ namespace Nektar
* evaluated
*/
MULTI_REGIONS_EXPORT void EvaluateBoundaryConditions(
const NekDouble time = 0.0);
const NekDouble time = 0.0,
int var = 0,
std::string varName = "");
inline const Array<OneD,const MultiRegions::ExpListSharedPtr>
&GetBndCondExpansions();
......@@ -271,6 +273,8 @@ namespace Nektar
virtual void v_EvaluateBoundaryConditions(