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

Fixed: Handling of chained multi-point constraints.

Browse Source 
For problems with time-varying dirichlet conditions,
the constraint equations were not properly updated
between the increments. This was due to the resolving
where we lost the connection between the actual master dof
beeing updated, and the slave in the other end of the chain.

The update is now done in a recursive manner instead, using
a pointer to the next MPC in the chain stored on the DOF level.

Notice that this fix also affected the FSI-slip regression test.

git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@2354 e10b68d5-8a6e-419e-a041-bce267b0401d
This commit is contained in:
kmo
2013-05-02 13:16:41 +00:00
committed by Knut Morten Okstad
parent b6b1877358
commit 8c0579dcdc
11 changed files with 215 additions and 54 deletions
Download Patch File Download Diff File Expand all files Collapse all files

53
src/ASM/ASMbase.C
View File

@@ -611,7 +611,18 @@ void ASMbase::mergeAndGetAllMPCs (const ASMVec& model, MPCSet& allMPCs)
}
void ASMbase::resolveMPCchains (const MPCSet& allMPCs)
/*!
If \a setPtrOnly is \e true, the MPC equations are not modified. Instead
the pointers to the next MPC in the chain is assigned for the master DOFs
which are slaves in other MPCs.
This is used in time-dependent/nonlinear simulators where the constraint
coefficients might be updated due to time variation. The resolving is then
done directly in the MMCEQ/TTCC arrays of the SAM data object.
\sa SAMpatch::updateConstraintEqs.
*/
void ASMbase::resolveMPCchains (const MPCSet& allMPCs, bool setPtrOnly)
{
#if SP_DEBUG > 1
std::cout <<"\nResolving MPC chains"<< std::endl;
@@ -620,16 +631,27 @@ void ASMbase::resolveMPCchains (const MPCSet& allMPCs)
int nresolved = 0;
for (MPCIter cit = allMPCs.begin(); cit != allMPCs.end(); cit++)
if (ASMbase::resolveMPCchain(allMPCs,*cit)) nresolved++;
if (setPtrOnly)
for (size_t i = 0; i < (*cit)->getNoMaster(); i++)
{
MPC master((*cit)->getMaster(i).node,(*cit)->getMaster(i).dof);
MPCIter chain = allMPCs.find(&master);
if (chain != allMPCs.end())
(*cit)->updateMaster(i,*chain); // Set pointer to next MPC in chain
}
else if (ASMbase::resolveMPCchain(allMPCs,*cit))
nresolved++;
if (nresolved > 0)
std::cout <<"Resolved "<< nresolved <<" MPC chains."<< std::endl;
}
// Recursive method to resolve (possibly multi-level) chaining in multi-point
// constraint equations (MPCs). If a master dof in one MPC is specified as a
// slave by another MPC, it is replaced by the master(s) of that other equation.
/*!
Recursive resolving of (possibly multi-level) chaining in multi-point
constraint equations (MPCs). If a master dof in one MPC is specified as a
slave by another MPC, it is replaced by the master(s) of that other equation.
*/
bool ASMbase::resolveMPCchain (const MPCSet& allMPCs, MPC* mpc)
{
@@ -672,6 +694,27 @@ bool ASMbase::resolveMPCchain (const MPCSet& allMPCs, MPC* mpc)
}
bool ASMbase::hasTimeDependentDirichlet (const std::map<int,RealFunc*>& func,
const std::map<int,VecFunc*>& vfunc)
{
std::map<int,RealFunc*>::const_iterator fit;
std::map<int,VecFunc*>::const_iterator vfit;
for (MPCMap::iterator cit = dCode.begin(); cit != dCode.end(); cit++)
if ((fit = func.find(cit->second)) != func.end())
{
if (!fit->second->isConstant())
return true;
}
else if ((vfit = vfunc.find(cit->second)) != vfunc.end())
{
if (!vfit->second->isConstant())
return true;
}
return false;
}
bool ASMbase::updateDirichlet (const std::map<int,RealFunc*>& func,
const std::map<int,VecFunc*>& vfunc, double time)
{

9
src/ASM/ASMbase.h
View File

@@ -266,16 +266,23 @@ public:
//! \brief Resolves (possibly multi-level) chaining in MPC equations.
//! \param[in] allMPCs All multi-point constraint equations in the model
//! \param[in] setPtrOnly If \e true, only set pointer to next MPC in chain
//!
//! \details If a master DOF in one MPC (multi-point constraint) equation
//! is specified as slave by another MPC, it is replaced by the master(s) of
//! that other equation. Since an MPC-equation may couple nodes belonging to
//! different patches, this method must have access to all patches.
static void resolveMPCchains(const MPCSet& allMPCs);
static void resolveMPCchains(const MPCSet& allMPCs, bool setPtrOnly = false);
//! \brief Initializes the multi-point constraint coefficients.
virtual bool initConstraints() { return true;}
//! \brief Checks for time-dependent in-homogeneous Dirichlet conditions.
//! \param[in] func Scalar property fields
//! \param[in] vfunc Vector property fields
bool hasTimeDependentDirichlet (const std::map<int,RealFunc*>& func,
const std::map<int,VecFunc*>& vfunc);
//! \brief Updates the time-dependent in-homogeneous Dirichlet coefficients.
//! \param[in] func Scalar property fields
//! \param[in] vfunc Vector property fields

115
src/ASM/SAMpatch.C
View File

@@ -13,7 +13,6 @@
#include "SAMpatch.h"
#include "ASMbase.h"
#include "MPC.h"
bool SAMpatch::init (const ASMVec& model, int numNod)
@@ -216,7 +215,7 @@ bool SAMpatch::initConstraintEqs (const ASMVec& model)
MPCIter cit;
for (j = 0; j < model.size(); j++)
for (cit = model[j]->begin_MPC(); cit != model[j]->end_MPC(); cit++)
nmmceq += 1 + (*cit)->getNoMaster();
nmmceq += 1 + (*cit)->getNoMaster(true);
// Initialize the constraint equation arrays
mpmceq = new int[nceq+1];
@@ -243,14 +242,16 @@ bool SAMpatch::initConstraintEqs (const ASMVec& model)
}
else if (msc[idof-1] < 0)
{
std::cerr <<"SAM: Ignoring constraint equation for dof "
<< idof <<" ("<< (*cit)->getSlave()
<<").\n This dof is already marked as SLAVE."<< std::endl;
std::cerr <<"SAM: Ignoring constraint equation for dof "<< idof
<<" ("<< (*cit)->getSlave() <<").\n"
<<" This dof is already marked as SLAVE."<< std::endl;
ip--;
nceq--;
continue;
}
(*cit)->iceq = ip-1; // index into mpmceq for this MPC equation
int ipslv = (mpmceq[ip]++) - 1;
mmceq[ipslv] = idof;
ttcc[ipslv] = (*cit)->getSlave().coeff;
@@ -258,29 +259,17 @@ bool SAMpatch::initConstraintEqs (const ASMVec& model)
// Master dofs ...
for (size_t i = 0; i < (*cit)->getNoMaster(); i++)
{
idof = madof[(*cit)->getMaster(i).node-1] + (*cit)->getMaster(i).dof-1;
if (msc[idof-1] > 0)
if (!this->initConstraintEqMaster((*cit)->getMaster(i),
(*cit)->getSlave(),
ttcc[ipslv],ip))
{
int ipmst = (mpmceq[ip]++) - 1;
mmceq[ipmst] = idof;
ttcc[ipmst] = (*cit)->getMaster(i).coeff;
}
else if (msc[idof-1] < 0)
{
// Master dof is constrained (unresolved chaining)
std::cerr <<"SAM: Chained MPCs detected"
<<", slave "<< (*cit)->getSlave()
<<", master "<< (*cit)->getMaster(i)
<<" (ignored)."<< std::endl;
// Something is wrong, ignore this constraint equation
(*cit)->iceq = -1;
mpmceq[ip] = mpmceq[ip-1];
ip--;
nceq--;
break;
}
}
(*cit)->iceq = ip-1; // index into mpmceq for this MPC equation
}
// Reset the negative values in msc before calling SYSEQ
@@ -291,6 +280,40 @@ bool SAMpatch::initConstraintEqs (const ASMVec& model)
}
bool SAMpatch::initConstraintEqMaster (const MPC::DOF& master,
const MPC::DOF& slave,
Real& offset, int ip, Real scale)
{
int idof = madof[master.node-1] + master.dof - 1;
if (msc[idof-1] > 0)
{
int ipmst = (mpmceq[ip]++) - 1;
mmceq[ipmst] = idof;
ttcc[ipmst] = master.coeff*scale;
}
else if (msc[idof-1] < 0)
// This master dof is constrained (unresolved chaining)
if (!master.nextc)
{
std::cerr <<" SAM: Chained MPCs detected, slave "<< slave
<<", master "<< master <<" (ignored)."<< std::endl;
return false;
}
else
{
scale *= master.coeff;
offset += master.nextc->getSlave().coeff*scale;
for (size_t i = 0; i < master.nextc->getNoMaster(); i++)
if (!this->initConstraintEqMaster(master.nextc->getMaster(i),
master.nextc->getSlave(),
offset,ip,scale))
return false;
}
return true;
}
bool SAMpatch::updateConstraintEqs (const ASMVec& model, const Vector* prevSol)
{
if (nceq < 1) return true; // No constraints in this model
@@ -306,24 +329,48 @@ bool SAMpatch::updateConstraintEqs (const ASMVec& model, const Vector* prevSol)
int ipeq = mpmceq[(*cit)->iceq] - 1;
if (msc[idof-1] > 0 || mmceq[ipeq] != idof)
{
std::cerr <<" *** Corrupted SAM arrays detected in update."<< std::endl;
std::cerr <<" *** SAM: Failed to update constraint equations from "
<< *cit << std::endl;
return false;
}
else if (!prevSol)
ttcc[ipeq] = 0.0;
else if (idof <= (int)prevSol->size())
ttcc[ipeq] = (*cit)->getSlave().coeff - (*prevSol)(idof);
else
ttcc[ipeq] = (*cit)->getSlave().coeff;
int jpeq = ipeq;
Real c0 = (*cit)->getSlave().coeff;
// Master dofs ...
for (size_t i = 0; prevSol && i < (*cit)->getNoMaster(); i++)
{
idof = madof[(*cit)->getMaster(i).node-1] + (*cit)->getMaster(i).dof-1;
if (msc[idof-1] > 0 && mmceq[++ipeq] == idof)
ttcc[ipeq] = (*cit)->getMaster(i).coeff;
}
this->updateConstraintEqMaster((*cit)->getMaster(i),c0,ipeq);
// Update the constant term of the constraint equation
if (!prevSol)
ttcc[jpeq] = 0.0;
else if (idof <= (int)prevSol->size())
ttcc[jpeq] = c0 - (*prevSol)(idof);
else
ttcc[jpeq] = c0;
#if SP_DEBUG > 1
std::cout <<"Updated constraint value for dof="<< idof <<": "<< ttcc[jpeq]
<<" from "<< **cit;
#endif
}
return true;
}
void SAMpatch::updateConstraintEqMaster (const MPC::DOF& master,
Real& offset, int& ipeq, Real scale)
{
int idof = madof[master.node-1] + master.dof - 1;
if (msc[idof-1] > 0 && mmceq[++ipeq] == idof)
ttcc[ipeq] = master.coeff*scale;
else if (master.nextc)
{
scale *= master.coeff;
offset += master.nextc->getSlave().coeff*scale;
for (size_t i = 0; i < master.nextc->getNoMaster(); i++)
this->updateConstraintEqMaster(master.nextc->getMaster(i),
offset,ipeq,scale);
}
}

