Skip to content

Commit aae1e9f

Browse files
kmokstadkmokstad
committed
Added: Input of general nodal couplings from file
1 parent adf2ab7 commit aae1e9f

3 files changed

Lines changed: 166 additions & 47 deletions

File tree

‎Elasticity/SIMElasticity.C‎

Lines changed: 33 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -366,13 +366,13 @@ void SIMElasticity<Dim>::preprocessA ()
366366

367367
/*!
368368
This method creates the multi-point constraint equations representing the
369-
rigid couplings in the model, if any.
369+
rigid and nodal couplings in the model, if any.
370370
*/
371371

372372
template<class Dim>
373373
bool SIMElasticity<Dim>::preprocessBeforeAsmInit (int& ngnod)
374374
{
375-
return this->addRigidMPCs(this,ngnod);
375+
return this->addRigidMPCs(this,ngnod) && this->addGeneralCouplings(this);
376376
}
377377

378378

@@ -642,26 +642,39 @@ bool SIMElasticity<Dim>::parse (char* keyWord, std::istream& is)
642642
template<class Dim>
643643
bool SIMElasticity<Dim>::parse (const tinyxml2::XMLElement* elem)
644644
{
645-
if (!strcasecmp(elem->Value(),"postprocessing"))
645+
if (strcasecmp(elem->Value(),myContext.c_str()))
646646
{
647-
const tinyxml2::XMLElement* child = elem->FirstChildElement();
648-
for (; child && !plotRgd; child = child->NextSiblingElement())
649-
if (!strcasecmp(child->Value(),"plot_rigid"))
650-
plotRgd = true;
651-
else if (!strcasecmp(child->Value(),"strain"))
652-
Elasticity::wantStrain = true;
653-
}
654-
else if (!strcasecmp(elem->Value(),"localsystem"))
655-
if (ElasticBase* elInt = this->getIntegrand(); elInt)
656-
static_cast<Elasticity*>(elInt)->parseLocalSystem(elem);
647+
if (!this->Dim::parse(elem))
648+
return false;
657649

658-
if (strcasecmp(elem->Value(),myContext.c_str()))
659-
return this->Dim::parse(elem);
650+
else if (!strcasecmp(elem->Value(),"geometry"))
651+
{
652+
for (const tinyxml2::XMLElement* child = elem->FirstChildElement();
653+
child; child = child->NextSiblingElement())
654+
if (!strcasecmp(child->Value(),"coupling"))
655+
this->parseCouplings(child);
656+
}
657+
else if (!strcasecmp(elem->Value(),"postprocessing"))
658+
{
659+
for (const tinyxml2::XMLElement* child = elem->FirstChildElement();
660+
child && !plotRgd; child = child->NextSiblingElement())
661+
if (!strcasecmp(child->Value(),"plot_rigid"))
662+
plotRgd = true;
663+
else if (!strcasecmp(child->Value(),"strain"))
664+
Elasticity::wantStrain = true;
665+
}
666+
else if (!strcasecmp(elem->Value(),"localsystem"))
667+
{
668+
if (ElasticBase* elInt = this->getIntegrand(); elInt)
669+
static_cast<Elasticity*>(elInt)->parseLocalSystem(elem);
670+
}
671+
return true;
672+
}
660673

661674
bool result = true;
662675

663-
const tinyxml2::XMLElement* child = elem->FirstChildElement();
664-
for (; child; child = child->NextSiblingElement())
676+
for (const tinyxml2::XMLElement* child = elem->FirstChildElement();
677+
child; child = child->NextSiblingElement())
665678

666679
if (!strcasecmp(child->Value(),"isotropic") ||
667680
!strcasecmp(child->Value(),"texturematerial"))
@@ -726,6 +739,9 @@ bool SIMElasticity<Dim>::parse (const tinyxml2::XMLElement* elem)
726739
else if (!strcasecmp(child->Value(),"rigid"))
727740
result &= this->parseRigid(child,this);
728741

