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 namespace Nektar
39 {
40 namespace GlobalMapping
41 {
42 
44 bool Mapping::m_init = false;
45 bool Mapping::m_isDefined = false;
46 
48 {
49  typedef Loki::SingletonHolder<MappingFactory,
50  Loki::CreateUsingNew,
51  Loki::NoDestroy > Type;
52  return Type::Instance();
53 }
54 
57  : m_session(pSession), m_fields(pFields)
58 {
59  switch (m_fields[0]->GetExpType())
60  {
61  case MultiRegions::e1D:
62  {
64  }
65  break;
66 
67  case MultiRegions::e2D:
68  {
70  }
71  break;
72 
73  case MultiRegions::e3D:
76  {
78  }
79  break;
80 
81  default:
82  ASSERTL0(0,"Dimension not supported");
83  break;
84  }
85 
87  (pSession->GetComm());
88 }
89 
90 /**
91  * This function initialises the Mapping object. It computes the coordinates
92  * and velocity coordinates, initialises the workspace variables, and calls
93  * UpdateGeomInfo, which will perform the calculations specific for each type
94  * of Mapping.
95  * @param pFields ExpList array used in the mapping
96  * @param pMapping xml element describing the mapping
97  */
100  const TiXmlElement *pMapping)
101 {
102  int phystot = m_fields[0]->GetTotPoints();
103  m_fromFunction = true;
104  // Initialise variables
108  for (int i = 0; i < 3; i++)
109  {
110  m_coords[i] = Array<OneD, NekDouble> (phystot);
111  m_coordsVel[i] = Array<OneD, NekDouble> (phystot);
112  coords[i] = Array<OneD, NekDouble> (phystot);
113  }
114 
115  // Check if mapping is defined as time-dependent
116  const TiXmlElement* timeDep = pMapping->
117  FirstChildElement("TIMEDEPENDENT");
118  if (timeDep)
119  {
120  string sTimeDep = timeDep->GetText();
121  m_timeDependent = ( boost::iequals(sTimeDep,"true")) ||
122  ( boost::iequals(sTimeDep,"yes"));
123  }
124  else
125  {
126  m_timeDependent = false;
127  }
128 
129  // Load coordinates
130  string fieldNames[3] = {"x", "y", "z"};
131  const TiXmlElement* funcNameElmt = pMapping->FirstChildElement("COORDS");
132  if (funcNameElmt)
133  {
134  m_funcName = funcNameElmt->GetText();
135  ASSERTL0(m_session->DefinesFunction(m_funcName),
136  "Function '" + m_funcName + "' not defined.");
137 
138  // Get coordinates in the domain
139  m_fields[0]->GetCoords(coords[0], coords[1], coords[2]);
140 
141  std::string s_FieldStr;
142  // Check if function from session file defines each component
143  // and evaluate them, otherwise use trivial transformation
144  for(int i = 0; i < 3; i++)
145  {
146  s_FieldStr = fieldNames[i];
147  if ( m_session->DefinesFunction(m_funcName, s_FieldStr))
148  {
149  EvaluateFunction(m_fields, m_session, s_FieldStr, m_coords[i],
150  m_funcName);
151  if ( i==2 && m_fields[0]->GetExpType() == MultiRegions::e3DH1D)
152  {
153  ASSERTL0 (false,
154  "3DH1D does not support mapping in the z-direction.");
155  }
156  }
157  else
158  {
159  // This coordinate is not defined, so use (x^i)' = x^i
160  Vmath::Vcopy(phystot, coords[i], 1, m_coords[i], 1);
161  }
162  }
163  }
164  else
165  {
166  m_fields[0]->GetCoords(coords[0], coords[1], coords[2]);
167  for(int i = 0; i < 3; i++)
168  {
169  // Use (x^i)' = x^i as default. This can be useful if we
170  // have a time-dependent mapping, and then only the
171  // initial mapping will be trivial
172  Vmath::Vcopy(phystot, coords[i], 1, m_coords[i], 1);
173  }
174  }
175 
176  // Load coordinate velocity if they are defined,
177  // otherwise use zero to make it general
178  string velFieldNames[3] = {"vx", "vy", "vz"};
179  const TiXmlElement* velFuncNameElmt = pMapping->FirstChildElement("VEL");
180  if (velFuncNameElmt)
181  {
182  m_velFuncName = velFuncNameElmt->GetText();
183  ASSERTL0(m_session->DefinesFunction(m_velFuncName),
184  "Function '" + m_velFuncName + "' not defined.");
185 
186  std::string s_FieldStr;
187  // Check if function from session file defines each component
188  // and evaluate them, otherwise use 0
189  for(int i = 0; i < 3; i++)
190  {
191  s_FieldStr = velFieldNames[i];
192  if ( m_session->DefinesFunction(m_velFuncName, s_FieldStr))
193  {
194  EvaluateFunction(m_fields, m_session, s_FieldStr,
195  m_coordsVel[i], m_velFuncName);
196  if ( i==2 && m_fields[0]->GetExpType() == MultiRegions::e3DH1D)
197  {
198  ASSERTL0 (false,
199  "3DH1D does not support mapping in the z-direction.");
200  }
201  }
202  else
203  {
204  // This coordinate velocity is not defined, so use 0
205  Vmath::Zero(phystot, m_coordsVel[i], 1);
206  }
207  }
208  }
209  else
210  {
211  for(int i = 0; i < 3; i++)
212  {
213  Vmath::Zero(phystot, m_coordsVel[i], 1);
214  }
215  }
216 
217  // Initialise workspace variables
218  int nvel = m_nConvectiveFields;
222  for (int i=0; i< nvel; i++)
223  {
224  m_tmp[i] = Array<OneD, NekDouble>(phystot,0.0);
225  for (int j=0; j< nvel; j++)
226  {
227  m_wk1[i*nvel+j] = Array<OneD, NekDouble>(phystot,0.0);
228  m_wk2[i*nvel+j] = Array<OneD, NekDouble>(phystot,0.0);
229  }
230  }
231 
232  // Calculate information required by the particular mapping
233  UpdateGeomInfo();
234 }
235 
236 /**
237  * This function replaces the expansion list contained in m_fields, and then
238  * proceeds reinitialising the Mapping with this new field.
239  * @param pFields New field to be used by the Mapping
240  */
243 {
244  m_fields = pFields;
245 
246  TiXmlElement* vMapping = NULL;
247 
248  if (m_session->DefinesElement("Nektar/Mapping"))
249  {
250  vMapping = m_session->GetElement("Nektar/Mapping");
251  }
252  InitObject(pFields, vMapping);
253 }
254 
255 /**
256  * This function is responsible for loading the Mapping, guaranteeing that a
257  * single instance of this class exists. When it is first called, it creates a
258  * Mapping and returns a pointer to it. On subsequent calls, it just returns
259  * the pointer.
260  * @param pSession Session reader object
261  * @param pFields Fields which will be used by the Mapping
262  * @return Pointer to the Mapping
263  */
265  const LibUtilities::SessionReaderSharedPtr& pSession,
267 {
268  if (!m_init)
269  {
270  TiXmlElement* vMapping = NULL;
271  string vType;
272  if (pSession->DefinesElement("Nektar/Mapping"))
273  {
274  vMapping = pSession->GetElement("Nektar/Mapping");
275  vType = vMapping->Attribute("TYPE");
276  m_isDefined = true;
277  }
278  else
279  {
280  vType = "Translation";
281  }
282 
284  vType, pSession, pFields,
285  vMapping);
286 
287  m_init = true;
288  }
289 
290  return m_mappingPtr;
291 }
292 
293 /**
294  * This function should be called before writing a fld or chk file. It updates
295  * the metadata with information from the Mapping. Also, if the mapping is not
296  * defined by a function, it writes a .map file containing the coordinates and
297  * velocity of the coordinates.
298  * @param fieldMetaDataMap Metadata of the output file
299  * @param outname Name of the output file
300  */
302  LibUtilities::FieldMetaDataMap &fieldMetaDataMap,
303  const std::string &outname)
304 {
305  // Only do anything if mapping exists
306  if (m_isDefined)
307  {
308  fieldMetaDataMap["MappingCartesianVel"] = std::string("False");
309  if (m_fromFunction)
310  {
311  // Add metadata
312  fieldMetaDataMap["MappingType"] = std::string("Expression");
313  fieldMetaDataMap["MappingExpression"] = m_funcName;
314  if (m_timeDependent)
315  {
316  fieldMetaDataMap["MappingVelExpression"] = m_velFuncName;
317  }
318  }
319  else
320  {
321  int expdim = m_fields[0]->GetGraph()->GetMeshDimension();
322  string fieldNames[3] = {"x", "y", "z"};
323  string velFieldNames[3] = {"vx", "vy", "vz"};
324 
325  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
326  = m_fields[0]->GetFieldDefinitions();
327  std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
328 
329  int ncoeffs = m_fields[0]->GetNcoeffs();
330  Array<OneD, NekDouble> fieldcoeffs(ncoeffs);
331 
332  bool wavespace = m_fields[0]->GetWaveSpace();
333  m_fields[0]->SetWaveSpace(false);
334  // copy coordinates Data into FieldData and set variable
335  for(int j = 0; j < expdim; ++j)
336  {
337  m_fields[0]->FwdTrans_IterPerExp(m_coords[j], fieldcoeffs);
338 
339  for(int i = 0; i < FieldDef.size(); ++i)
340  {
341  // Could do a search here to find correct variable
342  FieldDef[i]->m_fields.push_back(fieldNames[j]);
343  m_fields[0]->AppendFieldData(FieldDef[i], FieldData[i],
344  fieldcoeffs);
345  }
346  }
347  if (m_timeDependent)
348  {
349  //copy coordinates velocity Data into FieldData and set variable
350  for(int j = 0; j < expdim; ++j)
351  {
352  m_fields[0]->FwdTrans_IterPerExp(m_coordsVel[j],
353  fieldcoeffs);
354 
355  for(int i = 0; i < FieldDef.size(); ++i)
356  {
357  // Could do a search here to find correct variable
358  FieldDef[i]->m_fields.push_back(velFieldNames[j]);
359  m_fields[0]->AppendFieldData(FieldDef[i],
360  FieldData[i],
361  fieldcoeffs);
362  }
363  }
364  }
365 
366  std::string outfile = outname;
367  outfile.erase(outfile.end()-4, outfile.end());
368  outfile += ".map";
369 
370  m_fld->Write(outfile,FieldDef,FieldData,fieldMetaDataMap);
371 
372  // Write metadata to orginal output
373  fieldMetaDataMap["MappingType"] = std::string("File");
374  fieldMetaDataMap["FileName"] = outfile;
375 
376  m_fields[0]->SetWaveSpace(wavespace);
377  }
378  }
379 }
380 
383  std::string pFieldName,
384  Array<OneD, NekDouble>& pArray,
385  const std::string& pFunctionName,
386  NekDouble pTime)
387 {
388  ASSERTL0(pSession->DefinesFunction(pFunctionName),
389  "Function '" + pFunctionName + "' does not exist.");
390 
392  pSession->GetFunction(pFunctionName, pFieldName);
393 
394  Array<OneD, NekDouble> x0(1,0.0);
395  Array<OneD, NekDouble> x1(1,0.0);
396  Array<OneD, NekDouble> x2(1,0.0);
397 
398  ffunc->Evaluate(x0, x1, x2, pTime, pArray);
399 }
400 
401 
405  std::string pFieldName,
406  Array<OneD, NekDouble>& pArray,
407  const std::string& pFunctionName,
408  NekDouble pTime)
409 {
410  ASSERTL0(pSession->DefinesFunction(pFunctionName),
411  "Function '" + pFunctionName + "' does not exist.");
412 
413  unsigned int nq = pFields[0]->GetNpoints();
414  if (pArray.num_elements() != nq)
415  {
416  pArray = Array<OneD, NekDouble> (nq);
417  }
418 
420  vType = pSession->GetFunctionType(pFunctionName, pFieldName);
422  {
423  Array<OneD, NekDouble> x0(nq);
424  Array<OneD, NekDouble> x1(nq);
425  Array<OneD, NekDouble> x2(nq);
426 
427  pFields[0]->GetCoords(x0, x1, x2);
429  pSession->GetFunction(pFunctionName, pFieldName);
430 
431  ffunc->Evaluate(x0, x1, x2, pTime, pArray);
432  }
433  else if (vType == LibUtilities::eFunctionTypeFile)
434  {
435  std::string filename = pSession->GetFunctionFilename(
436  pFunctionName,
437  pFieldName);
438 
439  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
440  std::vector<std::vector<NekDouble> > FieldData;
441  Array<OneD, NekDouble> vCoeffs(pFields[0]->GetNcoeffs());
442  Vmath::Zero(vCoeffs.num_elements(), vCoeffs, 1);
443 
446  (m_session->GetComm());
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  const 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  switch(nvel)
1066  {
1067  case 2:
1068  {
1069  Array<OneD,NekDouble> Vx(physTot);
1070  Array<OneD,NekDouble> Uy(physTot);
1071  Array<OneD,NekDouble> Dummy(physTot);
1072 
1073  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0],
1074  inarray[1], Vx);
1075  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],
1076  inarray[0], Uy);
1077 
1078  Vmath::Vsub(physTot, Vx, 1, Uy, 1, Dummy, 1);
1079 
1080  m_fields[0]->PhysDeriv(Dummy,outarray[1],outarray[0]);
1081 
1082  Vmath::Smul(physTot, -1.0, outarray[1], 1, outarray[1], 1);
1083  }
1084  break;
1085 
1086  case 3:
1087  {
1088  // Declare variables
1089  Array<OneD,NekDouble> Ux(physTot);
1090  Array<OneD,NekDouble> Uy(physTot);
1091  Array<OneD,NekDouble> Uz(physTot);
1092 
1093  Array<OneD,NekDouble> Vx(physTot);
1094  Array<OneD,NekDouble> Vy(physTot);
1095  Array<OneD,NekDouble> Vz(physTot);
1096 
1097  Array<OneD,NekDouble> Wx(physTot);
1098  Array<OneD,NekDouble> Wy(physTot);
1099  Array<OneD,NekDouble> Wz(physTot);
1100 
1101  Array<OneD,NekDouble> Dummy1(physTot);
1102  Array<OneD,NekDouble> Dummy2(physTot);
1103 
1104  // Calculate gradient
1105  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0],
1106  inarray[0], Ux);
1107  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],
1108  inarray[0], Uy);
1109  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],
1110  inarray[0], Uz);
1111 
1112  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0],
1113  inarray[1], Vx);
1114  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],
1115  inarray[1], Vy);
1116  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],
1117  inarray[1], Vz);
1118 
1119  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0],
1120  inarray[2], Wx);
1121  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],
1122  inarray[2], Wy);
1123  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],
1124  inarray[2], Wz);
1125 
1126  // x-component
1127  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],
1128  Vx, Dummy1); //Vxy
1129  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],
1130  Uy, Dummy2); //Uyy
1131  Vmath::Vsub(physTot, Dummy1, 1, Dummy2, 1, outarray[0], 1);
1132 
1133  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],
1134  Uz, Dummy1); //Uzz
1135  Vmath::Vsub(physTot, outarray[0], 1, Dummy1, 1,
1136  outarray[0], 1);
1137 
1138  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0],
1139  Wz, Dummy2); //Wxz
1140  Vmath::Vadd(physTot, outarray[0], 1, Dummy2, 1,
1141  outarray[0], 1);
1142 
1143  // y-component
1144  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],
1145  Wz, Dummy1); //Wzy
1146  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],
1147  Vz, Dummy2); //Vzz
1148  Vmath::Vsub(physTot, Dummy1, 1, Dummy2, 1, outarray[1], 1);
1149 
1150  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0],
1151  Vx, Dummy1); //Vxx
1152  Vmath::Vsub(physTot, outarray[1], 1, Dummy1, 1,
1153  outarray[1], 1);
1154 
1155  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0],
1156  Uy, Dummy2); //Uyx
1157  Vmath::Vadd(physTot, outarray[1], 1, Dummy2, 1,
1158  outarray[1], 1);
1159 
1160  // z-component
1161  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0],
1162  Uz, Dummy1); //Uxz
1163  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],
1164  Wy, Dummy2); //Wyy
1165  Vmath::Vsub(physTot, Dummy1, 1, Dummy2, 1, outarray[2], 1);
1166 
1167  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0],
1168  Wx, Dummy1); //Wxx
1169  Vmath::Vsub(physTot, outarray[2], 1, Dummy1, 1,
1170  outarray[2], 1);
1171 
1172  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],
1173  Vz, Dummy2); //Vyz
1174  Vmath::Vadd(physTot, outarray[2], 1, Dummy2, 1,
1175  outarray[2], 1);
1176  }
1177  break;
1178  }
1179  }
1180 
1181  // Restore value of wavespace
1182  m_fields[0]->SetWaveSpace(wavespace);
1183 }
1184 
1186 {
1187  int physTot = m_fields[0]->GetTotPoints();
1188  int nvel = m_nConvectiveFields;
1189  int nfields = m_fields.num_elements();
1190  int nbnds = m_fields[0]->GetBndConditions().num_elements();
1191 
1192  // Declare variables
1193  Array<OneD, int> BCtoElmtID;
1194  Array<OneD, int> BCtoTraceID;
1199 
1200  Array<OneD, NekDouble> ElmtVal(physTot, 0.0);
1201  Array<OneD, NekDouble> BndVal(physTot, 0.0);
1202  Array<OneD, NekDouble> coordVelElmt(physTot, 0.0);
1203  Array<OneD, NekDouble> coordVelBnd(physTot, 0.0);
1204  Array<OneD, NekDouble> Vals(physTot, 0.0);
1205 
1206  Array<OneD, bool> isDirichlet(nfields);
1207 
1208  Array<OneD, Array<OneD, NekDouble> > values(nfields);
1209  for (int i=0; i < nfields; i++)
1210  {
1211  values[i] = Array<OneD, NekDouble> (physTot, 0.0);
1212  }
1213 
1216  Array<OneD, Array<OneD, NekDouble> > coordVel(nvel);
1217  for (int i = 0; i< nvel; i++)
1218  {
1219  tmp[i] = Array<OneD, NekDouble> (physTot, 0.0);
1220  tmp2[i] = Array<OneD, NekDouble> (physTot, 0.0);
1221  coordVel[i] = Array<OneD, NekDouble> (physTot, 0.0);
1222  }
1223 
1224  // Get coordinates velocity in transformed system (for MovingBody regions)
1225  GetCoordVelocity(tmp);
1226  ContravarFromCartesian(tmp, coordVel);
1227 
1228  // Get Cartesian coordinates for evaluating boundary conditions
1230  for (int dir=0; dir < 3; dir++)
1231  {
1232  coords[dir] = Array<OneD, NekDouble> (physTot, 0.0);
1233  }
1234  GetCartesianCoordinates(coords[0],coords[1],coords[2]);
1235 
1236  // Loop boundary conditions looking for Dirichlet bc's
1237  for(int n = 0 ; n < nbnds ; ++n)
1238  {
1239  // Evaluate original Dirichlet boundary conditions in whole domain
1240  for (int i = 0; i < nfields; ++i)
1241  {
1242  BndConds = m_fields[i]->GetBndConditions();
1243  BndExp = m_fields[i]->GetBndCondExpansions();
1244  if ( BndConds[n]->GetBoundaryConditionType() ==
1246  {
1247  isDirichlet[i] = true;
1248  // If we have the a velocity component
1249  // check if all vel bc's are also Dirichlet
1250  if ( i<nvel )
1251  {
1252  for (int j = 0; j < nvel; ++j)
1253  {
1254  ASSERTL0(m_fields[j]->GetBndConditions()[n]->
1255  GetBoundaryConditionType() ==
1257  "Mapping only supported when all velocity components have the same type of boundary conditions");
1258  }
1259  }
1260  // Check if bc is time-dependent
1261  ASSERTL0( !BndConds[n]->IsTimeDependent(),
1262  "Time-dependent Dirichlet boundary conditions not supported with mapping yet.");
1263 
1264  // Get boundary condition
1265  LibUtilities::Equation condition =
1266  boost::static_pointer_cast<
1268  (BndConds[n])->
1269  m_dirichletCondition;
1270  // Evaluate
1271  condition.Evaluate(coords[0], coords[1], coords[2],
1272  time, values[i]);
1273  }
1274  else
1275  {
1276  isDirichlet[i] = false;
1277  }
1278  }
1279  // Convert velocity vector to transformed system
1280  if ( isDirichlet[0])
1281  {
1282  for (int i = 0; i < nvel; ++i)
1283  {
1284  Vmath::Vcopy(physTot, values[i], 1, tmp[i], 1);
1285  }
1286  ContravarFromCartesian(tmp, tmp2);
1287  for (int i = 0; i < nvel; ++i)
1288  {
1289  Vmath::Vcopy(physTot, tmp2[i], 1, values[i], 1);
1290  }
1291  }
1292 
1293  // Now, project result to boundary
1294  for (int i = 0; i < nfields; ++i)
1295  {
1296  BndConds = m_fields[i]->GetBndConditions();
1297  BndExp = m_fields[i]->GetBndCondExpansions();
1298 
1299  // Loop boundary conditions again to get correct
1300  // values for cnt
1301  int cnt = 0;
1302  for(int m = 0 ; m < nbnds; ++m)
1303  {
1304  int exp_size = BndExp[m]->GetExpSize();
1305  if (m==n && isDirichlet[i])
1306  {
1307  for (int j = 0; j < exp_size; ++j, cnt++)
1308  {
1309  m_fields[i]->GetBoundaryToElmtMap(BCtoElmtID,
1310  BCtoTraceID);
1311  /// Casting the bnd exp to the specific case
1312  Bc = boost::dynamic_pointer_cast<
1314  (BndExp[n]->GetExp(j));
1315  // Get element expansion
1316  elmt = m_fields[i]->GetExp(BCtoElmtID[cnt]);
1317  // Get values on the element
1318  ElmtVal = values[i] +
1319  m_fields[i]->GetPhys_Offset(
1320  BCtoElmtID[cnt]);
1321  // Get values on boundary
1322  elmt->GetTracePhysVals(BCtoTraceID[cnt],
1323  Bc, ElmtVal, BndVal);
1324 
1325  // Pointer to value that should be updated
1326  Vals = BndExp[n]->UpdatePhys()
1327  + BndExp[n]->GetPhys_Offset(j);
1328 
1329  // Copy result
1330  Vmath::Vcopy(Bc->GetTotPoints(),
1331  BndVal, 1, Vals, 1);
1332 
1333  // Apply MovingBody correction
1334  if ( (i<nvel) &&
1335  BndConds[n]->GetUserDefined() ==
1336  "MovingBody" )
1337  {
1338  // get coordVel in the element
1339  coordVelElmt = coordVel[i] +
1340  m_fields[i]->GetPhys_Offset(
1341  BCtoElmtID[cnt]);
1342 
1343  // Get values on boundary
1344  elmt->GetTracePhysVals(
1345  BCtoTraceID[cnt], Bc,
1346  coordVelElmt, coordVelBnd);
1347 
1348  // Apply correction
1349  Vmath::Vadd(Bc->GetTotPoints(),
1350  coordVelBnd, 1,
1351  Vals, 1, Vals, 1);
1352  }
1353  }
1354  }
1355  else // setting if m!=n
1356  {
1357  cnt += exp_size;
1358  }
1359  }
1360  }
1361  }
1362 
1363  // Finally, perform FwdTrans in all fields
1364  for (int i = 0; i < m_fields.num_elements(); ++i)
1365  {
1366  // Get boundary condition information
1367  BndConds = m_fields[i]->GetBndConditions();
1368  BndExp = m_fields[i]->GetBndCondExpansions();
1369  for(int n = 0 ; n < BndConds.num_elements(); ++n)
1370  {
1371  if ( BndConds[n]->GetBoundaryConditionType() ==
1373  {
1374  BndExp[n]->FwdTrans_BndConstrained(BndExp[n]->GetPhys(),
1375  BndExp[n]->UpdateCoeffs());
1376  if (m_fields[i]->GetExpType() == MultiRegions::e3DH1D)
1377  {
1378  BndExp[n]->HomogeneousFwdTrans(BndExp[n]->GetCoeffs(),
1379  BndExp[n]->UpdateCoeffs());
1380  }
1381  }
1382  }
1383  }
1384 }
1385 
1387  const NekDouble time,
1388  const Array<OneD, Array<OneD, NekDouble> > &coords ,
1389  const Array<OneD, Array<OneD, NekDouble> > &coordsVel)
1390 {
1391  if (m_fromFunction)
1392  {
1393  std::string s_FieldStr;
1394  string fieldNames[3] = {"x", "y", "z"};
1395  string velFieldNames[3] = {"vx", "vy", "vz"};
1396  // Check if function from session file defines each component
1397  // and evaluate them, otherwise there is no need to update
1398  // coords
1399  for(int i = 0; i < 3; i++)
1400  {
1401  s_FieldStr = fieldNames[i];
1402  if ( m_session->DefinesFunction(m_funcName, s_FieldStr))
1403  {
1404  EvaluateFunction(m_fields, m_session, s_FieldStr, m_coords[i],
1405  m_funcName, time);
1406  if ( i==2 && m_fields[0]->GetExpType() == MultiRegions::e3DH1D)
1407  {
1408  ASSERTL0 (false,
1409  "3DH1D does not support mapping in the z-direction.");
1410  }
1411  }
1412  s_FieldStr = velFieldNames[i];
1413  if ( m_session->DefinesFunction(m_velFuncName, s_FieldStr))
1414  {
1415  EvaluateFunction(m_fields, m_session, s_FieldStr,
1416  m_coordsVel[i], m_velFuncName, time);
1417  if ( i==2 && m_fields[0]->GetExpType() == MultiRegions::e3DH1D)
1418  {
1419  ASSERTL0 (false,
1420  "3DH1D does not support mapping in the z-direction.");
1421  }
1422  }
1423  }
1424  }
1425  else
1426  {
1427  int physTot = m_fields[0]->GetTotPoints();
1428  int nvel = m_nConvectiveFields;
1429  // Copy coordinates
1430  for(int i = 0; i < 3; i++)
1431  {
1432  Vmath::Vcopy(physTot, coords[i], 1, m_coords[i], 1);
1433  }
1434 
1435  for(int i = 0; i < nvel; i++)
1436  {
1437  Vmath::Vcopy(physTot, coordsVel[i], 1, m_coordsVel[i], 1);
1438  }
1439  }
1440 
1441  // Update the information required by the specific mapping
1442  UpdateGeomInfo();
1443 }
1444 
1445 }
1446 }
static MappingSharedPtr m_mappingPtr
Definition: Mapping.h:434
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
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:47
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
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: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:428
Array< OneD, Array< OneD, NekDouble > > m_wk2
Definition: Mapping.h:440
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:301
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
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
string m_velFuncName
Name of the function containing the velocity of the coordinates.
Definition: Mapping.h:422
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:98
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:1009
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:1386
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:329
MultiRegions::Direction const DirCartesianMap[]
Definition: ExpList.h:86
string m_funcName
Name of the function containing the coordinates.
Definition: Mapping.h:420
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:381
GLOBAL_MAPPING_EXPORT Mapping(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
Constructor.
Definition: Mapping.cpp:55
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:402
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:241
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:264
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:1185