15
src/ASM/SAMpatch.h
View File

@@ -15,6 +15,7 @@
#define _SAM_PATCH_H
#include "SAM.h"
#include "MPC.h"
class ASMbase;
@@ -41,17 +42,25 @@ public:
//! \brief Updates the multi-point constraint array \a TTCC.
//! \param[in] model All spline patches in the model
//! \param[in] prevSol Previous primary solution vector in DOF-order
virtual bool updateConstraintEqs(const std::vector<ASMbase*>& model,
const Vector* prevSol = 0);
bool updateConstraintEqs(const std::vector<ASMbase*>& model,
const Vector* prevSol = 0);
protected:
//! \brief Initializes the nodal arrays \a MINEX, \a MADOF and \a MSC.
bool initNodeDofs(const std::vector<ASMbase*>& model);
//! \brief Initializes the element topology arrays \a MPMNPC and \a MMNPC.
bool initElementConn(const std::vector<ASMbase*>& model);
//! \brief Initializes the multi-point constraint arrays
//! \brief Initializes the multi-point constraint arrays.
//! \a MPMCEQ, \a MMCEQ and \a TTCC.
virtual bool initConstraintEqs(const std::vector<ASMbase*>& model);
private:
//! \brief Recursive helper method used by \a initConstraintEqs.
bool initConstraintEqMaster(const MPC::DOF& master, const MPC::DOF& slave,
Real& offset, int ip, Real scale = Real(1));
//! \brief Recursive helper method used by \a updateConstraintEqs.
void updateConstraintEqMaster(const MPC::DOF& master,
Real& offset, int& ipeq, Real scale = 1.0);
};
#endif

