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

Added: Output of extraction function value to VTF as element result

Browse Source
This commit is contained in:
Knut Morten Okstad 2019-05-15 19:08:52 +02:00
parent f08a50ebe9
commit e531d64d6e
3 changed files with 87 additions and 48 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

82
src/SIM/SIMoutput.C
View File

@ -1115,6 +1115,14 @@ bool SIMoutput::writeGlvM (const Mode& mode, bool freq, int& nBlock)
}
/*!
If \a dualPrefix is not null, it is assumed that the \a norms were
computed from the dual solution and therefore labelled as such.
In addition to the \a norms, also the \a dualField function value is written
as a scalar field to the VTF-file in that case, evaluated at the centre of
each visualization element since this function is typically discontinuous.
*/
bool SIMoutput::writeGlvN (const Matrix& norms, int iStep, int& nBlock,
const std::vector<std::string>& prefix,
int idBlock, const char* dualPrefix)
@ -1137,63 +1145,73 @@ bool SIMoutput::writeGlvN (const Matrix& norms, int iStep, int& nBlock,
return norm->hasElementContributions(iGroup,iNorm);
};
size_t idxW = 0;
size_t nrow = norms.rows();
if (dualField && dualPrefix) idxW = ++nrow;
std::vector<IntVec> sID(nrow);
Matrix field;
std::vector<IntVec> sID(norms.rows());
size_t i, j, k, l, m;
size_t i, j, k, l;
int geomID = myGeomID;
for (i = 0; i < myModel.size(); i++)
for (const ASMbase* pch : myModel)
{
if (myModel[i]->empty()) continue; // skip empty patches
if (pch->empty()) continue; // skip empty patches
if (msgLevel > 1)
IFEM::cout <<"Writing element norms for patch "<< i+1 << std::endl;
IFEM::cout <<"Writing element norms for patch "
<< pch->idx+1 << std::endl;
const ElementBlock* grid = myVtf->getBlock(++geomID);
myModel[i]->extractElmRes(norms,field);
pch->extractElmRes(norms,field);
if (grid->getNoElms() > field.cols())
size_t ncol = grid->getNoElms();
if (ncol > field.cols() || nrow > field.rows())
{
// Expand the element result array
Matrix efield(field);
field.resize(field.rows(),grid->getNoElms());
for (j = 1; j <= field.cols(); j++)
field.resize(nrow,ncol);
for (j = 1; j <= ncol; j++)
field.fillColumn(j,efield.getColumn(grid->getElmId(j)));
if (idxW && dualField->initPatch(pch->idx))
for (j = 1; j <= ncol; j++)
field(idxW,j) = dualField->getScalarValue(grid->getCenter(j));
}
j = l = 1;
for (k = m = 0; m < field.rows(); m++)
for (k = 0, i = j = l = 1; i <= nrow; i++, l++)
{
if (l > norm->getNoFields(j))
l = 1, ++j;
if (writeNorm(j,l++))
if (!myVtf->writeEres(field.getRow(1+m),++nBlock,geomID))
if (writeNorm(j,l) || i == idxW)
{
if (!myVtf->writeEres(field.getRow(i),++nBlock,geomID))
return false;
else
sID[k++].push_back(nBlock);
sID[k++].push_back(nBlock);
}
}
}
std::string normName;
j = l = 1;
for (k = 0; k < sID.size() && !sID[k].empty(); l++)
for (k = 0, i = j = l = 1; k < sID.size() && !sID[k].empty(); i++, l++)
{
if (l > norm->getNoFields(j))
l = 1, ++j;
if (!writeNorm(j,l))
continue;
if (writeNorm(j,l) || i == idxW)
{
if (i == idxW && dualPrefix)
normName = std::string(dualPrefix) + " extraction function";
else if (j == 1 && dualPrefix)
normName = norm->getName(j,l,dualPrefix);
else if (j > 1 && j-2 < prefix.size())
normName = norm->getName(j,l,prefix[j-2].c_str());
else
normName = norm->getName(j,l);
if (j == 1 && dualPrefix)
normName = norm->getName(j,l,dualPrefix);
else if (j > 1 && j-2 < prefix.size())
normName = norm->getName(j,l,prefix[j-2].c_str());
else
normName = norm->getName(j,l);
if (!myVtf->writeSblk(sID[k++],normName.c_str(),++idBlock,iStep,true))
return false;
if (!myVtf->writeSblk(sID[k++],normName.c_str(),++idBlock,iStep,true))
return false;
}
}
delete norm;
@ -1215,15 +1233,15 @@ bool SIMoutput::writeGlvE (const Vector& vec, int iStep, int& nBlock,
IntVec sID;
int geomID = myGeomID;
for (size_t i = 0; i < myModel.size(); i++)
for (const ASMbase* pch : myModel)
{
if (myModel[i]->empty()) continue; // skip empty patches
if (pch->empty()) continue; // skip empty patches
if (msgLevel > 1)
IFEM::cout <<"Writing element field "<< name
<<" for patch "<< i+1 << std::endl;
<<" for patch "<< pch->idx+1 << std::endl;
myModel[i]->extractElmRes(infield,field);
pch->extractElmRes(infield,field);
if (!myVtf->writeEres(field.getRow(1),++nBlock,++geomID))
return false;

42
src/Utility/ElementBlock.C
View File

@ -33,7 +33,7 @@ ElementBlock::ElementBlock (size_t nenod)
void ElementBlock::resize (size_t nI, size_t nJ, size_t nK)
{
coord.resize(nI*nJ*nK);
params.resize(nI*nJ*nK);
param.resize(nI*nJ*nK);
if (nen == 2 && nJ < 2 && nK < 2)
MMNPC.resize(2*(nI-1));
else if (nen == 3 && nK < 2)
@ -57,7 +57,7 @@ void ElementBlock::resize (size_t nI, size_t nJ, size_t nK)
void ElementBlock::unStructResize (size_t nEl, size_t nPts)
{
coord.resize(nPts);
params.resize(nPts);
param.resize(nPts);
MMNPC.resize(nen*nEl);
MINEX.resize(MMNPC.size()/nen,0);
std::iota(MINEX.begin(),MINEX.end(),1);
@ -84,9 +84,11 @@ bool ElementBlock::setCoor (size_t i, const Vec3& X)
bool ElementBlock::setParams (size_t i, Real u, Real v, Real w)
{
if (i >= params.size()) return false;
if (i >= param.size()) return false;
params[i] = {{u,v,w}};
param[i][0] = u;
param[i][1] = v;
param[i][2] = w;
return true;
}
@ -116,21 +118,35 @@ bool ElementBlock::addLine (Real x1, Real y1, Real z1,
void ElementBlock::merge (const ElementBlock* other, std::vector<int>& nodeNums)
{
nodeNums.resize(other->coord.size());
nodeNums.clear();
nodeNums.reserve(other->coord.size());
size_t i;
std::vector<Vec3>::const_iterator cit;
for (i = 0; i < nodeNums.size(); i++)
if ((cit = find(coord.begin(),coord.end(),other->coord[i])) == coord.end())
for (const Vec3& X : other->coord)
if ((cit = find(coord.begin(),coord.end(),X)) == coord.end())
{
nodeNums[i] = coord.size();
coord.push_back(other->coord[i]);
nodeNums.push_back(coord.size());
coord.push_back(X);
}
else
nodeNums[i] = cit - coord.begin();
nodeNums.push_back(cit - coord.begin());
for (i = 0; i < other->MMNPC.size(); i++)
MMNPC.push_back(nodeNums[other->MMNPC[i]]);
for (int inod : other->MMNPC)
MMNPC.push_back(nodeNums[inod]);
MINEX.insert(MINEX.end(),other->MINEX.begin(),other->MINEX.end());
}
Vec3 ElementBlock::getCenter (size_t i) const
{
if (i < 1 || i > MINEX.size())
return Vec3();
Vec3 XC;
for (size_t j = 0; j < nen; j++)
XC += coord[MMNPC[nen*(i-1)+j]];
XC /= nen;
return XC;
}

11
src/Utility/ElementBlock.h
View File

@ -77,21 +77,26 @@ public:
std::vector<Vec3>::const_iterator begin_XYZ() const { return coord.begin(); }
//! \brief Returns the end of the coord array.
std::vector<Vec3>::const_iterator end_XYZ() const { return coord.end(); }
//! \brief Returns the coordinate of a given node.
const Vec3& getCoord(size_t i) const { return coord[i]; }
//! \brief Returns a pointer to the parameter values of a given node.
const Real* getParam(size_t i) const { return params[i].data(); }
const Real* getParam(size_t i) const { return param[i].data(); }
//! \brief Returns a pointer to the element connectivity array.
const int* getElements() const { return MMNPC.data(); }
//! \brief Returns the coordinates of the center of the given elemment.
Vec3 getCenter(size_t i) const;
private:
typedef std::array<Real,3> Prm3; //!< Convenience type
std::vector<Vec3> coord; //!< Vector of nodal coordinates
std::vector<Prm3> param; //!< Vector of parameter values of the nodal points
std::vector<int> MMNPC; //!< Matrix of Matrices of Nodal Point Correspondance
std::vector<int> MINEX; //!< Matrix of Internal to External element numbers
size_t nen; //!< Number of Element Nodes
std::vector<std::array<Real,3>> params; //!< Parameter values where relevant
};
#endif