Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 using namespace std;
50 
51 namespace Nektar
52 {
53 
54 string IterativeElasticSystem::className = GetEquationSystemFactory().
55  RegisterCreatorFunction("IterativeElasticSystem",
56  IterativeElasticSystem::create);
57 
58 IterativeElasticSystem::IterativeElasticSystem(
60  : LinearElasticSystem(pSession)
61 {
62 }
63 
65 {
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.num_elements(); ++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 
100  m_toDeform.size());
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,
114  bndCondExp->GetCoeffs(), 1,
115  bndCondExp->UpdateCoeffs(), 1);
116  Vmath::Vcopy(nCoeffs, bndCondExp->GetCoeffs(), 1,
117  m_savedBCs[i][j], 1);
118  }
119  }
120  }
121 }
122 
124 {
126 }
127 
129 {
130  int i, j, k;
131 
132  // Write initial geometry for consistency/script purposes
133  WriteGeometry(0);
134 
135  // Now loop over desired number of steps
136  for (i = 1; i <= m_numSteps; ++i)
137  {
138  int invalidElmtId = -1;
139 
140  // Perform solve for this iteration and update geometry accordingly.
143  WriteGeometry(i);
144 
145  // Check for invalid elements.
146  for (j = 0; j < m_fields[0]->GetExpSize(); ++j)
147  {
149  m_fields[0]->GetExp(j)->GetGeom()->GetGeomFactors();
150 
151  if (!geomFac->IsValid())
152  {
153  invalidElmtId =
154  m_fields[0]->GetExp(j)->GetGeom()->GetGlobalID();
155  break;
156  }
157  }
158 
159  m_session->GetComm()->AllReduce(invalidElmtId, LibUtilities::ReduceMax);
160 
161  // If we found an invalid element, exit loop without writing output.
162  if (invalidElmtId >= 0)
163  {
164  if (m_session->GetComm()->GetRank() == 0)
165  {
166  cout << "- Detected negative Jacobian in element "
167  << invalidElmtId << "; terminating at"
168  " step: "<<i<< endl;
169  }
170 
171  break;
172  }
173 
174  if (m_session->GetComm()->GetRank() == 0)
175  {
176  cout << "Step: " << i << endl;
177  }
178 
179  // Update boundary conditions
180  if (m_repeatBCs)
181  {
182  for (j = 0; j < m_fields.num_elements(); ++j)
183  {
184  string varName = m_session->GetVariable(j);
185  m_fields[j]->EvaluateBoundaryConditions(m_time, varName);
186  }
187  }
188  else
189  {
190  for (j = 0; j < m_fields.num_elements(); ++j)
191  {
193  &bndCondExp = m_fields[j]->GetBndCondExpansions();
194 
195  for (k = 0; k < m_toDeform.size(); ++k)
196  {
197  const int id = m_toDeform[k];
198  const int nCoeffs = bndCondExp[id]->GetNcoeffs();
199  Vmath::Vcopy(nCoeffs,
200  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  if(!fs::is_directory(s.str()))
223  {
224  fs::create_directory(s.str());
225  }
226 
227  boost::format pad("P%1$07d.xml");
228  pad % m_session->GetComm()->GetRank();
229  filename = fs::path(s.str()) / fs::path(pad.str());
230  }
231  else
232  {
233  s << ".xml";
234  filename = fs::path(s.str());
235  }
236 
237  string fname = LibUtilities::PortablePath(filename);
238  m_fields[0]->GetGraph()->WriteGeometry(fname);
239 }
240 
241 }
#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
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:59
int m_numSteps
Number of steps to split deformation across.
STL namespace.
virtual void v_DoSolve()
Solve elliptic linear elastic system.
bool m_repeatBCs
Flag determining whether to repeat boundary conditions.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
LibUtilities::CommSharedPtr m_comm
Communicator.
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_savedBCs
Storage for boundary conditions.
Base class for linear elastic system.
std::vector< int > m_toDeform
Vector of boundary regions to deform.
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.
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:1047
void WriteGeometry(const int i)
Write out a file in serial or directory in parallel containing new mesh geometry. ...