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