Nektar++
IterativeElasticSystem.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File IterativeElasticSystem.cpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: IterativeElasticSystem solve routines
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <boost/format.hpp>
37 
40 #include <StdRegions/StdQuadExp.h>
41 #include <StdRegions/StdTriExp.h>
42 #include <LocalRegions/MatrixKey.h>
46 
48 
49 namespace Nektar
50 {
51 
53  RegisterCreatorFunction("IterativeElasticSystem",
55 
58  : LinearElasticSystem(pSession)
59 {
60 }
61 
63 {
65 
66  const int nVel = m_fields[0]->GetCoordim(0);
67 
68  // Read in number of steps to take.
69  m_session->LoadParameter("NumSteps", m_numSteps, 0);
70  ASSERTL0(m_numSteps > 0, "You must specify at least one step");
71 
72  // Read in whether to repeatedly apply boundary conditions (for e.g.
73  // rotation purposes).
74  string bcType;
75  m_session->LoadSolverInfo("BCType", bcType, "Normal");
76  m_repeatBCs = bcType != "Normal";
77 
78  if (!m_repeatBCs)
79  {
80  // Loop over BCs, identify which ones we need to deform.
82  &bndCond = m_fields[0]->GetBndConditions();
83 
84  for (int i = 0; i < bndCond.num_elements(); ++i)
85  {
86  if (boost::iequals(bndCond[i]->GetUserDefined(), "Wall"))
87  {
88  m_toDeform.push_back(i);
89  }
90  }
91 
92  int numDeform = m_toDeform.size();
93  m_comm->AllReduce(numDeform, LibUtilities::ReduceMax);
94  ASSERTL0(numDeform > 0, "You must specify at least one WALL tag on"
95  "a boundary condition");
96 
98  m_toDeform.size());
99 
100  for (int i = 0; i < m_toDeform.size(); ++i)
101  {
103  for (int j = 0; j < nVel; ++j)
104  {
105  const int id = m_toDeform[i];
106  MultiRegions::ExpListSharedPtr bndCondExp =
107  m_fields[j]->GetBndCondExpansions()[id];
108  int nCoeffs = bndCondExp->GetNcoeffs();
109 
110  m_savedBCs[i][j] = Array<OneD, NekDouble>(nCoeffs);
111  Vmath::Smul(nCoeffs, 1.0/m_numSteps,
112  bndCondExp->GetCoeffs(), 1,
113  bndCondExp->UpdateCoeffs(), 1);
114  Vmath::Vcopy(nCoeffs, bndCondExp->GetCoeffs(), 1,
115  m_savedBCs[i][j], 1);
116  }
117  }
118  }
119 }
120 
122 {
124 }
125 
127 {
128  int i, j, k;
129 
130  // Write initial geometry for consistency/script purposes
131  WriteGeometry(0);
132 
133  // Now loop over desired number of steps
134  for (i = 1; i <= m_numSteps; ++i)
135  {
136  int invalidElmtId = -1;
137 
138  // Perform solve for this iteration and update geometry accordingly.
141  WriteGeometry(i);
142 
143  // Check for invalid elements.
144  for (j = 0; j < m_fields[0]->GetExpSize(); ++j)
145  {
147  m_fields[0]->GetExp(j)->GetGeom()->GetGeomFactors();
148 
149  if (!geomFac->IsValid())
150  {
151  invalidElmtId =
152  m_fields[0]->GetExp(j)->GetGeom()->GetGlobalID();
153  break;
154  }
155  }
156 
157  m_session->GetComm()->AllReduce(invalidElmtId, LibUtilities::ReduceMax);
158 
159  // If we found an invalid element, exit loop without writing output.
160  if (invalidElmtId >= 0)
161  {
162  if (m_session->GetComm()->GetRank() == 0)
163  {
164  cout << "- Detected negative Jacobian in element "
165  << invalidElmtId << "; terminating at"
166  " step: "<<i<< endl;
167  }
168 
169  break;
170  }
171 
172  if (m_session->GetComm()->GetRank() == 0)
173  {
174  cout << "Step: " << i << endl;
175  }
176 
177  // Update boundary conditions
178  if (m_repeatBCs)
179  {
180  for (j = 0; j < m_fields.num_elements(); ++j)
181  {
182  string varName = m_session->GetVariable(j);
183  m_fields[j]->EvaluateBoundaryConditions(m_time, varName);
184  }
185  }
186  else
187  {
188  for (j = 0; j < m_fields.num_elements(); ++j)
189  {
191  &bndCondExp = m_fields[j]->GetBndCondExpansions();
192 
193  for (k = 0; k < m_toDeform.size(); ++k)
194  {
195  const int id = m_toDeform[k];
196  const int nCoeffs = bndCondExp[id]->GetNcoeffs();
197  Vmath::Vcopy(nCoeffs,
198  m_savedBCs[k][j], 1,
199  bndCondExp[id]->UpdateCoeffs(), 1);
200  }
201  }
202  }
203  }
204 }
205 
206 /**
207  * @brief Write out a file in serial or directory in parallel containing new
208  * mesh geometry.
209  */
211 {
212  fs::path filename;
213  stringstream s;
214  s << m_session->GetSessionName() << "-" << i;
215 
216  if (m_session->GetComm()->GetSize() > 1)
217  {
218  s << "_xml";
219 
220  if(!fs::is_directory(s.str()))
221  {
222  fs::create_directory(s.str());
223  }
224 
225  boost::format pad("P%1$07d.xml");
226  pad % m_session->GetComm()->GetRank();
227  filename = fs::path(s.str()) / fs::path(pad.str());
228  }
229  else
230  {
231  s << ".xml";
232  filename = fs::path(s.str());
233  }
234 
235  string fname = LibUtilities::PortablePath(filename);
236  m_fields[0]->GetGraph()->WriteGeometry(fname);
237 }
238 
239 }
static EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession)
Creates an instance of this class.
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
virtual void v_InitObject()
Set up the linear elasticity system.
virtual void v_GenerateSummary(SolverUtils::SummaryList &s)
Generate summary at runtime.
NekDouble m_time
Current time of simulation.
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: Misc.h:47
static std::string className
Name of class.
int m_numSteps
Number of steps to split deformation across.
virtual void v_DoSolve()
Solve elliptic linear elastic system.
vector< int > m_toDeform
Vector of boundary regions to deform.
bool m_repeatBCs
Flag determining whether to repeat boundary conditions.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:50
LibUtilities::CommSharedPtr m_comm
Communicator.
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_savedBCs
Storage for boundary conditions.
IterativeElasticSystem(const LibUtilities::SessionReaderSharedPtr &pSession)
Base class for linear elastic system.
virtual void v_InitObject()
Set up the linear elasticity system.
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:199
virtual void v_GenerateSummary(SolverUtils::SummaryList &s)
Generate summary at runtime.
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
Definition: ExpList.h:1340
void UpdateGeometry(SpatialDomains::MeshGraphSharedPtr graph, Array< OneD, MultiRegions::ExpListSharedPtr > &fields)
Update geometry according to displacement that is in current fields.
Definition: Deform.cpp:57
std::string PortablePath(const boost::filesystem::path &path)
create portable path on different platforms for boost::filesystem path
Definition: FileSystem.cpp:41
EquationSystemFactory & GetEquationSystemFactory()
boost::shared_ptr< GeomFactors > GeomFactorsSharedPtr
Pointer to a GeomFactors object.
Definition: GeomFactors.h:62
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
SpatialDomains::MeshGraphSharedPtr m_graph
Pointer to graph defining mesh.
virtual void v_DoSolve()
Solve elliptic linear elastic system.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1038
void WriteGeometry(const int i)
Write out a file in serial or directory in parallel containing new mesh geometry. ...