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

Added: Support for function-defined immersed boundaries

Browse Source
This commit is contained in:
Knut Morten Okstad
2019-09-17 07:56:47 +02:00
parent bcc7754b36
commit 53533aaa97
8 changed files with 146 additions and 47 deletions
Download Patch File Download Diff File Expand all files Collapse all files

14
src/ASM/ASMs2DIB.C
View File

@@ -89,6 +89,16 @@ void ASMs2DIB::addHole (double R, double X1, double Y1, double X2, double Y2)
}
bool ASMs2DIB::setGeometry (RealFunc* f, double power, double threshold)
{
if (myGeometry) return false;
myGeometry = new GeoFunc2D(f,power,threshold);
return true;
}
ElementBlock* ASMs2DIB::immersedGeometry () const
{
if (!myGeometry) return nullptr;
@@ -112,8 +122,8 @@ void ASMs2DIB::getNoIntPoints (size_t& nPt, size_t& nIPt)
if (myGeometry)
{
nPt = firstIp;
for (size_t e = 0; e < quadPoints.size(); e++)
nPt += quadPoints[e].size();
for (const Real2DMat& qp : quadPoints)
nPt += qp.size();
}
#ifdef SP_DEBUG

6
src/ASM/ASMs2DIB.h
View File

@@ -68,6 +68,12 @@ public:
//! \param[in] Y2 Y-coordinate of the second circle centre
virtual void addHole(double R, double X1, double Y1, double X2, double Y2);
//! \brief Defines the immersed geometry from a scalar function.
//! \param[in] f The function to use
//! \param[in] power Exponent to apply on the specified function
//! \param[in] threshold Inside/outside threshold
virtual bool setGeometry(RealFunc* f, double power, double threshold);
//! \brief Returns an additional geometry to visualize (immersed boundaries).
virtual ElementBlock* immersedGeometry() const;

40
src/ASM/LR/ASMu2DIB.C
View File

@@ -32,7 +32,7 @@ ASMu2DIB::ASMu2DIB (const ASMu2DIB& patch, unsigned char n_f)
: ASMu2D(patch,n_f), quadPoints(patch.quadPoints)
{
maxDepth = patch.maxDepth;
myGeometry = nullptr; // because we don't allow multiple patches (yet)
myGeometry = nullptr;
}
@@ -87,31 +87,42 @@ void ASMu2DIB::addHole (double R, double X1, double Y1, double X2, double Y2)
}
bool ASMu2DIB::setGeometry (RealFunc* f, double power, double threshold)
{
if (myGeometry) return false;
myGeometry = new GeoFunc2D(f,power,threshold);
return true;
}
ElementBlock* ASMu2DIB::immersedGeometry () const
{
return myGeometry ? myGeometry->tesselate() : nullptr;
}
void ASMu2DIB::getNoIntPoints (size_t& nPt, size_t&)
void ASMu2DIB::getNoIntPoints (size_t& nPt, size_t& nIPt)
{
nPt = 0;
for (size_t e = 0; e < quadPoints.size(); e++)
nPt += quadPoints[e].size();
firstIp = nPt;
this->ASMbase::getNoIntPoints(nPt,nIPt);
if (myGeometry)
{
nPt = firstIp;
for (const Real2DMat& qp : quadPoints)
nPt += qp.size();
}
}
bool ASMu2DIB::generateFEMTopology ()
{
if (!myGeometry)
{
std::cerr <<" *** ASMu2DIB::generateFEMTopology: No geometry description."
<< std::endl;
return false;
}
if (!this->ASMu2D::generateFEMTopology())
return false;
else if (!myGeometry)
return true;
size_t i, e, n;
const int nBasis = lrspline->nBasisFunctions();
@@ -182,7 +193,10 @@ bool ASMu2DIB::generateFEMTopology ()
bool ASMu2DIB::integrate (Integrand& integrand,
GlobalIntegral& glbInt, const TimeDomain& time)
{
return this->ASMu2D::integrate(integrand,glbInt,time,quadPoints);
if (!myGeometry)
return this->ASMu2D::integrate(integrand,glbInt,time);
return this->integrate(integrand,glbInt,time,quadPoints);
}

8
src/ASM/LR/ASMu2DIB.h
View File

@@ -46,11 +46,17 @@ public:
//! \param[in] Y2 Y-coordinate of the second circle centre
virtual void addHole(double R, double X1, double Y1, double X2, double Y2);
//! \brief Defines the immersed geometry from a scalar function.
//! \param[in] f The function to use
//! \param[in] power Exponent to apply on the specified function
//! \param[in] threshold Inside/outside threshold
virtual bool setGeometry(RealFunc* f, double power, double threshold);
//! \brief Returns an additional geometry to visualize (immersed boundaries).
virtual ElementBlock* immersedGeometry() const;
//! \brief Computes the total number of integration points in this patch.
virtual void getNoIntPoints(size_t& nPt, size_t&);
virtual void getNoIntPoints(size_t& nPt, size_t& nIPt);
//! \brief Generates the finite element topology data for the patch.
//! \details This method is overridden in this class, to include the

4
src/SIM/SIM1D.C
View File

@@ -163,7 +163,7 @@ bool SIM1D::parseGeometryTag (const TiXmlElement* elem)
}
IFEM::cout <<"\tConnecting P"<< slave <<" V"<< sVert
<<" to P"<< master <<" E"<< mVert << std::endl;
<<" to P"<< master <<" V"<< mVert << std::endl;
ASM1D* spch = dynamic_cast<ASM1D*>(myModel[slave-1]);
ASM1D* mpch = dynamic_cast<ASM1D*>(myModel[master-1]);
@@ -177,7 +177,6 @@ bool SIM1D::parseGeometryTag (const TiXmlElement* elem)
else if (!strcasecmp(elem->Value(),"projection") && !isRefined)
{
bool ok = true;
const TiXmlElement* child = elem->FirstChildElement();
if (child && !strncasecmp(child->Value(),"patch",5) && child->FirstChild())
{
@@ -189,6 +188,7 @@ bool SIM1D::parseGeometryTag (const TiXmlElement* elem)
for (ASMbase* pch : myModel)
pch->createProjectionBasis(false);
bool ok = true;
for (int pid = 1; isp->good() && ok; pid++)
{
IFEM::cout <<"\tReading projection basis for patch "<< pid << std::endl;

106
src/SIM/SIM2D.C
View File

@@ -15,6 +15,7 @@
#include "ModelGenerator.h"
#include "ASMs2DC1.h"
#include "ImmersedBoundaries.h"
#include "FieldFunctions.h"
#include "Functions.h"
#include "Utilities.h"
#include "Vec3Oper.h"
@@ -32,6 +33,7 @@ struct Interface
std::pair<ASMs2DC1*,int> master; //!< Patch and edge index of the master
std::pair<ASMs2DC1*,int> slave; //!< Patch and edge index of the slave
bool reversed; //!< Relative orientation toggle
//! \brief Constructor initializing an Interface instance.
Interface(ASMbase* m, int me, ASMbase* s, int se, bool r)
: master(std::make_pair(static_cast<ASMs2DC1*>(m),me)),
@@ -257,7 +259,19 @@ bool SIM2D::parseGeometryTag (const TiXmlElement* elem)
else if (!strcasecmp(elem->Value(),"immersedboundary"))
{
int patch = 0;
utl::getAttribute(elem,"patch",patch);
ASMbase* pch = nullptr;
if (utl::getAttribute(elem,"patch",patch))
pch = this->getPatch(patch,true);
else if (!myModel.empty())
pch = myModel.front();
if (!pch)
{
std::cerr <<" *** SIM2D::parse: Invalid patch index "
<< patch << std::endl;
return false;
}
utl::getAttribute(elem,"stabilization",Immersed::stabilization);
if (Immersed::stabilization != 0)
IFEM::cout <<"\tStabilization option: "<< Immersed::stabilization
@@ -272,10 +286,7 @@ bool SIM2D::parseGeometryTag (const TiXmlElement* elem)
utl::getAttribute(child,"R",R);
utl::getAttribute(child,"Xc",Xc);
utl::getAttribute(child,"Yc",Yc);
if (patch > 0 && patch <= (int)myModel.size())
myModel[patch-1]->addHole(R,Xc,Yc);
else
myModel.front()->addHole(R,Xc,Yc);
pch->addHole(R,Xc,Yc);
}
else if (!strcasecmp(child->Value(),"Oval"))
{
@@ -285,16 +296,52 @@ bool SIM2D::parseGeometryTag (const TiXmlElement* elem)
utl::getAttribute(child,"Y1",Y1);
utl::getAttribute(child,"X2",X2);
utl::getAttribute(child,"Y2",Y2);
if (patch > 0 && patch <= (int)myModel.size())
myModel[patch-1]->addHole(R,X1,Y1,X2,Y2);
pch->addHole(R,X1,Y1,X2,Y2);
}
else if (!strcasecmp(child->Value(),"levelset") && child->FirstChild())
{
double power = 1.0, threshold = 0.5;
std::string type("file");
utl::getAttribute(child,"power",power);
utl::getAttribute(child,"threshold",threshold);
utl::getAttribute(child,"type",type,true);
const char* val = child->FirstChild()->Value();
RealFunc* lfunc = nullptr;
IFEM::cout <<"\tLevel set ("<< type <<")";
if (type == "file")
{
IFEM::cout <<" \""<< val <<"\""<< std::endl;
std::ifstream ibs(val);
if (ibs.good())
lfunc = new FieldFuncStream({pch},ibs);
else
std::cerr <<" *** SIM2D::parse: Failed to open file \""<< val
<<"\""<< std::endl;
}
else
myModel.front()->addHole(R,X1,Y1,X2,Y2);
{
lfunc = utl::parseRealFunc(val,type);
IFEM::cout << std::endl;
}
if (!pch->setGeometry(lfunc,power,threshold))
return false;
if (type == "file")
myAddScalars[val] = lfunc;
else
{
static std::string fname = "level_set0";
++fname[9];
myAddScalars[fname] = lfunc;
}
}
}
else if (!strcasecmp(elem->Value(),"projection") && !isRefined)
{
bool ok = true;
const TiXmlElement* child = elem->FirstChildElement();
if (child && !strncasecmp(child->Value(),"patch",5) && child->FirstChild())
{
@@ -306,6 +353,7 @@ bool SIM2D::parseGeometryTag (const TiXmlElement* elem)
for (ASMbase* pch : myModel)
pch->createProjectionBasis(false);
bool ok = true;
for (int pid = 1; isp->good() && ok; pid++)
{
IFEM::cout <<"\tReading projection basis for patch "<< pid << std::endl;
@@ -642,22 +690,19 @@ bool SIM2D::parse (char* keyWord, std::istream& is)
}
/*!
\brief Local-scope convenience function for error message generation.
*/
static bool constrError (const char* lab, int idx)
{
std::cerr <<" *** SIM2D::addConstraint: Invalid "<< lab << idx << std::endl;
return false;
}
bool SIM2D::addConstraint (int patch, int lndx, int ldim, int dirs, int code,
int& ngnod, char basis)
{
// Lambda function for error message generation
auto&& error = [](const char* message, int idx)
{
std::cerr <<" *** SIM2D::addConstraint: Invalid "
<< message <<" ("<< idx <<")."<< std::endl;
return false;
};
if (patch < 1 || patch > (int)myModel.size())
return constrError("patch index ",patch);
return error("patch index",patch);
bool open = ldim < 0; // open means without its end points
bool project = lndx < -10;
@@ -676,7 +721,7 @@ bool SIM2D::addConstraint (int patch, int lndx, int ldim, int dirs, int code,
// Must dynamic cast here, since ASM2D is not derived from ASMbase
ASM2D* pch = dynamic_cast<ASM2D*>(myModel[patch-1]);
if (!pch) return constrError("2D patch ",patch);
if (!pch) return error("2D patch",patch);
switch (abs(ldim))
{
@@ -689,7 +734,7 @@ bool SIM2D::addConstraint (int patch, int lndx, int ldim, int dirs, int code,
case 4: pch->constrainCorner( 1, 1,dirs,abs(code),basis); break;
default:
IFEM::cout << std::endl;
return constrError("vertex index ",lndx);
return error("vertex index",lndx);
}
break;
@@ -714,7 +759,7 @@ bool SIM2D::addConstraint (int patch, int lndx, int ldim, int dirs, int code,
break;
default:
IFEM::cout << std::endl;
return constrError("edge index ",lndx);
return error("edge index",lndx);
}
break;
@@ -729,7 +774,7 @@ bool SIM2D::addConstraint (int patch, int lndx, int ldim, int dirs, int code,
default:
IFEM::cout << std::endl;
return constrError("local dimension switch ",ldim);
return error("local dimension switch",ldim);
}
return true;
@@ -835,3 +880,14 @@ Vector SIM2D::getSolution (const Vector& psol, double u, double v,
double par[2] = { u, v };
return this->SIMgeneric::getSolution(psol,par,deriv,patch);
}
bool SIM2D::writeAddFuncs (int iStep, int& nBlock, int idBlock, double time)
{
for (const std::pair<std::string,RealFunc*>& func : myAddScalars)
if (!this->writeGlvF(*func.second, func.first.c_str(), iStep,
nBlock, idBlock++, time))
return false;
return this->SIMgeneric::writeAddFuncs(iStep,nBlock,idBlock,time);
}

6
src/SIM/SIM2D.h
View File

@@ -133,6 +133,12 @@ protected:
virtual ASMbase* readPatch(std::istream& isp, int pchInd, const CharVec& unf,
const char* whiteSpace) const;
//! \brief Writes out the additional functions to VTF-file.
virtual bool writeAddFuncs(int iStep, int& nBlock, int idBlock, double time);
private:
std::map<std::string,RealFunc*> myAddScalars; //!< Scalar functions to output
protected:
CharVec nf; //!< Number of scalar fields
bool checkRHSys; //!< Check if all patches are in a right-hand system

9
src/SIM/SIMoutput.h
View File

@@ -310,6 +310,9 @@ protected:
//! \brief Preprocesses the result sampling points.
virtual void preprocessResultPoints();
//! \brief Writes out the additional functions to VTF-file.
virtual bool writeAddFuncs(int iStep, int& nBlock, int idBlock, double time);
private:
//! \brief Private helper to extract patch-level solution vectors.
bool extractNodeVec(const Vector& glbVec, Vector& locVec,
@@ -321,9 +324,6 @@ private:
int& nBlock, std::vector< std::vector<int> >& sID,
size_t* i = nullptr);
//! \brief Private helper to write additional functions to VTF-file.
bool writeAddFuncs(int iStep, int& nBlock, int idBlock, double time);
//! \brief Struct defining a result sampling point.
struct ResultPoint
{
@@ -332,7 +332,8 @@ private:
int inod; //!< Local node number of the closest node
double u[3]; //!< Parameters of the point (u,v,w)
Vec3 X; //!< Spatial coordinates of the point
// \brief Default constructor.
//! \brief Default constructor.
ResultPoint() : npar(0), patch(1), inod(0) { u[0] = u[1] = u[2] = 0.0; }
};