Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator 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  (pSession->GetComm());
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  int phystot = m_fields[0]->GetTotPoints();
105  m_fromFunction = true;
106  // Initialise variables
110  for (int i = 0; i < 3; i++)
111  {
112  m_coords[i] = Array<OneD, NekDouble> (phystot);
113  m_coordsVel[i] = Array<OneD, NekDouble> (phystot);
114  coords[i] = Array<OneD, NekDouble> (phystot);
115  }
116 
117  // Check if mapping is defined as time-dependent
118  const TiXmlElement* timeDep = pMapping->
119  FirstChildElement("TIMEDEPENDENT");
120  if (timeDep)
121  {
122  string sTimeDep = timeDep->GetText();
123  m_timeDependent = ( boost::iequals(sTimeDep,"true")) ||
124  ( boost::iequals(sTimeDep,"yes"));
125  }
126  else
127  {
128  m_timeDependent = false;
129  }
130 
131  // Load coordinates
132  string fieldNames[3] = {"x", "y", "z"};
133  const TiXmlElement* funcNameElmt = pMapping->FirstChildElement("COORDS");
134  if (funcNameElmt)
135  {
136  m_funcName = funcNameElmt->GetText();
137  ASSERTL0(m_session->DefinesFunction(m_funcName),
138  "Function '" + m_funcName + "' not defined.");
139 
140  // Get coordinates in the domain
141  m_fields[0]->GetCoords(coords[0], coords[1], coords[2]);
142 
143  std::string s_FieldStr;
144  // Check if function from session file defines each component
145  // and evaluate them, otherwise use trivial transformation
146  for(int i = 0; i < 3; i++)
147  {
148  s_FieldStr = fieldNames[i];
149  if ( m_session->DefinesFunction(m_funcName, s_FieldStr))
150  {
151  EvaluateFunction(m_fields, m_session, s_FieldStr, m_coords[i],
152  m_funcName);
153  if ( i==2 && m_fields[0]->GetExpType() == MultiRegions::e3DH1D)
154  {
155  ASSERTL0 (false,
156  "3DH1D does not support mapping in the z-direction.");
157  }
158  }
159  else
160  {
161  // This coordinate is not defined, so use (x^i)' = x^i
162  Vmath::Vcopy(phystot, coords[i], 1, m_coords[i], 1);
163  }
164  }
165  }
166  else
167  {
168  m_fields[0]->GetCoords(coords[0], coords[1], coords[2]);
169  for(int i = 0; i < 3; i++)
170  {
171  // Use (x^i)' = x^i as default. This can be useful if we
172  // have a time-dependent mapping, and then only the
173  // initial mapping will be trivial
174  Vmath::Vcopy(phystot, coords[i], 1, m_coords[i], 1);
175  }
176  }
177 
178  // Load coordinate velocity if they are defined,
179  // otherwise use zero to make it general
180  string velFieldNames[3] = {"vx", "vy", "vz"};
181  const TiXmlElement* velFuncNameElmt = pMapping->FirstChildElement("VEL");
182  if (velFuncNameElmt)
183  {
184  m_velFuncName = velFuncNameElmt->GetText();
185  ASSERTL0(m_session->DefinesFunction(m_velFuncName),
186  "Function '" + m_velFuncName + "' not defined.");
187 
188  std::string s_FieldStr;
189  // Check if function from session file defines each component
190  // and evaluate them, otherwise use 0
191  for(int i = 0; i < 3; i++)
192  {
193  s_FieldStr = velFieldNames[i];
194  if ( m_session->DefinesFunction(m_velFuncName, s_FieldStr))
195  {
196  EvaluateFunction(m_fields, m_session, s_FieldStr,
197  m_coordsVel[i], m_velFuncName);
198  if ( i==2 && m_fields[0]->GetExpType() == MultiRegions::e3DH1D)
199  {
200  ASSERTL0 (false,
201  "3DH1D does not support mapping in the z-direction.");
202  }
203  }
204  else
205  {
206  // This coordinate velocity is not defined, so use 0
207  Vmath::Zero(phystot, m_coordsVel[i], 1);
208  }
209  }
210  }
211  else
212  {
213  for(int i = 0; i < 3; i++)
214  {
215  Vmath::Zero(phystot, m_coordsVel[i], 1);
216  }
217  }
218 
219  // Initialise workspace variables
220  int nvel = m_nConvectiveFields;
224  for (int i=0; i< nvel; i++)
225  {
226  m_tmp[i] = Array<OneD, NekDouble>(phystot,0.0);
227  for (int j=0; j< nvel; j++)
228  {
229  m_wk1[i*nvel+j] = Array<OneD, NekDouble>(phystot,0.0);
230  m_wk2[i*nvel+j] = Array<OneD, NekDouble>(phystot,0.0);
231  }
232  }
233 
234  // Calculate information required by the particular mapping
235  UpdateGeomInfo();
236 }
237 
238 /**
239  * This function replaces the expansion list contained in m_fields, and then
240  * proceeds reinitialising the Mapping with this new field.
241  * @param pFields New field to be used by the Mapping
242  */
245 {
246  m_fields = pFields;
247 
248  TiXmlElement* vMapping = NULL;
249 
250  if (m_session->DefinesElement("Nektar/Mapping"))
251  {
252  vMapping = m_session->GetElement("Nektar/Mapping");
253  }
254  InitObject(pFields, vMapping);
255 }
256 
257 /**
258  * This function is responsible for loading the Mapping, guaranteeing that a
259  * single instance of this class exists. When it is first called, it creates a
260  * Mapping and returns a pointer to it. On subsequent calls, it just returns
261  * the pointer.
262  * @param pSession Session reader object
263  * @param pFields Fields which will be used by the Mapping
264  * @return Pointer to the Mapping
265  */
267  const LibUtilities::SessionReaderSharedPtr& pSession,
269 {
270  if (!m_init)
271  {
272  TiXmlElement* vMapping = NULL;
273  string vType;
274  if (pSession->DefinesElement("Nektar/Mapping"))
275  {
276  vMapping = pSession->GetElement("Nektar/Mapping");
277  vType = vMapping->Attribute("TYPE");
278  m_isDefined = true;
279  }
280  else
281  {
282  vType = "Translation";
283  }
284 
286  vType, pSession, pFields,
287  vMapping);
288 
289  m_init = true;
290  }
291 
292  return m_mappingPtr;
293 }
294 
295 /**
296  * This function should be called before writing a fld or chk file. It updates
297  * the metadata with information from the Mapping. Also, if the mapping is not
298  * defined by a function, it writes a .map file containing the coordinates and
299  * velocity of the coordinates.
300  * @param fieldMetaDataMap Metadata of the output file
301  * @param outname Name of the output file
302  */
304  LibUtilities::FieldMetaDataMap &fieldMetaDataMap,
305  const std::string &outname)
306 {
307  // Only do anything if mapping exists
308  if (m_isDefined)
309  {
310  fieldMetaDataMap["MappingCartesianVel"] = std::string("False");
311  if (m_fromFunction)
312  {
313  // Add metadata
314  fieldMetaDataMap["MappingType"] = std::string("Expression");
315  fieldMetaDataMap["MappingExpression"] = m_funcName;
316  if (m_timeDependent)
317  {
318  fieldMetaDataMap["MappingVelExpression"] = m_velFuncName;
319  }
320  }
321  else
322  {
323  int expdim = m_fields[0]->GetGraph()->GetMeshDimension();
324  string fieldNames[3] = {"x", "y", "z"};
325  string velFieldNames[3] = {"vx", "vy", "vz"};
326 
327  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
328  = m_fields[0]->GetFieldDefinitions();
329  std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
330 
331  int ncoeffs = m_fields[0]->GetNcoeffs();
332  Array<OneD, NekDouble> fieldcoeffs(ncoeffs);
333 
334  bool wavespace = m_fields[0]->GetWaveSpace();
335  m_fields[0]->SetWaveSpace(false);
336  // copy coordinates Data into FieldData and set variable
337  for(int j = 0; j < expdim; ++j)
338  {
339  m_fields[0]->FwdTrans_IterPerExp(m_coords[j], fieldcoeffs);
340 
341  for(int i = 0; i < FieldDef.size(); ++i)
342  {
343  // Could do a search here to find correct variable
344  FieldDef[i]->m_fields.push_back(fieldNames[j]);
345  m_fields[0]->AppendFieldData(FieldDef[i], FieldData[i],
346  fieldcoeffs);
347  }
348  }
349  if (m_timeDependent)
350  {
351  //copy coordinates velocity Data into FieldData and set variable
352  for(int j = 0; j < expdim; ++j)
353  {
354  m_fields[0]->FwdTrans_IterPerExp(m_coordsVel[j],
355  fieldcoeffs);
356 
357  for(int i = 0; i < FieldDef.size(); ++i)
358  {
359  // Could do a search here to find correct variable
360  FieldDef[i]->m_fields.push_back(velFieldNames[j]);
361  m_fields[0]->AppendFieldData(FieldDef[i],
362  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 
385  std::string pFieldName,
386  Array<OneD, NekDouble>& pArray,
387  const std::string& pFunctionName,
388  NekDouble pTime)
389 {
390  ASSERTL0(pSession->DefinesFunction(pFunctionName),
391  "Function '" + pFunctionName + "' does not exist.");
392 
394  pSession->GetFunction(pFunctionName, pFieldName);
395 
396  Array<OneD, NekDouble> x0(1,0.0);
397  Array<OneD, NekDouble> x1(1,0.0);
398  Array<OneD, NekDouble> x2(1,0.0);
399 
400  ffunc->Evaluate(x0, x1, x2, pTime, pArray);
401 }
402 
403 
407  std::string pFieldName,
408  Array<OneD, NekDouble>& pArray,
409  const std::string& pFunctionName,
410  NekDouble pTime)
411 {
412  ASSERTL0(pSession->DefinesFunction(pFunctionName),
413  "Function '" + pFunctionName + "' does not exist.");
414 
415  unsigned int nq = pFields[0]->GetNpoints();
416  if (pArray.num_elements() != nq)
417  {
418  pArray = Array<OneD, NekDouble> (nq);
419  }
420 
422  vType = pSession->GetFunctionType(pFunctionName, pFieldName);
424  {
425  Array<OneD, NekDouble> x0(nq);
426  Array<OneD, NekDouble> x1(nq);
427  Array<OneD, NekDouble> x2(nq);
428 
429  pFields[0]->GetCoords(x0, x1, x2);
431  pSession->GetFunction(pFunctionName, pFieldName);
432 
433  ffunc->Evaluate(x0, x1, x2, pTime, pArray);
434  }
435  else if (vType == LibUtilities::eFunctionTypeFile)
436  {
437  std::string filename = pSession->GetFunctionFilename(
438  pFunctionName,
439  pFieldName);
440 
441  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
442  std::vector<std::vector<NekDouble> > FieldData;
443  Array<OneD, NekDouble> vCoeffs(pFields[0]->GetNcoeffs());
444  Vmath::Zero(vCoeffs.num_elements(), vCoeffs, 1);
445 
448  (m_session->GetComm());
449  fld->Import(filename, FieldDef, FieldData);
450 
451  int idx = -1;
452  for (int i = 0; i < FieldDef.size(); ++i)
453  {
454  for (int j = 0; j < FieldDef[i]->m_fields.size(); ++j)
455  {
456  if (FieldDef[i]->m_fields[j] == pFieldName)
457  {
458  idx = j;
459  }
460  }
461 
462  if (idx >= 0)
463  {
464  pFields[0]->ExtractDataToCoeffs(
465  FieldDef[i],
466  FieldData[i],
467  FieldDef[i]->m_fields[idx],
468  vCoeffs);
469  }
470  else
471  {
472  cout << "Field " + pFieldName + " not found." << endl;
473  }
474  }
475  pFields[0]->BwdTrans_IterPerExp(vCoeffs, pArray);
476  }
477 }
478 
479 /**
480  * This function converts a contravariant vector in transformed space
481  * \f$v^{j}\f$ to the corresponding \f$\bar{v}^{i}\f$ in Cartesian (physical)
482  * space using the relation
483  * \f[\bar{v}^{i} = \frac{\partial \bar{x}^i}{\partial x^j}v^{j}\f]
484  *
485  * @param inarray Components of the vector in transformed space (\f$v^{j}\f$)
486  * @param outarray Components of the vector in Cartesian space (\f$\bar{v}^{i}\f$)
487  */
489  const Array<OneD, Array<OneD, NekDouble> > &inarray,
490  Array<OneD, Array<OneD, NekDouble> > &outarray)
491 {
492  if(inarray == outarray)
493  {
494  int physTot = m_fields[0]->GetTotPoints();
495  int nvel = m_nConvectiveFields;
496 
497  for (int i=0; i< nvel; i++)
498  {
499  Vmath::Vcopy(physTot, inarray[i], 1, m_tmp[i], 1);
500  }
501  v_ContravarToCartesian( m_tmp, outarray);
502  }
503  else
504  {
505  v_ContravarToCartesian( inarray, outarray);
506  }
507 }
508 
509 /**
510  * This function converts a covariant vector in transformed space
511  * \f$v_{j}\f$ to the corresponding \f$\bar{v}_{i}\f$ in Cartesian (physical)
512  * space using the relation
513  * \f[\bar{v}_{i} = \frac{\partial x^j}{\partial \bar{x}^i}v_{j}\f]
514  *
515  * @param inarray Components of the vector in transformed space (\f$v_{j}\f$)
516  * @param outarray Components of the vector in Cartesian space (\f$\bar{v}_{i}\f$)
517  */
519  const Array<OneD, Array<OneD, NekDouble> > &inarray,
520  Array<OneD, Array<OneD, NekDouble> > &outarray)
521 {
522  if(inarray == outarray)
523  {
524  int physTot = m_fields[0]->GetTotPoints();
525  int nvel = m_nConvectiveFields;
526 
527  for (int i=0; i< nvel; i++)
528  {
529  Vmath::Vcopy(physTot, inarray[i], 1, m_tmp[i], 1);
530  }
531  v_CovarToCartesian( m_tmp, outarray);
532  }
533  else
534  {
535  v_CovarToCartesian( inarray, outarray);
536  }
537 }
538 
539 /**
540  * This function converts a contravariant vector in Cartesian space
541  * \f$\bar{v}^{i}\f$ to the corresponding \f$v^{j}\f$ in transformed
542  * space using the relation
543  * \f[v^{j} = \frac{\partial x^j}{\partial \bar{x}^i}\bar{v}^{i}\f]
544  *
545  * @param inarray Components of the vector in Cartesian space (\f$\bar{v}^{i}\f$)
546  * @param outarray Components of the vector in transformed space (\f$v^{j}\f$)
547  */
549  const Array<OneD, Array<OneD, NekDouble> > &inarray,
550  Array<OneD, Array<OneD, NekDouble> > &outarray)
551 {
552  if(inarray == outarray)
553  {
554  int physTot = m_fields[0]->GetTotPoints();
555  int nvel = m_nConvectiveFields;
556 
557  for (int i=0; i< nvel; i++)
558  {
559  Vmath::Vcopy(physTot, inarray[i], 1, m_tmp[i], 1);
560  }
561  v_ContravarFromCartesian( m_tmp, outarray);
562  }
563  else
564  {
565  v_ContravarFromCartesian( inarray, outarray);
566  }
567 }
568 
569 /**
570  * This function converts a covariant vector in Cartesian space
571  * \f$\bar{v}_{i}\f$ to the corresponding \f$v_{j}\f$ in transformed
572  * space using the relation
573  * \f[\bar{v}_{j} = \frac{\partial \bar{x}^i}{\partial x^j}v_{i}\f]
574  *
575  * @param inarray Components of the vector in Cartesian space (\f$\bar{v}_{i}\f$)
576  * @param outarray Components of the vector in transformed space (\f$v_{j}\f$)
577  */
579  const Array<OneD, Array<OneD, NekDouble> > &inarray,
580  Array<OneD, Array<OneD, NekDouble> > &outarray)
581 {
582  if(inarray == outarray)
583  {
584  int physTot = m_fields[0]->GetTotPoints();
585  int nvel = m_nConvectiveFields;
586 
587  for (int i=0; i< nvel; i++)
588  {
589  Vmath::Vcopy(physTot, inarray[i], 1, m_tmp[i], 1);
590  }
591  v_CovarFromCartesian( m_tmp, outarray);
592  }
593  else
594  {
595  v_CovarFromCartesian( inarray, outarray);
596  }
597 }
598 
599 /**
600  * This function lowers the index of the contravariant vector \f$v^{i}\f$,
601  * transforming it into its associated covariant vector \f$v_{j}\f$
602  * according to the relation
603  * \f[v_{j} = g_{ij}v^{i}\f]
604  * where \f$g_{ij}\f$ is the metric tensor.
605  *
606  * @param inarray Components of the contravariant vector \f$v^{i}\f$
607  * @param outarray Components of the covariant vector \f$v_{j}\f$
608  */
610  const Array<OneD, Array<OneD, NekDouble> > &inarray,
611  Array<OneD, Array<OneD, NekDouble> > &outarray)
612 {
613  if(inarray == outarray)
614  {
615  int physTot = m_fields[0]->GetTotPoints();
616  int nvel = m_nConvectiveFields;
617 
618  for (int i=0; i< nvel; i++)
619  {
620  Vmath::Vcopy(physTot, inarray[i], 1, m_tmp[i], 1);
621  }
622  v_LowerIndex( m_tmp, outarray);
623  }
624  else
625  {
626  v_LowerIndex( inarray, outarray);
627  }
628 }
629 
630 /**
631  * This function raises the index of the covariant vector \f$v_{j}\f$,
632  * transforming it into its associated contravariant vector \f$v^{i}\f$
633  * according to the relation
634  * \f[v^{i} = g^{ij}v_{j}\f]
635  * where \f$g^{ij}\f$ is the inverse of the metric tensor.
636  *
637  * @param inarray Components of the contravariant vector \f$v^{i}\f$
638  * @param outarray Components of the covariant vector \f$v_{j}\f$
639  */
641  const Array<OneD, Array<OneD, NekDouble> > &inarray,
642  Array<OneD, Array<OneD, NekDouble> > &outarray)
643 {
644  if(inarray == outarray)
645  {
646  int physTot = m_fields[0]->GetTotPoints();
647  int nvel = m_nConvectiveFields;
648 
649  for (int i=0; i< nvel; i++)
650  {
651  Vmath::Vcopy(physTot, inarray[i], 1, m_tmp[i], 1);
652  }
653  v_RaiseIndex( m_tmp, outarray);
654  }
655  else
656  {
657  v_RaiseIndex( inarray, outarray);
658  }
659 }
660 
665 {
666  int physTot = m_fields[0]->GetTotPoints();
667 
668  out0 = Array<OneD, NekDouble>(physTot, 0.0);
669  out1 = Array<OneD, NekDouble>(physTot, 0.0);
670  out2 = Array<OneD, NekDouble>(physTot, 0.0);
671 
672  Vmath::Vcopy(physTot, m_coords[0], 1, out0, 1);
673  Vmath::Vcopy(physTot, m_coords[1], 1, out1, 1);
674  Vmath::Vcopy(physTot, m_coords[2], 1, out2, 1);
675 }
676 
678  Array<OneD, Array<OneD, NekDouble> > &outarray)
679 {
680  int physTot = m_fields[0]->GetTotPoints();
681 
682  for(int i = 0; i < m_nConvectiveFields; ++i)
683  {
684  outarray[i] = Array<OneD, NekDouble>(physTot, 0.0);
685  Vmath::Vcopy(physTot, m_coordsVel[i], 1, outarray[i], 1);
686  }
687 }
688 
690  const Array<OneD, Array<OneD, NekDouble> > &inarray,
691  Array<OneD, NekDouble> &outarray)
692 {
693  int physTot = m_fields[0]->GetTotPoints();
694 
695  outarray = Array<OneD, NekDouble>(physTot, 0.0);
696  if ( !HasConstantJacobian() )
697  {
698  // Set wavespace to false and store current value
699  bool wavespace = m_fields[0]->GetWaveSpace();
700  m_fields[0]->SetWaveSpace(false);
701 
702  // Get Mapping Jacobian
703  Array<OneD, NekDouble> Jac(physTot, 0.0);
704  GetJacobian(Jac);
705 
706  // Calculate inarray . grad(Jac)
707  Array<OneD, NekDouble> wk(physTot, 0.0);
708  for(int i = 0; i < m_nConvectiveFields; ++i)
709  {
710  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[i],
711  Jac, wk);
712  Vmath::Vvtvp(physTot, inarray[i], 1, wk, 1,
713  outarray, 1, outarray, 1);
714  }
715  m_fields[0]->SetWaveSpace(wavespace);
716  }
717 }
718 
720  const Array<OneD, Array<OneD, NekDouble> > &inarray,
721  Array<OneD, Array<OneD, NekDouble> > &outarray)
722 {
723  int physTot = m_fields[0]->GetTotPoints();
724  int nvel = m_nConvectiveFields;
725 
726  Array<OneD, Array<OneD, NekDouble> > g(nvel*nvel);
727 
728  GetMetricTensor(g);
729 
730  for (int i=0; i< nvel; i++)
731  {
732  outarray[i] = Array<OneD, NekDouble> (physTot, 0.0);
733  for (int j=0; j< nvel; j++)
734  {
735  Vmath::Vvtvp(physTot, g[i*nvel+j], 1, inarray[j], 1,
736  outarray[i], 1,
737  outarray[i], 1);
738  }
739  }
740 }
741 
743  const Array<OneD, Array<OneD, NekDouble> > &inarray,
744  Array<OneD, Array<OneD, NekDouble> > &outarray)
745 {
746  int physTot = m_fields[0]->GetTotPoints();
747  int nvel = m_nConvectiveFields;
748 
749  Array<OneD, Array<OneD, NekDouble> > g(nvel*nvel);
750 
752 
753  for (int i=0; i< nvel; i++)
754  {
755  outarray[i] = Array<OneD, NekDouble> (physTot, 0.0);
756  for (int j=0; j< nvel; j++)
757  {
758  Vmath::Vvtvp(physTot, g[i*nvel+j], 1, inarray[j], 1,
759  outarray[i], 1,
760  outarray[i], 1);
761  }
762  }
763 }
764 
766  const Array<OneD, Array<OneD, NekDouble> > &inarray,
767  Array<OneD, NekDouble> &outarray)
768 {
769  int physTot = m_fields[0]->GetTotPoints();
770  Array<OneD, NekDouble> wk(physTot, 0.0);
771 
772  Vmath::Zero(physTot, outarray, 1);
773 
774  // Set wavespace to false and store current value
775  bool wavespace = m_fields[0]->GetWaveSpace();
776  m_fields[0]->SetWaveSpace(false);
777 
778  // Get Mapping Jacobian
779  Array<OneD, NekDouble> Jac(physTot, 0.0);
780  GetJacobian(Jac);
781 
782  for(int i = 0; i < m_nConvectiveFields; ++i)
783  {
784  Vmath::Vmul(physTot,Jac, 1, inarray[i], 1, wk, 1); // J*Ui
785  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[i],
786  wk, wk); // (J*Ui)_i
787  Vmath::Vadd(physTot, wk, 1, outarray, 1, outarray, 1);
788  }
789  Vmath::Vdiv(physTot,outarray,1,Jac,1,outarray,1); //1/J*(J*Ui)_i
790 
791  m_fields[0]->SetWaveSpace(wavespace);
792 }
793 
795  const Array<OneD, Array<OneD, NekDouble> > &inarray,
796  Array<OneD, Array<OneD, NekDouble> > &outarray,
797  const NekDouble alpha)
798 {
799  int physTot = m_fields[0]->GetTotPoints();
800  int nvel = m_nConvectiveFields;
801 
803 
804  // Set wavespace to false and store current value
805  bool wavespace = m_fields[0]->GetWaveSpace();
806  m_fields[0]->SetWaveSpace(false);
807 
808  // Calculate vector gradient wk2 = u^i_(,k) = du^i/dx^k + {i,jk}*u^j
809  ApplyChristoffelContravar(inarray, m_wk1);
810  for (int i=0; i< nvel; i++)
811  {
812  if(nvel == 2)
813  {
814  m_fields[0]->PhysDeriv(inarray[i],
815  m_wk2[i*nvel+0],
816  m_wk2[i*nvel+1]);
817  }
818  else
819  {
820  m_fields[0]->PhysDeriv(inarray[i],
821  m_wk2[i*nvel+0],
822  m_wk2[i*nvel+1],
823  m_wk2[i*nvel+2]);
824  }
825  for (int k=0; k< nvel; k++)
826  {
827  Vmath::Vadd(physTot,m_wk1[i*nvel+k],1,m_wk2[i*nvel+k],1,
828  m_wk1[i*nvel+k], 1);
829  }
830  }
831  // Calculate wk1 = A^(ij) = g^(jk)*u^i_(,k)
832  for (int i=0; i< nvel; i++)
833  {
834  for (int k=0; k< nvel; k++)
835  {
836  tmp[k] = m_wk1[i*nvel+k];
837  }
838  RaiseIndex(tmp, m_tmp);
839  for (int j=0; j<nvel; j++)
840  {
841  Vmath::Vcopy(physTot, m_tmp[j], 1, m_wk1[i*nvel+j], 1);
842  }
843  }
844  //
845  // Calculate L(U)^i = (A^(ij))_(,j) - alpha*d^2(u^i)/dx^jdx^j
846  //
847 
848  // Step 1 :
849  // d(A^(ij) - alpha*du^i/dx^j)/d(x^j)
850  for (int i=0; i< nvel; i++)
851  {
852  outarray[i] = Array<OneD, NekDouble>(physTot,0.0);
853  for (int j=0; j< nvel; j++)
854  {
855  Vmath::Smul(physTot, alpha, m_wk2[i*nvel+j], 1,
856  m_tmp[0], 1);
857  Vmath::Vsub(physTot, m_wk1[i*nvel+j], 1, m_tmp[0], 1,
858  m_tmp[0], 1);
859 
860  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[j],
861  m_tmp[0],
862  m_tmp[1]);
863  Vmath::Vadd(physTot,outarray[i],1,m_tmp[1],1,outarray[i], 1);
864  }
865  }
866 
867  // Step 2: d(A^(ij))/d(x^j) + {j,pj}*A^(ip)
868  for (int i=0; i< nvel; i++)
869  {
870  for (int p=0; p< nvel; p++)
871  {
872  tmp[p] = m_wk1[i*nvel+p];
873  }
875  for (int j=0; j< nvel; j++)
876  {
877  Vmath::Vadd(physTot,outarray[i],1,m_wk2[j*nvel+j],1,
878  outarray[i], 1);
879  }
880  }
881 
882  // Step 3: d(A^(ij))/d(x^j) + {j,pj}*A^(ip) + {i,pj} A^(pj)
883  for (int j=0; j< nvel; j++)
884  {
885  for (int p=0; p< nvel; p++)
886  {
887  tmp[p] = m_wk1[p*nvel+j];
888  }
890  for (int i=0; i< nvel; i++)
891  {
892  Vmath::Vadd(physTot,outarray[i], 1, m_wk2[i*nvel+j], 1,
893  outarray[i], 1);
894  }
895  }
896 
897  // Restore value of wavespace
898  m_fields[0]->SetWaveSpace(wavespace);
899 }
900 
902  const Array<OneD, Array<OneD, NekDouble> > &inarray,
903  Array<OneD, Array<OneD, NekDouble> > &outarray)
904 {
905  int physTot = m_fields[0]->GetTotPoints();
906  int nvel = m_nConvectiveFields;
907 
908  // Declare variables
909  outarray = Array<OneD, Array<OneD, NekDouble> > (nvel*nvel*nvel);
911  for (int i=0; i< nvel*nvel*nvel; i++)
912  {
913  outarray[i] = Array<OneD, NekDouble>(physTot,0.0);
914  }
915 
916  // Set wavespace to false and store current value
917  bool wavespace = m_fields[0]->GetWaveSpace();
918  m_fields[0]->SetWaveSpace(false);
919 
920  // Calculate vector gradient u^i_(,j) = du^i/dx^j + {i,pj}*u^p
921  ApplyChristoffelContravar(inarray, m_wk1);
922  for (int i=0; i< nvel; i++)
923  {
924  if (nvel == 2)
925  {
926  m_fields[0]->PhysDeriv(inarray[i],
927  m_wk2[i*nvel+0],
928  m_wk2[i*nvel+1]);
929  }
930  else
931  {
932  m_fields[0]->PhysDeriv(inarray[i],
933  m_wk2[i*nvel+0],
934  m_wk2[i*nvel+1],
935  m_wk2[i*nvel+2]);
936  }
937  for (int j=0; j< nvel; j++)
938  {
939  Vmath::Vadd(physTot,m_wk1[i*nvel+j],1,m_wk2[i*nvel+j],1,
940  m_wk1[i*nvel+j], 1);
941  }
942  }
943 
944  //
945  // Calculate (u^i_,j),k
946  //
947 
948  // Step 1 : d(u^i_,j))/d(x^k)
949  for (int i=0; i< nvel; i++)
950  {
951  for (int j=0; j< nvel; j++)
952  {
953  if (nvel == 2)
954  {
955  m_fields[0]->PhysDeriv(m_wk1[i*nvel+j],
956  outarray[i*nvel*nvel+j*nvel+0],
957  outarray[i*nvel*nvel+j*nvel+1]);
958  }
959  else
960  {
961  m_fields[0]->PhysDeriv(m_wk1[i*nvel+j],
962  outarray[i*nvel*nvel+j*nvel+0],
963  outarray[i*nvel*nvel+j*nvel+1],
964  outarray[i*nvel*nvel+j*nvel+2]);
965  }
966  }
967  }
968 
969  // Step 2: d(u^i_,j)/d(x^k) - {p,jk}*u^i_,p
970  for (int i=0; i< nvel; i++)
971  {
972  for (int p=0; p< nvel; p++)
973  {
974  tmp[p] = m_wk1[i*nvel+p];
975  }
977  for (int j=0; j< nvel; j++)
978  {
979  for (int k=0; k< nvel; k++)
980  {
981  Vmath::Vsub(physTot,outarray[i*nvel*nvel+j*nvel+k],1,
982  m_wk2[j*nvel+k],1,
983  outarray[i*nvel*nvel+j*nvel+k], 1);
984  }
985  }
986  }
987 
988  // Step 3: d(u^i_,j)/d(x^k) - {p,jk}*u^i_,p + {i,pk} u^p_,j
989  for (int j=0; j< nvel; j++)
990  {
991  for (int p=0; p< nvel; p++)
992  {
993  tmp[p] = m_wk1[p*nvel+j];
994  }
996  for (int i=0; i< nvel; i++)
997  {
998  for (int k=0; k< nvel; k++)
999  {
1000  Vmath::Vadd(physTot,outarray[i*nvel*nvel+j*nvel+k],1,
1001  m_wk2[i*nvel+k],1,
1002  outarray[i*nvel*nvel+j*nvel+k], 1);
1003  }
1004  }
1005  }
1006 
1007  // Restore value of wavespace
1008  m_fields[0]->SetWaveSpace(wavespace);
1009 }
1010 
1012  const Array<OneD, Array<OneD, NekDouble> > &inarray,
1013  Array<OneD, Array<OneD, NekDouble> > &outarray,
1014  const bool generalized)
1015 {
1016  int physTot = m_fields[0]->GetTotPoints();
1017  int nvel = m_nConvectiveFields;
1018 
1019  // Set wavespace to false and store current value
1020  bool wavespace = m_fields[0]->GetWaveSpace();
1021  m_fields[0]->SetWaveSpace(false);
1022 
1023  // For implicit treatment of viscous terms, we want the generalized
1024  // curlcurl and for explicit treatment, we want the cartesian one.
1025  if (generalized)
1026  {
1027  // Get the second derivatives u^i_{,jk}
1029  Array<OneD, Array<OneD, NekDouble> > ddU(nvel*nvel*nvel);
1030  gradgradU(inarray, ddU);
1031 
1032  // Raise index to obtain A^{ip}_{k} = g^pj u^i_{,jk}
1033  for (int i = 0; i < nvel; ++i)
1034  {
1035  for (int k = 0; k < nvel; ++k)
1036  {
1037  // Copy to wk
1038  for (int j = 0; j < nvel; ++j)
1039  {
1040  tmp[j] = ddU[i*nvel*nvel+j*nvel+k];
1041  }
1042  RaiseIndex(tmp, m_tmp);
1043  for (int p=0; p<nvel; ++p)
1044  {
1045  Vmath::Vcopy(physTot, m_tmp[p], 1,
1046  ddU[i*nvel*nvel+p*nvel+k], 1);
1047  }
1048  }
1049  }
1050  // The curlcurl is g^ji u^k_{kj} - g^jk u^i_kj = A^{ki}_k - A^{ik}_k
1051  for (int i = 0; i < nvel; ++i)
1052  {
1053  outarray[i] = Array<OneD, NekDouble> (physTot, 0.0);
1054  for (int k = 0; k < nvel; ++k)
1055  {
1056  Vmath::Vadd(physTot, outarray[i], 1,
1057  ddU[k*nvel*nvel+i*nvel+k], 1,
1058  outarray[i], 1);
1059  Vmath::Vsub(physTot, outarray[i], 1,
1060  ddU[i*nvel*nvel+k*nvel+k], 1,
1061  outarray[i], 1);
1062  }
1063  }
1064  }
1065  else
1066  {
1067  switch(nvel)
1068  {
1069  case 2:
1070  {
1071  Array<OneD,NekDouble> Vx(physTot);
1072  Array<OneD,NekDouble> Uy(physTot);
1073  Array<OneD,NekDouble> Dummy(physTot);
1074 
1075  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0],
1076  inarray[1], Vx);
1077  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],
1078  inarray[0], Uy);
1079 
1080  Vmath::Vsub(physTot, Vx, 1, Uy, 1, Dummy, 1);
1081 
1082  m_fields[0]->PhysDeriv(Dummy,outarray[1],outarray[0]);
1083 
1084  Vmath::Smul(physTot, -1.0, outarray[1], 1, outarray[1], 1);
1085  }
1086  break;
1087 
1088  case 3:
1089  {
1090  // Declare variables
1091  Array<OneD,NekDouble> Ux(physTot);
1092  Array<OneD,NekDouble> Uy(physTot);
1093  Array<OneD,NekDouble> Uz(physTot);
1094 
1095  Array<OneD,NekDouble> Vx(physTot);
1096  Array<OneD,NekDouble> Vy(physTot);
1097  Array<OneD,NekDouble> Vz(physTot);
1098 
1099  Array<OneD,NekDouble> Wx(physTot);
1100  Array<OneD,NekDouble> Wy(physTot);
1101  Array<OneD,NekDouble> Wz(physTot);
1102 
1103  Array<OneD,NekDouble> Dummy1(physTot);
1104  Array<OneD,NekDouble> Dummy2(physTot);
1105 
1106  // Calculate gradient
1107  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0],
1108  inarray[0], Ux);
1109  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],
1110  inarray[0], Uy);
1111  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],
1112  inarray[0], Uz);
1113 
1114  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0],
1115  inarray[1], Vx);
1116  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],
1117  inarray[1], Vy);
1118  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],
1119  inarray[1], Vz);
1120 
1121  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0],
1122  inarray[2], Wx);
1123  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],
1124  inarray[2], Wy);
1125  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],
1126  inarray[2], Wz);
1127 
1128  // x-component
1129  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],
1130  Vx, Dummy1); //Vxy
1131  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],
1132  Uy, Dummy2); //Uyy
1133  Vmath::Vsub(physTot, Dummy1, 1, Dummy2, 1, outarray[0], 1);
1134 
1135  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],
1136  Uz, Dummy1); //Uzz
1137  Vmath::Vsub(physTot, outarray[0], 1, Dummy1, 1,
1138  outarray[0], 1);
1139 
1140  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0],
1141  Wz, Dummy2); //Wxz
1142  Vmath::Vadd(physTot, outarray[0], 1, Dummy2, 1,
1143  outarray[0], 1);
1144 
1145  // y-component
1146  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],
1147  Wz, Dummy1); //Wzy
1148  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],
1149  Vz, Dummy2); //Vzz
1150  Vmath::Vsub(physTot, Dummy1, 1, Dummy2, 1, outarray[1], 1);
1151 
1152  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0],
1153  Vx, Dummy1); //Vxx
1154  Vmath::Vsub(physTot, outarray[1], 1, Dummy1, 1,
1155  outarray[1], 1);
1156 
1157  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0],
1158  Uy, Dummy2); //Uyx
1159  Vmath::Vadd(physTot, outarray[1], 1, Dummy2, 1,
1160  outarray[1], 1);
1161 
1162  // z-component
1163  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0],
1164  Uz, Dummy1); //Uxz
1165  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],
1166  Wy, Dummy2); //Wyy
1167  Vmath::Vsub(physTot, Dummy1, 1, Dummy2, 1, outarray[2], 1);
1168 
1169  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0],
1170  Wx, Dummy1); //Wxx
1171  Vmath::Vsub(physTot, outarray[2], 1, Dummy1, 1,
1172  outarray[2], 1);
1173 
1174  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],
1175  Vz, Dummy2); //Vyz
1176  Vmath::Vadd(physTot, outarray[2], 1, Dummy2, 1,
1177  outarray[2], 1);
1178  }
1179  break;
1180  }
1181  }
1182 
1183  // Restore value of wavespace
1184  m_fields[0]->SetWaveSpace(wavespace);
1185 }
1186 
1188 {
1189  int physTot = m_fields[0]->GetTotPoints();
1190  int nvel = m_nConvectiveFields;
1191  int nfields = m_fields.num_elements();
1192  int nbnds = m_fields[0]->GetBndConditions().num_elements();
1193 
1194  // Declare variables
1195  Array<OneD, int> BCtoElmtID;
1196  Array<OneD, int> BCtoTraceID;
1201 
1202  Array<OneD, NekDouble> ElmtVal(physTot, 0.0);
1203  Array<OneD, NekDouble> BndVal(physTot, 0.0);
1204  Array<OneD, NekDouble> coordVelElmt(physTot, 0.0);
1205  Array<OneD, NekDouble> coordVelBnd(physTot, 0.0);
1206  Array<OneD, NekDouble> Vals(physTot, 0.0);
1207 
1208  Array<OneD, bool> isDirichlet(nfields);
1209 
1210  Array<OneD, Array<OneD, NekDouble> > values(nfields);
1211  for (int i=0; i < nfields; i++)
1212  {
1213  values[i] = Array<OneD, NekDouble> (physTot, 0.0);
1214  }
1215 
1218  Array<OneD, Array<OneD, NekDouble> > coordVel(nvel);
1219  for (int i = 0; i< nvel; i++)
1220  {
1221  tmp[i] = Array<OneD, NekDouble> (physTot, 0.0);
1222  tmp2[i] = Array<OneD, NekDouble> (physTot, 0.0);
1223  coordVel[i] = Array<OneD, NekDouble> (physTot, 0.0);
1224  }
1225 
1226  // Get coordinates velocity in transformed system (for MovingBody regions)
1227  GetCoordVelocity(tmp);
1228  ContravarFromCartesian(tmp, coordVel);
1229 
1230  // Get Cartesian coordinates for evaluating boundary conditions
1232  for (int dir=0; dir < 3; dir++)
1233  {
1234  coords[dir] = Array<OneD, NekDouble> (physTot, 0.0);
1235  }
1236  GetCartesianCoordinates(coords[0],coords[1],coords[2]);
1237 
1238  // Loop boundary conditions looking for Dirichlet bc's
1239  for(int n = 0 ; n < nbnds ; ++n)
1240  {
1241  // Evaluate original Dirichlet boundary conditions in whole domain
1242  for (int i = 0; i < nfields; ++i)
1243  {
1244  BndConds = m_fields[i]->GetBndConditions();
1245  BndExp = m_fields[i]->GetBndCondExpansions();
1246  if ( BndConds[n]->GetBoundaryConditionType() ==
1248  {
1249  isDirichlet[i] = true;
1250  // If we have the a velocity component
1251  // check if all vel bc's are also Dirichlet
1252  if ( i<nvel )
1253  {
1254  for (int j = 0; j < nvel; ++j)
1255  {
1256  ASSERTL0(m_fields[j]->GetBndConditions()[n]->
1257  GetBoundaryConditionType() ==
1259  "Mapping only supported when all velocity components have the same type of boundary conditions");
1260  }
1261  }
1262  // Check if bc is time-dependent
1263  ASSERTL0( !BndConds[n]->IsTimeDependent(),
1264  "Time-dependent Dirichlet boundary conditions not supported with mapping yet.");
1265 
1266  // Get boundary condition
1267  LibUtilities::Equation condition =
1268  boost::static_pointer_cast<
1270  (BndConds[n])->
1271  m_dirichletCondition;
1272  // Evaluate
1273  condition.Evaluate(coords[0], coords[1], coords[2],
1274  time, values[i]);
1275  }
1276  else
1277  {
1278  isDirichlet[i] = false;
1279  }
1280  }
1281  // Convert velocity vector to transformed system
1282  if ( isDirichlet[0])
1283  {
1284  for (int i = 0; i < nvel; ++i)
1285  {
1286  Vmath::Vcopy(physTot, values[i], 1, tmp[i], 1);
1287  }
1288  ContravarFromCartesian(tmp, tmp2);
1289  for (int i = 0; i < nvel; ++i)
1290  {
1291  Vmath::Vcopy(physTot, tmp2[i], 1, values[i], 1);
1292  }
1293  }
1294 
1295  // Now, project result to boundary
1296  for (int i = 0; i < nfields; ++i)
1297  {
1298  BndConds = m_fields[i]->GetBndConditions();
1299  BndExp = m_fields[i]->GetBndCondExpansions();
1300 
1301  // Loop boundary conditions again to get correct
1302  // values for cnt
1303  int cnt = 0;
1304  for(int m = 0 ; m < nbnds; ++m)
1305  {
1306  int exp_size = BndExp[m]->GetExpSize();
1307  if (m==n && isDirichlet[i])
1308  {
1309  for (int j = 0; j < exp_size; ++j, cnt++)
1310  {
1311  m_fields[i]->GetBoundaryToElmtMap(BCtoElmtID,
1312  BCtoTraceID);
1313  /// Casting the bnd exp to the specific case
1314  Bc = boost::dynamic_pointer_cast<
1316  (BndExp[n]->GetExp(j));
1317  // Get element expansion
1318  elmt = m_fields[i]->GetExp(BCtoElmtID[cnt]);
1319  // Get values on the element
1320  ElmtVal = values[i] +
1321  m_fields[i]->GetPhys_Offset(
1322  BCtoElmtID[cnt]);
1323  // Get values on boundary
1324  elmt->GetTracePhysVals(BCtoTraceID[cnt],
1325  Bc, ElmtVal, BndVal);
1326 
1327  // Pointer to value that should be updated
1328  Vals = BndExp[n]->UpdatePhys()
1329  + BndExp[n]->GetPhys_Offset(j);
1330 
1331  // Copy result
1332  Vmath::Vcopy(Bc->GetTotPoints(),
1333  BndVal, 1, Vals, 1);
1334 
1335  // Apply MovingBody correction
1336  if ( (i<nvel) &&
1337  BndConds[n]->GetUserDefined() ==
1338  "MovingBody" )
1339  {
1340  // get coordVel in the element
1341  coordVelElmt = coordVel[i] +
1342  m_fields[i]->GetPhys_Offset(
1343  BCtoElmtID[cnt]);
1344 
1345  // Get values on boundary
1346  elmt->GetTracePhysVals(
1347  BCtoTraceID[cnt], Bc,
1348  coordVelElmt, coordVelBnd);
1349 
1350  // Apply correction
1351  Vmath::Vadd(Bc->GetTotPoints(),
1352  coordVelBnd, 1,
1353  Vals, 1, Vals, 1);
1354  }
1355  }
1356  }
1357  else // setting if m!=n
1358  {
1359  cnt += exp_size;
1360  }
1361  }
1362  }
1363  }
1364 
1365  // Finally, perform FwdTrans in all fields
1366  for (int i = 0; i < m_fields.num_elements(); ++i)
1367  {
1368  // Get boundary condition information
1369  BndConds = m_fields[i]->GetBndConditions();
1370  BndExp = m_fields[i]->GetBndCondExpansions();
1371  for(int n = 0 ; n < BndConds.num_elements(); ++n)
1372  {
1373  if ( BndConds[n]->GetBoundaryConditionType() ==
1375  {
1376  BndExp[n]->FwdTrans_BndConstrained(BndExp[n]->GetPhys(),
1377  BndExp[n]->UpdateCoeffs());
1378  if (m_fields[i]->GetExpType() == MultiRegions::e3DH1D)
1379  {
1380  BndExp[n]->HomogeneousFwdTrans(BndExp[n]->GetCoeffs(),
1381  BndExp[n]->UpdateCoeffs());
1382  }
1383  }
1384  }
1385  }
1386 }
1387 
1389  const NekDouble time,
1390  const Array<OneD, Array<OneD, NekDouble> > &coords ,
1391  const Array<OneD, Array<OneD, NekDouble> > &coordsVel)
1392 {
1393  if (m_fromFunction)
1394  {
1395  std::string s_FieldStr;
1396  string fieldNames[3] = {"x", "y", "z"};
1397  string velFieldNames[3] = {"vx", "vy", "vz"};
1398  // Check if function from session file defines each component
1399  // and evaluate them, otherwise there is no need to update
1400  // coords
1401  for(int i = 0; i < 3; i++)
1402  {
1403  s_FieldStr = fieldNames[i];
1404  if ( m_session->DefinesFunction(m_funcName, s_FieldStr))
1405  {
1406  EvaluateFunction(m_fields, m_session, s_FieldStr, m_coords[i],
1407  m_funcName, time);
1408  if ( i==2 && m_fields[0]->GetExpType() == MultiRegions::e3DH1D)
1409  {
1410  ASSERTL0 (false,
1411  "3DH1D does not support mapping in the z-direction.");
1412  }
1413  }
1414  s_FieldStr = velFieldNames[i];
1415  if ( m_session->DefinesFunction(m_velFuncName, s_FieldStr))
1416  {
1417  EvaluateFunction(m_fields, m_session, s_FieldStr,
1418  m_coordsVel[i], m_velFuncName, time);
1419  if ( i==2 && m_fields[0]->GetExpType() == MultiRegions::e3DH1D)
1420  {
1421  ASSERTL0 (false,
1422  "3DH1D does not support mapping in the z-direction.");
1423  }
1424  }
1425  }
1426  }
1427  else
1428  {
1429  int physTot = m_fields[0]->GetTotPoints();
1430  int nvel = m_nConvectiveFields;
1431  // Copy coordinates
1432  for(int i = 0; i < 3; i++)
1433  {
1434  Vmath::Vcopy(physTot, coords[i], 1, m_coords[i], 1);
1435  }
1436 
1437  for(int i = 0; i < nvel; i++)
1438  {
1439  Vmath::Vcopy(physTot, coordsVel[i], 1, m_coordsVel[i], 1);
1440  }
1441  }
1442 
1443  // Update the information required by the specific mapping
1444  UpdateGeomInfo();
1445 }
1446 
1447 }
1448 }
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:161
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:794
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:578
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
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
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:677
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:488
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:428
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:227
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:53
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:303
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:640
virtual GLOBAL_MAPPING_EXPORT void v_DotGradJacobian(const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, NekDouble > &outarray)
Definition: Mapping.cpp:689
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
The base class for all shapes.
Definition: StdExpansion.h:69
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:199
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:661
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:100
boost::shared_ptr< FieldIO > FieldIOSharedPtr
Definition: FieldIO.h:225
double NekDouble
virtual GLOBAL_MAPPING_EXPORT void v_CurlCurlField(const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const bool generalized)
Definition: Mapping.cpp:1011
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:901
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:1388
virtual GLOBAL_MAPPING_EXPORT void v_Divergence(const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, NekDouble > &outarray)
Definition: Mapping.cpp:765
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:329
MultiRegions::Direction const DirCartesianMap[]
Definition: ExpList.h:86
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::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:404
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:609
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:548
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:518
virtual GLOBAL_MAPPING_EXPORT void v_RaiseIndex(const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: Mapping.cpp:742
virtual GLOBAL_MAPPING_EXPORT void v_LowerIndex(const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: Mapping.cpp:719
GLOBAL_MAPPING_EXPORT void ReplaceField(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
Replace the Expansion List used by the mapping.
Definition: Mapping.cpp:243
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:359
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
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:1047
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:266
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:285
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:169
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:1187