OilfieldToolsNet/IFEM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

added: support FieldFunctions for mixed spline fields

Browse Source 
If no dash is found in the basis name, we assume the field
names given are on separate bases, one per basis

this is used to instance fieldfunctions for div-compatible
solution fields
This commit is contained in:
Arne Morten Kvarving 2022-03-12 21:52:53 +01:00
parent 3dc01f07a4
commit 9bfe8af84e
8 changed files with 127 additions and 44 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

129
src/Utility/FieldFunctions.C
View File

@ -103,8 +103,12 @@ bool FieldFuncHDF5::load (const std::vector<std::string>& fieldNames,
size_t nOK = 0;
size_t nPatches = 0;
#ifdef HAS_HDF5
bool mixedField = basisName.find('-') == std::string::npos;
std::stringstream str;
str << level << "/" << basisName << "/fields/" << fieldNames.front();
str << level << "/" << basisName;
if (mixedField)
str << "-1";
str << "/fields/" << fieldNames.front();
nPatches = hdf5->getFieldSize(str.str());
#endif
if (nPatches == 0)
@ -120,39 +124,63 @@ bool FieldFuncHDF5::load (const std::vector<std::string>& fieldNames,
for (size_t ip = 0; ip < nPatches; ip++)
{
std::stringstream str;
str << level << "/" << basisName << "/basis";
str << level << "/" << basisName;
if (mixedField)
str << "-1";
str << "/basis";
if (hdf5->getFieldSize(str.str()))
{
if (patch[ip])
{
// We have an updated basis at this level, replace current one
if (patch[ip]->getNoParamDim() == 2) nFldC2D = patch[ip]->getNoFields();
if (patch[ip]->getNoParamDim() == 3) nFldC3D = patch[ip]->getNoFields();
delete patch[ip];
if (mixedField) {
std::vector<unsigned char> nFld(fieldNames.size(), 1);
for (size_t i = 1; i <= fieldNames.size(); ++i) {
std::string g2;
std::stringstream sbasis;
sbasis << level << "/" << basisName << '-' << i << "/basis/"<< ip+1;
hdf5->readString(sbasis.str(),g2);
if (i == 1) {
if (g2.compare(0,9,"200 1 0 0") == 0)
patch[ip] = ASM2D::create(ASM::Spline,2,nFld);
else if (g2.compare(0,9,"700 1 0 0") == 0)
patch[ip] = ASM3D::create(ASM::Spline,nFld);
else if (g2.compare(0,18,"# LRSPLINE SURFACE") == 0)
patch[ip] = ASM2D::create(ASM::LRSpline,2,nFld);
else if (g2.compare(0,17,"# LRSPLINE VOLUME") == 0)
patch[ip] = ASM3D::create(ASM::LRSpline,nFld);
}
std::stringstream strg2(g2);
patch[ip]->read(strg2, i);
}
} else {
if (patch[ip])
{
// We have an updated basis at this level, replace current one
if (patch[ip]->getNoParamDim() == 2) nFldC2D = patch[ip]->getNoFields();
if (patch[ip]->getNoParamDim() == 3) nFldC3D = patch[ip]->getNoFields();
delete patch[ip];
}
std::string g2;
std::stringstream sbasis;
sbasis << level << "/" << basisName << "/basis/"<< ip+1;
hdf5->readString(sbasis.str(),g2);
if (g2.compare(0,9,"200 1 0 0") == 0)
patch[ip] = ASM2D::create(ASM::Spline,nFldC2D);
else if (g2.compare(0,9,"700 1 0 0") == 0)
patch[ip] = ASM3D::create(ASM::Spline,nFldC3D);
else if (g2.compare(0,18,"# LRSPLINE SURFACE") == 0)
patch[ip] = ASM2D::create(ASM::LRSpline,nFldC2D);
else if (g2.compare(0,17,"# LRSPLINE VOLUME") == 0)
patch[ip] = ASM3D::create(ASM::LRSpline,nFldC3D);
else
patch[ip] = nullptr;
if (patch[ip])
{
std::stringstream strg2(g2);
patch[ip]->read(strg2);
}
else
std::cerr <<" *** FieldFuncHDF5::load: Undefined basis "<< sbasis.str()
<<" ("<< g2.substr(0,9) <<")"<< std::endl;
}
std::string g2;
std::stringstream sbasis;
sbasis << level << "/" << basisName << "/basis/"<< ip+1;
hdf5->readString(sbasis.str(),g2);
if (g2.compare(0,9,"200 1 0 0") == 0)
patch[ip] = ASM2D::create(ASM::Spline,nFldC2D);
else if (g2.compare(0,9,"700 1 0 0") == 0)
patch[ip] = ASM3D::create(ASM::Spline,nFldC3D);
else if (g2.compare(0,18,"# LRSPLINE SURFACE") == 0)
patch[ip] = ASM2D::create(ASM::LRSpline,nFldC2D);
else if (g2.compare(0,17,"# LRSPLINE VOLUME") == 0)
patch[ip] = ASM3D::create(ASM::LRSpline,nFldC3D);
else
patch[ip] = nullptr;
if (patch[ip])
{
std::stringstream strg2(g2);
patch[ip]->read(strg2);
}
else
std::cerr <<" *** FieldFuncHDF5::load: Undefined basis "<< sbasis.str()
<<" ("<< g2.substr(0,9) <<")"<< std::endl;
}
if (patch[ip])
@ -161,7 +189,10 @@ bool FieldFuncHDF5::load (const std::vector<std::string>& fieldNames,
for (size_t i = 0; i < nFldCmp; i++)
{
std::stringstream str;
str << level << "/" << basisName << "/fields/" << fieldNames[i] << "/" << ip+1;
str << level << "/" << basisName;
if (mixedField)
str << '-' << i+1;
str << "/fields/" << fieldNames[i] << "/" << ip+1;
hdf5->readVector(str.str(),coefs[i]);
#if SP_DEBUG > 1
std::cout <<"FieldFuncHDF5::load: Reading \""<< fieldNames[i]
@ -174,14 +205,25 @@ bool FieldFuncHDF5::load (const std::vector<std::string>& fieldNames,
if (nFldCmp > 1)
{
RealArray coef1;
coef1.reserve(nFldCmp*coefs.front().size());
for (size_t i = 0; i < coefs.front().size(); i++)
for (size_t j = 0; j < nFldCmp; j++)
coef1.push_back(coefs[j][i]);
this->addPatchField(patch[ip],coef1);
int basis, nf;
if (mixedField) {
for (const RealArray& v : coefs)
std::copy(v.begin(), v.end(), std::back_inserter(coef1));
basis = nFldCmp == 2 ? 12 : 123;
nf = coefs.size();
} else {
coef1.reserve(nFldCmp*coefs.front().size());
for (size_t i = 0; i < coefs.front().size(); i++)
for (size_t j = 0; j < nFldCmp; j++)
coef1.push_back(coefs[j][i]);
basis = 1;
nf = patch[ip]->getNoFields(1);
}
this->addPatchField(patch[ip],coef1,nf,basis);
}
else
this->addPatchField(patch[ip],coefs.front());
this->addPatchField(patch[ip],coefs.front(),
patch[ip]->getNoFields(1),1);
nOK++;
}
@ -215,7 +257,9 @@ void FieldFunction::clearField ()
}
void FieldFunction::addPatchField (ASMbase* pch, const RealArray& coefs)
void FieldFunction::addPatchField (ASMbase* pch,
const RealArray& coefs,
int, int)
{
field.push_back(Field::create(pch,coefs));
npch = field.size();
@ -298,9 +342,12 @@ void FieldsFuncBase::clearField ()
}
void FieldsFuncBase::addPatchField (ASMbase* pch, const RealArray& coefs)
void FieldsFuncBase::addPatchField (ASMbase* pch,
const RealArray& coefs,
int nf,
int basis)
{
field.push_back(Fields::create(pch,coefs,1,pch->getNoFields(1)));
field.push_back(Fields::create(pch,coefs,basis,nf));
npch = field.size();
}

16
src/Utility/FieldFunctions.h
View File

@ -79,7 +79,11 @@ protected:
//! \brief Adds a patch-wise field with the given coefficient values.
//! \param[in] pch The patch to define the field over
//! \param[in] coefs Field values
virtual void addPatchField(ASMbase* pch, const std::vector<Real>& coefs) = 0;
//! \param[in] nf Number of field components
//! \param[in] basis Basis to use
virtual void addPatchField(ASMbase* pch,
const std::vector<Real>& coefs,
int nf, int basis) = 0;
//! \brief Clears the field container.
virtual void clearField() = 0;
@ -121,7 +125,9 @@ protected:
//! \brief Adds a patch-wise field with the given coefficient values.
//! \param[in] pch The patch to define the field over
//! \param[in] coefs Field values
virtual void addPatchField(ASMbase* pch, const std::vector<Real>& coefs);
virtual void addPatchField(ASMbase* pch,
const std::vector<Real>& coefs,
int, int);
//! \brief Clears the field container.
virtual void clearField();
@ -184,7 +190,11 @@ protected:
//! \brief Adds a patch-wise field with the given coefficient values.
//! \param[in] pch The patch to define the field over
//! \param[in] coefs Field values
virtual void addPatchField(ASMbase* pch, const std::vector<Real>& coefs);
//! \param[in] nf Number of field components
//! \param[in] basis Basis to use
virtual void addPatchField(ASMbase* pch,
const std::vector<Real>& coefs,
int nf, int basis);
//! \brief Clears the field container.
virtual void clearField();

13
src/Utility/Test/TestFieldFunctions.C
View File

@ -23,6 +23,9 @@ TEST(TestFieldFunctions, 2D1P)
VecFieldFunction f2D_vec("src/Utility/Test/refdata/Field2D-1P",
"Stokes-1", "u", 0);
VecFieldFunction f2Dmx_vec("src/Utility/Test/refdata/Field2D-1Pmx",
"Stokes", "u_x|u_y", 0);
TensorFieldFunction f2D_ten("src/Utility/Test/refdata/Field2D-1P",
"Stokes-1", "u_x,x|u_x,y|u_y,x|u_y,y", 0);
@ -37,6 +40,9 @@ TEST(TestFieldFunctions, 2D1P)
Vec3 vec = f2D_vec(X);
EXPECT_NEAR(vec[0], x, 1e-14);
EXPECT_NEAR(vec[1], -y, 1e-14);
Vec3 vecmx = f2Dmx_vec(X);
EXPECT_NEAR(vecmx[0], x, 1e-14);
EXPECT_NEAR(vecmx[1], -y, 1e-14);
Tensor ten = f2D_ten(X); // see above
EXPECT_NEAR(ten(1,1), 1.0, 1e-14);
EXPECT_NEAR(ten(1,2), 0.0, 1e-14);
@ -89,6 +95,9 @@ TEST(TestFieldFunctions, 3D1P)
VecFieldFunction f3D_vec("src/Utility/Test/refdata/Field3D-1P",
"Stokes-1", "u", 0);
VecFieldFunction f3Dmx_vec("src/Utility/Test/refdata/Field3D-1Pmx",
"Stokes", "u_x|u_y|u_z", 0);
TensorFieldFunction f3D_ten("src/Utility/Test/refdata/Field3D-1P",
"Stokes-1", "u_x,x|u_y,x|u_z,x|u_x,y|u_y,y|u_z,y|u_x,z|u_x,y|u_z,z", 0);
@ -106,6 +115,10 @@ TEST(TestFieldFunctions, 3D1P)
EXPECT_NEAR(vec[0], y*(1.0-y), 1e-14); // y*(1-y)
EXPECT_NEAR(vec[1], z*(1.0-z), 1e-14); // z*(1-z)
EXPECT_NEAR(vec[2], x*(1.0-x), 1e-14); // x*(1-x)
Vec3 vecmx = f3Dmx_vec(X);
EXPECT_NEAR(vecmx[0], x, 1e-13);
EXPECT_NEAR(vecmx[1], y, 1e-13);
EXPECT_NEAR(vecmx[2], -2*z, 1e-13);
param[0] = param[1] = param[2] = 0.5;
Tensor ten = f3D_ten(X);
x = y = z = 1.0;

13
src/Utility/Test/TestFieldFunctionsLR.C
View File

@ -23,6 +23,9 @@ TEST(TestFieldFunctions, 2D1PLR)
VecFieldFunction f2D_vec("src/Utility/Test/refdata/Field2D-1PLR",
"Stokes-1", "u", 0);
VecFieldFunction f2Dmx_vec("src/Utility/Test/refdata/Field2D-1PLRmx",
"Stokes", "u_x|u_y", 0);
TensorFieldFunction f2D_ten("src/Utility/Test/refdata/Field2D-1PLR",
"Stokes-1", "u_x,x|u_x,y|u_y,x|u_y,y", 0);
@ -37,6 +40,9 @@ TEST(TestFieldFunctions, 2D1PLR)
Vec3 vec = f2D_vec(X);
EXPECT_NEAR(vec[0], x, 1e-14);
EXPECT_NEAR(vec[1], -y, 1e-14);
Vec3 vecmx = f2Dmx_vec(X);
EXPECT_NEAR(vecmx[0], x, 1e-14);
EXPECT_NEAR(vecmx[1], -y, 1e-14);
Tensor ten = f2D_ten(X); // see above
EXPECT_NEAR(ten(1,1), 1.0, 1e-14);
EXPECT_NEAR(ten(1,2), 0.0, 1e-13);
@ -89,6 +95,9 @@ TEST(TestFieldFunctions, 3D1PLR)
VecFieldFunction f3D_vec("src/Utility/Test/refdata/Field3D-1PLR",
"Stokes-1", "u", 0);
VecFieldFunction f3Dmx_vec("src/Utility/Test/refdata/Field3D-1PLRmx",
"Stokes", "u_x|u_y|u_z", 0);
TensorFieldFunction f3D_ten("src/Utility/Test/refdata/Field3D-1PLR",
"Stokes-1", "u_x,x|u_y,x|u_z,x|u_x,y|u_y,y|u_z,y|u_x,z|u_x,y|u_z,z", 0);
@ -106,6 +115,10 @@ TEST(TestFieldFunctions, 3D1PLR)
EXPECT_NEAR(vec[0], y*(1.0-y), 1e-14); // y*(1-y)
EXPECT_NEAR(vec[1], z*(1.0-z), 1e-14); // z*(1-z)
EXPECT_NEAR(vec[2], x*(1.0-x), 1e-14); // x*(1-x)
Vec3 vecmx = f3Dmx_vec(X);
EXPECT_NEAR(vecmx[0], x, 1e-13);
EXPECT_NEAR(vecmx[1], y, 1e-13);
EXPECT_NEAR(vecmx[2], -2*z, 1e-13);
param[0] = param[1] = param[2] = 0.5;
Tensor ten = f3D_ten(X);
x = y = z = 1.0;

BIN
src/Utility/Test/refdata/Field2D-1PLRmx.hdf5 Normal file
View File

Binary file not shown.

BIN
src/Utility/Test/refdata/Field2D-1Pmx.hdf5 Normal file
View File

Binary file not shown.

BIN
src/Utility/Test/refdata/Field3D-1PLRmx.hdf5 Normal file
View File

Binary file not shown.

BIN
src/Utility/Test/refdata/Field3D-1Pmx.hdf5 Normal file
View File

Binary file not shown.