13
src/SIM/SIMbase.C
View File

@@ -945,7 +945,8 @@ bool SIMbase::preprocess (const std::vector<int>& ignored, bool fixDup)
if (!this->initDirichlet()) return false;
// Resolve possibly chaining of the MPC equations
if (!allMPCs.empty()) ASMbase::resolveMPCchains(allMPCs);
if (!allMPCs.empty())
ASMbase::resolveMPCchains(allMPCs,this->hasTimeDependentDirichlet());
// Generate element groups for multi-threading
bool silence = msgLevel < 1 || (msgLevel < 2 && myModel.size() > 1);
@@ -1324,6 +1325,16 @@ size_t SIMbase::getNoSolutions () const
}
bool SIMbase::hasTimeDependentDirichlet () const
{
for (size_t i = 0; i < myModel.size(); i++)
if (myModel[i]->hasTimeDependentDirichlet(myScalars,myVectors))
return true;
return false;
}
bool SIMbase::initDirichlet (double time)
{
if (time == 0.0)

2
src/SIM/SIMbase.h
View File

@@ -224,6 +224,8 @@ public:
//! \brief Initializes time-dependent in-homogeneous Dirichlet coefficients.
//! \param[in] time Current time
bool initDirichlet(double time = 0.0);
//! \brief Checks for time-dependent in-homogeneous Dirichlet conditions.
bool hasTimeDependentDirichlet() const;
//! \brief Updates the time-dependent in-homogeneous Dirichlet coefficients.
//! \param[in] time Current time

12
src/Utility/ExprFunctions.h
View File

@@ -71,6 +71,7 @@ class EvalFunction : public RealFunc
Real* y; //!< Pointer to the Y-coordinate of the function argument
Real* z; //!< Pointer to the Z-coordinate of the function argument
Real* t; //!< Pointer to the time coordinate of the function argument
bool IAmConstant; //!< Indicates whether the time coordinate is given or not
public:
//! \brief The constructor parses the expression string.
@@ -78,6 +79,9 @@ public:
//! \brief The destructor frees the dynamically allocated objects.
virtual ~EvalFunction();
//! \brief Returns whether the function is time-independent or not.
virtual bool isConstant() const { return IAmConstant; }
protected:
//! \brief Evaluates the function expression.
virtual Real evaluate(const Vec3& X) const;
@@ -125,6 +129,14 @@ public:
delete p[i];
}
//! \brief Returns whether the function is time-independent or not.
virtual bool isConstant() const
{
for (size_t i = 0; i < p.size(); i++)
if (!p[i]->isConstant()) return false;
return true;
}
protected:
//! \brief Evaluates the function expressions.
virtual Ret evaluate(const Vec3& X) const;

