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