Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Forcing.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: Forcing.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: Abstract base class for forcing terms.
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
37 
38 namespace Nektar
39 {
40  namespace SolverUtils
41  {
43  {
44  typedef Loki::SingletonHolder<ForcingFactory,
45  Loki::CreateUsingNew,
46  Loki::NoDestroy > Type;
47  return Type::Instance();
48  }
49 
51  : m_session(pSession)
52  {
53 
54  }
55 
57  const Array<OneD, MultiRegions::ExpListSharedPtr>& pFields,
58  const unsigned int& pNumForcingFields,
59  const TiXmlElement* pForce)
60  {
61  v_InitObject(pFields, pNumForcingFields, pForce);
62  }
63 
64 
65  /**
66  * @param fields Expansion lists corresponding to input arrays
67  * @param inarray u^n from previous timestep
68  * @param outarray output array to append forcing to
69  */
71  const Array<OneD, MultiRegions::ExpListSharedPtr>& fields,
72  const Array<OneD, Array<OneD, NekDouble> >& inarray,
73  Array<OneD, Array<OneD, NekDouble> >& outarray,
74  const NekDouble& time)
75  {
76  v_Apply(fields, inarray, outarray, time);
77  }
78 
79 
80  /**
81  *
82  */
83  vector<ForcingSharedPtr> Forcing::Load(
85  const Array<OneD, MultiRegions::ExpListSharedPtr>& pFields,
86  const unsigned int& pNumForcingFields)
87  {
88  vector<ForcingSharedPtr> vForceList;
89 
90  if (!pSession->DefinesElement("Nektar/Forcing"))
91  {
92  return vForceList;
93  }
94 
95  TiXmlElement* vForcing = pSession->GetElement("Nektar/Forcing");
96  if (vForcing)
97  {
98  unsigned int vNumForcingFields = pNumForcingFields;
99  if (!pNumForcingFields)
100  {
101  vNumForcingFields = pFields.num_elements();
102  }
103 
104  TiXmlElement* vForce = vForcing->FirstChildElement("FORCE");
105  while (vForce)
106  {
107  string vType = vForce->Attribute("TYPE");
108 
109  vForceList.push_back(GetForcingFactory().CreateInstance(
110  vType, pSession, pFields,
111  vNumForcingFields, vForce));
112  vForce = vForce->NextSiblingElement("FORCE");
113  }
114  }
115  return vForceList;
116  }
117 
120  std::string pFieldName,
121  Array<OneD, NekDouble>& pArray,
122  const std::string& pFunctionName,
123  NekDouble pTime)
124  {
125  ASSERTL0(pSession->DefinesFunction(pFunctionName),
126  "Function '" + pFunctionName + "' does not exist.");
127 
129  pSession->GetFunction(pFunctionName, pFieldName);
130 
131  Array<OneD, NekDouble> x0(1,0.0);
132  Array<OneD, NekDouble> x1(1,0.0);
133  Array<OneD, NekDouble> x2(1,0.0);
134 
135  ffunc->Evaluate(x0, x1, x2, pTime, pArray);
136  }
137 
138 
140  Array<OneD, MultiRegions::ExpListSharedPtr> pFields,
142  std::string pFieldName,
143  Array<OneD, NekDouble>& pArray,
144  const std::string& pFunctionName,
145  NekDouble pTime)
146  {
147  ASSERTL0(pSession->DefinesFunction(pFunctionName),
148  "Function '" + pFunctionName + "' does not exist.");
149 
150  unsigned int nq = pFields[0]->GetNpoints();
151  if (pArray.num_elements() != nq)
152  {
153  pArray = Array<OneD, NekDouble> (nq);
154  }
155 
157  vType = pSession->GetFunctionType(pFunctionName, pFieldName);
159  {
160  Array<OneD, NekDouble> x0(nq);
161  Array<OneD, NekDouble> x1(nq);
162  Array<OneD, NekDouble> x2(nq);
163 
164  pFields[0]->GetCoords(x0, x1, x2);
166  pSession->GetFunction(pFunctionName, pFieldName);
167 
168  ffunc->Evaluate(x0, x1, x2, pTime, pArray);
169  }
170  else if (vType == LibUtilities::eFunctionTypeFile)
171  {
172  std::string filename = pSession->GetFunctionFilename(
173  pFunctionName,
174  pFieldName);
175 
176  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
177  std::vector<std::vector<NekDouble> > FieldData;
178  Array<OneD, NekDouble> vCoeffs(pFields[0]->GetNcoeffs());
179  Vmath::Zero(vCoeffs.num_elements(), vCoeffs, 1);
180 
183  fld->Import(filename, FieldDef, FieldData);
184 
185  int idx = -1;
186  for (int i = 0; i < FieldDef.size(); ++i)
187  {
188  for (int j = 0; j < FieldDef[i]->m_fields.size(); ++j)
189  {
190  if (FieldDef[i]->m_fields[j] == pFieldName)
191  {
192  idx = j;
193  }
194  }
195 
196  if (idx >= 0)
197  {
198  pFields[0]->ExtractDataToCoeffs(
199  FieldDef[i],
200  FieldData[i],
201  FieldDef[i]->m_fields[idx],
202  vCoeffs);
203  }
204  else
205  {
206  cout << "Field " + pFieldName + " not found." << endl;
207  }
208  }
209  pFields[0]->BwdTrans_IterPerExp(vCoeffs, pArray);
210  }
211  }
212 
213  }
214 }