2
src/Utility/Function.h
View File

@@ -36,6 +36,8 @@ namespace utl
//! \brief Returns whether the function is identically zero or not.
virtual bool isZero() const { return false; }
//! \brief Returns whether the function is time-independent or not.
virtual bool isConstant() const { return true; }
protected:
//! \brief Evaluates the function for the argument \a x.

7
src/Utility/Functions.h
View File

@@ -167,6 +167,8 @@ public:
//! \brief Returns whether the function is identically zero or not.
virtual bool isZero() const { return tfunc->isZero(); }
//! \brief Returns whether the function is time-independent or not.
virtual bool isConstant() const { return false; }
protected:
//! \brief Evaluates the time-varying function.
@@ -193,6 +195,8 @@ public:
//! \brief Returns whether the function is identically zero or not.
virtual bool isZero() const { return sfunc->isZero() || tfunc->isZero(); }
//! \brief Returns whether the function is time-independent or not.
virtual bool isConstant() const { return false; }
protected:
//! \brief Evaluates the space-time function.
@@ -359,6 +363,9 @@ public:
LinearRotZFunc(bool retX, Real a, Real x_0 = Real(0), Real y_0 = Real(0))
: rX(retX), A(a), x0(x_0), y0(y_0) {}
//! \brief Returns whether the function is time-independent or not.
virtual bool isConstant() const { return false; }
protected:
//! \brief Evaluates the rotation function.
virtual Real evaluate(const Vec3& X) const;