742+
else if (!strcasecmp(child->Value(),"coupling"))
743+
result &= this->parseCouplings(child);
744+
729745
else if (!strcasecmp(child->Value(),"anasol"))
730746
result &= this->parseAnaSol(child);
731747

‎Elasticity/SIMRigid.C‎

Lines changed: 115 additions & 26 deletions
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,7 @@
77
//!
88
//! \author Knut Morten Okstad / SINTEF
99
//!
10-
//! \brief Rigid Coupling handler for elasticity problems.
10+
//! \brief Rigid and nodal coupling handler for elasticity problems.
1111
//!
1212
//==============================================================================
1313

@@ -19,6 +19,9 @@
1919
#include "IFEM.h"
2020
#include "Vec3Oper.h"
2121
#include "tinyxml2.h"
22+
#include <iostream>
23+
#include <fstream>
24+
#include <sstream>
2225

2326

2427
bool SIMRigid::parseRigid (const tinyxml2::XMLElement* elem, SIMinput* mySim)
@@ -63,13 +66,69 @@ bool SIMRigid::parseRigid (const tinyxml2::XMLElement* elem, SIMinput* mySim)
6366
}
6467

6568

69+
bool SIMRigid::parseCouplings (const tinyxml2::XMLElement* elem)
70+
{
71+
// Lambda function parsing the nodal couplings for a patch.
72+
auto&& parseCpl = [](std::istream& is, Couplin& cpl, char delim = '\n')
73+
{
74+
std::string cline;
75+
int slave, master;
76+
double weight;
77+
while (std::getline(is,cline,delim))
78+
{
79+
std::istringstream ss(cline);
80+
ss >> slave >> master >> weight;
81+
while (ss)
82+
{
83+
cpl[slave].emplace_back(master,weight);
84+
ss >> master >> weight;
85+
}
86+
#if INT_DEBUG > 1
87+
std::cout <<"Slave ("<< slave <<") coupled to master(s):";
88+
for (const Master& mst : cpl[slave])
89+
std::cout <<" ("<< mst.first <<")*"<< mst.second;
90+
std::cout << std::endl;
91+
#endif
92+
}
93+
};
94+
95+
IFEM::cout <<" Parsing <"<< elem->Value() <<">"<< std::endl;
96+
97+
int patch = 1;
98+
utl::getAttribute(elem,"patch",patch);
99+
if (std::string fName; utl::getAttribute(elem,"file",fName))
100+
{
101+
std::ifstream fs(fName);
102+
parseCpl(fs,myCouplings[patch]);
103+
}
104+
else if (elem->FirstChild())
105+
{
106+
std::istringstream ss(elem->FirstChild()->Value());
107+
parseCpl(ss,myCouplings[patch],'\\');
108+
}
109+
else
110+
return false;
111+
112+
IFEM::cout <<"\tParsed "<< myCouplings[patch].size()
113+
<<" nodal couplings for Patch "<< patch;
114+
if (double xtol = 0.0; utl::getAttribute(elem,"xtol",xtol) && xtol > 0.0)
115+
{
116+
myGeomTols[patch] = xtol;
117+
IFEM::cout <<" (tolerance = "<< xtol <<")";
118+
}
119+
120+
IFEM::cout << std::endl;
121+
return true;
122+
}
123+
124+
66125
bool SIMRigid::addRigidMPCs (SIMinput* mySim, int& ngnod) const
67126
{
68-
PropertyVec::const_iterator pit;
69-
for (pit = mySim->begin_prop(); pit != mySim->end_prop(); ++pit)
127+
for (PropertyVec::const_iterator pit = mySim->begin_prop();
128+
pit != mySim->end_prop(); ++pit)
70129
if (pit->pcode == Property::RIGID)
71130
{
72-
std::map<int,TopItem>::const_iterator mit = myMasters.find(pit->pindx);
131+
PntMap::const_iterator mit = myMasters.find(pit->pindx);
73132
if (mit == myMasters.end())
74133
{
75134
std::cerr <<" *** SIMRigid::addRigidMPCs: Local error"
@@ -86,17 +145,51 @@ bool SIMRigid::addRigidMPCs (SIMinput* mySim, int& ngnod) const
86145

87146
// Find master point and the affected patch,
88147
// and create a rigid coupling between them
89-
SIMinput::IdxVec3* mst = mySim->getDiscretePoint(mit->second.item);
90-
ASMbase* pch = mySim->getPatch(pit->patch);
91-
if (mst && pch)
92-
if (pch->addRigidCpl(pit->lindx,pit->ldim,pit->basis,
93-
mst->first,mst->second)) ++ngnod;
148+
if (SIMinput::IdxVec3* mp = mySim->getDiscretePoint(mit->second.item); mp)
149+
if (ASMbase* pch = mySim->getPatch(pit->patch); pch)
150+
if (pch->addRigidCpl(pit->lindx,pit->ldim,pit->basis,
151+
mp->first,mp->second)) ++ngnod;
94152
}
95153

96154
return true;
97155
}
98156

99157

158+
bool SIMRigid::addGeneralCouplings (SIMinput* mySim) const
159+
{
160+
bool status = true;
161+
for (const CplMap::value_type& cpl : myCouplings)
162+
if (ASMbase* pch = mySim->getPatch(cpl.first); pch)
163+
{
164+
std::vector<Ipair> connectedNodes;
165+
for (const Couplin::value_type& mpc : cpl.second)
166+
if (mpc.second.size() == 1 && mpc.second.front().second == 1.0)
167+
// Direct coupling of matching nodes
168+
connectedNodes.emplace_back(mpc.first,mpc.second.front().first);
169+
else
170+
{
171+
// Create linear couplings
172+
IntVec masters;
173+
RealArray weights;
174+
for (const Master& mst : mpc.second)
175+
{
176+
masters.push_back(mst.first);
177+
weights.push_back(mst.second);
178+
}
179+
pch->addNodalCouplings(mpc.first,masters,weights);
180+
}
181+
182+
if (TolMap::const_iterator it = myGeomTols.find(cpl.first);
183+
it != myGeomTols.end())
184+
status &= pch->selfInterconnect(connectedNodes,it->second);
185+
else
186+
status &= pch->selfInterconnect(connectedNodes);
187+
}
188+
189+
return status;
190+
}
191+
192+
100193
ElementBlock* SIMRigid::rigidGeometry (SIMinput* mySim) const
101194
{
102195
if (myMasters.empty()) return nullptr; // no rigid couplings
@@ -105,31 +198,27 @@ ElementBlock* SIMRigid::rigidGeometry (SIMinput* mySim) const
105198
std::map<int,size_t> imap;
106199
ElementBlock* rgd = new ElementBlock(2);
107200
rgd->unStructResize(0,myMasters.size());
108-
for (const std::pair<const int,TopItem>& master : myMasters)
201+
for (const PntMap::value_type& master : myMasters)
109202
{
110-
SIMinput::IdxVec3* mst = mySim->getDiscretePoint(master.second.item);
111-
if (mst) rgd->setCoor(inod,mst->second);
203+
if (SIMinput::IdxVec3* mp = mySim->getDiscretePoint(master.second.item); mp)
204+
rgd->setCoor(inod,mp->second);
112205
imap[master.first] = inod++;
113206
}
114207

115208
int slvThick = mySim->opt.discretization == ASM::SplineC1 ? 2 : 1;
116209

117-
PropertyVec::const_iterator pit;
118-
std::map<int,TopItem>::const_iterator mit;
119-
for (pit = mySim->begin_prop(); pit != mySim->end_prop(); ++pit)
210+
for (PropertyVec::const_iterator pit = mySim->begin_prop();
211+
pit != mySim->end_prop(); ++pit)
120212
if (pit->pcode == Property::RIGID)
121-
if ((mit = myMasters.find(pit->pindx)) != myMasters.end())
122-
if (mit->second.patch == 0)
213+
if (PntMap::const_iterator mit = myMasters.find(pit->pindx);
214+
mit != myMasters.end() && mit->second.patch == 0)
215+
if (ASMbase* pch = mySim->getPatch(pit->patch); pch &&
216+
mySim->getDiscretePoint(mit->second.item))
123217
{
124-
SIMinput::IdxVec3* mst = mySim->getDiscretePoint(mit->second.item);
125-
ASMbase* pch = mySim->getPatch(pit->patch);
126-
if (mst && pch)
127-
{
128-
IntVec nodes;
129-
pch->getBoundaryNodes(pit->lindx,nodes,pit->basis,slvThick,0,true);
130-
for (int node : nodes)
131-
rgd->addLine(imap[mit->first],pch->getCoord(node));
132-
}
218+
IntVec nodes;
219+
pch->getBoundaryNodes(pit->lindx,nodes,pit->basis,slvThick,0,true);
220+
for (int node : nodes)
221+
rgd->addLine(imap[mit->first],pch->getCoord(node));
133222
}
134223

135224
return rgd;

‎Elasticity/SIMRigid.h‎

Lines changed: 18 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -7,22 +7,23 @@
77
//!
88
//! \author Knut Morten Okstad / SINTEF
99
//!
10-
//! \brief Rigid Coupling handler for elasticity problems.
10+
//! \brief Rigid and nodal coupling handler for elasticity problems.
1111
//!
1212
//=============================================================================
1313

1414
#ifndef _SIM_RIGID_H
1515
#define _SIM_RIGID_H
1616

1717
#include "TopologySet.h"
18+
#include <vector>
1819

1920
class SIMinput;
2021
class ElementBlock;
21-
namespace tinyxml2 { class XMLElement; }
22+
namespace tinyxml2 { class XMLElement; }
2223

2324

2425
/*!
25-
\brief Rigid coupling handler for elasticity problems.
26+
\brief Rigid and nodal coupling handler for elasticity problems.
2627
*/
2728

2829
class SIMRigid
@@ -35,15 +36,28 @@ class SIMRigid
3536

3637
//! \brief Parses a rigid coupling definition from an XML element.
3738
bool parseRigid(const tinyxml2::XMLElement* elem, SIMinput* mySim);
39+
//! \brief Parses general nodal couplings from an XML element.
40+
bool parseCouplings(const tinyxml2::XMLElement* elem);
3841

3942
//! \brief Creates multi-point constraint equations for the rigid couplings.
4043
bool addRigidMPCs(SIMinput* mySim, int& ngnod) const;
44+
//! \brief Creates multi-point constraint equations for the general couplings.
45+
bool addGeneralCouplings(SIMinput* mySim) const;
4146

4247
//! \brief Creates an element block visualizing the rigid couplings.
4348
ElementBlock* rigidGeometry(SIMinput* mySim) const;
4449

4550
private:
46-
std::map<int,TopItem> myMasters; //!< Discrete master points
51+
using PntMap = std::map<int,TopItem>; //!< Discrete master point definition
52+
using Master = std::pair<int,double>; //!< Independent node with weight
53+
using Masters = std::vector<Master>; //!< Independent node list
54+
using Couplin = std::map<int,Masters>; //!< General nodal coupling
55+
using CplMap = std::map<int,Couplin>; //!< Couplings for patches
56+
using TolMap = std::map<int,double>; //!< Geometry tolerance for patches
57+
58+
PntMap myMasters; //!< Discrete master point container
59+
CplMap myCouplings; //!< Nodal coupling container
60+
TolMap myGeomTols; //!< Geometry tolerance container
4761
};
4862

4963
#endif

0 commit comments

Comments
 (0)