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

Fixed: Processing of chained MPCs for time-dependent problems.

Browse Source 
Changed: Replaced private methods by recursive lambda functions.
Added: Use of std::setw to improve formatting of debug output.
This commit is contained in:
Knut Morten Okstad 2020-03-03 08:40:57 +01:00
parent ec5eb6ad25
commit 40e481a6fa
5 changed files with 213 additions and 189 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

165
src/ASM/ASMbase.C
View File

@ -22,6 +22,8 @@
#include "Function.h"
#include "Utilities.h"
#include <algorithm>
#include <functional>
#include <iomanip>
bool ASMbase::fixHomogeneousDirichlet = true;
@ -379,7 +381,8 @@ void ASMbase::printNodes (std::ostream& os) const
os <<"\n\nNodal coordinates for Patch "<< idx+1;
for (size_t inod = 1; inod <= X.cols(); inod++)
{
os <<'\n'<< inod <<' '<< MLGN[inod-1] <<':';
os <<'\n'<< std::setw(4) << inod
<< ' '<< std::setw(4) << MLGN[inod-1] <<':';
for (size_t i = 1; i <= X.rows(); i++)
os <<' '<< X(i,inod);
}
@ -394,13 +397,14 @@ void ASMbase::printElements (std::ostream& os) const
os <<"\n\nElement connectivities for Patch "<< idx+1;
for (size_t iel = 0; iel < MLGE.size(); iel++)
{
os <<'\n'<< iel+1 <<' '<< MLGE[iel] <<':';
os <<'\n'<< std::setw(4) << iel+1
<< ' '<< std::setw(4) << MLGE[iel] <<':';
if (iel < MNPC.size())
for (int node : MNPC[iel]) os <<" "<< node;
}
for (size_t ielx = MLGE.size(); ielx < MNPC.size(); ielx++)
{
os <<'\n'<< ielx+1 <<" ---:";
os <<'\n'<< std::setw(4) << ielx+1 <<" ----:";
for (int node : MNPC[ielx]) os <<" "<< node;
}
os << std::endl;
@ -644,7 +648,7 @@ int ASMbase::prescribe (size_t inod, int dirs, int code)
for (int dof : utl::getDigits(dirs))
if (dof <= nf)
{
MPC* mpc = new MPC(node,dof);
MPC* mpc = new MPC(node,dof,1.0);
if (!this->addMPC(mpc,code))
ignoredDirs = 10*ignoredDirs + dof;
}
@ -805,6 +809,12 @@ void ASMbase::mergeAndGetAllMPCs (const ASMVec& model, MPCSet& allMPCs)
/*!
Recursive resolving of (possibly multi-level) chaining in the 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.
Since an MPC-equation may couple nodes belonging to different patches,
this method must have access to all patches in the model.
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.
@ -815,24 +825,104 @@ void ASMbase::mergeAndGetAllMPCs (const ASMVec& model, MPCSet& allMPCs)
\sa SAMpatch::updateConstraintEqs.
*/
void ASMbase::resolveMPCchains (const MPCSet& allMPCs, bool setPtrOnly)
void ASMbase::resolveMPCchains (const MPCSet& allMPCs,
const ASMVec& model, bool setPtrOnly)
{
#if SP_DEBUG > 1
int count = 0;
std::cout <<"\nResolving MPC chains"<< std::endl;
for (MPC* mpc : allMPCs) std::cout << *mpc;
for (MPC* mpc : allMPCs) std::cout << std::setw(4) << ++count <<": "<< *mpc;
if (setPtrOnly) std::cout <<"\nResolved MPCs"<< std::endl;
count = 0;
#endif
// Recursive lambda function for resolving the chained MPC-equations.
std::function<bool(MPC*)> resolve = [&allMPCs,&resolve](MPC* mpc)
{
if (!mpc) return false;
bool resolved = false;
for (size_t i = 0; i < mpc->getNoMaster();)
{
MPC master(mpc->getMaster(i).node,mpc->getMaster(i).dof);
MPCIter cit = allMPCs.find(&master);
if (cit != allMPCs.end())
{
// We have a master dof which is a slave in another constraint equation.
// Invoke resolve() recursively to ensure that all master dofs of that
// equation are not slaves themselves.
resolve(*cit);
// Remove current master specification
double coeff = mpc->getMaster(i).coeff;
mpc->removeMaster(i);
// Add constant offset from the other equation
mpc->addOffset(coeff*(*cit)->getSlave().coeff);
// Add masters from the other equations
for (size_t j = 0; j < (*cit)->getNoMaster(); j++)
mpc->addMaster((*cit)->getMaster(j).node,
(*cit)->getMaster(j).dof,
(*cit)->getMaster(j).coeff*coeff);
resolved = true;
}
else
i++;
}
#if SP_DEBUG > 1
if (resolved) std::cout <<"Resolved constraint: "<< *mpc;
#endif
return resolved;
};
// Recursive lambda function for linking the chained MPC-equations.
std::function<bool(MPC*)> resolveLnk = [&allMPCs,&model,&resolveLnk](MPC* mpc)
{
if (!mpc) return false;
for (size_t i = 0; i < mpc->getNoMaster();)
{
int node = mpc->getMaster(i).node;
int ldof = mpc->getMaster(i).dof;
for (ASMbase* pch : model)
if (pch->isFixed(node,ldof))
{
ldof = -1;
break;
}
if (ldof > 0)
{
MPC master(node,ldof);
MPCIter chain = allMPCs.find(&master);
if (chain == allMPCs.end())
i++; // Free master DOF
else if (!resolveLnk(*chain))
mpc->removeMaster(i); // This master DOF is fixed, remove from link
else
mpc->updateMaster(i++,*chain); // Set pointer to next MPC in chain
}
else
mpc->removeMaster(i);
}
return mpc->getNoMaster() > 0 || mpc->getSlave().coeff != 0.0;
};
int nresolved = 0;
for (MPC* mpc : allMPCs)
if (setPtrOnly)
for (size_t i = 0; i < mpc->getNoMaster(); i++)
{
MPC master(mpc->getMaster(i).node,mpc->getMaster(i).dof);
MPCIter chain = allMPCs.find(&master);
if (chain != allMPCs.end())
mpc->updateMaster(i,*chain); // Set pointer to next MPC in chain
}
else if (ASMbase::resolveMPCchain(allMPCs,mpc))
{
resolveLnk(mpc);
#if SP_DEBUG > 1
count++;
if (mpc->isChained() || mpc->getNoMaster() == 0)
std::cout << std::setw(4) << count <<": "<< *mpc;
#endif
}
else if (resolve(mpc))
nresolved++;
if (nresolved > 0)
@ -840,53 +930,6 @@ void ASMbase::resolveMPCchains (const MPCSet& allMPCs, bool setPtrOnly)
}
/*!
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)
{
if (!mpc) return false;
bool resolved = false;
for (size_t i = 0; i < mpc->getNoMaster();)
{
MPC master(mpc->getMaster(i).node,mpc->getMaster(i).dof);
MPCIter cit = allMPCs.find(&master);
if (cit != allMPCs.end())
{
// We have a master dof which is a slave in another constraint equation.
// Invoke resolveMPCchain recursively to ensure that all master dofs
// of that equation are not slaves themselves.
ASMbase::resolveMPCchain(allMPCs,*cit);
// Remove current master specification
double coeff = mpc->getMaster(i).coeff;
mpc->removeMaster(i);
// Add constant offset from the other equation
mpc->addOffset(coeff*(*cit)->getSlave().coeff);
// Add masters from the other equations
for (size_t j = 0; j < (*cit)->getNoMaster(); j++)
mpc->addMaster((*cit)->getMaster(j).node,
(*cit)->getMaster(j).dof,
(*cit)->getMaster(j).coeff*coeff);
resolved = true;
}
else
i++;
}
#if SP_DEBUG > 1
if (resolved) std::cout <<"Resolved constraint: "<< *mpc;
#endif
return resolved;
}
void ASMbase::setNodeNumbers (const IntVec& nodes)
{
for (size_t i = 0; i < nodes.size() && i < myMLGN.size(); i++)

15
src/ASM/ASMbase.h
View File

@ -388,13 +388,10 @@ public:
//! \brief Resolves (possibly multi-level) chaining in MPC equations.
//! \param[in] allMPCs All multi-point constraint equations in the model
//! \param[in] model All spline patches 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, bool setPtrOnly = false);
static void resolveMPCchains(const MPCSet& allMPCs,
const ASMVec& model, bool setPtrOnly = false);
//! \brief Initializes the multi-point constraint coefficients.
virtual bool initConstraints() { return true; }
@ -814,12 +811,6 @@ protected:
//! is changed into the number of the other node.
static bool collapseNodes(ASMbase& pch1, int node1, ASMbase& pch2, int node2);
private:
//! \brief Recursive method used by resolveMPCchains().
//! \param[in] allMPCs All multi-point constraint equations in the model
//! \param mpc Pointer to the multi-point constraint equation to resolve
static bool resolveMPCchain(const MPCSet& allMPCs, MPC* mpc);
public:
static bool fixHomogeneousDirichlet; //!< If \e true, pre-eliminate fixed DOFs

134
src/ASM/SAMpatch.C
View File

@ -13,7 +13,9 @@
#include "SAMpatch.h"
#include "ASMbase.h"
#include "MPC.h"
#include "IFEM.h"
#include <functional>
bool SAMpatch::init (const std::vector<ASMbase*>& patches, int numNod,
@ -258,6 +260,31 @@ bool SAMpatch::initConstraintEqs ()
int ip = mpmceq[0] = 1;
if (nceq < 1) return true;
// Recursive lambda function for resolving master DOF dependencies.
std::function<void(const MPC::DOF&,const MPC::DOF&,Real&,int&,Real)> resolve;
resolve = [this,&resolve](const MPC::DOF& master, const MPC::DOF& slave,
Real& offset, int& ipeq, Real scale)
{
int idof = madof[master.node-1] + master.dof - 1;
if (master.nextc)
{
// This master DOF is constrained (unresolved chaining)
scale *= master.coeff;
offset += master.nextc->getSlave().coeff*scale;
for (size_t i = 0; i < master.nextc->getNoMaster(); i++)
resolve(master.nextc->getMaster(i),
master.nextc->getSlave(),
offset,ipeq,scale);
}
else if (msc[idof-1] > 0)
{
// This master DOF is free
int ipmst = (ipeq++) - 1;
mmceq[ipmst] = idof;
ttcc[ipmst] = master.coeff*scale;
}
};
mmceq = new int[nmmceq];
ttcc = new Real[nmmceq];
for (const ASMbase* pch : model)
@ -267,20 +294,12 @@ bool SAMpatch::initConstraintEqs ()
// Slave dof ...
int idof = madof[(*cit)->getSlave().node-1] + (*cit)->getSlave().dof - 1;
if (msc[idof-1] == 0)
if (msc[idof-1] <= 0)
{
std::cerr <<"SAM: Ignoring constraint equation for dof "
<< idof <<" ("<< (*cit)->getSlave()
<<").\n This dof is already marked as FIXED."<< std::endl;
ip--;
nceq--;
continue;
}
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 "
<< (msc[idof-1] == 0 ? "FIXED." : "SLAVE.") << std::endl;
ip--;
nceq--;
continue;
@ -293,19 +312,24 @@ bool SAMpatch::initConstraintEqs ()
ttcc[ipslv] = (*cit)->getSlave().coeff;
msc[idof-1] = -ip;
#if SP_DEBUG > 1
if ((*cit)->getNoMaster() > 0)
std::cout <<"Resolving "<< **cit;
#endif
// Master dofs ...
for (size_t i = 0; i < (*cit)->getNoMaster(); i++)
if (!this->initConstraintEqMaster((*cit)->getMaster(i),
(*cit)->getSlave(),
ttcc[ipslv],ip))
{
// Something is wrong, ignore this constraint equation
(*cit)->iceq = -1;
mpmceq[ip] = mpmceq[ip-1];
ip--;
nceq--;
break;
}
resolve((*cit)->getMaster(i),(*cit)->getSlave(),
ttcc[ipslv],mpmceq[ip],Real(1));
#if SP_DEBUG > 1
if ((*cit)->getNoMaster() > 0)
{
std::cout <<"Resolved CEQ: ";
this->printCEQ(std::cout,ip-1);
std::cout << std::endl;
}
#endif
}
#if SP_DEBUG > 1
@ -323,43 +347,26 @@ bool SAMpatch::initConstraintEqs ()
}
bool SAMpatch::initConstraintEqMaster (const MPC::DOF& master,
const MPC::DOF& slave,
Real& offset, int ip, Real scale)
bool SAMpatch::updateConstraintEqs (const Vector* prevSol)
{
int idof = madof[master.node-1] + master.dof - 1;
if (msc[idof-1] > 0)
if (nceq < 1) return true; // No constraints in this model
// Recursive lambda function for updating master DOFs.
std::function<void(const MPC::DOF&,Real&,int&,Real)> update;
update = [this,&update](const MPC::DOF& master,
Real& offset, int& ipeq, Real scale)
{
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
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++)
if (!this->initConstraintEqMaster(master.nextc->getMaster(i),
master.nextc->getSlave(),
offset,ip,scale))
return false;
update(master.nextc->getMaster(i),offset,ipeq,scale);
}
return true;
}
bool SAMpatch::updateConstraintEqs (const Vector* prevSol)
{
if (nceq < 1) return true; // No constraints in this model
};
MPCIter cit;
for (const ASMbase* pch : model)
@ -382,7 +389,7 @@ bool SAMpatch::updateConstraintEqs (const Vector* prevSol)
// Master dofs ...
for (size_t i = 0; prevSol && i < (*cit)->getNoMaster(); i++)
this->updateConstraintEqMaster((*cit)->getMaster(i),c0,ipeq);
update((*cit)->getMaster(i),c0,ipeq,Real(1));
// Update the constant term of the constraint equation
if (!prevSol || c0 == Real(0)) // Note: c0=0 means no constant offset
@ -400,20 +407,3 @@ bool SAMpatch::updateConstraintEqs (const Vector* prevSol)
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);
}
}

8
src/ASM/SAMpatch.h
View File

@ -15,7 +15,6 @@
#define _SAM_PATCH_H
#include "SAM.h"
#include "MPC.h"
class ASMbase;
@ -62,13 +61,6 @@ protected:
bool initConstraintEqs();
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);
std::vector<ASMbase*> model; //!< The spline patches of the model
};

80
src/SIM/SIMbase.C
View File

@ -203,46 +203,49 @@ bool SIMbase::preprocess (const IntVec& ignored, bool fixDup)
this->preprocessA();
// Create the classical FE data structures
if (!this->createFEMmodel('Y')) return false;
PatchVec::const_iterator mit;
ASMbase* pch;
size_t patch;
if (!this->createFEMmodel('Y'))
return false;
// Erase all patches that should be ignored in the analysis
for (int idx : ignored)
if ((pch = this->getPatch(idx)))
pch->clear();
{
ASMbase* pch = this->getPatch(idx);
if (pch) pch->clear();
}
// If material properties are specified for at least one patch, assign the
// property code 999999 to all patches with no material property code yet
PatchVec pchWthMat;
PatchVec patches;
for (const Property& p : myProps)
if (p.pcode == Property::MATERIAL && (pch = this->getPatch(p.patch)))
if (!pch->empty()) pchWthMat.push_back(pch);
if (p.pcode == Property::MATERIAL)
{
ASMbase* pch = this->getPatch(p.patch);
if (pch && !pch->empty())
patches.push_back(pch);
}
if (!pchWthMat.empty())
for (mit = myModel.begin(), patch = 1; mit != myModel.end(); ++mit, patch++)
if (std::find(pchWthMat.begin(),pchWthMat.end(),*mit) == pchWthMat.end())
myProps.push_back(Property(Property::MATERIAL,999999,patch,
(*mit)->getNoParamDim()));
if (!patches.empty())
for (size_t i = 0; i < myModel.size(); i++)
if (std::find(patches.begin(),patches.end(),myModel[i]) == patches.end())
myProps.push_back(Property(Property::MATERIAL,999999,i+1,
myModel[i]->getNoParamDim()));
if (fixDup)
{
// Check for duplicated nodes (missing topology)
int nDupl = 0;
std::map<Vec3,int> globalNodes;
for (mit = myModel.begin(), patch = 1; mit != myModel.end(); ++mit, patch++)
if (!(*mit)->empty())
for (ASMbase* pch : myModel)
if (!pch->empty())
{
IFEM::cout <<" * Checking Patch "<< patch << std::endl;
for (size_t node = 1; node <= (*mit)->getNoNodes(); node++)
IFEM::cout <<" * Checking Patch "<< pch->idx+1 << std::endl;
for (size_t node = 1; node <= pch->getNoNodes(); node++)
{
Vec3 X((*mit)->getCoord(node));
Vec3 X(pch->getCoord(node));
std::map<Vec3,int>::const_iterator xit = globalNodes.find(X);
if (xit == globalNodes.end())
globalNodes.insert(std::make_pair(X,(*mit)->getNodeID(node)));
else if ((*mit)->mergeNodes(node,xit->second))
globalNodes.insert(std::make_pair(X,pch->getNodeID(node)));
else if (pch->mergeNodes(node,xit->second))
nDupl++;
}
}
@ -267,8 +270,8 @@ bool SIMbase::preprocess (const IntVec& ignored, bool fixDup)
renum = ASMbase::renumberNodes(myModel,myGlb2Loc);
ngnod = g2l->size();
}
else for (mit = myModel.begin(); mit != myModel.end(); ++mit)
renum += (*mit)->renumberNodes(myGlb2Loc,ngnod);
else for (ASMbase* pch : myModel)
renum += pch->renumberNodes(myGlb2Loc,ngnod);
if (renum > 0)
IFEM::cout <<"\nRenumbered "<< renum <<" nodes."<< std::endl;
@ -340,19 +343,29 @@ bool SIMbase::preprocess (const IntVec& ignored, bool fixDup)
MPCSet allMPCs;
ASMbase::mergeAndGetAllMPCs(myModel,allMPCs);
// Set initial values for the inhomogeneous dirichlet conditions, if any,
// and compute coupling coefficients for the C1-continuity constraints
if (!this->initDirichlet()) return false;
// Compute coupling coefficients for C1-continuity constraints, if any
for (ASMbase* pch : myModel)
if (!pch->initConstraints())
return false;
// Set the inhomogeneous dirichlet conditions, in stationary case
bool timeDependent = this->hasTimeDependentDirichlet();
if (!timeDependent && !this->initDirichlet())
return false;
// Resolve possibly chaining of the MPC equations
if (!allMPCs.empty())
ASMbase::resolveMPCchains(allMPCs,this->hasTimeDependentDirichlet());
ASMbase::resolveMPCchains(allMPCs,myModel,timeDependent);
// Set initial values for the inhomogeneous dirichlet conditions, if any
if (timeDependent && !this->initDirichlet())
return false;
// Generate element groups for multi-threading
bool silence = msgLevel < 1 || (msgLevel < 3 && nGlPatches > 1);
for (mit = myModel.begin(); mit != myModel.end() && myProblem; ++mit)
if (!(*mit)->empty())
(*mit)->generateThreadGroups(*myProblem,silence,lagMTOK);
for (ASMbase* pch : myModel)
if (!pch->empty())
pch->generateThreadGroups(*myProblem,silence,lagMTOK);
for (const Property& p : myProps)
if (p.pcode == Property::NEUMANN ||
@ -685,11 +698,6 @@ bool SIMbase::hasTimeDependentDirichlet () const
bool SIMbase::initDirichlet (double time)
{
if (time == 0.0)
for (ASMbase* pch : myModel)
if (!pch->initConstraints())
return false;
Vector dummy;
return this->updateDirichlet(time,&dummy);
}