Nektar++
Mapping.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: Mapping.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: Abstract base class for mappings.
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
35 #include <boost/core/ignore_unused.hpp>
36 
37 #include <GlobalMapping/Mapping.h>
38 
39 #include <boost/algorithm/string/predicate.hpp>
40 
41 using namespace std;
42 
43 namespace Nektar
44 {
45 namespace GlobalMapping
46 {
47 
48 MappingSharedPtr Mapping::m_mappingPtr = MappingSharedPtr();
49 bool Mapping::m_init = false;
50 bool Mapping::m_isDefined = false;
51 
53 {
54  static MappingFactory instance;
55  return instance;
56 }
57 
58 Mapping::Mapping(const LibUtilities::SessionReaderSharedPtr &pSession,
60  : m_session(pSession), m_fields(pFields)
61 {
62  switch (m_fields[0]->GetExpType())
63  {
64  case MultiRegions::e1D:
65  {
67  }
68  break;
69 
70  case MultiRegions::e2D:
71  {
73  }
74  break;
75 
76  case MultiRegions::e3D:
79  {
81  }
82  break;
83 
84  default:
85  ASSERTL0(0, "Dimension not supported");
86  break;
87  }
88 
90 }
91 
92 /**
93  * This function initialises the Mapping object. It computes the coordinates
94  * and velocity coordinates, initialises the workspace variables, and calls
95  * UpdateGeomInfo, which will perform the calculations specific for each type
96  * of Mapping.
97  * @param pFields ExpList array used in the mapping
98  * @param pMapping xml element describing the mapping
99  */
102  const TiXmlElement *pMapping)
103 {
104  boost::ignore_unused(pFields);
105 
106  int phystot = m_fields[0]->GetTotPoints();
107  m_fromFunction = true;
108  // Initialise variables
112  for (int i = 0; i < 3; i++)
113  {
114  m_coords[i] = Array<OneD, NekDouble>(phystot);
115  m_coordsVel[i] = Array<OneD, NekDouble>(phystot);
116  coords[i] = Array<OneD, NekDouble>(phystot);
117  }
118 
119  // Check if mapping is defined as time-dependent
120  const TiXmlElement *timeDep = pMapping->FirstChildElement("TIMEDEPENDENT");
121  if (timeDep)
122  {
123  string sTimeDep = timeDep->GetText();
124  m_timeDependent = (boost::iequals(sTimeDep, "true")) ||
125  (boost::iequals(sTimeDep, "yes"));
126  }
127  else
128  {
129  m_timeDependent = false;
130  }
131 
132  // Load coordinates
133  string fieldNames[3] = {"x", "y", "z"};
134  const TiXmlElement *funcNameElmt = pMapping->FirstChildElement("COORDS");
135  if (funcNameElmt)
136  {
137  m_funcName = funcNameElmt->GetText();
138  ASSERTL0(m_session->DefinesFunction(m_funcName),
139  "Function '" + m_funcName + "' not defined.");
140 
141  // Get coordinates in the domain
142  m_fields[0]->GetCoords(coords[0], coords[1], coords[2]);
143 
144  std::string s_FieldStr;
145  // Check if function from session file defines each component
146  // and evaluate them, otherwise use trivial transformation
147  for (int i = 0; i < 3; i++)
148  {
149  s_FieldStr = fieldNames[i];
150  if (m_session->DefinesFunction(m_funcName, s_FieldStr))
151  {
152  EvaluateFunction(m_fields, m_session, s_FieldStr, m_coords[i],
153  m_funcName);
154  if (i == 2 && m_fields[0]->GetExpType() == MultiRegions::e3DH1D)
155  {
156  ASSERTL0(
157  false,
158  "3DH1D does not support mapping in the z-direction.");
159  }
160  }
161  else
162  {
163  // This coordinate is not defined, so use (x^i)' = x^i
164  Vmath::Vcopy(phystot, coords[i], 1, m_coords[i], 1);
165  }
166  }
167  }
168  else
169  {
170  m_fields[0]->GetCoords(coords[0], coords[1], coords[2]);
171  for (int i = 0; i < 3; i++)
172  {
173  // Use (x^i)' = x^i as default. This can be useful if we
174  // have a time-dependent mapping, and then only the
175  // initial mapping will be trivial
176  Vmath::Vcopy(phystot, coords[i], 1, m_coords[i], 1);
177  }
178  }
179 
180  // Load coordinate velocity if they are defined,
181  // otherwise use zero to make it general
182  string velFieldNames[3] = {"vx", "vy", "vz"};
183  const TiXmlElement *velFuncNameElmt = pMapping->FirstChildElement("VEL");
184  if (velFuncNameElmt)
185  {
186  m_velFuncName = velFuncNameElmt->GetText();
187  ASSERTL0(m_session->DefinesFunction(m_velFuncName),
188  "Function '" + m_velFuncName + "' not defined.");
189 
190  std::string s_FieldStr;
191  // Check if function from session file defines each component
192  // and evaluate them, otherwise use 0
193  for (int i = 0; i < 3; i++)
194  {
195  s_FieldStr = velFieldNames[i];
196  if (m_session->DefinesFunction(m_velFuncName, s_FieldStr))
197  {
198  EvaluateFunction(m_fields, m_session, s_FieldStr,
200  if (i == 2 && m_fields[0]->GetExpType() == MultiRegions::e3DH1D)
201  {
202  ASSERTL0(
203  false,
204  "3DH1D does not support mapping in the z-direction.");
205  }
206  }
207  else
208  {
209  // This coordinate velocity is not defined, so use 0
210  Vmath::Zero(phystot, m_coordsVel[i], 1);
211  }
212  }
213  }
214  else
215  {
216  for (int i = 0; i < 3; i++)
217  {
218  Vmath::Zero(phystot, m_coordsVel[i], 1);
219  }
220  }
221 
222  // Initialise workspace variables
223  int nvel = m_nConvectiveFields;
227  for (int i = 0; i < nvel; i++)
228  {
229  m_tmp[i] = Array<OneD, NekDouble>(phystot, 0.0);
230  for (int j = 0; j < nvel; j++)
231  {
232  m_wk1[i * nvel + j] = Array<OneD, NekDouble>(phystot, 0.0);
233  m_wk2[i * nvel + j] = Array<OneD, NekDouble>(phystot, 0.0);
234  }
235  }
236 
237  // Calculate information required by the particular mapping
238  UpdateGeomInfo();
239 }
240 
241 /**
242  * This function replaces the expansion list contained in m_fields, and then
243  * proceeds reinitialising the Mapping with this new field.
244  * @param pFields New field to be used by the Mapping
245  */
248 {
249  m_fields = pFields;
250 
251  TiXmlElement *vMapping = NULL;
252 
253  if (m_session->DefinesElement("Nektar/Mapping"))
254  {
255  vMapping = m_session->GetElement("Nektar/Mapping");
256  }
257  InitObject(pFields, vMapping);
258 }
259 
260 /**
261  * This function is responsible for loading the Mapping, guaranteeing that a
262  * single instance of this class exists. When it is first called, it creates a
263  * Mapping and returns a pointer to it. On subsequent calls, it just returns
264  * the pointer.
265  * @param pSession Session reader object
266  * @param pFields Fields which will be used by the Mapping
267  * @return Pointer to the Mapping
268  */
270  const LibUtilities::SessionReaderSharedPtr &pSession,
272 {
273  if (!m_init)
274  {
275  TiXmlElement *vMapping = NULL;
276  string vType;
277  if (pSession->DefinesElement("Nektar/Mapping"))
278  {
279  vMapping = pSession->GetElement("Nektar/Mapping");
280  vType = vMapping->Attribute("TYPE");
281  m_isDefined = true;
282  }
283  else
284  {
285  vType = "Translation";
286  }
287 
288  m_mappingPtr = GetMappingFactory().CreateInstance(vType, pSession,
289  pFields, vMapping);
290 
291  m_init = true;
292  }
293 
294  return m_mappingPtr;
295 }
296 
297 /**
298  * This function should be called before writing a fld or chk file. It updates
299  * the metadata with information from the Mapping. Also, if the mapping is not
300  * defined by a function, it writes a .map file containing the coordinates and
301  * velocity of the coordinates.
302  * @param fieldMetaDataMap Metadata of the output file
303  * @param outname Name of the output file
304  */
306  const std::string &outname)
307 {
308  // Only do anything if mapping exists
309  if (m_isDefined)
310  {
311  fieldMetaDataMap["MappingCartesianVel"] = std::string("False");
312  if (m_fromFunction)
313  {
314  // Add metadata
315  fieldMetaDataMap["MappingType"] = std::string("Expression");
316  fieldMetaDataMap["MappingExpression"] = m_funcName;
317  if (m_timeDependent)
318  {
319  fieldMetaDataMap["MappingVelExpression"] = m_velFuncName;
320  }
321  }
322  else
323  {
324  int expdim = m_fields[0]->GetGraph()->GetMeshDimension();
325  string fieldNames[3] = {"x", "y", "z"};
326  string velFieldNames[3] = {"vx", "vy", "vz"};
327 
328  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef =
329  m_fields[0]->GetFieldDefinitions();
330  std::vector<std::vector<NekDouble>> FieldData(FieldDef.size());
331 
332  int ncoeffs = m_fields[0]->GetNcoeffs();
333  Array<OneD, NekDouble> fieldcoeffs(ncoeffs);
334 
335  bool wavespace = m_fields[0]->GetWaveSpace();
336  m_fields[0]->SetWaveSpace(false);
337  // copy coordinates Data into FieldData and set variable
338  for (int j = 0; j < expdim; ++j)
339  {
340  m_fields[0]->FwdTransLocalElmt(m_coords[j], fieldcoeffs);
341 
342  for (int i = 0; i < FieldDef.size(); ++i)
343  {
344  // Could do a search here to find correct variable
345  FieldDef[i]->m_fields.push_back(fieldNames[j]);
346  m_fields[0]->AppendFieldData(FieldDef[i], FieldData[i],
347  fieldcoeffs);
348  }
349  }
350  if (m_timeDependent)
351  {
352  // copy coordinates velocity Data into FieldData and set
353  // variable
354  for (int j = 0; j < expdim; ++j)
355  {
356  m_fields[0]->FwdTrans(m_coordsVel[j], fieldcoeffs);
357 
358  for (int i = 0; i < FieldDef.size(); ++i)
359  {
360  // Could do a search here to find correct variable
361  FieldDef[i]->m_fields.push_back(velFieldNames[j]);
362  m_fields[0]->AppendFieldData(FieldDef[i], FieldData[i],
363  fieldcoeffs);
364  }
365  }
366  }
367 
368  std::string outfile = outname;
369  outfile.erase(outfile.end() - 4, outfile.end());
370  outfile += ".map";
371 
372  m_fld->Write(outfile, FieldDef, FieldData, fieldMetaDataMap);
373 
374  // Write metadata to orginal output
375  fieldMetaDataMap["MappingType"] = std::string("File");
376  fieldMetaDataMap["FileName"] = outfile;
377 
378  m_fields[0]->SetWaveSpace(wavespace);
379  }
380  }
381 }
382 
384  LibUtilities::SessionReaderSharedPtr pSession, std::string pFieldName,
385  Array<OneD, NekDouble> &pArray, const std::string &pFunctionName,
386  NekDouble pTime)
387 {
388  ASSERTL0(pSession->DefinesFunction(pFunctionName),
389  "Function '" + pFunctionName + "' does not exist.");
390 
392  pSession->GetFunction(pFunctionName, pFieldName);
393 
394  Array<OneD, NekDouble> x0(1, 0.0);
395  Array<OneD, NekDouble> x1(1, 0.0);
396  Array<OneD, NekDouble> x2(1, 0.0);
397 
398  ffunc->Evaluate(x0, x1, x2, pTime, pArray);
399 }
400 
403  LibUtilities::SessionReaderSharedPtr pSession, std::string pFieldName,
404  Array<OneD, NekDouble> &pArray, const std::string &pFunctionName,
405  NekDouble pTime)
406 {
407  ASSERTL0(pSession->DefinesFunction(pFunctionName),
408  "Function '" + pFunctionName + "' does not exist.");
409 
410  unsigned int nq = pFields[0]->GetNpoints();
411  if (pArray.size() != nq)
412  {
413  pArray = Array<OneD, NekDouble>(nq);
414  }
415 
417  vType = pSession->GetFunctionType(pFunctionName, pFieldName);
419  {
420  Array<OneD, NekDouble> x0(nq);
421  Array<OneD, NekDouble> x1(nq);
422  Array<OneD, NekDouble> x2(nq);
423 
424  pFields[0]->GetCoords(x0, x1, x2);
426  pSession->GetFunction(pFunctionName, pFieldName);
427 
428  ffunc->Evaluate(x0, x1, x2, pTime, pArray);
429  }
430  else if (vType == LibUtilities::eFunctionTypeFile)
431  {
432  std::string filename =
433  pSession->GetFunctionFilename(pFunctionName, pFieldName);
434 
435  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
436  std::vector<std::vector<NekDouble>> FieldData;
437  Array<OneD, NekDouble> vCoeffs(pFields[0]->GetNcoeffs());
438  Vmath::Zero(vCoeffs.size(), vCoeffs, 1);
439 
441  LibUtilities::FieldIO::CreateForFile(pSession, filename);
442  fld->Import(filename, FieldDef, FieldData);
443 
444  int idx = -1;
445  for (int i = 0; i < FieldDef.size(); ++i)
446  {
447  for (int j = 0; j < FieldDef[i]->m_fields.size(); ++j)
448  {
449  if (FieldDef[i]->m_fields[j] == pFieldName)
450  {
451  idx = j;
452  }
453  }
454 
455  if (idx >= 0)
456  {
457  pFields[0]->ExtractDataToCoeffs(FieldDef[i], FieldData[i],
458  FieldDef[i]->m_fields[idx],
459  vCoeffs);
460  }
461  else
462  {
463  cout << "Field " + pFieldName + " not found." << endl;
464  }
465  }
466  pFields[0]->BwdTrans(vCoeffs, pArray);
467  }
468 }
469 
470 /**
471  * This function converts a contravariant vector in transformed space
472  * \f$v^{j}\f$ to the corresponding \f$\bar{v}^{i}\f$ in Cartesian (physical)
473  * space using the relation
474  * \f[\bar{v}^{i} = \frac{\partial \bar{x}^i}{\partial x^j}v^{j}\f]
475  *
476  * @param inarray Components of the vector in transformed space (\f$v^{j}\f$)
477  * @param outarray Components of the vector in Cartesian space
478  * (\f$\bar{v}^{i}\f$)
479  */
481  const Array<OneD, Array<OneD, NekDouble>> &inarray,
482  Array<OneD, Array<OneD, NekDouble>> &outarray)
483 {
484  if (inarray == outarray)
485  {
486  int physTot = m_fields[0]->GetTotPoints();
487  int nvel = m_nConvectiveFields;
488 
489  for (int i = 0; i < nvel; i++)
490  {
491  Vmath::Vcopy(physTot, inarray[i], 1, m_tmp[i], 1);
492  }
493  v_ContravarToCartesian(m_tmp, outarray);
494  }
495  else
496  {
497  v_ContravarToCartesian(inarray, outarray);
498  }
499 }
500 
501 /**
502  * This function converts a covariant vector in transformed space
503  * \f$v_{j}\f$ to the corresponding \f$\bar{v}_{i}\f$ in Cartesian (physical)
504  * space using the relation
505  * \f[\bar{v}_{i} = \frac{\partial x^j}{\partial \bar{x}^i}v_{j}\f]
506  *
507  * @param inarray Components of the vector in transformed space (\f$v_{j}\f$)
508  * @param outarray Components of the vector in Cartesian space
509  * (\f$\bar{v}_{i}\f$)
510  */
512  const Array<OneD, Array<OneD, NekDouble>> &inarray,
513  Array<OneD, Array<OneD, NekDouble>> &outarray)
514 {
515  if (inarray == outarray)
516  {
517  int physTot = m_fields[0]->GetTotPoints();
518  int nvel = m_nConvectiveFields;
519 
520  for (int i = 0; i < nvel; i++)
521  {
522  Vmath::Vcopy(physTot, inarray[i], 1, m_tmp[i], 1);
523  }
524  v_CovarToCartesian(m_tmp, outarray);
525  }
526  else
527  {
528  v_CovarToCartesian(inarray, outarray);
529  }
530 }
531 
532 /**
533  * This function converts a contravariant vector in Cartesian space
534  * \f$\bar{v}^{i}\f$ to the corresponding \f$v^{j}\f$ in transformed
535  * space using the relation
536  * \f[v^{j} = \frac{\partial x^j}{\partial \bar{x}^i}\bar{v}^{i}\f]
537  *
538  * @param inarray Components of the vector in Cartesian space
539  * (\f$\bar{v}^{i}\f$)
540  * @param outarray Components of the vector in transformed space (\f$v^{j}\f$)
541  */
543  const Array<OneD, Array<OneD, NekDouble>> &inarray,
544  Array<OneD, Array<OneD, NekDouble>> &outarray)
545 {
546  if (inarray == outarray)
547  {
548  int physTot = m_fields[0]->GetTotPoints();
549  int nvel = m_nConvectiveFields;
550 
551  for (int i = 0; i < nvel; i++)
552  {
553  Vmath::Vcopy(physTot, inarray[i], 1, m_tmp[i], 1);
554  }
555  v_ContravarFromCartesian(m_tmp, outarray);
556  }
557  else
558  {
559  v_ContravarFromCartesian(inarray, outarray);
560  }
561 }
562 
563 /**
564  * This function converts a covariant vector in Cartesian space
565  * \f$\bar{v}_{i}\f$ to the corresponding \f$v_{j}\f$ in transformed
566  * space using the relation
567  * \f[\bar{v}_{j} = \frac{\partial \bar{x}^i}{\partial x^j}v_{i}\f]
568  *
569  * @param inarray Components of the vector in Cartesian space
570  * (\f$\bar{v}_{i}\f$)
571  * @param outarray Components of the vector in transformed space (\f$v_{j}\f$)
572  */
574  const Array<OneD, Array<OneD, NekDouble>> &inarray,
575  Array<OneD, Array<OneD, NekDouble>> &outarray)
576 {
577  if (inarray == outarray)
578  {
579  int physTot = m_fields[0]->GetTotPoints();
580  int nvel = m_nConvectiveFields;
581 
582  for (int i = 0; i < nvel; i++)
583  {
584  Vmath::Vcopy(physTot, inarray[i], 1, m_tmp[i], 1);
585  }
586  v_CovarFromCartesian(m_tmp, outarray);
587  }
588  else
589  {
590  v_CovarFromCartesian(inarray, outarray);
591  }
592 }
593 
594 /**
595  * This function lowers the index of the contravariant vector \f$v^{i}\f$,
596  * transforming it into its associated covariant vector \f$v_{j}\f$
597  * according to the relation
598  * \f[v_{j} = g_{ij}v^{i}\f]
599  * where \f$g_{ij}\f$ is the metric tensor.
600  *
601  * @param inarray Components of the contravariant vector \f$v^{i}\f$
602  * @param outarray Components of the covariant vector \f$v_{j}\f$
603  */
605  Array<OneD, Array<OneD, NekDouble>> &outarray)
606 {
607  if (inarray == outarray)
608  {
609  int physTot = m_fields[0]->GetTotPoints();
610  int nvel = m_nConvectiveFields;
611 
612  for (int i = 0; i < nvel; i++)
613  {
614  Vmath::Vcopy(physTot, inarray[i], 1, m_tmp[i], 1);
615  }
616  v_LowerIndex(m_tmp, outarray);
617  }
618  else
619  {
620  v_LowerIndex(inarray, outarray);
621  }
622 }
623 
624 /**
625  * This function raises the index of the covariant vector \f$v_{j}\f$,
626  * transforming it into its associated contravariant vector \f$v^{i}\f$
627  * according to the relation
628  * \f[v^{i} = g^{ij}v_{j}\f]
629  * where \f$g^{ij}\f$ is the inverse of the metric tensor.
630  *
631  * @param inarray Components of the contravariant vector \f$v^{i}\f$
632  * @param outarray Components of the covariant vector \f$v_{j}\f$
633  */
635  Array<OneD, Array<OneD, NekDouble>> &outarray)
636 {
637  if (inarray == outarray)
638  {
639  int physTot = m_fields[0]->GetTotPoints();
640  int nvel = m_nConvectiveFields;
641 
642  for (int i = 0; i < nvel; i++)
643  {
644  Vmath::Vcopy(physTot, inarray[i], 1, m_tmp[i], 1);
645  }
646  v_RaiseIndex(m_tmp, outarray);
647  }
648  else
649  {
650  v_RaiseIndex(inarray, outarray);
651  }
652 }
653 
657 {
658  int physTot = m_fields[0]->GetTotPoints();
659 
660  out0 = Array<OneD, NekDouble>(physTot, 0.0);
661  out1 = Array<OneD, NekDouble>(physTot, 0.0);
662  out2 = Array<OneD, NekDouble>(physTot, 0.0);
663 
664  Vmath::Vcopy(physTot, m_coords[0], 1, out0, 1);
665  Vmath::Vcopy(physTot, m_coords[1], 1, out1, 1);
666  Vmath::Vcopy(physTot, m_coords[2], 1, out2, 1);
667 }
668 
670 {
671  int physTot = m_fields[0]->GetTotPoints();
672 
673  for (int i = 0; i < m_nConvectiveFields; ++i)
674  {
675  outarray[i] = Array<OneD, NekDouble>(physTot, 0.0);
676  Vmath::Vcopy(physTot, m_coordsVel[i], 1, outarray[i], 1);
677  }
678 }
679 
681  const Array<OneD, Array<OneD, NekDouble>> &inarray,
682  Array<OneD, NekDouble> &outarray)
683 {
684  int physTot = m_fields[0]->GetTotPoints();
685 
686  outarray = Array<OneD, NekDouble>(physTot, 0.0);
687  if (!HasConstantJacobian())
688  {
689  // Set wavespace to false and store current value
690  bool wavespace = m_fields[0]->GetWaveSpace();
691  m_fields[0]->SetWaveSpace(false);
692 
693  // Get Mapping Jacobian
694  Array<OneD, NekDouble> Jac(physTot, 0.0);
695  GetJacobian(Jac);
696 
697  // Calculate inarray . grad(Jac)
698  Array<OneD, NekDouble> wk(physTot, 0.0);
699  for (int i = 0; i < m_nConvectiveFields; ++i)
700  {
701  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[i], Jac, wk);
702  Vmath::Vvtvp(physTot, inarray[i], 1, wk, 1, outarray, 1, outarray,
703  1);
704  }
705  m_fields[0]->SetWaveSpace(wavespace);
706  }
707 }
708 
710  Array<OneD, Array<OneD, NekDouble>> &outarray)
711 {
712  int physTot = m_fields[0]->GetTotPoints();
713  int nvel = m_nConvectiveFields;
714 
715  Array<OneD, Array<OneD, NekDouble>> g(nvel * nvel);
716 
717  GetMetricTensor(g);
718 
719  for (int i = 0; i < nvel; i++)
720  {
721  outarray[i] = Array<OneD, NekDouble>(physTot, 0.0);
722  for (int j = 0; j < nvel; j++)
723  {
724  Vmath::Vvtvp(physTot, g[i * nvel + j], 1, inarray[j], 1,
725  outarray[i], 1, outarray[i], 1);
726  }
727  }
728 }
729 
731  Array<OneD, Array<OneD, NekDouble>> &outarray)
732 {
733  int physTot = m_fields[0]->GetTotPoints();
734  int nvel = m_nConvectiveFields;
735 
736  Array<OneD, Array<OneD, NekDouble>> g(nvel * nvel);
737 
739 
740  for (int i = 0; i < nvel; i++)
741  {
742  outarray[i] = Array<OneD, NekDouble>(physTot, 0.0);
743  for (int j = 0; j < nvel; j++)
744  {
745  Vmath::Vvtvp(physTot, g[i * nvel + j], 1, inarray[j], 1,
746  outarray[i], 1, outarray[i], 1);
747  }
748  }
749 }
750 
752  Array<OneD, NekDouble> &outarray)
753 {
754  int physTot = m_fields[0]->GetTotPoints();
755  Array<OneD, NekDouble> wk(physTot, 0.0);
756 
757  Vmath::Zero(physTot, outarray, 1);
758 
759  // Set wavespace to false and store current value
760  bool wavespace = m_fields[0]->GetWaveSpace();
761  m_fields[0]->SetWaveSpace(false);
762 
763  // Get Mapping Jacobian
764  Array<OneD, NekDouble> Jac(physTot, 0.0);
765  GetJacobian(Jac);
766 
767  for (int i = 0; i < m_nConvectiveFields; ++i)
768  {
769  Vmath::Vmul(physTot, Jac, 1, inarray[i], 1, wk, 1); // J*Ui
770  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[i], wk,
771  wk); // (J*Ui)_i
772  Vmath::Vadd(physTot, wk, 1, outarray, 1, outarray, 1);
773  }
774  Vmath::Vdiv(physTot, outarray, 1, Jac, 1, outarray, 1); // 1/J*(J*Ui)_i
775 
776  m_fields[0]->SetWaveSpace(wavespace);
777 }
778 
780  const Array<OneD, Array<OneD, NekDouble>> &inarray,
781  Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble alpha)
782 {
783  int physTot = m_fields[0]->GetTotPoints();
784  int nvel = m_nConvectiveFields;
785 
787 
788  // Set wavespace to false and store current value
789  bool wavespace = m_fields[0]->GetWaveSpace();
790  m_fields[0]->SetWaveSpace(false);
791 
792  // Calculate vector gradient wk2 = u^i_(,k) = du^i/dx^k + {i,jk}*u^j
794  for (int i = 0; i < nvel; i++)
795  {
796  if (nvel == 2)
797  {
798  m_fields[0]->PhysDeriv(inarray[i], m_wk2[i * nvel + 0],
799  m_wk2[i * nvel + 1]);
800  }
801  else
802  {
803  m_fields[0]->PhysDeriv(inarray[i], m_wk2[i * nvel + 0],
804  m_wk2[i * nvel + 1], m_wk2[i * nvel + 2]);
805  }
806  for (int k = 0; k < nvel; k++)
807  {
808  Vmath::Vadd(physTot, m_wk1[i * nvel + k], 1, m_wk2[i * nvel + k], 1,
809  m_wk1[i * nvel + k], 1);
810  }
811  }
812  // Calculate wk1 = A^(ij) = g^(jk)*u^i_(,k)
813  for (int i = 0; i < nvel; i++)
814  {
815  for (int k = 0; k < nvel; k++)
816  {
817  tmp[k] = m_wk1[i * nvel + k];
818  }
819  RaiseIndex(tmp, m_tmp);
820  for (int j = 0; j < nvel; j++)
821  {
822  Vmath::Vcopy(physTot, m_tmp[j], 1, m_wk1[i * nvel + j], 1);
823  }
824  }
825  //
826  // Calculate L(U)^i = (A^(ij))_(,j) - alpha*d^2(u^i)/dx^jdx^j
827  //
828 
829  // Step 1 :
830  // d(A^(ij) - alpha*du^i/dx^j)/d(x^j)
831  for (int i = 0; i < nvel; i++)
832  {
833  outarray[i] = Array<OneD, NekDouble>(physTot, 0.0);
834  for (int j = 0; j < nvel; j++)
835  {
836  Vmath::Smul(physTot, alpha, m_wk2[i * nvel + j], 1, m_tmp[0], 1);
837  Vmath::Vsub(physTot, m_wk1[i * nvel + j], 1, m_tmp[0], 1, m_tmp[0],
838  1);
839 
840  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[j], m_tmp[0],
841  m_tmp[1]);
842  Vmath::Vadd(physTot, outarray[i], 1, m_tmp[1], 1, outarray[i], 1);
843  }
844  }
845 
846  // Step 2: d(A^(ij))/d(x^j) + {j,pj}*A^(ip)
847  for (int i = 0; i < nvel; i++)
848  {
849  for (int p = 0; p < nvel; p++)
850  {
851  tmp[p] = m_wk1[i * nvel + p];
852  }
854  for (int j = 0; j < nvel; j++)
855  {
856  Vmath::Vadd(physTot, outarray[i], 1, m_wk2[j * nvel + j], 1,
857  outarray[i], 1);
858  }
859  }
860 
861  // Step 3: d(A^(ij))/d(x^j) + {j,pj}*A^(ip) + {i,pj} A^(pj)
862  for (int j = 0; j < nvel; j++)
863  {
864  for (int p = 0; p < nvel; p++)
865  {
866  tmp[p] = m_wk1[p * nvel + j];
867  }
869  for (int i = 0; i < nvel; i++)
870  {
871  Vmath::Vadd(physTot, outarray[i], 1, m_wk2[i * nvel + j], 1,
872  outarray[i], 1);
873  }
874  }
875 
876  // Restore value of wavespace
877  m_fields[0]->SetWaveSpace(wavespace);
878 }
879 
881  Array<OneD, Array<OneD, NekDouble>> &outarray)
882 {
883  int physTot = m_fields[0]->GetTotPoints();
884  int nvel = m_nConvectiveFields;
885 
886  // Declare variables
887  outarray = Array<OneD, Array<OneD, NekDouble>>(nvel * nvel * nvel);
889  for (int i = 0; i < nvel * nvel * nvel; i++)
890  {
891  outarray[i] = Array<OneD, NekDouble>(physTot, 0.0);
892  }
893 
894  // Set wavespace to false and store current value
895  bool wavespace = m_fields[0]->GetWaveSpace();
896  m_fields[0]->SetWaveSpace(false);
897 
898  // Calculate vector gradient u^i_(,j) = du^i/dx^j + {i,pj}*u^p
900  for (int i = 0; i < nvel; i++)
901  {
902  if (nvel == 2)
903  {
904  m_fields[0]->PhysDeriv(inarray[i], m_wk2[i * nvel + 0],
905  m_wk2[i * nvel + 1]);
906  }
907  else
908  {
909  m_fields[0]->PhysDeriv(inarray[i], m_wk2[i * nvel + 0],
910  m_wk2[i * nvel + 1], m_wk2[i * nvel + 2]);
911  }
912  for (int j = 0; j < nvel; j++)
913  {
914  Vmath::Vadd(physTot, m_wk1[i * nvel + j], 1, m_wk2[i * nvel + j], 1,
915  m_wk1[i * nvel + j], 1);
916  }
917  }
918 
919  //
920  // Calculate (u^i_,j),k
921  //
922 
923  // Step 1 : d(u^i_,j))/d(x^k)
924  for (int i = 0; i < nvel; i++)
925  {
926  for (int j = 0; j < nvel; j++)
927  {
928  if (nvel == 2)
929  {
930  m_fields[0]->PhysDeriv(
931  m_wk1[i * nvel + j],
932  outarray[i * nvel * nvel + j * nvel + 0],
933  outarray[i * nvel * nvel + j * nvel + 1]);
934  }
935  else
936  {
937  m_fields[0]->PhysDeriv(
938  m_wk1[i * nvel + j],
939  outarray[i * nvel * nvel + j * nvel + 0],
940  outarray[i * nvel * nvel + j * nvel + 1],
941  outarray[i * nvel * nvel + j * nvel + 2]);
942  }
943  }
944  }
945 
946  // Step 2: d(u^i_,j)/d(x^k) - {p,jk}*u^i_,p
947  for (int i = 0; i < nvel; i++)
948  {
949  for (int p = 0; p < nvel; p++)
950  {
951  tmp[p] = m_wk1[i * nvel + p];
952  }
954  for (int j = 0; j < nvel; j++)
955  {
956  for (int k = 0; k < nvel; k++)
957  {
958  Vmath::Vsub(physTot, outarray[i * nvel * nvel + j * nvel + k],
959  1, m_wk2[j * nvel + k], 1,
960  outarray[i * nvel * nvel + j * nvel + k], 1);
961  }
962  }
963  }
964 
965  // Step 3: d(u^i_,j)/d(x^k) - {p,jk}*u^i_,p + {i,pk} u^p_,j
966  for (int j = 0; j < nvel; j++)
967  {
968  for (int p = 0; p < nvel; p++)
969  {
970  tmp[p] = m_wk1[p * nvel + j];
971  }
973  for (int i = 0; i < nvel; i++)
974  {
975  for (int k = 0; k < nvel; k++)
976  {
977  Vmath::Vadd(physTot, outarray[i * nvel * nvel + j * nvel + k],
978  1, m_wk2[i * nvel + k], 1,
979  outarray[i * nvel * nvel + j * nvel + k], 1);
980  }
981  }
982  }
983 
984  // Restore value of wavespace
985  m_fields[0]->SetWaveSpace(wavespace);
986 }
987 
989  Array<OneD, Array<OneD, NekDouble>> &outarray,
990  const bool generalized)
991 {
992  int physTot = m_fields[0]->GetTotPoints();
993  int nvel = m_nConvectiveFields;
994 
995  // Set wavespace to false and store current value
996  bool wavespace = m_fields[0]->GetWaveSpace();
997  m_fields[0]->SetWaveSpace(false);
998 
999  // For implicit treatment of viscous terms, we want the generalized
1000  // curlcurl and for explicit treatment, we want the cartesian one.
1001  if (generalized)
1002  {
1003  // Get the second derivatives u^i_{,jk}
1005  Array<OneD, Array<OneD, NekDouble>> ddU(nvel * nvel * nvel);
1006  gradgradU(inarray, ddU);
1007 
1008  // Raise index to obtain A^{ip}_{k} = g^pj u^i_{,jk}
1009  for (int i = 0; i < nvel; ++i)
1010  {
1011  for (int k = 0; k < nvel; ++k)
1012  {
1013  // Copy to wk
1014  for (int j = 0; j < nvel; ++j)
1015  {
1016  tmp[j] = ddU[i * nvel * nvel + j * nvel + k];
1017  }
1018  RaiseIndex(tmp, m_tmp);
1019  for (int p = 0; p < nvel; ++p)
1020  {
1021  Vmath::Vcopy(physTot, m_tmp[p], 1,
1022  ddU[i * nvel * nvel + p * nvel + k], 1);
1023  }
1024  }
1025  }
1026  // The curlcurl is g^ji u^k_{kj} - g^jk u^i_kj = A^{ki}_k - A^{ik}_k
1027  for (int i = 0; i < nvel; ++i)
1028  {
1029  outarray[i] = Array<OneD, NekDouble>(physTot, 0.0);
1030  for (int k = 0; k < nvel; ++k)
1031  {
1032  Vmath::Vadd(physTot, outarray[i], 1,
1033  ddU[k * nvel * nvel + i * nvel + k], 1, outarray[i],
1034  1);
1035  Vmath::Vsub(physTot, outarray[i], 1,
1036  ddU[i * nvel * nvel + k * nvel + k], 1, outarray[i],
1037  1);
1038  }
1039  }
1040  }
1041  else
1042  {
1043  m_fields[0]->CurlCurl(inarray, outarray);
1044  }
1045 
1046  // Restore value of wavespace
1047  m_fields[0]->SetWaveSpace(wavespace);
1048 }
1049 
1051 {
1052  int physTot = m_fields[0]->GetTotPoints();
1053  int nvel = m_nConvectiveFields;
1054  int nfields = m_fields.size();
1055  int nbnds = m_fields[0]->GetBndConditions().size();
1056 
1057  // Declare variables
1060 
1061  Array<OneD, bool> isDirichlet(nfields);
1062 
1063  Array<OneD, Array<OneD, NekDouble>> values(nfields);
1064  for (int i = 0; i < nfields; i++)
1065  {
1066  values[i] = Array<OneD, NekDouble>(physTot, 0.0);
1067  }
1068 
1071  Array<OneD, Array<OneD, NekDouble>> coordVel(nvel);
1072  for (int i = 0; i < nvel; i++)
1073  {
1074  tmp[i] = Array<OneD, NekDouble>(physTot, 0.0);
1075  tmp2[i] = Array<OneD, NekDouble>(physTot, 0.0);
1076  coordVel[i] = Array<OneD, NekDouble>(physTot, 0.0);
1077  }
1078 
1079  // Get coordinates velocity in transformed system (for MovingBody regions)
1080  GetCoordVelocity(tmp);
1081  ContravarFromCartesian(tmp, coordVel);
1082 
1083  // Get Cartesian coordinates for evaluating boundary conditions
1085  for (int dir = 0; dir < 3; dir++)
1086  {
1087  coords[dir] = Array<OneD, NekDouble>(physTot, 0.0);
1088  }
1089  GetCartesianCoordinates(coords[0], coords[1], coords[2]);
1090 
1091  // Loop boundary conditions looking for Dirichlet bc's
1092  for (int n = 0; n < nbnds; ++n)
1093  {
1094  // Evaluate original Dirichlet boundary conditions in whole domain
1095  for (int i = 0; i < nfields; ++i)
1096  {
1097  BndConds = m_fields[i]->GetBndConditions();
1098  BndExp = m_fields[i]->GetBndCondExpansions();
1099  if (BndConds[n]->GetBoundaryConditionType() ==
1101  {
1102  isDirichlet[i] = true;
1103  // If we have the a velocity component
1104  // check if all vel bc's are also Dirichlet
1105  if (i < nvel)
1106  {
1107  for (int j = 0; j < nvel; ++j)
1108  {
1109  ASSERTL0(m_fields[j]
1110  ->GetBndConditions()[n]
1111  ->GetBoundaryConditionType() ==
1113  "Mapping only supported when all velocity "
1114  "components have the same type of boundary "
1115  "conditions");
1116  }
1117  }
1118  // Check if bc is time-dependent
1119  ASSERTL0(!BndConds[n]->IsTimeDependent(),
1120  "Time-dependent Dirichlet boundary conditions not "
1121  "supported with mapping yet.");
1122 
1123  // Get boundary condition
1124  LibUtilities::Equation condition =
1125  std::static_pointer_cast<
1127  ->m_dirichletCondition;
1128  // Evaluate
1129  condition.Evaluate(coords[0], coords[1], coords[2], time,
1130  values[i]);
1131  }
1132  else
1133  {
1134  isDirichlet[i] = false;
1135  }
1136  }
1137  // Convert velocity vector to transformed system
1138  if (isDirichlet[0])
1139  {
1140  for (int i = 0; i < nvel; ++i)
1141  {
1142  Vmath::Vcopy(physTot, values[i], 1, tmp[i], 1);
1143  }
1144  ContravarFromCartesian(tmp, tmp2);
1145  for (int i = 0; i < nvel; ++i)
1146  {
1147  Vmath::Vcopy(physTot, tmp2[i], 1, values[i], 1);
1148  }
1149  }
1150 
1151  // Now, project result to boundary
1152  for (int i = 0; i < nfields; ++i)
1153  {
1154  BndConds = m_fields[i]->GetBndConditions();
1155  BndExp = m_fields[i]->GetBndCondExpansions();
1156  if (BndConds[n]->GetUserDefined() == "" ||
1157  BndConds[n]->GetUserDefined() == "MovingBody")
1158  {
1159  m_fields[i]->ExtractPhysToBnd(n, values[i],
1160  BndExp[n]->UpdatePhys());
1161 
1162  // Apply MovingBody correction
1163  if ((i < nvel) && BndConds[n]->GetUserDefined() == "MovingBody")
1164  {
1165  // Get coordinate velocity on boundary
1166  Array<OneD, NekDouble> coordVelBnd(
1167  BndExp[n]->GetTotPoints());
1168  m_fields[i]->ExtractPhysToBnd(n, coordVel[i], coordVelBnd);
1169 
1170  // Apply correction
1171  Vmath::Vadd(BndExp[n]->GetTotPoints(), coordVelBnd, 1,
1172  BndExp[n]->UpdatePhys(), 1,
1173  BndExp[n]->UpdatePhys(), 1);
1174  }
1175  }
1176  }
1177  }
1178 
1179  // Finally, perform FwdTrans in all fields
1180  for (int i = 0; i < m_fields.size(); ++i)
1181  {
1182  // Get boundary condition information
1183  BndConds = m_fields[i]->GetBndConditions();
1184  BndExp = m_fields[i]->GetBndCondExpansions();
1185  for (int n = 0; n < BndConds.size(); ++n)
1186  {
1187  if (BndConds[n]->GetBoundaryConditionType() ==
1189  {
1190  if (BndConds[n]->GetUserDefined() == "" ||
1191  BndConds[n]->GetUserDefined() == "MovingBody")
1192  {
1193  if (m_fields[i]->GetExpType() == MultiRegions::e3DH1D)
1194  {
1195  BndExp[n]->SetWaveSpace(false);
1196  }
1197 
1198  BndExp[n]->FwdTransBndConstrained(
1199  BndExp[n]->GetPhys(), BndExp[n]->UpdateCoeffs());
1200  }
1201  }
1202  }
1203  }
1204 }
1205 
1207  const NekDouble time, const Array<OneD, Array<OneD, NekDouble>> &coords,
1208  const Array<OneD, Array<OneD, NekDouble>> &coordsVel)
1209 {
1210  if (m_fromFunction)
1211  {
1212  std::string s_FieldStr;
1213  string fieldNames[3] = {"x", "y", "z"};
1214  string velFieldNames[3] = {"vx", "vy", "vz"};
1215  // Check if function from session file defines each component
1216  // and evaluate them, otherwise there is no need to update
1217  // coords
1218  for (int i = 0; i < 3; i++)
1219  {
1220  s_FieldStr = fieldNames[i];
1221  if (m_session->DefinesFunction(m_funcName, s_FieldStr))
1222  {
1223  EvaluateFunction(m_fields, m_session, s_FieldStr, m_coords[i],
1224  m_funcName, time);
1225  if (i == 2 && m_fields[0]->GetExpType() == MultiRegions::e3DH1D)
1226  {
1227  ASSERTL0(
1228  false,
1229  "3DH1D does not support mapping in the z-direction.");
1230  }
1231  }
1232  s_FieldStr = velFieldNames[i];
1233  if (m_session->DefinesFunction(m_velFuncName, s_FieldStr))
1234  {
1235  EvaluateFunction(m_fields, m_session, s_FieldStr,
1236  m_coordsVel[i], m_velFuncName, time);
1237  if (i == 2 && m_fields[0]->GetExpType() == MultiRegions::e3DH1D)
1238  {
1239  ASSERTL0(
1240  false,
1241  "3DH1D does not support mapping in the z-direction.");
1242  }
1243  }
1244  }
1245  }
1246  else
1247  {
1248  int physTot = m_fields[0]->GetTotPoints();
1249  int nvel = m_nConvectiveFields;
1250  // Copy coordinates
1251  for (int i = 0; i < 3; i++)
1252  {
1253  Vmath::Vcopy(physTot, coords[i], 1, m_coords[i], 1);
1254  }
1255 
1256  for (int i = 0; i < nvel; i++)
1257  {
1258  Vmath::Vcopy(physTot, coordsVel[i], 1, m_coordsVel[i], 1);
1259  }
1260  }
1261 
1262  // Update the information required by the specific mapping
1263  UpdateGeomInfo();
1264 }
1265 
1266 } // namespace GlobalMapping
1267 } // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
bool m_fromFunction
Flag defining if the Mapping is defined by a function.
Definition: Mapping.h:427
int m_nConvectiveFields
Number of velocity components.
Definition: Mapping.h:414
static MappingSharedPtr m_mappingPtr
Definition: Mapping.h:430
virtual GLOBAL_MAPPING_EXPORT void v_GetCartesianCoordinates(Array< OneD, NekDouble > &out0, Array< OneD, NekDouble > &out1, Array< OneD, NekDouble > &out2)
Definition: Mapping.cpp:654
GLOBAL_MAPPING_EXPORT void EvaluateTimeFunction(LibUtilities::SessionReaderSharedPtr pSession, std::string pFieldName, Array< OneD, NekDouble > &pArray, const std::string &pFunctionName, NekDouble pTime=NekDouble(0))
Definition: Mapping.cpp:383
LibUtilities::SessionReaderSharedPtr m_session
Session reader.
Definition: Mapping.h:402
virtual GLOBAL_MAPPING_EXPORT void v_UpdateBCs(const NekDouble time)
Definition: Mapping.cpp:1050
GLOBAL_MAPPING_EXPORT void ApplyChristoffelCovar(const Array< OneD, Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
Apply the Christoffel symbols to a covariant vector.
Definition: Mapping.h:229
GLOBAL_MAPPING_EXPORT void GetCartesianCoordinates(Array< OneD, NekDouble > &out0, Array< OneD, NekDouble > &out1, Array< OneD, NekDouble > &out2)
Get the Cartesian coordinates in the field.
Definition: Mapping.h:135
GLOBAL_MAPPING_EXPORT bool IsTimeDependent()
Get flag defining if mapping is time-dependent.
Definition: Mapping.h:338
virtual GLOBAL_MAPPING_EXPORT void v_VelocityLaplacian(const Array< OneD, Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble alpha)
Definition: Mapping.cpp:779
Array< OneD, Array< OneD, NekDouble > > m_wk2
Definition: Mapping.h:436
GLOBAL_MAPPING_EXPORT void ContravarToCartesian(const Array< OneD, Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
Convert a contravariant vector to the Cartesian system.
Definition: Mapping.cpp:480
GLOBAL_MAPPING_EXPORT void UpdateGeomInfo()
Recompute the metric terms of the Mapping.
Definition: Mapping.h:395
Array< OneD, Array< OneD, NekDouble > > m_coords
Array with the Cartesian coordinates.
Definition: Mapping.h:408
Array< OneD, Array< OneD, NekDouble > > m_coordsVel
Array with the velocity of the coordinates.
Definition: Mapping.h:410
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Definition: Mapping.h:406
std::string m_velFuncName
Name of the function containing the velocity of the coordinates.
Definition: Mapping.h:419
GLOBAL_MAPPING_EXPORT void InitObject(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const TiXmlElement *pMapping)
Initialise the mapping object.
Definition: Mapping.h:77
virtual GLOBAL_MAPPING_EXPORT void v_CovarToCartesian(const Array< OneD, Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray)=0
GLOBAL_MAPPING_EXPORT void Output(LibUtilities::FieldMetaDataMap &fieldMetaDataMap, const std::string &outname)
Output function called when a chk or fld file is written.
Definition: Mapping.cpp:305
GLOBAL_MAPPING_EXPORT void ContravarFromCartesian(const Array< OneD, Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
Convert a contravariant vector to the transformed system.
Definition: Mapping.cpp:542
GLOBAL_MAPPING_EXPORT void GetCoordVelocity(Array< OneD, Array< OneD, NekDouble >> &outarray)
Obtain the velocity of the coordinates.
Definition: Mapping.h:244
virtual GLOBAL_MAPPING_EXPORT void v_LowerIndex(const Array< OneD, Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
Definition: Mapping.cpp:709
Array< OneD, Array< OneD, NekDouble > > m_wk1
Definition: Mapping.h:435
GLOBAL_MAPPING_EXPORT void EvaluateFunction(Array< OneD, MultiRegions::ExpListSharedPtr > pFields, LibUtilities::SessionReaderSharedPtr pSession, std::string pFieldName, Array< OneD, NekDouble > &pArray, const std::string &pFunctionName, NekDouble pTime=NekDouble(0))
Definition: Mapping.cpp:401
virtual GLOBAL_MAPPING_EXPORT void v_CurlCurlField(Array< OneD, Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const bool generalized)
Definition: Mapping.cpp:988
virtual GLOBAL_MAPPING_EXPORT void v_CovarFromCartesian(const Array< OneD, Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray)=0
LibUtilities::FieldIOSharedPtr m_fld
Definition: Mapping.h:404
GLOBAL_MAPPING_EXPORT void GetMetricTensor(Array< OneD, Array< OneD, NekDouble >> &outarray)
Get the metric tensor .
Definition: Mapping.h:177
GLOBAL_MAPPING_EXPORT void CovarFromCartesian(const Array< OneD, Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
Convert a covariant vector to the transformed system.
Definition: Mapping.cpp:573
virtual GLOBAL_MAPPING_EXPORT void v_gradgradU(const Array< OneD, Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
Definition: Mapping.cpp:880
GLOBAL_MAPPING_EXPORT bool HasConstantJacobian()
Get flag defining if mapping has constant Jacobian.
Definition: Mapping.h:362
GLOBAL_MAPPING_EXPORT void gradgradU(const Array< OneD, Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
Second order covariant derivatives of a contravariant vector.
Definition: Mapping.h:306
GLOBAL_MAPPING_EXPORT void ReplaceField(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
Replace the Expansion List used by the mapping.
Definition: Mapping.cpp:246
GLOBAL_MAPPING_EXPORT void GetInvMetricTensor(Array< OneD, Array< OneD, NekDouble >> &outarray)
Get the inverse of metric tensor .
Definition: Mapping.h:184
virtual GLOBAL_MAPPING_EXPORT void v_DotGradJacobian(const Array< OneD, Array< OneD, NekDouble >> &inarray, Array< OneD, NekDouble > &outarray)
Definition: Mapping.cpp:680
virtual GLOBAL_MAPPING_EXPORT void v_Divergence(const Array< OneD, Array< OneD, NekDouble >> &inarray, Array< OneD, NekDouble > &outarray)
Definition: Mapping.cpp:751
virtual GLOBAL_MAPPING_EXPORT void v_ContravarFromCartesian(const Array< OneD, Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray)=0
GLOBAL_MAPPING_EXPORT void LowerIndex(const Array< OneD, Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
Lower index of vector: .
Definition: Mapping.cpp:604
GLOBAL_MAPPING_EXPORT void CovarToCartesian(const Array< OneD, Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
Convert a covariant vector to the Cartesian system.
Definition: Mapping.cpp:511
std::string m_funcName
Name of the function containing the coordinates.
Definition: Mapping.h:417
static GLOBAL_MAPPING_EXPORT MappingSharedPtr Load(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
Return a pointer to the mapping, creating it on first call.
Definition: Mapping.cpp:269
GLOBAL_MAPPING_EXPORT void GetJacobian(Array< OneD, NekDouble > &outarray)
Get the Jacobian of the transformation.
Definition: Mapping.h:155
virtual GLOBAL_MAPPING_EXPORT void v_ContravarToCartesian(const Array< OneD, Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray)=0
virtual GLOBAL_MAPPING_EXPORT void v_InitObject(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const TiXmlElement *pMapping)
Definition: Mapping.cpp:100
Array< OneD, Array< OneD, NekDouble > > m_tmp
Definition: Mapping.h:437
virtual GLOBAL_MAPPING_EXPORT void v_UpdateMapping(const NekDouble time, const Array< OneD, Array< OneD, NekDouble >> &coords=NullNekDoubleArrayOfArray, const Array< OneD, Array< OneD, NekDouble >> &coordsVel=NullNekDoubleArrayOfArray)
Definition: Mapping.cpp:1206
virtual GLOBAL_MAPPING_EXPORT void v_RaiseIndex(const Array< OneD, Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
Definition: Mapping.cpp:730
GLOBAL_MAPPING_EXPORT void ApplyChristoffelContravar(const Array< OneD, Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
Apply the Christoffel symbols to a contravariant vector.
Definition: Mapping.h:211
bool m_timeDependent
Flag defining if the Mapping is time-dependent.
Definition: Mapping.h:425
GLOBAL_MAPPING_EXPORT void RaiseIndex(const Array< OneD, Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
Raise index of vector: .
Definition: Mapping.cpp:634
virtual GLOBAL_MAPPING_EXPORT void v_GetCoordVelocity(Array< OneD, Array< OneD, NekDouble >> &outarray)
Definition: Mapping.cpp:669
NekDouble Evaluate() const
Definition: Equation.cpp:95
static std::shared_ptr< FieldIO > CreateForFile(const LibUtilities::SessionReaderSharedPtr session, const std::string &filename)
Construct a FieldIO object for a given input filename.
Definition: FieldIO.cpp:227
static std::shared_ptr< FieldIO > CreateDefault(const LibUtilities::SessionReaderSharedPtr session)
Returns an object for the default FieldIO method.
Definition: FieldIO.cpp:198
Provides a generic Factory class.
Definition: NekFactory.hpp:105
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:144
GLOBAL_MAPPING_EXPORT typedef std::shared_ptr< Mapping > MappingSharedPtr
A shared pointer to a Mapping object.
Definition: Mapping.h:50
MappingFactory & GetMappingFactory()
Declaration of the mapping factory singleton.
Definition: Mapping.cpp:52
std::shared_ptr< FieldIO > FieldIOSharedPtr
Definition: FieldIO.h:301
std::map< std::string, std::string > FieldMetaDataMap
Definition: FieldIO.h:52
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< Equation > EquationSharedPtr
Definition: Equation.h:130
MultiRegions::Direction const DirCartesianMap[]
Definition: ExpList.h:89
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
double NekDouble
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:209
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:574
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:359
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 Vdiv(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x/y.
Definition: Vmath.cpp:284
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:492
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1255
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:419