13
src/Utility/MPC.C
View File

@@ -57,6 +57,19 @@ bool MPCLess::operator() (const MPC* lhs, const MPC* rhs) const
}
size_t MPC::getNoMaster (bool recursive) const
{
size_t nMaster = master.size();
if (recursive)
for (size_t i = 0; i < master.size(); i++)
if (master[i].nextc)
nMaster += master[i].nextc->getNoMaster(true) - 1;
return nMaster;
}
void MPC::addMaster (const DOF& dof, Real tol)
{
if (dof.coeff >= -tol && dof.coeff <= tol) return; // ignore if zero coeff.

28
src/Utility/MPC.h
View File

@@ -20,7 +20,7 @@
/*!
\brief A class for representing a general multi-point constraint equation.
\brief A class representing a general multi-point constraint equation.
\details A multi-point constraint (MPC) equation is used to introduce
additional coupling between the degrees of freedom (DOF) in a FE grid.
@@ -44,9 +44,8 @@
prescribed (\f$c_0\ne0\f$) dof.
One or more of the master dofs may also be a slave in another constraint
equation (chained constrains). In this case the two equations are combined
to eliminate the master dof that is constrained. This is done while
preprocessing the model by the Grid::resolveMPCchains function.
equation (chained constraints). In this case the two equations are combined
during model preprocessing, to eliminate the master dof that is constrained.
*/
class MPC
@@ -54,19 +53,19 @@ class MPC
public:
/*!
\brief A struct for representing one term (DOF number and associated
coefficient) in a MPC equation.
\brief A struct representing one term in a multi-point constraint equation.
\details Each term consits of the DOF number and an associated coefficient.
*/
struct DOF
{
//! \brief Default constructor.
DOF() : node(0), dof(0), coeff(Real(0)) {}
DOF() : node(0), dof(0), coeff(Real(0)), nextc(0) {}
//! \brief Convenience constructor creating a valid DOF object.
//! \param[in] n Node number (1...NNOD)
//! \param[in] d The local DOF number (1...3)
//! \param[in] c Associated coefficient or constrained value
DOF(int n, int d, Real c = Real(0)) : node(n), dof(d), coeff(c) {}
DOF(int n, int d, Real c = Real(0)) : node(n), dof(d), coeff(c), nextc(0) {}
//! \brief Global stream operator printing a DOF instance.
friend std::ostream& operator<<(std::ostream& s, const DOF& dof)
@@ -83,6 +82,7 @@ public:
int node; //!< Node number identifying this DOF
int dof; //!< Local DOF number within \a node
Real coeff; //!< The constrained value, or master DOF scaling coefficient
MPC* nextc; //!< Points to another MPC in case of chained constraints
};
//! \brief Constructor creating a constraint for a specified slave DOF
@@ -115,6 +115,13 @@ public:
master[pos].coeff = c;
}
//! \brief Updates the next chain element of the \a pos'th master DOF.
void updateMaster(size_t pos, MPC* m)
{
if (pos < master.size())
master[pos].nextc = m;
}
//! \brief Removes the \a pos'th master DOF from the constraint equation.
void removeMaster(size_t pos)
{
@@ -143,12 +150,13 @@ public:
const DOF& getMaster(size_t i) const { return master[i]; }
//! \brief Returns the number of master DOFs.
size_t getNoMaster() const { return master.size(); }
size_t getNoMaster(bool recursive = false) const;
//! \brief Global stream operator printing a constraint equation.
friend std::ostream& operator<<(std::ostream& s, const MPC& mpc);
int iceq; //!< Global constraint equation identifier
int iceq; //!< Global constraint equation identifier
private:
DOF slave; //!< The slave DOF of this constraint equation
std::vector<DOF> master; //!< The master DOFs of this constraint equation